程 鳳,楊可明*,王 敏,2,李 燕,高 鵬,張 超
(1.中國礦業大學(北京)煤炭資源與安全開采國家重點實驗室,北京 100083;2.華北理工大學,河北唐山 063210)
重金屬的空間分布及其對生態系統結構功能的影響一直備受關注,其監測與治理已成為國內外科研工作的熱點[1-3]。根據2014年發布的“全國土壤污染狀況調查報告”[4],我國重金屬污染問題嚴峻,全國污染物的總超標率約為16%,其中銅污染物點位超標率達2.1%。重金屬污染具有累積性、隱蔽性、不可逆轉性等特點,其在土壤、水體中的富集會對生物有機體產生嚴重影響。部分重金屬離子(如Cu2+)是植物新陳代謝的必要元素,但過量的重金屬會對植物細胞產生毒害[1]。并且還會阻礙植物的光合作用,抑制光合產物的輸配,影響農作物產量[5-6]。一定濃度的重金屬會在植物的根、莖、葉、籽粒部位產生不同程度的累積,進而通過食物鏈進入人體,危害人類健康。近60年,學者已對植物重金屬的吸收和遷移做了大量研究,并提出了土壤-植物系統中元素解吸、遷移、吸收的聯合數學模型[7]。重金屬污染的監測與治理方面相繼發展出物理、化學、生物等方法。傳統的植物重金屬污染監測方法需使用專業儀器設備、定期觀測采樣來實現,成本高、過程繁瑣、費時費力。而高光譜遙感技術具有較高的光譜分辨率、生化成分反演能力強、便捷、可操作性高和視野寬等優點,為實現大范圍實時動態監測提供了可能[8]。因此運用高光譜遙感監測農作物重金屬污染是一項具有一定優勢的遙感前沿技術。
近年來,有很多學者采用光譜數據的多階微分、微分光譜角正切、小波變換、經驗模態分解(Empirical mode decomposition,EMD)等方法[8-11]對植物葉片原始光譜進行變換以提取反映重金屬污染的異常信息,探究不同濃度重金屬脅迫下植物的生長狀況。就光譜變換與分解方法而言,小波變換缺乏自適應性,必須根據光譜信號的特征選擇合適的小波基進行分析,而EMD能以信號本身的局部時間尺度來進行分解,具有良好的自適應性,但分解過程存在模態混疊現象[12]。集成經驗模態分解(Ensemble EMD,EEMD)方法克服了模態混疊這一缺點,已經廣泛應用于心電信號處理[12]、電力系統故障信號檢測[13]、混合信號分離識別[14]等不同領域,同時在土壤剖面反射光譜消噪方面也取得了一定進展[15]。Mallat算法是基于小波分析改進的信號分解和重構方法,已用于數字圖像處理[16]、雷達回波去噪[17]、圖像增強變換[18]、空氣質量指數預測[19]等方面。本文利用EEMD克服模態混疊這一優勢,結合MA算法和光譜的分形維數(Fractal dimension,FD)方法,構建不同濃度Cu2+脅迫下玉米葉片污染信息的EEMD-MA-FD光譜診斷模型,用于探測玉米葉片在不同Cu2+污染程度下的光譜弱差信息,為農作物重金屬污染監測提供一種新的高光譜遙感應用思路或方法。
1.1.1 植株培養
實驗選取“密糯8號”玉米作為研究對象,采用有底漏花盆培育玉米種子。對栽培所用顆粒土壤進行碾碎、去除石子及草根雜質、過篩等處理,以逐層噴灑、翻土混合的方式將不同濃度CuSO4·5H2O與土壤充分攪拌均勻并標注相應脅迫濃度。根據《土壤環境質量標準》(GB 15618—1995)中土壤污染標準設置Cu2+脅迫濃度分別為 0、200、400、600、800 μg·g-1,記為CK、Cu(200)、Cu(400)、Cu(600)、Cu(800),同一濃度設置3組平行樣,共15盆。其中CK組為空白對照組,土壤中本身含有Cu2+,3組平行樣Cu2+含量均值為19.6 μg·g-1。2017年5月10日對玉米種子進行催芽處理,5月12日栽培植株,且于生長期定期澆灌并添加適量NH4NO3、KNO3、KH2PO4營養液,培育期間保持各盆栽生長環境一致且通風良好。
1.1.2 光譜數據獲取
2017年7月19日,玉米穗期,在暗室內進行光譜數據采集。采用350~2500 nm波譜范圍的SVC HR-1024I高性能地物光譜儀進行測定,測量過程使用50 W的鹵素燈作為光源和25°視場角的探頭,探頭垂直于葉片表面5 cm。為了防止其他物體對玉米葉片光譜的影響,將玉米葉片裁剪后,在遮蓋黑布的實驗臺上進行測定。使用專用白板標準化光譜反射系數。實驗選取位于植株上、中、下部位具備代表性的新、中、老葉片的葉中位置進行光譜數據測量,各葉片每次測取3條光譜,取平均值作為最終結果。不同脅迫濃度下玉米葉片均值光譜數據如圖1所示。
1.1.3 玉米葉片Cu元素含量測定
各項實驗完成后,采集各脅迫濃度的玉米葉片,進行沖洗、烘干、粉碎等處理,對各個樣品封裝保存并進行標注、編號。在相同規范條件下,經微波處理后,采用電感耦合等離子發射光譜儀(ICP-OES)對已進行光譜數據采集的老、中、新葉片分別進行Cu元素含量測定,取均值作為每組脅迫濃度葉片Cu2+含量,所測結果如表1。
1.2.1 集成經驗模態分解(EEMD)

圖1 不同Cu2+脅迫濃度下玉米葉片光譜曲線Figure 1 Maize leaf spectra curves stressed by different Cu2+concentrations

表1 玉米葉片中Cu2+含量(μg·g-1)Table 1 Cu2+content in maize leaves(μg·g-1)
EEMD由Wu等[20]提出,主要用于解決EMD存在的虛假分量和模態混疊問題[21]。EEMD的核心是通過多次在原始信號中加入等長度的正態分布白噪聲序列,以保證不同尺度信號區域自動映射到與背景白噪聲相關的尺度上。將信號和噪聲重新組合為一個整體,視為待分解信號,對其進行多次EMD分解,對每次分解的本征模態函數(Intrinsic mode function,IMF)取均值,作為最終結果。需要注意的是,由于每次EMD分解時人為添加的白噪聲不同,且噪聲之間不具備相關性,因而每次分解所得的IMF分量均會存在微弱差異,當使用足夠測試的全體均值時,人為添加的噪聲就會被抵消,信號本身是穩固持久的[20]。設一條原始的光譜信號為x(t),則其EEMD處理步驟如下:
Ⅰ為原始光譜信號數據添加一個正態分布的白噪聲序列xm(t),即

式中:k為加入背景白噪聲的幅值系數;nm(t)為第m次加入的白噪聲。
Ⅱ將帶有白噪聲的光譜信號數據分解成一組IMF分量。
Ⅲ重復執行步驟Ⅰ和步驟Ⅱ,但每次加入不同的白噪聲序列。
Ⅳ將每次獲得的IMF分量各自取均值作為最終結果。

式中:ci,m為第m次EMD分解得到的第i個分量;N為EMD分解的次數。
1.2.2 Mallat算法(MA)
Mallat算法(記為MA)是由Mallat受塔式算法啟發于1988年提出的正交小波構造方法,是一種建立在多分辨率基礎上進行的分解和重構算法[22]。MA主要是把信號分解為離散平滑分量和離散細節分量,若一光譜信號f(t),Aj為f(t)在第j層低頻部分(近似部分)的小波系數,Dj為f(t)在第j層高頻部分(細節部分)的小波系數,則MA分解光譜信號的具體算法為[23]:

式中:t為波長;j為層數,j=1,2,…,log2N;k為離散平移位置,k=1,2,…,N;H、G為濾波器系數。可見,Aj是通過第2j-1尺度的小波系數Aj-1與分解濾波器H卷積獲得,Dj是通過第2j-1尺度的小波系數Aj-1與分解濾波器G卷積獲得[23]。
1.2.3 分形維數(FD)
分形理論是一種非線性科學理論,可根據某物質的自相似性描述復雜無序、標度不變的系統,探究混沌事物內部結構[24]。FD是刻畫分形體復雜結構的重要工具,反映了復雜形體占有空間的有效性。計算FD的方法有很多,例如量規法、盒維數法、結構函數法等[25-26]。本文采用盒維數法計算光譜曲線的FD,記為D,則[27]:

式中,設A是非空集合,在歐式距離下,用邊長為1/Tn的小盒子緊鄰地包含A,其中Nn(A)表示包含A所需要的最小盒子數,D為集合A的盒維數。
結合EEMD、MA和FD理論與方法,用于不同脅迫濃度下玉米葉片光譜的Cu污染診斷分析,稱為玉米葉片Cu2+污染信息的EEMD-MA-FD光譜診斷模型。
玉米葉片光譜信號可視為非線性非穩態信號,對不同脅迫濃度的光譜進行EEMD處理,為保證分解有效性,加入的白噪聲幅值通常設置在0.01~0.5之間,EMD分解次數設置在100~300之間[20]。本次添加白噪聲幅值系數k=0.01,執行EMD總次數為300次。分解后獲得各階IMF(CK時光譜分解的各階IMF如圖2所示),光譜信號可分解為10階IMF分量和1階殘余項r。EEMD分解的每一階IMF分量的振幅與頻率均不同,每階IMF分量在各自特征尺度上具有不同的變化特性。選擇可見光及近紅外部分波段波長范圍在410~780 nm進行探究,對原始光譜反射率與各階IMF分量之間進行相關性分析,相關系數如表2所示。根據表2,比較各IMF分量與原始光譜反射率的相關性,IMF4分量和IMF7分量與原始光譜反射率的相關系數較高,分別為0.755、0.811。并且原始光譜具有在0.54 μm附近形成反射峰,在0.76 μm附近反射率急劇上升的特征,IMF4分量較IMF7分量具有與原始光譜相似的特點,很大程度上保留了原始光譜的特征信息,因此選擇IMF4分量對監測重金屬Cu2+污染程度具有一定意義。

圖2 玉米葉片光譜曲線及各IMF分量Figure 2 Maize leaf spectra curves and various components of IMF

表2 各階IMF分量與原始光譜相關系數Table 2 Correlation coefficient between IMF components and original spectra
對不同Cu2+脅迫濃度的IMF4進行MA變換處理,IMF4分量被分解為逼近信號(sig1)和細節信號(sig2),分解結果如圖3所示。分析圖3發現,不同脅迫梯度sig1與sig2信號在不同波長所對應的振幅與頻率均不相同。信號高頻部分一般包含突變特征和細節信息等,不同脅迫梯度IMF4分量經過MA低通濾波處理后,分解所得的低頻部分(sig1)中細節信息消失,噪聲成分被有效去除;經MA高通濾波處理后,IMF4分量中低頻部分數據被零替代,高頻部分(sig2)保留。
基于IMF4分量分解的sig1信號和sig2信號計算曲線的分形維數,將結果與相應脅迫濃度玉米葉片Cu2+含量進行相關性分析。sig2分量曲線分形維數與玉米葉片Cu2+含量之間相關性較弱,而sig1信號具有一定的穩定性,其曲線分形維數與對應脅迫濃度玉米葉片Cu2+含量相關性較強,可用指標值刻畫其變化情況。將sig1信號分形維數結果記為EEMF,EEMF隨玉米葉片Cu2+含量增加而降低,如圖4所示,相關系數為-0.942 2。

圖3 不同脅迫濃度IMF4分量Mallat變換分解結果Figure 3 The results of Mallat transform decomposition of IMF4 components with different stress gradients
EEMD-MA-FD模型的EEMF隨著玉米葉片Cu2+含量增加而降低,為了驗證模型的優越性,將模型運用結果分別與傳統的監測參數藍邊最大值、紅邊最大值、紅肩最大值進行對比,傳統參數的具體定義如表3所示,分析結果見圖5。藍邊最大值與葉片中的Cu2+含量呈正相關,相關系數為0.907 5,紅肩最大值與葉片中Cu2+含量呈負相關,相關系數為-0.801 8,而紅邊最大值與葉片中Cu2+含量無明顯相關性。將經EMD分解、MA變換處理后的分形維數結果記為EMF,分析發現EMF與玉米葉片Cu2+含量無明顯相關性,如圖5d所示。由此可證明EEMD-MA-FD模型對監測重金屬污染程度具有優越性。

圖4 葉片Cu2+含量與EEMF擬合結果Figure 4 Fitting results on the EEMF and the Cu2+contents of leaves

表3 玉米葉片的光譜特征參數名稱及定義Table 3 The name and definition of spectral characteristic parameters of maize leaves′
為了驗證EEMD-MA-FD光譜診斷模型的穩定性,利用2016年培育的“中糯1號”玉米光譜數據進行對比驗證。隨機選擇Cu2+脅迫濃度為0、400、800 μg·g-1的玉米葉片原始光譜數據作為研究對象,計算各脅迫濃度光譜在EEMD-MA-FD模型的EEMF值,并選擇紅肩最大值與其作對比,結果如圖6所示。EEMF隨土壤重金屬脅迫濃度的增加呈下降趨勢,且與玉米葉片Cu2+含量呈負相關,相關系數為-0.993 7。紅肩最大值與葉片中的Cu2+含量無明顯相關性。由此可證明EEMD-MA-FD模型穩定性較強,對各品種玉米光譜均可進行良好監測。

圖5 葉片Cu2+含量與各監測方法計算值的擬合結果Figure 5 Fitting results on the computing values of the monitoring methods and the Cu2+contents of maize leaves

圖6 葉片Cu2+含量與EEMF和紅肩最大值的擬合結果Figure 6 Fitting results on the EEMF and the maximun of red shoulder and the Cu2+contents of maize leaves
本研究探索性地將EEMD應用于玉米葉片高光譜微弱信息提取中。結果表明,經EEMD分解后IMF4分量與原始光譜相關性較高,并有效抑制噪聲。同時結合Mallat算法及盒維數法分形維數構建EEMD-MA-FD光譜診斷模型,并與藍邊最大值、紅邊最大值等傳統監測方法及經EMD分解、MA變換后的分形維數結果分析比較。EEMD-MA-FD模型的EEMF值與玉米葉片Cu2+含量相關性較強,能有效甄別不同Cu2+脅迫濃度的光譜差異,而傳統監測方法及EMF值與葉片Cu2+含量均未達到此效果,驗證了EEMD-MA-FD模型的優越性。利用另一實驗數據對EEMD-MA-FD模型進行的普適性分析,表明EEMDMA-FD模型能夠在微弱的光譜差異中識別污染光譜信息,可作為玉米Cu2+污染監測的一種新方法。