1
數值地形模型特徵點選取的新方法
A New Approach of Feature Point Selection Method for DTMs
計畫編號：NSC 89-2211-E-151-006
執行期限：89 年8 月1 日至90 年7 月31 日
主持人：李良輝 國立高雄科學技術學院 土木系
一、中文摘要
本論文之研究內容即由小波理論之基 本概念出發，並經由實例分析之驗證，探 討利用小波理論為基礎進行數值地形模型 特徵點之選取。 實例分析中採用三種矩陣形式之小波 轉換：Haar 轉換、Daubechies D2 轉換， 以及 Quadratic Spline W-Matrix。分析結 果顯示，不同之小波因其自身特性的差 異，在地形起伏平緩抑或高差變化急遽之 區域分別有不同之適應性。其次，小波轉 換係數中隱含之高頻資訊，可作為選取地 形特徵點之依據；經與目前普遍使用之 VIP 法比較後，證實由小波轉換所選取之 特徵點具有更高的可靠性，可作為數值地 形模型分析、簡化與展現之重要參考。 關鍵詞：小波理論，數值地形模型，影像 處理ABSTRACT
In this study, we proceed from the fundamental concepts of wavelet theory first, and then take the evidence in the experiment to assess and investigate the feasibility of wavelet-based featue point selection for DTM data.
In the experiment we choose three wavelet in matrix form - Haar transform, Daubechies D2 transform, and quadratic spline W-matrix. The analysis result shows that different wavelet adapt to different sort of terrain, due to their inherent characteristics. In the second place, the high-frequency information in the wavelet transform coefficients can
be a secundum of selecting
characteristic points. Compared with the
method of Very Important Points (VIP) algorithm, we concluded that the points selected by the wavelet transform have more reliability and consequently the
method can be an considerable
reference of DTM analysis, compression, reduction, and presentation.
Keywords : Wavelet Theory, Digital Terrain Model, Image Processing
二、緣由與目的
隨著計算機軟硬體的快速發展，加上 高解析度衛星與航照影像的普遍使用，配 合數值地形模型（DTM）資料的使用，三 維空間的立體景觀模擬，便可栩栩如生地 展現在使用者眼前。然而，在許多實際應 用上，龐大的影像資料與地形資訊，對儲 存空間與計算時間造成相當的負荷；此 外，在某些應用中，高解析度與高精度的 資料乃非必要，多餘的資訊反而造成儲存 空間的浪費與處理效率的降低。因此，如 何針對不同應用的需求，應用多重解析度 分析（Multi-Resolution Analysis）來表達 必須且適當的資訊，便成為近年來一個重 要的研究課題。 在各種多重解析度分析模式中，小波 理論（Wavelet Theory）是近年來甚為流 行的一種分析工具，因為小波在時間域 （Time Domain）與頻率域（Frequency Domain）皆有良好的局部化特性，使其成 為極佳的多重解析度分析工具。目前在實 際應用領域中，諸如訊號與影像分析、地 震探勘、語音辨識與合成、電腦視覺及其 他 許 多 方 面 ， 皆 已 獲 得 突 破 性 的 成 果 [1-3]。本研究即擬以小波理論為基礎，進 行數值地形模型之多重解析度分析，並經 由實例分析探討其可行性。2
三、
小波理論基礎
連續小波轉換可說是為了要解決短時 傅利葉轉換（STFT）在解析度上之不足而 發展的，其基本概念與短時傅利葉轉換類 似，乃透過一分析窗對原函數進行分析。 關於小波理論背景的深入探討，可參閱 [4][5]。 定義函數 f(t)之連續小波轉換為∫
⋅ = Ψ = s f t t dt s CWTfψ(τ, ) ψf (τ, ) ( ) ψτ*,s( ) (1) 其 中 ψ(t) 稱 為 母 小 波 （ Mother Wavelet），而小波函數ψ 之定義則為：τ,s − = s t s s τ ψ ψτ 1 , (2) 在連續小波轉換中，透過改變尺度參 數（Scaling Parameter） s 與平移參數 （Translation Parameter）τ 來分析原函 數 中 所 存 在 的 頻 率 成 分 （ Frequency Component），而在實際之離散訊號分析 上，則使原始訊號通過一系列的高通與低 通 濾 波 器 （ High-Pass and Low-PassFilter）g[n]及h[n]來分析訊號的頻率。 m m n n n n g g h h g g g h h h g g g h h h D 2 2 2 1 2 1 2 1 2 1 2 1 2 1 . . .. .. .. .. × = （3） 濾波分析之過程，首先將轉換矩陣 D 乘上原訊號，所得之輸出訊號中包含高頻 與低頻資訊，經過次取樣(Subsampling) 的步驟將高低頻濾波訊號分別取出，以進 行後續的更高階轉換或分析。在 Haar 轉 換中，低通濾波向量h[n]（亦稱為尺度向 量 Scaling Vector）為[1,1]，而高通濾波 向 量 g[n] （ 亦 稱 為 小 波 向 量 Wavelet Vector ） 則 為 [1,-1] ， 長 度 均 為 2 ； 而 Daubechies D2 轉 換 之 h[n]=[-0.1294, 0.2241,0.8365, 0.4830]，g[n]=[-0.4830, 0.8365,-0.2241,-0.1294][6]，QS W-Matrix 之 h[n] 及 g[n] 則 分 別 為 [-1,3,3,-1] 與 [-1,3,-3,1]，並且在邊界值部分需加以修 正。在濾波分析過程中，每一次次取樣的 過程皆會將低頻訊號與高頻訊號加以分 離，亦即使樣本數減半，故在進行原訊號 之 重 建 時 ， 需 逆 向 進 行 上 取 樣 （Upsampling）以使樣本數加倍；而對二 維之網格式 DTM 資料而言，每一階完整 的轉換需經過下列步驟： (1)將轉換矩陣乘上 DTM 資料； (2)將乘積之奇數列（低通濾波輸出）與偶 數列（高通濾波輸出）加以分離而成為 二個新的矩陣； (3)步驟 2 所得之二矩陣加以轉置後，再分 別以轉換矩陣乘之； (4)對步驟 3 所得之二乘積，分別進行奇、 偶數列之分離而得到四個輸出訊號，分 別為整體低通濾波輸出及垂直方向、水 平 方 向 與 對 角 線 方 向 之 高 通 濾 波 輸 出，樣本數均為原 DTM 資料之 1/4。 由於每一階轉換均經過二階段之次取 樣，故在逆運算時，亦必須經過二次上取 樣，如此便可獲得與原訊號相同維度（或 樣本數）之一組訊號。若對完成一階小波 轉換所獲得之四種係數分別進行逆運算， 並 在 上 取 樣 的 步 驟 中 採 用 嵌 入 零 值 法 （Zero Padding），則得到的四個訊號皆 與轉換前訊號有相同之樣本數，其意義則 可解釋為原訊號之趨勢（低頻）分量、Y 方向高頻分量、X 方向高頻分量，以及對 角線高頻分量；其中趨勢分量保留了原 DTM 資料中的大致起伏資訊，而各高頻分 量中之數值較大者，則表示原地形在該處 之起伏較為劇烈。
四、多重解析度分析基本概念
現以一個有限頻寬的多頻訊號函數 x 為例 來扼要說明多重解析度分析的觀念。 x = T + D 圖 1 一個簡單多頻訊號函數 x 的分解3 舉例來說，該訊號可以是一維的影像 訊號或是一維的高程訊號。 x 可以分解成 一個趨勢訊號（Trend Component）T 與 一個剩餘訊號（或稱為細節分量，Detailed Component）D 的總和，亦即：x=T+D； 而 剩 餘 訊 號 可 再 分 解 為 一 個 低 頻 分 量 （Low-Frequency Component）D 與一lf 個 剩 餘 高 頻 分 量 D （ High-Frequency hf Component）的總和：D =Dlf +Dhf。 一般而言，趨勢訊號與細節低頻分量 的 總和T~=T+Dlf 具 有 比 T 還好的 解析 度，是一個比T 還好的訊號 x 之近似。根 據這樣簡單定性直觀的分析可簡明又有系 統的得出多重解析度分析。 趨勢訊號T 通常以如下的簡單線性函 數來表達： Z k a T k jk jk ∈ =
∑
φ , (4) 式中φ 為已知的基底函數，jk a 為其jk 適當的係數。所有依式（3-1）的形式所表 現的趨勢訊號之集合即構成一個所謂的近 似空間（Approximation Space）V ：j ∈ =∑
k jk jk jk j a a R V φ ; (5) T 是 x 在空間V 中的一個近似，用以j 下的符號表示之： j jx x A T = = (6) 其中A 為係數集合。同理，上述的定j 義也適用於另一個趨勢訊號T~與近似空間 1 + j V ，亦即： 1 1 , 1 , 1 ~ + + + +∑
= = = j k j k j k j A x x a T φ (7) ∈ =∑
+ + + + k k j k j k j j a a R V 1 1,φ 1, ; 1, (8) 由x 轉變到j x （亦即由T 轉變到Tj+1 ~ ） 定義了一個解析度變細的過程，例如由較 粗的空間V 變換到較細的空間j Vj+1。此二 個近似空間V 、j Vj+1之間的關係如式(9)所 示。亦即近似空間Vj+1可做如下的正交分 解： 1 + = ⊕ j j j W V V (9) 式 中W 稱 為 細 節 空 間 （ Detailed j Space），它是V 在j Vj+1裡的正交餘集/補 集。 數學家分別在此三個空間V 、j W 、j 1 + j V 定 義 了 正 交 基 底 函 數φ 、jk ψ 、jk k j+1, φ 。上述的三種信號 T 、D 、 Tlf ~可藉 著這些基底函數從原始訊號 x 粹取得出。 訊號函數 x 在兩相鄰近似空間V 、j Vj+1之 間的分解隨著 j 值變大而變得越來越細， 如此便表達了所謂多重解析度分析的根本 觀念。 如上所述，空間Vj+1可分解成兩空間 j V 、W ；相反的，空間j Vj+1可由空間V 、j j W 重組而得。與此相仿，T~可分解成兩訊 號T 與D ；反之，可由此二訊號 T 與lf D 重lf 組得出。這些關係在離散的網格式資料上 便演化得出吾人所謂的分解與重組的數學 關係式。五、實驗例分析
本研究的重點，乃在對網格式 DTM 資料，依小波理論的基本觀念進行轉換， 以轉換後所得之轉換係數進行多重解析度 分析。在粹取地形特徵點的方法上，則利 用 Daubechies D2 轉換，以其小波轉換係 數中之高通濾波係數為選點之依據；此 外，在分析中亦採用 VIP 法作為對照，以 進行選取地形特徵點之可靠度比較。 切取 40 米解析度台灣地區之原始 DTM 資 料 座 標 （ 255000,2550000 ） 至 （265240,2539760）之範圍為實驗區，共 256×256 個樣本點，其高程自 178m 至 1911m，該實驗區 DTM 之三維展現如圖 2 所示；由圖上可以看出，之所以選擇此一 範圍為實驗區，乃因此區同時包含西半部 之高山、河谷與東側之平原，亦即包含 DTM 資料之高頻部與低頻部，特徵明顯， 相當適合作為分析之用。4 圖 2. 實驗用 DTM 資料 定義地形特徵點為地表高程起伏較大 之點，則由小波轉換所獲得之高頻分量亦 可作為粹取地形特徵點之依據。理論上， D2 轉換之細節分量應可保存較多之高頻 資訊[9]，故在本節中僅就 D2 轉換與 VIP 選點之能力進行比較。首先由原 256×256 之整體細節分量矩陣中，以元素（3,3）為 左上角取出矩陣中間之 250×250 個元素 為新的高頻矩陣；如同 VIP 法選取地形特 徵點的步驟(2)，將一階小波轉換所得之三 種高頻分量相加而得到整體細節分量，再 經過預先設定的界定值篩選出其中元素值 較大者，記錄其網格點之座標(X,Y)與高程 值 Z。對 VIP 法，則依同樣方式取出顯著 性指標矩陣中之 250×250 個元素，與 D2 轉換同時進行選點。 對於兩種方法所篩選之地形特徵點， 為了評估其代表性與可靠度，另設計一實 驗以作為比較。首先設定主控制條件為壓 縮比（Compression Ratio），亦即所篩選 出之點數與原 62500 個點之比值。實驗之 操作程序如下： (1)對每一預設之壓縮比，分別計算所對應之 界定值與篩選點數； (2)紀錄所選之點的三維空間座標值，經 Kriging 方法內插還原為 250×250 個樣 本點； (3)計算每一樣本點與原 DTM 資料之差值； (4)計算差值之平均值與標準差。 計算結果如表 1 所示； 表 1 選取之特徵點經 Kriging 內插後與 DTM 差值 之比較（單位：m）
六、結論
小波轉換係數中的細節訊號，隱含空 間中的高頻變化資訊，利用此一資訊，可 將地形上高程變化較大的點位加以粹取得 出。經過與 VIP 法之實驗比較，證實在篩 選地形特徵點的能力上，小波理論確實優 於 VIP 法，可作為諸多應用中空間資訊分 析的有效工具。參考文獻
[1]柯正誼，何建邦，池天河，1993，數字 地面模型，中國科學技術出版社，第 1 ∼13 頁。 [2]李軍，1995，小波理論及其在攝影測量 與遙感中的應用，測繪通報第四期，第 8~13 頁。[3] Floriani, L.D., Magillo, P., Puppo, E., and Bertolotto, M., 1997, Variable
Resolution Operators on a
Multiresolution Terrain Model, 4th
ACM Workshop on Advances in
Geographic Information Systems
(ACM-GIS '96),RockVille, Maryland, November 1996, pp.121-128.
[4]Polikar, R., 1996, The Wavelet Tutorial, Part I ~ PartIV, p.1-60,
[5]Young, R.K., 1993, Wavelet Theory and Its Applications, Kluwer Academic Publishers, pp.1-38.
[6]Daubechies, I., 1991, Ten Lectures on Wavelets, CBMS-NSF Series Appl. Math., SIAM. VIP 法 D2 轉換係數 選點 方式 壓縮比 平均值 標準差 平均值 標準差 10％ -0.176 0.224 -0.219 0.148 20％ -0.025 0.131 0.292 0.106 40％ -0.189 0.070 -0.420 0.060 60％ 0.236 0.053 0.107 0.053 80％ 0.449 0.019 0.163 0.009
