基于無人機影像的無瓣海桑單木提取與地上生物量估算

余楚瀅，龔 輝，曹晶晶，劉燕君，劉 凱

[ 1. 中山大學 地理科學與規劃學院//廣東省公共安全與災害工程技術研究中心//廣東省城市化與地理環境空間模擬重點實驗室，廣州510006；2. 南方海洋科學與工程廣東省實驗室（珠海），廣東 珠海 519000；3. 廣州市城市規劃勘測設計研究院，廣州 510060 ]

紅樹林是熱帶亞熱帶地區海岸潮間帶的木本植物群落，具有顯著的防風消浪、促淤造陸和凈化海洋等效果（陳玉軍 等，2012；盧昌義 等，2019），其在保護海岸帶生態安全中發揮著重要作用。紅樹林擁有強大的固碳能力，準確的紅樹林生物量估算是定量研究其生態系統生產力和生態系統碳循環的重要內容（田義超 等，2019；王子予 等，2020）。無瓣海桑（Sonneratia apetala）是海桑科海桑屬喬木，其具有良好的環境適應性，且生物量積累快、固碳速率高，生產力水平處于中國紅樹林群落中的較高位置，是人工紅樹林營造應用中重要的優質紅樹樹種之一（朱可峰 等，2011；周元滿 等，2012）。相比于傳統的紅樹林生物量估算方式，遙感技術已成為紅樹林生物量估算研究的重要手段和數據源。然而，基于衛星遙感影像估算精度通常受限于較低的影像分辨率，且大多數的紅樹林生物量估算研究對所有樹種進行整體建模，較少考慮不同樹種的生物量差異。由于紅樹林群落冠層密集、結構復雜，目前，樹種尺度的紅樹林生物量估算研究仍較少。

單木尺度的森林結構信息有助于提高紅樹林生物量估算精度。基于遙感技術的單木提取包括單株立木探測和單木樹冠提取兩部分（董新宇 等，2019）。單株立木探測指對單木樹冠的所在位置進行探測，找出每個樹冠的中心點；單木樹冠提取指將探測到的樹冠中心點作為參照，找到樹木樹冠的邊界點，對樹冠輪廓進行描繪（劉曉雙 等，2010）。目前，國內外單木提取研究主要集中在針葉林（Dalponte et al., 2014; Yang et al., 2016），這是因為針葉林的樹冠呈尖塔形，樹冠高點明顯，而紅樹林由于樹冠比較密集，相鄰樹木之間高度差較小，在遙感數據中難以識別樹冠邊界，紅樹林的單木提取精度受限。近年來，廣受關注的低空無人機（Unmanned aerial vehicle, UAV）遙感技術為從高分辨率影像中提取單木樹冠信息提供新的機會，且已在紅樹林單木提取中得到應用。如Yin 等（2019）基于無人機載激光雷達（Light Detection And Ranging, LiDAR）數據實現了紅樹林的樹冠分割和參數估計，并分析了樹冠聚集度和點云密度對單木提取精度的影響；Navarro 等（2020）使用無人機影像生成高分辨率的冠層高度模型（Canopy Height Model, CHM），基于可變窗口濾波（Variable Window Filter, VWF）和控制分水嶺分割方法（Controlled Watershed Segmentation Method, CWSM）實現紅樹林單木提取。

以往研究主要利用地面樣方數據，采用皆伐法、平均木法和異速生長法等方法（薛立 等，2004）構建紅樹林地上生物量模型。其中，異速生長方程是一種通過部分單木參數得到單木生物量的方法，相比于皆伐法、平均木法，該法較為簡便，破壞性小，便于推廣（Castedo-Dorado et al.,2012； 金川 等，2012）。基于遙感技術的紅樹林地上生物量估算方法可分為兩類：一是提取光譜、指數、紋理等遙感特征，建立生物量反演模型（Lei et al., 2022）；二是通過高分辨率影像、LiDAR數據等獲取樹高、冠幅等林分因子，基于異速生長方程構建生物量估算模型（Pham et al., 2019）。前者結合多個特征變量，可達到較高的精度，但反演模型可推廣性較差。后者基于異速生長法理論，具有較強解釋能力。但衛星遙感影像通常受限于較低的影像分辨率（Zhu et al., 2017），難以實現單木尺度的紅樹林生物量估算。近年來，無人機遙感在精細尺度的植被生物量監測中備受關注，部分學者也開始探討其在紅樹林生物量估算中的應用潛力。如Qiu 等（2019）結合WorldView-2 高分辨率影像和無人機載LiDAR 數據進行單木和格網尺度的紅樹林地上生物量估算；Navarro 等（2020）利用無人機影像進行紅樹林單木檢測，并結合異速生長方程估算生物量。然而，無人機遙感在紅樹林生物量監測中的應用潛力仍有待進一步探究，且當前基于無人機遙感的單木尺度無瓣海桑生物量估算研究仍較少。

因此，本文以廣東省珠海市淇澳島紅樹林自然保護區為研究區，基于無人機影像和種子區域生長（Seed Region Growing, SRG）算法進行無瓣海桑單木提取，通過確定研究區無瓣海桑樹高和胸徑之間的回歸關系，構建基于樹高的無瓣海桑異速生長方程以準確估算其地上生物量。以期為紅樹林生態恢復、保護與管理提供決策支持。

1 研究區與數據

1.1 研究區介紹

廣東省珠海市淇澳島（22°23′40″-22°27′38″N，113°36′40″-113°23′40″ E）位于珠江口內伶仃洋西側，東望深圳、香港，北鄰中山、東莞，海島總面積約24 km2（王樹功 等，2005；Cao et al.,2021）。淇澳島屬熱帶海洋性季風氣候，雨量充沛、光照充足、雨熱同期（黃建榮 等，2011）。淇澳島紅樹林保護區內有秋茄（Kandelia candel）、桐花樹（Aegiceras corniculatum）、銀葉樹（Heritiera littoralis Dryand）、老鼠簕（Acanthus ilicifolius）和無瓣海桑（Sonneratia apetala）等主要紅樹植物。本研究區位于淇澳島紅樹林保護區的中、高潮區域（圖1），面積為2.89 hm2。研究區的植被群落中上層木主要為無瓣海桑，中層木主要為桐花樹和老鼠簕，其中無瓣海桑為研究區的優勢種。研究區的無瓣海桑為2001-2002年左右種植（Yu et al., 2020），2019年采集本文使用的無人機影像時無瓣海桑的樹齡約為17～18 a。

圖 1 廣東珠海淇澳島地理位置（a）與研究區域（b、c）及野外調查樣方點位置（c）Fig.1 Geographical location of Qi'ao Island, Zhuhai, Guangdong (a);Location of the study area (b, c); Quadrat positions in the study area (c)

1.2 無人機影像獲取與預處理

采用大疆DJI精靈4 PRO V2.0四旋翼消費級可見光無人機設備進行研究區數據采集。該無人機質量輕、體積小，搭載CMOS 傳感器，鏡頭焦距為8.8 mm，無人機平臺與RGB 傳感器集成一體，整機重量1 375 g。共有紅、綠、藍3個波段，傳感器分辨率為5 472×3 648（大疆創新，2021）。無人機影像的采集時間為2019-11-10，天氣狀況良好，飛行時段選擇正午T 10:00-14:00，風速約為1～4 m/s。使用Altizure生成無人機傾斜攝影航線，航向重疊率設置為80%，旁向重疊率設置為75%，航高設置為100 m，飛行速度約為9 m/s，傾斜航線的俯仰角度為45°。共獲取覆蓋研究區的無人機傾斜攝影影像566幅。

基于傾斜攝影建模軟件ContextCapture 生成研究區無人機正射影像圖（Digital Orthophoto Map, DOM） 和數字表面模型（Digital Surface Model, DSM）以及密度點云（圖2）。其中DOM空間分辨率約為0.04 m，DSM 空間分辨率為0.25 m，密度點云面積為0.343 km2，空間分辨率為0.03 m。

圖2 基于無人機航片生成的研究區正射影像圖DOM（a）、數字表面模型DSM（b）和密度點云（c）Fig.2 Digital Orthographic Model DOM (a), Digital Surface Model DSM (b) and point cloud (c) of the study area based on UAV Imagery

1.3 地面數據采集與預處理

1.3.1 實測單株無瓣海桑樹高和胸徑數據 使用Trimble SX10三維激光掃描儀、測高儀和卷尺對無瓣海桑進行調查，獲取無瓣海桑的樹高，并測量樹干1.3 m處的直徑（胸徑），以構建無瓣海桑樹高與胸徑的回歸關系。由于紅樹林生長在潮間帶且生境復雜，研究區內難以架設激光掃描儀設備，因而，采樣地點主要選取研究區邊緣的無瓣海桑群落。采集時間為2018-06-28，共采集了48棵無瓣海桑的樹高和胸徑數據，用于構建兩者的回歸關系。

1.3.2 樣方數據 針對研究區開展無瓣海桑地面樣方調查，數據采集時間為2020-11-03，在研究區內構建5 個10 m×10 m 的樣方（見圖1）。借助測高儀和近距離目視估測無瓣海桑的樹高；使用軟尺測量樹干1.3 m處的直徑（胸徑）。獲取的樹高和胸徑數據主要用于代入異速生長方程以獲取樣方地上生物量，驗證無人機獲取的結果。

1.4 技術路線

基于無人機傾斜攝影影像的無瓣海桑單木提取與地上生物量估算技術路線如圖3所示。首先，基于無人機影像生成冠層高度模型；然后，采用局部極大值算法和種子區域生長算法進行無瓣海桑單木提取；同時，利用地面調查數據確定樹高和胸徑之間的關系，構建基于樹高的無瓣海桑異速生長方程；最后，估算研究區的無瓣海桑地上生物量，并結合地面調查數據進行驗證分析。

圖3 無瓣海桑單木樹冠提取與地上生物量估算流程Fig.3 Flow chart of tree crown delineation and aboveground biomass estimation of Sonneratia apetala

2 基于無人機點云的無瓣海桑單木樹冠提取

2.1 點云濾波與冠層高度模型生成

基于無人機航片生成的研究區原始點云數據可分為植被點云和地面點云兩部分。為了精確獲取植被的冠層高度模型（Canopy Height Model, CHM），通常需先將地面點與非地面點分離開。已有研究證明，布料模擬濾波（Cloth Simulation Filtering,CSF）算法效果較好且所需參數較少（Zhang et al.,2016； 張昌賽 等，2018）。在CloudCompare 中采用CSF算法進行點云濾波，進而分離地面點與植被點。最終得到的點云濾波結果如圖4-a 所示，可以看出，地面點云主要包括林分外圍灘涂點云和林窗間隙內的地面點云。對CSF濾波處理得到的地面點云數據進行處理，采用最鄰近算法進行插值，并在ArcGIS中使用低通濾波進行平滑處理，生成數字高程模型（Digital Elevation Model, DEM），將DSM和DEM 進行差值運算得到研究區冠層高度模型（Canopy Height Model, CHM），重采樣后空間分辨率為0.25 m（圖4-b）。運用CHM 模型描述研究區林木從地表至冠層的高度分布。

圖4 研究區地面點云（a）和CHM（b）Fig.4 Ground point cloud (a) and CHM (b) of the study area

2.2 單木位置探測與樹冠分割

單木提取通常包括2個步驟：首先，進行單木位置探測，即將搜索到的樹冠頂點作為單木的空間位置；再以提取的樹冠頂點為標記點，進行單木的樹冠分割。采用經典的局部極大值算法（Yin et al.,2019）進行單木頂點提取。該算法將樹冠最高點作為單木位置，通過移動窗口逐步對柵格數據進行搜索，并判斷搜索窗口的中心點是否為局部極大值，若為局部極大值，則將此像素標記為單木頂點。單木探測結果的準確性通常與搜索窗口大小、影像分辨率和冠幅大小等因素有關。采用可變窗口濾波（Variable Window Filter, VWF） 算法（Popescu et al., 2004）確定搜索窗口大小。VWF算法基于樹高與冠幅的相關關系假設，即樹木越高則冠幅越大。該算法根據窗口中心的像素值確定移動窗口的大小，其中可變窗口函數的選擇非常關鍵。根據前期地面調查數據和多次試驗，進而確定函數y=0.06x為窗口變換函數（其中y為窗口大小，單位為m；x為樹高，單位為m）。以上步驟通過R 語言Forest Tools工具包實現，最終得到研究區紅樹林單木頂點的局部探測結果（圖5）。

圖5 無瓣海桑單木頂點探測結果局部（a. 局部位置示意圖；b. 局部區域正射影像；c. 局部區域單木頂點檢測結果）Fig.5 Partial result of individual tree detection of Sonneratia apetala (a. Location of the partial area; b. DOM of the partial area;c. Individual tree detection result of the partial area)

結合研究區CHM 影像，在SAGA GIS 中采用種子區域生長（Seed Region Growing,SRG）算法（Yin et al., 2019）進行單木樹冠分割。SRG 算法將圖像種子點標記為生長區域，進一步判斷相鄰像素與種子點的相似度，若相似度在閾值以內則并入生長區域。不斷循環上述操作直到滿足停止條件或遍歷結束。該算法利用特征閾值和距離閾值控制分割對象的生長。使用樹冠頂點作為種子點進行區域生長，以進行樹冠分割，得到的研究區域無瓣海桑群落分割結果及局部分割效果（圖6）。

圖6 無瓣海桑單木分割結果（a. 研究區總體單木分割結果與局部區域位置示意圖；b. 局部區域正射影像；c. 局部區域單木分割結果）Fig.6 The individual tree crown delineation result of Sonneratia apetala (a. Individual tree crown delineation result and the location of the partial area; b. DOM of the partial area; c. Individual tree crown delineation of the partial area)

2.3 單木提取精度

利用研究區的DOM影像、CHM影像，采用目視解譯，同時結合前期實地調研，人工隨機選取100 個無瓣海桑樣點，勾繪樣點的樹冠邊界，進而生成單木提取結果的驗證樣本，驗證效果如圖7所示。考慮到樹冠分割結果與樣本之間存在多種空間關系，參考已有研究，將分割結果與驗證樣本之間存在重合且重合部分同時占分割結果與驗證樣本樹冠50%以上（Zhen et al., 2013）定義為正確分割。采用探測精度（Detection Accuracy, DA）指標進行樹冠分割精度評價（Yin et al., 2019），公式為：

圖7 研究區局部無瓣海桑單木樹冠提取結果與驗證樣本（a. 局部位置示意圖；b. 局部區域正射影像；c. 局部區域單木分割結果與驗證樣本）Fig.7 Individual tree crown delineation result and reference samples of Sonneratia apetala (a. Location of the partial area;b. DOM of the partial area; c. Individual tree crown delineation result and reference crowns of the partial area)

式中：T為正確分割的樹冠數量（株）；R為所有樣本的數量；DA 為正確分割的樹冠數量T與所有樣本數量R的比例。

由于研究區內無瓣海桑冠層十分密集、群落內部冠層相互交疊、樹冠間縫隙較小，且樹冠之間的高度變化較為平滑，使得單木尺度的無瓣海桑提取工作具有挑戰性。使用SRG算法得到的無瓣海桑單木樹冠探測精度DA 為67%，基本能較為準確地將研究區的無瓣海桑單木樹冠分割出來。

3 單木尺度的無瓣海桑地上生物量估算

3.1 無瓣海桑生物量估算模型構建

異速生長方程是林木生物量估算的經典方法，其通過建立樹高、冠幅、胸徑等林木參數與生物量之間的經驗方程估算林木生物量（武高潔 等，2014）。已有研究表明，異速生長方程具有較高的生物量估算精度，且相比傳統的皆伐法和平均木法，其能夠大幅降低對森林的破壞性采樣，適用于自然保護區等特殊政策環境下的受保護樹種（Navarro et al., 2020）。同種樹種的異速生長關系往往有較高的相似性，區域差異較小。因此，對不同地域的同一樹種可以使用同一異速生長方程（Komiyama et al., 2008）。

參考Ren 等（2010）利用實地調查數據建立的無瓣海桑異速生長方程（表1）構建無瓣海桑單木地上質量估算模型。Ren 等的實驗區域為廣東省雷州灣紅樹林，該區域與廣東省珠海市淇澳島處于相近的緯度，且具有相似的紅樹林生境與樹種組成。同時，Zhu等（2021）也將該異速生長方程應用于深圳福田紅樹林保護區17～18年樹齡的無瓣海桑單木地上質量估算。基于表1中地上部分的干、枝、皮和葉生物量求和，得到無瓣海桑單木地上質量。

表1 基于胸徑和樹高建立無瓣海桑異速生長方程Table 1 The allometric growth equation of Sonneratia apetala based on DBH and tree height

已有研究表明，無瓣海桑的樹高和胸徑之間有較強的相關性（朱可峰 等，2011），因此，通過建立樹高和胸徑的回歸方程，利用樹高推算胸徑，進而代入異速生長方程以估算無瓣海桑單木地上質量。利用地面采集工作獲取的48株無瓣海桑單木樹高和胸徑，并采用最小二乘法進行回歸擬合，得到研究區無瓣海桑的樹高和胸徑的回歸方程，其公式為：

式中：D為胸徑，單位為cm；H為樹高，單位為m。

基于地面調查數據構建的無瓣海桑的樹高和胸徑線性回歸趨勢線（圖8），可以看出，樹高和胸徑具有較高的相關性，相關系數R2為0.871 3。

圖8 基于地面調查數據得到的無瓣海桑樹高-胸徑回歸線Fig.8 Regression line of height-DBH of Sonneratia apetala based on ground survey data

利用上述樹高-胸徑回歸方程，將異速生長方程（見表1）的胸徑參數用樹高參數表示，從而得到優化后的異速生長方程（表2）。

結合無瓣海桑單木提取結果，搜索每個樹冠多邊形內的最大值，將其作為無瓣海桑的單木樹高，代入優化后的異速生長方程（見表2），得到無瓣海桑單木地上質量。進而計算無瓣海桑地上生物量，公式為：

表2 優化后的無瓣海桑異速生長方程Table 2 The optimized allometric growth equation of Sonneratia apetala

式中：AGB是單位面積內無瓣海桑地上生物量，單位為t/hm2；W為單位面積內無瓣海桑單木地上質量總和，單位是kg；A為單位面積，單位為hm2。

3.2 地上生物量估算結果與精度評價

從研究區無瓣海桑單木地上質量分布（圖9）可以看出，無瓣海桑單木地上質量的范圍為29.60～388.44 kg，平均值為145.72 kg，地上質量總和為368.97 t。

總體上，研究區內無瓣海桑單木地上質量分布存在部分空間集聚現象。采用全局Moran'sI指數（Moran, 1950）進一步評估研究區內無瓣海桑單木地上質量分布的空間集聚程度，得到其Moran'sI指數值為0.594（顯著性檢驗結果p值為0，Z檢驗結果為50.934），表現為較強的空間聚集性。可以看出，林分邊緣和林窗部分的無瓣海桑單木地上質量較小，如圖10所示，主要可能有2個原因，一是部分研究區內的無瓣海桑正常死亡后留下林窗，其周圍的幼苗仍處于生物量累積階段，二是可能是由于林分邊緣的無瓣海桑易受自然和人為因素的干擾，其單木地上質量常常受到損失。

圖10 林分內部（a）與林分邊緣和林窗處（b）的無瓣海桑單木地上質量對比Fig.10 Comparison of individual tree mass of Sonneratia apetala inside the forest (a) with that at the edge of the forest and at the gap (b)

通過對比基于無人機影像估算的無瓣海桑單木地上生物量、他人估算的研究區地上生物量與實地樣方調查獲得的地上生物量，進而驗證利用無人機影像估算地上生物量的適用性和準確性。參考地面調查樣方大小，將無人機估算得到的研究區無瓣海桑地上生物量分布分割為399個10 m×10 m的格網。圖9-b 為按研究區10 m×10 m 樣方尺度統計的無瓣海桑地上生物量分布，計算得到每個樣方內無瓣海桑平均單木數為6.3棵，樣方內最多單木數為17棵，最少為1棵；得到的樣方平均地上生物量為92.14 t/hm2，最小值為2.99 t/hm2，最大值為247.24 t/hm2。根據前期實地樣方數據顯示，研究區內5 個10 m×10 m 無瓣海桑樣方的平均單木數為5.3 棵，平均地上生物量為130.89 t/hm2，最小值為52.69 t/hm2，最大值為269.09 t/hm2。可以看出，基于無人機影像估算得到的研究區無瓣海桑生物量總體上低于地面實測生物量數據。Zhu 等（2020）利用WorldView-2高分辨率影像和無人機DSM 數據，采用隨機森林算法對廣東珠海淇澳島紅樹林保護區地上生物量進行估算。首先將淇澳島內植被主要分為秋茄、無瓣海桑與其他植物3種類型，將遙感影像的波段特征和植被指數作為輸入變量，實地考察并通過異速生長方程得到的地上生物量作為輸出變量對模型進行訓練，進而得到整個保護區的地上生物量，RMSEr為30.48%。其得到的2016 年研究區平均地上生物量為92.53 t/hm2，最小值為39.24 t/hm2，最大值為249.99 t/hm2，總體上與本文估算的地上生物量結果一致。

圖9 研究區無瓣海桑單木地上質量分布（a）和樣方尺度地上生物量分布（b）Fig.9 Aboveground mass distribution (a) and quadrat scale biomass distribution (b) of Sonneratia apetala in the study area

4 結論

針對廣東珠海淇澳島紅樹林保護區的無瓣海桑群落，基于無人機影像和種子區域生長算法實現單木樹冠提取，構建了無瓣海桑樹高-胸徑回歸模型以優化生物量異速生長方程，進而實現了研究區無瓣海桑地上生物量的精確估算。得到的主要結論為：

1）基于無人機影像能夠實現無瓣海桑的單木分割，且單木樹冠探測精度為67%，并驗證其在單木尺度無瓣海桑地上生物量估算中的可用性，可為后續基于無人機遙感的紅樹林生物量監測提供技術參考。

2）利用地面生物量調查數據分析了研究區無瓣海桑的樹高和胸徑間較高的相關性，相關系數R2為0.871 3，進而構建了基于樹高估算無瓣海桑地上生物量異速生長方程。

3）基于無人機遙感估算得到研究區無瓣海桑的平均地上生物量為92.14 t/hm2，最低為2.99 t/hm2，最高為247.24 t/hm2，整體上略低于實地調查的地上生物量。

此外，本研究存在以下幾點不足：1）由于無人機未攜帶RTK 系統，可能會造成垂直定高上約1 m的誤差，同時無瓣海桑的樹種特點為同一植株樹冠可能存在多個樹頂點，這些因素會造成單木分割的誤差從而影響地上生物量計算；2）受到局地氣候條件與潮汐等因素的影響，研究區內有部分無瓣海桑處于倒伏狀態，樹冠被其他未倒伏無瓣海桑樹冠遮蓋，而這些情況在高空視角往往難以獲取，該部分無瓣海桑地上生物量在基于遙感的研究中通常不被考慮，但在研究區林下實地調查時卻可以被統計，因而，基于無人機影像與高分辨率衛星影像估算的無瓣海桑地上生物量相比地面實測數據總體上偏低；3）由于研究區的無瓣海桑屬于人工種植，相比灘涂區自然生長的無瓣海桑群落密度較低，因而本文單木分割算法中閾值的適用性在群落密度較高的區域仍有待探究。