基于Sentinel-2A的孙吴地区土壤有机质反演研究

2022-10-28陈超群戴慧敏冯雨林杨泽杨佳佳

物探与化探 2022年5期

陈超群，戴慧敏，冯雨林，杨泽，杨佳佳

(1.中国地质调查局沈阳地质调查中心，辽宁沈阳 110034；2.自然资源部黑土地演化与生态效应重点实验室，辽宁沈阳 110034；3.辽宁省黑土地演化与生态效应重点实验室，辽宁沈阳 110034)

0 引言

土壤有机质是土壤质量的重要参数，可为农作物提供各类养分，同时对元素表生地球化学特征有重要影响。黑土作为珍贵的土壤资源，其有机质含量是反映土壤质量的重要指标参数[1-2]。近年来，随着黑土逐渐退化，土壤中有机质成分明显减少，估算黑土有机质含量，扭转含量下降趋势，是黑土地保护的重要举措[3]。传统的土壤有机质监测主要是通过对监控区进行大量野外土壤样品采集和室内化学实验分析进行反演，这种方法周期较长，费时费力，精度受样品密度控制，难以满足现代农业快速发展的需求[4]。随着遥感技术日益成熟，通过有机质含量的光谱差异来测定土壤有机质的含量已成为一种有效手段。

土壤有机质含量遥感反演主要包含两个研究方向：光谱信息的处理与选择和反演模型的构建。常采用的光谱处理方法有倒数、对数、去包络线变换等，但选取的有机质特征波段因影像数据源不同而有所差异。屈冉等[5]选取Landsat TM 影像反演广西壮族自治区富川县的有机质含量，认为土壤有机质含量与Landsat TM 波段5和波段7 的DN值相关性最高。陈德宝等[6]借助Landsat 8遥感影像对农安县黑土区有机质进行建模反演，表明短波红外B6波段反射率所建模型拟合效果最好。陈思明等[7]对Landsat 7 土壤光谱进行线性波谱分离重建，认为重建光谱能显著增强与土壤有机质含量的相关性，提高土壤有机质反演精度。在以往研究中，线性回归和偏最小二乘回归模型(PLSR)常被用于土壤有机质含量反演。Dhawale 等[8]结合土壤样品有机质含量和相应的土壤反射率，选用PLSR建模，均方根误差不超过2.24%。马驰[9]对比Sentinel-2A遥感影像不同波段组合的多元回归模型，R2均大于0.7。目前针对有机质敏感波段的选择主要采用Pearson相关分析法，反演模型也多选择线性拟合。本次研究借助Sentinel-2A遥感影像，结合黑河市孙吴县实测土壤有机质含量，通过Peason 相关分析和随机森林(RF)选择不同特征波段作为模型输入量，采用PLSR和BP 神经网络建模，以期研究土壤表层有机质含量与遥感影像关系，并实现地面黑土区红旗林场的土壤有机质高精度快速反演。

1 研究区概况

研究区孙吴县地处黑龙江省黑河市中部，位于东经126°39′35″～128°1′6″，北纬48°59′00″～49°41′55″(图1)。东部紧靠逊克县，西边为嫩江县，南侧与五大连池市相挨，北方为黑河市爱辉区，总面积4 318.9 km2。孙吴县海拔110～755 m，属于低山丘陵区，地势总体呈西南高东北低趋势。地貌分界清晰，从西到东分别为低山沟谷区、丘陵河谷地区和沿江平原。土壤类型以暗棕壤和草甸黑土为主[10]。气候属于寒温带大陆性季风气候，年均气温-0.6 ℃，年均降雨约550 mm，冻结期较长，无霜期短[11]。本文选择孙吴县红旗林场地区进行遥感反演，红旗林场位于孙吴县西北方向，范围为东经126°41′25″～127°14′34″，北纬49°16′32″～49°30′58″，界内发育孙吴县最高山峰松木山。

图1 孙吴县遥感影像(a)及红旗林场位置(b)Fig.1 Remote sensing image of Sunwu County(a)and the location of Hongqi Forest Farm(b)

2 数据采集与处理

2.1 土壤采集与有机质测定

按照《土地质量地球化学评价规范》(DZ/T 0295—2016)采样要求，在孙吴县采集土壤时去除表面枯枝落叶等杂物，用刻槽法垂直采集地表至20 cm深土样，保证上下均匀采集，并去除动、植物残留体、砾石、肥料团块等。土壤有机质含量采用硫酸—重铬酸钾法测定。共计采集806个样品，其中564个土样作为建模集，242个样品为测试集，统计信息如表1所示。

表1 土壤样品中有机质含量统计信息Table 1 Statistical information of organic matter content in soil samples

2.2 遥感数据获取与处理

选取研究区内2018年11月7日裸土无雪时期的Sentinel-2A影像，云覆盖0%。影像的预处理包括几何校正、大气校正、图像镶嵌及图像剪裁等操作。所选影像为Level-1C 上层大气反射产品，已经过系统几何精校正处理，其精度在一个像元内，满足研究需求。借助SNAP软件中Sen2cor280工具箱实现大气校正，校正后丢失卷云波段B10。为提高土壤有机质与光谱反射率(R)相关性，对遥感影像进行倒数(1/R)、对数(lgR)、幂函数(Ra)、一阶微分(FDR)、二阶微分(SDR)及倒数对数一阶微分(FDLR)处理。

3 算法原理

3.1 Pearson相关分析

为获取土壤有机质光谱响应波段，研究中采用Pearson相关判断反射率与有机质含量之间的线性相关性。公式为：

(1)

3.2 随机森林

2001年Breiman提出随机森林(random forest，RF)算法，主要优势体现在处理多维数据集时无需降维，最大程度的保留数据集原始信息[12-13]。本次研究中以X.IncMSE指标为重要性选择依据，X.IncMSE值越大，表明该波段信息越能反应有机质含量。

3.3 偏最小二乘回归模型

偏最小二乘回归分析(partial least square regression，PLSR)常用于遥感光谱反演建模，优势在于建模过程中集中了主成分分析、典型相关分析和线性回归分析方法的优点[14]。其建模思路为：设自变量X=[x1，x2，…，xp]n×p，因变量为Y=[y]n×1，其中X为波谱曲线反射率及其变换形式，Y为土壤成分含量，n为采集样本数，p为特征波段数。从自变量信息中提取最大变异信息成分μ1，且与因变量呈最大相关性。在提取第一主成分μ1后，建立Y与μ1的回归方程，若模型精度满意，则算法停止；否则继续利用X和Y被成分解释后的残余信息进行多次迭代提取，直到回归方程达到满意精度。

为了检验反演模型的精度及稳定性，借助决定系数R2和均方根误差RMSE作为模型评价指标，公式如下：

(2)

(3)

3.4 BP神经网络

BP神经网络(误差反向传播)是人工神经网络的一种，由输入层—隐含层—输出层组成[16]。学习过程包含信号的正向传播与误差的反向传播，正向传播时，信号从输入层传入，经隐含层逐层处理后，最后达到输出层。若输出层的实际输出与期望输出不符，则转向误差的反向传播阶段[17]。BP神经网络具有较强的非线性处理能力和自适应特点，能够较好地拟合光谱反射率与土壤有机质含量的关系。通过经验公式(4)确定隐含层的节点数取值范围，结合训练结果的精度选择最佳隐含层的节点数：

(4)

式中：n为输入层节点数；m为输出层节点数；k为1～10之间的常数；N为隐藏层节点数。

4 黑土区土壤有机质反演

4.1 特征波段选择

4.1.1 相关性波段选择

在SPSS 26平台下计算土壤有机质含量与Sentinel-2A遥感影像反射率及其变换间的相关性。如图2所示，有机质含量与光谱反射率呈现负相关，但相关性不高，各波段不同数学变换的相关系数最高值绝大多数出现在FDLR变换，表明该预处理方法能有效提高Sentinel-2A反射率与土壤有机质的相关性。将不同数学变换中通过显著性检验的波段作为相应反演的特征波段，同时组合各波段中最高相关系数的变换形式作为一种响应波段参考，其中B8和B11无通过显著性检验的变换，故不予讨论(表2)。

图2 波段反射率及其变换与土壤有机质含量相关性Fig.2 Correlation between band reflectivity and transformations and soil organic matter content

表2 相关性选取特征波段Table 2 Feature bands selected by correlation analysis

4.1.2 RF重要性选择

采用R语言的randomForest包实现土壤有机质特征波段选取，其中默认生成500棵决策树，并进行5次重复十折交叉验证，结合最简原则选择不同光谱变换下的特征波段。以对数变换为例，图3交叉验证曲线展示了模型误差与用于拟合的自变量数量的关系，当波段数为6时，误差下降幅度基本保持不变，结合简约性原则，选择重要程度值从大到小排序前6的波段作为有机质反演建模的输入参数，实验中R、1/R、Ra、FDR、SDR以及FDLR均需要6个重要变量表示土壤有机质含量。为提高变量表达精度，将所有波段的变换作为RF的因变量，并需要X.IncMSE值前26个波段变换作为建模输入参数集，以精准表达有机质含量信息(图4)。RF重要波段选取结果如表3所示。

图3 对数变化交叉验证曲线Fig.3 Cross validation curve of lgR

图4 所有波段交叉验证曲线Fig.4 Cross validation curve of all bands

表3 RF重要波段Table 3 Important bands of RF

4.2 土壤有机质反演

4.2.1 PLSR模型反演

将Pearson相关分析(表2)和RF(表3)提取的特征波段作为自变量，土壤有机质含量作为因变量，建立有机质含量PLSR反演模型，如表4所示，结果显示相关-PLSR模型和RF-PLSR模型反演精度结果相近，FDLR变换和组合波段都能有效提高模型反演精度，其中针对传统的相关-PLSR模型，FDLR光谱变换的拟合程度最好，RF-PLSR模型中组合波段的效果更为显著。但PLSR模型下建模集和测试集的决定系数R2均未超过0.1。

4.2.2 BP神经网络模型反演

结合式(4)，根据相关性和RF重要程度选取获得的特征波段数，确定BP神经网络模型中隐含层层数。网络的训练函数为Trainlm，输入层和输出层传递函数分别为Tansig和Purelin，表5为相关分析和RF两种特征响应波段模拟结果，对比PLSR拟合结果(表4)，非线性拟合的多光谱遥感影像反射率与土壤有机质含量模型精度能得到显著提高。由于多光谱遥感的光谱分辨率较低，光谱包含的土壤信息较为复杂，因此无法类比高光谱土壤有机质遥感反演，线性回归拟合模型不能有效提取影像上土壤有机质含量信息[18-20]。相关性提取波段与RF提取的重要波段在进行BP神经网络建模时，建模集和测试集的R2主要集中在0.2～0.5，RMSE集中在1.3%～1.4%。相关-BP神经网络模型中FDLR建模拟合程度最高，建模集R2为0.623 7，RMSE为1.354 8%，测试集R2为0.444 6，RMSE为1.266 4%。RF-BP神经网络模型中组合波段建模拟合程度最高，建模集R2为0.724 5，RMSE为1.312 7%，测试集R2为0.541 8，RMSE为1.372 2%。

表4 基于PLSR模型的土壤有机质反演Table 4 Inversion of soil organic matter by PLSR

表5 基于BP神经网络模型的土壤有机质反演Table 5 Inversion of soil organic matter by BP neural network

结合表4中对比R、1/R等波段变换的不同提取方法，Pearson相关分析中选择的FDLR反射率变换在线性回归和非线性回归中都展现较高的拟合效果。因为Pearson相关性分析获取的是有机质含量与反射率间简单直线性相关的方向和密切程度，因此光谱处理的程度直接决定了与有机质含量的相关性，进而影响了模型反演精度。FDLR变换中光谱倒数对数计算可以有效放大相似光谱间的差别，再经过一阶微分处理后消除部分线性的背景，同时降低噪声光谱对目标光谱的影响程度。但相关性选择的组合波段包含了反射率变量处理的不同级别，当加入相关性低的波段变换，很有可能引入了土壤其他成分的特征信息，导致有机质反演精度降低。而在RF重要程度选择中存在误差验证，结合每次选择的特征集计算袋外误差率，最后选择袋外误差率最低的特征集作为回归模型的输入集。筛选出的波段可能与有机质含量相关性低，但叠加其他波段光谱特征反而提高了有机质估测精度。因此RF波段选取时，组合所有光谱反射率及变换信息后筛选出的特征响应波段更能充分反映有机质含量信息，建模精度显著提升。

4.3 红旗林场有机质空间分布

对比8种影像变换的相关-PLSR、相关-BP神经网络、RF-PLSR、RF-BP神经网络建模情况，选择模型拟合度最高、稳定性最好的RF-BP神经网络模型作为多光谱遥感数据的土壤有机质反演模型，并预测孙吴县红旗林场的有机质含量分布(图5)。红旗林场的有机质范围大致在0.1%～18.8%，主要集中在3.742 6%～12.455 2%，平均含量为7.939 9%，呈现出中间高、四周低的分布趋势。该地区土壤有机质地球化学分析结果显示，红旗林场有机质平均含量为8.51%，分布趋势为中西部含量偏高，向北侧逐渐降低，该结果与遥感反演程度几乎吻合。但由于在林场附近采集土壤样品较少，因此地球化学实测结果分布较为粗糙，而遥感反演获取的土壤有机质含量更为细致。