基于VMD-MSE的玉米銅污染信息提取與預測模型

楊可明 李 燕 程 鳳 高 鵬 張 超

(中國礦業大學(北京)煤炭資源與安全開采國家重點實驗室， 北京 100083)

0 引言

近年來，礦山開采、污水排放、礦物肥料使用等造成的土壤重金屬污染問題日益嚴重，給人類生產生活帶來嚴重影響[1]。銅(Cu)是土壤重金屬污染的主要元素之一，Cu在土壤中的累積對植被生長發育產生嚴重影響[2-3]。同時，通過食物鏈進入人體，將破壞人體神經系統和免疫系統[4]。因此，植被重金屬污染的實時、有效監測顯得尤為重要。研究表明，受重金屬污染的植被光譜曲線會發生畸變[5-7]，因此很多學者利用高光譜遙感技術獲取重金屬脅迫下的植被光譜信息，旨在通過監測植被光譜特征變化來診斷重金屬污染程度。

電磁波在傳播過程中不可避免受到大氣等影響而產生噪聲，噪聲對光譜數據中污染信息的解讀造成極大干擾[8]，尋找一種有效提取作物污染弱信息的方法成為高光譜遙感重金屬污染監測領域的研究熱點。近年來越來越多的特征提取方法被引入到高光譜遙感領域，主要有小波變換[9-10](Wavelet transform，WT)、分形特征提取法[11-12](Fractal dimension，FD)、經驗模態分解[13](Empirical mode decomposition，EMD)等。上述方法在一定程度上均可較好地提取出光譜變異信息，但仍存有一定局限性。WT難以選擇合適的小波基函數，FD受到信號長度的影響，EMD分解過程中存在模態混疊現象。變分模態分解(Variational mode decomposition，VMD)是近年來提出的一種新的自適應性信號分解方法，該方法解決了EMD方法的問題，并具有很好的噪聲魯棒性[14]，目前被廣泛應用于信號分析及故障特征提取等方面[15-16]，而在高光譜遙感領域鮮見報道。

本文引入VMD理論，并結合多尺度熵(Multiscale entropy，MSE)構建一種基于變分模態分解的MSE光譜弱信息提取及定量描述模型(VMD-MSE)。利用VMD對受Cu2+污染的玉米光譜數據進行分解獲得模態分量，以MSE方法對其進行定量描述獲得VMD-MSE模型值(VM)，探討VM值與玉米葉片中Cu2+含量的相關關系，并基于VM值構建玉米葉片中Cu2+含量預測模型，以實現玉米重金屬污染程度的診斷，為農作物重金屬污染監測研究提供參考。

1 理論與方法

1.1 變分模態分解(VMD)

VMD分解過程實質為變分問題的求解過程，主要分為變分約束問題的建立和求解兩部分，最終將信號分解為數個有限帶寬的固有模態分量[17-18]。該過程包含自適應濾波組，表現了較好的噪聲魯棒性。

(1)變分約束問題的建立

假設每個“模態”是具有中心頻率的有限帶寬，變分問題描述為尋求K個模態函數uk(t)(k=1，2，…，K)，使得每個模態的估計帶寬之和最小，約束條件為各模態之和等于輸入信號f。具體步驟如下：①每個模態通過Hilbert變換計算與之相關的解析信號。②對于每個模態，通過加入指數項調整各自估計的中心頻率，把模態的頻譜變換到基帶上。③通過對解調信號進行H1高斯平滑對帶寬進行估計。④得到一個變分約束問題

(1)

其中

(u1，u2，…，uK)=u

(ω1，ω2,…，ωK)=ω

式中δ(t)——單位脈沖函數

u——分解得到的K個模態分量

ω——各模態分量的中心頻率

f——輸入的原始信號

最后求解該問題。

(2)變分約束問題的求解

①引入二次懲罰因子α和拉格朗日算子λ(t)，將約束性變分問題變為非約束性變分問題，構成擴展的拉格朗日表達式為

L({uk},{ωk},λ)=

(2)

式中α——二次懲罰因子，在高斯噪聲存在的情況下可以保證信號的重構準確度

λ(t)——拉格朗日算子，用來保持約束條件的嚴格性

②利用ADMM算法迭代搜索求取上述擴展的拉格朗日函數的鞍點，即式(1)約束變分模型的最優解，其中解得模態分量uk及中心頻率ωk分別為

(3)

(4)

1.2 多尺度熵(MSE)

MSE是基于樣本熵(Sample entropy, SpEn)的一種時間序列復雜性的度量方法，用于從不同尺度度量時間序列的復雜程度[19]。假設原始數據為X={x1，x2，…，xN}，則MSE的具體計算步驟如下[20]：

(1)給定嵌入維數m、相似容限r及尺度因子τ=1，2，…，τmax，建立新的粗粒序列

(5)

對于每個τ，原序列被分為N/τ個長度為τ的粗粒序列。

(2)計算樣本熵

①由步驟(1)所得的粗粒化序列組成m維向量

Y(i)=(yi(τ),yi+1(τ),…，yi+m-1(τ))
(1≤i≤N-m)

(6)

②計算Y(i)、Y(j)間的距離

(7)

(8)

(9)

④將維數增至m+1，重復①～③，得到Bm+1(r)。

⑤此序列的樣本熵為

(10)

當N取有限值時，式(10)表示為

(11)

利用式(6)～(11)計算每個尺度序列的SampnEn，即可得到

MSE(X)=SampnEn(y(τ),m,r)

(12)

由以上計算過程可知，MSE由τ、m和r3個參數決定。經過多次實驗分析，選取τmax=5，m=2，r=0.15δ(其中δ為X的標準差)時，光譜奇異信息規律性最顯著，能夠較好地反映光譜信號中的弱信息特征。

1.3 VMD-MSE模型

基于不同濃度Cu2+脅迫梯度下的玉米葉片光譜形態仍極為相似的特點，難以利用其進行有效的污染診斷，提出了一種VMD-MSE光譜特征提取及定量描述模型：將各脅迫梯度下的玉米葉片光譜經VMD分解，選取能夠表征光譜奇異特征的模態函數uk，計算其多尺度熵，實現光譜變異特征的定量描述。算法流程如圖1所示。

圖1 VMD-MSE模型流程Fig.1 Flow chart of VMD-MSE model

2 實驗數據及模型應用

2.1 實驗數據獲取

(1)植株培養：選用有底漏的盆缽進行“中糯1號”玉米種子培育。設置脅迫濃度為0(ck)、100、200、300、400、700、800、900 μg/g的CuSO4·5H2O溶液，將其翻土加入到玉米實驗盆缽中，每梯度設置3組平行實驗，共24盆。2016年5月6日對玉米種子進行催芽處理，出苗后向盆栽中添加NH4NO3、KH2PO4和KNO3營養液。玉米培育期間定期通風與澆水，保持適宜的培育溫度與濕度。

(2)光譜數據采集：2016年7月17日對玉米葉片反射光譜進行測量。在50 W鹵素燈光源照射下，將光譜范圍為350～2 500 nm的SVC HR-1024I型地物光譜儀探頭視場角設置為4°并垂直于葉片進行光譜采集，采集的光譜使用平面白板進行標準化。選取每盆玉米植株的老、中、新3種代表性葉片進行光譜測試，每盆獲得3組數據。各脅迫梯度光譜由3組平行實驗的9條光譜數據求均值所得。各脅迫梯度下的玉米葉片光譜曲線如圖2所示。

圖2 不同濃度Cu2+脅迫下玉米葉片光譜數據Fig.2 Spectral data of corn leaves stressed by different Cu2+ concentrations

(3)Cu2+含量測定：2016年9月16日對采集過光譜數據的葉片進行沖洗、干燥、粉碎等樣品預處理，再經高純硝酸和高氯酸消化處理后用WFX-120型原子吸收分光光度計對葉片中Cu2+含量進行測定，每梯度測量3次后取平均值作為該梯度葉片中的Cu2+含量，測量結果如表1所示。隨著脅迫梯度的增大，葉片中Cu2+含量的總體趨勢表現為先增大后減小，原因可能是隨著土壤中Cu2+增加，玉米對其吸收逐漸增大；而當土壤中Cu2+濃度過高時，對植被根系產生毒害作用，造成根部受損，對Cu2+吸收逐漸降低。本文選取脅迫梯度為0、100、200、300、400 μg/g的5組Cu2+含量數據用于建立回歸模型，脅迫梯度為700、800、900 μg/g的3組數據用于模型驗證。

2.2 VMD-MSE光譜弱信息探測

以0、100 μg/g脅迫濃度下的玉米葉片光譜數據為例，分析VMD-MSE模型的特征提取與描述

表1 不同脅迫濃度下的玉米葉片中Cu2+含量Tab.1 Cu2+ contents in corn leaves stressed by different copper concentrations μg/g

過程：

(1)對葉片光譜進行分解模態數K為2～7的VMD處理，獲得第K模態函數uk，如圖3所示。由植被反射光譜特性曲線可知，在960 nm和1 900 nm附近，植被光譜受水分和大氣的影響，出現較大干擾信息。由圖3可知，光譜信號噪聲隨著分解次數的增大而增大。K=3時，960 nm和1900 nm處噪聲較弱，且全波段奇異信息顯著；從K=4開始，1 000 nm和1 900 nm處噪聲開始出現，且奇異信息逐漸弱化。因此選擇K=3既實現模態分離，又有效抑制噪聲，可有效提取玉米銅污染脅迫信息。

圖3 VMD分解結果Fig.3 Results of decomposition of VMD

(2)計算由VMD分解得到的第3模態函數u3的多尺度熵，即為VMD-MSE的模型值(記為VM)，不同尺度下的模型值記為VMi(i=1，2，…，τmax)。

3 結果與分析

3.1 VMD-MSE污染監測

對0、100、200、300、400 μg/g脅迫濃度下的玉米葉片光譜應用VMD-MSE模型進行污染信息提取與測度。各脅迫濃度下的玉米葉片光譜的VM值如圖4所示。不同脅迫濃度下，VM曲線變化趨勢有所差異，但總體而言，曲線變化趨勢相似，即隨著τ的增大而增大。在所有尺度因子上，VM值與作物受污染程度呈負相關關系，即受污染越嚴重，VM值越小。顯然，各脅迫濃度下的光譜細節特征通過不同尺度VM(即VMi)得到體現，能更好地揭示光譜曲線畸變特征。

圖4 不同脅迫濃度下玉米葉片光譜的VM值Fig.4 VM values of corn leaf spectrum under different stress concentrations

圖5 不同脅迫梯度下各尺度VM值與玉米葉片中Cu2+含量相關關系Fig.5 Relationships between VM on five scales and Cu2+ contents in corn leaves stressed by different Cu2+ concentrations

不同脅迫濃度下，VMi與葉片中Cu2+含量的相關變化如圖5所示，兩者呈現相反的變化趨勢。隨著脅迫濃度的增加，葉片中Cu2+含量先增加后減小，相應地VMi表現為先減小后增大的趨勢，兩者相關系數絕對值均大于0.90，其中VM1與Cu2+相關性最強，相關系數為-0.968。據此可認為，各尺度VM均可在一定程度上診斷玉米受銅脅迫程度，其中VM1效果最優。

3.2 VM值與玉米銅含量的回歸分析

為了進一步確定VMD-MSE模型值與玉米銅脅迫的定量關系，以不同脅迫濃度下玉米葉片中Cu2+含量為因變量y，以VMi為自變量x進行回歸分析，VMi值與銅含量的線性回歸方程、決定系數(R2)和均方根誤差(RMSE)如表2所示。VM1～VM5的回歸決定系數均大于0.88，且顯著性水平P<0.05。整體上，從VM1到VM4，R2逐漸減小，RMSE逐漸增大，這可能是由于隨時間序列復雜度的增加，樣本熵估計準確率下降所致。

為了驗證上述模型的可靠性與準確度，選用700、800、900 μg/g 3個脅迫濃度下的樣本數據對模型進行檢驗，圖6描述了玉米葉片中Cu2+含量實測值與模型預測值的相關變化。各預測模型的Cu2+含量預測值與實測數據具有很好一致性。VM1和VM2回歸模型的預測結果較優，相關系數大于0.70，其中VM1模型預測結果最優，相關系數為0.992；VM3～VM5回歸模型預測效果較差，相關系數小于0.650，可能是由于隨著污染程度的增加，光譜信號的奇異信息隨機性增加，且隨度量尺度的增大，熵估計準確率下降，使VM3～VM5模型穩定性及適用性變差。因此，通過VMD-MSE模型值VM1可以建立玉米Cu污染程度的預測模型，從而實現監測玉米Cu污染狀況的目的。

表2 基于不同尺度VM值的玉米葉片中Cu2+含量回歸模型Tab.2 Regression models for VM on different scales with Cu2+ content in corn leaves

圖6 玉米葉片中Cu2+含量實測值與預測值的相關關系Fig.6 Relationships between measured and predicted values of Cu2+ contents in corn leaves

4 結論

(1)將VMD運用到高光譜信息提取中，從相似光譜信號中挖掘有利于診斷重金屬Cu污染的畸變信息。結果表明，對原始光譜數據進行3次VMD分解后，可有效抑制噪聲，實現模態分離，從而提取出隱藏在光譜噪聲中的污染弱信息。

(2)將MSE方法應用于VMD分解結果，揭示了污染信息不同尺度的內在差異。

(3)基于VMD和MSE構建的VMD-MSE模型可有效診斷作物銅污染程度，模型值VM與玉米葉片中Cu2+含量呈現顯著負相關，兩者回歸精度達到0.88以上，其中VM1模型預測效果最優。本研究構建的預測模型能在一定程度上預測玉米受污染程度，但具有一定局限性，即預測精度隨污染程度增大而降低，這種情況可通過增加樣本數量對模型進行優化而得到改善。