幾種插值方法在構建南極冰蓋DEM中的比較*

楊元德 熊云琪 王澤民 鄂棟臣

(武漢大學中國南極測繪研究中心，武漢 400079)

幾種插值方法在構建南極冰蓋DEM中的比較*

楊元德 熊云琪 王澤民 鄂棟臣

(武漢大學中國南極測繪研究中心，武漢 400079)

在南極冰蓋選取不同實驗區域對常用的六種插值方法進行比較分析，采用均方根預測誤差和源數據均方根等評價指標對不同插值方法進行評價，結果表明克里金插值方法更適用于構建南極冰蓋DEM。

南極冰蓋;DEM;插值方法;均方根;克里金插值方法

1 引言

南極冰蓋數字高程模型(DEM，Digital Elevation Model)在南極研究中發揮著重要的作用。它可用于確定分冰嶺、冰流盆地、接地線的位置、計算冰流的大小及方向、平衡速度及底部剪應力等;它也是估算冰面溫度、降水、下降風的大小和方向的重要參數;精確的地形數據對于利用遙感方法進行南極測圖研究也有著十分重要的意義;南極冰蓋DEM與冰厚數據相結合，能夠計算冰體變形的速率及應變;此外，精確的DEM也是構建南極冰下地形模型的基礎［1－6］。表 1給出了基于測高衛星的南極冰蓋DEM的發展概況。

表1 南極冰蓋DEM的發展Tab.1 Development of Antarctic ice sheet DEM

測高衛星為南極科學研究提供了豐富的源數據，基于這些衛星測高數據構建南極冰蓋DEM時，不同研究者采用了不同的數據插值方法，因此數據插值效果直接影響南極冰蓋DEM的精度。本文主要在南極冰蓋選取不同地形特征的實驗區，通過比較幾種常用插值方法的插值效果，選擇最適宜構建南極冰蓋DEM的插值方法。

2 數據插值方法簡介

本文采用的數據插值方法由美國GOLDEN軟件公司的Surfer10.1軟件提供，且主要對六種插值方法進行比較分析。

反距離加權插值法的插值原理是待插值點鄰域內已知離散點屬性值的加權平均，權的大小與待插點的鄰域內離散點之間的距離有關，是距離k(0≤k≤2，k一般取2)次方的倒數。即:

克里金插值方法是在滿足無偏估計和最小方差條件的前提下求得估計值。設區域化變量f(x)滿足二階平穩條件，則待插值點P的估計值為

其中fi是n個已知點的函數值，wi是n個已知點的權系數。其中μ為拉格朗日算子，γ(xj－xi)為已知點間的變差函數值，γ(xp－xi)為已知點與待插值點間的變差函數值。

最小曲率法用公式表示為:

傳統的多項式插值方法在Surfer10.1中有三種形式:

最近鄰點插值法又稱泰森多邊形方法，其隱含假設條件是任一網格點P(x，y)的屬性值都使用距它最近位置點的屬性值，用每一個網格節點的最鄰點值作為該點的節點值。假若有一塊總面積為s的區域，共有N個網格點n(n=1，2，…，N)。各個網格點的屬性值分別為qi(i=1，2，…，N)，以各網格點之間連線的垂直平分線，把區域劃分為若干個多邊形，然后以各個多邊形的面積si(i=1，2，…，N)為權系數，計算各網格點的加權平均值，并將其作為區域的平均屬性值，即［8］:

假設三角片的3個頂點是A、B、C，其觀測值分別為f(A)、f(B)、f(C)，任何點相對于該三角形有唯一的重心，其坐標為:

線性插值可寫成

由式(6)和式(7)，點(x，f(x))可以表示為{(A，f(A));(B，f(B));(C，f(C))}的線性凸組合，也就是說這是一張平面;當0＜u、v、w＜1時，x位于三角形內部，(x，f(x))形成一個三角片［9］。由以上原理出發，需要求得系數u、v、w，本文采用三角形的3個頂點信息分別獲得系數值。

3 插值方法評價指標與實驗

3.1 評價指標

插值精度的評定指標概括起來主要分為:1)對于采樣點計算的均方根預測誤差(RMSPE，Root Mean Square Predictive Error)，其值越小說明插值結果越接近真實值，結果越可靠;2)計算源數據得到的殘差均方根(RMS)，源數據殘差均方根用于評價網格數據與源數據的一致性。因此本文主要以這二個指標對插值結果進行評價，并輔之以殘差絕對值最大、最小、平均值等參數。

均方根預測誤差的計算公式為:

式中k表示用于估計網格數據精度的樣點數，zj表示第j(1≤j≤k)個樣點插值后高程，z'j表示第j個樣點原始高程。

源數據殘差均方根計算公式為:

式中n表示整個實驗區的用于插值的源數據點數，zi表示第i個插值后高程，z'i表示第i(1≤i≤n)個原始觀測高程。

3.2 實驗數據和實驗區

實驗采用的ERS-1數據由歐洲太空局網站(http://www.esa.int/esaCP/index.html)下載。其數據采集時間為ERS-1在Phase E和Phase F工作段采集。ERS-1任務時間段如表2所示。

表2 ERS-1任務時間段表Tab.2 ERS-1 mission phase

為選擇最佳插值方法，選擇的三個地區的地形特征為:橫斷山脈的地區地形較為復雜，有些區域的高程變化比較急劇，是典型的山區地形;南極半島地區的地形非常復雜，基巖起伏不平，有大量山峰和山脈，海岸曲折，近海島嶼很多地勢陡峭;中山站至DomeA地區的地形比較簡單，變化較緩，地勢平坦。如表3所示。

表3 ERS-1數據統計的實驗區信息Tab.3 Statistics of ERS-1 data within experimental areas

3.3 實驗計算及結論

利用Surfer10.1軟件對實驗區進行隨機均勻采樣，在每個實驗區均勻隨機取約占總點數約5%的采樣點。橫斷山脈、南極半島、中山站至Dome A三個區域的采樣點分別為35 000、35 000和75 000。然后采用有效性評價指標對插值方法進行有效性評價。結果如表4～6。

表4 橫斷山脈地區不同插值方法評價指標計算結果(單位:m)Tab.4 Comparison of interpolation methods in Trans-Antarctic mountains(unit:m)

表5 南極半島地區不同插值方法評價指標計結果(單位:m)Tab.5 Comparison of interpolation methods in Antarctic peninsula(unit:m)

表6 中山站-Dome A地區不同插值方法評價指標計算結果(單位:m)Tab.6 Comparison of interpolation methods in the transit from Zhongshan station to Dome A area(unit:m)

對表4～6分析可以得到:

由RMSPE和殘差RMS統計結果可知，最近臨點方法在表5、6中較小，但在表4中較大，可見最近臨點方法依賴于數據或地形，其插值效果不穩定。而對于線性插值三角網、反距離加權插值、局部多項式插值在三個實驗區，無論是殘差均值平均值，還是RMSPE和RMS都比較大，這三種插值方法的插值效果較差。而克里金和最小曲率插值算法的RMSPE和RMS最小，說明這兩種方法的插值效果最好。為了對克里金插值算法和最小曲率插值效果進一步比較，圖1給出了兩種方法插值的中山站-Dome A地區等高線圖。由圖1可以看出，最小曲率的等高線上方出現明顯的錯誤。這是因為最小曲率主要考慮曲面的光滑性，使得插值結果容易失真，繪出的等值線與實際相差較大。綜上分析，在南極地區，插值方法效果最好的是克里金插值，其余的幾種插值方法都存在不足。

圖1 中山站-Dome A地區等高線圖Fig.1 Contour map of the transit from Zhongshan station to Dome A area

對比表4～6可以發現表6中殘差絕對值均值、RMS和RMSPE明顯小于其他兩個表格，該結果表明了地形平坦、坡度變化較緩的中山站-Dome A地區插值精度高于地形復雜、地勢陡峭的橫斷山脈和南極半島，即地形因素對于插值結果的精度具有很大影響，平坦地區插值精度較高，而坡度較大、地形劇烈變化地區插值精度較低。

4 結束語

對比分析結果表明，在眾多插值方法中，克里金插值方法有很好的適應性，對于不同復雜程度的地形均能實現有效的空間插值，精度較高，因此建議采用克里金算法構建南極冰蓋DEM。從在實驗地區的插值效果來看，大部分方法均適用于地勢相對平坦的區域，如中山站至DOME A區域，但對于坡度較大、地形劇烈變化地區，有些插值方法出現了較大誤差精度也較低，如南極半島。因此在選擇內插方法時，尤其是復雜區域，需要慎重。本文沒有對各方法的插值效率、數據分布均勻與否條件下的插值表現進行比較，這也是今后的重點研究方向。

1 張勝凱，等.南極數字高程模型研究進展［J］.極地研究，2006，18(4):301 －309(Zhang S k，et al.Progress in Antarctic DEM［J］.Polar Research，2006，18(4):301 －309)

2 Budd W F，Jenssen D and Smith I N.A three-dimensional time-dependent model of the West Antarctic ice sheet［J］.Ann Glaciol.，1984，5:29 －36.

3 Zwally H J，et al.Surface elevation contours of Greenland and Antarctic ice sheet［J］.J Geophys Res.，1983，88:1 589 －1 598.

4 Zwally H J，et al.Ice measurements by GEOSAT radar altimetry［R］.The Navy GEOSAT Mission，Johns Hopkins APL Tech Dig，1987，8(2):251－254.

5 Zwally H J，et al.Ice sheet topography，slopes and flow directions from ERS altimetry［A］.Florence，Italy:The 3rd ERS Symposium，Euro Space Agency［C］.1997.

6 吳宗敏.散亂數據擬合的模型、方法和理論［M］.北京:科學出版社，2008.(Wu Z M.Model，method and theories for scattered data fitting［M］.Beijing:Science Press，2008)

7 Smith W H F and Wessel P.Gridding with continuous curvature splines in tension［J］.Geophysics，1990，55，3:293 －305.

8 顧春雷，楊漾，朱志春.幾種建立DEM模型插值方法精度的交叉驗證［J］.測繪與空間地理信息，2010，34(5):99 －102.(Gu C L，Yang Y and Zhu Z C.Accuracy crossvalidation of several interpolation methods for DEM model［J］.Surveying＆ Spatial Information，2010，34(5):99－102)

9 武俊紅，汪云甲.基于Surfer的煤礦等值線空間插值方法有效性評價［J］.中國礦業，2007，16(1):108－120.(Wu J H and Wang Y J.Effectiveness evaluation of interpolation methods in coal contour based on Surfer［J］.China Mining，2007，16(1):108－120)

COMPARISON OF INTERPOLATION METHODS IN CONSTRUCTING ANTARCTIC ICE SHEET DEM

Yang Yuande，Xiong Yunqi，Wang Zeming and E Dongchen
(Chinese Antarctic Center of Surveying and Mapping，Wuhan University，Wuhan430079)

The paper compared and discussed six commonly used interpolation methods.The interpolation methods were evaluated using the parameters such as the root mean square prediction error and the source data root mean square，within different experimental areas in the Antarctica.The result showed that the Kriging interpolation method is more suitable for constructing the Antarctic ice sheet DEM.

Antarctic ice sheet;Digital Elevation Model(DEM);interpolation method;Root Mean Square(RMS);Kriging interpolation method

P207

A

1671-5942(2013)05-0063-04

2013-02-14

國家自然科學基金(41106163);國家海洋局極地科學重點實驗室開發研究基金(KP201102);南北極環境綜合考察與評估專項(CHINARE2013-03-03，CHINARE2013-01-03);中央高校基本科研業務費(121047);國家863計劃項目(2012AA12A304);國家測繪地理信息局科技項目(極地測繪技術試驗)

楊元德，男，副教授，主要從事衛星測高與衛星重力應用研究.E－mail:yuandeyang@whu.edu.cn