姚一飞 王 爽 张珺锐 黄小鱼 陈 策 张智韬
(1.西北农林科技大学水利与建筑工程学院， 陕西杨凌 712100；
2.西北农林科技大学旱区农业水土工程教育部重点实验室， 陕西杨凌 712100)
土壤水分是决定植被生长、农作物生产、区域生态环境稳定的关键因素，其动态变化趋势影响着各种生物的生存和发展。因此，准确、及时地反演土壤水分具有重要意义。遥感监测土壤含水率技术克服了传统田间实测法低效率、高成本、单点测量的缺点[1]，是大范围监测土壤含水率时空分布和动态变化的有效手段[2]，在精准农业中具有广阔应用前景[3]。
利用光学遥感监测裸土土壤含水率的可行性已被诸多学者验证[4-6]。但在现实情况下农田土壤常被不同作物覆盖，仅监测裸土土壤含水率无法满足农田灌溉的需要。因此，部分学者提出采用构造不同光谱指数可以定量监测植被生长状况及其覆盖下的土壤含水率。文献[7]利用地表温度与归一化植被指数(NDVI)之间的比例关系对地表的干旱程度开展动态监测；
土壤背景对用于计算植被指数的反射率具有不同程度的影响，HUETE[8]基于NDVI大量观测数据提出了土壤调节植被指数(SAVI)来减少土壤背景的影响[9]，而基于SAVI发展了多种新的指数，包括：土壤调节植被指数2(SAVI2)[10]、转换型土壤调节植被指数(TSAVI)、调整转换型土壤调节植被指数(ATSAVI)[11]等。夏权等[12]利用资源三号卫星的多光谱数据，建立增强型植被指数(EVI)反演土壤水分含量数学模型；
吴春雷等[13]针对垂直干旱指数(PDI)在农田植被区土壤湿度监测精度降低问题引入垂直植被指数(PVI)进行调整；
蔡亮红等[14]利用Landsat8 OLI遥感数据，将短波红外波段引入到可见光-近红外植被指数中，获得了较好的土壤水分反演结果；
WEST等[15]研究发现干旱条件下NDVI能较好反映土壤水分状况，并且Sentinel-2卫星反演效果优于Landsat8，在土壤水分监测中应用潜力大；
同样利用Landsat卫星，NATALIE等[16]使用归一化红外指数(NDII)和NDVI分析湿地土壤水分的时空趋势，REZA等[17]使用归一化差值含水指数(NDMI)对干旱半干旱地区的土壤表面水分进行了测量。GF-1作为我国重要的多光谱卫星，在土地利用分类[18]、船舶探测[19]、农作物光合有效辐射吸收比率估算[20]、气溶胶光学厚度反演[21]、冰川识别[22]方面已开展广泛研究，但目前用于植被覆盖条件下土壤水分的反演研究较少。
本文以河套灌区解放闸灌域沙壕渠为研究区，以GF-1卫星遥感影像为数据源，并同步采集不同深度土壤含水率，利用全子集筛选法确定不同深度下光谱指数的最优自变量组合，采用多元线性回归、BP神经网络和支持向量机3种算法，构建不同深度下土壤含水率反演模型，探究GF-1号卫星反演植被覆盖条件下农田土壤含水率的可行性，以期为全灌区干旱监测与预警提供技术依据。
1.1 研究区概况
内蒙古河套灌区是我国设计灌溉面积最大的灌区[23]。位于内蒙古自治区西部的巴彦淖尔盟，东西长约250 km，南北宽约50 km，总灌溉面积约为5 740 km2。灌区是典型的温带大陆性气候，降雨量少、蒸发量大，地处干旱的西北高原，但其区域内遍布人工水系，年引黄河水量高达50亿m3，形成了独特的干旱区“绿洲”生态，是国家和自治区重要的商品粮、油生产基地。保障河套灌区农业可持续发展对我国旱区农业生产生活意义重大，随节水改造工程的实施，河套灌区面临指令性的节水形势。
解放闸灌域是河套灌区的第二大灌域，位于河套灌区西部(图1a)，灌溉面积1 420 km2，灌溉水主要来源于过境的黄河水，年均降水量约为140 mm，年均蒸发量约为2 000 mm。解放闸灌域地形较为平缓，地势由东南向西微度倾斜[24]，平均坡度为0.2%，其平坦的地形地势使排水系统能力降低，加之地下水埋深较浅，为高效地利用水资源，优化灌溉制度，有必要进行大面积土壤含水率的监测。根据实际试验条件，本次选取河套灌区解放闸灌域沙壕渠试验站为研究区。
1.2 土壤样本采集与处理
7—9月是研究解放闸灌域植被覆盖区域的最佳时期，该时段内灌域的主要作物葵花、玉米(占所有作物面积80%)等都处于生长旺盛期，9月下旬灌域作物均已成熟，等待收割[25]，因此选择作物生长旺盛、植被信息丰富的2019年7月15—25日作为研究时段。采用手持GPS定位，利用打土钻取样，且记录采样点的植株高度、耕作状况、植被种类等，综合考虑了不同土地利用方式、种植结构、灌溉方式等情况，以及GF-1号卫星分辨率为16 m，在河套灌区沙壕渠共布置了120个16 m×16 m采样点，并尽量保证点的均匀性。采样点分布情况见图1b。在植被覆盖情况下，土壤水分胁迫将影响作物根系的生长，导致作物光谱特征及作物长势明显变化，卫星遥感监测是依据解译作物生长状况及作物长势进而判断根系主要活动层所在土壤深度的含水率，故采用五点取样法，根据研究区域该时段的作物主根系深度及相关研究成果[26-28]，在深度0～60 cm范围内每隔20 cm分层取样，将采取的土样充分混合后得到具有代表性的样品。
1.2.1土壤含水率测定
将野外采集到的各深度的120个土样快速装入已知质量的干燥铝盒中，贴好标签并密封后称量。采集工作完成后将封装好的土样托运回西北农林科技大学旱区农业水土工程教育部重点实验室进行后续处理，经105℃，24 h恒温处理后，取出铝盒称量。土壤质量含水率θm计算式为
(1)
式中W1——湿土加铝盒质量
W2——干土加铝盒质量
W3——空铝盒质量
结合地面植被情况，剔除小麦收割后、裸土以及存在明显粗误差的采样数据(由于试验测量过程中或记录相关数据时错读、错记造成含水率为负)，剩余100余个土壤样本用于本次7月植被覆盖条件下的灌区土壤含水率反演研究。
1.2.2土壤含水率统计特征
选取0～20 cm、20～40 cm、40～60 cm、0～40 cm和0～60 cm 5个深度土壤含水率为研究对象，其中0～40 cm和0～60 cm含水率为对应深度内实测数据的平均值。不同深度数据集按照土壤含水率降序或升序排列，采用每隔2个样本提取1个样本作为验证集的方法，使得建模集与验证集样本数比例近似为2∶1，以期实现建模样本和验证样本范围一致且分布均匀。土壤含水率统计特征见图2和表1。
表1 土壤含水率统计特征Tab.1 Statistical characteristics of soil moisture content
0～20 cm、20～40 cm和40～60 cm 3个深度土壤含水率平均值分别为24.46%、26.48%和29.21%，标准误差分别为0.27%、0.36%、0.28%，表明不同深度的土壤含水率存在显著差异，土壤含水率随着土壤深度而增加，与文献[29-30]研究结果一致。变异系数是表征土壤含水率变化幅度的指标，变异系数越小，土壤含水率越稳定，深层(40～60 cm)土壤含水率变异系数相对最小，主要原因是深层土壤受光照辐射、降雨量以及土壤蒸散发等外界环境因素的影响较小[31]。
各个深度下总体与建模集、验证集的土壤含水率变化范围和统计分布(平均值、中位数、标准差、变异系数)基本一致，确保了样本的代表性，同时也避免了建模和验证中的估计偏差。
1.3 遥感影像选取与处理
采用GF-1卫星遥感影像为数据源，空间分辨率16 m，共计4个波段，重访时间2 d，卫星影像的成像时间为2019年7月27日，与实测数据日期基本同步。卫星影像数据由中国资源卫星应用中心下载获得(http:∥www.cresda.com/CN/)。对下载的影像经过几何精校正、辐射校正、大气校正预处理后用于后续土壤含水率反演。卫星4个波段范围分别为0.45～0.52 μm、0.52～0.59 μm、0.63～0.69 μm、0.77～0.89 μm。
1.4 光谱指数提取与筛选
1.4.1光谱指数选取
选取空间分辨率为16 m×16 m的GF-1卫星遥感影像，仅用4个波段反射率提取土壤水分信息并分析其相关性有一定局限性，考虑到7月为灌区内作物生长发育的关键时期，植被覆盖率较高，可利用GF-1卫星数据的4个波段反射率计算相关光谱指数建立遥感图像与土壤含水率的定量关系。
本文共选取包括土壤亮度指数(BI)[32]、盐分指数(SI)[33]、盐分指数1(SI1)、盐分指数2(SI2)、盐分指数S1、盐分指数S2、盐分指数S6和EVI[34]、差值环境植被指数(DVI)[35]、PVI[36]、PDI[37]、SAVI2、土壤大气抗阻植被指数2(SARVI2)[38]、TSAVI、ATSAVI在内的光谱指数共15个。其中，EVI、DVI等是光学卫星反演土壤水分的常用植被指数，同时，由于研究区域为我国典型的盐渍化地区，与正常情况相比，当作物受到水分胁迫和盐分胁迫时，其状态的改变会在光谱特征上有所体现[39]。当研究区域内土壤水分在时空上具有一定动态变化时，土壤水分和盐分对土壤反射光谱有相似的影响[40]。故在此基础上，选取6个盐分指数共同用于本次植被覆盖条件下农田土壤含水率的反演。
1.4.2全子集筛选
全子集筛选法是遍历自变量全集中所有可能的组合方式，通过最小二乘法对每种变量组合进行拟合，并从中选取效果最佳的模型。主要计算过程如下：
1.5 回归模型构建
1.5.1多元线性回归模型
多元线性回归(Multiple linear regression，MLR)是一元线性回归的扩展，可分为单因变量和多因变量2种类型[48]。
当因变量受多个自变量影响时，可通过一般线性回归模型表示为
Y=βX+ε
(2)
式中X——自变量β——自变量系数
ε——满足正态独立同分布的随机误差项
与一元线性回归类似，参数β估计是在误差平方和最小的前提下，用普通最小二乘法求解[49]。利用Matlab软件进行多元线性回归的建模与验证。
1.5.2BP神经网络模型
BP神经网络(Back propagation neural network，BPNN)是目前应用较为广泛、相对通俗易懂的神经网络模型。该模型的非线性映射能力较强，能够根据输出结果与实测值之间的误差调节系统内部的参数，并依据精度要求反复进行信号前向传播和误差反向传播的操作，最终得到输出结果。该模型拓扑结构由输入层、隐藏层、输出层构成。输入层为不同深度下全变量(即4个波段反射率和15个光谱指数)及经过全子集筛选后的最优变量组合；
隐含层的网络层数设置为2，采用双曲正切激活函数；
输出层为土壤含水率，采用恒等激活函数；
目标误差为0.001。利用Matlab软件进行BP神经网络建模与验证。
1.5.3支持向量机模型
支持向量机(Support vector machines，SVM)是一种基于分类边界的方法[50]，可扩展用于回归分析。其能够根据有限的样本信息，平衡模型的复杂性和学习能力，找到最优解，在处理高维数据集、小样本和非线性问题中都具有良好表现。SVM在结构风险最小化原理的基础上，利用核函数映射的方法，通过非线性变换，将低维空间的点映射到高维空间，使它们成为线性可分，建立线性模型以此求解非线性问题，很大程度上克服了“离散值多”和“过学习”等问题。对建模精度影响最大的3个参数是核函数类型、惩罚参数C、核参量g[51]。
径向基核是应用最广泛的核函数，简单实用、普适性好，通过选择合适的参数，可以适用于任意分布的样本[52]。本次支持向量机采用RBF核函数进行计算，其解析式为
(3)
式中x——未知向量x′——支持向量
σ——函数宽度
利用Matlab软件进行支持向量机建模与验证。
1.5.4模型精度评价指标
2.1 土壤含水率与自变量相关性分析
选取GF-1卫星的4个波段反射率(B1、B2、B3、B4)和15个光谱指数，组成反演土壤含水率的自变量集。随机选取约2/3的样本数据组成估算数据集，对自变量集与土壤含水率的相关性进行分析，得到Pearson相关系数如表2所示。
根据相关系数检验临界值表对自变量集与土壤含水率的显著性进行检验。当自由度为70时，相关系数绝对值大于0.232即达到0.05显著性水平，大于0.302即达到0.01显著性水平。由表2可以看出，EVI、DVI、PVI、SARVI2在深度0～20 cm、20～40 cm、0～40 cm和0～60 cm下，与土壤含水率均未表现出明显的显著性，在深度40～60 cm 下与土壤含水率均达到0.05显著性水平；
B3在深度40～60 cm 下与土壤含水率未表现出明显的显著性，与其余深度达到0.05显著性水平；
SI、SI1在深度40～60 cm下和S6在深度0～20 cm下分别与土壤含水率达到0.05显著性，其余自变量与各深度土壤含水率均达到0.01显著性水平。
表2 波段及光谱指数与土壤含水率的Pearson相关系数Tab.2 Pearson correlation coefficient of bands and spectral indices with soil moisture content
2.2 基于全子集筛选的最优自变量组合
表3 基于全子集筛选的土壤含水率最佳自变量组合Tab.3 Best independent variable combination of soil moisture content based on total subset selection
选择SI、SI1、S1共3个自变量，0～60 cm处选择B1、SI1、S1共3个自变量，作为各深度下最优自变量组合方式。
2.3 土壤含水率反演模型建立与验证
2.3.1线性回归模型建立与验证
分别以不同深度下筛选前后的波段反射率及光谱指数为自变量，土壤含水率为因变量，运用多元线性回归模型进行不同深度下土壤含水率估算。其建模集和验证集精度如表4所示。
表4 基于多元线性回归全子集筛选前后的土壤含水率反演模型Tab.4 Inversion model of soil moisture content before and after screening of all subsets based on MLR
总体来说，在多元线性回归模型中，各深度筛选后模型拟合效果和泛化性能优于筛选前，且筛选后模型稳定性更强。全子集多元线性回归模型在不同深度土壤含水率的反演效果由好到差依次为：0～40 cm、0～60 cm、20～40 cm、0～20 cm、40～60 cm。
2.3.2BP神经网络模型建立与验证
分别以不同深度下筛选前后的波段反射率及光谱指数为自变量，土壤含水率为因变量，运用BP神经网络模型进行不同深度下土壤含水率估算。其建模集和验证集精度如表5所示。
表5 基于BP神经网络全子集筛选前后的土壤含水率反演模型Tab.5 Inversion model of soil moisture content before and after screening of all subsets based on BPNN
总体来说，在BP神经网络模型中，除深度20～40 cm外，其余各深度筛选后模型拟合效果和泛化性能优于筛选前，且筛选后模型稳定性更强。全子集BP神经网络模型在不同深度土壤含水率的反演效果由好到差依次为：0～40 cm、0～60 cm、20～40 cm、0～20 cm、40～60 cm。
2.3.3支持向量机模型建立与验证
分别以不同深度下筛选前后的波段反射率及光谱指数为自变量，土壤含水率为因变量，运用支持向量机模型进行不同深度下土壤含水率估算。其建模集和验证集精度如表6所示。
表6 基于支持向量机全子集筛选前后的土壤含水率反演模型Tab.6 Inversion model of soil moisture content before and after screening of all subsets based on SVM
总体来说，在支持向量机模型中，深度0～60 cm筛选后模型拟合效果略有下降，其余各深度筛选后模型拟合效果和泛化性能优于筛选前，过拟合现象得到一定改善，反演误差相当。全子集支持向量机模型在不同深度土壤含水率的反演效果由好到差依次为：0～60 cm、0～40 cm、0～20 cm、20～40 cm、40～60 cm。
2.4 反演效果综合评价
图3 土壤含水率反演深度综合对比图Fig.3 Comprehensive contrast diagrams of soil moisture content inversion depth
2.4.1反演深度评价
2.4.2反演模型评价
图4 土壤含水率反演模型综合对比Fig.4 Comprehensive contrast diagrams of soil moisture content inversion model
利用GF-1卫星4个波段和15个光谱指数，经过全子集筛选的方法确定不同土壤深度条件下敏感变量组合方式，个别模型和深度筛选后模型效果略有下降，其余模型的拟合效果和泛化能力均有提高，反演误差基本保持不变，过拟合现象有所减轻。3种模型模拟结果表明，筛选后的敏感变量组合方式仅用2～4个自变量即可得到筛选前同等精度的反演效果，且模型稳定性较大提高。罗湘华等[53]使用全子集有效地消除了共线性在土地利用变化统计分析中的影响，李长春等[44]在评估叶片表面指数时，也发现使用全子集筛选方法创建的回归模型与传统线性回归和机器学习模型相比更加贴合实际，符合全子集筛选有效的结论。
对比全子集筛选后的3种模型在不同深度下的反演结果，发现在植被覆盖条件下利用GF-1号卫星数据反演深度0～40 cm的土壤含水率效果最佳，深度0～60 cm次之。深层土壤(40～60 cm)反演效果较差，可能是由于深层含水率的空间变异性相对更小，不同样本点之间的光谱差异更多由浅层土壤造成，一定程度上掩盖了深层(40～60 cm)土壤含水率的光谱差异，并且深层(40～60 cm)土壤平均含水率更高，对植被的胁迫性也相对较小，所反映出来的差异也更小。张智韬等[54]利用TM5、TM7数据，应用归一化土壤湿度指数监测土壤含水率也得出反演效果最好的土层深度为0～40 cm的结论，屈创[55]基于多源遥感数据的白龙江流域水分反演研究里也表明回归模型对深度0～40 cm范围的土壤水分反演精度最高。其主要原因是：①表层土壤水分容易受到太阳辐射、降水、风速等多种因素的影响，并且卫星成像时间和野外实测土壤含水率数据很难达到完全同步，下层土壤处于一个被保护的状态，受到外界环境的影响较小，其土壤水分的变化在一段时间内能维持在一个相对稳定的状态。②作物耗水主要集中在0～40 cm的土层[56]，作物生长期根系对土壤水分反应敏感，受到土壤水分胁迫更加明显。在本次研究的河套灌区中，7月下旬的主要作物葵花处于苗期后期，生理需水量处于中等水平，主要吸水层集中在距地表0～50 cm土层之间[57]，最大吸水峰值在30 cm土层之间，与本文研究结果基本相符。
本研究表明利用GF-1卫星反演植被覆盖条件下农田土壤含水率具有可行性，评价了植被覆盖条件下GF-1卫星对不同深度土壤含水率的反演能力，研究不同深度下最优的土壤含水率反演模型，为全灌区干旱监测与预警提供了一定技术依据。CHEN等[58]率先对GF-1宽视场(WFV)传感器土壤水分反演观测能力做出了定量评价，NIE等[59]在进行光谱指数在土壤水分反演中的适宜性评价时也得出使用GF-1影像获得的整体反演结果优于使用Landsat8 OLI影像的结论；
GF-1光学数据与SAR数据结合能有效提高土壤水分的反演精度，WANG等[60]基于GF-1卫星计算的植被指数与Sentinel-1 SAR数据，结合水云模型(WMC)和高级积分方程模型(AIEM)反演土壤水分的R2为0.52～0.56，RMSE为0.052%～0.064%；
TAO等[61]使用多时相SAR数据和GF-1数据提出了改进的植被反向散射模型，结果表明该模型对于土壤水分估算准确。本文通过全子集筛选法确定了不同土壤深度的最优光谱指数组合，得到了与以上研究同等精度的土壤水分反演结果。
解放闸灌域盐渍化问题突出，地下水位偏高，会对土壤水分反演造成一定程度上的影响[62]，本文中采用的7个土壤指数(包括6个盐分指数)在以往多用于土壤含盐量的研究，首次将其利用在植被覆盖条件下土壤含水率的反演，结果表明这些盐分指数在各个深度下的土壤含水率反演模型中都表现良好，为未来在盐渍化地区的土壤含水率反演中提供了新的思路，未来可以进一步研究分析水盐之间的相互作用。此外，采用更高空间分辨率卫星遥感数据，考虑不同地区农田的种植结构差异，以及不同作物不同生育期根系影响土壤水分特点，本文结论是否具有更广泛的适用性还有待进一步研究。
(1)利用全子集筛选法确定了不同土壤深度下光谱指数的最优自变量组合，通过对比筛选前后各模型的反演效果，一定程度上证明了利用GF-1卫星反演植被覆盖条件下农田土壤含水率的可行性，结果表明筛选后模型反演效果有所提升，能减轻过拟合现象，全子集筛选有效。