張建紅，楊可明，韓倩倩，李艷茹，高 偉

中國礦業(yè)大學(xué)(北京)地球科學(xué)與測繪工程學(xué)院，北京 100083

引 言

近些年來環(huán)境污染，特別是重金屬污染在工業(yè)化和城市化快速發(fā)展地區(qū)尤為嚴(yán)重，因此重金屬污染監(jiān)測與治理技術(shù)等已成為當(dāng)今環(huán)境保護(hù)領(lǐng)域的研究熱點，其中農(nóng)業(yè)重金屬污染更是備受關(guān)注[1]。農(nóng)業(yè)重金屬污染主要來源于農(nóng)藥和近農(nóng)田區(qū)域采礦活動等[2]。農(nóng)田土壤受重金屬污染程度不易被察覺，而重金屬離子可以被農(nóng)作物吸收并富集，進(jìn)而能造成糧食減產(chǎn)、毒害食物鏈、威脅食品安全、危害人體健康。

如何實時快捷有效識別重金屬污染、更加準(zhǔn)確甄別重金屬污染的類別與程度是當(dāng)前亟待研究與突破的關(guān)鍵技術(shù)。傳統(tǒng)的重金屬污染監(jiān)測方法通常是基于野外取樣后，在室內(nèi)采用無焰原子吸收分光光度法(FAAS)、原子吸收光譜法(AAS)等進(jìn)行測定，雖然檢測精度高，但也具有耗費大量人力、物力和財力，損害被測對象，所檢測儀器昂貴、維護(hù)成本高、操作程序復(fù)雜和檢測時間長等缺點[3]；而今高光譜遙感技術(shù)具有圖譜合一、光譜連續(xù)，高效、無損、低成本等優(yōu)勢，正在被廣泛應(yīng)用于重金屬污染甄別與監(jiān)測的探索研究中[4]。Liu等[5]通過提取植被紅邊光譜曲線，構(gòu)造歸一化紅邊指數(shù)和紅邊葉綠素指數(shù)，用于識別植被受重金屬鎘(Cd)的脅迫；喬曉英等[6]以礦區(qū)周邊玉米、苦菜等作為研究對象，構(gòu)建植物中重金屬Cd、鉛(Pb)質(zhì)量比與植物光譜紅邊、藍(lán)谷等特征參數(shù)的關(guān)系模型，可間接反演植物中重金屬的含量；付萍杰等[7]對葉片光譜進(jìn)行經(jīng)驗?zāi)B(tài)分解并結(jié)合自相關(guān)函數(shù)、多重分形理論構(gòu)建Cu、Pb含量反演的多元線性回歸模型，成功識別了玉米葉片中重金屬類別及含量；朱葉青等[8]通過獲取不同生育期、不同Cu污染下的植被葉片光譜信息，采用7個特征波段結(jié)合光譜角方法研究Cu污染下的葉片光譜特征，表明Cu污染葉片與健康葉片的光譜存在明顯差異，同時證明了Cu污染后的葉片內(nèi)部結(jié)構(gòu)更加紊亂。

作為植被生長所必需的微量元素Cu，也屬于一種重金屬，當(dāng)濃度較低時有利于植被生長，濃度過高就會抑制植物的光合作用，對植被生長產(chǎn)生脅迫作用，輕則引起植物代謝紊亂，生長受阻，重則導(dǎo)致植物死亡。已有報道當(dāng)植被受到重金屬污染時，多數(shù)會長勢矮小，葉片失去綠色等癥狀[9]。本工作提出固有波長尺度分解(intrinsic wavelength-scale decomposition，IWD)的方法并運用于甄別農(nóng)作物重金屬污染的高光譜信息處理與分析中。依據(jù)受Cu污染的玉米葉片光譜數(shù)據(jù)，基于Hankel矩陣構(gòu)建、奇異值分解(singular value decomposition，SVD)等建立玉米葉片Cu污染信息預(yù)測的IWD-Hankel-SVD模型，該模型可有效預(yù)測玉米葉片中銅離子(Cu2+)含量，能為植物Cu污染信息預(yù)測、甄別提供一種新方法。

1 理論與算法

1.1 固有波長尺度分解

固有波長尺度分解(IWD)引自于固有時間尺度分解(intrinsic time-scale decomposition，ITD)。ITD是由Mark G Frei與Ivan Osorio首次提出[10]，相對于經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition，EMD)和小波變換等，ITD是一種新的時頻分析方法[11]。經(jīng)過ITD分解后的信號分量具有相對較為完整的時頻信息，能夠反映信號頻率變化等優(yōu)點。本工作用波長(wavelength，W)替代ITD中的時間(time，T)，探索形成光譜信號的IWD處理技術(shù)，并引入到高光譜數(shù)據(jù)的信息處理與分析中。IWD可以將波長序列(w)自適應(yīng)分解為幾個固有旋轉(zhuǎn)分量(proper rotation component，PRC)和一個趨勢分量(r)，每個旋轉(zhuǎn)分量只需要通過迭代一次就可以得到。對于一條準(zhǔn)備分解的光譜曲線Xw，有極值點Xk，其基線提取算子定義為L，可將基線分量Lw與準(zhǔn)備分解的光譜曲線Xw分開，再定義一個固有旋轉(zhuǎn)分離算子為H，于是有

Xw=LXw+(1-L)Xw=Lw+Hw

(1)

其中，Lw=LXw為基線分量，Hw=(1-L)Hw=HXw為PRC分量。

對于給定的光譜曲線，IWD的處理步聚如下：

(1)提取光譜曲線Xw上的極值點Xk(k=1,2,…,N，N表示極值點的個數(shù))，令τk是Xw局部極值點所在的波長集合，定義τ0=0作為Xw的起始波長，并且使Xw和Lw在τk處的值分別為Xτk和Lτk，為了表示方便，將Xτk和Lτk分別定義為Xk和Lk。

假設(shè)Lw和Hw被定義在波段[τ0,τk]上，Xw被定義在波段[τ0,τk+2]上，并且，假設(shè)Lw是Xw在區(qū)間(τk,τk+1]的仿射線性逼近，即Lw=mXw+nw，w∈(τk,τk+1]。因為Lw必須滿足局部極值點的上述邊界條件，所以Lw被表示為以下形式

(2)

假設(shè)光譜曲線的基本趨勢足夠平滑，則可以忽略各局部極值點間的變化，于是可以得出

(3)

因為Lk+1=Xk+1，于是有

Lk+1=αXk+1+(1+α)Xk+1=

(4)

式(4)中：α為線性縮放因子，可用于調(diào)整提取到的PRC幅度，α∈[0,1]，一般情況下，經(jīng)驗值取值0.5。

(2)將式(2)計算得出的基線分量Lw=LXw，按照式(1)計算提取固有旋轉(zhuǎn)分量Hw，即PRCi。

(3)將基線分量Lw作為一個輸入光譜再進(jìn)行下一次分解，重復(fù)第(1)和(2)步的分解，當(dāng)基線分量Lw變?yōu)閱握{(diào)或小于某個預(yù)設(shè)值的趨勢分量(r)時，分解終止。

IWD的總分解過程為

Xw=Hw+Lw=HXw+LXw=HXw+(H+L)LXw=

(5)

經(jīng)過IWD分解后，原始輸入的光譜曲線Xw即被分解為多個固有旋轉(zhuǎn)分量(PRCi)和一個單調(diào)趨勢分量(r)。該方法計算復(fù)雜度低，分解精度高，魯棒性強(qiáng)，有效的避免了波形疊加、模態(tài)混疊現(xiàn)象。

1.2 Hankel矩陣

漢克爾(Hankel)矩陣具有很多優(yōu)良的性能和特殊的性質(zhì)，因此，作為一種優(yōu)秀的數(shù)學(xué)工具被廣泛應(yīng)用于諸多領(lǐng)域的數(shù)值模擬計算中。Hankel矩陣的每一條逆對角線上的元素都相等?？蛇x擇經(jīng)IWD處理后較好保留了原始輸入光譜信息的PRC，構(gòu)造Hankel矩陣Za

(6)

式(6)中，1≤a≤J，且1

1.3 奇異值分解與奇異熵

奇異值分解(SVD)即對于任意一個維數(shù)的實矩陣A∈Rm×n，都可分解為三個矩陣的乘積，即

A=USVT

(7)

式(7)中，U為單位正交矩陣，即其列向量是單位向量并且相互正交，是奇異值對應(yīng)的特征向量[13]；V為單位正交矩陣，其行向量是單位向量并且相互正交，也是奇異值對應(yīng)的特征向量；S為對角矩陣，S=[diag(σ1,σ2,σ3,…,σm),O]，其中O表示零矩陣，它的對角線矩陣為奇異值。

奇異熵是用來度量奇異值所對應(yīng)用的信號分量信息量的多少，奇異熵的計算公式是

(8)

式(8)中，k為奇異熵所對應(yīng)的階次；ΔEi為奇異熵在階次i時得到的增量，其計算公式為

(9)

1.4 IWD-Hankel-SVD模型

基于IWD，Hankel矩陣和SVD可構(gòu)建IWD-Hankel-SVD模型，用于光譜的信息分析與特征參量反演等。IWD-Hankel-SVD模型的植物重金屬污染下葉片光譜特征信息提取步驟為：

(1)IWD分解。對需要甄別重金屬污染的植物光譜數(shù)據(jù)進(jìn)行IWD分解，得到PRC1，PRC2，…，以及r若干分量。

(2)分量篩選。對經(jīng)過IWD分解后得到的各PRC分量與原始輸入光譜進(jìn)行對比，選擇能夠較好保留原始輸入光譜信息的PRCi分量。

(3)特征波段選擇。在PRCi分量中選取與原始輸入光譜相似的波段作為特征波段進(jìn)行下一步的處理。

(4)構(gòu)造Hankel矩陣。利用從第(2)步得到的PRCi分量選擇的特征波段，構(gòu)造Hankel矩陣Z。

(5)SVD分解。對第(4)步構(gòu)造的Hankel矩陣進(jìn)行SVD分解，得到奇異值矩陣。

(6)求取奇異熵。通過第(5)步得到的奇異值計算對應(yīng)的奇異熵，作為表征重金屬污染下光譜信息分量的特征。

2 實驗部分

(1)實驗準(zhǔn)備。實驗采用盆栽的玉米植株土培法，玉米種子選用“中糯301”。重金屬Cu2+脅迫試劑選用純度較高、干擾雜質(zhì)較少的分析純級別CuSO4·5H2O。光譜數(shù)據(jù)獲取采用SVC HR-1024i便攜式地物光譜儀，其光譜測量范圍為350～2 500 nm，其中350～1 000 nm光譜分辨率為3.0 nm，1 000～1 900 nm光譜分辨率為9.5 nm，1 900～2 500 nm光譜分辨率為6.5 nm。葉綠素濃度相對值的測定采用SPAD-502手持便攜式葉綠素測定儀。

(2)植株培養(yǎng)?！爸信?01”玉米種植時采用底部有滲水孔的培土花盆，將分析純級別的CuSO4·5H2O配置成相應(yīng)梯度的Cu2+溶液拌入培土作為所種植玉米生長的脅迫試劑。培植實驗共設(shè)置三組平行實驗，每組設(shè)7個Cu2+脅迫梯度，即空白對照組(ck，不添加脅迫試劑)和100，200，400，600，800和1 000 μg·g-1等6級Cu2+脅迫梯度，分別標(biāo)注為ck(0)，Cu(100)，Cu(200)，Cu(400)，Cu(600)，Cu(800)和Cu(1 000)。為了提高玉米培養(yǎng)效率，事先對“中糯301”玉米種子進(jìn)行催芽處理。植株培養(yǎng)期間保持水量充沛，空氣流通。2個月后培植結(jié)束，測定玉米葉片光譜信息、葉片中Cu2+的含量和葉綠素濃度相對值。

(3)光譜數(shù)據(jù)采集。在室內(nèi)采用SVC HR-1024i光譜儀測定實驗盆栽玉米在不同濃度Cu2+脅迫下的葉片光譜數(shù)據(jù)，選用50W鹵素?zé)糇鳛閱我还庠础９庾V采集時分別測量每個平行組各脅迫濃度下玉米老(Old)、中(Middle)、新(New)葉片的光譜3次，每個脅迫濃度獲取9組數(shù)據(jù)，共計63組光譜數(shù)據(jù)。最終獲取的光譜數(shù)據(jù)是經(jīng)由去除異常值，求平均得到。

(4)葉綠素濃度相對值測定。同步利用SPAD-502葉綠素儀分別測定采集光譜的葉片葉綠素濃度(重復(fù)3次)，最終通過求平均得到了各Cu2+脅迫濃度下玉米葉片葉綠素濃度相對值。

(5)Cu2+含量的測定。在各類被測葉片光譜和葉綠素的數(shù)據(jù)采集結(jié)束后就即時進(jìn)行葉片樣本裁剪和保存處理，為葉片中Cu2+含量測定作準(zhǔn)備。最終測定的葉片中Cu2+含量是采用電感耦合等離子發(fā)射光譜儀(ICP-OES)進(jìn)行分析得到。

脅迫梯度以及對應(yīng)的玉米葉片葉綠素濃度相對值、Cu2+含量見表1所示。由表1可見，隨著Cu2+脅迫增大，葉綠素濃度相對值在逐漸減小，葉片中Cu2+含量在逐漸增大。各脅迫梯度的葉片葉綠素濃度相對值與Cu2+含量相關(guān)系數(shù)為-0.953，p值為0.001<0.01，可見兩者呈負(fù)相關(guān)，負(fù)相關(guān)性高度顯著，說明在重金屬Cu2+脅迫下，玉米葉片中葉綠素含量的積累受到影響[14]。

表1 不同濃度Cu2+脅迫下玉米葉片中葉綠素濃度相對值及Cu2+含量Table 1 Relative value of chlorophyll concentration and content of Cu2+ in corn leaves under different concentrations of Cu2+

3 結(jié)果與討論

3.1 IWD處理結(jié)果與分析

通過SVC HR-1024i光譜儀測定的不同Cu2+脅迫梯度下玉米葉片光譜曲線如圖1所示，受不同濃度的Cu2+脅迫，玉米葉片光譜曲線與對照組ck(0)的健康玉米葉片光譜曲線有所不同，但是不能準(zhǔn)確的甄別出受重金屬Cu2+脅迫的變化規(guī)律。為此，引入IWD對不同濃度Cu2+脅迫下的玉米葉片光譜曲線進(jìn)行變換，然后構(gòu)造Hankel矩陣，進(jìn)行SVD處理，最后求得奇異熵，通過線性擬合建立的IWD-Hankel-SVD模型監(jiān)測確定重金屬Cu2+脅迫下玉米葉片的污染程度。

采用IWD方法分解對照組ck(0)的健康葉片原始光譜，得到PRC1和PRC2的2個固有旋轉(zhuǎn)分量和1個趨勢分量r，如圖2所示，其中并非每一個分量都包含有豐富的光譜特征信息，有些分量對于光譜特征信息的提取貢獻(xiàn)很小。通過與原始輸入光譜信息比較，發(fā)現(xiàn)PRC1分量較好的保留了原始輸入光譜信息，并且對原始光譜中部分波段的反射率有所放大。同樣，通過對Cu(100)，Cu(200)，Cu(400)，Cu(600)，Cu(800)和Cu(1000)脅迫梯度下原始葉片光譜進(jìn)行IWD分解，發(fā)現(xiàn)不同脅迫濃度下的固有旋轉(zhuǎn)分量中PRC1分量均較好地保留了原始輸入光譜信息，如圖3(a—f)所示。因此，根據(jù)PRC1分量的光譜保留信息，從中選擇與原始輸入光譜相似的553～680，681～780，1 266～1 429，1 430～1 631，1 836～1 913和1 914～2 111 nm波段區(qū)間作為有效特征波段。

圖2 對照組ck(0)玉米葉片光譜及其IWD處理結(jié)果Fig.2 Spectrum and the results processed by IWD of corn leaves in ck(0) control group

圖3 不同濃度 Cu2+脅迫下玉米葉片光譜IWD處理的PRC1分量Fig.3 Spectral PRC1 components of the corn leaves stressed by the different concentrations of Cu2+ based on the IWD processing

3.2 單變量的IWD-Hankel-SVD模型應(yīng)用

通過對特征波段553～680，681～780，1 266～1 429，1 430～1 631，1 836～1 913和1 914～2 111 nm對應(yīng)的PRC1分量分別構(gòu)造Hankel矩陣，并對該矩陣進(jìn)行SVD分解，實現(xiàn)IWD-Hankel-SVD模型的應(yīng)用，最后求取模型的奇異熵用于分析監(jiān)測效果。各特征波段所對應(yīng)的奇異熵與玉米葉片中Cu2+含量的相關(guān)關(guān)系見表2。由表2可見，通過特征波段1 266～1 429和1 836～1 913 nm計算求得的奇異熵與玉米葉片中Cu2+含量的相關(guān)系數(shù)分別為0.950和0.913，說明經(jīng)過IWD分解后的波段1 266～1 429和1 836～1 913 nm對應(yīng)的PRC1分量對Cu污染較為敏感。選取相關(guān)系數(shù)大于0.9的特征波段計算得到的奇異熵分別與玉米葉片中Cu2+含量進(jìn)行線性擬合，擬合結(jié)果如圖4、圖5所示，可見：(1)根據(jù)特征波段1 266～1 429 nm計算得到的奇異熵與玉米葉片中Cu2+含量線性擬合，得到的模型決定系數(shù)R2為0.903，p值小于0.01，相關(guān)性高度顯著；(2)根據(jù)特征波段1 836～1 913 nm計算得到的奇異熵與葉片中Cu2+含量線性擬合，得到的模型決定系數(shù)R2為0.834，p值小于0.01，相關(guān)性高度顯著。因此得出，可通過IWD-Hankel-SVD模型來預(yù)測玉米植株受重金屬Cu污染的程度。

圖4 特征波段1 266～1 429 nm的奇異熵與玉米葉片中Cu2+含量變化關(guān)系Fig.4 Relationship between the singular entropy of 1 266～1 429 nm characteristic band and the content of Cu2+ in corn leaves

圖5 特征波段1 836～1 913 nm奇異熵與玉米葉片中Cu2+含量變化關(guān)系Fig.5 Change relationship between the singular entropy of 1 836～1 913 nm characteristic band and the content of Cu2+ in corn leaves

表2 特征波段的奇異熵與玉米葉片Cu2+含量的相關(guān)關(guān)系Table 2 The correlation between the singular entropy of characteristic bands and the content of Cu2+ in corn leaves

為了驗證依據(jù)所提取特征波段以及IWD-Hankel-SVD模型預(yù)測植物重金屬污染程度的優(yōu)越性，選擇傳統(tǒng)的藍(lán)邊最大值、紅肩最大值、綠峰高度等光譜特征參數(shù)(見表3)監(jiān)測方法以及與采用連續(xù)投影變換(SPA)選擇的特征波段如圖6、圖7所示，進(jìn)行應(yīng)用結(jié)果對比分析。通過藍(lán)邊最大值、紅肩最大值、綠峰高度的監(jiān)測結(jié)果與玉米葉片中Cu2+含量的線性擬合，得到擬合效果如圖8、圖9、圖10所示，通過連續(xù)投影變換(SPA)選擇的特征波長與玉米葉片中Cu2+含量的線性關(guān)系見表4，對比結(jié)果表明，IWD-Hankel-SVD模型的監(jiān)測能力明顯要優(yōu)于其他傳統(tǒng)光譜特征等參數(shù)所建立的模型。

表3 玉米葉片的光譜特征參數(shù)名稱與定義Table 3 Names and definitions of spectral characteristic parameters of corn leaves

表4 SPA選擇的特征波長與玉米葉片中Cu2+含量的預(yù)測表達(dá)式Table 4 The prediction expression of Cu2+ content in corn leaves and characteristic wavelengths selected by SPA

圖6 連續(xù)投影變換選擇特征變量的個數(shù)Fig.6 Number of characteristic bands selected by SPA

圖7 連續(xù)投影變換提取特征波長Fig.7 Characteristic bands extraction by SPA

圖8 藍(lán)邊最大值方法監(jiān)測結(jié)果與葉片中Cu2+含量的變化關(guān)系Fig.8 Relation between monitoring results based on the maximum of blue edge and Cu2+ content in corn leaves

圖9 紅肩最大值方法監(jiān)測結(jié)果與葉片中Cu2+含量的變化關(guān)系Fig.9 Relation between monitoring results based on the maximum of red shoulder and Cu2+ content in corn leaves

圖10 綠峰高度方法監(jiān)測結(jié)果與玉米葉片中Cu2+含量的變化關(guān)系Fig.10 Relation between monitoring results based on the green peak height and Cu2+ content in corn leaves

3.3 多變量的IWD-Hankel-SVD模型應(yīng)用

植被中的葉綠素含量是植物生長健康狀況的重要生物指標(biāo)，決定著綠色植物的光合作用效率。隨著Cu2+脅迫濃度的增大，玉米葉片中的Cu2+含量在逐漸增大，而葉綠濃度相對值在逐漸減少。這是因為重金屬脅迫對植物的生長具有負(fù)向影響，抑制植物進(jìn)行光合作用，造成葉綠素水解增加，含量降低，改變細(xì)胞結(jié)構(gòu)。植物葉片中葉綠素含量的變化會引起可見光和近紅外波段的反射率變化[15]。選擇數(shù)據(jù)采集時測得的葉綠素濃度相對值、特征波段1 266～1 429和1 836～1 913 nm分別計算得到的模型奇異熵作為參數(shù)，結(jié)合偏最小二乘回歸分析預(yù)測玉米葉片受重金屬污染下的Cu2+含量。經(jīng)驗證可得預(yù)測值與實測值的相關(guān)性更強(qiáng)，如圖11所示，其決定系數(shù)R2為0.947 6。與單變量的模型應(yīng)用結(jié)果相比，多變量的IWD-Hankel-SVD模型具有較強(qiáng)的魯棒性、穩(wěn)健性。

圖11 玉米葉片中Cu2+含量的多變量模型預(yù)測精度Fig.11 Accuracy of the multivariable model on predicting Cu2+ content in corn leaves

4 結(jié) 論

重金屬Cu2+脅迫引起的玉米葉片光譜特征變化難以依據(jù)光譜曲線形態(tài)甄別其污染程度。提出了IWD概念和方法，并通過構(gòu)造IWD-Hankel-SVD模型對玉米葉片高光譜數(shù)據(jù)的Cu污染信息進(jìn)行反演，甄別和預(yù)測玉米葉片受Cu污染的程度，認(rèn)為：(1)對Cu污染的玉米葉片光譜數(shù)據(jù)進(jìn)行IWD處理，分解提取到的固有旋轉(zhuǎn)分量PRC1能較好地保留了原始輸入光譜信息，有效避免了模態(tài)混疊現(xiàn)象；(2)通過特征波段1 266～1 429和1 836～1 913 nm求得的模型奇異熵與玉米葉片中的Cu2+含量存在很強(qiáng)的相關(guān)性，模型預(yù)測的決定系數(shù)R2均高達(dá)0.9左右。通過一些傳統(tǒng)預(yù)測方法和由SPA選擇的特征波段的預(yù)測結(jié)果比較，驗證了依據(jù)所提取特征波段這單一變量形式的IWD-Hankel-SVD模型在玉米葉片Cu污染信息和污染程度方面具有一定的優(yōu)越性和較好的解釋能力；(3)選擇玉米葉片中葉綠素濃度相對值以及葉片光譜特征波段1 266～1 429和1 836～1 913 nm分別計算得到的模型奇異熵作為參數(shù)，形成多變量的模型預(yù)測技術(shù)，可得到葉片中Cu含量的預(yù)測值與實測值之間相關(guān)性更強(qiáng)，其決定系數(shù)R2達(dá)到0.9476。與單變量模型相比，多變量模型具有較強(qiáng)的魯棒性和穩(wěn)健性。