基于克里金法的地磁數(shù)據(jù)插值研究＊

于運(yùn)治 田茂均

（海軍潛艇學(xué)院 青島 266044）

1 引言

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

地磁數(shù)據(jù)服務(wù)水下載體導(dǎo)航主要的方式就是構(gòu)成載體航行區(qū)域的先驗(yàn)地磁圖，達(dá)到這一目的最主要可靠的方法就是通過(guò)實(shí)測(cè)地磁數(shù)據(jù)法構(gòu)建。載體活動(dòng)區(qū)域通?？臻g尺度都非常大，要想實(shí)測(cè)如此大區(qū)域獲取地磁數(shù)據(jù)無(wú)疑是一項(xiàng)巨大而且長(zhǎng)遠(yuǎn)的工程。由于艦船等本身對(duì)磁場(chǎng)的影響和測(cè)量經(jīng)費(fèi)的制約，實(shí)測(cè)法不可能實(shí)現(xiàn)地磁數(shù)據(jù)很密集的測(cè)量，因此獲得地磁圖的精度和密度都有限。鑒于這樣的實(shí)際情況，目前許多學(xué)者加大了對(duì)地磁空間數(shù)據(jù)插值算法的研究，力求通過(guò)研究出較好的插值算法彌補(bǔ)實(shí)測(cè)數(shù)據(jù)稀疏的問(wèn)題實(shí)現(xiàn)地磁圖的加密重構(gòu)，從而提高地磁圖的匹配概率和匹配精度。在眾多的空間數(shù)據(jù)插值算法中，克里金法又成為了其中的佼佼者。

2 克里金插值法原理

Kriging插值法是一種用于空間插值的數(shù)理統(tǒng)計(jì)方法，以法國(guó)科學(xué)家Krige的名字命名，是一種空間自協(xié)方差最優(yōu)內(nèi)插算法，實(shí)質(zhì)是利用區(qū)域化變量的原始數(shù)據(jù)和變異函數(shù)的結(jié)構(gòu)特點(diǎn)，對(duì)未知樣點(diǎn)進(jìn)行線性無(wú)偏、最優(yōu)估計(jì)［2］。它通過(guò)對(duì)已知樣本點(diǎn)賦權(quán)重來(lái)計(jì)算未知點(diǎn)的值，不僅考慮距離還通過(guò)變異函數(shù)和結(jié)構(gòu)分析來(lái)考慮已知點(diǎn)與未知點(diǎn)的空間方位關(guān)系，這樣既反映了數(shù)據(jù)的空間結(jié)構(gòu)變化又反映了變量的隨機(jī)分布，還能描述誤差信息，因此有著廣泛的應(yīng)用［3］。

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

其中λi為已測(cè)量點(diǎn)Z（Xi）的權(quán)系數(shù)。上式中權(quán)系數(shù)的確定由克里金方程組決定，方程組可表示為

式中γ（Xi，Xj）為磁測(cè)量點(diǎn)的變異函數(shù)，μ為拉格朗日系數(shù)。變異函數(shù)的實(shí)質(zhì)是一個(gè)協(xié)方差函數(shù)，是同一個(gè)變量在一定相隔距離上差值平方的期望。函數(shù)值越小，說(shuō)明此段距離上該變量的相關(guān)性越好。因此變異函數(shù)可表示為

而當(dāng)測(cè)點(diǎn)間的間距為h時(shí)有：

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

為求估計(jì)誤差方差的極小值，將拉格朗日系數(shù)引入有：

求解使估計(jì)估計(jì)誤差方差最小有：

由此得到：

通過(guò)對(duì)上面方程的求解就可得到插值點(diǎn)的估值［6］。另外，克里金法還提供了一個(gè)檢驗(yàn)其可靠性的指標(biāo)，估計(jì)值與實(shí)測(cè)值之間的均方差σ2，它通過(guò)計(jì)算權(quán)系數(shù)和拉格朗日系數(shù)來(lái)得到：

綜上看出克里金插值法不但考慮估計(jì)點(diǎn)和已知點(diǎn)之間的距離和空間分布，還考慮了已知點(diǎn)與已知點(diǎn)之間的分布，因?yàn)橄噍^其它插值法而言有一定的優(yōu)越性。

3 實(shí)例分析

實(shí)例采用的數(shù)據(jù)是某海區(qū)5km×5km區(qū)域的實(shí)測(cè)磁異常數(shù)據(jù)，原始數(shù)據(jù)量為101×101，間距50m。為了驗(yàn)證插值方法的有效性，從中等間距的選取了11×11共計(jì)121個(gè)數(shù)據(jù)來(lái)作為“實(shí)驗(yàn)”數(shù)據(jù)。也就相當(dāng)于用目前間距為500m的數(shù)據(jù)經(jīng)過(guò)插值以后得到間距為50m甚至更密集的新數(shù)據(jù)。

圖1是原始數(shù)據(jù)未經(jīng)插值方法處理得到等值線圖。

3.1 三種常用空間數(shù)據(jù)插值算法效果對(duì)比

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

圖1 原始數(shù)據(jù)等值線圖

圖2 實(shí)驗(yàn)數(shù)據(jù)加權(quán)反距離插值圖

圖3 實(shí)驗(yàn)數(shù)據(jù)徑向基函數(shù)插值圖

圖4 實(shí)驗(yàn)數(shù)據(jù)普通克里金插值圖

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

3.2 數(shù)據(jù)空間結(jié)構(gòu)分析

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

圖5 實(shí)驗(yàn)數(shù)據(jù)變異函數(shù)線性模型

建模得出本例中數(shù)據(jù)的各項(xiàng)異性比為1.155，各項(xiàng)異性角度為142.6，傾率為0.1581。

3.3 交叉驗(yàn)證

空間插值法研究過(guò)程中通常都采用交叉驗(yàn)證法來(lái)評(píng)價(jià)插值算法的優(yōu)劣程度。交叉驗(yàn)證的基本思路是重復(fù)從測(cè)點(diǎn)集合中去掉一個(gè)或者多個(gè)測(cè)點(diǎn)，用剩余測(cè)點(diǎn)集合利用需要評(píng)價(jià)的插值算法估算得到被刪除點(diǎn)的值。然后用統(tǒng)計(jì)學(xué)的方法將測(cè)點(diǎn)上的估計(jì)值和實(shí)測(cè)值兩組數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析，評(píng)價(jià)插值方法的精度［8］；常用的交叉驗(yàn)證統(tǒng)計(jì)指標(biāo)［9～10］有：

1）平均誤差（ME）

2）標(biāo)準(zhǔn)差（SD）

3）均方根誤差（RMSE）

表1 插值方法交叉驗(yàn)證結(jié)果

4 結(jié)語(yǔ)

從上面的實(shí)例驗(yàn)證，可以得出結(jié)論。首先從直觀的角度可以看出，如果值用實(shí)測(cè)的數(shù)據(jù)得出等值線，那么得到的等值線圖不但不平滑，不利于匹配導(dǎo)航使用，而且“牛眼”現(xiàn)象也比較嚴(yán)重。另外也可以排除加權(quán)反距離插值法，因?yàn)橄啾容^于普通克里金法和徑向基函數(shù)法的插值效果，加權(quán)反距離法要遠(yuǎn)遜于這兩者，而且比較原始數(shù)據(jù)可以看出，其所得到的等值線圖不符合原始數(shù)據(jù)的性質(zhì)特征。其次，就交叉驗(yàn)證的結(jié)果可以看出，在各項(xiàng)指標(biāo)上，徑向基函數(shù)法是要優(yōu)于普通克里金法的。但是相比較與考慮了數(shù)據(jù)空間結(jié)構(gòu)的泛克里金，徑向基函數(shù)法又要稍遜于后者。

總結(jié)起來(lái)，通過(guò)插值的方法加密地磁數(shù)據(jù)，得到高密度的地磁數(shù)據(jù)，從而構(gòu)建高精度的地磁圖方法上是非常有效的。從交叉驗(yàn)證指標(biāo)看，總體來(lái)講誤差都比較小，可以使用載體匹配導(dǎo)航需要。但是也可以明顯看出原始數(shù)據(jù)在中偏右部有一個(gè)楔形峰值，幾種插值方法都沒(méi)有擬合出。問(wèn)題可能存在于一是測(cè)量值有存在的野值可能；其二文章采用等間距取點(diǎn)方式，勢(shì)必會(huì)再特殊峰值區(qū)存在插值信息不足的問(wèn)題，這些有待下一步研究的解決。

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

［2］侯景儒，黃競(jìng)先.地質(zhì)統(tǒng)計(jì)學(xué)的理論與方法［M］.北京：地質(zhì)出版社，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］張楊，康崇，呂金庫(kù)，等.區(qū)域地磁測(cè)量實(shí)驗(yàn)及水下載體對(duì)周圍磁場(chǎng)的影響分析［J］.中國(guó)慣性技術(shù)學(xué)報(bào)，2011，19（2）：205－208.

［5］靳國(guó)棟，劉衍聰，牛文杰.距離加權(quán)反比插值法和克里金插值法的比較［J］.長(zhǎng)春工業(yè)大學(xué)學(xué)報(bào)，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方法的空間數(shù)據(jù)插值研究［J］.2007，10（5）：53－57.5－8.

［8］武俊紅，汪云甲.基于Surfer的煤礦等值線空間插值方法有效性評(píng)價(jià)［J］.中國(guó)礦業(yè)，2007，16（1）：108－110.

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

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

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