北京市平谷区桃林土壤重金属元素含量反演估测

刘泓君，牛 腾，于 强*, 苏 凯，杨林哲，刘 维，王慧媛

1. 北京林业大学林学院，北京 100083 2. 广西大学林学院，广西 南宁 530005

引 言

近些年来，矿产资源的开采及工业污水的排放对矿区周边自然环境的可持续发展和人体健康造成了严重威胁。特别在种植经济作物的地区，当重金属含量严重超标后进入种植物体内，进而造成植物光合作用减弱、 叶片呈现棕色斑块等症状。人类直接或间接食用此类作物果实后，易引起重金属中毒，重则危害生命[1-2]。北京东北部平谷区，栽培桃树面积世界第一，但其矿石的冶炼及采选造成重金属物质大量残留在土壤中，对周围桃林造成了严重的污染[3]，因此了解平谷区重金属污染情况极为重要。传统以“点采样+实验室分析”、 电感耦合等离子体发射光谱法等方法分析重金属含量，但此类方法依赖于大量的外业数据，费时费力，且不具有推广性。多光谱遥感具有高分辨率、 高精度、 快速无损、 大面积监测的优势，为宏观获取植被覆盖区土壤重金属元素污染情况提供了新的途径[4]。

目前，基于遥感技术监测土壤重金属污染取得了丰硕的成果。Anne等基于高光谱影像，提取红树林植被反射光谱，构建差值、 比值植被指数，基于偏最小二乘模型，成功反演出土壤中的有机质、 黏土矿物等含量[5]。高伟等以玉米为研究对象，通过运用包络线去除处理叶片光谱，构建作物铜铅污染判别规则[6]；杨灵玉等基于高光谱影像植被光谱反射率间接估算出土壤Zn和Cd元素含量[7]。综上所述，前人研究多集中在运用几个敏感的光谱波段或者常用植被指数建立反演模型，而对于基于敏感波段构建植被指数来进行反演的研究方法较少，且以桃叶为研究对象的更少。

以北京平谷区为例，将桃林作为监测目标，以野外采集的土壤样本、 桃叶光谱数据、 Sential-2多光谱影像为基础，探究受重金属胁迫影响较为敏感的特征波段并构建植被指数，同时将其与红色归一化植被指数、 类胡萝卜素反射指数2、 结构不敏感色素指数、 归一化水指数四种常用植被指数做对比分析，验证该指数在监测重金属方面的有效性与优越性；通过植被指数与土壤重金属Cd，AS和Pb的拟合建立反演模型，最终获取果园重金属污染空间分布情况，为快速有效地监测经济作物区域土壤生态状况提供技术支持。

1 实验部分

1.1 研究区概况

北京市平谷区(116°55′—117°24′E，40°1′—40°22′N)，西距北京市区70 km，总面积1 075 km2，约有40万人。属于温带大陆性季风气候，年降水量为629.4 mm，降雨主要集中在6月—9月，仅28%的降雨量出现在其他月份, 平均最高气温为17.3 ℃。平谷是中国大桃第一区，大桃种植面积达16.8万亩，产量1.5亿公斤，其品种达到130多个。平谷区拥有大量丰富的金属矿和非金属矿等矿产资源，矿床、 矿点、 矿化点共79处。该地区矿区开采过程中，矿石矿渣中的重金属经雨水冲刷、 淋溶渗入周边土壤中或河流中，对周边生态环境造成污染。

图1 研究区概况及采样点分布位置Fig.1 Overview of the study area and distributionlocation of sampling points

1.2 样本采集

依据当地的地理环境特征，并遵循均匀性、 代表性、 规律性的原则来布置采样点，共87个采样点的空间分布如图1所示。选择无风晴朗天气，于2021年7月4日至25日10：00—14：00时间段，在每个布设点选取树龄为9年的京艳桃树(2株)采集相同数量的桃叶(共20片)、 使用GPS定位，详细记录采样点编号、 地理坐标和其所处空间类型等信息。同时，在矿区周围每个布设点处选取5～10个分样点采集土壤，采样深度为0～40 cm，等量多点采集后将其均匀混合在一起，从中抽取1 kg混合样品封装，编号带回实验室。

1.3 遥感数据及处理

遥感数据来自欧空局哥白尼数据中心(https://scihub.copernicus.eu/dhus/#/home)，下载2021年七月份云量在10%以下的Sential-2影像。Sential-2具有13个光谱波段，在红边范围内含有三个波段的数据，有益于监测指数健康信息。地面分辨率分别为10，20和60 m。利用ENVI对影像进行辐射定标、 大气校正、 几何校正、 拼接、 裁剪、 重采样，预处理后得到研究区地表反射率数据。

1.4 叶片光谱反射率与重金属含量测定

叶片光谱的采集使用ASD FieldSpec 4 便携式地物光谱仪，波长范围为350～2 500 nm，光谱分辨率为3 nm，光谱带宽为2.2 nm，采样时间为200 ms。测定前以白色参考板对光谱仪进行标准化，测量过程中，选择50 W的卤素灯为光源，光源照射方向与样品呈15°夹角，叶子置于水平台距光源30 cm，距探头5 cm处。采集每个样点20片桃叶，剔除350～399 nm信噪比低、 噪声大的边缘波段，最后取平均反射率作为该样点桃叶实际光谱反射数据。

将采集的土壤样品碾碎、 风干，过10目尼龙筛去除杂质后，置于60 ℃烘箱烘干后，再次过16目尼龙筛，密封装袋后用于实验室化学分析。土壤平均pH值为5.78，平均有机质含量为17.62 g·kg-1，黏粒含量为25.6%。其中，土壤As，Cr和Zn采用电感耦合等离子体质谱仪进行测试，使用xSPECTRAA-220型原子吸收光谱仪测量土壤Cd，Pb和Cu含量。采用陈同斌等[8]提出的北京市土壤重金属的背景值作为本研究的背景值。如表1所示，平谷区土壤中重金属平均含量与背景值比值大小依次是Cr>As>Cu>Cd>Pb>Zn。根据王威等[9]对平谷区重金属危害程度的研究，综合考虑比值与危害程度的大小，选取重金属Cd，AS和Pb元素作为反演对象。

表1 土壤重金属含量描述性统计Table 1 Descriptive statistics of heavy metal content in soil

1.5 方法

1.5.1 特征变量选取

当植物中的重金属浓度非常低时，反演难度大[10-11]。单一的影像光谱反射率很难突出其差异，为增强光谱响应特征波段，采用一阶微分(first deribative，FD)、 二阶微分(second derivative，SD)、 倒数对数(reciprocal logarithm，RL)、 连续统去除法(continum removal，CR)等变换对光谱数据进行预处理，增强有效波谱信息。为确定特征波段，采用相关性分析，选取在0.01水平下极显著相关的波段作为重金属特征波段。

1.5.2 基于光谱特征的植被指数计算

由于重金属污染对植被光谱特征有一定的影响，植被指数能较好的反应生长状况。依据前人研究，重金属对植被的胁迫主要体现在影响叶绿素的合成和阻碍植被中水分输送。固选择红色归一化植被指数(NDVI705)、 类胡萝卜素反射指数2(CRL2)、 结构不敏感色素指数(SIPI)、 归一化水指数(NDWI)以及基于特征波段构建的植被指数构建模型反演叶片的重金属含量。

2 结果与讨论

2.1 桃叶光谱对Cd元素的响应

将采集的叶片光谱曲线按照重金属含量分为两组，土壤中Cd，As和Pb含量小于背景值的组，代表正常叶片；其余样本的光谱曲线，代表受重金属胁迫的叶片。将两组光谱曲线分别取其平均值，图2显示了正常桃叶与受重金属胁迫桃叶光谱反射率曲线的差异。桃叶总体变化趋势基本一致，在450～750 nm范围内有两个明显的光谱吸收带和反射波峰与叶绿素吸收蓝光、 红光，反射绿光相关，由于植物吸收水分在1 450和1 900 nm形成两个吸收带，具有典型植被光谱曲线的特征。但受重金属胁迫叶片的平均光谱反射率高于正常叶片，且740～1 300 nm波长范围内差值明显，其原因是植被遭重金属胁迫时，叶绿素大量减少，叶黄素、 叶红素相对增加，导致绿光波段反射率变高。

图2 受重金属胁迫叶片与正常叶片反射光谱曲线(a) 及其一阶微分(b)

植被的蓝边、 黄边和红边位置可以反映植被的生长态势和健康状况，红边位置可以作为植物受污染后的重要指示波段。提取680～750 nm范围内反射率一阶导数的最大值作为红边的斜率，其最大值处所对应的波长代表红边位置，该波段范围内一阶导数所包围的面积作为红边面积。同理，用相同的方法确定蓝边(490～530 nm)及黄边(550～580 nm)的位置。由表2显示，受重金属胁迫的叶片相较于正常叶片红边位置发生了“蓝移”现象，蓝边和黄边位置变化不大，表现出了较强的抗干扰能力。红边斜率及红边面积随重金属污染程度加剧而减小。

表2 蓝边、 黄边、 红边参数Table 2 Parameter of blue edge, yellow edge and red edge

2.2 桃叶光谱特征波段的选取

为了更好增强信噪比、 提高分辨率、 突出叶片的光谱特性，提高叶片光谱数据与土壤各参数间的相关性，将原始光谱进行一阶微分、 二阶微分、 倒数对数、 连续统去除法等变换预处理，将不同预处理后的光谱分别与三种重金属含量进行Pearson相关性分析。图3依次表示为一阶微分、 二阶微分、 倒数对数、 连续统去除法与三种重金属含量的相关性曲线。结果表明，预处理后的光谱与三种重金属含量间的相关系数曲线均有明显峰值，而二阶微分处理后光谱与三种重金属含量相关系数曲线峰值分布较多且分散。选取每种变换中六个相关系数较高的峰值波段，相关系数如表3所示。

图3 四种变换形式的光谱数据Fig.3 Pearson correlations of spectral data using four transformations

表3 三种重金属含量与光谱变换特征吸收波段相关系数Table 3 Correlation coefficients between the contents of three heavy metals andthe characteristic absorption bands of four spectral transformation

相关系数筛选出的光谱特征波段虽与重金属含量有相关性，但还需要拟合回归方程确定对重金属含量预测起重要作用的光谱特征波段，通过剔除回归建模中贡献率不显著的波普，得到贡献率显著的波段，如表4所示。

将逐步回归筛选的重复特征波段进行有效整合，最终确定三种重金属的特征光谱变量为：780，945和1 375 nm。

2.3 植被指数的构建与筛选

植被指数可以很好反应植被物理参数，监测植被生长状况。为了充分利用780，945和1 375 nm三个特征波段的信息，便于直接分析特征波段与重金属含量，构建重金属胁迫植被指数(HMSVI)

HMSVI=(R1 375-R945)/(R1 375-R780)

(1)

并选取应用于植被胁迫性探测的指数：红色归一化植被指数(NDVI705)、 类胡萝卜素反射指数2(CRL2)、 结构不敏感色素指数(SIPI)、 归一化水指数(NDWI)。将以上五种植被指数分别与三种不同重金属含量做相关分析，相关系数如表5所示。NDVI705，SIPI和HMSVI与三种重金呈极显著正相关性；NDWI与三种重金属呈负相关性。四种植被指数均能较好反应叶片重金属污染程度，但NDVI705及SIPI与重金属的相关性较弱，CRL2与三种重金属都没有显著相关性；而HMSVI与三种重金属相关性系数均较高，固HMSVI用于三种重金属建模预测效果最佳。

2.4 重金属Cd，As和Pb反演模型

以植被指数为自变量，以重金属含量为因变量，分别建立线性、 二次多项式、 对数、 指数、 幂的回归模型。通过拟合精度(R2)、 均方根误差(RMSE)检验上述植被指数与重金属回归模型的精度，选取R2较高的公式作为各个重金属拟合模型，结果如表6所示。

表4 逐步回归确定的特征波段Table 4 Characteristic band determined bystepwise regression

表5 五种植被指数及其与三种重金属含量相关系数Table 5 Vegetation indexes and its correlation coefficientswith the contents of three heavy metals

表6 三种重金属含量的光谱模型Table 6 Spectral models of three heavy metals

2.5 重金属空间分布的反演与验证

以Sential-2影像为基础，利用三种重金属预测模型进行反演，并将重金属含量划分为5个等级，三种金属含量预测值的空间分布见图4，研究发现：

(1)由反演结果统计可知，Cd含量变化范围为0～0.85 mg·kg-1，As含量变化范围为0～85 mg·kg-1，Pb含量变化范围为0～98 mg·kg-1；Cd，As和Pb三种重金属污染较为严重的区域面积占比分别为: 28.06%，16.82%和19.35%；表明桃林土壤受污染程度较为严重。

(2)图4显示，Cd金属污染严重的区域主要分布在平谷西部刘家店镇、 东北部罗营镇、 东部黄松峪乡及金海湖镇、 东南部南独乐河镇等区域；As金属污染严重的区域主要分布在西部刘家店镇及峪口镇、 东部金海湖镇；Pb金属污染严重的区域主要分布在西部刘家店镇及峪口镇、 东部金海湖镇、 东南部南独乐河镇等区域。三种重金属高值区均大面积的分布在西部万庄尾矿库及刘店尾矿库、 东部黄松峪乡及金海湖镇附近，西部尾矿区比东部尾矿区重金属污染更为严重。

随机选取150组反演值和实测数据，利用均方根误差(RMSE)和决定系数R2对三类金属反演模型的结果进行精度分析。RMSE越接近0，R2越大，表示模型反演效果越好，精度越高。结果如图5所示，三种重金属预测模型RMSE分别为0.019，0.673和5.412，R2分别为0.668，0.794和0.834。表明三种重金属反演模型能准确反映北京平谷区桃林三种重金属空间分布。

3 结 论

以北京平谷桃林区为研究对象，通过光谱变换对桃叶片光谱数据进行处理，选取特征波段，构建植被指数，基于R2最大准则确定了三种金属元素的反演模型，并利用Sential-2影像对土壤重金属含量进行反演及空间分布成图，得到了以下结论:

(1)研究区桃叶遭重金属胁迫时，叶绿素大量减少，绿光波段反射率变高，其平均光谱反射率高于正常叶片，740～1 300 nm波长范围内差值明显。随重金属污染程度加剧，红边位置发生了“蓝移”现象，蓝边和黄边位置变化不大，红边面积及红边斜率减小。

(2)将桃叶光谱曲线进行光谱特征变换后，筛选确定780，945和1 375 nm为特征波段，构建植被指数HMSVI。HMSVI相比于NDVI705，CRL2，SIPI，NDWI与重金属Cd，As，Pb的相关性更高。

图4 三种重金属元素含量空间分布图Fig.4 Spatial distribution of three heavy metal elements

图5 Cd，As和Pb重金属含量预测模型精度验证Fig.5 Accuracy verification of heavy metal content prediction model for Cd, As and Pb

(3)以HMSVI为自变量，土壤Cd，As，Pb元素含量为因变量建立线性、 二次多项式、 对数和指数形式的预测模型，其最终预测模型分别为Y=0.44X+0.193，Y=7.436lnX+13.161，Y=-15.359X+13.583X2+23.541，R2均达到0.6以上，RMSE分别为0.22，1.394和2.403，模型稳定、 适应性强。

(4)利用实测数据验证了反演结果的合理性, 结果表明：北京平谷区Cd，As和Pb含量变化范围分别为0～0.85，0～85和0～98 mg·kg-1；三种重金属高值区均大面积的分布在西部万庄尾矿库及刘店尾矿库附近，东部黄松峪乡及金海湖镇；西部尾矿区比东部尾矿区重金属污染更为严重，且在东北部罗营镇存在较为严重的Cd污染。金属反演结果精度结果显示：RMSE分别为0.019，0.673和5.412，R2分别为0.668，0.794和0.834，三种模型均具有一定的预测能力。

结果表明构建的HMSVI指数，能较好的反演出桃林区土壤重金属Cd，As和Pb元素的分布，可用于监测桃林土壤重金属污染，并为植物覆盖区间接反演土壤重金属污染研究提供一定的科学支持。