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

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

引 言

近些年来环境污染，特别是重金属污染在工业化和城市化快速发展地区尤为严重，因此重金属污染监测与治理技术等已成为当今环境保护领域的研究热点，其中农业重金属污染更是备受关注[1]。农业重金属污染主要来源于农药和近农田区域采矿活动等[2]。农田土壤受重金属污染程度不易被察觉，而重金属离子可以被农作物吸收并富集，进而能造成粮食减产、毒害食物链、威胁食品安全、危害人体健康。

如何实时快捷有效识别重金属污染、更加准确甄别重金属污染的类别与程度是当前亟待研究与突破的关键技术。传统的重金属污染监测方法通常是基于野外取样后，在室内采用无焰原子吸收分光光度法(FAAS)、原子吸收光谱法(AAS)等进行测定，虽然检测精度高，但也具有耗费大量人力、物力和财力，损害被测对象，所检测仪器昂贵、维护成本高、操作程序复杂和检测时间长等缺点[3]；而今高光谱遥感技术具有图谱合一、光谱连续，高效、无损、低成本等优势，正在被广泛应用于重金属污染甄别与监测的探索研究中[4]。Liu等[5]通过提取植被红边光谱曲线，构造归一化红边指数和红边叶绿素指数，用于识别植被受重金属镉(Cd)的胁迫；乔晓英等[6]以矿区周边玉米、苦菜等作为研究对象，构建植物中重金属Cd、铅(Pb)质量比与植物光谱红边、蓝谷等特征参数的关系模型，可间接反演植物中重金属的含量；付萍杰等[7]对叶片光谱进行经验模态分解并结合自相关函数、多重分形理论构建Cu、Pb含量反演的多元线性回归模型，成功识别了玉米叶片中重金属类别及含量；朱叶青等[8]通过获取不同生育期、不同Cu污染下的植被叶片光谱信息，采用7个特征波段结合光谱角方法研究Cu污染下的叶片光谱特征，表明Cu污染叶片与健康叶片的光谱存在明显差异，同时证明了Cu污染后的叶片内部结构更加紊乱。

作为植被生长所必需的微量元素Cu，也属于一种重金属，当浓度较低时有利于植被生长，浓度过高就会抑制植物的光合作用，对植被生长产生胁迫作用，轻则引起植物代谢紊乱，生长受阻，重则导致植物死亡。已有报道当植被受到重金属污染时，多数会长势矮小，叶片失去绿色等症状[9]。本工作提出固有波长尺度分解(intrinsic wavelength-scale decomposition，IWD)的方法并运用于甄别农作物重金属污染的高光谱信息处理与分析中。依据受Cu污染的玉米叶片光谱数据，基于Hankel矩阵构建、奇异值分解(singular value decomposition，SVD)等建立玉米叶片Cu污染信息预测的IWD-Hankel-SVD模型，该模型可有效预测玉米叶片中铜离子(Cu2+)含量，能为植物Cu污染信息预测、甄别提供一种新方法。

1 理论与算法

1.1 固有波长尺度分解

固有波长尺度分解(IWD)引自于固有时间尺度分解(intrinsic time-scale decomposition，ITD)。ITD是由Mark G Frei与Ivan Osorio首次提出[10]，相对于经验模态分解(empirical mode decomposition，EMD)和小波变换等，ITD是一种新的时频分析方法[11]。经过ITD分解后的信号分量具有相对较为完整的时频信息，能够反映信号频率变化等优点。本工作用波长(wavelength，W)替代ITD中的时间(time，T)，探索形成光谱信号的IWD处理技术，并引入到高光谱数据的信息处理与分析中。IWD可以将波长序列(w)自适应分解为几个固有旋转分量(proper rotation component，PRC)和一个趋势分量(r)，每个旋转分量只需要通过迭代一次就可以得到。对于一条准备分解的光谱曲线Xw，有极值点Xk，其基线提取算子定义为L，可将基线分量Lw与准备分解的光谱曲线Xw分开，再定义一个固有旋转分离算子为H，于是有

Xw=LXw+(1-L)Xw=Lw+Hw

(1)

其中，Lw=LXw为基线分量，Hw=(1-L)Hw=HXw为PRC分量。

对于给定的光谱曲线，IWD的处理步聚如下：

(1)提取光谱曲线Xw上的极值点Xk(k=1,2,…,N，N表示极值点的个数)，令τk是Xw局部极值点所在的波长集合，定义τ0=0作为Xw的起始波长，并且使Xw和Lw在τk处的值分别为Xτk和Lτk，为了表示方便，将Xτk和Lτk分别定义为Xk和Lk。

假设Lw和Hw被定义在波段[τ0,τk]上，Xw被定义在波段[τ0,τk+2]上，并且，假设Lw是Xw在区间(τk,τk+1]的仿射线性逼近，即Lw=mXw+nw，w∈(τk,τk+1]。因为Lw必须满足局部极值点的上述边界条件，所以Lw被表示为以下形式

(2)

假设光谱曲线的基本趋势足够平滑，则可以忽略各局部极值点间的变化，于是可以得出

(3)

因为Lk+1=Xk+1，于是有

Lk+1=αXk+1+(1+α)Xk+1=

(4)

式(4)中：α为线性缩放因子，可用于调整提取到的PRC幅度，α∈[0,1]，一般情况下，经验值取值0.5。

(2)将式(2)计算得出的基线分量Lw=LXw，按照式(1)计算提取固有旋转分量Hw，即PRCi。

(3)将基线分量Lw作为一个输入光谱再进行下一次分解，重复第(1)和(2)步的分解，当基线分量Lw变为单调或小于某个预设值的趋势分量(r)时，分解终止。

IWD的总分解过程为

Xw=Hw+Lw=HXw+LXw=HXw+(H+L)LXw=

(5)

经过IWD分解后，原始输入的光谱曲线Xw即被分解为多个固有旋转分量(PRCi)和一个单调趋势分量(r)。该方法计算复杂度低，分解精度高，鲁棒性强，有效的避免了波形叠加、模态混叠现象。

1.2 Hankel矩阵

汉克尔(Hankel)矩阵具有很多优良的性能和特殊的性质，因此，作为一种优秀的数学工具被广泛应用于诸多领域的数值模拟计算中。Hankel矩阵的每一条逆对角线上的元素都相等。可选择经IWD处理后较好保留了原始输入光谱信息的PRC，构造Hankel矩阵Za

(6)

式(6)中，1≤a≤J，且1

1.3 奇异值分解与奇异熵

奇异值分解(SVD)即对于任意一个维数的实矩阵A∈Rm×n，都可分解为三个矩阵的乘积，即

A=USVT

(7)

式(7)中，U为单位正交矩阵，即其列向量是单位向量并且相互正交，是奇异值对应的特征向量[13]；V为单位正交矩阵，其行向量是单位向量并且相互正交，也是奇异值对应的特征向量；S为对角矩阵，S=[diag(σ1,σ2,σ3,…,σm),O]，其中O表示零矩阵，它的对角线矩阵为奇异值。

奇异熵是用来度量奇异值所对应用的信号分量信息量的多少，奇异熵的计算公式是

(8)

式(8)中，k为奇异熵所对应的阶次；ΔEi为奇异熵在阶次i时得到的增量，其计算公式为

(9)

1.4 IWD-Hankel-SVD模型

基于IWD，Hankel矩阵和SVD可构建IWD-Hankel-SVD模型，用于光谱的信息分析与特征参量反演等。IWD-Hankel-SVD模型的植物重金属污染下叶片光谱特征信息提取步骤为：

(1)IWD分解。对需要甄别重金属污染的植物光谱数据进行IWD分解，得到PRC1，PRC2，…，以及r若干分量。

(2)分量筛选。对经过IWD分解后得到的各PRC分量与原始输入光谱进行对比，选择能够较好保留原始输入光谱信息的PRCi分量。

(3)特征波段选择。在PRCi分量中选取与原始输入光谱相似的波段作为特征波段进行下一步的处理。

(4)构造Hankel矩阵。利用从第(2)步得到的PRCi分量选择的特征波段，构造Hankel矩阵Z。

(5)SVD分解。对第(4)步构造的Hankel矩阵进行SVD分解，得到奇异值矩阵。

(6)求取奇异熵。通过第(5)步得到的奇异值计算对应的奇异熵，作为表征重金属污染下光谱信息分量的特征。

2 实验部分

(1)实验准备。实验采用盆栽的玉米植株土培法，玉米种子选用“中糯301”。重金属Cu2+胁迫试剂选用纯度较高、干扰杂质较少的分析纯级别CuSO4·5H2O。光谱数据获取采用SVC HR-1024i便携式地物光谱仪，其光谱测量范围为350～2 500 nm，其中350～1 000 nm光谱分辨率为3.0 nm，1 000～1 900 nm光谱分辨率为9.5 nm，1 900～2 500 nm光谱分辨率为6.5 nm。叶绿素浓度相对值的测定采用SPAD-502手持便携式叶绿素测定仪。

(2)植株培养。“中糯301”玉米种植时采用底部有渗水孔的培土花盆，将分析纯级别的CuSO4·5H2O配置成相应梯度的Cu2+溶液拌入培土作为所种植玉米生长的胁迫试剂。培植实验共设置三组平行实验，每组设7个Cu2+胁迫梯度，即空白对照组(ck，不添加胁迫试剂)和100，200，400，600，800和1 000 μg·g-1等6级Cu2+胁迫梯度，分别标注为ck(0)，Cu(100)，Cu(200)，Cu(400)，Cu(600)，Cu(800)和Cu(1 000)。为了提高玉米培养效率，事先对“中糯301”玉米种子进行催芽处理。植株培养期间保持水量充沛，空气流通。2个月后培植结束，测定玉米叶片光谱信息、叶片中Cu2+的含量和叶绿素浓度相对值。

(3)光谱数据采集。在室内采用SVC HR-1024i光谱仪测定实验盆栽玉米在不同浓度Cu2+胁迫下的叶片光谱数据，选用50W卤素灯作为单一光源。光谱采集时分别测量每个平行组各胁迫浓度下玉米老(Old)、中(Middle)、新(New)叶片的光谱3次，每个胁迫浓度获取9组数据，共计63组光谱数据。最终获取的光谱数据是经由去除异常值，求平均得到。

(4)叶绿素浓度相对值测定。同步利用SPAD-502叶绿素仪分别测定采集光谱的叶片叶绿素浓度(重复3次)，最终通过求平均得到了各Cu2+胁迫浓度下玉米叶片叶绿素浓度相对值。

(5)Cu2+含量的测定。在各类被测叶片光谱和叶绿素的数据采集结束后就即时进行叶片样本裁剪和保存处理，为叶片中Cu2+含量测定作准备。最终测定的叶片中Cu2+含量是采用电感耦合等离子发射光谱仪(ICP-OES)进行分析得到。

胁迫梯度以及对应的玉米叶片叶绿素浓度相对值、Cu2+含量见表1所示。由表1可见，随着Cu2+胁迫增大，叶绿素浓度相对值在逐渐减小，叶片中Cu2+含量在逐渐增大。各胁迫梯度的叶片叶绿素浓度相对值与Cu2+含量相关系数为-0.953，p值为0.001<0.01，可见两者呈负相关，负相关性高度显著，说明在重金属Cu2+胁迫下，玉米叶片中叶绿素含量的积累受到影响[14]。

表1 不同浓度Cu2+胁迫下玉米叶片中叶绿素浓度相对值及Cu2+含量Table 1 Relative value of chlorophyll concentration and content of Cu2+ in corn leaves under different concentrations of Cu2+

3 结果与讨论

3.1 IWD处理结果与分析

通过SVC HR-1024i光谱仪测定的不同Cu2+胁迫梯度下玉米叶片光谱曲线如图1所示，受不同浓度的Cu2+胁迫，玉米叶片光谱曲线与对照组ck(0)的健康玉米叶片光谱曲线有所不同，但是不能准确的甄别出受重金属Cu2+胁迫的变化规律。为此，引入IWD对不同浓度Cu2+胁迫下的玉米叶片光谱曲线进行变换，然后构造Hankel矩阵，进行SVD处理，最后求得奇异熵，通过线性拟合建立的IWD-Hankel-SVD模型监测确定重金属Cu2+胁迫下玉米叶片的污染程度。

采用IWD方法分解对照组ck(0)的健康叶片原始光谱，得到PRC1和PRC2的2个固有旋转分量和1个趋势分量r，如图2所示，其中并非每一个分量都包含有丰富的光谱特征信息，有些分量对于光谱特征信息的提取贡献很小。通过与原始输入光谱信息比较，发现PRC1分量较好的保留了原始输入光谱信息，并且对原始光谱中部分波段的反射率有所放大。同样，通过对Cu(100)，Cu(200)，Cu(400)，Cu(600)，Cu(800)和Cu(1000)胁迫梯度下原始叶片光谱进行IWD分解，发现不同胁迫浓度下的固有旋转分量中PRC1分量均较好地保留了原始输入光谱信息，如图3(a—f)所示。因此，根据PRC1分量的光谱保留信息，从中选择与原始输入光谱相似的553～680，681～780，1 266～1 429，1 430～1 631，1 836～1 913和1 914～2 111 nm波段区间作为有效特征波段。

图2 对照组ck(0)玉米叶片光谱及其IWD处理结果Fig.2 Spectrum and the results processed by IWD of corn leaves in ck(0) control group

图3 不同浓度 Cu2+胁迫下玉米叶片光谱IWD处理的PRC1分量Fig.3 Spectral PRC1 components of the corn leaves stressed by the different concentrations of Cu2+ based on the IWD processing

3.2 单变量的IWD-Hankel-SVD模型应用

通过对特征波段553～680，681～780，1 266～1 429，1 430～1 631，1 836～1 913和1 914～2 111 nm对应的PRC1分量分别构造Hankel矩阵，并对该矩阵进行SVD分解，实现IWD-Hankel-SVD模型的应用，最后求取模型的奇异熵用于分析监测效果。各特征波段所对应的奇异熵与玉米叶片中Cu2+含量的相关关系见表2。由表2可见，通过特征波段1 266～1 429和1 836～1 913 nm计算求得的奇异熵与玉米叶片中Cu2+含量的相关系数分别为0.950和0.913，说明经过IWD分解后的波段1 266～1 429和1 836～1 913 nm对应的PRC1分量对Cu污染较为敏感。选取相关系数大于0.9的特征波段计算得到的奇异熵分别与玉米叶片中Cu2+含量进行线性拟合，拟合结果如图4、图5所示，可见：(1)根据特征波段1 266～1 429 nm计算得到的奇异熵与玉米叶片中Cu2+含量线性拟合，得到的模型决定系数R2为0.903，p值小于0.01，相关性高度显著；(2)根据特征波段1 836～1 913 nm计算得到的奇异熵与叶片中Cu2+含量线性拟合，得到的模型决定系数R2为0.834，p值小于0.01，相关性高度显著。因此得出，可通过IWD-Hankel-SVD模型来预测玉米植株受重金属Cu污染的程度。

图4 特征波段1 266～1 429 nm的奇异熵与玉米叶片中Cu2+含量变化关系Fig.4 Relationship between the singular entropy of 1 266～1 429 nm characteristic band and the content of Cu2+ in corn leaves

图5 特征波段1 836～1 913 nm奇异熵与玉米叶片中Cu2+含量变化关系Fig.5 Change relationship between the singular entropy of 1 836～1 913 nm characteristic band and the content of Cu2+ in corn leaves

表2 特征波段的奇异熵与玉米叶片Cu2+含量的相关关系Table 2 The correlation between the singular entropy of characteristic bands and the content of Cu2+ in corn leaves

为了验证依据所提取特征波段以及IWD-Hankel-SVD模型预测植物重金属污染程度的优越性，选择传统的蓝边最大值、红肩最大值、绿峰高度等光谱特征参数(见表3)监测方法以及与采用连续投影变换(SPA)选择的特征波段如图6、图7所示，进行应用结果对比分析。通过蓝边最大值、红肩最大值、绿峰高度的监测结果与玉米叶片中Cu2+含量的线性拟合，得到拟合效果如图8、图9、图10所示，通过连续投影变换(SPA)选择的特征波长与玉米叶片中Cu2+含量的线性关系见表4，对比结果表明，IWD-Hankel-SVD模型的监测能力明显要优于其他传统光谱特征等参数所建立的模型。

表3 玉米叶片的光谱特征参数名称与定义Table 3 Names and definitions of spectral characteristic parameters of corn leaves

表4 SPA选择的特征波长与玉米叶片中Cu2+含量的预测表达式Table 4 The prediction expression of Cu2+ content in corn leaves and characteristic wavelengths selected by SPA

图6 连续投影变换选择特征变量的个数Fig.6 Number of characteristic bands selected by SPA

图7 连续投影变换提取特征波长Fig.7 Characteristic bands extraction by SPA

图8 蓝边最大值方法监测结果与叶片中Cu2+含量的变化关系Fig.8 Relation between monitoring results based on the maximum of blue edge and Cu2+ content in corn leaves

图9 红肩最大值方法监测结果与叶片中Cu2+含量的变化关系Fig.9 Relation between monitoring results based on the maximum of red shoulder and Cu2+ content in corn leaves

图10 绿峰高度方法监测结果与玉米叶片中Cu2+含量的变化关系Fig.10 Relation between monitoring results based on the green peak height and Cu2+ content in corn leaves

3.3 多变量的IWD-Hankel-SVD模型应用

植被中的叶绿素含量是植物生长健康状况的重要生物指标，决定着绿色植物的光合作用效率。随着Cu2+胁迫浓度的增大，玉米叶片中的Cu2+含量在逐渐增大，而叶绿浓度相对值在逐渐减少。这是因为重金属胁迫对植物的生长具有负向影响，抑制植物进行光合作用，造成叶绿素水解增加，含量降低，改变细胞结构。植物叶片中叶绿素含量的变化会引起可见光和近红外波段的反射率变化[15]。选择数据采集时测得的叶绿素浓度相对值、特征波段1 266～1 429和1 836～1 913 nm分别计算得到的模型奇异熵作为参数，结合偏最小二乘回归分析预测玉米叶片受重金属污染下的Cu2+含量。经验证可得预测值与实测值的相关性更强，如图11所示，其决定系数R2为0.947 6。与单变量的模型应用结果相比，多变量的IWD-Hankel-SVD模型具有较强的鲁棒性、稳健性。

图11 玉米叶片中Cu2+含量的多变量模型预测精度Fig.11 Accuracy of the multivariable model on predicting Cu2+ content in corn leaves

4 结 论

重金属Cu2+胁迫引起的玉米叶片光谱特征变化难以依据光谱曲线形态甄别其污染程度。提出了IWD概念和方法，并通过构造IWD-Hankel-SVD模型对玉米叶片高光谱数据的Cu污染信息进行反演，甄别和预测玉米叶片受Cu污染的程度，认为：(1)对Cu污染的玉米叶片光谱数据进行IWD处理，分解提取到的固有旋转分量PRC1能较好地保留了原始输入光谱信息，有效避免了模态混叠现象；(2)通过特征波段1 266～1 429和1 836～1 913 nm求得的模型奇异熵与玉米叶片中的Cu2+含量存在很强的相关性，模型预测的决定系数R2均高达0.9左右。通过一些传统预测方法和由SPA选择的特征波段的预测结果比较，验证了依据所提取特征波段这单一变量形式的IWD-Hankel-SVD模型在玉米叶片Cu污染信息和污染程度方面具有一定的优越性和较好的解释能力；(3)选择玉米叶片中叶绿素浓度相对值以及叶片光谱特征波段1 266～1 429和1 836～1 913 nm分别计算得到的模型奇异熵作为参数，形成多变量的模型预测技术，可得到叶片中Cu含量的预测值与实测值之间相关性更强，其决定系数R2达到0.9476。与单变量模型相比，多变量模型具有较强的鲁棒性和稳健性。