賈文抖， 林春生， 孫玉繪， 翟國君

(1.海軍工程大學(xué) 兵器工程學(xué)院， 湖北 武漢 430033； 2.海軍工程大學(xué) 導(dǎo)航工程系， 湖北 武漢 430033；3.海軍海洋測繪研究所， 天津 300061)

0 引言

磁梯度張量具有能夠有效克服地磁場的干擾、突出目標(biāo)磁信號、更加全面地反映磁場細(xì)節(jié)等突出的優(yōu)勢，因此近年來基于磁梯度張量的探測技術(shù)逐步成為磁探測領(lǐng)域的研究熱點(diǎn)，德國、美國、澳大利亞等國已相繼研制出磁梯度張量探測系統(tǒng)，進(jìn)行了大量的野外試驗[1]，并逐步應(yīng)用于軍事和民用領(lǐng)域。

基于歐拉齊次方程的磁梯度張量定位方法是早期人們研究的一個重點(diǎn)，但該方法需要獲取磁目標(biāo)磁場中測量點(diǎn)上干凈的目標(biāo)三分量磁信號，這在地磁場環(huán)境中是難以實現(xiàn)的，實際應(yīng)用效果也不夠理想。磁梯度張量不變量是利用磁梯度張量數(shù)據(jù)計算出來的一些不隨測量坐標(biāo)系變化而變化的量，它反映了磁目標(biāo)的磁場特征，不易受地磁場的干擾，因此適用于地磁場中磁目標(biāo)的探測。另外，由于磁梯度張量不變量還具有不隨測量坐標(biāo)系的變化而改變的性質(zhì)，適合于搭載在移動探測平臺上。常見的不變量有磁梯度張量的特征值、跡、等效衰減磁矩μ等[2-3]。

江勝華等[4]依據(jù)磁梯度張量縮并理論對磁梯度張量的模Gt及參數(shù)k的空間分布特征進(jìn)行了分析，結(jié)果顯示Gt和參數(shù)k呈橢球形分布，并擬合得到k值分布的近似計算公式。陳謹(jǐn)飛等[5]建立了正六面體測量模型，利用磁梯度張量的模Gt進(jìn)行磁目標(biāo)定位，但由于參數(shù)k的存在使得定位結(jié)果有一定的偏差。若能對參數(shù)k的影響進(jìn)行補(bǔ)償，或者將其從磁梯度張量的模中完全消除掉，則可以進(jìn)一步提高磁定位的準(zhǔn)確性。Nara等[6]研究發(fā)現(xiàn)影響參數(shù)k分布變化的因素是磁矩與位置矢量之間的夾角，并給出了參數(shù)k隨夾角變化的計算公式。呂俊偉等[2]研究發(fā)現(xiàn)參數(shù)k可以用磁梯度張量特征值的組合公式進(jìn)行計算，從磁梯度張量的模Gt中消除參數(shù)k后，得到了一個僅與距離有關(guān)的不變量。該不變量就是等效衰減磁矩μ，而磁梯度張量的模Gt與等效衰減磁矩μ之間的比值就是參數(shù)k. 此外，Sui等[7]也對如何消除參數(shù)k帶來的橢圓誤差進(jìn)行了研究，并進(jìn)行了相應(yīng)的實驗，修正效果較好。實際上，關(guān)于不變量——等效衰減磁矩μ的研究在國外文獻(xiàn)中早有出現(xiàn)，在1985年Wilson[8]就提出過與磁目標(biāo)磁矩方向無關(guān)的不變量。Beiki等[9]研究發(fā)現(xiàn)當(dāng)直線上的測點(diǎn)距離磁目標(biāo)最近時，不變量等效衰減磁矩μ取得極大值這一特性。Clark[10]計算了磁梯度張量的3個特征值，發(fā)現(xiàn)可以用磁梯度張量的特征值直接計算出等效衰減磁矩μ，繼而可以根據(jù)μ的分布規(guī)律確定磁目標(biāo)的位置參數(shù)。Beiki等[11]進(jìn)一步研究發(fā)現(xiàn)等效衰減磁矩μ同樣滿足歐拉齊次方程，并可以利用多點(diǎn)處μ的梯度值反演得到磁目標(biāo)的位置。

在計算文獻(xiàn)[2，4]中正六面體模型的六面體側(cè)面磁梯度張量數(shù)據(jù)時，首先利用口字型布置的磁力儀測得的磁場數(shù)據(jù)近似計算磁梯度張量中的9個量，然后計算出相應(yīng)的磁梯度張量不變量；計算目標(biāo)源與測點(diǎn)的距離r時，根據(jù)上下2個側(cè)面的磁梯度張量不變量與距離之間的比例關(guān)系，借助麥克勞林公式簡化近似計算得到距離r的計算公式，本文稱這種方法為方法1. 為了進(jìn)一步降低上述計算誤差，得到更精確的磁梯度張量信息和距離r，本文對方法1進(jìn)行改進(jìn)，在計算正六面?zhèn)让嬷行狞c(diǎn)的磁梯度張量中的分量時，采用坐標(biāo)系旋轉(zhuǎn)的方法將口字型布置的磁力儀陣列變成十字型測量陣列，相比方法1，所用到的磁力儀數(shù)目由4減少到2. 由于正六面體中心點(diǎn)既不位于正六面體中某個正方形中心點(diǎn)也不位于十字型交叉點(diǎn)，該中心點(diǎn)的磁梯度張量既不能利用口字型磁力儀陣列也不能利用十字型磁力儀陣列差分計算得到。為此根據(jù)全微分概念，提出一種微分求解方法，利用正六面體上下對角磁力儀測量的磁梯度數(shù)據(jù)建立微分方程組，求解計算磁梯度張量中的9個元素量，解決了正六面體中心點(diǎn)磁梯度張量難以計算的問題。在計算目標(biāo)源與測點(diǎn)的距離r時，直接利用等效衰減磁矩μ的梯度公式推導(dǎo)得到r，消除了方法1中借助麥克勞林公式簡化計算帶來的誤差。最后根據(jù)μ的圓周分布特征進(jìn)行磁目標(biāo)位置參數(shù)的反演，實現(xiàn)磁目標(biāo)的定位。

1 基于等效衰減磁矩的定位方法

1.1 磁梯度張量特征值的旋轉(zhuǎn)不變特性

將磁偶極子視為研究對象，以磁偶極子位置為原點(diǎn)(0, 0, 0)建立直角坐標(biāo)系Oxyz，磁偶極子的磁矩M=(mx,my,mz)，位置矢量為r=(x,y,z)的測量點(diǎn)Q處磁感應(yīng)強(qiáng)度為

(1)

B的三分量形式為

(2)

磁梯度張量G為磁場三分量在空間x、y、z方向上的導(dǎo)數(shù)，其形式為

(3)

式中：

(4)

下標(biāo)i、j表示x、y、z方向中的任意2個，δij是克羅內(nèi)克函數(shù)，當(dāng)i=j時δij=1，當(dāng)i≠j時δij=0，ri、rj表示距離矢量r=(x,y,z)中的任意兩個分量。

根據(jù)特征方程det (G-λI)=0(I是與G維數(shù)相同的單位對角陣，λ是待求的磁梯度張量特征值)，計算得到G的特征值為

(5)

式中：

假設(shè)坐標(biāo)系Oxyz繞x軸旋轉(zhuǎn)角后得到坐標(biāo)系Oxαyαzα. 在坐標(biāo)系Oxαyαzα中，測量點(diǎn)Q的坐標(biāo)、磁矩M可分別表示為

(6)

(7)

將(6)式、(7)式代入(5)式，計算得到坐標(biāo)系Oxαyαzα中測量點(diǎn)Q處磁梯度張量對應(yīng)的特征值為

(8)

由(8)式可知，磁偶極子上坐標(biāo)系Oxyz繞x軸旋轉(zhuǎn)前后磁梯度張量G的特征值不變。當(dāng)坐標(biāo)系繞y軸、z軸旋轉(zhuǎn)時，同樣可得磁梯度張量特征值不隨測量坐標(biāo)系的旋轉(zhuǎn)而變化。由于位于磁偶極子處的任意直角坐標(biāo)系都可以由坐標(biāo)系Oxyz繞經(jīng)x軸、y軸、z軸旋轉(zhuǎn)得到，根據(jù)以上分析可知，當(dāng)固定了測量點(diǎn)與磁目標(biāo)點(diǎn)的相對位置，同一測量點(diǎn)的磁梯度張量特征值在不同直角坐標(biāo)系中都是相等的。因此，磁梯度張量的特征值是磁梯度張量的旋轉(zhuǎn)不變量。

1.2 等效衰減磁矩的定位原理

(9)

式中：

(10)

根據(jù)特征方程

det (G-λI)=-(μcosβ-λ)[(2μcosβ+
λ)(μcosβ-λ)+μ2sin2β]=0,

求得特征值為

(11)

由(11)式可得

(12)

當(dāng)矢量r與磁矩M同向(β=0)或反向(β=180)時G退化為對角矩陣，對角線上的元素即為特征值，磁梯度張量的特征值與等效衰減磁矩之間的關(guān)系同樣滿足(12)式。由于特征值是磁梯度張量的不變量，因此μ也是磁梯度張量的不變量，且與位置矢量r的模的4次方呈反比。

(13)

(14)

根據(jù)(13)式可得，測量點(diǎn)Q與磁偶極子之間的距離為

(15)

則位置矢量r為

r=-rμe.

(16)

2 正六面體定位模型

(17)

以計算x軸正向?qū)?yīng)面中心的μx-為例，磁力儀S1、S2、S3和S4測量的磁感應(yīng)強(qiáng)度分別為

(18)

由第1節(jié)的分析可知，μ不隨計算坐標(biāo)系的變化而變化。將磁力儀S1、S2、S3和S4測得的3個分量數(shù)據(jù)對應(yīng)到坐標(biāo)系Ox45°y45°z45°中(見圖3)，坐標(biāo)系Ox45°y45°z45°由坐標(biāo)系Oxyz繞x軸順時針旋轉(zhuǎn)45°得到。

測量數(shù)據(jù)由坐標(biāo)系Oxyz對應(yīng)到坐標(biāo)系Ox45°y45°z45°中的轉(zhuǎn)換公式為

(19)

式中：Bx、By、Bz為坐標(biāo)系Oxyz中磁力儀直接測得磁場的3個分量；Bx45°、By45°、Bz45°為對應(yīng)到坐標(biāo)系Ox45°y45°z45°中的3個分量。

在坐標(biāo)系Ox45°y45°z45°中，直接利用差分方法計算磁梯度張量：

G45°=

計算正六面體中心點(diǎn)P處的等效衰減磁矩μP時，由于正六面體中心點(diǎn)處的磁梯度張量無法按照差分方法計算得到，在此利用全微分的概念，建立微分方程組間接求解磁梯度張量各分量[12]。將圖2中的8個磁力儀按交叉線分成4組：S1和S8、S2和S7、S3和S6、S4和S5. 這4組磁力儀各自的連線均過正六面體中心點(diǎn)P. 以磁力儀S1和S8這一組為例，磁力儀S1和S8的測量數(shù)據(jù)差值為

ΔB81=(ΔB81x,ΔB81y,ΔB81z)=

(B8x-B1x,B8y-B1y,B8z-B1z).

(20)

根據(jù)全微分的概念，結(jié)合S1、S8的空間位置，有

(21)

結(jié)合其他3組磁力儀數(shù)據(jù)，則有方程組：

(22)

求方程組(22)式即可得到Gxx、Gxy和Gxz. 類似地，利用ΔBy和ΔBz的全微分表達(dá)式，可計算得到Gyx、Gyy、Gyz、Gzx、Gzy和Gzz. 得到正六面體中心點(diǎn)P處的磁梯度張量G后，可計算出等效衰減磁矩μP.

3 仿真分析

仿真條件設(shè)置為磁偶極子磁矩M=(3 000 A·m2，3 300 A·m2，4 200 A·m2)，磁目標(biāo)的真實位置坐標(biāo)r=(-4 m，-6 m，-11 m)，計算得到的位置參數(shù)用rP表示，位置參數(shù)的相對誤差為

下面分析地磁場、磁梯度測量噪聲、正六面體基線以及磁力儀測量精度等因素對位置參數(shù)反演結(jié)果的影響，并與方法1的反演結(jié)果進(jìn)行比較。

3.1 地磁場對定位結(jié)果的影響

假設(shè)研究環(huán)境中地磁場大小為50 000 nT，地磁場梯度0.02 nT/m[13]，地磁傾角I=42°，磁偏角T=3°(設(shè)北偏東為正)。正六面體模型中的坐標(biāo)系x軸正向朝東，y軸正向朝北，z軸正向朝上。測量數(shù)據(jù)分兩種：一組數(shù)據(jù)中包含地磁場；一組數(shù)據(jù)中不含地磁場。根據(jù)設(shè)定的磁目標(biāo)參數(shù)，利用仿真的兩組數(shù)據(jù)計算磁目標(biāo)的位置參數(shù)，計算結(jié)果見表1.

表1 有無地磁場時的定位結(jié)果

從表1中的計算結(jié)果可以看出，以磁偶極子作為磁目標(biāo)，利用本文方法進(jìn)行磁目標(biāo)位置參數(shù)的反演計算時，地磁場對定位結(jié)果的影響很微弱，在應(yīng)用中地磁場的影響完全可以忽略。

3.2 磁梯度噪聲對定位結(jié)果的影響

在測量過程中，由于周圍的磁場環(huán)境比較復(fù)雜，各種設(shè)備及人為干擾難以避免，且各磁力儀之間還存在匹配誤差，使得測量的目標(biāo)磁場梯度信號中不可避免地包含一定的磁干擾量，測得的磁梯度信號不純凈。本文將這些磁干擾等效為磁梯度噪聲，添加到目標(biāo)磁信號測量數(shù)據(jù)中，磁梯度測量噪聲對定位結(jié)果的影響見圖4.

圖4中藍(lán)線是磁梯度測量噪聲水平為0.353 6 nT/m、隨機(jī)進(jìn)行100次計算得到的定位誤差曲線，黑線是磁梯度測量噪聲水平為0.707 2 nT/m、隨機(jī)進(jìn)行100次計算得到的定位誤差曲線。由圖4可見：磁梯度噪聲水平為0.353 6 nT/m時，定位誤差不超過1 m，大部分誤差在0.5 m以內(nèi)；磁梯度噪聲為0.707 2 nT/m時，定位誤差幾乎全在1.5 m以內(nèi)，大部分不超過1 m. 由于初始設(shè)置磁目標(biāo)的磁矩不大，在測量點(diǎn)上目標(biāo)信號的磁梯度大小約為幾十納特斯拉每米，相當(dāng)于測量信號中的磁梯度噪聲不超過目標(biāo)信號強(qiáng)度的10%，定位相對誤差也基本維持在10%以內(nèi)。

3.3 基線長度對定位結(jié)果的影響

根據(jù)磁梯度張量G45°的表達(dá)式可知，采用差分替代微分的方法計算磁梯度張量時，基線長度會影響G45°的計算結(jié)果，并最終影響位置參數(shù)的反演結(jié)果。定位誤差隨基線長度變化的曲線見圖5.

從圖5可以看出，定位誤差隨正六面體邊長基線的增加而變大。主要是因為計算六面體側(cè)面的磁梯度張量時采用的是差分計算方法，與磁梯度張量的理論值之間存在一定的偏差，差分基線越長，偏差越大，定位誤差也越大。在本節(jié)設(shè)置的條件下，當(dāng)邊長基線不超過1.5 m時，定位誤差在0.7 m以內(nèi)。

3.4 磁力儀測量精度對定位結(jié)果的影響

當(dāng)磁力儀測量精度不同時，從表2中的定位結(jié)果可以看出，測量精度從0.01～1 nT逐漸降低時定位誤差逐漸加大，當(dāng)磁力儀的測量精度為1 nT時，定位誤差在1.4 m左右，相對誤差約為11%.

表2 不同測量精度時的定位誤差

3.5 本文方法與方法1的比較

表3 本文方法與方法1的計算結(jié)果比較

由表3可知，在磁力儀測量精度相同情況下，利用本文方法反演計算得到的位置參數(shù)與真實的位置參數(shù)更加接近，定位誤差為0.101 m，相對誤差為0.77%. 利用方法1反演計算磁目標(biāo)位置參數(shù)時，定位誤差為0.418 m，相對誤差為3.18%. 由此可見，本文方法的定位精度高于方法1. 另外，在磁梯度計算基線相等的情況下，本文正六面體模型的空間體積僅為方法1中模型體積的35.3%，顯著降低了模型所占的空間體積。

4 結(jié)論

本文研究了在地磁場中利用磁目標(biāo)的等效衰減磁矩信息實現(xiàn)磁目標(biāo)定位的方法。通過計算分析了等效衰減磁矩μ不隨測量坐標(biāo)系的改變而變化的特征。根據(jù)這一特征，利用了坐標(biāo)系旋轉(zhuǎn)的方法將口字型布置的磁力儀陣列轉(zhuǎn)變成十字型，在差分基線相等的情況下使正六面體的邊長縮減為原來的0.707倍。利用求解全微分方程組的方法計算得到正六面體中心點(diǎn)的磁梯度張量，解決了正六面體中心點(diǎn)磁梯度張量難以計算的問題。對磁目標(biāo)位置參數(shù)進(jìn)行了反演，結(jié)果表明在地磁場中利用等效衰減磁矩的定位方法可實現(xiàn)對磁目標(biāo)的有效定位。與之前的計算方法相比，在計算磁梯度的差分基線相等的情況下，本文方法正六面體模型的體積減小了64.7%，且定位效果更好。

)

[1] 于振濤，呂俊偉，畢波，等．四面體磁梯度張量系統(tǒng)的載體磁干擾補(bǔ)償方法[J]．物理學(xué)報，2014，63(11)：110702-1-110702-6．

YU Zhen-tao, LYU Jun-wei, BI Bo, et al. A vehicle magnetic noise compensation method for the tetrahedron magnetic gradiometer[J]. Acta Physica Sinica，2014, 63(11): 110702-1-110702-6. (in Chinese)

[2] 呂俊偉，遲鋮，于振濤，等．磁梯度張量不變量的橢圓誤差消除方法研究[J]．物理學(xué)報，2015，64(19)：190701-1-190701-8．

LYU Jun-wei, CHI Cheng, YU Zhen-tao, et al. Research on the asphericity error elimination of the invariant of magnetic gradient tensor[J]. Acta Physica Sinica，2015, 64(19): 190701-1-190701-8. (in Chinese)

[3] 萬成彪，潘孟春，張琦，等. 基于張量特征值和特征向量的磁性目標(biāo)定位[J]. 吉林大學(xué)學(xué)報：工學(xué)版，2017，47(2)：655-660．

WAN Cheng-biao, PAN Meng-chun, ZHANG Qi, et al. Magnetic object location with eigenvalue and eigenvector of tensor[J]. Journal of Jilin University: Engineering and Technology Edition, 2017, 47(2):655-660. (in Chinese)

[4] 江勝華，申宇，褚玉程．基于磁偶極子的磁場梯度張量縮并的試驗驗證及相關(guān)參數(shù)確定[J]．中國慣性技術(shù)學(xué)報，2015，23(1)：103-106．

JIANG Sheng-hua, SHEN Yu, CHU Yu-cheng. Experimental verification and related parameter's determination for magnetic gradient tensor contraction using magnetic dipole[J]. Journal of Chinese Inertial Technology, 2015, 23(1): 103-106. (in Chinese)

[5] 陳謹(jǐn)飛，張琦，潘孟春，等．基于正六面體結(jié)構(gòu)測量陣列的磁異常定位技術(shù)研究[J]．傳感技術(shù)學(xué)報，2012，25(8)：1088-1092．

CHEN Jin-fei, ZHANG Qi, PAN Meng-chun, et al. Research on geomagnetic anomaly localization based on cubic measurement array[J]. Chinese Journal of Sensors and Actuators, 2012, 25(8): 1088-1092. (in Chinese)

[6] Nara T, Ito W. Moore-penrose generalized inverse of the gradient tensor in Euler's equation for locating a magnetic dipole[J]. Journal of Applied Physics, 2014, 115(17): 17E504-1-17E504-3.

[7] Sui Y Y, Li G, Wang S L, et al. Asphericity errors correction of magnetic gradient tensor invariants method for magnetic dipole localization[J]. IEEE Transactions on Magnetics, 2012, 48(12):4701-4706.

[8] Wilson H. Analysis of the magnetic gradient tensor[R]. Canada: Defence Research Establishment Pacific, 1985: 5-13.

[9] Beiki M, Keating P, Clark D. Depth estimation of magnetic and gravity sources using normalized source strength calculated from gradient tensor[C]∥2012 SEG Annual Meeting.Las Vegas, CA，US: The Society of Exploration Geophysicists, 2012: 996.

[10] Clark D A. New methods for interpretation of magnetic vector and gradient tensor data I: eigenvector analysis and the normalised source strength[J]. Exploration Geophysics, 2012, 43(4):267-282.

[11] Beiki M, Clark D A, Austin J R, et al. Estimating source location using normalized magnetic source strength calculated from magnetic gradient tensor data[J]. Geophysics, 2012, 77(6):23-37.

[12] 賈文抖，林春生，孫玉繪，等．基于單個磁梯度計的磁目標(biāo)定位方法研究[J]．兵工學(xué)報，2017，38(8)：1572-1577．

JIA Wen-dou, LIN Chun-sheng, SUN Yu-hui, et al. Research on magnetic target location method based on a single magnetic gradiometer[J]. Acta Armamentarii, 2017, 38(8): 1572-1577. (in Chinese)

[13] 張光，張英堂，李志寧，等．載體平動條件下的磁梯度張量定位方法[J]．華中科技大學(xué)學(xué)報：自然科學(xué)版，2013，41(1)：21-24．

ZHANG Guang, ZHANG Ying-tang, LI Zhi-ning, et al. Localizing method of magnetic field gradient tensor under carriers moving parallelly[J]. Journal of Huazhong University of Science and Technology: Natural Science Edition, 2013, 41(1): 21-24. (in Chinese)