銅脅迫下玉米葉片的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。