應用地基激光雷達三維點云數據構建長白落葉松樹干削度方程1)

種雨絲 何培 張茲鵬 姜立春

(東北林業大學, 哈爾濱,150040)

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

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

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

1 研究方法

1.1 數據來源與處理

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

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

用FARO SCENE軟件完成同一樣地不同站點之間的拼接處理,獲得整個樣地的點云數據,并以“.xyz”格式將點云數據導出。將FARO SCENE軟件導出的數據(.xyz格式)加載到LiDAR360軟件中,完成重采樣、去噪、分離地面點,根據地面點歸一化以及點云分割的數據預處理流程,獲得樣地和單木三維點云。對單木點云提取胸徑、樹高、冠幅以及相對樹高(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塊樣地內共有落葉松75株,其中4株因樹木冠層遮擋嚴重,提取樹高與實際樹高存在較大差距,因此本研究最終選用71株落葉松,共計1 065對樹干相對高與直徑干形數據。樣木的胸徑為12.8～28.4 cm,平均胸徑為22.3 cm,標準差為4.2 cm;樹高為16.4～24.5 m,平均樹高為20.8 m,標準差為1.3 m。通過畫散點圖的方式剔除落葉松干形數據異常點(見圖1)。

圖1 數據處理前后相對高和相對直徑散點圖

1.2 基礎削度方程選擇

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

Kozak(1969)簡單削度方程：

MB(1976)分段削度方程：

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

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

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

Bi(2000)三角函數方程：

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

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

Lynch(2017)簡單削度方程：

1.3 分位數回歸

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

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

1.4 廣義加性模型

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

式中：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(·)表示不同的光滑樣條函數。

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

1.5 模型評價和檢驗的方法

模型評價指標包括確定系數(R2)、平均絕對誤差(MAE)、相對誤差(MPB)、均方根誤差(RMSE)和多重共線性指標條件數(CN)。當CN<30時,模型自變量不存在多重共線性;當30100時,模型自變量存在嚴重的共線性[22]。各統計量的計算表達式如下：

模型檢驗采用留一交叉驗證法進行驗證。該方法每次只選擇一個樣本數據作為驗證數據,其余數據作為擬合數據,依此循環,直至每個樣本數據都被作為一次驗證數據結束,擬合和驗證的次數等于樣本數。

1.6 模型相對排序值的計算方法

當多個模型的評價指標不統一時,使用相對排序法計算各模型的相對排序值[23],計算方法如下：

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

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

2 結果與分析

2.1 基礎模型的擬合結果

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

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

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

表1 落葉松削度方程的參數估計值

2.2 分位數模型的預測結果

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

式中：dτ為各分位點處的直徑,b1τ～b7τ為模型參數。

通過對分位數模型的9個分位點的擬合,9個分位點中,分位數τ=0.4、0.5、0.6分位點處的擬合精度較高。

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

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

表3 落葉松分位數回歸模型參數估計值及擬合效果

2.3 廣義加性模型的預測結果

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

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

由表4可知,模型4的擬合精度優于模型1、模型3;模型5的擬合精度高于模型2,與模型3的擬合精度基本相同;模型6與模型3的擬合精度相同。各模型的擬合效果從高到低可排序為：模型4、模型5、模型3(模型6)、模型2、模型1。

表4 落葉松非參數模型擬合效果

2.4 模型預測精度的檢驗

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

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

表5 落葉松削度方程檢驗結果

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

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

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

3 討論

選取了9個基礎削度方程進行比較,這些方程分別代表簡單、分段、三角函數和可變指數削度方程等4種類型。可變指數Bi(2000)模型的擬合效果最好,除Kozak(1969)模型外,其余7個模型的擬合統計量比較接近。但研究發現Kozak(1988)模型存在比較嚴重的多重共線性問題。鄒茂勝等[24]研究認為Kozak(1988)模型優于Kozak(2004)模型和曾偉生(1997)模型。

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

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

4 結論

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