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

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

0 引言

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

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

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

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

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

1.1 磁梯度張量特征值的旋轉不變特性

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

(1)

B的三分量形式為

(2)

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

(3)

式中：

(4)

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

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

(5)

式中：

假設坐標系Oxyz繞x軸旋轉角后得到坐標系Oxαyαzα. 在坐標系Oxαyαzα中，測量點Q的坐標、磁矩M可分別表示為

(6)

(7)

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

(8)

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

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

(9)

式中：

(10)

根據(jù)特征方程

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

求得特征值為

(11)

由(11)式可得

(12)

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

(13)

(14)

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

(15)

則位置矢量r為

r=-rμe.

(16)

2 正六面體定位模型

(17)

以計算x軸正向對應面中心的μx-為例，磁力儀S1、S2、S3和S4測量的磁感應強度分別為

(18)

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

測量數(shù)據(jù)由坐標系Oxyz對應到坐標系Ox45°y45°z45°中的轉換公式為

(19)

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

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

G45°=

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

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

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

(20)

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

(21)

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

(22)

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

3 仿真分析

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

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

3.1 地磁場對定位結果的影響

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

表1 有無地磁場時的定位結果

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

3.2 磁梯度噪聲對定位結果的影響

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

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

3.3 基線長度對定位結果的影響

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

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

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

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

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

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

表3 本文方法與方法1的計算結果比較

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

4 結論

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

)

[1] 于振濤，呂俊偉，畢波，等．四面體磁梯度張量系統(tǒng)的載體磁干擾補償方法[J]．物理學報，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]．物理學報，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] 萬成彪，潘孟春，張琦，等. 基于張量特征值和特征向量的磁性目標定位[J]. 吉林大學學報：工學版，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] 江勝華，申宇，褚玉程．基于磁偶極子的磁場梯度張量縮并的試驗驗證及相關參數(shù)確定[J]．中國慣性技術學報，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]．傳感技術學報，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] 賈文抖，林春生，孫玉繪，等．基于單個磁梯度計的磁目標定位方法研究[J]．兵工學報，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]．華中科技大學學報：自然科學版，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)