機載DEM特征點提取及壓縮方法研究

趙鴻森， 馮 琦， 周德云

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

0 引言

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

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

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

1 原始DEM數據去噪

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

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

表1 噪點判斷Table 1 Noise judgement

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

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

圖2 去噪后DEMFig.2 DEM excluding noises

2 特征點提取

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

2.1 去噪后特征點的提取

判斷矩陣形式DEM數據特征點的方法類似于噪點的判斷，也是利用3×3的滑動窗口，具體方法為假設有一個 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個數據點，可以對山頂點、山谷點、鞍點做出判斷。

2.2 山脊線和山谷線確定

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

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

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

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

確定山脊線算法如下。

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

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

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

終止條件如下。

段內出現區域最大值，即 Zmax=Zi，j+N1，0≤N1≤Δ，Δ 表示區域大小判斷 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≤Δ 不存在最大值，繼續判斷 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行繼續進行以上工作，找出第i+1行Zmax，直到找出N行。

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

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

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

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

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

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

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

3 徑向基函數神經網絡進行DEM數據壓縮

3.1 模型建立

構建具有兩個輸入和一個輸出的正則化RBF網絡來實現從R2到R1的映射，從而完成DEM圖像的壓縮。網絡的輸入節點分別對應平面坐標x，y值，輸出層節點對應高程H值。

3.2 隱層單元數的確定

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

3.3 樣本集的確定

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

4 仿真算例

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

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

圖6 訓練后網絡輸出地形Fig.6 Compressed DEM

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

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

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

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

5 小結

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

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

［2］ 劉愛利，閭國年.基于DCT域數字水印技術的DEM版權保護研究［J］.地球信息科學，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］.現代測繪，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］ 施松新，葉修梓，張三元，等.基于分塊的大規模地形實時渲染方法［J］.浙江大學學報:工學版，2007，41(12):2002-2006.