雷達地形測繪DEM空洞插補方法研究

龍四春，周威，文佳勝，陳鵬琦

（1.湖南科技大學煤炭資源清潔利用與礦山環(huán)境保護湖南省重點實驗室，湖南湘潭411201；2.湖南科技大學測量工程與形變監(jiān)測研究所，湖南湘潭411201；3.湖南科技大學能源學院，湖南湘潭411201）

雷達地形測繪DEM空洞插補方法研究

龍四春1，2，3，周威2，文佳勝2，陳鵬琦3

（1.湖南科技大學煤炭資源清潔利用與礦山環(huán)境保護湖南省重點實驗室，湖南湘潭411201；2.湖南科技大學測量工程與形變監(jiān)測研究所，湖南湘潭411201；3.湖南科技大學能源學院，湖南湘潭411201）

免費公開的SRTM DEM數(shù)據(jù)用圖廣泛，但存在大量的空洞。針對SRTM DEM空白處填補現(xiàn)狀，利用MATLAB編程特性，編寫線性插值、三次插值、最鄰近插值的MATLAB代碼，選擇湖南西部丘陵山地進行空白填補實驗。通過對比分析其最大值、最小值、標準差、中誤差，得出三次插值法填補的效果最佳，生產(chǎn)的等高線最光滑。該方法與結(jié)論能為類似地區(qū)DEM插補方法選擇提供參考。

MATLAB；SRTM DEM；線性插值；三次插值；最鄰近插值

0 引 言

航天飛機雷達地形測繪使命（Shuttle Radar Topography Mission，SRTM）數(shù)據(jù)是制作地形圖與地形分析的寶貴資源。但由于雷達側(cè)視的幾何特征、航天飛機軌道設計、信號干擾、雷達陰影或回波滯后等因素的影響，導致部分地區(qū)出現(xiàn)數(shù)據(jù)空洞，尤其在水體和高山峽谷地區(qū)。因此，要增強SRTM數(shù)據(jù)的實用性，必須對其數(shù)據(jù)空洞進行填補［1－5］。MATLAB作為專門的數(shù)值計算軟件，具有強大、專業(yè)的數(shù)據(jù)計算、分析和可視化等功能［6－9］，將之應用SRTM DEM空白插補可以提高效率。

國內(nèi)已有許多學者嘗試過各種方法對SRTM高程數(shù)據(jù)的空值區(qū)域進行了填補［1－5，10－12］，但由于填補的數(shù)據(jù)地理范圍、操作者所掌握的技術以及其他輔助性的數(shù)據(jù)不同，從而在過程和結(jié)果上就有所差異。2009年，左美蓉用等高線內(nèi)插生成分辨率與SRTM相同的DEM，然后將內(nèi)插出來的高程值取整來填補SRTM高程數(shù)據(jù)的空值區(qū)域，但填補后的數(shù)據(jù)顯得不連續(xù)，內(nèi)插出來的值與SRTM數(shù)據(jù)在空值邊界的值有較大差距［5］。2012年CSI（Consortium for Spatial Information）將等高線和SRTM數(shù)據(jù)整合起來再內(nèi)插，由于中間只有等高線上的高程值，而且相距較遠，對于大范圍空值區(qū)域，內(nèi)插結(jié)果很難消除帶狀效果［1］。同時，國外很多其他機構(gòu)和個人也在研究DEM生成的填補方法和填補工具［1－7，10，13－25］，如表1所示。

表1 SRTM數(shù)據(jù)空洞填補工具

1 SRTM DEM及MATLAB內(nèi)插原理

1.1 SRTM DEM數(shù)據(jù)

SRTM數(shù)據(jù)每經(jīng)緯度方格提供一個文件，精度有1arc－second和3arc－seconds兩種，稱作SRTM1和SRTM3，或者稱作30m和90m數(shù)據(jù)，SRTM1的文件里面包含3601×3601個采樣點的高度數(shù)據(jù)，SRTM3的文件里面包含1201×1201個采樣點的高度數(shù)據(jù)，每一個小像元代表邊長約90m的正方形區(qū)域，如圖1所示。

圖1 SRTM DEM文件示意圖

圖2 SRTM DEM空白示意圖

由于天線桿和姿態(tài)的量測精度、記時誤差、多路徑效應、相位量測誤差及雷達的熱噪聲等的影響，在SRTM中會出現(xiàn)一些空白值，如圖2中紅色方框所示。

1.2 SRTM DEM數(shù)據(jù)MATLAB內(nèi)插原理

通常，內(nèi)插貫穿在DEM的生產(chǎn)、質(zhì)量控制、精度評定和分析應用等各環(huán)節(jié)。DEM內(nèi)插實質(zhì)是根據(jù)分布在內(nèi)插點周圍的采樣點高程值求出待定點的高程值過程，根據(jù)內(nèi)插點分布范圍，可以將內(nèi)插分為整體內(nèi)插、分塊內(nèi)插和逐點內(nèi)插三類［10］。DEM內(nèi)插方法的選擇，要考慮諸多因素，不僅要滿足DEM內(nèi)插精度要求，還要盡可能顧及計算效率，也要考慮地面復雜函數(shù)和已知數(shù)據(jù)特點。

MATLAB語言是一個基于矩陣和矢量的高級語言，簡單易學，又具有面向?qū)ο蟮木幊烫攸c，編程效率高；有眾多的應用工具箱，具備很強的開放性，除內(nèi)部函數(shù)外，用戶可通過對源文件的修改或加入自己編寫的程序語句去構(gòu)成新的專用工具箱。

DEM內(nèi)插算法的空間相似性反映在由未知點附近已知點高程的加權平均值來確定其高程［10］，即任何一種DEM插值方法，待插點高程Z0都是已知點高程向量Z：Z1，Z2，Z3…，Zn（其中n為已知點個數(shù)）的函數(shù)，具體可由式（1）所示：

式（1）為DEM統(tǒng)一插值模型，Z0為待插點，i為參與內(nèi)插點Zi的點數(shù)，qi為分配給參與插值點i的權重。實際上，不同的內(nèi)插方法區(qū)別在于權重分配方式的不同。

Matlab中l(wèi)inear、cubic和nearest 3種內(nèi)插基本模型都是基于三角形，是按Delaunay方法先找出內(nèi)插點周圍的3個點，構(gòu)成包圍內(nèi)插點的三角形，再應用內(nèi)插模型進行插值。

2 實例驗證與MATLAB內(nèi)插

2.1 實驗區(qū)概況

本內(nèi)插實驗選取湖南西部山區(qū)的SRTM DEM數(shù)據(jù)，位于北緯28°東經(jīng)110°的區(qū)塊，文件名為N28E110.hgt.zip，大小約為2.75MB，具體位置為湖南省湘西土家族苗族自治州瀘溪縣六一村南300m的地方（圖3）。實驗區(qū)域是丘陵山地，該區(qū)域的SRTM DEM空白較多，用MATLAB讀取該區(qū)域SRTM DEM文件，顯示圖形如圖4所示，高海拔區(qū)黑點為高程空白區(qū)域。

圖3 實驗區(qū)域位置圖

使用Global Mapper 10打開實驗數(shù)據(jù)顯示的結(jié)果如圖5所示，可以清楚地觀察到空白區(qū)域主要是山脊、河流。使用Global Mapper 10選取空白區(qū)塊的一條剖面線，如圖6所示，可見空白對應處出現(xiàn)斷裂。

圖4 DEM MATLAB顯示

圖5 DEM Global Mapper顯示

圖6 空白區(qū)域剖面線

2.2 MATLAB內(nèi)插及數(shù)據(jù)處理流程

該hgt文件包含1弧度的區(qū)域，共有1201× 1201即144萬多個數(shù)據(jù)點，SRTM DEM空白點是用－32768表示，其最大高程為1410，排除空白點的－32768后，最小值為42。

普通PC機在處理這樣數(shù)量級的數(shù)據(jù)時內(nèi)存占用過大，內(nèi)插速度較慢，本實驗將數(shù)據(jù)分為3段，再逐段內(nèi)插處理，最后段合并，這樣MATLAB主程序需要進行一個大循環(huán)，即每次從SRTM數(shù)據(jù)中截取一塊進行內(nèi)插，完成內(nèi)插后再合并到另一個變量中。具體流程如圖7所示。

圖7 數(shù)據(jù)處理流程圖

2.3 內(nèi)插結(jié)果與精度分析

對SRTM DEM進行分段，采用linear、cubic和nearest內(nèi)插函數(shù)，將實驗數(shù)據(jù)N28E110.hgt進行插值處理，得到linear.hgt、cubic.hgt和nearest.hgt 3個文件。用global mapper軟件顯示這3個DEM文件和原始N28E110.hgt文件，具體如圖8所示。

用global mapper軟件導出同一區(qū)域的srtm v41.xyz和AST.xyz數(shù)據(jù)，基于這兩種數(shù)據(jù)經(jīng)過官方的外部數(shù)據(jù)修復和內(nèi)插填補處理［1，3，5］，具有較高的精度。將SRTM DEM V4＿1數(shù)據(jù)和ASTGTM GDEM數(shù)據(jù)與圖8數(shù)據(jù)進行對比，再使用MATLAB分別計算linear.hgt、cubic.hgt、nearest.hgt、N28E110.hgt、srtmv41.xyz以及AST.xyz的最大值、最小值、平均值、標準差，并分別以srtmv41.xyz和AST.xyz為基準計算中誤差。具體結(jié)果如表2、表3所示。

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

表2 以srtm v41為基準的統(tǒng)計數(shù)據(jù)

表3 以ASTGTM GDEM數(shù)據(jù)為基準的統(tǒng)計數(shù)據(jù)

從表2、表3看出，由于SRTM DEM數(shù)據(jù)中存在空白單元，原始數(shù)據(jù)中高程最小值為－32768，內(nèi)插前的高程平均值偏小，其標準差和中誤差均遠高于內(nèi)插處理后數(shù)據(jù)。而3種內(nèi)插方法填補空白后的數(shù)據(jù)平均值、標準差和中誤差幾乎一致，與srtm v41和ASTGTM GDEM的平均值、標準差與標準數(shù)據(jù)比較接近，三次內(nèi)插法的標準差最小?？梢?，通過內(nèi)插，數(shù)據(jù)質(zhì)量得到了很大提升。

根據(jù)以上插值函數(shù)，將內(nèi)插填補后的數(shù)據(jù)生成等高線，可得原始等高線與內(nèi)插等高線對比圖，如圖9所示。

圖9 原始等高線與內(nèi)插等高線對比圖

通常，可以通過對比圖9中等高線細節(jié)部分的光滑程度、自然程度，來判斷插值方法的好壞。由圖9（a）看出原始數(shù)據(jù)中存在空白空洞，導致等高線圖中出現(xiàn)了空白的塊；圖9（b）和圖9（c）比較接近，但是詳細比較，三次內(nèi)插結(jié)果更加光滑；如圖9（d）中等高線出現(xiàn)了不正常的疊加現(xiàn)象，可見，最鄰近插值法不能體現(xiàn)細微的變化，內(nèi)插后的數(shù)據(jù)不夠連續(xù)，而三次插值法內(nèi)插得到的等高線最光滑最自然。

3 結(jié)束語

SRTM DEM空白單元大多分布在海拔較高的山區(qū)，這些地方外部數(shù)據(jù)缺乏且現(xiàn)勢性也較差，利用其周圍數(shù)據(jù)直接內(nèi)插是最簡便的途徑。本文針對SRTM DEM空白處的填補，利用MATLAB編寫線性插值、三次插值、最鄰近插值代碼進行空白填補實驗，并對比分析3種插值法填補的效果，得出線性插值法計算量比較大，插值后的數(shù)據(jù)較光滑；最近鄰插值法輸出高程值就等于距離它映射到的位置最近點的高程值，算法簡單，但難以體現(xiàn)數(shù)據(jù)中細微的變化，插值后數(shù)據(jù)不連續(xù)；三次插值法內(nèi)插得到的等高線最光滑最自然，內(nèi)插填補數(shù)據(jù)質(zhì)量與官方最新SRTM v4＿1數(shù)據(jù)非常接近，也與ASTGTM GDEM結(jié)果較為接近，得出三次插值法是一種山區(qū)較合適的DEM內(nèi)插方法。但本文主要利用SRTM DEM周圍數(shù)據(jù)插值，插值填補法的精度有待結(jié)合外部數(shù)據(jù)進一步提高。

［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地形數(shù)據(jù)質(zhì)量的評估［J］.武漢大學學報，2005，30（11）：941－944.

［4］ 闞璦珂，朱利東，張瑞軍，等.基于數(shù)據(jù)融合的SRTM數(shù)據(jù)空洞填補方法［J］.地理空間信息，2007，（3）：67－69.

［5］ 左美蓉.SRTM高程數(shù)據(jù)及其應用研究［D］.長沙：中南大學，2009.

［6］ 陳本富，王貴武，沈慧，等.基于Matlab的數(shù)據(jù)處理方法在GPS高程擬合中的應用［J］.昆明理工大學學報（理工版），2009，（5）：7－10.

［7］ 劉衛(wèi)國.MATLAB程序設計教程［M］.北京：中國水利水電出版社，2006.

［8］ 徐建華.現(xiàn)代地理學中的數(shù)學方法［M］.北京：高等教育出版社，2002.

［9］ 張卡，盛業(yè)華，張書畢.MATLAB在測繪領域中的應用［J］.礦山測量，2004，17（1）：28－31.

［10］ 任志峰.DEM內(nèi)插評價模型與應用系統(tǒng)開發(fā)研究［D］.南京：南京師范大學，2008.

［11］ 游松財，孫朝陽.中國區(qū)域SRTM90m數(shù)字高程數(shù)據(jù)空值區(qū)域的填補方法比較［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］ 胡海，游漣，胡鵬，等.數(shù)字高程模型內(nèi)插方法的分析和選擇［J］.武漢大學學報（信息科學版），2011，（1）：86－89.

［18］ 蘭燕，王明華，劉珊紅，等.逐點內(nèi)插法建立DEM的研究［J］.測繪科學，2009，34（1）：214－216.

［19］ 李胤，楊武年，楊容浩，等.基于移動曲面擬合算法和加權平均算法的DEM內(nèi)插算法改進［J］.測繪，2010，33（4）：26－29.

［20］ 李志林，朱慶.數(shù)字高程模型［M］.武漢：武漢測繪科技大學出版社，2000.

［21］ 張蕾，陳曉宏.珠江三角洲網(wǎng)河區(qū)水位空間插值的Kriging方法［J］.中山大學學報（自然科學版），2004，（5）：112－114.

［22］ 周啟鳴，劉學軍.數(shù)字地形分析［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）；煤炭資源與環(huán)保湖南省重點實驗室基金（E21221）。

龍四春（1975—），男，副教授，博士后，研究方向為雷達遙感與大地測量。

E－mail：sclong＠hnust.edu.cn