基于向量角度差和擬合曲線融合的建筑物點云提取方法

劉茂華, 陳杰, 陳晗琳

(沈陽建筑大學(xué)交通與測繪工程學(xué)院, 沈陽 110168)

建筑物的提取在城市規(guī)劃、數(shù)字制圖以及更新地理信息系統(tǒng)數(shù)據(jù)庫等眾多領(lǐng)域中具有重要意義。近年來,航空激光雷達(dá)(light detection and ranging,LiDAR)技術(shù)因其具有高密度和高精度的點云數(shù)據(jù)而迅速發(fā)展,并成為提取建筑物的一種替代方法[1]。

目前,已經(jīng)有大量的研究致力于建筑物提取。一些方法采用了將LiDAR數(shù)據(jù)與圖像數(shù)據(jù)結(jié)合的策略[2-4]。鄧飛等[5]利用方向和梯度等信息以及基于活動輪廓的圖割算法,成功提取了建筑物的輪廓信息。王春林等[6]綜合利用LiDAR數(shù)據(jù)和影像特征,在復(fù)雜場景中實現(xiàn)了建筑物輪廓的提取。盡管LiDAR和圖像結(jié)合的方法可以提高準(zhǔn)確性,但也存在一些問題,如在數(shù)據(jù)融合過程中需要設(shè)置合適的閾值,以及融合后可能導(dǎo)致精度降低等挑戰(zhàn)。因此,為了克服上述困難,研究人員開始關(guān)注從單獨的LiDAR數(shù)據(jù)中提取建筑物信息。

點云領(lǐng)域中的建筑物提取方法通常分為有監(jiān)督方法[7-8]和無監(jiān)督方法[9-10]。有監(jiān)督方法利用算法如支持向量機(jī)[11]、隨機(jī)森林[12]和深度學(xué)習(xí)[13-15]等。這些方法在點云分類方面具有一定的優(yōu)勢,但也存在一些限制,如參數(shù)調(diào)整和重復(fù)計算特征會降低準(zhǔn)確性和計算效率。無監(jiān)督方法包括形態(tài)學(xué)、區(qū)域生長、擬合方法和掃描線等技術(shù)。李昂等[16]結(jié)合梯度約束濾波、標(biāo)記分水嶺變換和最大似然分類等方法提取建筑物。趙傳等[17]通過提取初始建筑物輪廓點、計算點云分布特征、聚類和區(qū)域增長等步驟實現(xiàn)建筑物提取。李強(qiáng)等[18]通過地面LiDAR點云數(shù)據(jù)的特征分析和構(gòu)建規(guī)則集提取建筑物。朱軍桃等[19]通過點云投影、邊緣點提取、角點提取和邊界線還原等步驟提取建筑物輪廓。呂富強(qiáng)等[20]通過去噪、聚類和點云分割提取建筑物,實現(xiàn)了自動化的建筑物點云提取。

處理大規(guī)模高密度LiDAR數(shù)據(jù)的主要挑戰(zhàn)之一是計算負(fù)擔(dān)。航空LiDAR系統(tǒng)中的掃描線提供了適合使用(graphics processing unit,GPU)進(jìn)行并行計算的一維數(shù)據(jù)結(jié)構(gòu),使得處理高密度區(qū)域的LiDAR數(shù)據(jù)更加高效。一些研究已經(jīng)在其方法中利用了掃描線。如在點云濾波[21-22],道路標(biāo)線提取等[23]。Han等[24]和Hu等[25]基于掃描線成功的提取了建筑物,然而,在這種方法中,當(dāng)植被具有非常平坦的冠層或建筑物具有粗糙和不規(guī)則的屋頂表面時,可能會產(chǎn)生錯誤的結(jié)果。這是因為該方法對建筑物屋頂表面的粗糙度和多樣性作了不適當(dāng)?shù)募僭O(shè)。

綜上所述,盡管以前的工作在某些方面取得了進(jìn)展,但在功能和性能方面仍存在重要的限制。為了解決這些限制,提出一種基于向量角度差的屋頂提取算法。該算法利用航空LiDAR獲取數(shù)據(jù)時的掃描線與x-o-y平面垂直,并且掃描線上相鄰點具有相同維度的空間向量的特點,將三維空間向量簡化為二維方向向量,然后計算方向向量與z軸之間的角度,并設(shè)置閾值角度來提取建筑物的屋頂平面。該方法能夠在有效提取建筑物平面屋頂?shù)耐瑫r,降低計算復(fù)雜度。提出一種基于歐氏距離和多項式曲線擬合的特征提取方法,用于提取曲線建筑物屋頂。首先,計算掃描線上相鄰點之間的歐氏距離。然后,執(zhí)行基于最小二乘法的曲線擬合。由于植被點分布是無序的,而建筑物屋頂表面是有序的,擬合曲線的平均殘差值可以有效地提取建筑物曲面屋頂。

1 方法提出

通過LiDAR掃描線數(shù)據(jù)建立索引,采用k-維樹(k-dimension tree,KD-tree)搜尋算法對點云數(shù)據(jù)進(jìn)行去噪,同時依靠布料模擬算法(cloth simulation filter, CSF)對點云數(shù)據(jù)進(jìn)行濾波處理,達(dá)到提取非地面點云的目的。首先,分別以每一條掃描線為研究對象計算非地面點的相鄰方向向量并求差,利用樣本確定差值閾值并提取出平面建筑物屋頂點云;其次,利用剩余點云的相鄰點歐氏距離進(jìn)行曲線擬合,得到曲面建筑物屋頂點云;最后,經(jīng)過細(xì)化后處理,對提取后數(shù)據(jù)中的植物點進(jìn)行剔除,得到完整建筑物點云數(shù)據(jù)的結(jié)果。算法流程如圖1所示。

圖1 本文算法流程圖Fig.1 The proposed algorithm flowchart

1.1 掃描線模型

目前,機(jī)載LiDAR點云數(shù)據(jù)的獲取掃描方式分為:有線掃描、圓錐掃描以及顯微光學(xué)掃描。本文算法針對的數(shù)據(jù)模型是應(yīng)用最為廣泛的有線掃描式,這種掃描方式的激光腳點投影在地面上形成“Z”形。圖2展示有效掃描的數(shù)據(jù)形式。圖3展示了一條掃描線的剖面圖。

藍(lán)色為一條掃描線圖2 掃描線模型Fig.2 Scan line model

1.2 基于方向向量的平面屋頂點云提取

1.2.1 向量角度差

機(jī)載LiDAR獲取點云數(shù)據(jù)時,運動方向平行于x-o-y平面,因此每條掃描線所在的平面與x-o-y平面垂直,可以使用該平面上所有方向向量與z軸之間的角度來區(qū)分平面和非平面。如圖4所示,x-o-y平面即三維點云沿著z軸后的投影平面,則投影面法向量nm為(0,0,1),θi為相鄰點的方向向量ni與nm的夾角。

紅色點為建筑物屋頂點;黃色點為建筑物立面點;藍(lán)色點為地面點;綠色點為高大植被點;棕色點為低矮植被點圖3 掃描線剖面圖Fig.3 Scan line profile

圖4 向量角模型Fig.4 Vector angle model

黑色框1～6為建筑物樣本;黑色框7為曲面屋頂樣本;紅色圈1～6為植物樣本圖5 樣本數(shù)據(jù)來源Fig.5 Source of sample data

得到非地面點后依次計算相鄰點之間的方向向量夾角θi可表示為

(1)

則相鄰?qiáng)A角的差值為

Δθi=θi+1-θi

(2)

1.2.2 角度差閾值確定

選取7個建筑物和6個植物樣本進(jìn)行實驗計算,樣本來源如圖5所示。通過統(tǒng)計相鄰?qiáng)A角差值來確定區(qū)分平面與非平面的角度差閾值。建筑物樣本中包含平頂,斜頂,高程變化復(fù)雜以及四周突起等建筑物,能充分表達(dá)城區(qū)內(nèi)平頂建筑物的種類。

如圖6所示,建筑物樣本的相鄰向量角度差高度集中在0°～15°,差值波動的主要原因在于:一是建筑物屋頂表面不是完全平滑的;二是掃描儀本身存在高程值與真實值存在誤差。大于30°以上的差值是因為建筑物屋頂由多個平面構(gòu)成而存在高度差。分別以10°、15°和20°作為閾值統(tǒng)計角度差的概率。

圖6 6個建筑物樣本直方圖Fig.6 Histogram of 6 building samples

各樣本數(shù)據(jù)如表1所示,平頂建筑物的相鄰向量角度差高度集中在20°以內(nèi),而閾值分別設(shè)置為10°、15°以及20°時,平均概率分別為83.91%、95.96%和97.69%。以10°為基準(zhǔn),閾值設(shè)置為15°時,概率增長率(增加的概率與原本概率比值)為14.4%;閾值為20°時概率增長率為16.4%。說明當(dāng)閾值設(shè)定為15°時的有效性提升遠(yuǎn)高于閾值設(shè)定為20°,這意味著提取建筑物時的準(zhǔn)確率提升明顯。

表1 建筑物樣本數(shù)據(jù)

此外,為驗證閾值普適性,選取6簇植物群樣本進(jìn)行相同計算。圖7為6簇植物樣本的相鄰向量角度差,與平面建筑物不同的是,植物群向量夾角的分布是比較分散的,樣本中角度差不均勻分布在角度的各個分段,其中更多的是集中在角度中段且跨度較大。

與建筑物樣本取相同閾值進(jìn)行比較,詳細(xì)樣本數(shù)據(jù)如表2所示。3個閾值下植物點分布概率分別為2.96%、4.69%以及8.47%,以10°為基準(zhǔn),15°的增長率為58.5%、20°的增長率為186.1%。這意味著,隨著角度差的增大,提取建筑物時植物引起的誤差逐漸增大。綜合建筑物樣本與植物樣本的閾值結(jié)果,閾值從10°增大到15°,提取建筑物的準(zhǔn)確率增加明顯,且植物引起的誤差較低;閾值從15°增加到20°時,建筑物的提取準(zhǔn)確率增長變緩且植物引起的誤差增加明顯,故將提取平面屋頂?shù)南蛄拷嵌炔铋撝翟O(shè)置為15°。對于分類為非平面點的建筑物點和分類為平面點的植物點可以利用后續(xù)的細(xì)化后處理進(jìn)行修正。

表2 植物群樣本數(shù)據(jù)

1.3 基于最小二乘法的曲線擬合

鑒于在點云數(shù)據(jù)中無法通過向量夾角區(qū)分植物點和彎曲屋頂點,引入曲線擬合方法以解決該問題。建筑物屋頂點云通常具有人為設(shè)計的特征,其分布呈現(xiàn)規(guī)律的連續(xù)曲線,并且相鄰點之間的歐氏距離存在一定規(guī)律。與此相反,植物點的點云分布通常不規(guī)則,并且相鄰點之間的歐氏距離變化較大。因此,可以根據(jù)這一特點計算點間歐氏距離并依靠最小二乘法擬合曲線,根據(jù)擬合曲線的殘差進(jìn)行建筑物曲面屋頂?shù)奶崛　?/p>

1.3.1 擬合曲線殘差

(3)

式(3)中:I為多元函數(shù);xi、yi為二維空間的橫縱坐標(biāo)值;min為多元函數(shù)的最小值。

由多元函數(shù)極值的必要條件,得

(4)

式(4)的系數(shù)矩陣是一個對稱正定矩陣,存在唯一解。解得ak(k=0,1,…,n)可得多項式為

(5)

式(5)為所求的擬合多項式,稱為最小二乘擬合多項式,殘差的平方計算公式為

(6)

1.3.2 曲線擬合階數(shù)及殘差閾值確定

根據(jù)圖8(a)所示,所研究的建筑物的屋頂呈現(xiàn)出具有曲率變化的曲面,并且呈現(xiàn)出不規(guī)則的特征,從左至右曲率逐漸增大。與此對應(yīng),圖8(b)描繪了建筑物中一條掃描線上相鄰點之間的歐氏距離變化情況?？梢钥闯?相鄰點之間的歐氏距離變化趨勢與曲面屋頂?shù)男螒B(tài)變化趨勢相一致,即相鄰點之間的距離逐漸增大。因此,以該建筑物作為目標(biāo)樣本,以探究通過歐氏距離擬合曲線的閾值。

圖8 建筑物樣本掃描線及相鄰點歐氏距離Fig.8 Scanning line of building sample and Euclidean distance of adjacent points

該建筑物樣本包含掃描線共224條,以該建筑物樣本的相鄰點歐氏距離為樣本數(shù)據(jù),分別對每條掃描線數(shù)據(jù)建立4階擬合曲線,5階擬合曲線以及6階擬合曲線,并計算每組曲線的殘差。最終通過所有掃描線的殘差平均值確定閾值。如圖9所示,以其中一條掃描線為例,展示擬合曲線及殘差的確定方式。

圖9 建筑物樣本擬合曲線及殘差Fig.9 Fitting curve and residual error of building samples

圖9(a)為建筑物取線擬合結(jié)果。圖9(b)為該條掃描線依靠相鄰點間距建立4階擬合曲線的殘差值,即曲線上的值與真實值的誤差。圖9(c)為5階曲線的殘差。圖9(d)為6階曲線的殘差。以該條樣本數(shù)據(jù)為例,4階曲線的殘差高達(dá)0.1,絕對值的平均值為0.031。5階最大值為-0.051,絕對值平均值為0.021。6階曲線殘差的最大值為-0.043,絕對值的平均值為0.017。

用相同方法對224條建筑物掃描線計算擬合曲線以及殘差,詳細(xì)數(shù)據(jù)如表3所示,隨著擬合曲線的階數(shù)提升,殘差值降低的趨勢逐漸減緩。而224條掃描線計算總時間隨著階數(shù)的增加而快速增加。綜合考慮時間成本與準(zhǔn)確率,確定算法的擬合曲線階數(shù)為5階。為了確定殘差的具體閾值,采用相同方法對植物樣本進(jìn)行實驗。

表3 建筑物樣本數(shù)據(jù)

植物樣本共174條掃描線,仍然以一條掃描線為例,植物樣本的結(jié)果如圖10所示。4階擬合曲線的殘差值極小值為-0.125,絕對值的平均值為1.836,5階最小值為-0.373,絕對值平均值為1.358,6階殘差最小值為0.064,絕對值的平均值為1.164。統(tǒng)計174條植物樣本掃描線的5階擬合曲線殘差值,絕對值的平均值為1.846。

圖10 植物樣本擬合曲線及殘差Fig.10 Fitting curve and residual error of plant samples

根據(jù)樣本實驗可以得出建筑物與植物的歐氏距離擬合曲線的殘差值差異明顯,植物樣本的5階極小值仍然明顯高于建筑物殘差平均值,利用正態(tài)分布原則,殘差閾值采用3倍平均值即0.057。詳細(xì)計算過程如圖11所示。

L為初始窗口;L0為側(cè)窗口尺度;L1為當(dāng)前尺寸窗口圖11 擬合曲線算法過程Fig.11 Algorithm process of curve fitting

利用擬合曲線對掃描線中的曲面提取時,需要建立窗口來提取有效的數(shù)據(jù)點。初始窗口L可以根據(jù)點的數(shù)量代替?zhèn)鹘y(tǒng)的距離,這意味著能夠非常有效的分離出前后不連續(xù)的掃描點數(shù)據(jù)。首先,當(dāng)確定L的大小時,計算當(dāng)前窗口所有點間距的擬合曲線以及殘差并記錄;之后從窗口L的開始處進(jìn)行縮減窗口并建立擬合曲線計算殘差,當(dāng)殘差值的平均值不再減小時,固定該側(cè)窗口尺度L0;其次,從L尾側(cè)逐漸縮減窗口并記錄擬合曲線的殘差值,當(dāng)殘差值不再減小時,記錄當(dāng)前尺寸窗口L1。

根據(jù)融合向量角度差和擬合曲線法提取建筑物時,可能存在誤差,其中包括少量混入建筑物點云數(shù)據(jù)的植物點和被誤刪的建筑物點?？梢岳肒D-tree搜尋法來剔除散落的植物點。建立KD-tree結(jié)構(gòu)后從點云數(shù)據(jù)中隨機(jī)抽取點,計算點之間平均距離d以及全局平均距離μ和標(biāo)準(zhǔn)差σ,其中d在μ±3σ之外的被認(rèn)為是離散的植物點而去除。為了恢復(fù)過刪除的建筑物點,一條掃描線的建筑物提取后,重新歷遍這條掃描線,當(dāng)兩部分建筑物的距離為3σ之內(nèi)時,重新提取這部分點并比較與周圍建筑物點的高程,若在范圍內(nèi)則可以認(rèn)定是過刪除的建筑物點并恢復(fù),進(jìn)而提高建筑物點云提取的準(zhǔn)確性和完整性。

2 實驗分析

2.1 實驗數(shù)據(jù)

實驗數(shù)據(jù)通過運5飛行器搭載Leica ALS07傳感器掃描得到。傳感器脈沖頻率159 kHz,平均飛行高度為1 080 m,總面積達(dá)1.78 km2,平均點密度為4.6個/m2。如圖12所示,該數(shù)據(jù)具有典型的高建筑密度特性的城區(qū),植物與建筑物之間距離較短,且部分植物被修剪成規(guī)則形狀。

圖12 原始數(shù)據(jù)Fig.12 Raw data

2.2 實驗結(jié)果

建筑提取實驗使用MATLAB2019b軟件平臺,實驗計算機(jī)具有英特爾酷睿i5-8300H,2.3 GHz CPU,16.00 GB內(nèi)存和64位Windows 10操作系統(tǒng)。

(1)去噪和濾波結(jié)果。圖13為去噪和濾波的結(jié)果,濾波后點云數(shù)量為738 957。

圖13 濾波后結(jié)果Fig.13 Filtered results

(2)定性結(jié)果。如圖14所示,將降噪和濾波后的點云作為輸入點云,進(jìn)行建筑物點云。圖14(a)為所提算法提取結(jié)果,圖14(b)為TIN算法提取結(jié)果。

藍(lán)色框1～4為同一區(qū)域內(nèi)的不同結(jié)果;紅色圓圈1和紅色圓圈2為一些手動修剪后的植被圖14 兩種算法的建筑物提取結(jié)果Fig.14 Building extraction results of two algorithms

如圖14所示,與TIN算法相比所提算法可以更加全面的提取建筑物。如圖14中藍(lán)色框1和藍(lán)色框3所示,所提算法能夠更加完整的提取邊界形態(tài)復(fù)雜的建筑物。如圖14中藍(lán)色框4所示,TIN算法在面對高度變化大的平面屋頂時表現(xiàn)出不穩(wěn)定性,不能完整提取建筑物。圖14中藍(lán)色框2為一個曲面屋頂,所提算法成功提取了整個屋頂,而TIN算法幾乎漏掉了整個曲面屋頂。如圖14(b)所示,紅色圓圈1和紅色圓圈2為一些手動修剪后的植被,在形狀上表現(xiàn)出一定規(guī)則性,TIN算法無法有效去除這些植被。綜上可知,所提算法相較于TIN算法能夠更加穩(wěn)健的提取平面和曲面屋頂?shù)耐瑫r能夠有效去除植被點云。

2.3 精度及效率分析

為了驗證所提算法對建筑物提取的有效性,需對算法進(jìn)行精度評定,采用經(jīng)典模型混淆矩陣對數(shù)據(jù)結(jié)果進(jìn)行評價[26]。

定義一類誤差T1為把非建筑物點錯誤地分類為建筑物,可表示為

(7)

式(7)中:b為算法的建筑物數(shù)據(jù)中非建筑物點的個數(shù);e為參考數(shù)據(jù)的建筑物點個數(shù)。

定義二類誤差T2為把建筑物點錯誤的分類為非建筑物,可表示為

(8)

式(8)中:c為算法的非建筑物數(shù)據(jù)中建筑物點的個數(shù);f為參考數(shù)據(jù)的非建筑物點個數(shù)。

則總誤差T3可表示為

(9)

式(9)中:n為濾波后非地面點個數(shù)。

于是總體精度(overall accuracy, OA)可定義為

(10)

同時,為增加算法精度評價的準(zhǔn)確性,使用kappa系數(shù)進(jìn)行對算法準(zhǔn)確性分析。

(11)

(12)

式中:p0為正確分類的建筑物點a和非建筑物點d的總和除以樣本總數(shù);pe為偶然一致性;ai為建筑物和非建筑物的真實點云數(shù);bi為由算法分類的建筑物和非建筑物點云。

此類混淆矩陣中kappa系數(shù)替換為

(13)

式(13)中:g為算法分類的建筑物點數(shù);h為算法分類的非建筑物點數(shù)。

表4為3種算法在提取建筑物二元分類時的精度結(jié)果。

表4 兩組實驗結(jié)果精度評定

本文算法的T1為2.77%,而TIN算法的T1為7.21%。這表明本文算法在將非建筑物正確分類為負(fù)類方面表現(xiàn)更好,能夠更加有效區(qū)分非建筑物點云。TIN算法在T2指標(biāo)上略好于本文算法,但是本文算法T2也保持在較低水平,說明兩種算法都能較好的識別出建筑物點云。兩種方法的總體精度都高于95%,說明兩種方法在整體上都具有較高的分類準(zhǔn)確性。本文算法kappa達(dá)到了91.48%,而TIN的為87.38%。說明本文算法在不同樣本上具有較高的一致性和可靠性。

與TIN算法相比,本文算法表現(xiàn)出更快的處理速度,TIN算法運行時間為294 s,本文算法的運行時間為118 s,僅為TIN算法的40.1%。TIN算法為了更好地提取建筑物,需要建立較多三角測量數(shù)據(jù)結(jié)構(gòu),耗時較長,而本文算法基于掃描線模型,將它將三維空間向量簡化為二維方向向量,減少了冗余計算,提高了算法效率。

3 結(jié)論

提出了一種機(jī)載LiDAR建筑物點云提取方法,通過實驗區(qū)數(shù)據(jù)處理,結(jié)果表明,能夠有效提取建筑物點云。該算法具有以下特點:通過構(gòu)建掃描線數(shù)據(jù)結(jié)構(gòu)計算相鄰向量角度差能夠快速準(zhǔn)確地提取平面建筑物屋頂;利用相鄰點歐氏距離擬合曲線,計算殘差可以準(zhǔn)確區(qū)分曲面建筑物屋頂與植物點;最后通過細(xì)化處理能夠更加完整地提取建筑物并消除植物的影響。該過程可以很好避免植物對結(jié)果的影響,同時能夠非常有效的提取曲面屋頂以及高差復(fù)雜的建筑物屋頂。該算法解決了面對大面積城區(qū)建筑物復(fù)雜且植物多而提取精度不夠的問題,同時能夠很好提取曲面屋頂建筑物。通過成熟的TIN算法進(jìn)行消融實驗,驗證了本文算法的有效性和穩(wěn)定性。