铜胁迫下玉米叶片的HHT包络谱变化与污染预测模型

杨可明 程 龙 郭 辉 王 敏

(1.中国矿业大学(北京)煤炭资源与安全开采国家重点实验室，北京 100083;2.中国矿业大学(北京)地球科学与测绘工程学院，北京 100083;3.安徽理工大学测绘学院，淮南 232001; 4.华北理工大学, 唐山 063210)

0 引言

农作物重金属污染引发的粮食安全问题备受各界关注[1-2]。对于植物而言，Cu是植物生长发育的必需微量元素，其在常量时对植被有益，如果在植物体内积累过多，会对植物产生Cu毒害作用，使植物体内的代谢紊乱，直接影响植物生长发育，乃至造成植株死亡[3-4]；更严重的是农作物受到Cu等重金属污染，不但影响农产品的品质，还可能通过食物链进入人体，对人类的生存和健康构成威胁[5-6]。已有研究表明，当植物受到重金属Cu胁迫时，Cu2+进入植物体叶片内在细胞壁中积累，并通过离子泵进入细胞，代替一些二价离子参与细胞代谢，从而影响细胞色素的形成，改变细胞的渗透压，最终造成细胞内部结构发生变化[7-8]。而且，植物叶片细胞内色素含量的变化将影响叶片对光的吸收与反射[9]，改变光在植物体内反射和散射的路径，使植被光谱特性发生变化[10]。通常，正常生长的绿色植物具有典型的反射光谱曲线，但是，受重金属污染后的植被光谱会有畸变发生，且这种畸变信息异常微弱，因此，在可见光和近红外波谱区间受污染与未受污染的植被光谱曲线形态相似度很高，但可利用具有数据量丰富、光谱分辨率高、对地面植被无破坏、能够实时动态监测等优势[11-12]的高光谱遥感技术实现重金属污染的光谱变异特征提取和污染监测。

近年来，高光谱遥感技术在农作物重金属污染检测中得到广泛应用。顾艳文等[13]研究了小白菜在不同镉胁迫梯度下叶片光谱反射特征，通过相关分析和逐步回归统计等方法发现基于红边面积的倒数模型能够较好地反演小白菜镉含量。王慧等[14]利用光谱分析方法研究不同浓度Cu、Zn胁迫下的小麦冠层光谱变化情况，发现重金属胁迫下的冠层光谱在可见光范围内反射率增加，而在近红外部分的反射率下降，“红边”和“红谷”均出现蓝移现象。关丽等[15]基于镉胁迫下水稻的生理生态表征构建光谱指数的二维模型，揭示了在不同浓度的镉处理下叶绿素、水分和细胞结构等响应因子在各光谱指数空间的分布规律。李婷等[16]运用BP人工神经网络模型构建水稻敏感光谱指数与叶绿素含量的关系，建立了水稻重金属污染胁迫的光谱分析模型。以上方法通过对光谱数据进行一阶微分、包络线去除等变换，可以增强光谱数据的相关性，筛选出敏感的特征波段和红边参数、红边“蓝移”指数等相关参数，以建立模型等进行污染监测和含量预测，但均是在光谱域和空间域中分析而并未从频率域中分析高光谱数据。因此，本文引入希尔伯特-黄变换(Hilbert-Huang transform,HHT)用于探索光谱信息特征与植被的重金属污染监测。基于重金属Cu2+胁迫下玉米盆栽实验及玉米叶片光谱测量等数据，利用HHT方法获取不同Cu胁迫梯度下光谱的包络谱，在比较分析所获得的Hilbert包络谱基础上，构建用于描述不同Cu污染程度的包络谱峰值指标(E1)和包络谱峭度系数(E2)两个谱特征参量，并基于HHT包络谱特征参量建立玉米叶片中Cu2+含量的预测模型，以实现对玉米叶片的Cu2+污染检测和污染程度判断。

1 理论与方法

HHT是由HUANG等[17]提出的新型信号时频联合分析方法，由经验模态分解(Empirical mode decomposition, EMD)方法和Hilbert变换组成，该方法能够在消除人为因素的情况下根据信号特点自适应分解非线性非平稳信号，得到分辨率高的谱结构[18-20]，适合于植物重金属污染的光谱畸变信息探测，主要利用HHT方法获取植被光谱在频域上分布的Hilbert包络谱(Envelope spectrum)[21-22]，通过研究包络谱的性质和规律，探索植物受重金属Cu胁迫下玉米叶片光谱响应特征。

1.1 经验模态分解

EMD可以把任意一个复杂的非平稳信号自适应地分解为有限个窄带信号，称为本征模态函数(Intrinsic mode function, IMF)，其中任何一个IMF都必须满足以下2个条件：整个信号中极值点数与零点数相等或最多相差1个；由极小值和极大值点确定的上下包络线均值为零。

对任意给定的信号进行经验模态分解的步骤如下：

(1)找出原始信号x(t)中所有的极大值点和极小值点，采用三次样条函数拟合原始数据数列的上下包络线。

(2)计算上下包络线的均值，记为u1(t)；并将原始数据序列x(t)减去该均值可得到一个去掉低频的新数据序列

y1(t)=x(t)-u1(t)

(1)

(3)判断y1(t)是否满足IMF条件，否则对y1(t)信号重复上述两过程，直到均值趋近于零为止，这样可得到第1个IMF分量c1(t)，它代表信号x(t)中最高频率的分量。

(4)将原始信号x(t)与第1个IMF分量c1(t)相减，得到原始信号中不包含最高频率分量的剩余信号

r1(t)=x(t)-c1(t)

(2)

(5)将r1(t)作为原始信号重复步骤(1)～(3)，得到其它的IMF分量ci(t)(i=2，3，…，n)，直到残余函数rn(t)为单调函数为止。

则原始信号可以表示为

(3)

各个分量包含了信号从高到低不同频率段的成分，每个频率段包含的频率分辨率都随信号本身变化，具有自适应多分辨分析特性。

1.2 Hilbert变换

通过对信号进行经验模态分解得到IMF后，就可以对每一个IMF作Hilbert变换得到有意义的瞬时频率，就可精确表达频率随时间的变化。

(4)

式中，H(ci(t))为第i阶IMF分量ci(t)的Hilbert变换函数。当ci(t)和H(ci(t))为共轭复数时，可构造解析信号

Zi(t)=ci(t)+jH(ci(t))=ai(t)ejφi(t)

(5)

(6)

(7)

式中ai(t)——解析信号的幅度或能量

φi(t)——解析信号的相位

由此，通过求出分量的瞬时频率wci(t)，可表达第i阶IMF分量ci(t)随时间的变化。

(8)

(9)

式中aci(t)——分量的幅度或能量

则原始信号x(t)可以表示为

(10)

式中n——IMF分量的个数

对解析信号Zi(t)求模即得到信号的包络

(11)

对式(10)的包络信号进行谱分析就可得到Hilbert包络谱。

1.3 HHT包络谱特征参量

为了定量地描述HHT包络谱的特征，参考文献[23]，定义HHT包络谱特征参量为：

(1)包络谱峰值指标(E1)：是指包络谱幅值在0～1 200 Hz频率内的最大值与其均方根之比，表达式为

(12)

式中 Max(xf)——包络谱信号x在频率f内峰值

ψf——包络谱信号x的均方根

(2)包络谱峭度系数(E2)：是归一化的四阶中心矩，可以用于反映包络谱信号的分布特性，表达式为

(13)

式中N——包络谱信号的长度

μ——包络谱信号x的均值

σ——信号x的标准差

xi——第i个频率对应幅值

1.4 精度评价指标

采用模型决定系数(R2)、标准差(S)和均方根误差(RMSE)对预测模型进行精度评定[24]，计算公式为

(14)

(15)

(16)

式中n——样本数量

yi、yj——实测值和预测值

2 实验设计与数据采集

2.1 实验设计

选用有底漏的花盆进行“蜜糯1号”玉米种子培植，用CuSO4·5H2O分析纯胁迫玉米发育生长。为保证玉米生长环境的统一性从播种到处理均在温室培养，2017年5月10日进行玉米种子催芽并于5月12日选取长势相同的玉米幼苗种植在盆栽土壤中。出苗后浇灌等量NH4NO3、KH2PO4和KNO3营养液。盆栽土壤分别设置了0(空白对照实验)、50、100、150、200、300、400、600、800、1 000、1 200 μg/g(以纯Cu2+计)共11种Cu2+胁迫梯度，每个浓度设置3组平行实验，共33盆盆栽，分别记为Cu(CK)-1、Cu(CK)-2、Cu(CK)-3，Cu(50)-1、Cu(50)-2、Cu(50)-3，…，Cu(1200)-1、Cu(1200)-2、Cu(1200)-3。在培植期定期浇水并保持每天通风换气；除了土壤中Cu2+胁迫梯度不同，所有盆栽的培养均在温室同一条件下，统一水肥管理。

2.2 光谱数据采集

采用350～2 500 nm波谱范围的SVCHR-1024I型高性能地物光谱仪测定“蜜糯1号”玉米叶片在不同Cu2+胁迫梯度下光谱数据。2017年7月19日在密闭室内进行玉米叶片的光谱采集，采集时将玉米叶片裁剪，平铺在方桌上测定，使用光谱仪配套功率为50 W的卤素灯光源和4°视场角的探头，探头垂直于叶片表面50 cm；所采集的光谱反射系数经专用平面白板进行标准化处理。为了防止其他物体对玉米叶片光谱的影响，全过程均用黑色塑料布盖住方桌。从而得到不同胁迫梯度下的玉米叶片光谱数据，由于Cu2+胁迫浓度过高，导致1 000、1 200 μg/g胁迫梯度下的盆栽玉米叶片在幼苗期开始渐渐枯萎，无法用于后续数据分析，200、300 μg/g植株长势较不理想，且200、300 μg/g胁迫梯度下的测定数据有异常，故不予采用。因此，本文选取0、100、150、400、600、800 μg/g为研究对象，其中Cu(CK)-1、Cu(100)-1、Cu(150)-1、Cu(400)-1、Cu(600)-1、Cu(800)-1记为第1组样本，Cu(CK)-2、Cu(100)-2、Cu(150)-2、Cu(400)-2、Cu(600)-2、Cu(800)-2记为第2组样本，Cu(CK)-3、Cu(100)-3、Cu(150)-3、Cu(400)-3、Cu(600)-3、Cu(800)-3记为第3组样本，所测得的光谱数据如图1所示，其中Cu(CK)可以近似认为是玉米叶片的无胁迫污染对照光谱。

图1 不同浓度Cu2+胁迫下玉米叶片的光谱采集数据Fig.1 Collected spectral data of corn leaves stressed by Cu2+ with different contents

2.3 Cu2+含量测定

将采集光谱的玉米叶片样品进行冲洗、干燥、粉碎等一系列预处理后，贴上标签，然后将样品袋带入实验室化验分析。在相同的规范条件下，经微波处理后，采用电感耦合等离子发射光谱仪(ICP-OES)测定各胁迫梯度下玉米叶片样品中的Cu2+含量，如表1所示。其中第1组和第2组样本的Cu2+含量数据用于建立预测模型，第3组样本数据用于预测模型验证。

表1 土壤中铜胁迫与玉米叶片中Cu2+对照表Tab.1 Comparison of copper stress gradients in soil and Cu2+ contents in corn leaves

3 HHT包络谱分析

3.1 玉米叶片光谱HHT分析

由信号理论可知，玉米叶片光谱为“非线性非稳定”信号，不能满足Hilbert变换的要求，因此，需要对叶片光谱信号进行EMD自适应地分解得到IMF后再做Hilbert变换和包络谱分析。随机选择无污染的玉米叶片光谱中的一组进行HHT分析，得到8个IMF分量和对应每个IMF及原始光谱的包络谱图如图2所示。

图2 玉米叶片光谱平行组HHT时频分析结果Fig.2 Results of HHT time frequency analysis on parallel group spectral data of corn leaves

3.2 不同Cu2+胁迫梯度下玉米叶片光谱HHT分析

玉米受到Cu2+污染必然会引起玉米叶片生理生态参数和光谱波形的变化。由于Cu2+胁迫梯度不同，这些参数的变化情况也有所差异，有时可能会出现某些参数变化过小而在光谱数据上没有产生显著响应，只是在局部波段发生微小的突变，从而导致了光谱的畸变性。因此，为了使玉米铜污染的诊断更加可靠和有效，一般情况下不采用原始光谱信号作为判别依据，而是利用原始光谱信号经HHT后得到的包络谱曲线进行诊断分析。通过研究包络谱的性质和规律，揭示玉米受重金属Cu污染的光谱特征信息。

实验对6个梯度玉米光谱进行HHT时频分析后，不同Cu2+胁迫梯度下玉米叶片光谱获得的各阶IMF数量不同，存在一定的差异，根据各阶IMF与原始信号的相关性，得到IMF4与原始信号的相关性最大，因而选取IMF4进行HHT分析。图3为第1组样本不同Cu2+胁迫梯度下玉米叶片光谱HHT的IMF4及包络谱，由IMF4的包络谱可以大致看出玉米受到不同程度Cu2+污染下的包络谱曲线在各个频段的差异且为连续谱，包络谱的频率主要分布在100 Hz以内，由IMF4前100 Hz的包络谱曲线中可以更加直观体现出不同Cu2+胁迫梯度下包络谱的变化信息，可以在辨别光谱整体形状相似性的前提下，增强对光谱局部特征差异性的分辨能力。虽然包络谱能较好地体现玉米受污染的变化信息，可以作为识别污染的一个重要特征，但也可以看出，并不能直接用于判别污染。因此需在包络谱的基础上，通过计算E1和E2两个包络谱特征参量值识别玉米叶片的Cu2+污染程度以及预测叶片中的Cu2+含量。

图3 不同Cu2+胁迫梯度下玉米叶片光谱的IMF4分析结果Fig.3 IMF4 analysis results on spectral of corn leaves with different Cu2+ stress gradients

3.3 包络谱特征参量的计算与分析

实验分别计算得到5个Cu2+胁迫梯度及对照组Cu(CK)-1、Cu(CK)-2共计12个样本IMF4的HHT包络谱特征参量E1和E2值及其平均值E1a、E2a。由表2可知，当Cu2+胁迫梯度为150 μg/g时，E1和E2值相比较Cu(CK)、Cu(100)有明显的增大趋势，且E1和E2值均达到最大值；而当Cu2+胁迫梯度高于150 μg/g后，E1和E2值停止增大开始呈现下降趋势，表现出Cu2+对玉米的生长具有“低促高抑”的效应，即：随着Cu2+浓度增加，玉米在低于150 μg/g浓度时仍可正常生长，而在高于浓度150 μg/g后生长受到抑制。

表2 不同Cu2+胁迫梯度下玉米叶片光谱平行组包络谱特征参量值及平均值Tab.2 Spectral parallel group envelope spectral parameters and mean value of corn leaves with different Cu2+ stress gradients

3.4 玉米叶片Cu2+污染的预测模型建立

通过对各胁迫梯度下光谱的HHT包络谱特征参量进行统计计算，发现其结果与玉米叶片中Cu2+含量变化规律相似。如图4所示，随机选取第1组样本数据处理得到E1、E2与叶片中Cu2+含量的相关变化图。从图中可以看出，E1、E2值与叶片中Cu2+含量呈现相同的变化趋势，即都与Cu2+含量呈正相关关系且E1和E2与叶片中Cu2+含量都具有良好的相关性，相关系数分别为0.58和0.76，因此可认为E1、E2都能在一定程度上区分玉米叶片的Cu2+污染程度并可预测其Cu2+含量，其中E2效果最优。

图4 不同胁迫梯度下包络谱特征参量值和叶片中Cu2+含量变化的相关性分析Fig.4 Correlation analysis of variations between envelope spectral characteristic parameter values and Cu2+ content in leaves with different Cu2+ stress gradients

根据以上相关性分析，以实验中的第1组和第2组不同Cu2+胁迫梯度下玉米叶片中的Cu2+含量为因变量，在E1、E2中选取单因素和双因素变量进行回归建模分析(表3)。依据预测模型决定系数R2最大和标准差S最小的原则，由表3可得，基于E1和E2双因素变量模型的R2最大，达到0.69(P<0.01),S最小，为6.992 2，具备最优的Cu2+含量预测潜力。E1模型的R2值最小，预测玉米叶片中Cu2+含量的能力相对较弱。

为了检验基于E1和E2单、双因素变量所建模型的准确性和可靠性，选用实验所测第3组不同胁迫梯度下玉米叶片光谱以及叶片中Cu2+含量数据进行模型验证，通过计算检验光谱的E1、E2包络谱特征参量值应用于相应的预测模型，其模型预测结果与实测的叶片中Cu2+含量线性拟合效果如图5所示。由图5可知，基于E1、E2双因素变量建立的模型预测结果最优，R达到了0.902 2，RMSE为4.348 3 μg/g，表明该模型预测的叶片Cu2+含量数值最接近真实值。E2模型预测能力较优，R为0.889 8，RMSE为6.967 0 μg/g，该模型也可用于预测玉米叶片中的Cu2+含量。E1的模型预测效果不理想。

表3 玉米叶片中Cu2+含量预测的包络谱特征参量模型Tab.3 Models on predicting Cu2+contents in corn leaves based on envelope spectral characteristic parameters

图5 玉米叶片中Cu2+含量实测值与预测值的拟合曲线Fig.5 Fitting curves between measured and predicted values of Cu2+ contents in corn leaves

4 结束语

采用HHT方法研究了不同Cu2+胁迫梯度下玉米叶片包络谱的变化，并计算了包络谱峰值指标(E1)和包络谱峭度系数(E2)等包络谱特征参量值。结果表明，玉米叶片光谱的包络谱是分布在100 Hz以内的连续谱；在不同Cu2+胁迫梯度下，E1和E2均与叶片中Cu2+含量呈现出相同的变化趋势，两个包络谱特征参量值均在Cu(150)位置达到最大值，体现出了不同浓度的Cu2+对玉米的生长有“低促高抑”的效果。同时对比分析得出包络谱特征参量值与叶片中Cu2+含量具有良好的相关性，因此，可把E1和E2作为玉米Cu污染与判断其污染程度预测指标。在构建包络谱特征参量E1、E2单因素和双因素变量的Cu2+含量预测模型的基础上，通过各种模型的应用比较和评价，得出基于E1、E2的双因素变量模型在预测玉米叶片中Cu2+含量时敏感性更优，其模型预测值与实测值的相关系数R达0.902 2，RMSE最小，为4.348 3 μg/g。