雷達地形測繪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