基于环境变量和机器学习的土壤水分反演模型研究

王思楠 李瑞平,2 吴英杰 赵水霞 王秀青

(1.内蒙古农业大学水利与土木建筑工程学院， 呼和浩特 010018；2.内蒙古自治区农牧业大数据研究与应用重点实验室， 呼和浩特 010018；3.中国水利水电科学研究院牧区水利科学研究所， 呼和浩特 010020；4.内蒙古自治区测绘地理信息中心， 呼和浩特 010050)

0 引言

土壤含水率是地表能量平衡的重要决定参数，在全球水循环中起着重要作用[1]。同时在农业应用中，土壤水分含量不仅是作物生长发育的基本条件，也是作物产量估算、干旱监测的关键参数[2]。因此，准确、及时地反演土壤水分具有重要意义。

土壤水分的分布受多个相互作用的因素影响，如土壤特性、植被覆盖和气候条件[3]。因此，用传统的单点测量方法，如烘干法、数字探头等，相对难以有效地获得大规模的土壤水分信息[4]。与地面单点测量方法相比，遥感技术因其覆盖面广、时效性强、成本低，逐渐得到应用[5-6]。光学遥感作为最早和最成熟的地球观测技术，一直发挥着重要的作用。学者们通过地物反射辐射特征变化来模拟地表覆盖类型、地表温度、土壤热惯量及地表蒸散发与土壤含水量的经验关系，从而实现土壤水分的反演[7-9]。然而，上述大多数方法都是经验性的，对于回归分析需要做出的一些统计假设可能会导致使用受限，如离群数据、非线性、异方差和多共线性，不能对土壤水分进行强有力的预测。机器学习方法能够发现数据之间有意义的关联、模式和规则，可以克服上述问题，如人工神经网络[10]、卷积神经网络[11]、支持向量机[12]、极限学习机(Extreme learning machine，ELM)[13]、随机森林(Random forest，RF)[14]，已经成功地应用于干旱半干旱地区土壤水分反演。研究证明基于机器学习方法能够构建综合多特征变量的土壤水分反演模型，可以提高土壤水分反演的时空分辨率以及反演精度[15-17]。研究表明，偏最小二乘回归法、极限学习机和随机森林是定量反演的可行方法[18-19]。总体而言，上述模型方法各具优缺点，对不同研究区域有其不同的要素要求和适用范围，有的土壤水分反演方法虽精度较高，但所需测定要素容量大，有的模型虽解决了大范围应用的难题，却在参数获取和模型建立方面要求很高。

本文以内蒙古自治区乌审旗作为研究区，对不同环境变量与表层土壤含水率进行相关分析，利用最优子集筛选环境变量作为模型输入变量，利用偏最小二乘(Partial least squares regression，PLSR)、极限学习机和随机森林等方法构建不同的土壤含水率反演模型，通过比较，探索毛乌素沙地腹部土壤含水率的高效反演模型及方法，在此基础上，反演研究区不同月份土壤水分并分析其时空变化。

1 数据与方法

1.1 研究区概况

乌审旗位于鄂尔多斯市西南部(图1)，地处毛乌素沙地腹部(37°38′～39°23′N、108°17′～109°40′E)，地势由西北向东南倾斜，平均海拔1 305 m。年平均降雨量350～400 mm，年平均气温6.8℃。主要土壤类型有栗钙土、草甸土、盐碱土、沼泽潜育土以及各类风沙土。主要土地利用类型为草地、沙地、林地、耕地、水体和建筑用地。

图1 研究区地理位置与采样区分布图Fig.1 Geographical location and sampling layout of study area

1.2 数据

1.2.1遥感数据

选取2016年4月21日、8月27日的Landsat8 OLI影像，数据来源于美国地质调查局网站(https:∥earthexplorer.usgs.gov)，通过ENVI 5.3软件完成影像预处理，具体包括：辐射定标、大气校正、影像镶嵌、裁剪和波段运算。从而进一步得到地表温度、地表反照率、缨帽变换要素、反射率、植被指数、水体指数、建筑指数和干旱指数等环境变量。其中植被指数和地表温度是表述地表特征的两个重要参数，以及二者构成的温度植被干旱指数均与地表土壤水分存在着紧密关系，是土壤水分监测的常用参数。反照率是影响地表辐射平衡的一个重要参数，其可以进一步影响土壤水分的变化。不同的光谱反射率对土壤水分的敏感性也不一样。缨帽变换要素的湿度对不同植被覆盖的土壤水分也有一定的表征。水体指数与建筑指数可以直接影响地表温度来间接影响土壤水分。哨兵1A数据是Level-1地距影像(Ground range detected，GRD)，成像方式为干涉宽幅 (Interferometric wide swath，IW)模式、极化方式(Vertical vertical(VV)和Vertical horizontal(VH))，经过辐射标定等得到所需要的后向散射系数。考虑地形条件对土壤水分的影响，本研究还通过地理空间数据云网站(http:∥www.gscloud.cn/)获取ASTER数字高程模型。

1.2.2野外实测数据

根据卫星过境时间在2016年4月20—22日、8月26—28日进行了两次地面采样实验。采样区域共24个，每个采样区域之间的最小间隔为1 km。采样时，首先记录每个采样区域5个点的GPS位置信息，各点间隔30 m，然后用土钻从土壤表面0～10 cm处分别采集土壤样品。最后，在现场对这5个点的土壤样品进行混合和称量。使用烘干称量法进行土壤含水率的测量，计算公式为

(1)

式中W——土壤质量含水率，%

W1——空土盒质量，g

W2——土盒和湿土质量，g

W3——土盒和干土质量，g

1.3 最优子集筛选

全子集回归是对所有预测变量的可能组合模型都进行拟合，然后根据贝叶斯信息准则(Bayesian information criterion，BIC)筛选出现有变量条件下的最佳模型，又叫最优子集筛选，计算公式为

BIC=mInc-2InL

(2)

式中m——模型参数个数

c——样本数量L——最大似然函数

1.4 建模方法及模型评价指标

偏最小二乘回归利用对系统中的数据信息进行分解和筛选的方式，提取对因变量的解释性最强的综合变量，辨识系统中的信息与噪声，从而更好地克服变量多重相关性在系统建模中的不良作用[20]。极限学习机是在单隐层前馈神经网络上发展起来的一种机器学习方法，可以随机初始化输入权重和偏置并得到相应的输出权重，弥补了传统神经网络中运行时间长等缺点[21]。随机森林是多重决策树的组合，输出类别由个别树输出的类别众数而决定的分类回归模型。不需要关于响应协变量关系的分布假设，该过程通过平均决策树进行统计上的可靠估算，降低过拟合风险[22]。

采用决定系数R2、均方根误差(RMSE)和平均绝对误差(MAE)验证反演模型的精度。RMSE、MAE越小，R2越大，表明模型反演精度越高。

2 结果与分析

2.1 描述性统计

图2 土壤含水率描述性统计Fig.2 Descriptive statistics of soil moisture content

为了对预测模型进行训练和验证，分别将2016年4月21日、8月27日土壤含水率数据(120个)随机分为两组：建模集(70%，84个样点)对模型进行训练，验证集(30%，36个样点)对模型进行检验。全集、建模集和验证集描述性统计结果如图2所示。4月，全集土壤含水率均值为4.84%，标准差为5.18%。建模集(0.45%～20.13%)和验证集(0.61%～15.89%)均值分别为4.73%和5.90%；8月，全集土壤含水率均值为7.51%，标准差为6.98%。建模集(0.82%～24.67%)和验证集(1.19%～24.00%)均值分别为7.20%和8.22%。表明建模集和验证集与全集土壤含水率保持类似的统计分布，在确保代表性样本的同时尽可能缩减了建模集和验证集中存在偏差的估计。

2.2 环境变量与表层土壤含水率的相关性分析

基于机器学习环境变量的数据来源与性质，包括17个变量：后向散射系数(σVV、σVH)、地表温度(Land surface temperature，LST)、地表反照率(Albedo)、亮度(Bright)、绿度(Green)、湿度(Wet)、波段反射率(红光、近红外、短波红外)、植被指数(Normalized difference vegetation index,NDVI)、水体指数(Normalized differential water index,NDWI)、建筑指数(Normalized differential building index，NDBI)、高程(DEM)、坡度(Slope)和温度植被干旱指数(Temperature vegetation drought index，TVDI)等要素(图3)。

图3 环境变量空间分布Fig.3 Surface biophysical characteristics and topographic parameters

通过对不同环境变量与表层土壤含水率进行相关分析，结果发现不同的环境变量与表层土壤含水率之间的相关程度不同(图4中*、**分别表示P<0.05、P<0.01)。首先，温度植被干旱指数与土壤含水率相关性最高，决定系数为0.64，主要原因是温度植被干旱指数综合植被和地表温度信息，从裸地到植被全覆盖可以准确体现土壤湿度状况[23]。其次为地表温度与土壤含水率相关性，决定系数为0.6，地表温度在土壤湿度和地、气相互循环过程中有重要作用，同时能间接反映土壤水分状况，表征旱情分布，因而对地表温度数据的研究是不可或缺的内容[24]，本研究地表温度与土壤含水率为高相关性，对土壤含水率高精准反演有重要价值。最后，σVV、植被指数、波段7反射率、σVH、波段6反射率、建筑指数、水体指数、反照率、亮度、绿度和湿度与土壤含水率之间的决定系数分别为0.38、0.38、0.35、0.32、0.30、0.18、0.13、0.26、0.14、0.15、0.06，其中NDVI与土壤含水率的相关性从4月到9月增加最明显，主要原因是NDVI在植被覆盖度高的月份更加敏感，反映的信息更为丰富；NDBI与土壤含水率的相关性从4月到9月减小最明显，主要原因是NDBI在植被覆盖度低的月份更加敏感。σVV、σVH在植被覆盖度低的月份对土壤含水率更加敏感。地形因素中DEM和坡度的正向相关性较好，这主要是地形会影响地表径流与该地所受的太阳光照、植被生长环境等[25]，从而对土壤水分产生影响，此外乌审旗的地形呈现西北高东南低的特征，这对DEM和坡度间相关性影响较大。反射率与土壤含水率也有很好的相关性，主要是由于土壤水分的变化会引起土壤颗粒物理性质和反射率也发生变化，通常随着土壤水分的减少，其反射率也会相应增大[25]；而植被的生长状态能够间接反映出土壤水分，通常来说，植被密度越大，生长态势越好的地方的土壤含水率越高。其中，波段反射率中B6、B7波段反射率与表层土壤含水率的相关性高于B4、B5波段，主要原因是短波红外光谱域包括水分吸收带，影响植被和土壤中水分含量的反射率敏感性[26]。

图4 不同月份环境变量与土壤含水率的相关性分析Fig.4 Correlation analysis of environmental variables and soil moisture in different months

图5 4月土壤含水率实测值与预测值比较Fig.5 Comparison of measured and predicted soil moisture content in April

总体而言，2个月间表层土壤含水率与后向散射系数均值R2(0.35)、波段反射率均值R2(0.26)、地表温度、地表反照率、缨帽变换要素、植被指数、水体指数、建筑指数均值的R2(0.24)以及干旱指数均值的R2(0.6)高于表层土壤含水率与地形参数的均值R2(0.19)。这一结果表明，表层土壤含水率对研究区后向散射系数、干旱指数、波段反射率和地表生物物理特征的依赖性较高。

2.3 最优子集筛选特征变量

基于BIC最小的原则，利用最优子集筛选法(Best subset selection，BSS)选取不同月份最优的环境变量组合，如表1所示。

表1 全子集筛选结果统计Tab.1 Statistics of full subset selection

2.4 土壤水分反演模型的精度评价与分析

图6 8月土壤含水率实测值与预测值比较Fig.6 Comparisons of measured and predicted soil moisture content in August

为了验证3类模型的可靠性和适用性，使用2期验证集中的36个样点实测数据分别进行精度检验。图5、6分别为2016年4月21日和2016年8月27日PLSR、ELM、RF 3种模型预测值与实测数据的散点图。由图5、6可知，预测值与实测值都具有很高的相关性。由表2可知，4月，建模集中3种模型方法的精度差异较大，相比于未筛选变量所构建的模型，利用BSS筛选变量所构建的模型精度均有所提升，且在RF模型的R2最高，为0.88，而RMSE和MAE均为最小值，分别为8.25%、5.17%；PLSR模型的R2最低，为0.66，RMSE和MAE分别为10.07%、9.43%，均为建模集中最大值；ELM模型性能介于二者间。由表3可知，8月，建模集中3种模型方法的精度差异较大，相比于未筛选变量所构建的模型，利用BSS筛选变量所构建的模型精度均有所提升，且在RF模型的R2最高，为0.89，而RMSE和MAE均为最小值，分别为7.36%、4.61%；PLSR模型的R2最低，为0.69，RMSE和MAE分别为9.74%、9.37%，均为建模集中最大值；ELM模型性能介于二者间。综上所述，建模集中，RF模型性能明显优于PLSR和ELM方法，反演效果由优到劣依次为RF、ELM、PLSR。在验证集中，3种模型的各指标相较于建模集均无明显下降，表明模型比较稳定，但三者相比较而言，仍然是RF模型的各项指标明显优于ELM和PLSR，而ELM和PLSR相比，ELM模型的R2明显高于PLSR模型，而RMSE和MAE低于PLSR模型。综合考察各模型建模集与验证集的评价指标，3种模型的预测性能和稳定性从高到低排序依次为RF、ELM、PLSR，基于PLSR的模型精度最低，这是由于土壤的组分非常复杂，功能团多样，偏最小二乘回归方法仅仅将环境变量与土壤含水率进行线性回归，同时没有将部分与土壤含水率非线性相关的关系考虑进来，导致所构建的模型具有一定的缺陷。而ELM、RF模型充分考虑了环境变量对表层土壤含水率的影响，并且具有对非线性问题的强解析能力和模型的稳健性。研究区表层土壤含水率模型预测精度结果进一步表明了RF模型自身的优越性和利用环境变量预测研究区表层土壤水分空间分布的可行性。

表2 4月土壤含水率反演模型与精度Tab.2 Soil moisture content retrieval models and accuracies in April

图7 4月表层土壤水分空间分布特征Fig.7 Spatial distribution characteristics of surface soil moisture in April

表3 8月土壤含水率反演模型与精度Tab.3 Soil moisture content retrieval models and accuracies in August

2.5 土壤水分的空间分布特征

图8 8月表层土壤水分空间分布特征Fig.8 Spatial distribution characteristics of surface soil moisture in August

由图7、8可知，不同月份ELM和RF模型预测的土壤水分分布图与PLSR模型预测的土壤水分分布图比较相似。高土壤含水率出现在研究区的北部和东南部，主要由于这些地区有密集的植被覆盖，土壤持水性好。中北部平坦地区的土壤含水率较低，主要由于这些地区植被覆盖度低，地表蒸散强烈。此外， ELM反演的土壤水分分布(4月，1.06%～22.47%；8月，2.48%～30.57%)与RF模型(4月，1.56%～23.71%；8月，2.63%～34.93%)相似。而PLSR模型(4月，0.72%～22.03%；8月，2.16%～29.92%)的土壤水分分布图更强烈地表现出所有部分的低土壤含水率。一般情况下，PLSR方法反演土壤水分的空间变异性低于ELM和RF方法，其主要原因是ELM和RF方法更能识别和反演土壤水分中的局部微小变化，能够更详细地表达地表信息和空间异质性。

3 讨论

土壤水分是了解地表过程、陆-气相互作用、干旱预测、作物生长模式等的基础，是一个动态变量，受多种因素的影响，如植被覆盖、地表粗糙度、土壤类型、地形等，在不同的时空尺度上，即使在较小的区域内也会发生显著的变化。本文在建模方法上选择了PLSR、ELM和RF构建研究区土壤含水率反演模型，经过对比分析发现，在相同条件下ELM、RF算法效果优于PLSR算法，因为土壤含水率与环境变量之间并非简单的线性关系，PLSR模型在处理土壤含水率与环境之间复杂关系时具有一定的局限性[27]，而ELM、RF算法在非线性问题中具有较强的解析力和较高的模型鲁棒性，这与葛翔宇等[18]和蔡亮红等[13]研究结果一致。但在ELM、RF土壤含水率反演模型中，本文建模集和验证集R2均不小于0.65，其中RF模型的精度最高，这与王浩等[28]利用RF方法考虑温度、蒸散发、地形等因子对土壤水分影响的结果相一致。土壤水分不仅与这3个变量有关，还与坡度、波段反射率、地表反照率和植被指数等其他变量有关[29]。本研究与传统的经验模型相比，基于机器学习的检索算法避免了复杂的公式，提高了土壤水分的检索效率[30]。但各环境变量对研究区土壤水分的影响不同，在红色、近红外、短波红外1和短波红外2波长范围内，土壤含水率和土壤光谱反射率之间存在不显著的线性关系，主要由于水在非饱和砂中的水力特性，土壤水分低时光谱反射率的降低非线性[26]，本研究发现短波红外1和短波红外2波长与土壤水分之间线性相关。主要原因是短波红外波段裸土和植被土壤的反射率与表面粗糙度的相关性最好。但是光谱反射率会受到土壤固有因素的影响，例如有机质含量、粒径分布、矿物成分、表面粗糙度和土壤元素的颜色[31]。温度植被干旱指数对土壤水分的影响最大，它可以将植被指数和地表温度结合起来，从水分胁迫开始到植被指数出现变化的时刻之间的时间延迟是最小的，从而避免了红波段和近红外波段光谱响应的延迟[32]，这一结果与王思楠等[33]利用不同干旱指数研究土壤水分的结果较为相似。地表温度对土壤水分的影响也很大，因为它直接影响能量平衡分量和地表蒸散发。随着地表温度的增大，土壤水分含量减小。BABAEIAN等[34]和CARLSON等[35]也表明，影响土壤水分空间变异性的最重要参数是地表温度。因此，在基于遥感数据常用的TVDI[36]、条件温度植被指数(Vegetation temperature condition index，VTCI)[37]、植被供水指数(Vegetation supply water index，VSWI)[38]等反演土壤水分时，都使用了地表温度参数。因此，有研究表明，土壤水分也受到这些环境变量因素的直接和间接影响。在环境变量中，研究区NDVI和绿度的增加表明植被覆盖度增加，蒸腾作用增加，因此土壤水分增加。增加NDBI和亮度会增加地表反照率，从而降低土壤水分。地形参数也是一个控制地表温度的重要影响因素，高程的增加导致气温下降，地表蒸散发减少，从而土壤水分增加。在该区域，表面生物物理性质对土壤水分的影响大于地形的影响。这些结果与其他一些研究结果相似[30]。

土壤水分受不同环境因素的显著影响，不同模型在不同月份的反演结果相似。在毛乌素沙地腹部，4月降水量变化小，而且农作物处于生长初期，蒸发量小，土壤水分变化量不大；8月降水高度集中，沙壤土入渗率高，降水很快能渗入土壤，故很少能形成径流，土壤水分出现明显增加。土地利用在确定土壤水分变化的空间变异性方面非常重要，因为它影响植被覆盖、入渗和径流速率、蒸散过程、土壤表面特征。土地利用甚至可以消除地形相关参数对土壤含水率的影响。一些研究提供了通过反射图像和辅助地理空间数据估计土壤湿度的解决方案。乌审旗南部地区的高土壤水分可能是由于土壤表面的蒸发潜力相对较低，降水多，植被对降水的截留作用比较大。北部地区土壤含水率低可能是降水稀少，植被覆盖率低，基本上为沙地，水分下渗非常快。当沙丘被梭梭稳定后，粘土和淤泥的数量迅速增加，因为悬浮颗粒堆积和细颗粒是由沙子的风化机制产生的。粘土和淤泥含量的增加也与土壤含水率有轻微的正相关关系[39]。沙地的表层土壤含水率低是因为由于高太阳辐射，土壤表面蒸发潜力高，植物的生长主要受土壤水分的控制，而土壤水分是光合作用所必需的。此外，土壤湿度调节植物蒸腾和蒸发的速率，进而影响近地表温度、湿度和大气水蒸气[40]。本研究由于其他土壤特性(例如土壤孔隙度、堆积密度和土壤有机质含量)很难获取，仅使用了17个环境变量，可能会导致反演结果的不确定性。除此之外，光学遥感受天气的影响较大，本研究选取的遥感影像获取时间范围内研究区处于晴空状态下，因此未考虑降雨条件下的不同模型的土壤水分估算表现。

4 结论

(1)通过对不同环境变量与表层土壤含水率进行相关分析，发现不同的环境变量与表层土壤含水率之间的相关程度不同。环境变量温度植被干旱指数和地表温度均与土壤含水率相关性较高。

(2)对比PLSR、ELM和RF 3种模型评价指标发现，4月，RF模型的R2为0.74，RMSE为8.85%；8月，RF模型的R2为0.75，RMSE为8.66%，表明RF为土壤水分反演最优模型。

(3)不同月份ELM和RF模型预测的土壤水分分布图与PLSR模型预测的土壤水分分布图比较相似。高土壤含水率出现在研究区的北部和东南部，中北部平坦地区的土壤含水率较低。此外， ELM与RF模型反演的土壤水分分布，突出了研究区地理位置的高低土壤含水率。而PLSR模型更强烈地表现出所有的低土壤含水率。