青藏高原，是世界上海拔最高、面积大的高原，其独特的地形使得其成为全球气候变化的敏感区(Zeng et al., 2015a).许多的研究均表明，青藏高原通过大气和水文过程极大地影响着全球的能量和水循环(Yang et al., 2013).在青藏高原气候和水循环中，土壤水分扮演着重要角色，是陆气相互作用中的主要物理量之一.它通过影响陆地表面蒸散、水的运移和碳循环与地球气候系统相互作用，从而对气候系统及其变化产生重要影响(Yang et al., 2013; Su et al., 2013).因此，获取青藏高原地区高精度的大尺度范围的土壤水分信息对于研究全球的气候变化有着至关重要的作用.

传统的土壤水分测量方法，如重量法测量和时域反射仪法都是基于点的测量，该方法操作简单、精度高，但时间和人力成本高，且无法应用于大范围的地表土壤水分测量(Huisman et al., 2001; Noborio, 2001).随着遥感技术的发展，光学和微波遥感技术成为估测土壤水分的重要手段(陈书林等, 2012; 康悦等, 2014).光学遥感的温度植被干旱指数(Temperature Vegetation Dryness Index, TVDI)是研究土壤水分的重要方法(Sandholt et al., 2002; Patel et al., 2009).但在传统的地表温度(Surface Temperature, Ts)与NDVI构成的TVDI模型中，由于NDVI在植被生长后期会出现饱和现象的问题，因此，利用NDVI-Ts反演TVDI模型将产生一定的误差(Girolimetto and Venturini, 2013)；Gao等(2011)利用修正型土壤调整植被指数(Modified Soil Adjusted Vegetation Index, MSAVI)修正了TVDI模型，其经过修正的TVDI模型精度更高.同时，研究区内地形起伏对地表温度也将产生一定影响，使TVDI模型在估算大区域的土壤水分时反演精度会有所降低(冉琼等, 2005; Howell et al., 2010).赵跃环等(2015)利用数字高程模型(Digital Elevation Model, DEM)对鄂尔多斯盆地地表温度进行了订正，结果表明经过DEM校正获取的TVDI能更好地反映土壤水分.

微波遥感监测土壤水分的物理基础是土壤的介电特性与土壤水分之间存在较密切的关系.被动微波遥感在地表土壤水分的反演方面，以其较高的时间分辨率和较低的成本，成为区域土壤水分监测的重要手段(Liu et al., 2012; Mladenova et al., 2014).但是，由于被动微波传感器的空间分辨率通常较低，极大地限制了其产品的适用范围.如何联合光学数据及被动微波数据，发挥两者的优势，是当今研究的热点之一.赵杰鹏等(2012)结合TVDI与AMSR-E数据，通过复合反演模型反演了新疆地区的土壤水分，反演结果的精度相比于单一数据有一定的提高，相关系数由0.36提高到0.86.

基于以上分析，本文选取对气候变化最为敏感地区之一的青藏高原地区为研究区，建立一种基于光学遥感及微波遥感数据融合反演土壤水分的新方法.首先，基于使用MSAVI参数代替NDVI，以减少NDVI易饱和的缺点；同时利用Aster GDEM数据对地表温度进行修正；然后辅以AMSR-E/AMSR2与降雨数据，建立研究区的BP神经网络模型，并应用该模型反演了观测网地区的土壤水分分布，利用实测数据对土壤水分反演结果进行了验证.结果表明，该算法在精度上优于传统方法.与此同时，本文应用该算法得到了青藏高原地区连续的土壤水分分布图，并将土壤水分的月变化趋势与实测降水量月变化趋势比较，结果表明二者月变化量的正负关系一致，该算法具有良好的实用性.

本研究的AMSR-E亮温数据来自NSIDC (National Snow and Ice Data Center)的AMSR-E的三级产品升轨的土壤水分数据.由于AMSR-E传感器于2011年9月停止传输数据，本文使用AMSR-E传感器的第二代产品AMSR2传感器对2012年AMSR2那曲观测网的实测数据进行后续研究，使用的AMSR2数据产品为三级产品升轨的土壤水分数据.

用于反演TVDI的数据为从美国USGS下载的地表温度数据MOD11A2和陆面地表反射率数据MOD09Q1.MOD11A2为MODIS/Terra卫星8天合成的地表温度数据，空间分辨率为1 km (Chu et al., 2003).MOD09Q1数据为MODIS/Terra卫星8天合成的地表反射率数据，空间分辨率为250 m (de Beurs et al., 2008).MODIS数据处理主要包括数据格式转换、影像拼接、重采样、融合等几个步骤(陈顺云等, 2013; 杨婷等, 2014).

本研究使用与实测土壤湿度观测网数据时间同步的AMSR-E/AMSR2及MODIS数据进行后续研究(表 1).

表 1(Table 1)

表 1 遥感数据集及获取时间 Table 1 Remote sensing datasets and acquisition time

表 1 遥感数据集及获取时间 Table 1 Remote sensing datasets and acquisition time

青藏高原地形、植被变化存在显著的空间差异，因此，需考虑不同的气候及地表覆盖类型的实测数据.CAMP/Tibet、玛曲和那曲县为青藏高原地区典型研究区，具有不同的气候与地表覆盖类型(Lu et al., 2012; Chen et al., 2013; Zeng et al., 2015b)(图 1).本文采用部署在青藏高原的CAMP/Tibet、玛曲和那曲三个观测网土壤湿度的实测数据，为后续模型反演及验证本文反演的土壤湿度产品精度提供数据支持，数据详细介绍见表 2.考虑到目前光学遥感的穿透能力，本研究利用的是三个观测网所有站点表层4~5 cm的日观测土壤水分数据.

图 1

Fig. 1

图 1 观测网分布 Fig. 1 Distribution of in-situ networks

表 2(Table 2)

表 2 观测网数据及介绍 Table 2 Introduction of the dataset of the in-situ networks

表 2 观测网数据及介绍 Table 2 Introduction of the dataset of the in-situ networks

本研究首先对Ts-NDVI模型改进，利用MSAVI代替NDVI，与地表温度(Ts)建立TVDI模型；然后利用Aster GDEM数据获取地面高程，对地表温度(Ts)的差异进行校正；最后，辅以被动微波AMSR-E/AMSR2数据，以及降雨TRMM 3B43数据，建立神经网络模型以反演三个观测网(CAMP/Tibet、玛曲及那曲)的土壤水分，演生成流程见图 2.

图 2

Fig. 2

图 2 土壤水分反演流程图 Fig. 2 Process flow diagram of soil moisture estimation

目前国内外利用遥感数据进行TVDI模型反演中，NDVI是使用最广泛的植被指数.但是，研究表明，气温和降水的持续性和滞后性及植被覆盖度的影响可导致NDVI反演精度的降低；同时，当植被覆盖度较低时，NDVI对土壤背景的变化较为敏感(Girolimetto and Venturini, 2013; Eswar et al., 2016).为了研究精确反演植被覆盖较少的地区土壤水分的变化，必须利用消除了土壤背景影响的植被指数进行监测.MSAVI考虑了裸土土壤线，在描述植被覆盖度和土壤背景的方面比NDVI具有较高的优势(Gao et al., 2011).因此，本文选用MSAVI作为TVDI的植被指数，以消除土壤背景的影响，更准确地反映地表土壤水分.其计算公式为：

(1)

其中，DN_band1为MOD09Q1数据中可见光通道数据DN值，DN_band2为近红外通道数据DN值.

为比较Ts-MSAVI模型及Ts-NDVI模型的反演精度，本文选择样本点，分别对二者反演的TVDI值与对应的实测土壤水分数据的拟合程度进行对比.由于青藏高原地区7—8月降水量相对较多，因此，本文选取了三个观测网不同年份的7—8月的实测土壤水分数据与TVDI进行拟合，实验结果如图 3—5所示.从拟合图的相关系数(R²)中可以看出，相比于Ts-NDVI的模型，经Ts-MSAVI的模型得到的TVDI与实测土壤水分数据有着更好的拟合效果.对于CAMP/Tibet观测网，Ts-NDVI模型得到的TVDI与实测土壤水分数据拟合的R²为0.5236；经Ts-MSAVI模型得到的TVDI与实测土壤水分数据拟合的R²提高到0.5937；玛曲观测网中，R²从0.5013提高到0.6141；那曲观测网中，Ts-NDVI模型的R²为0.5253，而在Ts-MSAVI模型的R²提高到0.6013.

图 3

Fig. 3

图 3 CAMP/Tibet观测网TVDI与实测土壤水分拟合(2002年7—8月) (a) Ts-NDVI模型; (b) Ts-MSAVI模型. Fig. 3 Correlation analysis of the CAMP/Tibet in-situ measurements with TVDI (July and August, 2002) (a) Ts-NDVI model; (b) Ts-MSAVI model.

图 4

Fig. 4

图 4 玛曲观测网实测土壤水分与TVDI拟合效果比较(2008年7—8月) (a) Ts-NDVI模型; (b) Ts-MSAVI模型. Fig. 4 Correlation analysis of Maqu in-situ measurements with TVDI (July and August, 2008) (a) Ts-NDVI model; (b) Ts-MSAVI model.

图 5

Fig. 5

图 5 那曲观测网实测土壤水分与TVDI拟合效果比较(2010年7—8月) (a) Ts-NDVI模型; (b) Ts-MSAVI模型. Fig. 5 Correlation analysis of Nagqu in-situ measurements with TVDI (July and August, 2010) (a) Ts-NDVI model; (b) Ts-MSAVI model.

青藏高原地区平均海拔超过4500 m，最高海拔超过8800 m，海拔跨度较大.Sandholt等研究发现，气温、太阳辐射等因素都会对地温产生影响，而地面高程则是影响地表温度的重要因素之一(Sandholt et al., 2002; Garcia et al., 2014).因此，随着高程的变化，TVDI模型精度也受到一定影响.

冉琼等利用原始的未经地表温度修正的TVDI模型得到全国的TVDI分布图.研究发现，虽根据实测数据得到的青藏高原的土壤水分小于我国南部丘陵区，但反演结果显示青藏高原土壤水分明显高于南部丘陵区.这种不合理的分布是由于青藏高原地区高程远远高于南部地区，气温等因素对地表温度的降低作用显著.因此，需针对青藏高原地区高海拔的特性，将青藏高原地区的Ts修正到大地水准面的Ts值.本文利用Aster GDEM数据获取地面高程，对地表温度的差异进行校正，从而达到对大区域地表温度差异进行校正的目的.利用高程的校正公式如下：

(2)

式中，Ts₂为修正后的地表温度；Ts₁为修正前的地表温度；H为高程值；a为高程对土壤水分影响的程度.结合相关参考文献，本文采用的陆面温度随高程变化的系数a为0.6 ℃·100 m^-1.

此部分将基于文中3.1节的样本点，选择经过Ts-MSAVI的模型的构建的TVDI数据，进行Ts校正的实验.

如图 6所示，从拟合图的相关系数中可以看出，经过Ts校正的TVDI数据，与实测土壤水分数据有着更好的拟合效果.

图 6

Fig. 6

图 6 经过Ts校正的TVDI与实测土壤水分拟合效果比较 (a) CAMP/Tibet观测网；(b)玛曲观测网；(c)那曲观测网. Fig. 6 Correlation analysis of the in-situ measurements with TVDI (a) The CAMP/Tibet network; (b) The Maqu network; (c) The Nagqu network.

对于CAMP/Tibet观测网，经过Ts校正的TVDI与实测土壤水分数据拟合的R²由0.5937提高到0.6495；玛曲观测网中，R²从0.6141提高到0.7062；那曲观测网中，未经过Ts校正的R²为0.6013，而经过Ts校正的R²提高到0.6518.

土壤水分反映了在气候、气象和水文等因素共同作用下的当前地表状况，而气象数据提供了丰富的土壤状态信息.降雨对土壤渗透、径流和蒸发蒸腾有强烈的影响，是影响土壤水分的一个主要控制因子之一(Brocca et al., 2012).因此，本文从地气交互作用的普遍性出发，将降雨作为辅助信息，与土壤水分之间的作用关系进行研究.

本文使用的降雨数据来源于热带降雨观测(Tropical Rainfall Measurement Mission, TRMM)的测雨产品3B43，该数据集由TRMM全球3 h降雨估测数据3B42、NOAA气候预测中心气候异常监测系统的全球地面雨量计测量资料和全球降水气候中心的全球降水资料综合制作而成(Duan et al., 2013).

BP神经网络神经元网络算法具有极强的自组织性、非线性模拟的优点，非常适合解决非线性的复杂问题(余凡等, 2012).本文采用BP神经网络建立TVDI、AMSR-E/AMSR2及TRMM 3B43数据与实测土壤水分数据之间的关系模型，实现大区域连续土壤水分的遥感反演.输入层为TVDI、AMSR-E/AMSR2及TRMM 3B43数据三个要素；输出层为土壤水分(图 7).

图 7

Fig. 7

图 7 BP神经网络结构图 Fig. 7 The neural network used for soil moisture estimation

基于以上的算法和流程，本研究使用MODIS数据、AMSR-E/AMSR2数据、TRMM 3B43以及地面观测网实测数据，通过对原始的TVDI模型进行NDVI的改进、Ts校正等步骤，分别利用BP神经网络建立CAMP/Tibet、玛曲和那曲观测网的土壤水分反演模型.剔除无效值后，本文在CAMP/Tibet观测网选择了335个实测土壤水分样本点，玛曲观测网选取了804个样本点，那曲地区选取了2167个样本点.同时采用未参与建模的样本点对该BP神经网络模型进行结果拟合和验证(其中CAMP/Tibet地区93个，玛曲地区268个，那曲地区782个).

三个观测网经过模型预测的土壤水分和实测水分的拟合结果见图 8—10.如图所示，预测土壤水分与实测水分两者数据比较吻合，在数值上也较贴近.

图 8

Fig. 8

图 8 CAMP/Tibet观测网模型预测的土壤水分与实测土壤水分的对比 Fig. 8 Site-average and estimated soil moisture in CAMP/Tibet

图 9

Fig. 9

图 9 玛曲观测网模型预测的土壤水分与实测土壤水分的对比 Fig. 9 Site-average and estimated soil moisture in Maqu

图 10

Fig. 10

图 10 那曲观测网模型预测的土壤水分与实测土壤水分的对比 Fig. 10 Site-average and estimated soil moisture in Nagqu

对没有参与建模的样本点验证结果显示，CAMP/Tibet观测网模型平均误差(ME)为0.051 m³·m^-3，均方根误差(RMSE)为0.031 m³·m^-3；玛曲观测网模型ME为0.050 m³·m^-3，RMSE为0.035 m³·m^-3；那曲观测网模型ME为0.063 m³·m^-3，RMSE为0.041 m³m^-3；研究结果表明利用神经网络模型反演土壤水分结果相对可靠.相对于前人的研究成果，本研究中土壤水分的估算精度有了进一步提高.ME及RMSE的计算公式为

(3)

(4)

表 3(Table 3)

表 3 观测网实测数据验证结果 Table 3 Validation Results by the in-situ measurements

表 3 观测网实测数据验证结果 Table 3 Validation Results by the in-situ measurements

基于以上的算法，本研究得到了2010年青藏高原土壤水分的连续空间分布图(图 11).从图中可发现，青藏高原地降水随季节的变化明显，湿季降水较多，干季降水相对较少.其中，6—9月土壤水分较高，在8月和9月达到最高值；1—2月及11—12月的土壤水分为全年最低，这与气象资料的结果表现为一致的相关性.

图 1

Fig. 1

图 1 2010年青藏高原月平均土壤水分 Fig. 1 Monthly average soil moisture in the Tibetan Plateau

同时，东部和南部地表土壤湿度较高；藏北地区土壤水分相对较低.整体趋势为，土壤水分自东南向西北逐渐降低，这也符合青藏高原的不同地区气候区的空间分布(Yang et al., 2013; 刘强等, 2013).

因青藏高原地区地面土壤水分实测较为数据缺乏，实测站点对区域尺度范围的代表性较差，而土壤水分和降水的强相关性，因此，本文基于中国气象局提供的《中国地面气候资料月值数据集》的73个气象站点的月降水量，与土壤水分月变化率进行比较.

同时，为了减少基于点的站观测与基于像素的卫星观测在空间尺度上不匹配带来误差，本文采用聚类的方法，按照地理坐标，将高原的东西南北部分布的站点各聚集为一类，即将青藏高原地区的73个气象站分为四类，用每一类的月降水量均值代表这一区域的降水情况(图 12).

图 12

Fig. 12

图 12 站点聚类空间分布图 蓝色：高原北部站点聚类；绿色：高原南部站点聚类；黄色：高原西部站点聚类；红色：高原东部站点聚类. Fig. 12 Site points clustering spatial distribution map Blue: site clustering of the North plateau; Green: site clustering of the South plateau; Yellow: site clustering of the Western plateau; Red: site clustering of the Eastern plateau.

为了验证土壤水分和降水量的空间变化趋势，本文对高原月降水变化量空间分布和月土壤水分变化率空间分布情况进行对比(图 13).几乎所有的月份降水量月变化量与所在格网点的土壤水分月变化量的正负关系一致，表明本文产品的土壤水分2010年的月变化趋势基本上响应了该年降雨的月变化.

图 13

Fig. 13

图 13 四个聚类平均土壤水分月份差值与降水月份差值时间序列图对比(2010) (a)高原东部站点聚类序列图对比; (b)高原南部站点聚类序列图对比; (c)高原西部站点聚类序列图对比; (d)高原北部站点聚类序列图对比. Fig. 13 Comparison of the monthly difference between the average soil moisture and precipitation in 2010 (a) Difference of the cluster of the eastern part of the plateau; (b) Difference of the cluster of the south part of the plateau; (c) Difference of the cluster of the western part of the plateau; (d) Difference of the cluster of the northern part of the plateau.

本文提出了一种多特征变量融合的土壤水分反演新方法，以适应青藏高原地区地表覆盖类型、高程差异较大的特点，该方法首先针对TVDI模型在大范围土壤水分遥感反演中存在的不足，进行了两个方面的改进：以MSAVI代替NDVI进行TVDI的计算，以及对TVDI模型进行Ts修正.

同时，本文利用BP神经网络的方法，结合改进后的TVDI数据、AMSR-E/AMSR2数据、降雨数据以及地表实测土壤水分数据，建立融合多种变量的土壤水分反演模型.结合地面实测数据的验证结果表明，本文提出的基于多源遥感数据融合的土壤水分反演算法可以较高精度的反演土壤水分(CAMP/Tibet：RMSE=0.031 m³·m^-3；玛曲：RMSE=0.035 m³·m^-3；那曲：RMSE=0.041 m³·m^-3).相比于传统的方法，该算法即考虑了微波穿透性强的优势，又能发挥MODIS数据空间分辨率高的优点，更精确地反演较高空间与时间分辨率的土壤水分.

本文主要结论如下：

在TVDI模型中，NDVI指数对TVDI模型精度的影响较大，基于MSAVI-Ts构建的TVDI模型反演精度高于基于NDVI-Ts的TVDI模型.同时，由于地表高程变化引起的地表温度变化也会对TVDI模型精度产生影响.

提供了青藏高原月土壤水分月平均土壤水分的空间分布数据产品，并将其与降水量空间分布进行比较分析，二者在空间变化趋势上基本一致，因此，该算法在青藏高原地区具有良好的适用性.

联合光学MODIS数据与被动微波AMSR-E/AMSR2数据共同反演土壤水分是可行的方法，且反演精度较高，为陆气相互作用研究提供更为精确的数据.

本文融合了光学与被动微波数据，同时，也可与其他类型的遥感数据融合，比如主动微波数据或GNSS-R数据，取各自的优势，联合反演.如何融合的精度更高的遥感数据，是今后需要进一步研究的内容.

感谢中国科学院青藏高原研究所阳坤研究员提供实测数据，感谢美国家冰雪数据中心NSIDC提供MODIS，AMSE-R及AMSE2数据，感谢审稿专家提出宝贵意见.