基于偏最小二乘的土壤重金属镉间接反演模型

贺军亮，韩超山，韦 锐，周智勇，东启亮

(1.石家庄学院资源与环境科学学院，石家庄 050035；2.河北省水文工程地质勘查院，石家庄 050021)

0 引言

土壤是不可再生的重要资源，是人类赖以生存的物质基础，保护土壤资源，防治土壤环境污染，是推进生态文明建设和维护国家生态安全的重要内容[1-2]。随着工业的发展，过量重金属通过污水排放、大气降尘等途径进入土壤，造成土壤重金属污染[3]。土壤重金属不能被土壤微生物分解，具有残留时间长、难迁移、易累积的特点，而且可通过食物链以有害的浓度在人体内蓄积，严重危害人体健康和环境安全[4-5]。

土壤重金属含量的调查监测是进行土壤污染有效防治的前提。高光谱遥感具有快速、无损、光谱分辨率高等特点，可以在土壤定量遥感监测中发挥重要作用[6-7]。Kemper等[8]利用多元线性逐步回归等方法基于土壤反射光谱与重金属的相关分析构建模型，反演了Aznalcollar矿区土壤Pb，As，Fe和Hg的含量；郭云开等[9]基于水稻冠层光谱变化特征，利用最小二乘方法拟合建立了Zn，Pb和Cd等土壤重金属含量的反演模型；程先锋等[10]在兰坪铅锌矿区利用逐步回归方法构建了土壤重金属含量的高光谱估算模型。

一般情况下，土壤重金属含量较低，属于痕量级，即使在重污染状态下，重金属元素的直接光谱响应也非常微弱[11]。但是，当土壤有机质含量高时，对重金属的富集具有一定的吸附作用，使得利用有机质光谱特征间接反演重金属含量成为可能[12]。贺军亮等[13]利用水稻土有机质的敏感波段构建了有机质光谱诊断指数，用于重金属Cu和Pb含量间接反演模型的构建；兰泽英等[14]在乐安河流域采用土壤高光谱数据间接反演了Cu，Zn和Pb的含量，为该区域土壤生态环境监测提供了相关参考。除原始光谱反射率外，微分、倒数和对数等多种光谱变换指标能够从不同角度和程度上突出反演对象的光谱特征，降低背景噪声的影响[13-15]。前人所构建的重金属最优反演模型大多采用其中一种指标作为自变量建模，多种光谱变换指标的集成建模有可能会提高模型的估算精度和稳定度。

本文以石家庄市地表水源保护区土壤重金属Cd为研究对象，基于间接反演的研究机理，尝试通过多种光谱变换和相关性分析方法提取有机质光谱诊断特征，并建立多指标集成估算模型来间接反演重金属Cd的含量，进一步丰富土壤重金属高光谱反演研究的案例，以期为土壤定量遥感研究提供科学参考。

1 研究区及数据源

1.1 研究区概况

石家庄市地表水源保护区位于石家庄市西部的平山县和井陉县境内，保护区范围为N38°01′～38°45′，E113°34′～114°18′。一级保护区包括滹沱河干流、岗南水库和黄壁庄水库。一级保护区之外按缓冲距离和行政区划设二级和三级保护区。整个保护区地势西高东低，高差悬殊。虽然该保护区是石家庄市生活用水的主要供应地，也是北京市备用水源地，但是该区矿产开发历史悠久，人地矛盾突出，铁矿和石灰矿等矿产资源开采存在尾矿威胁、植被破坏、水土流失和土壤污染等状况，对保护区生态安全造成了一定影响。考虑到土壤空间分布的不均匀性及研究区地形的影响，参考研究区内主要采矿点和主要冶炼企业的区位，在研究区共采集69个土壤样品，土壤类型以褐土为主，采集深度为0～20 cm，同时利用GPS对采样点进行了定位(图1)。

图1 研究区采样点分布Fig.1 Sampling point distribution in the study area

所有样点采集的土壤样品带回室内后经风干、研磨、过筛(100目)后备用。土壤有机质含量采用重铬酸钾法测定，重金属Cd全量利用电感耦合等离子体质谱仪ICP-MS检测方法测定。

1.2 土壤光谱采集

土壤反射光谱测量在南京师范大学地理科学学院实验室完成。土壤反射率采用美国ASD公司生产的FieldSpec Pro便携式地物光谱仪测定，光谱采集范围为350～2 500 nm。光谱采集工作在室内进行，把卤素灯作为唯一光源。将土壤样品置于放在黑色绒布上的直径10 cm、深2 cm玻璃培养皿内，用工具将表面刮平。光源入射角为45°，距离土壤样品30 cm，采集枪垂直于土壤样品并保持15 cm的距离。每次测定之前都需要进行白板校正，每个样品测10次光谱曲线，取其平均值[16]。受检测周期长、土样残留水分以及光谱测量背景环境等影响，全光谱范围两端噪声较大，剔除光谱两端和水汽吸收影响较大的波段，最终保留1 640个有效波段进行后续建模研究。

1.3 光谱数据预处理

光谱变换可以很大程度上消除土壤背景的影响，进一步提高光谱信噪比，突出光谱的吸收和反射特征[17]。本文光谱实测值为土壤光谱反射率，采用以下方法对原始光谱反射率进行变换：一阶微分(first derivative，FD)、二阶微分(second derivative，SD)、倒数变换(reciprocal transformation，RT)、倒数的一阶微分(reciprocal transformation and first derivative，RTFD)、倒数的二阶微分(reciprocal transformation and second derivative，RTSD)、倒数的对数变换(absorbance transformation，AT)、倒数对数的一阶微分(absorbance transformation and first derivative，ATFD)、倒数对数的二阶微分(absorbance transformation and second derivative，ATSD)、连续统去除(continuum removal，CR)。将有机质含量实测值与以上各光谱数据进行相关性分析，在每项变换通过0.01水平显著性检验的波段中，找出相关系数绝对值最大处所对应的波段，作为该变换下有机质的敏感波段。光谱变换工作在Origin2017和ENVI5.3软件中进行，相关性分析在SPSS24软件中进行。

2 研究方法

2.1 建模方法选择

FD，SD，RT和AT等多种光谱变换指标能够从不同角度和程度上突出反演对象的光谱特征，考虑间接反演误差累积的可能影响，有必要对以上各光谱变换指标进行筛选。筛选方法采用多元线性逐步回归方法(multiple linear stepwise regression，MLSR)，保留对模型精度影响较大的变换指标。

反演模型的构建采用偏最小二乘回归方法。偏最小二乘法是一种数学优化技术，通过最小化误差的平方和找到一组数据的最佳函数匹配[18]。前人所构建的重金属最优反演模型大多只采用一种光谱指标作为自变量建模，即单光谱变换指标偏最小二乘(univariate partial least squares regression，U-PLSR)模型。多种光谱变换数据集成建模能够综合反映各种光谱变换形式的特点，通过尝试建立多光谱变换指标偏最小二乘(multivariate partial least squares regression，M-PLSR)模型，并与U-PLSR模型进行对比分析。

2.2 模型精度评价

Rank-KS方法可有效提升建模样本与验证样本选择的合理性[19]。采用Rank-KS方法将69个样本分为2组，其中46个样本用来建模，23个样本用来验证模型精度。在模型精度分析过程中，主要参考拟合优度调节参数R2，并结合均方根误差(root mean square error，RMSE)和相对分析误差(relative percent deviation，RPD)等参数进行模型对比。R2和RPD越高，RMSE越小，说明模型的拟合程度越好。RPD为样本标准差和RMSE的比值，可以用来判断模型的预测能力。当RPD<1.4时，认为该模型不具有预测能力；当1.4≤RPD<2时，可以对样本进行粗略预测；当RPD≥2时，模型具有极好的预测能力[20]。

3 结果分析

3.1 Cd特征量统计与污染现状分析

表1为Cd的特征量统计值。所采集土壤样本中，Cd的含量分布范围为0.138～0.359 mg/kg，平均值达0.220 mg/kg。总体样本、建模样本和验证样本的各统计量值都较为接近，标准差较小，说明Rank-KS选取样本较为合理，分类数据具有一定的代表性。

表1 Cd特征量统计Tab.1 Statistical characteristics of heavy metal Cd in soil

采用单因子污染指数评价法[21]对研究区Cd的污染现状进行评价分析，其计算公式为

Pi=Ci/Si，

(1)

式中：Pi为土壤中污染物的环境质量指数；Ci为实测值；Si为背景值。

石家庄地区Cd的背景值为0.08 mg/kg，当污染指数P≥1时，可以认定该地区存在Cd污染，指数越大污染状况越严重[22]。根据统计可以得出所有样本的污染指数P>1，92.8%的样本的污染指数P>2。因此，土壤样本采集区域存在较为严重的Cd污染，需要加强对该地区Cd含量的全面监测。

表2 土壤重金属Cd的污染指数统计Tab.2 Statistical analysis of pollution index of heavy metal Cd in soil

3.2 MLSR模型分析

经计算，总体样本的重金属Cd和有机质含量的相关系数达到0.7，并通过了0.01水平的显著性检验，说明两者之间存在一定的吸附赋存关系[12,14]。根据吸附机理，将有机质含量与各光谱变换指标进行相关性分析，分析结果如表3所示。

表3 土壤有机质含量与光谱变量的最大相关系数Tab.3 Maximum correlation coefficients of soil organic matter content and spectral variables

①R为光谱反射率；②**表示通过0.01水平的显著性检验。

从表3可以看出，原始光谱反射率在797 nm波段处与有机质含量存在最大相关性，与前人研究结果一致[23]，ATFD与有机质含量的相关性最大，且呈负相关，FD与有机质含量存在最大的正相关关系。根据有机质含量与光谱变量的相关分析，将各敏感波段对应的光谱变换指标作为Cd估算模型的自变量因子，Cd含量实测值作为因变量，构建MLSR模型。经MLSR分析，保留对模型精度影响较大的变换指标(ATSD1409和FD1396)，模型散点图及验证结果如图2和表4所示，具体模型为

Y=0.432+1 389.565XATSD-66.253XFD。

(2)

(a)建模样本 (b)验证样本

图2 MLSR模型散点图

Fig.2ScatterplotsofthetwosetsofsamplesoffittingandtestingforMLSR

表4 MLSR模型结果Tab.4 Results of MLSR model

建模样本MLSR模型R2为0.81，说明该模型对数据具有良好的解释能力。其RMSE仅为0.02，RPD为2.31，说明模型误差较小并具有良好的预测能力。而验证样本在R2和RMSE与建模样本相近的情况下，RPD却只有1.23，模型预测能力较差。通过图2可以看出，建模样本MLSR模型的趋势线与1∶1线夹角较小，模型预测值与实测值较为接近；验证样本的趋势线与1∶1线夹角较大，说明模型存在一定误差，模型稳定度有待进一步提高。

3.3 U-PLSR模型分析

经过3.2节MLSR分析，筛选出对模型精度影响较大的变换指标(ATSD和FD)。分别采用ATSD和FD光谱指标作为自变量，构建U-PLSR模型。ATSD-U-PLSR模型散点图及验证结果如图3和表5所示，具体模型为

Y=0.445 2+1 271.535XATSD。

(3)

(a)建模样本 (b)验证样本

图3 ATSD-U-PLSR模型散点图

Fig.3ScatterplotsofthetwosetsofsamplesoffittingandtestingforATSD-U-PLSR

表5 ATSD-U-PLSR模型结果Tab.5 Results of ATSD-U-PLSR model

与MLSR模型结果对比，ATSD-U-PLSR建模样本的R2降低了0.03，RMSE提高了0.017，重金属Cd实测值和预测值的拟合度下降，但RPD提高了12.09；验证样本R2，RMSE和RPD分别提高了0.05，0.024和7.87，相比MLSR模型数据预测能力有大幅提升。

FD-U-PLSR模型散点图及验证结果如图4和表6所示，具体模型为

Y=0.287 9+151.316 6XFD。

(4)

(a)建模样本 (b)验证样本

图4 FD-U-PLSR模型散点图

Fig.4ScatterplotsofthetwosetsofsamplesoffittingandtestingforFD-U-PLSR

表6 FD-U-PLSR模型结果Tab.6 Results of FD-U-PLSR model

从图4和表6可以看出，FD-U-PLSR模型建模样本与验证样本的拟合度均过低，R2仅分别为0.24和0.42。虽然总体样本FD光谱变量与有机质含量存在最大的正相关关系，但是相对于集成多光谱变换指标构建的MLSR模型，FD-U-PLSR模型估算误差较大，反映出单光谱变换指标估算模型的不稳定性。

3.4 M-PLSR模型分析

在MLSR和U-PLSR模型对比分析的基础上，考虑不同光谱变换形式提高光谱信噪比的能力差异，分别基于建模样本和验证样本数据集，进一步建立M-PLSR模型。所建模型解释变量与MLSR模型一致，仍为ATSD与FD。M-PLSR模型散点图及验证结果如图5和表7所示，具体模型为

Y=0.452 1+1 379.77XATSD-34.186 1XFD。

(5)

(a)建模样本 (b)验证样本

图5 M-PLSR模型散点图

Fig.5ScatterplotsofthetwosetsofsamplesoffittingandtestingforM-PLSR

表7 M-PLSR模型结果Tab.7 Results of M-PLSR model

M-PLSR模型建模样本与验证样本R2都超过了0.8，RPD都超过了2，RMSE也较低，并且验证样本与建模样本的预测精度相近。与MLSR模型和U-PLSR模型相比，M-PLSR模型拟合优度有所提高，模型稳定度得到明显增强。

4 结论与讨论

以石家庄市水源保护区褐土为研究对象，采用偏最小二乘法构建土壤重金属Cd的高光谱间接反演模型，具体结论如下：

1)石家庄市水源保护区内土壤样本采集区重金属Cd含量的平均值为0.220 mg/kg，高于背景值，并且总体样本的污染指数均大于1，存在较为严重的Cd污染。

2)研究区土壤重金属Cd和有机质含量在0.01水平上显著相关，相关系数达到0.7，两者之间存在一定的吸附赋存关系，基于土壤有机质的光谱特征间接反演重金属Cd含量具有一定的合理性。

3)有机质的原始光谱反射率对应的敏感波段为797 nm，各种光谱变换中倒数对数的一阶微分与有机质含量的相关性最大，可达到-0.766，一阶微分与有机质含量存在最大的正相关关系，达到0.721。

4)采用多元线性逐步回归方法筛选确定的模型解释变量为ATSD1409和FD1396。对于所建多光谱变换指标偏最小二乘模型，其建模样本和验证样本的R2分别为0.83和0.80，模型预测值与重金属Cd含量实测值较为接近，整体拟合效果优于单光谱变换指标模型。

与以往基于单光谱变换指标的建模方法研究不同，采用多种光谱变换指标进行集成建模，具有更好的建模效果和预测能力，可以对重金属Cd含量起到较好的预测作用。然而，由于影响土壤重金属反射光谱特征的因素不是单一的，仅仅考虑重金属与有机质的光谱响应特征还不全面。因此，在进一步研究工作中，可以尝试考虑其他因素如铁锰氧化物的综合影响，提取多因素特征波段构建多光谱变换指标估算模型。

此外，由于气候、地质地貌等条件的影响，不同地区土壤具有其独特的区域性。选取石家庄市水源保护区褐土为研究样本，土样经过研磨、风干处理后，虽然基本消除了土壤质地、土壤湿度等对土壤光谱的影响，但是所建立的土壤重金属Cd含量的估算模型是否适用于其他地区，有待进一步研究。