Kriging的兩步極小法在GPS高程擬合中的應(yīng)用*

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

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

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

1 二次多項式擬合法

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

(1)

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

(2)

即

ζ=AX+Δ

(3)

誤差方程為：

v=ζ-AX

(4)

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

2 普通Kriging法原理

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

(5)

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

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

(6)

由二階平穩(wěn)和本征假設(shè)進一步可得：

(7)

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

(8)

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

(9)

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

實際應(yīng)用中，變異函數(shù)可由下式計算：

(10)

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

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

(11)

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

γ(hm)=C0+k1h+k2h3

(12)

3 Cook距離的原理

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

(13)

(14)

(15)

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

4 算例分析

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

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

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

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

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

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

表 2 已知數(shù)據(jù)的Cook距離

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

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

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

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

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

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

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

5 結(jié) 論

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