基于葉片光譜反射率的毛竹氮元素含量估測(cè)研究

林靈辰，余坤勇，曾 琪，姚 雄，鄧洋波，范華棟，劉 健

（1.福建農(nóng)林大學(xué)林學(xué)院，福建 福州 350002； 2.3S 技術(shù)與資源優(yōu)化利用福建省高校重點(diǎn)實(shí)驗(yàn)室，福建 福州 350002）

毛竹是我國(guó)重要竹類樹(shù)種，同時(shí)也是經(jīng)濟(jì)價(jià)值很高的筍、材兩用竹種[1]，其廣泛地分布于我國(guó)南方低山丘陵區(qū)[2]。根據(jù)我國(guó)第八次森林資源清查統(tǒng)計(jì)，中國(guó)毛竹面積已達(dá)430 萬(wàn)hm2，占所有竹林總面積的74%，在我國(guó)所有竹類資源中占有絕對(duì)優(yōu)勢(shì)[3]。氮元素是毛竹生長(zhǎng)發(fā)育過(guò)程中的必要營(yíng)養(yǎng)元素之一[4]，毛竹葉片中的氮元素含量在一定程度上能夠反映毛竹生長(zhǎng)狀況以及毛竹林地土壤肥力，所以估測(cè)毛竹林氮元素含量對(duì)毛竹生長(zhǎng)脅迫、氮元素缺乏情況、土壤肥力的間接估測(cè)等具有重要作用。

通過(guò)建立植被的生化參數(shù)與植被氮元素含量的模型進(jìn)行反演是當(dāng)前氮元素含量估測(cè)的主要方式，而光譜數(shù)據(jù)能夠獲取十分精細(xì)的植被信息，因此被廣泛應(yīng)用于植被氮元素含量的反演的研究中[5]，利用光譜數(shù)據(jù)對(duì)植物葉片生化組分進(jìn)行快速、無(wú)損的監(jiān)測(cè)已成為植物生長(zhǎng)狀況評(píng)價(jià)的重要內(nèi)容[6]。1997年宮鵬和浦瑞良率先使用光譜數(shù)據(jù)完成了森林氮元素含量的反演[7]，隨后薛紅利[8]、Clark[9]、Mutang[10]等采用簡(jiǎn)單的逐步回歸方法估測(cè)植被氮元素含量，其估測(cè)值與實(shí)測(cè)值具有相當(dāng)高的相關(guān)性。邢東興等[11]利用葉片原始反射率及其一階導(dǎo)數(shù)光譜、一階導(dǎo)數(shù)的對(duì)數(shù)光譜、二階微分光譜等基于原始光譜的變換形式與果樹(shù)葉片的全氮元素、全磷等含量的相關(guān)性進(jìn)行分析從中選擇最優(yōu)的光譜數(shù)據(jù)形式構(gòu)建逐步多元線性回歸方程。隨著研究的深入，有學(xué)者發(fā)現(xiàn)逐步多元線性回歸方法存在過(guò)擬合[12]、參與波段的多重共線性且與生化組分無(wú)關(guān)[13]等問(wèn)題。隨機(jī)森林算法（RF）具有很強(qiáng)的抗過(guò)擬合能力，因而被廣泛地運(yùn)用于分類與回歸分析中，李旭青 等[14]基于隨機(jī)森林算法估測(cè)的水稻冠層氮元素含量研究結(jié)果表明隨機(jī)森林算法具有普適性、不會(huì)過(guò)擬合等優(yōu)點(diǎn)。支持向量機(jī)（SVM）方法是由Vapnik 基于內(nèi)核統(tǒng)計(jì)理論于1995年提出，大量研究表明它可以很好地運(yùn)用于非線性回歸問(wèn)題，梁亮等[15]基于一階微分敏感波長(zhǎng)采用差值、比值以及歸一化的方法構(gòu)建了12 種光譜指數(shù)以及常用的22種植被指數(shù)，以此為參數(shù)采用最小二乘支持向量機(jī)回歸算法對(duì)小麥冠層氮元素含量進(jìn)行估測(cè)，同時(shí)優(yōu)選出最佳的植被指數(shù)。

毛竹氮元素含量對(duì)研究毛竹生長(zhǎng)狀況以及估測(cè)毛竹林地土壤肥力具有重要意義，因此探索如何快速并精準(zhǔn)地獲取毛竹氮元素含量具有重要意義。為此，基于毛竹葉片光譜反射率構(gòu)建了3 種毛竹葉片氮元素含量估測(cè)模型，并從中篩選出毛竹葉片氮元素含量的最優(yōu)估測(cè)模型，為土壤肥力的間接估測(cè)以及進(jìn)一步深入研究基于冠層尺度光譜估測(cè)毛竹氮元素含量奠定基礎(chǔ)。

1 研究區(qū)概況

試驗(yàn)點(diǎn)分布于福建省順昌縣大干鎮(zhèn)的武坊村、土壟村、干山村、良坊村以及1 個(gè)鄉(xiāng)林場(chǎng)。其 經(jīng) 緯 度 位 置 為26°50′30″ ～26°58′0″N，117°33′30″～117°47′30″E（圖1），屬于中亞熱帶季風(fēng)氣候，氣候溫和、雨量充沛，年平均降水量約為1 752 mm，年平均日照數(shù)約為1 700 h，年平均氣溫為16.3℃，海拔范圍在118 ～1 282 m 之間。研究區(qū)地處山區(qū)，森林資源、水利資源十分豐富。大干鎮(zhèn)總面積20 942 hm2，森林覆蓋率高達(dá)78.8%，其中毛竹面積3 533.33 hm2，占林業(yè)用地面積的21.41%。

圖1 研究區(qū)采樣點(diǎn)Fig.1 Sampling points in the study area

2 數(shù)據(jù)采集及處理

2.1 葉片光譜數(shù)據(jù)的采集與預(yù)處理

2017年7月8—24日在福建省順昌縣大干鎮(zhèn)布設(shè)59 個(gè)半徑為3.26 m 的樣圓串，并在3 個(gè)樣圓內(nèi)各選取1 株生長(zhǎng)良好的毛竹作為樣本樹(shù)，將樣本樹(shù)伐倒后摘取冠層的上中下3 層葉子各2 片，每片葉子在黑色背景下測(cè)量5 次其正面光譜，剔除異常值后將3 個(gè)樣本樹(shù)的光譜數(shù)據(jù)取均值作為樣地實(shí)測(cè)光譜數(shù)據(jù)，采用ISI921VF-256 野外地物光譜輻射計(jì)測(cè)量進(jìn)行光譜的測(cè)量，光譜范圍為可見(jiàn)-近紅外（380 ～1 080 nm），光譜分辨率約為4 nm。采用The Unscrambler X10.4 軟件中的Savizky-Golay 卷積平滑算法對(duì)原始光譜進(jìn)行預(yù)處理，并用Origin 9.0 軟件對(duì)其進(jìn)行線性插值，光譜間隔為1 nm。

2.2 葉片氮元素含量的測(cè)定

每個(gè)樣圓串的樣本樹(shù)冠層的上中下三層共摘取6 袋葉子，裝入信封中后放入烘箱中進(jìn)行105 ℃ 殺青20 min，之后75 ℃恒溫進(jìn)行烘干，直到葉片樣品恒質(zhì)量為止。用粉碎機(jī)粉碎后過(guò)篩形成待測(cè)樣品，取0.1 g 待測(cè)樣，使用德國(guó)VARIO MAX 碳氮元素分析儀進(jìn)行測(cè)定，其工作原理是將待測(cè)樣品完全燃燒后測(cè)定樣品的全氮、全碳含量，實(shí)驗(yàn)對(duì)同一樣品進(jìn)行3 次重復(fù)平行試驗(yàn)，取3 個(gè)數(shù)據(jù)中的2 個(gè)誤差較小的值取均值作為葉片全氮元素的實(shí)測(cè)數(shù)據(jù)。共59 個(gè)樣地的實(shí)測(cè)數(shù)據(jù)，剔除6 個(gè)異常樣本后隨機(jī)選取43 個(gè)樣本作為建模集，10 個(gè)樣本作為驗(yàn)證集。

2.3 模型變量的選取

大量研究表明一階微分光譜能夠降低環(huán)境因素對(duì)目標(biāo)光譜反射率的影響[16]，在構(gòu)建氮元素估測(cè)模型時(shí)可以作為一個(gè)重要的變量類型。光譜植被指數(shù)是依據(jù)目標(biāo)地物的光譜特性將某幾個(gè)波段利用數(shù)學(xué)方式組合成的光譜指數(shù)，其具有指示植物體內(nèi)某種生化組分的情況、減少背景干擾和突出目標(biāo)信息等作用[17]，在氮元素含量估算研究中，常見(jiàn)的光譜植被指數(shù)有光化學(xué)植被指數(shù)（Photochemical reflectance index，PRI）、比值植被指數(shù)（Ratio vegetation index，RVI）、歸一化植被指數(shù)（Normalized vegetation index，NDVI）等，參考國(guó)內(nèi)外學(xué)者的研究，共選取了25 種與氮元素相關(guān)的光譜植被指數(shù)（表1）與毛竹葉片氮元素含量進(jìn)行相關(guān)性分析并篩選模型變量。

2.4 模型構(gòu)建原理

2.4.1 多元線性模型

多元線性回歸是使用多個(gè)自變量對(duì)因變量進(jìn)行解釋的回歸方法[18]，利用多元逐步回歸方法能夠自動(dòng)篩選與目標(biāo)變量相關(guān)關(guān)系最大的自變量進(jìn)入模型從而去除無(wú)關(guān)變量的優(yōu)點(diǎn)[19]，逐步篩選與毛竹葉片氮元素含量最密切的參數(shù)構(gòu)建模型。

2.4.2 隨機(jī)森林

隨機(jī)森林算法（RF）是Breiman 等于2001年提出一種機(jī)器學(xué)習(xí)算法[20],由多個(gè)單獨(dú)決策回歸樹(shù)組成，其工作原理為：用boostrasp 算法有放回地從N 個(gè)原始數(shù)據(jù)袋中隨機(jī)抽取n 個(gè)樣本形成新的建模集，未抽到的樣本作為預(yù)測(cè)集樣本估計(jì)模型的性能。若原始數(shù)據(jù)集有P 個(gè)變量，則每個(gè)決策樹(shù)的每個(gè)節(jié)點(diǎn)抽取mtry 個(gè)變量，一般mtry=P/3，并根據(jù)這些特征變量計(jì)算最佳的分裂方式并對(duì)節(jié)點(diǎn)進(jìn)行分裂，在運(yùn)行不斷生成大量決策樹(shù)直到達(dá)到預(yù)先設(shè)定的決定樹(shù)數(shù)目從而達(dá)到回歸分析的目的。

表1 選取的植被指數(shù)Table 1 Selected vegetation index

2.4.3 支持向量機(jī)

支持向量機(jī)（SVM）方法是由Vapnik 基于內(nèi)核統(tǒng)計(jì)理論于1995年提出，其原理是通過(guò)一個(gè)非線性映射P 將樣本空間映射到一個(gè)高維空間，并在這個(gè)空間將非線性可分問(wèn)題轉(zhuǎn)化為特征空間中的線性可分問(wèn)題[21-23]，其實(shí)現(xiàn)關(guān)鍵在于核函數(shù)k(xi,xj)=φ( xi) ·φ(xj)，常見(jiàn)的核函數(shù)有4種：線性核函數(shù)、多項(xiàng)式核函數(shù)、徑向基核函數(shù)、Sigmoid核函數(shù)。

3 結(jié)果與分析

3.1 不同氮元素含量的毛竹葉片光譜特征

氮元素對(duì)可見(jiàn)光和近紅外波段不存在明顯吸收、反射與透射特征，其敏感的吸收波段主要為短波紅外，但是吸收特征相對(duì)較弱。氮元素的吸收特征主要是由于氮與葉綠素存在較強(qiáng)的相關(guān)關(guān)系，進(jìn)而通過(guò)葉綠素的吸收特征進(jìn)行表達(dá)[24]。圖2 為毛竹不同氮含量的葉片光譜反射率。不同氮含量的毛竹葉片光譜反射率的規(guī)律具有很好的一致性，在可見(jiàn)光波段范圍（380 ～760 nm）內(nèi)，由于葉綠素的強(qiáng)反射有一綠色反射峰，在600 ～760 nm 的波段范圍內(nèi)，由于葉綠素強(qiáng)烈地吸收輻射能而形成吸收谷，其光譜反射率小于10%，約在680 ～760 nm 的波段范圍，反射率急劇上升，在760 ～1 080 nm，形成一個(gè)高反射平臺(tái)，該平臺(tái)光譜反射率范圍在40%～60%之間，隨著氮含量的增加，此波段范圍的反射率先降低后升高[25]。

圖2 不同氮元素含量的毛竹葉片反射率光譜Fig.2 Reflectance spectra of leaves of moso bamboo with different nitrogen content

3.2 氮元素估測(cè)模型變量的選取

將毛竹葉片氮元素含量與原始光譜、一階微分光譜進(jìn)行相關(guān)性分析，結(jié)果如圖3 所示，在原始光譜386 ～389 nm、394 ～396 nm 和411 ～415 nm的波段范圍與氮元素含量處于極顯著相關(guān)，其中在波長(zhǎng)R387處相關(guān)性達(dá)到最大，相關(guān)系數(shù)為0.445 4，與一階微分光譜的相關(guān)系數(shù)在380 ～1 080 nm 范圍內(nèi)的變化幅度比與原始光譜相關(guān)系數(shù)大，在波長(zhǎng)385 nm、523 ～524 nm、600 ～601 nm、663 ～ 664 nm、666 nm、691 ～692 nm、773 nm、761 nm、 794 ～795 nm 處與葉片氮元素含量呈極顯著相關(guān)，最大相關(guān)系數(shù)為0.490 0，對(duì)應(yīng)波長(zhǎng)為DR663。由于毛竹葉片氮元素含量與原始光譜、一階微分光譜的相關(guān)性浮動(dòng)較大，相關(guān)波長(zhǎng)比較接近，因而選取R387、DR663作為構(gòu)建氮元素估測(cè)模型的變量。

將毛竹葉片氮元素含量與25 種植被指數(shù)進(jìn)行相關(guān)性分析，根據(jù)兩者之間的相關(guān)性以及顯著性進(jìn)行敏感參數(shù)的選取，結(jié)果如表2，達(dá)到極顯著相關(guān)的植被指數(shù)有歸一化植被指數(shù)NDVIg-b(R575、R440)、結(jié)構(gòu)不敏感色素SIPI、光化學(xué)反射指數(shù)PRI 和色素比值植被指數(shù)PPR，相關(guān)系數(shù)分別為 0.472 1、0.382 4、0.401 5 和0.380 3，故選取這4種植被指數(shù)作為變量進(jìn)行模型的構(gòu)建。

3.3 氮元素含量估測(cè)模型的構(gòu)建

3.3.1 多元線性模型

將與毛竹葉片氮元素含量相關(guān)性較高的6 個(gè)敏感參數(shù)作為輸入自變量，利用建模集中的43 個(gè)樣本進(jìn)行逐步回歸擬合，結(jié)果如表2所示，模型共進(jìn)入了DR663、PRI 2 個(gè)變量，模型方程為Y=0.333- 25.631*DR663-0.939*PRI，擬合系數(shù)R2為0.363。

圖3 毛竹葉片氮元素含量與葉片原始光譜、一階微分光譜的相關(guān)性Fig.3 The correlation between leaf nitrogen content and leaf original spectrum and first-order differential spectrum of moso bamboo

3.3.2 非線性模型

為了保證隨機(jī)森林模型以及支持向量機(jī)模型擬合得到的結(jié)果具有對(duì)比性，均采用篩選得到的6 個(gè)敏感參數(shù)組合作為隨機(jī)森林和支持向量機(jī)的輸入變量，同樣利用建模集中的43 個(gè)樣本進(jìn)行擬合。隨機(jī)森林的mtry 與決策樹(shù)數(shù)量設(shè)置分別為2和1 000，從多次試驗(yàn)的結(jié)果中取得的隨機(jī)森林最優(yōu)擬合結(jié)果的擬合系數(shù)R2為0.353 2。SVM 算法采用徑向基函數(shù)進(jìn)行葉片氮元素含量估測(cè)模型的構(gòu)建，經(jīng)多次反復(fù)試驗(yàn)發(fā)現(xiàn)當(dāng)懲罰因子C 和核參數(shù)Sigma 分別設(shè)置為3 和0.1 時(shí)的擬合效果最好，其估測(cè)結(jié)果的擬合系數(shù)R2達(dá)到了0.799 7。

表2 毛竹氮元素含量估測(cè)模型Table 2 Estimation model of nitrogen content in moso bamboo

3.4 模型精度檢驗(yàn)分析

用驗(yàn)證集中的10 個(gè)樣本對(duì)3 個(gè)模型進(jìn)行精度檢驗(yàn)，采用決定系數(shù)R2、均方根誤差RMSE、預(yù)測(cè)偏差值bias 和總體平均精度來(lái)評(píng)價(jià)模型的估測(cè)能力。結(jié)果見(jiàn)圖4，三個(gè)模型的估測(cè)值與實(shí)測(cè)值都在Y=X 附近波動(dòng)，支持向量機(jī)的估測(cè)值與實(shí)測(cè)值的擬合度較高，而多元線性模型和隨機(jī)森林模型的擬合點(diǎn)相對(duì)分散。從精度分析表3 可以看出，3 個(gè)估測(cè)模型的總體精度均達(dá)到了90%以上，其中支持向量機(jī)的R2是三個(gè)模型中最高的，達(dá)到了0.803 1，其預(yù)測(cè)偏差值與均方根誤差相當(dāng)小，綜合均方根誤差值、R2以及總體精度來(lái)看隨機(jī)森林模型估測(cè)效果要略優(yōu)于多元線性模型，二者與支持向量機(jī)模型估測(cè)結(jié)果相比效果較差。綜合三個(gè)模型的精度檢驗(yàn)結(jié)果來(lái)看，懲罰因子C 和核參數(shù)Sigma 分別設(shè)置為3 和0.1 的支持向量機(jī)模型的估測(cè)效果最佳。

4 結(jié)論與討論

氮元素含量是影響毛竹生產(chǎn)的主要因素之一，同時(shí)也是反映毛竹生長(zhǎng)狀況的重要指標(biāo)，因此研究如何快速獲取毛竹氮元素含量具有重要意義，本研究利用地面實(shí)測(cè)光譜數(shù)據(jù)，從葉片角度探索毛竹林氮元素含量估測(cè)模型并取得了一定成果。

圖4 模型精度比較Fig.4 Comparison of model accuracy

表3 精度分析Table 3 Precision analysis

1）結(jié)合毛竹葉片氮元素含量與葉片原始光譜、一階微分光譜和植被指數(shù)三者之間的相關(guān)性分析結(jié)果發(fā)現(xiàn)：葉片原始光譜在波長(zhǎng)387 nm 處相關(guān)性達(dá)到最大，一階微分光譜在波長(zhǎng)663 nm 處相關(guān)性達(dá)到最大，NDVIg-b、SIPI、PRI 和PPR4種植被指數(shù)與葉片氮元素含量達(dá)到極顯著相關(guān)?？梢钥闯雠c毛竹葉片氮元素含量達(dá)到極顯著相關(guān)的原始光譜波段相對(duì)較少，這可能是因?yàn)樵诳梢?jiàn)-近紅外（380 ～1 080 nm）范圍內(nèi)其他生化組分的吸收波段與氮元素相近而產(chǎn)生干擾[26]，并且在實(shí)際的光譜操作過(guò)程中，不同的光譜分辨率、不同的光譜處理方法等也會(huì)對(duì)植被生化參數(shù)估測(cè)造成一定的影響。一階微分光譜與氮元素含量達(dá)到極顯著相關(guān)的波段遠(yuǎn)遠(yuǎn)多于原始光譜，這一結(jié)果與李哲等[27]、 丁雅等[28]的研究結(jié)果相一致，說(shuō)明原始光譜的一階微分變化形式在研究光譜與植被氮元素含量之間的關(guān)系時(shí)具有一定的優(yōu)勢(shì)。

2）利用6 個(gè)氮元素敏感特征參數(shù)作為輸入變量構(gòu)建毛竹氮元素含量多元線性估測(cè)模型、隨機(jī)森林模型以及支持向量機(jī)模型，綜合模型精度分析結(jié)果來(lái)看，支持向量機(jī)模型的R2、均方根誤差、總體精度是3 個(gè)模型中最好的，這一結(jié)果與楊曦光[29]、李百超[30]的研究結(jié)果相一致。多元線性模型的決定系數(shù)以及整體精度略低于隨機(jī)森林模型，但其預(yù)測(cè)偏差值過(guò)高，因此多元線性模型在用于估測(cè)毛竹氮元素含量時(shí)不具有優(yōu)勢(shì)，支持向量機(jī)模型在3 種模型中的估測(cè)結(jié)果效果是最好的，這與模型本身的特性有關(guān)，支持向量機(jī)在用于擬合樣本較少的數(shù)據(jù)時(shí)具有更大優(yōu)勢(shì)[30]。

基于葉片光譜數(shù)據(jù)構(gòu)建毛竹氮元素含量估測(cè)模型，對(duì)比精度檢驗(yàn)結(jié)果發(fā)現(xiàn)支持向量機(jī)在估測(cè)葉片氮元素含量方面具有明顯優(yōu)勢(shì)，其實(shí)測(cè)值與預(yù)測(cè)值的擬合系數(shù)與總體精度分別為0.803 1、94.02%，因此能夠作為毛竹葉片氮元素含量的較佳估測(cè)方法。實(shí)驗(yàn)過(guò)程中存在幾點(diǎn)不足，在之后的研究中將加以深化：可以增加多種原始光譜變化形式以及不同輸入變量組合對(duì)毛竹氮元素含量的影響，同時(shí)增加樣本數(shù)量來(lái)避免由于在實(shí)際光譜測(cè)量過(guò)程以及氮元素含量的測(cè)定過(guò)程中產(chǎn)生的實(shí)驗(yàn)誤差而導(dǎo)致模型擬合精度下降的情況。