基于机器学习的绿洲土壤盐渍化尺度效应研究

陈香月 丁建丽 葛翔宇 王 飞 王敬哲

(1.新疆大学智慧城市与环境建模自治区普通高校重点实验室， 乌鲁木齐 830046；2.新疆大学绿洲生态教育部重点实验室， 乌鲁木齐 830046;3.深圳大学自然资源部大湾区地理环境监测重点实验室， 深圳 518060)

0 引言

土壤盐渍化指水溶性盐分向土壤表面迁移和聚集的过程，是自然因素和不合理的人类活动共同作用影响下的产物[1-3]。在人口快速增长背景下，现有土地资源难以继续承载更多的人口生产生活所需，合理开发利用后备耕地资源，对缓解土地资源稀缺绿洲地区的人地矛盾尤为关键，也成为近年来新的研究热点。如何准确、快速、廉价获取盐渍化信息是科学合理开发盐渍土地资源的重要依据。近年来，遥感技术因其具有回访周期快、覆盖面积广、成本低、获取方便，且非侵入土壤等优势，已成为土壤盐渍化信息数字制图研究中广泛应用的有效工具[4]。

地表的异质性意味着就单一尺度研究地学问题具有一定局限性，通常基于特定尺度的研究结果较难直接推广至其他尺度，尺度上的细微变化都将引出新的科学问题[5-6]。尺度问题研究的根本目的是挖掘时空尺度特征信息，从而找到与研究对象的观测尺度相匹配的本征时空尺度。土壤盐渍化往往具有极强的时空变异性，受各类环境要素影响极大，对土壤盐渍化现象发生发展的深刻理解都离不开尺度问题这一环节[7-8]。由于遥感影像具有多分辨率(多尺度)特性，不同尺度在遥感数据中所呈现的地物信息丰富度存在显著差异，通常某些规律往往在特定的尺度下才被发现，尺度的不适配将会降低土壤盐度模拟效果，故探寻最优研究尺度对于干旱区绿洲土壤盐渍化研究具有重要意义。

在土壤研究中，基于遥感数据计算的环境变量常被用于预测土壤属性的空间分布。盐分含量及其组成成分是影响盐渍化土壤光谱特征的主要因素，已有研究表明盐分含量增加会引起反射率增加[9-10]。植被作为土壤盐渍化程度的间接指示器，其光谱特性是各种生态环境因子的综合表达，而为了有效理解这种表达，需要在合适的尺度下开展研究[11-12]。目前，已有部分学者对土壤属性与遥感数据尺度(栅格单元大小)之间的关系进行了初步分析[13-14]，并在此基础上对其衍生的生物物理过程开展相关研究[15-16]，但基于植被光谱特征信息间接推理土壤属性研究多停留在单一尺度上，对不同尺度下环境变量对土壤属性的解析尚不够深入。而忽略尺度效应可能会导致难以充分理解“土壤-环境”间的复杂关系，更会限制土壤盐渍化遥感模型的性能[17]。

鉴于此，本研究针对土壤盐渍化这一典型土壤退化过程开展尺度效应研究，以新疆北部典型绿洲——奇台绿洲为研究区，利用遥感衍生数据与土壤表层实测电导率数据，定量描述土壤盐度与环境变量关系之间的尺度影响，结合机器学习算法通过多环境变量对土壤盐度分尺度进行模拟，并获取不同尺度下各环境变量对土壤盐度的解析能力，以期找出最优响应尺度，充分理解不同尺度下环境变量与土壤盐度之间的复杂关系，对绿洲土壤盐渍化状况进行有效评估，为干旱区绿洲农业可持续发展提供科学参考与数据支持。

1 数据与方法

1.1 研究区概况

奇台绿洲位于新疆维吾尔自治区昌吉回族自治州奇台县境内，地处新疆维吾尔自治区东北部，博格达山北麓，天山山脉东侧。研究区(43.95°～44.20°N，89.20°～90.00°E，图1)位于奇台县平原地区，属温带大陆性干旱半干旱气候，夏季高温少雨，日温差、年温差较大，多年平均降水量为211.9 mm，多年平均蒸发量为1 735.7 mm，蒸降比高达8∶1，强烈的蒸散发作用及水资源的不合理利用导致奇台县内广泛存在土壤盐碱化问题。全区约有31%农业土地存在盐碱化现象，严重危害到当地的生态环境与农业生产生活[18]。

1.2 采样数据

本研究利用的样本采样时间为2009年6月，样本量为58个，样本主要集中在绿洲平原区，主要的土地利用与覆盖类型包括：农田、荒漠、灌丛、稀疏草地、草地等。以梅花五点混合采样方法采取表层土壤(0～20 cm)约500 g，同时记录样品相应的地理坐标，标记后用自封袋封装带回实验室。采集得到的样本在室内剔除杂质晾干后进行研磨与过筛处理(2 mm)，采用土水质量比1∶5制备土壤悬浊液，利用土壤溶液电导率仪对制备的浸提取液进行电导率(dS/m)测定[17]。

1.3 遥感数据

Landsat TM是搭载于Landsat 4/5卫星上的光学成像传感器，提供了大量中分辨率和长时间序列的遥感数据。其光谱范围从可见光波长(0.45 μm)到中红外波长(12.5 μm)间共7个波段，重访周期为16 d，空间分辨率为30 m。从美国USGS官网下载2009年6月过境的Landsat TM影像，行列号为141/29，以完全覆盖所有采样点且云量覆盖低于10%为准。通过ENVI软件对数据进行辐射定标、大气校正，最后将辐射率转为研究需要的地表反射率以进行后续分析。

1.4 尺度转换

地理空间数据多具有尺度依赖性，尺度效应被认为是对地观测中不可忽略的问题，严重制约了遥感、地信等学科的发展，故在研究中考虑尺度问题非常必要[19-20]。魏阳等[21]研究表明，就干旱区绿洲土壤而言，在地面分辨率1 km左右开展土壤盐度尺度效应研究是有效的。基于此，本研究采用栅格重采样和邻域滤波两种升尺度方法，将最大分辨率控制在1 km左右，以此为阈值对影像进行升尺度转换。栅格重采样方法选用最邻近法，基于像元对影像进行重采样计算，邻域滤波方法主要是基于邻域窗口的选择，在特定空间分辨率下按照n×n移动窗口分别计算3×3、5×5、…、31×31多个邻域滤波，其中窗口尺度为奇数，数值间隔为2，见图2。具体而言，利用ENVI软件中的重采样模块(Resize data)将空间分辨率为30 m的原始遥感数据进行栅格重采样处理，研究以30 m为步长将原始遥感数据依次重采样至990 m。在邻域滤波计算中，选择了7种原始分辨率，分别为30、60、90、120、150、180、210 m，邻域滤波窗口尺度选择了3×3、5×5、… 31×31共计15个级别，共对应39个尺度，利用ENVI中的邻域滤波模块(Occurrence measures)实现。

1.5 环境变量

基于尺度变换后的遥感影像，计算应用于土壤盐度评价的植被指数和归一化红外指数[17]以及图像衍生变换参量(表1)。植被与土壤关系密不可分，作为地表变化的直接观测因子，其长势能够间接反映土壤盐渍化程度[22-23]。且已有研究表明，相较于单纯利用土壤信息进行盐度预测，植被光谱信息的引入可以有效改善预测效果[24]。其中，归一化植被指数(NDVI)、扩展归一化植被指数(ENDVI)以及扩展增强性植被指数(EEVI)已经被证实可以间接用以判断土壤盐度[25-26]。归一化红外指数(NDII)，亦称土壤湿度/水分含量指数，通过融合Landsat TM的短波红外波段信息，使得NDII对冠层湿度与土壤盐度具有一定指示作用。图像变换衍生因子某种程度上可以对原始波段信息进行高度凝练，对海量数据进行有效降维，避免维数灾难。其中主成分分析(Principal components analysis, PCA)通过生成互不相关的输出波段，研究表明主成分分析的前3个波段(PCA1、PCA2、PCA3)占数据整体信息高达90%以上，达到隔离噪声和减少数据集的维数目的；缨帽变换(Tasseled cap transformation，TC)是一种特殊的主成分分析，该变换将地表信息有效地转换成3个主题因子，分别为亮度(TCB)、植被绿度(TCG)和土壤湿度(TCW)，可近似作为土壤属性的相关指示因子。通过上述两种尺度转换方法共产生720个环境变量用于耦合多尺度和机器学习算法的土壤盐度估算。

表1 选取的环境变量Tab.1 Environmental variables derived from pretreated remote sensed imagery

1.6 梯度提升决策树

梯度提升决策树(Gradient boosting decision tree，GBDT)[27-28]是一种由多棵决策树组成的集成学习方法，它采用Boosting思想将多个弱分类器组成一个强分类器。该算法的基本原理是通过迭代计算，不断减小模型的残差，在梯度方向上建立一个新的模型来代替旧的模型[29]，并在构建新模型的过程可以迭代地生成基本分类器的组合，从而将损失函数最小化。作为典型的集成算法，GBDT减少模型在建模过程中梯度方向上的残差，相较于传统的机器学习模型具有更好的预测能力[30]，不仅能够灵活地处理各类离散数据，同时无需过分调节超参数，更为重要的是，GBDT对输入变量的数量及其自相关不敏感，还可以使用一些较为稳健的损失函数，对异常值具有较好的鲁棒性[31-32]。针对遥感影像这一高维数据源，GBDT算法可最大程度实现信息的准确提取。因此，本研究基于GBDT，在不同尺度下利用环境变量进行土壤盐度模拟，并获取不同环境变量对土壤盐度的贡献度。随机抽取建模集，并固定验证集，以此优化模型参数，经多次训练后，GBDT模型具体参数设置为：迭代回归树数量为300，学习速率为0.01，叶子节点的最小个数为7，交叉验证折数为3，算法依托R语言平台实现，流程如图3所示。

2 结果与分析

2.1 GBDT模型构建及验证

2.2 栅格重采样变换后土壤盐度估算

基于GBDT算法，对不同栅格重采样尺度下土壤盐度进行预测，并计算相应尺度下不同环境变量对土壤盐度的贡献权重及其总体解析能力。为可视化表达不同尺度下环境变量对土壤盐度的解析效果，本文引入泰勒图进行相关表述(图4，图中黑色虚线表示调整R2，蓝色弧线表示RMSE(dS/m)，当调整R2越高、RMSE越小时模型效果越好)。图4a为不同栅格重采样尺度环境变量对土壤盐度的解析能力，尺度范围为30～990 m，在不同尺度条件下，不同环境变量对土壤盐度的响应也不尽相同。

表2 土壤盐度预测效果Tab.2 Statistical results of estimated soil salinity based on various scale conversion modes

由图4a可知，上述栅格重采样尺度的调整R2在0.56～0.75之间，RMSE在1.10～1.47 dS/m之间，其中在栅格重采样尺度60 m下，对于解析土壤表层盐度-环境变量效果最好，解析能力达到了75.31%(R2>0.75)，为单一栅格重采样尺度下研究土壤表层盐度-环境变量关系的最佳尺度。结合不同尺度下各环境变量对土壤盐度贡献权重(图5)发现，对应栅格重采样尺度60 m下，对土壤表层盐度贡献率从高到低依次为NDII、EEVI、PC2、ENDVI，其贡献率总和达到了77.32%，解释了77.32%的土壤盐渍化变异，这也佐证了本研究使用植被指数对表层土壤盐度进行表征是可行的。随着尺度的上升，GBDT算法对于环境变量和土壤盐度的关系描述效果越来越不理想。其中，起主导作用的环境变量从原有与植被相关的TCG、植被指数等变成了图像衍生变换参量，这可能归因于随着分辨率的下降，原始数据所包含的信息丰富度损失，植被特征被弱化，而PC、TC等参数以其独有的信息优势过渡为主导因素。

2.3 邻域滤波变化后土壤盐度估算

对30～210 m栅格重采样分辨率下对应计算邻域窗口3×3、5×5、…、31×31的邻域滤波进行土壤盐度的预测，并获取各环境变量对土壤盐度影响权重，结果见图4、6。由图4b可知，在30 m分辨率多种邻域滤波尺度下，土壤盐度-环境变量之间的调整R2最大值为0.79，RMSE为1.02 dS/m，可解释78.55%的土壤盐度变化，对应邻域窗口11×11，地面分辨率为330 m、响应的环境变量依次为TCG、ENDVI、PC3、NDVI、NDII，能够解释近78.57%的盐度变异性。图4c、4d为60～210 m栅格重采样分辨率下多种邻域滤波尺度对土壤盐度与环境变量关系的表征，调整R2为0.66～0.78之间。当基础栅格重采样尺度为60 m时，R2在0.68～0.78之间，RMSE为1.05～1.25 dS/m，其中当地面分辨率为300 m、邻域窗口为5×5时效果最好；90 m时，R2在0.69～0.74之间，RMSE为1.05～1.24 dS/m，其中当地面分辨率为270 m、邻域窗口为3×3时效果最好；120～210 m时，R2在0.66～0.76之间，其中当窗口为5×5时效果最好；对土壤盐度的解析力分别达到76.01%、68.11%、69.41%、67.97%。可以发现，当初始分辨率为30、60、90 m和120～210 m时，对应最佳邻域窗口分别为11×11、5×5、3×3、5×5，其中当初始空间分辨率为30～90 m时，最优土壤盐度表征尺度为300 m左右。单一栅格重采样尺度下，伴随着邻域窗口的变大，土壤盐度对环境变量的解析程度也在下降；单一邻域尺度下，随着栅格重采样尺度的增大其对环境变量的响应也在弱化。即随着尺度上升，调整R2总体呈下降趋势，RMSE呈上升趋势，其中与土壤盐度相关的环境变量集中为与植被相关的TCG和植被指数。

2.4 多尺度协同作用下土壤盐度估算

通过GBDT算法结合多尺度环境变量对土壤盐度进行模拟，并获取多尺度下环境变量指征土壤盐度的最优组合。图7为多尺度下基于原始分辨率数据和邻域滤波尺度，考量环境变量对土壤盐度-环境变量的协同影响。在参与分析的30～990 m所有栅格重采样尺度时，其调整R2可达到0.91，RMSE为0.68 dS/m，对土壤盐度的解析度达到90.63%。对邻域滤波，分别以30、60、90、120、150、180、210 m为基础分辨率，讨论在上述7个基础栅格重采样分辨率下随着邻域窗口的变化环境变量与土壤盐度的响应情况，结果如图7所示，在空间分辨率为30、60、90、120、150、180、210 m时，其对土壤盐度预测的调整R2分别为0.86、0.85、0.82、0.80、0.77、0.73、0.67，RMSE分别为0.81、0.85、0.94、0.99、1.06、1.14、1.27 dS/m，随着基础空间分辨率的下降，调整R2降低，环境变量对土壤盐度的解析力度逐渐降低。当考虑30～210 m所有邻域滤波的情况下，对土壤盐度的调整R2可达到0.88，RMSE为0.78 dS/m，由此可见，相较于邻域滤波变换，栅格重采样可以更为显著提高土壤盐渍化建模效果。当考虑全尺度情况下，调整R2达到0.91，RMSE低至0.68 dS/m，对土壤盐度的解析度达到90.66%，与栅格重采样多尺度下性能接近，故当考虑全尺度时，增加邻域滤波尺度对环境变量与土壤盐度关系的研究改善有限。综上所述，多尺度协同极大地提高了环境变量对土壤盐度的解释力，其中栅格重采样尺度对环境变量与土壤盐度关系的研究起到主导作用。

3 讨论

在植被覆盖度高的区域，土壤盐渍化会造成植被红光波段反射率增加、近红外波段反射率下降，故植被指数的变化有助于推测盐渍化的程度和发展趋势[3]。基于此，在绿洲地区，已有许多学者开展了相关研究，结果表明，因干旱区地理环境存在显著的差异性，环境变量在某个区域表现出的高解析力，其泛化能力尚不可知；不同盐度等级下的土壤盐度光谱信息存在明显差异，遥感指数在盐渍化表征上具有复杂的不确定性；受多种因素影响，环境变量的数量和类别上还有待进一步优化[21]。在本研究中，随着栅格分辨率的上升，环境变量对土壤盐度的解析能力从与植被相关的TCG、植被指数等变成了图像衍生变换参量，如PC1、PC2。对邻域滤波而言，TCG、PC3、NDII、TCW和ENDVI可解析69.83%的土壤盐度变异，对土壤盐度的预测具有一定的应用优势。不同栅格重采样和邻域滤波窗口尺度下，环境变量与土壤盐度的关系存在较为明显的差异，这说明了环境变量与土壤盐度之间存在强尺度依赖性。MULLER等[33]在多种分辨率下对比了植被指数与土壤盐度之间的耦合关系，发现随着分辨率的降低其响应尺度也在降低，这与本研究得到的结论基本一致，如在30 m分辨率时，最佳邻域窗口为11×11，随着分辨率变为60、90 m，最佳邻域窗口分别为5×5、3×3。然而，值得注意的是当分辨率在120～210 m时，最佳邻域窗口稳定在5×5，这可能与研究区独特的地理特征及土壤属性的空间变异性有关。

尺度问题是生态学研究中最具挑战的问题之一，也是遥感生态学需要面对的主要问题，反演的准确性和遥感数据的空间精度有很大的关系[34-35]。本研究以国内外研究者广泛选用的Landsat 影像为数据源，针对土壤盐渍化现象，分析了土壤盐度-环境变量关系的尺度效应。但该数据的空间分辨率(观测尺度)为30 m，仅依据单次采样数据结合遥感对土壤盐度进行模拟，仍存在一定的不确定性。随着遥感技术的发展及传感器分辨率的不断提升，更优质的遥感数据的出现(更高空间分辨率、更宽光谱范围等，如欧洲空间局的Sentinel-2多光谱遥感数据)，为进一步研究尺度对土壤盐度-环境变量之间耦合关系提供了新的选择。同时，本研究使用固定观测尺度，栅格重采样的步长为30 m(像元逐个加入)，但这可能造成某些尺度被遗漏。而高空间分辨率数据的引入则可以缩短研究间隔，使得更为细致地研究环境变量-土壤盐度的尺度效应成为可能。需要注意的是，由于地表空间异质性的存在会影响土壤属性或者过程模拟的精确程度，故合适的栅格分辨率和邻域滤波尺度的选择最终依赖于研究区的地域特点，即被模拟属性的空间变异性和分布规律。

4 结论

(1)单一尺度下，邻域滤波尺度对绿洲土壤盐度的模拟总体优于重采样模式。在单一栅格重采样尺度下，60 m栅格重采样尺度对土壤表层盐度-环境变量解析效果最好，解析能力达到75.31%；单一邻域滤波尺度下，初始分辨率为30 m、邻域窗口为11×11是土壤表层盐度-环境变量最佳表征尺度，解析能力为78.55%，随着邻域窗口尺度的外推，其对土壤盐度的响应效果逐渐减弱。

(2)栅格重采样模式相对邻域滤波而言，其调整R2波动范围更为宽泛，分布范围为0.56～0.75，而后者变化平缓。一定程度上说明，栅格重采样方法对土壤盐度-环境变量关系的变化较为敏感，邻域滤波对土壤盐度-环境变量关系的敏感程度相对较弱。

(3)当考虑融合多种尺度时，环境变量对土壤盐度的模拟精度相比仅利用单一尺度而言，得到极大程度提升，具体表现在R2从0.78提高到0.91，RMSE从1.02 dS/m降低到0.68 dS/m，对土壤盐度的解析度最高可达到90.66%，这也体现了不同环境变量协同作用可最大程度发挥其最大响应效力，同时变量之间的信息进行了有效互补。