基于多重分形克里金的逐步插值校正法構建局部地磁基準圖

趙玉新, 邢文

哈爾濱工程大學 自動化學院，黑龍江 哈爾濱 150001

?

基于多重分形克里金的逐步插值校正法構建局部地磁基準圖

趙玉新, 邢文

哈爾濱工程大學 自動化學院，黑龍江 哈爾濱 150001

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

關鍵詞：地磁導航；克里金法；逐步插值校正；分形理論；多重分形克里金插值法；地磁基準圖

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

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

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

針對上述問題，文中提出一種基于多重分形克里金的逐步插值校正法，實驗表明，該方法構建的局部地磁基準圖精度較高，為后續將其應用于地磁匹配導航作保障。

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

1.1克里金法基本原理

克里金法是插值精度較高的空間插值方法之一。它是在平均殘差接近于零、估計誤差方差最小的前提下，將已有數據加權線性結合得到估計值，是一種最優、線性且無偏的估計。它通過以距離為自變量的變異函數計算權值，相對于其他空間插值方法而言，更加符合實際。克里金插值法在估計待插值點屬性值的過程中，不僅考慮了待插值點與鄰近已知點的空間位置關系，還考慮了各相鄰已知點之間的位置關系，充分利用了已知數據點的空間分布結構特征。但由于克里金插值的推導過程是建立在空間相關系數矩陣極小意義下的，這決定了它是一種滑動加權平均或低通濾波過程，因此不可避免地削弱了高頻、局部及弱信號，無法重現局部區域隆起或凹陷的奇異特征。

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

式中：Z*為待插值點x0處的屬性值；Zi為估計區域內1～n個采樣點的實測地磁異常場強度值；λi為待求權重系數。

1.2 多重分形理論

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

(1)

根據分維定義，擬合最佳逼近時，尺度與測度的關系可表示為

(2)

式中：a、b為常數。a為奇異系數，表征了待插值點小鄰域內的凹凸特性。在測度—尺度的雙對數坐標中，a是在最小方差意義下所擬合直線的斜率。

2多重分形克里金插值法

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

多重分形克里金插值法的原理為利用多重分形的理論對待插值點小鄰域內克里金插值結果進行奇異性校正，即用較大鄰域內地磁異常場強度的均值去估計中心小鄰域內的均值，作為待插值點的插值結果。

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

(3)

(4)

式(3)與式(4)聯立得

式中：Z*為N·δ鄰域內地磁異常場強度的均值，也是克里金法對基準數據的估值結果；Zδ為δ鄰域內地磁異常場強度的均值，也是待插值點的地磁異常場估計值；N為最大尺度與最小尺度的比值；N2-a為校正系數。

圖1　二維空間測度計算示意

2.2基于多重分形克里金插值法的仿真實驗

本文選取3個實驗區域對插值算法進行仿真實驗，并選取4個有效性評價指標對插值后的網格數據進行評估，如表1所示。由于在研究區域的邊界缺少參考點，以至于邊界處插值精度較低，因此文中在構建地磁基準圖的仿真過程中均是在考慮邊界效應的情況下進行的。經實驗驗證，第1周和第2周的參考點對邊界區域的構圖精度有著重要的影響，不可缺少。為了方便，將待插值數據的最外圍5周作為邊界區域。

表1　 有效性評價指標

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

表2　 區域4145、4053、4109B的經緯度范圍

對編號為4 053、4 145、4 109B數據的中心區域進行多重分形克里金插值，流程如圖2所示，比較多重分形克里金插值法與克里金插值法的構圖精度。

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

改變盒子數N的值，觀察其對插值精度的影響。以4 053數據為例，實驗結果見表3。

表3　 4 053不同盒子數下的多重分形克里金插值精度驗證

由表3的數據可以看出，多重分形克里金插值法的構圖精度不高，原因可能有2點：

1)現有數據較稀疏，在用局部去估計整體的過程中，不可避免地導致部分分形特征丟失。

2)在某些點上克里金插值的結果已經與基準數據很接近，再用校正系數去校正，就會增大結果的誤差。

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

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

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

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

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

表4　 文中方法的構圖精度

(a) 基準數據中心區域三維圖

(b) 克里金插值構建的地磁基準圖

(c) 文中方法所構地磁基準圖圖4　區域4053三維比較圖

分析表4的數據可知，文中方法的構圖精度優于克里金插值法的構圖精度。由圖4也可看出，基于逐步插值校正的多重分形克里金插值法對局部位置的刻畫較細致，能夠表現出更多的隆起或下陷等奇異特征，與基準數據的三維圖形態更加相近，具有良好的插值效果。

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

圖5　區域4 053中的剖線x=82

圖6　區域4 145中的剖線x=166

圖7　區域4 109B中的剖線x=106

從剖線圖中更能清晰地看出，文中方法所構地磁基準圖在細節上與基準數據曲線更加相似。文中方法是在克里金插值法描繪整體趨勢的基礎上，采用多重分形插值對局部位置逐步進行奇異性校正。

實驗表明，在數據非常稀疏的情況下，文中方法相對于常用插值算法而言仍具有較高的插值精度，如表5所示。表中數據的分辨率為0.016°，中心區域共7×7個采樣點。

表5　 編號4 053稀疏數據的插值精度比較

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

5結束語

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

參考文獻:

[1]喬玉坤, 王仕成, 張金生, 等. 采用矩諧分析和支持向量機的地磁導航基準圖構建方法[J]. 西安交通大學學報, 2010, 44(10): 47-51.

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

[3]喬玉坤, 王仕成, 張金生, 等. 基于BP網絡的地磁基準圖制備及其精度評價[J]. 中國慣性技術學報, 2009, 17(1): 53-58.

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

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

[6]王哲, 王仕成, 張金生, 等. 一種地磁匹配制導基準圖制備方法及其有效性評價[J]. 系統工程與電子技術, 2008, 30(11): 2207-2211.

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

[8]楊功流, 張桂敏, 李士心. 泛克里金插值法在地磁圖中的應用[J]. 中國慣性技術學報, 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]孫霞, 吳自勤. 規則表面形貌的分形和多重分形描述[J]. 物理學報, 2011, 50(11): 2126-2131.

[11]楊娟, 卞保民, 閆振綱, 等. 典型隨機信號特征參數統計分布的分形特性[J]. 物理學報, 2011, 60(10): 86-92.

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

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

網絡出版地址：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-)，男，教授，博士生導師；

基金項目：國家自然科學基金資助項目(51379049).

收稿日期：2015-03-25.網絡出版日期：2015-12-06.

中圖分類號：TP391

文獻標志碼：A

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

doi:10.11991/yykj.201503027