以湖南省郴州市东河流域铅锌矿矿区为研究区, 海拔335~451 m, 当地多种植茶树, 植被覆盖率高, 土壤呈酸性, 主要土壤类型为红壤, 黄壤和水稻土^[15]。 郴州开采了众多矿山, 包括煤矿、 铁矿、 铅锌矿、 石墨矿等。 其中, 清江铅锌矿即为郴州东河流域的重要矿山之一。 1985年郴州的东河流域发生溃坝事件, 郴州铅锌矿的尾矿库垮塌, 虽然后续做了大量的清理工作, 但是该流域土壤仍存在严重的重金属 (As、 Zn、 Pb等)污染。

1.1.1 土壤样本采集与Zn含量测定

2015年沿东河西岸在 113°17'18″E—113°17'36″E, 25°45'39″N—25°46'1″N 采集铅锌矿矿区土壤表层(0~20 cm)样本共73个, 取样点沿道路布设, 其中黄壤样本38个, 红壤样本35个, 采样点分布如图1所示。 采样过程中, 采用五点采样法在10 m²范围内采集表层20 cm深度范围内的土壤, 混合提取约500 g样本装入密封的土壤袋中, 将采集的土壤在实验室通风处晾干, 用有机玻璃压碎去杂, 然后用玛瑙研钵研磨, 最后过100目的筛子得到分选良好的土壤样本, 将土壤样本分为两部分, 一部分用于测量重金属Zn的浓度, 另一部分用于土壤样本的实验室光谱测量, 均存储在专用的容器中。 称取处理后的土壤样本0.1 g, 采用HNO₃-HF-HClO₄消煮法在微波消解仪进行消解, 用电感耦合等离子体质谱ICP-MS测定Zn元素的含量, 对采集的每个土壤样本平行测定三次, 取均值记为改样本的重金属含量。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 研究区及土壤采样点分布[15]Fig.1 Study area and distribution of soil sampling points[15]

73个土壤样本的重金属Zn含量统计值如表1所示, 其中黄壤样本38个, 红壤样本35个。 由表1可知研究区内总样本Zn含量的平均值为704.70 mg·kg^-1, 黄壤样本和红壤样本的Zn含量平均值分别为535.34、 888.57 mg·kg^-1, Zn含量远高于国家《土壤环境质量农用地土壤污染风险管控标准》(GB15618—2018)风险筛选值200 mg·kg^-1, 研究区内Zn含量严重超标。

表1

Table 1

表1(Table 1)

表1 土壤样本中Zn含量统计值 Table 1 Statistical value of Zn content in soil samples 样本类别 Min/ (mg·kg-1) Max/ (mg·kg-1) Mean/ (mg·kg-1) Median/ (mg·kg-1) Std/ (mg·kg-1) CV 总样本 60.44 4946.59 704.70 334.74 872.35 1.24 黄壤样本 60.44 2141.69 535.34 292.84 615.58 1.15 红壤样本 87.97 4946.59 888.57 364.60 1 064.08 1.12

表1 土壤样本中Zn含量统计值 Table 1 Statistical value of Zn content in soil samples

1.1.2 土壤反射率光谱测量

光谱测定在一个暗室内进行, 以一个50W的卤素灯作为光源, 利用PSR-3500便携式光谱仪(Spectral Evolution Inc.)和一个视场角为25°的光纤测量土壤样本的反射率光谱。 将光钎探头固定在土壤样本上方, 为减少实验台的背景反射光, 在实验台上铺一块黑布。 将过100目筛的土壤样本以培养皿的形状依次、 等距排列在黑布上。 测量过程中使用一个BaSO₄板作为定标板, 用于土壤反射率定标。 对每个土壤样品测量10条光谱, 并记录平均光谱作为土壤样本的最终光谱。

1.1.3 光谱预处理

光谱预处理主要去除噪声对光谱的影响。 首先, 由于光谱仪在波长区间两端的信噪较低, 导致土壤反射光谱在350~400和2 400~2 500 nm噪声非常严重, 故将这两个区间的光谱去除; 其次, 在余下的400~2 400 nm范围内仍存在不同程度的噪声, 分为三个噪声区间: 900~1 000和1 800~2 400 nm为严重噪声区间; 400~500和1 000~1 800 nm为中度噪声区间; 剩余的500~900 nm为轻度噪声区间。 采用分段Savitzky-Golay(SG)滤波进行噪声去除^[16], 分别采用窗口大小为21个、 15个和7个波段的二次多项式对严重噪声、 中度噪声和轻度噪声区间进行滤波处理。 经过去噪后的73个土壤样本的光谱反射率曲线图如图2所示, 黄壤和红壤样本的光谱反射率曲线图如图3所示。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 湖南郴州73个土壤样本光谱曲线Fig.2 Spectral reflectance curves of 73 soil samples in Chenzhou, Hunan

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 黄壤和红壤样本光谱反射率曲线图 (a): 黄壤样本光谱反射率曲线图; (b): 红壤样本光谱反射率曲线图Fig.3 Spectral reflectance curves of red soil and yellow soil samples (a): Spectral reflectance curves of yellow soil samples; (b): Spectral reflectance curves of red soil samples

为减小高光谱波段之间的相关性, 增强光谱的反射特征, 对光谱进行一阶微分处理, 一阶微分公式如式(1)

R'i=Ri-Ri-1, i=2, 3, …, n(1)

式(1)中, R'为光谱的一阶微分结果, R_i为第i个波段的反射特征, n为光谱波段总数。

技术流程图如图4所示。 首先, 获取土壤样本的实验室光谱并对其进行SG滤波和包络线去除; 其次, 根据土壤光谱活性物质对Zn的吸附作用, 从土壤光谱中提取有机质和黏土矿物的特征谱段; 最后, 分别对黄壤样本和红壤样本进行基于IGA-PLSR算法的Zn含量估算, 并对结果进行精度评价。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 研究方法技术流程图Fig.4 Flowchart of the proposed method

1.2.1 土壤光谱活性物质特征谱段提取

基于土壤反射光谱重金属含量估算机理研究表明, 有机质和黏土矿物对Zn的吸附作用是土壤反射光谱对Zn含量估算的主要机理^[17], 有机质对土壤反射光谱影响最大的范围在600~800 nm之间^[6], 对于风干后的土壤黏土矿物(主要指硅酸盐类黏土矿物, 包括高岭石、 蒙脱石、 蛭石等), 受羟基官能团和矿物富含的层间水分影响, 在1 400、 1 900和2 200 nm附近具有较强的光谱吸收特征^[18]。 因此, 土壤反射率光谱曲线600~800 nm附近的吸收峰可认为是土壤有机质的吸收特征, 以1 400、 1 900和2 200 nm为中心的吸收峰为土壤黏土矿物的吸收特征。 为获取有机质和黏土矿物的特征谱段, 采用包络线去除的方法增强土壤光谱的吸收特征。

1.2.2 训练集和验证集划分

训练样本和验证样本的划分对重金属反演建模的精度有很大影响^[19]。 选择Rank Kennard-Stone(Rank-KS)算法进行样本划分, 主要原因在于该算法同时考虑到光谱空间和性质空间(如重金属含量)对样本划分的影响, 使得训练集和验证集的划分更加合理^[19]。 Rank-KS首先将样本按照性质空间的大小进行排序, 将排序后的样本n等分, 然后在等分的n个小区间内分别利用Kennard-Stone (KS)算法选出在光谱空间具有代表性的样本作为训练集^[20], 从而提高训练集划分时在性质空间分布的均匀性。 n取反演精度最好时对应的值。

1.2.3 基于遗传算法改进的波段选择算法

在应用土壤光谱活性物质特征谱段进行反演建模前, 采用改进的遗传算法(improved genetic algorithm, IGA)^[21]进一步减少参与建模的波段的个数, 降低波段间相关性, 以提高模型估算精度。 IGA将模拟退火算法中的Metropolis准则和传统的GA算法相结合, 在生成新种群时, 当新个体目标函数值小于父个体目标函数值时, 用新个体取代父个体, 否则以一定的概率接受该新个体取代父个体, 从而解决GA算法经常会存在的过早收敛问题。 IGA能够在运行前期跳出局部解区域, 寻找更优解。 本研究采用IGA优选能够获取适应度函数最大值的波段组合, 将该组合作为自变量集合, 进行Zn含量的反演建模。

参考现有土壤重金属含量估算研究中对GA参数的设置^[12], IGA参数的设置如下: 染色体个数为20, 代际间隙为90%, 迭代次数为1 000次, 基因变异概率为10%, 目标函数为交叉验证的均方根误差。

1.2.4 偏最小二乘回归算法及其评价指标

野外土壤采样受限于采样、 制样与化验操作环节多、 成本高的限制, 往往仅限于少数几个野外站点, 很难获得大量的土壤样本^{[13, 22]}。 偏最小二乘回归算法(partial least squares regression- PLSR)兼具主成分分析、 典型相关分析和线性回归的优点, 在样本数量较少, 甚至样本点个数小于自变量个数、 自变量之间相关性比较严重的情况下依然可以比较好的进行建模分析^[23]。 基于此, 采用PLSR建立Zn含量的反演模型。

此外, 为降低土壤样本数量较少造成的过拟合的风险, 通过确定最大PLSR变量个数来限制PLSR的变量个数。 根据PLSR变量对因变量信息量的解释能力, 设置一个阈值来确定最大PLSR变量个数。 如果新的PLSR变量的加入对因变量信息量的解释能力不能超过这个阈值时, 就停止加入新的PLSR变量, 将当前的变量个数认为是最大的PLSR变量个数。 在最大的PLSR变量个数内, 采用交叉验证的最小均方根误差确定建模所用的PLSR变量个数。 结合PLSR变量个数对因变量信息量的解释能力, 将阈值设置为4%。 即, 当新加入的PLSR对土壤属性变量信息量的解释不足4%时, 停止新加入变量, 当前已经加入的变量个数为PLSR的最大变量个数。

对PLSR建立的土壤重金属反演模型, 采用预测均方根误差(RMSEP)、 相对偏差(RPD)、 决定系数(R²)进行评价。 三者的计算公式如式(2)

RMSEP=∑i=1n(yi-ŷi)2n(2)

RPD=stdRMSEP(3)

R2=∑i=1m(ŷi-y¯i)2∑i=1m(yi-ŷi)2(4)

其中, y_i和ŷi分别是第i个土壤样本的重金属含量分析值和经过PLSR建模之后的估算值, m是验证样本的个数, std是验证样本数据的标准差,y¯是验证样本数据的平均值。

模型的评价参考现有的土壤属性含量高光谱估算的评价标准: 出色模型, R²>0.9, RPD>3.0; 良好模型, 0.9>R²>0.82, 3.0>RPD>2.5; 近似模型, 0.82>R²>0.65, 2.5>RPD>2.0; 具有一定估算能力, 0.65>R²>0.50, 2.0>RPD>1.5; 不具备估算能力, 0.50>R², RPD<1.5^{[24, 25]}。

根据1.2.1节方法, 研究区土壤光谱经包络线去除后, 提取得到有机质相关特征谱段为638.2~727 nm, 共66个光谱波段, 黏土矿物相关特征谱段为1 343.3~1 473.7、 1 890~1 961.9及2 175.6~2 247.9 nm, 共96个光谱波段(如图5所示)。 由图6可知, 黄壤和红壤的有机质、 黏土矿物特征谱段光谱曲线形状相似, 但是红壤的有机质特征谱段反射率比黄壤高10%左右, 黏土矿物特征谱段反射率比黄壤高20%左右; 此外, 将黄壤的有机质和黏土矿物特征谱段曲线向上平移后可知, 在有机质特征谱段部分, 混壤光谱的吸收峰斜率^[26]大于黄壤, 在1 400和2 200 nm附近, 红壤光谱的吸收峰深度^[26]大于黄壤。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 黄壤和红壤有机质、 黏土矿物特征谱段Fig.5 Characteristic wavebands of organic matter and clay minerals in yellow soil and red soil

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 黄壤和红壤平均光谱的有机质、 黏土矿物特征谱段Fig.6 The characteristic wavebands of organic matter and clay minerals selected from the mean spectra of yellow soil and red soil

将上述提取的有机质和黏土矿物的特征谱段用于重金属Zn含量的反演建模。

不考虑土壤类别的情况下, 将全部73个土壤样本的光谱用于重金属Zn含量的估算。 采用Rank-KS算法划分为50个样本的训练集和23个样本的验证集, 统计特征如表2所示。

表2

Table 2

表2(Table 2)

表2 样本集的Zn含量数据统计 Table 2 Statistics of Zn content in sample set 样本集 样本 个数 Min/ (mg·kg-1) Max/ (mg·kg-1) Mean/ (mg·kg-1) Median/ (mg·kg-1) Std/ (mg·kg-1) CV 训练集 50 71.19 4 946.59 721.24 334.74 938.53 1.30 验证集 23 60.44 2 434.68 666.36 313.03 713.92 1.07

表2 样本集的Zn含量数据统计 Table 2 Statistics of Zn content in sample set

表3为全谱段和特征谱段估算精度的对比, 全谱段的估算精度R²和RPD分别为0.624和1.668, 而基于有机质煌黏土矿物特征谱段的估算精度R²和RPD分别为0.755和2.069, 根据2.3节的评价指标, 该模型达到近似模型的标准, 说明特征谱段的提取有利于提高土壤光谱速算Zn含量的能力。

表3

Table 3

表3(Table 3)

表3 不考虑土壤类型时重金属Zn含量估算精度 Table 3 Estimation accuracy of Zn content without considering soil type 光谱集 R2 RPD RMSEP/(mg·kg-1) 全谱段 0.624 1.668 385.670 特征谱段 0.755 2.069 345.079

表3 不考虑土壤类型时重金属Zn含量估算精度 Table 3 Estimation accuracy of Zn content without considering soil type

利用Rank-KS算法将38个黄壤样本划分为26个训练集和12个验证集。 针对黄壤样本, 进行基于有机质和黏土矿物特征谱段的Zn含量反演建模; 在保证与黄壤样本验证集相同的情况下, 针对全部土壤样本, 进行基于有机质和黏土矿物特征谱段和基于全谱段的Zn含量反演建模, 并与黄壤样本的建模结果作比较。

实验结果如表4所示, 在相同的12个黄壤样本的验证集上, 特征谱段的建模结果优于全谱段的建模结果, R²由0.732提升到0.761, RPD由2.018提升到2.137, 从RPD的角度来看, 特征谱段的建模结果满足了近似模型的条件; 黄壤样本的建模精度明显优于全部样本特征谱段的建模结果, R²由0.761提升到0.879, RPD由2.137提升到3.001, 从RPD的角度来看, 黄壤样本特征谱段的建模结果满足出色模型的条件, 且RMSEP的值在三个建模样本里最小。

表4

Table 4

表4(Table 4)

表4 12个黄壤样本的Zn含量验证精度对比 Table 4 Comparison of Zn content verification accuracy of 12 yellow soil samples 不同样本及 光谱集 R2 RPD RMSEP/ (mg·kg-1) 选择的波 段数量 有机质、黏土矿物特征谱段(黄壤) 0.879 3.001 184.725 73 有机质、黏土矿物特征谱段(全部样本) 0.761 2.137 259.462 92 全谱段(全部样本) 0.732 2.018 274.789 460

表4 12个黄壤样本的Zn含量验证精度对比 Table 4 Comparison of Zn content verification accuracy of 12 yellow soil samples

图7分别展示了三种样本集在运行5次取估算精度最高的验证集散点图, 可以看出黄壤样本特征谱段的估算结果明显优于另外两种估算结果, 其散点图更加集中于1∶1线附近。

图7

Fig.7

Figure Option ViewDownloadNew Window

图7 12个黄壤样本的Zn含量验证精度散点图 (a): 全部样本的全谱段估算结果; (b): 全部样本的特征谱段估算结果; (c): 黄壤样本的特征谱段估算结果Fig.7 Scatter plots of Zn content verification accuracy of 12 yellow soil samples (a): Estimation result obtained by using the entire spectral range of all samples; (b): Estimation result obtained by using the characteristic wavebands of all samples; (c): Estimation result obtined by using the characteristic wavebands of yellow soil samples

利用Rank-KS算法将35个红壤样本划分为23个训练集和12个验证集。 针对红壤样本, 进行基于有机质和黏土矿物特征谱段的Zn含量反演建模; 在保证与红壤样本验证集相同的情况下, 针对全部土壤样本, 进行基于有机质和黏土矿物特征谱段和基于全谱段的Zn含量反演建模, 并与红壤样本的建模结果作比较。

实验结果如表5所示, 在相同的12个红壤样本作为验证集的条件下, 特征谱段的建模结果优于全谱段的建模结果, R²由0.814提升到0.866, RPD由2.420提升到2.848; 红壤特征谱段的估算精度明显优于全部样本特征谱段的估算结果, R²由0.866提升到0.939, RPD由2.848提升到4.212, 从RPD和R²的角度来看, 红壤特征谱段的估算结果都满足出色模型的条件, 且RMSEP的值在三个建模样本里最小。

表5

Table 5

表5(Table 5)

表5 12个红壤样本的Zn含量验证精度对比 Table 5 Comparison of Zn content verification accuracy of twelve red soil samples 不同样本及光谱集 R2 RPD RMSEP/ (mg·kg-1) 选择的波 段数量 有机质、黏土矿物特征谱段(红壤样本) 0.939 4.212 186.623 73 有机质、黏土矿物特征谱段(全部样本) 0.866 2.848 275.981 66 全谱段(全部样本) 0.814 2.420 324.780 445

表5 12个红壤样本的Zn含量验证精度对比 Table 5 Comparison of Zn content verification accuracy of twelve red soil samples

图8分别展示了三种样本集及特征在运行5次估算精度最高的验证集散点图, 可以看出红壤样本的特征谱段估算精度明显优于另外两种的估算精度, 其散点图更加集中于1∶1线附近。

图8

Fig.8

Figure Option ViewDownloadNew Window

图8 12个红壤样本的Zn含量验证精度散点图 (a): 全部样本的全谱段估算结果; (b): 全部样本的特征谱段估算结果; (c): 红壤样本的特征谱段估算结果Fig.8 Scatter plot of Zn content verification accuracy of 12 red soil samples (a): Estimation result obtained by using the entire spectral range of all samples; (b): Estimation result obtained by using the characteristic wavebands of all samples; (c): Estimation result obtained by using the characteristic wavebands of red soil samples