趙鴻森， 馮 琦， 周德云

(西北工業(yè)大學(xué)電子信息學(xué)院，西安 710129)

0 引言

機載數(shù)字高程模型(Digital Elevation Model，DEM)數(shù)據(jù)信息量龐大［1］，100 km2的3 m分辨率DEM數(shù)據(jù)就包含11億個數(shù)據(jù)，以32位字節(jié)的浮點數(shù)表示其高程，將至少需要8.5 GB存儲空間。另一方面，這些數(shù)據(jù)因為受到目前各種通訊網(wǎng)絡(luò)帶寬的限制，無法高效地進行網(wǎng)絡(luò)傳輸。未來DEM數(shù)據(jù)要顯示在手持終端(PDA)上，關(guān)鍵的一個問題就是DEM數(shù)據(jù)的壓縮與解壓算法的硬件實現(xiàn)。為了提高數(shù)據(jù)傳輸?shù)男?，以及解決這一類海量數(shù)據(jù)的存儲問題，必須進行DEM數(shù)據(jù)的壓縮。研究數(shù)字高程模型數(shù)據(jù)的壓縮技術(shù)有重要的意義。

機載DEM數(shù)據(jù)無損壓縮可以采用普通的柵格數(shù)據(jù)壓縮方式，如游程編碼、塊碼等，但是由于DEM數(shù)據(jù)反映了地形的連續(xù)起伏變化，通常比較“破碎”，普通壓縮方式難以達到很好的效果。因此對于網(wǎng)格DEM數(shù)據(jù)，可以采用哈夫曼編碼進行無損壓縮。有時，在犧牲細節(jié)信息的前提下，可以對網(wǎng)格DEM進行有損壓縮，通常的壓縮大都是基于離散余弦變換或小波變換的［1－2］，由于小波變換具有能較好保持細節(jié)的特性，因此，將小波變換應(yīng)用于DEM數(shù)據(jù)處理的研究較多。

已有壓縮方法都沒有充分考慮DEM數(shù)據(jù)不同于其他海量數(shù)據(jù)集的特點:地形的連續(xù)變化特征，鮮有高程的劇烈變化部分;幾乎所有分塊地形都由一些特征點或特征線確定，如山脊線、山谷線等，因此，壓縮效率(比率)還有潛力可挖。提取DEM特征點(線)，以這些點集作為樣本集采用RBF神經(jīng)網(wǎng)絡(luò)來完成壓縮，將存儲的網(wǎng)格點高程轉(zhuǎn)換為存儲神經(jīng)網(wǎng)絡(luò)權(quán)值，可以進一步提高壓縮比率。

1 原始DEM數(shù)據(jù)去噪

由于測量設(shè)備和傳輸介質(zhì)等的不完善，DEM數(shù)據(jù)在其形成、傳輸記錄過程中往往會引入多種噪聲［3］，噪聲往往表現(xiàn)為或大或小的極值。濾波去噪，即在盡量保留地形細節(jié)特征的條件下對噪聲進行抑制，其處理效果將直接影響到后續(xù)特征點提取的有效性。

噪點大部分是隨機性的，它們對其他位置的點不會產(chǎn)生影響，和鄰近各點相比，該點的某一坐標分量值將會有明顯的不同?；谶@一分析，可以用鄰域平均的方法來判斷每一點是否為噪聲點，并用適當插值方法消除。在一個3×3的柵格窗口中，如Zi，j表示該點的高程值，直接利用中心格網(wǎng)高程值與8個鄰域格網(wǎng)點的高程數(shù)值大小來進行判斷，如表1所示。

表1 噪點判斷Table 1 Noise judgement

式(1)中，ε為門限值，它可以根據(jù)對誤差允許范圍的程度由用戶設(shè)定選擇合適的值。這種方法可以去除DEM數(shù)據(jù)中的噪點，起到濾波的作用，如圖1、圖2所示。

圖1 含噪點的原始DEMFig.1 Original DEM

圖2 去噪后DEMFig.2 DEM excluding noises

2 特征點提取

地形特征點主要包括山頂點(peak)、凹陷點(pit)、脊點(ridge)、谷點(channel)、鞍點(pass)、平地點(plane)等［4］。利用DEM提取地形特征點，可通過一個3×3或者更大的柵格窗口，利用中心網(wǎng)格點與領(lǐng)域8個點的高程關(guān)系來進行判斷與獲取。即在一個局部內(nèi)，用X方向和Y方向上關(guān)于高程Z的二階導(dǎo)數(shù)的正負組合關(guān)系來判斷。

2.1 去噪后特征點的提取

判斷矩陣形式DEM數(shù)據(jù)特征點的方法類似于噪點的判斷，也是利用3×3的滑動窗口，具體方法為假設(shè)有一個 3 ×3 的窗口，如果(Zi，j－1－ Zi，j)(Zi，j+1－ Zi，j) ＞0，

1) 當 Zi，j+1＞ Zi，j則 VR(i，j)= － 1;

2) 當 Zi，j+1＜ Zi，j則 VR(i，j)=1;

如果(Zi－1，j－ Zi，j)(Zi+1，j－ Zi，j) ＜ 0，

3) 當 Zi+1，j＞ Zi，j則 VR(i，j)=1;

4) 當 Zi+1，j＜ Zi，j則 VR(i，j)= － 1。

如果1)和4)或者2)和3)同時成立，則VR(i，j)=2，如果以上條件都不成立，則VR(i，j)=0。

其中:VR(i，j)= －1，表示谷點;VR(i，j)=1，表示脊點;VR(i，j)=2，表示鞍點;VR(i，j)=0，表示其他點。

這種判斷特征點的方法利用了判斷點Zi，j周圍4個數(shù)據(jù)點，可以對山頂點、山谷點、鞍點做出判斷。

2.2 山脊線和山谷線確定

在地貌特征描述中，山脊線和山谷線發(fā)揮著巨大作用，可以在提取特征點之前確定下來。

山頂點通常是在山脊線上，山谷點通常存在于山谷線上。DEM數(shù)據(jù)特征點的求取可以先求出山脊線和山谷線，然后進行局部判斷，在山脊線上求出局部最大值為山頂點，在山谷線上求取局部最小值作為山谷點。

定義設(shè):L為選擇判斷山脊線的段長;N為選擇段長內(nèi)數(shù)據(jù)個數(shù);n為DEM數(shù)據(jù)間距。則N=L/n，取整數(shù)部分，N≥3。為了避免求取山脊線過程中無法求取次山脊線的情況，必須滿足L≤d，其中，d表示相鄰不相交山脊線之間的最小距離，也可根據(jù)經(jīng)驗取值。

設(shè):Zi，j為 DEM 數(shù)據(jù)點的高程值;Di，j為山脊線上的點;Gi，j為山谷線上的點;i，j為該數(shù)據(jù)點的位置。獲取d，取段長L=d，由N=L/n求取數(shù)據(jù)點最大個數(shù)N。第 i行上對選定段長數(shù)據(jù) Zi，j到 Zi，j+N進行比較，設(shè)段內(nèi)最大值為Zmax。

確定山脊線算法如下。

1) 若{|Zmax－ Zi，j+h|≤ε|1≤h≤N，0≤ε}，該段地形較平坦，不選取特征點。

2) 若存在唯一 Zmax且，則將Zmax保存為 Di，j，i，j表示 Zmax對應(yīng)位置。

若 Zmax=Zi，j并且 Zmax＞ Zi，j－1，或者 Zmax=Zi，j+N并且Zmax＞Zi，j+N+1，記錄 Zmax;若出現(xiàn)區(qū)域最大值，即 Zmax=Zi，j+N1，0≤N1≤Δ，Δ 表示區(qū)域大小。無法判斷段內(nèi)最大值時，判斷山脊線是否終止。

終止條件如下。

段內(nèi)出現(xiàn)區(qū)域最大值，即 Zmax=Zi，j+N1，0≤N1≤Δ，Δ 表示區(qū)域大小判斷 Zi+1，j+N1，0≤N1≤Δ:若 Zi+1，j+N1，0≤N1≤Δ 存在最大值 Z′max，Z′max坐標位置為 i′，j′，則Zmax的坐標選取為 I，j′;若 Zi+1，j+N1，0≤N1≤Δ 不存在最大值，繼續(xù)判斷 Zi+2，j+N1，0≤N1≤Δ，Zi+2，j+N1，0≤N1≤Δ存在最大值 Z′max，Z′max坐標位置為 i′，j′，則 Zmax的坐標選取為 i，j′;若 Zi+2，j+N1，0≤N1≤Δ 不存在最大值，則山脊線終止于Zmax。

3)在第i+1行繼續(xù)進行以上工作，找出第i+1行Zmax，直到找出N行。

4)將記錄的i～i+N個Zmax相互連接得到山脊線，若相鄰兩個最值的列坐標相差超過3則斷開山脊線，進行山脊線延伸。

在山脊線上相鄰3個數(shù)據(jù)中求出局部最大值Di，j。

山脊線延伸:若按行選取的最大值相鄰兩行之間兩個最大值分別為 Zi，j，Zi+1，j，若 |j′－ j|≥3，則山脊線在該行分為兩條，在斷開點山脊線不一定終止，需要在Zi，j和 Zi+1，j的領(lǐng)域內(nèi)進行判斷，在 i+1 行判斷距 Zi+1，j最近的3個點數(shù)據(jù)，取最大值作為該條山脊線的下一節(jié)點，繼續(xù)重復(fù)進行，直到相鄰相距最近的3個點不存在最值終止，對于Zi+1，j則是向上進行領(lǐng)域判斷，直到相鄰相距最近的3個點不存在最值終止。

5)如果{Zi，j＞ Zi+1，j+h|0≤h≤N}，Zi，j到 Zi，j+N可能為坡面上的點或者山脊點，這兩種情況在行掃瞄中無法確定，要求出完整的山脊線還需要在垂直方向上進行以上4步的工作。

這樣在設(shè)定的N×N的范圍內(nèi)就得到了山脊線、頂點，同樣的方法也可以用來求取山谷線。圖3為高程表上確定的山脊線，圖4為50×50范圍內(nèi)選取的山脊線和山谷線點。

圖3 DEM表上山脊線Fig.3 Ridge line on DEM

圖4 50×50山脊線Fig.4 50 ×50 ridge line

3 徑向基函數(shù)神經(jīng)網(wǎng)絡(luò)進行DEM數(shù)據(jù)壓縮

3.1 模型建立

構(gòu)建具有兩個輸入和一個輸出的正則化RBF網(wǎng)絡(luò)來實現(xiàn)從R2到R1的映射，從而完成DEM圖像的壓縮。網(wǎng)絡(luò)的輸入節(jié)點分別對應(yīng)平面坐標x，y值，輸出層節(jié)點對應(yīng)高程H值。

3.2 隱層單元數(shù)的確定

DEM間距一般都是固定不變的，為了減少存儲量，首先需要提取壓縮范圍內(nèi)的特征點，對于特征點較少的區(qū)域，如果計算結(jié)果誤差較大，可通過增大特征點選取過程中的閥值的方法增加特征點的個數(shù)，直至誤差滿足實際要求。即在特征點之間，特征點與區(qū)域邊界之間增加新的特征點，新的特征點初值選為其對應(yīng)位置的高程數(shù)據(jù)。

3.3 樣本集的確定

RBF網(wǎng)絡(luò)的優(yōu)點［5］是收斂速度快、逼近精度高、訓(xùn)練結(jié)果唯一，缺點是泛化性能較差。所以當其輸入偏離訓(xùn)練樣本集時，系統(tǒng)表現(xiàn)出的魯棒性就會很差。為提高網(wǎng)絡(luò)泛化能力，在DEM表中等間距選取樣本集。間距大小的確定與DEM的范圍和精度有密切的關(guān)系，還要考慮地面計算設(shè)備性能。因此，對于較大范圍的DEM數(shù)據(jù)，初始樣本點不宜選得過多，先按照粗略的訓(xùn)練樣本集進行訓(xùn)練，若不能滿足總樣本集的誤差要求，則應(yīng)當逐步減小采樣間距以細化樣本集，同時增加特征點個數(shù);如果是確定山脊線和山谷線求取特征點，則需要減小段長以獲得更多的特征點，基于DEM相關(guān)性獲取特征點的方法則需要增大誤差閥值，重新進行學(xué)習(xí)訓(xùn)練。

4 仿真算例

DEM數(shù)據(jù)取自ASTER GDEM，右上角為北緯25°東經(jīng)95°的一分度經(jīng)度×一分度緯度的實際DEM數(shù)據(jù)，按照30 m的間隔選取，總共有50×50個DEM數(shù)據(jù)。全局相對誤差精度選為2%，實際地形如圖5所示，樣本點數(shù)量為1250，特征點個數(shù)為150，因此隱層單元數(shù)選為150，初值均取1。壓縮后輸出地形圖如圖6所示。

圖5 實際地形圖Fig.5 Actual DEM

圖6 訓(xùn)練后網(wǎng)絡(luò)輸出地形Fig.6 Compressed DEM

為分析壓縮效率，分別選取50×50、60×60和70×70依次進行仿真實驗。

對比以上不同大小樣本集的仿真實驗，可看出隨著樣本集的擴大，數(shù)據(jù)最大的壓縮比率也在顯著增大，對于70×70大小的樣本集，就可以達到20∶1的壓縮比。

在實施大范圍DEM壓縮時，需要先通過分塊策略對原始DEM進行區(qū)域分割［6］，并建立分塊與相應(yīng)的RBF網(wǎng)絡(luò)權(quán)值及隱層單元之間對應(yīng)關(guān)系的檢索表，然后在各子塊中應(yīng)用該壓縮方法。

表2 各種大小的樣本集的對比仿真結(jié)果Table 2 Comparison of emulational outcome

5 小結(jié)

充分考慮DEM數(shù)據(jù)自身特點，通過提取山脊線、山谷線等地貌特征，以這些地形特征點作為樣本點訓(xùn)練集，從而具有根據(jù)地形自適應(yīng)確定網(wǎng)絡(luò)結(jié)構(gòu)的特點。采用RBF神經(jīng)網(wǎng)絡(luò)進行DEM數(shù)據(jù)壓縮，對于70×70大小的樣本集，即可以達到20∶1的壓縮比，并且隨著樣本點數(shù)目的增加壓縮比逐漸增大，滿足機載DEM壓縮和解壓要求。

［1］ 羅永，成禮智，陳波，等.數(shù)字高程模型數(shù)據(jù)小波壓縮算法［J］.國防科技大學(xué)學(xué)報，2005，27(2):118-123.

［2］ 劉愛利，閭國年.基于DCT域數(shù)字水印技術(shù)的DEM版權(quán)保護研究［J］.地球信息科學(xué)，2008，10(2):214-223.

［3］ DARNELL A R，TATE N J，BRUNSDON C.Improving user assessment of error implications in digital elevation models［J］.Computers，Environment and Urban Systems，2008，32:268-277.

［4］ 袁江紅，歐建良，查正軍.等值線DEM地形特征點提取與分類［J］.現(xiàn)代測繪，2009，32(3):3-6.

［5］ NAMPHOL A，CHIN S H，AROZULLAH M.Image compression with a hierarchical neural network［J］.IEEE Transactions on Aerospace and Electronic Systems，2008，32(1):327-337.

［6］ 施松新，葉修梓，張三元，等.基于分塊的大規(guī)模地形實時渲染方法［J］.浙江大學(xué)學(xué)報:工學(xué)版，2007，41(12):2002-2006.