基于機(jī)載LiDAR點(diǎn)云和建筑物輪廓線構(gòu)建DSM的方法

李 遷，肖春蕾，陳 潔，楊達(dá)昌

(1．中國國土資源航空物探遙感中心，北京 100083;2．中國地質(zhì)大學(xué)(北京)，北京 100083)

0 引言

LiDAR(light detection and ranging，LiDAR)系統(tǒng)包含激光雷達(dá)、全球定位系統(tǒng)(GPS)和慣性導(dǎo)航系統(tǒng)(INS)。激光脈沖不受陰影和太陽角度影響，能夠快速、直接且連續(xù)自動地獲取地面三維數(shù)據(jù)。這些數(shù)據(jù)經(jīng)過簡單的處理(如粗差剔除、格網(wǎng)化)，便可以得到一種重要的空間數(shù)據(jù)——數(shù)字表面模型(digital surface model，DSM)。

激光點(diǎn)云為不規(guī)則的三維離散點(diǎn)數(shù)據(jù)，可通過采用逐點(diǎn)內(nèi)插的方法建立DSM。常用的內(nèi)插方法有鄰近距離內(nèi)插和三角網(wǎng)線性插值等。鄰近距離內(nèi)插算法能保留建筑物和周圍地面的差異，得到高精度的建筑物信息［1］，但得到的DSM會出現(xiàn)鋸齒化的建筑物邊緣［2］，不能應(yīng)用于正射影像、城市三維建模等的生產(chǎn);三角網(wǎng)線性插值算法雖然在大部分情況下能滿足幾何精度的要求，但因沒有顧及點(diǎn)云所表達(dá)的地物之間存在的高程關(guān)系［3］，存在同一個(gè)三角網(wǎng)同時(shí)“穿越”了地面和建筑物，會有明顯的三角面的出現(xiàn)，造成定位精度不準(zhǔn)［4］?；谏鲜鰡栴}，本文提出了基于建筑物輪廓線構(gòu)建DSM的方法。構(gòu)造Delaunay三角網(wǎng)(以下簡稱D－三角網(wǎng))時(shí)嵌入多邊形約束條件，即利用三維激光地面點(diǎn)和建筑物輪廓線構(gòu)建不規(guī)則三角網(wǎng)(triangulated irregular network，TIN)。實(shí)驗(yàn)證明，基于該方法構(gòu)建的DSM能較為精細(xì)地表達(dá)建筑物邊緣，可用于高精度DSM的生產(chǎn)。

1 基于建筑物輪廓線的DSM構(gòu)建方法

建立不規(guī)則三角網(wǎng)的基本過程是將最鄰近的3個(gè)離散點(diǎn)連接成三角形，同時(shí)考慮地性線(如建筑物輪廓矢量線)、地物等特征線對格網(wǎng)的影響。為了保證DSM格網(wǎng)最大限度地符合實(shí)際地形，應(yīng)用中通常把地性線等地形特征線作為TIN中三角形的邊［5］。機(jī)載LiDAR系統(tǒng)獲取的離散三維點(diǎn)包括地面點(diǎn)、人工建筑物(房子、煙囪、塔、輸電線等)及自然植被(樹、灌木、草)等，對點(diǎn)云進(jìn)行濾波處理可分離出地面點(diǎn)和非地面點(diǎn)，建筑物點(diǎn)云可以從非地面點(diǎn)中檢測出來，繼而可以提取建筑物點(diǎn)云深度影像中的建筑物邊緣。在構(gòu)建不規(guī)則三角網(wǎng)時(shí)加入建筑物矢量邊緣線，能使DSM中的建筑物信息表達(dá)得更精確。本文基于建筑物輪廓線構(gòu)建DSM的流程如圖1所示。

圖1 基于建筑物輪廓線構(gòu)建DSM流程Fig．1 Flow chart of constructing DSM based on building contour lines

該流程主要包括以下4個(gè)步驟:

1)點(diǎn)云的濾波處理。該步驟可為下一步建筑物點(diǎn)云的檢測做準(zhǔn)備。本文采用Terra Scan軟件對數(shù)據(jù)進(jìn)行濾波，將點(diǎn)云分為地面點(diǎn)和非地面點(diǎn)2類。該軟件采用的濾波算法為不規(guī)則三角網(wǎng)漸進(jìn)濾波算法，具有很強(qiáng)的斷線檢測能力，適用于地物比較復(fù)雜的城區(qū)，能成功地濾除大多數(shù)的建筑物信息。

2)建筑物點(diǎn)云的檢測。在上述濾波處理得到非地面點(diǎn)的基礎(chǔ)上，基于點(diǎn)云的高程紋理信息檢測建筑物腳點(diǎn)。從激光點(diǎn)云的空間分布特征出發(fā)，本文設(shè)定3個(gè)參數(shù)閾值(平面距離閾值R、高程閾值H以及角度閾值?)，用以提取建筑物的平面屋頂信息，以確保下一步提取建筑物邊緣的準(zhǔn)確性。

3)建筑物邊緣的提取及規(guī)則化?；谝越ㄖ稂c(diǎn)云生成的深度影像圖，利用Canny算子提取建筑物邊緣，用方位角聚類規(guī)則法對建筑物邊緣規(guī)則化，得到規(guī)則化的建筑物輪廓線。

4)構(gòu)建約束D－三角網(wǎng)，重采樣生成DSM。首先不考慮約束多邊形的影響，由激光點(diǎn)云地面點(diǎn)數(shù)據(jù)構(gòu)建非約束D－三角網(wǎng)，然后將建筑物輪廓線作為約束數(shù)據(jù)入網(wǎng)，對建筑物輪廓線內(nèi)部的三角形進(jìn)行清空處理，生成DSM。

2 數(shù)據(jù)處理中的主要算法

2．1 基于高程紋理的建筑物檢測

在機(jī)載LiDAR數(shù)據(jù)中，紋理可定義為因局部區(qū)域高程變化而產(chǎn)生的對比度、均勻性等物理特征，即高程紋理［6］。由圖2可以看出，建筑物與鄰近點(diǎn)的高差比較固定，沒有明顯的高程突變(圖2的P1，P2點(diǎn))，可以近似看作連續(xù)的水平平面;植被間的高差較大，且不規(guī)則(圖2的P3，P4，P5點(diǎn))，其與鄰近點(diǎn)之間的斜率很大，很難看作近似的水平面。

圖2 LiDAR點(diǎn)云空間分布Fig．2 Spatial characteristics of LiDAR point clouds

為了提高建筑物點(diǎn)云提取的精確性，本文在提取建筑物點(diǎn)云的過程中加入鄰近點(diǎn)的斜率信息，通過設(shè)定一定的參數(shù)來擬合這種近似平面(圖3)［7］。

圖3 基于高程紋理提取建筑物點(diǎn)示意圖Fig．3 Detecting building points based on height texture

圖3上P點(diǎn)為建筑物中心點(diǎn);P1為相鄰點(diǎn);為P1到中心點(diǎn)P所在平面的投影點(diǎn);S為P1P之間的空間距離;d為P1P之間的平面距離;h為P1P之間的高差;θ為P1點(diǎn)與P點(diǎn)所在平面的夾角。設(shè)定高差閾值H，平面距離閾值R和角度閾值?，判定P1與P是否在同一個(gè)平面上需要滿足3個(gè)條件，即d＜R且|h|＜H，θ＜?。

一般高差閾值H和平面距離閾值R分別設(shè)為平均高差和平均平面距離的1～2倍;傾角閾值?設(shè)為10°～20°。擬合傾斜平面所需要的參數(shù)閾值一般大于平面的參數(shù)閾值。過小的閾值會增大平面點(diǎn)的遺漏誤差，過大的閾值將增大平面的侵蝕作用。為了避免某些少量滿足條件的點(diǎn)被歸為平面，需設(shè)定閾值N，只有滿足條件的P點(diǎn)所有鄰近點(diǎn)的數(shù)目大于N，才能歸為平面，認(rèn)為P點(diǎn)為建筑物點(diǎn);反之，為非平面，即P點(diǎn)為其他地物點(diǎn)。

2．2 建筑物輪廓線提取及其規(guī)則化

建筑物邊緣特征是構(gòu)成建筑物模型的重要信息，能夠描述建筑物結(jié)構(gòu)和幾何特征，大多數(shù)建筑物的特征線均為直線，檢測到建筑物的各個(gè)直線特征就意味著檢測到了建筑物的主要結(jié)構(gòu)信息。因此，建筑物直線特征的檢測是建筑物提取的重要工作。一般常用的邊緣檢測算子有梯度算子、Laplace算子、Sobel算子、Canny算子以及方向算子等非線性算子，此外還有曲面擬合法等等。其中，Canny算子具有方向性，能更好地應(yīng)用于邊緣強(qiáng)度估計(jì)，能產(chǎn)生梯度方向和強(qiáng)度2個(gè)信息［8］。因此，本文利用Canny算子檢測建筑物腳點(diǎn)深度影像中的建筑物邊緣。

本文采用Ruijin［9］提出的方位角聚類比較法規(guī)則化建筑物輪廓。該方法結(jié)合了聚類和調(diào)整2個(gè)過程，聚類方法類似于K－means方法。其算法描述如下:

1)計(jì)算建筑物線段各自的方位角，根據(jù)方位角將線段分為2類。求這2類方位角的平均值a和b，比較候選線段與這2個(gè)平均值的差別，差異小的為該線段所屬的類別。經(jīng)處理后，這2類線段間應(yīng)為互相垂直關(guān)系。

2)在每一類中計(jì)算方位角的權(quán)重均值，因?yàn)檩^長線段比較短線段方位角精度高，所以權(quán)重取決于線段的長度，即

式中:li為其中一類中第i條線段的長度(i=1，2，3，…，n);αazimi為第i條線段的方位角;αazim為其中一類的方位角。

3)使用Gauss－Markov模型調(diào)整2類權(quán)重，使得2類線段嚴(yán)格垂直。這里的權(quán)重為這2條線段的總長度。

4)將每一條線段繞該線段中點(diǎn)旋轉(zhuǎn)至該類方位角。到此為止，建筑物的輪廓線段之間不是垂直，就是平行。

5)設(shè)定平行線段間距離閾值d(可取1 m)，當(dāng)2條線段間距離小于該值，則融成一條新線段。通過計(jì)算各線段中點(diǎn)的權(quán)重均值獲取新線段的中心點(diǎn)。比如，當(dāng)2條線段融合成新的線段時(shí)，x坐標(biāo)求取公式為

式中:l1和l2分別為2條線段的長度;x1和x2為2條線段的中心點(diǎn)坐標(biāo)。

該算法的優(yōu)勢在于無需較好的建筑物初始輪廓，就能獲取規(guī)則的建筑物邊界。也就是說，不需要非常準(zhǔn)確的邊緣檢測算子即可提取出建筑物輪廓。

2．3 TIN 的構(gòu)造

對于約束D－三角網(wǎng)剖分算法，根據(jù)約束邊嵌入時(shí)機(jī)的不同可分為2類:第1類是指在構(gòu)網(wǎng)的同時(shí)考慮約束邊的影響，直接構(gòu)建約束D－三角網(wǎng);第2類是指首先不考慮約束邊的影響，構(gòu)建數(shù)據(jù)集的非約束D－三角網(wǎng)，在已構(gòu)好的三角網(wǎng)中強(qiáng)行嵌入約束邊以構(gòu)建約束D－三角網(wǎng)，即2步法［10］。本文采用2步法進(jìn)行構(gòu)網(wǎng)，嵌入算法的步驟如下:

1)按照一定規(guī)則對數(shù)據(jù)區(qū)域進(jìn)行格網(wǎng)劃分，基于逐點(diǎn)內(nèi)插算法對離散地面點(diǎn)集生成D－三角網(wǎng)。

2)確定線段的影響區(qū)域。如圖4所示，首先確定線段AB的端點(diǎn)A所在的任一△A，然后由△A來確定線段方向上的首三角形(記為△S);在△A(如△APB)中，確定點(diǎn)A的對應(yīng)邊b(如PB)與有向線段AB相交否，若相交，則首三角形為△A(△APB=△S);如果判斷邊c(如DE)在有向線段AB的右側(cè)，則以端點(diǎn)A為圓心逆時(shí)針方向搜索與△A相鄰的下一個(gè)三角形;若為左側(cè)，則以順時(shí)針方向搜索。

圖4 搜索線段影響區(qū)域Fig．4 Detecting affecting area of the line

圖5 影響區(qū)域的三角剖分Fig．5 Delaunay triangulation by incremental insertion algorithm for affecting area

3)影響區(qū)域的三角剖分。如圖5所示，以有向線段AB作為擴(kuò)展邊，在擴(kuò)展邊右側(cè)影響區(qū)域(圖5中陰影部分)的點(diǎn)集中取一點(diǎn)C，使得C點(diǎn)與擴(kuò)展邊的兩端點(diǎn)的連線組成的夾角最大(即最大角準(zhǔn)則)，生成新的△ABC。同時(shí)，在有向線段AB左側(cè)構(gòu)網(wǎng)時(shí)，則需要將有向線段AB改為BA。將多邊形的每條線段都重復(fù)上述過程，最終實(shí)現(xiàn)多邊形入網(wǎng)。

4)多邊形內(nèi)部清空處理。從多邊形內(nèi)部出發(fā)，根據(jù)拓?fù)潢P(guān)系，向八方向輻射，搜索位于多邊形內(nèi)部的三角網(wǎng)并將其移除。圖6(a)為待移除位于內(nèi)部多邊形的三角網(wǎng);圖6(b)為重新構(gòu)造的D－三角網(wǎng)。

圖6 嵌入約束多邊形前(左)后(右)D－三角網(wǎng)的變化示意圖Fig．6 Comparison of constrained delaunay triangulation(left)and delaunay triangulation(right)

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

3．1 實(shí)驗(yàn)結(jié)果

實(shí)驗(yàn)區(qū)位于廣西壯族自治區(qū)柳州城區(qū)，其正射影像如圖7所示。區(qū)內(nèi)有多種地物，如植被、道路、汽車及大小高度形狀各異的建筑物等。實(shí)驗(yàn)數(shù)據(jù)由ALS50－II系統(tǒng)獲得，LiDAR點(diǎn)云密度為8點(diǎn)/m2左右，數(shù)據(jù)的垂直精度優(yōu)于15 cm，水平精度優(yōu)于0．5 m，總點(diǎn)數(shù)為1 976 040。圖8為實(shí)驗(yàn)區(qū)激光點(diǎn)云高程設(shè)色圖。

圖7 正射影像Fig．7 Orthophoto

圖8 高程設(shè)色圖Fig．8 Display at height

首先，對點(diǎn)云進(jìn)行預(yù)處理。利用TerraScan濾波得到地面點(diǎn)數(shù)據(jù)，實(shí)驗(yàn)區(qū)屬于丘陵地區(qū)，角度閾值設(shè)為6．0°;根據(jù)測區(qū)的實(shí)際建筑物面積，將建筑物最大邊長設(shè)為160 m;構(gòu)建三角形過程中的高差閾值設(shè)為1．4 m。當(dāng)所加點(diǎn)構(gòu)成的三角形每條邊短于5 m時(shí)，阻止向三角形內(nèi)部加點(diǎn)。濾波后獲取的地面點(diǎn)如圖9所示?？梢钥闯?，對于植被和建筑物等高程突變比較明顯的地物來說，TerraScan濾波方法具有比較好的效果。

圖9 地面點(diǎn)Fig．9 Ground points

圖10 濾除地面點(diǎn)的nDSMFig．10 nDSM after filting ground points

其次，進(jìn)行建筑物點(diǎn)的檢測。首先獲取去除地面點(diǎn)的nDSM(normalized DSM)(圖10)，利用3 m的高度閾值濾掉低矮灌木叢，為下一步的建筑物平面擬合減少誤差;然后基于高程紋理提取建筑物點(diǎn)，令閾值R=4 m，H=1 m，?=15°，N=40。建筑物點(diǎn)的提取結(jié)果如圖11所示。

圖11 建筑物點(diǎn)云Fig．11 Building points

然后，提取建筑物矢量輪廓線。首先，基于Canny邊緣檢測算法提取由建筑物點(diǎn)云所生成的nDSM深度圖像中建筑物邊緣(圖12);然后，利用方位角聚類比較法規(guī)則化建筑物輪廓(圖13);最后，在規(guī)則化的建筑物二維輪廓線上任取一點(diǎn)，由離散三維點(diǎn)云內(nèi)插出其高程，將高程值賦與此建筑物。應(yīng)用該方法得到的建筑物輪廓線與本區(qū)域的正射影像疊合如圖14所示?？梢钥闯?，提取的建筑物邊緣基本與正射影像上的建筑物邊緣重疊。

圖12 Canny算子提取的建筑物邊緣Fig．12 Building edges detected by Canny

圖13 建筑物邊緣規(guī)則化Fig．13 Building edges normalized

圖14 建筑物輪廓與DOM疊合Fig．14 Combination of building contour and DOM

最后，利用濾波得到的離散三維地面點(diǎn)與建筑物矢量輪廓線構(gòu)建TIN格網(wǎng)，這里格網(wǎng)間距設(shè)為1 m，得到DSM的深度圖像(圖15)。

圖15 DSM深度圖像Fig．15 Depth image of DSM

3．2 誤差分析

如圖14所示，實(shí)驗(yàn)區(qū)域內(nèi)共有71棟建筑物，其中有66棟被檢測出來，對于建筑物頂部為平面的規(guī)則建筑物檢測效果比較明顯，而對于不規(guī)則和尖頂房有漏檢的現(xiàn)象。將提取的結(jié)果與通過TerraScan軟件(鄰近點(diǎn)高差閾值設(shè)為0．2 m)自動提取的建筑物結(jié)果相比較(圖16，17)，可以看出，本文方法對于建筑物的提取更為精確，建筑物邊緣有很好的保留，幾乎沒有建筑物邊緣點(diǎn)被誤分為高植被點(diǎn)的現(xiàn)象。

圖16 TerraScan提取的建筑物點(diǎn)和高植被點(diǎn)Fig．16 Building points and high vegetation points detected by TerraScan

圖17 本文方法提取的建筑物點(diǎn)和高植被點(diǎn)Fig．17 Building points and high vegetation points detected by the paper method

為了定量描述建筑物的分類誤差，本文以TerraScan的半自動、半手工分類結(jié)果為參考數(shù)據(jù)，計(jì)算本文方法的分類誤差(表1)。實(shí)驗(yàn)結(jié)果表明，本文方法的分類誤差小于10%，說明方法比較精確有效。

表1 分類誤差分析Tab．1 Error analysis of classification

比較本文方法構(gòu)建的DSM(圖18)及采用逐點(diǎn)內(nèi)插得到的規(guī)則格網(wǎng)(格網(wǎng)間距也為1 m)所構(gòu)建的DSM(圖19)建筑物細(xì)節(jié)，可以得出，本文構(gòu)建的DSM建筑物邊緣更準(zhǔn)確和規(guī)則。

圖18 本文方法構(gòu)建的DSM建筑物邊緣Fig．18 Building edge of DSM extracting by proposed method

圖19 規(guī)則格網(wǎng)構(gòu)建的DSM建筑物邊緣Fig．19 Building edge of DSM based on regular grid

建筑物矢量輪廓線的高程信息是由二維邊緣線上任意一點(diǎn)內(nèi)插得到的，所以構(gòu)建的DSM存在一定的誤差。本文利用抽稀建筑物激光點(diǎn)云內(nèi)插得到的DSM高程誤差DZ如表2所示。高程平均誤差DZ平均=0．018 5 m，均方根誤差RMSE=0．069 5 m。表明利用本文方法所構(gòu)建的DSM建筑物邊界信息比較精確，高程精度比較高。

表2 DSM高程誤差Tab．2 Evaluation of DSM accuracy

4 結(jié)論

本文提出了基于建筑物輪廓線構(gòu)建DSM的方法。在TIN構(gòu)造時(shí)，采用2步法嵌入建筑物輪廓線約束多邊形，從而得到格網(wǎng)間距為1m的DSM。主要結(jié)論如下:

1)通過定性誤差分析，本文方法構(gòu)建的DSM較規(guī)則格網(wǎng)構(gòu)建的DSM更能準(zhǔn)確表達(dá)規(guī)則的建筑物邊緣信息;通過定量高程誤差分析，本文方法構(gòu)建的DSM很好地避免了無約束條件構(gòu)建D－三角網(wǎng)引起的定位不準(zhǔn)，所構(gòu)建的DSM建筑物邊緣更為精確。

2)為了獲取準(zhǔn)確的建筑物邊緣信息，本文基于高程紋理檢測三維建筑物腳點(diǎn)，以TerraScan軟件進(jìn)行的半自動、半手工分類結(jié)果作為參考數(shù)據(jù)進(jìn)行了定量的誤差分析，誤分點(diǎn)比例為8．09%，說明結(jié)果較為理想。

3)建筑物輪廓規(guī)則化是在建筑物的特征線為直線的前提下進(jìn)行的，對于復(fù)雜建筑物邊緣的細(xì)化，還需要進(jìn)一步地研究。

［1］ 熊俊華，方源敏，付亞梁，等．機(jī)載LiDAR數(shù)據(jù)的建筑物三維重建技術(shù)［J］．科學(xué)技術(shù)與工程，2011，11(1):189－192．Xiong JH，F(xiàn)ang Y M，F(xiàn)u Y L，et al．The research on 3D reconstruction of buildings based on airborne LiDAR data［J］．Science Technology and Engineering，2011，11(1):189－192．

［2］ Shen W．Building boundary extraction based on LiDAR point clouds data［C］//Li Z L．Anthology of International Society for Photogrammetry and Remote Sensing．England:CRC Press，2008:153－157．

［3］ 鄔建偉，馬洪超，李 奇．顧及語義的機(jī)載LiDAR點(diǎn)云網(wǎng)格化方法［J］．測繪科學(xué)技術(shù)學(xué)報(bào)，2008，25(2):237－240．Wu JW，Ma H C，Li Q．LiDAR data cloud gridding based on semantics［J］．Journal of Geomatics Science and Technology，2008，25(2):237－240．

［4］ Yang B，ShiW，Li Q．An integrated TIN and grid method for constructing multi－resolution digital terrain models［J］．International Journal of Geographical Information Science，2005，19(10):1019－1038．

［5］ 古林玉．機(jī)載LiDAR點(diǎn)云構(gòu)建高精度DSM的關(guān)鍵技術(shù)研究［D］．焦作:河南理工大學(xué)，2010．Gu L Y．Research of constructing high－accuracy DSM using airborne LiDAR points［D］．Jiaozuo:Henan Polytechnic University，2010．

［6］ Elberink S O，Mass H G．The use of anisotropic height texture measures for the segmentation of airborne laser scanner data［C］//International Archives of Photogrammetry and Remote Sensing．Amsterdam，Netherlands:Science Press，2000:678－684．

［7］ 徐花芝．基于航空LiDAR點(diǎn)云數(shù)據(jù)的建筑物提取研究［D］．西安:長安大學(xué)，2008．Xu H Z．Study on the building extraction from airborn LiDAR points cloud data［D］．Xi’an:Chang’an University，2008．

［8］ 龔 亮，李正國，包全福．融合航空影像的LiDAR地物點(diǎn)云分類［J］．測繪工程，2012，21(2):34－39．Gong L，Li ZG，Bao Q F．Classification of LiDAR object points by fusing aerial image［J］．Engineering of Surveying and Mapping，2012，21(2):34－39．

［9］ Ruijin M．Building model reconstruction from LiDAR data and aerial photographs［D］．Ohio:Ohio State University，2004．

［10］ 劉學(xué)軍，龔健雅．約束數(shù)據(jù)域的Delaunay三角剖分與修改算法［J］．測繪學(xué)報(bào)，2001，30(1):82－88．Liu X J，Gong JY．Delaunay triangulation of constrained data set［J］．Acta Geodaetica et Cartographica Sinica，2001，30(1):82－88．