2.4　常用DEM插值算法

縱觀DEM插值算法研究的發展歷程，反距離加權插值算法、改進謝別德插值算法、徑向基函數插值算法、克里格插值算法是常用的DEM插值算法。

? 2.4.1　反距離加權插值算法

反距離加權插值算法（Inverse Distance Weighted，IDW）最早是由氣象學家和地質工作者提出的（王建等，2004），是空間數據插值最常見的算法之一。

反距離加權插值算法基于相近、相似的原理（盧華興，2008），每個采樣點都對插值點有一定的影響，即權重。權重隨著采樣點和插值點之間距離的增加而減弱，距離插值點越近的采樣點的權重越大；當采樣點在距離插值點一定范圍以外時，權重可以忽略不計。任意插值點的值是各采樣點的權重之和（王家華等，1999），如式（2.3）所示。

式中，zp為插值點的高程值，λi為第i個點的權重，di為第i個采樣點到插值點的距離，d-u為距離衰減函數，冪指數u具有隨著距離的增加減小其他位置影響的作用（de Smith et al.，2007）。當u=0時，距離沒有影響；當u=1時，距離的影響是線性的；當u?1時，快速地減小遙遠位置的影響。冪指數u通常取值為1或2（Lam，1983），但是大多數學者認為，冪指數采用2將取得更好的實驗效果（Declercq，1996）。

反距離加權插值算法的計算易受采樣點集群的影響，導致在采樣點附近局部出現明顯的隆起或凹陷的“牛眼”效應（Johns，1998）。

另外，由于反距離加權插值算法是一種精確性插值算法，因此插值生成的最大值和最小值只會出現在采樣點處。這會直接導致出現山頂高程被降低、山谷高程被抬高的局部細節湮沒。

基于反距離加權插值算法的缺點，許多學者提出了反距離加權插值算法的多種改進形式，用于克服上述缺點。

（1）給距離衰減函數增加一個平滑參數：一個微小的距離增量t，即，這樣可能導致平滑，而不是精確插值。

（2）給距離衰減函數增加一個調和參數：基于最遠點的距離R調整權重值，即

（3）使用其他的距離衰減函數，如高斯函數，即λi=e-（d/m）2。

? 2.4.2　改進謝別德插值算法

謝別德插值算法（Modified Shepards Method，SPD）由南非地質學家 Shepard最早提出，本質上是一種標準的導數距離加權過程（王金玲等，2010），權函數如式（2.5）所示。

式中，wi為權重，di為第i點距待插值點的距離，r為調整距離。

改進謝別德插值算法一般存在兩種變化形式。

一是基于最遠點的距離（在整個數據集中或在給定搜索半徑范圍內）調整權重。假設最遠距離為r，那么修正的距離倒數加權函數如式（2.6）所示。

二是使用擬合的局部二次多項式調整權重，即參與距離倒數加權函數的高程值并不是原始采樣點的高程值，而是使用擬合的局部二次多項式修正的高程值，如式（2.7）所示。

式中，zj為待插值點的高程值；，為插值點至采樣點的距離；δ為平滑因子，當δ=0時為精確性插值，當δ≠0時為非精確性插值；Qi為二次多項式函數；u為權指數。

? 2.4.3　徑向基函數插值算法

徑向基函數插值（Radial Basis Functions，RBF）算法是一系列用于精確插值算子的統稱（de Smith et al.，2007）。它來源于Hardy的多面函數法，其插值原理是任何一個表面都可以使用多個曲面的線性組合逼近（盧華興，2008）。在多數情況下，徑向基函數插值算法與地統計插值算法相似，但具有不需要分析半變異函數模型的優點，而且不需要有關采樣點的任何假設（除了非共線性）。

在通常情況下，徑向基函數插值算法可以表述為兩個部分之和（Mitasova and Mitas，1993），即

式中，zp為插值點的高程值；λi為第i個點的權重；di為第i個采樣點到插值點的距離；φ（di）為徑向基函數，它代表第j個核函數對多層疊加曲面的貢獻；fj（x）為“趨勢”函數，是次數小于m的基本多項式函數，由于fj（x）并不能提高插值的精度，因此在插值過程中不考慮“趨勢”函數的影響（de Smith et al.，2007）。徑向基函數插值算法的解算過程可以使用矩陣符號表示為如下步驟。

步驟1：計算源數據中所有（x，y）點對的點間距離構成的n×n階矩陣D，有

步驟2：對D中的每個矩陣值應用選擇的徑向基函數φ（·），從而產生一個新的矩陣Φ，有

步驟3：用單位列矢量和單位行矢量增大矩陣Φ，并且在位置（n+1，n+1）處插入零值，稱這個增廣矩陣為A，有

步驟4：計算從格網點p到用來創建D的每個源數據點間的距離構成的列矢量r，有

步驟5：將選擇的徑向基函數應用于r中的每個距離產生一個列矢量φ，然后生成一個n+1階的列矢量c，它由φ加上元素1構成，即

步驟6：計算矩陣積b=A-1c。這樣就給出了用于計算p點的估計值的n個權重，即

徑向基函數插值算法可以選用許多不同的徑向基函數（見表2.3）。

表中，c為光滑因子，一般由用戶指定。c的值取決于對插值結果產生重要影響的采樣點的數目、高程、空間分布等因素（Rippa，1999）。

對于如何確定c，沒有普遍認可的方法，但也有一些學者提出了各種方法。Hardy（1971）使用c=0.815d，其中，，di為第i個點到其最近鄰的距離。Franke（1982）使用替換d，其中，D是數據集最小外接圓的直徑，于是建議使用。Foley（1987）做出了和Franke類似的建議，不過使用數據集的最小外接矩形的邊長代替最小外接圓的直徑。Rippa（1999）提出了使用遞歸算法尋找使得插值表面全局誤差最小的參數c的方法。Aguilar等（2005）認為，在MQF插值算法和MLF插值算法中應當使用接近于零的光滑因子；在IMQF插值算法、NCSF插值算法、TPSF插值算法中則應當使用非常大的光滑因子，因為在IMQF插值算法、NCSF插值算法、TPSF插值算法中如果使用較小的光滑因子，將會產生顯著的數值不穩定性。

徑向基函數插值算法作為一種精確的插值算法，不同于局部多項式插值算法。局部多項式插值算法作為一種非精確的插值算法，并不要求表面經過所有的采樣點。徑向基函數插值算法和同為精確插值算法的反距離加權插值算法的不同之處在于，反距離加權插值算法不能計算出高于或低于采樣點的插值點的值，而徑向基函數插值算法則可以計算出高于或低于采樣點的插值點的值。

? 2.4.4　克里格插值算法

克里格插值算法也稱局部估計插值或空間局部插值，是地統計學的兩大主要內容之一（張景雄，2008）。地統計學源于20世紀50年代Krige在地質和采礦業方面的工作，1963年法國學者Matheron發表了專著《應用地質統計學》，提出了區域化變量理論，并且給出了地統計學的概念：以區域化變量理論為基礎，以變異函數為主要工具，研究在空間分布上既有隨機性又有結構性的自然現象的科學（侯景儒，1998）。

1.區域化變量

區域化變量是以空間采樣點x的三維直角坐標（xu，xv，xw）為自變量的隨機場函數Z（xu，xv，xw）=Z（x），在對其進行一次觀測后，就得到隨機場函數Z（x）的一個具體實現z（x）。在空間的每個點取某個確定的數值后，當由一個點移動到下一個點時，函數具體實現z（x）是變化的。

區域化變量具有隨機性和結構性的雙重特征。隨機性是指區域化變量在具體實現時表現出一定的不規則特征；結構性是指區域化變量在不同的空間方位具有某種程度的空間自相關性。

地形表面作為一個連續的隨機場表面，符合區域化變量的雙重特征，因此以區域化變量理論為基礎的“地統計學”在地形建模、空間分析等方面的應用方興未艾。

2.半變異函數

半變異函數是一種空間變量相關性的定量化描述模型。當空間采樣點在一維軸x上變化時，區域化變量在x和x+h處的值為Z（x）和Z（x+h），兩者之差的方差的一半定義為區域化變量在x軸上的半變異函數，記為γ（x，h）。

在二階平穩假設下，有

那么γ（x，h）可以改寫成

從式（2.17）可知，半變異函數依賴于兩個自變量x和h，當半變異函數γ（x，h）與位置x無關時，它僅依賴于分隔兩個采樣點之間的距離，那么γ（x，h）可以改寫成γ（h），即

在通常情況下，半變異函數值隨著采樣點間距的增加而增大，并在到達某一個間距值后趨于穩定（見圖2.4）。半變異函數具有3個重要的參數，分別是塊金值（Nugget）、基臺值（Sill）和變程（Range），它們表示區域化變量在一定尺度上的空間變異性和相關性。

（1）塊金值。根據半變異函數的定義，理論上當h=0時，半變異函數值應等于0。但是，由于采樣誤差等原因，即使兩個采樣點之間的距離h很小，其變量依然存在差異，這表示區域化變量在小于觀測尺度時的非連續性變異。

（2）基臺值。基臺值表示半變異函數隨著間距遞增到一定程度時出現的平穩值，即C0+C。其中，C稱為結構方差（或拱高），在數值上等于基臺值與塊金值之間的差值，代表由于樣本數據中存在空間相關性而引起的方差變化的范圍。

（3）變程。變程表示半變異函數達到基臺值時的距離，反映了空間采樣點的自相關距離尺度。在變程距離之內，空間上越近的點之間的相關性越大；當h大于變程時，空間采樣點之間不具備自相關性，除非半變異函數具有周期性變化特征。更為重要的是，變程表示了空間插值的極限距離，選擇在變程范圍內的采樣點參與插值才有意義（張仁鐸，2005）。

在實際插值估計中，由于空間采樣點是離散的，無法獲取半變異函數γ（h）的理論值，所以需要通過實驗方法獲得實驗半變異函數值γ*（h），即

式中，N（h）是近似相隔h的采樣點對的數目。

之后，根據實驗半變異函數值選擇合適的理論半變異函數模型，并且擬合半變異函數模型的基本參數。

理論半變異函數模型包括幾種簡單的模型，如圖2.5所示。

（1）線性模型（LINE）：

（2）球形模型（SPHERE）：

當C0=0、C=1時，球形模型稱為標準球形模型。

（3）指數模型（EXP）：

當C0=0、C=1時，指數模型稱為標準指數模型。

（4）高斯模型（GAUSS）：

當C0=0、C=1時，高斯模型稱為標準高斯模型。

3.普通克里格插值算法

克里格插值算法包括簡單克里格插值算法、普通克里格插值算法、通用克里格插值算法、指標克里格插值算法、概率克里格插值算法、分離克里格插值算法、分層克里格插值算法、聯合克里格插值算法、因子克里格插值算法等20多種不同的變形形式（de Smith et al.，2007），但是所有的克里格插值算法都是基于式（2.24）的微小變異。

式中，m為整個區域內所有采樣數據的均值，λi是克里格權重，n是以x0為中心的、指定搜索區域內的參與克里格插值的采樣點個數，m（x0）是指定搜索區域內的采樣點均值。

當m為已知參數時，克里格插值稱為簡單克里格插值；當m為未知參數時，克里格插值稱為普通克里格插值。

從式（2.24）可以看出，克里格插值的關鍵在于求解克里格權重λi，并且克里格權重λi必須滿足無偏條件，使估計方差最小。

其中，無偏條件的數學表達式為

估計方差表示為

式中，為協方差函數，協方差函數和半變異函數之間具有如下關系：

式中，C（0）為區域化變量的Z（x）的方差。

要使估計方差在無偏條件下最小，則問題變為一個求解條件極值的方程，可以采用標準拉格朗日乘數法求解。依據拉格朗日原理構造函數F，即

式中，μ為拉格朗日乘數法。

分別對F求λi和μ的偏導，可得

式（2.29）是一個n+1階線性方程組，有n個未知數λi和1個未知數μ。因此，可以建立n+1維線性方程組，即

將已知采樣點數據代入式（2.30），即可解得λi。