基于磁梯度張量不變量的橢圓誤差消除方法*

袁 鵬,祗會(huì)強(qiáng),靳 鴻,段曉倩,孫赫軒

(中北大學(xué)儀器科學(xué)與動(dòng)態(tài)測(cè)試教育部重點(diǎn)實(shí)驗(yàn)室,太原 030051)

磁異常探測(cè)技術(shù)是未知爆炸物探測(cè)的重要手段,具有極高的軍事意義和民用價(jià)值[1-2]。磁探測(cè)初始階段測(cè)量磁場總量,后來又測(cè)量磁場分量和梯度,到現(xiàn)在發(fā)展到測(cè)量磁梯度張量。磁梯度張量測(cè)量技術(shù)和傳統(tǒng)的磁場測(cè)量相比有巨大的優(yōu)勢(shì),它的最大優(yōu)點(diǎn)是能夠有效的克服地磁場的影響,提高磁性目標(biāo)的定位準(zhǔn)確度,不足之處就是搭載平臺(tái)對(duì)其定位結(jié)果的影響較大。

磁梯度張量不變量是由磁梯度張量通過某些變換獲得的標(biāo)量,這些標(biāo)量不隨坐標(biāo)系的變化而變化,它適合于機(jī)動(dòng)平臺(tái)對(duì)磁性目標(biāo)的定位,目前已經(jīng)成為國內(nèi)外的研究熱點(diǎn)[3-5]。

磁梯度張量不變量定位方法,由于橢圓系數(shù)的存在,導(dǎo)致求解單位距離矢量和距離矢量模的過程中都存在一定的誤差,單位距離矢量是測(cè)量點(diǎn)指向磁源的單位方向矢量,距離矢量模是測(cè)量點(diǎn)到磁源的距離。針對(duì)該問題,提出了一種利用磁梯度張量不變量消除橢圓誤差的改進(jìn)定位方法。該方法利用單位定位矢量和定位矢量的模相乘,來求出定位矢量,本文的創(chuàng)新點(diǎn)是:①提出一種新的單位定位矢量計(jì)算方法,該方法主要利用測(cè)量點(diǎn)和磁性目標(biāo)形成的距離矢量與磁梯度張量矩陣的絕對(duì)值最小的特征值對(duì)應(yīng)的特征向量的垂直關(guān)系。②提出一種距離矢量模的修正方法,該方法利用磁梯度張量的特征值不變量來消除橢圓系數(shù),從而修正定位矢量的模。

本文所提出的基于磁梯度張量不變量的改進(jìn)定位方法,可以有效地避免地球磁場的影響,實(shí)現(xiàn)隱蔽磁性目標(biāo)的定位,而且適用于移動(dòng)平臺(tái),有很高的應(yīng)用價(jià)值。

1 磁梯度張量不變量定位方法

當(dāng)探測(cè)距離大于2.5倍磁性目標(biāo)自身長度時(shí),磁性目標(biāo)就可以看作一個(gè)磁偶極子[6-7]。距磁偶極子r處的磁場大小可以表示成

(1)

式中:μ0是真空磁導(dǎo)率;M是磁性目標(biāo)的磁矩;r是測(cè)量點(diǎn)到磁性目標(biāo)的距離;r0是沿r方向的單位方向向量。磁梯度張量表示磁場三分量在3個(gè)方向的空間變化率,由9個(gè)元素組成,其表達(dá)式為

(2)

式中:G為磁梯度張量。磁性目標(biāo)在其周圍產(chǎn)生的磁場我們將它看作靜磁場,由麥克斯韋方程組可知,B的旋度和散度都為0,所以G是對(duì)稱矩陣,只需要計(jì)算G的5個(gè)元素。定義CT=|G|,在磁偶極子模型場源條件下其表達(dá)式為

(3)

式中:參數(shù)k是在空間中為非球面分布的常數(shù),在定位矢量正交于磁矩的方向上取得最小值,在定位矢量平行于磁矩的方向上取得最大值。CT的等值線是以磁源為中心的近似橢圓[8]。

地球磁場約為5萬nT,分量有幾千到幾萬nT,但地磁場的梯度[9-10]小于等于0.02 nT/m,所以測(cè)量磁梯度張量可以消除地磁場的干擾。

磁梯度張量系統(tǒng)的測(cè)量陣列采用正六面體結(jié)構(gòu),如圖1所示,坐標(biāo)原點(diǎn)位于正六面體的中心,三個(gè)坐標(biāo)軸分別垂直穿過對(duì)應(yīng)正六面體的三個(gè)平面的中心點(diǎn),8個(gè)磁傳感器分別位于正六面體的8個(gè)頂點(diǎn),傳感器的三軸與坐標(biāo)系的三軸方向保持一致。正六面體每個(gè)平面內(nèi)都有4個(gè)傳感器,可以測(cè)得六個(gè)面中心的磁梯度張量,以Z軸正方向的平面為例:

圖1 張量系統(tǒng)測(cè)量陣列

(4)

(5)

式(4)中G為磁梯度張量,式(5)中S為每個(gè)磁傳感器測(cè)得的磁場3分量,d為正六面體測(cè)量陣列的邊長。

由空間梯度定義可知,標(biāo)量函數(shù)某一點(diǎn)的空間梯度指向函數(shù)在該點(diǎn)處空間變化率最大的方向。對(duì)測(cè)量陣列而言,其測(cè)量基點(diǎn)磁梯度張量的?？臻g梯度為

(6)

式中:i、j、k為正交坐標(biāo)系的基向量,ΔCT,X為X方向的兩個(gè)測(cè)量平面磁梯度張量的模的差值,即 ΔCT,X=CT,+X-CT,-X,其余兩個(gè)方向同理.CT的方向由測(cè)量陣列基點(diǎn)指向磁性目標(biāo),將其單位化即可得到單位定位矢量:

(7)

(8)

同理

(9)

(10)

聯(lián)立式(8)和式(10)得

(11)

求出單位定位矢量和定位矢量的模,最終得到定位矢量為

(12)

2 橢圓誤差消除方法

把Wiegert提出的STAR定位方法[11-12]稱為原始方法,原始方法由于存在橢圓系數(shù)而有一定的誤差。誤差來源于式(3)中的橢圓系數(shù)k,k值的大小因?yàn)榇啪嗯c距離矢量方向的夾角的變化而變化,而原始方法把k看作一個(gè)不變的常數(shù),所以造成了誤差。由于橢圓系數(shù)k的變化造成的計(jì)算誤差主要有三處:①式(6)對(duì)CT中的變量求導(dǎo),其中有兩個(gè)變量r和k,原始方法只對(duì)r進(jìn)行求導(dǎo),造成單位距離矢量的誤差。②式(10)中正六面體的X軸正負(fù)方向平面的CT相除,原始方法把k看作一個(gè)常數(shù),相互抵消,實(shí)際上兩個(gè)不同點(diǎn)的k有微小差別,造成距離矢量模的誤差。③式(11)求解定位矢量模R過程中,使用了式(6)計(jì)算得到的有誤差的單位距離矢量,也會(huì)造成距離矢量模的誤差.

通過上述分析可知,因?yàn)闄E圓系數(shù)的存在,原始方法計(jì)算得到的單位距離矢量及距離矢量的模都存在誤差。

2.1 新的單位距離矢量計(jì)算方法

本文提出了一種新的單位定位矢量計(jì)算方法,所提方法利用測(cè)量點(diǎn)和磁性目標(biāo)形成的位置矢量與磁梯度張量矩陣絕對(duì)值最小的特征值對(duì)應(yīng)的特征向量之間的垂直關(guān)系[13],這個(gè)垂直關(guān)系由磁梯度張量不變量推導(dǎo)得出。

該方法使用如圖1所示的正六面體結(jié)構(gòu)的測(cè)量陣列,坐標(biāo)原點(diǎn)位于正六面體的中心,三個(gè)坐標(biāo)軸分別垂直穿過對(duì)應(yīng)正六面體的兩個(gè)平面的中心點(diǎn),測(cè)量陣列每個(gè)面都有4個(gè)三軸磁傳感器,可以測(cè)出每個(gè)正六面體每個(gè)面中心處的磁梯度張量。

設(shè)X軸正方向平面中心的點(diǎn)是P1,P1點(diǎn)的磁梯度張量為G1,矩陣G1的特征值為λ1、λ2、λ3,且|λ1|>|λ3|,|λ2|>|λ3|,λ3對(duì)應(yīng)的特征向量是V3,Y軸正方向平面中心的點(diǎn)是P2,P2點(diǎn)的磁梯度張量為G2。矩陣G2的特征值為λ4、λ5、λ6,且|λ4|>|λ6|,|λ5|>|λ6|,λ6對(duì)應(yīng)的特征向量是V6,Z軸正方向平面中心的點(diǎn)是P3,P3點(diǎn)的磁梯度張量為G3。矩陣G3的特征值為λ7、λ8、λ9,且|λ7|>|λ9|,|λ8|>|λ9|,λ9對(duì)應(yīng)的特征向量是V9。V3、V6和V9是正六面體3個(gè)面磁梯度張量矩陣絕對(duì)值最小的特征值對(duì)應(yīng)的特征向量。

正六面體的中心O到磁性目標(biāo)之間的位置矢量為r,P1點(diǎn)到磁性目標(biāo)之間的位置矢量為r1,P2點(diǎn)到磁性目標(biāo)之間的位置矢量為r2,P3點(diǎn)到磁性目標(biāo)之間的位置矢量為r3,由磁梯度張量不變量推導(dǎo)出的垂直[13]關(guān)系可得

(13)

由O、P1、P2的位置關(guān)系,利用向量運(yùn)算可以得到:

(14)

將式(14)中的等式代入式(13)可以得到

(15)

f=max(|f1|,|f2|,|f3|)

(16)

利用粒子群優(yōu)化算法對(duì)式(16)中x,y和z進(jìn)行優(yōu)化,使得f最小,即可求解出x,y,z。

粒子群優(yōu)化算法是一種啟發(fā)式優(yōu)化算法,它是從隨機(jī)解出發(fā),通過迭代來尋找最優(yōu)解,通過適應(yīng)度函數(shù)來判斷解是否最優(yōu),通過追隨當(dāng)前最優(yōu)解來尋找全局最優(yōu)解。這種算法的優(yōu)點(diǎn)是實(shí)現(xiàn)容易、精度高、收斂快,在解決實(shí)際問題的過程有巨大的優(yōu)越性。所以單位距離矢量為

(17)

2.2 距離矢量模的修正方法

正六面體的X軸正負(fù)方向平面中心的橢圓系數(shù)k有微小的差別,所以兩個(gè)平面中心點(diǎn)的橢圓系數(shù)分別設(shè)為k1,k2,聯(lián)立式(3)、式(10)和式(17)可得

(18)

聯(lián)立式(8)和式(17)得

(19)

由文獻(xiàn)[14]可知橢圓系數(shù)的公式

(20)

由文獻(xiàn)[13]可知磁梯度張量的特征值不變量公式

(21)

聯(lián)立式(20)和式(21)得

(22)

所以有

(23)

將式(23)代入式(19)可以得到修正的距離矢量的模

(24)

求出單位距離矢量和距離矢量的模,得到修正后的距離矢量為

(25)

3 仿真實(shí)驗(yàn)分析

針對(duì)本文提出的橢圓誤差修正方法,對(duì)可能造成誤差的因素進(jìn)行仿真,并分析仿真結(jié)果。仿真實(shí)驗(yàn)工具有MATLAB軟件和COMSOL多物理場仿真軟件?？赡苡绊懚ㄎ痪鹊囊蛩赜?①測(cè)量系統(tǒng)的基線d;②高斯噪聲;③磁傳感器的測(cè)量精度。

仿真時(shí)通過式(12)計(jì)算原始方法的距離矢量,通過式(25)計(jì)算改進(jìn)方法的距離矢量。

仿真實(shí)驗(yàn)的初始條件:測(cè)量陣列是圖1所示正六面體,以正六面體的中心為坐標(biāo)原點(diǎn)建立笛卡爾坐標(biāo)系,設(shè)磁源坐標(biāo)為(10,20,30)單位m;磁性目標(biāo)的磁矩m=(5 000 000,5 000 000,7 071 100)Am2;當(dāng)?shù)氐牡卮艌稣鎸?shí)值為(30 000,2 200,44 000)nT。

3.1 正六面體測(cè)量系統(tǒng)的基線

仿真方法:仿真時(shí)正六面體測(cè)量系統(tǒng)的基線d從0.1 m增加到2.2 m,其他量保持初始條件不變,加入地磁場.比較理論值與仿真結(jié)果,得出結(jié)論。

仿真結(jié)果如圖2,圖2橫坐標(biāo)表示正六面體基線d從0.1 m增加到2.2 m,圖2縱坐標(biāo)為測(cè)量點(diǎn)與理論位置的絕對(duì)誤差,單位是m。三種線型的曲線表示三種傳感器精度對(duì)應(yīng)的仿真結(jié)果。

圖2 測(cè)量系統(tǒng)基線d變化時(shí)的定位誤差

從圖2可知,正六面體基線d小于等于0.4 m時(shí),定位誤差較小,比較穩(wěn)定,當(dāng)d大于0.4 m時(shí),系統(tǒng)的定位誤差逐漸增大。

3.2 高斯噪聲仿真實(shí)驗(yàn)

仿真方法:正六面體的基線d=0.4,測(cè)量系統(tǒng)仿真時(shí),加入高斯噪聲,噪聲的信噪比從10增加到60,其他量保持初始條件不變,加入地磁場,每個(gè)信噪比都仿真100次,然后求解100次仿真結(jié)果的平均值。比較理論值與仿真結(jié)果,得出結(jié)論。

仿真結(jié)果如圖3,圖3橫坐標(biāo)表示高斯噪聲的信噪比從10增加到60,圖3縱坐標(biāo)為100次仿真結(jié)果的平均值與理論位置的絕對(duì)誤差,單位是m。三種線型的曲線表示三種傳感器精度對(duì)應(yīng)的仿真結(jié)果。

圖3 高斯噪聲信噪比變化時(shí)的定位誤差

從圖3可知,隨著高斯噪聲信噪比的增大,測(cè)量系統(tǒng)的定位誤差迅速減小,當(dāng)信噪比大于等于50時(shí),定位誤差達(dá)到最小并趨于穩(wěn)定。

3.3 磁傳感器的精度

仿真方法:正六面體的基線d=0.4,磁源沿一條直線運(yùn)動(dòng),在不同傳感器精度條件下,其他量保持初始條件不變,分別測(cè)量出對(duì)應(yīng)磁源的坐標(biāo)。具體方法是,磁源從坐標(biāo)(10,20,0)沿Z軸正方向移動(dòng)到(10,20,70),單位是m,每次移動(dòng)2.5 m,選擇1 nT、0.1 nT和0.01 nT三種傳感器精度,在其他量保持初始條件不變的前提下,加入地磁場,分別測(cè)量這三種精度下每次磁源移動(dòng)后的坐標(biāo)。

仿真結(jié)果如圖4,圖4橫坐標(biāo)表示磁源的Z坐標(biāo)從0增加到70,圖4縱坐標(biāo)為三種傳感器精度測(cè)量的絕對(duì)誤差,單位是m。三種線型的曲線表示三種傳感器精度對(duì)應(yīng)的仿真結(jié)果。

圖4 傳感器精度變化時(shí)的定位誤差

在圖2到圖4中三種線型的曲線表示1 nT、0.1 nT和0.01 nT三種傳感器精度對(duì)應(yīng)的仿真結(jié)果,可以看出不論是在系統(tǒng)基線的仿真實(shí)驗(yàn),還是高斯噪聲實(shí)驗(yàn)中,還是在磁源直線運(yùn)動(dòng)的傳感器精度實(shí)驗(yàn)中,都有以下結(jié)論:傳感器精度越高,定位誤差越小。

由圖4還可以看出當(dāng)傳感器精度為1 nT,磁源在點(diǎn)(10,20,20)時(shí),測(cè)量誤差開始增大;當(dāng)傳感器精度為0.1 nT,磁源在點(diǎn)(10,20,30)時(shí),測(cè)量誤差開始增大;當(dāng)傳感器精度為0.01 nT,磁源在點(diǎn)(10,20,55)時(shí),測(cè)量誤差開始增大;由此可以得出以下結(jié)論:本文所提改進(jìn)定位方法,傳感器精度為1 nT時(shí)的測(cè)量半徑是30 m,傳感器精度為0.1 nT時(shí)的測(cè)量半徑是37 m,傳感器精度為0.01 nT時(shí)的測(cè)量半徑是60 m左右。

3.4 磁源的圓形軌跡運(yùn)動(dòng)驗(yàn)證實(shí)驗(yàn)

為了測(cè)試系統(tǒng)的整體性能,本文進(jìn)行了磁源圓形軌跡運(yùn)動(dòng)定位實(shí)驗(yàn)。在加入地磁場基礎(chǔ)上,正六面體的基線d=0.4,噪聲的信噪比為50,傳感器精度為0.1 nT,磁源沿圓的軌跡移動(dòng),使用原始方法和本文所提改進(jìn)方法分別對(duì)磁源進(jìn)行定位,觀察定位結(jié)果,分析兩種定位方法的軌跡,得出結(jié)論。

正六面體測(cè)量陣列以正六面體的中心為坐標(biāo)原點(diǎn)建立笛卡爾坐標(biāo)系。在地磁場和高斯噪聲環(huán)境下,其他量保持初始條件不變,磁源沿圓的軌跡從初始位置(10,20,30)繞Z軸逆時(shí)針旋轉(zhuǎn)一周,每次旋轉(zhuǎn)10°,共旋轉(zhuǎn)36次,每次磁源位置更新后,分別使用原始方法和本文所提改進(jìn)方法分別對(duì)磁源位置進(jìn)行定位,記錄定位結(jié)果。

仿真結(jié)果如圖5、圖6和圖7,三個(gè)圖中的實(shí)線是圓形軌跡的理論值。

圖5 磁源圓形軌跡運(yùn)動(dòng)原始方法橫縱坐標(biāo)仿真圖

圖6 磁源圓形軌跡運(yùn)動(dòng)改進(jìn)方法橫縱坐標(biāo)仿真圖

圖7 磁源圓形軌跡運(yùn)動(dòng)兩種方法z坐標(biāo)坐標(biāo)仿真圖

圖5中的虛線是使用原始方法測(cè)量出的磁源軌跡橫縱坐標(biāo),圖6中的虛線是使用本文所提改進(jìn)方法測(cè)量出的磁源軌跡橫縱坐標(biāo),圖5和圖6的橫坐標(biāo)都是磁源軌跡的橫坐標(biāo),單位是m,圖5和圖6的縱坐標(biāo)都是磁源軌跡的縱坐標(biāo),單位是m,圖7中的點(diǎn)線是使用原始方法測(cè)量出的磁源軌跡z坐標(biāo),圖7中的虛線是使用本文所提改進(jìn)方法測(cè)量出的磁源軌跡z坐標(biāo),圖7的橫坐標(biāo)表示磁源的第n次旋轉(zhuǎn),圖7的縱坐標(biāo)表示兩種方法測(cè)量出的磁源軌跡z坐標(biāo),單位是m。

從圖5～圖7可知,在磁源的圓形軌跡運(yùn)動(dòng)定位實(shí)驗(yàn)中,本文所提改進(jìn)定位方法在x,y和z軸的誤差都在1 m以內(nèi),相對(duì)于原始方法的誤差有很大程度的減小,定位精度大幅度提高。

4 總結(jié)

本文針對(duì)現(xiàn)有的磁梯度張量不變量定位算法存在橢圓誤差的問題,提出了一種利用磁梯度張量不變量消除橢圓誤差的改進(jìn)定位方法。根據(jù)仿真實(shí)驗(yàn)的結(jié)果可得,本文所提的改進(jìn)方法消除了橢圓誤差,定位誤差相比于原始方法有大幅度的減小,可以實(shí)現(xiàn)地磁環(huán)境和一定信噪比的高斯噪聲情況下的實(shí)時(shí)精確定位,而且該方法適用于移動(dòng)平臺(tái)。在改進(jìn)方法中,測(cè)量系統(tǒng)的基線d=0.4 m時(shí),可以取得較小的定位誤差;信噪比大于50的高斯噪聲對(duì)測(cè)量系統(tǒng)的影響較小;傳感器的精度越高,定位誤差越小。所以本文所提的改進(jìn)方法有很好的定位效果,不足的地方就是對(duì)傳感器的精度要求較高。