Volume 12, No.3, September 2007, pp. 197-216
小波演算法重建點雲覆蓋面
蔡展榮
1嚴晟瑋
2摘要
本文提出一個重建點雲覆蓋面的演算法，它應用具備碎形表達能力的二維三階 Daubechies 小波之尺 度函數來組成每一個點的觀測方程式，運用由粗而細的求解策略，並使用二元格點位置上的虛擬觀測量 PHO 及 POI 來解決最小二乘平差求解過程常出現的劣態問題。同時，本文也提出了一個全自動重新給權 的模式，來降低改正數絕對值大於兩倍先驗高程精度的光達點高程觀測值之權值；空載光達點雲實驗結 果驗證了此法已有效降低「吉布斯效應」之影響，而且當二元格點的點間隔(1.25m)與點雲之平均點間距(約 1m)相當時，後驗單位權中誤差為±20~23cm，已相當於點雲的先驗高程精度±25cm。
關鍵詞：物面重建、空載光達點雲、小波、劣態條件、吉布斯效應
1. 前言
近年來，光達(LIght Detection And Ranging, LIDAR)的資料廣泛地被應用於不同的場合，例如 空載光達產製高解析力數值地形模型、三維城市模 型之萃取與重建等；雖然光達測量具有”快速”取得 大量三維點雲的優勢，可是，目前的光達測量仍有 一些嚴重的缺點，例如：可靠度為零（光達點無多 餘觀測、缺乏自我檢核錯誤的能力）、光達點未落 在物面特徵上、（對高品質應用而言）解析力仍嫌 不足，所以它必須和其它測量方法合作以確保成果 的品質。另外，在實務上，欲重建高精度的點雲覆 蓋面（包含地面、房屋面、植被表面、路面等），
則所使用的物面模式要能夠表達碎形的真實面函 數。
直至今日，表達一個物面的方式可概略地分為 兩 大 類: (1) 以 一 些 簡 單 的 幾 何 元 件 (geometry primitives)組成一個複雜的物面、(2)以歐幾里德幾 何學、非歐幾里德幾何學、或碎形幾何學為基礎的 物面表達法。在第一大類中，最常見的例子有: 以
三角面或規則網格面或立方體元(voxel)來組成物 面、CSG(Constructive Solid Geometry)模型法等；
其中 CSG 模型法雖可藉由高階曲面(high-order surface primitive)之方式來描述較複雜的物面(You et al., 2003)，然而它所能表達的複雜度仍會受限於 所使用的幾何元件。相對地，在第二大類中，由於 物面是由一些所給定的參數以及特定的數學模式 產生，此類型的物面即可具有較高的複雜度與不連 續度(Hill, 2000)。如：在歐幾里德幾何學中，有 Bezier 曲面、B-spline 曲面等表達方式(Busé and Galligo, 2004; Xie et al., 2003; Anca et al., 2004)；在 非歐幾里德幾何學中，可藉由非正交的歐幾里德空 間來表達物面(Hill, 2000)；在碎形幾何學中，物面 可由所給定的隨機模式與碎形維度來產生(Peitgen and Saupe, 1988)。
然而，上述兩大類之物面表達法於實務應用上 仍有一些限制與問題。實務上，現有的物面表達法 僅能夠描述（局部）較平滑、連續、且不具碎形特 徵之物面，使得依此等表達法所設計之物面重建模 式，無法精確地重建物面的局部細節特徵。例如:
1國立成功大學測量及空間資訊學系副教授
2國立成功大學測量及空間資訊學系碩士
收到日期:民國 95 年 09 月 25 日 修改日期:民國 96 年 08 月 29 日 接受日期:民國 96 年 08 月 30 日
Hoppe et al.(1992)提出了一個整合型的物面重建演 算 法 ， 應 用 等 值 面(iso-surface) 及 Lorensen and Cline(1987)提出的移動方塊法之原理，將未知的物 面重建出。Park et al.(1999, 2000)以不規則三角面 (Triangulated Irregular Network, TIN)之方式重建初 始模型，進而將重建模型進行分割，將分割後各小 塊之物面模型以 NURBS (Non-Uniform Rational B-Spline)來模擬，得到一個整體的 NURBS 物面模 型。Carr et al.(2003)運用 Radial Basis Function(RBF) 演算法和點雲來重建出平滑趨勢物面，它能夠對不 規則、不平滑、或具有破洞的點雲進行內插計算，
將內插所得的平滑變化之點群，以 TIN 組成重建 模 型 。Tsay(2000, 2002) 提 出 一 個 以 小 波 理 論 (Wavelet Theory)為基礎之物面重建演算架構，運用 最小二乘法，進行不規則的離散取樣點與近似函數 面(由 Haar 小波函數所構成)間的套合計算，重建多 解析力物面。
檢視現今的物面重建法可知，它們受限於所能 夠使用的物面表達法，導致無法表達局部物面之不 連續、不平滑、以及碎形的特徵，使得重建的結果 往往僅只是較平滑且連續的物面趨勢而已。此外，
一些演算法並非基於測量領域之應用面來設計，因 此點雲所包含的各式誤差往往是被忽略的，使得各 重建模式缺乏可供品管使用的後驗精度指標。因 此，本文欲將現有的小波訊號重建法(Tsay, 2000;
Tsay, 2002)予以推廣，使用具碎形表達能力之三階 Daubechies 尺度函數，提出一個能夠表達連續的與 不連續的、平滑的與不平滑的、規則的與不規則 的、以及碎形的與非碎形的三維物面之重建演算 法，且將點雲之偶然誤差特性考量於此模式中，以 克服現有的物面重建法之缺點和限制，使得重建後 的物面得以具備一些可供參考的後驗精度指標，並 以此做為真實的三維複雜物面重建之初步研究。
2. 小波表達之面函數
令 f(x,y)是一個有限能量的(square integrable) 連 續單值 實數 函數(continuous single-valued real function)，根據小波理論，可以使用(1)式來表達面
函數 f 的近似函數Ajf( yx, )。
] 1 , 0 [ , ] 1 , 0 [ , ) , ( ) , ( 2( 1)1
1 1 ) 1 ( 2
1 ,⋅ ,,, ∀∈ − ∀∈ −
= ∑−− ∑
−
=
−
−
−
=
N y M x y x c y
x f A
M
k N
l
l k j p j
l k j
j j
φ (1)
) 2 ( 2 ) 2 ( 2 ) ,
,(
,
,jkl x y j p jx k j p jy l
p = φ − φ −
φ (2)
∫ ∫ ( )
∞
∞
−
∞
∞
−
⋅
=
= f f x y x y dxdy
ckj,l ,φp,j,k,l ( , ) φp,j,k,l , (3)
其中，j 是尺度函數的尺度參數，∀j∈Z。k、l 代 表尺度函數的空間平移參數，∀ ,kl∈Z。φp,j,k,l代表 p 階的二維 Daubechies 尺度函數φp( yx, )經參數 j 的 尺度縮放及參數 k、l 的平移調整後之基底函數，
由(2)式計算得到。ckj,l代表兩實數函數 f 與φp,j,k,l的 內積值，稱為尺度係數(scaling coefficient)，由(3) 式計算得到。同時，我們將所欲重建之訊號函數 f 之計算區範圍令為[0, M-1]x[0, N-1]，於此條件下，
根據Daubechies 尺度函數的承載區(support)可推導 得到 k、l 值的累加範圍分別為[2p-1, 2j(M-1)]、[2p-1, 2j(N-1)]；而又為了避免尺度函數之承載區落在計 算區之外圍，本文將 k、l 值的累加範圍分別定為[-1, 2j(M-1)-1]、[-1, 2j(N-1)-1] (嚴晟瑋，2006)。
由(1)式可知，若各解析力空間下之各項尺度 係數可求得，則可得到函數 f 於各解析力下之近似 函數。本文考量到所欲重建的覆蓋面之特性(如：
連續度、平滑度、不規則性、碎形特性等)與計算 量(階數 p 越大，尺度函數值的計算量越大)，因此 選擇以具碎形描述能力(Daubechies, 1992; Kaiser, 1994)之三階 Daubechies 尺度函數，做為後續物面 重建之基底函數。
3. 小波演算法重建覆蓋面
3.1 基本假設
一般的傳統測量會施測重要的地面及物面特 徵點，相較之下，光達掃描像是打散彈槍一樣（王 蜀嘉、等，2006），（幾乎完全）沒有光達點落在重 要的地面和物面特徵上。在這種取樣方式下，重建 原始物面的局部細節特徵的必要資料無法直接被 掃瞄記錄下來，因此，本文重建點雲覆蓋面的小波
法有以下的基本假設。
基本假設一：令未包含於點雲的隱性特徵集合中的 原物面細節分量是可被忽略的。重建的覆蓋面應產 生良好套合的曲面，並且應避免在鄰近點之間產生 巨大擺盪的曲面。其中，所謂的「點雲的隱性特徵」
是 指 從 局 部 點 雲 推 導 出 的 物 面 特 徵(derived features)，譬如李姝儀(2005)使用的物面角特徵就 是隱性特徵的一個例子。
基本假設二：令欲重建的覆蓋面是一個單值函數 面。
3.2 小波法表達單值函數面
由於二維的尺度函數是連續的單値函數，因此 (1)式僅能夠用來描述一個單值函數面，無法描述 多值函數面。然而，如同一般常見的物面表達法，
如:Bezier、B-spline 曲面等，可使用面貼片(surface patch)之方式，將許多的單值函數面組合銜接在一 起，來產生真實的三維物面模型。因此，本文僅將 研究重心放在以小波重建單值函數面的課題上。
將點雲的參考坐標系統之 xy 平面線性投影到 二維尺度函數的定義域，使得散佈於此空間中之點 雲為單値連續函數 z(x, y)之取樣點，以滿足 3.1 節 所述之基本假設二。
3.3 函數面套合於點雲上
令(1)式中之Ajf( yx, )與ckj,l分別是離散的點 雲構成的面函數及未知數向量元素，則可得到每一 點的觀測方程式，如(4)式：
∑−− ∑
−
=
−
−
−
=
⋅
=
=
+ 2( 1)1
1 1 ) 1 ( 2
1 , , ,
, ( , )
)
, ( )
,
( M
k N
l
i i l k j p j
l k i
i j i i i i
j j
y x c y
x f A v y x
z φ (4)
其中，
(
xi,
yi,
zi)
代表第 i 個光達點的三維坐標，vi代表高程値zi之改正數。將全部點的觀測方程 式(4)式以矩陣表達，即可得如下的線性方程式系 統：
AX V
L+ = (5)
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎣
⎡
× =
+ n n
z z z
L M
1 0
1 ) 1
( ;
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎣
⎡
× =
+ n n
v v v
V M
1 0
1 ) 1
( ;
⎥⎥
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎢⎢
⎣
⎡
×=
j le ke
j le ks
j ls ks
u
c c c X
, , ,
1
M M
(6)
⎥⎥
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢⎢
⎢
⎣
⎡
× =
+
) , ( )
, ( )
, (
) , ( )
, ( )
, (
) , ( )
, ( )
, (
, , , ,
, , ,
, ,
1 1 , , , 1 1 , , , 1 1 , , ,
0 0 , , , 0 0 , , , 0 0 , , ,
) 1 (
n n le ke j p n n le ks j p n n ls ks j p
le ke j p le
ks j p ls
ks j p
le ke j p le
ks j p ls
ks j p
u n
y x y
x y
x
y x y
x y
x
y x y
x y
x
A
φ φ
φ
φ φ
φ
φ φ
φ
L L
M L L
L L
(7)
) 1 ( ) 1 (
1 ) 1 ( 2
; 1 ) 1 ( 2
; 1
+
−
× +
−
=
−
−
=
−
−
=
−
=
=
ls le ks ke u
N le M
ke ls
ks j j
(8) 其中，由於尺度函數值可使用 Strang 法(蔡展榮，
2004)計算得到，故係數矩陣 A 中之各元素皆可以 數值計算之方式得到。因此，當選定了某一解析力 階層 j 後，即可根據最小二乘平差法，將某一空間 涵蓋範圍內的每一個點的觀測方程式係數，逐筆乘 積累加組成法矩陣，進而求解每個未知數(即尺度 係數ckj,l)。當所有的未知數皆求解得到時，即可利 用(1)式將該解析力階層 j 下之近似函數面繪出，得 到點雲覆蓋之物面。
3.4 劣態問題
在前述之重建模式下，由於(4)式可藉由不同 的尺度參數值 j 來改變尺度函數之承載區大小，隨 著 j 值提高，尺度函數之承載區即愈小，愈能重建 出局部點雲描述的細節面特徵。因此可以很容易地 得到在不同的解析力下之點雲覆蓋面。
圖 1、點雲(圓點)與尺度函數之承載區(實線)的相對位置例
Case A
Case B
Case C
Case D
在實際的應用上，由於點雲的點位分布不是等 間距的，當某一個尺度函數的承載區近似等於(圖 1 的Case A)或小於(圖 1 的 Case B)該區兩相鄰點的 間距時，或當點落在尺度函數的承載區邊緣時(圖 1 的Case C 與 Case D)，則該尺度係數所對應的法矩 陣行(列)向量會等於或近似等於零向量，導致法方 程式系統產生所謂的劣態問題(ill-posed problem)。
為了解決劣態問題，本文採取由粗而細的求解 策略。由粗解析力的計算階層開始求解，直到最細 的解析力階層為止；其中，粗解析力的計算階層沒 有劣態問題，可以得到粗略的覆蓋面。隨著解析力 提高，於產生劣態問題的細解析力階層加入合宜的 虛擬觀測量，以求得該階層的解。
3.5 引入兩種虛擬觀測量
為克服細解析力階層求解覆蓋面經常出現的 劣態問題，本文設計了兩種虛擬觀測量。
a. 虛擬的點雲高程觀測量 PHO
如果某一個粗解析力空間Vj下之各項尺度係 數及其後驗中誤差皆已求定，則由於所得之近似函 數面Ajf
( y
x, )
是一個連續的單值函數面，因此，我們可引入(9)式產生的虛擬的點雲高程觀測量，
使得下一個細解析力階層可以求得最小二乘套合 面 。 以 下 ， 本 文 簡 稱 這 樣 的 虛 擬 觀 測 量 為 PHO(Pseudo Height Observation)，並將它參考的解 析力階層 j 使用下標註記為 PHOj。
∑
− −∑
−
=
−
−
−
=
⋅
=
= 2( 1)1
1
1 ) 1 ( 2
1
, , ,
, ( , )
)
, (
M
k N
l
l k j p j
l k j
j
j j
y x c
y x f A
PHO φ (9) b. 虛擬的點雲高程觀測量 POI
然而，依據小波理論的多解析力分析原理可 知，PHOj將局部的細節面分量忽略不計，所以在 具有局部細節面特徵的地方（例如房屋牆面），
PHOj較不適用。以圖2 為例，其中圓點(○)代表原 始的點雲點位，實線代表於階層 j 重建所得之面剖 線，很明顯的，此時若以該重建面來產生下一階層 計算所需之虛擬觀測量，則所得的PHO 點就是圖 2 中的三角形(△)標示的點位。於原覆蓋面的局部 低頻區（例如平坦地表或路面），PHO 點雲與原始 點雲兩者的覆蓋面是相當一致的，然而於原覆蓋面 的斷線(或斷面)處，PHO 點雲就無法正確反映這些 局部的面特徵。
考量到光達測量之原理，以本文使用的空載光 達系統為例，平均航高約2000m，它測得的光達點 之 平 面 精 度 約 為 2000m×1/2000 1.0m(Optech, 2004)。因此，本文在每一個二元格點(dyadic point) 上，以格點為圓心，距離格點1.0m 半徑的圓形範 圍內，進行最鄰近內插，如圖3 所示，得到二元格 點上的另一群新的虛擬觀測量，本文稱之為 POI (Pseudo Observation by Interpolation)。圖 2 中之方 型(□)即為對原始的點雲逐點進行最鄰近內插所得 之POI 點雲。很明顯地，POI 點雲能夠更有效地保 留原始物面的局部細節特徵。
圖 2、PHO、POI、和原始點雲在高程剖面 xz(或 yz) 平面上的位置圖
圖 3、最鄰近內插 POI 點與原始點雲的平面位置圖
由於最鄰近內插法僅在前述的圓形搜尋範圍 內，找出距離興趣格點最鄰近的一點，令它的高程 值等於該格點的 POI 虛擬高程觀測值；如果在某 一個格點的搜尋範圍內無點（局部點雲分佈較稀疏 之區域），則該格點無POI 虛擬高程觀測值，此時，
本文使用該格點位置上的PHO 值，做為它的虛擬 觀測量。
4. 實驗測試與分析
本文使用中興測量有限公司提供的點雲資 料，它們是由空載Optech ALTM 3070 系統掃瞄台
中交流道鄰近區域而得。表1 說明了 ALTM 3070 雷射掃描儀的一些掃瞄規格參數。圖4 與表 2 顯示 了本文使用的三組測試區之相關資料。
表 1、Optech ALTM 3070 規格表(摘錄自 Optech ALTM 使用手冊)
Optech ALTM 3070 Specification 雷射波長 1064 nm
取樣點之高程精度 ±15 cm (1σ, 1200 m 航高) ±30 cm (1σ, 3000 m 航高) 取樣點之平面精度 ±1/2000×航高 (1σ)
圖 4、三組測試區的空載點雲點位分布圖
表 2、各測試區的點雲分布、衛星影像圖、點數量及測試區大小
測試區 I 測試區 II 測試區 III
點雲分布
測試區之 衛星影像 (影像擷取
自 Google Earth)
資料點數 10921 9578 10047
涵蓋範圍
Δ
X≅ Δ
Y≅
100 mΔ
X≅ Δ
Y≅
100 mΔ
X≅ Δ
Y≅
100 m 高程範圍ΔZ ≅
14.5 mΔZ ≅
28.1 mΔZ ≅
30.5 m測試區 II
測試區 測試區 I
本文的各項實驗都是使用二維的三階不對稱 Daubechies 尺度函數做為物面重建的基底函數，並 以尺度參數j=0 做為起始的解析力階層，同時將每 一個測試區的平面位置範圍（約 100m×100m）線性 投 影 到 該 尺 度 函 數 (j=0) 的 定 義 域 裡 的 區 域 [0, 10]×[0, 10]。
以本文使用的測試資料為例，表 2 所列的各測 試區之平均點間距約為 1m，因此在以下各小節的 實驗中，本文皆以尺度參數j=3 做為重建計算之終 止階層，此時，尺度函數的二元格點間格為 10m/23
=1.25m。
4.1 測試區 II 於階層 0~3 之重建 成果與分析
圖 5a~5d 分別顯示了測試區 II 於階層j=0~3 重 建的覆蓋面成果，其中，圖 5c、5d 之成果分別是 由 22、23分點位置之 POI 結合 PHO1、PHO2做為虛 擬觀測量重建而得；圖 5e、5f 分別顯示了圖 5c、
5d 兩階層 j=2、3 計算時所使用的虛擬觀測量之 POI、PHO 點位。表 3 列出了以上各階層解算成果 的相關數據，其中，
σ ˆ
0為後驗單位權中誤差。表 3、測試區 II 於階層 j=0~3 重建成果之相關數據
階層 二元格點間格(m) 未知數個數 虛擬觀測量總數 POI 點數 PHO 點數
± σ ˆ
0(m)0 10 121 × × × 3.73
1 5 441 × × × 2.64
2 2.5 1681 1681 1645 36 1.97
3 1.25 6561 6561 6439 122 1.50
圖 5、測試區 II 於階層 j=0~3 重建的覆蓋面成果(a)~(d)；階層 j=2、3 計算時所使用之虛擬觀測量 POI、PHO 點位(e)、(f)
在以上之成果中，二元格點位置之 POI 和 PHO 虛擬觀測值以及原始點雲的觀測值都是以彼此等 權之方式加入平差計算中。此時，絕大多數的二元 格點都可以在它們各自的搜尋區裡找到鄰近點，兩 階層j=2、3 分別有 97.86% (=1645/1681)、98.14% (=
6439/6561)的格點有虛擬觀測值 POI，只有很少量的 格點找不到鄰近點而使用 PHO 虛擬觀測值，分別 佔 2.14%、1.86%。
另一方面，隨著尺度函數的解析力提高，重建出 的覆蓋面愈能顯示局部細節特徵；引入二元格點位 置上的虛擬觀測量 POI 或 PHO 後，每一個階層不 僅能得到解，不再出現劣態問題，而且都能得到套 合精度越來越好的覆蓋面，重建而得的覆蓋面於絕 大多數之區域也都能夠滿足基本假設一。然而，本 文也發現，在階層j=3 中，即二元格點間隔(1.25m) 相當於原始點間距(約 1m)的情況下，點雲和覆蓋面 的整體套合精度為±1.5m，顯著大於點雲的先驗高 程精度±25cm；同時，在部分（近似）零階不連續 面的局部區域裡，如斷面或牆面處，重建的覆蓋面 呈現不合理的擺盪面，本文將這一個典型的吉布斯
現象(Gibbs phenomenon)(蔡展榮，1999)之影響效應 稱為「吉布斯效應」(Gibbs effect)。造成小波級數 和傅立葉級數的吉布斯現象之原因已分別在(蔡展 榮，1999)和(Farge et al., 1993)做一分析與說明。
4.2 吉布斯效應的實務解決辦法
本文以圖 5c 的ㄧ個斷面區域的局部放大圖(圖 6) 來扼要說明產生此等吉布斯效應之原因，其中，黑 色點代表原始點雲點位。圖 7 顯示了該斷面之重建 面及點雲間之關係剖面圖。
觀察圖 5c 及圖 5d 可發現，並非所有斷面處皆 會產生上述之吉布斯效應。依據小波級數的吉布斯 現象(Gibbs phenomenon)（蔡展榮，1999）之原理，
只有在原始點雲於(圖 7)斷線區域呈現的高頻局部 細節面超過了該解析力下的小波級數所能夠表達 之詳細程度時，於最小二乘法之計算中，圖 7 中之 理想的重建面將受到該斷線區的點雲影響，使得實 際所得的重建面呈現如圖 7 中灰色實線之趨勢，呈 現不合理的擺盪面。
依據吉布斯現象之原理可知，在斷線、斷點鄰 域處，覆蓋面函數(1)式會有所謂的「逾越誤差 (overshoot value)」（蔡展榮，1999），它將造成此區 的點雲有大改正數值。例如，圖 8 分別顯示了圖 5c 與圖 5d 的計算成果中的改正數絕對値|vi|大於兩
倍先驗的點雲高程精度 (即 )之點位，本文 測試使用的點雲之 =±25cm。由圖 8 可知，幾乎 所有改正數絕對值大於 之點雲皆落在（近似）
零階不連續的斷線鄰域處（例如牆面、樹木區邊界 線），而且解析力越細， 的點雲越集中於覆
:點雲位置 :理想的重建面 :實際的重建面
:尺度函數之 2 元格點位置 圖 7、斷面區域之剖面示意圖 圖 6、斷面區域的局部放大圖
蓋面的局部斷線處，這一個現象也普遍出現於其它 的測試區。
這些測試結果說明了: (1)由於改正數絕對值較大的 點位幾乎全部皆落在斷面(如:房屋之牆面)、獨立物 (如:電力桿、樹木、行道樹)等處，可知此法有很好 的潛力將點雲資料中的隱性特徵點、特徵線、特徵 物偵測出來，詳細的偵測方法尚待後續的研究；(2)
解析力越細，此法能重建越詳細的覆蓋面特徵；(3) 重建的覆蓋面與點雲的套合精度顯著大於點雲的 先驗高程精度(由表 3 可知)；(4)為了抑制、甚至消 除吉布斯效應，同時使得覆蓋面的後驗單位權中誤 差可達與點雲之先驗高程精度相當的等級，本文提 出一個自動重新給權之模式。
圖 8、測試區 II 分別於階層j=2（左）、j=3（右）之重建成果的 的點雲位置
4.3 重新給權之模式
產生這一個重新給權模式的靈感來自以下的 原則：「重新給權後，後驗的單位權中誤差應能達 到與先驗的單位權中誤差相當之等級」，詳細的公 式推導請參閱（嚴晟瑋，2006）第 120 頁。(10)式 說明了每一個改正數絕對值大於 的點位之新 權值的計算公式。
n u n v
p r e
r i
) ( ]
ˆ2 [ 2 −
⋅
⋅
=
σ
(10)其中，r 代表 之點雲總點數，
σ ˆ
i代表 理想的後驗單位權中誤差（即點雲的先驗高程精度），
[
v ]2 er代表由 的r 個點計算所得的改正數平方總和，p 代表需重新給權的 r 個點之新權 值，n 代表點雲觀測量及虛擬觀測量之總點數，u 代表尺度係數未知數之個數。於後續之計算中，此 重新給權之模式皆根據(10)式，以全自動化之方式 由程式來決定該等點位之新權值。
因為 的點雲幾乎全部落在斷面和斷
點鄰域區裡，所以調權處理只調降它們的權值，只 影響斷面和斷點鄰域區的覆蓋面重建結果，對於其 它區的覆蓋面重建結果不會產生顯著之影響，因為 其它區的點雲之 ，所以它們的權值不變，
且 Daubechies 小 波 是 有 限 承 載 函 數 (compactly supported functions)，所以此區域的局部面函數重建 結果無顯著改變。
換言之，重新給權處理只調降 的點雲 高程觀測值之權值， 的點雲之權值以及虛 擬觀測量 POI 及 PHO 之權値仍然維持不變，它們 的權值仍為 1。
4.4 調權抑制吉布斯效應之效果
使用(10)式對 的點雲重新給權後，得 到階層j=2、3 之重建成果，如圖 9a、9b 所示；兩 者的 的點位分別如圖 9c、9d 所示。表 4 顯 示測試區 II 於階層j=2、3 重建計算的一些統計數 據，其中，*代表重新給權後之計算階層；P_O 與 P_R 分別代表不需要及需要重新給權之點群的權
値。
在以上的計算中，虛擬觀測量與 的點 雲高程觀測值均一律以等權 P_O=1 處理之。相較 於圖 5，重新給權後，斷面處的吉布斯效應已明顯 地被有效抑制掉了。同時，圖 9c、9d 也顯示「重
新給權後，|vi|較大的點群亦幾乎皆落在斷面處」。
此外，當 的點雲之權值降低後，如表 4 所 示，於階層j=3 時，重新給權後之後驗單位權中誤 差 ( ±23cm ) 已 相 當 於 先 驗 的 點 雲 高 程 精 度 (±25cm)。
(c) (d)
圖9、重新給權後，測試區 II 重建的覆蓋面：(a)j=2、(b)j=3 ;兩者的 的點位分布圖：(c)j=2、(d)j=3
表 4、測試區 II 於階層j=2、3 重建計算的一些統計數據
階層 未知數個數 虛擬觀測量總數 POI 點數 PHO 點數 P_O P_R
± σ ˆ
0(m)2 1 1 1.97
2* 1681 1681 1645 36
1 0.0076 0.30
3 1 1 1.50
3* 6561 6561 6439 122
1 0.0072 0.23
4.5 測試區 I、III 於階層 0~3 之重 建成果與分析
根據測試區 II 的經驗，本節使用二元格點上 的POI 及 PHO 虛擬觀測量，配合前述的吉布斯效 應之解決辦法，分別對兩測試區I、III 進行各階層 的重建計算，包括等權處理、與重新給權處理，兩
者的重建成果分別顯示並分析如下。
圖10a~10d 分別顯示了測試區 I 於階層 j=0~3 的等權處理之重建成果。其中，圖10c、10d 之成 果分別是由22、23分點位置上的POI 結合 PHO1、 PHO2做為虛擬觀測量重建而得。圖11a、11b 分別 顯示了圖 10c、10d 的計算成果中之 的點
: LIDAR (b) : LIDAR
(a)
位分布圖。
依據(10)式調降 的點雲高程觀測值之 權值，求得的測試區I 於階層 j=2、3 之覆蓋面重 建成果，分別如圖12a、12b 所示。圖 12c、12d 分 別為 12a、12b 兩圖裡的 的點位分布圖。
表5 列出了測試區 I 於階層 j=0~3 重建成果之相關 數據。
觀察圖10c、10d 可發現，相較於測試區 II，
在相同的高程比例尺下，測試區I 的吉布斯效應較 不顯著，究其原因，這是由於測試區I 的點雲高程
值範圍 約為測試區II 者 之一
半，且此區的斷面起伏較緩和之故。另外，由圖 11a、11b 可發現， 的點位也是落在斷面及 獨立物等處，且表5 顯示「等權處理所得的後驗單 位權中誤差±0.57cm 也是顯著大於點雲的先驗高 程精度 =±25cm」，因此，在此測試區中，本文仍 進行了重新給權之處理，以使得在二元格點間格 (1.25m)約等於點雲的平均點間格(1m)時，覆蓋面 重建的後驗單位權中誤差(±20cm)得以反映出點雲 的先驗高程精度( =±25cm)。
圖 10、測試區 I 的覆蓋面重建成果（等權處理）
(a) j=2 (b) j=3
圖 11、在測試區 I 的重建成果中， 的點位分布圖（等權處理）
: LIDAR : LIDAR
: LIDAR : LIDAR
(a) j=0 (b)j=1
(c) j=2 (d) j=3
(c) j=2 (d) j=3
圖 12、測試區 I 重新給權後的覆蓋面重建成果(a)、(b)及其 的點位分布圖(c)、(d)
表 5、測試區 I 於階層 j=0~3 重建成果之相關數據
階層 二元格點間格(m) 未知數個數 POI 點數 PHO 點數 P_O P_R
± σ ˆ
0(m)0 10 121 × × × × 1.47
1 5 441 × × × × 1.11
2 1 1 0.85
2* 2.5 1681 1657 24
1 0.028 0.24
3 1 1 0.57
3* 1.25 6561 6514 47
1 0.032 0.20
圖13a~13d 分別顯示了測試區 III 於階層 0~3 以等權處理而得的重建成果，其中，圖 13c、13d 分別是由22、23分點位置之POI 結合 PHO1、PHO2
做為虛擬觀測量重建所得的覆蓋面。在圖13c、13d
的計算成果中， 的點位分布，分別如圖
14a、14b 所示。
對測試區 III 的 的點雲進行重新給權 後，重建而得的覆蓋面如圖15a、15b 所示，其中，
階層 j=2、3 下的 的點位分布圖，如圖15c、
15d。表 6 列出了以上各階層於解算時之相關數據。
在三個測試區中，測試區III 的點雲高程值分 佈範圍(約為 30.5m)最廣，有高建築物、電力桿、
和高聳分立的路邊行道樹，因此，圖13c、13d 的 部分斷面處，亦可發現吉布斯效應產生的局部不合 理擺盪面。同時，圖 14a、14b 顯示 的點 雲也都落在斷點與斷面的鄰近區域內。
: LIDAR : LIDAR
(a) j=2 (b) j=3
圖 13、測試區 III 於階層 0~3 之重建成果（等權處理）
(a) j=2 (b) j=3
圖 14、測試區 III 之重建成果的 的點位分布圖（等權處理）
為降低吉布斯效應之影響，俾能重建精確的覆 蓋面，並使得覆蓋面重建的後驗單位權中誤差得以 反映點雲的先驗高程精度，本文亦進行了重新給權 之計算。此時，觀察圖 15a 與 15b 所示的重新給權
後之成果可知，斷面處的吉布斯效應皆已完全地被 消除了。另外，重新給權之計算後，部分獨立物（例 如獨立樹）之局部細節面仍能重建良好。
: LIDAR : LIDAR
(a) j=0 (b) j=1
: LIDAR : LIDAR
(c) j=2 (d) j=3
(c) j=2 (d) j=3
圖 15、重新給權後的測試區 III 覆蓋面重建成果(a)、(b)及 的點位分布圖(c)、(d)
表 6、測試區 III 的重建成果之相關數據
階層 二元格點間格(m) 未知數個數 POI 點數 PHO 點數 P_O P_R
± σ ˆ
0(m)0 10 121 × × × × 2.97
1 5 441 × × × × 2.37
2 1 1 2.06
2* 2.5 1681 1529 152
1 0.0062 0.28
3 1 1 1.84
3* 1.25 6561 6010 551
1 0.0053 0.22
4.6 測試區 I~III 重建成果之總結
從前述的三個測試區 I~III 使用空載光達點雲 進行實驗之成果可知，引入二元格點位置上的虛擬 觀測量POI 及 PHO，配合本文提出的重新給權之 模式，依照由粗而細的解算步驟，確實可以有效克 服吉布斯效應之問題，無論是（近似）零階不連續 的斷面、斷點鄰近區域，或是其它的零階連續區
域，都能得到良好的覆蓋面。此外，由於重新給權 之模式僅調整 點雲高程值之權値，因此，
在POI 及 PHO 的輔助下，覆蓋面的局部獨立物（例 如獨立樹、電桿）的面細節特徵仍然可以在重建後 完整的被保留下來。當二元格點間距(1.25m)相當 於點雲的平均點間距(約 1m)時，測試區 I、II、III 求 定 出 的 覆 蓋 面 之 後 驗 單 位 權 中 誤 差 分 別 為
±20cm、±23cm、±22cm，它們反映了零階連續區
: LIDAR : LIDAR
(a) j=2 (b) j=3
域的覆蓋面與點雲的套合精度已相當於點雲的先 驗高程精度±25cm，由圖 9、12、15 可知，（近似）
零階不連續的斷面、斷點區域的吉布斯效應已被有 效抑制，不再產生顯著不合理的擺盪面。
圖16 顯示了三個測試區於階層 j=1~3 重建所 得之的高程影像圖，其中，圖a、f、k 分別是測試 區I、II、III 的點雲高程影像，圖 b~d、圖 g~i、圖 l~n 分別是測試區 I、II、III 於 j=1~3 三個階層求定 的覆蓋面高程影像，圖e、j、o 分別是圖 a~d、圖 f~i、圖 k~n 之高程值所對應的色階圖。在圖 16 中，
圖a、f、k 明顯繪出空載雷射掃描線的軌跡。由階 層 j=1~3 重建的覆蓋面函數(1)式可計算任意格點 的高程值，藉此產生圖16 各階層的高程影像圖。
很明顯地，解析力越細，重建而得的覆蓋面越能夠 呈現更細節的面特徵，而且它和點雲之間的套合精 度越佳。以圖16g~16i 為例，圖 16g 只能顯示出房 屋的粗略結構，此時房屋的牆面區仍呈現較平滑的 變化。隨著解析力的提高，由圖16h 可發現，部分 房屋的牆面已能夠清楚的與地表面分離，並且部分 較細微的結構物亦開始出現在該高程影像中。最 後，當二元格點間格(1.25m)與點雲的平均點間距 (約 1m)相當時（即 j=3），高程影像圖 16i 已能夠 精確重建原始點雲呈現的局部細節特徵，例如圖 16i 的口字型房屋之天井、ㄇ字型房屋之屋頂排氣 孔等。
4.7 目視檢查測試區貼圖產生 的覆蓋面模型
為目視檢查前述的覆蓋面成果之正確性，本文 使用表 2 所示的衛星影像敷貼到重建計算得到的 覆蓋面上，分別得到測試區I~III 於階層 3 之覆蓋 面模型，如圖17 所示。由圖 17 可知，本文提出的 小波法確實能夠重建及表達出連續的(如圖 17 中之 球形及山形屋頂面)、不連續的(如房屋之牆面處、
圖17 中之電力桿)、規則的(如房屋之平坦屋頂面、
道路面、及地表面)、不規則的(如圖 17 的屋頂之 排氣孔)、以及碎形的(如圖 17 的樹叢)覆蓋面。
4.8 小波演算法的自動化程度 及計算複雜度
本文以Borland C++ Builder 6.0 做為小波演算 法的開發平台，配合Open GL 之應用程式介面做 為重建模型之繪圖工具。圖18 顯示了此系統之圖 形操作介面。
在本文開發的重建系統介面中，使用者僅需於 各階層解算時，進行必要的設定，系統就可以全自 動地進行重建計算，例如：1.資料檔輸入/輸出之檔 案路徑設定，2.是否重新給權、3.是否要產生虛擬 觀測量 PHO/POI。故在此設計下，除了前述的系 統初始值設定需人工輔助外，其餘的系統計算處理 是全自動化的。
表 7 列出了此重建系統的測試平台之規格資 料，表8 顯示了各階層之平均解算時間。其中，平 均計算時間包含了: 讀取各項輸入資料(如點雲資 料、虛擬觀測量等)、法矩陣之產生及求逆、未知 數之求解計算、輸出各項重建資料(如 PHO、變方- 協變方矩陣、未知數矩陣、重建模型等)等。
由表8 可知，隨著未知數個數的增加，此重建 系統所需的解算時間即明顯提升。這是由於在此系 統之開發中，本文並未使用任何商用的或共享的矩 陣運算程式集，本文是以 Cholesky 矩陣分解法來 進行線性系統的求解及求逆計算。因此，即使已盡 可能地對此系統進行程式最佳化之處理(包含程式 碼及編譯器之最佳化)，程式的計算時間仍會受到 所使用之矩陣求解及求逆之計算複雜度影響，使得 此系統的運算時間隨著未知數之增加而提高。因 此，若能配合更有效率的法矩陣求解程式集，或結 合現有的商用矩陣處理之程式集，此求解系統之執 行效率必能夠獲得顯著的改善。
5. 結論
本文提出一個小波演算法來重建點雲覆蓋之 物面（包含地面、房屋面、植被表面、路面等），
它應用具備碎形表達能力的二維三階 Daubechies 小波之尺度函數，做為描述覆蓋面的基底函數，並
據以組成每一個點的觀測方程式，進而以最小二乘 平差法來重建點雲覆蓋面。為了克服平差求解過程 常出現的劣態問題，本文運用由粗而細的求解策 略，並使用二元格點位置上的虛擬觀測量 POI 及 PHO；同時，又為了抑制所謂的「吉布斯效應」，
本文也提出了一個全自動重新給權的模式，來降低 改正數絕對值大於兩倍先驗高程精度之(光達)點 高程觀測值的權值。
空載光達點雲實驗結果顯示，此法可以有效抑 制「吉布斯效應」並克服劣態問題，且重建物面的 後驗高程精度已相當於光達點的先驗高程精度。以 本文的三個測試區I、II、III 為例，當二元格點的 點間隔(1.25m)與點雲之平均點間距(約 1m)相當 時，後驗單位權中誤差分別為±20cm、±23cm、
±22cm，已相當於點雲的先驗高程精度±25cm。
(a) (f) (k)
(b) (g) (l)
(c) (h) (m)
(d) (i) (n)
(e) (j) (o) 圖 16、各測試區於階層 j=1~3 重建所得之高程影像圖(詳請參閱本文之說明)
圖 17、不同視角敷貼影像產生的三個測試區於階層 3 重建的覆蓋面模型
圖 18、點雲覆蓋面重建系統之圖形操作介面
表 7、測試平台之規格 Motherboard ASUS P5AD2-E Premium CPU Intel Pentium 4 640 oc. 3.8MHz.
RAM Transcend DDRII 533 1GB oc. 622 MHz.
Hard Disk WD WD1600JD 160GB
VGA ASUS EN6600 GT/TOP/TD oc. 1220MHz.
Others CPU and VGA Watercooling System
表 8、各階層重建解算所需之平均計算時間
階層 未知數個數 平均計算時間 CPU 使用率
0 121 約 0.7 秒 約50 % 1 441 約 1.03 秒 約51 % 2 1681 約 12.03 秒 約52 % 3 6561 約 661.43 秒 約51 %
本文實驗成果驗證了本文提出的小波演算法 之可行性。現有的物面重建法容易受到它們所使用 的物面表達模式之影響，而產生重建後的物面遺失 了局部細節面特徵之結果。相較於此，本文提出的 小波演算法能夠表達連續的與不連續的、平滑的與 不平滑的、規則的與不規則的、以及碎形的與非碎 形的三維物面，將不同性質的面特徵一併考量在 內，使得重建後的光達點雲覆蓋面仍可保留局部細 節面特徵，例如連續的球形和山形屋頂、不連續的 房屋牆面、規則的道路面和地表面、不規則的屋頂 排氣孔特徵、以及具碎形特徵之樹叢表面，同時將 點雲之偶然誤差特性考量於此模式中，以克服現有 的物面重建法之缺點和限制，使得重建後的物面得 以具備一些可供參考之後驗精度指標，並以此做為 真實的三維複雜物面重建之初步研究。
此外，它也能夠在完全沒有任何資料預處理之 情況下，以相當微量的人工輔助，全自動的進行物 面重建計算，並能夠精確重建各式不同的面特徵。
因為產生的覆蓋面函數是單值的面函數，所以此法 可對各階層重建之覆蓋面以常見的規則面網或三 角面網之格式輸出，可供各種常見的三維繪圖平台 展繪多解析力物面模型之用。
由於吉布斯現象造成斷面和斷點區的點雲一 定有顯著不合理的改正數(遠大於這些點的高程觀 測誤差)，如(蔡展榮，1999)的分析與圖 8、9、11、
12 所示，隨著解析力越細，吉布斯現象的影響範 圍越小，大改正數的點雲越集中於斷面和斷點的小 區裡，本文所提的方法能重建越精確的此區覆蓋 面，換言之，它有偵測地物表面上越精確的斷面和 斷點中心位置的應用潛力，是一個尚待研究的課 題。
誌謝
非 常 感 謝 國 科 會 提 供 補 助 經 費 ( 計 畫 編 號 NSC94-2211-E-006-072)、中興測量有限公司提供空 載光達實驗資料，讓本文的研究工作得以順利進 行，兩位審查委員提供非常珍貴的修正意見讓本文 更臻完善，特此表達由衷的萬分謝意。
參考文獻
王蜀嘉、張祖勛，2006。航測數位像機對空載雷射 掃瞄帶來的衝擊，第25屆測量及空間資訊研討
會論文集，2006年9月7~8日，清雲科技大學，
中壢，pp. 179~188。
李姝儀，2005。從地面雷射點雲萃取物面角特徵供 多測站連結之研究，國立成功大學測量及空間 資訊學系碩士論文，台南。
蔡展榮，1999。吉布斯問題的一個實務解法，中國 土木水利工程學刊，第11卷，第3期，民國88 年9月，第175~182頁。
蔡展榮，2004。小波基礎理論與應用，授課講義，
國立成功大學，台南。
嚴晟瑋，2006。以小波為基礎之最小二乘物面重建 演算法，國立成功大學測量及空間資訊學系碩 士論文，台南。
Anca, A., Gaildrat, V., Barthe, L., 2004. Modeling and representing surfaces: Interactive modeling from sketches using spherical implicit functions, Proceedings of the 3rd international conference on Computer Graphics, Virtual Reality, Visualization and Interaction in Africa, pp.
25-34.
Busé, L. and Galligo, A., 2004. Using semi-implicit representation of algebraic surfaces, IEEE Transactions on Shape Modeling Applications 2004, pp. 342-345.
Carr, J. C., Beatson, R. K., McCallum, B. C., Fright, W. R., McLennan, T. J., and Mitchell, T. J., 2003.
Representation: smooth surface reconstruction from noisy range data, Proceedings of the 1st international conference on Computer Graphics and Interactive Techniques in Australasia and South East Asia, pp. 119-126.
Daubechies, I., 1992. Ten lectures on wavelets, SIAM, Philadelphia, PA.
Farge, M., Hunt, J. C. R., Vassilicos, J.C., 1993.
Wavelets, Fractals, and Fourier Transforms, Clarendon Press, Oxford, ISBN 0-19-853647-X, pp.123-142.
Hill, Francis S., 2000. Computer graphics using
OpenGL, 2nd edition, Prentice Hall, ISBN:
0023548568, USA.
Hoppe, H., DeRose, T., Duchamp, T., McDonald, J., and Stuetzle, W., 1992. Surface reconstruction from unorganized points, ACM SIGGRAPH 1992, pp. 71-78.
Kaiser, G., 1994. A Friendly Guide to Wavelets, Birkhäuser Boston.
Lorensen, W. E. and Cline, H. E., 1987. Marching cubes: A high resolution 3D surface construction algorithm, ACM SIGGRAPH 1987, pp.
163-169.
Optech Incorporated, 2004. ALTM 3100 System Specifications, http://www.optech.ca/
pdf/Specs/specs_altm_3100.pdf.
Park, I. K., Yun, I. D., and Lee, S. U., 1999.
Constructing NURBS surface model from scattered and unorganized range data, IEEE Transactions on 3-D Digital Imaging and Modeling 1999, pp. 312-320.
Park, I. K., Yun, I. D., and Lee, S. U., 2000.
Automatic 3-D model synthesis from measured range data, IEEE Transactions on Circuits and Systems for Video Technology 2000, Vol. 10, Issue 2, pp. 293-301.
Peitgen, H. O. and Saupe, D., 1988,. The Science of Fractal Images, Springer-Verlag, New York.
Tsay, Jaan-Rong, 2000. A new algorithm for automatic precise reconstruction of a real object surface, the XIXth Congress of the International Society for Photogrammetry and Remote Sensing (ISPRS), Amsterdam, The Netherlands, 16-23 July 2000, International Archives of Photogrammetry and Remote Sensing (IAPRS), Vol. XXXIII, Supplement B3, pp. 59-66.
Tsay, Jaan-Rong, 2002. A concept and algorithm for 3D city surface modeling, proceedings of the International Society for Photogrammetry and
Remote Sensing (ISPRS) Commission III symposium, Graz, Austria, September 9-13, 2002, “Photogrammetric Computer Vision (PCV’02)”, Vol. 34 (ISSN 1682-1777), Part 3B, pp. B-283-B-286.
Xie, H., Wang, J., Hua, J., Qin, H., and Kaufman, A., 2003. Piecewise C/sup 1/ continuous surface reconstruction of noisy point clouds via local implicit quadric regression, IEEE Transactions on Visualization 2003, pp. 91-98.
You, S., Hu, J., and Fox, P., 2003. Urban site modeling from Lidar, Proceedings of Second International Workshop on Computer Graphics and Geometric Modeling, Vol. 2668, pp.
579-588.
A Wavelet Algorithm for DSM Reconstruction by Using 3D Point Cloud
Jaan-Rong Tsay
1Sheng-Wei Yen
2ABSTRACT
This paper proposes an algorithm for reconstruction of an object surface covered with 3D points.
It utilizes the 2D Daubechies scaling functions of 3rd order, which can describe fractal geometry, to derive the observation equation for each point. The linear system is then solved by the least-squares adjustment (LSA) and the reconstructed surface can then be generated. To overcome the ill-posed problem which often emerges in LSA, we employ a from-coarse-to-fine strategy and use the pseudo observations on dyadic points, called POI (Pseudo Observations by Interpolation) and PHO (Pseudo Height Observations). Moreover, a full-automated weighting model is proposed in order to eliminate the so-called Gibbs effect. It reduces the weights of the points whose absolute residuals are larger than twice the a priori height accuracy of the LIDAR point. Some tests are done by using airborne LIDAR points. They verify that the artifacts caused by the Gibbs phenomenon can be eliminated to the large extent by combining the pseudo observations and the weighting model. While the dyadic points have approximately the point interval of LIDAR points, the a posteriori standard deviations of unit weight in our tests are about ±20~23cm which are all to the extents of the a priori height accuracy, ±25cm.
Key Words: object surface reconstruction, airborne LIDAR point cloud, wavelets, ill-posed
condition, Gibbs effect
1 Associate Prof. Dr.-Ing., Department of Geomatics, National Cheng Kung University
2 M.Sc., Department of Geomatics, National Cheng Kung University
Received Date: Sep. 25, 2006 Revised Date: Aug. 29, 2007 Accepted Date: Aug. 30, 2007
