基于Landsat 8数据的荒漠土壤水分遥感反演

高 琪, 彭 杰, 冯春晖, 宋 奇, 马自强, 王玉珍

(1.塔里木大学 植物科学学院， 新疆 阿拉尔 843300;2.北京大学 地球与空间科学学院， 遥感与地理信息系统研究所， 北京 100000)

水资源是地球生态系统的重要组分，并在全球水汽循环中发挥着积极作用，土壤水分作为水资源的重要分支，是连接植被、土壤和大气的关键因子，随着全球气候变暖，三者间物质和能量循环加剧，从而土壤水分蒸散严重，干旱现象频繁发生，如何合理的监测土壤水分运移，改善旱情是当今全球的热点问题[1-4]。传统土壤水分监测是基于地面有限点的测定，不但耗损大量经济成本和人力资源，而且空间连续性缺乏，时空信息获取困难，导致传统方法难以满足大范围土壤水分监测的需求[5]。遥感技术因具有宏观性、高效性、时效性、经济成本低和地面限制因素少等独特优势，已逐渐发展为现今土壤水分研究的关键方法[6]。

20世纪60年代至今，以遥感监测土壤水分的技术日渐成熟，并逐渐形成了微波遥感法、热惯量法、温度—植被干旱指数法和光谱指数法等[7-8]。其中，微波遥感在监测过程中易受地形、植被覆盖度和地表温度的影响[9-10]；热惯量法存在卫星数据获取时效性差、遥感图像解译效果低等劣势[11]；温度-植被干旱指数法不适用于高植被覆盖区域土壤水分的监测；光谱指数法在监测过程中虽对遥感数据要求较多，但在植被和土地类型多样且覆盖程度复杂区域有较好的监测效果，因此，以光谱指数法监测地貌和植被类型复杂的空台里克地区土壤水分变化，其适应度更高，监测效果更为显著。鉴于光谱指数监测土壤水分的特性，国内外诸多学者以光谱指数法对土壤水分监测进行了大量研究，如Ghulam等[12]通过构建Nir-Red特征空间，描述植被覆盖信息，发现改进垂直干旱指数(modified perpendicular drought index， MPDI)比传统指数能更好的反演旱情；张学艺等[13]结合MODIS-EVI得到改进型干旱指数(modified temperature vegetation dry ness index， MTVDI)，通过构建各指数与土壤含水量的相关函数，表明经改进后MTVDI的监测效果显著；姚云军等[14]发现并提出了短波红外土壤湿度指数(shortwave infrared soil moisture index， SIMI)，以SIMI反映土壤水分信息有较好的效果。诸多研究结果表明基于单一光谱指数研究均能反映土壤水分表征旱情，但面对植被覆盖复杂、地理因素多样，同时外延到研究区以外大范围监测研究时，其适用性尚未验证。为此，有学者在光谱指数法的基础上构建综合干旱监测模型，如孙灏等[15]比较13种典型光谱指数的适用特点，构建了农业干旱遥感监测的分类体系；葛少青等[16]在综合NDPI，TVDI和PDI共3种指数监测干旱区土壤水分取得了良好效果；Du等[17]考虑植被、土壤和降水亏缺等干旱因素，以主成分变换构建综合干旱监测指数(spatial drought index，SDI)，表明综合光谱指数监测土壤水分的精度和适用性均优于单一指数，同时以光谱指数监测土壤水分的发展逐步呈现为单指标向多指标、单因素向多因素发展的整体趋势。

迄今为止，以多种光谱指数、多因素为模型因子对土地利用类型和地表景观复杂多样区域土壤水分监测研究的相关报道甚少，其监测效果有待进一步验证。本文以南疆干旱区阿克苏地区的空台里克冲积扇为研究区，该研究区以沙漠—绿洲过渡带为主要土地利用类型，同时地表景观类型多样且交错分布，其中，包括高覆盖区域植被—中覆盖区域植被—低覆盖区域植被、裸地和农田；区域内土壤质地空间分布变异性强，从而地表粗糙度差异较大，特别是在盐结皮区域和非盐结皮区域之间差异较为显著。为更好的探究综合光谱指数法对地表景观复杂多样的空台里克区域荒漠土壤水分监测的效果及适用性，因此，本研究以Landsat 8数据拟构建一种基于多个光谱指数，同时考虑地形(digital elevation model， DEM)和环境数据等因子的土壤水分综合监测模型，为干旱区旱情监测和荒漠土壤水分空间分布研究提供一定理论依据和方法支持。

1 材料与方法

1.1 研究区概况

研究区位于南疆阿克苏地区温宿县境内的空台里克冲积扇(40°40′10″—41°32′36″N，80°45′09″—81°12′45″E)，北靠天山山脉，西临塔里木盆地边缘，地势北高南低，西高东低[18]，平均海拔1 200 m，面积约2 176 km2。研究区属于典型沙漠—绿洲过渡带，光热资源充足，气候干燥，年平均积温4 200 ℃，无霜期212 d，年平均降水55.1 mm，而年平均蒸发却高达2 428 mm；区域交通单一，以南北方向长110 km和东西方向长80 km的省级公路为主且贯穿整个研究区，沿省级公路两侧由北向南土地利用类型依次为沙漠、沙漠—绿洲过渡带、河流和农田；土壤质地以沙土、沙壤土为主；地表植被类型分布多样，以柽柳、骆驼刺，盐惠木等荒漠植被为主且呈高(植被全覆盖)—中(植被半覆盖)—低(裸地)的趋势集中分布于大面积沙漠—绿洲过渡带，小面积农田以棉花种植为主。

1.2 遥感数据获取及预处理

通过美国地质调查局USGS官网(https:∥earth explorer.usgs.gov/)下载过境时间为2019年6月15日，轨道号为146/31,146/32，云量小于10%，空间分辨率为30 m的2景Landsat 8遥感影像。通过ENVI 5.3软件完成影像预处理，具体包括：辐射定标、FLAASH大气校正、影像镶嵌、裁剪和波段运算，经预处理后对研究区遥感影像进行支持向量机分类(SVM Classical)；以遥感影像提取比辐射率和黑体辐射亮度等数据信息用于地表温度计算。通过ArcGIS 10.2软件获取各光谱指数信息。

通过地理空间数据云(http:∥www.gscl oud.cn/)下载空间分辨率为30 m的DEM数据，以ArcGIS 10.2软件完成洼地填充、影像裁剪、几何精度校正和高程数据提取。

1.3 土壤样品采集与测定

本研究于2019年6月15日Landsat 8卫星过境当日完成土壤样品采集，共计160个，为使野外采样充分反映研究区土壤质地及植被覆盖，同时考虑到研究区交通较为单一的情况，沿东西走向长110 km的省级公路两侧以近似等间距布设20个采样样方，各样方内分布8个采样点且遵循如下采样原则：①远离省级公路300 m以上；②避开河流和农业用水渠等；③依据土地利用分类结果，在不同覆盖类型内均布设样方；④以五点法采集各样方表层0—20 cm的土样并均匀混合。为防止土样水分蒸发而导致含水量数据误差，将采集的土样快速装入自封袋中密封保存，带回实验室称取50 g左右，采用质量烘干法测定土壤含水率，每个样品设置3次重复，各重复间相对误差控制在5%以内，取3个重复样的平均值为该样品的土壤水分实际测定值。剩余土壤样品放置于室内风干，经充分研磨后统一过筛(2 mm)，使用稀释溶液法测定pH值和电导率(EC)。

1.4 建模方法与模型精度评价

本研究选取经典线性模型偏最小二乘法(partial least squares regression， PLSR)、二分类模型支持向量机(support vector machine， SVM)和机器学习模型随机森林(random forest， RF)3种方法，分别通过The Unscrambler10.5和R语言软件构建模型。将160个样品实测数据在Excel软件中依次由小到大排列，以2∶1比例分成107个建模集样品，53个预测集样品，通过决定系数(R2)、均方根误差(RMSE)、平均绝对误差(MAE)、相对分析误差(RPD)对模型精度验证；当R2值越大、RMSE值越小，说明模型精度越高，预测效果越好；MAE表示模型误差范围，其值越低，模型精度越高；相对RPD值而言，当RPD<1.5时，表明模型没有预测能力；当1.5≤RPD<2.0，RPD≥3.0，表明模型预测能力较差，只能进行粗略估算；当2.0≤RPD<2.5时，表明模型预测能力较好；当2.5≤RPD<3.0时，表明模型预测能力很好；当RPD≥3时，表明模型预测能力极出色，能进行高精度估算[19]。其中，各精度评价指标通过如下(1)—(5) 公式计算：

(1)

(2)

(3)

(4)

(5)

1.5 研究方法

1.5.1 光谱指数 不同光谱指数描述的地表干旱特征有所差异，进而反映土壤水分表征旱情的结果也不一致[15]，由此造成单一光谱指数监测土壤水分精度较低，本研究在此基础上针对空台里克地区植被覆盖和土地利用等复杂情况，选取TVDI，NR和NDVI等31个光谱指数并充分考虑了不同指数对土壤水分反映的基本特征，以多个光谱指数为建模因子构建反演模型，探究研究区土壤水分变化规律。

1.5.2 地表温度(Ts) 地表温度是区域尺度地表物理参量的一个关键因子，作为地球表面能量平衡的重要指示指标，地表温度能间接的反映土壤水分状态表达旱情特征。其计算方法见公式(6)—(8)：

热红外辐射亮度值的表达式如下：

Lλ=〔εB(Ts)+(1-ε)L↓〕ε+L↑

(6)

式中ε为地面比辐射率数据;Ts为地面实际温度;B(Ts)为黑体在地表实际温度的热辐射亮度;τ为大气;L↑为大气向上的辐射亮度;L↓为大气向下由地表反射的辐射亮度。

黑体辐射亮度B(Ts)的计算公式如下：

(7)

式中：Lλ表示大气校正后的热红外波段辐射定标值。

得到黑体辐射亮度数据后通过普朗克函数计算地表温度(Ts)。

(8)

式中：Ts为地表实际温度;B(Ts)为黑体在Ts的热辐射亮度;K1，K2为卫星预设常数；Landsat 8预设K1=774.89 W/(m2·sr·μm)，K2=1 321.08 K。

1.5.3 地形数据(DEM) 研究区地势北高南低，自西北向东南倾斜，考虑到地势高低起伏的特征是影响土壤水分、植被分布和气候特征等地域分异的关键因子[20]，因此，本研究以DEM数据作为建模因子，分析因地形变化引起土壤水分分异的规律，更进一步了解空台里克地区土壤水分分布的总体特征。

2 结果与分析

2.1 土壤属性统计性描述

表1对实地采集的160个表层土壤样品基础属性进行统计描述。由表1可知，表层土壤水分最小值为4.30%，而最大值为44.29%，极差高达39.99%，平均值为21.58%，说明表层土壤水分(0—20 cm)分布差异较大，引起这种差异的可能原因为不同土壤质地和植被覆盖，从而造成研究区表层各部分土壤的持水能力差异。EC，pH值的极差分别为78.53 dS/m，1.8，平均值分别为24.01 dS/m，8.16，变异系数分别为44.57%和5.12%，同时表层土壤水分变异系数为30.74%。而变异系数一般分为强变异(>50%)、中等变异(20%～50%)和弱变异(0～20%)[21]，引起土壤水分和EC值中等变异一是植被覆盖差异较大，高植被覆盖区域植被对土壤水分和养分元素需求较高，低覆盖区域直接受太阳高温辐射作用较强；二是研究区人为活动现象明显，新开垦农耕田、撂荒地不断增加以及人为设立保护区存在，导致土壤结构人为破坏严重，使土壤水分和盐分较强变异。土壤属性的偏度和峰度值分别描述数据偏斜程度和各特征值的高低分布情况，表中所有数据偏度值均大于0，峰度均与3有一定差距，表明数据总体呈现偏态分布且均不符合正态分布规律。

表1 土壤样品基础属性

2.2 建模因子与土壤水分的相关性分析

为探究各建模因子与土壤水分相关性，将地表温度(Ts)和地形数据(DEM)以及经波段运算的31个光谱指数与实测土壤水分在Excel 2019软件中进行相关分析并附各因子计算公式(具体见表2)。通过表2可知，26个光谱指数与土壤水分为极显著相关，其中，温度植被干旱指数(TVDI)与土壤水分相关性最高，为0.53，表明26个光谱指数均能在不同程度内指示土壤水分，但TVDI指示效果优于其余光谱指数；绿度差值植被指数(green difference vegetation index， GDVI)、土壤改良植被指数2(modified soil adjustment vegetation， MSAVI2)和非线性植被指数(nonlinear vegetation index， NLI)的相关系数分别为0.06，0.08,0.15，与土壤水分均未达到显著或极显著相关，表明3个光谱被指数均不能有效的指示土壤水分信息，不适宜作为建模因子，地表温度和地形与土壤水分相关系数分别为-0.41和-0.43，均为极显著负相关，表明两者均能影响土壤水分分布且相关性越高，土壤含水量则越低。因此，本研究在相关系数的基础上筛选出26个光谱指数，地表温度和地形为反演模型建模因子。

表2 建模因子与土壤水分相关系数

2.3 不同建模方法对比分析

表3为PLSR，SVM，RF3种模型方法构建土壤水分建模集和预测集的精度统计表。通过表3可知，建模集中3种模型方法的精度差异较大，RF模型的R2最高，为0.93，而RMSE和MAE均为最小值，分别为1.49，1.09；SVM模型的R2最低，为0.43，RMSE和MAE分别为3.66,2.82，均为建模集中最大值；PLSR模型性能介于二者间。综上所述，建模集中，RF模型性能明显优于PLSR和SVM方法，反演效果优劣排序依次为RF>PLSR>SVM。在预测集中，3种模型的各指标相较于建模集均无明显下降，表明模型比较稳定，但三者相比较而言，仍然是RF模型的各项指标明显优于PLSR和SVM，而PLSR和SVM二者间相比，PLSR模型的R2明显高于SVM模型，而RMSE和MAE要低于SVM模型。就RPD指标而言，PLSR模型和SVM模型的1.5≤RPD<2.0，表明这两种模型只能对土壤水分进行粗略估计。RF模型的RPD值为3.90(3.90>3.0)，表明RF模型有出色预测能力，能精准反映研究区土壤水分。综合考察各模型建模集与预测集的评价指标，3种模型的预测性能和稳定性从高到低排序依次为RF，PLSR和SVM，表明神经网络模型比线性模型和二分类模型监测反演土壤水分更理想，模型稳定性更强，精度更高。

表3 不同模型下土壤水分反演精度

2.4 不同建模方法可信度检验

分析土壤水分在PLSR，SVM和RF模型预测结果下与实测土壤水分的相关性，进一步确定预测土壤水分数据的精准性(见图1)。通过图1可以看出，RF模型(c)实测水分与预测水分R2最高，为0.91，表明基于机器学习方法的RF模型反演土壤水分最优，PLSR(a)和SVM(b)的R2为0.63,0.32，相比较RF模型，对土壤水分反演效果较差，但都具有一定的可信度。

图1 不同模型反演数据与实测数据相关性

2.5 对比不同模型的土壤水分空间分布

图2为最优模型(RF)反演表层(0—20 cm)荒漠土壤水分的分布情况，通过图2可知，土壤含水量以0～5%和5%～12%集中分布于研究区东北部；15%～20%或>20%大面积分布于北部、中部和南部区域。

图2 最优模型反演土壤水分分布结果

表4对土壤水分反演结果的分布面积进行统计，从表4可知，0～5%土壤含水量面积最小，为104 km2；>20%的分布面积最大，为1 041.18 km2；0～15%之间的分布面积总占比为18.21%；而15%～20%和>20%为81.79%，表明研究区土壤含水量低值较少，总体呈偏高的分布特征。为更进一步探究研究区土壤水分变化特征，结合土地利用分类结果和最优模型土壤水分分布结果(图2)可知，东北部为沙漠区，土壤水分含量集中为0～12%，南部多为高作物覆盖农田区，土壤含水量多为15%～20%和>20%，同时农田内存在无作物覆盖的新开垦田和撂荒地，土壤含水量为0～15%，中部、北部为典型沙漠—绿洲过渡带区，土壤含水量多集中为15%～20%或者>20%，表明在不同土地利用类型下，土壤水分含量存在较大差异。

表4 土壤水分分布面积及比例

3 讨 论

3.1 多因素建模可行性

本研究以26个光谱指数、环境数据因子和地形因子(DEM)构建PLSR，SVM和RF3种模型反演土壤水分。研究区地势北高南低，由1 400 m向980 m依次降低，而土壤含水量与地形相关系数为-0.44，达极显著相关，表明DEM变化对土壤水分分布有较大影响，而本研究引入DEM数据作为模型因子，校正因地形引起土壤水分反演精度的偏差，该结果与蔡亮红等[20]在TVDI基础上引入DEM数据校正，校正后TVDI能更好的反映土壤水分状况的研究结论较为一致。地表温度在土壤湿度和地、气相互循环过程中有重要作用，同时能间接反映土壤水分状况表征旱情分布，因而对地表温度数据的研究是不可或缺的内容[21-22]，本研究地表温度与土壤水分为极显著负相关，表明土壤水分随地表温度的增大而迅速减少，以地表温度为环境数据因子，构建土壤水分综合模型，对土壤水分高精准反演有重要价值，与蔡国印等[23]以热惯量法研究土壤水分的结果较为相似。结果表明在模型中引入地形和地表温度数据对土壤水分反演模型的精度有较大提升。

3.2 土壤水分空间分布特征

通过最优模型RF反演空台里克地区土壤水分，结合图2和表4可知，东北部沙漠土壤含水量最低，集中为0～12%，中部沙漠—绿洲过渡带土壤含水量最高，集中在15～20%和>20%且分布面积为81.79%，由于沙漠和沙漠—绿洲过渡带土壤质地相差较大，土壤结构疏密度和土壤持水能力不同，因而在不同土地类型下土壤含水量差异较大，这与姜红等[24]研究结论较为一致。相比较，沙漠—绿洲过渡带土壤含水量普遍较高，由于卫星遥感在监测中具有不确定性，空台里克是典型的高盐渍化区域，盐分表聚现象严重，地表上层形成大量盐结皮，而下层土壤水分受盐结皮封阻蒸发困难，导致下湿上干的现象，与Peng等[25]对空台里克地区盐渍化程度研究结论相符；南部农耕田土壤含水量呈0～12%，15～20%和>20%的空间分布特征，而农耕田内存在大量新开垦农田和撂荒地，在不同作物覆盖区域，土壤水分的运移速率有所差异，导致变异性较高(44.57%)，而含水量低值区(新开垦田、撂荒地)受人为因素干扰，土壤表层结构破坏严重，下层土壤水分受6月高温蒸发上移现象明显，与王思楠等[26]研究结论较为相似。以光谱指数法研究土壤水分，受植被覆盖、大气、土壤背景粗糙度等自然因素、遥感传感器等自身因素影响较为严重，后续研究需进一步探索不同影响因素对土壤水分的贡献率，筛选出一套适合研究区甚至新疆南部类似干旱区域的高精度土壤水分反演方法。

4 结 论

本研究有效结合卫星遥感和地面实测数据，通过PLSR，SVM和RF 3种方法构建多指数、多因素的综合监测模型反演土壤水分，选取最优模型对南疆阿克苏地区空台里克表层荒漠土壤(0—20 cm)水分预测分析。得到以下基本结论：

(1) 各因子与土壤水分相关性分析表明TVDI，NR，NDVI，EVI，GOSAVI等26个光谱指数，环境因子地表温度(Ts)、地形(DEM)均与土壤水分极显著相关，适合作为南疆荒漠土壤水分遥感监测的指示因子。

(2) 对比PLSR，SVM和RF3种模型各评价指标精度，RF模型的R2为0.91，RMSE为1.49，同时RPD值高达3.90，表明RF为土壤水分反演最优模型。

(3) 以最优模型RF反演研究区表层土壤水分，呈现沙漠0～12%集中分布，沙漠—绿洲过渡带15%～20%和>20%广泛分布，农田区域交错分布的土壤水分总体分布特征，其结果与土地利用分类特征和表1描述性统计结果较为一致。