依據(jù)機器學(xué)習(xí)算法的杉木干形模擬1)

梁瑞婷 周來 謝運鴻 丁志丹 孫玉軍

(森林資源和環(huán)境管理國家林業(yè)和草原局重點開放性實驗室(北京林業(yè)大學(xué))，北京,100083)

林木干形是描述樹干不同高度處的截面直徑隨樹干縱向位置的變化，是研究林木最重要的因子之一，它既與樹干的出材量有關(guān)，又影響木材的質(zhì)量[1]。通過削度模型可以精確估算單株林木的樹干材積，對重建樹干輪廓、估算林分蓄積和森林碳儲量、評價森林經(jīng)營水平和森林規(guī)劃設(shè)計等方面具有重要意義[2]。確定干形，解析木截段測量是最精確的方法，但需要耗費大量人力、物力，而且破壞性大，因此，一般采用削度方程。削度方程即是利用樹木胸徑、樹高、截面高等易測量的林木特征因子作為自變量預(yù)測截面直徑的回歸方程[3]。盡管過去已有多種統(tǒng)計學(xué)方法建立林木的削度方程，如：線性與非線性回歸[4]、分位數(shù)回歸[5-6]和混合效應(yīng)模型[7]等，但這些方法在應(yīng)用時對數(shù)據(jù)分布有很強的要求，需要滿足一定的統(tǒng)計學(xué)假設(shè)前提，如數(shù)據(jù)的獨立、正態(tài)、等方差等[8]。而林木生長是由多個因素交互作用、共同影響[9]，其數(shù)據(jù)通常難以滿足這些假設(shè)條件，需要探索新的建模方法作為補充。

機器學(xué)習(xí)是新近興起的一門涉及多領(lǐng)域的交叉學(xué)科，作為一種數(shù)據(jù)驅(qū)動的算法，機器學(xué)習(xí)能夠自動從數(shù)據(jù)中分析規(guī)律，并應(yīng)用于新的樣本，非常適合處理復(fù)雜的非線性問題[10-11]。與傳統(tǒng)回歸相比，機器學(xué)習(xí)對數(shù)據(jù)分布沒有要求，能高效處理變量之間非正態(tài)、非線性和共線性關(guān)系[12]。而且機器學(xué)習(xí)可以通過一定措施改進模型，進一步提高模型精度，如增加訓(xùn)練集數(shù)據(jù)、調(diào)整學(xué)習(xí)參數(shù)或改變模型結(jié)構(gòu)。近年來，機器學(xué)習(xí)逐漸開始應(yīng)用于林學(xué)領(lǐng)域的研究，如：Martins Silva J P et al.[13]利用機器學(xué)習(xí)和混合模型估算了巴西大草原上單株林木的材積和生物量，Kilham P et al.[14]利用隨機森林估測了德國西南部多種森林類型的林分蓄積，Ozcelik R et al.[15]人利用人工神經(jīng)網(wǎng)絡(luò)和非線性回歸建立了紅松樹高曲線模型[16]。這些研究都取得了良好的結(jié)果，表明機器學(xué)習(xí)在森林生長收獲和森林資源監(jiān)測方面具有很大的應(yīng)用潛力，也為研究林木干形提供了一種新思路。國外有學(xué)者展開相關(guān)研究，如Nunes et al.[16]采用隨機森林和人工神經(jīng)網(wǎng)絡(luò)對巴西復(fù)雜雨林內(nèi)多種林木的干形進行估測，發(fā)現(xiàn)人工神經(jīng)網(wǎng)絡(luò)的精度高于傳統(tǒng)削度方程。Sanquetta et al.[17]采用K-近鄰法和人工神經(jīng)網(wǎng)絡(luò)對巴西地區(qū)柳杉的干形進行模擬，發(fā)現(xiàn)最近鄰法能夠提高干形估測的精度，是一種簡單明了的估測方法。Schikowski et al.[18]應(yīng)用3種機器學(xué)習(xí)算法對黑松的樹干削度和材積建立模型，表明機器學(xué)習(xí)一種估測森林參數(shù)可行的方法。迄今為止，國內(nèi)應(yīng)用機器學(xué)習(xí)對林木干形模擬的研究較少。

本文以福建將樂國有林場的杉木(Cunninghamialanceolata)為研究對象，利用最近鄰算法、隨機森林和人工神經(jīng)網(wǎng)絡(luò)3種機器學(xué)習(xí)模擬其干形，估測樹干不同高度處對應(yīng)的直徑，并與傳統(tǒng)削度模型進行比較。通過4個模型評價指標(biāo)，結(jié)合殘差圖和相對偏差圖進行排序，以得到高效、低偏的干形估計方法，為杉木干形的精確估測提供參考。

1 研究區(qū)概況

本研究的解析木數(shù)據(jù)來源于福建省三明市將樂縣(117°5′～117°40′E，26°26′～27°4′N)的將樂國有林場，該地區(qū)的年均氣溫20°左右，年均降水量約2 700 mm，氣候溫暖濕潤，雨熱同期，以中、低山為主，海拔多在400～800 m范圍內(nèi)，土壤肥沃深厚。主要樹種有馬尾松(PinusmassonianaLamb.)、毛竹(Phyllostachysheterocycla(Carr.) Mitford cv.Pubescens)等。

2 材料與方法

選擇不同齡組且有充分代表性的杉木人工純林，并在林分內(nèi)設(shè)置標(biāo)準(zhǔn)地(標(biāo)準(zhǔn)地面積為20 m×20 m或20 m×30 m)，共29塊。根據(jù)每木檢尺的結(jié)果計算林分的平均胸徑和平均樹高，并選取生長健康的平均木，盡量使樣木均勻分布于各個徑階，共得46株杉木解析木。伐倒前測量樹高(H)和胸徑(D)。伐倒后按1 m區(qū)分段截取圓盤，不足1 m的區(qū)分段作為梢頭。同時，在樹干基部和胸徑1.3 m處也分別截取圓盤。然后測量每個圓盤東西、南北兩個方向的直徑，計算其平均值作為對應(yīng)高度處的直徑。共計獲得793組數(shù)據(jù)，解析木的數(shù)據(jù)統(tǒng)計見表1。

表1 杉木調(diào)查因子統(tǒng)計

2.1 應(yīng)用非線性回歸的杉木削度方程

本研究從前人對杉木干形的研究成果中，選取了精度最高的削度模型[3-4]，該模型由曾偉生提出，其具體表達見公式(1)，方程的參數(shù)通過非線性最小二乘法確定，利用R4.00的nls()函數(shù)擬合求得。

(1)

式中:a1、a2、a3和a4為模型參數(shù)，D、H、di、hi分別表示胸徑、樹高、樹干直徑和樹干距離地面的高度。

2.2 機器學(xué)習(xí)算法

本研究采用的3種機器學(xué)習(xí)算法都是利用R4.0.0軟件，為了避免單一測試集的偶然性和隨機性，采用十折交叉驗證，即把所有數(shù)據(jù)10等分，每次用9/10的數(shù)據(jù)作為訓(xùn)練樣本，剩余1/10作為驗證樣本。這樣重復(fù)10次，每次采用不同的訓(xùn)練樣本和檢驗樣本，最后以10次結(jié)果的平均值為評價標(biāo)準(zhǔn)。

randomForest包實現(xiàn)隨機森林算法，隨機森林算法的精度通常與“樹”的個數(shù)有關(guān)，盡管設(shè)置的樣本(樹)數(shù)量越多，結(jié)果越穩(wěn)定[19]，樣本(樹)的設(shè)置需要依據(jù)具體的研究數(shù)據(jù)而定，默認(rèn)的樣本量為500，一般而言，樣本(樹)在500以后誤差趨于穩(wěn)定[20-21]，但是考慮本研究的數(shù)據(jù)量較少，取100≤樣本(樹)≤500以50為間隔進行參數(shù)調(diào)優(yōu)，最終確定樣本(樹)為300，以往的采用隨機森林研究樹干削度也采用樣本(樹)=300，模型表現(xiàn)最好[17]。根據(jù)周志華建議[19]，枝(mtry)設(shè)置為輸入變量數(shù)目的1/3左右，參考前人研究[22]，枝(mtry)分別取1、2和3，結(jié)果表明不同的枝(mtry)，模型擬合精度沒有明顯差異。

網(wǎng)絡(luò)包實現(xiàn)人工神經(jīng)網(wǎng)絡(luò)，在學(xué)習(xí)過程中，通常采用人工神經(jīng)網(wǎng)絡(luò)建模時為了防止過擬合，一般都將權(quán)值衰減參數(shù)(decay參數(shù))設(shè)為0.001，最大迭代次數(shù)(maxit)設(shè)為1 000次，如果達到目標(biāo)，就會提前終止迭代[19，23],隱蔽單元個數(shù)則是(size參數(shù))設(shè)置為15[18]。訓(xùn)練過程如圖1所示，由輸入層，2個隱藏層和一個輸出層組成，上一層的輸出通過權(quán)重轉(zhuǎn)化后作為下一層的輸入，不斷迭代，通過梯度下降法尋找最小誤差平方和。輸入因子為樹高(H)、胸徑(D)和不同高度(hi)，輸出因子為沿樹干不同高度處對應(yīng)的直徑(di)。

圖1 ANN模型的結(jié)構(gòu)圖

2.3 模型評價指標(biāo)

模型的評價包括擬合精度和擬合優(yōu)度兩方面，本研究采用決定系數(shù)(R2)、均方根誤差(ERMS)、平均誤差(EM)和平均絕對相對誤差(EMA)4個模型評價指標(biāo)，依據(jù)十折交叉驗證，各指標(biāo)的平均值來進行模型評價，具體計算公式如下：

(2)

(3)

(4)

(5)

3 結(jié)果與分析

為了比較3種機器學(xué)習(xí)算法和傳統(tǒng)削度模型對杉木干形的擬合能力，采用評分法對模型進行排名[19]，指標(biāo)最優(yōu)的模型得分最高，最后每個模型把所有指標(biāo)的得分相加就是該模型的最終得分，結(jié)果見表2。各模型排名為人工神經(jīng)網(wǎng)絡(luò)、最近鄰法、傳統(tǒng)削度模型、隨機森林，4個模型的R2在0.97～0.99，ERMS在1左右。最近鄰法與人工神經(jīng)網(wǎng)絡(luò)的表現(xiàn)優(yōu)于其他模型，其R2均是0.99，其余指標(biāo)也都較小。人工神經(jīng)網(wǎng)絡(luò)模型的EM最小，但其EMA比最近鄰法的高了65.3%。最近鄰法是以K個最鄰近的目標(biāo)的平均值作為參考，因此，最近鄰法生成的不是擬合曲線，而是類似實際數(shù)據(jù)分布的散亂的點云圖，所以與平均回歸相比，其誤差往往更小。

表2 基于非線性回歸與機器學(xué)習(xí)算法的擬合結(jié)果(訓(xùn)練集)

為了檢驗各模型的泛化預(yù)測能力，利用不參與模型構(gòu)建的檢驗集對模型進行獨立性檢驗(表3)，模型在檢驗集上的表現(xiàn)略差于訓(xùn)練集，各模型的R2都有減小，但減小幅度很小，這與前人的結(jié)論一致[19]。人工神經(jīng)網(wǎng)絡(luò)模型的變化最小，說明人工神經(jīng)網(wǎng)絡(luò)模型的泛化預(yù)測能力更好，穩(wěn)定性更強。人工神經(jīng)網(wǎng)絡(luò)模型的R2最大，為0.98，其ERMS、EM和EMA均是各模型間最小。最近鄰法模型的EM比傳統(tǒng)削度模型大0.046，但其EMA與ERMS分別比傳統(tǒng)削度模型小了16.52%和2.53%。由于傳統(tǒng)削度方程是采用平滑的數(shù)學(xué)函數(shù)描述干形輪廓，而人工神經(jīng)網(wǎng)絡(luò)與最近鄰法算法直接受數(shù)據(jù)驅(qū)動，從數(shù)據(jù)中尋找規(guī)律，能最大限度逼近復(fù)雜的函數(shù)，所以精度一般會更高。

表3 基于非線性回歸與機器學(xué)習(xí)算法的獨立性檢驗(檢驗集)

各模型的殘差分布情況如圖2所示，整體來看，各個模型的殘差均勻分布在y=0直線的兩側(cè)，且都沒有出現(xiàn)極端殘差，殘差的波動范圍正常，檢驗集的殘差與訓(xùn)練集的殘差分布非常接近，說明沒有出現(xiàn)過擬合，這一現(xiàn)象在機器學(xué)習(xí)過程中很容易出現(xiàn)，通過合理地調(diào)參可以有效避免[17]。最近鄰法與隨機森林的殘差散點圖趨勢表現(xiàn)相似，當(dāng)di<2 cm時，殘差大都小于0，而當(dāng)di>25 cm時，殘差都大于0。說明這兩個模型對樹干基部有預(yù)測偏小的趨勢，而對樹干上部，有預(yù)測偏大的趨勢。但很明顯最近鄰法的殘差很小，精度高于隨機森林。盡管人工神經(jīng)網(wǎng)絡(luò)模型的擬合殘差表現(xiàn)不是最好，但它的預(yù)測殘差表現(xiàn)最優(yōu)，有最好的預(yù)測泛化能力，這在實際應(yīng)用中非常重要。

圖2 不同模型的殘差圖

各模型的相對偏差分布情況如圖3所示，可以看出，無論訓(xùn)練集還是檢驗集，各個模型的相對偏差分布都呈現(xiàn)單峰趨勢，在x=0處為峰值，說明絕大多數(shù)預(yù)測的相對偏差都為0，或是非常接近0。最近鄰法模型的峰值相差最大，相差0.15，隨機森林模型的差值為0.10，說明這兩個模型穩(wěn)定性不好，能較好地描述已知數(shù)據(jù)，但不能準(zhǔn)確地預(yù)測未知數(shù)據(jù)。人工神經(jīng)網(wǎng)絡(luò)模型的訓(xùn)練集和檢驗集的相對偏差分布幾乎重合，其峰值相差極小，模型的擬合結(jié)果和預(yù)測結(jié)果較理想。人工神經(jīng)網(wǎng)絡(luò)模型由于具有大量的“神經(jīng)元”，通過調(diào)整權(quán)重而不斷逼近[24]，采用誤差反向傳播算法訓(xùn)練的，能夠?qū)θ魏晤愋偷妮斎氡M可能映射到輸出。研究采用提前終止，沒有出現(xiàn)過擬合現(xiàn)象。

圖3 相對偏差的密度圖

各模型對干形的模擬精度隨著相對高度的變化而變化(表4)，在相對高度為0.1～0.6時，各模型的EMA較小。說明削度模型對樹干中部的預(yù)測精度更高，對樹干基部的和樹干上部的直徑的預(yù)測誤差較大，這與前人結(jié)論一致[18]，符合林木樹干的形態(tài)特征，樹干基部與上部的削度變化較大，而樹干中部類似圓柱體，更為規(guī)律。在相對高度≤0.1時，傳統(tǒng)削度模型對干形的模擬精度最高，而在相對高度為0.2～1.0時，人工神經(jīng)網(wǎng)絡(luò)模型與最近鄰法模型對干形的模擬精度較高。說明不同模型對樹干不同部位的模擬效果不同，如果多個模型互為補充，配合使用，其預(yù)測精度會更高。人工神經(jīng)網(wǎng)絡(luò)模型對樹干絕大部分的預(yù)測精度最高，可以較精確地模擬杉木干形，但對樹干基部的預(yù)測效果略差于傳統(tǒng)削度模型，若采用傳統(tǒng)削度模型作為補充可能會進一步提高預(yù)測精度。

表4 各模型對不同樹干部位預(yù)測的平均絕對誤差(EMA)

4 結(jié)論與討論

傳統(tǒng)回歸在處理復(fù)雜的非線性問題時，尤其當(dāng)變量個數(shù)較多、變量間關(guān)系復(fù)雜時，往往出現(xiàn)不收斂的情況[14，25-27]。樹木干形的影響因素很多，采用傳統(tǒng)的模型增加自變量會使模型結(jié)構(gòu)復(fù)雜，實際應(yīng)用不便。機器學(xué)習(xí)是一種能夠自主發(fā)現(xiàn)變量間的關(guān)系、給出模型結(jié)果的算法，而且可以處理有噪聲、缺失等低質(zhì)量數(shù)據(jù)。

本研究以我國重要的用材樹種的杉木為對象，引入機器學(xué)習(xí)，分別采用最近鄰法(KNN)、隨機森林(RF)和人工神經(jīng)網(wǎng)絡(luò)(ANN)3種算法對干形進行模擬，同時構(gòu)建了傳統(tǒng)非線性削度方程作為比較。結(jié)果表明：人工神經(jīng)網(wǎng)絡(luò)模型是4個模型間表現(xiàn)最優(yōu)的，R2最大，殘差范圍和殘差密度范圍最小，對樹干絕大部分的干形模擬都較為精確，特別是其預(yù)測泛化能力，優(yōu)于其他模型，這與Schikowski et al.的對巴西草原的黑松干形的研究結(jié)果一致[18]。人工神經(jīng)網(wǎng)絡(luò)模型之所以表現(xiàn)如此出色，是因為該算法能夠自動發(fā)掘“真實數(shù)據(jù)”之間的復(fù)雜關(guān)系，考慮變量之間的相互依賴性，這對傳統(tǒng)非線性回歸模型來說，是不容易實現(xiàn)的。其次表現(xiàn)較好的是最近鄰法模型，最近鄰法模型的擬合能力是4個模型間最好的，但是其預(yù)測能力不如人工神經(jīng)網(wǎng)絡(luò)，在Sanquetta et al.利用最近鄰法模型估測材積時也發(fā)現(xiàn)這一特點[17]。最近鄰法模型與傳統(tǒng)回歸模型之間的誤差非常相近，甚至在某些情況下，最近鄰法的誤差更小。最近鄰法模型能夠通過若干鄰近值進行預(yù)測[20]，不需要任何回歸擬合過程，具有簡單通用和較強的可塑性等優(yōu)點。隨機森林模型能夠有效處理生態(tài)領(lǐng)域的分類和回歸問題，具有很好的泛化能力和統(tǒng)計可靠性[22]。但本研究中，隨機森林對杉木不同高度處直徑的預(yù)測不高，尤其在檢驗集內(nèi)，該模型有對較小直徑的預(yù)測偏小，對較大直徑的預(yù)測偏大的趨勢。這種特殊趨勢是基于回歸樹模型特有的，在其他的研究中也有類似發(fā)現(xiàn)[17]。隨機森林表現(xiàn)不如人工神經(jīng)網(wǎng)絡(luò)，是內(nèi)在算法不同導(dǎo)致，隨機森林的每個節(jié)點由平均值給出，可能會低估較低值而高估較高值。

機器學(xué)習(xí)模型和傳統(tǒng)削度模型都表現(xiàn)出良好的擬合能力，訓(xùn)練集的決定系數(shù)(R2)均高于0.97，均方根誤差(ERMS)在1 cm左右波動。盡管傳統(tǒng)回歸能夠滿足干形預(yù)測精度，但是在模型選型和參數(shù)估計上仍是比較復(fù)雜。此外，最近鄰法模型與人工神經(jīng)網(wǎng)絡(luò)模型在沒有引入更多變量的前提下，擬合精度與預(yù)測精度均高于傳統(tǒng)削度模型，比傳統(tǒng)回歸能更精確地預(yù)測杉木沿樹干不同高度處的直徑，最近鄰法與人工神經(jīng)網(wǎng)絡(luò)更適合模擬杉木干形，具有更好的適用性。未來的研究可以基于機器學(xué)習(xí)算法研究其他森林參數(shù)，擴展其應(yīng)用范圍，以提高森林參數(shù)的估測精度。