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

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

雷達(dá)地形測(cè)繪DEM空洞插補(bǔ)方法研究

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

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

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

MATLAB；SRTM DEM；線(xiàn)性插值；三次插值；最鄰近插值

0 引 言

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

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

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

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

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

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

圖1 SRTM DEM文件示意圖

圖2 SRTM DEM空白示意圖

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

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

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

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

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

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

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

2 實(shí)例驗(yàn)證與MATLAB內(nèi)插

2.1 實(shí)驗(yàn)區(qū)概況

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

圖3 實(shí)驗(yàn)區(qū)域位置圖

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

圖4 DEM MATLAB顯示

圖5 DEM Global Mapper顯示

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)束語(yǔ)

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

［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］ 陳俊勇.對(duì)SRTM3和GTOPO30地形數(shù)據(jù)質(zhì)量的評(píng)估［J］.武漢大學(xué)學(xué)報(bào)，2005，30（11）：941－944.

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

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

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

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

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

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

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

［11］ 游松財(cái)，孫朝陽(yáng).中國(guó)區(qū)域SRTM90m數(shù)字高程數(shù)據(jù)空值區(qū)域的填補(bǔ)方法比較［J］.地理科學(xué)進(jìn)展，2005，24（6）：88－92.

［12］ 詹蕾.SRTM DEM的精度評(píng)價(jià)及其適用性研究［D］.南京：南京師范大學(xué)，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］.武漢大學(xué)學(xué)報(bào)（信息科學(xué)版），2011，（1）：86－89.

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

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

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

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

［22］ 周啟鳴，劉學(xué)軍.數(shù)字地形分析［M］.北京：科學(xué)出版社，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

國(guó)家自然科學(xué)基金（41474014）；湖南省教育廳科研重點(diǎn)項(xiàng)目（15A060）；大地測(cè)量與地球動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室基金（SKLGED2014－5－3－E）；桂科能基金（1207115－21）；煤炭資源與環(huán)保湖南省重點(diǎn)實(shí)驗(yàn)室基金（E21221）。

龍四春（1975—），男，副教授，博士后，研究方向?yàn)槔走_(dá)遙感與大地測(cè)量。

E－mail：sclong＠hnust.edu.cn