種雨絲 何培 張茲鵬 姜立春

(東北林業(yè)大學(xué), 哈爾濱,150040)

立木材積的準(zhǔn)確估計(jì)是森林有效管理和經(jīng)營的重要環(huán)節(jié)之一。目前在國內(nèi)外林業(yè)生產(chǎn)實(shí)踐中,都使用削度方程估計(jì)立木材積和商品材積。通過伐倒木分段測量獲取干形數(shù)據(jù)構(gòu)建削度方程的方法耗時(shí)費(fèi)力,具有破壞性,浪費(fèi)木材資源[1]。隨著測量技術(shù)的不斷發(fā)展,采用非破壞性測量的地基激光雷達(dá)(TLS)技術(shù)獲取干形數(shù)據(jù)[2-3],利用TLS技術(shù)調(diào)查森林資源,從確定單木位置、提取單木胸徑和樹高到預(yù)測樹干直徑、估算材積、監(jiān)測森林的動(dòng)態(tài)變化等獲得高質(zhì)量林木影像信息,提取林木豐富的水平和垂直信息[4-8]。

削度方程是以樹木胸徑、樹高等單木特征因子作為自變量,以單木任意高度處直徑作為因變量的函數(shù)[9]。削度方程在擬合時(shí),通常使用非線性最小二乘法得到參數(shù)估計(jì),該方法需要假定模型誤差服從獨(dú)立正態(tài)分布、非共線性和無自相關(guān)等[10],由于干形數(shù)據(jù)存在自相關(guān)性的問題,采用非線性回歸需要消除自相關(guān)性,增加了模型參數(shù)估計(jì)的復(fù)雜程度。近年來,應(yīng)用分位數(shù)回歸和廣義加性模型方法建立了樹高-胸徑模型、削度方程和枝下高模型等[11],由于分位數(shù)回歸方法由多個(gè)分位點(diǎn)構(gòu)建多個(gè)模型,靈活地反映林木數(shù)據(jù)的變化規(guī)律。廣義加性模型以加性項(xiàng)組合的形式表達(dá)自變量與因變量之間的關(guān)系,在國內(nèi)外林業(yè)生產(chǎn)研究領(lǐng)域得到廣泛應(yīng)用,并取得了較好的預(yù)測效果[12-14]。

長白落葉松(Larixolgensis)是我國東北林區(qū)的主要造林樹種。以吉林省東豐縣一面山林場和楊木林林場長白落葉松人工林為研究對象,使用TLS點(diǎn)云數(shù)據(jù)提取落葉松人工林單木干形數(shù)據(jù)。根據(jù)提取的數(shù)據(jù),選取9個(gè)基礎(chǔ)削度方程分別代表4種類型(簡單、可變指數(shù)、三角函數(shù)和分段),利用非線性回歸確定基礎(chǔ)最優(yōu)模型,構(gòu)建分位數(shù)回歸模型和廣義加性模型,綜合分析削度方程構(gòu)建方法優(yōu)劣,為落葉松干形精準(zhǔn)預(yù)測提供技術(shù)支持。

1 研究方法

1.1 數(shù)據(jù)來源與處理

本研究落葉松干形數(shù)據(jù)采集地點(diǎn)為吉林省東豐縣一面山林場和楊木林林場,地貌屬于低山丘陵區(qū),坡度平緩。該區(qū)域?qū)儆诩撅L(fēng)區(qū)中溫帶濕潤氣候,主要樹種為落葉松、樟子松和紅松等。于2021年7月在該區(qū)域內(nèi)共選取了3塊落葉松人工林樣地,樣地面積為0.04 hm2,樣地內(nèi)落葉松年齡為48～51 a,對樣地進(jìn)行每木檢尺,測量胸徑、樹高、冠幅等單木特征因子。

TLS數(shù)據(jù)由FARO Focus 3D X 330掃描儀對樣地進(jìn)行掃描獲得。依據(jù)各樣地的坡度、郁閉度等環(huán)境因子調(diào)查情況,在各樣地外圍及內(nèi)部分別布設(shè)5個(gè)站點(diǎn),并在站點(diǎn)范圍內(nèi)設(shè)置5個(gè)靶球,保證每兩個(gè)站點(diǎn)間可匹配到最少兩對靶球,便于后期每個(gè)樣地內(nèi)不同站點(diǎn)的點(diǎn)云數(shù)據(jù)拼接。

用FARO SCENE軟件完成同一樣地不同站點(diǎn)之間的拼接處理,獲得整個(gè)樣地的點(diǎn)云數(shù)據(jù),并以“.xyz”格式將點(diǎn)云數(shù)據(jù)導(dǎo)出。將FARO SCENE軟件導(dǎo)出的數(shù)據(jù)(.xyz格式)加載到LiDAR360軟件中,完成重采樣、去噪、分離地面點(diǎn),根據(jù)地面點(diǎn)歸一化以及點(diǎn)云分割的數(shù)據(jù)預(yù)處理流程,獲得樣地和單木三維點(diǎn)云。對單木點(diǎn)云提取胸徑、樹高、冠幅以及相對樹高(0、0.02、0.04、0.06、0.08、0.10、0.15、0.20、0.30、0.40、0.50、0.60、0.70、0.80、0.90)處的樹木直徑。

3塊樣地內(nèi)共有落葉松75株,其中4株因樹木冠層遮擋嚴(yán)重,提取樹高與實(shí)際樹高存在較大差距,因此本研究最終選用71株落葉松,共計(jì)1 065對樹干相對高與直徑干形數(shù)據(jù)。樣木的胸徑為12.8～28.4 cm,平均胸徑為22.3 cm,標(biāo)準(zhǔn)差為4.2 cm;樹高為16.4～24.5 m,平均樹高為20.8 m,標(biāo)準(zhǔn)差為1.3 m。通過畫散點(diǎn)圖的方式剔除落葉松干形數(shù)據(jù)異常點(diǎn)(見圖1)。

圖1 數(shù)據(jù)處理前后相對高和相對直徑散點(diǎn)圖

1.2 基礎(chǔ)削度方程選擇

削度方程可大致分為簡單模型、分段模型、三角函數(shù)和可變指數(shù)模型。簡單模型通過使用基本初等函數(shù)模擬樹干形狀,模型容易構(gòu)建,計(jì)算簡便;分段模型借助拐點(diǎn)將樹干劃分不同的形狀,相較于簡單模型可提高干形預(yù)測精度;三角函數(shù)和可變指數(shù)模型則可通過變化的指數(shù)描述同一樹干不同高度的樹干干形變化[9]。本研究選用了Kozak(1969)、Lynch(2017)簡單削度方程、Max-Burkhart(1976)分段削度方程(簡稱MB(1976))、Bi(2000)三角函數(shù)方程和5個(gè)可變指數(shù)削度方程作為本研究的基礎(chǔ)方程[15-19],模型具體形式如下：

Kozak(1969)簡單削度方程：

MB(1976)分段削度方程：

Kozak(1988)可變指數(shù)削度方程：

Kozak(1994)可變指數(shù)削度方程：

曾偉生(1997)可變指數(shù)削度方程：

Bi(2000)三角函數(shù)方程：

Kozak(2001)可變指數(shù)削度方程：

Kozak(2004)可變指數(shù)削度方程：

Lynch(2017)簡單削度方程：

1.3 分位數(shù)回歸

與傳統(tǒng)回歸法相比,分位數(shù)回歸靈活性強(qiáng),對數(shù)據(jù)分布沒有要求。通過擬合任一分位點(diǎn)的數(shù)據(jù),可估計(jì)不同分位點(diǎn)下的回歸參數(shù)。主要研究不同分位點(diǎn)處自變量對因變量的變化趨勢或研究不同分位點(diǎn)處哪些自變量是主要影響因素[20]。在樹木干形預(yù)測方面,可以通過對比不同分位點(diǎn)回歸模型對樹干不同高度處直徑的預(yù)測結(jié)果,分析不同分位點(diǎn)模型的干形預(yù)測精度。對分位數(shù)回歸模型進(jìn)行參數(shù)估計(jì)時(shí),損失函數(shù)(Lmin)采用加權(quán)的最小二乘法,即通過最小化加權(quán)離差的絕對值之和獲得方程參數(shù)。

本研究采用quantreg包的nlrq()函數(shù)擬合各分位數(shù)模型?；A(chǔ)削度方程利用qnls()函數(shù)進(jìn)行擬合,由于干形數(shù)據(jù)存在自相關(guān)特征,在基礎(chǔ)削度方程擬合時(shí)引入CAR(l)函數(shù)來消除變量的自相關(guān)性[21]。

1.4 廣義加性模型

廣義加性模型是一種多元回歸的非參數(shù)化平滑回歸形式,用光滑函數(shù)代替線性函數(shù),在建模之前不需要分析因變量與自變量之間的關(guān)系。廣義加性模型由不同的光滑樣條函數(shù)組成,模型中自變量和因變量之間可以為任意線性或非線性的關(guān)系,引入的光滑樣條函數(shù)可增加對因變量預(yù)測的準(zhǔn)確性,靈活地探測數(shù)據(jù)間的復(fù)雜關(guān)系,解釋因變量隨自變量變化規(guī)律。構(gòu)建削度廣義加性模型,可選取樹干直徑或樹干相對直徑作為因變量,選取胸徑、樹高和相對樹高作為自變量,同時(shí)對自變量進(jìn)行平方、根號(hào)等變換。廣義加削度模型如下：

式中：dij表示第i株樹第j位置處的直徑;Di表示第i株樹的胸徑;dij/Di表示第i株樹的相對直徑;Hi表示第i株樹的樹高;hij表示第i株樹第j位置處的高度;hij/Hi表示第i株樹的相對高度;i=1、2、3、…、n,j=1、2、3、…、n。α表示截距;f(·)表示不同的光滑樣條函數(shù)。

利用R軟件mgcv包的gamm()函數(shù)對廣義加模型進(jìn)行擬合,對不同的加性項(xiàng)選擇三次回歸樣條函數(shù)(CR)、B-樣條函數(shù)(BS)、薄板回歸樣條函數(shù)(TP)、P-樣條函數(shù)(PS)、Duchon樣條函數(shù)(DS)和高斯過程平滑樣條函數(shù)(GP)。

1.5 模型評價(jià)和檢驗(yàn)的方法

模型評價(jià)指標(biāo)包括確定系數(shù)(R2)、平均絕對誤差(MAE)、相對誤差(MPB)、均方根誤差(RMSE)和多重共線性指標(biāo)條件數(shù)(CN)。當(dāng)CN<30時(shí),模型自變量不存在多重共線性;當(dāng)30100時(shí),模型自變量存在嚴(yán)重的共線性[22]。各統(tǒng)計(jì)量的計(jì)算表達(dá)式如下：

模型檢驗(yàn)采用留一交叉驗(yàn)證法進(jìn)行驗(yàn)證。該方法每次只選擇一個(gè)樣本數(shù)據(jù)作為驗(yàn)證數(shù)據(jù),其余數(shù)據(jù)作為擬合數(shù)據(jù),依此循環(huán),直至每個(gè)樣本數(shù)據(jù)都被作為一次驗(yàn)證數(shù)據(jù)結(jié)束,擬合和驗(yàn)證的次數(shù)等于樣本數(shù)。

1.6 模型相對排序值的計(jì)算方法

當(dāng)多個(gè)模型的評價(jià)指標(biāo)不統(tǒng)一時(shí),使用相對排序法計(jì)算各模型的相對排序值[23],計(jì)算方法如下：

式中：Ri為模型i的相對排序值(i=1、2、…、a),Ri值越小表示該模型的預(yù)測效果越好;a為參與排序的模型數(shù);Si為模型i的評價(jià)指標(biāo)值;Sx和Sn分別為Si的最大值和最小值。

在此排序方法中,最優(yōu)模型排序值為1,最差模型的排序值為a,其余模型排序值為1～a。確定系數(shù)(R2)值越大,表示模型的擬合度越好,而平均絕對誤差(MAE)、相對誤差(MPB)和均方根誤差(RMSE)越小模型精度越高。計(jì)算確定系數(shù)(R2)的相對排序值時(shí),Sx和Sn分別表示R2的最小值和最大值。最終以確定系數(shù)(R2)、平均絕對誤差(MAE)、相對誤差(MPB)、均方根誤差(RMSE)等4個(gè)指標(biāo)的平均值作為該模型的相對排序值,排序值最小的模型即為最優(yōu)模型。

2 結(jié)果與分析

2.1 基礎(chǔ)模型的擬合結(jié)果

基于TLS落葉松提取的干形數(shù)據(jù),利用R軟件擬合9個(gè)基礎(chǔ)削度方程,在擬合過程中,若模型中存在不顯著的參數(shù)(P>0.05),就刪除不顯著參數(shù)后再次擬合,直至模型中所有參數(shù)顯著。

由表1可知,MB(1976)模型的參數(shù)估計(jì)值(b5和b6)表明了影響落葉松干形變化的兩個(gè)拐點(diǎn)分別位于相對高度0.73和0.07處。

由表2可知,Bi(2000)模型的擬合效果最好,其次是Kozak(1988)和Kozak(2004)模型,而Kozak(1969)模型的擬合效果相對較差。MB(1976)、Kozak(1988)、Kozak(1994)、曾偉生(1997)和Lynch(2017)模型的多重共線性指標(biāo)條件數(shù)(CN)大于100,說明各模型自變量存在嚴(yán)重的多重共線性;Bi(2000)、Kozak(2001)和Kozak(2004)模型的CN值為30～100,各模型自變量多重共線性在可接受范圍內(nèi);Kozak(1969)模型的CN值小于30,模型自變量不存在共線性。因此,基礎(chǔ)削度方程中Bi(2000)模型的擬合效果最好,CN值也在可接受范圍內(nèi)。

表1 落葉松削度方程的參數(shù)估計(jì)值

2.2 分位數(shù)模型的預(yù)測結(jié)果

根據(jù)擬合精度最高的Bi(2000)模型,構(gòu)建9個(gè)分位數(shù)(τ=0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9)削度方程形式如下：

式中：dτ為各分位點(diǎn)處的直徑,b1τ～b7τ為模型參數(shù)。

通過對分位數(shù)模型的9個(gè)分位點(diǎn)的擬合,9個(gè)分位點(diǎn)中,分位數(shù)τ=0.4、0.5、0.6分位點(diǎn)處的擬合精度較高。

由表3可知,不同分位點(diǎn)下模型的參數(shù)估計(jì)值均不相同,說明在不同分位點(diǎn)下,落葉松干形的擬合結(jié)果存在一定的差異。隨著分位點(diǎn)數(shù)的增加,分位數(shù)模型的平均絕對誤差、相對誤差、均方根誤差值呈現(xiàn)先減小后增加的特征,確定系數(shù)(R2)呈現(xiàn)先增加后減小的特征,在分位點(diǎn)τ=0.5處,模型擬合效果最好(平均絕對誤差為0.847、相對誤差為4.510、均方根誤差為1.211、確定系數(shù)(R2)為0.963),且與非線性回歸擬合效果基本相近。

表2 落葉松削度模型擬合效果

表3 落葉松分位數(shù)回歸模型參數(shù)估計(jì)值及擬合效果

2.3 廣義加性模型的預(yù)測結(jié)果

對廣義加性模型中的非參數(shù)平滑項(xiàng)f(·),使用6種樣條函數(shù)(TP、CR、BS、PS、DS和GP)擬合。經(jīng)對比分析各模型的擬合統(tǒng)計(jì)量,當(dāng)不同加性項(xiàng)使用同一樣條函數(shù)時(shí),6個(gè)樣條函數(shù)中,CR、TP和DS樣條函數(shù)的擬合效果較好,分別記作模型1、模型2和模型3。

根據(jù)擬合精度最高的DS模型(模型3),胸徑和樹高所在加性項(xiàng)依然使用DS樣條函數(shù),相對樹高所在加性項(xiàng)依次使用其余5個(gè)樣條函數(shù)(TP、CR、BS、PS和GP)重新進(jìn)行擬合。經(jīng)比較分析,DS樣條函數(shù)分別與CR、TP和GP樣條函數(shù)組合后的模型擬合效果較好,分別記作模型4、模型5和模型6。

由表4可知,模型4的擬合精度優(yōu)于模型1、模型3;模型5的擬合精度高于模型2,與模型3的擬合精度基本相同;模型6與模型3的擬合精度相同。各模型的擬合效果從高到低可排序?yàn)椋耗Ｐ?、模型5、模型3(模型6)、模型2、模型1。

表4 落葉松非參數(shù)模型擬合效果

2.4 模型預(yù)測精度的檢驗(yàn)

根據(jù)3種建模方法(基礎(chǔ)模型、分位數(shù)模型和廣義加性模型)取得的最優(yōu)模型,使用留一交叉驗(yàn)證法對Bi(2000)非線性模型、Bi(2000)分位數(shù)模型和廣義加性模型進(jìn)行檢驗(yàn)。

由表5可知,根據(jù)分位數(shù)回歸構(gòu)建的Bi(2000)模型的預(yù)測精度略高于傳統(tǒng)非線性回歸模型,廣義加性模型的平均絕對誤差、相對誤差、均方根誤差、排序值均為最小,確定系數(shù)(R2)為最大,預(yù)測精度最高。

表5 落葉松削度方程檢驗(yàn)結(jié)果

為進(jìn)一步比較各模型(ONLS、QR和GAM)在樹干不同高度處的預(yù)測能力,計(jì)算了各模型在9個(gè)相對樹高處的平均絕對誤差和均方根誤差值。

由圖2可知,根據(jù)平均絕對誤差和均方根誤差值隨相對樹高的變化情況,非線性回歸模型和分位數(shù)回歸模型在不同相對樹高處的預(yù)測精度基本一致。在樹干相對樹高10%以下,廣義加性模型的平均絕對誤差值略大于非線性回歸模型和分位數(shù)回歸模型,均方根誤差值與非線性回歸模型和分位數(shù)回歸模型基本相同;在相對樹高10%～80%范圍內(nèi),廣義加性模型的平均絕對誤差和均方根誤差值均小于非線性回歸模型和分位數(shù)回歸模型;在相對樹高80%～90%處,廣義加性模型的平均絕對誤差值處于非線性回歸模型和分位數(shù)回歸模型之間,均方根誤差值小于非線性回歸模型和分位數(shù)回歸模型。因此,廣義加性模型在樹干大部分的預(yù)測精度都優(yōu)于非線性回歸模型和分位數(shù)回歸模型。

圖2 落葉松削度方程在樹干不同高度處的預(yù)測能力

3 討論

選取了9個(gè)基礎(chǔ)削度方程進(jìn)行比較,這些方程分別代表簡單、分段、三角函數(shù)和可變指數(shù)削度方程等4種類型?？勺冎笖?shù)Bi(2000)模型的擬合效果最好,除Kozak(1969)模型外,其余7個(gè)模型的擬合統(tǒng)計(jì)量比較接近。但研究發(fā)現(xiàn)Kozak(1988)模型存在比較嚴(yán)重的多重共線性問題。鄒茂勝等[24]研究認(rèn)為Kozak(1988)模型優(yōu)于Kozak(2004)模型和曾偉生(1997)模型。

在最優(yōu)基礎(chǔ)模型的基礎(chǔ)上構(gòu)建了9個(gè)分位點(diǎn)(τ=0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9)的分位數(shù)回歸模型。中位數(shù)(τ=0.5)模型的擬合效果最好。通過比較2個(gè)模型的預(yù)測精度發(fā)現(xiàn)非線性模型和中位數(shù)模型都可以較好的描述干形變化,但中位數(shù)回歸模型沒有假設(shè)要求,建議使用分位數(shù)回歸模型。但兩個(gè)參數(shù)模型在樹干大部分范圍內(nèi)的預(yù)測精度都略低于廣義加性模型,因此,3種建模方法中,推薦使用廣義加性模型預(yù)測樹干直徑。

地基激光雷達(dá)技術(shù)作為干形數(shù)據(jù)獲取的一種手段,由于受林分密度和枝條生長狀況等因素的影響,該技術(shù)也存在一定的局限性。對樹冠體積較大的單木,在樹干相對高80%以上,樹干直徑信息量低于中下部,需手動(dòng)提取樹干直徑,提取精度相對較低,會(huì)對干形預(yù)測精度有影響。在設(shè)置多站點(diǎn)掃描的情況下,可通過多數(shù)據(jù)源融合的方法補(bǔ)充樹干上部點(diǎn)云信息,提高樹干直徑提取精度,進(jìn)而提高干形預(yù)測精度。

4 結(jié)論

本研究基于地基激光雷達(dá)技術(shù)獲得了落葉松人工林樣地的三維點(diǎn)云信息,提取了落葉松的樹高、胸徑和不同相對樹高處的樹干直徑等單木信息,并對比分析了傳統(tǒng)削度方程、分位數(shù)回歸削度方程和廣義加性削度方程的預(yù)測精度。廣義加性模型的預(yù)測精度優(yōu)于傳統(tǒng)削度模型和分位數(shù)回歸模型,其中以相對直徑為因變量,以胸徑的平方、樹高和相對樹高的平方根為自變量,使用組合樣條函數(shù)的廣義加性削度方程更適合描述該區(qū)域的落葉松干形變化。