基于全子集-分位数回归的土壤含盐量反演研究

张智韬 韩 佳 王新涛 陈皓锐 魏广飞 姚志华

(1.西北农林科技大学旱区农业水土工程教育部重点实验室， 陕西杨凌 712100； 2.西北农林科技大学水利与建筑工程学院， 陕西杨凌 712100； 3.中国水利水电科学研究院水利研究所， 北京 100038)

0 引言

土壤盐渍化程度是影响灌区土地资源经济效益和作物产量的关键因素，及时、准确地监测盐渍化动态是进行科学防治的重要前提。与传统的土壤盐渍化监测方法相比，卫星遥感具有无创、动态、综合、快速、宏观等优点，已成为现阶段大面积监测土壤盐渍化的重要研究方向。

国内外学者在利用卫星遥感监测土壤盐渍化程度方面取得了一定的进展。目前，常用方法主要有遥感图像分类及解译法[1]、光谱指数法[2-3]、特征空间分析法[4]。其中，光谱指数法是一种便捷高效的方法。陈红艳等[2]在传统的植被指数基础上加入短波红外波段，明显提高了遥感反演土壤盐渍化的精度。EL HARTI等[3]在盐分指数(SI)的基础上加入蓝波段，构建了新的盐分指数OLI-SR，提高了Morocco灌区盐渍化反演的精度。但是，光谱指数易受到下垫面条件和大气状况等多种因素的影响，具有明显的地域性和时效性。因此，筛选敏感光谱指数对于土壤含盐量监测具有重要意义。目前，常用的变量筛选方法有灰度关联法[5]、岭回归[6]、LASSO回归[7]、变量投影重要性分析[8]等，这些筛选方法均为局部最优筛选。全子集筛选法列举全部可能组合方式，建立全局最优模型，以包含最少自由变量的模型解释因变量，进而消除共线性的影响。李长春等[9]在估测叶面积指数时发现，采用全子集筛选法建立的回归模型优于PLSR、SVM和RF模型。冯克伟[10]经过全子集筛选，对野生二粒小麦和硬粒小麦的耐盐指数分别拟合，预测决定系数达到0.95以上，进一步说明全子集筛选法对于筛选光谱指数具有一定的可行性。

在以往的土壤盐渍化监测中，大多采用经典线性回归模型[11-12](如最小二乘法和偏最小二乘法)和机器学习模型[13-14]。经典线性回归对随机误差的分布特征要求严格，对于一些实际问题，难以得到无偏、有效的参数估计值[15]。人工神经网络模型和支持向量机模型是常用的机器学习算法。人工神经网络模型基于梯度学习方法，在层间传播的非参数非线性模型，因其具有较强的自学习和自适应能力、非线性映射能力和容错能力[16]，而被广泛用于遥感反演中。支持向量机模型主要用于对高维度、小样本的分类和回归，将有限样本信息在模型的复杂性和学习能力之间寻求最佳折中，在保证泛化能力的前提下达到最优学习效果[17]。与机器学习方法相比，分位数回归模型所需参数少，训练时间短，不易产生过拟合现象。分位数回归可以在任意一个分位点准确地描述自变量对于因变量的变化范围以及条件分布形状的影响，适合处理波动性大、异质性强的数据[18]，广泛应用于经济、医学等领域，但在监测盐渍化方面的研究还很少[19]。

本文利用河套灌区解放闸灌域GF-1卫星遥感影像数据，计算多种光谱指数，将通过全变量及全子集筛选法得到的最优光谱指数组合方式作为自变量，分别构建不同深度下的人工神经网络、支持向量机和分位数回归3种模型，进而对比分析得到土壤含盐量反演的最优模型，以期提高植被覆盖条件下的土壤含盐量反演精度。

1 材料与方法

1.1 研究区概况

研究区设在内蒙古自治区河套灌区解放闸灌域(40°30′～41°17′N，106°33′～107°31′E)，是河套灌区的第二大灌域，位于河套灌区西部，西接一干灌域，东邻永济灌域，灌溉面积1 420 km2，灌溉水主要来源于过境的黄河水。主要农作物包括玉米、小麦、葵花。该灌域地处干旱半干旱气候区，年均降水量约为140 mm，年均蒸发量约为2 000 mm[20]。研究区内平均坡度为0.2%，地势平坦使得排水系统的排水能力降低，地下水埋深较浅。加之降雨少、蒸发大、高灌低排，使得研究区内土壤盐渍化问题较为严重。

1.2 遥感数据获取及预处理

以国产高分一号卫星影像(GF-1 WFV相机)为数据源，通过中国资源卫星应用中心进行下载(www.cresda.com)。卫星影像的成像时间为2018年7月27日，与实测数据日期基本同步，高分一号卫星数据的重访周期为4 d，空间分辨率为16 m，包括4个波段，分别为蓝(0.45～0.52 μm)、绿(0.52～0.59 μm)、红(0.63～0.69 μm)和近红外波段(0.77～0.89 μm)。对下载的影像经过几何精校正、辐射校正、大气校正、镶嵌、裁剪预处理过程。并对影像进行低通滤波处理[21]，保存图像的低频部分，使图像平滑，在一定程度上消除了高频随机噪声，提高数据的信噪比。

1.3 土壤样本采集及化学分析

采样点的选定考虑到灌域内盐分分布、植被覆盖种类(包括葵花35个、玉米22个、小麦8个、果蔬8个和荒地7个)，并顾及点位分布的均匀性，本次试验在解放闸灌域共设置了80个不同土壤盐渍化程度的采样点，采样点分布情况见图1。采样时间为2018年7月12—16日，根据该时期灌区内主要作物的根系活动层所在深度，每个采样点按照0～20 cm、20～40 cm、40～60 cm分层采样，并通过手持GPS记录采样点位置信息及周围环境信息。采用五点法采样，采样单元为16 m×16 m，每个采样点取5个土样，编号后带回实验室。

将野外收集的土样干燥研磨处理后，配置土水比为1∶5的土壤溶液，经搅拌、静置、沉淀、过滤后，采用电导率仪(DDS-307A型，上海佑科仪器公司)测定土壤溶液电导率EC1∶5，对每个采样点的5个土样电导率取平均值作为该样本样点处的电导率，并通过经验公式计算土壤含盐量S(S=(0.288 2EC1∶5+0.018 3)×100%，剔除图像中小麦收割后的8个样本，剩余72个土壤样本用于本次植被覆盖条件下的灌区土壤含盐量反演。将测得的土壤含盐量从大到小排序，每隔2个样本取出1个作为验证集样本。可保证建模样本和验证样本范围一致且分布均匀。土壤含盐量数据的统计特征如表1所示。土壤盐渍化等级划分参照文献[22]。表1中0～40 cm土壤含盐量为0～20 cm和20～40 cm深度实测含盐量的平均值；0～60 cm土壤含盐量为0～20 cm、20～40 cm、40～60 cm深度实测含盐量的平均值。

图1 研究区地理示意图及采样点分布Fig.1 Schematic of study area and distribution of sampling points

表1 采样点的土壤含盐量统计Tab.1 Summary statistics of soil salinity of sampling points

1.4 光谱指数计算及全子集筛选

高分一号遥感影像自身的光谱分辨率不高，对于植被覆盖条件下的土壤盐渍化反演模型，仅用多个波段的反射率对比分析提取土壤盐渍化信息有明显的局限性，本研究利用高分一号卫星数据的4个波段反射率和12个光谱指数(包括5个土壤指数、3个盐分指数和4个植被指数，见表2)建立遥感图像与土壤含盐量的定量关系。

表2 光谱指数Tab.2 Spectral index summary

注：RB、RG、RR、和RNIR分别为高分一号卫星数据对应4个波段的光谱反射率。

全子集筛选是基于不同自变量的所有可能的组合方式，对缩减后的变量组合通过最小二乘法进行拟合，在所有可能的模型中选择一个最优模型。该方法简单直观，适合自变量较少的条件下使用。其主要计算步骤如下：记K为子变量数目(K=1,2，…，P)；拟合1～P个预测变量的模型；在1～P个模型中，根据“调整后R2最大”准则来选择P个最优模型；根据验证集调整后决定系数(Coefficient of determination，R2)、均方根误差(Root mean squared error，RMSE)从P个模型中选择一个最佳自变量组合。全子集筛选的模型建立和预测通过Matlab R2017b完成。

1.5 模型建立

1.5.1分位数回归模型

分位数回归研究自变量与因变量的条件分位数之间的关系，通过最小化离差绝对值的加权和，依据因变量的条件分位数对自变量进行回归，进而得到所有分位数下的回归模型[27]。相较于传统的回归方法，其优点主要表现在：①无需对模型中变量进行随机扰动或进行正态变化，在非正态的干扰下，分位数模型精度更高。②模型中异常点对于模型整体精度影响较小，模型稳定性强。③可以给出任意分位点的参数估计，解释自变量对不同水平因变量的影响[28]。分位数回归模型的建立与预测在SPSS 23.0软件中完成。

对于因变量Y的一个随机样本{y1,y2,…,yn}，通常τ分位数的样本分位数线性回归要求满足

其中，yi为因变量，β为自变量的回归系数，β(τ)为特定分位数τ的最终回归系统，x′i为自变量。

对于任意的0<τ<1，求解参数估计值的公式为

(1)

求得参数x′iβ(τ)即为唯一的第τ回归分位数。

1.5.2人工神经网络模型

人工神经网络模型从信息处理的角度对人脑的神经元进行抽象，建立某种简单模型，按不同的连接方式组合不同的网络。该模型拓扑结构由输入层、隐含层、输出层构成。经大量试算，本文输入层为不同深度下全变量及经过全子集筛选后的敏感光谱指数组合，输出层为土壤含盐量，隐含层中的网络层数为1。其中隐含层采用双曲正切激活函数，输出层采用恒等激活函数，可以在一定程度上防止梯度消失且结果可靠；训练拟合目标误差拟定为0.1%；培训类型采用适用于小数据集且精度高的批处理，并选用对应的标度共轭梯度估计突出权重进行网络训练；为了消除量纲、数量级等对数据分析造成的影响，将光谱指数标准化处理。本文采用SPSS 23.0软件中的多层感知器神经网络进行模型的建立与预测。

1.5.3支持向量机模型

支持向量机模型是在结构风险最小化原理的基础上，通过一个非线性变换将数据变换到一个高维特征空间，并在该特征空间建立线性模型来拟合回归函数，很大程度上克服了“离散值多”和“过学习”等问题[29]。径向基核(RBF)是应用最广泛的核函数。惩罚参数C和RBF核参数δ是影响SVM的主要参数，其中参数C直接影响模型的稳定性，避免模型在学习和推广过程中产生欠学习和过学习问题，决定了适应误差的最小化和平滑程度；参数δ反映了支持向量之间的相关程度，直接影响支持向量之间联系的松弛度，避免产生欠拟合和过拟合问题，决定模型预测的推广能力和泛化性[30]。本文核函数采用RBF核函数进行计算，利用R3.5.1软件进行支持向量机模型的建立与预测。

1.6 统计分析

本文通过决定系数R2、均方根误差RMSE、赤池信息准则AIC[31](Akaike information criterion)和施瓦茨信息准则SIC(Schwarz information criterion)来综合评价全子集筛选的效果。通过决定系数R2、均方根误差RMSE综合评价土壤含盐量拟合效果。其中R2越接近于1，表示土壤含盐量拟合效果越好，RMSE通过预测值和实测值的偏差度判断模型的准确性，值越小，表示预测值和实测值越接近，AIC和SIC是建立在熵的概念基础上，衡量统计模型拟合优良性的一种标准，有效避免过拟合现象的发生。AIC和SIC参数值越小表示该模型能够以最少自由变量最好地解释因变量[32]。其计算公式为

(2)

(3)

(4)

式中n——样本数k——自变量数目

L——似然函数

RSS——残差平方和

2 结果与分析

2.1 自变量与土壤含盐量的相关性分析

从遥感图像提取4个波段反射率和12个光谱指数，组成各个深度土壤含盐量模型构建的样本数据集，随机选择2/3的样本数据组成估算数据集(48个样本)，与土壤含盐量进行相关性分析，其结果如图2所示。参考相关系数检验临界值进行变量的显著性检验，当自由度为48，相关系数的绝对值大于0.365时，达到0.01显著性水平。从图2可知，B1、B2、B3、SI、SI1、SI3、S1、S2、RVI、NDVI、MSAVI、ARVI在各个深度下与土壤含盐量的相关系数的绝对值均大于0.365，达到0.01显著性水平。

图2 自变量与土壤含盐量的Pearson相关系数分析结果Fig.2 Results of Pearson correlation coefficient analysis between independent variables and soil salinity

2.2 基于全子集的最优自变量组合方式筛选确定

利用全子集筛选法列举不同数目自变量的随机组合方式；通过验证集R2确定不同深度下不同数目自变量的最佳组合方式见表3；根据R2、RMSE、AIC、SIC4种不同的验证指标确定不同深度下最优自变量组合方式。

从表3可以看出，在同一深度下随着自变量数目的增加，R2呈先增大后减小的趋势，AIC、SIC、RMSE呈先减小再增大的趋势。在0～20 cm、20～40 cm下，RMSE随着深度增加而减小，这是由于上层土壤容易受大气、人为等外界环境影响造成土壤含盐量变异性大，随着深度增加外部环境的影响基本消除精度逐渐升高。总体来说，R2呈现先增大后减小的趋势，在20～40 cm处达到最大值。结合表2、3可以看出，在相关性分析中B4和BI与不同深度土壤含盐量没有呈现极显著的相关关系，经过全子集筛选后B4和BI的组合方式在不同深度下与土壤含盐量均呈现良好的相关性，并成为各深度下自变量数目为2的最优自变量组合。在植被遥感中应用最为广泛的归一化植被指数NDVI与不同深度的土壤含盐量均呈现极显著的相关性，在20～40 cm、0～40 cm、40～60 cm、0～60 cm土壤深度下与B4、BI组合后成为该深度下自变量数目为3的最优自变量组合。在0～20 cm处，其最优组合方式为B4、BI、MSAVI，这是由于在0～20 cm处多是土壤和主要作物的侧根，土壤背景和植被的相互作用通过改进型土壤调整植被指数MSAVI减小土壤亮度的影响[33]。当自变量数目增加到4，主要增加了由B2和B3计算得出的SI1、SI3两种光谱指数，这是由于土壤一般都有很高的溶解性盐分，在520～770 nm平均反射率最高[34]，且属于确定不同盐渍化过程中积盐状态特征的6个光谱区间之一[35]。随着模型复杂程度的增加，AIC和SIC逐渐增大，但R2减小、RMSE增大，模型灵敏性降低。综合分析全子集筛选的各个评价指标，确定在0～20 cm、0～40 cm处选择B4、BI、SI1、SI3共4个自变量，在20～40 cm、40～60 cm、0～60 cm处选择B4、BI、NDVI共3个自变量作为各深度下最优自变量组合方式。

表3 全子集筛选最佳组合方式结果统计Tab.3 Best combination result after total subset selection

注：** 表示0.01显著性水平，*表示0.05显著性水平。

2.3 人工神经网络模型建立与分析

分别以不同深度下筛选前后的反射率及光谱指数为自变量、土壤含盐量为因变量，运用人工神经网络模型进行筛选前后不同深度下土壤含盐量估算，其建模及验证结果如表4所示。

表4 基于不同深度土壤含盐量的人工神经网络模型Tab.4 ANN model of soil salinity at different depths

2.4 支持向量机模型建立与分析

以不同深度下筛选前后的反射率及光谱指数为自变量，土壤含盐量为因变量，运用支持向量机模型进行筛选前后不同深度下土壤含盐量估算。其建模和验证结果如表5所示。

表5 基于不同深度土壤含盐量的支持向量机模型Tab.5 SVM model of soil salinity at different depths

2.5 分位数回归模型建立与分析

对于不同深度选定的敏感植被指数建立分位数回归模型，由于采样点土壤含盐量变异性较大，选取τ= 0.5分位点可以较好地解决最小二乘法中某些“离群值”影响回归显著性的问题。同时，由于0.5分位点处于因变量的中间位置，在对所有的数据进行拟合时较为适宜[36]。故本文以不同深度下筛选前后的反射率及光谱指数为自变量，以土壤含盐量为因变量，运用分位数回归模型中的0.5分位点进行筛选前后不同深度下土壤含盐量估算，其建模及验证结果如表6所示。

表6 基于不同深度土壤含盐量的分位数回归模型Tab.6 QR model of soil salinity at different depths

2.6 模型综合评价

综上所述，人工神经网络模型的反演效果最差，分位数回归模型和支持向量机模型在本次土壤含盐量估算中就模型效果而言，二者基本相同。但是分位数回归模型仅采用3个自变量，支持向量机模型采用16个自变量。因此，本文选择拟合效果好、验证精度高、模型稳定性强、简洁高效的20～40 cm深度的全子集-分位数回归模型作为本次最佳土壤盐渍化估算模型。

图3 土壤含盐量空间分布图Fig.3 Soil salinity distribution map

基于全子集-分位数回归方法建立模型，进行解放闸灌域土壤含盐量估测，结果如图3所示。图中白色部分为城镇所在区域不参与本次含盐量反演。在该时期，灌域内盐渍化程度较轻，有利于作物的生长发育。灌域内轻度盐渍化土所占比例最大为54%，非盐土所占比重次之，为26%。这是由于该时期处于作物生长关键期，灌水次数多且灌水量大，土壤盐分淋洗到下层土壤，土壤处于夏季脱盐状态。重盐渍化土所占比重较少，为13%，对于灌水次数较少的地方，土壤仍产生积盐现象。盐土占总面积的比例最少，为7%，主要集中在灌域西部(多为盐荒地)，其余盐土零散分布于灌域内，主要包括部分盐荒地和弃耕地。

3 讨论

土壤盐渍化问题在中国北方干旱灌区较为突出，已成为制约灌区农业可持续发展的重要因素。黄河流域上游的河套灌区是受盐渍化影响的典型区，土壤盐渍化面积约占总面积的69%[37]。利用卫星遥感对土壤含盐量进行大面积监测，符合未来精准农业发展的要求。在本次植被覆盖条件下，作物长势可间接识别土壤根域盐渍化水平，研究结果显示20～40 cm深度下的土壤含盐量与高分一号卫星数据相关性最高，土壤盐分含量对作物主根系的影响明显。史海滨等[38]对比各生育期内不同含盐处理的各层主根重的变化情况发现，向日葵的苗期后期(7月27日)主根深度在30 cm之内，土壤盐分含量对作物主根系层的生长影响明显。与本文研究成果基本一致。由此可见，应用卫星遥感大范围监测植被覆盖条件下不同土层的盐渍化水平具有一定的可行性。

采用相关性分析只能获得单一自变量与含盐量的相关关系，本研究采用全子集法列举全部组合方式，得到不同深度下对于土壤含盐量最优自变量数目的敏感变量组合，筛选后的结果与未筛选的结果在拟合优度和模型精度方面略有提升。表明全子集筛选过程简单、高效、灵敏度强，适用于本次盐分含量的估测。OKCU等[39]采用全子集法生成一个新的适用于河流和水槽数据的泥沙运移公式，精准预测泥沙运移的非线性关系。但是由于自变量数目的增加，全子集筛选的运行效率减弱，HOFMANN等[40]改进和扩展全子集回归的计算方法，采用基于回归树的RadiusBBA和HBBA方法解决大规模模型选择问题的穷举和启发式策略。全子集法和相关性分析法的筛选结果可作为自变量寻优的一种参考，而筛选结果与因变量是否有内在联系，需要做更严谨的物理机制和数学推理工作。

本文对于不同深度选定的敏感指数建立的ANN、SVM、QR 3种模型进行土壤含盐量反演，发现QR模型的反演精度最高。这是由于QR采用非对称权重的方法估计参数，适用于离群值多、变量间关系较弱和需要了解因变量分布的情况，且无需对模型进行任何分布的假定，对于强变异性(变异系数约100%)的实测土壤盐分具有抗耐性。AMAKOR[19]进行钙质盐碱土电导率模拟时，实测电导率变异系数在100%左右，采用QR模型决定系数R2达到了0.88。MUELLER等[41]、GERBER等[42]和CARSLAW等[43]采用QR模型，均很好地解决了实测数据重尾分布、变异性强的问题。τ=0.5分位数处于土壤含盐量的中间位置，采用该分位数进行各个深度的土壤盐渍化反演可以保证轻、重度盐渍化土权重较大，非盐土和盐土权重较小，从而减小离群值对于模型精度的影响。前人也得到类似的结果，如王蕾等[36]发现在0.5分位数下进行冬小麦估测结果最为可靠。

研究植被覆盖条件下不同深度的解放闸灌域土壤盐渍化模型，对改善灌区土壤盐渍化情况具有现实意义。然而最优土壤含盐量反演模型会因作物的种类、生育期、天气状况、灌溉情况、甚至所使用的遥感平台而异。本文所得的最佳含盐量反演模型也仅限于本次测量结果，研究区内其他生育阶段和其他地区还可尝试更多图像特征(如颜色特征[44]、地形特征[45])，寻找更优的含盐量反演指标，从而建立精度更高且实用性更强的土壤含盐量反演模型。

4 结论

(1)全子集筛选具有全局性，经过全子集筛选后敏感变量组合方式不尽相同。B4、BI、SI1、SI3是0～20 cm、0～40 cm处土壤含盐量的敏感变量组合，B4、BI、NDVI为20～40 cm、40～60 cm、0～60 cm处土壤含盐量的敏感变量组合。筛选后的敏感变量组合方式仅用3～4个自变量，与筛选之前的效果相差不多，因而全子集筛选具有灵敏度强、简单、高效的优势。

(2)在不同深度下的土壤含盐量反演模型中，分位数回归、人工神经网络、支持向量机模型均取得良好的建模验证精度。分位数回归模型由于其抗异性较强，其精度与支持向量机模型相差不多，优于人工神经网络模型。土壤深度对含盐量反演模型精度具有一定的影响。在使用相同的建模方法时，20～40 cm深度下反演模型效果优于其他深度。

(3)在20～40 cm深度下建立的模型效果最为理想，且分位数回归抗异性、稳健性强，在保证模型估测精度的前提下，比支持向量机模型和人工神经网络模型操作简单、结果可靠。因此，全子集-分位数回归模型是本次估算植被覆盖条件下土壤含盐量的最优模型。