土壤盐渍化是严重威胁农业生产的全球性问题，我国受盐渍化影响的土地面积(3.6×10⁷ hm²)占全国可利用土地资源总量的4.88%^[1]。绿洲区域作为新疆重要的粮食储备区，受盐渍化的影响，制约了农业生产的可持续发展。随着现代数学发展以及遥感数据类型日益广泛和易于获取，一种新兴、高效的技术方法——数字化制图(digital soil mapping, DSM)在土壤空间分布制图研究中得以应用，并得到了飞速发展^[2]。目前主流的数字土壤制图方法主要有两种，第一种是基于要素相关性，包括统计学方法、机器学习方法和基于专家知识的方法，第二种是基于空间自相关，包括全局自相关和局域自相关分析，其中克里格插值法是应用最为广泛的一种方法。

构建反演模型定量化的区域土壤盐分动态监测，有助于盐渍化区域制定改良措施和优化政策，已有大量研究尝试不同方法对土壤盐渍化进行反演。Vermeulen和van Niekerk^[3]运用k邻近、支持向量机和随机森林等方法，结合DEM预测土壤盐害影响范围，证明高程数据在监测灌区盐分积累方面具有很大潜力。Zhang等^[4]运用MODIS时间序列植被数据，通过内插法模拟年际植被指数变化，结合季度植被指数对黄河三角洲的土壤盐分进行数字制图，结果表明结合长时间序列植被信息较大地提升了预测精度。Wang等^[5]结合环境变量预测新疆渭库绿洲土壤盐度信息比较了13种算法及其性能，证明随机森林方法具有较好效果。

以往大多数研究证明植被指数在盐分预测方面表现良好^[6-7]，但多数研究仅针对特定时期进行，MODIS长时间序列数据能满足长期观测的条件，而植被物候作为连续的环境参量在之前的盐分制图工作中被忽略或低估^[8]。此外对于土壤盐渍化监测，大都集中在田间尺度和流域尺度^[9-11]。随着土壤环境问题涉及尺度愈发扩大，当前土地规划、生态文明及粮食安全需要决策者从宏观角度进行把握^[12]。因此，本研究以新疆绿洲区域作为研究靶区，借助随机森林(random forest, RF)模型，结合研究区域的典型采样点，耦合物候参数和地表参数反演土壤盐分信息，通过绘制新疆绿洲地区土壤盐渍化空间分布图，深入理解绿洲区域土壤盐渍化的形成及分布，以期为物候信息在大尺度的盐渍化特征监测应用中提供参考。

1 材料与方法 1.1 研究区概况

新疆维吾尔自治区(73°40′~ 96°18′E，34°25′~ 48°10′N)，深居欧亚大陆腹地，远离海洋，四周为高山环抱，海洋气流通达性不强，形成明显的温带大陆性气候。年均降水量50 ~ 200 mm，蒸发量高达1 500 ~ 3 000 mm，蒸发强烈造成地表土壤无机盐含量升高，形成大面积盐碱土。除气候因素外，人为因素加速土壤次生盐渍化，导致约32% 的耕地出现次生盐渍化。依据第二次新疆地区土壤普查结果，新疆土壤可划分为7个土纲及32个土类，其中风沙土、棕漠土、棕钙土、寒冻土及石质土占比较大。新疆主要植被类型为旱生灌木、小半灌木和乔木等，其稀疏根系发达且较耐干旱，农田植被主要由玉米、小麦和棉花等构成，大多分布在绿洲内部水资源丰富的区域。新疆绿洲面积约占总面积的5% ~ 8%，集中了新疆90% 的耕地、人口和国内生产总值。由于水资源的不合理利用及不科学的耕作方式，致使绿洲生态环境恶化，盐渍化现象依然普遍，严重影响了绿洲的可持续发展。

1.2 土壤盐分数据采集

根据土壤盐渍化状况的不同，本研究在新疆地区选择两个不同的研究靶区采集土样，分别是位于天山以南的渭干河–库车河绿洲(简称渭库绿洲)和天山以北的艾比湖流域。基于相似性的土壤制图方法无需要求采样点的数量或分布，它根据每个地点与采样点的环境相似性来预测土壤盐度值。在这两个地区选择代表性采样区，保证其包括不同地理环境和典型气候区，融合两个地区的实测数据，使之具备与新疆其他区域近似的地理和气候环境。从研究区中心到外围按均匀分层采样，采用“S”型布点避开城镇用地，选择大量的典型覆被类型，例如绿洲荒漠交错带、盐渍地、沙地、荒漠进行采样，尽量保证样点均匀地分布在样区内部，如图 1所示。采用网格法布设样方，在地势较为平坦之处选取土壤样方(30 m×30 m)，保持土壤性质相对一致，环境因素相似。采用五点法进行表层土壤(0 ~ 10 cm)样品的采集，将测试的数据求平均值作为本样点的实际观测值。采集的土样自然风干清除石块和杂草后，研磨过2 mm孔筛，制备水土质量比为5∶1的浸提液，采用Cond7310土壤测试仪测定土壤上清液电导率，采用残渣烘干法测定土壤全盐量。采样时间为2019年7月，剔除无效样点，选取209个采样点(渭库绿洲113个，艾比湖96个)进行分析。

1.3 地表覆盖数据获取

本研究将新疆绿洲区域作为研究靶区，采用王智等^[13]的分区方法，首先按照新疆行政区边界将全疆划分为北疆、南疆、东疆和伊犁4个大区。其次，结合DEM数据在四大分区基础上对北疆和伊犁地区以海拔小于1 300 m为基础划分出平原区域，而南疆和东疆地区则以海拔小于2 000 m为阈值同样划分出平原地区。最后，运用2009年MODIS 500分辨率的土地利用产品数据MCD12Q1，提取出平原地区的耕地和城镇居民用地，以深入了解绿洲边缘与部分交错带区域的土壤盐分分布，在消除破碎斑块后，对提取出的两类土地利用数据向外做3 km缓冲区划为绿洲区域。

1.4 地表参数获取

土壤盐渍化在不同尺度下受土壤、气候、地形和生物等因素的综合影响，对于不同研究区，影响土壤形成发展最主要的环境因子可能有所不同^{[12, 14]}。本文选择红波段(R)、蓝波段(B)和近红外波段(NIR)作为反射率参数，并以蓝、红、近红外波段组合出23种地表遥感参数，见表 1。

1.5 物候参数计算

在高植被生物量水平下NDVI容易受到土壤背景变化、噪声、气溶胶消光系数和饱和度的影响^[23]，因此选取EVI进行植被物候参量的提取^[15]。本文所选用的MODIS EVI植被指数产品MOD13A2数据虽然已通过16 d的MVC预处理，但仍不可避免大气和噪声等的影响，直接提取对物候参量影响较大。通过TIMESAT 3.2的Data plotting模块对比Asymmetric Gaussian、Double Logistic和Savitzky-Golay三种平滑拟合方法，发现Savitzky-Golay滤波法对时间序列包络数据保留更加完整且对年生长季原始植被指数拟合更精准，因此，选用Savitzky-Golay滤波法重建物候曲线。

新疆地区横跨不同气候带且植被类型多样，不同土地覆盖类别其物候曲线也会有所差异。采用阈值法通过对植被指数设定阈值条件来确定植被物候的始末期指标EOS和SOS，该方法可以忽略地表的覆盖类型和人为活动的影响，根据研究区植被不同的生理特点，形象地描述植被的生长发育过程^[24-25]。根据相关研究^[26]，本文将提取生长季始末日期的系数设为30%，其计算公式如下：

$ \text{EVI}（\text{SOS}）=（{\text{EVI}}_{\text{max}}-{\text{EVI}}_{\text{min1}}）\times 30\% $

(1)

$ \text{EVI}（\text{EOS}）=（{\text{EVI}}_{\text{max}}-{\text{EVI}}_{\text{min2}}）\times 30\% $

(2)

式中：EVI(SOS)和EVI(EOS)分别表示最后确定的阈值，EVI_max为植被指数最大值，EVI_min1为植被指数曲线上升阶段最小值，EVI_min2为曲线下降阶段最小值。

根据平滑和参数设置，用TIMESAT软件计算出SOS、EOS和LEN等10个物候参量，其代表的生态学意义见表 2。

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

随机森林(RF)方法是一种利用多棵决策树对样本进行训练并预测的机器学习算法，该算法的优势在于具备非线性挖掘能力，数据的分布不需要符合任何假设^[27]。RF模型需要用户定义3个参数：终端节点树(ntree)、作为每棵树的预测特征的特征个数(mtry)，以及每个终端节点的最小值^[28]。ntree的默认值为500，但一般认为它不足以产生可靠的结果^[29]，本研究中通过遍历确定ntree的最优值为1 000。在使用mtry确定单一树与模型中其他树之间的相关性时，随着mtry值的增加每棵树和树之间的相关强度逐渐增加^[30]。本研究以预测误差最小化为目标函数，通过遍历比较确定mtry的最优值为9，终端节点大小为5。同时采用较常用的交叉验证法将样点数据随机分为5个子集，其中4/5作为训练集，1/5用于模型验证。选择决定系数(R²)、均方根误差(RMSE)和平均绝对误差(MAE)来评估预测模型性能，经过反复迭代获得最终的精度。R²范围在0 ~ 1之间，越接近1，模型的稳定性越高，RMSE和MAE越小，模型的预测能力越强、精度也越高。

$ {R^2} = 1 - \frac{{\sum\limits_{i = 1}^n {{{\left( {{V_{pi}} - {V_{oi}}} \right)}^2}} }}{{\sum\limits_{i = 1}^n {{{\left( {{{\bar V}_{oi}} - {V_{oi}}} \right)}^2}} }} $

(3)

$ {\text{RMSE}} = \sqrt {\frac{1}{n}{{\sum\limits_{i = 1}^n {\left( {{V_{oi}} - {V_{pi}}} \right)} }^2}} $

(4)

$ {\text{MAE}} = \frac{1}{n}\sum\limits_{i = 1}^n {\left| {{V_{oi}} - {V_{pi}}} \right|} $

(5)

式中：V_oi为样本点i处的实测值；V_pi为点i处的预测值；n为总样点个数。

2 结果 2.1 物候参数提取结果

对绿洲区域的10种物候参数进行提取，根据其对于盐渍化预测的相对重要性，仅对生长季大生物累积量(LSI)和生长季小生物累积量(SSI)2个物候参数进行展示。从图 2可以看出，两类物候参量分布情况高度一致，均呈现由绿洲外围向绿洲内部递增状态，从侧面可以反映出土壤盐渍化信息，即在绿洲外围盐渍化情况一般较重，植被长势不佳，生物量较低，值域在1.32 ~ 3.14。对比发现LSI和SSI在北疆区域分布几乎一致，整个伊犁地区生物积累量最大集中在4.58 ~ 11.04，说明其植被长势良好，而在南疆绿洲区域存在LSI明显高于SSI的情况，说明南疆绿洲区域大生物累积量高而小生物累积量较低，这可能与该地区本身的气候干燥有关。纵观LSI可以发现，其在北疆绿洲区域和南疆绿洲区域的分布情况基本一致，大部分区域集中于4.61 ~ 9.07。一般正常生长的植被生物量较高，而受盐渍化或病虫害影响的植被长势受限，积累生物量会呈现较低水平，因此，LSI和SSI参考在监测土壤盐渍化水平中具有极高的应用潜力。

2.2 模型精度验证及变量遴选

由表 3可得，耦合物候参数的预测模型精度(R² = 0.61, RMSE = 9.14, MAE = 7.76)高于仅用陆表参数的预测模型(R² = 0.53, RMSE = 10.75, MAE = 9.44)，说明在融合土壤湿度、植被和地形等因素的同时，加入受不同程度盐渍化影响的植被长势情况所对应的物候参数，能最大化地提升模型预测精度。将物候信息这类长时间序列反映年际变化的数据源加入到环境因子中，能更好地阐明复杂景观异质性条件下土壤盐渍化与环境之间的响应关系。但其均方根误差RMSE相差不大，可能是由于这33个参数之间存在多重共线性。Keskin等^[31]认为，当变量增加到一定程度时其精度将不再继续提升，继续增加变量会与仅有少量变量一样使模型运作变差甚至出现精度“负增长”的情况。因此，有必要对参数进行筛选来确定更合理的预测变量，本研究结合变量相对重要性，对比反演精度，经模型反复筛选最终选出23个取得较好预测结果的参数。

通过RF模型迭代100次获得环境协变量的相对重要性(relative importance, RI)，并绘制重要性玫瑰图。如图 3所示，在全部参量(图 3A)中，S1、S2、TCA、MOS和SOS 5个参数重要性较高，说明盐分指数、典型物候期和地形特征是监测土壤盐渍化的重要指标。DVI、MSAVI和SAVI等植被指数在相对重要性的排名中都接近中间和偏后的位置，说明单纯靠植被覆盖指数并不能很好地作为反映土壤盐渍化的指标，因为在盐渍化较重区域，盐生植被一般长势较好。在仅地表参量(图 3B)中，植被指数相对重要性有了一定提升，遥感反射率参数也有一定重要性，进一步说明地表可见光的反射率能有效地反映土壤盐渍化的具体分布情况。在筛选参数(图 3C)中，作为物候参数代表生物光合作用的SSI和LSI起重要作用，因为不同类型植被对土壤盐分含量响应程度不同，最终会导致生物积累量的不同，所以LSI和SSI对土壤盐渍化的表征能力总体优于其他几种物候参数。

2.3 盐渍化空间分布特征

将筛选出的23个参数作为预测变量输入RF模型，根据乔木等^[32]拟定的新疆盐渍土含盐量划分标准(非盐渍化 < 8 g/kg；8 g/kg < 轻盐渍化 < 10 g/kg；10 g/kg < 中盐渍化 < 15 g/kg；15 g/kg < 重盐渍化 < 20 g/kg；极重度盐渍化 > 20 g/kg)，得到新疆绿洲区域2019年土壤盐渍化空间分布如图 4所示。总体上看，新疆绿洲大部分区域分布的是非盐渍化土和轻盐渍化土，普遍分布在绿洲内部，含盐量大多在8 g/kg以下，中度及以上盐渍化土多分布在绿洲外围，尤其是东疆绿洲和南疆绿洲表现较为明显，即含盐量都在10 g/kg以上。位于南疆，在塔里木河流经的区域周边普遍是盐渍化土，特别是渭库南面区域为重度盐渍化土，主要是由于该区地下水位较高，加之蒸发作用强烈，盐分易累积于地表形成盐渍化土，这与野外实地调研结果一致。而在伊犁地区沿伊犁河的绝大部分区域都是非盐渍化土，据相关研究^[33]，伊犁河沿岸有大量地表水，在自然的洗盐作用下降低了土壤含盐量。

3 讨论

不同盐度的土壤具备不同光谱特征，土壤盐分随土壤样品光谱反射率增加，这使得利用多光谱波段及其衍生的光谱指数进行表土盐分监测和制图成为可能。盐度指数与植被指数在盐分预测方面表现良好，然而，每种植被类型的生境和耐受性有很大差异，为了监测土壤盐渍化，必须考虑盐生植被的光谱反射率，特别是在严重盐渍化的干旱和半干旱地区^[34]。本研究除了加入常用的植被指数，还挖掘了物候参数对土壤盐分的预测能力，发现代表生物光合作用的SSI和LSI在土壤盐分预测中起重要作用。此外，新疆地区气候干旱，蒸发强烈，土壤中的可溶性盐随水分运输积聚表层，土壤水分与盐分关系密切，而温度影响水热分布从而影响着土壤盐分水平梯度变化，反演地表温度LST和温度植被干旱指数TVDI对盐分的预测具有较大潜力。Masoud等^[35]利用Landsat影像反演地表温度，预测土壤盐渍化取得较好效果并推理土壤盐分未来可能的变化。

研究区地质环境复杂，气候环境跨度较大，导致土壤盐碱类型多样。研究中通过采样区域的物候与非物候特征得出的建模关系，运用到全疆其他区域时，通过预测发现其含盐量的总体分布特征符合实际，说明全疆的绿洲区域在总体的环境变量上存在某些内在与外在的相似性，使得通过小区域来推广到大尺度区域的研究具有一定的可行性和现实意义。考虑到具体环境的差异性，今后可以从以下两方面进行深入挖掘研究：

1) 新疆地区的物候验证数据杂乱且难以获取，对植被物候节点始末的预测也略显不足。因此，在今后的野外采样过程中应该对典型样区的植被物候数据予以阶段性记录并观测，以寻找最佳物候始末期，提高盐渍化推理制图精度。

2) 若能将研究区划分成具有不同等级的区域系统，针对具体环境筛选合适的变量，将进一步提高大尺度土壤盐分信息的反演精度。探讨不同分区的土壤盐分信息，也有助于决策者因地制宜地进行生态安全的保护和农业生产的建设。

4 结论

本文基于机器学习算法，利用少量典型样点，结合关键环境协变量预测新疆绿洲盐渍化土壤的空间分布情况，得出全疆绿洲区域盐渍化水平是伊犁平原 < 北疆绿洲 < 南疆绿洲 < 东疆绿洲。充分挖掘不同物候参数对土壤盐渍化程度的响应，发现代表生物积累量的LSI和SSI能较好地反映土壤盐渍化信息。在耦合物候参数后模型对土壤盐渍化的预测精度明显提升，决定系数从0.53提升到0.61，经模型反复迭代，最终遴选出适合该研究区的23个环境协变量，模型决定系数提升至0.73。