基于CWT-HHT的玉米葉片銅離子污染信息探測

郭輝，石海*，張全旺

（1.安徽理工大學空間信息與測繪工程學院，安徽 淮南 232001；2.安徽理工大學礦山環境與災害協同監測煤炭行業工程研究中心，安徽 淮南 232001）

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

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

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

1 材料與方法

1.1 實驗設計

1.1.1 玉米植株培養

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

1.1.2 數據采集

玉米出苗后使用美國SVC 公司生產的SVC HR-1024I 高性能地物光譜儀采集玉米葉片光譜，光譜范圍為350～2 500 nm。選取位于玉米中部的葉片在室內密閉的環境下進行光譜采集。以50 W 鹵素燈為光源，使用4°視場角的探頭垂直于玉米葉片表面，探頭距離葉片表面40 cm。為了避免因光源強度分布不均勻導致暗電流噪聲影響光譜數據質量，每次測量玉米葉片光譜前，先用白板進行標準化處理。光譜儀連續測量玉米葉片反射光譜3 次，隨后平均3 次測量值作為該脅迫濃度下玉米葉片光譜數據。

將采集完光譜的玉米葉片洗凈后放到烘箱中烘烤至質量恒定，用剪刀剪成小塊裝入樣品袋中并貼上標簽。經過預處理后，在化學分析實驗室中采用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 連續小波變換

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

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

1.3 希爾伯特-黃變換

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

1.3.1 經驗模態分解

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

其具體算法如下：

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

1.3.2 Hilbert變換

信號經過經驗模態分解后，就可以對每一個IMF作Hilbert變換，公式為：

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

根據歐拉公式可得其解析信號為：

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

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

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

式中：n為IMF個數；Re為殘差項。

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

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

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

1.4 數據處理

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

2 結果與討論

2.1 光譜波段選擇

實驗采集的一組玉米葉片光譜如圖1 所示。從光譜曲線可以看出，不同濃度Cu2+脅迫下葉片光譜均表現出高度的相似性，在可見光范圍，550 nm 附近因光合色素對綠光的強烈反射形成綠峰，650 nm 附近因葉綠素對紅光的吸收形成紅谷，750 nm 附近因反射率急劇上升形成葉片光譜最明顯的紅邊特征；在近紅外波段范圍，760～1 250 nm 附近受葉片細胞內部結構影響形成一個高反射平臺，在1 450 nm 與1 930 nm呈現出兩個明顯的水吸收帶，難以直接進行重金屬污染探測[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污染監測

2.2.1 小波系數提取

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

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

2.2.2 CWT-HHT銅污染信息探測

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

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

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

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

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

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

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

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

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

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

3 結論

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

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

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