Volume10, 1, March 2005, pp. 27-46
整合光達點雲與空照影像重建三維建物模型
賴彥中
1陳良健
2饒見有
3摘要
本研究從資訊融合角度出發,結合光達資料及數位化彩色空照影像進行三維建物模型之重建。研 究重點分成兩個部分,第一部份是房屋區塊偵測,第二部份是房屋模型重建。房屋區塊偵測部分,依 原始光達資料之高度資訊,判別出地上物區域。利用空照影像之色彩分析,將地上物區域之樹木部分 去除。殘餘之區塊設定一面積門檻,將面積過小之區塊濾除。之後進行區塊形狀以及表面分析,去除 不合理之區塊,房屋區塊即建立完成。房屋模型重建部分,首先進行光達資料前處理,預估出概略屋 緣位置,同時萃取出位於屋頂平面之光達點。將預估之概略屋緣位置投影至空照影像上,設定工作區 域,在此區域內進行精確之直線偵測。所得之直線線段組合成房屋之候選屋頂面,搭配前處理中萃取 出之光達點進行屋頂面之判別與高度計算後,投影回物空間獲得三維之屋頂面,最終結合 SMS 方法建 立房屋模型。研究中使用新竹科學園區之光達資料及大比例尺數位空照影像進行測試,並使用一千分 之一數值地形圖以及人工量測航測立體對產生之房屋模型做為檢核。
關鍵詞︰光達點雲、空照影像、資料融合、房屋偵測,房屋重建,房屋模型
1. 前言
三維都市模型已被廣泛運用在各種領域上,其
中房屋模型更為必要的單元。傳統上以遙測的方法 來重建三維房屋模型時,是採用立體製圖的方式,透過空照立體像對,經由人工量測的方式從點到 線,再用線組成面建構房屋模型。雖然經由立體製 圖可得到高品質之房屋模型,但此過程須耗費大量 人力與時間,因此以電腦運算代替人工量測是目前 研究值得研究的方向。然而電腦對於影像辨識能力 遠遠未及人工辨識程度,只應用影像資訊全自動的 產生三維房屋模型仍有許多困難。因此在利用影像 重建房屋模型之應用中,有加入人工的判別以半自
動之方式補足電腦辨識能力之不足,如(Gruen &
Dan, 1997)。
自 從 空 載 雷 射 掃 瞄 系 統 (Airborne Laser Scanning System) 或 稱 光 達 系 統 (LIDAR, LIght Detecting And Ranging)問世之後,對於地表之描述 更為詳盡,且在房屋資訊之萃取方面提供了一個新 的方向。光達系統之發展可回溯到一九七零年代 (Ackermann, 1999),藉著發射雷射波以及精準的量 測回波所經過的時間進而得知感測器與目標物之 距離。光達系統使用連續雷射波,經由相位差獲得 目標物與感測器間的距離,此類系統籍由全球定位 系統(GPS,Global Positioning System)以及慣性導航 系統(INS,Inertial Navigation System)之輔助可得目
1國立中央大學土木工程所空間資訊組碩士
2國立中央大學太空及遙測研究中心暨土木工程學系教授
3國立中央大學太空及遙測研究中心航太技術師
收到日期:民國 93 年 07 月 08 日 修改日期:民國 93 年 11 月 22 日 接受日期:民國 93 年 11 月 24 日
標之精確三維座標(Morgan & Tempfli, 2000),加上 密集之掃瞄頻率,對於地表之描述在精確度與解析 力上均有重大之突破(Wehr & Lohr, 1999),且對於 三維房屋模型重建具有自動化之潛力。近年來應用 光達資料於模型重建之研究文獻有許多,如以網格 化之光達資料重建房屋模型(Alharthy & Bethel, 2002),或是結合光達與其他輔助資料,如結合空 照影像(Rottensteiner & Briese, 2002)等。
空照影像或是光達資料上,對於地表之描述各 有優勢。如光達資料具有高程精度佳以及平面萃取 容易之特性等,空照影像具有平面解析力佳邊界萃 取容易及含有光譜資訊之特性等(Schenk & Csatho, 2002)。藉由此特性,本研究將研擬一套流程,以 資訊融合(Data fusion)為出發點,結合數位空照影 像與光達資料之優勢。期能以最少人工之介入,進 行三維房屋模型的重建,測試將以可靠度及精度為 主要考慮。
本研究中執行之步驟上,首先以分割之方法,
區分出原始的資料中房屋所在位置。再將範圍縮小 至房屋之區域,搭配空照影像與原始光達資料,進 行三維房屋模型重建。目前之研究針對房屋分佈較 離散之地區進行實驗,假設房屋之輪廓皆由直線線 段構成(Haala & Brenner, 1998),且限制屋頂為平頂 構造之房屋,在此設定之處理範圍內,期望能重建 房屋之概略輪廓以及屋頂之突出物。研究之測試區 位於新竹科學園區一帶,包含光達資料以及高解析 力之數位化彩色空照影像,測試將以可靠度及精度 為主要考慮。
2. 房屋區塊建立
本研究使用之資料中,光達資料為一描述地表 起伏之三維點雲,而空照影像則為具有地表覆蓋物 之光譜反射資訊之影像。所以在進行房屋模型重建
的工作時,首先要在資料中確立出進行模型重建之 地區,也就是偵測出資料中屬於房屋的區塊,再進 而將目標鎖定於房屋區塊內進行模型重建。
在進行房屋區塊偵測是以網格化後之光達資 料進行分析。實際分析採取下列六步驟,分別是(1) 光達資料網格化、(2)高度分析、(3)色彩分析、(4) 面積分析、(5)形狀分析以及(6)表面分析。經過此 分析步驟,可將光達資料中屬於房屋的區塊切割出 來,在之後房屋模型重建之步驟便可將目標鎖定於 房屋區塊內。
2.1 光達資料網格化
光達之原始資料為三維點雲,在進行房屋區塊 偵測處理時將之轉換成為規則之網格化資料,利用 影像處理之觀念進行分析。光達原始點雲(all points) 利用濾波器過濾雜訊後,以其多重回波(multiple returns)之特性(Wehr & Lohr, 1999),可初步分類之 地面點(ground points)與地表點(surface points),便 可網格化得兩份資料,分別為數值地表模型(DSM) 與數值高程模型(DTM) (Briese et al., 2002)。
在進行網格化時,由於地面之起伏較平緩,不 同於房屋之牆面,針對地面點以及地表點進行網格 化時採用不同之方法。為求效率,DTM 之產生是 採 用 簡 易 之 反 距 離 權 重 法 (Inverse Distance Weighting, IDW)(Shepard, 1968),將最靠近網格中 心之 16 個點,以其與網格中心之距離倒數為權,
進行加權平均可得 DTM。此加權平均後之 DTM 符合地表起伏較平緩之特性。對於產生 DSM,則 是採用最近點內插法(Nearest neighbor interpolation) (Richards,1986)。以此方法進行內插時並沒有對於 光達點之高度進行平均,而是取最鄰近該網格中心 之光達點高度取代網格中心之高度,因此網格化後 之 DSM 可保留房屋邊緣之高差資訊。
2.2 高度分析
獲得 DSM 與 DTM 之後,DSM 代表的是所有 地表覆蓋物之高度資訊,例如房屋、樹木以及地面 等,而 DTM 代表的純粹是地面之資訊。可依兩組 資料之高度差異,判斷出光達資料上屬於地上物之 區域(above ground area)。
進行地上物區域判別時,是以網格為基礎 (pixel-based)。將 DSM 與 DTM 中相對應之網格高 度相減,並設一高度差門檻值。若相對應之網格相 減後高度差大於此門檻值,則該網格點視為地上物 之區域。
由於樹木於資料中也具有一定之高度,高度分 析之結果除了將人工建物萃取出之外,同時也會將 樹木之區域保留下來。在進行完此步驟時,將這些 殘留之區域稱為地上物區域。
2.3 色彩分析
在進行光達資料之高度分析後,所得之地上物 區域中尚包含大量之樹木區域。由於樹木在空照影 像中含有特定之顏色,且與房屋之顏色有一定差 異,在本步驟加入空照影像中之光譜反射資訊來濾 除地上物區域中之樹木區域。
色彩分析之步驟,首先在彩色空照影像上,依 照光譜反射之特性,進行影像分類,判別出屬於綠 色樹木之區域。進行影像分類時,採用監督式分類 (Supervised classification)(Richards, 1986)之方法,
目的是要將樹木區從影像中區分出。在進行監督式 分類時,樹木區域之灰度值通常較低,若只是針對 樹木以及非樹木區域進行分類,則在陰影區容易造 成分類錯誤。為了避免此錯誤,在進行分類時,是 將影像分成四類,分別是樹木區域、樹木陰影區 域、非樹木區域、非樹木陰影區域四類。再將前兩 類結合成樹木區域,將後兩類結合成非樹木區域,
得一張分類之影像。
已偵測出之地上物區域,是以物空間(X,Y,Z) 之座標表示,而分類影像是以像空間之像座標(x,y) 表示,影像之外方位參數已知,可經由共線條件式 完成兩類資料之對應。將光達前處理判別出之地上 物區域,以網格中心之地面座標,搭配該網格點之 高程,透過共線條件式投影至分類影像上,將地上 物區域之網格與分類影像套合。若地上物區域之網 格對應於分類影像上為樹木之區域,則將該點於地 上物區域中刪除,最後得一個濾除樹木區域之光達 高度影像。
2.4 面積分析
在高度分析以及色彩分析時,是以單一網格為 基礎,針對每個網格進行分析與判斷。由此步驟開 始,將色彩分析後殘餘之網格群聚(clumping)成為 獨立區塊。以獨立區塊為單位(region-based),進行 之後步驟之分析。
進行色彩分析時,因光達資料與數位空照影像 資料取樣時間不同、套合誤差之影響或是分類之誤 差等,使得無法針對樹木區完全濾除。真實世界中 房屋之面積具有一定之大小,可藉由區塊之面積大 小來確立屬於房屋之區塊。在此步驟設一面積門 檻,區塊面積大於此門檻,則該區塊為候選之房屋 區塊,小於面積門檻之區塊即在光達資料中刪除。
2.5 形狀分析
此步驟是針對區塊之形狀進行分析,以濾除不 合乎房屋形狀特性之區塊。對於形狀判斷之指標為 面積周長比,若房屋區塊為狹長之形狀,則周長除 以面積之值會較大。若房屋區塊之形狀較為緊密 (compactness),則周長除以面積之值便來的較小。
因此便利用此指標對於候選之房屋區塊進行判斷。
將通過面積門檻之區塊,計算其邊緣之網格數
(周長),以及區塊網格總數(面積),以周長除以面 積為指標來判別房屋區塊之形狀是否符合房屋之 外型。此步驟必須針對周長除以面積之值設定一門 檻,若通過門檻之區塊,則在下一步驟進行表面分 析,作為最終判斷房屋之依據。
2.6 表面分析
一般房屋之屋頂多為平面組合而成,與樹木表 面高低起伏明顯之特徵有顯著差異。表面分析是針 對區塊之表面之粗糙度進行評估,判斷候選之房屋 區塊是否由平面所構成,作為最終確定房屋區塊之 依據。
分析時以一個 3*3 之移動視窗,如下圖 1,對 於候選房屋區塊內網格之高度進行運算,計算該視 窗中心點位四個方向之斜率差並取最小值。若運算 過後之值大於所設定門檻之網格數,佔整個候選房 屋區塊之網格某比例以上時,代表此區塊之高低起 伏明顯,並不符合房屋區塊為平面構成之特徵,因 此將此區塊濾除(邵怡誠及陳良健,2003)。計算如 式(1)。
S1 = abs(H1-H0) – abs(H5-H0) S2 = abs(H2-H0) – abs(H6-H0) S3 = abs(H3-H0) – abs(H7-H0) S4 = abs(H4-H0) – abs(H8-H0) dS = min(abs(Si))
(1)
圖 1. 房屋模型重建流程圖
3. 房屋模型重建
從資料中獲得房屋之位置以及概略形狀後,將 目標鎖定於房屋區域內,於小範圍內進行模型重
建。本研究對於房屋模型重建之概念,首先於影像 空間中偵測屋緣線段,利用線段構成房屋屋頂面。
屋頂面加入光達之高度資訊計算出其高度後,將屋 頂面投影回物空間組成房屋模型。
在實際操作流程上,以前步驟之房屋區塊偵測 成果為起始,將相對應之空照影像、DTM 以及光 達原始點雲切割出。利用光達原始點雲中之高度資 訊判斷出概略屋緣之位置後,將概略屋緣位置投影 至影像上,設定直線偵測工作區域。運用影像對於 房屋邊界萃取容易之特性,在影像空間中偵測得屋 緣線段,並利用線段之幾何關係組合成房屋之候選 屋頂面。候選屋頂面加入光達資料之高度資訊後,
對於候選屋頂面進行判別並計算屋頂面之高度。最 終將像空間之屋頂面投影回物空間,利用 SMS 法 (Chen & Rau, 2003)組合房屋模型。模型重建流程 圖如下圖 2。
真實世界之房屋之外型,大多具有規則且明顯 的幾何特性,如房屋屋緣為直線且且相鄰邊界成正 交。在本研究中對於模型重建之處理範圍為房屋輪 廓為正交之直線線段組成,且為平頂構造之房屋。
在此設定之處理範圍內,期望能重建房屋之概略輪 廓以及屋頂之突出物。
圖 2. 房屋模型重建流程圖
房屋區塊 空照影像
光達原始 點雲
光達資料前處理
屋緣直線線段偵測
候選屋頂面建立
屋頂面判別與高度計算 像空間
物空間 SMS法房屋模型重建 房屋區塊
DTM 物空間
3.1 模型重建輸入資料
房屋模型的重建流程,以一例作為說明。如圖 3 所示,流程初始以房屋區塊偵測之成果出發。利 用已偵測出房屋區塊之平面位置,如圖 3(a),切割
出由光達資料中之地面點網格化後之數值高程模 型,如圖 3(b),以及光達原始點雲,如圖 3(c)。在 影像空間則依照共線條件式將該區塊對應之空照 影像切割出,如圖 3(d)。運用這四筆相對應之資料 當作模型重建之初始資料。
(a) (b)
(c) (d) 圖 3. 模型重建輸入資料示意圖
(a)房屋區塊,(b)房屋區塊 DTM,(c)房屋區域光達點雲,(d)房屋空照影像
3.2 光達資料前處理
在本步驟中,將光達原始點雲進行初步之處 理。目標是要利用光達之高度資訊預估出房屋概略 屋緣之位置以及萃取出位於屋頂平面上之光達 點。概略屋緣區域可於之後直線偵測時設定直線偵 測之工作區域,而屋頂面之光達點則可用來判別候 選屋頂面之正確性以及計算屋頂面之高度。
在進行前處理時,是將房屋區塊附近之光達點 雲切割出後,先在二維空間組合成三角網群。獲得 三角網群後,給予每個三角網頂點高度,便可獲得
具三維資訊之三角網群。前處理時,以每個三角網 為單位,利用三角形網頂點之高度資訊以及房屋區 塊之平面位置預估出房屋概略屋緣之位置,並加入 房屋區塊 DTM 以萃取出位於屋頂平面之光達點。
3.2.1 概略屋緣預估
房屋屋緣在真實世界中具有明顯之高度差 異,此差異可由光達資料中明顯觀察出。在本步驟 中由光達點雲組合後之三角網群出發,進行高度判 斷後,預估出概略屋緣位置,在之後步驟投影至影 像上將直線偵測之範圍所小。
概略屋緣預估之方法,是將房屋區塊附近之光 達點雲切割出後,組合成不規則三角網群。之後以 每個獨立之三角網為單位,設定兩判斷準則如下:
(1) 三角網之三頂點最高與最低點之差大於某門 檻值。
(2) 三角網之三頂點必須位於房屋區塊以及設定之 房屋區塊緩衝區(buffer)內。
通過上述兩準則之三角網群所涵蓋之區域,則視為 概略屋緣位置。房屋區域原始點雲組合成之三角網 群如圖 4(a),預估出之概略屋緣位置如圖 4(b)。
(a) (b)
圖 4. 概略屋緣預估示意圖
(a)房屋區塊三角網(顏色代表三角網高度),(b)概略屋緣位置(含顏色之三角網群)
3.2.2 屋頂平面光達點萃取
本步驟之目的為萃取出位於屋頂平面之光達 點,於之後投影至影像上,對候選屋頂面進行判別 以及高度計算。在萃取具有屋頂平面特性之光達 點,利用組合後之三角網群,對於每個三角網設定 三判斷準則如下:
(1) 三角網之三頂點最高與最低點之差小於某門
檻值。
(2) 組成三角網之三頂點必須位於房屋區塊以及 設定之房屋區塊緩衝區內。
(3) 組成三角網之三頂點必須高於地表所設定門 檻之上。
若通過此準則之三角網群之頂點便視為具有 屋頂平面特性之光達點。圖 5(a)為房屋區域原始光 達點雲,圖 5(b)為萃取出屋頂平面之光達點。
(a) (b) 圖 5. 屋頂面光達點萃取示意圖
(a)萃取前之光達點雲,(b)萃取後位於屋頂平面之光達點
3.3 屋緣直線線段偵測
本研究在模型重建的必要過程為偵測影像上 屋緣之直線線段。若是在整張空照影像上偵測直
線,則如道路邊緣以及房屋陰影等不屬於屋緣之線 段都會被偵測出。利用前步驟預估出之概略屋緣位 置,在影像空間中設立直線偵測之工作範圍,可減 少直線偵測所產生之雜訊。
光達資料預估出概略邊緣位置後,由於光達資 料之平面解析力較空照影像低,無法獲得精確之屋 緣平面位置。在此步驟中,將概略位置投影至影像 空間中,運用空照影像平面解析力佳以及直線特徵 萃取容易之特性,在影像空間中進行精確之屋緣線 段偵測。
本研究中精確屋緣直線偵測,主要包含五步 驟,分別為(1)工作區域建立、(2)影像平滑處理、
(3)梯度強度影像偵測、(4)細化及線追蹤以及(5)直 線偵測。經過直線偵測之步驟,可在影像空間上獲 得房屋精確之屋緣直線線段。
3.3.1 工作區域建立
此步驟之目的,是在影像上設定一工作區域,
將直線偵測之範圍鎖定於此區域內,以減少其他雜 訊之干擾。光達資料前處理中預估出之概略屋緣位 置,是以三角網群所涵蓋之區域。在建立工作區域 時,以構成概略屋緣位置之每個三角網為單位,投 影至影像空間中,在影像空間中同時形成由多個三 角形構成之區域。精確之直線偵測步驟,便在此工 作區域內進行。
將三角網群投影至像空間時,是將每個三角網 頂點之平面座標(X,Y),搭配三頂點內最高點之高 程(Zmax)投影至像空間。目的為去除高差移位之影 響,使得投影至影像空間之三角網群能涵蓋房屋屋 緣。
3.3.2 影像之平滑處理
影像上之雜訊可能源自影像品質不佳或非目 標之細小地物之影響,易受到此雜訊之干擾。進行 房屋邊界偵測時,容易得到含有雜訊之特徵,影響 後續房屋特徵萃取與辨識工作之進行。
最簡單之平滑雜訊之方法是以平均值(mean) 來濾除雜訊,但此法往往也將邊界資訊平滑化,以 致特徵弱化,造成邊界損失。針對抑制雜訊以及邊
緣資訊保留兩種考量因素,本研究採用最鄰近勻化 法(Symmetric Nearest Neighbor Mean, SNN-Mean) (Harwood et al, 1987)。採用此方法之優點在濾除雜 訊的同時又可保持建築物邊緣資訊的反差,以供後 續進行直線特徵萃取之用。
SNN-Mean 原理是以一個 n × n 之移動視窗(n 為奇數),在影像中由左至右、由上而下對影像進 行過濾,其演算法如下:
(1) 以中心像元為原點,在 n × n 之視窗中可建構 出(n2-1)/2 組對稱像元。
(2) 找出對稱像元中與中心像元之灰階值最接近 者,可得(n2-1)/2 個灰階值。
(3) 計算此(n2-1)/2 個像元之平均灰階值,並取代中 心像元。
以一個 5*5 之視窗為例,如圖 6 所示,可取得 12 組對稱像元對。將對稱像元對中灰度值較接近 中心像元者取出並平均後,取代中心像元,即可對 影像進行過濾。假設一個邊界經過此視窗,從對稱 像元對中挑選出之點與中心像元皆位於邊緣之同 一側,因此在進行平均時並不會破壞邊緣之反差資 訊。
圖 6. SNN-Mean 示意圖
3.3.3 梯度強度影像偵測
此步驟之目的為偵測影像上灰度值之反差強 度。由於房屋的邊界之灰度值在影像上大多會有明 顯之變化,先將反差大之地區偵測出,為具有屋緣 特徵之位置。為求效率,實際運算採用 Sobel edge detection (Parker, 1997),以 3×3 之移動視窗,由左
邊界位置
對稱像元對
至右、由上至下對影像進行計算。演算法如下:
(1) 計算影像中灰度值 Pn 與 Gx、Gy 兩個 filter 的 乘積,相加後取絕對值,如圖 7。
(2) 將所計算得之灰度值取代 3 × 3 視窗之中心像 元。
圖 7. Sobel edge detection 示意圖
3.3.4 細化及線追蹤
經由 Sobel edge detection 所偵測到之梯度強度 影像,代表著空照影像中灰度值之變化量,具有強 度且具有寬度。在梯度強度影像上,必須進行線條 細化之動作,才能提供之後直線偵測之用。本研究 中 採 用 局 部 最 大 值 追 蹤 法 ( 饒 見 有 及 陳 良 健 , 1997),在已偵測出之梯度強度影像上,沿著局部 最大值進行追蹤,可將梯度強度影像細化成只有單 個像元之寬度,以供之後直線偵測之用。
局部最大值追蹤法演算法如下:
(1) 選取種子點:
影像上由左至右、由上而下選取種子點,種子 點必須具有一定灰度值,並且至少在兩對稱像元對 上之強度為最大。
(2) 探詢新種子點:
由種子點向八方向探詢,進行局部最大值之追 蹤,若搜尋到之點比相鄰兩點之強度大,且通過設 定之灰階強度,則設定此點為邊緣特徵點,同時記 錄在影像上之位置,並繼續追蹤下個點,進而完成 整條線段之細化。
(3) 記錄分叉點:
若有一點具有兩個以上之邊界點,也就是分叉 點(branch point),則記錄此點作為下一線段之起始 點。
(4) 線段長度檢核:
尋找完所有符合條件之特徵點,且相連結之特 徵點像元數必須通過給定之最大像元數,若大於則 保留線段,反之則刪除。
3.3.5 直線偵測
Hough(1962)提出了一種方法,可用來辨識影 像中特定之圖形,這些圖形可以特定的數學模式來 表現出來,例如直線或是圓弧。在本研究中主要是 針對直線之偵測,將影像上所有經過細化線追蹤後 之特徵點,經由數學模式轉換至參數空間後,將所 對應之參數值進行群集偵測,可得到擁有直線特徵 之特徵點之參數,再轉換為影像空間,得到影像上 直線之線段。Hough 轉換之另個優點為抗干擾能力 強,如果帶偵測之線段有小的擾動或是斷裂,甚至 是需線之偵測,進行轉換後在參數空間仍能得明顯 之峰點。
於直線偵測時,是將所得之直線之參數於房屋 區域之影像內延伸,目的為利用線段將影像分割成 許多區域,於下步驟建立房屋之候選屋頂面。且在 直線偵測之過程中,為了避免過多之雜訊產生,因 此對於所偵測出之直線線段設定三項約制條件。約 制內容如下:
(1) 直線偵測出之線段必須具有邊界特徵點數某 門檻之上。
(2) 直線偵測出之線段必須與最明銳之線段角度 差在某門檻內。
(3) 直線與直線之間相距在某距離以上。
直線偵測步驟示意圖,利用三角網群建立之直 線偵測工作區域如圖 8(a)。Sobel edge detection 之 梯度強度影像如圖 8(b)。細化及線追蹤成果如圖 8(c)。直線偵測之成果如圖 8(d)。
(P1 2 P2 P3) (P7 2 P8 P9) (P3 2 P6 P9) (P1 2 P4 P7)
G= + × + − + × + + + × + − + × +
(a) (b)
(c) (d) 圖 8. 屋緣直線偵測示意圖
(a)直線偵測工作區域,(b)梯度強度影像,(c)細化線追蹤成果,(d)直線偵測成果
3.4 候選屋頂面建立
此步驟之目的,為建立屋頂面之封閉多邊形。
由前步驟偵測出之屋緣直線線段,利用線與線之關 係將房屋區域之影像分割成數個多邊形。將封閉之 多邊形保留後,視為候選房屋屋頂面。
影像空間中之候選屋頂面,為二維影像上之多 邊形所構成,其中並無高度資訊。直線偵測過程容 易偵測出不屬於屋緣之線段,且在影像上進行切割 時,是將線段構成之所有多邊形保留下來。因此在 建立候選屋頂面時會出現不屬於屋頂面之多邊 形,必須在下個步驟中加入光達資訊判別候選屋頂 面之正確性並且計算候選屋頂面之高度。
3.5 候選屋頂面判別與高度計算
此步驟之目的,是對於候選屋頂面之正確性進 行判別,以及計算屋頂面於物空間之高度。建立出 候選屋頂面後,假設候選屋頂面為具有單一高度之 水平面。依照此假設,加入房屋區塊之前處理後位 於屋頂平面之光達點,以每個候選屋頂面為單位,
判斷其正確性,將不正確之候選屋頂面刪除。並利
用光達之三度空間資訊計算出屋頂面之高度後,將 高度相近之候選屋頂面合併。
3.5.1 候選屋頂面之判別
此步驟將前處理過後位於屋頂平面之光達 點,投影至影像空間。利用影像上之光達點對於候 選房屋屋頂面進行判別,刪除不屬於房屋之屋頂 面。
將屬於屋頂平面之光達點投影至影像後,理論 上影像空間中錯誤的候選屋頂面內是不會有光達 點存在的。因前處理過程不完美、直線偵測以及光 達點位之誤差等因素,導致錯誤的房屋屋頂面內仍 然會有些許光達點的存在,但相對於正確的屋頂面 內之光達點數,顯然少了許多。
假設光達點於物空間之分佈是均勻的,在固定 面積之內之光達點數會相當接近。在前處理過程中 會將不屬於屋頂平面之點濾除,將前處理過後之光 達點投影至影像後,在錯誤的房屋屋頂面中固定面 積之光達點數便會減少,在正確的候選屋頂面中固 定面積之光達點數就會相對較高。
本步驟利用原始光達點雲分佈均勻之特性,將
經過前處理後之光達點投影至影像上,計算出候選 屋頂面內之光達點與該屋頂面面積之比值,利用此 比值為約制,對於候選屋頂面進行判別。此比值必 須大於設定之門檻,將比值過低之候選房屋刪除,
剩餘的屋頂面視為正確之房屋屋頂面。
3.5.2 候選屋頂面之高度計算
此步驟以屋頂面為單位,搭配影像上屋頂面內 之光達點,採用迭代之方法計算出屋頂面在物空間 之高度。計算時於屋頂面內求得光達點高度之平均 值以及標準差,若通過收斂門檻,則高度平均值為 屋頂面高度。反之則依照標準偏差將信賴區間以外
之光達點濾除,並重新計算高度以及標準差,直到 通過收斂門檻。
進行上述迭代計算後,可獲得每個屋頂面於物 空間之高度。可依照每個屋頂面之高度資訊與相鄰 近之屋頂面合併,並將合併屋頂面之高度平均,為 合併後之屋頂面高度。屋頂面合併之條件如下:
(1) 屋頂面必須相鄰。
(2) 相鄰屋頂面之高度足夠接近。
屋頂面建立示意如圖 9,利用直線偵測成果建 立之候選屋頂面如圖 9(a),圖 9(b)為投影至影像空 間上之平面之光達點。經過判別後正確之屋頂面如 圖 9(c)。合併後之屋頂面如圖 9(d)所示。
3.6 房屋模型重建
建構房屋模型,是以面為單位,將影像空間中 房屋屋頂面多邊形之節點影像座標,搭配屋頂面之 高 度 , 便 可 透 過 共 線 條 件 式 , 以 光 線 追 蹤 法 (Ray-Tracing),將影像中二維的平面投影回物空間
產生三維屋頂面,三維屋頂面示意圖如圖 10(a)。
將影像空間之二維平面投影至地面三度空間 時,因為屋頂面間高度不同,以致於投影回物空間 之三維屋頂面會有高差移位之情形產生。此高差移 位使得原本於影像空間中相連之屋頂面,投影回物 空間中之平面位置會有偏移之情形發生。此偏差可
(a) (b)
(c) (d) 圖 9. 屋頂面建立示意圖
(a)候選屋頂面,(b)影像上光達資訊,(c)判別後之屋頂面,(d)合併後之屋頂面
用 SMS 法修正,並同時組合三維屋頂面,成為最 終之房屋模型。房屋模型如圖 10(b)。
4. 研究成果討論
研究成果討論中可分為房屋區塊偵測與以及 房屋模型重建兩部份。房屋區塊偵測之成果,與一 千分之一向量圖進行比對分析。房屋模型重建之成 果,與人工量測空照立體對且重建完成之房屋模型 為檢核用之模型進行成果分析。
4.1 測試資料簡介
研究中之測試資料位於新竹科學園區,包含光 達資料以及高解析力之數位化彩色空照影像。其中 光達資料為農委會所提供,為 2002 年 4 月間由 Leica 公司之 ALS 40 之系統所收集而得。光達原 始點雲解析力大約為 1.6 point/m2,高程精度大約 15cm,平面精度大約 30cm,且包含經過分類之地 表點與地面點。空照影像為中華顧問工程司所提 供,拍攝日期為 2002 年 3 月 20 日,採用大比例尺 之數位化彩色影像,為 10 公分解析力之影像,經 空中三角測量評估,平面精度大約 21 公分。相關 資料簡介如表 1。
4.2 房屋區塊偵測成果評估
在房屋區塊偵測之部分,挑選兩測試區,分別
包含光達資料以及高解析力之數位化彩色空照影 像。測試區一以及測試區二之大小分別是 400m × 320m 以及 260m × 400m。兩測試區之局部放大空 照影像如圖 11。
進行房屋區塊偵測,分別將兩測試區之光達資 料中已分類完成之地表點與地面點,網格化成一公 尺解析力之 DSM 與 DTM。測區一之 DSM 如圖 12(a),測區二之 DSM 如圖 12(b)。測區一之 DTM 如圖 13(a),測區二之 DTM 如圖 13(b)。空照影像 中進行監督式分類,將影像中樹木部分萃取出,分 類影像如圖 14。偵測出之房屋區塊如圖 15。檢核 用之向量圖網格化後如圖 16。
房屋區塊成果以房屋棟數為單位進行評估,測 試區一之檢核用向量圖中共有房屋 40 棟,本研究 偵測出正確之房屋區塊共 32 棟,遺漏 8 棟房屋。
測試區二中檢核用向量圖中共有房屋 18 棟,本研 究偵測出正確之房屋區塊共 14 棟,遺漏 4 棟房屋。
以房屋棟數為單位進行評估,成果如表 2.。
房屋區塊成果評估時以網格為單位,在測試區 一之房屋偵測成功率為 77.6%。測試區二房屋偵測 成功率為 81.2%。測試區一之偵測成果如表 3。
表 1. 測試資料簡介 系統名稱 Leica ALS
40
航照形 式
數位化彩 色影像 航高 800m FOV 57°
FOV 35° 航高 1543.43m 雷射波長 1.064 μm 焦距 305.11mm 掃描頻率 29.4 Hz 空間解
析力
10cm
雷射脈衝 率
38 KHz 平面精 度
約 0.211m 水平精度 約 30cm
高程精度 約 15cm
(a) (b) 圖 10. 房屋模型重建示意圖
(a)物空間三維屋頂面,(b)SMS 法組合完成之房屋模型
(a) (b) 圖 11. 測試區局部空照影像
測試區一,(b)測試區二
(a) (b)
圖 12. 測試區 DSM (a)測試區一,(b)測試區二
(a) (b)
圖 13. 測試區 DTM (a)測試區一,(b)測試區二
(a) (b) 圖 15. 表面分析後最終房屋區塊
(a)測試區一,(b)測試區二
(a) (b)
圖 14. 分類影像(白色:樹木區域、黑色:非樹木區域) (a)測試區一,(b)測試區二
(a) (b)
圖 16. 一千分之一向量圖網格化 (a)測試區一,(b)測試區二
表 2. 房屋區塊偵測成果分析(以棟數為單位) 單位:棟 測試區一 測試區二 共有房屋 40 18 偵測房屋數 33 14 偵測正確 32 14
偵測錯誤 1 0
遺漏 8 4
偵測成功率(%) 80.0% 77.8%
表 3. 房屋區塊偵測成果分析(以網格為單位) 單位:網格 測試區一 測試區二
共有房屋 30565 11735 偵測房屋數 27444 16392 偵測正確 23729 14401 偵測錯誤 3715 1991
遺漏 6836 3334 偵測成功率(%) 77.6% 81.2%
4.3 房屋模型重建成果
模型重建之成果部分,對於房屋外型輪廓較為 簡單,且符合研究設定處理範圍之地區,即測試區 一,進行模型重建正確率以及精度評估。對於房屋 外型較複雜之且不於研究所設定範圍內之地區,即 測試區二,則因模型重建能力尚不足則未加以評 估。
在房屋外型輪廓較為簡單部分,即測試區一中 之部分地區,共偵測出 21 個房屋區塊。其中 19 個房屋區塊之房屋符合本研究中設置之處理範 圍,有兩棟不符合,針對 19 棟符合處理範圍之房 屋進行重建成功率以及模型精度評估。
4.4.1 測試區一房屋模型重建成果
測試區一房屋模型重建之成果,在偵測出之 19 個房屋區塊中,依本研究提出之方法,可重建
出 18 棟房屋模型。測試區一重建出之房屋模型成 果如圖 17(a),經由人工量測且重建完成之模型如 圖 17(b)。
4.4.2 測試區一房屋模型正確率評估
房屋模型正確率評估方面,針對測試區一內之
重建成果,利用人工於空照立體對中量測且重建完 成之房屋模型作為檢核,進行重建成功率以及模型 精度之評估。在重建成功判斷之準則方面,以每棟房屋模型 為單位進行評估。若重建出之模型輪廓與檢核用之 模型偏差經由人工量測在 2 公尺以上,則視為重建 失敗。重建正確率分析結果,在測試區內共有 19 個房屋區塊,由本研究之方法可重建出之 18 棟房 屋模型。由檢核用之模型評估之結果,共有 15 棟 房屋重建成功,重建成功率為 79%。在屋頂突出物 方面,檢核用之模型中共有 45 個屋頂突出物,本 研究之方法共重建出 24 個,其中 23 個重建成功,
重建正確率為 51%。房屋模型重建正確率評估如表 4。
4.4.3 測試區一房屋模型精度評估 與討論
在房屋模型精度評估,針對測試區一內重建正 確之房屋模型中,挑選出 102 個角點與檢核用之模 型進行分析。挑選時模型之角點必須在檢核用之模 型中為明確之點(如角點於檢核用之模型中為圓弧 等則不列入精度評估)。挑選出之 102 個點經分析 後,其均方根誤差在 X、Y、Z 方向分別為為 0.53、
0.47、0.43 公尺,精度評估如表 5。平面誤差向量 圖如圖 18(a),高程誤差向量圖模圖 18(b)。
表 4. 測試區一房屋模型重建正確率評估
單位:個數 房屋 屋頂突出物
總數(Total) 19 45
重建數 18 24
重建成功 15 23
重建失敗 3 1
遺漏 4 22
正確率 79% 51%
表 5. 測試區一房屋模型精度評估
Mean_X(m) 0.24 Rmse_X(m) 0.53 Mean_Y(m) 0.21 Rmse_Y(m) 0.47 Mean Z(m) 0.23 Rmse Z(m) 0.43
(a) (b)
圖 18. 房屋模型誤差向量圖 (a)平面誤差,(b)高程誤差 依誤差向量圖所示,平面方向誤差有局部之系
統誤差存在,產生之地方主要為屋頂突出物。此誤 差產生之主要原因,是因為屋頂突出物之邊界特徵 不明顯,使該邊緣被附近其他較明顯之特徵線段取
代。如圖 19(a)所示,箭頭所指之突出物之邊界,
被兩旁房屋主要屋頂面之邊緣輪廓線段所取代,使 得建立屋頂突出物之多邊形與實際之邊緣會有一 段空隙。在屋頂面計算得高度且投影回物空間後,
(a) (b)
圖 17. 測試區一房屋模型重建成果 (a)重建之模型,(b)檢核用之房屋模型
利用 SMS 法組合出之房屋模型如圖 19(b)。在箭頭 所指之屋頂突出物與房屋之邊緣是緊連在一起,而 檢核用之房屋模型中之突出物是與房屋之邊緣有
一定之距離,如圖 19(c)箭頭所指處,因此在屋頂 突出物處會有局部之系統誤差存在。
(a) (b) (c) 圖 19. 局部系統誤差產生示意圖
(a)重建之房屋屋頂面,(b)本文重建之模型,(c)檢核用之房屋模型 在實驗流程中直線偵測之門檻設定對於模型
之成果具有一定之敏感度,舉兩例說明。如約制線 段與線段間之平行或正交之角度差在於 4 度內,使 得屋緣線段未能偵測出,如圖 20(a)箭頭所指處。
若於此例子中將角度差門檻調整為 5 度,則遺漏之 屋緣線段便可在此放寬之門檻下被偵測出,如圖 20(b)箭頭所指處。另外在直線偵測之約制條件為
線段與線段間必須相距 30 個像元以上,使得圖 20(c)箭頭所指處之屋緣線段在此限制之下未能被 偵測出。若將此例中直線偵測之線段與線段間之約 制門檻調整為 20 個像元,則未能偵測出之屋緣線 段便能在此約制條件下被偵測出,如圖 20(d)箭頭 所指處。
(a) (b)
(c) (d)
圖 20. 門檻值敏感度示意圖
(a)角度門檻 4 度之直線偵測成果,(b)角度門檻 5 度之直線偵測成果,
(c)距離門檻 30 像元之直線偵測成果,(d)距離門檻 20 像元之直線偵測成果
4.4.4 測試區二房屋模型成果與討論
房屋外型較複雜之測試區二模型重建成果,在 偵測出之 14 個房屋區塊中,可重建出 11 棟房屋模
型。測試區二依本研究之方法重建出之模型成果如 圖 21(a),經由人工量測且重建完成之模型如圖 21(b)。
(a) (b)
圖 21. 測試區二房屋模型重建成果 (a)本研究重建之房屋模型,(b)檢核用之房屋模型
在此測區中,房屋之大致輪廓可重建出。本研 究對房屋輪廓之限制,為必須由直線線段所構成。
在對房屋輪廓特徵萃取時,是對於直線特徵進行偵 測。因此針對外型複雜之房屋,如房屋輪廓非直線
所構成,則無法偵測出,抑或被直線所取代。如圖 22(a)標示處為不具直線特徵之房屋屋緣,在直線 偵測時無法偵測出該屋緣線。重建完成之模型如圖 22(b),標示處為重建錯誤之處。
(a) (b)
圖 22. 複雜房屋之模型重建
(a)複雜房屋之空照影響,(b)本研究重建之房屋模型
5. 結論與建議
本研究以資訊融合為概念,提出完整之程序,
結合光達資料以及數位空照影像進行房屋區塊偵 測以及房屋模型重建。以研究結果來看,本法具有
自動化之潛力。
房屋區塊偵測之成果,運用房屋在光達資料中 具有之高度及形狀外觀等特性,以及空照影像之色 彩資訊,將房屋區塊偵測出。房屋區塊偵測真成果 部分對於樹木所造成之影響尚無法完全克服。如分
類錯誤造成之影響,使房屋區塊被錯誤濾除。在房 屋與樹木相連結之部分,若樹木之色彩資訊不足,
則無法將房屋與樹木區分開,使房屋區塊邊緣有許 多不平整之情形發生。
房屋模型重建之成果,以偵測出之房屋區塊出 發,搭配空照影像平面精度佳以及邊界萃取易,以 及光達資料高程精度佳之優勢進行模型重建。房屋 模型重建之部分,對於房屋外型做了諸多限制,如 屋緣需正交以及屋頂為平頂等特性。在此限制條件 之下,則於外觀輪廓較複雜之房屋無法完整重建,
但偵測出之房屋區塊仍有可用之處。在未來可依區 塊偵測成果為出發點,嘗試不同之方法於複雜房屋 之重建。
在建議方面,可加入其他輔助資料以獲得更完 整之資訊。如加入多光譜影像,以近紅外波段對於 植物區分容易之特性更完整的濾除樹木區域,建立 更完整之房屋區塊。在影像特徵萃取部分,加入空 照影像立體對,可減少遮蔽之情形發生。對於偵測 出之屋緣線段於兩張影像上進行線段特徵匹配,以 獲得更可靠之屋緣線段。
模型重建中直線偵測之成果決定房屋模型之 外觀與輪廓,若影像上房屋之屋緣不明銳,易導致 模型重建失敗。在直線偵測之門檻設定對於模型之 成果具有一定之敏感度,關於門檻給定之部分,需 以經驗法則獲得較佳之成果。如何降低門檻之敏感 度,有待進一步之研究。
在方法上之建議,在影像空間中除了偵測直線 特徵外,可加入圓弧偵測,對於建物之輪廓描述更 為完備。在光達資料中可組合三角網後,以區塊增 長(region growing)之方式萃取出完整之平面,並加 入斜面分析,對於面之判別更為完備。
本研究所產生之房屋模型,可視為初始之房屋 模型。在後續之工作,對於此初始房屋模型更進一 步琢磨,以獲得更精確之房屋模型。
致謝
本研究承蒙行政院國家科學委員會研究計畫 (NSC92-2211-E008-052)支持得以順利完成。另感 謝行政院農委會農林航空測量所提供光達資料,以 及中華顧問工程司提供大比例尺空照影像。
參考文獻
邵怡誠、陳良健,「利用光達資料於 DTM 生產及 房屋偵測」,第二十二屆測量學術及應用研討 會論文集,第 87-94 頁,2003。
饒見有、陳良健,「局部最大值追蹤法於邊緣影像 細化及向量化之研究」,八十六年電子計算機 於土木水利工程應用論文研討會論文集,第 471-476 頁,1997。
Ackermann, F.,1999, “Airborne Laser Scanning – present status and future expectations”, ISPRS Journal of Photogrammetry & Remote Sensing, Vol. 54, pp. 64-67.
Alharthy, A., and Bethel, J., 2002, “Heuristic filtering and 3D feature extraction from LIDAR data”, IAPRS, vol. XXXIII, pp. 29-35, Graz, Austria.
Briese, C., Pfeifer, N., and Dorninger, P., 2002,
“Application of the Robust Interpolation for DTM Determination”, IAPRS, vol. XXXIII, pp.55-61, Graz, Austria.
Chen, L. C., and Rau, J. Y., 2003, “Robust Reconstruction of Building Models from Three-Dimensional Line Segments”, Photogrammetry Engineering & Remote Sensing, Vol. 69, No .2, pp. 181-188
Gruen, A., and Dan, H., 1997, “TOBAGO – a topology builder for the automated generation of
building models”, Automatic Extraction of Man-Made Objects from Aerial and Space Images (II) Monte Vertia, pp. 149-160.
Halla, N., and Brenner, C., 1998, “Interpretation of Urban Surface Models Using 2D Building Information”, Computer Vision and Image Understanding, Vol. 72, No.2, pp.204-214.
Harwood, D., Subbarao, M., Hakalahti, H., and Davis, L. S., 1987, “A new class of edge-preserving smoothing filters”, Pattern Recog. Lett., vol. 6, pp. 155-162, August 1987.
Morgan, M., and Tempfli, L., 2000, “Automatic Building Extraction from Airborne Laser Scanning Data”, IAPRS, vol. XXXIII. Part B3.
pp. 616-623 Amsterdam.
Parker, J. R., 1997, “Algorithms for Image Processing and Computer Vision”, John Wiley
& Sons, Inc, U.S. pp16-19.
Hough, P.V.C., 1962, “Methodsand Means for Recognising Complex Patterns”, U.S. patent No.306954
Richards, J. A., 1986, “Remote sensing digital image analysis”, Springer-Verlag, Berlin, 281 pages.
Rottensteiner, F., and Briese, Ch., 2002. “A New Method For Extraction In Urban Areas From High-Resolution LIDAR Data”, ISPRS, vol.
XXXIII, pp. 295-301, Graz, Austria.
Schenk, T., and Csathó, B., 2002, “Fusion of LIDAR Data and Aerial Imagery for a More Complete Surface Description”, IAPRS, vol. XXXIII, pp.
310-317, Graz, Austria.
Shepard, D., 1968, “A two-dimensional interpolation function for irregularly-spaced data”, Proc. 23rd National Conference ACM, ACM, 517-524.
Wehr, A., and Lohr, U.,1999, “Airborne Laser Scanning – An Introduction and Overview”, ISPRS Journal of Photogrammetry & Remote Sensing, Vol. 54, pp. 68-82.
Integration of LIDAR Point Clouds and Aerial Imagery for Building Reconstruction
Yen-Chung Lai
1Liang-Chien Chen
2Jiann-Yeou Rau
3ABSTRACT
From the viewpoint of data fusion, we integrate LIDAR point clouds and digital aerial image to perform 3D building modeling in this study. The proposed scheme comprises two major parts: (1) building block detection and (2) building model reconstruction.
In the first step, height differences from LIDAR point clouds are analyzed to detect the above ground areas. Color analysis is then performed for the exclusion of tree areas. Remaining regions are restricted by area threshold. After the refinement by shape and texture analysis, building blocks are established.
In the second step, we analyze the height information of LIDAR point clouds to estimate the approximate roof-edges and extract the LIDAR points located on the rooftop. Based on the approximate roof-edges, accurate 2D edges are calculated in aerial image and then transformed into candidate roof-planes. Combined with the extracted LIDAR points, the correct roof-planes are selected and projected back to object space as 3D roof-planes. At last, a patented method SMS (Split-Merge-Shape) is employed to generate building models using the 3D roof-planes. LIDAR point clouds and large scale aerial image in Hsin-Chu Science-based Industrial Park were used in the test. A 1/1000 scale topographic map and reconstructed building models measured manually from stereo pairs were used to evaluate the results.
Key Words: LIDAR point clouds, aerial image, building detection, building reconstruction, building model
Received Date: July 08, 2004 Revised Date: Nov. 22, 2004 Accepted Date: Nov. 24, 2004
1 Graduate Student, Department of Civil Engineering, National Central University
2 Professor, Center for Space and Remote Sensing Research and Department of Civil Engineering, National Central University
3 Specialist, Center for Space and Remote Sensing Research, National Central University