郭輝，石海*，張全旺

（1.安徽理工大學空間信息與測繪工程學院，安徽 淮南 232001；2.安徽理工大學礦山環(huán)境與災害協(xié)同監(jiān)測煤炭行業(yè)工程研究中心，安徽 淮南 232001）

近年來，由于礦山開采、污水灌溉、大氣沉降等因素導致農(nóng)田土壤中重金屬含量普遍升高[1-2]。重金屬污染下農(nóng)作物的產(chǎn)量和質(zhì)量均受到嚴重影響，并且積累在農(nóng)作物中的重金屬一旦進入食物鏈，將最終危害人類健康[3-4]。因此，快速、有效地監(jiān)測農(nóng)作物重金屬污染對于農(nóng)業(yè)生產(chǎn)、糧食安全具有重要意義。傳統(tǒng)的光學檢測法、電化學檢測法和生物檢測法需要對樣品進行破環(huán)性采集，檢測過程繁瑣且容易造成二次污染，不僅費時費力且難以滿足快速、實時、大面積檢測的需求[5]。

高光譜遙感技術具有光譜連續(xù)精細、波段數(shù)多、光譜分辨率高等特點，為農(nóng)作物重金屬污染監(jiān)測提供了技術保障[6]。目前，相關學者通常采用簡單的光譜變換或提取植被指數(shù)的方法，開展農(nóng)作物重金屬污染監(jiān)測研究。劉來等[7]利用光譜微分技術篩選鎘污染下油菜葉片光譜敏感波段；Yu 等[8]基于光譜一、二階微分技術選取重金屬污染下冬小麥光譜特征；程鳳等[9]構建一種植被指數(shù)探測玉米重金屬銅污染程度；Martinez等[10]將鋰脅迫下農(nóng)作物光譜反射率提取的植被指數(shù)作為預測因子。以上研究取得了一定成果，但在重金屬污染監(jiān)測方面缺乏針對性，且易受農(nóng)作物中色素、水分、有機質(zhì)等生化參數(shù)的影響[11]。近年來，希爾伯特-黃變換（Hilbert-Huang transform，HHT）時頻分析技術被廣泛用于非線性非平穩(wěn)信號的處理與分析，如張世文等[12]將HHT 算法運用于探地雷達數(shù)據(jù)分析，有效探測出覆土層厚度；王英樂等[13]利用HHT算法提取微震信號的特征參數(shù)，實現(xiàn)微震監(jiān)測與預警；Zao 等[14]采用HHT 算法分析噪聲語音信號，有效估計噪聲語音的基本頻率；Yin 等[15]運用HHT 算法處理電磁輻射信號，提高了煤層破壞前特征識別程度。以上研究表明HHT 算法在非線性非平穩(wěn)信號的特征提取上有一定優(yōu)勢，但HHT 算法在高光譜數(shù)據(jù)處理與分析的應用較少。因此，本研究引入HHT 算法用于玉米葉片高光譜數(shù)據(jù)分析，以探測玉米重金屬銅離子污染信息。

考慮到HHT 算法容易受到信號噪聲的影響[16]，本研究結合連續(xù)小波變換（Continuous wavelet transform，CWT）與HHT 算法的優(yōu)勢，首先對玉米葉片原始光譜進行CWT 處理去除信號噪聲，然后進行HHT時頻分析，通過研究不同濃度Cu2+脅迫下玉米葉片光譜信號的瞬時能量特征進行污染信息探測。同時與常見的紅邊位置（REP）、紅邊歸一化指數(shù)（NDVI705）和紅邊植被脅迫指數(shù)（RVSI）等植被指數(shù)監(jiān)測農(nóng)作物重金屬污染方法作對比分析，進一步驗證CWT-HHT 方法在探測玉米葉片Cu2+污染信息方面的可行性與有效性，以期為農(nóng)作物重金屬污染信息探測提供一種新的方法，并為農(nóng)作物重金屬含量反演奠定基礎。

1 材料與方法

1.1 實驗設計

1.1.1 玉米植株培養(yǎng)

以“中糯1 號”品種玉米作為研究對象，優(yōu)選籽粒飽滿、大小均勻的種子進行重金屬Cu2+脅迫土培實驗。實驗采用CuSO4·5H2O分析純?nèi)芤好{迫玉米發(fā)育生長，設置了4 種脅迫梯度，分別為0、100、300、500μg·g-1，分別標記為CK0、Cu100、Cu300、Cu500。每個濃度梯度設置3 組平行實驗，共計12 盆玉米盆栽。2019年5月對玉米種子進行催芽處理，出苗后定期澆灌NH4NO3、KH2PO4和KNO3營養(yǎng)液。所有盆栽玉米的培育條件除了土壤中的Cu2+濃度不同，均在同一條件下進行。

1.1.2 數(shù)據(jù)采集

玉米出苗后使用美國SVC 公司生產(chǎn)的SVC HR-1024I 高性能地物光譜儀采集玉米葉片光譜，光譜范圍為350～2 500 nm。選取位于玉米中部的葉片在室內(nèi)密閉的環(huán)境下進行光譜采集。以50 W 鹵素燈為光源，使用4°視場角的探頭垂直于玉米葉片表面，探頭距離葉片表面40 cm。為了避免因光源強度分布不均勻?qū)е掳惦娏髟肼曈绊懝庾V數(shù)據(jù)質(zhì)量，每次測量玉米葉片光譜前，先用白板進行標準化處理。光譜儀連續(xù)測量玉米葉片反射光譜3 次，隨后平均3 次測量值作為該脅迫濃度下玉米葉片光譜數(shù)據(jù)。

將采集完光譜的玉米葉片洗凈后放到烘箱中烘烤至質(zhì)量恒定，用剪刀剪成小塊裝入樣品袋中并貼上標簽。經(jīng)過預處理后，在化學分析實驗室中采用WFX-120 原子吸收分光光度計進行玉米葉片樣品Cu2+含量測定，每份葉片樣品平均分成3份，將實驗測出的Cu2+含量值取平均得到該脅迫濃度下玉米葉片樣品中Cu2+的含量。實驗中測得0、100、300、500μg·g-1脅迫濃度下的玉米葉片Cu2+含量分別為：9.77、29.62、75.78、54.51μg·g-1。

1.2 連續(xù)小波變換

連續(xù)小波變換通過小波基函數(shù)將光譜信號分解為低頻和高頻兩個部分，在低頻部分具有較低的時間分辨率和較高的頻率分辨率，而在高頻部分具有較高的時間分辨率和較低的頻率分辨率，適合于各類信號噪聲的處理。連續(xù)小波變換公式為[17]：

式中：Wf（a，b）為小波變換結果，f（λ）為玉米葉片高光譜反射率，ψa，b（λ）為小波基函數(shù)，a為尺度因子，b為平移因子。

1.3 希爾伯特-黃變換

希爾伯特-黃變換（HHT）主要用于非線性非平穩(wěn)信號的處理，包含經(jīng)驗模態(tài)分解（Empirical mode decomposition，EMD）和Hilbert變換兩部分[18]。

1.3.1 經(jīng)驗模態(tài)分解

EMD 可以把任意一個復雜的非平穩(wěn)信號自適應地分解為有限個窄帶信號，稱為本征模態(tài)函數(shù)（Intrinsic mode function，IMF），每個IMF 需要滿足以下兩個條件：（1）零點數(shù)目與極值點數(shù)目相同或至多相差1 個；（2）函數(shù)關于局部平均對稱。最終將光譜信號x（t）分解為若干個IMF 分量ci（t）與單調(diào)殘差rn（t）之和：

其具體算法如下：

（1）求取信號x（t）的極值點；（2）采用三次樣條函數(shù)擬合信號的極小值點與極大值點得到上下包絡線；（3）計算上下包絡的平均值；（4）將原始信號減去上下包絡的平均值，提取IMF分量；（5）重復以上步驟至均值趨近于零；（6）計算殘差；（7）重復運算直到殘差中不含IMF函數(shù)。

1.3.2 Hilbert變換

信號經(jīng)過經(jīng)驗模態(tài)分解后，就可以對每一個IMF作Hilbert變換，公式為：

式中：H[ci（t）]為第i階IMF 分量ci（t）的Hilbert 變換函數(shù)。

根據(jù)歐拉公式可得其解析信號為：

式中：ai（t）為解析信號的幅度或能量；θi（t）為解析信號的相位。

由此可以求出IMF分量的瞬時頻率wci（t）：

則原始信號x（t）的希爾伯特譜表示為：

式中：n為IMF個數(shù)；Re為殘差項。

則Hilbert瞬時能量譜E（t）為：

1.3.3 CWT-HHT銅污染信息探測流程

基于CWT 與HHT 所構建的CWT-HHT 探測方法處理具體過程如下：（1）對原始玉米葉片光譜數(shù)據(jù)進行CWT 分解，有利于光譜信號降噪；（2）提取連續(xù)小波分解后的第5 層小波系數(shù)進行EMD 分解，得到有限個單一頻率的IMF 分量；（3）用Hilbert 變換對IMF分量進行處理，獲取解析信號，以提取信號的瞬時頻率信息；（4）將解析信號表示為希爾伯特譜形式；（5）對各IMF 分量的希爾伯特譜的平方在頻率域進行積分，得到光譜信號的瞬時能量；（6）根據(jù)波長、頻率和能量構建三維信息圖，選取瞬時能量峰值E作為特征參數(shù)，以探測不同濃度銅污染下玉米葉片光譜重金屬污染信息。

1.4 數(shù)據(jù)處理

本研究采用Matlab 2018b 編程對玉米葉片光譜數(shù)據(jù)進行連續(xù)小波變換與希爾伯特-黃變換處理分析，利用SPSS 進行相關性分析，在Excel 2013 中進行顯著性檢驗，采用Origin 2018進行制圖。

2 結果與討論

2.1 光譜波段選擇

實驗采集的一組玉米葉片光譜如圖1 所示。從光譜曲線可以看出，不同濃度Cu2+脅迫下葉片光譜均表現(xiàn)出高度的相似性，在可見光范圍，550 nm 附近因光合色素對綠光的強烈反射形成綠峰，650 nm 附近因葉綠素對紅光的吸收形成紅谷，750 nm 附近因反射率急劇上升形成葉片光譜最明顯的紅邊特征；在近紅外波段范圍，760～1 250 nm 附近受葉片細胞內(nèi)部結構影響形成一個高反射平臺，在1 450 nm 與1 930 nm呈現(xiàn)出兩個明顯的水吸收帶，難以直接進行重金屬污染探測[19]。因此，需要選取合適的光譜波段以提取玉米葉片重金屬污染光譜特征。為了提高探測精度，去除噪聲較大的350～400 nm 波段以及受水分強吸收影響的1 300～2 500 nm 波段光譜，選取400～1 300 nm 波段進行分析。

圖1 不同Cu2+脅迫濃度的玉米葉片光譜曲線Figure 1 Spectral curves of corn leaves under different Cu2+stress gradients

2.2 基于CWT-HHT污染監(jiān)測

2.2.1 小波系數(shù)提取

選取不同的小波基取得的結果不盡相同。經(jīng)過多次對比分析，采用Daubechies 小波系中的“Db5”小波函數(shù)對重金屬Cu2+脅迫下的玉米葉片光譜進行連續(xù)小波變換，提取分解至第5層的小波系數(shù)，如圖2所示，可以看出，在400～1 300 nm 波段范圍內(nèi)小波信號較原始光譜曲線更加平滑，并且在指示植物重金屬污染的紅邊（680～750 nm）這一敏感位置小波系數(shù)變化劇烈，可見連續(xù)小波變換既有效降低了光譜噪聲又增強了玉米受重金屬Cu2+污染信息。

圖2 不同Cu2+脅迫濃度的玉米葉片光譜小波系數(shù)Figure 2 Spectral wavelet coefficients of corn leaves under different Cu2+stress gradients

2.2.2 CWT-HHT銅污染信息探測

將連續(xù)小波變換后的光譜數(shù)據(jù)進一步作HHT 處理得到希爾伯特譜，并將希爾伯特譜的平方在頻率域上進行積分，構建葉片光譜信號的波長-頻率-能量三維信息圖。如圖3 所示，橫坐標表示波長，縱坐標表示對應波長信號的頻率，能量以顏色深淺顯示，顏色越深對應的波段能量就越高。

圖3 不同Cu2+脅迫濃度下的希爾伯特瞬時能量Figure 3 Hilbert instantaneous energy under different Cu2+stress gradients

由圖3 分析可知，隨著Cu2+脅迫濃度的不斷增加，玉米葉片光譜的Hilbert 瞬時能量峰值呈現(xiàn)先增加后降低的趨勢。不同濃度Cu2+污染下，瞬時能量峰值均在700 nm 波段附近出現(xiàn)，并且在1 100 nm 波段附近出現(xiàn)第2 個能量峰，表明葉片光譜在這兩個波段存在能量波動，并且在700 nm 波段附近能量波動異常明顯，這與提取的小波系數(shù)劇烈變化位置一致。而信號頻率集中分布在0～0.5 Hz較窄的部分，這是因為原始光譜經(jīng)連續(xù)小波變換后提取的第5 層小波系數(shù)的頻率單一，去除了不同頻率信號的干擾。

為了定量描述玉米受重金屬Cu2+污染程度，提取瞬時能量峰值E作為探測指標，圖4 為能量峰值E與玉米葉片Cu2+含量變化趨勢?？梢钥闯?，不同濃度Cu2+脅迫下，能量峰值E與葉片Cu2+含量變化趨勢一致。因此，可以認為瞬時能量峰值E作為玉米重金屬Cu2+污染監(jiān)測的指標具有可行性。

圖4 能量峰E與葉片Cu2+含量變化趨勢Figure 4 Energy peak E and transformation trend of Cu2+content in leaves

2.3 CWT-HHT探測方法應用與驗證

由于玉米受到不同濃度重金屬Cu2+污染后其葉片光譜信號差異微弱，難以直接進行污染程度判別。實驗首先對采集的玉米葉片光譜進行CWT 處理，獲得單一頻率的小波系數(shù)。葉片光譜信號經(jīng)CWT 處理后，能夠降低光譜信號噪聲并增強玉米重金屬污染信息，可以提高HHT算法探測玉米受Cu2+污染程度信息的精度。隨后利用HHT 算法獲取光譜信號的三維信息圖并提取信號的特征值，以描述在不同濃度Cu2+污染下玉米污染程度信息。因HHT 算法能夠自適應分解光譜信號，在消除人為因素干擾的同時可以進一步突出信號本身的局部特征，適合玉米在重金屬Cu2+污染下光譜畸變特性的探測研究。

為了驗證CWT-HHT 方法在玉米Cu2+污染信息探測方面的有效性與優(yōu)越性，將其與紅邊位置（REP）、紅邊歸一化指數(shù)（NDVI705）和紅邊植被脅迫指數(shù)（RVSI）3 個常規(guī)的植被指數(shù)監(jiān)測重金屬污染方法進行對比分析。其中，REP 為反射光譜一階微分在670～780 nm 最大值對應波長，NDVI705計算公式為（R750-R705）/（R750+R705），RVSI 計算公式為[（R712+R752）/2]-R732。

不同濃度Cu2+脅迫下葉片光譜的CWT-HHT 探測方法結果E值和REP、NDVI705、RVSI 與葉片中Cu2+含量的相關系數(shù)以及顯著性檢驗結果如表1 所示。分析發(fā)現(xiàn)，E值與玉米葉片中Cu2+含量的相關系數(shù)最高達到0.980 3，表明玉米受重金屬Cu2+污染程度越大，葉片光譜的瞬時能量特征越顯著。同時，將實驗獲取的玉米葉片中Cu2+含量與各光譜參數(shù)進行線性擬合分析，擬合結果如圖5 所示。結果表明CWTHHT 探測方法提取的E值與葉片中Cu2+含量擬合程度最優(yōu)，其擬合決定系數(shù)R2最高，達到了0.949 8，遠高于常規(guī)植被指數(shù)監(jiān)測方法的擬合系數(shù)。綜上所述，CWT-HHT方法可以有效探測出Cu2+污染下玉米葉片光譜微弱的差異信息，也驗證了其在農(nóng)作物重金屬污染監(jiān)測方面的可行性。

圖5 不同方法與玉米葉片中Cu2+含量的擬合結果Figure 5 Fitting results of different methods and Cu2+content in corn leaves

3 結論

（1）隨著Cu2+脅迫濃度的增加，玉米葉片光譜瞬時能量峰值參數(shù)呈現(xiàn)先升高后降低的趨勢，與相同濃度下葉片Cu2+含量變化趨勢一致，因此可將其作為描述玉米葉片受Cu2+污染程度的指標。

（2）基于CWT-HHT 方法提取的瞬時能量峰值參數(shù)能夠有效甄別相似光譜間微弱的差異信息，且與葉片Cu2+含量的相關系數(shù)達到了0.980 3，擬合決定系數(shù)達到了0.949 8。

（3）通過與常見的紅邊位置（REP）、紅邊歸一化指數(shù)（NDVI705）和紅邊植被脅迫指數(shù)（RVSI）等植被指數(shù)監(jiān)測方法對比，進一步驗證了CWT-HHT 方法在玉米葉片Cu2+污染信息探測方面具有較好的有效性與可行性。