馬 驤，姜文龍，閆保銀，金 悅，崔 立

(1.常州市金壇區(qū)林業(yè)技術(shù)指導(dǎo)站，江蘇 金壇 213200； 2.南京國(guó)圖信息產(chǎn)業(yè)有限公司，江蘇 南京 210037)

森林是地球生態(tài)系統(tǒng)的重要組成部分，在給人類帶來(lái)林副產(chǎn)品產(chǎn)生經(jīng)濟(jì)效益的同時(shí)，還具備涵養(yǎng)水源、調(diào)節(jié)區(qū)域氣候條件，實(shí)現(xiàn)碳平衡、碳循環(huán)等功能，產(chǎn)生生態(tài)效益。在IPCC發(fā)布的《全球溫控1.5 ℃特別報(bào)告》中明確指出，“林業(yè)是實(shí)現(xiàn)1.5 ℃溫控目標(biāo)的重要路徑之一，通過(guò)林業(yè)固碳除碳等技術(shù)，可實(shí)現(xiàn)溫室氣體的大幅減少，在促進(jìn)森林生態(tài)系統(tǒng)恢復(fù)，提高基于生態(tài)系統(tǒng)和社區(qū)的適應(yīng)性以及濕地管理方面潛力巨大?！痹诖吮尘跋?，如何快速準(zhǔn)確地獲取森林資源信息，實(shí)現(xiàn)動(dòng)態(tài)監(jiān)測(cè)，為森林保護(hù)、森林經(jīng)營(yíng)等提供依據(jù)，成為當(dāng)前森林資源管理的重要任務(wù)之一。

激光雷達(dá)(LiDAR)的出現(xiàn)為森林資源調(diào)查方式的轉(zhuǎn)變提供契機(jī)。作為一種主動(dòng)遙感技術(shù)，具有精度高、成本低等優(yōu)勢(shì)，早在20世紀(jì)80年代已用于林業(yè)研究[1]。1984年，Nelson等[2]通過(guò)研究激光雷達(dá)與樹木冠層郁閉度關(guān)系，指出兩者間相關(guān)性顯著；MacLean等[3]通過(guò)研究樹冠垂直剖面與森林蓄積量相關(guān)性，指出2者的對(duì)數(shù)呈線性相關(guān)。進(jìn)入21世紀(jì)后，越來(lái)越多學(xué)者也相繼證明激光雷達(dá)對(duì)于估測(cè)樹高、胸高斷面積、生物量等森林參數(shù)準(zhǔn)確性較高[4-6]。Morsdorf等[7]采用林木冠幅分割算法對(duì)瑞士國(guó)家森林公園針葉林進(jìn)行分割，再生成種子點(diǎn)進(jìn)行聚類分析，獲得樹冠直徑后按圓面積公式反演了林木胸徑；劉清旺等[8-9]基于機(jī)載激光雷達(dá)技術(shù)，使用多元線性回歸、隨機(jī)森林模型等方法，對(duì)森林平均高、生物量、郁閉度等參數(shù)進(jìn)行提取分析，其反演結(jié)果較光學(xué)遙感而言精度更高；Zhao等[10]采用2種尺度不變模型估測(cè)生物量，結(jié)果表明該2種方法R2范圍在0.80—0.95之間，可靠性較高。本文以江蘇省常州市金壇區(qū)碳匯計(jì)量山區(qū)監(jiān)測(cè)點(diǎn)中的國(guó)外松片林作為研究對(duì)象，進(jìn)行機(jī)載激光雷達(dá)蓄積量反演試驗(yàn)，評(píng)價(jià)這一技術(shù)在江蘇省蘇南地區(qū)國(guó)外松林的適用性，以及在蓄積量反演方面的可能性，為今后森林資源調(diào)查新技術(shù)應(yīng)用提供參考。

1 研究區(qū)概況與數(shù)據(jù)獲取

1.1 研究區(qū)概況

研究區(qū)位于江蘇省常州市金壇區(qū)4 km×4 km碳匯計(jì)量山區(qū)監(jiān)測(cè)點(diǎn)范圍內(nèi)，國(guó)家2000大地坐標(biāo)系(China Geodetic Coordinate System 2000, CGCS2000)范圍坐標(biāo)：西北至(20 718 118.46，3 517 999.08)，西南至(20 718 118.46，3 513 999.07)，東南至(20 722 118.48，3 513 999.07)，東北至(2 0722 118.48，3 517 999.08)。區(qū)內(nèi)地貌以丘陵為主，為北亞熱帶季風(fēng)區(qū)，年均氣溫15.3 ℃四季分明，雨量充沛，年降水量1 063.5 mm。日照充足，日照率46%，無(wú)霜期228 d，年平均濕度為78%。

研究區(qū)面積為16 km2，根據(jù)2018年碳匯計(jì)量監(jiān)測(cè)成果，區(qū)內(nèi)森林總蓄積量為48 486.04 m3，分布樹種主要有濕地松(PinuselliottiiEngelmann)、杉木(Cunninghamialanceolata)、香樟(Cinnamomumcamphora)、女貞(Ligustrumlucidum)、柳杉(Cryptomeriafortunei)等。本次研究對(duì)象位于4 km×4 km碳匯計(jì)量山區(qū)監(jiān)測(cè)點(diǎn)(見(jiàn)圖1)，77，250號(hào)國(guó)外松小班范圍內(nèi)。其中77號(hào)小班用于激光雷達(dá)森林參數(shù)反演，250號(hào)小班用于反演數(shù)據(jù)驗(yàn)證。

圖1 試驗(yàn)區(qū)范圍

1.2 LiDAR數(shù)據(jù)獲取

分別于2021年4月16日、5月12日，利用激光雷達(dá)在4 km×4 km碳匯計(jì)量山區(qū)監(jiān)測(cè)點(diǎn)，77，250號(hào)小班范圍內(nèi)，進(jìn)行機(jī)載激光雷達(dá)進(jìn)行試驗(yàn)數(shù)據(jù)采集。使用大疆經(jīng)緯M300 RTK作為飛行平臺(tái)，飛行航速為6 m/s，總覆蓋面積約為10 hm2，飛行時(shí)天氣晴朗，少量高云分布，對(duì)于LiDAR飛行影響較少。

在本次飛行中，激光雷達(dá)系統(tǒng)采用大疆禪思L1(DJIL1)，集成Livox激光雷達(dá)模塊、高精度慣導(dǎo)、測(cè)繪相機(jī)、三軸云臺(tái)等模塊。飛行時(shí)設(shè)置航高80 m，旁向重疊度50%。鏡頭角度-90°，激光掃描方式為三回波重復(fù)掃描，類似于傳統(tǒng)線掃激光雷達(dá)，能獲取更均勻、精度更高的掃描結(jié)果，也更穩(wěn)定。

1.3 地面樣地?cái)?shù)據(jù)獲取

2021年5月對(duì)77，250號(hào)小班進(jìn)行外業(yè)調(diào)查。

本次地面實(shí)測(cè)調(diào)查分為大樣地與小樣地。大樣地以77號(hào)小班范圍為界，在其中按照20 m×20 m劃分小樣地，共計(jì)34塊。在樣地范圍內(nèi)記錄樹種、胸徑、樹高等參數(shù)，所得數(shù)據(jù)將以小樣地為單位進(jìn)行林分參數(shù)反演。

在250號(hào)小班內(nèi)，按照《森林資源規(guī)劃設(shè)計(jì)調(diào)查技術(shù)規(guī)程》[11]中小班測(cè)樹因子調(diào)查的相關(guān)要求，深入小班內(nèi)部，選擇具有代表性的調(diào)查點(diǎn)，以調(diào)查點(diǎn)為中心布設(shè)樣地并每木檢尺，記錄樹種、胸徑、樹高、郁閉度等測(cè)樹因子，按照《江蘇省主要樹種一元材積表》計(jì)算公頃蓄積，所得數(shù)據(jù)用于林分參數(shù)反演的數(shù)據(jù)可靠性驗(yàn)證。

圖2 樣地設(shè)置

2 研究方法

2.1 樣地?cái)?shù)據(jù)處理

根據(jù)樣地內(nèi)每木檢尺獲得的胸徑，對(duì)照《江蘇省主要樹種一元材積表》，查詢樹種為“馬1”一列，以及胸徑對(duì)應(yīng)的徑階，填寫單株立木的材積，最后再統(tǒng)一求和，獲得樣地內(nèi)總蓄積量。計(jì)算公式為

M樹種=V樹種/ ΣAi；

M全小班=∑(M樹種×A)。

式中，A：小班面積；Ai:第i塊小樣地面積；V樹種：某樹種各樣方蓄積之和；M樹種：某樹種小班每公頃蓄積；M全小班：全小班蓄積。

對(duì)于77號(hào)小班，其范圍內(nèi)的34塊小樣地需逐塊計(jì)算平均胸徑，用于后期反演模型的構(gòu)建，250號(hào)小班需求得小班平均胸徑、每公頃蓄積量，用于反演結(jié)果驗(yàn)證。

2.2 LiDAR數(shù)據(jù)處理

在林分反演模型構(gòu)建之前，首先需將點(diǎn)云去噪，其次對(duì)點(diǎn)云進(jìn)行地面點(diǎn)分類，以此分離地面點(diǎn)與地物點(diǎn)，然后分別生成數(shù)字地面模型(DTM)、數(shù)字表面模型(DSM)和冠層高度模型(CHM)，并基于CHM生成種子點(diǎn)，在點(diǎn)云歸一化處理(NPC)，消除地形影響后，最后進(jìn)行單木分割。本次對(duì)于LiDAR數(shù)據(jù)的處理均在LiDAR 360軟件中完成。

地面點(diǎn)分類是點(diǎn)云數(shù)據(jù)處理的基礎(chǔ)操作，本文采用改進(jìn)的漸進(jìn)加密三角網(wǎng)濾波算法(IPTD)[12]分類地面點(diǎn)。其原理是：對(duì)于不含建筑物的點(diǎn)云，以默認(rèn)值作為格網(wǎng)大小，取格網(wǎng)內(nèi)最低點(diǎn)作為起始種子點(diǎn)，構(gòu)建初始三角網(wǎng)。遍歷所有待分類點(diǎn)后，查詢各點(diǎn)水平投影所落入的三角形，并計(jì)算點(diǎn)到三角形的距離d及點(diǎn)到三角形3個(gè)頂點(diǎn)與三角形所在平面所成角度的最大值，分別與迭代距離與迭代角度進(jìn)行比較，如果小于對(duì)應(yīng)閾值，則將此點(diǎn)判定為地面點(diǎn)，并加入三角網(wǎng)中。重復(fù)此過(guò)程，直至所有地面點(diǎn)分類完畢。

地面點(diǎn)分類后，需經(jīng)插值處理生成DSM，DTM。LiDAR 360軟件中提供了反距離權(quán)重插值(IDW)、克里金插值(Kriging)和不規(guī)則三角網(wǎng)插值(TIN)方法。本文采用TIN插值法，分別生成DSM和DTM，兩者作差后生成CHM。利用生成的DTM高程對(duì)點(diǎn)云作歸一化處理，消除地形影響。設(shè)置地面高度、最小樹高再進(jìn)行單木分割。

2.3 基于單木分割的反演模型構(gòu)建

胸徑是計(jì)算林木蓄積量的重要參數(shù)之一，若能準(zhǔn)確獲取胸徑值，則可根據(jù)地方相應(yīng)的一元材積表求算蓄積量。

樹高、冠幅等參數(shù)可通過(guò)機(jī)載激光雷達(dá)直接測(cè)得，而胸徑則需通過(guò)建立回歸方程的方式間接獲取[13]。本文擬參考劉清旺[14]、何祺勝等[15]的研究方法，以34塊小樣地為單位，統(tǒng)計(jì)各樣地內(nèi)LiDAR估測(cè)的單木平均樹高、平均冠幅等參數(shù)，作為自變量，實(shí)測(cè)平均胸徑作為因變量，建立實(shí)測(cè)胸徑與LiDAR估測(cè)參數(shù)的一元回歸方程，間接估測(cè)胸徑。共設(shè)置3組：

(1)實(shí)測(cè)胸徑與LiDAR樹高

(2)實(shí)測(cè)胸徑與LiDAR冠幅

(3)實(shí)測(cè)胸徑的自然對(duì)數(shù)與LiDAR樹高的自然對(duì)數(shù)

一元回歸方程建立后，根據(jù)相關(guān)系數(shù)，確定最優(yōu)估測(cè)方程。

2.4 反演模型可靠性評(píng)價(jià)

反演模型采用擬合優(yōu)度(R2)、均方根誤差(RMSE)和相對(duì)均方根誤差(rRMSE)進(jìn)行可靠性評(píng)價(jià)。

R2表示模型的擬合優(yōu)度，一般R2值越高，表明該模型的擬合度就越高。計(jì)算公式如下：

RMSE表示實(shí)際值與預(yù)測(cè)值之間的均方根,誤差越小，則方程擬合效果越好，計(jì)算公式為

rRMSE為相對(duì)均方根誤差，其值應(yīng)當(dāng)越小越好，計(jì)算公式為

3 結(jié)果與分析

3.1 樣地實(shí)測(cè)數(shù)據(jù)分析

77號(hào)小班的LiDAR數(shù)據(jù)經(jīng)數(shù)據(jù)處理后共獲得單木309株，與實(shí)測(cè)株數(shù)381株相比單木識(shí)別率達(dá)81%。在小班內(nèi)按20 m×20 m大小劃分小樣地34塊，以小樣地為單位，分別統(tǒng)計(jì)LiDAR估測(cè)數(shù)據(jù)(含樹高、冠幅、冠幅面積和冠幅體積)和實(shí)測(cè)胸徑的平均數(shù)、中位數(shù)、最大值、最小值以及標(biāo)準(zhǔn)差，結(jié)果如表1所示。

表1 34塊小樣地林分因子統(tǒng)計(jì)

在LiDAR估測(cè)數(shù)據(jù)中，平均樹高為18.7 m，其最大值為21.1 m，最小值為14.2 m，標(biāo)準(zhǔn)差達(dá)到1.5；平均冠幅為7.1 m，其最大值為12.8 m，最小值為3.9 m，標(biāo)準(zhǔn)差2.0；平均冠幅面積與體積標(biāo)準(zhǔn)差分別為24.6和128.3，數(shù)據(jù)離散程度較高。

在實(shí)測(cè)胸徑方面，平均胸徑為29.4 cm，最大值為33.8 cm，最小值為20 cm，標(biāo)準(zhǔn)差為2.6。

3.2 胸徑反演模型構(gòu)建

建模前需對(duì)34塊小樣地?cái)?shù)據(jù)進(jìn)行樣本異常值診斷和處理。經(jīng)分析后發(fā)現(xiàn)，第16號(hào)小樣地的LiDAR估測(cè)平均樹高較低，其原因是由于該樣地范圍內(nèi)林下分布植被，導(dǎo)致單木分割時(shí)錯(cuò)將林下植被作為國(guó)外松進(jìn)行分割，降低了樣地平均樹高值。因此需重新設(shè)置最小樹高，將林下植被進(jìn)行過(guò)濾，修正平均樹高后進(jìn)行胸徑反演。

基于LiDAR估測(cè)與樣地實(shí)測(cè)數(shù)據(jù)結(jié)果，按照基于單木分割的反演模型構(gòu)建中的方法，設(shè)置3組反演模型：

(1)實(shí)測(cè)胸徑與LiDAR樹高

(2)實(shí)測(cè)胸徑與LiDAR冠幅

(3)實(shí)測(cè)胸徑的自然對(duì)數(shù)與LiDAR樹高的自然對(duì)數(shù)

反演結(jié)果如圖3、表2所示。

表2 實(shí)測(cè)胸徑與估測(cè)胸徑回歸方程

在第(1)組模型中，以LiDAR估測(cè)樹高作為自變量，實(shí)測(cè)胸徑作為因變量建立反演模型。結(jié)果顯示：獲得方程y=26.7 lnx-48.7，擬合優(yōu)度R2=0.741 1，均方根誤差(RMSE)=0.004。

在第(2)組模型中，以LiDAR估測(cè)冠幅作為自變量，實(shí)測(cè)胸徑作為因變量建立反演模型。結(jié)果顯示：獲得方程y=- 0.233 9x2+ 3.142 3x+ 19.797，擬合優(yōu)度R2=0.603 1，均方根誤差(RMSE)=0.008。

在第(3)組模型中，首先將LiDAR估測(cè)樹高與實(shí)測(cè)胸徑分別進(jìn)行自然對(duì)數(shù)轉(zhuǎn)換，再以LiDAR估測(cè)樹高的自然對(duì)數(shù)作為自變量，實(shí)測(cè)胸徑的自然對(duì)數(shù)為因變量，建立反演模型。結(jié)果顯示：獲得方程y=- 2.184 9x2+13.569x-17.605，擬合優(yōu)度R2=0.806 3，均方根誤差(RMSE)=0.001。

圖3 3組胸徑反演模型

綜上所述，在3組反演模型中，分別將LiDAR估測(cè)樹高和實(shí)測(cè)胸徑取自然對(duì)數(shù)后，再建立反演模型的效果最佳，此時(shí)擬合優(yōu)度(R2)達(dá)到最高，均方根誤差(RMSE)最小。因此，結(jié)合本文樣地實(shí)測(cè)數(shù)據(jù)分析，當(dāng)LiDAR估測(cè)樹高范圍在14—22 m時(shí)，推薦使用第3組方程y=- 2.184 9x2+13.569x-17.605作為胸徑的最優(yōu)估測(cè)方程。

3.3 LiDAR估測(cè)蓄積量驗(yàn)證

250號(hào)小班面積為3.269 6 hm2，在其內(nèi)部選擇2個(gè)具有代表性的點(diǎn)位布設(shè)樣地，并每木檢尺，測(cè)得該小班平均胸徑為25.8 cm，平均樹高12.2 m，公頃蓄積量為183.6 m3/hm2。

250號(hào)小班經(jīng)LiDAR飛行后，獲得估測(cè)平均樹高為12.7 m；將該小班單木估測(cè)樹高代入胸徑反演模型y=- 2.184 9x2+13.569x-17.605進(jìn)行胸徑反算，得LiDAR估測(cè)平均胸徑為26.6 cm。查《江蘇省主要樹種一元材積表》中“馬1”列，根據(jù)LiDAR估測(cè)的單木胸徑統(tǒng)計(jì)對(duì)應(yīng)徑階的材積，得估測(cè)公頃蓄積量為157.2 m3/hm2。

根據(jù)《森林資源規(guī)劃設(shè)計(jì)調(diào)查技術(shù)規(guī)程(GB/T 26424-2010)》“13 調(diào)查精度”中允許誤差的要求，主要小班調(diào)查因子允許誤差分經(jīng)營(yíng)單位性質(zhì)、小班森林類別等分為A，B，C 3個(gè)等級(jí)。250號(hào)小班位于金壇區(qū)薛埠鎮(zhèn)東進(jìn)村金沙灣高爾夫球場(chǎng)內(nèi)，屬金壇區(qū)重點(diǎn)林區(qū)，因此小班調(diào)查因子允許誤差適用等級(jí)“A”。

經(jīng)LiDAR估測(cè)數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)對(duì)比(如表3)，平均胸徑誤差為3.1%，平均樹高誤差為4.1%，公頃蓄積量誤差為14.4%，均在《森林資源規(guī)劃設(shè)計(jì)調(diào)查技術(shù)規(guī)程(GB/T 26424-2010)》中允許誤差范圍內(nèi)(平均胸徑、平均樹高、每公頃蓄積量誤差需控制在5%，5%，15%范圍以內(nèi))，符合精度要求。

表3 LiDAR估測(cè)與實(shí)測(cè)數(shù)據(jù)對(duì)比

4 結(jié)論與討論

4.1 LiDAR估測(cè)受多因素制約

LiDAR技術(shù)較傳統(tǒng)的人工調(diào)查而言可節(jié)約大量人力物力，提升工作效率[16]。然而，LiDAR技術(shù)在實(shí)際應(yīng)用時(shí)易受多因素制約，影響精度。制約因素主要包括：LiDAR采樣密度、航帶匹配誤差和數(shù)據(jù)處理方法[17]。

在本研究中，選擇的250，77號(hào)小班均為國(guó)外松片林，周圍無(wú)建筑物遮擋，從單木結(jié)構(gòu)角度來(lái)說(shuō)確保了樹冠上部的回波信號(hào)，減少了信號(hào)的滅失，同時(shí)點(diǎn)密度為100個(gè)/m2，較好地反映了試驗(yàn)區(qū)現(xiàn)狀；在航帶匹配方面，對(duì)LiDAR飛行小班均進(jìn)行了航帶平差處理，保證了定位精度；在數(shù)據(jù)處理方面，由于試驗(yàn)區(qū)地貌均為丘陵，地形趨于平緩，因此在地面點(diǎn)分離方面采用改進(jìn)的漸進(jìn)加密三角網(wǎng)濾波算法(IPTD)，將點(diǎn)云分為地面點(diǎn)與植被點(diǎn)，然后分別生成DTM，DSM，CHM和NPC，對(duì)點(diǎn)云進(jìn)行單木分割。總體而言本次LiDAR估測(cè)盡量減少了試驗(yàn)誤差，估測(cè)精度較高。

4.2 LiDAR估測(cè)蓄積量可能性探討

一般而言，林分蓄積量與生物量隨平均胸徑的增加而增加[18]。研究表明，LiDAR估測(cè)樹高與實(shí)測(cè)樹高關(guān)系密切，而樹高也與胸徑、蓄積量有較高的相關(guān)性[19-20]。

本研究通過(guò)建立3組胸徑反演模型，結(jié)合擬合優(yōu)度(R2)挑選出優(yōu)選方程，反算了250號(hào)小班的林木胸徑。最后查《江蘇省主要樹種一元材積表》，得公頃蓄積量為157.2 m3/hm2。雖然其結(jié)果符合精度要求，但與實(shí)測(cè)公頃蓄積量183.6 m3/hm2仍有較大差距。分析原因，筆者認(rèn)為部分區(qū)域樹冠較為茂密，同時(shí)林下植被較為復(fù)雜，單木分割時(shí)對(duì)于樹冠邊緣難以區(qū)分，最終導(dǎo)致獲得的單木數(shù)量較實(shí)際偏少，估測(cè)公頃蓄積量較實(shí)測(cè)數(shù)據(jù)偏小。

4.3 LiDAR技術(shù)應(yīng)用綜合展望

本次采用LiDAR技術(shù)對(duì)金壇區(qū)國(guó)外松片林小班進(jìn)行了估測(cè)，并獲得了良好的效果，對(duì)于森林資源調(diào)查領(lǐng)域新技術(shù)應(yīng)用提供了技術(shù)參考。

目前對(duì)于機(jī)載激光雷達(dá)用于森林資源調(diào)查更偏向于理論研究，實(shí)際應(yīng)用較少。其主要原因在于LiDAR飛行時(shí)客觀因素制約較多，影響數(shù)據(jù)的誤差，而且飛行成本偏高，若進(jìn)行區(qū)縣范圍內(nèi)的全域飛行，將投入大量物力。此外，飛行效果視林分結(jié)構(gòu)而定，若林分結(jié)構(gòu)、樹種結(jié)構(gòu)、林下植被較為簡(jiǎn)單，則LiDAR估測(cè)精度較高，若林分結(jié)構(gòu)復(fù)雜區(qū)域，其精度則會(huì)降低[16]。

綜上所述，LiDAR技術(shù)是未來(lái)發(fā)展的大趨勢(shì)，前景應(yīng)用廣闊。在精度符合要求的情況下，可在樣地水平成果的基礎(chǔ)上大范圍應(yīng)用，節(jié)約人力物力，提高工作效率。