基于小波分析的地形多分辨率建模方法

王江云 曾 瓊

(北京航空航天大學 自動化科學與電氣工程學院，北京100191)

王會霞

(北京航天自動控制研究所，北京100854)

早在20世紀50年代，針對地形建模，研究者們就提出了數字地形模型(DTM，Digital Terrain Model)的概念[1－2].DTM 按照模型各自功能的不同可以分為數字地面模型、數字高程模型(DEM，Digital Elevation Model)和數字地形(地貌)模型[3].DEM具有數據結構簡單，便于計算機管理、分析與繪制，數據可直接測繪獲取等優點，在實際工程中得到了廣泛的應用.但隨著DEM數據量的日益增長，對數據處理、繪制和存儲提出的要求也越來越高.另一方面，實際應用中并不是一直需要最高分辨率的原始DEM數據，在很多情況下可能需要的是較為宏觀的多分辨率層次的DEM數據，在不影響數據建模質量的同時，對DEM數據進行適當的簡化處理并實現高效的存儲和繪制成為相關領域的一個研究熱點.

采用小波分析技術能夠有效地解決以上問題.小波分析是一種有效、快速的時頻分析方法，作為對傅里葉分析的突破性發展，它既是一種快速的分析技術，又是一種快速的計算工具[4－5].小波分析技術所表現出的高壓縮比、高速簡化以及在變換中不改變原有信號特征的特點，使其在信號處理中被廣泛的應用.借助于小波變換技術，可實現對數據的時、頻域局部特征的分析和提取，以及對數據的多分辨率分析和表示.

1 小波變換原理

一個函數φ(t)成為小波函數的首要條件，稱為小波的容許條件，表示為

滿足小波的容許條件并具有帶通性、波動性和時頻局部化特性的平方可積函數φ(t)為一個基小波或小波母函數[6].由小波母函數φ(t)派生出的連續小波函數的表達式為

式中，a為伸縮因子或尺度因子;b為平移因子.

對于函數f(t)∈L2(R)，定義其小波變換為

小波變換根據數據的不同類型，又可以分為連續小波變換和離散小波變換.連續小波變換實現對數據的連續參數變化，一般適用于理論分析[7].離散小波變化實現了參數離散化的變換運算，更加適用于對實際的數值問題的解決.

小波函數的尺度及位移因子a和b的離散化表達式為

式中，a0和b0為任取的滿足a0＞1，b0＞0的變換因子;m和n為平移因子，均為整數.

將連續小波函數中的a，b用離散化后的形式表達即得到離散小波函數，表示為

相應的離散小波變換[8]表示為

在實際應用中，上述一維離散小波變換可以推廣到二維空間分析處理二維信號.文獻[9]以二進制小波變換為例，闡述了二維離散小波變化的原理，如下所示:

式中，φ(x)和φ(y)為定義的2個一維尺度函數;ψ(x)和ψ(y)為定義的2個一維小波函數;φ(x，y)為構造得到的二維尺度函數;ψH，ψD和 ψV為構造得到的3個二維小波函數.

以上3個小波函數將沿著不同方向的響應信號值變化:ψD對應對角線方向上信號的變化，ψH對應水平方向即列方向上信號的變化，ψV對應垂直方向即行方向上的信號的變化.

這樣尺寸為M×N的二維信號f(x，y)的二維離散小波變換表達式為

式中i∈{H，D，V}代表小波變換的3個方向.

2 DEM的小波變換方法

基于小波的多分辨率分析方法可用于對數據的多尺度近似、建模，以及對圖像的表現分析和可視化等領域[10－11].從多分辨率分析的角度講，小波分析能夠使數據中任意的局部和細節特性被充分地放大和暴露，因而小波被譽為“數學顯微鏡”.小波分析可應用于DEM數據的處理，從而建立多分辨率的DEM數據模型.

DEM是一種描述地形起伏變化的數字模型，通過對地形區域中一組離散點的平面坐標和高程值的測定建立.通過使用特定的格網間隔對地面高程值進行采集，并投影為在二維平面上符合規則化的格網分布規律的空間點的方法建立的DEM，稱為規則格網DEM.由此可知DEM數據在結構上與圖像數據具有一定的相似性，它等價于一個二維的灰度圖像，可以作為二維信號處理.

另一方面，地形所表現出的起伏特性包括了整體的起伏、變化規律，以及疊加于整體起伏、變化規律之上的局部、細微的起伏、變化規律.因此，如果使用小波分析方法對規則格網DEM進行多分辨率分析與表示的話，則可以實現對變化相對平緩的地形基本骨架和變化相對劇烈的起伏細節特征的分解，最終得到的“低頻分量”可作為低分辨率的 DEM 使用[12].

文獻[13]依據多分辨率分析的理論提出了快速小波變換算法(Mallat算法)，大大降低了小波變換的計算量和時間復雜度.它的基本思路是:通過離散濾波來實現在獨立的兩級分辨率之間的基變換.從最高分辨率開始，遞歸地進行這種濾波;在每一步濾波后，都有一些高頻細節被分離出來，剩下一個更為粗糙、頻率更低的表示.這一個過程被稱為進行小波分解或分析，而它的逆過程則被稱為重構或合成.

Mallat算法處理一維信號的原理如圖1和圖2所示.

圖1 基于Mallat算法的小波分解Fig.1 Mallat algorithm based wavelet decomposition

圖2 基于Mallat算法的小波重構Fig.2 Mallat algorithm based wavelet reconstruction

本文以Mallat算法為基礎，提出對DEM數據的快速處理及參數選取方法.其原理是，通過對DEM進行行元素與列元素的多分辨率分析，實現對DEM的二維小波變換.在變換中，一種策略是首先對各行(或各列)進行一次多分辨率分析，再輪換維度，對各列(或各行)進行一次多分辨率分析;另一種策略是對行和列交替地進行多分辨率分析.經過這樣的一次二維小波變換，可將當前尺度的DEM數據分解成為下一尺度的4個部分，即一個低頻部分(總體地形)，以及水平、垂直和對角線方向上的高頻部分(細節地形)，基本的變換過程如圖3所示.

圖3 DEM的離散小波分解示意圖Fig.3 Discrete wavelet decomposition of DEM

以DEM的二級離散小波分解為例闡述DEM小波分解的基本過程:設DEM的尺寸為M×N，且該DEM數據的尺度為j=0.對其進行一級小波分解，得到尺度為j=1的分量.考慮分解時所進行的“下二抽取”操作，則原始DEM相應地分離成為4幅尺寸相等且為M/2×N/2的子圖，如圖4所示.圖4中LL1表示的是尺度為j=1的概貌分量，表示地形的骨架特征，HL1，LH1和HH1為3個細節分量，分別表示水平、垂直和對角方向的地形的細節特征，并且根據下二抽取原理，每幅圖像的數據量均為原始DEM數據量的1/4.再將低頻部分即LL1視為原始模型，并對其執行與上一步相同的分解操作，可得到4個尺寸相同且均為M/4×N/4的子圖，其中LL2為尺度j=2的概貌分量，表示地形的骨架特征.

上述過程迭代進行，逐級地對 LLj(j=1，2，…)進行分解，就可得到一組分辨率不同的DEM序列.

DEM數據被變換到小波空間后，小波基及小波系數就形成了對表面的逼近.在實際應用中，對DEM進行小波變換后，往往并不是直接將低頻部分的平滑骨架信息作為DEM的低分辨率模型，而是需要根據不同的仿真需要，適當地保留部分高頻的細節地形信息.這就需要尋找一個合適的準則來控制逼近的精度，對此本文采取的方法是根據仿真需要選擇恰當的閾值或濾波門限來對小波系數進行濾波.濾波完成后將保留下來的細節信息加入數據集中進行小波逆變換進行重構，即生成了多分辨率DEM數據.

基于小波的DEM數據多分辨率建模的算法流程如圖5所示.

圖5 基于小波的DEM數據多分辨率建模算法Fig.5 Wavelet based multi-resolution modeling of DEM

3 小波函數及邊界延拓方式的選取

大量學者對基本的小波函數(系)的構造進行了研究，并提出一些經典、常用的小波系可供分析時選擇，包括Haar小波、Daubechies小波等.不同類型的小波函數體現出不同的時、頻域特性，因此對于特定的應用領域，應當選取與研究目的相適應的小波.DEM數據量龐大，對小波變換的選取有比較明確的需求:

1)必須是二維離散小波變換;

2)能夠較好地實現對原始數據的逼近;

3)能夠支持對數據的完全重構;

4)變換后的數據量小，提高實現效率.

對DEM數據多分辨率建模而言，應當選取對稱、正交并且具有緊支撐性、平滑性和高階消失矩的小波系，以滿足分析的需要[14].但是實際上兼具這些特性的小波函數是不存在的，只能選擇最為適當的小波函數.

另一方面，多分辨率分析中的Mallat快速算法是針對無窮區間信號的.DEM數據是一個行和列均有限的二維矩陣，直接對DEM數據進行小波變換會在矩陣邊界產生數據失真，且這種失真效應會隨著變換尺度的增加而加劇[15].解決這一問題的方法是在小波變換之前，首先進行DEM數據邊界的延拓處理[16].常用的邊界延拓處理方法主要有補零延拓、平滑延拓、周期延拓和對稱延拓[17－18].

設DEM數據矩陣的各行或各列為f(n)，邊界延拓后為g(n).4種延拓方式的表達式如下:

平滑延拓和對稱延拓在邊界上是連續的，而補零延拓和周期延拓在邊界上是不連續的.4種延拓方式各有特點，在進行DEM小波變換時，必須通過仿真實驗選擇最適合的邊界延拓方式.

4 實驗結果與分析

在基于小波的DEM數據多分辨率建模研究中，多分辨率分析的精度受小波函數類型的選擇，以及原始數據邊界延拓方式的選取等主觀因素的影響.基于本文提出的小波處理方法，利用Matlab工具，對DEM數據進行多分辨率分解與重構，通過對重構結果與原始DEM數據的比較分析，計算出誤差信息作為評價指標，從而得出最佳的小波函數及邊界延拓方式.

設原始DEM的高程矩陣值為 f(i，j)，重構DEM 的高程矩陣值為 g(i，j)，其中 i=1，2，…，m，j=1，2，…，n，則有

其中，eRMS能夠從整體意義上描述重構的DEM數據偏離原始數據的程度，是目前最為廣泛接受的精度評價指標.

本文選取的DEM仿真實驗數據來源于美國USGS官方網站，選取的是美國西部山地類型DEM數據.Matlab仿真實驗的結果如表1所示.

由表1中的仿真結果可以看出，coif1小波和周期延拓方式為基于小波的DEM數據多分辨率建模的最佳參數選擇，得到的重構后的DEM精度最高.

表1 仿真實驗結果Table 1 Simulation result

本文選取的仿真數據為山地類型DEM，最大高程3390m，最小高程2094m，高差1296m，分辨率30 m.采用coifl小波和周期延拓模式對DEM進行多尺度變換，可獲得較優的變換結果.對所選DEM進行多尺度變換，效果如圖6所示.

圖6 原始、一尺度近似、二尺度近似DEM對比圖Fig.6 Original DEM and its one-scale and two-scale approximation

實驗結果表明，小波分析方法適用于大規模數據建模，并能快速地完成對數據的處理，并且在基于小波的DEM數據多分辨率建模研究中，通過小波分析方法處理后的各層次DEM數據集在數據總量上不會增加，低分辨率模型數據僅為原始數據總量的25%，大大提高了存儲設備的利用效率.

5 結論

本文在經典的Mallat多分辨率分析算法的基礎上，對DEM建模問題進行了研究:

1)提出了基于小波的DEM數據多分辨率建模方法，得到的地形多分辨率模型精度高，并且具有較好的運算效率和存儲效率;

2)提出了小波函數的選取與DEM數據的邊界延拓方法，實現對DEM建模所需小波函數的選取和通過邊界延拓方式解決數據失真問題;

3)提出了重構DEM的精度評估的指標和方法，并通過仿真實驗和比較分析給出最適合DEM多分辨率建模的小波參數.

References)

[1]Moore I D，Grayson R B，Ladson A R.Digital terrain modelling:a review of hydrological，geomorphological，and biological applications[J].Hydrological Processes，1991，5(1):3 －30

[2]Weibel R，Heller M.Digital terrain modelling[M].Oxford:Oxford University Press，1993

[3]Wang G X，Zhang，Y B，Li J.Research and practice of accuracy evaluation method in DEM[J].Science of Surveying and Mapping，2006，31(3):73 －75

[4]Torrence C，Compo G P.A practical guide to wavelet analysis[J].Bulletin of the American Meteorological Society，1998，79(1):61－78

[5]Wu J，Amaratunga K.Wavelet triangulated irregular networks[J].International Journal of Geographical Information Science，2003，17(3):273 －289

[6]Burrus C S，Gopinath R A，Guo Haitao，et al.Introduction to wavelets and wavelet transforms:a primer[M].Upper Saddle River:Prentice Hall，1998

[7]Torrence C，Compo G P.A practical guide to wavelet analysis[J].Bulletin of the American Meteorological Society，1998，79(1):61－78

[8]Daubechies I.The wavelet transform，time-frequency localization and signal analysis[J].IEEE Transactions on Information Theory，1990，36(5):961 －1005

[9] 黃為，魏迎梅，宋漢辰，等.基于小波分析的 DEM 數據多分辨率表達研究[J].計算機仿真，2010，27(5):80 －83

Huang Wei，Wei Yingmei，Song Hanchen，et al.Multi-resolution representation of DEM based on wavelet analysis[J].Journal of Computer Simulation，2010，27(5):80 －83(in Chinese)

[10]Yang J S，Zhou L P.Iterative algorithm of wavelet moments method based on the multi-resolution analysis characteristics of wavelet[C]//2011 2nd IEEE International Conference on Microwave Technology and Computational Electromagnetics.Piscataway，NJ:IEEE，2011:462 －465

[11]Gaouda A M，Salama M M A，Sultan M R，et al.Power quality detection and classification using wavelet-multiresolution signal decomposition[J].IEEE Transactions on Power Delivery，1999，14(4):1469 －1476

[12]Mallat S.A wavelet tour of signal processing[M].Salt Lake City:Academic Press，1999

[13]Mallat S G.A theory for multiresolution signal decomposition:the wavelet representation[J].IEEE Transactions on Pattern Analysis and Machine Intelligence，1989，11(7):674 －693

[14]Danovaro E，De Floriani L，Papaleo L，et al.A multi-resolution representation for terrain morphology[C]//Geographic Information Science.Berlin Heidelberg:Springer，2006:33 － 46

[15]Danovaro E，De Floriani L，Magillo P，et al.Morphology-driven simplification and multiresolution modeling of terrains[C]//Proceedings of the 11th ACM International Symposium on Advancesin Geographic Information Systems.New York:ACM，2003:63－70

[16]袁禮海，宋建社.小波變換中的信號邊界延拓方法研究[J].計算機應用研究，2006，23(3):25 －27

Yuan Lihai，Song Jianshe.Research on signal extended methods in wavelet transform [J].Application Research of Computers，2006，23(3):25 －27(in Chinese)

[17]Herley C.Boundary filters for finite-length signals and time-varying filter banks[J].IEEE Transactions on Circuits and Systems II:Analog and Digital Signal Processing，1995，42(2):102－114

[18]Kharitonenko I，Zhang X，Twelves S.A wavelet transform with point-symmetric extension at tile boundaries[J].IEEE Transactions on Image Processing，2002，11(12):1357 －1364