1、三種空間插值算法在地磁匹配基準(zhǔn)圖構(gòu)建中的應(yīng)用于運(yùn)治1田茂均1（1：海軍潛艇學(xué)院，青島，266044）摘要：針對地磁匹配導(dǎo)航對構(gòu)建高密度高精度地磁基準(zhǔn)圖的需求，結(jié)合目前應(yīng)用較為廣泛效果較好的三種空間插值方法，用某地區(qū)的實(shí)測地磁數(shù)據(jù)進(jìn)行了插值效果分析。對于插值效果的評價(jià)引入了交叉驗(yàn)證法，并借助誤差最大值、最小值、平均誤差、標(biāo)準(zhǔn)差和運(yùn)算時(shí)間等六項(xiàng)評價(jià)指標(biāo)對三種插值算法進(jìn)行了綜合評價(jià)。通過研究說明了三種插值算法各自的特點(diǎn)和應(yīng)用范圍。關(guān)鍵詞：地磁圖構(gòu)建；插值分析；交叉驗(yàn)證Theapplicationofthreespatialinterpolationalgorithmingeomagneticmatc

2、hingreferencemapbuildingYuYunzhiTianMaojuni(1:NavySubmrineAcademyQindao,266044)Abstract：Accordingtogeomagneticmatchingnavigationdemandforhigh-densityandhigh-precisiongeomagneticreferencemap,combinationthreebetterandwidelyusedspatialinterpolationmethod,withameasuredgeomagneticdatatoanalysisinterpolat

3、ioneffect.Introductionofcross-validationmethodtoevaluatetheinterpolationeffect,andwitherrormaximum,minimum,averageerror,standarddeviation,andcomputationtimesixfactorsevaluationofthreeinterpolationalgorithm.Bystudyingdescriptionofthefeaturesandapplicationsofeachofthethreeinterpolationalgorithmsingeom

4、agneticreferencemapbuliding.Keywords：geomagneticdata,interpolationanalysis,cross-validation1引言地磁場的場強(qiáng)和方向均是位置的函數(shù)，載體理論上可以通過測量地磁場要素，與地磁圖匹配得到唯一的載體位置。地磁場的特點(diǎn)使得地磁匹配導(dǎo)航具有了抗干擾、隱蔽、完全自主等優(yōu)點(diǎn)，也成為了目前各國重點(diǎn)發(fā)展的載體導(dǎo)航手段之一。地磁匹配導(dǎo)航與地形匹配等匹配導(dǎo)航手段的基本思路類似，它是通過載體實(shí)測的地磁數(shù)據(jù)信息與先驗(yàn)的地磁基準(zhǔn)圖進(jìn)行匹配來得到載體的位置信息。由此看出先驗(yàn)地磁圖的優(yōu)劣很大程度上決定了地磁匹配導(dǎo)航的位置精度。獲得先驗(yàn)地

5、磁基準(zhǔn)圖的方法主要有兩種，一是利用地磁模型法；二是通過實(shí)測地磁數(shù)據(jù)法。前一種方法建立的地磁圖尺度范圍大，往往是通過一個或多個數(shù)學(xué)物理模型建立較大區(qū)域乃至全球的地磁圖。這樣的地磁圖通常在實(shí)驗(yàn)室基于衛(wèi)星測磁數(shù)據(jù)，并通過大量的計(jì)算得到，但是由于建立的地磁圖空間尺度大，而且衛(wèi)星高度過高得到的地磁數(shù)據(jù)往往濾除了許多的磁異常信息，這樣的地磁圖并不適合于載體導(dǎo)航定位使用。后一種方法也是目前獲得地磁匹配基準(zhǔn)圖的主要方法，它使用艦船或者飛機(jī)搭載磁敏感器實(shí)測地磁數(shù)據(jù)來構(gòu)建地磁圖。但是由于艦船本身對磁場的影響和測量經(jīng)費(fèi)的制約，實(shí)測法不可能很密集的測量，因此獲得的地磁圖精度和密度都有限。目前許多學(xué)者加大了對地磁空間數(shù)

6、據(jù)插值算法的研究，力求通過較好的插值算法彌補(bǔ)實(shí)測數(shù)據(jù)稀疏的問題，實(shí)現(xiàn)地磁圖的插值加密重構(gòu)，從而提高地磁圖的匹配概率和精度。2三種主要的空間插值算法插值算法在許多工程領(lǐng)域都有極大的應(yīng)用，發(fā)展歷史也較長，按分類和性質(zhì)的不同也形成了為數(shù)眾多的方法。特別是20世紀(jì)初以來,將時(shí)間序列方法擴(kuò)展到二維與多維分析，對空間插值問題進(jìn)行了大量卓有成效的研究，建立了如反距離加權(quán)插值法、克里金法、徑向基函數(shù)法、樣條函數(shù)法、多元回歸法等一系列模型方法。其中反距離加權(quán)插值法、克里金法、徑向基函數(shù)法在地磁數(shù)據(jù)插值中的應(yīng)用較為普遍。2.1反距離加權(quán)插值法反距離加權(quán)插值法或是叫做加權(quán)反距離插值法是一種加權(quán)平均的插值方法。它認(rèn)為

7、未知點(diǎn)的值受較近點(diǎn)的影響要大于較遠(yuǎn)的點(diǎn)，并且這種影響與點(diǎn)之間距離的乘方成反比，也就是說離未知點(diǎn)越近的點(diǎn)所占的權(quán)重越大，越遠(yuǎn)的點(diǎn)所占的權(quán)重越小，如果乘方數(shù)為1說明變化率是恒定的，該方法就成為了線性插值法；如果乘方數(shù)為2或者更高階則說明較近點(diǎn)的影響遠(yuǎn)大于較遠(yuǎn)的點(diǎn)，而影響位置點(diǎn)值的所有有效點(diǎn)權(quán)重之和為1.0。其基本公式可表示為：dk其中Z0為未知點(diǎn)的值；Zi為影響Z0的第i個點(diǎn)值；k表示乘方數(shù)，可是1、2或者更大的數(shù)；d=；（x-x）2+（y一y）2，表示點(diǎn)Z）與Z。iii之間的距離；n為影響未知點(diǎn)的有效作用點(diǎn)數(shù)。反距離加權(quán)插值法的適用不僅取決于乘方數(shù)的選擇，還取決于確定有效作用點(diǎn)的方法。其次反距離

8、加權(quán)插值法是一種精確的插值法，插值得到的面最大值和最小值只會出現(xiàn)在采樣點(diǎn)，因?yàn)槠溆?jì)算易受到數(shù)據(jù)點(diǎn)集群的影響，因此計(jì)算結(jié)果有時(shí)會出現(xiàn)孤立點(diǎn)數(shù)據(jù)明顯高于周圍數(shù)據(jù)點(diǎn)的“鴨蛋”分布。而且，當(dāng)增加、刪除或者改變一個點(diǎn)時(shí)都需要從新計(jì)算權(quán)函數(shù),所以它存在著一定的缺陷。但是其優(yōu)點(diǎn)也非常明顯，因算法的結(jié)構(gòu)簡單，因此運(yùn)算非常迅速。2.2克里金（Kiging）法Kriging插值法是一種用于空間插值的數(shù)理統(tǒng)計(jì)方法，以法國科學(xué)家Krige的名字命名，是一種空間自協(xié)方差最優(yōu)內(nèi)插算法，實(shí)質(zhì)是利用區(qū)域化變量的原始數(shù)據(jù)和變異函數(shù)的結(jié)構(gòu)特點(diǎn)，對未知樣點(diǎn)進(jìn)行線性無偏、最優(yōu)估計(jì)。算法認(rèn)為數(shù)據(jù)的空間變化既不是完全隨機(jī)的也不是完全變化

9、的，而是包括了三個部分：一是具有恒定均值的確定性部分；二是與空間變化有關(guān)的隨機(jī)變量；三是與空間變化無關(guān)的隨機(jī)噪聲。因此運(yùn)用該方法首先要分析數(shù)據(jù)的空間結(jié)構(gòu)，提出相應(yīng)的變異函數(shù)，然后才能在該模型基礎(chǔ)上進(jìn)行插值運(yùn)算。它與反距離加權(quán)法的共同點(diǎn)在于都通過對已知樣本點(diǎn)賦權(quán)重來計(jì)算未知點(diǎn)的值;所不同的是在賦權(quán)重時(shí)，反距離加權(quán)法只考慮已知點(diǎn)與未知點(diǎn)的距離，而克里金插值法不僅考慮距離還通過變異函數(shù)和結(jié)構(gòu)分析來考慮已知點(diǎn)與未知點(diǎn)的空間方位關(guān)系，這樣既反映了數(shù)據(jù)的空間結(jié)構(gòu)變化又反映了變量的隨機(jī)分布，還能描述誤差信息，因此有著廣泛的應(yīng)用。克里金法常采用的變異函數(shù)有常數(shù)模型、指數(shù)模型、球狀模型、高斯模型等，由于采用的模

10、型不同和優(yōu)化時(shí)又引入的趨勢、殘差、協(xié)變量、協(xié)方差等因子，克里金法又產(chǎn)生了許多變種，這里主要以普通克里金法進(jìn)行分析，其原理如下4,5:對于研究區(qū)域A，設(shè)在區(qū)域內(nèi)采樣的位置坐標(biāo)為Xi，變量的觀測值為Z(Xj),i=1,2,3,.,n，則預(yù)測點(diǎn)X0的估計(jì)值7(X0)可用這n個樣本點(diǎn)的線性組合來表示：Z(x)=Xxz(x)0iii=1其中X.為已測量點(diǎn)Z(X)的權(quán)系數(shù)。上式中權(quán)系數(shù)的確定由克里金方程組決定，方程組可表示為：,X)+卩=y(X,X)ijj0/、(3)XX=1ii=1)為磁測量點(diǎn)的變異函數(shù)，卩為拉格朗日系數(shù)。變異函數(shù)的實(shí)質(zhì)是一個協(xié)方差函數(shù)，是同一個變量在一定相隔距離上差值平方的期望。函數(shù)值

11、越小，說明此段距離上該變量的相關(guān)性越好。因此變異函數(shù)可表示為：y(X,X)=y(XX)=1Et(X)Z(X)】(4)ijij2ij而當(dāng)是測點(diǎn)間的間距為h時(shí)有：y(h)=1Ez(X)Z(X,h)2(5)2式中，h通常被稱為滯后矢量；Z(x)是位置X處的值；Z(X,h)為與位置X偏離量h處的實(shí)測值。同時(shí)，克里金法還給出了估值誤差的方差:S=2乙Xy(XX)乙乙XXy(XX)(6)ii0ijiji=1i=1j=1為求估計(jì)誤差方差的極小值，將拉格朗日系數(shù)引入有：S=2XXy(XX)XXXXy(XX)2y(XX1)TOC o 1-5 h zii0ijijii=1i=1j=1i=1(7)求解使估計(jì)估計(jì)誤差

12、方差最小有：asas=o,=0a九a九12由此得到:y(x,x).11.y(x,x)11nX1y(x,x)01y(x,x).y(x,x)1X=y(x,x)n1nnn0n1.111通過對上面方程的求解就可得到插值點(diǎn)的估值。另外，克里金法還提供了一個檢驗(yàn)其可靠性的指標(biāo)，估計(jì)值與實(shí)測值之間的均方差a2，它通過計(jì)算權(quán)系數(shù)和拉格朗日系數(shù)來得到：a2=XXy(XX)+卩i0ii=1綜上看出克里金插值法不但考慮估計(jì)點(diǎn)和已知點(diǎn)之間的距離和空間分布，還考慮的已知點(diǎn)與已知點(diǎn)之間的分布，因?yàn)橄噍^其它插值法而言有一定的優(yōu)越性。2.3徑向基函數(shù)法徑向基函數(shù)法(RBF,RadialBasisFunction)是一種精確的

13、插值方法，它有一個很大的優(yōu)點(diǎn)是解決了系數(shù)矩陣奇異的問題。它主要適用于對大量點(diǎn)數(shù)據(jù)進(jìn)行插值計(jì)算，從而獲得平滑表面，且能預(yù)測比樣點(diǎn)高或低的未知點(diǎn)的值，因而具有較高的預(yù)測精度，能較好地反應(yīng)數(shù)據(jù)變化情況。其基本過程為：首先將插值函數(shù)構(gòu)造成下面形式的空間點(diǎn)函數(shù)：TOC o 1-5 h zF(x)=為九Q(|x一x|)(1)iii二1式中F(x)是已知點(diǎn)(i=1,2,.,n)；九是待定系ii數(shù)；0是指定的基函數(shù)；|是歐氏范數(shù)。當(dāng)F(x)滿足插值條件F(x)=f時(shí)，有下述i線性方程組：F(x)=九(|xx|)=“iii=1fi,i=1,2,.,n(11)令A(yù)=(|xx.|)；i尢=(九，1，九)Tnf=(f

14、,/.)T則有A九=f(12)由式(12)解出待定系數(shù)九可以定出F(x)。i基函數(shù)的選擇對徑向基函數(shù)法的精度至關(guān)重要，常用的有逆二次曲面、二次曲面、自然三次樣條、Gauss函數(shù)、薄板樣條函數(shù)。其中，復(fù)二次函數(shù)在大地測量、地質(zhì)及采礦、水文測量、地球物理等領(lǐng)域得到了較好的應(yīng)用效果，在數(shù)據(jù)點(diǎn)數(shù)量不太大的情況下計(jì)算量也不算大。因此本文采用的基函數(shù)是復(fù)二次函數(shù)(Multiquadratic)函數(shù)12:F(h)=(h2+R2)(13)其中h表示待插值點(diǎn)(x,y)到第i點(diǎn)(x,y)的距ii離；R參數(shù)是平滑因子。3數(shù)據(jù)驗(yàn)證數(shù)據(jù)來源采用的是某單位在某海區(qū)開展科研課題時(shí)所采集。數(shù)據(jù)是去掉了地磁主磁場的地磁異常場，

15、另外涉及到數(shù)據(jù)保密性等原因還去掉了經(jīng)緯度而以絕對距離表示。原始區(qū)域?yàn)?0kmX10km；數(shù)據(jù)量為201X201，精度50米。為了驗(yàn)證插值方法的有效性，從中等間距的選取了21X21共計(jì)441個數(shù)據(jù)來作為“采樣”數(shù)據(jù)。驗(yàn)證思路為使用“采樣”數(shù)據(jù)進(jìn)行插值處理得到與原始數(shù)據(jù)相同的數(shù)據(jù)量(201X201個)，然后通過圖形和交叉驗(yàn)證指標(biāo)數(shù)值兩個方面對插值算法進(jìn)行評價(jià)。31交叉驗(yàn)證空間插值法研究過程中常采用交叉驗(yàn)證法來評價(jià)插值算法的優(yōu)劣程度。交叉驗(yàn)證的基本思路是重復(fù)從測點(diǎn)集合中去掉一個測點(diǎn)，用剩余測點(diǎn)集合利用需要評價(jià)的插值算法估算得到被刪除點(diǎn)的值。然后用統(tǒng)計(jì)學(xué)的方法將測點(diǎn)上的估計(jì)值和實(shí)測值兩組數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分

16、析，評價(jià)插值方法的精度吻；常用的交叉驗(yàn)證統(tǒng)計(jì)指標(biāo)I0】有：誤差最大/小值(Max/Min)最大值和最小值一定程度上表征了算法的穩(wěn)定性，如果兩者相差較大，說明算法誤差波動范圍大。誤差均值(ME)ME:n；式中，g為真值與估計(jì)值ii之間的誤差亍也稱為殘差；它總體概括表征了預(yù)測數(shù)據(jù)和實(shí)測數(shù)據(jù)相近程度，但是由于存在正負(fù)抵消問題，表征能力有限。平均絕對離差(AAD)AAD=1XM它是殘差的絕對值均niei=1值；與均值不同的是它取了絕對值，因此能杜絕誤差中正負(fù)值抵消的問題，且與均值一樣對數(shù)據(jù)反應(yīng)敏銳、意義明確，能較好的反映數(shù)據(jù)的離散程度。只因?yàn)槿〗^對值不利于計(jì)算而已。標(biāo)準(zhǔn)差(SD)SD=.,Z(gME)

17、2.；n；式中，標(biāo)準(zhǔn)差越小,i=1i說明越多的估計(jì)值接近于實(shí)測值，同時(shí)說明插值算法的誤差比較穩(wěn)定。運(yùn)行時(shí)間(T)實(shí)驗(yàn)使用的數(shù)據(jù)量是較小的，但實(shí)際運(yùn)用過程中往往需要處理較大的數(shù)據(jù)量，因此運(yùn)行時(shí)間必須作為重要算法評價(jià)的重要指標(biāo)。3.2驗(yàn)證結(jié)果原始數(shù)據(jù)及實(shí)驗(yàn)數(shù)據(jù)經(jīng)插值后的等值線圖如下：4.2584.2564.2544.2524.253.523.543.563.583.6東西距離/m(a)原始數(shù)據(jù)等值線圖x1054.2483.5東西距離/mx105m儘距北南4.2584.2564.2544.2524.253.523.543.563.583.6東西距離/mx105(b)反距離加權(quán)法等值線圖(c)克里金法

18、等值線圖4.2584.2564.2544.2524.25(d)徑向基函數(shù)法等值線圖4243.53.523.543.563.58東西距離/m3.6105圖3.1原始數(shù)據(jù)與各插值算法結(jié)果等值線圖交叉驗(yàn)證各個指標(biāo)的計(jì)算結(jié)果如下：表3.1插值方法交叉驗(yàn)證結(jié)果統(tǒng)計(jì)比較空間插值方法MaxMinMEAADSDT反距離加權(quán)插值83.2243-137.18531.110925.059236.08922.42克里金插值281.4521-44.8306-1.74694.27657.974627.3徑向基函數(shù)插值25.9883-39.33671.12734.79137.885112.54結(jié)論從上述的等值線圖來看，加權(quán)

19、反距離法雖然大體上反映出了區(qū)域的等值線走向，但是“牛眼”眾多，與原始數(shù)據(jù)等值線圖相去甚遠(yuǎn)；而克里金法和徑向基函數(shù)法與原始數(shù)據(jù)結(jié)果都很接近，無法僅通過圖形得出答案。而從交叉驗(yàn)證結(jié)果看則反映出三種插值方法的這樣一些特點(diǎn)：一是三種方法在某一項(xiàng)指標(biāo)上都有優(yōu)于其它方法；二是基于不同的目的得到的最優(yōu)算法結(jié)論不同；三是反映出算法間可以彼此互補(bǔ)，說明改進(jìn)的空間很大。統(tǒng)計(jì)結(jié)果說明，總體來說插值算法效果基本滿足要求，但是進(jìn)一步分析可以對三種算法做出這樣的評價(jià)：加權(quán)反距離法運(yùn)算速度最快，非常適合于磁場的簡單反演插值，可以在極短時(shí)間內(nèi)得到區(qū)域磁場的大致特點(diǎn)；克里金法誤差最小，但是相對徑向基函數(shù)法來說極差較大，且運(yùn)算時(shí)

20、間也較長，更適合于離線計(jì)算；徑向基函數(shù)法誤差相對較小，能很好的反映磁場結(jié)構(gòu)，運(yùn)算時(shí)間也適中，可以適應(yīng)實(shí)時(shí)性要求不是特別高但又要保證一定插值精度的地磁圖構(gòu)建。通過該研究筆者認(rèn)為算法還可以從以下方面進(jìn)一步提咼效果：改變系數(shù)或者基函數(shù)等使之最大限度的符合數(shù)據(jù)的物理特性；加深數(shù)據(jù)插值前的分析，提高數(shù)據(jù)的完整性并多次驗(yàn)證找到最合適的模型參數(shù)。地磁匹配導(dǎo)航關(guān)鍵技術(shù)之一就是高密度高精度的地磁基準(zhǔn)圖構(gòu)建。隨著裝備制造和材料工藝的進(jìn)步以及算法研究的不斷深入，地磁匹配導(dǎo)航將會取得更大的發(fā)展。參考文獻(xiàn)：GoldenbergF.GeomagneticnavigationbeyondthemagneticcompassC/PositionLocationandNavigationSymposium.Washington:IEEE,2006:684-694.朱會義，劉述林，賈紹鳳自然地理要素空間插值的幾個問題J.地理研究，2004,23(4)：425-432.黃學(xué)功，房建成，劉剛，蔣顏瑋地磁圖制備方法及其有效性評估J.北京航空航天大學(xué)學(xué)報(bào)，2009,35(7)：891-894.張楊，康崇，呂金庫，楊月芳，張悅區(qū)域地磁測量實(shí)驗(yàn)及水下載體對周圍磁場的影響分析J.中國慣性技術(shù)學(xué)報(bào)，2011,19(2)：205-208.靳國棟，劉衍聰,牛文杰距離加權(quán)反比插值法和克里金插值法的比較J.長春工