銅脅迫下玉米葉片光譜STFT分析與葉片銅離子濃度反演

孟 飛 崔 宇 付萍杰

(山東建筑大學(xué)測繪地理信息學(xué)院， 濟(jì)南 250101)

0 引言

土壤重金屬污染是影響農(nóng)作物正常生長的重要因素之一，其中銅(Cu)、鉛(Pb)、鎘(Cd)等重金屬污染較為突出[1]。土壤重金屬污染隱蔽性強、毒性大、治理難且周期長[2]。以Cu為例，常量時Cu是植物生長發(fā)育的必需微量元素[3]，而土壤Cu2+含量的超標(biāo)會破壞農(nóng)作物細(xì)胞組織結(jié)構(gòu)，引起植物代謝紊亂[4]。受重金屬脅迫的農(nóng)作物不但生長會受到抑制，而且土壤生理化參數(shù)也會發(fā)生改變[5]，這些變化都會反映在光譜曲線上。因此，利用光譜技術(shù)對農(nóng)作物重金屬污染脅迫等弱信息進(jìn)行提取與判別已成為高光譜遙感的研究熱點，快速有效地監(jiān)測農(nóng)作物重金屬污染對生態(tài)系統(tǒng)的安全與穩(wěn)定具有重要意義。研究利用葉片光譜判別玉米銅脅迫程度的最佳生長階段同樣具有重要意義[6]。

高光譜數(shù)據(jù)具有波段眾多、光譜分辨率高的優(yōu)勢，可實現(xiàn)農(nóng)作物快速、無損、實時的動態(tài)監(jiān)測[7-8]，也為監(jiān)測和防治土壤總金屬污染提供了技術(shù)保障[9]。劉厚田等[10]基于一、二階微分分析了土壤中過量銅對水稻葉片光譜反射特性的影響；李苑溪等[11]研究發(fā)現(xiàn)，銅脅迫下玉米葉片光譜紅邊波長與葉片中葉綠素含量存在一定的相關(guān)關(guān)系；高鵬等[12]基于銅脅迫下玉米葉片光譜的MRSVD特征構(gòu)建奇異特征參數(shù)，進(jìn)行了光譜奇異性分析，并采用逐步回歸分析法建立奇異特征參數(shù)與玉米葉片Cu2+含量之間的定量關(guān)系模型。光譜時頻分析在前人研究中同樣有了長足的發(fā)展。楊可明等[13]采用結(jié)合EMD算法的HHT時頻分析方法，對不同濃度Cu2+污染下土壤光譜進(jìn)行信息挖掘及診斷，并依據(jù)信息進(jìn)行Cu2+濃度的反演；姜雪芹等[14]將諧波分析與BP神經(jīng)網(wǎng)絡(luò)結(jié)合，利用土壤高光譜遙感的傳統(tǒng)處理數(shù)據(jù)進(jìn)行土壤含水率反演，取得較好的效果；李旭青等[15]基于小波變換和BP神經(jīng)網(wǎng)絡(luò)建立水稻冠層重金屬含量反演模型，該模型精度較高，具有一定的實用性。信號的STFT時頻分析方法是針對非平穩(wěn)信號非常有效的分析工具[16-17]。晁嬌[18]對不同類型復(fù)雜雷達(dá)信號經(jīng)STFT變換后的結(jié)果進(jìn)行統(tǒng)計，形成了基于短時傅里葉變換的雷達(dá)信號脈內(nèi)特征自動識別流程，這對電子偵察情報的獲取及應(yīng)用具有重要意義；同曉榮[19]提出一種基于STFT的聲信號盲分選算法，優(yōu)化了聲信號盲分選仿真系統(tǒng)的實時性，提高了分選精度；張博等[20]基于軌道高低不平順數(shù)據(jù)提出一種基于STFT時頻分析模型的高速鐵路路基局部變形識別方法；LIU等[21]提出了一種基于STFT的端點檢測方法，應(yīng)用于光纖干涉振動傳感系統(tǒng)；NARASIMHAN等[22]提出一種新的STFT控制算法，以減少基礎(chǔ)隔震建筑物的響應(yīng)；CHIKKERUR等[23]提出了一種基于STFT分析的指紋增強新方法，提高了傳統(tǒng)指紋識別精度。盡管高光譜遙感已廣泛應(yīng)用于農(nóng)業(yè)檢測領(lǐng)域，但傳統(tǒng)的農(nóng)作物重金屬污染監(jiān)測多數(shù)圍繞植被原始光譜獲取較為復(fù)雜的光譜特征參數(shù)，而重金屬脅迫下農(nóng)作物葉片光譜信號產(chǎn)生的變化微小，根據(jù)傳統(tǒng)監(jiān)測方法和常規(guī)的光譜特征參數(shù)很難區(qū)別光譜之間的微弱差異[24]，因此引入在雷達(dá)信號、聲音信號、光線信號、電信號等異常檢測中廣泛應(yīng)用的STFT時頻分析方法，研究重金屬脅迫下農(nóng)作物葉片光譜響應(yīng)。

傳統(tǒng)的STFT對較為平滑的信號反應(yīng)并不敏感，本文結(jié)合一階微分和STFT方法的優(yōu)勢，對農(nóng)作物原始平滑光譜做一階微分處理后進(jìn)行STFT時頻分析，研究重金屬脅迫下的農(nóng)作物葉片光譜信號能量隨波長和頻率聚集形態(tài)的眾多參數(shù)，通過分析揭示重金屬脅迫下農(nóng)作物葉片光譜的時頻能量變化規(guī)律，同時依據(jù)其峰值能量參數(shù)建立玉米葉片Cu2+濃度反演的多變量偏最小二乘回歸模型，以期實現(xiàn)基于光譜時頻分析法的農(nóng)作物重金屬脅迫的快速、簡單、有效監(jiān)測。

1 實驗設(shè)計與數(shù)據(jù)采集

1.1 實驗設(shè)計

選用含有Cu2+的CuSO4·5H2O溶液作為重金屬脅迫源，植物樣品選用“中糯1號”玉米進(jìn)行培養(yǎng)，實驗分別設(shè)置0、50、200、300、400、800、1 200 μg/g 7個濃度梯度，分別標(biāo)記為CK、Cu(50)、Cu(200)、Cu(300)、Cu(400)、Cu(800)和Cu(1200)，其中CK(0 μg/g設(shè)計濃度的植物樣品)為對照組，未添加CuSO4·5H2O，每個濃度梯度設(shè)置3組同期培養(yǎng)的平行實驗，組1為實驗組，組2、3為驗證組，共計21盆盆栽，作為研究用培養(yǎng)植物，用于獲取玉米葉片光譜。

1.2 玉米葉片光譜采集

采用美國SVC公司生產(chǎn)的SVC HR-1024I型高性能地物光譜儀測定Cu2+脅迫下玉米葉片光譜。光譜采集使用功率為50 W的鹵素?zé)艄庠春?°視場角探頭，探頭垂直于葉片表面40 cm，每盆玉米葉片光譜測量3次，由原始掃描光譜自動平均得到輸出結(jié)果，光譜反射系數(shù)經(jīng)專用平面白板標(biāo)準(zhǔn)化處理，由此得到不同Cu2+濃度脅迫下的玉米葉片光譜。不同濃度梯度Cu2+脅迫下玉米葉片光譜如圖1所示。

1.3 葉片Cu2+濃度測定

對采集的玉米葉片樣品經(jīng)過沖洗、干燥、粉碎等預(yù)處理后封裝在樣品袋中，樣品經(jīng)高純硝酸、高氯酸消化后，用WFX-120型原子吸收分光光度計測定Cu2+濃度。各濃度梯度Cu2+污染脅迫下3組玉米葉片樣本中Cu2+濃度見表1。

表1 不同濃度梯度Cu2+脅迫下玉米葉片Cu2+實測濃度Tab.1 Measured concentration of Cu2+ in maize leaves under different concentrations gradient Cu2+ stress μg/g

2 研究方法

2.1 光譜信號處理

2.1.1光譜一階微分

光譜信號的預(yù)處理包括歸一化、一階微分、統(tǒng)去除等，其中對光譜的一階微分處理是最為常用的方法。光譜數(shù)據(jù)進(jìn)行一階微分處理可以去除線性或接近線性的背景噪聲對目標(biāo)光譜的影響[25],也能提供比原光譜更高分辨率和更清晰的光譜輪廓變換[26]，且相比原始平滑光譜曲線，光譜一階微分曲線更加不平穩(wěn)，進(jìn)一步有利于信號的STFT變換處理。光譜一階微分計算式為

SD(λi)=(ρ(λi+1)-ρ(λi-1))/Δλ

(1)

式中λi——光譜波長

ρ(λi+1)、ρ(λi-1)——光譜反射率

SD(λi)——光譜一階微分

Δλ——波長λi+1與λi-1的波段間隔

2.1.2光譜短時傅里葉變換(STFT)

STFT變換的基本原理是把原始信號和窗函數(shù)進(jìn)行相乘并進(jìn)行一維的傅里葉變換，然后通過窗函數(shù)的滑動得到一系列頻譜函數(shù)，將頻譜函數(shù)隨時域依次展開可得到該信號的三維時頻圖，即信號隨時間和頻率變化的復(fù)數(shù)幅度。STFT變換的時間分辨率和頻率分辨率之間不能兼得，窗長越長，頻率分辨率越高，時間分辨率越差[27]。STFT常用窗函數(shù)有Hamming窗、Hanning窗、矩形窗、高斯窗等。STFT變換公式為

(2)

式中Z(u)——u時刻的原信號

g(u-t)——窗函數(shù)f——角頻率

2.2 回歸模型

2.2.1多元線性回歸

多元線性回歸模型(MLR)通常用于研究一個因變量與多個自變量的變化關(guān)系，若二者的關(guān)系可以用線性形式來描述，則可以建立多元線性模型進(jìn)行分析。本文分別以7組受不同濃度Cu2+脅迫下葉片光譜一階微分參數(shù)及光譜一階微分STFT分析參數(shù)為因變量，實測葉片Cu2+濃度為自變量，在SPSS中建立葉片Cu2+濃度反演的多元線性回歸模型。由于光譜一階微分及光譜一階微分STFT分析參數(shù)較多，SPSS軟件自動剔除對因變量無顯著影響的自變量，保留6個對因變量有顯著影響參數(shù)，分別建立六變量多元線性回歸模型，反演葉片Cu2+濃度。MLR模型為

y=β0+β1x1+…+β6x6

(3)

式中x1、x2、…、x6——非隨機變量，葉片光譜一階微分參數(shù)及光譜一階微分STFT能量幅值

y——因變量，模型預(yù)測Cu2+濃度

β0、β1、…、β6——回歸系數(shù)

2.2.2偏最小二乘回歸

偏最小二乘回歸(PLSR)是一種統(tǒng)計學(xué)方法，與主成分回歸有關(guān)，利用多因變量對多自變量進(jìn)行回歸建模。PLSR模型通過投影分別將頻譜和相關(guān)響應(yīng)變量投影到一個二維空間，來尋找一個線性回歸模型，從而降低維數(shù)并避免噪聲[28]。相比傳統(tǒng)回歸分析，在自變量多、樣本數(shù)量較少的情況下，PLSR有更明顯的優(yōu)勢[29]。本文分別以7組受不同濃度Cu2+脅迫下葉片光譜一階微分參數(shù)及光譜STFT時域特征波段上隨頻域變化的STFT時頻能量幅值(每個特征波段包含65個幅值)為因變量，以實測葉片Cu2+濃度為自變量分別建立多變量PLSR模型，反演葉片Cu2+濃度。PLSR模型流程圖如圖2所示(以STFT能量參數(shù)為例)。

2.3 模型精度驗證

采用決定系數(shù)R2評價模型預(yù)測值與實測值的擬合精度。決定系數(shù)R2是回歸平方和占總誤差平方和的比例，反映回歸直線的擬合程度，其取值范圍介于0～1之間。R2越趨近于1，說明回歸方程擬合精度越高；R2越趨近于0，說明回歸方程擬合精度越差。同時采用均方根誤差RMSE評價模型精度，反映模型預(yù)測值與實測值之間的離散程度，RMSE越小，說明模型預(yù)測值與實測值之間離散程度越小，模型精度越高。

3 結(jié)果與分析

3.1 光譜一階微分分析

用式(1)對獲取的波段340～1 000 nm玉米葉片光譜信號進(jìn)行一階微分處理，發(fā)現(xiàn)不同濃度Cu2+脅迫下玉米葉片光譜一階微分處理后在波段480～670 nm、700～750 nm處產(chǎn)生較為明顯的峰值及低谷變化，原因是該波段范圍是較易受重金屬脅迫影響的玉米葉片光譜響應(yīng)特征的“紫谷-綠峰-紅谷”所對應(yīng)的波段范圍[30]。通過信號的STFT變換，能夠在該明顯“峰值-低谷”變化范圍獲取其對應(yīng)時頻能量的變化規(guī)律，并針對相應(yīng)規(guī)律及參數(shù)做進(jìn)一步研究。不同濃度梯度Cu2+脅迫下玉米葉片光譜一階微分處理結(jié)果如圖3所示。

3.2 光譜的STFT分析

以不同濃度梯度Cu2+脅迫下波段340～1 000 nm玉米葉片光譜的一階微分?jǐn)?shù)據(jù)為原始信號，為突出光譜一階微分STFT時頻能量在一階微分光譜“峰值-低谷”位置的能量值及其變化情況，設(shè)置采樣間隔為1，F(xiàn)FT長度為128，窗口長度為32，采樣重疊度為31，窗口選用Hamming窗進(jìn)行STFT變換，獲取不同濃度梯度Cu2+脅迫下波段340～1 000 nm玉米葉片光譜響應(yīng)的一階微分處理數(shù)據(jù)的時頻圖，將時域轉(zhuǎn)換為虛擬光譜，如圖4所示。能量振幅峰值所處波段隨Cu2+濃度梯度變化情況見表2。

由圖4可看出，隨著Cu2+濃度的不斷增加，能量振幅峰值呈先降低后升高趨勢。結(jié)合表2，無重金屬Cu2+脅迫的正常植株的光譜一階微分STFT峰值所處波長為722 nm附近。隨著Cu2+濃度不斷升高，實驗組與驗證組的光譜一階微分STFT時頻能量峰值均在樣品Cu2+濃度設(shè)置為400 μg/g時為最小值，峰值能量所處波段均出現(xiàn)在690～730 nm之間，且沿波段軸向左，即短波方向遷移，僅驗證組組2與組3在Cu2+濃度800 μg/g與1 200 μg/g時振幅與所處波段出現(xiàn)極小反常情況，但未影響整體規(guī)律。紅邊(REP)是綠色植物葉子光譜曲線在波段680～740 nm之間變化率最快的點，也是一階導(dǎo)數(shù)光譜在該區(qū)間內(nèi)的拐點[31]，而重金屬Cu脅迫程度增加會造成農(nóng)作物紅邊“藍(lán)移”程度增強[32]，因此玉米葉片光譜受重金屬脅迫影響而產(chǎn)生的紅邊“藍(lán)移”是其STFT時頻峰值能量遷移的重要原因。

3.3 玉米葉片重金屬Cu2+濃度反演

3.3.1基于光譜一階微分變量的MLR及PLSR模型反演

選取17組不同Cu2+濃度梯度脅迫下玉米葉片光譜一階微分處理參數(shù)為自變量，葉片實測Cu2+濃度為因變量，在SPSS軟件中做多元回歸線性擬合，建立Cu2+濃度反演的MLR模型。光譜一階微分前共661個波段，一階微分后得659個參數(shù)，將參數(shù)編號為Fod1、Fod2、…、Fod659，軟件模型參數(shù)選取及對應(yīng)系數(shù)如表3所示，模型R2=0.98(P<0.01)。

分別選取組2與組3光譜一階微分參數(shù)對應(yīng)參數(shù)代入模型進(jìn)行Cu2+濃度反演，Cu2+濃度反演值與實測值對比見圖5。由圖5可得，光譜一階微分參數(shù)建立的線性預(yù)測模型精度低，組2實測值與預(yù)測值R2=0.114，RMSE為5.696 μg/g；組3R2=0.120 1，RMSE為4.159 μg/g，且組2反演結(jié)果存在一個負(fù)值，表明簡單的光譜一階微分參數(shù)建立的六變量參數(shù)MLR模型盡管本身精度較高，但存在一定的局限，不能實現(xiàn)對驗證組玉米葉片Cu2+濃度的準(zhǔn)確預(yù)測。

選取17組不同Cu2+濃度梯度脅迫下玉米葉片光譜一階微分全部參數(shù)為自變量，實測Cu2+濃度為因變量，在Matlab中建立一階微分參數(shù)PLSR Cu2+濃度反演模型。模型自變量系數(shù)見圖6a，精度驗證見圖6b，模型精度R2=1(P<0.01)，由于模型本身精度過高，存在過擬合現(xiàn)象，因此僅將此模型用于對比，而不能對其進(jìn)行參考。

表3 一階微分參數(shù)反演MLR模型選取編號及對應(yīng)系數(shù)Tab.3 First-order differential parameter inversion MLR model selection number and corresponding coefficient

依據(jù)組1光譜一階微分參數(shù)建立的PLSR模型，對組2與組3光譜一階微分參數(shù)進(jìn)行Cu2+濃度反演，模型反演值與實測值擬合見圖7。由圖7可見，組2與組3反演精度R2分別為0.619 2(P<0.05)和0.451 2，RMSE分別為6.05、5.429 μg/g，盡管模型對驗證組反演精度有所提高，但總體離散程度較大，因此光譜一階微分STFT參數(shù)建立的PLSR模型雖發(fā)揮了自變量多的優(yōu)勢，且模型本身精度極高，但其反演結(jié)果僅能用作對比，而不予對其參考。

3.3.2基于光譜一階微分STFT參數(shù)的MLR及PLSR模型反演

同樣選取17組不同Cu2+濃度梯度脅迫下玉米葉片光譜一階微分STFT能量峰值所在時域波段上的能量參數(shù)為自變量，實測葉片Cu2+濃度為因變量，在SPSS軟件中做多元回歸線性擬合，建立Cu2+濃度反演模型。光譜一階微分STFT處理后在頻域上得65個參數(shù)，將參數(shù)編號為St1、St2、…、St65，軟件模型參數(shù)選取及對應(yīng)系數(shù)見表4，模型精度R2=0.97(P<0.01)，精度較高。忽略組2、3真實峰值所處波段，仍以組1建模波段為參數(shù)選取波段，將組2與組3光譜一階微分STFT數(shù)據(jù)對應(yīng)參數(shù)代入模型進(jìn)行葉片Cu2+濃度反演，葉片Cu2+濃度模型反演值與實測值對比見圖8。由圖8可得，光譜一階微分STFT參數(shù)建立的MLR預(yù)測模型組2實測值與預(yù)測值R2為0.509 5，RMSE為6.152 μg/g，組3R2為0.275 1，RMSE為4.921 μg/g，與光譜一階微分參數(shù)建立的MLR模型相比，精度略有提高，但組3反演效果仍不理想，僅不存在負(fù)值反演結(jié)果。

為進(jìn)一步提高STFT參數(shù)模型反演精度，發(fā)揮PLSR模型的優(yōu)勢，選取組1 Cu2+濃度梯度脅迫下玉米葉片一階微分處理光譜的STFT變換虛擬光譜波段上特定波長隨頻域變化的65個全部能量幅值用作7組自變量，7個實測葉片Cu2+濃度用作因變量，在Matlab中建立PLSR Cu2+濃度預(yù)測模型。建立的PLSR線性回歸模型自變量系數(shù)見圖9a，模型精度驗證見圖9b。忽略驗證組實際振幅峰值波段，仍以實驗組建模參數(shù)波段為特征波段，依次選取組2與組3光譜響應(yīng)一階微分?jǐn)?shù)據(jù)的STFT變換結(jié)果的對應(yīng)波段全部能量幅值作為驗證參數(shù)，代入PLSR模型進(jìn)行葉片Cu2+濃度含量反演，并將反演結(jié)果與驗證組實測Cu2+濃度進(jìn)行線性擬合，擬合結(jié)果如圖10所示。組2與組3R2分別為0.880 6與0.733 1(P<0.01)，RMSE分別為1.563、2.619 μg/g，相比SPSS軟件建立的6自變量MLR模型，STFT參數(shù)的PLSR多自變量模型反演精度進(jìn)一步提高，且反演結(jié)果離散程度較其他參數(shù)和模型的反演結(jié)果，有了大幅降低。

表4 一階微分STFT參數(shù)反演線性模型選取編號及對應(yīng)系數(shù)Tab.4 First-order differential STFT parameter inversion linear model selection number and corresponding coefficient

4 結(jié)論

(1)隨著Cu2+脅迫濃度升高，受農(nóng)作物“紅邊藍(lán)移”現(xiàn)象的影響，光譜一階微分STFT能量峰值參數(shù)呈先降低、后升高趨勢，且峰值參數(shù)頻譜隨Cu2+濃度升高沿波段軸不斷向短波方向遷移。峰值參數(shù)頻譜的遷移狀況與峰值參數(shù)可以作為判別玉米葉片光譜重金屬脅迫程度的標(biāo)志。

(2)相比一階微分參數(shù)、一階微分STFT分析參數(shù)的6自變量MLR模型，以及一階微分參數(shù)PLSR模型，以峰值參數(shù)頻譜的全部能量幅值為自變量、以實測玉米葉片Cu2+濃度為因變量建立的多自變量光譜一階微分STFT參數(shù)的PLSR模型Cu2+濃度預(yù)測精度最高，且在對驗證組兩組參數(shù)進(jìn)行反演時，預(yù)測值與實測值R2分別為0.880 6和0.733 1(P<0.01)，RMSE分別為1.563、2.619 μg/g。由此可見，若驗證樣本組數(shù)量足夠多，且本文建立的模型能夠利用相同方式獲取參數(shù)對玉米葉片Cu2+濃度進(jìn)行準(zhǔn)確預(yù)測，則表明本文通過光譜STFT時頻分析方法獲取參數(shù)建立的Cu2+濃度預(yù)測模型存在一定的普適性。