土壤水分是农业过程研究和环境因子评价的重要组成部分，是区域乃至全球尺度水循环研究中一个必不可少的参数，在水文、气候变化及土地退化研究等方面具有重要意义。传统的土壤水分监测主要依靠人工手动测量或站点监测，数据的精确度高，但是获取成本过高，无法满足长时间序列、大面积的土壤水分监测。随着遥感技术的发展，实现了大面积、动态获取土壤水分的目标。目前，土壤水分遥感监测方法主要有可见光-近红外遥感、高光谱遥感、热红外遥感、微波遥感等方法^[1]。其中，以光学遥感监测应用相对广泛和比较成熟，具体方法主要包括热惯量法、植被指数法、蒸散模型等^[2-4]。另外，微波遥感由于具有高时相、坚实的物理基础以及不受云干扰的优势，在全球土壤水分监测中也具有很好的应用前景^[5]。近年来，许多携带微波传感器的遥感卫星相继发射，进一步推动了微波遥感，尤其是被动微波遥感反演土壤水分的发展^[6-8]。我国第二代极轨气象卫星风云三号B星搭载的微波成像仪(microwave radiation imager，FY3B/MWRI)，在反演陆表参数中得到了较好的应用。王国杰等^[9]基于风云三号卫星微波亮温数据反演了土壤湿度资料。然而，不管是光学遥感还是微波遥感，都因存在自身的不足而在区域土壤水分反演中受到一定的制约^[10]。因此，不少学者尝试联合多源遥感数据进行较大区域范围的土壤水分估算，取得了较好的效果。Zhang等^[11]基于高级合成孔径雷达(advanced synthetic aperture radar，ASAR)数据和光学数据，以归一化植被指数(normalized difference vegetation index, NDVI)参数化植被冠层微波辐射传输特性，并结合水云模型构建了植被覆盖下地表土壤含水量的反演算法。马红章等^[12]基于SEMX02数据集，分析了L波段土壤发射率与地表土壤水分之间的关系，将NDVI引入土壤水分的反演算法中。张显峰等^[13]利用热惯量法和温度植被干旱指数法估算土壤水分基准值，利用高级微波扫描辐射计(advanced microwave scanning radiometer，AMSR-E)传感器的X波段反演土壤水分的日变化量，然后集二者建立土壤水分协同反演模型，得到较高时空分辨率的土壤水分结果。

藏北地区地处青藏高原腹地，海拔高，很多地方人迹罕至，土壤水分观测资料十分缺乏。遥感技术无疑是获取藏北地区土壤水分数据的最佳手段。然而，由于地表覆被类型复杂多样，常规的土壤水分遥感方法在藏北地区的监测效果不佳。一方面，植被光谱与土壤水分的关系易受植被类型和植被盖度的影响，单一的遥感模型不能很好适用于整个区域；另一方面，地理环境比较复杂，土壤水分具有较强烈的时空差异，遥感指标与土壤水分关系的不确定性增加。因此，本文拟在藏北地区植被覆盖类型划分的基础上，耦合适用于不同地表类型的光学遥感模型，并协同被动微波遥感反演表层土壤水分，以期提高研究区遥感反演土壤水分的精度，为藏北地区生态环境评价、水资源调查以及全球气候变化研究提供重要的土壤水分信息，同时也为其他陆表参数的定量遥感反演提供相应的参考。

1 研究区概况与研究方法 1.1 研究区概况

研究区范围包括西藏那曲地区和阿里地区，平均海拔4 500 m以上，地理位置大概在78° ~ 95°E，29°~ 36°N(图 1)。该区是西藏自治区重要的畜牧业生产基地，草地是区内主要的生态系统类型。研究区东西跨度大，地形非常复杂。东南部受构造作用、冰川作用及强烈的流水切割作用，为典型的高山峡谷地貌；中部高原地形保持完整，为高原宽谷区；西部地区为高原湖盆区。

区内属于典型的温带大陆性气候，年平均气温为-2.8 ~ 1.6 ℃，春冬多大风，降水主要集中在6—8月，年降水量247.3 ~ 513.6 mm。受大气环流和地形的影响，降水量总体呈现出由东向西、由东南向西北递减的趋势。此外，植被覆盖从东南向西北依次出现亚高山疏林、灌丛草甸、高寒草甸、高寒草原、高寒荒漠草原等基带^[14]。

1.2 数据来源

综合交通条件和地表覆被情况，本研究于2015年7月25日至8月6日在研究区范围内进行了土壤水分野外调查，共获得36个采样点数据(图 1)。调查中利用土壤水分速测仪TDR200获取10 cm深度的土壤水分数据，同时利用GPS记录采样点位置信息。考虑到光学遥感影像的空间分辨率，实地采样时在1 km × 1 km的范围内大致按对角线每隔200 ~ 300 m距离进行土壤水分值的测定，依据拉依达准则^[15]剔除其中粗大误差的异常值，最后对剩余数据求取平均值得到样点土壤含水量。选取26个样点数据作为建模样本，剩余10个数据进行模型验证和精度分析。

微波遥感数据选取2015年7月20日至8月6日的风云三号B星MWRI数据。风云三号MWRI可全天候、全天时地监测云中液态水含量、大气可降水量、地面降水量等。MWRI在10.65 ~ 89 GHz频段内设有水平和垂直两种极化方式的5个微波频点，共有10个探测频道。其中，低频10.65 GHz通道具有穿透云雨大气的能力，并对地表粗糙度和介电常数比较敏感，可以用于全天候获取风速、土壤水分含量等地球物理参数。研究选用FY3B/MWRI的一级(L1)产品降轨资料(10.65 GHz通道)，数据来源于国家卫星气象中心(http://www.nsmc.org.cn)。获取的对地观测亮温数据首先经地理位置查找表法(geographic lookup table, GLT)进行几何校正，然后利用ENVI5.0的band math工具进行辐射定标，最后经裁剪后得到研究区MWRI数据。

光学遥感数据选用NASA的Land Processes DAAC数据中心(lpdaac.usgs.gov/data access/data pool)的MODIS数据，包括Aqua卫星(PM 14:30过境)的8 d合成地表温度产品(MYD11 A2)、地表反照率产品(MCD43 B4)及16 d合成植被指数产品(MYD13 A2)，空间分辨率为1 km。以上产品数据首先经美国地质调查局(USGS)提供的MODIS数据处理工具(MODIS Reprojection Tool, MRT)进行影像拼接和投影转换处理。

其他辅助数据包括空间分辨率为30 m的全球数字高程数据(GDEMV2)，数据来源于地理空间数据云(http://www.gscloud.cn)，经重采样生成空间分辨率为1 km的数字高程数据(digital elevation model，DEM)，用于地表温度校正。此外，还有由西藏气象局提供的2015年7月和8月藏北地区15个气象监测站点的日累积降水量和气温数据，用于对土壤水分估算结果的精度分析。

1.3 研究方法 1.3.1 光学遥感和被动微波遥感协同反演地表土壤水分的方法

根据美国学者生产土壤水分产品使用的半经验回归模型^[16]，Zhang等^[17]对该模型作了进一步的改进，使之更适用于干旱区表层土壤水分的遥感监测与估算，改进的土壤水分遥感估算半经验回归模型可以表示为：

$ s{{m}_{ij}}\left( t \right)={{m}_{ij}}+\Delta {{m}_{ij}}\left( t \right) $

(1)

式中：m_ij为观测周期T内遥感影像的第i行第j列像元的土壤水分基准；Δm_ij(t)为t时刻第i行第j列土壤水分与基准m_ij的差异变化量；sm_ij(t)为该像元通过遥感反演得到的土壤水分。公式(1)表示假定在不长的一个观测时段T内，表层土壤水分变化可以分解成基准值和日变化量两部分。其中，基准值(m_ij)反映这段时期内土壤水分的一个最低状况，与区域自然条件和气候相关；变化量Δm_ij(t)反映了气象因子如降水、蒸散等引起的土壤水分短周期波动。正常情况下，在较短的时间范围内植被和地表粗糙度没有显著变化，那么，估算表层土壤水分变化量时只需消除地表温度的差异即可。因此，利用光学遥感数据可获得土壤水分基准值的反演，利用MWRI数据可获得土壤水分日变化量，两者协同能够实现表层土壤水分估算。

1.3.2 光学遥感反演地表土壤水分基准的方法

1) 裸地或低植被覆盖下的热惯量法。热惯量法具有较好的重复性、准确性和同一性，可以通过土壤水分含量与其热惯量间的良好线性关系直接反演土壤水分，但是只适合于裸地或低植被覆盖区。同时，由于求解真实热惯量所需参数比较复杂难解，实际应用中常用表观热惯量代替真实热惯量建立表观热惯量与土壤水分之间的遥感统计模型^[18]。表观热惯量的计算公式为：

$ \text{ATI}=\left( 1-\alpha \right)/\Delta T $

(2)

式中：ATI为表观热惯量；ΔT为昼夜温差；α为全波段反照率。α采用Liang^[19]针对MODIS数据提出的公式计算：

$ \begin{align} & \alpha ={{1.160}_{\rho }}_{1}+{{0.291}_{\rho }}_{2}+{{0.243}_{\rho }}_{3}+{{0.116}_{\rho }}_{4}+{{0.112}_{\rho }}_{5}+ \\ & {{0.081}_{\rho }}_{7}-0.0015 \\ \end{align} $

(3)

式中：ρ_i(i=1，2，3，4，5，7)为MODIS产品的各波段地物反射率。

2) 中等植被覆盖下的温度植被干旱指数(tem perature vegetation dryness Index, TVDI)模型。TVDI是利用Ts-NDVI特征空间提取的水分胁迫指标来估算陆面表层土壤水分的一种方法，可以表示为^[20]：

$ \text{TVDI=}\frac{{{T}_{\text{s}}}-{{T}_{\min }}}{{{T}_{\max }}-{{T}_{\min }}}=\frac{{{T}_{\text{s}}}-({{a}_{1}}+{{b}_{1}}\text{NDVI})}{({{a}_{2}}+{{b}_{2}}\text{NDVI})-({{a}_{1}}+{{b}_{1}}\text{NDVI})} $

(4)

式中：T_s和NDVI分别为给定像元的地表温度和归一化植被指数。a₁、b₁和a₂、b₂分别是湿边(T_min)和干边(T_max)的拟合系数。由此看出，TVDI的计算由T_s和NDVI决定，但稀疏植被区的T_s和NDVI受到土壤和植被的影响变得更加复杂，直接影响参数的估计。因此，TVDI模型比较适合用于中等植被覆盖区域。研究采用MODIS地表温度产品进行云掩膜和16 d平均值合成(利用连续两期的MODIS 8 d合成地表温度产品进行均值合成得到)以降低云对表层土壤水分反演的影响，针对地势起伏引起的地表温度差异问题则通过地面高程进行校正^[21]，利用重采样成1 km DEM数据对其进行高程校正，最后结合NDVI产品数据构建TVDI模型。

3) 密集植被覆盖下的植被供水指数(vegetation supply water index, VSWI)模型。在密集植被覆盖区域，NDVI值越大，TVDI的不确定性增大，因此在较高的植被覆盖区监测土壤水分不宜采用TVDI法。VSWI适合于高植被覆盖情况下的土壤水分监测，该方法获取资料容易、时效性强、物理意义明确，是进行大范围旱情监测的有效和常用方法^[22]。VSWI值越大，表明植被蒸腾旺盛，土壤水分含量较高，反之表明植被供水不足，土壤水分含量较低。同时，本文估算土壤水分基准值的前提假设是较短时期内植被和地表粗糙度基本维持不变，因此，密集植被覆盖下的土壤水分采用植被供水指数法进行估算^[23]。以MODIS 31通道的亮温数据反映植被冠层温度，计算NDVI与之植被冠层温度的比值得到VSWI。

1.3.3 被动微波遥感估算地表土壤水分日变化量的方法

一般情况下，亮温极化率P_r是指某一频率亮温的水平、垂直极化差与两者之和的比值，用其可以有效降低大气与地表温度对土壤水分的影响。因此，利用微波P_r和半经验模型可计算土壤水分日变化量^[24]：

$ \Delta {{m}_{v}}={{m}_{1}}({{P}_{r}}-{{P}_{r}}_{_{\min }}){{P}_{r}}_{_{\min }}^{{{m}_{2}}} $

(5)

式中：Δm_v为土壤水分日变化量；m₁、m₂为通过地面实测值得到的经验系数；P_r为亮温极化率；${{P}_{r}}_{_{\min }}$为10.65 GHz通道计算的观测周期内的P_r最小值。从公式(5)可以看出，在一定的时间尺度内，土壤水分的变化量主要由${{P}_{r}}_{_{\min }}$与${{P}_{r}}-{{P}_{r}}_{_{\min }}$共同决定，在较短时间内植被含水量和地表粗糙度的变化相对较小，降雨、蒸发等气象条件的改变会导致土壤水分发生较大变化。

2 结果与分析 2.1 土壤水分基准值反演

参考相关文献^[25-26]，依据2015年第209天到217天16 d合成的NDVI数据将研究区地表类型划分为3类：裸地或低植被覆盖区(0<NDVI<0.2)、中等植被覆盖区(0.2≤NDVI<0.4)和密集植被覆盖区(NDVI≥0.4) (图 2)。分别采用本文1.3中所介绍的模型方法计算相应地表类型区的ATI、TVDI和VSWI值，利用Origin 9.1软件，建立ATI、TVDI和VSWI值与实测土壤水分的拟合模型(表 1)。

由表 1可知，在裸地或低植被覆盖和密集植被覆盖下，ATI值、VSWI值与实测土壤水分拟合的幂函数模型效果均优于线性模型的结果，决定系数R²分别为0.62、0.85，而中等植被覆盖的TVDI值与地面实测土壤水分存在良好的一元线性相关，决定系数R²为0.63。3种地表类型区域的拟合模型的相关性均通过了置信度为0.05的显著性检验，说明依据地表类型划分的结果选择相应的模型对土壤水分进行估算可以取得较好的效果。

综合表 1中的拟合方程，估算得到2015年7月21日至8月5日16 d合成的土壤水分基准值(图 3)。由图 3可以看出，藏北地区土壤水分干湿状况存在明显的空间分异：东南、西北部地区土壤水分含量较高，冰雪覆盖的较高海拔的高山区土壤水分也较充足，土壤含水量在10%以上；中部地区的尼玛县大部分地区、改则县南部地区及班戈县等地，土壤水分含量较低，土壤水分含量在0 ~ 10%，部分区域土壤处于缺水状态，土壤含水量低于5%。

利用地面实测数据对3种遥感模型联合反演的土壤水分基准值进行精度分析(图 4)，可以看出，土壤水分基准值的估算结果与地面实测数据的决定系数R²为0.65，均方根误差RMSE为1.19，说明依据地表覆盖类型构建的光学遥感估算模型在一定程度上可以实现对土壤水分的估算，但估算精度仍然不高，主要原因一方面可能是NDVI最大值法虽然对去云是有效的，但却忽略了地表方向性反射，造成由NDVI反演的地表参数的误差；另一方面，由于采样沿线地表特征差异较大，通过TDR测量获得的真实土壤水分数据本身也有一定的误差，导致土壤水分基准值估算结果的整体精度偏低，最高90.21%，最低54.73%，平均为74.69%，反映出复合模型反演的土壤水分基准值误差的波动性较大。另外，实测数据是10 cm的土壤水分含量，而卫星遥感反演结果主要反映表层土壤水分，这也会对估算结果产生影响。因此，还需要多源遥感数据协同反演提高土壤水分精度。

2.2 基于微波遥感的土壤水分变化值估算

选择2015年7月29日和8月5两天的FY3B/ MWRI辐射亮温数据，采用公式(5)计算土壤水分日变化量，运用卷积算法重采样生成空间分辨率为1 km的土壤水分日变化量结果图(图 5)。由图 5可以看出，两期的土壤水分日变化量估算结果存在较大差异。8月5日土壤水分日变化量(图 5B)明显高于7月29日(图 5A)，平均变化量为18.5%，最小仅为0.16%，最大达47.88%。土壤水分日变化量较高的区域主要分布在纳木错湖、昆仑山南坡、双湖县、班戈县、改则和尼玛县北部，土壤水分变化量多在20%以上。结合气象资料分析，土壤水分日变化量较高的主要原因可能与当日研究区大部分地区有不同强度的降雨有关，导致土壤水分变化量明显增大。研究区中部土壤水分日变化量的空间差异较明显，这可能受降水时空分布不均、陆表蒸散发量较大等因素的影响，气温升高引起地表蒸发加快，导致土壤水分变化量减小^[27]。另外，土壤水分日变化量在5%以下的地区主要分布在西部、东部，且这些区域土壤水分变化量的值比较稳定。分析原因，主要是这部分地区植被覆盖度较高，植被对土壤水分的调节作用使得土壤水分变化波动较小^[28]。

2.3 光学和微波遥感协同反演土壤水分

根据公式(1)，联合土壤水分基准和日变化量得到藏北地区表层土壤水分结果(图 6)。总体上看，整个研究区土壤水分空间分异明显，在高山区、林区以及湖泊周边地区土壤水分含量较高，这与高山积雪融化补给以及降水量高有关。同时，结合气象条件分析，气温和降水引起土壤水分产生时空分布变化。研究区2015年7月29日天气状况为晴天，而8月5日在安多县北部、双湖县、改则县中部地区、尼玛县中部、革吉县东部等地均有降雨过程，因而8月5日当天土壤水分含量明显增加，尤其对于裸土及稀疏植被覆盖区域，土壤水分基准值偏低，降水过程使该区域的表层土壤水分迅速增加^[29]，因而土壤水分变化量较大。在中等植被覆盖区或密集植被覆盖区，一定强度的降水易引起土壤水分出现饱和而使其随地表径流流走，因而土壤水分的变化量相对较小。

利用预留的10个实测土壤含水量数据验证协同反演模型的精度与合理性(图 7)。由图 7可知，协同反演值与实测值表现出显著的线性相关，决定系数R²为0.89，均方根误差RMSE为0.97。结合图 4的结果可以看出，相比于利用单一的光学遥感方法，微波与光学遥感协同反演的方法明显提高了对土壤水分含量的估算精度。进一步利用2015年7月25日至8月6日藏北地区15个气象站点记录的地面实测降水数据对土壤水分协同反演估算结果的合理性进行分析(图 8)。由图 8可以看出，降水量与土壤水分含量估算结果之间有较好的对应关系，两者的变化趋势基本一致，相关系数达到0.79，这也反映出估算结果的可信性。

综上所述，针对藏北地区土壤水分实测数据获取困难，结合MODIS数据的适中空间分辨率特性和微波数据对土壤水分定量估算的良好效果，协同光学与微波遥感反演土壤水分的方法简单易行且效果较好。由于在利用光学和微波遥感协同反演土壤水分时考虑了植被覆盖类型、地形对估算结果的影响，加之仅对比了较短时期土壤水分估算结果的动态变化，因此，在本研究中地表参数条件对土壤水分估算结果的影响较小。而实际上，环境条件的变化，比如空气湿度、风速、局部地形差异、土壤质地等，都会引起估算结果产生偏差。因此，下一步工作需要考虑更多的影响因素对藏北地区土壤水分遥感估算结果的影响。另外，受实测数据的时间限制，对降水引起的土壤水分变化，本研究仅仅分析了降水当日的土壤水分变化量，土壤水分随降水过程的变化情况还有待更深入的研究。

3 结论

1) 利用NDVI阈值法划分地表植被覆盖类型，并在不同覆盖类型区选用适用的光学土壤水分估算模型，可以在一定程度上提高地表类型复杂区域的光学遥感估算土壤水分的精度。

2) 区分不同的地表植被覆盖类型，协同光学遥感与微波遥感，可以取得更好的土壤水分估算结果。本研究光学与微波遥感协同反演的土壤水分含量与实测值之间具有显著的线性相关性，决定系数R²为0.89，均方根误差RMSE为0.97，表明在藏北地区利用该方法进行土壤水分的反演是可行的。