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

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

0 引言

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

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

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

目前缺少針對無人機影像匹配點云數(shù)據(jù)的單木識別算法，因此我們借鑒基于LiDAR數(shù)據(jù)的“局部最大值”算法進(jìn)行單木識別研究?！熬植孔畲笾怠彼惴ǜ鶕?jù)搜索窗口的不同分為兩種：一種是“固定窗口”的局部最大值；另一種是“變化窗口”的局部最大值。變化窗口的局部最大值通常需要精確的林下地形信息，而無人機影像點云在郁閉林區(qū)難以獲得精確的林下地形信息，因而難以采用變化窗口的局部最大值法進(jìn)行單木識別[15]。

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

為此，本文提出了一種聯(lián)合“局部最大值”與“樹冠結(jié)構(gòu)分析”的單木識別算法，該算法可解決 “局部最大值”算法的窗口大小設(shè)置的困境。

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

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)

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

1.2 無人機匹配點云數(shù)據(jù)

無人機影像的獲取時間為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信息，則需要進(jìn)行攝像機定標(biāo)。

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

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

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

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

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

1.3 樣地實測數(shù)據(jù)

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

表1 各樣地實測信息

2 研究方法

2.1 單木識別

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

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

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

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

第3步，逐層進(jìn)行擬合分析，具體如下：

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

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

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

(1)

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

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

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

圖2 真?zhèn)螁文九c過識別單木示例

2.2 完整度參數(shù)

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

3 實驗結(jié)果與分析

3.1 真?zhèn)螁文緲涔诮Y(jié)構(gòu)展示

本研究提出的單木識別算法涉及到的參數(shù)設(shè)置：“局部最大值搜索”使用的搜索窗口大小d1=0.5 m，小于當(dāng)?shù)刈钚〉墓诜担粏文军c云選取中用到的點云到候選單木頂點的距離閾值d2=5 m大于當(dāng)?shù)刈畲蟮墓诜霃?；根?jù)點云密度及當(dāng)?shù)貥浞N樹冠特點，點云水平切片厚度t=0.3 m，第1層點云到候選單木頂點距離閾值d3=1.5 m，冠幅外擴參數(shù)c=1.0，正確單木冠層點云厚度閾值n0=3，其單調(diào)遞增連續(xù)層數(shù)最小值nc0=3，單木間最小距離閾值d4=1.0 m。

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

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

圖3 真單木樹冠結(jié)構(gòu)注：黑圈為點云擬合圓，R為擬合圓半徑；藍(lán)色點為單木5 m范圍內(nèi)的點，紅色點為本株單木的點；藍(lán)色射線為圓環(huán)完整度間隔線；CI為完整度。

圖4 偽單木樹冠結(jié)構(gòu)注：黑圈為點云擬合的圓，R為擬合圓半徑；藍(lán)色點為單木5 m范圍內(nèi)的點，紅色點為本株單木的點；藍(lán)色射線為圓環(huán)完整度間隔線；CI為完整度。

3.2 單木識別分步結(jié)果

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

圖5 單木識別分步結(jié)果

3.3 各樣地的單木識別結(jié)果

為了驗證本算法的效果，選取窗口大小分別為1.0 m、2.0 m的“局部最大值”算法，對每個樣地的DSM數(shù)據(jù)進(jìn)行單木識別，其中當(dāng)?shù)氐淖钚」诜睆綖?.0 m，DSM柵格大小為0.1 m。圖6至圖13為8個樣地的單木識別結(jié)果，從左到右，分別是：實測可見單木(黑色點)、基于DSM數(shù)據(jù)窗口大小為1.0 m的局部最大值單木識別結(jié)果(亮綠色點)、基于DSM數(shù)據(jù)的窗口大小為2.0 m的局部最大值單木識別結(jié)果(紫色點)，本文算法得到的單木識別結(jié)果(紅色點)，其中藍(lán)色大圓為樣地邊界。

圖6 P1樣地結(jié)果

圖7 P2樣地結(jié)果

圖9 P4樣地結(jié)果

圖10 P5樣地結(jié)果

圖11 P6樣地結(jié)果

圖12 P7樣地結(jié)果

圖13 P8樣地結(jié)果

3.4 單木識別精度評價標(biāo)準(zhǔn)

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

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

表2 單木識別結(jié)果類別

本文的單木識別結(jié)果使用信息檢索與統(tǒng)計學(xué)中的準(zhǔn)確率、召回率和F測度進(jìn)行評價，其表達(dá)式分別如下：

(2)

(3)

(4)

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

3.5 單木識別結(jié)果的精度評價

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

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

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

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

表3 各樣地單木識別結(jié)果統(tǒng)計對比

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

4 結(jié)束語

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

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