趙玉新, 邢文

哈爾濱工程大學(xué) 自動(dòng)化學(xué)院，黑龍江 哈爾濱 150001

?

基于多重分形克里金的逐步插值校正法構(gòu)建局部地磁基準(zhǔn)圖

趙玉新, 邢文

哈爾濱工程大學(xué) 自動(dòng)化學(xué)院，黑龍江 哈爾濱 150001

摘要：為滿足地磁導(dǎo)航技術(shù)對(duì)高精度地磁基準(zhǔn)圖的需求，對(duì)局部地磁基準(zhǔn)圖構(gòu)建技術(shù)進(jìn)行深入研究。針對(duì)克里金插值法具有低通濾波性的缺陷及多重分形克里金法對(duì)稀疏實(shí)測的地磁異常場數(shù)據(jù)進(jìn)行插值時(shí)分形特征丟失的問題，提出一種基于多重分形克里金的逐步插值校正法，該方法從增加采樣點(diǎn)密度的角度出發(fā)，考慮了邊界效應(yīng)問題，利用克里金法將稀疏數(shù)據(jù)插值為過渡分辨率，再用多重分形克里金法進(jìn)行奇異性校正。仿真實(shí)驗(yàn)表明，該方法的構(gòu)圖精度具有明顯優(yōu)勢，適用于構(gòu)建高精度、高分辨率的局部地磁基準(zhǔn)圖。

關(guān)鍵詞：地磁導(dǎo)航；克里金法；逐步插值校正；分形理論；多重分形克里金插值法；地磁基準(zhǔn)圖

邢文(1990-)，女，碩士.

地磁導(dǎo)航技術(shù)具有自主性、全天時(shí)、全天候、高精度等優(yōu)點(diǎn)，是現(xiàn)有導(dǎo)航方式的重要補(bǔ)充，而高精度的局部地磁基準(zhǔn)圖是實(shí)現(xiàn)精確地磁導(dǎo)航的重要基礎(chǔ)。地磁基準(zhǔn)圖的制備以實(shí)測地磁數(shù)據(jù)為基礎(chǔ)，并可分為2種方式，一是局部地磁場建模方法[1-4]，二是基于插值估計(jì)的數(shù)據(jù)加密方法。前者在描述小尺度的磁場特征，尤其是磁異常信息時(shí)，需要很高的模型階數(shù)，計(jì)算量大，而且很難充分利用高密度的局部磁場測量數(shù)據(jù)，造成浪費(fèi)；后者主要是根據(jù)未知位置周圍的已知數(shù)據(jù)對(duì)其進(jìn)行估計(jì)，能夠充分利用所有測量數(shù)據(jù)資源，而且計(jì)算效率高，因而在精確地磁導(dǎo)航研究中，得到了更為廣泛的應(yīng)用。

黃學(xué)功[5]、王哲[6]、張曉明[7]、楊功流[8]等對(duì)制備地磁圖的幾種網(wǎng)格插值方法進(jìn)行評(píng)估，實(shí)驗(yàn)結(jié)果顯示徑向基函數(shù)、克里金法具有較高的插值精度。但是對(duì)奇異性較強(qiáng)的局部區(qū)域的刻畫誤差較大。Spector[9]等在對(duì)航磁數(shù)據(jù)譜進(jìn)行統(tǒng)計(jì)分析的過程中證實(shí)了地磁異常具有分形特征。大多研究[10-12]還局限于對(duì)整體分形特征的描述，在小尺度上，例如對(duì)文中涉及的局部地磁異常場分形特征的研究較少。文獻(xiàn)[13]利用逐步插值校正法刻畫了磁異常場在小尺度上的奇異特征，但該方法在實(shí)測數(shù)據(jù)較豐富的情況下插值效果較好。

針對(duì)上述問題，文中提出一種基于多重分形克里金的逐步插值校正法，實(shí)驗(yàn)表明，該方法構(gòu)建的局部地磁基準(zhǔn)圖精度較高，為后續(xù)將其應(yīng)用于地磁匹配導(dǎo)航作保障。

1克里金法及多重分形理論

1.1克里金法基本原理

克里金法是插值精度較高的空間插值方法之一。它是在平均殘差接近于零、估計(jì)誤差方差最小的前提下，將已有數(shù)據(jù)加權(quán)線性結(jié)合得到估計(jì)值，是一種最優(yōu)、線性且無偏的估計(jì)。它通過以距離為自變量的變異函數(shù)計(jì)算權(quán)值，相對(duì)于其他空間插值方法而言，更加符合實(shí)際?？死锝鸩逯捣ㄔ诠烙?jì)待插值點(diǎn)屬性值的過程中，不僅考慮了待插值點(diǎn)與鄰近已知點(diǎn)的空間位置關(guān)系，還考慮了各相鄰已知點(diǎn)之間的位置關(guān)系，充分利用了已知數(shù)據(jù)點(diǎn)的空間分布結(jié)構(gòu)特征。但由于克里金插值的推導(dǎo)過程是建立在空間相關(guān)系數(shù)矩陣極小意義下的，這決定了它是一種滑動(dòng)加權(quán)平均或低通濾波過程，因此不可避免地削弱了高頻、局部及弱信號(hào)，無法重現(xiàn)局部區(qū)域隆起或凹陷的奇異特征。

普通克里金插值法的估值公式為

式中：Z*為待插值點(diǎn)x0處的屬性值；Zi為估計(jì)區(qū)域內(nèi)1～n個(gè)采樣點(diǎn)的實(shí)測地磁異常場強(qiáng)度值；λi為待求權(quán)重系數(shù)。

1.2 多重分形理論

多重分形研究的關(guān)鍵問題之一就是測度與尺度的關(guān)系。對(duì)于曲面而言，δ尺度下的測度是以s點(diǎn)為中心、δ為邊長的正方形鄰域的質(zhì)量。表示為

(1)

根據(jù)分維定義，擬合最佳逼近時(shí)，尺度與測度的關(guān)系可表示為

(2)

式中：a、b為常數(shù)。a為奇異系數(shù)，表征了待插值點(diǎn)小鄰域內(nèi)的凹凸特性。在測度—尺度的雙對(duì)數(shù)坐標(biāo)中，a是在最小方差意義下所擬合直線的斜率。

2多重分形克里金插值法

2.1多重分形克里金插值法基本原理

多重分形克里金插值法的原理為利用多重分形的理論對(duì)待插值點(diǎn)小鄰域內(nèi)克里金插值結(jié)果進(jìn)行奇異性校正，即用較大鄰域內(nèi)地磁異常場強(qiáng)度的均值去估計(jì)中心小鄰域內(nèi)的均值，作為待插值點(diǎn)的插值結(jié)果。

據(jù)式(1)、(2)，將尺度為N·δ鄰域及尺度為δ鄰域內(nèi)的測度表示為

(3)

(4)

式(3)與式(4)聯(lián)立得

式中：Z*為N·δ鄰域內(nèi)地磁異常場強(qiáng)度的均值，也是克里金法對(duì)基準(zhǔn)數(shù)據(jù)的估值結(jié)果；Zδ為δ鄰域內(nèi)地磁異常場強(qiáng)度的均值，也是待插值點(diǎn)的地磁異常場估計(jì)值；N為最大尺度與最小尺度的比值；N2-a為校正系數(shù)。

圖1　二維空間測度計(jì)算示意

2.2基于多重分形克里金插值法的仿真實(shí)驗(yàn)

本文選取3個(gè)實(shí)驗(yàn)區(qū)域?qū)Σ逯邓惴ㄟM(jìn)行仿真實(shí)驗(yàn)，并選取4個(gè)有效性評(píng)價(jià)指標(biāo)對(duì)插值后的網(wǎng)格數(shù)據(jù)進(jìn)行評(píng)估，如表1所示。由于在研究區(qū)域的邊界缺少參考點(diǎn)，以至于邊界處插值精度較低，因此文中在構(gòu)建地磁基準(zhǔn)圖的仿真過程中均是在考慮邊界效應(yīng)的情況下進(jìn)行的。經(jīng)實(shí)驗(yàn)驗(yàn)證，第1周和第2周的參考點(diǎn)對(duì)邊界區(qū)域的構(gòu)圖精度有著重要的影響，不可缺少。為了方便，將待插值數(shù)據(jù)的最外圍5周作為邊界區(qū)域。

表1　 有效性評(píng)價(jià)指標(biāo)

研究區(qū)域的編號(hào)分別為4 145、4 053、4 109B，它們的地磁異常場強(qiáng)度數(shù)據(jù)來源于美國國家海洋和大氣管理局(national oceanic and atmospheric administration，NOAA)發(fā)布的航空磁測數(shù)據(jù)。表2列出了原始數(shù)據(jù)、基準(zhǔn)數(shù)據(jù)及中心區(qū)域數(shù)據(jù)的經(jīng)緯度范圍。基準(zhǔn)數(shù)據(jù)由原始數(shù)據(jù)網(wǎng)格化得到，共257×257個(gè)點(diǎn)，經(jīng)、緯度分辨率均為0.001°。從基準(zhǔn)數(shù)據(jù)中每隔3個(gè)點(diǎn)提取1個(gè)點(diǎn)作為待插值點(diǎn)，則一共可提取出65×65個(gè)點(diǎn)組成待插值區(qū)域(即中心區(qū)域)，其經(jīng)、緯度分辨率均為0.004°。

表2　 區(qū)域4145、4053、4109B的經(jīng)緯度范圍

對(duì)編號(hào)為4 053、4 145、4 109B數(shù)據(jù)的中心區(qū)域進(jìn)行多重分形克里金插值，流程如圖2所示，比較多重分形克里金插值法與克里金插值法的構(gòu)圖精度。

圖2　多重分形克里金插值法示意

改變盒子數(shù)N的值，觀察其對(duì)插值精度的影響。以4 053數(shù)據(jù)為例，實(shí)驗(yàn)結(jié)果見表3。

表3　 4 053不同盒子數(shù)下的多重分形克里金插值精度驗(yàn)證

由表3的數(shù)據(jù)可以看出，多重分形克里金插值法的構(gòu)圖精度不高，原因可能有2點(diǎn)：

1)現(xiàn)有數(shù)據(jù)較稀疏，在用局部去估計(jì)整體的過程中，不可避免地導(dǎo)致部分分形特征丟失。

2)在某些點(diǎn)上克里金插值的結(jié)果已經(jīng)與基準(zhǔn)數(shù)據(jù)很接近，再用校正系數(shù)去校正，就會(huì)增大結(jié)果的誤差。

為解決上述不足，提出基于多重分形克里金的逐步插值校正法。

3基于多重分形克里金法的逐步插值校正構(gòu)圖方法

為解決分形特征丟失的問題，現(xiàn)從增加采樣點(diǎn)密度的思路出發(fā)，采用逐步插值校正的方法完成多重分形克里金插值過程，提出基于多重分形克里金的逐步插值校正法。其原理為先對(duì)分辨率為0.004°×0.004°，65×65個(gè)已知待插值數(shù)據(jù)點(diǎn)進(jìn)行克里金插值，插值成分辨率為0.002°×0.002°的數(shù)據(jù)，共129×129個(gè)點(diǎn)，完成加密重構(gòu)過程。再利用插值后的數(shù)據(jù)計(jì)算校正系數(shù)，完成多重分形克里金插值，如圖3所示，提取中心區(qū)域數(shù)據(jù)進(jìn)行精度驗(yàn)證。

圖3　基于多重分形克里金的逐步插值校正法示意

觀察表3可知，對(duì)于4 053這組數(shù)據(jù)，多重分形克里金法在N=2時(shí)插值精度最高，因此仿真過程中盒子數(shù)N均取為2。同理4 145、4 109B這兩組數(shù)據(jù)的盒子數(shù)分別取為2和3。盒子數(shù)是通過反復(fù)實(shí)驗(yàn)測得的，文中假設(shè)已知數(shù)據(jù)兩點(diǎn)之間的距離為單位1，得到盒子數(shù)的最佳取值。仿真結(jié)果見表4和圖4。

表4　 文中方法的構(gòu)圖精度

(a) 基準(zhǔn)數(shù)據(jù)中心區(qū)域三維圖

(b) 克里金插值構(gòu)建的地磁基準(zhǔn)圖

(c) 文中方法所構(gòu)地磁基準(zhǔn)圖圖4　區(qū)域4053三維比較圖

分析表4的數(shù)據(jù)可知，文中方法的構(gòu)圖精度優(yōu)于克里金插值法的構(gòu)圖精度。由圖4也可看出，基于逐步插值校正的多重分形克里金插值法對(duì)局部位置的刻畫較細(xì)致，能夠表現(xiàn)出更多的隆起或下陷等奇異特征，與基準(zhǔn)數(shù)據(jù)的三維圖形態(tài)更加相近，具有良好的插值效果。

為了更清楚地觀察文中方法所構(gòu)地磁基準(zhǔn)圖的準(zhǔn)確性，將觀測角度細(xì)化到每一條剖線上。圖5為4 053數(shù)據(jù)中x=82的剖線，圖6為4145數(shù)據(jù)中x=166的剖線，圖7為4 019B數(shù)據(jù)中x=106的剖線。

圖5　區(qū)域4 053中的剖線x=82

圖6　區(qū)域4 145中的剖線x=166

圖7　區(qū)域4 109B中的剖線x=106

從剖線圖中更能清晰地看出，文中方法所構(gòu)地磁基準(zhǔn)圖在細(xì)節(jié)上與基準(zhǔn)數(shù)據(jù)曲線更加相似。文中方法是在克里金插值法描繪整體趨勢的基礎(chǔ)上，采用多重分形插值對(duì)局部位置逐步進(jìn)行奇異性校正。

實(shí)驗(yàn)表明，在數(shù)據(jù)非常稀疏的情況下，文中方法相對(duì)于常用插值算法而言仍具有較高的插值精度，如表5所示。表中數(shù)據(jù)的分辨率為0.016°，中心區(qū)域共7×7個(gè)采樣點(diǎn)。

表5　 編號(hào)4 053稀疏數(shù)據(jù)的插值精度比較

綜上所述，基于多重分形克里金的逐步插值校正法的構(gòu)圖精度較高，局部位置奇異特征較明顯，在局部地磁基準(zhǔn)圖構(gòu)建方面具有顯著的優(yōu)勢。

5結(jié)束語

文中為改善克里金插值法的低通濾波性，引入多重分形理論，增加對(duì)局部位置的奇異性描述，還原輸入信號(hào)中更多的高頻信息，研究了多重分形克里金插值方法。由于稀疏的已知采樣點(diǎn)直接插值會(huì)導(dǎo)致部分分形特征丟失，插值精度降低，因此文中從增加參考數(shù)據(jù)密度的角度出發(fā)，在考慮邊界效應(yīng)問題的條件下，提出基于多重分形克里金的逐步插值校正法。實(shí)驗(yàn)證明該方法相對(duì)于克里金插值法得到了較大改善，適合于構(gòu)建高精度的局部地磁基準(zhǔn)圖。

參考文獻(xiàn):

[1]喬玉坤, 王仕成, 張金生, 等. 采用矩諧分析和支持向量機(jī)的地磁導(dǎo)航基準(zhǔn)圖構(gòu)建方法[J]. 西安交通大學(xué)學(xué)報(bào), 2010, 44(10): 47-51.

[2]李明明, 黃顯林, 盧鴻謙, 等. 基于矩諧分析的高精度局部地磁場建模研究[J]. 宇航學(xué)報(bào), 2010, 31(7): 1730-1736.

[3]喬玉坤, 王仕成, 張金生, 等. 基于BP網(wǎng)絡(luò)的地磁基準(zhǔn)圖制備及其精度評(píng)價(jià)[J]. 中國慣性技術(shù)學(xué)報(bào), 2009, 17(1): 53-58.

[4]孫淵, 張金生, 王仕成, 等. 基于主磁場長期變化模型的地磁導(dǎo)航基準(zhǔn)圖時(shí)變修正[J].中國慣性技術(shù)學(xué)報(bào), 2011, 19(5): 543-548, 552.

[5]黃學(xué)功, 房建成, 劉剛, 等. 地磁圖制備方法及其有效性評(píng)估[J]. 北京航空航天大學(xué)學(xué)報(bào), 2009, 35(7): 891-894.

[6]王哲, 王仕成, 張金生, 等. 一種地磁匹配制導(dǎo)基準(zhǔn)圖制備方法及其有效性評(píng)價(jià)[J]. 系統(tǒng)工程與電子技術(shù), 2008, 30(11): 2207-2211.

[7]張曉明, 趙剡. 基于克里金法插值的局部地磁圖的構(gòu)建[J]. 電子測量技術(shù), 2009, 32(4): 122-125.

[8]楊功流, 張桂敏, 李士心. 泛克里金插值法在地磁圖中的應(yīng)用[J]. 中國慣性技術(shù)學(xué)報(bào), 2008, 16(2): 162-166.

[9]SPECTOR A, GRANT F S. Statistical models for interpreting aeromagnetic data[J]. Geophysics, 1970, 35(2): 293-302.

[10]孫霞, 吳自勤. 規(guī)則表面形貌的分形和多重分形描述[J]. 物理學(xué)報(bào), 2011, 50(11): 2126-2131.

[11]楊娟, 卞保民, 閆振綱, 等. 典型隨機(jī)信號(hào)特征參數(shù)統(tǒng)計(jì)分布的分形特性[J]. 物理學(xué)報(bào), 2011, 60(10): 86-92.

[12]馮馳, 胡楊, 王兆豐. 基于分形理論的渦輪葉片特征提取[J]. 應(yīng)用科技, 2015, 42(4): 64-69.

[13]趙玉新, 常帥, 張振興. 地磁異常場的多重分形譜分析及構(gòu)圖法[J]. 測繪學(xué)報(bào), 2014, 43(5): 529-536. doi: 10.11991/yykj.201503032

網(wǎng)絡(luò)出版地址：http://www.cnki.net/kcms/detail/23.1191.U.20151206.1009.008.html

Building a local geomagnetic reference map by the multifractal Kriging

based step-by-step interpolation correction method

ZHAO Yuxin, XING Wen

College of Automation, Harbin Engineering University, Harbin 150001, China

Abstract：To meet the demand of geomagnetic navigation technology for high precision geomagnetic reference maps, the method of building a local geomagnetic reference map is studied in this paper. In order to improve the low-pass filtering characteristics of the Kriging method and reduce the loss of fractal characteristics when interpolating with sparse geomagnetic anomaly field data by multifractal Kriging, this paper proposes a step-by-step interpolation correction method based on the multifractal Kriging. This method increases the sampling point density and considers the boundary effect, firstly interpolating sparse data into a transitional resolution by Kriging and then conducting singularity correction by multifractal Kriging. Experimental result shows that this method has obvious advantages for building a high-precision high-resolution local geomagnetic reference map.

Keywords：geomagnetic navigation; Kriging; step-by-step interpolation correction; fractal theory; multifractal Kriging interpolation; geomagnetic reference map

通信作者：邢文，E-mail：544166264@qq.com.

作者簡介：趙玉新(1980-)，男，教授，博士生導(dǎo)師；

基金項(xiàng)目：國家自然科學(xué)基金資助項(xiàng)目(51379049).

收稿日期：2015-03-25.網(wǎng)絡(luò)出版日期：2015-12-06.

中圖分類號(hào)：TP391

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1009-671X(2015)06-001-05

doi:10.11991/yykj.201503027