王 亮，劉 鵬，賈旭斌

（1.青海省柴達(dá)木綜合地質(zhì)礦產(chǎn)勘查院，青海 格爾木 816099；2.青海省柴達(dá)木盆地鹽湖資源勘探研究重點(diǎn)實(shí)驗(yàn)室，青海 格爾木 816099）

區(qū)域增長(zhǎng)法是一種簡(jiǎn)單高效的分割方法，它通過選擇種子點(diǎn)將鄰近點(diǎn)聚集在一起，但容易受到種子點(diǎn)選取和增長(zhǎng)方法的影響[1]。本文采用一種分步提取屋頂面的方法，首先通過Alpha Shapes 算法確定屋頂面的邊界，并利用法向量和曲率特征選擇合適的屋頂面種子點(diǎn)，然后通過區(qū)域增長(zhǎng)算法提取出屋頂面，最后利用張量投票算法進(jìn)行面片后處理，從而實(shí)現(xiàn)從點(diǎn)云中提取屋頂面。

1 方法原理

本文分步提取屋頂面方法的流程如圖1所示。

圖1 本文流程圖

1.1 屋頂面輪廓線提取

Alpha Shapes[2]是一種可以有效提取離散點(diǎn)邊界的算法，該算法只需要設(shè)置滾動(dòng)圓半徑r，即可提取邊界點(diǎn)。如圖2 中，設(shè)置滾動(dòng)圓半徑為r，以點(diǎn)P1(x1,y1)和P2(x2,y2)為例，根據(jù)2點(diǎn)坐標(biāo)和半徑r計(jì)算圓心坐標(biāo)，i代表圓心R1和R2，則圓心坐標(biāo)表示為式（1）、（2）、（3）。

圖2 Alpha Shapes示意圖

計(jì)算點(diǎn)P1、P2周圍的點(diǎn)到圓心R1或R2距離，若到其中一個(gè)圓心的距離小于r的點(diǎn)不存在，到另一個(gè)圓心的距離小于r 的點(diǎn)存在，則表明此處為邊界。通過滾動(dòng)圓一次遍歷所有點(diǎn)，從而篩選出邊界點(diǎn)。

1.2 屋頂面種子點(diǎn)選取

法向量和曲率是提取建筑物屋頂面的重要依據(jù)。在三維空間中，垂直于平面的向量稱為該平面的法向量，由于點(diǎn)云在空間中是離散分布的，所以點(diǎn)云的法向量需要借助周圍鄰域點(diǎn)組成的局部平面來(lái)確定，鄰域點(diǎn)的選擇借助K-d Tree 實(shí)現(xiàn)，則點(diǎn)的法向量可以定義為該點(diǎn)和k個(gè)鄰域點(diǎn)構(gòu)成的最佳擬合平面確定的法向量[3-5]。設(shè)Pi為點(diǎn)云中一點(diǎn)，{Pi1,Pi2,…,Pik}為k個(gè)鄰域點(diǎn)，通過鄰域點(diǎn)的最小二乘擬合平面記為S，使所有鄰域點(diǎn)到平面S的擬合殘差平方和最小，平面S可表示為：

式中，為平面S 的法向量，d為坐標(biāo)原點(diǎn)到平面S的距離。

擬合平面會(huì)經(jīng)過鄰域點(diǎn)的重心，所以，擬合平面的求解可以轉(zhuǎn)換為鄰域點(diǎn)協(xié)方差矩陣特征值的求解。假設(shè)其對(duì)應(yīng)的協(xié)方差矩陣為C，鄰域重心為，則協(xié)方差矩陣C和鄰域重心表示為：

由公式（5）、（6）可求得協(xié)方差矩陣的特征值λi1、λi2、λi3，且λi1≤λi2≤λi3，特征值對(duì)應(yīng)的特征向量分別為，最小特征值對(duì)應(yīng)的特征向量即為點(diǎn)Pi的法向量，即

曲率則可以表示為：

屋頂面種子點(diǎn)的選取依據(jù)法向量和曲率值。種子點(diǎn)的單位法向量與鄰域點(diǎn)的單位法向量點(diǎn)積為1，與鄰近點(diǎn)的曲率足夠小，則將其作為屋頂面種子點(diǎn)用于屋頂面點(diǎn)云的聚類。

1.3 屋頂面面片提取

區(qū)域生長(zhǎng)算法是一種利用定義準(zhǔn)則由種子點(diǎn)向鄰域擴(kuò)張生成更大區(qū)域的算法。在種子點(diǎn)的基礎(chǔ)上，用法向量作為約束，向周圍鄰域點(diǎn)進(jìn)行擴(kuò)張，再以鄰域點(diǎn)為中心繼續(xù)向外生長(zhǎng)。圖3顯示以點(diǎn)P1為種子點(diǎn)向外生長(zhǎng)的示意圖，搜索8 連通方向的鄰域點(diǎn)，判斷鄰域點(diǎn)是否滿足約束條件，若滿足，則對(duì)鄰域點(diǎn)進(jìn)行標(biāo)記，再以標(biāo)記的鄰域點(diǎn)為中心，繼續(xù)搜索周圍未標(biāo)記的鄰域點(diǎn)，以此類推，直到所有的點(diǎn)檢索完成。

圖3 區(qū)域生長(zhǎng)示意圖

1.4 面片后處理

上述步驟中，我們已經(jīng)對(duì)提取的建筑物面片點(diǎn)和非面片點(diǎn)進(jìn)行了標(biāo)記，但是屋頂面交界處點(diǎn)云存在屋頂面競(jìng)爭(zhēng)的問題。因交界處點(diǎn)的法向量容易存在偏差，影響屋頂面交界處的點(diǎn)云提取，面片后處理主要是通過張量投票算法[6-7]實(shí)現(xiàn)，假設(shè)標(biāo)記為非面片點(diǎn)的點(diǎn)集為F={q1,q2,q3,…qi,…,qm}，首先檢索qi的鄰域點(diǎn)，利用鄰域點(diǎn)組成向量矩陣，并利用法向量張量對(duì)qi進(jìn)行投票，然后判斷qi與鄰域點(diǎn)構(gòu)成面片的可能性。當(dāng)qi與鄰域點(diǎn)構(gòu)成面片且qi的鄰域點(diǎn)為面片點(diǎn)時(shí)，則認(rèn)為點(diǎn)qi為面片點(diǎn)。

2 實(shí)驗(yàn)結(jié)果與精度分析

2.1 實(shí)驗(yàn)數(shù)據(jù)

為驗(yàn)證本文方法的可行性，選取國(guó)際攝影測(cè)量與遙感學(xué)會(huì)提供的4 組城區(qū)點(diǎn)云數(shù)據(jù)csite1、csite2、csite3和csite4進(jìn)行實(shí)驗(yàn)。其中，csite1數(shù)據(jù)有缺口，地形起伏較大，主要包含大量建筑物、樹木、道路等；csite2數(shù)據(jù)有缺口，包含建筑物、植被、道路、橋梁；csite3數(shù)據(jù)有缺口，包含建筑物、道路、植被；csite4數(shù)據(jù)有缺口，主要包含建筑物、道路和植被。因本文主要研究建筑物屋頂面點(diǎn)云的提取，所以對(duì)數(shù)據(jù)進(jìn)行了濾波處理，采用CSF濾波算法得到實(shí)驗(yàn)數(shù)據(jù)的非地面點(diǎn)云。

2.2 評(píng)價(jià)指標(biāo)

為了定量分析點(diǎn)云分割精度，選用完整率Cp、正確率Cr和提取質(zhì)量Cq對(duì)屋頂面提取結(jié)果進(jìn)行評(píng)價(jià)。

式中，TP為正確提取的屋頂面點(diǎn)個(gè)數(shù)；FP為錯(cuò)誤提取的屋頂面點(diǎn)個(gè)數(shù)；FN為未被提取的屋頂面點(diǎn)個(gè)數(shù)。

2.3 實(shí)驗(yàn)結(jié)果

本文采用上述方法對(duì)4 組實(shí)驗(yàn)數(shù)據(jù)濾波后的非地面點(diǎn)進(jìn)行處理，濾波的效果會(huì)對(duì)屋頂面種子點(diǎn)的篩選產(chǎn)生影響，如果有較多的地面點(diǎn)被錯(cuò)分為非地面點(diǎn)，那么篩選的屋頂面種子點(diǎn)中就會(huì)存在一部分錯(cuò)分的地面點(diǎn)，將地面點(diǎn)誤認(rèn)為是屋頂面點(diǎn)。

本文需要設(shè)置的參數(shù)主要有滾動(dòng)圓半徑r、法向量閾值ε和鄰域點(diǎn)K。滾動(dòng)圓半徑r 表示離散點(diǎn)提取邊界的精細(xì)程度，r越小提取的離散點(diǎn)邊界越精細(xì)。本文中依據(jù)點(diǎn)與點(diǎn)距離的平均值設(shè)置滾動(dòng)圓半徑r。法向量是以點(diǎn)所在局部平面的法向量作為該點(diǎn)的法向量，在計(jì)算法向量時(shí)一般以該點(diǎn)8鄰域或16鄰域內(nèi)的點(diǎn)所構(gòu)成的局部平面進(jìn)行計(jì)算，從而篩選出與鄰域點(diǎn)的夾角滿足條件的種子點(diǎn)，本文中設(shè)置的夾角閾值ε為15°，當(dāng)該點(diǎn)與鄰域點(diǎn)的夾角小于15°時(shí)認(rèn)為時(shí)平面點(diǎn)，可作為種子點(diǎn)。利用篩選的種子點(diǎn)進(jìn)行生長(zhǎng)向外擴(kuò)張，遍歷所有種子點(diǎn)鄰域內(nèi)的激光腳點(diǎn)，對(duì)滿足條件的點(diǎn)和不符合條件的點(diǎn)分別進(jìn)行標(biāo)記，直到所有的點(diǎn)檢索完成。鄰域點(diǎn)K通常設(shè)置為8。

如圖4中所示csite1局部區(qū)域的處理，計(jì)算點(diǎn)云的法向量和曲率特征，統(tǒng)計(jì)鄰域點(diǎn)法向量夾角大小并設(shè)置閾值條件篩選種子點(diǎn)，以此作為區(qū)域生長(zhǎng)的種子點(diǎn)，并以Alpha Shapes 提取的邊界點(diǎn)限制生長(zhǎng)的屋頂面面片范圍，最后采用張量投票對(duì)同一屋頂面的不同面片進(jìn)行合并，得到屋頂面點(diǎn)云。

圖4 csite1局部區(qū)域

4組實(shí)驗(yàn)數(shù)據(jù)的提取結(jié)果如圖5所示。

圖5 數(shù)據(jù)csite1、csite2、csite3、csite4的建筑物屋頂面提取結(jié)果

圖5中屋頂面的提取結(jié)果表明，本文方法能很好地從點(diǎn)云數(shù)據(jù)中提取出建筑物屋頂面。為了進(jìn)一步驗(yàn)證本文算法提取屋頂面的效果，利用完整率、正確率和提取質(zhì)量作為評(píng)價(jià)因子，與RANSAC算法、區(qū)域增長(zhǎng)算法提取屋頂面的效果進(jìn)行對(duì)比，可知，本文方法可以在完整率Cp上達(dá)到與RANSAC算法、區(qū)域增長(zhǎng)算法相近的精度，但在正確率Cr和提取質(zhì)量Cq上，本文方法表現(xiàn)出明顯的優(yōu)勢(shì)，正確率Cr達(dá)到95%以上，提取質(zhì)量Cq在94%以上，主要在于2種算法在提取屋頂面點(diǎn)云時(shí)存在過度分割或欠分割的現(xiàn)象，如圖6所示，以數(shù)據(jù)2中部分區(qū)域?yàn)槔故玖薘ANSAC方法、區(qū)域增長(zhǎng)方法與本文方法提取的屋頂面效果對(duì)比圖，從圖中可以看出RANSAC方法和區(qū)域增長(zhǎng)方法在進(jìn)行屋頂面提取時(shí)，存在相鄰屋頂面欠分割的情況，RANSAC方法在提取小面積屋頂面邊緣時(shí)存在錯(cuò)誤提取的問題。

圖6 屋頂面提取效果對(duì)比圖

通過定量、定性的對(duì)比本文方法與其他2 種方法的提取效果，說(shuō)明本文方法可以較好地實(shí)現(xiàn)屋頂面的提取，且能達(dá)到很好地提取效果。

3 結(jié) 語(yǔ)

本文提出的分步提取屋頂面點(diǎn)云的方法，能較好地實(shí)現(xiàn)屋頂面地提取，滿足實(shí)際工程的需要。但本文方法也存在步驟較多、數(shù)據(jù)處理連貫性不強(qiáng)的缺點(diǎn)，且實(shí)驗(yàn)未對(duì)本文算法處理數(shù)據(jù)的效率進(jìn)行評(píng)估，計(jì)算效率方面的研究還需要繼續(xù)跟進(jìn)，在今后的研究學(xué)習(xí)中將繼續(xù)改進(jìn)和完善其中的不足之處。