1 引言

土壤湿度作为表示土层干湿程度的物理量，是水文、气候、生物、生态过程、环境、农业等领域的重要参数，是陆—气系统中一个关键的物理因子(郭维栋等，2007；Seneviratne et al., 2010)。它作为陆地和大气相互作用过程中的关键变量，在地表—大气间的水分和能量交换中起到调控和支配作用(卓嘎等，2015)。近年来，土壤湿度研究受到国内外学者广泛关注(Walker and Rowntree, 1977；De Lannoy et al., 2009；朱智和师春香，2014；沈润平等，2016；韩帅等，2017)。

土壤湿度时空变异较大(Moiwo et al., 2012)，地面观测只能获取站点的土壤湿度，在空间上难以达到连续。随着卫星遥感技术的发展与完善，已经研发了基于可见光遥感、热红外遥感和微波遥感监测土壤湿度的方法，使得获取时空连续高精度的土壤湿度数据成为可能(赖欣等，2014；席家驹等，2014)。土壤水分主要遥感反演方法有：热惯量法、植被指数法和微波遥感法，其中微波遥感法由于能直接捕获土壤的介电性以及能够全天候观测的能力被认为是目前及未来监测土壤湿度最具有潜力的手段，目前微波遥感土壤湿度产品主要有ECV CCI、AMSR-E、ASCAT及SMOS等(Schmugge, 1985；乔平林等，2007；王磊等，2009)。Dorigo et al.(2012, 2015)实验发现ECV CCI数据产品与观测数据的平均均方根误差在0.05 mm^{3 mm −3}左右，并利用ECV CCI、再分析资料和降水数据分析全球土壤湿度变化，ECV CCI的土壤湿度变化规律与其余数据集一致。Su et al. (2013)对AMSR-E、ASCAT及SMOS三种微波遥感土壤湿度产品在澳大利亚南部的质量表现进行了研究，发现ASCAT产品优于AMSR-E和SMOS产品。但微波遥感反演土壤湿度仍存在着空间分辨率较粗、易受积雪、高密度植被覆盖和地形起伏等因素的影响，且目前只能较好的反应表层土壤(5 cm)的湿度变化(Kerr, 2007；马思源等，2016)。利用陆面模式进行数值模拟是研究土壤湿度的另外一个比较有效的途径，通用陆面模式(CLM)模拟土壤湿度在中国区域研究已经得到了较为广泛的应用(李明星等，2011)。李明星等(2011)发现CLM3.5模式模拟的土壤湿度结果能较好地体现我国土壤湿度时空分布特征。李震坤等(2011)对CLM3.0的冻土过程参数化进行了改进，发现改进后的模式提高了青藏高原地区土壤水热过程模拟能力，土壤温度的模拟值更接近实测。朱晨等(2013)利用NCEP/NCAR的再分析数据作为驱动数据对中国区域不同深度的土壤湿度模拟和评估，发现0~10 cm和10~20 cm深度的模拟结果能够很好地表现出土壤湿度在时间和空间上的变化，在70~100 cm深度上的模拟结果空间分布有较好的一致性，而时间变率较小。熊建胜等(2014)利用CLM4.0对玛曲高原草原陆面过程进行了数值模拟，发现模式模拟结果能较好地表现出观测站土壤温度、湿度、地表辐射、湍流通量等的变化趋势，但土壤温度模拟偏低，感热通量模拟偏大。Fang et al.(2016)通过改进CLM4.5中土壤水力性质参数化方案来研究青藏高原地区土壤冻融过程，研究发现默认的参数化方案实验和改进的参数化方案实验均能够较好地模拟出土壤湿度的变化，其中改进参数化方案的实验对于高原地区土壤冻融循环的起始时间模拟结果更为合理。Wang and Zeng(2011)讨论了大气驱动场中降水和气温的质量对陆面模式模拟水循环变量的影响。CLDAS大气驱动数据是国家气象信息中心自主研制的高分辨率大气驱动数据，其研究组成员对CLDAS产品质量进行了验证与评估，并将其作为陆面模式的大气驱动场进行数值模拟，分析评估模拟结果(龚伟伟，2014；朱智和师春香，2014；朱智，2016；韩帅等，2017；孟现勇等，2017；孙帅等，2017)。以上研究表明土壤湿度的模拟准确性除了与模式参数化方案有关，还与大气驱动场的质量有关，存在不确定性。

青藏高原地区海拔高，环境条件恶劣，广泛分布着冰川和冻土下垫面，主要为多年冻土和季节性冻土，土壤湿度观测困难，遥感技术因空间分辨率、卫星过境时间、反演精度和积雪覆盖的限制，获取的土壤湿度数据有较大的局限性(郭东林和杨梅学，2010；万国宁等，2012；刘川等，2015；石磊等，2016)。因此，利用高时空分辨率的大气驱动数据开展高原地区土壤湿度模拟研究则尤为重要。本文利用最新的高时空精度(1 km、1 h)的中国陆面数据同化系统(HRCLDAS-V1.0)大气近地面强迫资料，驱动通用陆面模式(CLM)，开展青藏高原地区2015年1~9月土壤湿度模拟和分析评估研究，以期得到高分辨率高精度的土壤湿度模拟结果，为青藏高原地区土壤湿度研究提供资料支持。

2 研究数据与方法 2.1 研究区概况

青藏高原(25°N~40°N，70°E~105°E)平均海拔在4 km以上，是我国最大、世界海拔最高的高原，有“世界屋脊”和“第三极”之称(叶笃正和高由禧，1979；陈宇航等，2016)。由于海拔高，青藏高原的空气比较干燥、稀薄，太阳辐射年总量较高。青藏高原地域辽阔，高原东岸南部地区年平均气温可达到20 ℃，年降水量为2000 mm，植被种类十分丰富；而高原的西北部地区年平均气温在－6 ℃以下，年降水量只有50 mm。高原腹部地区年平均气温在0 ℃以下，广泛分布着冰川和冻土下垫面，主要为多年冻土和季节性冻土，其地表最显著的特征是土壤的冻结和消融(陈渤黎等，2012)。

2.2 研究数据 2.2.1 大气驱动数据

使用国家气象信息中心研制的最新HRCLDAS-V1.0(High Resolution CMA Land Data Assimilation System Version1.0)产品(http://10.1.64.154/idata/web/common/common/cldas?funItemId=124&level1NodeID=124010000&level3NodeID=124010101 [2017-09-26])，包括地面气压、2 m气温、2 m比湿、降水、10 m风速和太阳短波辐射等信息(龚伟伟，2014；韩帅等，2015；朱智，2016)。气温、气压、湿度和风速格点产品主要通过LAPS/STMAS系统(张涛等，2014)将ECMWF预报场数据与全国4万多个站点观测数据进行融合实现；短波辐射采用风云二号静止气象卫星可见光信息通道计算值作为背景场与使用Hybrid模型(Yang et al., 2001)计算得到的站点辐射值利用STMAS(the Space and Time Multiscale Analysis System)方法(Stamnes et al., 1988)进行融合获取；降水驱动数据则采用多卫星集成降水产品EMSIP(徐宾等，2015)作为背景场，与全国4万多个站点观测数据进行融合，获取0.01°分辨率降水融合产品。与之前的CLDAS系统数据相比，HRCLDAS-V1.0产品具有高时空分辨率特征，空间分辨率由6 km提高到1 km，时间分辨率为1 h，为土壤湿度高精度模拟提供了基础。

2.2.2 观测数据

由第三次青藏高原科学考察计划组提供观测数据，利用CS616土壤水分传感器每隔0.5 ms获取一次数据，由传感器通过其探头测得被测物体的介电常数，从而计算出含水量。本次试验利用青藏高原的中部地区6个站点(安多、那曲、班戈、聂荣、比如、嘉梨)的土壤湿度观测数据(图 1)。

2.3 研究方法

CLDAS (CMA Land Data Assimilation System Version)是国家气象信息中心研制的高分辨率多源土壤温湿度融合分析系统(师春香等，2013)。目前共有3个版本，CLDAS-V1.0产品主要包括气温、气压、比湿、风速、降水和太阳短波辐射等陆面驱动产品和土壤湿度产品，它利用数据同化技术，对地面观测、数值模式和卫星观测等产品进行融合，获得该数据集中的气压、温度、湿度、降水、风速、和太阳辐射等信息，土壤湿度采用耦合的CLM3.5陆面模式模拟制作。CLDAS-V2.0在CLDAS-V1.0的基础上耦合了多个陆面模式，包括CLM3.5、CoLM和Noah-MP三个模式。HRCLDAS-V1.0在CLDAS-V2.0和CLDAS-V1.0的基础上，引入高分辨率DEM和Albedo参数，结合FY卫星1 km通道信息研发高分辨率(1 km)大气强迫数据。同时，HRCLDAS-V1.0对模式地表土壤参数进行了优化，增加了产品实时质量检验评估等功能，及时检验评估系统产品质量。

CLM(Community Land Model)从CoLM(Common Land Model)的基础上发展而来，是由NACR发布建立的通用陆面模式，是目前国际上较为完善可靠的陆面模式之一(Dai et al., 2003；Oleson et al., 2004)。CLM3.5模式中空间异质性是由嵌套的次网格来描述的，每个网格划分为多个小网格，包括陆地单元、柱单元和植被功能类型。第一层为陆地单元，分为植被、湿地、冰川、湖泊和城市5类；第二层是柱单元，柱单元由10层土壤、地表覆盖层和最多5层积雪(积雪层数由雪深决定)来表示；第三层为植被功能类型，每个柱单元中可容纳4~17种不同植被(包括裸土)。相较于CLM3.0模式，CLM3.5版本采用了新的土壤冻融参数化方案(Niu and Yang, 2006)，新的土壤冻融参数化方案有以下几点改进：1)考虑到低于0 ℃时，土壤液态水和冰共存，并用零点抑制方程表示；2)在计算垂直水通量时引进一个部分渗透的概念，即一个网格被分为渗透部分和非渗透部分；3)将土壤液态水和冰作为土壤总体水分含量来计算土壤水势和渗透系数；其对高原地区土壤冻融过程模拟具有一定的优势。

CLM3.5通过一个10层模型来模拟土壤湿度的分布，土壤层为不均匀分布，第i层土层深度z_i计算公式如下：

$ {z_i} = \frac{{\exp \left( {\frac{{i - 0.5}}{2}} \right) - 1}}{4}{\rm{ }}, $

(1)

不透水率${f_{{\rm{frz}},i}}$计算公式(Oleson et al., 2007)为

$ {f_{{\rm{frz}},i}} = \exp \left[ { - \alpha \left( {1 - \frac{{{w_{{\rm{ice}},i}}}}{{{w_{{\rm{ice}},i}} + {w_{{\rm{liq}},i}}}}} \right)} \right] - {\rm{exp}}( - \alpha ){\rm{ }}, $

(2)

其中，可调整参数$\alpha = 3$，${w_{{\rm{ice}},i}}$和${w_{{\rm{liq}},i}}$分别为第i层土壤中冰和液态水的含量。

相邻两层土壤水渗透系数$k({z_{{\rm{h}},i}})$计算公式(Oleson et al., 2007)为

$ k({z_{{\rm{h}},i}}) = \left\{ \begin{array}{l} [1 - 0.5({f_{{\rm{frz}},i}} + {f_{{\rm{frz}},i + 1}})]{k_{{\rm{sat}}}}({z_{{\rm{h}},i}}) \cdot \\ {\left[ {\frac{{0.5({\theta _i} + {\theta _{i + 1}})}}{{0.5({\theta _{{\rm{sat}},i}} + {\theta _{{\rm{sat}},i + 1}})}}} \right]^{2{B_i} + 3}}{\rm{, }}1 \le i \le 9\\ \left( {1 - {f_{{\rm{frz}},i}}} \right){k_{{\rm{sat}}}}({z_{h,i}}){\left( {\frac{{{\theta _i}}}{{{\theta _{{\rm{sat}},i}}}}} \right)^{2{B_i} + 3}}{\rm{, }}i = 10 \end{array} \right. $

(3)

其中，${z_{{\rm{h}},i}}$为第i层土壤层厚度，${f_{{\rm{frz}},i}}$和${f_{{\rm{frz}},i + 1}}$分别为第i层和第i+1层土壤的不透水率，${k_{{\rm{sat}}}}({z_{{\rm{h}},i}})$为相邻两层土壤水饱和渗透系数，$\theta $为整体土壤体积含水量(即冰水混合)，${\theta _{{\rm{sat}}}}$为最大土壤体积含水量(即此时冰的含量达到最大)，B为可调整参数。

第i层土壤水势${\mathit{\Psi} _i}$计算公式为

$ {\mathit{\Psi} _i} = {\mathit{\Psi} _{{\rm{sat}},i}}{\left( {\frac{{{\theta _i}}}{{{\theta _{{\rm{sat}},i}}}}} \right)^{ - {B_i}}}, $

(4)

其中${\mathit{\Psi} _{{\rm{sat}},i}}$为第i层饱和土壤水势。

土壤水分通量根据Darcy定律(Darcy, 1856)计算得到：

$ q = - k\left( {\frac{{\delta (\mathit{\Psi} + z)}}{{\delta z}}} \right){\rm{ }}, $

(5)

其中，k为土壤水渗透系数，$\psi $为土壤水势，z为土层深度。

根据公式(3)，第i层土壤水分通量的计算公式(Oleson et al., 2004)为

$ {q_i} = - k({z_{{\rm{h}},i}})\left[ {\frac{{\left( {{\mathit{\Psi} _i} - {\mathit{\Psi} _{i + 1}}} \right) + ({z_{i + 1}} - {z_i})}}{{{z_{i + 1}} - {z_i}}}} \right]. $

(6)

本研究利用HRCLDAS-V1.0的大气驱动数据驱动CLM3.5陆面模式对青藏高原地区的土壤湿度进行模拟，得到的不同土层深度高时空分辨率(1 km、1 h)土壤湿度分布，考虑到研究区冬季大多被积雪覆盖，本次试验对2015年1月1日至9月30日的土壤湿度开展模拟与验证研究，试验每6小时输出一次模拟结果，即每日00:00(协调世界时，下同)、06:00、12:00和18:00四个时次输出土壤湿度值。为了与第三次青藏科学考察观测站点的土壤湿度测定的验证数据相匹配，本研究选取0~5 cm、0~10 cm和10~40 cm 3个土层深度的土壤湿度进行验证，并将他们分别称为第一层、第二层和第三层，将模式模拟的逐6小时各层土壤湿度双线性插值到站点，并与观测值进行对比分析。

土壤湿度模拟结果主要利用平均偏差(B_ias)、均方根误差(R_MSE)和相关系数(C_orr)进行评估，计算公式如下：

$ {B_{{\rm{ias}}}} = \frac{1}{N}\mathop \sum \limits_{i = 1}^N ({S_i} - {G_i}){\rm{ ,}} $

(7)

$ {R_{{\rm{MSE}}}} = \sqrt {\frac{{\mathop \sum \nolimits_{i = 1}^N {{({S_i} - {G_i})}^2}}}{N}} {\rm{ ,}} $

(8)

$ {C_{{\rm{orr}}}} = \frac{{\mathop \sum \nolimits_{i = 1}^N ({S_i} - \bar s)({G_i} - \bar G)}}{{\sqrt {\mathop \sum \nolimits_{i = 1}^N {{({S_i} - \bar s)}^2}} \sqrt {\mathop \sum \nolimits_{i = 1}^N {{({G_i} - \bar G)}^2}} }}{\rm{ }}, $

(9)

其中，N为样本数量，S_i为被模拟数据，G_i为观测站点数据，s和G分别为模拟数据和观测站点数据的平均值。

3 结果与分析 3.1 逐6小时土壤湿度模拟

逐6小时的土壤湿度变化能体现出更多的细节，从土壤湿度逐6小时变化(图 2)中可以看出，与观测值比较，安多、那曲、比如和聂荣等站点模式模拟的土壤湿度较好地表现出各层土壤湿度随时间的变化趋势。从不同站点表层土壤湿度测定值来看，安多站和那曲站第一层和第二层土壤湿度在0~600时次内变化较小，随后较快速上升和下降，整个研究时段变化显现两个峰值三个谷值态势；比如站和聂荣站点土壤湿度0~300时次变化平缓，随后土壤湿度变化波动显著加强，比如站整个时间变化显现两峰三谷，而聂荣站显现三峰四谷，模式模拟值基本上反映了变化规律，但是对于谷值普遍存在高估现象。从相同站点不同土层来看，与第一层和第二层相比，各个站点第三层(聂荣站第三层土壤湿度缺值)土壤湿度随时间变化波动明显变小，模式模拟值反映了这一特点，但存在高估现象。

从均值逐6小时时间序列(图 3)来看，模拟的土壤湿度能够较好地反映土壤湿度的逐6小时变化，第一层和第二层土壤湿度随着时间变化的波动较大，第三层波动较为平缓，这与土壤表层气象要素和积雪融化对较深层影响较少有关，因此，深层土壤湿度变化较平缓。第一层和第二层土壤湿度从300时次开始缓慢上升，650时次土壤湿度快速增加，600~1000时次内土壤湿度变化较为剧烈，模拟的土壤湿度能够很好地表现春季积雪融化，夏季降水增加和高温蒸发对土壤湿度的影响，但高估了第一层和第二层土壤湿度谷值，低估了峰值，这是由于模式中输入的叶面积指数是MODIS多年月平均叶面积指数，其空间分辨率较低，叶面积指数偏大，较多的植被阻挡了降水的下渗以及土壤水分的蒸发，模式模拟的蒸发量以及径流量较小，导致模式低估了土壤湿度峰值，而高估了土壤湿度的谷值，使得模拟值随时间的变化幅度小于观测值。模式对于第二层0~300时次土壤湿度略有低估，整体高估了第三层土壤湿度。从相关系数(图 4)来看，3个土层的模式模拟值与观测值的相关系数均在0.8以上，且均通过了0.01显著性水平检验，这表明模式能够较好地模拟出土壤湿度的逐6小时变化趋势，第一层相关系数达到0.88，随着土层的加深，相关系数略有下降。

3.2 土壤湿度月平均变化的模拟

从土壤湿度模拟值月变化(图 5)可以看出在1~9月，土壤湿度随着月份整体呈现上升的趋势，这与青藏高原地区夏季和秋季降水较多，冬季和春季土壤降水少的实际相符合。第一层和第二层土壤湿度在1~2月增长较少，3月土壤湿度开始显著增长，这是由于高原地区随着气温的升高，冰雪冻土开始消融，6~9月土壤湿度变化较为平缓，相对于前一个月份，7月土壤湿度略有下降，这可能与7月气温较高，土壤水分蒸发较强有关。第三层土壤湿度1~3月变化较缓，4月土壤湿度显著上升，5~9月土壤湿度基本没有变化。第一层土壤湿度模拟值与观测值的偏差在4~6月和9月小于0.02 mm³ mm^-3，7月和8月偏差在0.025~0.028 mm³ mm^-3之间，1~3月相对较大。第二层土壤湿度偏差在1月、4月、5月和9月较小，3月和8月偏差相对较大，达到0.04 mm³ mm^-3以上。第三层土壤湿度模拟值与观测值的偏差整体相对较大，在0.04 mm³ mm^-3左右，8月最大达到0.08 mm³ mm^-3。从相同月份不同土层来看，4~9月第三层土壤湿度模拟值与观测值的偏差较大，第二层次之，第一层最小；在1~2月，第一层土壤湿度偏差较大，第三层最小，3月各层土壤湿度的偏差最大(达到0.055 mm³ mm^-3)。

3.3 土壤湿度整体模拟结果差异分析

研究统计了2015年1~9月土壤湿度观测值与模拟值的整体相关系数、均方根误差及偏差，从表 1可以看出，第一层土壤湿度模拟值与观测值的相关系数最高，达到0.92；第二层次之；第三层土壤湿度模拟值与观测值的相关系数为0.83，相对较小，这表明在3个不同土层，模型对土壤湿度均具有较好的模拟效果。土壤湿度模拟值与观测值的均方根误差，第一层最小，第二层次之，第三层相对较大，但3个不同土层均小于0.045 mm³ mm^-3。土壤湿度模拟值与观测值的偏差统计表明，第二层最小，第一层次之，第三层偏差为0.039 mm³ mm^-3，相对较大。从整体上来看，第一层土壤湿度模拟结果最好，第二层次之，第三层模拟结果相对较差，说明模式在高原地区对于深层土壤湿度模拟效果有待改进。

3.4 土壤湿度模拟结果空间特征分析

从1~9月模拟的1 km分辨率土壤湿度平均值空间分布(图 6)来看，本地区土壤湿度呈现出北部低、南部较高、西边低、东边较高的分布特点，与高原地区年降水量从东南向西北逐渐减少实际相符合。高原地区植被覆盖度低、辐射强度大、蒸发强、降水量较少，深层土壤缺少水分补给，第三层土壤湿度整体上小于第一层和第二层土壤湿度。东南地区土壤湿度能够达到0.3~0.4 mm³ mm^-3，在西北地区三层土壤湿度均在0.1 mm³ mm^-3左右，土壤湿度较低，这与高原西北地区降水较少，地表覆盖多为荒漠和裸土有关(刘振元等，2017)。在中部地区有明显的土壤湿度较高的似圆形区域，该区域内部土壤湿度相对较大，外围土壤湿度低，这可能与高原中部有村落(如那曲、林芝等)存在以及高原中部以荒漠为主，草甸植被点缀为辅的地表覆盖有关。

4 小结与讨论

利用国家气象信息中心最新研发的HRCLDAS/CLM系统，模拟了2015年青藏高原地区不同深度土层1~9月土壤湿度变化，并利用第三次青藏高原科学考察的6个观测站点数据进行了验证，结果表明：

(1) 从站点逐6小时、月均和整个研究时段的土壤湿度模拟结果分析，模拟结果较好地反映了土壤湿度随时间的变化趋势和规律，即0~5 cm和0~10 cm土壤湿度变化趋势较为一致，变化波动较大，10~40 cm土壤湿度随着时间变化较为平缓。0~5 cm和0~10 cm土壤湿度从3月中旬开始加速上升，4~5月为次高峰值，6月和8月出现了土壤湿度的两个高峰值，7月为土壤湿度的低谷期；而10~40 cm土壤湿度从4月才开始上升，5~6月土壤湿度维持在0.1 mm³ mm^-3左右，7月土壤湿度有所下降，8月又开始缓慢上升。但是对于研究时段土壤湿度随时间的变化模拟值幅度略小于观测值，且模拟能力随着土层深度的加深而减弱。

(2) 从各层土壤湿度模拟结果来看，0~5 cm深度模式模拟的土壤湿度与观测值的相关系数最大，达到0.92，10~40 cm较小为0.83，均方根误差则相反，0~5 cm和0~10 cm土壤湿度模拟值与观测值的偏差均小于0.02 mm³ mm^-3，随着土层的加深，模式对于土壤湿度的模拟能力有所减弱。

(3) 从各层土壤湿度模拟结果空间分布来看，高原地区各层土壤湿度整体上从东南到西北逐渐变干，东南地区0~5 cm和0~10 cm土壤湿度高于10~40 cm土壤湿度，且随着土层的加深，较为湿润地区向东移动；西北地区三层土壤湿度均较低在0.1 mm³ mm^-3左右。

(4) 在前人(陈渤黎等，2012；陈晓磊等，2013；卓嘎等，2017)对青藏高原土壤湿度研究基础上，首次利用最新的高时空分辨率(1 km、1 h)的HRCLDAS-V1.0大气强迫数据驱动CLM模式，获得了更高空间分辨率(1 km)逐小时高原地区土壤湿度资料，更能揭示出该地区土壤湿度的时空变化的细节特征，揭示了青藏高原地区土壤湿度的时空分布特征，模拟的各层土壤湿度均能较好地反映出土壤湿度的实际变化规律，提高了土壤湿度的模拟精度，0~5 cm和0~10 cm深度土壤湿度模拟值与实测值的偏差减小到0.016 mm³ mm^-3以下。

(5) 本文利用观测站点6个站点9个月的观测资料进行分析，初步取得了较好的模拟结果，模式长期趋势的模拟能力还有待于更多站点观测长期资料的进一步验证。另外，陆面模式的模拟效果不仅与大气强迫驱动数据有关，还依赖于地表多种参数的精确度、模式的改进、模式物理过程参数化方案以及大量参数的设定，随着较高版本的陆面通用模式(如CLM4.0或CLM4.5)发展，应加强研究，以改进模拟的不足。