真實三維碎形物面表達和求定之新方法以及它們在空載和地面雷射掃瞄點雲重建物表三維實景模型之應用(II)

行政院國家科學委員會專題研究計畫 成果報告

真實三維碎形物面表達和求定之新方法以及它們在空載和

地面雷射掃瞄點雲重建物表三維實景模型之應用(II)

研究成果報告(精簡版)

計 畫 類 別 ： 個別型

計 畫 編 號 ： NSC 94-2211-E-006-072-

執 行 期 間 ： 94 年 08 月 01 日至 95 年 10 月 31 日

執 行 單 位 ： 國立成功大學測量及空間資訊學系（所）

計 畫 主 持 人 ： 蔡展榮

計畫參與人員： 碩士班研究生-兼任助理：鄭邦寧、嚴晟瑋、葉佳欣、陳英

煥、劉金燁、陳仕堯、江師榮

報 告 附 件 ： 出席國際會議研究心得報告及發表論文

處 理 方 式 ： 本計畫可公開查詢

中 華 民 國 96 年 01 月 31 日

行政院國家科學委員會補助專題研究計畫成果報告

真實三維碎形物面表達和求定之新方法以及它們在空載和地

面雷射掃瞄點雲重建物表三維實景模型之應用(II)

New Methods for Representing and Determining 3D Fractal

Object Surfaces and their Applications on Reconstruction of 3D

Virtual Reality of Object Surfaces Based on Airborne and

Terrestrial LIDAR Point Clouds (II)

計畫類別：; 個別型計畫 □ 整合型計畫

計畫編號：NSC94－2211－E－006－072－

執行期間：94 年 08 月 01 日至 95 年 10 月 31 日

計畫主持人：蔡展榮

計畫參與人員：鄭邦寧、嚴晟瑋、葉佳欣、陳英煥、劉金燁、陳仕堯、江師榮

成果報告類型(依經費核定清單規定繳交)：;精簡報告 □完整報告

本成果報告包括以下應繳交之附件：

□赴國外出差或研習心得報告一份

□赴大陸地區出差或研習心得報告一份

;出席國際學術會議心得報告及發表之論文各一份

□國際合作研究計畫國外研究報告書一份

處理方式：除產學合作研究計畫、提升產業技術及人才培育研究計畫、列

管計畫及下列情形者外，得立即公開查詢

□涉及專利或其他智慧財產權，□一年□二年後可公開查詢

執行單位：國立成功大學測量及空間資訊學系

中 華 民 國 九 十 六 年 一 月 三 十 日

一、中英文摘要

此研究案使用不同地區的空載光達點雲來測試本研究提出的「加權最小二乘小波演算

法」

，它使用多布吉斯小波來定義一個單值的碎形面函數，運用金字塔演算法來取得虛擬觀

測值，配合一個自行研發的自動化給權模式，施行架橋功能以解決劣態問題，藉以求定真

實三維碎形物面。在全部的實驗區重建出擬合精度±14~23cm 的數值覆蓋面，已和輸入的空

載光達點雲先驗高程精度±15~35cm 相當。同時，本研究也提出一些新構想和新演算法，擬

從離散的光達點群來求定數值地形模型（含地形斷線資料）

，應用吉布斯效應來偵測與重建

房屋、樹木、橋樑等地上物，萃取可靠的物面特徵，它結合了浮測模型法、三維坐標變換、

最小二乘擬合法，進行整體平差，期能萃取與重建可靠的物面特徵供光達測量之各式應用。

關鍵詞：多布吉斯小波、碎形面、光達測量法、吉布斯效應

ABSTRACT

This study uses diverse sets of airborne LiDAR data to validate and evaluate our wavelet-based

least squares algorithm which fits all 3D points on a surface of interest object onto a fractal

surface function defined by Daubechies wavelets. It also utilizes a pyramidal approach to

automatically generate pseudo-observations with proper weighs given by our weighting model.

Thus, a bridging function is utilized to solve the ill-posed problem. Then, the corresponding

digital surface model (DSM) is reconstructed. More test results are analyzed. They concluded that

the DSM has the fitting accuracy of ±14~23cm, which corresponds to their a priori standard

deviations ±15~35 cm of LIDAR points in heights. Furthermore, this study also presents some

novel ideas and algorithms which try, on the one hand, to determine digital terrain models (DTMs)

inclusive of breakline and terrain feature line data, and to detect and reconstruct different objects

such as buildings, trees, bridges and so on. On the other hand, reliable features on object surfaces

are to be extracted and reconstructed by utilizing the floating model approach, 3D coordinate

transformation and least squares fitting method as well, and integrating them into a combined

least square adjustment. The thus derived features should suit for later diverse applications of

LIDAR-grammetry.

Key words: Daubechies Wavelets, Fractal Surfaces, LIDAR-grammetry, Gibbs effect

二、前言

近年來，測量暨空間資訊科技(geomatics)蓬勃發展，繼所謂的 3S[遙感探測 RS(remote

Sensing）、大地測量及衛星定位測量 GPS、空間暨地理資訊系統 GIS]之後，雷達測量

(radargrammetry)、光達(LiDAR，Light Detection And Ranging)測量也跟踵而至，其中，光達

測量是一個典型的直接定位定向(direct georeferencing)的測量法，它可以在短時間獲取被探

測物面的大量密集分佈的高精度三維點雲，是直接獲取物空間資料與資訊的新利器；其中，

空載光達測量以高精度、高解析力、高度自動化且高效率的優勢，已成為世界各國進行大

面積數值地表資料測製的另一個新興主流與趨勢，其多重反射回波的特性，可同時獲取地

面及其覆蓋物(房屋、植被、電力線等)之三維坐標，而透水光達系統更可穿透水體而量測

水底的地形起伏。光達所產製之高精度、高解析力數值高程模型 DEM(Digital Elevation

Model)，可做為土地利用、工程建設規劃、都市計畫管理，河海地形、潮間帶、集水區、

山坡地監測、地理資訊系統、防災、礦業、農業、林業、公共管線等方面數值化、自動化

等應用之基礎[內政部，2005]。

三、國內外相關研究概況簡述

內政部「高精度及高解析度數值地形模型建置計畫」引進空載光達測量技術，分別於

2004 年及 2005 年針對國內中、低海拔及海岸地區各種地形進行測試及分析，以快速提供

地形、地表資訊，提高防災與救災應變機能。同時，內政部也委託工研院辦理「LiDAR 測

區之高精度及高解析度數值地形測繪、資料庫建置與應用推廣工作」案，分別在桃園、新

竹、南投草嶺潭、外傘頂洲、台南縣市、高雄、屏東等地區，針對不同地形及地貌進行測

量、分析及應用推廣等工作，以評估其精度及適用之區域，做為將來應用之依據[王定平、

王成機、陳思仁，2004] 。

在國外，應用空載光達測量法來生產高精度的 DEM、萃取建物、道路、樹木等地表物

供景觀規劃、災害評估等等之應用，已從研究測試階段逐漸進展到實務應用層面[Axelsson,

1999; Haala and Brenner, 1997; Priestnall et al., 2000;Cobby et al., 2001;Steinel et al., 2001]。

四、真實物面之特性

目前為止，使用數值航空攝影測量或空載光達測量等相關的測量科技可以高自動化的

方式來快速獲取任意一個目標物的表面上的三維點雲資料，進而產生真實環境的「三維虛

擬實景地圖(3D virtual reality maps)」

，而真實環境的某一個表面各處的局部表面性質常不相

同，某些地方呈連續的或不連續的、平滑的或凹凸不平的、規則的或不規則的、碎形的或

非碎形的特性，如圖 1 和[Tsay, 2002]所示的各種實例。

圖 1. 三維虛擬實景地圖(3D virtual reality maps)例

圖 2. 表達複雜的真實碎形物面之合宜方法

五、表達真實碎形物面之合宜方法

如圖 2 所示，使用歐幾里德幾何學、或一些簡易的幾何元件、或 Bezier 面、或仿樣面

(spline surfaces) 等 均 無 法 有 效 表 達 精 確 的 真 實 碎 形 物 面 [Ebert, 2003; Musgrave, 1993;

Prusinkiewicz and Hammel, 1993; Deussen et al., 1998; Musgrave et al., 1989; Peitgen and

Saupe, 1988]，而碎形幾何可以表達真實環境裡具碎形性質的物面，所以有學者稱讚碎形幾

何才是正確的數學[Jähne, 1991]；其中，多布吉斯小波能表達碎形幾何[Kaiser, 1994]。

六、小波表達之物面函數

]

1

,

0

[

,

]

1

,

0

[

,

)

,

(

)

,

(

1 ) 1 ( 2 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)

其中，

φ

_p,j,k,l

代表在

V

_j

空間中平移了(k,l)單位後之 p 階多布吉斯尺度函數

φ

_p

(

x

,

y

)

；

代表

對應的尺度係數。

j l k

c

_,

七、觀測方程式

(

)

(

)

(

_{∑ ∑}

−)− ( )

(

− = − − − =

⋅

=

=

+

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

φ

)

(2)

其中，假設高程面函數z(x,y)是一個單值函數(single-valued function)，(x

i

, y

i

, z

i

)是某一個物面

點的三維坐標，未知數是尺度係數

，"k,l。每一個觀測方程式(2)都是線性的，它們組成

的法方程式系統是一個線性系統，所以不需要給定未知數的起始近似值、也不需要迭代計

算，這也是加權最小二乘小波演算法在實務應用上的另一個優點。

j l k

c

_,

八、劣態問題 IPP 及其解決辦法

圖 3. 產生劣態問題之原因

如圖 3 所示，每一個圓點代表某一個物面點，A~D 四例的水平橫線代表相應的多布吉斯尺

度函數

φ

_p,j,k,l

的承載區(support)，其中，B 例將造成該

φ

_p,j,k,l

的尺度係數未知數在法方程式矩

陣對應的行（列）向量是一個零向量，而 A、C、D 三例將造成該

φ

_p,j,k,l

的尺度係數未知數

在法方程式矩陣對應的行（列）向量趨近於一個零向量，它們將造成所謂的劣態問題

(ill-posed problem, IPP)。

為了解決前述的法方程式線性系統求解的劣態問題，本研究案提出 PHO (Pseudo Height

Observation)與 POI (Pseudo Observation derived by Interpolation)兩個演算法，兩者的主要構

想是考量空載光達點的平面精度較高程精度差之特性，配合金字塔由粗而細的演算策略，

使能自動重建物表模型。相對於金字塔的最上面的幾層之重建計算而言，空載光達點雲分

佈很密集，可直接求解得出物面未知數，不需給予任何外加的架橋功能。它們求定的物面

函數可計算下一層格點的高程估值，PHO 法就是使用這些高程估值做為虛擬的格點高程觀

測值，以解決 IPP 問題。而 POI 法是使用在格點一定距離（例如：原始光達平均的點間距）

範圍裡的原始光達點，依據最鄰近法來給定該格點的高程估值，如果無光達點落在該範圍

裡，則由 PHO 在該格點上的高程估值給定之，聯合這些格點的虛擬高程值之觀測方程式以

及原始的光達點觀測方程式，配合適當的給權，以解算物面未知數。

九、演算流程

加權最小二乘小波演算流程，如圖 4 所示。

1. 讀取尺度函數資料和 LiDAR 點雲資料。

式累加到法方程式矩陣中。

權方式累

4.

為零矩陣。

測方程式並將其係數依此權值累加到法方程式矩

7.

架橋功能，則組成每一個虛擬觀測量的觀測方程式並將其係數依等權方式累

8.

未知數之解。

由粗而細的各個解析力階層之最小二乘計算：

2. 組成每一點的觀測方程式並將其係數依等權方

3. 如果需要架橋功能，則組成每一個虛擬觀測量的觀測方程式並將其係數依等

加到法方程式矩陣中。

求解法方程式系統。

5. 重新設定法方程式矩陣

6. 由改正數計算權值，組成每一點的觀

陣中。

如果需要

加到法方程式矩陣中。

求解法方程式系統並輸出

9. 如果需要品質指標，則計算協變方矩陣。

10. 如果需要，可輸出數值高程模型 DEM。

圖 4. 加權最小二乘小波演算流程

十、實驗成果分析

目前已使用中興測量公司提供的台中地區、台南地區、和南一中操場的 Optech

ALT

蓋面上的七條剖線和相應的局部點雲分布圖，圖 9 由上而下

分別

M3070 空載光達資料進行相關的各種測試，全部的實驗區（包括 Tsay and Yen(2006)發

表的成果）重建出擬合精度±14~23cm 的數值覆蓋面，已和輸入的空載光達點雲先驗高程精

度±15~35cm 相當。例如圖 5 顯示南一中操場司令台附近的實驗區影像及其空載光達點雲三

維分布圖，此區也包含樹表面、司令台屋頂面、操場跑道、操場草地區、籃球場、籃球架、

網球場、夜間照明燈，面積為 50m ä 50m，共有 15576 個光達點，平均點間距為 40cm，平

面和高程的先驗精度分別為≤25cm 和≤10cm。圖 6 顯示加權最小二乘小波演算法求定的實

驗區幾何覆蓋面套疊空照彩色影像產生的四種不同視點和視角下的實景模型，而圖 7 顯示

地面攝影取得的實驗區不同角落的影像，目視檢查這些成果圖，顯示加權最小二乘小波演

算法確實能夠在一次計算中將不同性質的各處局部物表面同時求定，且求定的物表面和實

地調查的真實物面有高一致性，其中，由於司令台屋頂面是一傾斜面而且有 3 x 2 個下凹

的矩形區塊，加上所謂的吉布斯效應(Gibbs effect) [蔡展榮，1998、1999]導致求定的此區覆

蓋面產生明顯的彎曲起伏。

圖 8 顯示實驗區求定的覆

是粗解析力到細解析力的第 0~3 階層的大改正數點位分布圖，由左而右依序是

0

ˆ

σ

<|v|§2

σ

ˆ

₀

（紅色）、2

σ

ˆ

₀

<|v|§3

σ

ˆ

₀

（黃色）

、3

σ

ˆ

₀

<|v|（藍色）的點位，它們清楚顯示了

不連續

局部物表面

出現吉

斯效應；同時

大改正數幾乎都出現在零階不連續的

零階

的

區

布

，

局部物表面區，例如房屋的牆面、樹木的邊緣輪廓、籃球架。換言之，它們顯示了應用吉

布斯效應來偵測與重建房屋、樹木、橋樑等地上物之潛力。另外，從這些不同地區的實驗

成果也顯示了加權最小二乘小波演算法從離散的物表覆蓋面的點群來求定數值地形模型

（含地形斷線資料）具有良好的應用潛力。

圖 5. 南一中操場司令台附近的實驗區影像（左）及其空載光達點雲三維分布圖（右）

圖 6. 小波法求定實驗區的幾何面套疊空照影像產生四種不同視點和視角下的實景模型

表 1. 不同解析力階層的小波面和點雲的套合精度

σ

ˆ

₀

、不同等級的改正數之點數和百分比

解析力

統計以及最大改正數 max |v|

階層

0

ˆ

σ

(m)

σ

ˆ

₀

<|v|§2

σ

ˆ

₀

2

σ

ˆ

₀

<|v|§3

σ

ˆ

₀

3

σ

ˆ

₀

<|v|

max |v| (m)

0

2.

23

293

1 (18.8

%)

936 (6.0%)

171 (1.1%)

9.932

1

1.90

1687 (10.8%)

1162 (7.5%)

324 (2.1%)

9.376

2

1.80

1455 (9.3%)

1090 (7.0%)

370 (2.4%)

9.042

3

1.80

1231 (7.9%)

1090 (7.0%)

331 (2.1%)

9.148

3(調權)

0.14

445 (2.9%)

191 (1.2%)

3795 (24.4%)

9.563

圖 8. 實驗區求定的覆蓋面上的七條剖線和相應的局部點雲位置分布圖

圖 9. 由上而下分別是粗解析力到細解析力的第 0~3 階層的大改正數點位分布圖，由左而右

依序是

σ

ˆ

₀

<|v|§2

σ

ˆ

₀

（紅色）

、2

σ

ˆ

₀

<|v|§3

σ

ˆ

₀

（黃色）

、3

σ

ˆ

₀

<|v|（藍色）的點位

十一、結論

實驗成果驗證了加權最小二乘小波演算法重建真實物表覆蓋面的可行性，它也提供了

相關的品質指標，同時具有過濾點位高程的量測誤差之能力；因為每一個觀測方程式都是

線性的，它們組成的法方程式系統也直接就是一個線性系統，所以不需要給定未知數的起

始近似值、也不需要迭代計算，這也是加權最小二乘小波演算法在實務應用上的另一個優

點。吉布斯效應的影響範圍也確實會隨著小波法使用的解析力越細而縮小，使得重建精確

的零階不連續的物表面(包括房屋牆面、樹木邊緣輪廓)成為可行。加權最小二乘小波演算

法可以同時求定一個各個局部表面性質不同的覆蓋面，例如它包括了圓球型的屋頂面、山

形屋頂面、碎形的樹木表面、高聳的建築物牆面、圍牆、平坦的地面等。全部的實驗區重

建出擬合精度±14~23cm 的數值覆蓋面，已和輸入的空載光達點雲先驗高程精度±15~35cm

相當。

零階不連續的局部物表面區出現吉布斯效應，同時，大改正數的點群幾乎都出現在零

階不連續的局部物表面區，例如房屋的牆面、樹木的邊緣輪廓、籃球架。換言之，它們顯

示了應用吉布斯效應來自動偵測與重建房屋、樹木、橋樑等地上物之應用潛力。另外，從

這些不同地區的實驗成果也顯示了加權最小二乘小波演算法從離散的物表覆蓋面的點群來

求定數值地形模型（含地形斷線資料）具有良好的應用潛力，相關的課題將繼續深入研究

之。

另一方面，本研究也同時提出一些新構想和新演算法，擬從離散的光達點群來求定數

值地形模型（含地形斷線資料）

，應用吉布斯效應來偵測與重建房屋、樹木、橋樑等地上物，

萃取可靠的物面特徵，它結合了浮測模型法、三維坐標變換、最小二乘擬合法，進行整體

平差，期能萃取與重建可靠的物面特徵供光達測量之各式應用。

十二、參考文獻

內政部，2005。內政部『辦理LIDAR之高精度及高解析度數值地形測繪、資料庫建置與應

用推廣工作案』成果發表暨應用研討會邀稿序文。

王定平、王成機、陳思仁，2003。發展國家基本測量，第6 屆衛星定位測量研討會論文集。

蔡展榮，1998。取樣定理內的正交級數之訊號表達精度與收斂性，測量工程，第四十卷，

第三期，民國87年9月，第45~70頁。

蔡展榮，1999。吉布斯問題的一個實務解法，中國土木水利工程學刊，第11卷，第3期，民

國88年9月，第175~182頁。

Axelsson, P., 1999, Processing of laser scanner data－algorithms and applications, ISPRS

Journal of Photogrammetry & Remote Sensing, vol.54, pp.138-147, 1999.

Cobby, D.M., Mason, D.C. and Davenport, I.J., 2001, Image Prcessing of Airborne Scanning

Laser Altimetry Data for Improved River Flood Modeling, ISPRS Journal of Photogrammetry

& Remote Sensing, vol. 56, pp. 121-138.

Deussen, O., Hanrahan, P., Lintermann, B., Mech, R., Pharr, M., and Prusinkiewicz, P., 1998.

Realistic modeling and rendering of plant ecosystems, Proceedings of SIGGRAPH '98, pages

275 - 286, ACM SIGGRAPH, New York.

Ebert, D. S., 2003. Texturing and Modeling: A Procedural Approach. Morgan Kaufmann, USA.

Haala, N., Brenner, C., 1997a“Interpretation of Urban Surface Models using 2D Building

Information”, Automatic Extraction of Man-Made Objects from Aerial and Space Images(II),

Birkhauser Verlag, Berlin, pp.213-222,1997.

Jaan-Rong Tsay and Cheng Wei Yen, 2006. A Wavelet Approach for Determining A Surface

Covered with Airborne LiDAR Points. The CD of the proceedings of the 27

th

Asian

Conference on Remote Sensing (ACRS2006), 9-13 October 2006, Chinggis Khaan Hotel,

Ulaanbaatar, Mongolia, oral presentation on October 12, at the Technical Session M-2, room

1, Theme: Photogrammetry/Lidar. (sponsored by the National Science Council (NSC) under

the grants of NSC94-2211-E-006-072)

Jaan-Rong Tsay, 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.

Jaan-Rong Tsay, 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.

Jähne, B., 1991. Digital Bildverarbeitung. Springer Verlag, Berlin.

Kaiser, G., 1994,

A Friendly Guide to Wavelets, Birkhäuser, Boston.

Musgrave, F. K., Kolb, C. E., and Mace, R. S., 1989. The synthesis and rendering of eroded

fractal terrain, Proceedings of SIGGRAPH' 89, in Computer Graphics 23, 3, pages 41–50,

ACM SIGGRAPH, New York.

Musgrave, F. K., 1993. Methods for Realistic Landscape Imaging. PhD. thesis, Yale.

Peitgen, H.-O., and Saupe, D., 1988. The Science of Fractal Images. Springer-Verlag, New York.

Priestnall, G., J. Jaafar, and A. Duncan, 2000. Extracting Urban Feature from LiDAR Digital

Surface Models, Computers, Environment and Urban Systems, vol. 24, pp. 65-78, 2000.

Prusinkiewicz, P., and Hammel, M., 1993. A Fractal Model of Mountains with Rivers.

Proceedings of Graphics Interface '93, pages 174-180, Toronto.

Steinel , E., Kiema, J., Leebmann, J. And Bähr, L.H., 2001,” Laserscanning for Analysis of

Damages Caused by Earthquake Hazard,“ OEEPE Workshop on Airborne Laserscanning and

Interferometric SAR for Detailed Digital Elevation Models, Stockholm, pp.88-99.

誌謝

衷心感謝中興測量公司及林志交先生提供台中地區空載光達資料與說明資料參數，王

淼先生提供讀取和顯示光達資料之程式，同時，謝謝國科會補助研究經費（計畫編號NSC

94-2211-E-006-072與NSC93-2211-E-006-067）

，讓本研究各項實驗測試與研究工作得以順利

行政院國家科學委員會補助國內專家學者出席國際學術會議報告

95 年 10 月 22 日

報告人姓名 蔡展榮

服務機構

及職稱

國立成功大學測量及空間資訊學

系副教授

時間

會議

地點

95 年 10 月 9-13 日

蒙古烏蘭巴托城

本會核定

補助文號

NSC95-2221-E-006-339

會議

名稱

(中文)第 27 屆亞洲遙感探測研討會

(英文)27

th

Asian Conference on Remote Sensing (27

th

ACRS)

發表

論文

題目

(中文)求定空載光達點雲覆蓋面的一個小波演算法

(英文)

A Wavelet Approach for Determining a Surface Covered with Airborne LiDAR Points

報告內容應包括下列各項：

一、參加會議經過

亞洲遙感探測學會(Asian Association on Remote Sensing, AARS)成立於 1980 年，經

過相關教授們(例如中大陳哲俊教授、成大王蜀嘉教授及廖揚清教授)和國際大師們（例

如德國 Ackermann 教授、日本 Shunji Murai 教授）的努力，我國和大陸同時是創會會員，

分別名為 China Taipei 和 China Beijing，近年來，大陸團未經 AARS 同意逕行將他們的

會員名稱從原先的 China Beijing 改為 People Republic of China，不時對本國團會員名稱

和資格提出反對和抗議，幸 AARS 和本國代表們都具有共識，AARS 是以亞洲地區和國

家為單元，不是以國家為會員，且以促進學術合作交流為目標，不涉及政治，因此，每

每能在大陸團的各種動作下，化險為夷，讓台灣團繼續是 AARS 的創始會員。

今(95)年的 10 月 9~13 日在蒙古烏蘭巴托城成吉思汗旅館(Chinggis Khaan Hotel,

Ulaanbaatar, Mongolia)舉行第 27 屆亞洲遙感探測研討會(ACRS)，共有 20 個國家和地區

的會員參加，參加 ACRS 的人數達 389 人。其中，台灣團的註冊人數達 83 人，僅次於

日本 90 人，居第二位，比大陸團 28 人多。共計發表 294 篇論文，包括 168 篇口頭宣讀

的論文、和 126 篇海報型論文，其中包括農業、災害、地理資訊系統、森林、資料處理、

都市、製圖、合成孔徑雷達、海洋大氣氣象、環境生態系統、衛星系統、攝影測量與光

達、地質、土地使用、水文等共計 15 大主題。由於提供網路上傳論文的伺服器問題，

造成 40 篇摘要和 37 篇全文遺失，主辦單位公開表達歉意。另外，共計有 19 家公司和

單位參加這一次的 ACRS 大會的商業展示(commercial presentation)。

今年，台灣團表現依舊搶眼，文化表演(Culture Show)得第一名，並在日本攝影測量

及遙感探測學會 JSPRS 提供給 AARS 頒發給 35 歲以下年輕學者的論文獎中，台灣團中

央大學張智安博士生得到最佳海報型論文獎，讓與會外籍人士相當刮目相看。

為了參加這屆的 ACRS2006 研討會，我搭乘 10 月 7 日 03:40 時的統聯客運從台南

火車站前站出發，前往桃園中正機場，提早在台灣團的集合時間(10:00)前的兩小時（約

8:00）到達中正機場韓航櫃台前，由於航班的關係，幾經努力請兩家旅行社協助確定機

附件三

亞洲地區

其它地區

地區或國家

註冊人數

地區或國家

註冊人數

1

日本

90

法國

5

2

台灣

83

德國

5

3

蒙古

71

荷蘭

5

4

泰國

40

美國

4

5

中國北京

28

俄羅斯

1

6

韓國

27

瑞士

1

7

馬來西亞

10

英國

1

8

澳洲

6

小計

22

9

新加坡

4

10

尼泊爾

4

11

斯里蘭卡

2

12

印度

1

13

越南

1

小計

367

共計

389

口頭宣讀研討主題

場次

口頭宣讀研討主題

場次

1

農業

3 2

災害

3

3

地理資訊系統

3 4

森林

3

5

資料處理

3 6

都市

2

7

製圖

2 8

合成孔徑雷達

2

9

海洋大氣氣象

2 10

環境生態系統

2

11

衛星系統

2 12

攝影測量與光達

2

13

地質

1 14

土地使用

1

15

水文

1

共計 32 場次的口頭宣讀研討。另有 1 個學生場次(student session)以及 1 個

特別場次(special session)。

位後，我們成大共六人終於加入中大團，由於購買團票之故，所以到達機場後，我才拿

到登機證和護照。整團約 82 人搭乘 12:15 時從中正機場起飛的韓航 KE692 航班，15:40

時抵韓國首爾機場，在那裡，等候約 4 小時後，搭乘 19:50 時從韓國首爾機場起飛的韓

航 KE867 航班，22:15 時抵達蒙古烏蘭巴托機場，約 23:30 時，住進成吉思汗旅館，經

過一整天的行程，終於到達 ACRS2006 會場，因為打包行李，整夜未眠，所以到達旅館

後，很快就入眠了。

10 月 8 日星期日，也就是 ACRS2006 研討會的前一天，我先到會場辦理現金繳交註

冊費 USD100 等相關的註冊手續，並領取大會資料，包括論文 CD、大會手冊、蒙古傳

統音樂 CD、資料袋等。並在旅館房間裡，閱讀大會手冊，勾選擬前往聆聽的各場次宣

讀的論文。

10 月 9 日星期一 10:00~12:00，於另一個場地”文化劇院(cultural palace)”舉行開幕

式，大會主席 M. Saandar 博士分別邀請貴賓演講，包括蒙古教育文化科學部部長 Oe.

表 Y04

1 CNCRS (Chinese National Committee for Remote Sensing, China)

2

INFORMATICS AND REMOTE SENSING INSTITUTE OF MAS

(Mongolia)

3 MAP PRODUCT OF MONGOLIA LTD. (Mongolia)

4 NATIONAL UNIVERSITY OF MONGOLIA (Mongolia)

5 ERSDAC

(Japan)

6 MonMap

(Mongolia)

7 UNIVERSITY OF AGRICULTURE (Mongolia)

8 GISTDA

(Thailand)

9 SKY C&C (Mongolia)

10 JAXA (Japan)

11 RESTEC (Japan)

12 DIGITAL GLOBE (USA)

13 GAIA 3D (South Korea)

14 INFORMATION AND COMPUTER CENTER, MNE (Mongolia)

15 ITT VISUAL INFORMATION SOLUTIONS (Japan)

16 MACRES (Malaysia)

17 ER MAPPER (Australia)

18 BEIJING SPOT IMAGE (France)

19

MONGOLIAN UNIVERSITY OF SCIENCE AND TECHNOLOGY

(Mongolia)

共 19 個參展的公司和單位

Enkhtuvshin 博士、亞洲遙測協會秘書長 Shunji Murai 教授、國際攝影測量及遙感探測協

會 ISPRS 第一副主席 John Trinder 教授、日本太空署(Japan Aerospace Exploration Agency,

JAXA)副主席 Kaoru Mamiya 先生、蒙古國會議員 J. Gurragchaa 博士。在約 1 小時的貴

賓演講後，隨即進行非常精采的蒙古傳統文化歌舞表演。12:00 開幕式結束後，全部人

員再搭乘巴士回到會場成吉思汗旅館五樓用餐。下午有兩個場次的全員參加的研討

(plenary sessions)，分別討論蒙古的遙測發展、遙測對全球暖化效應研究之貢獻、日本

ALOS 衛星、大陸的風雲衛星和海岸衛星、台灣的福爾摩沙衛星、AVHRR 替代衛星

VIIRS、地球觀測系列衛星 MetOP 到 NPOESS，讓與會人士認識各國衛星發展與服務現

況。晚上，在另一個場地”Bayangol Hotel”舉行歡迎晚宴暨文化之夜；歷年來的 ACRS

研討會常將歡迎晚宴和文化之夜分別於不同的兩天舉行，這次大會將兩者合併舉行，因

此，韓國團因班機航班限制，來不及參加第一天的議程及文化之夜，儘管如此，還是讓

與會人士感受到強烈的文化洗禮，約二三十位身材高挑的蒙古小姐穿著不同的蒙古服飾

穿梭會場，讓各國人士與她們拍照留念，緊接著各國代表團的表演，我們台灣團由中大

八位女同學和八位男同學穿著台灣原住民服飾上台載歌載舞表演，台下的團員也隨著音

樂旋律自動一個接著一個手牽手地跳舞，一些老外也自動加入行列，帶動現場熱烈氣

氛，終場由瑞士 A. Gruen 教授宣布文化表演前三名的安靜時刻，首先是第三名由泰國隊

得獎，就在 A. Gruen 教授宣布第二名是日本隊的時候，現場台灣隊歡聲雷動，因為大家

都確定第一名就是非臺灣團莫屬，果然，A. Gruen 教授終於宣布第一名就是台灣中大團

的表演，不僅表演的學生很興奮，全團的老師們以及老外們都鼓掌肯定，就這樣又再一

次增加台灣在國際社會上的能見度並建立許多國際友誼，是一場很好的國民外交。

由於這次的 ACRS2006 共有四個研討會場，都位於成吉思汗旅館三樓，四個研討會

場彼此就在隔壁，距離很近，因此，我在緊接的 10 月 10 日到 10 月 13 日的每天研討活

動中，根據大會排定的時程場地，我預先挑選自己感興趣的論文，穿梭於不同的研討會

場，並逐一做一筆記，記錄我自己感興趣的新課題，頗覺收獲不少。

10 月 11 日星期三晚上，中華民國航空測量及遙感探測學會理事長宴請台灣團全體

團員，晚宴氣氛異常熱烈，拉近了台灣遙測界各單位人員之間的距離，達到一個很好的

聯誼效果。

我在 10 月 12 日星期四下午 15:30~17:30 時於第一個研討會場room 1 舉行的第M-2

場次（主題photogrammetry/LiDAR）發表自己的論文，刊載於Proceeding CD of the 27

th

ACRS, 9-13 Oct., 2006, Ulaanbaatar, Mongolia。口頭宣讀完成後，一些國外教授們都對我

們的研究內容很感興趣，例如韓國Kim教授，這場次的主席B. Enkhtuvshin也予以贊賞，

可惜我的學生嚴晟瑋先生因服兵役中而無法親自上場報告，否則他會有很高的機會得到

JARS頒發的青年論文獎。

我也利用這幾天的議程空檔，到一樓的廠商展示場地拜訪各家廠商攤位，並收集新

資料。

10 月 13 日星期五 11:00~12:00 舉行畢幕式，亞洲遙測協會秘書長 Shunji Murai 教授

報告了這次 ACRS2006 期間，各個會員國代表們於 10 月 10 日與 10 月 12 日兩天會議的

具體結論，包括確定第 28 屆 ACRS 研討會將於 2007 年 11 月 12-16 日於馬來西亞首都

吉 隆 坡 (Kuala Lumpur, Malaysia) 的 世 貿 中 心 (World Putra Trade Center) 舉 行 （ 網 址

http://www.macres.gov.my/acrs2007）

、第 29 屆 ACRS 研討會將於 2008 年在斯里蘭卡首都

可倫坡舉行。另外，也決定了 AARS 的學術期刊 Asian Journal of Geoinformatics (AJG)

的期刊主編日本本洲千葉大學(Chiba University, Japan)的 Ryutaro Tatershi 教授將任職到

2007 年一月底，2007 年二月起改由澳洲新南威爾斯大學教授 Prof. Bruce Forster 接任。

再者，也決定成立研討會籌辦支援委員會，以維持國際合作關係、維持 ACRS 研討會的

品質、維持標準化論文集、並組織特殊活動，例如教育課程、學生專屬的交流研討場次

(student session)等。另外，也建議未來的 ACRS 活動鼓勵各會員國的政府和民間公司的

決策負責人及主管能積極參加 ACRS、研討議題也能涵括全球議題、自動化產生數值高

程模型、正射影像、特徵物萃取，建議規劃出版 ACRS 屆滿 30 年的紀念刊物，包括製

作全亞洲無雲鑲嵌像片圖、亞洲各國首都的影像地圖、和出版介紹太空計畫與應用的一

本專書。

10 月 14 日星期六，我參加中華民國航空測量及遙感探測學會安排的參觀蒙古國家

野生動物園的活動，見識了蒙古黃色大草原的雄偉，於當天晚上，全體團員搭車前往烏

蘭巴托機場，搭乘 10 月 15 日星期日 00:20 時起飛的韓航 KE868 班機，04:20 抵韓國首

爾機場，等候約 5 小時後，搭乘 09:30 時起飛的韓航 KE691 班機，11:05 抵達中正機場，

返抵國門，正式結束了這次的 ACRS2006 活動，搭車返回台南。

二、與會心得

非常感謝國科會和中華民國航空測量及遙感探測學會分別補助這次活動的機票

費、簽證費、註冊費、保險費、和生活費各計新台幣 40000 元和 21000 元，使我能在較

少的自費負擔之下，有機會參與這次的國際性學術研討會，不但能吸收許多國外最新的

研究成果和研究現況，發表自己的研究論文，增加自己在這樣的國際會議上的歷練和經

驗，認識知名的國際教授們，增加未來的國際合作研究對象和機會，讓國際人士認識臺

灣在這領域的研究狀況，也能將攜回的資料和經驗分享並提供給有興趣的國內各校教授

和研究生們，對國內相關領域的研究提供最新的資訊暨資料，在此，特向國科會和中華

民國航空測量及遙感探測學會補助此行的相關費用，致以由衷萬分謝意！

每參加一次國際會議，就有一次的收穫和進步，提問題和回答問題更漸自然，在這

一次活動中，我在感興趣的幾篇論文研討活動中，提出問題，吸收許多別人的研究經驗

和心得，更獲得外國教授友人的友誼，其中，得以有機會和一些知名的教授們親自交談

聯誼，尤其難得，這些人物平時均相當忙碌，很難有機會和他們當面多談幾分鐘，例如：

ISPRS 前任主席暨現任第一副主席 J. Trinder 教授、大會主席 M. Saandar 博士、大陸武漢

大學李校長德仁教授等人，再者，有不少人對我在論文裡提出的構想和演算法感到興

趣，也因此得以在問題討論交流之中彼此互相認識，最讓我印象深刻。

三、考察參觀活動(無是項活動者省略)

略。

四、建議

雖然近年的國家財政狀況不佳，導致各項預算支出均必須刪減，若能客觀規劃積極支持

有意義的國際學術活動，增加臺灣在國際舞台上的能見度，也提攜國內年輕一代增加他

們的國際活動歷練，對國內科技生根紮實迎頭趕上甚至超越國際水準，都是一項值得投

資的多效益規劃，因此，建議在國家財政能力許可下，多鼓勵支持並補助國內學者參與

國際學術活動。以今年為例，國立中央大學太空及遙測研究中心的學生從教授們的研究

計劃結餘款獲得補助，得以出國參加是項大會，上台發表論文，也為台灣團拿到 culture

show 第一名，表現台灣年輕一代的活力、友善、以及傑出的學術研究表現，讓外國友人

對台灣團表現稱許，提高臺灣在國際舞台上的能見度，的確是一件相當可喜的現象，盼

成大也能仿傚中大的做法，讓自己的學生能增加國際歷練，邁向國際一流大學。

五、攜回資料名稱及內容

1. 第 27 屆亞洲遙感探測研討會論文光碟一片。

2. 蒙古傳統音樂光碟一片。

3. JERS-1 SAR Global Rain Forest Mapping Project, Insular South-East Asia & Papua New

Guinea, 100m L-band SAR image mosaics CD.

4. Sample Products Dataset for the Advanced Land Observing Satellite (ALOS), Earth

Observation Research Center, Japan Aerospace Exploration Agency.

六、其他

略

A WAVELET APPROACH FOR DETERMINING A SURFACE COVERED WITH AIRBORNE LIDAR POINTS

Jaan-Rong TSAY

Associate Professor, Department of Geomatics National Cheng Kung University 1 University Road, Tainan, Taiwan

Tel: (886)-6-2370876 ext. 838 Fax: (886)-6-2375764 E-mail: [email protected]

Cheng Wei Yen

M.Sc., Department of Geomatics, National Cheng Kung University 1 University Road, Tainan, Taiwan

KEY WORDS: Wavelet, Surface Reconstruction, Fractal, LiDAR

ABSTRACT: This paper presents a new approach for reconstructing object surface covered with 3D points. It utilizes

the 2D Daubechies scaling functions of 3rd order, which can describe fractal geometry, to formulate 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 designed on dyadic points, called PHO (Pseudo Height Observations) and POI (Pseudo Observations by Interpolation). Moreover, a full-automated weighting model is proposed 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. Tests are done by using airborne LiDAR points. They verify that the artifacts can be completely eliminated by adopting 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 of our tests are about ±20~23cm which are all to the extents of the a priori height accuracy, ±25cm.

1. INTRODUCTION

Especially in the past decade, the technique of light detection and ranging (LiDAR) becomes a well-known method for fast acquisition of precise coordinates of densely distributed 3D points on surfaces of interest objects. It is widely utilized in diverse applications, e.g. 3D cyber city modeling or generation of virtual environment. A key advantage of the so-called LiDAR-grammetry is its ability on fast acquisition of precise point cloud positions in 3D object space. Nevertheless, LiDAR-grammetry still has some very critical disadvantages:

(1) Unavoidable data errors cannot be detected and corrected without the corresponding a priori information or overlaying other data sources in the laser scanned area such as high-resolution (satellite or aerial) images or digital thematic maps or the overlapped LiDAR data scanned from other strips or other stations. Briefly to say, only LiDAR itself has completely no reliability.

(2) LiDAR points are almost always not located on surface features of interest objects. (3) Theresolution of LiDAR datais still limited until now.

Therefore, LiDAR should and must be integrated with other geomatic sciences and techniques in order to become a really powerful and applicable method for fast acquisition of GIS and space data with acceptable data quality (inclusive of accuracy and reliability).

This paper aims at another basic processing in LiDAR-grammetry, namely on object surface reconstruction by using the dense 3D point cloud covered on an interest surface. The dense 3D point cloud data can be acquired e.g. by laser scanning or digital photogrammetry or GPS.

For the present, the techniques for object surface reconstruction by using the dense 3D point cloud could be roughly categorized into two fields: (1) generation of a complex surface by connecting simple basic geometry primitives, and (2) surface reconstruction based on Euclidean geometry or non-Euclidean geometry such as fractal geometry. In the first field, the often used simple primitives are a triangle or a grid or a volume element (voxel) in 3D space. The constructive solid geometry (CSG) belongs to the first field. CSG can present a more complex surface by using high-order surface primitives, but the degree of complexity of the surface is still restricted by the adopted primitives (You et al., 2003).

The techniques in the second category can generate a surface with higher degree of complexity and discontinuity by using some surface functions and mathematical parameters such as Bezier surface and B-spline surface (Hill, 2000). For instance, Busé and Galligo (2004) adopted semi-implicit representation of algebraic surfaces. Anca et al. (2004) used spherical implicit functions for interactive modeling a surface. Xie et al. (2003) presented a method for piecewise C/sup 1/ continuous surface reconstruction of noisy point clouds via local implicit quadric regression. The non-Euclidean geometry formulates surfaces in non-orthogonal Euclidean spaces (Hill, 2000). Object surfaces are generated by some stochastic models and fractal dimension (Peitgen and Saupe, 1988).

Most methods in the afore-mentioned two fields still have some limitations. They only describe a piecewise smooth and continuous surface. Some of them cannot present a fractal surface. For example, Hoppe et al. (1992) reconstructed surfaces from unorganized points by utilizing iso-surfaces and the principle of marching cubes (Lorensen and Cline, 1987). Carr et al. (2003) applied the radial basis functions (RBF) to reconstruct smooth surface from noisy range data. This paper presents a wavelet approach for determining a complex 3D surface covered with airborne LiDAR points. The surface can be locally continuous, discontinuous, smooth, rough, regular, irregular, fractal, or non-fractal. Briefly to say, it is a real surface. We utilize the 2D Daubechies scaling function of 3rd order to define a surface function. In fact, the Daubechies scaling functions and wavelets display a fractal geometry (Kaiser, 1994).

2. SURFACE FUNCTION FORMULATED BY WAVELETS

Let f(x, y) be a surface function of two variables x and y. f(x, y) is a continuous single-valued real function. It can be approximately described in a multi-resolution manner according to the wavelet theory by the following equation:

)

,

(

)

,

(

1 ) 1 ( 2 1 1 ) 1 ( 2 1 , , , ,

∑

− −

∑

− = − − − =

⋅

=

M k N l l k j p j l k j j j

y

x

c

y

x

f

A

φ

(1) where j is a dilation parameter of the Daubechies scaling function fp of p-th order. In this paper, p=3 is adopted. In the

computation area, 3D object points with horizontal coordinates (Xi, Yi) are mapped linearly into the interval of 0§xi

§M-1 and 0§yi§N-1, "i. k and j are translation parameters in the x- and y-direction, respectively. is the scaling

coefficient for the scaling function f

j l k c, p,j,k,l. 3. OBSERVATION EQUATIONS

The airborne LiDAR provides 3D points with the coordinates (xi, yi, zi), where zi is equal to the original height Zi

determined by LiDAR. These 3D points are assumed to cover a surface, on which a point at (xi, yi) has only one height

value zi. The well-known least-squares adjustment (LSA) is applied here to determine a best fitting surface. In that

processing step, the following observation equation is given for each point.

∑

− −

∑

− = − − − =

⋅

=

=

+

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

φ

(2)

where vi is the residual of the height observation zi. The function value fp,j,k,l(xi, yi) can be computed. Therefore, the

linear equation only has the unknowns c_kj_,l. They can be determined by LSA.

4. ILL-POSED PROBLEM AND ITS SOLUTION

The airborne LiDAR produces very dense 3D points on different scanning lines. These points are not uniformly distributed in 3D space, i.e. the point interval is not constant in the entire computation area. The Daubechies scaling functions fp,j,k,l(x, y) are compactly supported. That means the function value fp,j,k,l(x, y) is equal to or very near to zero

for any point which is located outside the support or near the support boundary, respectively. On a resolution level, if no point is located inside the support or points are only located near the support boundary, the corresponding row (and column) vector for that unknown in the normal matrix is then equal to or near a zero vector. The so-called ill-posed problem emerges.

j l k

c_,

In order to solve the before-mentioned and often emerged ill-posed problem, a from-coarse-to-fine strategy is adopted. The approach starts from the coarser level, namely a smaller j and a larger interval between neighboring dyadic points.

In that case, the normal matrix is positive definite and all unknowns can be solved. While all coefficients are known, the approximation values on all dyadic points at the next finer level can be computed by the equation (3). They are called the pseudo height observations (PHOs) in this paper.

j l k c_, c_kj_,l

∑

− −

∑

− = − − − =

⋅

=

=

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

φ

(3) On the other hand, another pseudo observations are computed for all dyadic points at that next finer level. The nearest neighborhood interpolation is computed for each dyadic point. The neighborhood is a circle region centered at a dyadic point and with a radius sxy, where sxy denotes the a priori horizontal accuracy of LiDAR points. If there is no LiDAR

point inside the neighborhood, the pseudo observation PHO is adopted for that dyadic point. Otherwise, the interpolated value named POI is used as the pseudo observation for that dyadic point. In other word, each dyadic point has a either POI or PHO as its pseudo observation. Both are called POI and PHO points, respectively. They are used together with all original LiDAR points in LSA at that level. In this manner, a precise surface at the finest level can then be determined.

5. TEST DATA

In order to test the efficiency of our approach, airborne LiDAR data are used. They are acquired by the Optech ALTM 3070 system. The scanning density is about one points per meter. These LiDAR points have the a priori horizontal and vertical accuracy of about ≤1m and ±25 cm, respectively. Figure 1 illustrates our three test areas A, B, and C. Table 1 shows they contain different types of surfaces such as buildings, trees, semi-spherical roof, etc.

Figure 1. LiDAR points in the test areas A, B, and C. Table 1. Illustration of the test areas A, B, and C

Test area A Test area B Test area C

LiDAR point cloud

Google Earth image in the test area

Number of LiDAR points 10921 9578 10047

Horizontal dimension

∆

X

≅

∆

Y

≅

100 m

∆

X

≅

∆

Y

≅

100 m

∆

X

≅

∆

Y

≅

100 m

Table 2. Statistic figures of the computations in the test area B for the levels j=0 to j=3.

j Groundel size (m) Number of unknowns Number of POI points Number of PHO points ±

σ

ˆ₀(m)

0 10 121 × × 3.73

1 5 441 × × 2.64

2 2.5 1681 1645 36 1.97

3 1.25 6561 6439 122 1.50

Table 3. Some statistic figures of the computations in the test area B for the levels j=2 and j=3 j Number of unknowns Number of POI points Number of PHO points W1 W2 ±

σ

ˆ₀(m)

2 1 1 1.97

2* 1681 1645 36 1 0.0076 0.30

3 1 1 1.50

3* 6561 6439 122 1 0.0072 0.23

6. GIBBS EFFECT

Table 2 expresses some statistic figures of the computations in the test area B from the coarse level j=0 to the fine level j=3. The groundel size means the distance between neighboring dyadic points. The finer the resolution is, the better the fitting accuracy of the determined surface is. In case of j=3, the a posteriori standard deviation is ≤1.5m which is significantly much larger than the a priori height accuracy of ≤25cm. Moreover, The determined surfaces also show clearly that more detailed surface structures are reconstructed with finer resolution. Nevertheless, unreasonably large undulation surfaces are reconstructed near local discontinuous surface such as walls, boundary of high trees, electric poles, etc. It is the well-known Gibbs phenomenon.

7. WEIGHTING MODEL

In order to eliminate the Gibbs effect, we propose a full-automated weighting model to reduce the weights of the LiDAR points whose absolute residuals are larger than twice the a priori height accuracy of the LiDAR data, namely |vi|>2sZ. The equation (4) shows their weights, where σˆi denotes the ideal a posteriori standard deviation, e.g.

i

σˆ=sZ=≤25 cm in our tests. The number of such LiDAR points with larger residual is r and the sum of their square

residuals is denoted by . The total number of LiDAR points and POI points and PHO points as well is expressed by n. These LiDAR points with larger residuals own the weight p. The other points have the unit weights.

e r v ] [ 2 n u n v r p _e r i ) ( ] [ ˆ2⋅ ₂ ⋅ − =

σ

(4)

8. TEST RESULTS AND ANALYSIS

Table 3 shows some statistic figures of the computations in the test area B for the levels j=2 and j=3, where * denotes the reconstructions using the weighting function (4); W1 and W2 are the weights of observations with |v|§2sZ and |v|>2sZ, respectively; sZ denotes the a priori height accuracy (sZ = ≤25cm in our cases). Apparently, the a posteriori standard deviation of unit weight is successfully reduced to

σ

ˆ

₀=≤23cm after the afore-mentioned weighting model is

utilized. Figure 2 illustrates the location of the LiDAR points with larger residuals, namely with the weight W2, in test area B for j=2 and j=3. Apparently, almost all of those LiDAR points with larger residuals |v|>2sZ are more nearby located at break-lines than the case at coarser level. Similar results are also obtained in another two test areas A and C. They indicate that our approach has the potential for detecting local discontinuous structures on object surfaces. Moreover, the test results verify that the afore-mentioned weighting model really can eliminate the Gibbs effect. Figure 3 shows the reconstructed surfaces with j=2 and j=3 after the weighting model is adopted. Similar results are also obtained in the test areas A and C. They demonstrate the good ability of our approach for simultaneously reconstructing different types of surfaces such as local plane, crown canopy of trees, semi-spheric roof, etc., in a computation.

Table 4 and Table 5 show the statistic figures of the computations in the test area A and C for the levels j=0 to j=3. They demonstrate similar results as to the case in the test area B as shown in Table 3. While the interval of dyadic

points reaches the level of LiDAR point interval, the a posteriori standard deviations of unit weight in the test areas are about ≤20cm, ≤23cm, ≤22cm. They are all to the extents of the a priori height accuracy, ±25cm.

Figure 2. Black points denote the observation points with the weight W2 in the test area B for j=2 (left) and j=3 (right)

Figure 3. LiDAR points draped on the reconstructed surface of the level j=2 (left) and j=3 (right) in the test area B Table 4. Statistic figures of the computations in the test area A for the levels j=0 to j=3

j Groundel size (m) Number of unknowns Number of POI points Number of PHO points _{W1 W2}

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 Table 5. Statistic figures of the reconstructions in test area C for the levels j=0 to j=3

j Groundel size (m) Number of unknowns Number of POI points Number of PHO points W1 W2 ±

σ

ˆ₀(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 Figure 4 illustrates the reconstructed surfaces draped with the corresponding Google Earth images in these three test areas A, B, and C. The afore-mentioned concluding remarks are verified that a finer surface can be reconstructed by our approach with finer resolution, namely larger j.

9. CONCLUDING REMARKS

This paper presents a wavelet approach for determining a complex 3D surface covered with airborne LiDAR points. The surface can be locally continuous, discontinuous, smooth, rough, fractal, or non-fractal. We utilize the 2D Daubechies scaling function of 3rd order, which can describe fractal geometry, to write the observation equations of the point cloud. Furthermore, the linear system is solved by the least-squares adjustment and the surface can then be determined. To overcome the often emerged ill-posed problem, we employ a from-coarse-to-fine strategy and use the pseudo observations on dyadic points, namely PHO (Pseudo Height Observations) and POI (Pseudo Observations by Interpolation), to stabilize the linear system and to get the solutions.

Our experimental results show that with assistance of the pseudo observations on dyadic points, the algorithm can yield a stable solution and can meet with our basic hypotheses. Besides, we find: ① Gibbs effect: some irregular artifacts emerge around the abrupt areas of the reconstructed surface. ②points with larger residuals are largely located on abrupt areas such as walls, poles, or isolated trees. To eliminate the artifact effect, we propose a full-automated weighting model to reduce the weights of the points whose absolute residuals are higher than twice the a priori height accuracy of the LiDAR data. The results reveal that by combining the pseudo observations and the weighting model, artifacts can be completely eliminated and the a posteriori standard deviation of unit weight of the reconstructed surface can reach the same level of the height accuracy of LiDAR data points. For instance, while the interval of dyadic points reaches the level of LiDAR point interval, the a posteriori standard deviations of unit weight in our tests are about ± 20~23cm and are all to the extents of the a priori height accuracy, ±25cm.

By comparing to the diverse currently available surface reconstruction algorithms, the proposed approach can handle irregular, non-smooth, and fractal signals quite well and significant surface features registered in the original discrete sample points can be clearly expressed in the reconstructed surface. After some computation parameters are manually given, without the need on any other data preprocessing, our reconstruction system can automatically reconstruct a precise, highly complex, and multi-resolution surface model from discrete LiDAR points.

ACKNOWLEDGEMENT

We sincerely appreciate the National Science Council at Executive Yuen, R.O.C. for supporting this work under the project numbers NSC94-2211-E-006-072, and the Chung Hsing Surveying Co., Ltd. for providing the airborne LiDAR data for our tests.

REFERENCES

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. & 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.

Hill, Francis S., 2000. Computer graphics using OpenGL, 2nd edition. ISBN: 0023548568, Prentice Hall, 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. & 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 Peitgen, H. O. & Saupe, D., 1988. The Science of Fractal Images. Springer-Verlag, New York.

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.

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.