嚴(yán)羽靈，唐清善，白 創(chuàng)

(長(zhǎng)沙理工大學(xué) 物理與電子科學(xué)學(xué)院，湖南 長(zhǎng)沙 410114)

地磁場(chǎng)是屬于地球的天然資源[1]。大部分動(dòng)物能夠利用地磁場(chǎng)提供的地磁信息進(jìn)行導(dǎo)航[2-3]。包括鯨魚(yú)、海龜、候鳥(niǎo)在內(nèi)的遷徙動(dòng)物[4]受各種氣候環(huán)境的影響還可測(cè)定精確的位置。地磁導(dǎo)航包括地磁場(chǎng)模型[5]建立、地磁場(chǎng)實(shí)時(shí)測(cè)量和地磁導(dǎo)航算法[6]。地磁導(dǎo)航算法可以分為地磁匹配算法[7]和地磁濾波算法。地磁濾波算法雖然實(shí)時(shí)性強(qiáng)，但是地磁場(chǎng)相對(duì)地理位置具有較強(qiáng)的非線性，導(dǎo)致濾波容易發(fā)散。地磁匹配算法[8-9]原理簡(jiǎn)單，且匹配精度較高，在地磁導(dǎo)航中應(yīng)用較多。

目前比較常見(jiàn)的算法是地磁場(chǎng)等值線匹配(Magnetic Contour Matching，MAGCOM)算法和迭代最近等值線匹配(Iterative Closest Contour Point，ICCP)算法[10]。MAGCOM 算法只能修正位移誤差，不能校正角度誤差[11]。而 ICCP 算法[12-13]只能修正角度誤差和位移誤差，不能修正伸縮誤差，因此精度不高。本文提出了一種基于仿射變換的地磁匹配定位算法。仿射變換[14-15]比較常見(jiàn)的有旋轉(zhuǎn)、平移、伸縮變換、反射和剪切。由于慣導(dǎo)系統(tǒng)輸出的軌跡與真實(shí)軌跡之間只存在位移、角度和伸縮誤差，所以只考慮前面3種變換。為了實(shí)現(xiàn)該算法，本文改進(jìn)了第12代國(guó)際地磁參考場(chǎng)模型，并輸出地磁數(shù)據(jù)庫(kù)。再通過(guò)MATLAB對(duì)參考軌跡、真實(shí)軌跡以及該算法進(jìn)行了仿真。然后通過(guò)Mathematica中的迭代和數(shù)值求解函數(shù)解非線性方程組，并將解輸入到MATLAB中進(jìn)行匹配定位，將各組解的匹配效果進(jìn)行了比較。

1 仿射變換算法原理

1.1 問(wèn)題描述

慣導(dǎo)系統(tǒng)存在隨時(shí)間累加的定位誤差，導(dǎo)致慣導(dǎo)系統(tǒng)所輸出的軌跡與真實(shí)軌跡之間存在位移、角度和伸縮誤差。該情況可以用如圖1所示的軌跡示意圖來(lái)描述。

圖1 軌跡示意圖

其中，曲線Rf表示捷聯(lián)慣導(dǎo)系統(tǒng)輸出的參考軌跡；曲線M是經(jīng)過(guò)仿射變換算法得到的匹配軌跡，其位于真實(shí)軌跡Re附近；作曲線M的平行線M′。從圖1可以看出匹配軌跡與參考軌跡存在角度誤差α。除了角度誤差，還存在初始位置誤差和伸縮誤差，Δx和Δy分別表示初始經(jīng)度、緯度誤差；k1和k2表示經(jīng)度、緯度方向的伸縮誤差。令(a,b)T是參考軌跡上的一點(diǎn)(a表示經(jīng)度，b表示緯度)，參考軌跡的初始點(diǎn)坐標(biāo)為(a1,b1)T，那么(a,b)T對(duì)應(yīng)的匹配軌跡上的點(diǎn)是(u,w)T。

1.2 仿射變換模型

仿射變換是二維坐標(biāo)到二維坐標(biāo)之間的線性變換，能夠保持二維圖形的“平直性”和“平行性”，即直線經(jīng)變換后還是直線，圓弧經(jīng)變換后還是圓弧，且二維圖形之間的相對(duì)位置關(guān)系不會(huì)變。它的變化包括旋轉(zhuǎn)、平移、伸縮、反射和剪切。本文針對(duì)慣導(dǎo)系統(tǒng)匹配定位，根據(jù)具體問(wèn)題只考慮前面3種。圖2所示為仿射變換-旋轉(zhuǎn)變化圖。

圖2 仿射變換-旋轉(zhuǎn)變化圖

從圖2可知，原始坐標(biāo)系下的一個(gè)點(diǎn)是(a,b)，要對(duì)其進(jìn)行旋轉(zhuǎn)操作。在基于原點(diǎn)的情況下，通過(guò)旋轉(zhuǎn)坐標(biāo)軸就能得到旋轉(zhuǎn)之后的點(diǎn)的坐標(biāo)，并得到新的坐標(biāo)系。在這個(gè)坐標(biāo)系中通過(guò)立體幾何關(guān)系確定點(diǎn)(a,b)在新坐標(biāo)系中的坐標(biāo)，在新坐標(biāo)系中它的X、Y坐標(biāo)分別為(acosα+bsinα)和(bcosα-asinα)。然后在這個(gè)位置的基礎(chǔ)上加上其在X軸和Y軸的偏移量Δx、Δy。又因點(diǎn)(a,b)只是參考軌跡上的一個(gè)點(diǎn)，除了旋轉(zhuǎn)和平移，參考軌跡與匹配軌跡之間還存在伸縮變換關(guān)系。由此可以建立如下所示的仿射變換模型

(1)

1.3 地磁匹配

地磁匹配問(wèn)題可表示為平面上兩條曲線的相關(guān)性計(jì)算，采用平方差(Squared Difference，SD)算法計(jì)算曲線之間的相關(guān)性，即計(jì)算匹配軌跡上各點(diǎn)對(duì)應(yīng)的地磁場(chǎng)特征值與實(shí)際測(cè)量得到的地磁場(chǎng)特征值之差的平方和。因?yàn)樵撝翟叫t相差度越小，所以該值達(dá)到最小時(shí)對(duì)應(yīng)的軌跡序列就是最優(yōu)的匹配位置。那么SD相關(guān)性指標(biāo)函數(shù)表示如下

(2)

式中，I(ui,wi)表示點(diǎn)(ui,wi)T所在位置對(duì)應(yīng)的地磁數(shù)據(jù)庫(kù)中地磁特征值；Ir(ai,bi)表示載體在參考點(diǎn)(ai,bi)T時(shí)，磁強(qiáng)計(jì)實(shí)時(shí)測(cè)量得到的地磁特征值，N表示參考軌跡采樣點(diǎn)的總個(gè)數(shù)。

由于匹配軌跡位于真實(shí)軌跡附近，所以匹配結(jié)果點(diǎn)所在的地磁特征值約等于對(duì)應(yīng)的真實(shí)軌跡上的點(diǎn)的地磁特征值[16]，即存在式(3)。

I(ui,wi)≈Ir(ai,bi)

(3)

為了實(shí)現(xiàn)匹配，就要求參考軌跡在匹配軌跡附近，即可將I(ui,wi)在參考點(diǎn)(ai,bi)T處泰勒展開(kāi)，稱I(ui,wi)在點(diǎn)(ai,bi)T的一階泰勒公式，如式(4)所示。

(4)

其中，I(ai,bi)表示在點(diǎn)(ai,bi)T時(shí)所對(duì)應(yīng)的地磁數(shù)據(jù)庫(kù)中地磁特征值；?I(ai,bi)/?x表示地磁特征值對(duì)經(jīng)度方向的梯度在點(diǎn)(ai,bi)T上的取值；?I(ai,bi)/?y表示地磁特征值對(duì)緯度方向的梯度在點(diǎn)(ai,bi)T上的取值；O2表示的是泰勒展開(kāi)后所有大于一階的高階小項(xiàng)。

聯(lián)立式(3)和式(4)，并且忽略高階小項(xiàng)O2,就可得如下式

(5)

將上式簡(jiǎn)寫(xiě)成如下式

Ix,i(ui-ai)+Iy,i(wi-bi)+It,i≈0

(6)

其中，Ix,i=?I(ai,bi)/?x,Iy,i=?I(ai,bi)/?y,It,i=?I(ai,bi)-Ir(ai,bi)。

將式(1)代入式(6)可得

(7)

那么平方差相關(guān)性指標(biāo)函數(shù)如下

(8)

其中，a′i=ai-a1；b′i=bi-b1。

因此地磁匹配問(wèn)題可以轉(zhuǎn)化為求取仿射模型參數(shù)Δx、Δy、α、k1、k2的值，以便使得式(8)達(dá)到極小值?？梢詫⑵鋵?duì)仿射模型參數(shù)Δx、Δy、α、k1和k2求一階偏導(dǎo)數(shù)，并使其結(jié)果為零，如下所示。

(9)

把式(8)代入式(9)可得如下計(jì)算式。

(10)

Δy)+Iy,iIt,i]=0

(11)

(-a′icosα-b′isinα)]×[Ix,i(k1(a′icosα+b′isinα)-

a′i+Δx)+Iy,i(k2×(-a′isinα+b′icosα)-b′i+Δy)+

It,i]=0

(12)

(13)

(14)

以上5個(gè)計(jì)算式組成一個(gè)以仿射參數(shù)Δx、Δy、α、k1和k2為變量的非線性方程組

(15)

因此問(wèn)題又轉(zhuǎn)化為解非線性方程組(15)，在本文中通過(guò)Mathematica對(duì)它進(jìn)行求解，然后將計(jì)算所得的5個(gè)誤差變量Δx、Δy、α、k1和k2代入仿射變換模型式(1)即可得到匹配結(jié)果。

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

2.1 地磁數(shù)據(jù)庫(kù)的建立

為了建立匹配區(qū)域的地磁數(shù)據(jù)庫(kù)，本文改進(jìn)了第12代國(guó)際地磁參考場(chǎng)模型。給定經(jīng)緯度的精度、起始點(diǎn)和終點(diǎn)的經(jīng)緯度、年份、海拔高度、地點(diǎn)，可以輸出起點(diǎn)和終點(diǎn)之間所有點(diǎn)的地磁特征值。地磁特征值包括磁偏角、磁傾角、總磁場(chǎng)強(qiáng)度以及北向、東向、垂直方向的磁場(chǎng)強(qiáng)度。本文選取總磁場(chǎng)強(qiáng)度作為地磁特征值。

2.2 仿射變換算法驗(yàn)證

通過(guò)MATLAB模擬參考軌跡和真實(shí)軌跡，首先對(duì)地磁數(shù)據(jù)庫(kù)進(jìn)行網(wǎng)格化處理，并計(jì)算出數(shù)據(jù)庫(kù)里面各個(gè)點(diǎn)的梯度。然后設(shè)定參考軌跡的初始點(diǎn)為(110°,30°)，沿經(jīng)度和緯度方向的初始速度分別為120 m·s-1和160 m·s-1，每隔50 s取點(diǎn)并給速度加幅值為1 m·s-1粉噪聲，再經(jīng)過(guò)雙線性插值，將會(huì)得到軌跡上各點(diǎn)的經(jīng)緯度、梯度及相應(yīng)的點(diǎn)在地磁數(shù)據(jù)庫(kù)中檢索得到的地磁特征值。最后令真實(shí)軌跡的初始點(diǎn)和終點(diǎn)分別為(110°,30°)和(111°,31°)，總速度為200 m·s-1，對(duì)其正交分解得到分速度，每隔50 s取點(diǎn)，再使用雙線性插值輸出軌跡序列點(diǎn)和這些點(diǎn)在庫(kù)中的地磁特征值，并給地磁特征值加幅值為3 nT的白噪聲作為測(cè)量值。

將以上相關(guān)數(shù)據(jù)作為輸入，通過(guò)Mathematica解非線性方程組(15)，采用了兩種方法求解：迭代求解和數(shù)值求解。表1是關(guān)于使用這兩種方法計(jì)算所得的仿射參數(shù)，序號(hào)1是迭代法求得的解，序號(hào)2～13是數(shù)值求解的結(jié)果。

表1 仿射參數(shù)

將表1的仿射參數(shù)Δx、Δy、α、k1、k2輸入到MATLAB中，輸出軌跡圖和誤差圖。最終實(shí)驗(yàn)結(jié)果表明，迭代法所求解的匹配效果最好，圖3～圖5分別是它的軌跡圖、經(jīng)度誤差圖和緯度誤差圖。由圖得出，慣導(dǎo)輸出參考軌跡的最終定位點(diǎn)與真實(shí)軌跡的終點(diǎn)相差6 328.5 m，參考軌跡與真實(shí)軌跡之間的經(jīng)緯度誤差分別為-0.034°和0.052°；匹配軌跡的最終的定位點(diǎn)與真實(shí)軌跡的終點(diǎn)相差293.3 m，而且匹配軌跡與真實(shí)軌跡之間的經(jīng)緯度誤差分別為0.003°、0.000 4°。由此可知，相對(duì)于慣導(dǎo)系統(tǒng)，匹配軌跡的定位精度明顯提高。

圖3 軌跡圖

圖4 經(jīng)度誤差圖

圖5 緯度誤差圖

3 結(jié)束語(yǔ)

本文針對(duì)慣導(dǎo)系統(tǒng)的定位誤差隨時(shí)間累積的問(wèn)題，提出了一種基于仿射變換的地磁匹配定位算法。該算法將地磁匹配問(wèn)題表示為平面上兩條曲線的相關(guān)性計(jì)算，采用了SD算法計(jì)算曲線之間的相關(guān)性，即計(jì)算匹配軌跡上各點(diǎn)對(duì)應(yīng)的地磁場(chǎng)特征值與真實(shí)測(cè)量得到的地磁場(chǎng)特征值之差的平方和。依據(jù)相關(guān)性準(zhǔn)則，該值達(dá)到最小時(shí)，其對(duì)應(yīng)的軌跡序列就是最優(yōu)的匹配位置。實(shí)驗(yàn)表明，匹配軌跡的最終定位點(diǎn)與真實(shí)軌跡終點(diǎn)相比誤差較小，為慣導(dǎo)系統(tǒng)輸出的軌跡誤差的4.63%。因此本文提出的基于仿射變換的地磁匹配定位算法可以有效修正慣導(dǎo)系統(tǒng)的初始位置誤差、初始航向誤差和初始速度誤差，實(shí)現(xiàn)精準(zhǔn)的匹配定位。