基于SfM的針葉林無人機影像樹冠分割算法

楊全月 董澤宇 馬振宇 吳 悠 崔 琪 盧 昊

(1.北京農學院計算機與信息工程學院， 北京 102206； 2.北京林業大學信息學院， 北京 100083；3.弗萊堡大學環境與自然資源學院，弗萊堡 79106)

0 引言

樹冠是樹木獲取光能并進行能量轉換的主要場所，冠幅反映了樹木的生長活力及林木空間屬性[1]，準確獲取樹冠信息有助于監測樹木長勢、估算樹木生物量及小班蓄積量、預防樹木病蟲害[2]。利用傳統的光學遙感手段，可進行面向對象的圖像分割樹冠提取，一般使用光譜和紋理信息[3]、圖像梯度[1]或圖像分水嶺[4-5]等方法。隨著無人機(Unmanned aerial vehicle，UAV)技術的迅速發展，無人機航拍影像廣泛應用于高精度的森林資源調查等領域[6]。利用無人機影像可快速獲取林區高分辨率影像數據，運用傳統遙感圖像處理手段提取單木冠幅信息[7]。一般認為，激光雷達(Light detection and ranging，LiDAR)點云具有很高的幾何精度和細節反映能力，可精確提取單木樹高、冠幅、樹冠體積、林隙等信息[8-11]，基于機載LiDAR數據的單木分割算法中，具有代表性的分別為基于冠層高度模型(Canopy height model，CHM)[12-13]的分割算法和基于原始點云的分割算法[14]。

隨著攝影測量技術和計算機視覺技術的發展，基于雙目視覺原理，可對圖像構建立體像對、進行前方交會，求得林下樹木幾何結構信息[15-17]。由于計算性能的大幅度提高，傳統攝影測量的空三加密發展為密集匹配算法[18]，出現了快速有效的運動恢復結構(Structure from motion，SfM)技術。通過密集匹配點云提取林分參數已取得了初步進展，有效地解決了高分辨率圖像的過分割問題。曾健等[19]運用傾斜攝影測量技術提取了落葉樹人工林地形，發現無人機影像匹配點云與機載LiDAR點云具有很高的一致性。LI等[20]對有人機(Charge-coupled device，CCD)航片進行密集匹配，有效反演了東北地區地上生物量和葉面積指數。陳崇成等[21]利用密集匹配點云CHM模型分割了苗圃單木。劉見禮等[22]提出一種點云分層切片方法，實現了單木樹冠提取。許子乾等[23]結合LiDAR數字高程模型(Digital elevation model，DEM)研究亞熱帶森林林分特征變量，發現Lorey’s高等林分指標與UAV點云指標具有較高敏感性。

分析這些方法發現，無人機拍攝一般在人工林[19]、稀疏林地[22]、苗圃[21]等郁閉度較低的林分中才能直接觀測地表，得到DEM；在密林[23]、天然林[18, 20]中，影像密集匹配點云往往需要借助LiDAR等主動遙感手段獲取的DEM進行高度歸一化，再進行林木參數提取。這種依賴性顯然制約了無人機影像在大范圍天然林區自動提取森林參數的實用性。本文利用無人機快速機動航拍的優勢獲取地勢陡峭、坡度較大的生長季針葉林高重疊度影像，基于SfM技術進行三維表面重建，借鑒激光雷達點云分水嶺分割單木的思想，提出一種無需DEM數據支持的自適應鄰域分水嶺算法，旨在解決在茂密林區高分辨率無人機影像三維重建無法得到地形信息以及樹冠過分割的問題，從而提升無人機影像在森林資源調查中的實用性。

1 實驗區概況與數據獲取

1.1 實驗區概況

本研究區位于北京市西部百花山國家級自然保護區，地處北京市門頭溝區清水鎮境內，是以保護暖溫帶華北石質山地次生落葉闊葉林生態系統為主的自然保護區。森林植被垂直分異明顯，海拔為1 000～1 900 m。選擇區內一處北向坡面為實驗林地，海拔為1 100～1 200 m，面積約為8 000 m2，位置分布如圖1所示。林地坡度約為40°，地勢陡峭，生長植被主要為華北落葉松(Larixprincipis-rupprechtiiMayr)。

圖1 實驗區位置Fig.1 Overage of experiment region

1.2 無人機數據獲取

數據獲取時間為2018年9月，實驗區樹木處于生長季末期，部分林冠葉片顏色開始發黃，區分度較大。使用大疆精靈4型(DJI Phantom 4)無人機，搭載FC330型相機，內置Sony Exmor R CMOS傳感器，有效像素為1 200萬，感光單元尺寸為1.58 μm，焦距為3.61 mm，因此，在相對飛行高度100 m時地面采樣分辨率(Ground sample distance，GSD)約4.4 cm。由于目標區域坡度較大，坡面上下高差超50 m，為保證對不同高度處的林木均有近似分辨率成像，依據山勢設計了5種不同高度(1 235、1 244、1 265、1 276、1 318 m)航線進行重疊覆蓋。無人機作業采用自動巡航模式，使用DJI GS Pro型地面站進行控制。通過指定成像范圍和影像重疊度的方式自動計算飛行軌跡和曝光點，以定點懸停拍攝的方式完成區域成像。通過云臺控制鏡頭始終垂直向下以保證成像姿態角的相對穩定。無人機內置全球定位系統(Global positioning system，GPS)模塊記錄曝光點位置并以(Exchangeable image file format，EXIF)信息的形式存儲在圖像中。最終5個不同高度共獲得高分辨率無人機圖像509幅，實驗林地坡面區域內的影像曝光點如圖1所示。

2 原理與方法

本文方法總體思路是從二維影像構建三維模型，再從三維模型回歸林木二維結構信息，技術流程如圖2所示。

圖2 總體技術流程圖Fig.2 Overview technical flowchart

2.1 無人機影像SfM三維重建

2.1.1密集影像SfM處理

為提取林地三維表面結構，首先對原始無人機影像進行SfM三維重建。SfM算法步驟為：密集影像同名特征提取和匹配。利用同名特征點進行影像相對定向并計算外方位元素。空三處理生成稀疏點云及迭代光束法平差。基于稀疏點云進行加密得到密集點云。本研究使用Pix4Dmapper軟件處理無人機影像。特征點搜索使用結合特征匹配和像素匹配的Daisy算子，無人機內置GPS提供的曝光時刻機身坐標作為影像外方位元素的坐標初始值進行迭代優化，最后使用面片多視角立體視覺(Patch-based multi-view stereo，PMVS)算法搜索圖像中的像素以得到更多的匹配點，加密點云數據。SfM技術進行同名點搜索、匹配和點云加密的方式如圖3所示，圖3a下方背景為SfM生成彩色點云，空中藍色圓點為平差后像主點位置，藍色矩形表示相機投影方向，一組綠色射線表示地表一個同名點被多個影像同時觀測，前方交會于地物空間。圖3b表示該同名點在一系列不同影像中的搜索位置。初始相機內參數修正誤差為2.48%，聯合光束法平差的平均誤差為0.2像元，符合測繪要求。最終得到實驗區點云密度為600～1 000個/m2。同時，對原始無人機影像進行拼接處理和地理編碼，平均采樣分辨率為2.78 cm，利用點云高程信息進行了正射糾正并生成空間分辨率為0.05 m的彩色影像作為實驗區底圖(圖1)。

圖3 SfM技術密集匹配彩色點云示意圖Fig.3 Illustrations of densely matched true color point cloud from SfM

2.1.2三維表面模型處理

基于3D點云分割的單木提取方法依賴樹冠的點云包絡或較完整的點云冠層模型[14,22]，但由于研究區生長季樹葉茂密，林地郁閉度較高，無人機影像難以拍攝到林冠內部細節、林窗及林下地表，或無法在多張影像上構成同名點觀測。由圖3a可知，絕大部分點集中在樹冠頂層，導致上述3D點云分割方法適用性較差。但由于冠層表面三維點云密度極高，可利用表面模型處理方法進行樹冠提取。針對該特點，使用點云數字表面模型(Digital surface model，DSM)內插方法生成分辨率0.05 m的DSM模型，并進行了表面平滑處理以消除表面噪聲。

2.2 kNN自適應鄰域分水嶺分割

2.2.1單木檢測

單木檢測是為了找出樹木中心位置，作為后續樹冠分割的起始位置。為了充分利用無人機影像高分辨率紋理和SfM三維結構信息，設計了一種結合圖像信息和高程信息的檢測策略。

(1)基于圖像的單木檢測

首先，進行圖像預處理以增強樹冠邊緣和邊界差異。利用Sobel算子和Laplace算子增強圖像中樹冠的邊緣，使樹與樹之間、樹與地面之間的邊界更加清晰明確。利用中值濾波可以在一定程度上去除增強時出現的噪聲點。在預處理過程中，兩種算子和中值濾波交替使用。形態學開運算和閉運算可以使樹冠內部緊密、體積縮小，并且樹冠之間的間距增大，直方圖均衡化使圖像的樹冠部分和非樹冠部分的灰度差距增大。進行上述處理后，使用滑動窗口的大津閾值法將圖像分為多個小塊，對每一塊使用大津閾值法得到二值圖像。再對該二值圖像使用像素之間的歐氏距離變換計算初始樹冠中心[24]，最后使用圖像分水嶺法分割出樹冠。由于圖像進行了形態學開運算和閉運算，樹冠圖像不能準確代表真實樹冠，因此僅保留樹木位置信息。

(2)基于高程信息的單木檢測

由于樹木本身有向上生長的特點，在地形高差變化較大的區域，樹冠仍然符合樹頂高于其他部分的特點。基于高度信息進行樹冠檢測，一般利用樹冠幾何形態特征作為先驗知識，樹頂位置往往在高程上處于局部最大值點。參照文獻[22,25]，設定一個最小樹冠半徑R，對于表面模型柵格中的每一個點，判斷與距離小于R的所有點是否全部低于該點，將符合條件的像素視為初始樹頂點。分析發現，一般針葉林林分都較為均勻，定值R足夠進行檢測，自適應調整的影響不大。實際上，文獻中一般使用CHM進行檢測，但在坡度較大情況下，樹木冠幅內部不同像素對應的高程可能不一樣，導致CHM樹冠形狀發生扭曲，直接應用DSM可在很大程度避免該問題。

圖4 自適應kNN鄰域示意圖Fig.4 Illustrations of adaptive kNN neighborhood

(3)單木樹冠檢測結果合并

樹冠最高點一般是較為準確的樹頂點，而圖像識別的中心點不一定是準確的樹頂點，依據高程信息檢測的結果往往更加準確，但圖像信息可以在三維重建不精確的位置補充識別部分漏檢單木。由于三維模型的準確度有限且針葉樹木的特殊冠型(樹冠凹凸不平滑)，需設定較大的樹冠半徑R才能較好地避免過分割，但導致有些樹木不能被識別。對未被識別的樹，使用圖像信息作為補充，二者合并的具體策略是：對每個高程檢出點，去除其周圍半徑為R圓周內的圖像檢出點；再把剩余圖像檢出的樹頂點加入到高程檢出的樹頂點，形成完整的樹頂點檢測結果。

2.2.2自適應kNN鄰域檢測

(1)

其中dii=0(i=1,2,…,m)。待分割的目標樹為T，它的kNN鄰域定義為

(2)

(3)

2.2.3鄰域分水嶺分割

圖5 單木檢測結果Fig.5 Results of individual tree detection

分水嶺分割是一種常用于圖像分割的數學形態學方法，基本思想是將輸入的灰度圖像視為拓撲地貌，基于流溢模型進行分水嶺變換。對于圖像分割而言，通常在梯度圖像上進行分水嶺變換，而對于表示冠層高度的柵格模型，分水嶺可直接應用，不同之處在于樹冠高度模型經過“倒置”，一般由最高點開始“浸水”。由于已檢測了單木初始位置，本方法是一種標記控制的分水嶺算法，目的是基于kNN(T)和Ttree求得Ctree集合。Ttree為樹頂檢測得到的頂點集合，表示分水嶺算法的初始狀態。在每個kNN鄰域內，從k+1個單木初始點開始進行分水嶺處理，設kNN鄰域的高程范圍是[Hmin,Hmax]，隨著水位從Hmin到Hmax不斷增加，鄰域內的地形會被水漫過，每個被漫過像素將被分配給不同樹冠，當這k+1個樹冠兩兩之間出現連通時，在連通像素上標記阻斷。最終，鄰域內存在一些點屬于兩個或多個樹冠，這些點構成“分水線”，進而確定了樹冠像素集合，如圖4c所示。由于鄰域是圍繞目標樹構建的，只有包圍目標樹冠的分水線最接近自然邊界，故只取鄰域分水嶺結果中目標樹的分割區域作為該鄰域分割結果，其他鄰接樹冠予以舍棄，如圖4d所示。最后將所有單木作為目標樹的分割結果合并成完整的林分樹冠分割結果。綜上可知，鄰域分水嶺與傳統分水嶺最大區別在于鄰域內高程范圍的不同。傳統分水嶺作用于全局CHM模型，范圍一般為從0到最大樹高，而鄰域分水嶺由于不需高程歸一化，隨著不同目標樹鄰域在地形表面移動，自適應地獲得較小的高程范圍，避免高程范圍過大提取錯誤的分水線。

3 結果驗證與分析

本文使用于旭宅等[5]和陳崇成等[21]研究的驗證方式，利用地理信息系統(Geographic information system，GIS)在高分辨率無人機正射影像上進行目視解譯勾繪樹冠輪廓作為驗證。從圖1可看出，實驗區坡面西南側樹木冠幅相對較小，以2～4 m居多，東北側樹木冠幅相對較大，分布在4～6 m不等。以此為依據劃分了2塊驗證區域：區域A約為66 m×50 m，包括目視解譯的樹木246棵；區域B約為54 m×38 m，包括目視解譯的樹木112棵。區域A可加強表現分割算法的總體適用性，區域B樹木可分析局部效果。對GIS勾繪的冠幅多邊形計算幾何中心作為參考的單木位置，多邊形面積作為參考冠幅面積。

3.1 單木檢出精度

樣地樹頂檢測結果如圖5a所示，算法識別單木位置為圖中黑點，綜合了圖像檢出點和高程信息檢出點。單木檢測精度以檢測率表示，檢測結果有正檢、漏檢和過檢3種情況。正檢定義為在一個參考樹冠上，有且僅有一個檢測出的樹頂點(圖5b)；漏檢定義為一個參考樹冠上沒有檢測到樹頂點(圖5c)；過檢定義為在一個參考樹冠上檢測到兩個或多個樹頂點或不是樹冠的位置檢測到樹頂點(圖5d)。圖5b～5d中，黃點為目視解譯得到樹頂點，紅圈為算法識別樹頂點。正檢率D是評價結果的主要指標，用戶精度U突出過檢(錯分)，生產者精度P突出漏檢(漏分)，精度定義為

(4)

(5)

(6)

式中Nd——正檢數

Nc——樹冠內檢測出樹頂的參考樹冠數量

Nr——參考樹冠數量

Nall——算法檢測出的樹冠數量

統計2個區域各項檢測指標如表1所示。區域A檢出率為91.46%，區域B檢出率為92.86%，略高于區域A。可發現總體檢出率均較高，對于2～6 m冠幅的樹木而言，冠幅變化對檢出率的影響不大。但由于算法主要依賴高程信息進行樹頂檢測，樹冠越大，DSM冠型越趨于實際形態，因此檢測率略有提高。

表1 實驗區單木檢測結果Tab.1 Individual tree detection results of experiment regions

3.2 樹冠面積提取精度

檢出率反映了算法對初始單木位置的識別率，但檢出的樹頂點在樹冠分割過程中未必與參考冠幅有較好的匹配。為了進一步評價分割算法的總體性能，定義分割冠幅面積相對誤差為

(7)

式中S——算法提取的樹冠面積

Sr——對應位置參考樹冠面積

基于相對誤差，可使用位置和樹冠面積綜合判定的策略：在樹木中心位置符合正檢判定的前提下，冠幅面積相對誤差低于判定閾值(δt)時視為正確提取，故正確率定義為

(8)

兩個區域的單木樹冠提取結果如圖6所示，彩色多邊形分塊表示為本文算法分割樹冠結果。可看出本文方法對樹冠的提取無論在位置上還是形態上均基本符合實際樹冠情況。少數樹冠存在“直邊”分水線，發現是由于三維重建時DSM模型本身對茂密冠層沒有較深的透入，不能準確反映樹冠邊界輪廓所致，這一現象對面積提取影響不大。

圖6 提取樹冠結果Fig.6 Results of individual tree crown segmentation

圖7 相對誤差閾值對提取結果的影響Fig.7 Influence of relative error threshold on extraction results

對人工勾繪樹冠面積和算法提取面積進行線性回歸并計算決定系數R2，由于相對誤差閾值δt對正確率存在影響，計算不同δt下正確提取樹冠數量和面積R2的關系如圖7所示。由圖7可看出，R2隨δt增加而減小，因相對誤差閾值增加時，保留了更多面積誤差較大的樹木，從而降低了面積提取精度，同時樹木檢測數量(右側縱軸)隨之增加。

當δt取50%時，區域A正確率為76.02%，R2=0.834 9；區域B正確率為82.14%，R2=0.817 2。線性回歸如圖8所示。可以發現，主要為小樹的區域A樹冠面積精度略高于大樹居多的區域B，區域B內部算法提取樹冠面積總體偏大，但區域B正確率高于區域A。面積精度的差異主要來自DSM模型：樹冠之間存在連通性。對于常規DEM支持的分割算法，一般剔除一定高度閾值以下的CHM像素，以避免低矮植被和地物的影響，從而很大程度消除了樹冠之間的連通性。本研究由于不使用高度歸一化，故有少數樹冠之間的像素誤分為樹冠邊緣(圖4a～4c)。通過密集匹配點云提取單木，因難以得到冠層以下及地表觀測，SfM三維重建時僅在最頂層以拉伸填充的方式表示冠層，在茂密林區不可忽視。由于區域B樹冠較大，株數密度較低，導致樹冠邊緣像素引起的誤差也更明顯。對于正確率的差異，由于區域B大樹的檢出率較高，且面積分布更為均勻，相對誤差閾值變化引起的正確率變化較小。

圖8 樹冠面積回歸關系Fig.8 Regression relations of tree crown areas

林地樹冠參數提取結果如表2所示。由表2可知，本文結果滿足林業資源調查對樹冠提取的需求。

表2 單木樹冠提取結果與對比Tab.2 Results and comparisons of individual tree crown extraction

4 討論

4.1 算法比較

利用多種方法進行了樹冠分割提取的比較。實現了傳統的應用于全圖的分水嶺分割算法[12](以下稱為全局分水嶺)和SILVA等[13]提出的針對LiDAR數據CHM的分割算法(以下稱為SILVA2016)，正確率和面積R2比較如圖9所示。可見，不同算法在區域A的正確率差異不大(76.00%～79.00%)，區域B本文算法正確率為82.14%，略高于全局分水嶺(80.36%)，明顯高于SILVA2016(75.00%)。因不同算法使用的DSM模型相同，且均主要基于高程信息檢測初始單木，同時受到相對誤差剔除的影響較小。對于面積精度，在區域A本文算法(R2為0.834 9)明顯高于全局分水嶺(R2=0.685 8)和SILVA2016(R2=0.509 4)，區域B內本文算法與全局分水嶺接近(R2為0.817 2和0.799 0)，遠高于SILVA2016(R2=0.485 2)。由比較可知，本文算法避免了在未高程歸一化條件下的分割錯誤問題，在區域A體現最為明顯。區域A地形高差變化較大、樹冠為中等尺寸，全局分水嶺在沒有高程歸一化的條件下很難準確分割。

圖9 算法精度橫向比較Fig.9 Comparison of algorithm accuracy

圖10為全局分水嶺和SILVA2016方法的樹冠面積回歸情況，可以看出全局分水嶺同樣存在對大樹樹冠面積的高估(圖10c)。與圖8對比可發現，本文算法對面積提取的準確度更高，集中度更好。

圖10 算法的面積精度比較Fig.10 Area accuracy compare of algorithms

4.2 精度影響因素分析

本文單木樹冠分割算法結合了三維和二維數據處理流程，因此各環節處理效果對最終樹冠提取結果均有直接或間接影響。研究發現，SfM三維重建對林木參數提取影響最大，其精度直接決定樹冠形態是否真實。SfM重建受到幾個關鍵因素影響：①圖像本身拍攝質量，反映為圖像的清晰度、畸變和空間分辨率等特征。②圖像重疊度。③GPS初始坐標精度。較差的影像質量對同名點搜索匹配質量造成不可逆的降低，很難在算法處理環節恢復。因此，無人機作業應盡量在光照穩定、無風條件下進行，避免光照陰影和過度曝光，防止枝葉的搖晃造成同名點配準誤差。盡量提高影像重疊度以增加冗余觀測，可有效消除影像質量不佳對SfM總體結果的影響，同時可在不同航高疊加拍攝避免地形高差變化帶來航拍漏洞。在SfM處理中，利用區域網平差可有效調整曝光點坐標，因此無人機GPS精度對結果的影響相對可控。

同名點搜索與匹配直接影響點云坐標的精度，錯誤的冠層點云會使樹冠形態發生扭曲和形變，降低冠層的天然連續性，最有可能形成過檢(圖5d)。由于生長季樹葉茂密，無人機無法觀測到林冠下層，表面模型在樹冠邊緣和間隙普遍存在一定的模糊性，無法像激光雷達點云一樣準確地描述樹梢細節，會導致分水嶺分割算法出現“直邊”現象，也使分水嶺提取樹冠面積略高于真實面積。為了解決該問題，需提高無人機成像的清晰度和空間分辨率，避免在樹冠間隙進行錯誤的點云加密。本文算法選取kNN鄰元參數為4，可有效進行分割，對更復雜林分可嘗試增大k或進行二次分水嶺對邊緣進行優化。

復雜地形條件對無人機遙感單木識別具有較大的挑戰性，本研究區坡面地形條件復雜，坡度陡峭，代表了典型山地環境下的林分地形條件。一方面，復雜地形增大了無人機航線規劃難度，為了獲取完整坡面數據，使用了多個航高疊加的成像方式。另一方面，復雜地形條件下，地形坡度大、高程變化劇烈，無人機航拍獲取林下DEM存在很大的精度問題，采用DSM和DEM求差獲取CHM的方法受到了制約，逐像元的差值運算甚至可能扭曲原始冠型。本研究針對性地采用了DSM分割方法解決了這一問題，避免了地形變化對冠層提取的精度造成影響，可達到與傳統CHM分割方法相當的精度。

5 結束語

針對目前無人機技術在森林資源調查中的應用現狀和實際林區的特點，提出了一種基于SfM三維重建的kNN自適應鄰域分水嶺分割算法，充分利用密集航片的數據冗余特點進行三維表面重建，利用高程和紋理信息結合的方式實現了較高的樹木檢出率。自適應鄰域分水嶺分割方法無需無人機拍攝林下地表，無需第三方DEM數據進行高度歸一化，即可準確分割樹冠，同時，避免了高分辨率圖像的過分割問題，有效提升了無人機影像在森林資源調查中的實用性。