雷达地形测绘DEM空洞插补方法研究

龙四春，周威，文佳胜，陈鹏琦

（1.湖南科技大学煤炭资源清洁利用与矿山环境保护湖南省重点实验室，湖南湘潭411201；2.湖南科技大学测量工程与形变监测研究所，湖南湘潭411201；3.湖南科技大学能源学院，湖南湘潭411201）

雷达地形测绘DEM空洞插补方法研究

龙四春1，2，3，周威2，文佳胜2，陈鹏琦3

（1.湖南科技大学煤炭资源清洁利用与矿山环境保护湖南省重点实验室，湖南湘潭411201；2.湖南科技大学测量工程与形变监测研究所，湖南湘潭411201；3.湖南科技大学能源学院，湖南湘潭411201）

免费公开的SRTM DEM数据用图广泛，但存在大量的空洞。针对SRTM DEM空白处填补现状，利用MATLAB编程特性，编写线性插值、三次插值、最邻近插值的MATLAB代码，选择湖南西部丘陵山地进行空白填补实验。通过对比分析其最大值、最小值、标准差、中误差，得出三次插值法填补的效果最佳，生产的等高线最光滑。该方法与结论能为类似地区DEM插补方法选择提供参考。

MATLAB；SRTM DEM；线性插值；三次插值；最邻近插值

0 引 言

航天飞机雷达地形测绘使命（Shuttle Radar Topography Mission，SRTM）数据是制作地形图与地形分析的宝贵资源。但由于雷达侧视的几何特征、航天飞机轨道设计、信号干扰、雷达阴影或回波滞后等因素的影响，导致部分地区出现数据空洞，尤其在水体和高山峡谷地区。因此，要增强SRTM数据的实用性，必须对其数据空洞进行填补［1－5］。MATLAB作为专门的数值计算软件，具有强大、专业的数据计算、分析和可视化等功能［6－9］，将之应用SRTM DEM空白插补可以提高效率。

国内已有许多学者尝试过各种方法对SRTM高程数据的空值区域进行了填补［1－5，10－12］，但由于填补的数据地理范围、操作者所掌握的技术以及其他辅助性的数据不同，从而在过程和结果上就有所差异。2009年，左美蓉用等高线内插生成分辨率与SRTM相同的DEM，然后将内插出来的高程值取整来填补SRTM高程数据的空值区域，但填补后的数据显得不连续，内插出来的值与SRTM数据在空值边界的值有较大差距［5］。2012年CSI（Consortium for Spatial Information）将等高线和SRTM数据整合起来再内插，由于中间只有等高线上的高程值，而且相距较远，对于大范围空值区域，内插结果很难消除带状效果［1］。同时，国外很多其他机构和个人也在研究DEM生成的填补方法和填补工具［1－7，10，13－25］，如表1所示。

表1 SRTM数据空洞填补工具

1 SRTM DEM及MATLAB内插原理

1.1 SRTM DEM数据

SRTM数据每经纬度方格提供一个文件，精度有1arc－second和3arc－seconds两种，称作SRTM1和SRTM3，或者称作30m和90m数据，SRTM1的文件里面包含3601×3601个采样点的高度数据，SRTM3的文件里面包含1201×1201个采样点的高度数据，每一个小像元代表边长约90m的正方形区域，如图1所示。

图1 SRTM DEM文件示意图

图2 SRTM DEM空白示意图

由于天线杆和姿态的量测精度、记时误差、多路径效应、相位量测误差及雷达的热噪声等的影响，在SRTM中会出现一些空白值，如图2中红色方框所示。

1.2 SRTM DEM数据MATLAB内插原理

通常，内插贯穿在DEM的生产、质量控制、精度评定和分析应用等各环节。DEM内插实质是根据分布在内插点周围的采样点高程值求出待定点的高程值过程，根据内插点分布范围，可以将内插分为整体内插、分块内插和逐点内插三类［10］。DEM内插方法的选择，要考虑诸多因素，不仅要满足DEM内插精度要求，还要尽可能顾及计算效率，也要考虑地面复杂函数和已知数据特点。

MATLAB语言是一个基于矩阵和矢量的高级语言，简单易学，又具有面向对象的编程特点，编程效率高；有众多的应用工具箱，具备很强的开放性，除内部函数外，用户可通过对源文件的修改或加入自己编写的程序语句去构成新的专用工具箱。

DEM内插算法的空间相似性反映在由未知点附近已知点高程的加权平均值来确定其高程［10］，即任何一种DEM插值方法，待插点高程Z0都是已知点高程向量Z：Z1，Z2，Z3…，Zn（其中n为已知点个数）的函数，具体可由式（1）所示：

式（1）为DEM统一插值模型，Z0为待插点，i为参与内插点Zi的点数，qi为分配给参与插值点i的权重。实际上，不同的内插方法区别在于权重分配方式的不同。

Matlab中linear、cubic和nearest 3种内插基本模型都是基于三角形，是按Delaunay方法先找出内插点周围的3个点，构成包围内插点的三角形，再应用内插模型进行插值。

2 实例验证与MATLAB内插

2.1 实验区概况

本内插实验选取湖南西部山区的SRTM DEM数据，位于北纬28°东经110°的区块，文件名为N28E110.hgt.zip，大小约为2.75MB，具体位置为湖南省湘西土家族苗族自治州泸溪县六一村南300m的地方（图3）。实验区域是丘陵山地，该区域的SRTM DEM空白较多，用MATLAB读取该区域SRTM DEM文件，显示图形如图4所示，高海拔区黑点为高程空白区域。

图3 实验区域位置图

使用Global Mapper 10打开实验数据显示的结果如图5所示，可以清楚地观察到空白区域主要是山脊、河流。使用Global Mapper 10选取空白区块的一条剖面线，如图6所示，可见空白对应处出现断裂。

图4 DEM MATLAB显示

图5 DEM Global Mapper显示

图6 空白区域剖面线

2.2 MATLAB内插及数据处理流程

该hgt文件包含1弧度的区域，共有1201× 1201即144万多个数据点，SRTM DEM空白点是用－32768表示，其最大高程为1410，排除空白点的－32768后，最小值为42。

普通PC机在处理这样数量级的数据时内存占用过大，内插速度较慢，本实验将数据分为3段，再逐段内插处理，最后段合并，这样MATLAB主程序需要进行一个大循环，即每次从SRTM数据中截取一块进行内插，完成内插后再合并到另一个变量中。具体流程如图7所示。

图7 数据处理流程图

2.3 内插结果与精度分析

对SRTM DEM进行分段，采用linear、cubic和nearest内插函数，将实验数据N28E110.hgt进行插值处理，得到linear.hgt、cubic.hgt和nearest.hgt 3个文件。用global mapper软件显示这3个DEM文件和原始N28E110.hgt文件，具体如图8所示。

用global mapper软件导出同一区域的srtm v41.xyz和AST.xyz数据，基于这两种数据经过官方的外部数据修复和内插填补处理［1，3，5］，具有较高的精度。将SRTM DEM V4＿1数据和ASTGTM GDEM数据与图8数据进行对比，再使用MATLAB分别计算linear.hgt、cubic.hgt、nearest.hgt、N28E110.hgt、srtmv41.xyz以及AST.xyz的最大值、最小值、平均值、标准差，并分别以srtmv41.xyz和AST.xyz为基准计算中误差。具体结果如表2、表3所示。

图8 linear.hgt、cubic.hgt、nearest.hgt及N28E110.hgt图

表2 以srtm v41为基准的统计数据

表3 以ASTGTM GDEM数据为基准的统计数据

从表2、表3看出，由于SRTM DEM数据中存在空白单元，原始数据中高程最小值为－32768，内插前的高程平均值偏小，其标准差和中误差均远高于内插处理后数据。而3种内插方法填补空白后的数据平均值、标准差和中误差几乎一致，与srtm v41和ASTGTM GDEM的平均值、标准差与标准数据比较接近，三次内插法的标准差最小。可见，通过内插，数据质量得到了很大提升。

根据以上插值函数，将内插填补后的数据生成等高线，可得原始等高线与内插等高线对比图，如图9所示。

图9 原始等高线与内插等高线对比图

通常，可以通过对比图9中等高线细节部分的光滑程度、自然程度，来判断插值方法的好坏。由图9（a）看出原始数据中存在空白空洞，导致等高线图中出现了空白的块；图9（b）和图9（c）比较接近，但是详细比较，三次内插结果更加光滑；如图9（d）中等高线出现了不正常的叠加现象，可见，最邻近插值法不能体现细微的变化，内插后的数据不够连续，而三次插值法内插得到的等高线最光滑最自然。

3 结束语

SRTM DEM空白单元大多分布在海拔较高的山区，这些地方外部数据缺乏且现势性也较差，利用其周围数据直接内插是最简便的途径。本文针对SRTM DEM空白处的填补，利用MATLAB编写线性插值、三次插值、最邻近插值代码进行空白填补实验，并对比分析3种插值法填补的效果，得出线性插值法计算量比较大，插值后的数据较光滑；最近邻插值法输出高程值就等于距离它映射到的位置最近点的高程值，算法简单，但难以体现数据中细微的变化，插值后数据不连续；三次插值法内插得到的等高线最光滑最自然，内插填补数据质量与官方最新SRTM v4＿1数据非常接近，也与ASTGTM GDEM结果较为接近，得出三次插值法是一种山区较合适的DEM内插方法。但本文主要利用SRTM DEM周围数据插值，插值填补法的精度有待结合外部数据进一步提高。

［1］ CGIAR－CSI.SRTM 90mDEM Digital Elevation Database［EB／OL］.http：／／srtm.csi.cgiar.org，2012.

［2］ DENKER H.Evaluation of SRTM3and GOTOP30terrain data in Germany［D］.Proceeding of GGSM 2004，IAG，Porto，Portugal，2004.

［3］ 陈俊勇.对SRTM3和GTOPO30地形数据质量的评估［J］.武汉大学学报，2005，30（11）：941－944.

［4］ 阚瑷珂，朱利东，张瑞军，等.基于数据融合的SRTM数据空洞填补方法［J］.地理空间信息，2007，（3）：67－69.

［5］ 左美蓉.SRTM高程数据及其应用研究［D］.长沙：中南大学，2009.

［6］ 陈本富，王贵武，沈慧，等.基于Matlab的数据处理方法在GPS高程拟合中的应用［J］.昆明理工大学学报（理工版），2009，（5）：7－10.

［7］ 刘卫国.MATLAB程序设计教程［M］.北京：中国水利水电出版社，2006.

［8］ 徐建华.现代地理学中的数学方法［M］.北京：高等教育出版社，2002.

［9］ 张卡，盛业华，张书毕.MATLAB在测绘领域中的应用［J］.矿山测量，2004，17（1）：28－31.

［10］ 任志峰.DEM内插评价模型与应用系统开发研究［D］.南京：南京师范大学，2008.

［11］ 游松财，孙朝阳.中国区域SRTM90m数字高程数据空值区域的填补方法比较［J］.地理科学进展，2005，24（6）：88－92.

［12］ 詹蕾.SRTM DEM的精度评价及其适用性研究［D］.南京：南京师范大学，2008，23－24.

［13］ DEIACO S，MYERSD E，POSAD N.On separable space－time covariance models：Some parametric families［J］.Mathematical Geology，2002，（34）：23－42.

［14］ JENSONS K，DOMINGUE J.Extracting topographic structure from digital elevation data for geographical information system analysis［J］.Photogrammetric Engineering and Remote Sensing，1988，54（11）：1593－1600.

［15］ MA C.Families of spatio－temporal stationary covariance models［J］.J Stat Plan Infer，2003，（116）：489－501.

［16］ MA C.Spatio temporal covariance functions generated by mixtures［J］.Mathematical Geology，2002，34（8）：965－975.

［17］ 胡海，游涟，胡鹏，等.数字高程模型内插方法的分析和选择［J］.武汉大学学报（信息科学版），2011，（1）：86－89.

［18］ 兰燕，王明华，刘珊红，等.逐点内插法建立DEM的研究［J］.测绘科学，2009，34（1）：214－216.

［19］ 李胤，杨武年，杨容浩，等.基于移动曲面拟合算法和加权平均算法的DEM内插算法改进［J］.测绘，2010，33（4）：26－29.

［20］ 李志林，朱庆.数字高程模型［M］.武汉：武汉测绘科技大学出版社，2000.

［21］ 张蕾，陈晓宏.珠江三角洲网河区水位空间插值的Kriging方法［J］.中山大学学报（自然科学版），2004，（5）：112－114.

［22］ 周启鸣，刘学军.数字地形分析［M］.北京：科学出版社，2006.

［23］ BRUNRTTI M，MAUGERI M，NANNI T.High－resolution temperature climatology for Italy：Interpolation method intercomparison［J］.International Journal of Climatology，2014，34（4）：1278－1296.

［24］ AHN J S，CHUNG W J，JUNG C D.Realization of orientation interpolation of 6－axis articulated robot using quaternion［J］.Journal of Central South University，2012（19）：3407－3414.

［25］ NADAV C，AMIR B.A cartesian non－uniform grid interpolation method for fast field evaluation on elongated domains［J］.International Journal of Numerical Modelling：Electronic Networks，Devices and Fields，2012，25（5－6）：645－655.

SRTM DEM Voids Interpolation Method Based on MATLAB

LONG Si－chun1，2，3，ZHOU Wei2，WEN Jia－sheng2，CHEN Peng－qi3
（1.Hunan Key Laboratory of Coal Resources Clean－utilization and Mine Environment Protection，Hunan University of Science and Technology，Xiangtan411201；2.Institute of Geomatics and Geformation Monitoring，Hunan University of Science and Technology，Xiangtan411201；3.Mining and Safety Engineering，Hunan University of Science and Technology，Xiangtan411201）

SRTM DEM，which is free to the public，contains a lot of blank cells.This paper introduced the characteristics of SRTM DEM data and the research status of the voids filling at present.Taking advantage of object oriented programming features of MATLAB，we wrote codes which contain the calculation methods of linear interpolation，third order polynomial interpolation，nearest neighbor interpolation to fill SRTM DEM blank cells for Hunan west hills，and made comparative analysis of the processing results of three methods.We found that the third order polynomial interpolation method works best and produces the most smooth contours，which offers a new reference in selection of interpolation method of SRTM DEM for similar DEM areas.

MATLAB；SRTM DEM；triangle－based linear interpolation；third order polynomial interpolation；nearest neighbor interpolation

10.3969／j.issn.1000－3177.2015.04.004

TP31

A

1000－3177（2015）140－0020－05

2014－07－29

2014－10－09

国家自然科学基金（41474014）；湖南省教育厅科研重点项目（15A060）；大地测量与地球动力学国家重点实验室基金（SKLGED2014－5－3－E）；桂科能基金（1207115－21）；煤炭资源与环保湖南省重点实验室基金（E21221）。

龙四春（1975—），男，副教授，博士后，研究方向为雷达遥感与大地测量。

E－mail：sclong＠hnust.edu.cn