銅脅迫下玉米光譜變化的奇異性診斷指數(shù)與污染甄別

李 燕，楊可明*，王 敏，2，程 鳳，高 鵬，張 超

（1.中國礦業(yè)大學（北京）地球科學與測繪工程學院，北京 100083；2.華北理工大學，河北 唐山 063210）

近年來環(huán)境污染問題日趨嚴重，其中土壤中銅（Cu）、鉛（Pb）、鎘（Cd）、鉻（Cr）、汞（Hg）和砷（As）等重金屬污染較為突出，植被吸收土壤中高標重金屬后，會在體內產生富集效應，當重金屬含量達到一定限度后，會造成植被過氧化脅迫、影響細胞酶活性、破壞細胞構造，進而影響植被正常代謝[1]，甚至危害人類健康[2]，因此重金屬污染及其實時有效監(jiān)測越來越受到社會關注。在重金屬污染監(jiān)測方面，傳統(tǒng)的生化等監(jiān)測方法費時費力且不能大面積應用，而高光譜遙感因其豐富的光譜信息、無破壞和實時動態(tài)大面積監(jiān)測等優(yōu)勢彌補了傳統(tǒng)方法的不足[3]，已成為當今重金屬污染遙感監(jiān)測的一個重要研究內容，其中，農作物重金屬污染更是其研究的熱點之一。受重金屬污染的植被反射光譜特性會發(fā)生改變[4-6]，因而植被光譜變化的奇異信息就成為了重金屬污染監(jiān)測的重要依據。任艷紅等[7]、王圓圓等[8]研究發(fā)現(xiàn)水稻、鹽膚木等作物的光譜特征與Pb、Zn等含量之間有顯著相關性；杜華強等[9]驗證植物光譜曲線的分形測量可用于診斷植被健康狀況且分形維數(shù)能客觀地反映植株生長狀態(tài)；Noh等[10]建立神經網絡模型提取光譜變化特征對植被脅迫狀況進行分析。上述研究都是直接利用植被光譜曲線特征進行污染監(jiān)測分析，但由于植被生化組分的復雜性，其微弱量變不會引起光譜曲線的顯著變化，此時直接分析光譜特征診斷植被健康狀態(tài)達不到最佳效果，因此尋找一種更為靈敏的監(jiān)測方法是植被重金屬遙感監(jiān)測領域的熱點。

近年來，具有多分辨率特點的小波變換（Wavelet transform，WT）技術被引進到高光譜領域，在高光譜遙感植被重金屬污染監(jiān)測和有機質含量反演研究中具有很好的應用效果[11-13]。WT具有時頻局部特性，可應用于信號細節(jié)特征的提取，但由于植被光譜反射率變化尺度分布在一個較大范圍內，并且WT分解過程不具有自適應性以及能量泄露等局限性，提取光譜突變信息并不能達到很好效果。經驗模態(tài)分解（Em?pirical mode decomposition，EMD）是近年來出現(xiàn)的一種新的時頻分析方法，該方法從信號自身出發(fā)，在保留原始信號特征的同時具有很高分辨率，可以清晰地提取出信號的細微變化，在分析處理非平穩(wěn)非線性信號方面優(yōu)勢明顯，目前在超聲信號去噪[14]、機械振動故障[15]和地震信號分析[16]等領域得到廣泛應用，而在高光譜遙感分析處理領域很少涉及。針對這些問題，本文將結合EMD與WT構建時頻分布特征提取的EMD-WT模型，對受不同程度污染的玉米光譜信息奇異特征進行提取，并結合歐氏距離（Euclidean dis?tance，ED）、光譜角（Spectral angle，SA）等構建奇異性診斷指數(shù)（SI），對玉米奇異性變化進行定性分析，從而實現(xiàn)玉米受重金屬銅脅迫程度的甄別。同時，通過與常規(guī)監(jiān)測方法所得結果進行比較分析，驗證本文所應用方法的有效性與優(yōu)越性，為植被重金屬污染監(jiān)測方面的研究提供參考依據。

1 材料與方法

1.1 脅迫實驗與數(shù)據采集

1.1.1 玉米培育

設置脅迫梯度為 0（CK）、100、300、500 μg·g-1的CuSO4·5H2O溶液，翻土將其加入到玉米實驗盆缽中，每梯度設置3組平行試驗，共12盆。2016年5月6日對玉米種子進行催芽處理，出苗后向盆栽中添加NH4NO3、KH2PO4和KNO3營養(yǎng)液。玉米培育期間定期通風與澆水，保持適宜的培育溫度與濕度。

1.1.2 光譜數(shù)據采集

2016年7月17日對玉米葉片反射光譜進行測量。玉米光譜采集在50 W鹵素燈光源照射下進行，將光譜范圍為350～2500 nm的SVC HR-1024I型地物光譜儀探頭視場角設置為4°并垂直于葉片25 cm處進行光譜采集，采集的光譜使用平面白板進行標準化。選取每盆玉米植株的老、中、新3種代表性葉片進行光譜測試，獲得3組光譜數(shù)據（老、中、新），計算3組數(shù)據平均值作為各盆玉米光譜值。最后分別計算各脅迫梯度下3組平行試驗玉米光譜數(shù)據均值作為各梯度的光譜值。

1.1.3 玉米葉片中Cu2+含量測定

2016年9月16日對采集過光譜數(shù)據的葉片進行沖洗、烘干、粉碎等樣品預處理，再經高純硝酸和高氯酸消化處理后用WFX-120型原子吸收分光光度計進行葉片中Cu2+含量的測定，每梯度測量3次后取平均值作為該梯度葉片中的Cu2+含量，測量結果如表1。葉片中Cu2+含量隨脅迫濃度的增大表現(xiàn)為先增加后降低的趨勢，根據Cu（500）的3盆玉米植株長勢（出現(xiàn)枯黃葉片）推斷，這可能是由于向土壤中施加較高濃度的Cu2+時，植株根部受到較嚴重損傷，對Cu2+吸收降低所致。

1.2 經驗模態(tài)分解

經驗模態(tài)分解（EMD）法是由美國學者N.E.Huang于1996年提出的一種時頻分析方法，該方法不基于任何數(shù)學函數(shù)，而是根據信號自身的特點，將其分解為不同時間尺度特征的本征模函數(shù)（Intrinsic mode function，IMF）分量和殘差余項r（t）。對給定的一段信號序列V（x），EMD過程為[17]：

（1）計算包絡線均值M（t）。找出信號所有極大、極小值點后采用3次樣條函數(shù)擬合信號的上、下包絡線，計算上、下包絡線的均值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的變換可以實現(xiàn)時間窗口的伸縮與平移；φ（t）為母小波，φa，b（t）為基函數(shù)，該基函數(shù)由a、b變換得到

分解過程中，基函數(shù)沿時間軸移動去逼近或表示原始信號，通過基函數(shù)的伸縮平移，在高頻處有較高時間分辨率，低頻處有較高頻率分辨率。該過程實現(xiàn)了信號時頻局部特征提取，能夠用于探測信號中的突變信息。但由于小波分析過程中只能使用一個基函數(shù)，該基函數(shù)對整條光譜并不具有普適性。此外由于小波基長度有限，窗口內的能量會泄漏到其他波段，其他波段能量也會滲透到函數(shù)窗內[19]，造成分解結果中存在較大冗余，對分辨率造成影響，因此，采用EMD來改善其中的不足。

1.4 EMD-WT的時頻分布特征提取模型

EMD-WT模型的主要思想是：對原始信號進行EMD處理，提取IMF1分量后對其進行小波分解，IMF1分量包含信號的高頻率信息，信號的突變信息包含在該成分中，對IMF1分量進行WT的局部特征分析，賦予IMF1分量具體的物理意義，降低噪音干擾，從而得到信號精確的時頻分布特征。

采用Daubechies小波系列中的“Db5”小波函數(shù)（圖1）對信號進行分解，“Db5”小波具有正交性、雙正交性、非對稱性等特點[12]，而且時域長度較長，具有較高時間分辨率，對光譜奇異特征有較好的探測效果；小波系數(shù)表示小波基函數(shù)與信號序列的逼近程度，利用小波系數(shù)曲線來描述信號序列的突變特征，如圖2所示，當信號頻率發(fā)生突變時，小波函數(shù)與信號的逼近程度增加，所得小波系數(shù)為一個非零值x；在信號頻率平穩(wěn)處，所得小波系數(shù)為0，因而根據小波系數(shù)模值大小能夠準確判別信號序列頻率的突變程度。

1.5 光譜奇異性診斷指數(shù)

圖1 Db5小波函數(shù)Figure 1 Db5 wavelet function

圖2 小波變換擬合突變信號原理Figure 2 The principle of wavelet transform fitting the mutant signal

對原始光譜信號的IMF1分量進行Db5小波變換，提取第五層小波系數(shù)，以CK（0）小波系數(shù)曲線作為基準曲線，分別計算100、300、500 μg·g-1脅迫濃度下的小波系數(shù)曲線與基準曲線之間的光譜角（SA）及歐氏距離（ED）[20-22]。由于歐式距離可以表征兩條相似曲線在變化幅度上的差別，光譜角可以表征兩條曲線整體相似程度，兩者結合構造奇異性診斷指數(shù)（SI），可以有效地反映光譜信號奇異性特征的變化，以此來表征玉米受污染程度。SI可表示為：

式中：ED、SA分別為各脅迫濃度系數(shù)曲線與CK（0）組系數(shù)曲線之間的歐式距離及光譜角。

2 結果與討論

2.1 玉米光譜奇異信息提取

由于1300～2500 nm波長范圍的玉米反射率受水分和大氣的影響很大，會對結果的分析產生干擾，因此本文選取350～1300 nm之間的玉米光譜進行研究。圖3是采集的12盆玉米光譜原始數(shù)據，分析發(fā)現(xiàn)，不同濃度脅迫的玉米光譜曲線間并沒有顯著差異。

對Cu（300）光譜進行EMD處理，提取IMF1分量，所得結果如圖4所示。從圖4和表2分析可知，IMF1分量的極值點處，受銅脅迫的玉米光譜反射率增量[ΔCu（300）]比未受銅脅迫的玉米光譜反射率增量[ΔCK（0）]大，表明受銅污染后的光譜反射率發(fā)生突變；且IMF1極值點模值越大，兩個脅迫濃度的光譜反射率增量差值（Δ）越大，光譜頻率突變越大，該點奇異性越強，表征玉米受銅脅迫越嚴重。因此根據光譜信號IMF1分量極值點與原始光譜奇異點的對應關系，可以對信號奇異點出現(xiàn)的位置做定性分析。對IMF1分量的極值點及相應光譜突變點進行統(tǒng)計分析，表2結果顯示光譜曲線的突變主要集中在590～720 nm波段。

基于EMD提取所得的各脅迫濃度光譜IMF1分量，對其分別進行Db5小波5層分解，得到1～5層的小波系數(shù)，相應結果如圖5。由圖5（a）可知，各脅迫濃度光譜奇異變化信息主要在590～720 nm波段，在400 nm和1000 nm波段附近其出現(xiàn)大量的噪聲信息。由圖5（b）～圖5（f）可以看出，隨著小波變換分解尺度的增加，曲線噪聲模值越來越小，奇異信息模值將越來越大，其中d5小波系數(shù)能夠很好地反映奇異波段的幅值與位置信息，且能夠很好地消除噪聲信息的干擾。

綜上所述，對原始光譜進行EMD處理，IMF1分量能有效地監(jiān)測到奇異波段，再利用“Db5”小波的d5尺度對IMF1分量進行分解，其中第5層小波系數(shù)能有效地消除噪聲，突出奇異信息，進而實現(xiàn)不同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小波系數(shù)曲線如圖6所示，分析發(fā)現(xiàn)，隨著玉米葉片中Cu2+含量的增加，曲線極值點的幅度與數(shù)目不斷增加，說明玉米光譜奇異性不斷增強。不同脅迫濃度下玉米葉片中Cu2+含量及不同Cu2+脅迫濃度相對于對照組CK（0）的奇異性變化計算結果如表3，奇異性診斷指數(shù)隨著葉片中重金屬Cu含量的增加而增加，其相關系數(shù)達到0.972 4。

圖5 Db5小波函數(shù)分解的不同尺度小波系數(shù)Figure 5 Wavelet coefficients from different scales using Db5 wavelet function

2.3 方法應用結果比較與擬合分析

為了驗證奇異性診斷指數(shù)在診斷玉米葉片光譜變異信息與污染探測方面的有效性和優(yōu)越性，將其與常規(guī)的綠峰高度（GH）、紅邊最大值（MR）、紅邊一階微分包圍面積（FAR）3個參數(shù)進行污染監(jiān)測應用結果比較分析，每個參數(shù)的定義及計算方法如表4所示。

不同Cu2+脅迫濃度下葉片光譜的GH、MR、FAR、SI及其與葉片中Cu2+含量的相關系數(shù)計算結果如表5所示。分析發(fā)現(xiàn)，SI與玉米葉片中Cu2+含量的相關系數(shù)達到0.972 4，優(yōu)于其他監(jiān)測參數(shù)。同時將各參數(shù)與葉片中Cu2+含量進行線性擬合分析，擬合結果如圖7，其中SI與葉片中Cu2+含量擬合程度較高，其擬合判定系數(shù)達到0.945 4，高于其他監(jiān)測參數(shù)的判定系數(shù)。由此可知，光譜奇異性診斷指數(shù)具有更好的監(jiān)測效果。

表3 不同Cu2+脅迫濃度下玉米葉片奇異性診斷指數(shù)Table 3 The singularity diagnostic index of corn leaves under Cu2+different stress concentration

圖6 不同污染程度下玉米葉片光譜d5系數(shù)曲線Figure 6 Curves of d5 coefficient of corn leaves under different pollution levels

表4 玉米葉片光譜特征參數(shù)名稱及定義Table 4 The name and definition of corn leaves′differential spectral characteristic parameters

3 結論

（1）EMD-WT時頻分布特征提取模型能有效提取玉米光譜曲線銅脅迫弱信息。EMD分解的IMF1分量可準確地探測到光譜奇異范圍，結果表明，銅脅迫下玉米光譜奇異信息主要在590～720 nm波段范圍內。Db5小波函數(shù)對IMF1分量進行5層分解，其中第5層系數(shù)能有效地提取出玉米銅脅迫的奇異信息。

（2）奇異性診斷指數(shù)能有效地診斷玉米銅脅迫程度，SI與玉米葉片中重金屬銅含量顯著相關，SI越大，重金屬Cu2+含量越高，其相關系數(shù)達到0.972 4。同時與GH、MR、FAR監(jiān)測方法的應用結果進行比較分析，驗證了奇異性診斷指數(shù)在玉米銅污染信息監(jiān)測中具有更好的效果。

表5 玉米葉片Cu2+污染監(jiān)測奇異性診斷指數(shù)與常規(guī)方法應用結果Table 5 Application result on the singularity diagnostic index and some convention methods for monitoring Cu2+pollution of corn leaves

圖7 玉米葉片中Cu2+含量與各監(jiān)測方法計算值擬合結果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]孫彤彤,楊可明,張 偉,等.玉米光譜小波奇異熵的銅污染監(jiān)測研究[J].環(huán)境科學學報,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]王圓圓,陳云浩,李 京,等.指示冬小麥條銹病嚴重度的兩個新的紅邊參數(shù)[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]杜華強,金 偉,葛宏立,等.用高光譜曲線分形維數(shù)分析植被健康狀況[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].農業(yè)工程學報,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]于 雷,洪永勝,周 勇,等.連續(xù)小波變換高光譜數(shù)據的土壤有機質含量反演模型構建[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].農業(yè)工程學報,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]羅玉昆,羅詩途,羅飛路,等.激光超聲信號去噪的經驗模態(tài)分解實現(xiàn)及改進[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]張 猛,王華忠,隋志強,等.基于經驗模態(tài)分解和小波變換的地震瞬時頻率提取方法及應用[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]龔志強,鄒明瑋,高新全,等.基于非線性時間序列分析經驗模態(tài)分解和小波分解異同性的研究[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]張浚哲,朱文泉,鄭周濤,等.高光譜數(shù)據的相似性測度對比研究[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]張修寶,袁 艷,景娟娟,等.信息散度與梯度角正切相結合的光譜區(qū)分方法[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.