基于改进型光谱指数的荒漠土壤水分遥感反演

时间：2024-07-28

高琪，王玉珍，冯春晖，马自强，柳维扬，彭杰，季彦桢

(1.塔里木大学植物科学学院，阿拉尔 843300； 2.地质环境监测站，昌吉 831100；3.北京大学地球与空间科学学院遥感与地理信息系统研究所，北京 100000)

0 引言

土壤水分作为地球水资源的重要组成部分，不仅是植被生长，群落演替和生态防护的关键因子，而且在大气、植被和土壤的物质、能量循环中起着重要作用。西北南疆地区位处内陆，受远洋气候和高山环绕等多方面因素影响，导致该地区干旱现象尤为严重，不但造成生态环境退化，土壤质地破碎，而且制约了区域经济的发展，因此，合理、高效地监测和规划土壤水分对区域旱情防控和荒漠化治理具有重要价值[1-2]。传统土壤水分研究，虽然测定精度高，但样点数量少，代表性差，难以准确地反映区域土壤水分变化规律。而遥感技术具有时效性强、范围性广和成本价值低等优势，不仅实现了区域土壤水分的大范围监测研究，而且在干旱监测和水土防治中有着极其重要的意义[3]。

20世纪中期，部分研究者通过卫星遥感监测天气学视角特征，以可见光—近红外、热红外和短波红外光谱特征空间，并结合植被、土壤、地表温度和地形数据等与土壤水分变化的相关特征，以光谱指数法对区域土壤水分特征开展研究[4-5]。詹志明等[6]通过反射率信息构建近红外-红光光谱特征空间，以SMMRS模型监测土壤水分，发现地表实测土壤水分与预测值曲线趋势一致，表明以反射率特征光谱空间监测土壤水分的可行性；李喆等[7]以特征光谱信息构建近红外-红光特征空间监测土壤水分，发现不同特征空间的适应范围有一定差异，但对土壤水分监测均具备较高精度，这类研究多基于可见光—近红外光谱，缺乏短波红外和热红外波段对土壤水分监测的研究。为此，相关学者如Thenkabail等[8]通过研究发现短波红外波段对土壤水分变化具有重要指示作用且对土壤水分反演精度提升有一定价值； Hain等[9]通过ALEXI模型估算土壤含水量，发现热红外波段对土壤水分变化敏感度更高，能有效反映出下垫面土壤水分状态。以上研究表明不同类型特征光谱信息对土壤水分监测有一定效果，但多基于特征光谱所构建的传统指数，在此基础上，本研究在传统光谱指数内加入对土壤水分信息敏感性更出色的热红外波段(b10)探索改进后光谱指数监测土壤水分的效果。

传统光谱指数监测土壤水分的研究多集中于单一的裸土区、植被覆盖区之间，如葛少青等[10]发现温度植被干旱指数(temperature vegetation dryness index，TVDI)、垂直干旱指数(perpendicular dryness index，PDI)和归一化干旱监测指数(normalized perpendicular drought index，NPDI)3种干旱指数均能很好地在不同土层深度处进行监测；蔡亮红等[11]选取4种传统植被指数并在此基础上加入短波红外波段进行改进，表明改进后的优势明显，但面对裸土与植被覆盖区交错分布的区域很难适用，因此本研究选取3种类型(土壤指数、植被指数和干旱指数)9个光谱指数并加入热红外波段(b10)进行改进，以期为南疆干旱/极端干旱区域荒漠土壤水分反演选取最佳的遥感指示因子。

综上所述，本研究选取空台里克为研究区，该区域属于典型荒漠-绿洲过渡带，地表景观类型复杂、植被覆盖类型多样；区域盐渍化现象严重，盐结皮地带土壤水分分布尤为特殊。为更好地探索不同类型的传统型和改进型光谱指数对荒漠区域土壤水分监测的适应性和模型精度，通过相关系数和方差膨胀因子筛选最优光谱指数，同时考虑地形数据作为建模因子构建线性和非线性土壤水分综合反演模型，以期为荒漠土壤水分定量遥感快速精准监测和水土治理提供一定的科学依据。

1 研究区概况及数据源

1.1 研究区概况

研究区位于南疆阿克苏地区空台里克自然保护区(E80°48′～81°12′，N40°41′～41°18′)，北靠天山山脉，西邻塔里木盆地，地势北高南低，见图1(a)。该区域年蒸发比高达40∶1，年均气温为11.2 ℃，年均日照时数为2 929 h，为典型干旱地带，同时又属于荒漠-绿洲过渡带，土地利用类型多为盐碱地荒漠，面积约占80%，其余为新开垦的农田、河流和沙漠(图1(b))；土壤质地多为沙壤土；植被类型主要包括骆驼刺、柽柳和甘草等耐旱植被；区域交通以一条贯穿研究区南北走向(110 km)的S215省级公路为主，且由北向南植被覆盖呈现低-中-高的总体趋势。近年来，人口数量不断上升，为满足区域经济，南部新开垦农田、农业灌水量不断增加，制约了生态环境的稳定与发展。

(a) 研究区位置及采样点(b) 土地利用分类

1.2 数据源及预处理

本研究选取的遥感数据主要包括Landsat8 OLI/TIRS影像和地形数据。通过USGS官网下载2019年6月15日轨道号分别为LC146_31和LC146_32云覆盖量低于20%的2景影像，以ENVI5.3软件完成辐射定标、大气校正和光谱指数波段运算等预处理过程，同时采用支持向量机(support vector machines,SVM)分类解译遥感影像地物信息；地形影像数据(30 m)通过地理数据空间云(http: //www.gscloud.cn/)下载，并以ArcGIS10.7 软件进行洼地填充和数字高程模型(digital elevation model,DEM)数据提取等过程。

1.3 土壤样品采集与室内分析

于2019年6月15日采集160个0～20 cm的表层土壤样品，并定位记录各点经度和纬度信息。根据研究区基础交通和植被覆盖状况，沿南北走向长110 km的S215省级公路布设20个不同植被覆盖度的采样样方(图1(a)所示)，各样方的面积为300 m×300 m，单个样方内采集8个表层土壤混合样，每个表层土壤混合样采用梅花形五点法取样，面积为(30 m×30 m)。当日内以质量烘干法测定各样品土壤质量含水量且测定一组重复，将重复相对误差控制在5%以内，取各组平均值为单一土壤样品的实际质量含水量。剩余土壤样品经自然风干、研磨和过筛(2 mm)后称取50 g土壤样品，以稀释溶液法(水土比5∶1)测定电导率和pH值。

2 研究方法

2.1 光谱指数计算

本研究在土壤指数、植被指数和干旱指数的基础上引入对土壤水分敏感性更高的热红外波段(b10)进行改进，其中，土壤水分遥感监测领域常用的光谱指数为裸土指数(bare soil index，BSI)、黏土指数(clay index，CI)、碳酸盐指数(carbonated index，Cal)、比值植被指数(ratio vegetation index，RVI)、归一化差值植被指数(normalized difference vegetation index，NDVI)、差值植被指数(difference vegetation index，DVI)、归一化差值水体指数(normalized difference water index，NDWI)、PDI和全球植被水分指数(global vegetation moisture index，GVMI)，计算公式见表1。

表1 常用传统和改进光谱指数及其计算公式

2.2 方差膨胀因子分析

VIF=1/(1-r2)，

(1)

式中r为各光谱指数间的相关系数。

2.3 反演模型构建

为确保建模集与预测集样本土壤含水量数据间的均衡性，将160个土壤样品的含水量数据依次从小到大排序，按照2∶1比例原则选取107个样本作为建模集，53个样本作为预测集，以优选出的光谱指数为建模因子，并结合地形数据，通过The Unscrambler X 10.5和R软件分别构建多元线性回归(multivariable linear regression，MLR)和随机森林(random forest，RF)综合模型反演土壤水分。

MLR模型通过系数判断解释变量所对应因变量的影响程度，优点为: ①简化了数据结构； ②分清各影响因素的贡献率。RF模型是一种新颖的机器学习模型，其预测精度高，但运算量较大[21]，主要分为： ①训练集的生成； ②随机特征选取； ③决策树构建。训练集生成的特点是随机、有放回，且新的训练样本数量与原始训练样本数量相同；随机特征的选取过程中引入了随机性，即随机选择输入变量，可提高模型的预测精度；决策树构建后，利用所有树预测值的平均值来代替最终预测结果。

2.4 模型精度评价

反演模型精度和稳定性采用决定系数R2、均方根误差RMSE和相对分析误差RPD进行评价。其中，R2表示实测值与预测值间的拟合程度，其值越大，表明实测值与预测值越接近。对于RPD而言，当RPD<1.5，说明模型不具备预测能力；当1.5≤RPD<2.0，说明模型只能粗略的预测；当2.0≤RPD<2.5，说明模型预测能力较好；当2.5≤RPD<3.0，说明模型能较精确的预测，当RPD≥3.0，说明模型预测能力极为出色[22]。计算公式为：

(2)

(3)

(4)

3 结果与分析

3.1 土壤属性统计性描述

为掌握空台里克研究区土壤水分特征和土壤基础属性，对160个表层土壤样品基础属性值进行统计，结果见表2。由表2可知，土壤水分、电导率和pH值的平均值分别为21.58%，24.01 dS·m-1和8.16，极差分别为39.99%，78.53 dS·m-1和1.8，标准差依次为44.10，113.64和0.18。其中，标准差描述数据的分布程度，表明研究区土壤水分和电导率整体较高，总体分布趋势明显。而变异系数可反映数据间离散趋势，三者变异系数为电导率(44.57%)>土壤水分(30.74%)>pH值(5.12%)，表明土壤水分和电导率变异性较强[23]。引起土壤水分和土壤电导率的变异性较强的因素为： ①土壤质地、植被覆盖程度和盐分特质性等区域生态环境； ②人工自然保护区、地下水开采和人工灌溉因素。

式中，x1、x2、x3、x4分别代表指标下限、上限、最优值下限、最优值上限，要求x1

表2 土壤样品基础特征描述性统计

3.2 建模因子筛选

3.2.1 光谱指数相关性分析

表3为传统光谱指数和改进光谱指数与土壤水分相关系数及显著性检验结果。

表3 光谱指数相关系数分析

①**为P<0.01水平下极显著相关，样本数n为160。

通过表3可知，光谱指数BSI，CI，Cal，NDVI和PDI在加入热红外波段(b10)改进后，相关系数提升了0.02～0.11，其中，ECI的相关系数最高，为0.73，ECal相关系数提升幅度最大，为0.11，表明加入热红外波段(b10)能有效地增强以上传统光谱指数与土壤水分间的相关性，更利于指示土壤水分变化。而RVI，DVI，NDWI和GVMI改进后的相关系数下降了0.01～0.16，其中，ENDWI相关系数降为-0.07，与土壤水分间无显著相关，说明引入热红外波段降低了这些指数与土壤水分间的响应程度，因而改进光谱指数ERVI，EDVI，ENDWI和EGVMI不适应作为本研究荒漠土壤水分遥感指示因子。综合考虑下，采用改进后有提升的光谱指数同原本达到极显著水平但改进效果不显著的因子来建模，即EBSI，ECI，ECal，ENDVI，EPDI，RVI，DVI，NDWI，GVMI；同时，选取传统光谱指数BSI，CI和Cal等9个因子构建模型进行精度对比。

3.2.2 各指数间方差膨胀因子分析

为更进一步筛选建模因子，以VIF检验传统和改进光谱指数因子多重共线性，结果见表4—5。由表4可知，传统型光谱指数DVI，RVI和NDVI间VIF分别为11.48，16.10和10.51，均大于10，相关系数大小关系依次为NDVI>RVI>DVI，其余各光谱指数VIF均小于10，因此仅保留光谱指数BSI，CI，Cal，NDVI，PDI，GVM和NDWI作为传统建模因子。由表5可知，DVI与RVI间VIF为11.48，其余各光谱指数VIF均小于10，因此，选取光谱指数EBSI，ECI，ECal，ENDVI，EPDI，GVMI，RVI和NDWI为改进建模因子。

表4 传统光谱指数间方差膨胀因子

表5 改进光谱指数间方差膨胀因子

3.3 反演模型精度验证

表6为线性MLR和非线性RF方法构建典型干旱地区荒漠土壤水分传统和改进综合反演模型。由表6可知，MLR-I和RF-I模型建模集的R2分别为0.64和0.88，预测集R2为0.59和0.86，RPD为1.48和2.63，相比较，改进后MLR-II和RF-II模型建模集的R2为 0.75和0.92，分别提升了0.11和0.04；而预测集R2和RPD分别提升了0.12,0.05和0.35,0.49；依据RPD指标，MLR-I模型没有预测能力，RF-I模型能准确地预测土壤水分，而改进后MLR-II模型能粗略地预测土壤水分，RF-II模型能出色地预测土壤水分；综合分析表明改进模型土壤水分预测效果整体优于传统模型。

表6 线性和非线性反演模型精度验证

基于改进后线性和非线性模型，通过各指标间精度变化，以筛选出荒漠土壤水分反演的最优模型。其中，MLR-II建模集的R2为0.75，预测集的R2为0.71，RPD为1.83；而RF-II建模集的R2为0.92，预测集的R2为0.91，RDP高达3.12，RF-II各指标均优于MLR-II模型，表明非线性模型预测精度相比线性模型效果更好，精度更高。

3.4 土壤水分空间分布特征

本研究选取改进后MLR-II模型和RF-II模型预测荒漠土壤水分含量，并进行制图分析，根据干旱分级标准[24]和图1(b)描述土壤水分空间分布特征。通过图2可知，土壤水分空间分布趋势呈现一定的相似性，但不同土地利用类型区域含水量却具有差异，其中，[15,20)%和[20,40)%的土壤含水量大面积集中在荒漠-绿洲过渡带，而该区域土壤含水量整体偏高； [0,5)%和[5,12)%的土壤含水量主要集聚在东北部沙漠区，且在南部农田区域交错分布，中部条带状则多为[12,15)%。

(a) MLR-Ⅱ(b) RF-Ⅱ

图2 土壤水分空间分布特征

图3分别为不同改进模型土壤含水量面积统计结果，表明各等级土壤含水量分布面积呈现相似的增长或降低趋势，但依旧存在较明显差异。为定量化分析土壤水分，根据最优模型预测结果以图2(b)和图3(b)进行分析，[0，5)%和[5，12)%的土壤含水量集中分布于东北部沙漠区域，零星交错分布于南部农田区域，总面积为314.2 km2，南部区域人为干扰现象明显，撂荒地普遍存在，导致地表盐结皮破坏，从而在极端干旱条件下使土壤水分大量流失； [15，20)%和[20，40)%的土壤含水量主要分布于荒漠-绿洲过渡带，同时交错分布于农田区域，总面积为1 689.3 km2，荒漠-绿洲过渡带地处沙漠与农田中部，植被覆盖差异较大，土壤类型以裸地为主，裸露地表受盐结皮保水和隔离的作用使含水量普遍偏高，而农田高含水量区域受灌溉影响造成土壤水分交错分布的空间特征。

(a) MLR-Ⅱ(b) RF-Ⅱ

4 讨论

4.1 模型精度影响因子

在9个传统光谱指数的基础上加入热红外波段(b10)改进，选取相关性和多重共线性双重检验后的光谱指数为建模因子，并结合地形因子构建荒漠土壤水分综合反演模型。改进后MLR-II和RF-II模型预测集的R2和RPD分别提升了0.12，0.05和0.35，0.49，引起模型精度提升的可能原因包括：一是根据热红外波段的特性，在3.5～14 μm处波段信息能有效探测地表热变化且得到的信息量更大[25]，对土壤水分的反应更为敏感，表明加入热红外波段不仅能进一步丰富光谱指数信息，而且能有效增强对土壤水分的指示作用；二是不同波段或光谱指数间具有较强相关性[26]，在传统光谱指数内加入热红外波段，很大程度上降低波段数据间的冗余度，更好地解译光谱指数信息，使指示效果更为清晰，这与蔡亮红等[11]针对渭河-绿洲以短波红外改进传统植被指数的研究结果一致。

4.2 土壤水分空间分布驱动因素

选取最优模型RF-II(图2(b))描述土壤水分空间分异格局并进行驱动因素分析，其中，土壤含水量整体呈现东北部沙漠地区多为[0,5)%和[5,12)%，北部、中部等荒漠-绿洲过渡带地区多为[15,20)%和[20,40)%，南部区域呈现交错分布的趋势，结果与描述性统计值一致。东北部沙漠区域土壤含水量多为[0,5)%和[5,12)%，而温度、土壤质地和地形地貌与土壤水分相关性强，是土壤水分分布不同的主要驱动因素[27]，研究区年均温度为11.2 ℃，日照时数约2 900 h，且土壤质地为沙土，从而造成地表上层土壤水分滞留时间较短，土壤孔隙度较大，保水能力差，形成大面积干旱现象；其次，研究区地形呈北高南低趋势，而地形是土壤水分分异格局的主要驱动因素之一，北高南低造成区域高低含水量分布格局，这与蔡亮红等[28]发现地形校正后能更好地反映土壤水分状况相一致。

荒漠-绿洲过渡带土壤含水量多呈现[15，20)%和[20，40)%，引起该区域土壤含水量较高是由于研究区低植被覆盖区分布面积远大于高植被覆盖区，而遥感监测存在异物同谱现象，将低植被覆盖区域误分为裸地，其中，裸地区域位处高盐渍化地带，盐分表聚现象突出，该结果与Peng等[29]对空台里克地区土壤盐分含量及分布制图结论相一致。为此，对盐结皮区域进行遥感监测不确定分析从而探讨盐结皮作用对土壤水分分布的影响，在极端干旱地区研究发现，盐结皮对土壤水分蒸发抑制效果显著[30]，而盐结皮具有表层干底层湿的现象，干土层的厚度与盐结皮的厚度一致。本研究区盐结皮的厚度大致在0～10 cm，而卫星遥感监测的仅是地表0～3 cm的信息[31]，地表0～3 cm受降水、温度等因素影响导致过度干旱，若仅使用单一卫星遥感数据进行土壤墒情判读，很容易将盐结皮区域误判为严重干旱区域。而根据地面调查结果，盐结皮层土壤含水量普遍在10%以下，但盐结皮以下的土壤含水量相对高于20%，造成上干下湿的土壤水分分布特征，因而以多源数据(微波遥感、光学遥感)协同反演的研究具有很大的发展潜力。

南部农田区域土壤含水量呈现为交错分布的趋势，这是由于农田区域内部存在一定数量的新开垦农田和撂荒地，人为开垦破坏地表盐结皮及土壤孔隙度，从而使土壤水分蒸发速率加快，导致部分地区土壤含水量为[0，5)%或[5，12)%；同时耕地人为灌溉痕迹较为明显，容易造成判断误读现象，这与龚新梅[32]以遥感手段对新疆荒漠化治理研究中人为驱动因素对土壤水分产生影响的结果相一致。研究结果表明空台里克地区土壤水分分布受到自然(土壤质地、温度、地形、植被作物覆盖度和盐结皮机理)和人为(人类经济活动、人为水资源开采)2种驱动因素的影响，形成不同土壤水分分布特征。