土壤在地球表面与大气圈、水圈、岩石圈及生物圈相交并与其进行着物质循环和能量交换，是多种物理、化学和生物过程发生的场所^[1]。土壤水分作为土壤的重要组成部分，极大地影响着土壤性质，并且作为重要的溶剂，充当了土壤中各种过程的水环境。干旱区水资源匮乏，土壤水分含量是影响干旱区水分和能量平衡的关键变量，是区域水文和植被过程中最为有效的驱动因子，同时也是干旱区生态系统演变最为重要的制约因素^[2-3]。因此，土壤水分含量的准确监测在干旱区陆地表面过程的科学研究中具有重要意义^[4-5]。

由于时空变化较大，传统的监测方法无法及时准确地获取大面积的土壤水分含量信息，加之耗时费力，具有很大局限性。随着遥感技术的快速发展，利用卫星遥感进行大范围土壤水分含量监测已成为最为有效的研究手段。其中，基于光学遥感数据构建干旱指数模型以综合反映区域土壤水分含量状况、实现旱情监测的研究在国内外取得了较大进展。Sandholt等^[6]在分析了归一化差值植被指数（Normalized Difference Vegetation Index，NDVI）与地表温度（Land Surface Temperature，LST）的关系后提出了温度植被干旱指数（Temperature Vegetation Dryness Index，TVDI），以表征土壤水分含量状况；随后，Wan等^[7]与Jain等^[8]分别利用TVDI指数监测了美国南部平原与印度西部的土壤水分含量状况和干旱情况，并进行了干旱等级划分；Patel^[9]与Chen等^[10]等利用TVDI分别对印度北部和我国黄淮海平原的半湿润地区土壤水分含量的空间分布状况进行了评估；曹雷等^[11]、刘婕等^[12]与王娟等^[13]基于TVDI分别反演了新疆艾比湖地区、精河流域绿洲和四川绵竹山地平原过渡带的土壤水分含量，并探讨了区域土壤水分含量的时空变化特征，取得了较好的反演效果；詹志明等^[14]、阿布都瓦斯提·吾拉木^[15]提出了基于光谱特征空间的垂直干旱指数（Perpendicular Drought Index，PDI），并对裸露地表土壤水分含量进行了监测；冯海霞等^[16]在PDI基础上进行了改进，提出了归一化干旱监测指数（Normalized Perpendicular Drought Index，NPDI），并对宁夏永宁县的土壤水分含量进行了反演；为摆脱土壤线对模型的影响，刘英等^[17]提出了不依赖于土壤线参数的土壤湿度监测指数（Soil Moisture Monitoring Index，SMMI），评估了陕西省榆林市神东矿区的土壤水分含量分布情况；葛少青等^[18]基于TVDI、PDI和NPDI三种干旱监测方法，对甘肃省疏勒河中下游沼泽湿地进行了土壤水分含量遥感反演，并结合野外数据进行了比较和验证。上述光学遥感指数模型通过构建不同的特征空间，进而利用特征空间的形状、土壤线或者二维光谱空间的距离反演土壤水分含量，简单实用，在裸地和植被覆盖地区，在山地和平原地区、在干旱和半湿润地区，在国内和国外均已取得较好的监测结果。居延泽地区位于黑河下游，气候干旱，降水稀少，季节分配不均匀，蒸发强烈，生态环境极端脆弱，土地荒漠化问题突出，水是制约当地生态环境和经济可持续发展最为关键的因子。目前，有关区域土壤水分含量遥感监测的研究非常有限。本文基于Landsat-8 OLI/TIRS遥感数据，构建该地区多种干旱指数模型，结合不同深度的野外实测土壤水分含量数据，进行拟合精度验证，反演研究区多时相的土壤水分含量，并结合随机森林分类法探究土壤水分含量在不同地类的分布规律，以期为了解区域土壤水分含量动态及其与植物生长的关系、为区域生态环境保护与可持续发展提供科学依据和决策参考。

1.   研究区概况

居延泽位于内蒙古自治区阿拉善盟额济纳旗东南部（图1），介于41°47′ ~ 42°01′ N，101°31′ ~ 101°51′ E之间，面积约800 km²，北邻阿尔泰山脉，南接巴丹吉林沙漠。地势中间高，两侧低平，由东、西居延泽两个子盆地组成，西居延泽北部有居延泽干涸后残留的湖沼——天鹅湖^[19]。地质历史时期，古居延泽与苏古淖尔、嘎顺淖尔相连，共同组成黑河下游终端湖。居延泽深居欧亚大陆腹地，气候干燥，为典型的大陆型气候，年平均降水量为37.0 mm，年平均蒸发量为3841.51 mm，年平均温度为8.9 ℃。该地区生态环境极端脆弱，地表覆盖以沙地、裸地和盐碱地为主，兼有少量滩涂分布，植被稀疏，以旱芦苇和白刺等沙生植被为主，土地荒漠化、盐碱化问题突出^[20]。

图  1  研究区地理位置

Figure  1.  Location of the study area

下载: 全尺寸图片 幻灯片

2.   数据源及预处理

2.1   土壤水分含量实测数据

2017年8月14日至8月17日在研究区进行了野外采样。鉴于研究区自然条件极端恶劣，在考虑通达性的基础上，采用均匀布点的原则，并且尽可能涵盖研究区各种土地利用类型。对该研究区共设置50个土壤采样点，样点位置如图2所示，对各采样点人工开挖长宽深均约为50 cm的探坑，采用分层取样法沿垂直方向分别对0 ~ 10 cm、10 ~ 20 cm和20 ~ 30 cm的土层采集土壤样品，每个采样点采集3个土壤样品，总计150个土壤样品。样品采集后用铝盒密封，在实验室采用烘干法^[21]测定土壤体积含水量。50个样点中，0 ~ 10 cm、10 ~ 20 cm与20 ~ 30 cm的土壤水分含量均值分别为9.57%、14.60%和18.28%，最大值分别为38.07%、36.34%和48.36%，最小值分别为0.23%、0.31%和0.42%。

图  2  采样点分布图

Figure  2.  Map of sampling sites

下载: 全尺寸图片 幻灯片

2.2   遥感数据

采用2015年8月28日、2016年7月29日和2017年9月2日获取的Landsat-8 OLI/TIRS影像，轨道号/行号为133/31。Landsat-8由美国国家航空航天局于2013年2月11日成功发射，该卫星携带了两个主要载荷，即运行陆地成像仪（Operational Land Imager，OLI）和热红外传感器（Thermal Infrared Sensor，TIRS）^[22]。OLI共有9个波谱段，其中全色波段空间分辨率为15 m，其他波段空间分辨率为30 m。TIRS包含两个空间分辨率为100 m的热红外波段。对三期影像进行了辐射定标和大气校正等预处理。

3.   研究方法

3.1   温度植被干旱指数（TVDI）

Sandholt等^[6]发现LST与NDVI特征空间近似为三角形，因而提出了温度植被干旱指数的概念。任意一个NDVI值对应一组干边和湿边，因而LST-NDVI特征空间的拟合线可以反映区域土壤水分含量状况。LST-NDVI特征空间的计算表达式为：

式中，LST表示任意像元的地表温度（℃）；$ {LST}_{\mathrm{m}\mathrm{a}\mathrm{x}} $和$ {LST}_{\mathrm{m}\mathrm{i}\mathrm{n}} $分别表示相同NDVI条件下地表温度的最大值（℃）和最小值（℃）；$ {a}_{1} $、$ {b}_{1} $分别表示将LST-NDVI特征空间线性回归后干边的截距与斜率，$ {a}_{2} $、$ {b}_{2} $表示湿边的截距与斜率。

将公式（2）、（3）代入到（1）中可以得到TVDI的计算公式：

TVDI的值随着地表温度的升高而增加，土壤水分含量随TVDI值增大而减小，因而土壤水分含量与温度、TVDI值呈负相关关系。

其中NDVI与LST的计算公式如下：

式中，$ {\rho }_{NIR}{\text{、}}{\rho }_{R} $分别为近红外波段与红光波段的反射率；$ B\left(LST\right) $表示地表辐射亮度；$ {k}_{1} $、$ {k}_{2} $为Landsat-8 Band 10的反演常数，分别为774.89 $ W{\mathrm{m}}^{-2}{\mu \mathrm{m}}^{-1}{\mathrm{s}\mathrm{r}}^{-1} $、1321.08 K^[23]；$ \tau $、$ {L}^{\uparrow } $和$ {L}_{\downarrow } $分别表示大气透过率、大气向上辐射亮度和大气向下辐射亮度，三者均可通过大气校正参数计算器（http://atmcorr.gsfc.nasa.gov）查询得到；$ {L}_{\lambda } $表示热红外波段辐射亮度值；$ \varepsilon $表示比辐射率，参考前人文献[24]，水体的比辐射率取0.995，其他地表采用如下公式计算：

式中，$ {\varepsilon }_{\text{ surfare }} $为自然表面像元的比辐射率；$ Pv $为植被覆盖度；$ {NDV}{{I}}_{{veg}} $、$ {NDV}{{I}}_{soil} $分别为结合野外考察选定的纯植被像元与纯裸土像元对应的NDVI值。

3.2   垂直干旱指数（PDI）

詹志明等^[14]研究发现，使用Landsat-7 ETM + 遥感影像的红光与近红外波段反射率数据建立的光谱特征空间呈典型三角形。将该三角形下边线上的点线性回归即可得到土壤线^[25]。该三角形可以表现土壤水分含量状况和环境干旱程度，沿土壤线向上干旱程度加重。为了定量描述该规律，阿布都瓦斯提·吾拉木^[15]构建了PDI指数，其表达式为：

式中，$ {\rho }_{NIR}{\text{、}}{\rho }_{R} $分别为近红外波段与红光波段的反射率；m表示红光－近红外二维光谱特征空间中对三角形下边线线性回归得到的土壤线斜率。PDI指数定量描述了在红光－近红外二维光谱空间中土壤水分含量的分布规律，即PDI值越大，对应土壤水分含量越低的区域，干旱程度越严重；反之，PDI值越低，对应土壤水分含量越高的区域，环境状况越湿润。

3.3   归一化干旱监测指数（NPDI）

在PDI基础上，冯海霞等^[16]选择对水分敏感的红光和短波红外波段，以短波红外与红光的差值和短波红外与红光之和构建光谱空间，并且考虑了土壤和植被的光谱特征分异规律，提出了NPDI指数，其表达式为：

式中，$ {R}_{d} $和$ {R}_{s} $分别表示短波红外和红光波段反射率之差与反射率之和；m表示$ {R}_{d} $-$ {R}_{s} $二维光谱特征空间中对三角形下边线线性回归得到的土壤线斜率。

3.4   土壤湿度监测指数（SMMI）

刘英等^[17]发现在近红外和红光波段的二维光谱特征空间中，任意一点到原点的距离可以表征土壤水分含量的高低，在此基础上提出了SMMI指数。该指数的表达式为：

式中，$ {\rho }_{NIR}{\text{、}}{\rho }_{R} $分别为近红外波段与红光波段的反射率。

4.   结果与分析

4.1   四种干旱指数精度验证与分析

将2017年实测土壤水分含量与TVDI、PDI、NPDI和SMMI四种指数进行线性拟合，结果如图3所示。由图3可见，四种干旱指数在不同深度土壤中均随着土壤水分含量的增加而呈现下降趋势。从不同土层看，四种干旱指数均随着土层深度的增加拟合效果由优变劣，0 ~ 10 cm拟合效果最优，10 ~ 20 cm次之，20 ~ 30 cm最劣。从不同干旱指数看，在0 ~ 10 cm处，干旱指数与土壤水分含量的拟合效果从优到劣依次为TVDI、NPDI、SMMI和PDI；在10 ~ 20 cm处，拟合效果从优到劣的顺序同上；在20 ~ 30 cm处，干旱指数与土壤水分含量的拟合效果略有不同，从优到劣依次为NPDI、TVDI、SMMI和PDI。

图  3  TVDI、PDI、NPDI和SMMI与不同土层土壤实测水分含量之间的关系

Figure  3.  Relationship between TVDI, PDI, NPDI, SMMI and measured soil moisture content at different depths

下载: 全尺寸图片 幻灯片

综上，TVDI的综合表现最优，尤其在0 ~ 10 cm土壤水分含量的拟合精度上表现突出，与实测土壤水分含量值相关性显著，R²达到了0.76。推测原因是TVDI在反演土壤水分含量时，耦合了地表温度与植被指数这两个描述地表特征的重要参数，通过地表温度与植被指数特征空间分析土壤水分含量状况。干旱区地表蒸发旺盛，地表温度直接影响土壤水分含量的蒸发速率，由于干燥缺水，土壤水分含量成为植被生长最为重要的胁迫因素，因而TVDI在研究区的应用中具有明显优势；PDI和NPDI的表达式是基于反射率以及通过二维光谱空间拟合的土壤线确定的。土壤线的直线特征仅在一定范围内存在，会随土壤类型、有机质和矿物组成等的不同而动态变化，并非一条固定的直线。因此，对于地表状况存在差异的区域，PDI和NPDI中的m值是动态变化的。在本研究区中存在沙地、裸地和盐碱地等不同土壤类型，将m值设定为固定值可能是引起误差的原因^[26]。NPDI在PDI的基础上进行了改进，加入了对土壤水分含量变化敏感的短波红外波段，这可能是NPDI综合表现要优于PDI的原因。SMMI虽然避免了土壤线造成的监测精度的不确定性，但把水分对反射率的复杂影响简单抽象为光谱空间的点到原点的距离，可能是导致拟合精度较低的原因。因此，本研究选用TVDI与0 ~ 10 cm土壤实测水分建立的线性经验模型反演研究区土壤水分含量，以探求土壤水分含量在研究区的分布规律。

4.2   基于TVDI的0 ~ 10 cm土壤水分含量反演结果与分析

基于上述经验模型，得到研究区2015年8月28日、2016年7月29日与2017年9月2日的土壤水分含量反演结果（图4），采用随机森林分类法^[27]将研究区的土地分为裸地、沙地、盐碱地、植被和滩涂五类，三年分类精度为91.27%、92.74%和90.38%。分别统计三年各地类平均土壤水分含量，结果如表1所示。研究区天鹅湖附近与北部的铁路两侧水分较高，由于天鹅湖附近以滩涂为主，并附生植被，北部的铁路两侧建有防风带，对水分的保持有一定作用。研究区南部的东西两侧分别为东居延泽与西居延泽古湖区，以古湖岸堤相隔。由于东、西居延泽古湖区表面覆盖数厘米厚的盐壳，而光学遥感由于不具有穿透能力，所以居延泽古湖区内土壤水分含量反演值表现略低。西居延泽古湖区内部分区域存在稀疏分布的低矮植被，在反演结果的对应区域出现点状分布的相对高值，与实际相符。根据表1，同一年份的土壤水分含量由低到高均为沙地 < 盐碱地 < 裸地 < 植被 < 滩涂。滩涂与植被多分布于天鹅湖附近，地势较低，且植被可以遮挡阳光直射，减少水分蒸发，因此植被覆盖区域的土壤水分含量较高；沙地、盐碱地与裸地三年间土壤水分含量一直维持在较低水平，且变化幅度很小，推测原因可能是上述地类大多平坦无遮挡，蒸发强烈，且相较于植被与滩涂水分涵养能力差，无法有效保持水分。裸地与沙地的土壤构成存在差异，裸地的土壤持水性能要优于沙地^[29]，所以裸地土壤水分含量较沙地高。模型反演结果与研究区实际情况保持了较高的一致性。

图  4  TVDI模型土壤水分含量反演结果

Figure  4.  Retrieved soil moisture content using TVDI model

下载: 全尺寸图片 幻灯片

表  1  土壤体积含水量分类统计

Table  1.  Statistics of soil moisture content

地类Land type	2015年8月28日August 28, 2015			2016年7月29日July 29, 2016			2017年9月2日September 2, 2017
	像元比例（%）Pixel ratio	平均土壤水分含量（%）Average soil moisture content		像元比例（%）Pixel ratio	平均土壤水分含量（%）Average soil moisture content		像元比例（%）Pixel ratio	平均土壤水分含量（%）Average soil moisture content
裸地	12.73	4.65		12.33	5.07		12.47	5.57
沙地	63.59	2.13		68.59	2.19		65.69	2.88
盐碱地	18.36	4.10		15.63	4.37		12.53	4.15
滩涂	1.21	31.76		0.52	21.11		2.32	47.70
植被	1.42	9.24		0.61	15.75		2.72	18.33
水体	2.69	—		2.32	—		4.26	—

下载: 导出CSV 

| 显示表格

5.   结论

（1）在不同土层中，TVDI、PDI、NPDI与SMMI均与土壤实测水分呈负相关，即TVDI、PDI、NPDI与SMMI值越大，土壤水分含量越低，土壤越干燥。

（2）在不同土层中，TVDI、PDI、NPDI与SMMI均与0 ~ 10 cm土壤水分含量具有较好的拟合效果，随着土层深度的增加，拟合精度逐渐变劣。TVDI的综合表现最优，尤其在0 ~ 10 cm土壤水分含量的拟合精度上表现突出，与实测土壤水分含量值相关性显著，决定系数达到0.76，拟合精度最高。

（3）基于TVDI指数模型的反演结果表明，不同地类土壤水分含量存在差异，呈现出沙地、盐碱地、裸地、植被、滩涂依次升高的规律，可能与土壤持水性能、植被覆盖、地形和地势等因素有关。其中，沙地、盐碱地和裸地的土壤水分含量一直维持在较低水平，变化幅度不大。