北京市平谷区桃林土壤重金属元素含量反演估测
刘泓君1, 牛腾1, 于强1,*, 苏凯2, 杨林哲1, 刘维1, 王慧媛1
1.北京林业大学林学院, 北京 100083
2.广西大学林学院, 广西 南宁 530005
*通讯作者 e-mail: yuqiang@bjfu.edu.cn

作者简介: 刘泓君, 1998年生, 北京林业大学林学院硕士研究生 e-mail: liuhongjun_keai9@bjfu.edu.cn

摘要

矿产资源开采中产生的废渣废液长期堆存后产生的渗滤液向土壤中扩散易造成周围土壤的重金属污染, 影响作物生长; 人类通过食物链食用含重金属元素的果实后, 会引起神经系统的神经衰弱、 手足麻木, 消化系统的消化不良, 血液中毒和肾损伤等症状; 这种对生态环境和人身安全的污染和损害是十分严重的。 因此如何快速有效摸清矿区周围农作物土壤污染情况尤为重要。 多光谱遥感由于具备光谱分辨率高、 实时无损、 大面积监测等优势, 在突破植被屏障监测土壤重金属上具有巨大的潜力。 以平谷区主要的农作物桃树为研究对象, 利用桃叶的高光谱数据、 土壤采样数据, 分析桃叶光谱曲线的响应特性, 对桃叶反射光谱进行一阶/二阶微分、 标准正态、 连续去统等四种变换, 结合相关分析及多元线性回归模型确定光谱特征变量, 构建植被指数HMSVI; 结果表明HMSVI与土壤中Cd, AS和Pb含量的相关性较常用植被指数高。 运用线型回归方法进行元素含量与植被指数HMSVI建模后, 选取拟合较好的模型, 实现了叶片高光谱与土壤重金属含量的统计建模, 最后利用Sentinel-2遥感影像反演三种重金属含量空间分布, 并对结果进行精度验证。 结果表明: 受重金属胁迫叶片的平均光谱反射率高于正常叶片且红边位置发生了“蓝移”现象。 780, 945和1 375 nm三个波段对三种重金属污染都较为敏感, 利用三个波段构建的植被指数建立的反演模型能较好的用于桃林土壤重金属元素含量预测, 其预测模型分别为 Y=0.44 X+0.193, Y=7.436ln X+13.161, Y=-15.359 X+13.583 X2+23.541, 且具有较好稳定性和适宜性。 空间反演结果表示, 三种重金属高值区均大面积的分布在平谷区刘家店尾矿库、 万庄尾矿库、 金海湖尾矿库附近, 西部相比东部矿区重金属污染更为严重。 研究结果可以为北京平谷区桃林重金属污染的预防与治理提供基础数据支持。

关键词: 土壤重金属; 特征波段; 多光谱遥感影像; 反演; 空间分布
中图分类号:S127 文献标志码:A
Inversion and Estimation of Heavy Metal Element Content in Peach Forest Soil in Pinggu District of Beijing
LIU Hong-jun1, NIU Teng1, YU Qiang1,*, SU Kai2, YANG Lin-zhe1, LIU Wei1, WANG Hui-yuan1
1. College of Forestry, Beijing Forestry University, Beijing 100083, China
2. Forestry College, Guangxi University, Nanning 530005, China
*Corresponding author
Abstract

The leachate generated from the long-term storage of waste residue and waste liquid produced in mining mineral resources diffuses into the soil, which is easy to cause the surrounding soil to be polluted by heavy metals and affects its crop growth. After human beings eat fruits containing heavy metals through the food chain, they will cause neurasthenia of the nervous system, numbness of hands and feet, indigestion of the digestive system, blood poisoning, kidney injury and other symptoms. Then it will pollute and damage the ecological environment and personal safety. therefore, how to quickly find out the situation of soil pollution is particularly critical. With the development of remote sensing technology, multispectral remote sensing has great potential in breaking through the vegetation barrier to monitor heavy soil metals because of its high spectral resolution and real-time non-destructive and large-area monitoring advantages. This study takes peach trees, the main crop in Pinggu District, as the research object. Using hyperspectral data of peach leaves and field soil sampling data, the response characteristics of peach leaf spectral curves were analyzed. The reflectance spectra of peach leaves were transformed by first-order/second-order differentiation, standard normal transformation and continuous de unification. The characteristic variables are determined by correlation analysis and multiple linear regression model, construction of vegetation index HMSVI, the correlation between HMSVI and Cd, As and Pb content is higher than that of common vegetation index. After modeling element content and vegetation index HMSVI by linear regression method, selecting the model with good fitting, the statistical modeling of leaf hyperspectral reflectance spectrum and soil heavy metal content was realized, and the spatial distribution of heavy metal content was retrieved from sentinel-2 remote sensing image, and the results were verified. The results show that: the average spectral reflectance of leaves under heavy metal stress was higher than that of normal leaves, and the phenomenon of “blue shift” occurred. 780, 945 and 1 375 are the most sensitive to heavy metal pollution. The inversion model established by using the vegetation index constructed in three bands can be better used to predict the content of heavy metal elements in peach forest soil. The prediction models are y=0.44 x+0.193, y=7.436ln x+13.161, y=-15.359 x+13.583 x2+23.541 respectively. The spatial inversion results show that the high-value areas of the three heavy metals are widely distributed near the liujiadian tailings pond, Wanzhuang tailings pond and Jinhai Lake tailings pond in Pinggu District. Heavy metal pollution is more serious in the West than in the East. The mapping results can provide basic data support for preventing and treating heavy metal pollution in Taolin, Pinggu District, Beijing.

Keyword: Soil heavy metals; Characteristic band; Multispectral remote sensing image; Inversion; Spatial distribution
引言

近些年来, 矿产资源的开采及工业污水的排放对矿区周边自然环境的可持续发展和人体健康造成了严重威胁。 特别在种植经济作物的地区, 当重金属含量严重超标后进入种植物体内, 进而造成植物光合作用减弱、 叶片呈现棕色斑块等症状。 人类直接或间接食用此类作物果实后, 易引起重金属中毒, 重则危害生命[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.2 样本采集

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

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

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)Fig.2 Reflectance spectra (a) and first derivative transformation (b) of leaves under heavy metal stress and normal leaves

植被的蓝边、 黄边和红边位置可以反映植被的生长态势和健康状况, 红边位置可以作为植物受污染后的重要指示波段。 提取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 and the characteristic absorption bands of four spectral transformation

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

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

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

2.3 植被指数的构建与筛选

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

HMSVI=(R1375-R945)/(R1375-R780)(1)

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

表5 五种植被指数及其与三种重金属含量相关系数 Table 5 Vegetation indexes and its correlation coefficients with the contents of three heavy metals
2.4 重金属Cd, As和Pb反演模型

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

表6 三种重金属含量的光谱模型 Table 6 Spectral models of three heavy metals
2.5 重金属空间分布的反演与验证

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

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

(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。 表明三种重金属反演模型能准确反映北京平谷区桃林三种重金属空间分布。

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

3 结论

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

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

(2)将桃叶光谱曲线进行光谱特征变换后, 筛选确定780, 945和1 375 nm为特征波段, 构建植被指数HMSVI。 HMSVI相比于NDVI705, CRL2, SIPI, NDWI与重金属Cd, As, 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元素的分布, 可用于监测桃林土壤重金属污染, 并为植物覆盖区间接反演土壤重金属污染研究提供一定的科学支持。

参考文献
[1] Mamattursun Eziz, Ajigul Mamut, Anwar Mohammad, et al(麦麦提吐尔逊·艾则孜, 阿吉古丽·马木提, 艾尼瓦尔·买买提, ). Acta Geographica Sinica(地理学报), 2017, 72(9): 1680. [本文引用:1]
[2] LIANG Ya-ya, YI Xiao-yun, DANG Zhi, et al. (梁雅雅, 易筱筠, 党志, ). Journal of Agro-Environment Science(农业环境科学学报), 2019, 38(1): 103. [本文引用:1]
[3] JIANG Yu-zhuo, WANG Xue-mei, ZHOU Ying, et al(蒋玉琢, 王雪梅, 周颖, ). Environmental Pollution & Control(环境污染与防治), 2020, 42(11): 1392. [本文引用:1]
[4] CHENG Yong-sheng, ZHOU Yao(成永生, 周瑶). The Chinese Journal of Nonferrous Metals(中国有色金属学报), 2019, 39(12): 3880. [本文引用:1]
[5] Anne N J P, Abd-Elrahman A H, Lewis D B, et al. International Journal of Applied Earth Observations & Geoinformation, 2014, 33: 47. [本文引用:1]
[6] LIU Cong, YANG Ke-ming, XIA Tian, et al(刘聪, 杨可明, 夏天, ). China Environmental Science(中国环境科学), 2017, 37(10): 3952. [本文引用:1]
[7] YANG Ling-yu, GAO Xiao-hong, ZHANG Wei, et al(杨灵玉, 高小红, 张威, ). Chinese Journal of Applied Ecology(应用生态学报), 2016, 27(6): 1775. [本文引用:1]
[8] CHEN Tong-bin, ZHENG Yuan-ming, CHEN Huang, et al(陈同斌, 郑袁明, 陈煌, ). Environmental Science(环境科学), 2004, 25(1): 117. [本文引用:1]
[9] WANG Wei, YAN Guang-xin, WANG Li-fa, et al(王威, 闫广新, 王立发, ). Mineral Exploration(矿产勘查), 2019, 10(2): 344. [本文引用:1]
[10] HE Jun-liang, ZHANG Shu-yuan, ZHA Yong, et al(贺军亮, 张淑媛, 查勇, ). Remote Sensing Technology and Application(遥感技术与应用), 2015, 30(3): 407. [本文引用:1]
[11] WANG Jin-feng, WANG Shi-jie, BAI Xiao-yong, et al(王金凤, 王世杰, 白晓永, ). Spectroscopy and Spectral Analysis(光谱学与光谱分析), 2019, 39(12): 3873. [本文引用:1]