赵鸿森， 冯 琦， 周德云

(西北工业大学电子信息学院，西安 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.