玉米叶片铜铅污染元素种类光谱判别的EC-PB规则

吴 兵，杨可明，高 伟，李艳茹，韩倩倩，张建红

中国矿业大学(北京)地球科学与测绘工程学院，北京 100083

引 言

工业化和城市化发展以及一些自然或人为干扰，使人类赖以生存的土壤中重金属浓度逐渐增大，当土壤中重金属含量超过了当地环境背景值时，其累积效应会影响土地承载力，最终导致土壤重金属污染[1-2]，引发生态环境恶化和粮食安全危机。重金属污染物带有毒性并且不可降解，一旦通过生态系统或食物链进入人体，就会对人类身体健康构成严重威胁[3-4]。传统的重金属污染检测方法通常是对野外采集的样品开展实验室化学分析，该方法测量精度高、准确性强，但由于检测环节多、耗时长、成本高，很难快速获取大区域的污染物含量及分布信息[5]。快速判别重金属污染的元素种类以及污染程度，为污染区的针对性治理计划及时提供依据，成为当前亟待解决的关键步骤[6-7]。

由于高光谱遥感具有光谱信息量丰富、技术成本低、理化反演能力强，分析操作实时便捷、可非接触式无损检测等优点[8]，从而如何将光谱分析技术应用于重金属污染检测已成为了现今遥感领域研究的一个热点，并取得了一定进展。Hou等利用可见-近红外光谱数据分析山东省邹城市某煤矿周边的土壤重金属浓度，发现利用光谱数据建立的偏最小二乘回归模型预测土壤中Cu和Pb重金属浓度精度较高[9]。也有一些学者研究表明，土壤中的重金属能对植物的生理结构特征产生影响，当土壤中的重金属含量过高时，会导致植物叶片变黄、分枝、落叶甚至死亡[10]，从而改变植物的光谱特性。在农作物重金属污染方面，杨可明等利用玉米叶片光谱的分数阶微分(fractional order differential，FOD)信息构建玉米叶片光谱红边位置的铜铅敏感指数集群，实现铜铅胁迫识别[11]；黄钟霆等以湖南省某典型关停锰矿区为研究对象，采集矿区周边的农作物及其对应的土壤样品，探讨土壤及对应农作物间重金属迁移规律，并进行生态风险评价[12]。上述研究在重金属元素含量预测以及识别方面均得到了较好的结果，但光谱处理繁琐，算法较为复杂，计算量大，不易复现。

以不同铜铅离子浓度梯度胁迫生长下的玉米叶片光谱为研究对象，使用ASD FieldSpec 4便携式地物光谱仪获取不同胁迫梯度下穗期玉米叶片光谱数据，基于常规的包络线去除(continuum removed，CR)，光谱比值(spectral ratio，SR)以及分数阶微分(FOD)光谱预处理方法构建铜铅元素识别指数(copper and lead identification index，CLI)，联合欧式聚类(euclidean cluster，EC)与垂直平分线(perpendicular bisector，PB)法建立二维(2D)与三维(3D)铜、铅污染元素判别规则，实现依据被污染玉米叶片光谱的铜铅元素准确识别，为重金属元素快速识别，污染区域治理提供新思路。

1 实验部分

1.1 试验设计

玉米是我国重要的农作物之一，也是我国种植面积最大的农作物。为了探究重金属对农作物的影响，试验以盆栽种植的玉米植株为研究对象，设置CuSO4·5H2O和Pb(NO3)2含量分别为0，50，100，150，200，300，400，600，800和1 000 μg·g-1的污染胁迫梯度土壤用于培育植株生长，培育场所室温可控，水气环境适宜，营养元素充足，不同重金属胁迫梯度设置3个平行试验样本组。光谱数据采集在暗室进行，使用波段范围为350～2 500 nm的ASD FieldSpec 4便携式地物光谱仪获取不同胁迫梯度下穗期玉米叶片光谱数据，每个样本测量3次光谱后剔除异常信息取平均光谱，使用电感耦合等离子发射光谱仪测定玉米叶片Cu2+和Pb2+含量(如表1所示)并用于研究结果验证。

表1 不同浓度梯度下玉米叶片中Cu2+，Pb2+含量

本试验共设置9个不同铜铅浓度梯度胁迫下生长的玉米植株，每个浓度梯度设置3个平行组，共54组样本数据，由Python中的sklearn库随机划分样本数据的70%(38组)为训练集数据，30%(16组)为验证集数据。

1.2 光谱变换

(1)包络线去除(CR)。CR也称连续统去除，是一种突出光谱曲线的吸收和反射特征，并将反射率归一化的光谱预处理方法。根据Clark提出的外壳系数法[13]，通过计算光谱曲线上的最大极大值点作为包络线的一个端点，计算该点与波长增长方向各个极大值点连线的斜率，以斜率最大点作为包络线的下一个端点，以此循环，将所有端点连接，形成光谱曲线的包络线，用实际光谱曲线的反射率值除以包络线上相应波段的反射率值即可得到原始光谱曲线去除包络线后的结果，如式(1)和式(2)所示。

(1)

(2)

式中，Si为波段i的包络线去除值，K为极大值始点和端点之间的斜率，ρi为波段i的原始光谱反射率，ρs和ρe为极大值始点和端点的原始光谱反射率，λs和λe为ρs和ρe对应的波段值。本研究中所有玉米叶片光谱数据均进行了CR处理。

(2)光谱比值(SR)。SR最早应用于消除遥感图像中的阴影。随着高光谱技术的成熟，比值技术在混合高光谱信号解混中表现出巨大潜力，其核心思想是通过相同波段不同光谱数据反射率的比值进行SR变换，SR技术可以抑制作为分母的光谱特征，而突出分子光谱的影响。本研究中为了突出重金属元素铜铅的光谱特征，以相同环境下受不同程度铜铅污染的玉米叶片光谱作为分子，无污染的正常的玉米叶片光谱作为分母进行SR处理。

(3)分数阶微分(FOD)。光谱微分可以消除光谱数据之间的系统误差，减弱大气辐射、散射和大气吸收对待测物体光谱的影响，突出待测光谱曲线特征吸收峰和斜率的微小变化，光谱的FOD是整数阶微分的推广，根据Grünwald-Letnikov定义形式[14]，FOD计算方式如式(3)所示。

(3)

式(3)中，λ表示光谱曲线某一波段值，[m,n]为波段区间，λ∈[m,n]，h为光谱采样间隔，a为常数参数，q为任意阶数，Г为Gamma函数，其积分定义如式(4)所示。

(4)

本研究中对经过CR、SR处理后的玉米叶片光谱按照0.1的间隔进行0～2的各阶次FOD处理。

1.3 铜铅污染元素判别特征点(CLDFP)

作物中所含的重金属元素属于微量元素，在光谱曲线中即使经过了增强处理也很难仅通过光谱曲线进行重金属元素识别，因此需要建立一种判别规则对作物样本中的重金属元素进行准确的识别。本研究基于改进红边比值指数(modified red edge simple ratio index，MSR)[15]建立一种铜铅元素识别指数(CLI)如式(5)所示。

(5)

式(5)中，ρ(λa)和ρ(λb)分别表示波长为λa和λb处的光谱值。本研究中将经过CR-SR-FOD处理后的各组玉米叶片光谱代入CLI中，挑选出与铜铅元素种类相关性最强的三个阶数的CLI，分别为CLI1，CLI2和CLI3。以CLI1和CLI2为坐标分量构建2D坐标系，(CLI1为x轴，CLI2为y轴)；以CLI1，CLI2和CLI3为坐标分量构建3D坐标系，(CLI1为x轴，CLI2为y轴，CLI3为z轴)构建铜铅元素判别特征点(copper and lead discriminantfeature points，CLDFP)，如式(6)和式(7)所示。

CLDFP2D=(CLI1，CLI2)

(6)

CLDFP3D=(CLI1，CLI2，CLI3)

(7)

1.4 铜铅污染元素判别规则

欧式聚类(EC)[16]即基于欧几里得距离的一种聚类方法，针对CLDFP需要确定不同元素各样本点内最大的两个欧式距离d1和d2，2D与3D坐标系下的欧式距离公式如式(8)和式(9)所示。

(8)

式(8)中，x1，y1，x2和y2分别为2D坐标系下不同污染元素各样本点内最大欧式距离的CLI值。

(9)

式(9)中，x1，y1，z1，x2，y2和z2分别为3D坐标系下不同污染元素各样本点内最大欧式距离的CLI值。

在2D坐标系中，以最大欧式距离d2D1和d2D2为直径作圆(圆心为两点的中点)，将训练集样本分为铜污染和铅污染两类，连接两圆圆心，以圆心连线的垂直平分线(PB)作为铜铅元素的判别规则线(copper and lead discriminant rule lines，CLDRL2D)；同理，在3D坐标系中，以最大欧式距离d3D1和d3D2为直径作球(球心为两点的中点)，将训练集样本分为铜污染和铅污染两类，连接两球球心，在3D坐标系中将垂直平分线进一步推广为垂直平分面，以球心连线的垂直平分面作为铜铅元素的判别规则面(copper and lead discriminant rule planes，CLDRP3D)实现玉米叶片光谱重金属铜铅元素种类的准确识别。

2 结果与讨论

2.1 光谱变换处理与相关性分析

不同胁迫梯度下的玉米叶片原始光谱数据如图1所示。通过图1可以看出不同胁迫梯度的玉米叶片原始光谱曲线均表现出典型植被光谱特征，不同Cu2+和Pb2+浓度胁迫下仅在反射率中表现出部分差异，规律性不强，很难通过原始数据实现铜铅元素类别的区分，因此需要进行不同的光谱变换处理增强光谱特征以及不同污染元素光谱曲线之间的差异性。

图1 不同胁迫梯度下的玉米叶片原始光谱数据

以土壤中CuSO4·5H2O含量为50 μg·g-1下培育的玉米植株叶片光谱数据为例进行不同光谱变换处理，结果如图2所示。通过图2可以看出经过包络线去除(CR)处理后的光谱曲线峰谷特征表现的更加明显；光谱比值(SR)处理后的光谱曲线已无植被的光谱特征，进一步突出了铜元素的光谱特征；1.2阶分数阶微分(FOD)处理后的光谱曲线消除了光谱数据中的无用信息，使得光谱数据中的有用信息能够更有效的被利用。

图2 不同光谱变换处理的光谱

相关性可以反映数据间的关联程度。通过计算不同光谱变换处理后各波段光谱与叶片中重金属铜铅元素种类之间的相关系数，反映光谱信息与铜铅元素种类的关联性。分析得出：不作任何变形处理的原始光谱数据与叶片中铜铅元素种类之间相关系数绝对值的最大值为0.21，平均值为0.09；而CR处理后的光谱信息与叶片中铜铅元素种类之间相关系数绝对值的最大值为0.51，平均值为0.17；SR是在CR基础上的变换，其相关性与CR相同；0.1～2.0各阶次FOD处理后的光谱信息与叶片中铜铅元素种类之间的相关系数如表2所示。

表2 各阶次FOD与重金属铜铅元素种类相关系数

通过观察不同光谱变换处理后各波段光谱数据与叶片中铜铅元素种类之间的相关系数，可以发现经过光谱变换处理后相关系数最大值、平均值与原始数据相比均有较为明显的增加，说明光谱变换处理增加了光谱数据与铜铅元素种类之间的关联程度；不同阶次FOD处理后相关系数各不相同，为了更好地区分铜铅元素种类，需挑选出相关性较高的阶次用于后续研究。

2.2 建立与挑选CLI

训练集样本数据经过CR、SR以及各阶次FOD处理后代入式(5)计算铜铅元素识别指数CLI并计算0.1～2.0各阶次FOD对应的CLI与重金属铜铅元素种类之间的相关性，相关系数随波长变化情况如图3所示。

图3 各阶次CLI与铜铅元素种类相关性

其中，0.1—0.4阶FOD对应的CLI与铜铅元素种类相关系数值在-0.60～0.63之间；0.5—0.7阶FOD对应的CLI与铜铅元素种类相关系数值在-0.70～0.69之间；0.8—0.9阶FOD对应的CLI与铜铅元素种类相关系数值在-0.67～0.66之间；1.0—1.2阶FOD对应的CLI与铜铅元素种类相关系数值在-0.72～0.70之间；1.3—1.5阶FOD对应的CLI与铜铅元素种类相关系数值在-0.67～0.69之间；1.6—2.0阶FOD对应的CLI与铜铅元素种类相关系数值在-0.66～0.67之间。通过各阶次FOD对应的CLI与铜铅元素种类相关系数可以发现，随着阶次的增加，相关性呈现先递增，后递减的趋势，其中相关系数最高的三个阶次分别为1.2阶、0.7阶、1.0阶。

以1.2阶CLI作为CLI1，其构成波段为639 nm(λa)和1 150 nm(λb)，与铜铅元素种类相关系数为-0.719；以0.7阶CLI作为CLI2，其构成波段为567 nm(λa)和393 nm(λb)，与铜铅元素种类相关系数为-0.701；以1.0阶CLI作为CLI3，其构成波段为2 009 nm(λa)和636 nm(λb)，与铜铅元素种类相关系数为0.700。

2.3 识别效果与验证

(1)铜铅污染元素种类的判别规则线(CLDRL)

分别将CLI1和CLI2对应的训练集样本值带入2D坐标系，形成二维铜铅元素判别特征点CLDFP2D，逐点寻找不同元素各样本点内最大两点的欧式距离d2D1和d2D2进行铜污染与铅污染的分类，其中铜污染判别特征点间最大距离点为(1.267，0.701)与(-0.527，0.329)，两点间欧式距离d2D1为1.833，即以圆心为(0.370，0.515)(最大距离两点中点)，半径为0.916(最大距离的二分之一)作铜污染判别特征类；铅污染判别特征点间最大距离点为(-0.424，1.665)与(2.970，1.457)，两点间欧式距离d2D2为3.400，即以圆心为(1.273，1.561)(最大距离两点中点)，半径为1.700(最大距离的二分之一)作铅污染判别特征类，连接两圆圆心，圆心连线的垂直平分线即为CLDRL2D，其表达式如式(10)所示，2D坐标系下铜铅判别特征点与判别规则线分布情况如图4所示。

CLDRL2D：y=-0.863 0x+1.746 6

(10)

通过图4可以明显看出CLDFP2D根据两个判别特征类分成了不同的区域，大部分CLDFP2D可以由CLDRL2D正确区分，总训练集样本为38个，判别正确样本为30个，判别正确率达到了78.95%。

图4 训练集中2D铜铅判别特征点与判别规则线分布

将CLI1和CLI2对应的验证集样本值带入2D坐标系，以相同的CLDRL2D进行判别，结果如图5所示，总验证集样本为16个，判别正确样本为12个，判别正确率为75.0%。

图5 验证集中2D铜铅判别特征点与判别规则线分布

(2)铜铅污染元素种类的判别规则面(CLDRP)

为了增加样本的直观性，对2D判别规则进一步推广，将CLI1，CLI2和CLI3对应的训练集样本值带入3D坐标系，构建三维铜铅元素判别特征点CLDFP3D，逐点寻找不同元素各样本点内最大两点的欧式距离d3D1和d3D2进行铜污染与铅污染的分类，其中铜污染判别特征点间最大距离点为(0.169，1.689，-0.831)与(0.447，0.517，2.266)，两点间欧式距离d3D1为3.323，即以球心为(0.308，1.103，0.717)(最大距离两点中点)，半径为1.662(最大距离的二分之一)作铜污染判别特征类；铅污染判别特征点间最大距离点为(0.659，1.578，3.290)与(2.114，0.670，-0.904)，两点间欧式距离d3D2为4.531，即以球心为(1.387，1.124，1.193)(最大距离两点中点)，半径为2.266(最大距离的二分之一)作铅污染判别特征类，连接两球球心，球心连线的垂直平分面即为CLDRP3D，其表达式如式(11)所示，3D坐标系下铜铅判别特征点与判别规则面分布情况如图6所示。

图6 训练集中3D铜铅判别特征点与判别规则面分布

CLDRP3D：z=-2.268 0x-0.043 0y+2.924 6

(11)

通过图6可以发现CLDFP3D同样根据两个判别特征类形成了不同的区域，判别效果更为直观，由于新增加了一个维度，误差也在相应累加，总训练集样本为38个，判别正确样本为29个，判别正确率为76.32%。

将CLI1，CLI2和CLI3对应的验证集样本值带入3D坐标系，以相同的CLDRP3D进行判别，结果如图7所示，总验证集样本为16个，判别正确样本为12个，判别正确率为75.0%。

图7 验证集中3D铜铅判别特征点与判别规则面分布

经2D、3D坐标系下训练集与验证集的判别结果可知，通过光谱变换与CLI形成的CLDFP包含了其所受铜铅元素种类相关的信息，使得EC-PB法建立的CLDRL2D与CLDRP3D能够较好的区分铜铅污染下的玉米叶片光谱数据，其中2D数据判别精度略高于3D数据；3D数据可从多个角度查看，判别效果更为直观、易懂。

3 结 论

(1)玉米叶片光谱数据中的铜铅元素污染信息难以通过原始光谱曲线进行准确判别，原始光谱数据与铜铅元素种类之间的相关性较低；CR，SR以及FOD光谱变换处理均增加了玉米叶片光谱数据与铜铅元素种类之间的相关性；不同阶次FOD处理后的光谱数据与铜铅元素种类之间的相关性各不相同。

(2)各阶次FOD对应的CLI与重金属铜铅元素种类相关系数各不相同，随着阶次的增加，相关性呈现先递增，后递减的趋势，其中相关系数最高的三个阶次分别为1.2阶，0.7阶，1.0阶，其对应的相关系数为-0.719，-0.701，0.700。

(3)基于EC-PB建立的2D、3D坐标系下的铜铅元素判别规则线CLDRL2D与铜铅元素判别规则面CLDRP3D能够较好地区分玉米叶片光谱数据中的铜铅污染元素种类，在2D坐标系下训练集判别正确率为78.95%，验证集判别正确率为75.0%；在3D坐标系下训练集判别正确率为76.32%，验证集判别正确率为75.0%。其中2D数据判别精度略高于3D数据；3D数据可从多个角度查看，判别效果更为直观、易懂。

利用EC-PB能较为理想地取得玉米叶片污染的铜铅元素判别与识别效果，同时也存在一定的局限性，如模型的普适性、除铜铅外其他重金属元素的污染识别以及野外大面积农作物重金属污染元素种类识别还需要做进一步的研究。