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

孟 飛 崔 宇 付萍杰

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

0 引言

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

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

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

1 實驗設計與數據采集

1.1 實驗設計

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

1.2 玉米葉片光譜采集

采用美國SVC公司生產的SVC HR-1024I型高性能地物光譜儀測定Cu2+脅迫下玉米葉片光譜。光譜采集使用功率為50 W的鹵素燈光源和4°視場角探頭，探頭垂直于葉片表面40 cm，每盆玉米葉片光譜測量3次，由原始掃描光譜自動平均得到輸出結果，光譜反射系數經專用平面白板標準化處理，由此得到不同Cu2+濃度脅迫下的玉米葉片光譜。不同濃度梯度Cu2+脅迫下玉米葉片光譜如圖1所示。

1.3 葉片Cu2+濃度測定

對采集的玉米葉片樣品經過沖洗、干燥、粉碎等預處理后封裝在樣品袋中，樣品經高純硝酸、高氯酸消化后，用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光譜一階微分

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

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

(1)

式中λi——光譜波長

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

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

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

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

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

(2)

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

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

2.2 回歸模型

2.2.1多元線性回歸

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

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

(3)

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

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

β0、β1、…、β6——回歸系數

2.2.2偏最小二乘回歸

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

2.3 模型精度驗證

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

3 結果與分析

3.1 光譜一階微分分析

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

3.2 光譜的STFT分析

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結論

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

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