解偉男，李清華，屈楨深，奚伯齊

（哈爾濱工業(yè)大學(xué)空間控制與慣性技術(shù)研究中心，哈爾濱150001）

地磁場(chǎng)為近地空間提供了一個(gè)天然的導(dǎo)航坐標(biāo)系，許多動(dòng)物能夠利用地磁信息進(jìn)行導(dǎo)航，如海龜[1]、鳥類[2]等。地磁匹配作為地磁導(dǎo)航中一類主要方式，由于其不需要對(duì)地球局部地磁場(chǎng)建模，而是在地磁數(shù)據(jù)庫(kù)中進(jìn)行線圖匹配，因此特別適用于低空的導(dǎo)航定位。此外，地磁導(dǎo)航具有無源、無輻射、誤差不隨時(shí)間累積等優(yōu)點(diǎn)，與慣性導(dǎo)航組合可以實(shí)現(xiàn)兩種導(dǎo)航方式的優(yōu)勢(shì)互補(bǔ)，達(dá)到長(zhǎng)時(shí)間高精度自主導(dǎo)航的目的，具有重要的軍事意義[3-5]。

地磁日變是地磁場(chǎng)的一個(gè)基本組成部分[6-7]，其起源于大氣電離層中一直存在的電流體系，在地磁導(dǎo)航中相當(dāng)于地磁測(cè)量時(shí)的干擾。由于地磁日變的平均變化幅度約為幾十納特，因此會(huì)對(duì)地磁匹配導(dǎo)航定位精度和可靠性產(chǎn)生嚴(yán)重的影響[8]。然而地磁日變并未形成較為統(tǒng)一的數(shù)學(xué)模型，不同地點(diǎn)的地磁日變曲線有較大差異，同一地點(diǎn)的地磁日變曲線隨著季節(jié)不同也不盡相同，在地磁匹配過程中難以進(jìn)行有效的補(bǔ)償或剔除，因此只有少數(shù)學(xué)者考慮了地磁匹配中的地磁日變影響。文獻(xiàn)[9]指出了地磁日變是地磁匹配中的一個(gè)主要噪聲源，并通過數(shù)值仿真比較了不同的相關(guān)準(zhǔn)則對(duì)噪聲的適用性。文獻(xiàn)[7]采用FMI方法擬合地磁日變場(chǎng)，在地磁測(cè)量值中減去地磁日變場(chǎng)得到實(shí)時(shí)圖，該方法不需選取基線值，減小了人為誤差，所提取的日變更客觀、真實(shí)，然而該方法需以擬合當(dāng)天及前后共三天的數(shù)據(jù)為計(jì)算基礎(chǔ)，在實(shí)際應(yīng)用中實(shí)現(xiàn)困難。

本文在文獻(xiàn)[10]的基礎(chǔ)上，提出一種地磁日變影響下的地磁匹配算法。該算法根據(jù)地磁日變的特性，在每一段地磁匹配過程中，把地磁日變看成一個(gè)恒定的磁場(chǎng)偏差，在此基礎(chǔ)上構(gòu)建地磁日變影響下的均方差相關(guān)性準(zhǔn)則，引入匹配曲線的參數(shù)化模型，將地磁匹配問題轉(zhuǎn)化為非線性方程組的求解，最終通過快速迭代實(shí)現(xiàn)地磁日變影響下的匹配定位。該算法可以有效地消除地磁日變對(duì)匹配精度的影響，并具有較高的實(shí)時(shí)性。

1 地磁場(chǎng)的組成和地磁日變

地磁場(chǎng)根據(jù)其性質(zhì)不同可以劃分為穩(wěn)定磁場(chǎng)和變化磁場(chǎng)兩部分：

其中：T表示地磁總磁場(chǎng)；0T表示穩(wěn)定磁場(chǎng)；TΔ表示變化磁場(chǎng)。變化磁場(chǎng)相對(duì)于總磁場(chǎng)來說很弱，在特殊的磁暴條件下也只有總磁場(chǎng)的2%～4%，而多數(shù)情況僅有總磁場(chǎng)的1%以下。

穩(wěn)定磁場(chǎng)根據(jù)其性質(zhì)不同可以劃分為主磁場(chǎng)和異常場(chǎng)兩部分：

其中：mT表示主磁場(chǎng)；aT表示異常場(chǎng)。主磁場(chǎng)主要是由地核內(nèi)電流的對(duì)流形成，其強(qiáng)度約為穩(wěn)定磁場(chǎng)的94%；異常場(chǎng)主要是由地殼淺部具有磁性的巖石或礦石引起的，其強(qiáng)度約為穩(wěn)定磁場(chǎng)的6%。

變化磁場(chǎng)根據(jù)其變化特點(diǎn)又可分為平靜變化和擾動(dòng)變化。平靜變化起源于電離層中一直存在著的穩(wěn)定電流體系，其特點(diǎn)是變化連續(xù)，持續(xù)存在，并具有一定周期特性，平靜變化包括太陽靜日變化Sq和太陰日變化 L。擾動(dòng)變化起源于太陽噴射出來的帶電粒子流所產(chǎn)生的多種短暫的電流體系的變化而產(chǎn)生的，其特點(diǎn)是缺乏長(zhǎng)期連續(xù)性，出現(xiàn)時(shí)間不規(guī)則，變化形態(tài)復(fù)雜，擾動(dòng)變化包括磁暴、亞磁暴、地磁脈動(dòng)等。

地磁匹配導(dǎo)航依據(jù)不同地點(diǎn)的地磁場(chǎng)強(qiáng)度差異進(jìn)行匹配定位。雖然變化磁場(chǎng)強(qiáng)度只占地磁總場(chǎng)強(qiáng)度很小一部分，但是變化磁場(chǎng)強(qiáng)度相對(duì)于不同地點(diǎn)的地磁變化梯度來說卻很巨大，會(huì)嚴(yán)重影響匹配定位精度或產(chǎn)生誤匹配，因此在匹配定位時(shí)必須予以考慮。擾動(dòng)變化是一種復(fù)雜、不規(guī)則的短期磁干擾，應(yīng)避免在擾動(dòng)變化條件下進(jìn)行地磁導(dǎo)航，因此本文不考慮擾動(dòng)變化，只考慮平靜變化中的太陽靜日變化Sq和太陰日變化L。

太陽靜日變化Sq簡(jiǎn)稱靜日變化或日變，變化周期為1個(gè)太陽日（24 h），其幅度約為幾十nT。太陰日變化L變化周期為1個(gè)太陰日（約為25 h），其幅度很微弱，僅有1～2nT，很難將太陰日變化從太陽靜日變化中分離出來，因此通常將太陰日變化包括在太陽靜日變化之內(nèi)。

太陽靜日變化Sq具有如下特點(diǎn)：

1）Sq場(chǎng)隨著地方時(shí)變化而周期改變；

2）Sq場(chǎng)白天變化劇烈，夜晚變化平緩；

3）Sq場(chǎng)隨著緯度變化而改變；

4）Sq場(chǎng)隨著季節(jié)變化而改變，夏季變化幅度較大，冬季變化幅度較??；

5）Sq場(chǎng)的不同分量在赤道兩側(cè)對(duì)稱分布或反對(duì)稱分布。

2 地磁日變影響下的地磁匹配算法

地磁日變是地磁場(chǎng)的一個(gè)基本組成部分。地磁日變并未形成較為統(tǒng)一的數(shù)學(xué)模型，不同地點(diǎn)的地磁日變曲線有較大的差異，同一地點(diǎn)的地磁日變曲線隨著季節(jié)不同也不盡相同，因此如何實(shí)現(xiàn)地磁日變影響下的地磁匹配定位具有重要的意義。

2.1 問題描述

地磁匹配可以表示為平面上兩條曲線的相關(guān)性計(jì)算，如圖1所示的地磁匹配示意圖：曲線Rf表示慣導(dǎo)系統(tǒng)輸出的參考軌跡；曲線M為匹配算法求得的匹配軌跡；曲線Rf′與匹配軌跡M 平行，參考軌跡Rf與曲線Rf′之間的夾角為α，即航向誤差角。從圖中可以看出參考軌跡與匹配軌跡間存在著平移和旋轉(zhuǎn)，對(duì)應(yīng)著慣導(dǎo)系統(tǒng)的初始位置誤差和初始航向誤差。令點(diǎn)為參考軌跡上的第i個(gè)采樣點(diǎn)，其對(duì)應(yīng)的匹配軌跡點(diǎn)為(u ,v )T，其中 a 、u表示點(diǎn)的經(jīng)度，b、v表示點(diǎn)的緯度。

圖1 地磁匹配示意圖Fig.1 Schematic diagram of geomagnetic matching

2.2 性能指標(biāo)

若采用均方差準(zhǔn)則（MSD）計(jì)算曲線的相關(guān)性，需計(jì)算匹配軌跡曲線上各點(diǎn)所對(duì)應(yīng)的地球磁場(chǎng)特征值與真實(shí)磁場(chǎng)測(cè)量值之差的平方和的均值。然而如前文所述，若存在地磁日變時(shí)，該相關(guān)性計(jì)算方法并不能真實(shí)反映曲線的相關(guān)性，需要對(duì)地磁日變進(jìn)行補(bǔ)償。根據(jù)地磁日變特性可知，變化磁場(chǎng)TΔ是時(shí)間的函數(shù)。在不考慮擾動(dòng)變化的條件下，相對(duì)于地磁匹配過程中的每一段匹配周期而言，地磁日變變化緩慢。因此，在每一段軌跡的匹配過程中，地磁日變相當(dāng)于在匹配局部區(qū)域疊加了一個(gè)近似恒定的磁場(chǎng)偏差，由此可得地磁日變影響下的MSD相關(guān)性準(zhǔn)則：

考慮到參考軌跡在匹配軌跡附近，即可將 I(ui, vi)泰勒展開：

將泰勒展開式(4)帶入式(3)，并忽略高階小項(xiàng) O2，可以得到如式(5)的相關(guān)性約束：

在不引起混淆的情況下，可將相關(guān)性約束(5)寫成如下簡(jiǎn)化形式：

其中，

與文獻(xiàn)[10]相比，相關(guān)性約束(6)中不僅存在 ui和vi兩組未知變量，而且存在未知的地磁日變補(bǔ)償場(chǎng)δ，因此式(6)為地磁日變影響下的相關(guān)性約束。只要找到相應(yīng)的 ui、vi、δ使相關(guān)性約束(6)取極小值即可求得地磁日變影響下的匹配結(jié)果。

為了實(shí)現(xiàn)匹配結(jié)果的快速求解，建立了匹配軌跡的參數(shù)化模型，并將該模型引入相關(guān)性約束，從而將遍歷搜索方法轉(zhuǎn)化為快速迭代算法。如圖1所示，當(dāng)慣導(dǎo)系統(tǒng)存在初始位置誤差和初始航向誤差時(shí)，匹配曲線M的參數(shù)化模型為

將式(8)帶入相關(guān)性約束(6)，可以得到地磁日變影響下的均方差約束指標(biāo)函數(shù)：

因此地磁匹配問題可以轉(zhuǎn)化為尋求Δx、Δy、α和δ，使指標(biāo)函數(shù)(9)取極小值。

2.3 匹配求解

為使指標(biāo)函數(shù)(9)取極小值，可分別將指標(biāo)函數(shù)對(duì)自變量xΔ、yΔ、α、δ求一階偏導(dǎo)數(shù)，并令其為零，即：

將均方差約束指標(biāo)函數(shù)(9)帶入式(10)可得：

因此地磁匹配可以轉(zhuǎn)化為由xΔ、yΔ、α、δ組成的非線性方程組(15)的求解問題，該非線性方程組可通過Broyden迭代法求解。匹配結(jié)果可以通過式(8)計(jì)算得到。

3 迭代計(jì)算

地磁日變影響下的地磁匹配算法的實(shí)現(xiàn)步驟如下：

第1步：將參考軌跡離散化成N個(gè)點(diǎn)，對(duì)應(yīng)離散點(diǎn)的坐標(biāo)為 ( ai, bi)T， i = 1 ,2,…,N ；讀取載體在離散點(diǎn)(a,b )T時(shí)磁傳感器實(shí)時(shí)測(cè)量得到的地磁特征值

第2步：根據(jù)參考軌跡離散化的N個(gè)位置，從預(yù)先存儲(chǔ)的地磁數(shù)據(jù)庫(kù)中讀取該位置的地磁特征值，并根據(jù)地磁數(shù)據(jù)庫(kù)信息計(jì)算該位置的地磁特征值梯度信息 Ix,i和 Iy,i。

第3步：初始化迭代變量

第4步：根據(jù)式(17)～(19)解算迭代參數(shù)G、F和H：

其中，

第5步：計(jì)算位置誤差的增量、航向誤差的增量和地磁日變誤差的增量δM：

第6步：更新位置誤差、航向誤差和地磁日變誤差M：

第7步：判斷是否滿足終止迭代條件，若滿足則停止迭代并跳到第10步，否則跳到第8步。

終止迭代的條件有兩個(gè)，滿足任意一個(gè)即終止迭代：1）迭代次數(shù)達(dá)到預(yù)設(shè)的最大迭代次數(shù)；2）迭代誤差增量δM 的2范數(shù)小于設(shè)定值，即

其中，0ε＞為預(yù)先給定的迭代最小誤差。

第8步：根據(jù)更新后的M 計(jì)算參數(shù)K和δK：

第9步：根據(jù)公式(37)和公式(38)更新迭代變量H和G，然后跳到第5步。

第10步：根據(jù)迭代計(jì)算得到位置誤差、航向誤差和地磁日變誤差xΔ、yΔ、α、δ，將所得到結(jié)果代入匹配曲線的參數(shù)化模型(8)即得到匹配軌跡。

根據(jù)上述匹配步驟可知，迭代計(jì)算過程中地磁特征值梯度信息Ix,i和Iy,i并不改變，因此可以將匹配區(qū)域地磁特征值梯度信息事先存入載體計(jì)算機(jī)中，匹配計(jì)算前直接在載體計(jì)算機(jī)中檢索得到，也可以通過檢索載體計(jì)算機(jī)中的地磁特征值信息后計(jì)算得到。

需要指出的是，地磁日變影響下的地磁匹配算法假設(shè)參考軌跡位于匹配軌跡附近，在忽略高階小項(xiàng)2O的基礎(chǔ)上，以一階差分的形式近似相關(guān)性約束。然而當(dāng)慣導(dǎo)系統(tǒng)的初始誤差較大時(shí)，忽略的高階小項(xiàng)2O會(huì)嚴(yán)重影響算法的精度。為了克服上述問題，可以采用多次迭代的方式提高匹配算法的精度。在多次迭代時(shí)，可以將當(dāng)前的迭代結(jié)果作為下一次迭代的參考軌跡。

4 仿真實(shí)驗(yàn)

為了驗(yàn)證上述算法的有效性，本節(jié)對(duì)算法進(jìn)行仿真驗(yàn)證。根據(jù)地球磁場(chǎng)異常網(wǎng)格（EMAG2）繪制某局部區(qū)域地磁異常圖作為背景磁圖，其中地球磁場(chǎng)異常網(wǎng)格間距為2′，分辨率為0.01 nT，主要仿真參數(shù)設(shè)定如表1所示。此外設(shè)定載體初始經(jīng)度為127.00°，初始緯度為45.05°，以50 m/s的速度向北運(yùn)動(dòng)，慣性導(dǎo)航系統(tǒng)的解算頻率為100 Hz，地磁場(chǎng)參數(shù)的測(cè)量頻率為0.025 Hz。初始誤差假設(shè)為：經(jīng)度誤差0.014°，緯度誤差0.014°，航向誤差3°。圖2給出了真實(shí)軌跡和慣導(dǎo)系統(tǒng)輸出的參考軌跡。

表1 仿真參數(shù)Tab.1 Main parameters of the simulation

圖2 真實(shí)軌跡和參考軌跡Fig.2 The real track and the reference track

仿真時(shí)首先假設(shè)不存在地磁日變，圖3和圖4分別給出了仿真結(jié)果的經(jīng)度誤差曲線和緯度誤差曲線，其中算法 1為文獻(xiàn)[10]中二次迭代的基于迭代計(jì)算地磁輪廓線匹配算法，算法2為本文所提出的地磁日變影響下的地磁匹配算法，算法3為本文所提出的二次迭代的地磁日變影響下的地磁匹配算法。為了驗(yàn)證本文所提出算法對(duì)地磁日變的有效性，分別在磁場(chǎng)測(cè)量數(shù)據(jù)中加入10 nT和50 nT地磁日變。圖5～6和圖7～8分別給出不同地磁日變影響下的誤差曲線。表2給出了不同仿真條件下的匹配結(jié)果。根據(jù)匹配結(jié)果可以看出：在沒有日變影響時(shí)，本文所提出的算法與文獻(xiàn)[10]算法的匹配精度相當(dāng)；當(dāng)?shù)卮湃兆優(yōu)?0 nT時(shí)，本文所提出算法的匹配誤差僅為文獻(xiàn)[10]算法匹配誤差的27.3%；當(dāng)?shù)卮湃兆優(yōu)?0 nT時(shí)，文獻(xiàn)[10]所提出的算法發(fā)生了誤匹配，而本文所提出的算法由于在相關(guān)性準(zhǔn)則中考慮了地磁日變的影響，依然能夠保證較高的匹配定位精度。

圖3 無地磁日變時(shí)匹配經(jīng)度誤差曲線Fig.3 Longitude errors of matching results without geomagnetic diurnal variation

圖4 無地磁日變時(shí)匹配緯度誤差曲線Fig.4 Latitude errors of matching results without geomagnetic diurnal variation

圖5 10nT地磁日變時(shí)匹配經(jīng)度誤差曲線Fig.5 Longitude errors of matching results with geomagnetic diurn al variation10 nT

圖6 10 nT地磁日變時(shí)匹配緯度誤差曲線Fig.6 Latitude errors of matching results with geomagnetic diurnal variation10nT

圖7 50nT地磁日變時(shí)匹配經(jīng)度誤差曲線Fig.7 Longitude errors of matching results with geomagnetic diurnal variation50nT

圖8 50nT地磁日變時(shí)匹配緯度誤差曲線Fig.8 Latitude errors of matching results with geomagnetic diurnal variation50nT

由匹配結(jié)果可以看出，地磁日變影響下的地磁匹配算法可以有效補(bǔ)償?shù)卮湃兆儗?duì)匹配精度的影響，實(shí)現(xiàn)地磁日變影響下的高精度地磁匹配。此外，該算法采用多次迭代，匹配精度可以進(jìn)一步提高。利用Matlab中的tic和toc語句計(jì)算匹配算法的運(yùn)行時(shí)間，采用二次迭代的地磁日變影響下的地磁匹配算法在仿真計(jì)算機(jī)中總耗時(shí)為 16 ms。因此所提出的算法具有較高的實(shí)時(shí)性，可以保證地磁匹配導(dǎo)航的在線計(jì)算。

表2 匹配結(jié)果Tab.2 Matching results

5 結(jié) 論

本文提出了一種地磁日變影響下的地磁匹配新方法。該方法通過建立地磁日變影響下的相關(guān)性準(zhǔn)則函數(shù)，補(bǔ)償?shù)卮湃兆儗?duì)匹配結(jié)果的影響。仿真結(jié)果表明：在無地磁日變時(shí)，所提出算法的匹配誤差與文獻(xiàn)[10]相當(dāng)；當(dāng)?shù)卮湃兆優(yōu)?0 nT時(shí)，所提出算法的最大匹配誤差為文獻(xiàn)[10]的27.3%；當(dāng)?shù)卮湃兆優(yōu)?0 nT時(shí)，文獻(xiàn)[10]發(fā)生誤匹配，而本文所提出算法仍然可以實(shí)現(xiàn)高精度的匹配。此外，所提出的算法在Matlab仿真中的總耗時(shí)為16 ms。因此提出的地磁日變影響下的地磁匹配算法可以有效地補(bǔ)償?shù)卮湃兆儗?duì)匹配結(jié)果的影響，實(shí)現(xiàn)高精度的實(shí)時(shí)地磁匹配。

（

）：

[1] Brothers J R,Lohmann K J.Evidence for geomagnetic imprinting and magnetic navigation in the natal homing of sea turtles[J].Current Biology,2015,25(3):392-396.

[2] Qi X K,Ye D X,Sun Y Z,et al.Simulations to true animals’long-distance geomagnetic navigation[J].IEEE Transactions on Mag-netics,2017,53(1):1-8.

[3] Guo C F,Cai H,van der Heijden G H M.Feature extraction and geo-magnetic matching[J].The Journal of Navigation,2013,66(6):799-811.

[4] Wang P,Wu Y X,Hu X P,et al.Geomagnetic aided navigation suitability evaluation based on principal component analysis[C]//Inter-national Conference on Industrial Control and Electronics Engineering.2012:324-329.

[5] Xie W N,Qu Z S,Li Q H.A fast algorithm of the geomagnetic correlation matching based on MSD[C]//The3rd International Conference on Control,Automation and Systems Engineering.2013:59-62.

[6] Nizam M,Abdullah E J.Nanosatellite navigation with the WMM2005geomagnetic field model[J].Turkish Journal of Engineering&Envi-ronmental Sciences,2006,30:43-55.

[7] 鄧翠婷.地磁匹配中的地磁日變效應(yīng)與修正[D].南京:南京航空航天大學(xué),2013.Deng C T.Research on affection and modification of magnetic diurnal variation for geomagnetic matching[D].Nanjing:Nanjing University of Aeronautics and Astronautics,2013.

[8] 黃朝艷,鄧翠婷,趙華.地磁匹配中的地磁日變效應(yīng)[J].測(cè)繪工程,2013,22(4):17-21.Huang C Y,Deng C T,Zhao H.The affection of magnetic diurnal variation in geomagnetic matching[J].Engineering of Surveying and Mapping,2013,22(4):17-21.

[9] 喬玉坤,王仕成,張金生,等.基于相關(guān)分析的飛行器地磁匹配算法仿真研究[J].地球物理學(xué)進(jìn)展,2009,24(2):721-727.Qiao Y K,Wang S C,Zhang J S,et al.Simulation of geomagnetic matching algorithm based on correlation analysis for flying objects[J].Progress in Geophysics,2009,24(2):721-727.

[10] 解偉男,李清華,奚伯齊,等.基于迭代計(jì)算的地磁輪廓線匹配算法[J].中國(guó)慣性技術(shù)學(xué)報(bào),2015,23(5):631-635.Xie W N,Li Q H,Xi B Q,et al.Geomagnetic contour matching algorithm based on iterative method[J].Journal of Chinese Inertial Technology,2015,23(5):631-635.