基于国产资源一号02D卫星和机器学习算法的钦州湾滨海土壤盐分反演

田义超,郑丹琳,张 强,卢 芳,黄友菊,陶 进,张亚丽,林俊良,姚贵钊,姚媛元(.北部湾大学资源与环境学院,北部湾海洋发展研究中心,广西 钦州 55000；2.北部湾大学,广西北部湾海洋环境变化与灾害研究重点实验室,海洋地理信息资源开发利用重点实验室,广西 钦州 55000；.广西壮族自治区自然资源遥感院,广西 南宁 50028)

土壤盐渍化又称土壤盐碱化,指在自然和人为作用下土壤表层盐分含量不断增加,以至超过某一限度的现象和过程.近年来随着人为破坏、自然环境变化以及海平面不断上升,土壤盐碱化成为了粮食产量不断减少的诱因之一.据估计,全球约有2.3 亿hm2的灌溉土地和4500 万hm2的农田受到盐碱化的影响[1].因此,为了农田的可持续发展,对土壤盐分进行估算、监察,实现大范围土壤盐分制图刻不容缓.

国内外利用卫星遥感进行土壤盐渍化监测研究始于20世纪70年代[1-2].1992年,Dwivedi 等[3]通过 TM 影像研究了印度恒河冲积平原土壤盐渍化问题,他们发现不同波段组合会对土壤盐渍化的估算结果产生影像.随着科技的进步,研究方法日趋成熟,研究者们开始尝试构建不同的土壤盐分反演模型,主要分为两类,一是利用统计方法进行建模,主要的统计方法有偏最小二乘回归(PLSR)[4]、逐步多元线性回归(SMLR)[5]等,其中偏小二乘回归被证明相较于其他回归方法能够提供更好的预测精度[6-9].二是利用机器学习方法进行建模,主要有神经网络(NN)[10]、支持向量机(SVM)[11]、随机森林回归(RFR)等[12],其中神经网络算法[13]由于其收敛速度快等优势被广泛使用.Wang 等[14]基于不同的数据源创建Bootstrap-BP 土壤盐分反演模型,认为与以往研究中光谱指标的线性回归模型相比,Bootstrap-BP 模型的预测精度有很大提高.2017年Vermeulen 等[15]分别建立神经网络、支持向量机、决策树(DT)和随机森林4 个模型对南非Vaalharts 和Breede 区域的盐分进行遥感监测,结果表明利用高程及其导数,结合地统计学和最大似然法,能够很好的识别盐渍土.Wang 等[16]采用偏最小二乘回归、卷积神经网络(CNN)、支持向量机和随机森林四种预测建模方法对土壤盐度进行建模和估计.2020年Wang 等[17]以新疆干旱荒漠地区土壤盐分为研究对象,比较了两种遥感数据源(OLI 和MSI)在土壤盐度监测中的性能,发现MSI 数据集与Cubist 模型的组合总体上具有更好的模型性能和精度指标.

随着科学技术的发展,越来越多的机器学习方法被引入到土壤盐分的反演中,但较少涉及到集成学习.集成学习主要有Bagging 算法[18]、Boosting 算法[19]和Stacking 算法等,其中Boosting 算法可以将弱学习者转换为强学习者以提高模型的准确性,该算法有极端梯度Boosting(XGBoost)、分类Boosting(CatBoost)、自适应Boosting(AdaBoost)、轻梯度Boosting(LightGBM)和梯度提升决策树(GBDT).王飞等[20]以我国新疆为研究区,比较了5 种机器学习算法构建的土壤盐分含量反演模型,认为随机梯度增进算法(SGT)预测精度相对最高.2021年Ma 等[21]人通过对新疆塔里木盆地中北部渭-库绿洲的土壤盐渍化进行研究,发现XGBoost 模型的预测精度优于分类与回归树和随机森林.这类算法具有更准确的预测性能,在模型预测能力优于参数回归模型,被广泛应用于空气质量监测[22]、灾害预测[23]、水质反演[24]等,但目前在研究土壤盐分方面,利用集成学习方法对土壤盐分反演的相关研究较少.

滨海盐渍土主要是由沿海地区盐渍淤泥发育而成,积盐过程的盐分补给方式主要是海水浸渍以及海水倒灌[25].钦州湾位于广西沿海中部,属于我国北部湾的一部分,是广西重要的海湾,其作为南亚热带滨海平原区,自然环境复杂多变,生态脆弱性较高[26-27].同时,随着城市化、港口建设和填海等人为活动的加剧,钦州湾滨海区域环境压力逐渐增加,土壤质量、农业生产和区域生态系统的安全受到严重影响.但目前关于该地区土壤盐分空间分布状况及分布规律还尚无文献可循.有鉴于此,本研究以广西钦州湾滨海区域土壤为研究对象,以国产资源一号02D 卫星作为数据源并从中提取特征变量.利用AdaBooost、LightGBM、XGBoost、RFR 以及CatBoost 五种不同机器学习算法,设置了5 种不同变量组合来研究不同机器学习算法及不同变量组合对模型精度的影响,并基于最优算法的最优变量组合对研究区域内的土壤盐分进行预测.本研究的3个主要目标如下:(1)比较不同机器学习模型在估算土壤盐分方面的性能;(2)确定不同变量组中特征变量对机器学习模型精度的影响;(3)验证不同特征变量对模型重要性的贡献.

1 材料与方法

1.1 研究区概况

研究区选取位于钦州湾滨海地区的一处水田,介于21°38′10″～21°38′40″N,108°46′34″～108°50′25″E 之间,东西贯长约6.621 km,南北贯长约1.051 km,面积约352.4 hm2.该区域属于南亚热带季风气候,夏季湿热,冬季干冷,一年中,气温最高在35℃左右,最低在11℃左右,全年平均气温在22℃左右,年均降水量1658 mm,年总日照为1673 h[28].该区域位于广西沿海中部的钦州湾,钦州湾是广西海岸带中段,其西北为台地、丘陵,茅岭江出海口,东北为饮江平原,东南为犀牛角平原.钦州湾滨海地区的土壤成土母质主要是花岗岩和砂页岩,在地带性气候的作用下,陆地部分形成了典型浅海赤红壤,土壤有机质较为丰富[29].研究区具体位置如图1所示.

图1 研究区及采样点空间分布Fig.1 Distribution of the study area and sampling points

1.2 地面数据采集与处理

本文于2021年6月中旬对研究区进行实地调查与土壤采样.在采样过程中,设置了10 m×10 m 的样点格网,以格网交叉点作为采样点,利用GPS 确定每个采样点的经纬度坐标,通过手持式多参数土壤速测仪对每个样点的土壤电导率进行多次测量,待稳定后记录,取测量值均值作为各样点电导率,保证了结果的可靠性,最后共采集94 份土壤样本.参考前人的研究[30],将各样方点实测电导率转换为土壤盐分含量,并导入ArcGIS10.2 中进行矢量化.电导率转换为土壤盐分含量公式如下:

式中:SS 为土壤盐分含量,g/kg;EC 为电导率,µS/cm.

1.3 遥感数据与预处理

资源一号02D 卫星于2019年9月12日发射,由自然资源部主持建造,是资源一号02C 星的接续星,属于中等分辨率遥感业务卫星[31].该卫星具有幅宽大、光谱分辨率高的突出特点,海岸蓝、黄、红边等多光谱谱段及高光谱影像可提供丰富的地物光谱信息,拓展了资源调查与监测的手段,为及时掌控自然资源质量、分布情况以及变化趋势提供了支撑.资源一号02D 卫星采用太阳同步轨道,轨道高度778km,轨道倾角为98.5°,回归周期为55d,重访周期为3d.卫星配置了9 谱段可见近红外相机和166 谱段高光谱相机,其中,全色谱段与多光谱波段空间分辨率分别为2.5 和10m;高光谱相机的光谱范围400～2500nm,空间分辨率30m[32].

结合野外采样时间以及遥感影像云量等因素,本研究所选取的数据源为2021年6月16日获取的ZY1-02D VNIC 影像,包括2.5 m 全色波段与10 m 多光谱波段,数据来源于中国资源卫星应用中心(https://data.cresda.cn/).为使影像尽可能反映地表实际信息,对该数据进行预处理,过程包括辐射定标、大气校正、正射校正、几何校正、图像融合以及裁剪等.其中大气校正采用FLAASH 大气校正法,影像全色波段和多光谱波段的融合采用Gram-Schmidt Pan Sharpening 法.对于对应土壤采样点的光谱反射率本研究提取了所有土壤样点在ZY1-02D VNIC 影像上该采样点的平均反射率,作为遥感反射率来反演土壤盐度.以上步骤在ENVI 5.3 中完成.

1.4 特征提取

在ENVI5.3 以及ArcGIS10.2 等遥感软件的支持下,本研究共选择了37 个特征变量构建盐分反演模型,其中包括遥感影像的4 个原始波段、9个植被指数、13 个盐分指数、3 个主成分波段以及8 个纹理特征.表1 显示了37 个特征变量的计算方式.

表1 反演参数计算公式Table 1 Calculation formula of inversion parameters

(1)原始波段.遥感影像原始波段承载大量光谱信息,是地物识别的基础,本研究提取经预处理后的ZY1-02D 遥感影像单波段光谱信息,分别为蓝波段(B)、绿波段(G)、红波段(R)、近红外波段(NIR).

(2)植被指数.植被指数是对地表植被状况的有效、简单和经验的度量[33],可以有效地反映植被信息.本研究依据研究区植被覆盖特点及前人研究成果提取了9 种植被指数,分别为归一化植被指数(NDVI)、增强植被指数(EVI)、差值植被指数(DVI)、比值植被指数(RVI)、大气阻抗植被指数(ARVI)、广义差分植被指数(GDVI)、绿波段比值植被指数(GRVI)、联合光谱响应指数(COSRI)、绿色大气阻抗指数(GARI)、土壤调节植被指数(SAVI).

(3)盐分指数.盐分指数可以有效反映土壤盐分信息,表征土壤盐渍化信息.根据前人研究成果,提取了11 种盐分指数,分别为:盐分指数(S1)、盐分指数(S2)、盐分指数(S3)、盐分指数(S4)、盐分指数(S5)、盐分指数(S6)、盐分指数(SI-T)、盐分指数(SI1)、盐分指数(SI2)、盐分指数(SI3)、盐分指数(INT1)以及归一化盐分指数(NDSI).

(4)主成分波段.由于遥感影像波段数量较多,造成一定的信息冗余,于是本研究使用主成分分析变换特征,提取影像主成分变换的前3 个主成分作为特征变量,分别为主成分分析波段 PC1、PC2、PC3.

(5)纹理特征.纹理特征包含了丰富的影像信息,被广泛运用于高分辨率影像信息识别中.本研究在ENVI5.3 软件支持下对ZY1-02D 数据进行主成分变换后,通过Co-occurrence Measures 工具采用灰度共生矩阵(GLCM)的方法对第一主成分波段进行纹理特征提取,共提取ZY1-02D 影像PCA1 的8 个纹理特征值,即均值(MEA)、方差(VAR)、同质性(HOM)、对比度(CON)、相异性(DIS)、熵(ENT)、二阶距(SEC)、相关性(COR).

1.5 研究方法

1.5.1 研究路线 在资源一号卫星数据与实地采样数据的基础上,本研究通过提取ZY1-02D遥感影像的波段信息、植被指数、盐分指数、主成分波段以及纹理特征,构建了能够反演土壤盐分的特征变量数据集.通过比较AdaBooost、LightGBM、XGBoost、RFR以及CatBoost 五种机器学习算法以及不同变量组合的预测精度,得出最优反演模型来估算研究区土壤的盐分含量.具体技术路线如图2所示.

图2 土壤盐分反演技术路线Fig.2 Technical route for soil salinity inversion

1.5.2 机器学习方法 Adaboost 是Freund 等在1997年提出的,他们将加权投票的相关研究成果[49]与在线分配问题结合,并在Boosting 问题框架下进行对应推广,得到了著名的AdaBoost[50]算法.该算法的优点在于使用简单的弱分类器,通过不断的进行权值分配与迭代最终得到强分类器,可以深挖分类器的能力,也不存在过拟合的现象,因此能显著提高学习精度.

RFR 是Breiman[51]提出的一种机器学习算法,该算法可以有效处理分类和回归问题.随机森林集成了所用的分类投票结果,投票最高的结果作为最后输出的结果.目前,大量的理论研究已经证明,RFR具有较高的预测精度,能够很好地容忍异常值和噪声,并且不容易过拟合.

Xgboost 是2016年开发的一种在GBDT 算法基础上加以改进机器学习算法[52],现已广泛用于解决数据挖掘问题.XGBoost 通过二阶的泰勒展开表达式优化损失函数,并且在目标函数之外加入了正则项,整体求最优解,避免了过拟合,从而提高了预测的精度[53].

LGBM 是微软在2017年提出的一种梯度提升算法,是基于传统GBDT 算法改进的模型[54].LGBM通过梯度单边采样技术去除梯度小的数据,用剩余的数据对信息增益进行估计,并且利用独立特征合并技术减少特征的数量.

CatBoost 是俄罗斯Yandex 学者于2017年提出的一种新型开源机器学习库,它由Categorical 和Boosting 组成[55],是一种以对称决策树基学习器实现的模型[56].

1.5.3 模型精确度评估 本研究采用了决定系数R²以及均方根误差RMSE 来评估机器学习模型在土壤盐分估算中的精度.其中决定系数R²指的是预测值与实测值的拟合程度,当R²的值越接近1 即代表模型的精度越高.RMSE 主要用于模型验证,反映预测值与真实值之间的偏差.数值越小,模型精度越高[57].决定系数R²和RMSE 的计算公式如下所示:

2 结果与分析

2.1 研究区土壤盐度特征

本研究所调查的94 个采样点均均匀分布在研究区范围内.根据前人研究成果[58],将区域土壤盐渍化程度划分为非盐渍化、轻度盐渍化、中度盐渍化、重度盐渍化以及盐土5个等级,分级标准见表2.由表3 可知,本研究94 个采样点的土壤盐度范围为0.740～10.352 g/kg,74 个建模样本与20 个检验样本均值分别为1.739 和1.465 g/kg,标准差为1.927 和0.958 g/kg,变异系数(变异系数反映样点值的离散程度,变异系数<0.1 为弱变异性;0.1<变异系数<1 为中等变异;变异系数>1 为强变异)分别为1.10 和0.65,属于强变异与中等变异.根据土壤盐渍化程度分级标准,统计实测样点土壤含量数据,建模集与检验集中非盐渍化样点个数分别占37.838%和40%,轻度盐渍化样点个数占48.649%和45%,中度盐渍化样点个数占5.405%和10%,重度盐渍化样点个数2.703%和5%,盐土样点个数5.405%及0%,占比最大均为轻度盐渍化.

表2 土壤盐渍化程度分级标准Table 2 Grading standard of soil salinization degree

表3 采样点统计分析Table 3 Statistical analysis of sampling points

2.2 特征变量与盐分的相关性

根据采样点的坐标信息,在 ArcGIS10.2 及RStudio 软件的支持下,从反演参数中提取对应实测样本点的各特征变量值,使用R 语言中的Cor 函数对实测样本土壤盐分含量与对应的特征变量值进行皮尔森相关系数计算.实测点土壤盐分含量与反演参数的相关性见表4.从相关性结果来看:(1)通过计算植被指数变量组(r=0.062),原始波段变量组(r= 0.019),盐分指数变量组(r=0.078),遥感影像数据指数变量组(r=0.012),纹理特征变量组(r=0.068)共5 组变量相关性的平均绝对值,可以得到盐分指数变量组相较于其他特征变量组而言,与土壤盐度之间存在较高的相关性,而遥感影像数据指数及原始波段与土壤盐分含量相关性较低.(2)在37 个特征变量中,盐分指数S2、S1、S6 与土壤盐分含量关系较为密切,其中S2 的相关性系数r=0.201,在37 个特征变量中与盐分相关性最高,其次为S1 相关性系数 r=0.196,两者均为正相关;S6 的相关性系数r=0.192,与土壤盐分含量呈负相关.这3 个指数是适合作为ZY1-02D 影像在滨海区域进行盐渍化信息评估的盐分指数.相对而言,盐分指数中的SI-T、SI2则与土壤盐分含量相关性较低.(3)除盐分指数外,植被指数中的RVI(r=0.121),纹理特征中的MEA(r=-0.160)、COR(r=0.146)也与研究区土壤盐分含量存在较高的相关性.

2.3 不同模型以及不同方案结果与对比

将所提取的37 个特征变量作为全变量组分别构建AdaBooost、LightGBM、XGBoost、RFR、CatBoost 五种机器学习算法的土壤盐分含量估算模型,模型精度如图3 和表5所示.通过计算训练集和测试集的R2以及RMSE,可以看出在五种机器学习算法中,CatBoost 模型在测试阶段中精度最高,R2达到了 0.8317,RMSE 达到了 0.3957g/kg.LGBM、XGBR、RFR 在测试阶段中的模型拟合精度较为接近,三者的R2分别为0.5575,0.5768,0.5914,RMSE分别为0.6288,0.6243,0.6317g/kg,但在预测精度上与CatBoost 还有一定的差距.在5 个机器学习模型中AdaBoost 的预测精度最低,其中R2仅为0.2247,RMSE 为0.8244 g/kg,远远低于其他4 个机器学习模型,由此可见AdaBoost 不适合于本研究区土壤盐度含量反演.

表5 基于全变量组的不同模型精度Table 5 Different model accuracies based on full variable groups

图3 基于全变量组的不同模型精度Fig.3 Different model accuracies based on full variable groups

为探究不同特征变量在同一机器学习算法下对土壤盐分含量反演的准确性的影响,本研究设计了5 种不同变量组合场景,基于CatBoost 机器学习算法对5 个场景进行建模与训练.其中场景一为仅使用原始波段组,场景二为原始波段与植被指数13个特征变量,场景三为原始波段、植被指数与盐分指数26 个特征变量,场景四为原始波段、植被指数、盐分指数与主成分波段29 个特征变量,场景五为全变量组共37 个特征变量.不同场景测试结果如表6所示.

表6 基于CatBoost 模型的不同变量组精度Table 6 Accuracy of different variable groups based on CatBoost model

可以发现,拟合效果最好的变量组合为SC5(R2=0.8317,RMSE=0.3957g/kg),其次是SC4,该场景测试集R2为0.7348,RMSE 为0.6199g/kg,R2值较全变量组相差0.086,RMSE 值较全变量组相差0.2242g/kg.相比之下,仅使用原始波段作为特征变量的模型拟合程度最差,R2值为0.2016,RMSE 值为2.2587g/kg,并且出现一定程度的过拟合现象;其次使用原始波段与植被指数作为输入变量的模型拟合效果也不理想,R2值为0.3505,RMSE 值为1.2702g/kg,也存在一定程度的过拟合现象.综上所述,可以认为仅将原始波段作为输入变量以及原始波段、植被指数作为输入变量的反演模型不适于研究区土壤盐分含量反演.

2.4 特征变量重要性分析

CatBoost 算法可以识别特征变量的重要性,本研究将37 个特征变量放入模型进行重要性排序,检验特征变量对于土壤盐分含量反演的重要性.结果重要性值越高,则代表该特征变量在土壤盐分反演中越重要.重要性排序结果如图4所示,其中可以发现对于滨海地区土壤盐分反演重要性最大的特征是MEA,该特征表示的是遥感影像第一主成分纹理特征的均值,次重要的特征是植被指数ARVI,原始波段蓝波段(B)、植被指数COSRI、盐分指数SI1,近红外波段(NIR).在原始波段变量组中,蓝波段与近红外波段对模型贡献度较高,而红波段贡献相对较低;植被指数中,由近红外波段、红波段、蓝波段计算得到的ARVI,由绿波段、蓝波段、近红外波段计算得到的COSRI 以及由近红外波段、绿波段计算得到的GRVI 对模型有重要贡献,但SAVI在模型中贡献较低;盐分指数中,SI1 以及S3 对模型有重要贡献,但NDSI 以及S5 对模型精度影响不大;主成分波段中,第三主成分波段与第一主成分波段也对模型精度有重要影响;特征变量MEA 与COR对模型有较大贡献,但其中的DIS 以及ENT 对模型贡献度较低.

图4 基于CatBoost 模型全变量组的变量重要性Fig.4 Variable importance based on CatBoost model full variable group

2.5 研究区盐分反演制图

通过以上结果,可以发现利用CatBoost 机器学习算法,并且选用全变量组作为输入变量的模型拟合效果最好.因此本研究采用这个模型对研究区的土壤盐分含量进行反演,利用ArcGIS10.2 对反演结果进行可视化,得到的研究区土壤盐分含量反演图.结果如图5所示,研究区土壤盐分含量预测值在0～8.784 g/kg 范围内,均值为2.478g/kg.研究区域内土壤均呈现不同程度的盐渍化,盐渍化土以轻度等级为主,与实测样本的盐分含量结果相符合,表明反演结果准确率较高.通过统计各盐分等级土壤所占的像元,结果表明,轻度盐渍土分布广泛,面积占比最大为50.841%,主要集中分布在研究区西部,在东部地区分布较零散;其次为中度盐渍土,面积占比为34.792%,在研究区区域内分布较为零散;重度盐渍土面积占比为12.118%,分布主要集中在西南角以及中心区域;盐土面积占比为1.634%,主要集中在南边边界线附近;非盐渍化土最少面积仅占0.616%.

图5 研究区土壤盐分反演结果Fig.5 Inversion results of soil salinity in study area

3 讨论

3.1 不同机器学习算法的选择

近年来,越来越多的机器学习算法被运用在土壤盐分的研究中,也取得较好的结果.本研究中利用国产资源一号卫星数据,选用了5 个不同的机器学习算法对钦州湾滨海土壤盐分含量进行估算.通过比较5 个模型的结果,得到了CatBoost 模型反演土壤盐分精度最高,其中R2达到了0.8317,RMSE 为0.3957g/kg.该值高于2019年奚雪等[59]在黄河三角洲垦利县使用BP和SVM机器学习模型对当地滨海土壤进行盐分反演的精度(R2分别为 0.436 和0.387;RMSE 分别为1.297 和1.006g/kg);也高于2018年王多多等[60]以上海崇明东滩为研究区利用多元样条自回归模型(MARS)建立的滨海土壤盐分的估算模型(R2=0.74,RMSE=0.96 ds/m).由此可以发现本研究所采用的CatBoost 新型机器学习算法在滨海地区土壤盐分含量估算方面相较于传统机器学习模型以及线性回归模型有更好的拟合效果.

但本研究所采用的CatBoost 模型预测精度低于冯雪力等[61]利用BP 人工神经网络对内蒙古河套灌区进行盐分反演的模型精度(R2=0.8908),主要的原因可能在于所选择的数据源以及相关变量有关.冯雪力等[61]的研究选用的数据源是RADARSAT-2精细四极化单视复型(SLC)雷达影像,并且模型的输入变量是基于多元回归筛选的原始光谱反射率二阶导数、雷达后向散射特性以及地表组合粗糙度,这在一定程度上大大提高模型精度,该数据与本研究采用的国产资源一号卫星光学传感器数据相比,属于主动传感器成像方式,这种数据的雷达后向散射系数和土壤含水量之间具有较好的相关性,这对土壤盐渍化建模非常有用,因此冯雪力等[61]采用的雷达极化数据反演精度更高,但考虑到本研究主要探讨的是国产资源一号卫星在滨海土壤盐渍化反演中的适用性,因此本研究尚未加入雷达数据到训练模型中,未来可以考虑将SAR、PolSAR 及InSAR主动传感器卫星数据加入到模型中来提高模型训练的精度.但是,本研究所采用的CatBoost 模型精度显著高于张雅莉等[62]基于 landsat8OLI 影像光谱采用多元逐步回归的方法建立的土壤盐分预测模型(R2=0.7776),说明CatBoost 新型机器学习算法相较于传统的线性回归模型能够大大提升土壤盐渍化的模型拟合精度.

3.2 不同特征变量的选择

本研究中设置了不同的输入变量场景,发现当输入变量仅为4 个原始波段(SC1)以及由原始波段、植被指数(SC2)所构建模型拟合效果较差(R2分别为0.2016、0.3505).当引入盐分指数后,模型的预测值拟合能力得到较大的提升,预测精确度提高了34.53%.这可能与本研究中原始波段以及植被指数与盐分含量相关系数较低有关,并且盐分指数作为反映土壤盐分含量的直接参数,相较于其他间接参数对模型精度的提高起到较大的作用.这与边玲玲等[63]研究所得的结果一致,他们认为盐分指数与土壤含量之间存在一定的相关性,在估算应用方面具有很大贡献.Aldakheel 等[64]的研究也表明盐度指数对预测盐度有重要作用.在本研究中,当模型输入变量为全变量(SC5)时,得到了最好的土壤盐分反演模型精度(R2=0.8317),拟合效果最好.这可能是由于全变量组中所包含的信息较多,可以在一定程度上提高模型精度.全变量组中相较于SC3 引入了主成分波段以及由PCA1 计算得到的纹理特征,模型精度提高了12.45%.因此,可以证明在本研究中主成分波段以及纹理特征的引入可以有效提高机器学习模型的反演精度.此前,已有学者对主成分波段在土壤盐分含量反演过程中的表现进行探究,Metternicht 等[65]认为从理论上讲,主成分分析可以应用于盐分的识别和变化检测,因为PCA1 的稳定亮度和PCA2 的稳定绿度可以使得盐渍土壤从非盐渍土壤中分离出来,而PCA3 的差异亮度和PCA4 的差异绿度则解释了地表盐度发生的变化.Judkins 等[66]的研究也发现PCA1 与表层土壤盐分具有较高的相关性,对表层土壤电导率的变化具有较高的光谱响应性.

本研究对全变量组的特征变量重要性进行排序,排序结果如图4所示.可以发现纹理特征中的均值(MEA)对土壤盐分含量最为敏感,其次纹理特征中的相关性(COR)也对土壤盐分较为敏感.这可能是由于纹理信息是物体的固有属性,其反映的是图像灰度性质及其空间关系,能够在一定程度上增强模型的预测性能[67].对于纹理特征在土壤盐渍化方面的表现,许多学者都进行了相关研究,如Ren 等[68]的研究表明利用遥感光谱结合纹理特征可以提高滨海地区盐渍土分类精度,为准确掌握土壤盐渍化提供了新途径.台翔等[69]通过建立植被覆盖下的土壤含盐量极限学习机反演模型,发现加入纹理特征参数能有效提高模型精度.Hoa 等[70]通过构建越南湄公河三角洲地区土壤盐分反演模型,认为从图像波段提取的纹理特征对土壤盐度建模有重要作用.在本研究中纹理特征贡献了重要作用,这与此前的研究结论大抵相符.

除纹理特征外,在本研究中蓝波段(B1)与近红外波段(B4)对土壤盐分含量也较为重要.此前已有大量研究证明光谱指数对于土壤盐度估算有所贡献.在这些研究中,可见的红色、中红外和近红外波段被认为是土壤盐度估计最重要的预测变量[71],Fan等[72]指出350～549nm(含OLI 传感器深蓝、蓝光和绿光波段),6695～1001nm(含OLI 传感器近红外波段),都是比较重要的辨别土壤盐碱化的光谱信息,能较好的提取土壤盐渍化的信息.Metternicet 等[73]发现,受盐分影响的土壤具有较高的光谱反射率,尤其是在蓝色波段.本研究大体上与这些结论一致.

在前人的研究中,植被指数被认为是估算土壤盐度的关键变量,许多学者通过构建土壤盐分与植被指数的模型也很好地对研究区进行了反演.本研究的结果也表明,植被指数对于土壤盐分估算具有一定的重要性.在本研究中,重要性前十个特征变量中,有三个是植被指数,分别是ARVI、CORSI 以及GRVI.其原因可能是这些植被指数主要由蓝波段、近红外波段以及红波段计算得到.边慧芹等[74]的研究表明了植被指数ARVI 与土壤含盐量相关性较高,能在一定程度上反映土壤盐分含量,可用于土壤盐分估算模型的构建.Vermeulen等[15]的研究证明由蓝色、红色和近红外波段计算得到的COSRI 是土壤盐分反演中要考虑的重要变量.本研究与以上研究结论较为一致,ARVI、CORSI 在模型中重要性排名第二、第四位,在估算研究区土壤盐分方面起到重要作用.

3.3 模型的不确定以及未来研究方向

本研究利用Catboost 模型对研究区土壤盐分进行反演,通过模型预测研究区的土壤盐分范围为0～8.784g/kg 范围内,均值为:2.478g/kg,得到模型精度R2=0.8317.利用机器学习模型拟合时,我们也注意到拟合参数与样本数的关系,本研究中土壤样本点数量有97 个,其中建模样本有74 个,检验样本有20 个,实际投入模型训练的样本数74 个远远大于模型的拟合变量达34 个,但是从模型的拟合效果中确实可以发现,模型的训练精度较高,而验证的精度稍低,这具有一定的局限性.但是,本研究在验证阶段的精度(R2=0.8317)也高于王多多等[60]在上海崇明岛东滩滩涂的土壤盐分反演精度(R2=0.74)以及张雅莉等[62]在艾比湖流域盐渍化土壤盐分的反演精度(R2=0.7776),因此研究结果可以应用于钦州湾滨海土壤盐分反演的研究中,但是应该注意到该模型对于不同的时空尺度土壤盐分的反演还有待进一步探索.

机器学习算法相对于传统的回归分析方法具有一定的优势,但机器学习方法的选择、参数的确定均对预测模型的精度有一定的影响,这可能会使得模型的预测结果存在一定的不确定性.其次,本研究所使用的土壤盐分含量是由土壤电导率利用公式进行转换,也可能会影响模型的估算结果;最后,卫星自身传感器误差以及遥感图像在进行处理的过程中存在的误差也可能会对模型的拟合结果产生一定的影响.在未来的研究中可以进行改进:适当增加土壤样本数量,提供更加真实可靠的样本数据;采用Python 语言Scikit-learn 提供的GridSearch 算法对机器学习算法进行超参数寻优;引入Sentinel-1A、Radarsat-2、Landsat 和SRTM DEM 等数据,提取相关的环境变量如极化雷达参数、地表温度和地形因子等构建土壤盐分含量模型;将CatBoost 与粒子群算法相结合,采用优选变量的方式,减少变量间的冗余程度以提高土壤盐分估算的精度.

近年来,遥感技术作为大面积土壤盐分的反演提供了可能性.以往的研究多采用Sentinel 卫星[59]、Landsat[60]以及RADARSAT-2[61]等主被动卫星数据对不同地区的土壤盐渍化进行反演,但是这些研究多采用国外的卫星数据源,而针对国内的众多卫星数据探讨较少.国产卫星数据像高分系列中的众多卫星如高分1 号和高分2 号均提供了光学卫星数据、而高分3 号也提供了SAR 数据,这些国产主被动卫星数据源均可以提取本研究所用到的各种土壤特征参数和变量,融合国产多源卫星影像数据的不同优势,这可为国产卫星数据在土壤盐渍化中的反演提供了可能,未来均可成为中国土壤盐渍化反演的重要数据来源.

本研究所采用的资源一号02D 卫星数据能够衍生大量的盐分指数、植被指数以及纹理指数等变量,提取优势特征成为提高土壤盐度估算精度的关键问题.为此,本研究基于样地实测数据和资源一号02D 卫星数据衍生的各种特征变量,首次对比了AdaBooost、LightGBM、XGBoost、RFR 以及CatBoost 五种不同机器学习算法在滨海土壤盐分估算中的优劣特性,同时探讨了不同变量选择对模型精度的影响.本文所提供的技术方法虽然应用于我国的南亚热带地区,但是该方法论同样的适用于北方的干旱半干旱地区,这为国际同行利用机器学习算法和国产卫星数据研究土壤盐渍化问题提供了技术经验和方法参考,其研究结果也开创了利用机器学习算法和国产卫星数据研究我国土壤盐渍化问题的先河.

4 结论

4.1 研究区实测土壤盐分含量范围为 0.727～10.352g/kg,均值为1.739g/kg,实测样点中轻度盐渍化占比最大.

4.2 对比5 种机器学习算法及不同输入变量,CatBoost 相较于AdaBooost、LightGBM、XGBoost、RFR 有更好的预测性能,CatBoost 结合全变量组在预测阶段取得了最高精度,R2达到0.8317,RMSE达到0.3957g/kg,适用于研究区土壤盐分含量反演研究.

4.3 纹理特征中的均值对土壤盐分含量最为敏感;第一主成分纹理特征表现最好,其次为植被指数与原始波段.

4.4 研究区土壤盐分含量预测范围为0～8.784g/kg,均值为2.478g/kg,轻度盐渍土分布广泛,主要集中分布在研究区西部,在东部地区分布较零散;其次为中度盐渍土,在研究区分布较为零散;重度盐渍土分布主要集中在西南角以及中心区域;盐土主要集中在南边边界线附近.