空间关联随机森林模型结合Sentinel-2 影像估算潍北地区裸土期土壤盐分

2023-12-26彭远新王泽强周忠科宋晓宁徐夕博

中国环境监测 2023年5期

彭远新,王泽强,周忠科,宋晓宁,徐夕博

1.枣庄学院旅游与资源环境学院,山东枣庄 277160

2.山东师范大学地理与环境学院,山东济南 250358

土壤是地表生物的各项生命活动正常运行的基础,在生态系统的物质转化和能量传递中起到关键的媒介作用[1-2]。土壤盐分含量(SSC,本文均指质量分数)是评价土壤质量的一项重要指标,其对绿色植物生长异常敏感。土壤在发生盐化后,钠离子会发生富集进而限制植物细胞中的物质交换,导致植物生长受限甚至死亡[3-4]。在全球变暖和人类活动的作用下,全球范围内约有7%的土地受到不同程度的盐化影响,1/3 的世界耕地存在盐渍化的风险,土壤盐化对干旱区域乃至世界的生态环境安全带来极大的威胁[5-6]。所以,为实现土壤盐化的风险预警和治理恢复,快速精准地进行区域SSC 识别与分异状态监测尤为重要。

近年来,遥感作为一种新型对地观测技术,因其具有覆盖广、数据更新快和获取成本低等优势,被广泛用于土壤制图与环境监测领域。彭杰等[7-10]在实验室控制条件下利用土壤盐分光谱定量分析技术,分析出SSC 在450 ～680 nm 和短波红外的波谱范围内存在确定的响应波段,并利用其作为模型输入变量完成SSC 的光谱定量估算,证实SSC 光谱定量估算方式的可行性。在此基础上,LIU 等[11-14]分析了卫星影像波段反射率(Landsat-5、Landsat-8、Sentinel-2 和HJ-1 等)与土壤盐分间的相关关系,探寻光谱波段吸收特征与石膏、碳酸钙、硫酸钠等盐分矿物间微妙且稳定的联系性,最后建立遥感估算模型,分析得出区域尺度SSC 分布规律。

盐分在表层土壤中的积累是地貌、降水和人类活动等因素共同作用下的一种综合表现,因而在利用统计模型拟合SSC 与卫星接收的光谱发射率间的关系时,仅通过简单线性模型是很难取得较高精度的,而基于非线性回归技术的机器学习算法在构建SSC 遥感估算模型中具有更大的优势。例如,王飞等[15-17]在卫星遥感数据基础上采用机器学习技术(例如人工神经网络、支持向量机和随机森林等)完成了土壤盐分在区域尺度上的反演制图。通常来说,机器学习过程中通过增多隐藏层数量、增加结构参数设置和决策树节点数目,可以实现对土壤属性到光谱信号间的信息转换。但是,输入样本数据在土壤环境、人类活动和采样过程等因素影响下表现出的特异性和空间非关联性增加了目标地物与光谱特征间映射关系的复杂度,使机器学习模型的泛化能力减弱和出现局部极值问题,限制了估算精度的提升和模型的可应用性。针对此问题,方利民等[18-22]探索在遥感估算模型构建之前通过预先对样本数据进行随机选择、含量分布、光谱特征和土壤类型划分等方式减弱或消除特异性样本的影响,选取最佳训练样本集,在确保验证集精度的同时有效提升训练效率和准确度。上述研究对特异样本数据的识别均是基于样本的含量特征信息,然而,采集到的地理样本除带有描述统计的含量特征数据,还具有显著的地理位置特征信息[23-24]。根据地理学第一定律[25],地理事物之间距离越近,表现出的空间相关性越强,而特异性地理数据则呈现出弱的地理相关性并在空间上表现出非关联性。所以,在立足于决策树集成学习的随机森林中加入空间关联函数将空间位置信息引入,用以探寻样本数据在连续地理空间上的空间关联度,并完成对空间特异数据的识别和最优训练集的构建,减少样本量和提升计算效率,在区域尺度下的SSC状态识别与分布规律监测中具有一定的潜力。

潍北平原位于莱州湾南岸地区,土壤的母质及发育过程兼具海洋和陆地的特征,因此在此特定景观环境下构建准确高效的盐分估算制图模型,一直是资源环境领域的重点和难点工作。本文选取潍北平原为研究区,野外采集233 个土壤样本并获取同时相的Sentinel-2 多光谱影像,分别构建两波段差值、比值、归一值指数、相邻三波段正切指数和单波段共135 个光谱参量;进一步采用随机森林变量重要度评估技术,选择SSC 敏感光谱参量为输入自变量,测试得到的盐分含量值为因变量,分别建立基于随机森林和空间关联随机森林算法的遥感估算;最后完成区域尺度下的SSC 含量估算和数字制图,以期为实现土地资源优化管理和土壤环境政策制定提供理论支持和方法依据。

1 研究方法

1.1 研究区概况

研究区位于莱州湾南岸潍北平原地区(图1),地理坐标为北纬37°8'24″～37°16'36″和东经118°37'40″～118°46'24″,占地总面积约为59.14 km2。作为典型的温带大陆性季风气候区,该地受海洋环境影响明显且干湿季分明。据历史气象资料记录,多年平均气温和降水量分别为12.7 ℃和608.5 mm。研究区内整体地形平坦,以平原为主,适合农业种植,主要种植作物类型为玉米和小麦。地质构造上属辽冀台向斜第四系沉积层,当前地质状况较为稳定。因距海较近,地下水受到海咸水入侵,在冬季强烈的蒸发作用下土壤表现为轻微的盐碱化(1 g/kg ＜SSC ＜ 2 g/kg)[26]。

图1 研究区及采样点示意图Fig.1 The study area and samplings sites

1.2 样点数据与测试分析

在室内首先以ArcGIS 10.2 软件的数字底图为基础,综合考虑土地利用、水文地质和道路交通可达性等因素,完成233 处土壤样点的预设。采样过程中,将预设样点周边10 m 内设定为采样范围,采用多点混合法将土壤综合至1 kg 左右;其后装入聚乙烯密封袋中,送往实验室待测。同时,采用手持式GPS 确定并记录采样点的真实地理坐标,采样过程在2019 年1 月完成。在实验室内,首先去除土壤中小石块、木棒和草根等明显异质物;接下来,土壤样品在室温条件(25 ℃)下风干并过1 mm 筛,完成待测前处理;最后,运用质量法测定土壤样品的全盐含量,即吸取一定量基质水浸出液,并蒸干除去有机质,再烘干称量测得全盐量[27]。

1.3 影像获取及预处理

研究所使用的Sentinel-2(哨兵二号)多光谱影像数据(2019 年1 月17 日)由美国地质调查局网站提供,数据获取条件为晴朗无云天气及地面无明显积雪。下载得到的影像数据为L1C 级的大气表观反射率产品,在盐分估算中,要将L1C级数据转换为L2A 的地表真实反射率产品。需对数据进行辐射校正和大气校正,具体实现过程由欧空局推荐的SNAP 软件(已安装Sen2Cor 大气校正插件)完成[28]。哨兵二号影像数据中的波段1、波段9 和波段10 设计用于大气和水分特征的检测,不能用于土壤盐分含量的反演,因此在接下来的计算过程中将其排除。为确保各个光谱波段与地面信息间的匹配性,所有波段均进行重采样操作(空间分辨率为10 m)。并且影像与地面采样数据进行地理校正,确保像元偏差不大于0.5 个像元。此外,在土地覆盖以裸地为主,部分存在植被的区域,参照裸土的判别标准(NDVI ＜0.2)[29-30],剔除存在光谱干扰的植被区域。地面样品和影像数据获取时间正值冬季,降水较为稀少,因此土壤水分对光谱反射率的影响最小。

1.4 光谱指数构建

在土壤有机质、铁和黏土的重叠吸收作用下,土壤盐分在可见光-近红外范围内的光谱吸收特征通常具有非特异性,表现出覆盖范围广、光谱信号较弱的特点[31]。因此,除对多光谱影像的10个波段进行光谱分析以外,本研究还选择光谱差值指数(D)、光谱比值指数(RI)、光谱归一化指数(ND)和相邻三波段正切指数(tan)共4 种方式用于挖掘土壤盐分的特征光谱信号,计算过程:

式中:Bi和Bj分别表示哨兵数据的第i、j 波段的光谱反射率;tanpmq为相邻的3 个波段(波段p、波段m 和波段q)之间夹角的正切值,用于表征光谱曲线形状对土壤盐分的响应能力;kpm是波段p 和波段m 之间连线的斜率;kmq是波段m 和波段q之间连线的斜率。波段运算在ENVI 5.3.1 软件中的IDL 二次开发平台中完成。

1.5 空间关联随机森林模型

空间关联随机森林模型的构建是在随机森林算法[32]的基础之上(流程见图2),在样本数据输入侧引入空间权重函数[式(5)～式(7)],用以评估样本数据间的关联度和空间特异性,在完成输入样本数据的最优组合聚类后,进行回归模型的拟合计算。模型计算过程如图2 所示。

图2 空间关联随机森林模型流程图Fig.2 Work flowcharts for the spatial random forest model

式中:F 表示因变量yi、自变量xi和位置信息(μi,νi)间的函数关系; Fi(φ) 为对应位置i 处的空间权重高斯函数; δ 和L 分别代表随机误差和土壤盐分真实值;t 表示样本数据的空间关联度值,数值越小,样本数据空间关联性越强,依据此数据实现对输入数据集的组合聚类。

决策树数量(ntree)、总特征集上的特征子集(m try)、变量重要度值(VIS)和地理要素信息协同变量(Xi)是空间关联随机森林模型构建过程中的4 项重要参数。随机森林算法的核心思想是构建若干棵决策树(ntree),从决策树上的总特征集中,随机选择m 个特征子集(m try)用作预测变量,最后投票汇总结果;VIS 的大小用于衡量输入自变量对因变量的影响强度和解释能力的高低,在模型构建中将根据VIS 的大小选择最佳的光谱参数用于SSC 估算;此外,为增加空间权重评估函数对样本数据的位置信息解释能力和稳定性,选择典型且易于计算的5 个环境要素特征作为辅助信息,分别是地面高程、植被归一化指数、土壤湿度指数、地表温度和离海距离。地面高程数据免费下载自地理数据共享平台(http:/ /www.gscloud.cn/),NDVI、土壤湿度指数和地表温度数据的求算参照文献[33-34],获取时间与地面采样时间相同,样点离海岸线边界的距离的计算在ArcGIS 10.2 软件中的欧式距离计算模块中完成。在本研究中,ntree 和m try 分别设置为1 000 和3,模型的实现与计算在Python 3.8 中进行。

1.6 精度评估

采集到的土壤样本总集随机划分为建模集(n=209)和验证集(n=24),建模样本集用于估算模型的训练,估算模型精度的验证在验证集进行。决定系数(R2)和相对分析误差(RPD)是评价模型精度的重要统计指标。 R2表示自变量对因变量的信息解释能力的强弱,RPD 表征反演模型比只使用均值模型性能提升的能力。 RPD 大于2.0 时,说明模型取得满意的精度;RPD 处在1.4 ～2.0 之间时表明反演结果尚可接受,有待进一步提升;而RPD 小于1.40 时,意味着模型尚不具备反演能力[17]。通常,R2和RPD 越大,表明估算模型的精度和稳定性越强。

2 结果与讨论

2.1 土壤盐分含量描述性统计特征

潍北平原表层土壤(0 ～20 cm)SSC 的描述性统计特征如表1 所示。从表1 可以看出,土壤样品总集的SCC 为0.49 ～21.22 g/kg 内,均值为2.29 g/kg,覆盖范围较大且分布不均匀。 SSC 较高的采样点主要分布在中部和北部盐田附近,其余大部分采样点位于耕地,SSC 较低。在本研究中,变异系数主要用来评估数据的离散程度[35],3个样品数据集的变异系数均大于0.36,已处于高度变异,表明区域内SSC 分异较大。 SSC 数据存在空间变异会增加地表反演模型的训练难度,但一定程度的数据变异会有利于模型收敛并对模型的计算效率产生积极影响,所以随机抽取到的建模集和验证集用于土壤盐分的遥感建模反演是可靠的[36]。峰度和偏度用于描述数据分布的集中和双尾特性[37],3 个数据集的偏度值均大于2,表现出正偏,表明外部活动已经干扰了成土之初SSC 正态分布状态。土壤样本在总集合上的峰度偏高,说明盐分含量分布呈现出非集中的均匀分布态势。 SSC 的统计特征变异性,一定程度上也反映了土壤盐碱化的状态趋势。潍北平原地区工农业生产发展造成地下水过度开采,引发海咸水倒灌,加之旱季蒸发强烈,可溶性盐运移至表层土壤中产生结晶富集;土壤表现出轻微的盐碱化,并带来土地肥力下降等一系列的环境问题,这也是SSC 的数据分布表现显著变异的潜在驱动因素[38]。

表1 土壤SSC 描述性统计特征Table 1 Statistical characteristics of soil salt contents in the surface soils

2.2 选择光谱特征作为模型输入

通过式(1)至式(4)共计算得到125 个光谱指数,包括39 个D、39 个RI、39 个ND 和8 个tan。对SSC 最敏感的前10 个单波段和光谱指数的变量重要度(VIS)如图3 所示。 B11 是最重要的光谱波段,VIS 达到35.1%,其次是B3 和B8,其余波段的VIS 均小于10%。

图3 光谱参数变量重要度Fig.3 Variable importance scores of the spectral parameters

通过两两波段变换得到的光谱指数,对SSC信息的响应能力提升明显,波段比值指数的光谱信息提升能力最强;例如,B3 和B4 的变量重要度(VIS)分别为11.6%和4.4%,但经过波段比值变换后,RI34的VIS 达到27.5%,有效增强了土壤盐分的光谱特征。 B11 是重要度最高的单波段,在RF 变量重要度评估算法输出中的10 个最优光谱指数中,B11 参与了6 个光谱指数的构建,在盐分遥感反演模型和光谱指数的构建中应给予重点关注。基于此,选择3 个波段(B11、B3 和B8)和4个光谱指数(RI34、RI711、ND611和D45)作为SSC 的特征光谱参数,并在SSC 估算模型构建中充当输入自变量。

土壤盐化通常是盐分离子在蒸发作用下发生运移,在土壤表层生成白色晶体聚合物和盐壳的过程。盐分分子及其基团在晶格结构上振动的合频与倍频,会产生特异性的吸收特征峰[39]。WANG 等[3]和SIDIKE 等[40]在对SSC 的光谱吸收特征的室内控制实验分析时发现,波长在560 nm 和1 640 nm 附近波段与SSC 具有显著的相关性。 FOURATI 等[41]在利用卫星影像构建估算模型对土壤盐化水平进行识别时,得出了同样的结果,即处在短波红外范围内的2 个波段(波长在1 600 nm 和2 200 nm 附近)对SSC 表现出极显著的相关性。而本研究筛选得到的特征波段B3 和B11 的中心波长分别处在560 nm 和1 610 nm 附近,证实了本研究所选波段的典型性和估算模型构建的可靠性。此外,张俊华等[8]在分析宁夏银北地区盐化土壤的原位光谱分析时指出,600 nm附近是SSC 的一个明显吸收峰,并且640 ～870 nm光谱区间内波段的光谱反射率与碳酸根盐分离子显著相关,这与哨兵影像的波段8 的中心波长范围高度一致。总之,SSC 估算模型的输入特征波段有着明确的物理意义,能够反映地表土壤盐分的光谱反射率特征。

2.3 估算模型构建与精度对比

不同的特征光谱参数(特征波段和最优光谱指数) 构建的估算模型的预测结果如图4所示。

图4 不同估算模型的土壤盐分含量估算值与实测值对比散点图Fig.4 The scatter plots of estimated SSC values versus measured SSC values based on the different model

图4(a)和图4(b)分别是以3 个特征波段(B11 、B3 和B8)和4 个光谱指数(RI34、RI711、ND611和D45)作为输入自变量,测得的土壤SSC作为因变量,训练得到的基于RF 算法的反演模型的结果。可以看出,2 个估算模型的验证精度RPD 均小于1.40,实测值和估算值偏离1 ∶1 线较大,表明构建模型尚不能够进行土壤盐分的准确估算。图4(d)和图4(e)表示的估算模型是以空间关联随机森林算法为基础构建,对比RF 算法构型的模型[图4(a)和图4(b)],以空间关联随机森林算法为基础构建的估算模型[图4(d)和图4(e)] 的验证精度RPD 分别提升至1.17 和1.60,以4 个光谱指数(RI34、RI711、ND611和D45)为输入变量的空间关联随机森林模型,基本能够较准确地实现土壤盐分的反演[RPD ＞1.40,图4(e)]。

图4(f)是组合3 个特征波段和4 个光谱指数作为输入自变量,得到的空间关联随机森林模型估算效果,土壤盐分实测值与估算值的对比的散点基本在1 ∶1 线两侧,估算值与实测值偏差较小,反演结果最佳(R2=0.89 和RPD=2.04);而采用相同输入变量的RF 估算模型[图4(c)],则尚不能准确进行SSC 的遥感估算。总之,仅用特征波段或最佳光谱指数来构建SSC 估算模型均不能取得满意精度,空间关联随机森林模型的精度和稳定度要高于随机森林模型,再将特征波段和最佳光谱指数组合,输入到空间关联随机森林算法中去,建立得到的反演模型能够较准确地完成土壤盐分的遥感估算。

2.4 土壤盐分空间分布

以特征波段和光谱指数作为输入自变量的空间关联随机森林模型在不同空间关联度(t)下的精度变化如图5 所示。当t 为4 时,共有197 个数据样本(n 表示输入样本数)参与估算模型的建立,R2为0.77;t 继续减少,选取得到的数据样本间的空间关联度增强,t 为2 时,共有185 个土壤样本参与模型构建,取得了最高的估算精度(R2=0.89);随着样本空间关联度要求提升,样本量减少,模型逐渐不够完全训练,精度随之下降。

图5 不同t 时空间关联随机森林模型精度Fig.5 The perform ance changes of spatial random forest model with the different t value

为进一步验证构建的土壤盐分反演模型的稳定性和可靠度,将构建得到的最优反演模型推广至整个研究区,计算得到的盐分空间分布如图6所示。

图6 土壤盐分含量分布图Fig.6 The spatial distribution of the soil salt contents

从图6 可以看出,土壤盐分含量在中部地区略高于北部和南部地区,盐分含量高值区与盐田的分布位置密切相关。此外,研究区土地利用主要为农用耕地,地质地貌类型也较为单一,土壤盐分含量的变异并没有在大的空间尺度上表现得很明显。但是,SSC 在田块尺度上表现出显著变化,这与当地的微地貌和海水倒灌密切相关。研究区是国内知名的蔬菜生产基地,生产蔬菜的耕地需要进行大量抽水灌溉,地下水储备不足时引发海水倒灌,可溶性盐分运移至表层土壤产生一定程度富集,所以不同田块间的SSC 存在明显差异。

在遥感估算中,地理环境过程时刻在变化,土壤盐分的光谱特征仅强调瞬时信号是远远不够的。下一步要将遥感反演的经验统计模型与土壤盐分动态和环境的交互作用结合起来,从而最大化地减弱水汽、粗糙度和植被残积物对反演的影响,最终建立起稳健的遥感反演模型,实现动态大尺度层面上的土壤属性监测。其次,计算得到的光谱参数的重要度存在一定差异(即光谱参数对SSC 的响应能力不同),但在模型输入中并未考虑不同光谱参数的权重,若能将其对土壤盐分的不同光谱参数的响应强度考虑进估算模型,将有助于选取最佳分裂节点和决策树数目,提升模型的运算精度和效率。最后,混合像元问题是遥感估算中难以忽略的问题[12,42]。混合像元的存在,会减弱土壤盐分的目标光谱特征,增大信号噪声和模型误差。因此,在研究中选用空间分辨率(10 m)较高的Sentinel-2 影像作为数据输入源,土壤样点的采集也均在纯像元中进行,并且影像处理中掩膜掉道路、植被和建设用地等非裸土区域,以此最小化混合像元对反演结果的影像;并利用光谱指数构建法凸显SSC 的光谱信息,同时减弱环境弱信号。但是,实现像元解混合采用超高空间分辨率影像可能是解决混合像元问题的最佳方案[43-44],在下一步研究中需要参考改进。