鄧小波，趙軍瑞，宋 勝

(1.中國人民解放軍61711部隊，新疆 喀什 844000；2.中國電子科技集團(tuán)公司第五十四研究所，河北 石家莊 050081)

0 引言

地磁場是地球本身固有的基本物理場[1]。在主地磁場基礎(chǔ)上，不同位置的地形和組成物質(zhì)不同，導(dǎo)致局部的異常場[2-4]，而且異常場與空間地理位置存在一定的對應(yīng)關(guān)系?；诋惓銎ヅ涞牡卮牌ヅ鋵?dǎo)航是一種無源、自主的導(dǎo)航方式，具有良好的應(yīng)用前景[5-8]。

在地磁導(dǎo)航算法的研究中，穆華等人[9]采用了擴(kuò)展卡爾曼濾波的方法研究了慣性/地磁組合導(dǎo)航技術(shù)在船舶導(dǎo)航方面的應(yīng)用，仿真結(jié)果證明慣性/地磁導(dǎo)航魯棒性好，可實現(xiàn)長航時工作；楊功流等[10]研究了基于Sigma點采樣的UKF方法構(gòu)建慣性/地磁卡爾曼濾波器，相比于使用EKF方法避免了非線性系統(tǒng)的線性化誤差；朱占龍等[11]研究了使用新息正交法的自適應(yīng)濾波方法，仿真結(jié)果證明該方法提高了組合導(dǎo)航系統(tǒng)對野值的適應(yīng)能力，提高了系統(tǒng)的可靠性；劉睿[12]研究了慣性/地磁/天文的組合導(dǎo)航方法，綜合了多種導(dǎo)航方式的優(yōu)點，構(gòu)造了長航時高精度的組合導(dǎo)航系統(tǒng)。上述研究主要針對地面、水下以及低空域進(jìn)行地磁匹配算法研究，然而高空高速飛行器由于地磁場環(huán)境和運動特性有所不同，需要對其導(dǎo)航方法進(jìn)行針對性研究。

地磁圖的準(zhǔn)確性和包含可用信息量的大小是影響地磁匹配導(dǎo)航精度的重要因素。而不同位置地磁異常場特性不同，能用于匹配導(dǎo)航的信息量也不同。因此，進(jìn)行地磁匹配導(dǎo)航前需要選擇適合匹配的區(qū)域[13-16]。同時，穿越地磁區(qū)域的航跡方向也會影響到地磁匹配的效果，這也是進(jìn)行載體航跡規(guī)劃的重要依據(jù)[17-18]。

本文根據(jù)對高空地磁圖特性的分析，利用高空大范圍地磁梯度呈現(xiàn)趨勢性的特點，提出了地磁方向熵的概念，由此給出了一種高空高速飛行器的地磁匹配軌跡選取方法。

1 高空地磁圖匹配的特性分析

目前，高空地磁異常場測量在實現(xiàn)上面臨較大困難，通常通過采集地面或低空地磁異常場數(shù)據(jù)，然后通過延拓方法得到高空地磁場的方案[19]。這使得高空地磁異常場的特性與地面或低空地磁異常場有所不同[20]。取1 000 m高度地磁圖和30 000 m高度地磁圖繪制三維曲面圖，如圖1所示。

(a) 1 000 m

(b) 30 000 m圖1 同一區(qū)域不同高度的地磁場三維曲面圖Fig.1 3D surface map of the geomagnetic field at different altitudes in the same area

通過對比圖1不同高度的地磁三維圖可以看出，高空地磁異常場的變化相比低空更加平滑。因此，高空地磁匹配算法需要適應(yīng)高空地磁圖的特點。地磁異常場的削弱意味著對地磁傳感器測量精度的要求比地面地磁匹配要高。而高空地磁場這種平滑有明顯趨勢性變化的特性也會對地磁匹配方法提出新的要求和需要解決的問題，即飛行器飛入方向和飛行軌跡問題。高空地磁異常圖的波動較小，起伏變化具有明顯的與地理位置有關(guān)的趨勢性，這使得飛行軌跡的航向和等值線變化趨勢的角度會影響到匹配的性能。飛行軌跡與地磁等值線處于近平行狀態(tài)時，匹配性能較好，而處于近垂直狀態(tài)時匹配誤差增大，且匹配概率大大降低?；谝陨霞夹g(shù)背景，提出地了磁梯度網(wǎng)絡(luò)和地磁方向熵結(jié)合的地磁匹配航跡選擇方法。

2 數(shù)字地磁圖的地磁梯度向量場

然而，實測的地磁圖一般都是離散的，以二維數(shù)組的形式存儲在計算機(jī)中。二維數(shù)組可以以更加簡單的方式計算地磁梯度，即使用地磁矩陣的數(shù)值梯度(Numerical Gradient)來計算數(shù)字地磁圖沿著x方向和y方向的地磁場梯度。

矩陣的數(shù)值梯度的定義為二維函數(shù)在x方向和y方向的偏導(dǎo)，其值為：

(1)

對于矩陣這種離散數(shù)據(jù)，可以通過計算差分的方式來取代求偏導(dǎo)。例如對于單位間距的矩陣A，其內(nèi)部的數(shù)據(jù)可以采用中心差分的方式計算梯度。內(nèi)部梯度值G(:,j)為：

G(:,j)=0.5×(A(:,j+1)-A(:,j-1))，

(2)

式中，j∈(2～N-1)，對于矩陣邊緣的數(shù)據(jù)，可以采用單邊差分的方法計算梯度：

(3)

通過差分計算的方法可以得到離散的數(shù)字地磁矩陣關(guān)于x方向和y方向的梯度分量。其向量和為當(dāng)前位置的地磁梯度向量，該向量的方向即為地磁場強度增大的方向。該向量的角度可利用x梯度與y梯度按照四象限反正切的方法計算出來：

(4)

式中，Y和X表示地磁場y方向和x方向地磁場梯度值;θ表示該向量按照北偏東為正方向的角度。將地磁梯度向量以箭頭的形式，在地磁圖內(nèi)繪制出來，就稱為向量場地磁圖，或向量網(wǎng)格地磁圖。以高空30 000 m地磁圖的一個區(qū)域為例，其梯度向量場地磁圖如圖2所示，向量場地磁圖中的網(wǎng)格數(shù)為100×100，橫縱坐標(biāo)分別表示網(wǎng)格的序號。因較為密集，圖中給出了部分區(qū)域的放大。

圖2 高空30 000 m區(qū)域的梯度向量場地磁圖Fig.2 Gradient vector magnetic field map at 30 000 m altitude

可以看出，向量場中的向量箭頭指向的是地磁場強度增大的反向，即由深藍(lán)色等值線變化到黃色等值線，與右側(cè)色標(biāo)相對應(yīng)。通過將式(1)中的四象限反正切公式改為：

(5)

由此畫出的地磁梯度網(wǎng)格如圖3所示。

圖3 高空30 000 m區(qū)域的平行等值線向量場地磁圖Fig.3 Parallel isoline vector magnetic field map at 30 000 m altitude

相比式(4)，式(5)求解出的θ結(jié)果相當(dāng)于逆時針旋轉(zhuǎn)90°，變?yōu)榕c地磁場強度梯度垂直、與地磁等值線平行的曲線；同樣，可以從圖3放大部分看出，由式(5)繪制出的向量場地磁圖中的箭頭與等值線方向平行。

隨著高空地磁場隨高度衰減，地磁等值線的變化趨勢也逐漸趨于平緩。換言之，地磁場等值線就是高空地磁圖中逐漸在大范圍內(nèi)趨于平行的直線。高空30 000 m的最大地磁信息熵區(qū)塊中地磁圖區(qū)塊的地磁圖等值線就符合這種特性，整個圖內(nèi)的等值線幾乎都沿著一個方向呈現(xiàn)平行分布的趨勢。因此，在該地磁區(qū)域內(nèi)，地磁匹配方法應(yīng)當(dāng)適應(yīng)這種特性。常用的MSD方法在航跡與等值線呈現(xiàn)過大角度，例如接近垂直于等值線時，平行于等值線方向的軌跡的慣導(dǎo)誤差較大，但這些誤差無法通過MSD方法匹配到真實位置附近，地磁匹配幾乎失效。而當(dāng)航跡平行于等值線方向時，慣導(dǎo)誤差就能通過MSD方法匹配到，匹配的精度就能大大提高。

3 地磁方向熵

向量場地磁圖表現(xiàn)了地磁數(shù)據(jù)中每個位置地磁場等值線的方向情況。實際上，為了衡量一個地磁圖區(qū)域內(nèi)整體性的地磁場等值線的復(fù)雜程度，可以構(gòu)建出地磁方向熵的概念，其公式為：

(6)

該式的物理意義為：通過式(4)和式(5)可以計算出向量場中每個地磁梯度或等值線平行矢量的方向(該方向以北為起點，北偏東順時針為正，范圍為-180°～180°)。將這些矢量的角度排序后成為矢量角度數(shù)據(jù)集，將-180°～180°的范圍劃分為K個等跨度區(qū)間。矢量角度數(shù)據(jù)集落在對應(yīng)K個等跨度區(qū)間內(nèi)，每個區(qū)間的數(shù)量為ni。按照式(6)計算得出的Hdir就是對應(yīng)地磁場區(qū)域的地磁方向熵。

從地磁方向熵的物理意義來看，地磁方向熵反映出地磁場內(nèi)等值線的變化復(fù)雜程度。地磁方向熵越小，地磁矢量角的分布越集中，這意味著該區(qū)域內(nèi)的地磁場等值線在大范圍內(nèi)呈現(xiàn)出平行直線的形狀，其角度方向統(tǒng)一。地磁方向熵越大，說明地磁場等值線在大范圍內(nèi)呈現(xiàn)出較為豐富的形狀變化，或大范圍內(nèi)地磁等值線形成閉環(huán)。對高空地磁圖來說，其本身伴隨著地磁等值線峰值的削弱，地磁等值線越平直，越適合降低地磁匹配算法的匹配誤差。

4 仿真與分析

仿真實驗在圖2所示的地磁區(qū)域內(nèi)進(jìn)行。通過地磁方向熵計算方法，將地磁區(qū)域內(nèi)的所有地磁矢量角劃分為36個等跨度區(qū)間，計算得到地磁方向熵值為0.126 6。這個值意味著地磁等值線中的地磁矢量角高度集中于跨度區(qū)間[40°，50°]，即和等值線近似平行的角度。根據(jù)預(yù)測，航向角屬于該區(qū)間內(nèi)的軌跡的匹配效果更好。為了證明這個預(yù)測結(jié)果，設(shè)計3條航跡，分別為平行等值線、呈45°和垂直等值線來驗證仿真結(jié)果。根據(jù)地磁方向上設(shè)定3條運動軌跡的飛入點都在(38.39°N,115.45°E)，飛行速度為3 400 m/s。地磁等值線的方向大約為北偏東45°，本文選擇的3條運動軌跡由飛入點進(jìn)入后的飛行方向與地磁等值線間的初始夾角分別選為：航跡① 0°(近平行狀態(tài))、航跡② 45°(中間狀態(tài))和航跡③ 90°(近垂直狀態(tài))；對應(yīng)的航向分別為：航跡① 北偏東50°、航跡② 北偏東90°和航跡③ 北偏東130°。3條運動軌跡如圖4所示。

圖4 3條預(yù)設(shè)仿真軌跡Fig.4 Three preset simulation tracks

與3條運動軌跡相對應(yīng)，設(shè)定的慣導(dǎo)系統(tǒng)初始誤差、慣性器件誤差和地磁誤差如表1所示。

表1 仿真誤差設(shè)置Tab.1 Simulation error setting

飛行器的仿真導(dǎo)航算法采用慣性/地磁匹配松組合，將MSD法計算得到的匹配位置輸入松組合卡爾曼濾波器，用以修正慣導(dǎo)軌跡。

飛行器分別按照以上3條運動軌跡運動，運動過程中匹配定位結(jié)果分別如圖5、圖6和圖7所示。

(a) 原始大小

(b) 部分放大圖5 軌跡1的匹配效果圖Fig.5 Matching result of track 1

(a) 原始大小

(b) 部分放大

(a) 原始大小

(b) 部分放大

圖5、圖6和圖7的(a)圖為匹配區(qū)域內(nèi)的航跡的匹配定位圖，(b)圖為相應(yīng)3幅圖中黑色細(xì)實線框內(nèi)小段軌跡的放大示意圖。各幅圖中，紅色圓圈為軌跡真值位置，綠色實線為純慣導(dǎo)輸出位置，黑色星號點為地磁匹配定位結(jié)果，藍(lán)色細(xì)實線為地磁匹配算法校正慣性導(dǎo)航后輸出的軌跡。

圖5中，軌跡①匹配效果良好，定位結(jié)果基本和真值重疊，匹配誤差為530.19m(1σ)。軌跡②匹配結(jié)果接近真實位置，但其定位誤差較大，匹配誤差為2 051.11m(1σ)。軌跡③地磁匹配定位結(jié)果逐漸偏離真值，呈發(fā)散趨勢，仿真中匹配定位誤差達(dá)到4 162.34m(1σ)。

5 結(jié)束語

針對高空地磁圖等值線呈現(xiàn)趨勢性和全局隨機(jī)性下降的現(xiàn)象，本文結(jié)合地磁信息熵與地磁方向熵2種量化指標(biāo)，研究了高空地磁匹配區(qū)域飛行航跡選取方法。由仿真可見，在軌跡與等值線呈近平行狀態(tài)時，其匹配定位誤差較小，隨著運動軌跡與等值線夾角增加，匹配定位誤差加大，且在近垂直狀態(tài)時其匹配結(jié)果可信度很差。數(shù)據(jù)處理結(jié)果說明了應(yīng)用地磁方向熵選擇飛行軌跡方向?qū)μ岣叩卮牌ヅ涠ㄎ痪染哂兄匾饔煤托Ч?/p>