無人機影像匹配點云單木識別算法

劉見禮，張志玉，倪文儉，汪垚，張大鳳

(1.中國科學院遙感與數字地球研究所 遙感科學國家重點實驗室，北京 100101；2. 中國科學院大學，北京 100049)

0 引言

森林是陸地生態系統的主體，也是人類賴以生存的最重要的自然資源之一[1-3]。森林的單木信息在森林的生態研究、資源監測等方面具有重要作用[4-5]。

單木信息提取的前提是單木的準確識別，目前單木識別主要利用機載或星載的高分辨率光學影像數據的光譜與空間特征信息，以及LiDAR點云數據的三維幾何信息[6-7]。由于直觀的三維幾何信息更易于單木識別，因而成為目前的研究熱點[8-9]。

常用的三維幾何信息通常來自LiDAR點云數據或其衍生的數字表面模型(digital surface model，DSM)、冠層高度模型(canopy hight model，CHM)數據[10-11]。單木識別的方法主要為“局部最大值”法[12-13]。由于LiDAR數據獲取成本較高，限制了其推廣應用，而無人機影像匹配點云數據的獲取成本較低，在單木識別中具有較高的應用潛力[14]。

目前缺少針對無人機影像匹配點云數據的單木識別算法，因此我們借鑒基于LiDAR數據的“局部最大值”算法進行單木識別研究?！熬植孔畲笾怠彼惴ǜ鶕阉鞔翱诘牟煌譃閮煞N：一種是“固定窗口”的局部最大值；另一種是“變化窗口”的局部最大值。變化窗口的局部最大值通常需要精確的林下地形信息，而無人機影像點云在郁閉林區難以獲得精確的林下地形信息，因而難以采用變化窗口的局部最大值法進行單木識別[15]。

因此使用固定窗口大小的局部最大值法進行無人機影像點云數據的單木識別。已有的研究表明，固定窗口的局部最大值單木識別中，當窗口設置較大時，識別出的單木準確性較高，但會發生遺漏；窗口設置較小時，識別單木的數量會增多，但易引入了錯識別單木[16]。目前基于局部最大值的單木識別缺乏針對錯分單木的判別算法。發展有效的單木真偽判別算法，解決“小窗口”局部最大值算法產生的錯識別單木問題，是提高單木的識別精度的關鍵。

為此，本文提出了一種聯合“局部最大值”與“樹冠結構分析”的單木識別算法，該算法可解決 “局部最大值”算法的窗口大小設置的困境。

主要思路：(1)基于無人機匹配點云數據，通過小窗口的“局部最大值”，獲得盡可能多的候選單木；(2)對每個候選單木樹冠進行連續水平切片分析，獲得單木樹冠結構曲線，通過單木樹冠結構曲線對候選單木進行真偽判別，濾去偽單木；(3)根據樹冠切片中心對單木初始位置進行改正，對位置改正后相互聚攏的單木進行合并，最終實現單木準確識別。

1 研究區與數據

1.1 研究區

研究區位于大興安嶺北段西坡，內蒙古呼倫貝爾市東北部，屬于根河林業局。研究區氣候為寒溫帶濕潤型森林氣候，全年平均氣溫-5 ℃左右，年平均降水量294 mm。區域內優勢樹種為興安落葉松(larix gmelini)，伴生樹種有樟子松(pinus sylvestris var)、白樺(betula platyphylla)、山楊(populus davidiana dode)等。

1.2 無人機匹配點云數據

無人機影像的獲取時間為2015年8月17日至24日，采用大疆六旋翼無人機，搭載Sony NEX-5T相機，影像大小4 912×3 264，飛行高度300 m，影像分辨率8 cm，航向、旁向重疊度分別為90%、60%。

無人機影像處理采用Agisoft Photoscan軟件，主要步驟如下：

(1)加載包含可交換圖像文件格式(exchangeable image file format，EXIF)信息的無人機影像，若影像不包含EXIF信息，則需要進行攝像機定標。

(2)影像外方位元素計算，首先獲取影像上的特征點，然后進行影像間特征點匹配，通過影像間的匹配點使用“8點法”計算得到影像間的外方位元素信息。

(3)無人機影像匹配點云的生成，無人機影像匹配點云是利用上一步獲得的影像外方位元素，逐像素進行特征點的提取與同名點的匹配，然后使用光束法平差進行誤匹配噪聲點判別與去除。

(4)控制點的放置，通過放置控制點的方法，對影像匹配點云進行地理編碼。

(5)DSM的生成，為后續的單木識別算法評價，將無人機影像匹配點云進行柵格化，通過插值生成柵格大小為0.1 m的DSM數據。

(6)結果的導出，最終得到的平均點密度約為40個/m2的無人機影像匹配點云數據，柵格大小為0.1 m的DSM數據，數據坐標采用WGS84坐標系，投影為UTM投影北半球6°分帶的第51帶。

1.3 樣地實測數據

實測樣地數據獲取于2016年8月15日至24日。使用2015年獲取的無人機正射影像數據作為樣地布設參考，按樹種組成、株密度、高度等布設了8個直徑為30 m的圓形樣地。以正射影像作為底圖參考，通過手持GPS導航樣地位置，利用單木特征(如樹種、枯木、倒木等)對正射影像中的單木進行現場解譯，將解譯的單木進行編號并標注于正射影像圖中。實測的單木位置作為參考數據用于單木識別算法識別的精度評價，為了便于對單木識別算法進行評價，將實測單木中的不被遮擋，能夠目視識別的單木定義為“可見單木”，表1為各樣地實測信息。

表1 各樣地實測信息

2 研究方法

2.1 單木識別

單木識別有兩種常見錯誤：漏識別、錯識別。漏識別是指單木因為遮擋或靠近大樹等原因無法被識別；錯識別分為兩種，一種是因為噪聲點等原因將偽單木誤識別為單木；另一種是由于單木有較大分支或多個頂點而被識別為多株。本文提出的基于無人機影像匹配點云的單木識別算法，通過“較小的搜索窗口”獲取盡可能多的候選單木，然后對候選單木進行“樹冠結構分析”降低錯識別率。具體過程如下：

第1步，使用“局部最大值搜索”算法確定候選單木頂點。設定搜索窗口大小為d1，遍歷點云，判斷每個點在以其為中心的搜索窗口內是否為最高點，若為最高點則標記為候選單木頂點。

第2步，候選單木冠層點云水平分層切片。在原始點云中選取到候選單木頂點的水平距離小于d2的點云作為待分析點云。d2應大于研究區最大的單木冠幅半徑。對選取的單木待分析點云以厚度t自上而下進行水平分層切片，直到沒有點云為止，總層數記為n。對各層自上而下按1到n的順序進行編號(圖1)。

圖1 單木點云分層切片示例

第3步，逐層進行擬合分析，具體如下：

(1)第1層點云分析。在第1層點云中選取到候選單木頂點距離為小于d3的點進行圓形擬合，得到圓心坐標記為(x1，y1)、半徑為r1。

(2)第2至n層點云分析?？紤]到樹冠通常呈現一定的幾何體的形狀(如圓錐體，橢球體或球體)，樹冠點云會隨著層號的增加而逐漸外擴，因此點云的水平選取范圍參數d3需要根據層號和樹冠形狀進行調整，調整方法如下：

d3，i=r1+t×c×(i-1)

(1)

式中：i為層號，取值范圍為[2，n]；d3，i為d3在第i層的取值，r1為第一層擬合圓的半徑，t為層厚，c為冠幅外擴參數(即樹冠表面坡度的余切值)；每層進行圓形擬合得到圓心坐標為(xi，yi)，半徑為ri。

第4步，偽單木頂點剔除(圖2(a))。由于匹配點云不能穿透冠層，而只能刻畫冠層表面的特征，所以對于理想的樹冠而言，其冠幅會隨著層號的增大而逐步增大，直到達到最大值為止。當某一層擬合圓圓心位置較頂點位置發生明顯偏離，偏移距離超過擬合圓半徑或冠幅突然變小時，則可認定點云出現異常，記錄其上一層號為nc。對于偽單木而言，通常難以出現圓心保持基本不變而冠幅單調遞增的規律性現象。因而，可以通過冠層點云總厚度(即層數n)和點云異常層號nc兩個參數進行偽單木頂點識別，當n

第5步，臨近過識別單木頂點合并(圖2(b))。某些樹種會由于叢生、伴生或多個一級枝的存在等因素而造成過識別現象。因此在對所有候選單木頂點完成第2步至第4步的處理后，仍有可能會存在過識別的問題。但可以想見，對于由叢生、伴生或多個一級枝引起的過識別候選單木頂點，其擬合圓心會隨著層數的增加而逐漸聚攏?？梢愿鶕c層擬合圓心對初始單木位置進行改正，對于位置改正后距離小于d4的單木進行合并。d4為設定的閾值。

圖2 真偽單木與過識別單木示例

2.2 完整度參數

前一部分介紹了單木識別的方法，其中的關鍵在于擬合圓的圓心、半徑的應用，考慮到單木冠層水平切片得到的點云會圍繞著圓形或圓弧形存在一定的變動，且隨著層數的增加圓形或圓弧可能出現缺口，為了提高擬合圓的可信度，增加一個對點云分布圓形完整度的描述參數：完整度指數(completeness index，CI)。將平面按10度的間隔劃分為36個區間，如果某區間內有點云落入則取值為1，否則取值為0。36個區間的總和作為完整度指標。因此完整度最低為1，最高為36。當完整度低于18時，意味著點云所占據的圓弧長度低于擬合圓周長的一半，此時得到的擬合圓可信度降低。如前一部分第2步和第3步所述，對每木點云從第2層開始自上而下逐層分析，取2至nc層間最后一層完整度不低于18的層數作為最終的nc層。

3 實驗結果與分析

3.1 真偽單木樹冠結構展示

本研究提出的單木識別算法涉及到的參數設置：“局部最大值搜索”使用的搜索窗口大小d1=0.5 m，小于當地最小的冠幅值；單木點云選取中用到的點云到候選單木頂點的距離閾值d2=5 m大于當地最大的冠幅半徑；根據點云密度及當地樹種樹冠特點，點云水平切片厚度t=0.3 m，第1層點云到候選單木頂點距離閾值d3=1.5 m，冠幅外擴參數c=1.0，正確單木冠層點云厚度閾值n0=3，其單調遞增連續層數最小值nc0=3，單木間最小距離閾值d4=1.0 m。

圖3和圖4分別展示了真、偽候選單木冠層點云切片結果(圖3和圖4所展示的單木分別為圖5(a)中的編號為NO.21和NO.16單木)。圖3(a)與圖3(b)為該單木第4層和第7層結果，藍色點為候選單木中d2范圍內的點，紅色點為通過d3或d3，i截取的單木點云，擬合圓的半徑為“R”，藍色射線為點云圓環的完整度指數(completeness index，CI)，即為每層圓環的完整度。圖3(c)為該單木樹冠結構曲線圖，可以看到隨著層數的增加半徑先持續單調增加，直至由于點云不完整導致半徑發生突變；圖3(d)展示的候選單木的冠幅中心偏離初始候選單木頂點的距離很平穩，隨著深度增加，點云的不完整導致發生距離突變。

圖4(a)至圖4(b)分別對應第1層與第3層。由圖4(c)可知，偽單木的冠幅沒有隨著層數的增加而增加，相反呈現減少的趨勢，原因在于該偽單木為噪音點，不具有普通單木樹冠的基本特征，因而半徑變化無規律；同時，由圖4(d)可以看到，偽單木的冠幅中心偏離初始單木頂點的距離隨著層數的增加而增加，這不符合真實單木的特點；通過分析可以看出真實單木的冠層結構與偽單木的冠層結構存在明顯的差異，因而通過對候選單木冠層結構的分析，能夠對偽單木與真實單木進行判別。

圖3 真單木樹冠結構注：黑圈為點云擬合圓，R為擬合圓半徑；藍色點為單木5 m范圍內的點，紅色點為本株單木的點；藍色射線為圓環完整度間隔線；CI為完整度。

圖4 偽單木樹冠結構注：黑圈為點云擬合的圓，R為擬合圓半徑；藍色點為單木5 m范圍內的點，紅色點為本株單木的點；藍色射線為圓環完整度間隔線；CI為完整度。

3.2 單木識別分步結果

由上述方法對大興安嶺根河林區的8個樣地的無人機影像匹配點云數據進行處理，獲得每個樣地的單木識別結果。圖5所示為P1樣地單木識別的分步結果，底圖為P1樣地正射影像圖。圖5(a)將搜索窗口d1設為0.5 m，得到的候選單木頂點(紫色點)，可以看到，由于窗口設置較小，引入了較多的偽單木，如NO.9、NO.16、NO.40、NO.50、NO.53等；以及過識別單木，如NO.14與NO.19、NO.58與NO.59以及NO.67與NO.69等，NO.表示單木的編號。圖5(b)是根據候選單木樹冠結構分析，剔除掉偽單木后的結果(綠色點)，可以看到圖5(a)中錯識別的偽單木單木，如NO.9、NO.16、NO.40、NO.50、NO.53等被剔除，但也存在NO.18、NO.49這種小冠幅單木由于冠形特征不明顯被誤刪除的現象。圖5(c)對剔除掉偽單木后的單木進行單木位置改正后合并的結果(紅色點)，看以可到，圖5(b)的錯識別中的過識別單木，如NO.11與NO.14、NO.21與NO.22、NO.43與NO.44、NO.51與NO.53等，進行了合并。

圖5 單木識別分步結果

3.3 各樣地的單木識別結果

為了驗證本算法的效果，選取窗口大小分別為1.0 m、2.0 m的“局部最大值”算法，對每個樣地的DSM數據進行單木識別，其中當地的最小冠幅直徑為1.0 m，DSM柵格大小為0.1 m。圖6至圖13為8個樣地的單木識別結果，從左到右，分別是：實測可見單木(黑色點)、基于DSM數據窗口大小為1.0 m的局部最大值單木識別結果(亮綠色點)、基于DSM數據的窗口大小為2.0 m的局部最大值單木識別結果(紫色點)，本文算法得到的單木識別結果(紅色點)，其中藍色大圓為樣地邊界。

圖6 P1樣地結果

圖7 P2樣地結果

圖9 P4樣地結果

圖10 P5樣地結果

圖11 P6樣地結果

圖12 P7樣地結果

圖13 P8樣地結果

3.4 單木識別精度評價標準

單木既可以由單木樹冠頂點來表示，也可以由單木樹冠邊界輪廓來表示，因此存在“點”與“多邊形”兩種不同形式的結果。同樣，參考數據也存在點與多邊形兩種形式：一般通過GPS或全站儀測量單木樹干的坐標，形式為點；由于識別結果與參考結果的形式不同，其匹配的方法也不同。識別與參考結果的匹配主要存在4種方式：“點對點”“點對多邊形”“多邊形對點”“多邊形對多邊形”。

由于實測的可見單木位于單木樹冠的頂點，本文的3種單木識別算法得到的也是單木的頂點。因此本文采用“點對點”的匹配方式對單木的識別結果進行評價。將樣地內的實測可見單木頂點作為參考，將3種算法得到的單木頂點與參考單木頂點進行匹配。匹配的標準是：若1個實測可見單木頂點存在1個識別單木頂點與其對應，則為正確識別，存在0個識別單木頂點則為漏識別，存在2個以上頂點或0個實測單木頂點存在1個以上的識別單木與其對應，則為錯識別，具體如表2所示。

表2 單木識別結果類別

本文的單木識別結果使用信息檢索與統計學中的準確率、召回率和F測度進行評價，其表達式分別如下：

(2)

(3)

(4)

式中：Pd稱為準確率；Pr稱為召回率(或查全率)；F稱為F測度；Nc為正確識別單木總數；Nr為實測可見單木總數；Nd為識別出的單木總數。準確率Pd的含義為：在識別出來的所有單木中，真實單木所占的比例。Pr的含義為：識別出來的真實單木占總真實單木數的比例。F測度是對準確率和查全率的綜合描述，當能正確檢出全部真實單木時，F=100%；當檢出的單木全部為偽單木時，F=0%；F測度越高代表結果越好。

3.5 單木識別結果的精度評價

表3為8個樣地的單木識別結果統計對表，為便于敘述，“1 m窗口大小的局部最大值識別算法”稱為“A”算法，“2 m窗口大小的局部最大值識別算法”稱為“B”算法，“本文算法”稱為“C”算法。

從表3可知，8個樣地整體上的A、B、C 3種算法的準確率上A低于B接近20%，而B與C結果接近，說明隨著窗口的縮小，錯識別單木增加，而由于C成功實現偽單木的剔除，使得識別單木的準確性最高98.61%。在查全率上A高于B接近10%，C高于A接近9%，說明隨著窗口的縮小，單木漏識別問題得到抑制，這與前人的研究結論相符。A、B、C 3種方法的F測度分別為74.82%、77.35%與90.45%，即本文算法相比1 m、2 m窗口大小的局部最大值算法分別高近15%與13%，說明了本算法相比傳統的單木識別算法，能夠提高單木的識別精度。

針對稀疏樣地P2、P3、P5與P7樣地，除P5樣地，3種算法的查全率都比較高，均超過85%，原因在于稀疏樣地中單木較為獨立，樹冠相互遮擋的問題較輕，因而查全率較高，P5樣地雖整體較為稀疏，然而單木樹冠相連較多，使得有些單木樹頂點不明顯，導致查全率較低；3種算法的準確率差別較大，其中B與C的準確率都較高，而A算法的準確率較低，原因在于稀疏樣地存在較多的背景噪音被誤識別為偽單木，因而使得小窗口的A算法的精度偏低，而算法C能夠剔除偽單木，所以精度受影響較小。3種算法的F測度與準確率相似，由于小窗口的局部最大值算法，即A，引入了大量的偽單木使得其F測度低于B與C的F測度。

對于密集的樣地，P1、P4、P6、P8，3種算法的查全率都不高，都不大于80%，說明了在密集林區局部最大值單木識別算法存在明顯的局限性，原因在于密集林區的部分單木由于樹冠之間的相互遮擋，使得單木頂點不明顯；而3種算法的單木準確性都較高，原因在于密集林區的背景噪聲少。P6樣地比較特殊，是由于其樣地內部單木多伴生，樹冠相連，使得單木識別存在較嚴重的漏識別問題P1與P6樣地的內均存在大量的白樺(表1)，然而P1樣地的查全率與F測度均高于P6樣地，主要原因在于單木的空間分布特征，P1樣地雖然白樺較多，但都較為獨立，樹冠頂點較為明顯，因而易識別，而P6樣地內樹冠相連較多，使得單木樹冠頂點不明顯，因而降低了單木識別的結果。

表3 各樣地單木識別結果統計對比

注：A、B方法分別為1 m、2 m窗口大小的局部最大值法，C方法為本文的方法。

4 結束語

本文利用無人機影像匹配點云數據，進行了單木識別研究，針對單木漏識別與錯識別的問題，提出了聯合“局部最大值”與“單木樹冠結構分析”的單木識別算法。得出以下結論：(1)本算法相比于傳統的“局部最大值”算法，能夠獲得更精準的單木位置信息，為其他單木參數的獲取提供準確的初始參數；(2)對于郁閉林區(缺乏林下地形信息)，利用影像匹配點云獲取準確的單木位置與株數等信息，拓寬了固定窗口的“局部最大值”單木識別算法的應用范圍。

本文的研究還存在以下的問題：(1)由于本算法涉及多個閾值，而閾值的確定以實測數據的統計為基礎，這是本算法的不足之處，下一步考慮通過對無人機影像匹配點云的分析，實現閾值的自適應；(2)單木漏識別的問題，由于部分單木沒有明確的單木樹冠頂點導致局部最大值算法也存在飽和的現象，接下來將考慮綜合使用無人機影像點云的光譜與幾何信息，進行單木識別研究。