Kriging的兩步極小法在GPS高程擬合中的應用*

王東東，郭 戩，趙旭坤，劉 兵，王 靜

(1.浙江建設職業技術學院，浙江 杭州 311231; 2.陜西鐵路工程職業技術學院，陜西 渭南 714000)

工程實踐過程中，利用GPS手段求取正常高，關鍵要先得到高程異常值。高程異常的獲取可以通過地球重力場模型求取，重力場方面的數據由于保密和工程成本的原因，很少使用；通過函數擬合的方法求取高程異常方法簡單，可操作性比較強，因此，這種方法有廣闊的應用前景。由于隨機場往往包含趨勢性和隨機性，因此采用兩步極小法進行處理更加合理[1]。Kriging方法是法國學者G.Materon于1962年在“應用地質統計學論”中首次提出和發展的。它是地質統計學的核心，是對局部區域隨機過程的一種預測方法,是建立在變異函數和結構分析的基礎之上，具有求線性、無偏和最優的內插估計值，統計特性是方差最小和平滑效應。普通Kriging法是以待測點為中心的鄰域內一些已測點數據、已測點之間的距離和變異函數提供的結構信息，得到對每個已測點的權系數，然后進行加權平均計算未測點的估值[2-5]。所謂兩步極小法，即首先利用二次多項式處理趨勢性部分[6-8]，然后Kriging處理隨機性部分，該方法兼顧了趨勢和隨機兩個部分，使構建模型更加符合實際情況。在變異函數的擬合過程中，對擬合點的選取很少考慮擬合點對擬合的貢獻大小，采用Cook距離篩選擬合點，找出一些對擬合貢獻較大的點，同時剔除一些對擬合產生干擾的點，可以提高變異函數的參數求解精度。

1 二次多項式擬合法

采用二次多項式函數作為曲面擬合的函數模型：

(1)

式中，ζi為高程異常值；a0、a1、a2、a3、a4、a5為待求參數；xi、yi為平面坐標；Δi為誤差。利用該函數模型擬合一個曲面，作為大地水準面，對待求點進行插值。該函數模型的矩陣形式為：

(2)

即

ζ=AX+Δ

(3)

誤差方程為：

v=ζ-AX

(4)

通過n(n≥6)個水準重合點，利用最小二乘法解算待求參數，帶入式(3)，解算出插值點的高程異常。

2 普通Kriging法原理

Kriging是基于變異函數的局部區域的最優無偏估計。假設區域變量滿足本征假設或二階平穩，其屬性值為Z(xi)。未測點的估計值是由該區域已測點的觀測值加權求和得到的，即為：

(5)

式中，Z*(x0)為未測點的估計值；λi為已測點屬性值的權系數。

由無偏性和方差最小可得：

(6)

由二階平穩和本征假設進一步可得：

(7)

式中，j=1,2,3,…,n;γ(xi-xj)為變異函數；μ為拉格朗日乘子。上式的矩陣形式為：

(8)

式中，γij為已測點xi和xj之間的變異函數；γi0為已測點xi和未測點x0之間的變異函數。由變異函數及式(8)解得λi和μ。變異函數的模型有指數模型、高斯模型和球狀模型[4]等，GPS高程擬合一般采用球面模型，其模型為：

(9)

式中，C0為塊金值；C為拱高；a為變程。

實際應用中，變異函數可由下式計算：

(10)

式中，N為已測點中距離為hm的點的對數。從模型上可以看出，它是一個距離的函數。

若要求取球面模型的參數，需要利用式(10)計算各個間隔上的變異函數值。計算出所有已測點之間的距離并按間距d進行分組。由于距離比較離散，點對之間的距離若在我們設定的步長范圍之內時，就歸為一組，此范圍內所有點對距離的平均值為本組變異函數的距離h，步長值可設置為d/2。分組模型為：

(11)

式中，maxh和minh為所有已測點之間距離的最大值和最小值；Nh為所有已測點之間距離的分組數。將求得的變異函數的值帶入式(9)，利用最小二乘法求解未知參數C0、C和a。由于球面模型是非線性函數，可將其線性化為：

γ(hm)=C0+k1h+k2h3

(12)

3 Cook距離的原理

1977年，Cook距離從參數置信區域提出，經過多年的發展已經應用在各種統計模型中。本文主要應用Cook距離篩選變異函數的擬合點，提高參數求解的精度，進而提高Kriging插值的精度。一般線性模型為：

(13)

(14)

(15)

Cook距離越小，表明剔除i行數據后，其對回歸參數的影響比較小，該點適合參與變異函數的擬合；Cook距離越大，表明剔除i行數據后，其對回歸參數的影響比較大，則該點不利于變異函數的參數的求解，擬合時應該將其剔除，從而達到提高變異函數參數求解的精度。

4 算例分析

為了驗證Cook距離對Kriging法中變異函數的參數求解中的有效性，選取了某地區地形平緩，無粗差、同精度的GPS點和水準重合點17個，其GPS網按照國家B級網施測，正常高按照二等水準測量施測。高程異常通過大地高和正常高容易得到，相關數據如表1所示。

表1 已知點坐標和高程異常[11]/m

本文設置三組實驗進行計算，三組實驗都采用兩步極小法進行處理。先進行二次多項式擬合，對數據的趨勢項部分進行處理；然后利用普通Kriging法對隨機部分進行處理。實驗過程中，采用球面模型作為變異函數。

方案Ⅰ：令點5、26、4、28為檢核點，其余13個點作為變異函數的擬合點。

方案Ⅱ：令點5、26、4、28、22、6、9為檢核點，其余10個點作為變異函數的擬合點。

方案Ⅲ：從17個點中選取4個Cook距離最大的值作為檢核點，如表2所示，分別為26、22、6、9，其余13個點作為變異函數的擬合點。

表 2 已知數據的Cook距離

從表2可以看出，6、9、22、26點的距離數值比較大，對回歸參數的影響較大，不利于模型的構建。

通過對三個實驗方案進行解算，得到二次多項式法和基于二次多項式的Kriging法的外符合精度，如表3所示。

表3 三個方案的外符合精度/mm

由實驗方案和表2可知，在擬合點的選擇上，方案Ⅰ是隨機的，方案Ⅱ和方案Ⅲ是經過Cook距離選擇的。三種方案的解算結果如表3所示，從計算結果比較中可以得出：

(1)方案Ⅲ中采用二次多項式法和基于二次多項式的Kriging法的外符合精度分別為4.2 mm、3.2 mm，均優于方案Ⅰ，表明經過Cook距離篩選擬合點是可行的，可以剔除部分強擾動點，進一步提高模型精度。

(2)方案Ⅰ、方案Ⅱ和方案Ⅲ中，基于二次多項式法的外符合精度分別為6.0 mm、3.1 mm、4.2 mm；基于二次多項式的Kriging法的外符合精度分別為4.9 mm、2.9 mm、3.2 mm。可以得出，基于二次多項式的Kriging法的外符合精度相比二次多項式法的外符合精度有一定的改善，驗證了兩步極小法能夠提高GPS高程擬合的精度。

(3)由實驗方案可知，方案Ⅱ采用的擬合點數量要比方案Ⅰ、方案Ⅲ少，方案Ⅱ采用二次多項式法和基于二次多項式的Kriging法的外符合精度均優于方案Ⅰ、方案Ⅲ，這是由于方案Ⅱ采用的擬合點是Cook距離較小的點，提高了模型解算精度。可以得出：在模型構建的過程中，不僅要顧及擬合點的數量，也要考慮擬合點擾動性的強度。

5 結 論

Kriging方法的關鍵之處是變異函數選擇和參數的求解，GPS高程擬合往往采用球面模型，此時只有得到高精度的參數解，才能得到很好的插值精度。變異函數參數求解過程中，由于一些點的強擾動性影響了參數求解的精度。利用Cook距離剔除了一些對求解過程中會產生較大干擾的點，可以提高變異函數參數求解的精度。Cook距離較大時表明該點離趨勢面較遠，屬于離群數據，選擇擬合點時要剔除這些離群數據，進一步提高建模的精度。兩步極小法可以兼顧隨機場的趨勢性和隨機性，實驗也驗證了這一理論的正確性。由于本文采用的數據點比較少，下一步將會對更大范圍的數據進行相關研究。