陳周娟,程 光,卜元坤,黃 維,陳佳卉,李衛(wèi)忠

(1.西北農(nóng)林科技大學(xué) 林學(xué)院,陜西 楊凌 712100；2.陜西省林業(yè)科學(xué)院,西安 710082)

單木結(jié)構(gòu)參數(shù)作為森林資源調(diào)查最主要的測樹因子,不僅能反映林木個體的生長情況,還能夠反映森林生態(tài)系統(tǒng)對地上部分資源的利用能力,對森林生物量的估測、森林碳庫的估算等方面起著重要作用。因此,高效獲取林木參數(shù)是掌握林木生長狀態(tài)及改進森林經(jīng)營管理工作的重要前提[1]。傳統(tǒng)的每木檢尺調(diào)查方法不僅耗時長、成本高,而且精度往往難以滿足現(xiàn)代精準(zhǔn)林業(yè)的需求。近年來,無人機傾斜攝影測量技術(shù)(Unmanned Aerial Vehicles Tilt Photogrammetric Technology)的發(fā)展與完善為高效獲取林木信息提供了新思路。與傳統(tǒng)森林資源調(diào)查方式相比,無人機傾斜攝影技術(shù)在采集數(shù)據(jù)方面具有經(jīng)濟、高效、方便、靈活等特點,其數(shù)據(jù)可以滿足人們對三維信息的需求[2-4]。無人機傾斜攝影技術(shù)在林木精準(zhǔn)測量、森林科學(xué)經(jīng)營方面具有廣闊的應(yīng)用前景,在森林資源調(diào)查上正逐步成為不可或缺的高新技術(shù)手段。

國內(nèi)外學(xué)者基于無人機傾斜影像數(shù)據(jù)的單木參數(shù)提取策略,通常將無人機影像數(shù)據(jù)經(jīng)過運動恢復(fù)結(jié)構(gòu)(Structure from Motion,SfM)獲取影像匹配點云[5-7],再將影像匹配點云轉(zhuǎn)化為柵格化模型——冠層高度模型(Canopy Height Model,CHM),單木參數(shù)提取的過程隨即轉(zhuǎn)變成二維影像的檢測與分割問題。識別單木的一般過程,包括樹頂探測和樹冠邊界勾繪過程,其中,局部最大值法[8-9]常用于樹頂點識別,區(qū)域增長算法[10]、分水嶺算法[11-12]、模板匹配法[13]等則常見于樹冠提取。Mochan等[14]以混交針葉林為實驗樣地,基于無人機衍生的CHM,利用局部最大值法進行單木識別,準(zhǔn)確率高達85%(F得分為0.86)；劉江俊等[15]研究證明了基于無人機影像利用局部最大值法可以較好地提取密集針闊混交林的樹頂點和樹高,樹頂點識別的F得分最高可達到0.77,樹高提取的均方根誤差為1.41m；袁梓健等[16]研究發(fā)現(xiàn),利用CHM最大值法分割華山松時,其精度比點云分割法高,林木株樹分割精度分別為87.17%和80.79%；陳崇成等[17]以2個分別種植羅漢松(Podocarpusmacrophyllus(Thunb.)D.Don)和桂花樹(Osmanthusfragrans(Thunb.)Lour.)的苗圃為樣地,運用種子區(qū)域增長算法分割單木樹冠,結(jié)果顯示,兩個樣地中樹冠估測的F得分均達到了0.89以上,冠幅誤差在0.14m以內(nèi)；杜意鴻等[18]研究表明分水嶺算法和均值漂移分割算法在提取單木油松冠幅上具有較大的潛力,冠幅估測的F得分均達到了0.8以上。以上研究表明,當(dāng)前基于無人機傾斜影像的林木參數(shù)提取技術(shù)主要是面向針葉林,少數(shù)基于無人機影像的闊葉樹相關(guān)研究也主要聚焦于樹種識別等方面。相較于針葉樹樹冠的錐形特征,闊葉樹常存在單個闊葉冠內(nèi)有多個局部峰值,樹冠輪廓更為復(fù)雜的問題,如何從無人機影像數(shù)據(jù)中準(zhǔn)確提取闊葉樹的頂點及樹冠輪廓是目前無人機影像數(shù)據(jù)提取技術(shù)需要攻克的重點和難點。

基于此,本研究以闊葉樹種銀杏(GinkgobilobaL.)為研究對象,基于影像匹配點云,采用局部最大值算法、種子區(qū)域增長算法和標(biāo)記控制分水嶺算法,對研究區(qū)內(nèi)單木進行定位識別及樹冠輪廓分割,旨在探索基于無人機傾斜影像快速精確地提取闊葉樹種單木參數(shù)的可行性與精度,為無人機傾斜攝影技術(shù)在森林資源調(diào)查中的推廣以及基于無人機影像的闊葉樹種單木參數(shù)快速提取提供理論和方法支撐。

1 研究區(qū)概況及數(shù)據(jù)來源

1.1 研究區(qū)概況

研究區(qū)位于陜西省咸陽市楊凌農(nóng)業(yè)高新技術(shù)產(chǎn)業(yè)示范區(qū)(34°14′～34°20′N,107°59～108°08′E),海拔在418.0～540.1m之間,屬東亞暖溫帶半濕潤半干旱氣候區(qū),具有春暖多風(fēng),夏熱多雨、秋熱涼爽多連陰雨、冬寒干燥等明顯的大陸性季風(fēng)氣候特征。年均氣溫12.9℃,無霜期211d。年均日照時數(shù)2 163.8h,年總太陽輻射量為114.86kcal/cm2,年均降水量635.1mm。

研究區(qū)內(nèi)主要的喬木樹種為銀杏。銀杏是銀杏科、銀杏屬落葉喬木,為陽性樹種,深根性,對氣候環(huán)境及土壤條件的適應(yīng)范圍廣,旺盛生長于冬春溫寒干燥或者溫涼濕潤、夏秋溫暖多雨,土壤肥沃、排水良好的地區(qū)。銀杏適應(yīng)能力強,是主要的造林樹種、防護樹種、抗病蟲樹種及耐污染樹種,在改善生態(tài)環(huán)境、固土保水、防治蟲害和保護生物多樣性等方面具有重要意義。

圖1 研究區(qū)位置示意圖

1.2 數(shù)據(jù)來源

1.2.1無人機數(shù)據(jù)采集

無人機數(shù)據(jù)采集于2020年9月11日中午進行,天氣晴朗,風(fēng)力小,適合無人機航拍作業(yè)。飛行利用大疆四旋翼無人機經(jīng)緯M300 RTK,搭載賽爾五目鏡頭PSDK 102S,為獲取地面分辨率達到1cm的無人機傾斜影像數(shù)據(jù),設(shè)置無人機飛行高度為64m,飛行速度為7.9m/s,旁向重疊度75%,航向重疊度80%,外擴邊距64m(表1),獲取5組研究區(qū)影像照片,共2 510張相片。

表1 PSDK 102S相機參數(shù)及M300RTK無人機飛行參數(shù)

1.2.2地面實測數(shù)據(jù)采集

地面實測數(shù)據(jù)采集于2020年9月25日開始。在研究區(qū)銀杏林分內(nèi),選取地形平坦(坡度<5°)且株數(shù)密度大于200株/hm2的典型樣地1塊,如圖1所示。樣地為梯形,樣地面積為8 160m2,樣地內(nèi)共有銀杏林木184株。使用載波相位差分技術(shù)(Real-Time Kinematic,RTK)記錄每株林木的GPS坐標(biāo)與樣地四角點坐標(biāo),以便在ArcGIS中準(zhǔn)確定位,定位精度在3～5cm之間。利用傳統(tǒng)方法對胸徑大于5cm的銀杏進行每木調(diào)查,調(diào)查包含胸徑、樹高、東西和南北方向的冠幅長度等因子(表2)。由于實地采集樹冠面積“真值”具有一定的難度,且精度難以保證,因此,本研究將采用研究區(qū)的數(shù)字正射影像(Digital Orthophoto Map,DOM)運用eCongnition軟件進行多尺度分割,再由ArcGIS軟件對分割結(jié)果進行手動修改和描繪,得到的冠幅面積作為實冠幅面積“真值”。

表2 實測單木參數(shù)統(tǒng)計

2 研究方法

2.1 數(shù)據(jù)處理與分析

2.1.1無人機影像預(yù)處理

為提取研究區(qū)三維結(jié)構(gòu)信息,需先對無人機影像進行運動恢復(fù)結(jié)構(gòu)(SfM)三維重建。SfM算法首先根據(jù)攝影相機參數(shù)及無人機拍攝時記錄的位姿數(shù)據(jù)(POS數(shù)據(jù)文件)對初始無人機影像數(shù)據(jù)進行特征點匹配,在稀疏重建處理之后得到稀疏結(jié)構(gòu)點云,經(jīng)空中三角測量加密進行多視角影像密集匹配得到稠密點云。點云數(shù)據(jù)經(jīng)插值生成數(shù)字表面模型(Digital Surface Model,DSM),經(jīng)地面點分類、歸一化等處理,生成數(shù)字高程模型(Digital Elevation Model,DEM),DSM和DEM之差,得到CHM；同時,由實景三維模型可獲取數(shù)字正射影像(Digital Orthophoto Map,DOM)(圖2)。本文采用Benltly公司的Context Capture軟件和R語言中的LidR包處理無人機影像數(shù)據(jù)。

圖2 無人機影像預(yù)處理流程圖

2.1.2單木識別和分割

1)基于局部最大值法的單木識別

CHM是表達樹冠上表面距離地面高度和樹冠水平分布的模型,樹冠邊緣與樹冠頂部有較大區(qū)別,一般而言,樹冠頂部亮度值較高。本研究使用基于CHM探測單木頂點位置的方法是局部最大值法(Local-maximum Based Algorithm)。局部最大值法可以通過設(shè)置一個移動窗口探測CHM的局部最大值為樹冠中心點位置。獲取準(zhǔn)確局部最大值的關(guān)鍵在于設(shè)置一個準(zhǔn)確合理的移動窗口范圍[19]。結(jié)合研究區(qū)的林分結(jié)構(gòu)與地面實測數(shù)據(jù),本文設(shè)置3種移動窗口大小(3m×3m,5m×5m,7m×7m)進行單木頂點探測。

2)基于種子區(qū)域增長算法的樹冠提取

種子區(qū)域生長算法(Seeded Region Growing Algorithm)的基本思想是將有相似性質(zhì)的像素點合并到一起,對每一個區(qū)域要先指定一個種子點作為生長的起點,然后將種子點周圍領(lǐng)域的像素點和種子點進行對比,將具有相似性質(zhì)的點合并起來繼續(xù)向外生長,直到?jīng)]有滿足條件的像素被包括進來為止。本研究應(yīng)用的種子區(qū)域增長算法,在LidR軟件包的dalponte2016函數(shù)[20]中實現(xiàn),用探測的單木頂點位置作為區(qū)域生長法的起始點,并根據(jù)單木實際生長特征,設(shè)定算法的判斷閾值,相關(guān)參數(shù)設(shè)置為th_seed=0.2,th_cr=0.75,其余參數(shù)值均為默認。

3)基于標(biāo)記控制分水嶺算法的樹冠提取

標(biāo)記分水嶺分割算法(Marked-controlled Watershed Algorithm)是根據(jù)分水嶺的結(jié)構(gòu)來進行圖像的分割,是一種模擬浸水過程的數(shù)學(xué)形態(tài)學(xué)分割方法。將CHM翻轉(zhuǎn)之后,假設(shè)模型局部中存在極小值點,在極小值點處開始浸水,隨著水的上升,圖形逐漸被浸沒,最終形成積水盆。過程結(jié)束后,若積水盆中的極小值點為單木頂點,則相鄰積水盆之間合并形成的分水嶺即為樹冠邊界。本研究使用的標(biāo)記控制分水嶺算法與種子區(qū)域算法相似,為避免傳統(tǒng)分水嶺分割造成的過分割現(xiàn)象,本研究首先在判斷單木頂點位置后,再根據(jù)林木實測因子參考值調(diào)試分水嶺算法中的參數(shù)閾值。標(biāo)記控制分水嶺算法使用LidR軟件包中的segment_trees函數(shù),相關(guān)參數(shù)閾值設(shè)置為th_tree=6,其余參數(shù)值均為默認。

2.2 精度驗證

2.2.1樹頂點識別驗證

樹冠頂點與單木根部平面位置通常存在偏差,因此,根據(jù)Hirschmugl等[21]提出的檢測方法,在實測單木根部位置建立1m緩沖區(qū),由預(yù)測頂點與緩沖區(qū)的相對位置判斷預(yù)測頂點位置的準(zhǔn)確性。3種常見單木樹頂檢測結(jié)果如圖3所示,在緩沖區(qū)內(nèi)若只有唯一提取值則作為正檢林木(圖3(a)-(b)),若有多株提取木則取距真實位置最近值作為正檢值,緩沖區(qū)外且樹冠范圍內(nèi)的作為誤檢林木(圖3(c)),若樹冠范圍內(nèi)沒有被檢測到林木頂點則記為漏檢林木(圖3(d))。

統(tǒng)計3種不同窗口下進行局部最大值法的單木探測結(jié)果,并進行精度分析與評價。本研究采用以下3個指標(biāo)進行單木識別衡量[22]:

注:紅點表示實測單木根部位置,藍點表示預(yù)測單木頂點位置,圓表示以實測單木根部位置為中心建立的1m緩沖區(qū)。

(1)

(2)

(3)

式中:r為探測率,指正確預(yù)測樹頂點的數(shù)量和研究區(qū)內(nèi)所有樹頂點數(shù)量的比率,衡量預(yù)測系統(tǒng)的查全率；p為探測準(zhǔn)確率,指正確預(yù)測樹頂點的數(shù)量與檢測出的所有樹頂點之間的比率,衡量模型的查準(zhǔn)率；一般情況下r和p指標(biāo)往往是相互制約的,F得分則為r與p的加權(quán)調(diào)和平均值,當(dāng)F得分越高時,則說明單木識別的方法比較有效；TP,FN,FP分別代表正檢、漏檢和誤檢的數(shù)量。

對預(yù)測樹高與實測樹高、預(yù)測冠幅與實測冠幅的擬合效果采用決定系數(shù)(R2)、均方根誤差(RMSE)2個指標(biāo)進行評價。其中,決定系數(shù)(R2)表示模型的擬合程度,均方根誤差(RMSE)表示模型的回歸效果。如下:

(4)

(5)

2.2.2樹冠分割驗證

判斷樹冠分割的準(zhǔn)確性,一方面需要考慮預(yù)測冠幅與實測冠幅之間的相對誤差,另一方面則需要判斷樹頂點是否與分割冠幅有較好的匹配。根據(jù)楊全月等[23]提出的方法,可使用單木頂點位置與冠幅面積綜合判定策略:樹木屬于正檢且冠幅面積小于判斷閾值時判定為正確率,則正確率的定義為:

(6)

式中:δ為預(yù)測冠幅與實測冠幅的相對誤差,δt為判定閾值,Nr為實測樹冠數(shù)量。預(yù)測冠幅與實測冠幅的相對誤差計算如下:

(7)

式中:S為算法提取的樹冠面積；Sr為對應(yīng)位置參考樹冠面積。

3 結(jié)果與分析

3.1 單木識別結(jié)果

利用局部最大值法進行樹頂點探測時,不同窗口大小的識別結(jié)果都存在一定的漏檢、誤檢現(xiàn)象(表3)。研究中發(fā)現(xiàn),采用5m×5m窗口大小進行識別效果最好,準(zhǔn)確率達到了80.11%,F得分為0.88；采用3m×3m窗口時,出現(xiàn)較多誤檢現(xiàn)象,而采用7m×7m窗口時出現(xiàn)較多漏檢現(xiàn)象(圖4)。

樹高預(yù)測由樹頂點所在位置的相對高程決定,圖5列出了研究區(qū)內(nèi)預(yù)測樹高與實測樹高之間的回歸方程,從中可以看出,研究區(qū)的預(yù)測樹高與實測樹高表現(xiàn)出很強的線性關(guān)系,線性方程的斜率接近于1,R2為0.99。無人機傾斜影像數(shù)據(jù)在重建三維模型時,可能會造成樹頂數(shù)據(jù)的缺失,以及受到點云密度的影響,預(yù)測樹高值總體上偏低。

表3 樹頂探測精度驗證結(jié)果

圖4 5m×5m窗口單木樹頂點識別結(jié)果

圖5 單木樹高提取結(jié)果精度評定

3.2 單木樹冠分割結(jié)果

兩種算法的樹冠分割結(jié)果如圖6所示。從圖6中可以看出,分水嶺算法在提取部分密集林木樹冠時表現(xiàn)出一定優(yōu)勢,對于樹冠邊緣復(fù)雜、細節(jié)多的影像輪廓勾繪較敏感,而種子區(qū)域生長算法對林隙大的銀杏冠層中勾繪更為精細。本研究對于樹冠分割的精度驗證方式是在相對誤差閾值梯度下,分別統(tǒng)計了兩種算法中樹冠分割結(jié)果正檢數(shù)量和正確率,并根據(jù)建立預(yù)測樹冠面積與實測樹冠面積參考值之間的回歸關(guān)系計算出了相應(yīng)的R2。兩種算法對于相對誤差閾值、正確率與R2的關(guān)系如圖7所示。由圖7可以看出,R2隨著相對誤差閾值的增大而減小,正確率則相反,這是因為相對誤差閾值增加時,樹冠投影面積預(yù)測誤差較大的樹被保留了下來,從而降低冠幅面積的提取精度,同時被判斷為正確探測的樹木數(shù)量也隨之增加。

圖6 單木樹冠分割結(jié)果

圖7 相對誤差閾值對提取結(jié)果的影響

當(dāng)相對誤差閾值取30%時,種子區(qū)域增長算法提取的樹冠正確率為73.14%,R2為0.98；標(biāo)記控制分水嶺算法提取的樹冠正確率為63.43%,R2為0.97。兩種算法中,預(yù)測樹冠面積與實測樹冠面積參考值建立的線性回歸關(guān)系如圖8所示。由圖可知,種子區(qū)域增長算法和標(biāo)記控制分水嶺算法的樹冠預(yù)測結(jié)果均與實測參考值顯著相關(guān),但種子區(qū)域算法的預(yù)測結(jié)果(RMSE=1.79m2)相對于標(biāo)記控制分水嶺算法預(yù)測結(jié)果(RMSE=2.44m2)誤差更小。本研究的單木樹冠提取結(jié)果,整體上較為準(zhǔn)確,達到了半自動化、快速精準(zhǔn)的樹冠分割的目標(biāo),能夠滿足林業(yè)精準(zhǔn)化、智能化的發(fā)展需求。

圖8 樹冠面積提取結(jié)果精度評定

4 結(jié)論與討論

4.1 結(jié)論

本研究以銀杏闊葉林為研究對象,基于無人機傾斜影像數(shù)據(jù),運用三維重建技術(shù)提取了單木三維點云數(shù)據(jù),通過不同探測窗口下的局部最大值法識別了單木頂點位置并提取了樹高,同時運用種子區(qū)域增長算法和標(biāo)記控制分水嶺算法實現(xiàn)了銀杏樹冠的提取,并以地面實測數(shù)據(jù)作為“真值”對單木參數(shù)提取結(jié)果進行了精度驗證。研究結(jié)果表明:采用5m×5m固定窗口進行單木識別的效果最好,銀杏的樹頂點識別的F得分達到了0.88,樹高預(yù)測值總體偏低；在樹冠預(yù)測結(jié)果相對誤差為30%的情況下,種子區(qū)域增長算法的樹冠提取精度高于標(biāo)記控制分水嶺算法,標(biāo)記控制分水嶺算法預(yù)測面積值相比之下較離散,說明在林分密度小的銀杏闊葉林中,種子區(qū)域生長算法比分水嶺算法能更好地處理樹冠之間的林隙。本研究證實了運用無人機傾斜攝影數(shù)據(jù)可以對闊葉樹進行單木識別與分割,達到了半自動化、快速精準(zhǔn)的樹冠分割的目標(biāo),可為無人機傾斜攝影技術(shù)在單木參數(shù)提取提供有效的理論與方法支撐。

4.2 討論

本研究中,當(dāng)移動窗口大小為3m×3m,5m×5m時,單木分割效果較好,7m×7m時漏檢數(shù)目較多,因為研究區(qū)中的銀杏闊葉林平均冠幅為5.8m,過大的窗口會使得多個樹冠被誤判為一個樹冠,而造成頂點的漏檢。同時,5m×5m窗口樹頂點檢測優(yōu)于3m×3m窗口,因為銀杏為闊葉樹種,不同于針葉樹種的頂端優(yōu)勢,觀測位置不同時可能存在多個樹冠頂點,3m×3m窗口下,漏檢有所降低但誤檢樹木明顯增多。本研究中的樣地為銀杏闊葉林,其樹冠結(jié)構(gòu)相比于針葉樹更為復(fù)雜,因此提取難度更大,這與Guerra-Hernández等[24]研究結(jié)果相似。因此,單木分割窗口選擇應(yīng)考慮樹種與單木實際冠幅。此外,通過可變窗口大小[2]、改變CHM像素大小[25]、基于冠層形態(tài)的局部最大值識別[26]等方式,可在一定程度上提高單木識別的準(zhǔn)確率,故單木識別算法需要不斷地改進,以適應(yīng)不同喬木類型的調(diào)查需求。

森林資源調(diào)查中很多重要的單木結(jié)構(gòu)參數(shù),例如,胸徑、樹高、冠幅、郁閉度等,都可以從CHM中直接或間接地獲得。但從CHM獲取的這些參數(shù)通常會比實際值低。在前人的研究中發(fā)現(xiàn),由于LiDAR數(shù)據(jù)中樹頂周圍的點云錯失或?qū)Ⅻc云數(shù)據(jù)轉(zhuǎn)換為CHM過程中存在的一些無效值,樹高通常會被低估1～3m[27-28],這與本文從無人機傾斜攝影測量技術(shù)獲取的CHM中估測樹高的結(jié)果相一致。同時,由于林木間的競爭關(guān)系,使得樹冠相互遮擋與覆蓋,冠幅往往也被低估[29-30]。本研究中的銀杏闊葉林林分密度中等,林木之間的相互遮擋相對較少,因此,樹冠的估測效果較好。

在以往大多數(shù)研究中,樹冠分割算法都是基于樹冠形狀為錐形這一設(shè)定來研發(fā)的,適合運用于針葉樹種。本研究針對銀杏闊葉林,運用標(biāo)記分水嶺算法與種子區(qū)域增長算法提取樹冠,取得了較好的預(yù)測效果。從總體來看,標(biāo)記分水嶺算法對與誤差閾值的變化更為敏感,種子區(qū)域增長算法在相對誤差閾值為50%正確率趨于穩(wěn)定,且R2的損失較小,這是由于種子區(qū)域增長算法能將具有相同特征的聯(lián)通區(qū)域分割出來,避免由于闊葉樹樹冠結(jié)構(gòu)復(fù)雜對分水嶺的干擾[31]。本研究中銀杏闊葉林分密度為中等水平,受陰影干擾較小,種子區(qū)域增長算法表現(xiàn)更佳。當(dāng)前,大多數(shù)的單木樹冠提取研究集中在比較簡單、規(guī)則或郁閉度不大的林分中,在未來的研究中,可進一步研究和探索復(fù)雜林分中林木參數(shù)提取的方法與精度,明晰林木快速自動檢測與分割機制,滿足未來更高效、更精確的森林資源調(diào)查需求。