基于PlanetScope 影像的典型绿洲土壤盐渍化数字制图

李 科， 丁建丽， 韩礼敬， 葛翔宇，顾永昇， 周 倩， 吕阳霞

（1.新疆大学地理与遥感科学学院/智慧城市与环境建模自治区普通高校重点实验室，新疆 乌鲁木齐 830046；2.新疆大学绿洲生态教育部重点实验室，新疆 乌鲁木齐 830046；3.新疆阿克苏地区渭干河流域管理局，新疆 阿克苏 842000）

20世纪70年代，新疆地区通过第二次全国土地调查获取了土壤盐度数据，使用矢量空间数据结构记录其空间信息，这种数据记录方式难以准确地反映土壤属性的空间异质性，也难以捕获盐度信息的动态变化，无法满足新疆地区对土壤盐渍化管理和预防的需求［7］。研究证明，光谱学知识结合遥感技术可以快速且定量地获取土壤盐度信息，并成功用于绘制高分辨率的土壤盐度图，提供了更全面的土壤盐渍化空间信息［8］。盐渍化土壤相较于非盐渍化土壤，在可见光到短波红外波长范围内具有更高的光谱反射率，随着土壤盐度增加，光谱反射率进一步提高，在近红外波长附近出现具有诊断性的吸收带［9-11］。因此，归一化盐度指数（Normalized difference salinity index，NDSI）和多种盐度指数（Salinity index，SI）被广泛用于评估土壤盐度。植物对土壤盐碱变化的耐受，通过生长状况进行响应，长势不好的植被光合作用减弱，叶绿素含量变少，可见光范围内反射率增强，近红外范围内反射率降低［12］。因此，归一化植被指数（Normalized differential vegetation index，NDVI）和土壤调节植被指数（Soil adjusted vegetation index，SAVI）等多种植被指数（Vegetation indices，VIs）被广泛用于评估土壤盐度的间接指标。Gorji 等［13］基于Landsat8 OLI 和Sentinel-2A MSI影像，成功提取了与土壤电导率相关的SI，并结合多元线性回归技术成功绘制了乌尔米耶湖流域土 壤 盐 度 图。Fernandez-Buces 等［14］基 于Lansat ETM、机载近红外和实测光谱数据获取了不同植被和裸地的反射率，提取并调整了NDVI，并计算了综合光谱指数（Combined spectral response index，COSRI），以制作出高精度的原墨西哥特斯科科湖区域盐度图。Ding 等［15］使用近端遥感装置EM38 测量表观电导率数据，基于Landsat TM 影像提取了SI、NDSI、亮度指数（Britness index，BI）、NDVI 和差值植被指数（Difference vegetation index，DVI），将表观电导率数据和遥感指数信息应用于回归克里金插值法，成功评估了土壤盐度并获得了满意的结果。

20 世纪末，机器学习算法被广泛使用，在土壤科学领域，基于机器学习算法绘制了多种数字土壤地图，如Farifteh 等［16-17］基于随机森林方法（Random forest，RF）开发了绘制土壤有机碳浓度和储量数字图的框架。Ma 等［18-19］基于RF 和极端梯度提升等机器学习方法绘制了渭库绿洲土壤盐度图，有效解决了土壤和环境因子间的非线性问题。Wang 等［20］在绘制不同区域的土壤盐度图时，比较了5 种机器学习算法，发现基于机器学习算法进行土壤盐度预测研究时需要用户定义的参数更少，拥有更高的计算效率，能够处理数值序数和离散预测变量等优点。然而，在构建高效、准确和泛化程度高的土壤盐度预测模型时，特征变量的选择至关重要。Allbed等［19-21］的研究结果表明植被指数和盐分指数等特征变量可以有效用于土壤盐渍化监测研究，但是这些特征变量均是通过波段间的运算得到，造成了特征变量间不同程度的信息冗余。最相关最小冗余（Max-relevance and min-redundancy，mRMR）特征选择算法具有较低的复杂度，与当前常用的特征选择算法相比具有原理简单和对数据要求低的优点，mRMR 具有的优势提升了土壤盐度预测模型的预测效率［22］。Peng 等［23］使用Pearson 相关分析筛选相关性高的特征变量进行土壤盐渍化监测，并没有考虑选择特征间存在的信息冗余。马慧琴等［24］在筛选特征用于构建遥感监测模型时，对比了mRMR 与常用相关分析方法，研究结果表明mRMR 方法比常用相关分析方法更具优势。

目前，遥感影像的中低空间分辨率仍然是影响土壤盐度图精度的主要因素［25］。当受盐影响区域小于像素大小时，会限制土壤盐度图对土壤盐度进行准确评估的能力，空间分辨率低于5 m 的多光谱传感器正成为绘制土壤盐度图不可或缺的装置［26］。Bouaziz 等［27］证实了基于MODIS terra 数据提取的光谱指数与实测数据能有效绘制土壤盐度图，但较低的空间分辨率是造成二者相关性较差的主要原因。Eldeiry 等［28］比较了IKONOS 与Landsat 影像在农业地区评估土壤盐度的性能，结果表明高空间分辨率能获取更准确的评估结果。自2016 年以来，Planet公司运营的地球成像小型卫星星座PlanetScope，能够以较高的空间分辨率（＜5 m）和时间分辨率（1～2 d）获取影像数据。在2019 年4 月PlanetScope 携带了第三代传感器PSB.SD 能获取可见光近红外范围内的8 个波段［29］。PlanetScope 获取的影像与Landsat8 OLI［30］和Sentinel-2A MSI［31］相比，具有更高的空间分辨率和时间分辨率，而且国内利用PlanetScope小型星座卫星影像进行土壤盐度制图的研究并未广泛展开。促进高空间和高时间分辨率的遥感影像用于土壤盐度监测研究，能获取更翔实的空间信息、动态变化的盐度信息和更精确的评估结果。

本文以渭干河-库车河流域三角洲绿洲（简称渭库绿洲）为研究区，基于PlanetScope 卫星的原始波段数据和波段运算所得的光谱指数，共21个特征变量。采用mRMR 方法筛选重要特征作为自变量输入Bagging 回归算法，建立了土壤盐度预测模型，成功绘制了研究区高分辨率土壤盐度分布图。这为新疆地区制定管理改善土壤盐碱策略，维持绿洲可持续发展和保障粮食安全提供重要参考。

1 材料与方法

1.1 研究区概况

渭库绿洲，属于南疆绿洲区。地势北高南低、地貌复杂多样，气候干燥、蒸发强烈、降水稀少，属于典型的温带大陆性干旱气候［15］。复杂的地形地貌、极端的气候条件、农业上排水系统不健全、灌溉技术粗放等问题［32］，导致研究区存在不同程度的盐渍化。绿洲外围植被稀疏，生态系统脆弱，主要包括芦苇（Phragmites australis）、柽柳（Tamarix chinensis）、碱蓬（Suaeda glauca）和骆驼刺（Alhagi sparsifolia）等植物。实际研究区位于库车市绿洲外围农田荒漠交错带，地理位置介于83°23′～83°29′E、41°23′～41°27′N之间（图1）。

图1 研究区概况图Fig.1 Overview map of the study area

1.2 样品采集与处理

实验团队于2021 年6 月23—24 日在研究区进行了实地调查，结合以往的实地调查成果和采样经验，在农田荒漠交错带选择了84个具有代表性的采样点。用便携式GPS 记录采样点的地理位置，以采样点为中心，在距离采样点1 m 的范围内，以0～10 cm 的深度使用土壤取样器均匀的采集5 个样品，然后将5 个样品混合在一起，取混合后的样品约500 g，装入有标签的防水密封袋。在实验室里，所有的土壤样品均经过风干、研磨和过筛（2.0 mm）等一系列处理。取20 g 处理过的土样和100 mL 纯水充分振荡，静置24 h。在25 ℃的温度下，用多参数测量仪（WTWinoLab®Multi3420 set B，WTW GmbH，Germany）测量静置后的浸出液，得到电导率（EC）值。由于土壤中的盐分溶解在浸出液中，浸出液的含盐量与25 ℃下电极的导电能力呈正相关关系，通过测量土壤浸出液的EC 值，可以作为评估土壤盐度的参考指标［33］。

1.3 遥感影像预处理与光谱指数提取

PlanetScope 是通过多组发射的单独卫星组成的卫星星座，目前约有130 颗CubeSat 3U（10×10×30 cm）卫星［34］。目前已经搭载了第3代传感器PSB.SD，幅宽为32.5×19.6 km，能采集8 波段的多光谱影像，空间分辨率为3 m。这些卫星大多位于太阳同步轨道，所有卫星组成的卫星星座PlanetScope 拥有每天对地球成像的能力。本研究使用的Planet-Scope 波段参数见表1，使用蓝波段（B2）、绿波段Ⅱ（B4）、红波段（B6）、红边波段（B7）和近红外波段（B8）进行了后续研究。本文使用的是PlanetScope Level-3B数据，过境时间为2021年7月24日，Level-3B 产品使用6S 辐射传输模型和MODIS 数据作为辅助对PlatnetScope获取的遥感影像进行大气校正，并经过了传感器校正和正射校正。

HPLC切换波长法同时测定健脾止泻宁颗粒中盐酸小檗碱和黄芩苷的含量 ……………………………… 黄传俊等（10）：1324

表1 PlanetScope影像光谱波段参数Tab.1 Spectral band parameters of PlanetScope images

将光谱指数作为预测土壤特性模型的输入变量时，能有效降低光谱反射率误差［35］。通过波段运算得到的盐度指数与植被指数被证明能准确反演土壤盐分［36］。本研究使用PlanetScope 原始波段数据计算光谱指数，计算公式如表2所示。

表2 基于PlanetScope影像光谱指数计算Tab.2 Spectral indices calculation based on PlanetScope images

1.4 mRMR特征变量选择算法

mRMR 是一种过滤型特征选择算法，其算法原理简单且对数据的要求低，在土壤盐度监测领域还没有相关的研究探讨该算法的应用。因此，本研究结合mRMR 方法构建了土壤盐度预测模型，为土壤盐度监测领域提供了参考。mRMR 考虑了特征之间相似性造成的数据冗余问题，通过计算特征与离散化电导率数据类别之间的互信息，特征与特征之间的互信息，来度量它们之间的相关性。然后，保留与离散化电导率数据类别之间互信息值较高的特征，同时去除特征之间互信息值较高导致数据冗余的特征［44］。任意2 个特征x和y的互信息计算公式为：

式中：I(x,y)为特征x和y的互信息；P(x,y)为特征x和y的联合概率密度函数；P(x)为特征x的边缘概率密度函数；P(y)为特征y的边缘概率密度函数。特征与离散化电导率数据类别之间的互信息计算公式为：

式中：D(S,p)为特征与类别之间的平均互信息；S为变量集合；I(xi,p)为特征i与类别p之间的互信息；xi为特征i；p为离散化电导率数据类别。特征与特征之间的互信息计算公式为：

式中：R(S)为特征与特征间的平均互信息；I(xi,xj)为特征i与特征j之间的互信息；xi为特征i；xj为特征j。结合公式（2）与公式（3）得到最大相关最小冗余的判断条件商标准：

式中：φ(D,R)为互信息商标准；D为特征与类别之间的互信息；R为特征与特征间的互信息。mRMR 算法通过Matlab 2019b运行fscmrmr函数实现，fscmrmr函数自动将电导率数据离散化处理。

1.5 建模策略及评价

1.5.1 Bagging回归算法装袋算法（Bootstrap aggregating, Bagging）是Breiman［45］提出的一种并行集成学习方法，将多个单一的预测器聚合在一起，对结果进行平均处理，进行分类或者回归预测相关的研究。Bagging 在训练集中进行Bootstrap 采样，即从m个样本中进行随机的、有放回的采样，得到n个与原始训练集样本个数相等的采样集。单一预测器基于这些采样集进行独立的训练，在进行回归分析时，将n个预测器的均值作为综合预测结果，过程如图2 所示。Bagging 回归算法有效降低了训练集的方差，避免预测结果出现过拟合。通过Python 3.9导入机器学习库sklearn构建Bagging回归模型。

图2 Bagging回归算法预测过程Fig.2 Prediction process of Bagging regression algorithm

1.5.2 模型构建及评价并非所有特征都能提升Bagging 回归算法预测性能，不相关特征和冗余特征往往会降低模型预测的精度和效率。为了降低负面影响，探索PlanetScope 多光谱数据预测土壤盐分的潜力，所有特征构成集合Ⅰ，mRMR 算法选择前30%的特征构成集合Ⅱ，基于集合Ⅰ和集合Ⅱ，开发了两种模型Model-Ⅰ和Model-Ⅱ（表3）。将采集的84 个土壤样品进行随机划分（在Python 3.9 中test_size=0.3 用于划分数据集）。所有样品的70%，约58 个土样被划分为训练集，用于训练模型；所有样品的30%，约26 个土样被划分为验证集，用于模型的验证。为了确保模型的稳健性，提高模型的适用性，降低随机划分样本对模型预测性能的影响，本研究随机划分土壤样本3 次，获得3 组训练集和验证集（A、B 和C）。为了确保模型的可比性，每次划分训练集和验证集时，对于土壤样本与集合Ⅰ和集合Ⅱ构成的数据集，划分模式保持一致。

表3 PlanetScope数据建模策略细节Tab.3 Details of the PlanetScope modeling strategies

为了评估Bagging 回归模型预测EC 的能力，使用决定系数（R2）、均方根误差（RMSE）和四分位数的相对预测误差（RPIQ）3个指标进行衡量。模型预测的准确性与R2、RPIQ 成正比，与RMSE 成反比，即R2与RPIQ 越大，RMSE 越小模型预测性能越好，计算公式如下：

式中：p^i为预测值；pi为实际测量值；pˉ为测量值的平均值；n为采样点数量；IQR 为验证集中测量值的四分位间距；RMSEP为验证集的均方根误差；Q3 为上四分位值；Q1为下四分位值。

2 结果与分析

2.1 土壤电导率描述性统计

EC 作为土壤盐分的参考，EC 的描述性统计结果如图3 和表4 所示。整个EC 的范围为0.35～124.90 dS·m-1，平 均 值 为38.04 dS·m-1，标 准 差 为33.13 dS·m-1，变异系数为0.87，为中等变异。EC 描述性统计结果显示，土壤样本随机划分3 次得到A、B 和C 数据集，它们的训练集和验证集分布与整个数据集的分布相似。因此，训练集和验证集的样本可以代替整个数据集的样本进行建模和验证。

表4 数据集A、B和C与全集EC描述性统计Tab.4 Descriptive statistics of datasets A,B,and C compared to the entire dataset EC

图3 电导率统计图Fig.3 Statistical diagram of electrical conductivity

2.2 特征变量优选

集合Ⅰ由5 个波段反射率数据B2、B4、B6、B7 和B8，4 个植被光谱指数RVI、EVI、SAVI 和NDVI，12个土壤盐度光谱指数S1、S2、S3、S4、S5、S6、SI、SI1、SI2、SI3、BI和NDSI组成。从集合Ⅰ中选择最能反映土壤盐渍化发生发展状况的特征变量构建土壤盐渍化预测模型，能有效提升模型效率和预测精度。目前mRMR 特征选择方法还未广泛用于土壤盐渍化监测领域，本研究使用Matlab 2019b 统计和机器学习工具箱执行fscmrmr 函数，通过mRMR 算法计算互信息，并依据互信息商标准保留最相关特征，过滤冗余特征。经过互信息商标准选择前30%特征作为集合Ⅱ，包括1 个波段反射率数据B4，1 个植被光谱指数EVI，4 个土壤盐度光谱指数S6、S1、SI2和NDSI，图4 为波段反射率与光谱指数特征重要性。

图4 PlanetScope波段反射率和光谱指数的特征重要性Fig.4 Importance of PlanetScope band reflectance and spectral indices

2.3 预测模型和性能评估

为了提高模型的预测精度和泛化程度，使用基于步长为1 的网格搜索，对模型影响较大的基础学习器数量进行了1000次连续迭代调整，根据模型R2值确定模型参数最优值。在进行参数优化时，其余参数为默认值。图5显示了模型R2与模型参数在连续迭代过程中的变化关系，当R2最大时确定基础学习器数量数值，重复此过程直到确定参数最优值，Model-Ⅰ模型基础学习器数量为4，Model-Ⅱ模型基础学习器数量为9。

图5 模型参数调整Fig.5 Model parameters adjustment

基于优化后的参数，将集合Ⅰ和集合Ⅱ随机划分的训练集作为自变量输入Bagging 模型，构建EC预测模型Model-Ⅰ和Model-Ⅱ，利用随机划分的验证数据集对模型的预测性能进行验证。模型Model-Ⅱ的R2与RPIQ均高于Model-Ⅰ，RMSE均低于Model-Ⅰ，这说明mRMR 算法可以选择更重要的特征代替全部特征，降低数据的复杂性和冗余性，建立更准确和更高效的模型（表5）。在随机划分的三组数据集A、B 和C 中，Model-Ⅱ与Model-Ⅰ相比，R2平均提高了27%，RPIQ 平均提高了20%，RMSE 平均降低了16%。mRMR 在选择特征数量较少的条件下，能保留最相关最小冗余特征，能最大化特征包含的有用信息。在集合Ⅱ中光谱指数EVI 和S6 重要性最高，表明EVI对土壤盐渍化引起的植被状况与光谱反射率变化敏感，S6 对盐渍土表面状况与光谱反射率变化敏感。光谱指数EVI 与S6 在土壤盐渍化监测中具有巨大潜力。

表5 模型预测性能比较Tab.5 Comparison of model prediction performance

2.4 土壤盐度空间分布

基于模型Model-Ⅰ和Model-Ⅱ，使用随机划分的数据集A、B 和C 绘制了3 张高分辨率的土壤EC 分布图，对模型Model-Ⅰ和Model-Ⅱ的3张预测图进行平均，得到图6a 和图6b。图6a 和图6b 产生了相似的土壤EC 空间分布，Model-Ⅱ模型产生了最理想的结果，预测的土壤EC范围（0.70～100.83 dS·m-1）更符合实际测量的土壤EC 值。图6b 相较图6a，更准确的显示了农田区域与荒漠区域土壤EC 分布的细微变化，更符合该区域土壤盐度分布。受地形影响，土壤盐度高值主要分布在荒漠区域；受地形和灌溉的影响，土壤盐度低值主要分布在农田区域。

图6 土壤电导率空间分布Fig.6 Spatial distributions of soil electrical conductivity

3 讨论

PlanetScope 数据结合mRMR 特征选择方法首次用于绘制高空间分辨率的土壤EC 图，土壤盐渍化监测领域缺乏类似的研究作为参考，因此借鉴了马慧琴等［24］和Jing 等［44］的研究，通过比较模型Model-Ⅰ和Model-Ⅱ平均预测性能（表5），得出与其相似的研究结果，mRMR 方法能优化建模特征并改善模型预测性能。基于PlanetScope 数据建立的模型Model-Ⅰ和Model-Ⅱ获取的预测值如图7 所示，Model-Ⅱ获取的预测值与实测值更接近1:1 拟合线，表明Model-Ⅱ能获取更准确的预测结果。基于PlanetScope影像低于5 m 的空间分辨率，能避免如Scudiero 等［46］使用Landsat7 多光谱影像进行土壤盐度评估时，结果难以表征土壤盐度的空间变异性等情况。也能避免使用MODIS系列、Landsat系列或Sentinel 系列数据进行土壤盐度监测时，预测结果难以具备指导意义等情况［47］。基于此，使用PlanetScope数据结合mRMR特征选择方法建立的模型Model-Ⅱ对研究区土壤EC 进行预测并绘制成图，避免了由于空间分辨率不足对土壤盐度监测和绘制带来的限制。

图7 模型Model-Ⅰ和Model-Ⅱ预测性能比较Fig.7 Comparison of prediction performance of model Model-Ⅰand Model-Ⅱ

本研究将光谱指数作为自变量输入Bagging 模型，改善模型的预测性能，提高模型定量描述土壤盐度的能力，丁建丽等［48］与曹雷等［49］的研究结果证实这是可行的。并非所有光谱指数都能用于优化模型，其中不相关和冗余的特征变量反而会起反作用。本文通过mRMR 方法对所有特征变量进行了重要性排序，将排序后的特征变量作为自变量输入模型Model-Ⅱ中。当选择的特征变量数小于30%时，模型的预测精度较低；当选择的特征变量数超过30%时，模型的预测性能略低基本保持稳定。

重要性排序中EVI 重要性最高，远远高于其他波段与光谱指数（图4），对土壤盐渍化具备更多的解释能力。这可能是因为植物叶面积和冠层结构的发育受到土壤盐渍化抑制更明显，而EVI 对于植物冠层结构变化更敏感，故EVI 在建模时更重要。Rao 等［50］的研究表明，盐渍土具备更高的光谱反射率，且在可见光和近红外波长范围，盐渍化程度与光谱反射率成正比，在盐渍化程度较高的农田荒漠接壤处与荒漠区域可以直接利用PlanetScope原始波段反射率对土壤盐度进行建模分析。mRMR选择的特征集合Ⅱ中仅包含原始波段B4，故B4相较于其他原始波段，对于土壤盐渍化更敏感，更适合用于描述土壤盐度信息。在本研究中除植被覆盖度较高的农田以外，均能使用对土壤盐度敏感的盐度指数定量描述土壤盐渍化情况，mRMR 选择的前30%特征变量中，S6、B4、S1 和NDSI 均为盐度指数，故盐度指数对于建立土壤盐度预测模型必不可少。冯娟等［51］、Metternicht 等［52］和丁建丽等［53］的研究结果证明了这些观点是合理的。

土壤样本的随机划分和模型参数的设定使模型的误差出现累积。随机划分的训练集和验证集与整个数据集的分布相似，可以代替整个数据集的样本进行建模和验证。为避免样本随机划分为模型带来的不确定性，将随机划分3 次得到的数据集A、B 和C 用于训练模型和验证模型，最终对评价指标取均值来评估模型，增加了模型预测的稳健性，Jing 等［44］的研究结果证明这是有效的。虽然多次随机划分样本能有效降低模型预测的不确定性，但是本研究并没有确定最有效的随机划分次数。在设定模型参数时，为了降低模型参数对预测性能的影响，基于步长的网格搜索，对模型影响较大的基础学习器数量进行了连续迭代调整，结果如图5所示。基于步长的网格搜索是否是参数优化的最佳方式，这也存在不确定性。除上述两点，人类活动也增加了模型预测的不确定性，图6a和图6b左上角区域显示，通过Model-Ⅰ和Model-Ⅱ预测的农田区域土壤盐度被高估，这与实际情况不符。通过野外调查验证，该区域为了达到节水增产的目的，推广了膜下滴灌技术。覆膜区域光谱反射率不同于裸土区域，整体偏高，由于多光谱影像的光谱分辨率不足，使得盐碱地与覆膜地难以进行区分，在未来的研究中将尝试利用高光谱技术来消除农田区域覆膜给土壤盐度监测带来的影响。此外，本研究仅考虑了Bagging 回归模型，马国林等［54］的研究结果表明，仅通过一个模型难以取得最佳预测效果。因此，在未来的研究中将构建多种模型对土壤盐度进行预测，并对比模型性能选择最适合本研究数据集的模型。

4 结论

本研究基于PlanetScope 影像原始波段数据、植被光谱指数和土壤盐度指数共21个特征变量，使用mRMR 方法筛选前30%特征变量并输入Bagging 回归模型，建立Model-Ⅰ和Model-Ⅱ土壤盐度预测模型，成功绘制了研究区3 m 分辨率的EC 空间分布图，为准确监测随时空动态快速变化的土壤盐度信息提供参考。得出以下结论：

（1）将PlanetScope 影像原始波段数据、植被光谱指数和土壤盐度指数共21 个特征变量作为自变量建立的Model-Ⅰ模型，基本可以实现对EC 的预测（验证集平均R2=0.52，平均RMSE=21.32 dS·m-1，平均RPIQ=2.68）。

（2）将mRMR 筛选后的特征变量作为自变量建立的Model-Ⅱ模型，降低了预测EC 的特征维度，提高了模型的预测效率和准确性（验证集平均R2=0.66，平均RMSE=18.00 dS·m-1，平均RPIQ=3.21）。Model-Ⅱ与Model-Ⅰ相 比，EVI、S6、B4、S1、SI2 和NDSI 对农田荒漠交错带土壤盐渍化监测具有重要作用。