铜胁迫下玉米光谱变化的奇异性诊断指数与污染甄别

李 燕，杨可明*，王 敏，2，程 凤，高 鹏，张 超

（1.中国矿业大学（北京）地球科学与测绘工程学院，北京 100083；2.华北理工大学，河北 唐山 063210）

近年来环境污染问题日趋严重，其中土壤中铜（Cu）、铅（Pb）、镉（Cd）、铬（Cr）、汞（Hg）和砷（As）等重金属污染较为突出，植被吸收土壤中高标重金属后，会在体内产生富集效应，当重金属含量达到一定限度后，会造成植被过氧化胁迫、影响细胞酶活性、破坏细胞构造，进而影响植被正常代谢[1]，甚至危害人类健康[2]，因此重金属污染及其实时有效监测越来越受到社会关注。在重金属污染监测方面，传统的生化等监测方法费时费力且不能大面积应用，而高光谱遥感因其丰富的光谱信息、无破坏和实时动态大面积监测等优势弥补了传统方法的不足[3]，已成为当今重金属污染遥感监测的一个重要研究内容，其中，农作物重金属污染更是其研究的热点之一。受重金属污染的植被反射光谱特性会发生改变[4-6]，因而植被光谱变化的奇异信息就成为了重金属污染监测的重要依据。任艳红等[7]、王圆圆等[8]研究发现水稻、盐肤木等作物的光谱特征与Pb、Zn等含量之间有显著相关性；杜华强等[9]验证植物光谱曲线的分形测量可用于诊断植被健康状况且分形维数能客观地反映植株生长状态；Noh等[10]建立神经网络模型提取光谱变化特征对植被胁迫状况进行分析。上述研究都是直接利用植被光谱曲线特征进行污染监测分析，但由于植被生化组分的复杂性，其微弱量变不会引起光谱曲线的显著变化，此时直接分析光谱特征诊断植被健康状态达不到最佳效果，因此寻找一种更为灵敏的监测方法是植被重金属遥感监测领域的热点。

近年来，具有多分辨率特点的小波变换（Wavelet transform，WT）技术被引进到高光谱领域，在高光谱遥感植被重金属污染监测和有机质含量反演研究中具有很好的应用效果[11-13]。WT具有时频局部特性，可应用于信号细节特征的提取，但由于植被光谱反射率变化尺度分布在一个较大范围内，并且WT分解过程不具有自适应性以及能量泄露等局限性，提取光谱突变信息并不能达到很好效果。经验模态分解（Em⁃pirical mode decomposition，EMD）是近年来出现的一种新的时频分析方法，该方法从信号自身出发，在保留原始信号特征的同时具有很高分辨率，可以清晰地提取出信号的细微变化，在分析处理非平稳非线性信号方面优势明显，目前在超声信号去噪[14]、机械振动故障[15]和地震信号分析[16]等领域得到广泛应用，而在高光谱遥感分析处理领域很少涉及。针对这些问题，本文将结合EMD与WT构建时频分布特征提取的EMD-WT模型，对受不同程度污染的玉米光谱信息奇异特征进行提取，并结合欧氏距离（Euclidean dis⁃tance，ED）、光谱角（Spectral angle，SA）等构建奇异性诊断指数（SI），对玉米奇异性变化进行定性分析，从而实现玉米受重金属铜胁迫程度的甄别。同时，通过与常规监测方法所得结果进行比较分析，验证本文所应用方法的有效性与优越性，为植被重金属污染监测方面的研究提供参考依据。

1 材料与方法

1.1 胁迫实验与数据采集

1.1.1 玉米培育

设置胁迫梯度为 0（CK）、100、300、500 μg·g-1的CuSO4·5H2O溶液，翻土将其加入到玉米实验盆钵中，每梯度设置3组平行试验，共12盆。2016年5月6日对玉米种子进行催芽处理，出苗后向盆栽中添加NH4NO3、KH2PO4和KNO3营养液。玉米培育期间定期通风与浇水，保持适宜的培育温度与湿度。

1.1.2 光谱数据采集

2016年7月17日对玉米叶片反射光谱进行测量。玉米光谱采集在50 W卤素灯光源照射下进行，将光谱范围为350～2500 nm的SVC HR-1024I型地物光谱仪探头视场角设置为4°并垂直于叶片25 cm处进行光谱采集，采集的光谱使用平面白板进行标准化。选取每盆玉米植株的老、中、新3种代表性叶片进行光谱测试，获得3组光谱数据（老、中、新），计算3组数据平均值作为各盆玉米光谱值。最后分别计算各胁迫梯度下3组平行试验玉米光谱数据均值作为各梯度的光谱值。

1.1.3 玉米叶片中Cu2+含量测定

2016年9月16日对采集过光谱数据的叶片进行冲洗、烘干、粉碎等样品预处理，再经高纯硝酸和高氯酸消化处理后用WFX-120型原子吸收分光光度计进行叶片中Cu2+含量的测定，每梯度测量3次后取平均值作为该梯度叶片中的Cu2+含量，测量结果如表1。叶片中Cu2+含量随胁迫浓度的增大表现为先增加后降低的趋势，根据Cu（500）的3盆玉米植株长势（出现枯黄叶片）推断，这可能是由于向土壤中施加较高浓度的Cu2+时，植株根部受到较严重损伤，对Cu2+吸收降低所致。

1.2 经验模态分解

经验模态分解（EMD）法是由美国学者N.E.Huang于1996年提出的一种时频分析方法，该方法不基于任何数学函数，而是根据信号自身的特点，将其分解为不同时间尺度特征的本征模函数（Intrinsic mode function，IMF）分量和残差余项r（t）。对给定的一段信号序列V（x），EMD过程为[17]：

（1）计算包络线均值M（t）。找出信号所有极大、极小值点后采用3次样条函数拟合信号的上、下包络线，计算上、下包络线的均值M（t）；

表1 不同胁迫浓度下玉米叶片中Cu2+含量（μg·g-1）Table 1 The Cu2+contents in corn leaves under different stress concentration（μg·g-1）

（2）获取IMF1分量。原信息V（x）减去M（t）得到新信息序列H（x），将H（x）视为新的原始信号V（x），若还存在负的局部极大值和正的局部极小值，重复以上步骤，直到获取基本IMF分量，记为IMF1；

（3）获取IMFn及r（t）。把IMF1分量从原始信号中分离出来，得到新序列Q（x），将其视为新的V（x），重复步骤（1）、（2）直到分解出第n个IMF分量及残差r（t）为止。其中，r（t）代表了信号的平均趋势，IMF1～IMFn包含原始信号不同时间尺度的信息。

从以上分解过程可知，EMD是一种基于极值点的特征时间尺度的分解方法，反映信号随时间变化的局部特征[18]，在分解过程中首先把最主要的信息分解出来，其中IMF1包含了信号的最主要能量。

1.3 小波分析

小波变换（WT）是1982年法国物理学家Morlet提出的一种信号变换方法。信号f（t）的小波变换定义为：

式中：a和b分别为尺度因子和平移因子，a和b的变换可以实现时间窗口的伸缩与平移；φ（t）为母小波，φa，b（t）为基函数，该基函数由a、b变换得到

分解过程中，基函数沿时间轴移动去逼近或表示原始信号，通过基函数的伸缩平移，在高频处有较高时间分辨率，低频处有较高频率分辨率。该过程实现了信号时频局部特征提取，能够用于探测信号中的突变信息。但由于小波分析过程中只能使用一个基函数，该基函数对整条光谱并不具有普适性。此外由于小波基长度有限，窗口内的能量会泄漏到其他波段，其他波段能量也会渗透到函数窗内[19]，造成分解结果中存在较大冗余，对分辨率造成影响，因此，采用EMD来改善其中的不足。

1.4 EMD-WT的时频分布特征提取模型

EMD-WT模型的主要思想是：对原始信号进行EMD处理，提取IMF1分量后对其进行小波分解，IMF1分量包含信号的高频率信息，信号的突变信息包含在该成分中，对IMF1分量进行WT的局部特征分析，赋予IMF1分量具体的物理意义，降低噪音干扰，从而得到信号精确的时频分布特征。

采用Daubechies小波系列中的“Db5”小波函数（图1）对信号进行分解，“Db5”小波具有正交性、双正交性、非对称性等特点[12]，而且时域长度较长，具有较高时间分辨率，对光谱奇异特征有较好的探测效果；小波系数表示小波基函数与信号序列的逼近程度，利用小波系数曲线来描述信号序列的突变特征，如图2所示，当信号频率发生突变时，小波函数与信号的逼近程度增加，所得小波系数为一个非零值x；在信号频率平稳处，所得小波系数为0，因而根据小波系数模值大小能够准确判别信号序列频率的突变程度。

1.5 光谱奇异性诊断指数

图1 Db5小波函数Figure 1 Db5 wavelet function

图2 小波变换拟合突变信号原理Figure 2 The principle of wavelet transform fitting the mutant signal

对原始光谱信号的IMF1分量进行Db5小波变换，提取第五层小波系数，以CK（0）小波系数曲线作为基准曲线，分别计算100、300、500 μg·g-1胁迫浓度下的小波系数曲线与基准曲线之间的光谱角（SA）及欧氏距离（ED）[20-22]。由于欧式距离可以表征两条相似曲线在变化幅度上的差别，光谱角可以表征两条曲线整体相似程度，两者结合构造奇异性诊断指数（SI），可以有效地反映光谱信号奇异性特征的变化，以此来表征玉米受污染程度。SI可表示为：

式中：ED、SA分别为各胁迫浓度系数曲线与CK（0）组系数曲线之间的欧式距离及光谱角。

2 结果与讨论

2.1 玉米光谱奇异信息提取

由于1300～2500 nm波长范围的玉米反射率受水分和大气的影响很大，会对结果的分析产生干扰，因此本文选取350～1300 nm之间的玉米光谱进行研究。图3是采集的12盆玉米光谱原始数据，分析发现，不同浓度胁迫的玉米光谱曲线间并没有显著差异。

对Cu（300）光谱进行EMD处理，提取IMF1分量，所得结果如图4所示。从图4和表2分析可知，IMF1分量的极值点处，受铜胁迫的玉米光谱反射率增量[ΔCu（300）]比未受铜胁迫的玉米光谱反射率增量[ΔCK（0）]大，表明受铜污染后的光谱反射率发生突变；且IMF1极值点模值越大，两个胁迫浓度的光谱反射率增量差值（Δ）越大，光谱频率突变越大，该点奇异性越强，表征玉米受铜胁迫越严重。因此根据光谱信号IMF1分量极值点与原始光谱奇异点的对应关系，可以对信号奇异点出现的位置做定性分析。对IMF1分量的极值点及相应光谱突变点进行统计分析，表2结果显示光谱曲线的突变主要集中在590～720 nm波段。

基于EMD提取所得的各胁迫浓度光谱IMF1分量，对其分别进行Db5小波5层分解，得到1～5层的小波系数，相应结果如图5。由图5（a）可知，各胁迫浓度光谱奇异变化信息主要在590～720 nm波段，在400 nm和1000 nm波段附近其出现大量的噪声信息。由图5（b）～图5（f）可以看出，随着小波变换分解尺度的增加，曲线噪声模值越来越小，奇异信息模值将越来越大，其中d5小波系数能够很好地反映奇异波段的幅值与位置信息，且能够很好地消除噪声信息的干扰。

综上所述，对原始光谱进行EMD处理，IMF1分量能有效地监测到奇异波段，再利用“Db5”小波的d5尺度对IMF1分量进行分解，其中第5层小波系数能有效地消除噪声，突出奇异信息，进而实现不同Cu2+胁迫梯度下玉米光谱奇异性信息的提取。

图4 光谱信号突变点与IMF1分量Figure 4 Abrupt points of spectral signal and IMF1 components

表2 光谱信号突变点与IMF1值极值点计算结果Table 2 Statistical result of abrupt points of spectral signal and extreme values of IMF1

2.2 不同Cu2+胁迫浓度下玉米光谱奇异性分析

不同Cu2+胁迫浓度下玉米光谱IMF1分量的d5小波系数曲线如图6所示，分析发现，随着玉米叶片中Cu2+含量的增加，曲线极值点的幅度与数目不断增加，说明玉米光谱奇异性不断增强。不同胁迫浓度下玉米叶片中Cu2+含量及不同Cu2+胁迫浓度相对于对照组CK（0）的奇异性变化计算结果如表3，奇异性诊断指数随着叶片中重金属Cu含量的增加而增加，其相关系数达到0.972 4。

图5 Db5小波函数分解的不同尺度小波系数Figure 5 Wavelet coefficients from different scales using Db5 wavelet function

2.3 方法应用结果比较与拟合分析

为了验证奇异性诊断指数在诊断玉米叶片光谱变异信息与污染探测方面的有效性和优越性，将其与常规的绿峰高度（GH）、红边最大值（MR）、红边一阶微分包围面积（FAR）3个参数进行污染监测应用结果比较分析，每个参数的定义及计算方法如表4所示。

不同Cu2+胁迫浓度下叶片光谱的GH、MR、FAR、SI及其与叶片中Cu2+含量的相关系数计算结果如表5所示。分析发现，SI与玉米叶片中Cu2+含量的相关系数达到0.972 4，优于其他监测参数。同时将各参数与叶片中Cu2+含量进行线性拟合分析，拟合结果如图7，其中SI与叶片中Cu2+含量拟合程度较高，其拟合判定系数达到0.945 4，高于其他监测参数的判定系数。由此可知，光谱奇异性诊断指数具有更好的监测效果。

表3 不同Cu2+胁迫浓度下玉米叶片奇异性诊断指数Table 3 The singularity diagnostic index of corn leaves under Cu2+different stress concentration

图6 不同污染程度下玉米叶片光谱d5系数曲线Figure 6 Curves of d5 coefficient of corn leaves under different pollution levels

表4 玉米叶片光谱特征参数名称及定义Table 4 The name and definition of corn leaves′differential spectral characteristic parameters

3 结论

（1）EMD-WT时频分布特征提取模型能有效提取玉米光谱曲线铜胁迫弱信息。EMD分解的IMF1分量可准确地探测到光谱奇异范围，结果表明，铜胁迫下玉米光谱奇异信息主要在590～720 nm波段范围内。Db5小波函数对IMF1分量进行5层分解，其中第5层系数能有效地提取出玉米铜胁迫的奇异信息。

（2）奇异性诊断指数能有效地诊断玉米铜胁迫程度，SI与玉米叶片中重金属铜含量显著相关，SI越大，重金属Cu2+含量越高，其相关系数达到0.972 4。同时与GH、MR、FAR监测方法的应用结果进行比较分析，验证了奇异性诊断指数在玉米铜污染信息监测中具有更好的效果。

表5 玉米叶片Cu2+污染监测奇异性诊断指数与常规方法应用结果Table 5 Application result on the singularity diagnostic index and some convention methods for monitoring Cu2+pollution of corn leaves

图7 玉米叶片中Cu2+含量与各监测方法计算值拟合结果Figure 7 Fitting results on the computing values of the monitoring methods and the Cu2+contents of corn leaves

[2]Song B,Lei M,Chen T B,et al.Assessing the health risk of heavy met⁃als in vegetables to the general population in Beijing,China[J].Journal of Environmental Sciences,2009,21：1702-1709.

[3]Kloiber S M,Brezonik P L,Olmanson L G,et al.A procedure for region⁃al lake water clarity assessment using Landsat multispectral data[J].Re⁃mote Sinsing of Environment,2002,82（1）：38-47.

[4]朱叶青,屈永华,刘素红,等.重金属铜污染植被光谱响应特征研究[J].遥感学报,2014,18（2）：335-352.ZHU Ye-qing,QU Yong-hua,LIU Su-hong,et al.Spectral response of wheat and lettuce to copper pollution[J].Journal of Remote Sensing,2014,18（2）：335-352.

[5]李庆亭.重金属污染胁迫下盐肤木的生化效应及波谱特征[J].遥感学报,2008,12（2）：284-290.LI Qing-ting.Biogeochemistry responses and spectral characteristics of Rhus chinensis mill under heavy metal contamination stress[J].Journal of Remote Sensing,2008,12（2）：284-290.

[6]孙彤彤,杨可明,张 伟,等.玉米光谱小波奇异熵的铜污染监测研究[J].环境科学学报,2017,37（11）：4360-4365.SUN Tong-tong,YANG Ke-ming,ZHANG Wei,et al.Monitoring cop⁃per pollution based on wave singular entropy of corn leaves[J].Acta Sci⁃entiae Circumstantiae,2017,37（11）：4360-4365.

[7]任红艳,庄大方,潘剑君,等.重金属污染水稻的冠层反射光谱特征研究[J].光谱学与光谱分析,2010,30（2）：430-434.REN Hong-yan,ZHUANG Da-fang,PAN Jian-jun,et al.Study on can⁃opy spectral characteristics of paddy polluted by heavy metals[J].Spec⁃troscopy and Spectral Analysis,2010,30（2）：430-434.

[8]王圆圆,陈云浩,李 京,等.指示冬小麦条锈病严重度的两个新的红边参数[J].遥感学报,2007,11（6）：875-881.WANG Yuan-yuan,CHEN Yun-hao,LI Jing,et al.Two new red edge indices as indicators for stripe rust disease severity of winter wheat[J].Journal of Remote Sensing,2007,11（6）：875-881.

[9]杜华强,金 伟,葛宏立,等.用高光谱曲线分形维数分析植被健康状况[J].光谱学与光谱分析,2009,29（8）：2136-2140.DU Hua-qiang,JIN Wei,GE Hong-li,et al.Using fractal dimensions of hyperspectral curves to analyze the healthy status of vegetation[J].Spectroscopy and Spectral Analysis,2009,29（8）：2136-2140.

[10]Noh H,Zhang Q,Shin B,et al.A neural network model of maize crop nitrogen stress assessment for a multi-spectral imaging sensor[J].Bio⁃systems Engineering,2006,94（4）：477-485.

[11]刘美玲,刘湘南,李 婷,等.水稻锌污染胁迫的光谱奇异性分析[J].农业工程学报,2010,26（3）：191-197.LIU Mei-ling,LIU Xiang-nan,LI Ting,et al.Analysis of hyperspec⁃tral singularity of rice under Zn pollution stress[J].Transactions of the Chinese Society of Agricultural Engineering,2010,26（3）：191-197.

[12]于 雷,洪永胜,周 勇,等.连续小波变换高光谱数据的土壤有机质含量反演模型构建[J].光谱学与光谱分析,2016,36（5）：1428-1433.YU Lei,HONG Yong-sheng,ZHOU Yong,et al.Inversion of soil or⁃ganic matter content using hyperspectral data based on continuous wavelet transformation[J].Spectroscopy and Spectral Analysis,2016,36（5）：1428-1433.

[13]何汝艳,乔小军,蒋金豹,等.小波法反演条锈病胁迫下冬小麦冠层叶片全氮含量[J].农业工程学报,2015,31（2）：141-146.HE Ru-yan,QIAO Xiao-jun,JIANG Jin-bao,et al.Retrieving cano⁃py leaf total nitrogen content of winter wheat by continuous wavelet transform[J].Transactions of the Chinese Society of Agricultural Engi⁃neering,2015,31（2）：141-146.

[14]罗玉昆,罗诗途,罗飞路,等.激光超声信号去噪的经验模态分解实现及改进[J].光学精密工程,2013,21（2）：479-487.LUO Yu-kun,LUO Shi-tu,LUO Fei-lu,et al.Realization and im⁃provement of laser ultrasonic signal denoising based on empirical mode decomposition[J].Optics and Precision Engineering,2013,21（2）：479-487.

[15]李富香,张 韡.EMD技术在机械震动故障中的诊断方法[J].舰船科学技术,2016,38（20）：154-156.LI Fu-xiang,ZHANG Wei.The mechanical vibration fault diagnosis method based on EMD technology[J].Ship Science and Technology,2016,38（20）：154-156.

[16]张 猛,王华忠,隋志强,等.基于经验模态分解和小波变换的地震瞬时频率提取方法及应用[J].石油地球物理勘探,2016,51（3）：565-571.ZHANG Meng,WANG Hua-zhong,SUI Zhi-qiang,et al.Seismic in⁃stantaneous frequency extraction based on empirical mode decomposi⁃tion transform[J].Oil Geophysical Prospecting,2016,51（3）：565-571.

[17]Huang N E,Shen Z,Long S R,et al.The empirical mode decomposi⁃tion and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J].Proceedings：Mathematical,Physical and Engineer⁃ing Sciences,1998,454（1971）：903-995.

[18]李 卿.EMD时频分析方法的研究与应用[D].武汉：华中师范大学,2008.LI Qing.Research and application of time-frequency analysis method based on empirical mode decomposition[D].Wuhan：Central China Normal University,2008.

[19]龚志强,邹明玮,高新全,等.基于非线性时间序列分析经验模态分解和小波分解异同性的研究[J].物理学报,2005,54（8）：3947-3957.GONG Zhi-qiang,ZOU Ming-wei,GAO Xin-quan,et al.On the dif⁃ference between empirical mode decomposition and wavelet decompo⁃sition in the nonlinear time series[J].Acta Physica Senica,2005,54（8）：3947-3957.

[20]张浚哲,朱文泉,郑周涛,等.高光谱数据的相似性测度对比研究[J].测绘科学,2013,38（6）：33-36.ZHANG Jun-zhe,ZHU Wen-quan,ZHENG Zhou-tao,et al.Compar⁃ative study on similarity measures of hyperspectral remote sensing data[J].Science of Surveying and Mapping,2013,38（6）：33-36.

[21]闻兵工,冯伍法,刘 伟,等.基于光谱曲线整体相似性测度的匹配分类[J].测绘科学技术学报,2009,26（2）：128-131.WEN Bing-gong,FENG Wu-fa,LIU Wei,et al.Matching and classifi⁃cation based on the whole comparability measure of spectral curve[J].Journal of Geomatics Science and Technology,2009,26（2）：128-131

[22]张修宝,袁 艳,景娟娟,等.信息散度与梯度角正切相结合的光谱区分方法[J].光谱学与光谱分析,2011,31（3）：853-857.ZHANG Xiu-bao,YUAN Yan,JING Juan-juan,et al.Spectral dis⁃crimination method information divergence combined with gradient an⁃gle[J].Spectroscopy and Spectral Analysis,2011,31（3）：853-857.