研究区位于重庆市南川区(东经106° 54'— 107° 27', 北纬28° 46'— 29° 30'), 幅员面积2 602 km², 境内属亚热带湿润季风气候, 气候温和, 雨量充沛, 既无严寒, 又无酷暑, 四季分明, 霜雪稀少, 无霜期, 热量丰富, 年均气温16.6 ℃。 南川区境内多山, 地形以山地为主, 地势呈东南向西北倾斜, 最高点金佛山风吹岭海拔2 251 m, 最低点骑龙鱼跳岩海拔340 m。 重庆市南川区主要土壤类型为紫色土、 黄壤和石灰(岩)土, 有机质含量低, 磷、 钾丰富, 石灰(岩)土质地粘重, 偏碱性, 土壤阳离子交换量和盐基饱和度均高。

基于重庆地理数据库随机选取694个样点, 用梅花五点采样法采集土样, 均匀混合后作为该点的土壤样品, 提取样点的土壤锰含量数据, 土壤锰含量区间为20~4 462 mg· kg^-1。 此外, 在几何校正后的Landsat影像上, 提取锰含量有关的光谱数据, 获取样点的植被指数。 选取其中的555个样点数据进行建模与测试, 用随机抽样法按3: 1的比例对样点划分了建模集[图1(a)绿点]和测试集[图1(b)黄点]。 为了验证基于植被光谱响应指标的土壤锰含量制图精度, 将剩余的139个样点作为制图精度验证点[图1(c)红点]。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 研究区及样本点分布Fig.1 Study area and sampling distribution

影像数据来源于地理空间数据云(https://www.gscloud.cn)的Landsat8 OLI影像。 Landsat-8卫星上携带两个传感器, 分别是OLI陆地成像仪(Operational Land Imager)和TIRS热红外传感器(Thermal Infrared Sensor), 卫星一共有11个波段, 波段1~7, 9~11的空间分辨率为30 m, 波段8为15 m分辨率的全色波段, 成像宽幅为185 km× 185 km。 本研究获取的数据成像时间为2019年8月13日, 云量为0.09%, 影像质量较高, 符合使用要求, 为消除传感器和大气影响, 利用ENVI 5.3软件对Landsat8影像进行辐射定标和大气校正等预处理。 此外, 土地覆盖数据来源于GlobeLand30(http://www.globallandcover.com/)。

利用Origin软件对694个土壤样点中土壤锰含量和植被样点进行描述性分析, 包括样点数据的取值范围、 平均值、 标准差和平均绝对偏差等能够反映数据基本情况的各项指标, 结果如表1。

表1

Table 1

表1(Table 1)

表1 土壤锰含量及植被样点描述性统计 Table 1 Descriptive statistics of soil Mn content and vegetation sampling sites 土壤属性 样点数 范围 平均值 标准差STD 平均绝对偏差 偏度 峰度 土壤锰/(mg· k g-1-1) 694 20~4 462 1 021.45 785.34 611.42 1.33 1.65 NDVI(gr) 694 -0.41~0.73 0.53 0.15 0.10 -2.11 6.71 RVI 694 1.01~16.88 8.94 3.36 2.68 0.05 -0.41 VARI 694 -0.19~0.45 0.28 0.13 0.09 -1.38 1.51

表1 土壤锰含量及植被样点描述性统计 Table 1 Descriptive statistics of soil Mn content and vegetation sampling sites

1.3.1 植被光谱指标的筛选

土壤中锰含量过少, 植被叶绿体结构将受到破坏, 叶绿素含量随之下降, 将会导致叶片灰绿色或出现斑点等^[4]。 锰含量过高将阻碍植物对铁、 钙、 钼的吸收, 经常出现缺钼症状。 叶片出现褐色斑点, 叶缘白化或变紫, 幼叶卷曲等。 大量研究表明植被光谱指标能够表征植被生长健康状况^[5]。 因此选用了11个植被光谱指标, 其计算公式分别如下:

增强植被指数(enhanced vegetation index, EVI)

2.5×(B5-B4)/(B5+6B4-7.5B2-1)(1)

比值植被指数(ratio vegetation index, RVI)^[6]

B5/B4(2)

绿度总和指数(sum green index, SG)

B5/B3(3)

三角植被指数(triangular vegetation index, TVI)^[7]

60×(B5-B3)-100×(B4-B3)(4)

可见光大气阻抗植被指数(visible atmospherically resistant index green, VARI)^[6]

(B3-B4)/(B3+B4-B2)(5)

归一化植被指数(normailized difference vegetation index, NDVI)^[8]

(B5-B4)/(B5+B4)(6)

绿度归一化植被指数(green NDVI, GNDVI)^[9]

(B5-B3)/(B5+B3)(7)

蓝度归一化植被指数(blue NDVI, BNDVI)

(B5-B2)/(B5+B2)(8)

蓝-绿通道归一化植被指数(green-blue NDVI)^[10]

(B5-B3-B2)/(B5+B3+B2)(9)

红-绿归一化植被指数(green-red NDVI)^[11]

(B5-B3-B4)/(B5+B3+B4)(10)

叶绿素指数(modified chlorophyll reflectance index Green, mCRIG)^[12]

[(B2)-1-(B3)-1]×B5(11)

式中B₂, B₃, B₄, B₅分别为蓝波段、 绿波段、 红波段和近红外波段反射率。

为了获取最优植被光谱响应指标, 利用皮尔逊相关系数对11个植被指数进行筛选, 皮尔逊相关系数表达式如式(12)

R=∑i=1n(Xi-X̅)(Yi-Y̅)∑i=1n(Xi-X̅)2∑i=1n(Yi-Y̅)2(12)

式(12)中: R为相关系数; X与Y分别表示2个用于比较的变量, X̅, Y̅分别表示2个变量的平均值。 |R|值越大, 变量之间的线性相关程度越高; |R|值越接近0, 变量之间的线性相关程度越低。 在此基础上, 利用方差膨胀因子(variance inflation factor, VIF)对相关系数筛选的结果进行进一步筛选, 消除指标之间的冗余, 确定反演土壤锰的最佳植被光谱响应指标。 VIF判断规则如下: 当0< VIF< 10时, 表明指标之间不存在多重共线性; 当10< VIF< 100时, 表明指标之间存在较强的多重共线性; 当VIF> 100时, 表明指标之间存在严重的多重共线性^[13]。

1.3.2 土壤锰含量最佳反演模型的筛选

为了获取土壤锰含量最佳反演模型, 基于最佳植被光谱响应指标, 利用偏最小二乘法(partial least-squares regression, PLSR)、 多元逐步回归(multiple stepwise regression, MSR)和BP神经网络(back propagation neural network, BPNN)构建植被光谱响应指标和土壤锰元素之间的关系模型, 比较分析三个模型的精度, 确定最佳反演模型。

(1)偏最小二乘法

偏最小二乘回归(PLSR)综合了典型相关分析、 多元线性回归和主成分分析三种方法, 可以用于解决很多普通多元回归无法解决的问题, 具有计算简便、 预测能力强、 模型稳健等特点^[14]。 该方法的最简形式是用线性模型来描述预测变量x和独立变量y的关系, 如式(13)所示。

Y=Xβ+ε(13)

式(13)中, Y为经过归一化的因变量(土壤锰含量); X为经过归一化处理的自变量(最佳植被光谱指标); β 为系数矩阵; ε 为残差矩阵。

(2)多元逐步回归

多元逐步回归(MSR)是回归分析中一种筛选变量的过程, 结合了前向选择和后向剔除的优点, 其本质是建立反映变量之间变化关系的最优多元回归模型。 采用逐步回归法来建立植被指数与土壤锰元素含量的统计回归模型^[15]。 回归方程如式(14)

Y=b0+b1x1+b2x2+…+bpxp+ε(14)

式(14)中, Y为经过归一化的因变量(土壤锰含量); X_p为经过归一化处理的自变量(最佳植被光谱指标); b₀为截距; b_i为回归系数; ε 为误差项。

(3)BP神经网络

BP神经网络是一种由输入层(input layer)、 隐含层(hide layer)和输出层(output) 构成的前馈式误差反向传播神经网络, 该网络的特点是神经信号向前传递, 误差沿着网络反向传播, 该网络包含了神经网络的精髓内容, 可以实现由任意的m维到n维的映射关系。 BP算法的中心思想是调整权值使网络总误差最小, 通过把学习的结果反馈到中间层次的隐含层单元, 改变它们的权系数矩阵, 从而达到预期的学习目的^[16]。 典型的BP神经网络结构如图2。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 典型BP神经网络结构图Fig.2 Typical BP neural network structure

1.3.3 精度验证

参考现有的土壤属性含量光谱估算的评价标准, 采用均方根误差(RMSE)、 相对分析误差(RPD)和决定系数(R²)3个评价指标评估模型反演效果^{[17, 18]}。 决定系数R²用于检验回归方程对样本预测值的拟合程度[式(15)], RMSE用于对模型的稳定性进行评价[式(16)], RPD能够反映回归模型的可靠性[式(17)、 式(18)]。

R2=SSRSST=∑i=1n(y^i-y̅)2∑i=1n(yi-y̅)2(15)

RMSE=∑i=1n(yi-y^i)2n(16)

STD=∑i=1n(yi-y̅)2n(17)

RPD=STDRMSE(18)

其中, n为样品数目, y_i为实测值, y^i为模型预测值, y̅为实测值的平均值, STD是实测值的标准差。 RMSE值越小, RPD值越大, R²值越接近1, 说明反演模型的精度越高; 反之, RMSE值越大, RPD值越小, R²值越小, 说明反演模型的精度越低。

利用Landsat8数据, 按照表1中的公式获取11种植被光谱响应指标(RVI、 GVI、 mCRIG、 GBNDVI、 GRNDVI、 GNDVI、 BNDVI、 NDVI、 VARI、 EVI、 TVI), 结果如图3。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 植被光谱指标Fig.3 Vegetation spectral indices

利用IBM SPSS Statistics26对南川区森林覆盖区域进行土壤锰和11种植被光谱响应指标的相关性分析, 结果如图4所示, 选取相关性大于0.5且相关性显著的植被光谱响应指标进行多重共线性诊断。 获得3个与土壤锰元素有较高相关性且VIF均小于10的植被光谱响应指标(如表2)。 因此, 本文选取RVI、 GRNDVI和VARI为最佳光谱响应指标。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 土壤锰与植被光谱指标相关性Fig.4 Correlation between soil manganese and vegetation spectral indexes

表2

Table 2

表2(Table 2)

表2 最佳植被光谱响应指标 Table 2 Best vegetation spectral response index 植被指数 RVI NDVI(gr) VARI Pearson相关性 0.753* * 0.591* * 0.531* * VIF 7.310 8.559 4.481 注: * * 在0.01 级别(双尾), 相关性显著 Note: * * Significance at 0.01 level (two-tailed), indicating significant correlation

表2 最佳植被光谱响应指标 Table 2 Best vegetation spectral response index

将最佳植被光谱响应指标作为自变量, 南川区土壤锰含量作为因变量, 利用PLSR、 MSR和BP神经网络构建两者的关系模型。 基于PLSR的锰含量反演模型为: y=282.69+419.07RVI-5263.27GRNDVI-752.27VARI(R²=0.66), 结果如图5(a)。 基于MSR的锰含量反演模型为: y=272.776+417.174RVI-5216.026GRNDVI-748.28VARI(R²=0.66), 结果如图5(b)。 基于BPNN的锰含量反演模型中设置了3个隐含层和10个神经元节点网络迭代次数设置为2 000, 学习率为0.01, 学习目标为0.01; 输入层设置共3个, 包括GRNDVI、 RVI和VARI; 输出层为土壤锰含量, 这里设为1; 模型的建立和运行在MATLAB R2019a上完成, 结果如图5(c)。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 建模集实测值与估测值散点图Fig.5 Scatter plots of measured and estimated values of modeling set

由图5可知, PLSR模型和SMR模型的建模结果较为相近, 在土壤锰含量0~1 500 mg· kg^-1区间具有较好的拟合能力, 预测点能够较为均匀分布在1∶ 1线, 而对于锰含量1 500以上的预测偏差值较大。 BPNN模型相对于前两种模型拟合程度更高, 预测效果最优, 在土壤锰含量0~4 000 mg· kg^-1区间都具有较好的拟合能力。

三种模型的测试精度如图6所示。 在测试过程中, PLSR模型和MSR模型在土壤锰含量0~1 500 mg· kg^-1区间具有较好的拟合能力, 预测点能够较为均匀分布在1∶ 1线, 而对于锰含量1 500以上的预测偏差值较大。 总体而言, 基于BPNN模型所建立的反演模型表现出最优的预测效果, BPNN模拟的重金属含量与实测值比较更接近于1∶ 1直线。 其中, BPNN模型建模的R²为0.78; 测试的R²为0.71; PLSR模型和MSR模型的拟合精度较为接近, 建模R²均小于0.7。 已有研究证明, BPNN在处理复杂的非线性建模问题中有较好的拟合能力, 而PLSR和MSR通常被用于确定自变量和因变量之间的线性关系。 本研究中, 非线性模型(BPNN)的验证精度明显高于其他两种线性模型(PLSR和MSR), 其原因可能是土壤锰含量和部分植被指数之间存在明显的非线性的关系。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 测试集实测值与估测值散点图Fig.6 Scatter plots of measured and estimated values of test set

由上述结果可知, BPNN模型明显优于其他模型, 可以较好地反映土壤锰与最佳植被指数之间的映射关系。 因此选用BPNN模型进行土壤锰空间制图, 其结果见图7。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 南川区土壤锰含量空间分布Fig.7 Spatial distribution of soil manganese content in Nanchuan District

土壤锰含量反演能够充分体现出区域每个像元位置的土壤锰含量的变化特征。 从土壤锰含量空间分布来看, 研究区的东北— 西南方向出现锰含量偏高的条带, 土壤锰含量相对偏高的区域, 是金佛山国家级自然保护区, 该地区可能存在较为原始的未开发锰矿区, 这也与样本点实测值的分布情况基本一致。

为了验证BPNN在区域尺度上估算土壤锰含量的可行性, 使用139个样点作为验证集, 对比分析其实测值与估算值的离散程度, 结果如图8所示, 可以看出, 锰含量的实测值与预测值的R²、 RMSE和RPD分别为0.69, 567.64, 1.30, 制图效果较佳, 可以基本反映土壤锰含量的分布特征。 说明利用BPNN模型进行土壤锰空间制图具有一定的可行性。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 验证集实测值与估测值散点图Fig.8 Scatter plots of measured and estimated values of validation set

利用遥感影像数据进行土壤锰含量的反演有助于高效监测区域尺度土壤锰含量。 当前已有研究利用土壤光谱数据直接对土壤锰含量进行反演, 并且反演精度较高(R²约为0.70~0.78), 可以满足土壤锰的监测需求。 但在植被常年覆盖的南方地区, 难以从卫星影像中获取土壤光谱, 此方法难以执行。

本研究引入植被光谱响应指标, 利用线性(PLSR和MSR)和非线性(BPNN)方法构建土壤锰含量反演模型, 与PLSR和MSR相比, BPNN能显著提高土壤锰含量的估算精度, 非线性回归模型在大面积反演土壤锰含量方面具有很大的潜力, 这表明土壤锰和植被指数之间的关系是较为复杂的, 不能用简单的线性关系来表达。 有学者研究土壤光谱和锰含量关联时, 发现在400~900 nm波段对土壤锰含量有明显响应, 会随着锰浓度增加出现依次降低的趋势, 特别是红波段响应最为明显, 而在绿波段其反射率变化不太大。 这与我们利用植被光谱进行土壤锰含量反演研究比较一致。

在利用植被光谱响应指标反演土壤锰含量时, 样点建模和空间制图均表现出不错的精度, 说明在区域尺度进行大面积反演土壤锰含量是可行的。 但受限于Landsat8数据具有光谱范围较宽、 波段数量较少等缺点, 反演精度偏低。 在未来研究中, 将尝试利用高光谱影像(如Hyperion影像等)反演土壤锰含量。 此外, 仅利用BPNN一种非线性方法反演土壤锰含量, 这可能会导致植被光谱响应指标与土壤锰含量响应关系未能被准确刻画, 未来引入更多非线性算法(随机森林回归算法、 决策树回归算法等)以期提高反演精度。