张家港市农田土壤重金属含量高光谱遥感监测模型构建

钱家炜，刘晓青，张静静，周卫红，2，李建龙，*

(1.南京大学 生命科学学院 应用生态研究所，江苏 南京 210093； 2.江苏科技大学 苏州理工学院，江苏 张家港 215600)

近年来，随着经济和城市化的发展，重金属污染问题日益严重。重金属污染会严重影响土壤的物理、化学性质和生物学特性，还会抑制生物酶的活性[1]，并通过土壤-作物系统危害农产品质量安全与人体健康[2]。重金属含量预测是预防、评估和治理重金属污染的关键环节。目前，常规的重金属含量预测技术多是通过采集样本和实验室化学分析，借助于ArcGIS软件进行空间插值完成的。当采集的样本数量有限时，其结果并不能代表整个研究区；但大规模采集样本又不切实际，耗时，且费用昂贵，而且实验室使用强酸等物质进行分析时也存在安全隐患[3]。遥感作为一种新型的探测手段，可以方便、快捷、动态地获得空间上连续的地物光谱信息[4]，其中，可见光和近红外反射光谱已被广泛用于测量土壤的物理和化学性质[5]。

国内外许多学者利用高光谱遥感监测对土壤重金属含量进行研究，并取得了很好的成果。Ma等[6]利用超限学习机、支持向量机等方法建模，估测土壤中5种重金属的含量。刘华等[7]根据EO-1(Earth Observing-1，地球观察者1号)卫星Hyperion高光谱仪的波段设置，结合实测数据，基于光谱反射率运用偏最小二乘法很好地估测了土壤重金属含量。袁中强等[8]分析了环境一号卫星上搭载的超光谱成像仪的地表反射率数据与实测的土壤重金属元素含量的关系，建立了预测土壤中5种重金属含量的原始反射率、一阶导数、倒数对数3个回归模型。St Luce等[9]利用偏最小二乘法建立可见近红外光谱与土壤中总锌、总镉的模型。Choe等[10]研究发现，可见和近红外范围内610 nm和500 nm光谱反射率的比值、2 200 nm处的吸收面积，以及2 200 nm处吸收特征的不对称性分别与Pb、Zn、As浓度存在显著相关性。综上所述，利用高光谱与土壤重金属含量建立的反演模型可以预测重金属含量。但前人的研究多是用不同变换方式的光谱反射率建立多个模型，然后通过比较精度筛选出最优模型。与此不同，本文拟在建立模型的过程中，分析5种变换方式与8种重金属含量的相关性，先筛选出每种重金属对应的最优变换方式，再建立模型。为此，特以长三角地区张家港市为研究对象，选取对人体有伤害、且存在范围比较广的8种重金属做总结性研究，在室内用ASD地物光谱仪(美国Analytica Spectra Devices公司)获取土壤高光谱数据，与经电感耦合等离子体原子发射光谱(ICP-AES)法测量的土壤As、Cd、Cr、Cu、Zn、Ni、Pb、Hg含量相结合，光谱经过各种形式的变换，根据相关性分析选出每种重金属对应的敏感特征波段，利用逐步回归法建立土壤重金属的高光谱反演模型，预测土壤重金属含量。

1 材料与方法

1.1 研究区概况

张家港市位于江苏省苏州市北部，长江下游南岸，是沿海开放地区和长三角经济带交汇区较早开放的工业城市，地处120°21′～120°52′E、31°43′～32°02′N，年平均气温17 ℃，年降水量1 200 mm以上，年日照时数在1 600 h以上，属北亚热带南部湿润性气候。张家港市现有耕地面积30 646.7 hm2，常年种植水稻18 133.3 hm2、小麦19 000.0 hm2、油菜1 600.0 hm2、蔬菜3 666.7 hm2、果品1 400 hm2。张家港市下辖8镇2区，分别为杨舍镇、塘桥镇、金港镇、锦丰镇、乐余镇、凤凰镇、南丰镇、大新镇，及常阴沙现代农业示范园区、双山岛旅游度假区。

图1 采样点分布图Fig.1 Distribution of soil samples

1.2 样品采集与制备

在地理空间数据云上下载张家港市的TM数据，并将数据在ENVI 5.3软件中进行监督分类，提取农田土壤，以混淆矩阵法算出分类结果精度为0.87。在张家港市8个镇和1个现代农业示范园区的农田土壤中均匀地设置23个采样点。依据样点分布(图1)，用GPS精准定位后去野外采集0～20 cm的农田土壤样品。将收集的土壤样本带回实验室自然风干，粗略去除砂砾和植物残体，研磨、过100目塑料尼龙筛后混合均匀。采用四分法取样，一式2份，一份用于测定重金属含量，另一份用于光谱测定。土壤As、Cd、Cr、Cu、Zn、Ni、Pb用硝酸-盐酸-高氯酸消解，Hg经硝酸-盐酸混合试剂在沸水浴中加热消解，消解后均采用ICP-AES法[11]进行测定。

1.3 光谱测定

在室内条件下，采用ASD光谱仪测定经过处理的土壤样品的光谱反射率，光谱波长范围在350～2 500 nm。将土壤盛装在直径10 cm、深3 cm的黑色器皿中。在进行光谱测定之前，先对土壤表面做刮平处理[12]。光谱测定在暗室中进行，以功率50 W的卤素灯作为唯一光源，测量方法严格按照有关规范操作：光源入射角为45°，光源距离土样表面中心30 cm，探头距离土样15 cm。设置40 cm×40 cm的参考白板，用于获取绝对反射率。每个土样测量10次光谱曲线，取平均值，以减少误差，提高精度。利用View Spec Pro软件剔除异常曲线后，取光谱反射率平均值作为初始反射率光谱值[13]。

1.4 光谱变换

在对目标进行光谱测定的过程中，背景、大气、仪器、光照条件等因素影响都可能覆盖目标自身的光谱信息特征。因此，尝试进行多种形式的光谱变换，以减弱背景噪声的影响[4]，增强目标光谱信息，提高信噪比。为描述方便，本文将经过平滑处理后的反射光谱称为原始反射光谱。对原始反射光谱分别进行一阶导数、倒数一阶导数、倒数的对数一阶导数、平方根一阶导数和连续统去除[14]共5种形式的变换处理。

1.5 模型建立与精度分析

1.5.1 相关性分析

在高光谱数据与农田土壤重金属含量建模的过程中，筛选与重金属含量相关性高的波段作为特征敏感波段，相关性越高，波段响应越敏感[15]。将经过变换后的光谱反射率与土壤重金属含量进行皮尔逊(Pearson)相关分析，考虑不同重金属光谱反演的需要，根据每种重金属的光谱特征，筛选出相关系数通过P<0.05水平显著性检验的波段，作为建立高光谱定量估算模型的自变量。

1.5.2 建模方法

逐步回归的基本思想简述如下：(1)将变量逐个引入模型；(2)每引入一个解释变量后都要进行F检验，并对已选入的解释变量逐个进行t检验；(3)当原变量因为新解释变量的引入变得不再显著时，将其剔除，以确保每次引入新的变量之前回归方程中只包含主动变量；(4)重复进行这一过程，直到既没有显著的解释变量选入回归方程，也没有不显著的解释变量从回归方程中剔除为止，以保证最后所得到的解释变量集最优[16]。

依据上述思想，利用逐步回归筛选并剔除引起多重共线性的变量[17]，步骤如下：先用被解释变量对每一个所考虑的解释变量做简单回归，然后以对被解释变量贡献最大的解释变量所对应的回归方程为基础，再逐步引入其余解释变量[18]。经过逐步回归，使得最后保留在模型中的解释变量既是重要的，同时又可减小多重共线性[19]。

2 结果与分析

2.1 土壤样品统计分析

土壤重金属含量测定结果见表1。从标准差来看，Hg、Cd、Cr的标准差都小于1 mg·kg-1，其余5种重金属元素的标准差为1.173～5.568 mg·kg-1。从变异系数来看，Cd、Hg的变异系数分别为0.342、0.315，其余的变异系数均在0.1左右，变异程度为中等变异性[20]。从平均值来看，Cd、Hg均超过背景值，其余重金属元素的平均含量则均低于背景值，说明该区域土壤环境质量总体良好；但从最大值看，有部分样点的土壤重金属含量与背景值相近，部分样点的Cd、Hg、Cu、Zn含量已超出背景值。根据土壤重金属污染的单因子指数法测算，As、Cr、Pb、Ni的污染指数均小于0.7，含量在安全范围内；Zn、Cu的污染指数接近1，存在潜在危害；Cd的污染指数为1.324，属于轻度污染；Hg的污染指数为2.724，属于中度污染。因此，张家港市应该加强土壤质量调查与动态监测，以便及时发现并控制耕地土壤的重金属污染问题。

表1 土壤重金属含量统计特征

2.2 光谱变换分析

如图2所示，将原始光谱进行一阶导数、连续统去除、平方根一阶导数、倒数一阶导数、倒数的对数一阶导数变换后发现，除了倒数一阶导数减少了原来的光谱特征信息外，一阶导数、连续统去除、平方根一阶导数、倒数的对数一阶导数变换均能够将原始光谱波峰、波谷等信号特征放大，产生更大幅度的波动，增大与背景光谱的差异。其中，一阶导数和连续统去除变换方法都在相同的几个波段处变化剧烈，放大了原始光谱信息，在消除或减弱土壤背景噪声、提高信噪比、增强目标光谱信息方面作用明显，更容易从光谱数据中提取出能够反映土壤特性的微弱信号。

a，原始光谱；b，连续统去除；c，一阶导数；d，平方根一阶导数；e，倒数一阶导数；f，导数的对数一阶导数。a， Raw spectrum; b， Continuum removal; c， First derivative; d， Square root of first derivative; e， Reciprocal of first derivative; f， Logarithm of reciprocal of first derivative.图2 光谱变换结果Fig.2 Spectral transformation result

2.3 土壤重金属含量与光谱反射率的相关性分析

分别计算每种重金属元素含量和对应原始光谱的一阶导数、连续统变换的相关系数，比较不同光谱变换形式对该相关程度的影响，找出特征波段。

将土壤重金属含量与原始光谱的一阶导数进行相关性分析，得到每种重金属与对应光谱反射率的相关系数曲线，并对其进行显著性检验(显著性水平设定为P<0.05)。从图3可以看出，As、Cd、Cr、Pb、Zn、Ni、Hg含量与光谱反射率相关系数曲线中的许多波段都通过了P<0.05水平的显著性检验。其中，As、Cr、Pb、Zn、Ni、Hg的相关系数在900～2 500 nm通过显著性检验的波段比较多；Cd相关系数通过显著性检验的波段数少于除Cu以外的其他6种重金属；Cu与光谱反射率的相关系数较低，未能通过P<0.05水平的显著性检验。

从图4可以看出，Cr、Zn、Hg含量与反射率连续统去除的相关系数通过显著性检验的波段很少；Ni相关系数通过显著性检验的波段集中在360～600 nm，且与反射率呈现负相关；As相关系数通过显著性检验的波段集中在355～585、620～780 nm；Pb相关系数通过显著性检验的波段集中在995～1 305、1 490～1 790 nm；Cd、Cu通过显著性检验的波段较多，Cd相关系数通过显著性检验的波段集中在360～800、850～935、1 800～1 840 nm，Cu相关系数通过显著性检验的波段集中在860～1 340、1 450～1 660、2 180～2 485 nm。综上所述，As、Cr、Zn、Pb、Ni、Hg含量与原始反射率一阶导数的相关性高，Cd、Cu含量与原始反射率连续统去除的相关性较高。

2.4 高光谱模型建立与验证

将23个样本分成2组，从第1个样本开始，每5个抽取1个样本用于模型检验；因此，模型构建样本数为18，模型检验样本数为5。建立模型时，选择与重金属含量相关性高的光谱变换形式的反射率作为自变量(x)，以重金属含量作为因变量(y)。在SPSS 20.0中进行逐步回归分析，默认设置F值概率(P)小于0.05进入，大于0.1删除，得到高光谱重金属含量的定量估算模型。

图3 反射率的一阶导数与土壤重金属含量的相关系数Fig.3 Correlation coefficient between first derivative of reflectance and soil heavy metal content

图4 反射率的连续统去除与土壤重金属含量的相关系数Fig.4 Correlation coefficient between continuum removal of reflectance and soil heavy metal content

基于最优光谱形式组合和重金属含量，用逐步回归法建立模型。由表2可以看出，各种重金属和光谱建立的线性模型中P值均小于0.01，说明所建模型可以反映光谱与重金属含量间的关系。As、Cd、Cr和Hg定量估算模型的拟合度(R2)均达到0.7以上，其中Hg的估算模型的拟合度最高(R2=0.777)，Zn、Ni的拟合度达到0.6以上，Cu、Pb的拟合度在0.5～0.6。8种重金属估算模型的F统计量均较大，除了Pb的估算模型F值为9.654以外，其他模型的F值都大于10，其中Cd的估算模型的F值最大，为27.710。综上，建立的8种重金属估算模型可以在一定程度上反映研究区的土壤重金属含量。

如图5所示，用之前的5个验证样本对所构建的模型进行验证，得到模型实际值与验证值的拟合曲线。从拟合曲线可以看出，Cd、Hg估算模型的实际值与验证值的拟合度大于0.8，分别为0.874、0.879，Cr估算模型的拟合度为0.800，As、Cu、Zn、Ni、Pb估算模型的拟合度也都大于0.5，说明8种重金属的定量估算模型验证精度较高，具有预测张家港市农田土壤重金属的能力，能进一步用于大面积遥感估算。

表2 所建立的光谱与重金属含量模型

3 讨论

基于本研究试验数据，张家港市农田土壤重金属As、Cd、Cr、Cu、Pb、Zn、Ni、Hg的含量在0.053～99.258 mg·kg-1，部分样点的Cd、Hg、Cu、Zn含量超出背景值。依土壤重金属污染的单因子指数法测算：Hg的污染指数为2.724，属于中度污染；Cd的污染指数为1.324，属于轻度污染；Zn、Cu的污染指数都接近1，存在潜在危害；As、Cr、Pb、Ni的污染指数小于0.7，重金属含量在安全范围内。

通过土壤重金属与光谱反射率的相关性分析，选取As的敏感波段为2 069、2 201、2 267 nm，Cd的敏感波段为523、1 809 nm，Cr的敏感波段为2 065、2 413、1 634 nm，Cu的敏感波段为2 225、2 431 nm，Zn的敏感波段为2 332、2 060 nm，Ni的敏感波段为2 332、1 922 nm，Pb的敏感波段为2 332、2 038 nm，Hg的敏感波段为2 312、1 573、997 nm。选其用作张家港市土壤重金属含量反演模型的自变量。

8种重金属和光谱建立的线性模型，P值均小于0.05，说明构建的逐步回归模型整体显著。在土壤重金属含量和高光谱波段建立的模型中，Cu、Pb的拟合优度达到0.5以上，Zn、Ni的拟合优度达到0.6以上，As、Cd、Cr、Hg的拟合优度达到0.7以上，拟合优度较高。从模型实际值与验证值的拟合曲线得出，Cd、Hg实际值与验证值的拟合度分别为0.874、0.879，As、Cr、Cu、Zn、Ni、Pb的拟合度也都大于0.5。由此可知，8种重金属估算模型验证精度良好。以上结果说明，对As、Cr、Zn、Pb、Ni和Hg光谱反射率进行一阶求导处理，对Cd和Cu光谱反射率进行连续统去除处理后，所建立的8种重金属光谱模型经验证效果理想，可有效预测土壤重金属含量，是替代传统重金属预测方法的优化选择，可在实际生产中加以应用。

本研究采用多元线性逐步回归方法对8种重金属含量进行了建模分析，但由于只设置了23个采样点，建模集和验证集数量均有一定不足，因此，今后应进一步扩大采样范围以优化模型性能。与已有的基于高光谱方法预测土壤重金属含量的研究相比，李琼琼等[21]对土壤Cu、Pb、Zn含量进行反演研究的结果表明，最小二乘回归方法的建模精度优于多元线性逐步回归；王金凤等[22]运用随机森林、支持向量机、偏最小二乘3种算法对土壤Zn含量进行建模，发现基于二阶微分变换光谱的随机森林算法准确度最高；许吉仁等[23]用支持向量机方法对土壤Cd含量建立高光谱监测模型，R2为0.947。为了探究高光谱方法预测土壤重金属含量的更优方法，未来可考虑引入多元散射校正、标准正态化等光谱预处理方法和偏最小二乘、随机森林的建模方法来进行比选优化。