基于克里金法的地磁數據插值研究＊

于運治 田茂均

（海軍潛艇學院 青島 266044）

1 引言

對于水下載體，由于其不同于水上載體可以利用GPS等手段不間斷提供高質量的位置信息，它所能利用的定位手段極其有限。地磁場不像無線電波那樣在穿越水層時有衰減，它存在于大洋的任何深度。同時地磁場的場強和方向均是位置的函數，水下載體理論上可以通過測量地磁場要素與地磁圖匹配得到唯一的載體位置。地磁場特點使得其有望成為水下載體導航定位的又一重要手段，因此也成為了目前各國重點發展的載體導航手段之一［1］。

地磁數據服務水下載體導航主要的方式就是構成載體航行區域的先驗地磁圖，達到這一目的最主要可靠的方法就是通過實測地磁數據法構建。載體活動區域通常空間尺度都非常大，要想實測如此大區域獲取地磁數據無疑是一項巨大而且長遠的工程。由于艦船等本身對磁場的影響和測量經費的制約，實測法不可能實現地磁數據很密集的測量，因此獲得地磁圖的精度和密度都有限。鑒于這樣的實際情況，目前許多學者加大了對地磁空間數據插值算法的研究，力求通過研究出較好的插值算法彌補實測數據稀疏的問題實現地磁圖的加密重構，從而提高地磁圖的匹配概率和匹配精度。在眾多的空間數據插值算法中，克里金法又成為了其中的佼佼者。

2 克里金插值法原理

Kriging插值法是一種用于空間插值的數理統計方法，以法國科學家Krige的名字命名，是一種空間自協方差最優內插算法，實質是利用區域化變量的原始數據和變異函數的結構特點，對未知樣點進行線性無偏、最優估計［2］。它通過對已知樣本點賦權重來計算未知點的值，不僅考慮距離還通過變異函數和結構分析來考慮已知點與未知點的空間方位關系，這樣既反映了數據的空間結構變化又反映了變量的隨機分布，還能描述誤差信息，因此有著廣泛的應用［3］。

克里金法由于采用的模型不同和優化時又引入的趨勢、殘差、協變量、協方差等因子，克里金法又產生了許多變種，這里主要以普通克里金法進行分析，其原理如下［4～5］：

其中λi為已測量點Z（Xi）的權系數。上式中權系數的確定由克里金方程組決定，方程組可表示為

式中γ（Xi，Xj）為磁測量點的變異函數，μ為拉格朗日系數。變異函數的實質是一個協方差函數，是同一個變量在一定相隔距離上差值平方的期望。函數值越小，說明此段距離上該變量的相關性越好。因此變異函數可表示為

而當測點間的間距為h時有：

式中，h通常被稱為滯后矢量；Z（X）是位置X處的值；Z（X，h）為與位置X偏離量h處的實測值。同時，克里金法還給出了估值誤差的方差：

為求估計誤差方差的極小值，將拉格朗日系數引入有：

求解使估計估計誤差方差最小有：

由此得到：

通過對上面方程的求解就可得到插值點的估值［6］。另外，克里金法還提供了一個檢驗其可靠性的指標，估計值與實測值之間的均方差σ2，它通過計算權系數和拉格朗日系數來得到：

綜上看出克里金插值法不但考慮估計點和已知點之間的距離和空間分布，還考慮了已知點與已知點之間的分布，因為相較其它插值法而言有一定的優越性。

3 實例分析

實例采用的數據是某海區5km×5km區域的實測磁異常數據，原始數據量為101×101，間距50m。為了驗證插值方法的有效性，從中等間距的選取了11×11共計121個數據來作為“實驗”數據。也就相當于用目前間距為500m的數據經過插值以后得到間距為50m甚至更密集的新數據。

圖1是原始數據未經插值方法處理得到等值線圖。

3.1 三種常用空間數據插值算法效果對比

在空間數據插值方法中主要有三種，一是最早也是最簡單的加權反距離法，二是徑向基函數法，三是克里金法。其中基于加權反距離法的改進又有了改進的謝別得法。而克里金法主要有普通克里金法和泛克里金，當然如之前所述引入的趨勢、殘差、協變量等因子還有許多變種。這里首先應用加權反距離法、徑向基函數法和普通克里金法生成實驗數據的等值線圖。

圖1 原始數據等值線圖

圖2 實驗數據加權反距離插值圖

圖3 實驗數據徑向基函數插值圖

圖4 實驗數據普通克里金插值圖

如上圖我們可以看出加權反距離插值是不適用本例中地磁數據的，而徑向基函數法和普通克里金法在直觀上較為貼近實測數據。這里所采用的普通克里金法有一個前提就是默認數據是符合二階平穩的，而且沒有結算漂移。但是我們都知道，實測數據由于各種誤差是肯定存在一定漂移的，所以普通克里金通常都需要改進，泛克里金法就考慮了數據的非平穩性質和數據存在的漂移。

3.2 數據空間結構分析

泛克里金法通過變異函數分析空間數據的結構特點。變異函數在前一節克里金法原理已有提到，這里不贅述。這里主要從實驗變異函數角度實現磁異常數據這一區域化變量的變異函數建模。任何空間數據，在分析其結構時都要重點考慮其各向異性［11］。變異函數通過一些特點模型實現對數據各向異性的最大擬合，實驗變異函數模型有指數模型、對數模型、線性模型、高斯模型、球狀模型等等［7］。本例分析得出使用線性模型最能擬合數據特性，如圖5所示。

圖5 實驗數據變異函數線性模型

建模得出本例中數據的各項異性比為1.155，各項異性角度為142.6，傾率為0.1581。

3.3 交叉驗證

空間插值法研究過程中通常都采用交叉驗證法來評價插值算法的優劣程度。交叉驗證的基本思路是重復從測點集合中去掉一個或者多個測點，用剩余測點集合利用需要評價的插值算法估算得到被刪除點的值。然后用統計學的方法將測點上的估計值和實測值兩組數據進行統計分析，評價插值方法的精度［8］；常用的交叉驗證統計指標［9～10］有：

1）平均誤差（ME）

2）標準差（SD）

3）均方根誤差（RMSE）

表1 插值方法交叉驗證結果

4 結語

從上面的實例驗證，可以得出結論。首先從直觀的角度可以看出，如果值用實測的數據得出等值線，那么得到的等值線圖不但不平滑，不利于匹配導航使用，而且“牛眼”現象也比較嚴重。另外也可以排除加權反距離插值法，因為相比較于普通克里金法和徑向基函數法的插值效果，加權反距離法要遠遜于這兩者，而且比較原始數據可以看出，其所得到的等值線圖不符合原始數據的性質特征。其次，就交叉驗證的結果可以看出，在各項指標上，徑向基函數法是要優于普通克里金法的。但是相比較與考慮了數據空間結構的泛克里金，徑向基函數法又要稍遜于后者。

總結起來，通過插值的方法加密地磁數據，得到高密度的地磁數據，從而構建高精度的地磁圖方法上是非常有效的。從交叉驗證指標看，總體來講誤差都比較小，可以使用載體匹配導航需要。但是也可以明顯看出原始數據在中偏右部有一個楔形峰值，幾種插值方法都沒有擬合出。問題可能存在于一是測量值有存在的野值可能；其二文章采用等間距取點方式，勢必會再特殊峰值區存在插值信息不足的問題，這些有待下一步研究的解決。

［1］Goldenberg F.Geomagnetic navigation beyond the magnetic compass［C］／／Position Location and Navigation Symposium.Washington：IEEE，2006：684－694.

［2］侯景儒，黃競先.地質統計學的理論與方法［M］.北京：地質出版社，1990：263－280.

［3］Soo－ChangPei，Ching－Min Cheng.Color image processing by using binary quaternion－moment－preserving theres holding technique［J］.IEEE TANSACTIONS ON IMAGE PROCESSING，1999，8（5）：614－628.

［4］張楊，康崇，呂金庫，等.區域地磁測量實驗及水下載體對周圍磁場的影響分析［J］.中國慣性技術學報，2011，19（2）：205－208.

［5］靳國棟，劉衍聰，牛文杰.距離加權反比插值法和克里金插值法的比較［J］.長春工業大學學報，2003，24（3）：53－57.

［6］Rice H，Kelmenson S.Mendelsohn L.Geophysical navigation technologies and applications［C］／／Position Location and Navigation Symposium.April 26－29，2004：618－624.

［7］曾懷恩，黃聲響.基于Kriging方法的空間數據插值研究［J］.2007，10（5）：53－57.5－8.

［8］武俊紅，汪云甲.基于Surfer的煤礦等值線空間插值方法有效性評價［J］.中國礦業，2007，16（1）：108－110.

［9］湯剛安，劉學軍，閭國軍.數字高程模型及地學分析的原理與方法［M］.北京：科學出版社，2005，180－211.

［10］張曉明，趙剡.基于克里金插值的局部地磁圖的構建［J］.電子測量技術，2009，32（4）：122－125.

［11］張仁鐸.空間變異理論及應用［M］.北京：科學出版社，2005，273－284.