冷杉葉片氮含量高光譜反演技術研究

宋雪蓮,王志偉,張文，2，張威，丁磊磊，柳嘉佳，阮璽睿,王普昶

(1.貴州省農業科學院草業研究所，貴州 貴陽 550006；2.貴州陽光草業科技有限責任公司，貴州 貴陽 550006；3.貴州省水利水電勘測設計研究院有限公司,貴州 貴陽 550002；4.貴州省農業科學院,貴州 貴陽 550006)

氮是植物生長的基本營養元素，是陸地生態系統植物生長的主要限制因子，與植物的光合作用和細胞生長分類等重要生理活動相關[1]。缺氮會影響植物的光合作用能力，導致葉片生理及形態結構的相應變化，因而引起葉片光譜反射特性的變化[2]。傳統的化學方法能夠精確地監測出植物氮含量，但具有高損耗、復雜、時滯性等缺點，不能進行全面、大范圍、快速的營養診斷。高光譜光譜分辨率高，能夠獲得連續的光譜信息，已成為監測植物氮含量的有效手段。

目前有不少學者采用敏感波段篩選和植被指數進行高光譜植物氮含量檢測的相關研究。朱艷等[4]研究表明，單波段光譜在610 nm和680 nm處的水稻冠層反射率與葉片氮含量具有較好的相關性，提出采用回歸系數來提高水稻葉片氮含量估測的準確性。劉冰峰等[5]指出720 nm處的反射光譜以及550，720 nm處的一階光譜導數是夏玉米氮含量的敏感參數。Tarpley等[6]研究認為，利用紅邊位置和短波近紅外波段光譜反射率比值可預測棉花葉片氮濃度。Pacheco-Labrador等[7]采用歸一化指數NDIS和三波段指數TBIs有效估算Holm oak葉片氮含量；Abdel-Rahman等[8]采用一階導數的SR指數(波段743，1 316 nm)(波段743，1 317 nm)(波段741，1 323 nm)估算的葉片氮含量的決定系數在0.75左右；Yao等[9]采用歸一化指數NDSI和比值指數RSI一階導數對葉片氮的累積量進行了估算，其決定系數高達0.81；Ullah 等[10]利用MERIS數據采用植被指數(NDVI，NBDI，SAVI，TSAVI，REIP，MTCI)以及波段深度參數對草地生物量及氮含量進行了估算和對比。Sanches等[11]采用反射率，吸光度及其衍生形式和偏最小二乘方法對植物氮磷鉀進行估算，結果表明一階導數形式的估算精度最高。另有學者著重從算法角度進行了植物氮含量高光譜反演技術研究，Zhang等[12]利用單變量線性回歸、逐步多元線性回歸、偏最小二乘對濕地植物蘆葦氮含量進行估算，并采用留一交叉驗證比較3種模型的精度，認為氮的敏感波段集中在紅波段和綠波段，單變量線性回歸對665 nm和680 nm處的歸一化植被指數效果最好，3種方法中逐步線性回歸的精度最高。Yao等[13]采用原始反射率、一階導數、植被指數等變量利用SMLR、PLSR、ANNS、SVMs等反演了冠層葉片氮含量，結果表明一階導數的支持向量機方法精度更高。Wang 等[3]驗證了多核支持向量機在估算小麥葉片氮含量的有效性，并與多元線性回歸，偏最小二乘，人工神經網絡，單核支持向量機進行了對比。總體來說，現有的關于高光譜植物氮含量監測的大部分研究集中在植被指數反演氮含量以及敏感波段篩選方面，關于建模算法的研究較少[3]。且已有的研究或集中于建模之前敏感特征的篩選或集中于建模算法的比較，本文嘗試將不同特征提取方法與建模方法相結合，通過不同方法間的組合，篩選出一套完整的能有效估算葉片氮含量的流程。

1 材料和方法

1.1 光譜數據與氮含量數據獲取

數據集采用Accelerated canopy chemistry program(ACCP)(https://doi.org/10.3334/ORNLDAAC)[14]，該數據集旨在研究不同生態系統中植被冠層氮含量及木質素含量遙感反演的理論基礎，包含野外樣品的實驗室化學分析數據、光譜數據、小型冠層試驗的化學分析數據及冠層建模數據。本研究選擇ACCP中91組道格拉斯冷杉幼苗期新鮮葉片光譜數據及相應的實驗室化學分析數據，采集于1992年11月。光譜數據為400～2 498 nm處的吸光度[log(1/反射率)]，光譜間隔為2 nm，分辨率為10 nm。在采集光譜數據后的6 h，采集葉片測量葉片氮及葉綠素含量。采用常規的實驗室方法，對每棵樹幼苗期的葉片樣本進行干燥、研磨和全氮、全葉綠素(a+b)分析。樣品經硫酸-氧化汞催化劑(Perstorp分析)在塊狀消化器中消化后，用Alpkem連續流自動分析儀測定總氮,氮含量的單位為單位干重百分比。

1.2 數據分析

利用Matlab對試驗數據進行了異常數據剔除、計算了原始反射率、一階導數R′、反射率倒數的對數log(1/R)及其一階導數log(1/R)′，以及這些參數與氮含量的相關系數，采用連續投影算法、LASSO、相關系數篩選以上4種參數形式的特征波段，并采用隨機森林、偏最小二乘以及支持向量機建立氮含量反演模型。其中，相關系數、連續投影算法篩選波段組合，偏最小二乘建模在Matlab中完成，LASSO算法篩選變量組合、隨機森林和支持向量機建模在R中完成。

1.2.1 相關系數法 相關系數法是最簡單直接的能夠幫助理解特征與相應變量之間關系的方法，衡量的是變量之間的線性相關性，其值絕對值越大表示相關性越強。很多研究通過選擇最大相關性的變量作為特征變量進行建模研究。λ1-λ2，λ1/λ2表示兩變量間的差異性，其與因變量間的相關系數表示兩變量組合對因變量的解釋量，選取相關系數最大的兩個波長即為特征波長。

1.2.2 連續投影算法 連續投影算法(SPA)是一種矢量空間共線性最小化的前向變量選擇算法，能夠有效從全波段中提取特征波段，消除原始光譜矩陣中的冗余信息，降低模型的復雜度，在波長選取中取得了較好的效果[15]。連續投影算法是一種前向循環選擇算法，設樣本集M和波段數K組成一個光譜矩陣XMXK，分別記Xk(0)和N為初始的迭代向量和需要提取的波段個數。從一個波長開始，每次循環計算它在未選入的波長上的投影，將投影向量最大的波長引入波長組合，直到循環N次，每一個新入選的波段，都與前一個線性關系最小。通過循環會得到N×K對波段組合，對每一對XK(0)和N所決定的組合分別建立多元回歸模型，并用預測均方根誤差RMSE來決定所建模型的優劣，最小的RMSE對應的Xk(0)為最佳的波段組合。

1.2.4 隨機森林 隨機森林是一種并行式集成學習法，以決策樹為基礎，并在決策樹訓練過程中引入隨機屬性選擇。它通過自助法重采樣技術，從原始訓練樣本集N中有放回地重復隨機抽取K個樣本生成新的訓練樣本集合，然后根據自助樣本集生成K個分類樹組成隨機森林[17-18]。

1.2.5 支持向量機 支持向量機是一種基于統計學習理論的機器學習算法，它建立在VC(Vapnik-Chervonenkis Dimension)維理論和結構風險最小原理基礎上，能較好地解決小樣本、非線性、高維數和局部極小點等實際問題，同時能獲得較好的泛化能力[19]。

1.2.6 偏最小二乘 偏最小二乘回歸是一種新型的多變量回歸分析方法，可以實現回歸建模、簡化數據結果和分析兩組變量見得相關性，給多元數據統計分析帶來極大便利，它能在回歸建模過程中采用數據將為、信息綜合和篩選技術，提取對系統最佳解釋能力的新綜合成分[20]。

1.3 建模方法與精度評價

分別采用偏最小二乘，隨機森林與支持向量機建立葉片氮含量反演模型，并采用決定系數R2、均方根誤差來評價各個模型的精度。R2表示相關密切程度，RMSE用來衡量估測值與真實值之間的偏差程度。R2越大，RMSE的絕對值越小，表示模擬結果精度更高。

2 結果與分析

2.1 不同光譜形式與氮含量的相關性分析

原始反射率R與氮含量呈負相關關系，在500～730 nm波段的相關系數超過-0.6，在550、702 nm左右相關性達到最大，且差異顯著(P<0.05)，在550 nm左右相關性最大是由于葉綠素對綠光的反射作用，702 nm位于紅邊區域，有相關研究發現紅邊波段與植物葉綠素含量、生物量等參數存在顯著相關性；氮含量與原始光譜的相關性在近紅波段較小(表1)。

一階導數R′能夠有效消除葉面積變化的影響，反射率的一階微分與氮含量的相關性變化較為復雜，在490～544、556～710和736～800 nm超過0.6，其中R′在530、574、626，694 nm達到最大，近紅波段的相關性較小，但在1 200、1 672、2 164、2 310 nm相關性分別達到峰值0.6以上(圖1)。

根據工程的復雜程度，對設計資質提出具體要求，如對提灌站、壓力輸水管道、100 m以上的深井、規模較大的建筑物等技術含量較高的工程必須聘請有相應資質的單位進行設計，從源頭抓好工程質量。

圖1 R，R′，log(1/R)，log(1/R)′與氮含量的相關性

Log(1/R)相關性曲線與原始反射率的相關性曲線線型一致，但相關性相反。log(1/R)′與氮含量的相關性變化也較為復雜，在490～546，620～670，672～694，702～804 nm相關性在0.6以上，在紅外波段1 204，1 672，2 158 nm的相關性形成峰值，達到0.6以上。

應用λ1-λ2，λ1/λ2與氮含量的相關系數篩選出的波段組合如表1所示，R與log(1/R)篩選出的波段集中在可見光波段，兩種導數形式篩選出的波段集中在近紅波段；原始波段反射率及其變換形式篩選出的波段組合都包含522，526 nm；導數形式篩選出的波段組合均包含2 158，2 044和2 070 nm 3個波段，雖然這幾個單獨的波段與氮含量的相關性并不強，但其波段組合與氮含量的相關性達到0.86以上，說明該波段組合包含了氮含量的大部分信息。

表1 相關系數法篩選的波段組合

2.2 連續投影算法篩選敏感波段特征篩選

ACCP葉片數據集共91組鮮葉片數據，隨機選取其中65組數據進行連續投影算法。連續投影方法將65組數據分為訓練集和校正集，以校正集的預測均方根誤差來確定最佳的光譜變量總數，均方根誤差和決定系數來評價最佳的波段變量組合。對以上4種形式的光譜參數分別采用連續投影算法，確定氮含量反演的最佳光譜變量(表2)。

表2 連續投影法篩選的波段組合

原始反射率所選的波段在波段范圍內分布較為均勻(圖2)，位于可見光的波段有6個，2個波段位于紅邊區域，位于波峰或波谷處的波段有6個。一階導數只選擇了3個變量，其中一個位于可見光波段。Log(1/R)投影算法選擇了兩個紅光波段，其余均為近紅波段。log(1/R)′僅選擇了3個紅外波段。校正數據集RMSE最低的是Log(1/R)形式，RMSE最高的是一階導數形式。

圖2 連續投影算法篩選出的波段位置

2.3 LASSO算法

采用Lasso算法對4種光譜形式進行變量篩選，各光譜形式篩選出的變量數均比連續投影篩選出的變量數多。原始光譜反射率篩選出22個波段，其中有12個波段位于可見光區。Log(1/R)篩選出的變量數為16個，其中9個可見光波段。兩種導數形式均篩選出23個變量數，一階導數篩選出的可見光波段8個，Log(1/R)′篩選出的可見光波段6個(表3)。

表3 LASSO篩選的波段組合

2.4 模型預測結果與對比

分別將相關系數篩選法、連續投影算法、LASSO算法篩選出的波段作為偏最小二乘、隨機森林、支持向量機的輸入來估算葉片氮含量，以65組數據為建模數據，26組數據為驗證數據。不同方法組合的反演結果如表4及圖3所示。

表4 不同方法組合的反演結果

3種變量篩選方法和3種建模方法的組合中，RMSE在0.19～0.38，R2在0.6～0.89，其中，隨機森林算法對相關系數法和LASS0算法篩選出的Log(1./R)′形式變量的反演誤差最小，RMSE在0.19～0.20，能夠解釋氮含量的變化，但該算法對篩選出的其他形式的變量反演誤差較大，表現并不穩定;不論采用何種建模方法，對于相關系數法，其log(1/R)′形式的變量反演效果較好；對于連續投影法，其篩選出的R′變量，反演誤差較小；對于LASSO算法，其篩選出的R'和log(1/R)′形式的變量反演效果要明顯優于其余兩種變量形式。同時，分別采用偏最小二乘，隨機森林、支持向量機進行全波段的建模，反演結果分別為R2=0.605、RMSE=0.4542，R2=0.659、RMSE=0.3352，R2=0.730 4、RMSE=0.348 4。表明，相關系數&Log(1./R)′，LASSO&Log(1./R)′，SPA&R′能夠有效篩選出冷杉葉片氮含量的敏感波段組合，且前兩種方法組合效果更優。

3 討論

3.1 不同變量篩選方法和建模方法的組合效果比較

從圖3可以看出兩種導數形式估算出的氮含量與實測氮含量1:1的對應關系更好，散點基本上沿著y=x線分布，具有較好的線性關系。對于R及log(1/R)形式的反演，三種建模方法的結果與實測值相比普遍偏高，大部分點分布在y=x線以上。兩種導數形式的反演結果較為均勻緊致地分布在y=x線兩側，反演結果與實測數據更接近。從實測數據與反演結果的擬合圖可以看出，相關系數法&log(1/R)′以及LASSO&log(1/R)′中實測數據與反演結果的擬合函數具有較大的斜率，較小的截距，且R2較大，RMSE較小，具有較高的反演精度(圖3)。

圖3 實測氮含量與不同方法反演氮含量散點圖

3.2 3種篩選方法的敏感波段分析

4 結論

采用相關系數法、連續投影法、LASSO算法對4種形式的光譜變量R，R′，log(1/R)，log(1/R)′進行敏感波段篩選，對篩選后的敏感波段分別采用偏最小二乘、隨機森林、支持向量機對氮含量進行建模反演，得出以下結論：

1)兩種導數形式變量的反演誤差最小。

2)相關系數&Log(1/R)′，LASSO&Log(1/R)′能夠有效篩選出冷杉葉片氮含量的敏感波段組合，無論采取何種建模方法，其估算效果在幾種篩選方法組合中最好，R2>0.84，RMSE在0.19～0.24，估算效果明顯優于采用3種建模方法進行全波段建模。

3)隨機森林算法對相關系數法和LASS0算法篩選出的Log(1/R)′形式變量的反演誤差最小，但對其他形式篩選出的變量反演結果的誤差變化范圍較大，表現并不穩定。

4)3種變量篩選方法篩選出的R形式的變量與前人研究相符，連續投影算法能篩選出更多與葉片其他化學含量相關的波段，證明了3種變量篩選方法的有效性。