李青竹 李志寧 張英堂 范紅波

(陸軍工程大學(xué)車輛與電氣工程系，河北石家莊 050003)

0 引言

磁梯度張量(Magnetic gradient tensor，MGT)探測技術(shù)及其數(shù)據(jù)解釋與應(yīng)用可視為磁法勘探領(lǐng)域的一項突破[1]。相較于磁總場強(qiáng)度和磁場矢量，磁梯度張量作為磁場矢量三正交方向上的空間變化率，對測量點的空間取向和旋轉(zhuǎn)噪聲不敏感，梯度測量與背景勻強(qiáng)場環(huán)境無關(guān)。磁梯度張量技術(shù)可提取關(guān)于目標(biāo)體更豐富的姿態(tài)信息和磁源信息，分辨率高且抗干擾能力強(qiáng)，對復(fù)雜環(huán)境的適應(yīng)性較強(qiáng)[2-3]。利用磁梯度張量數(shù)據(jù)可實現(xiàn)航磁探測、礦產(chǎn)勘探與土壤黑色金屬搜索、探尋未爆彈、排雷、潛艇偵查或水下金屬目標(biāo)定位及反演識別等[4-6]，應(yīng)用前景十分廣闊。

當(dāng)觀測距離超過目標(biāo)體尺度的2.5倍時，磁性目標(biāo)可簡化為磁偶極子[7]。在實測過程中，對多點根據(jù)磁偶極子場源理論進(jìn)行目標(biāo)定位時，測量結(jié)果如果相對穩(wěn)定，則表明此時目標(biāo)體可以當(dāng)做磁偶極子。Nara等[8]提出基于歐拉反褶積的磁偶極子單點定位算法，通過計算張量矩陣逆與磁場矢量的乘積來估算磁性目標(biāo)的位置。然而，實際應(yīng)用中很難對背景地磁場與目標(biāo)體的磁異常場進(jìn)行分離，因而對目標(biāo)體的定位精度受到限制。若想在地磁背景場中有效提取磁異常場源信息，可利用不同測點的場數(shù)據(jù)構(gòu)建線性方程組，求解位置矢量，即多點定位[9-10]，但是定位精度會受到航線選擇的限制，且對系統(tǒng)姿態(tài)及位置的精度要求很高。

為在單一測量點處實現(xiàn)磁異常目標(biāo)體的精確定位，本文對磁偶極子張量場和矢量場均滿足的歐拉反褶積方程求導(dǎo)，進(jìn)而構(gòu)建一、二階張量矩陣與位置矢量的關(guān)系方程，消除背景磁場，利用張量衍生不變關(guān)系得到目標(biāo)位置信息的限定方程；同時提出二階張量系統(tǒng)的概念以及測量方法，以此對磁性目標(biāo)進(jìn)行單點定位。由于差分測量二階張量數(shù)據(jù)對傳感器輸出誤差十分敏感，故需配套磁梯度張量系統(tǒng)校正技術(shù)以提升張量數(shù)據(jù)的解釋精度。目前已有較為成熟的校正方法，如李青竹等[11-12]提出的平面十字形磁梯度張量系統(tǒng)的兩步線性校正和基于橢球擬合的磁梯度張量系統(tǒng)集成校正方法、尹剛等[13]提出的磁梯度張量系統(tǒng)非線性校正方法等，均能較準(zhǔn)確地估計系統(tǒng)誤差參數(shù)，使單點定位方法可實現(xiàn)較高的目標(biāo)定位精度。

1 磁梯度張量要素與張量不變量

磁梯度張量為磁場分量在正交方向上的空間變化率，共有9個元素，其中無源靜態(tài)磁場中某點磁場矢量旋度和散度為零[14]，故張量矩陣對稱且無跡，可表示為

(1)

式中：Bx、By、Bz為磁場矢量的三軸分量；Bij(i,j=x,y,z)表示磁梯度張量分量，共5個獨立分量，一般由差分磁梯度張量系統(tǒng)測得，如由四個磁通門傳感器組成的平面十字形結(jié)構(gòu)張量系統(tǒng)[11-13]。張量運算時不隨坐標(biāo)系旋轉(zhuǎn)而改變的量，稱為張量不變量。一些基本不變量[15]如下

(2)

式中： tr(·)代表求矩陣的跡；det(·)表示求行列式的值；λ1、λ2和λ3為G的特征值，滿足特征方程λ3-I0λ2+I1λ-I2=0。

(3)

(4)

式中μ為介質(zhì)磁導(dǎo)率， 在空氣中μ≈μ0，μ0=4π×10-7N·A-2為真空磁導(dǎo)率。

2 磁偶極子場張量衍生不變關(guān)系

2.1 特征值分析

特征值λ1、λ2、λ3及特征方程系數(shù)I0、I1、I2與場源磁偶極矩有關(guān)。若m、r已知，聯(lián)立式(2)和式(4)，可推導(dǎo)張量系統(tǒng)測量點處磁梯度張量矩陣的特征值

(5)

式中： |λ1|≥|λ3|、|λ2|≥|λ3|且λ2≥λ3≥λ1。令

(6)

根據(jù)G的對稱性與無跡性，可推導(dǎo)張量矩陣特征值λ1、λ2、λ3對應(yīng)的特征向量

(7)

2.2 張量衍生不變關(guān)系

2.2.1 磁矩夾角不變關(guān)系

磁偶極子場源相對固定，當(dāng)張量系統(tǒng)姿態(tài)變換時，同一測量點各姿態(tài)下不同的張量測量數(shù)據(jù)必然對應(yīng)于恒定的磁矩矢量m和位置矢量r，因此m與r間的夾角θ0是不變的，將這個衍生不變關(guān)系稱為磁矩夾角不變關(guān)系，如圖1所示。假設(shè)該點處特征值λ1、λ2、λ3由實際張量系統(tǒng)測量并求解得到，則通過矢量的點積運算，聯(lián)立式(5)特征值解算式，可得到θ0的表達(dá)式

(8)

式中u為歸一化磁源強(qiáng)度(normalized source strength,NSS)。

圖1 磁矩矢量與位置矢量夾角不變關(guān)系

2.2.2 特征向量空間不變關(guān)系

磁梯度張量系統(tǒng)空間姿態(tài)方向可以是任意的，選擇空間坐標(biāo)系以磁偶極子為原點，r矢量方向為x軸正向，m與r的平面為坐標(biāo)系xoy平面，則磁矩矢量與位置矢量的夾角不變關(guān)系可用如圖2所示的空間直角坐標(biāo)系表示。

根據(jù)文獻(xiàn)[16]給出的偶極子位置表示方法，即磁矩矢量為m=(Mcosθ0,Msinθ0,0)T，位置矢量可表示為r=(r,0,0)T。聯(lián)立式(4)～式(8)，可推導(dǎo)出特征值λ1、λ2、λ3及其對應(yīng)的特征向量v1、v2、v3

(9)

(10)

由式(10)可知，在建立的空間直角坐標(biāo)系中，特征向量v3垂直于矢量m、r所在平面xoy，而特征向量v1、v2則與矢量m、r所在的平面xoy共面(圖2)。由于該空間坐標(biāo)系沒有對磁偶極子的磁矩信息、張量系統(tǒng)測量姿態(tài)施加任何約束，而只與坐標(biāo)系的選擇有關(guān)，因此可衍生為一般性結(jié)論：絕對值最小的特征值λ3對應(yīng)的特征向量v3垂直于磁偶極子的磁矩矢量m和位置矢量r所在平面；最大和最小的兩個特征值λ1、λ2對應(yīng)的特征向量v1、v2與磁偶極子的磁矩矢量m和位置矢量r所在平面共面。將上述兩個衍生不變關(guān)系稱為特征向量空間不變關(guān)系。

圖2 特征向量空間不變關(guān)系在直角坐標(biāo)系中的表示

以上的張量衍生不變關(guān)系均與系統(tǒng)姿態(tài)無關(guān)，而僅與目標(biāo)源磁矩矢量和位置矢量有關(guān)，故必然能為單點定位提供關(guān)于目標(biāo)位置信息的限定條件。

3 二階張量歐拉反褶積定位方法

3.1 二階磁梯度張量矩陣與系統(tǒng)

磁梯度張量矩陣G各分量是磁矢量分量對于三軸坐標(biāo)的一階偏導(dǎo)數(shù)，因此又稱為一階磁梯度張量。通過求解磁矢量在三軸方向上的二階偏導(dǎo)數(shù)，可得到二階磁梯度張量矩陣GII，共由27個元素組成

(11)

二階張量較一階張量可提取出更淺表的附加磁源信息，具有更強(qiáng)的磁源分辨能力。在無源靜態(tài)磁場中，GII的27個元素中僅有7個相互獨立，例如一組獨立元素為Bxxx、Bxyx、Bxzx、Byxy、Byyy、Byzy、Bzzz。以平面十字形磁梯度張量系統(tǒng)為原型，將其擴(kuò)展為二階張量系統(tǒng)，能夠差分測量出磁矢量分量的二階偏導(dǎo)數(shù)。設(shè)計的該二階系統(tǒng)結(jié)構(gòu)如圖3所示。

圖3 平面十字形二階磁梯度張量系統(tǒng)結(jié)構(gòu)示意圖

由圖3可知，在原有平面十字形結(jié)構(gòu)基礎(chǔ)上，在中心o點處添增加了一個傳感器5。該系統(tǒng)能夠測量得到二階張量矩陣27個元素中的6個和一階張量矩陣G

(12)

(13)

3.2 三維歐拉反褶積公式

利用歐拉反褶積法可以獲得場源位置信息。若n階齊次方程f(x,y,z)滿足[17]

(14)

該方程稱為歐拉方程。設(shè)地質(zhì)異常體的坐標(biāo)為(x0,y0,z0)，則在測量點(x,y,z)處的位場函數(shù)

(15)

=-Nf(x-x0,y-y0,z-z0)

(16)

當(dāng)空間坐標(biāo)系中磁異常場源中心位置坐標(biāo)為(0,0,0)，則磁異常場源位置矢量r=(-x,-y,-z)T，該場源產(chǎn)生的磁場矢量三分量為

Bi(x,y,z)=f(x-0,y-0,z-0)

(17)

式(17)滿足式(16)的位場歐拉方程形式

(18)

式(18)可寫為矩陣形式

(19)

式(19)即為三維歐拉反褶積公式。Nara定位法[8]推導(dǎo)了磁偶極子梯度差，得到了式(19)中N=3的特例。當(dāng)位場中僅存在磁異常場，對于給定待測場源的構(gòu)造指數(shù)N，測得點(x,y,z)處張量矩陣的9個元素及磁場矢量的3個元素，就能反解出場源位置矢量r=(-x,-y,-z)T。然而，實際位場不會僅僅包括磁異常場。由于此處磁場矢量(Bx、Bx、Bz)僅由場源產(chǎn)生，而張量系統(tǒng)測得的磁場矢量為磁異常場與背景磁場的矢量疊加，故Nara法對實測數(shù)據(jù)的定位精度必然受到影響。

3.3 磁性目標(biāo)位置的二階張量歐拉反褶積求解

為了從疊加場中有效分離出磁性體異常場，對式(18)中各等式兩邊分別對x、y、z求偏導(dǎo)

(20)

磁偶極子的構(gòu)造指數(shù)為3， 用于磁性目標(biāo)定位的式(20)可寫為矩陣乘積的形式

(21)

本文定義上式為二階張量歐拉反褶積公式，描述了磁偶極子在測量點處產(chǎn)生的一階、二階張量與其位置矢量之間的關(guān)系?；?.1節(jié)中的二階張量系統(tǒng)可測得式(21)中二階張量矩陣的前6個元素和一階張量分量Bxx、Byy、Bzz。然而，位置矢量未知數(shù)為3， 則已知的兩個方程構(gòu)成的線性方程組是欠定的。

由張量衍生不變關(guān)系可知，磁偶極子場源下張量矩陣G的絕對值最小特征值λ3對應(yīng)的特征向量v3垂直于磁偶極子磁矩矢量m和位置矢量r所在平面，有

v3·r=[v3xv3yv3z][xyz]T=0

(22)

則式(22)可拓展為

(23)

式中Gv、gv均可通過二階張量系統(tǒng)測量后求解得到，據(jù)此可解得磁偶極子的位置矢量

(24)

再利用式(4)可解得磁矩矢量m。該方法僅需利用單一測量點的一階、二階張量數(shù)據(jù)便能實現(xiàn)磁偶極子的單點定位。

4 仿真分析

將磁偶極子置于真空空間(8m,5m,-4m)處，磁矩模為8000A·m-2，磁偶極子磁偏角為30°，磁傾角為40°。假設(shè)地磁背景場為勻強(qiáng)磁場，地磁總場強(qiáng)度(Total magnetic intensity，TMI)為55000nT，磁傾角為60°，磁偏角為-7°(西偏)。模擬搭建一個二階平面十字形磁梯度張量系統(tǒng)[11]進(jìn)行測量，基線距離d設(shè)為0.4m。單航線測量從(-20m,0,0)到(40m,0,0)，測量試驗如圖4。

圖4 磁偶極子場源張量系統(tǒng)單航線測量試驗示意圖

航線運動過程保持系統(tǒng)姿態(tài)不變，采樣間隔為1m。在測量航線上各測量點處測得的磁場矢量三分量值如圖5所示?？梢娫邶嫶蟮牡卮艌霰尘跋拢纱排紭O子產(chǎn)生的部分磁場矢量分量場基本被淹沒，磁異常場難以有效分離。由文獻(xiàn)[11]知，磁傳感器實際輸出可表示為

(25)

式中：B1為傳感器實際輸出；ix、iy、iz為零位偏差；cx、cy、cz為靈敏度標(biāo)度因子；φ、θ、ψ為非正交角；α、β、γ為非對準(zhǔn)誤差角；w為測量噪聲；B2為理想輸出。利用圖3中的5個傳感器的位置仿真出該點處的航線測量數(shù)據(jù)，共三組作對比。

圖5 單航線測量中磁場分量

①理想組：誤差為零、測量噪聲為零；

②實際組：隨機(jī)預(yù)設(shè)各傳感器的12個誤差參數(shù)，并加入均值為0、方差為1nT的高斯噪聲，預(yù)設(shè)參數(shù)值列于表1中；

③校正組：使用兩步線性校正方法[11]對②中預(yù)設(shè)誤差后的各傳感器的12種誤差參數(shù)進(jìn)行估計，并以此對實際組中的測量值進(jìn)行校正。

由式(1)可知G為對稱矩陣，但式(13)中Bxy與Byx測量值存在差異，可視其為兩個獨立分量，此時G中共有6個獨立分量。利用本文方法對航線上所有測量點的①、②、③組測量數(shù)據(jù)分別計算磁偶極子位置坐標(biāo)，則該測量航線三組測量數(shù)據(jù)下該二階張量系統(tǒng)測得的各測量點的一階張量的6個獨立分量、二階張量的6個獨立分量、估計得到的磁性目標(biāo)空間位置圖、磁性目標(biāo)位置的坐標(biāo)分別如圖6、圖7和圖8所示。

圖6 各傳感器不加誤差、噪聲時的單航線張量測量結(jié)果及目標(biāo)定位結(jié)果

圖7 各傳感器加入誤差、測量噪聲時的單航線張量測量結(jié)果及目標(biāo)定位結(jié)果

圖8 各傳感器進(jìn)行兩步線性校正后的單航線張量測量結(jié)果及目標(biāo)定位結(jié)果

表1 預(yù)設(shè)5個傳感器系統(tǒng)誤差與非對準(zhǔn)誤差

參考文獻(xiàn)[11]中的仿真結(jié)果，在相同測量工況的12個誤差參數(shù)的預(yù)設(shè)條件下，都引入了均值為0、方差為1nT的高斯噪聲，由第③步中估計的各個傳感器的誤差接近于0。由三組測量數(shù)據(jù)的定位結(jié)果可知： ①當(dāng)未加入噪聲和誤差時，該方法對磁偶極子場源的位置估計誤差接近0，且在航線上所有點位均能實現(xiàn)有效探測并估計出準(zhǔn)確的場源坐標(biāo)；②各傳感器引入12個誤差和預(yù)設(shè)噪聲時，在該航線上無法有效探測到磁源目標(biāo)；③加入均值為0、方差為1nT的高斯噪聲后，該二階張量系統(tǒng)在(-5m,0,0)～(20m,0,0)范圍內(nèi)基本實現(xiàn)了有效探測及精確定位，超出該范圍的點無法被有效探測。

5 實驗驗證

為驗證方法對實際目標(biāo)體的定位準(zhǔn)確性及張量系統(tǒng)校正對定位精度的提升能力，設(shè)計如下實驗。

①搭建一個一階平面十字形磁梯度張量系統(tǒng)(圖9)，基線距離d為0.5m。二階張量系統(tǒng)在中心點o處增加傳感器5。為最大限度避免結(jié)構(gòu)誤差，在單航線測量中可利用同軸其他傳感器對傳感器5讀數(shù)進(jìn)行間接測量，從而獲得二階張量數(shù)據(jù)。間接測量示意圖見圖10。

②使用標(biāo)量質(zhì)子磁強(qiáng)計選擇地磁場穩(wěn)定的某測區(qū)，面積為2.1m×2.1m。系統(tǒng)置于滑動臺架上，測量平面高0.65m。測得該測區(qū)內(nèi)平均TMI標(biāo)量為53911.48nT。盡管臺架盡量選擇無磁材質(zhì)，但局部連接處仍存在磁干擾。為盡量避免磁干擾影響定位精度，利用橢球擬合集成補(bǔ)償方法[12]對測區(qū)內(nèi)的張量系統(tǒng)進(jìn)行校正參數(shù)估計，得到一組適用于該區(qū)域的系統(tǒng)各傳感器誤差補(bǔ)償參數(shù)(集成誤差補(bǔ)償系數(shù)F和集成零偏向量Iw)和對準(zhǔn)參數(shù)(旋轉(zhuǎn)矩陣T)，見表2。該方法屬于間接校正法，不需要估計傳感器具體參數(shù)便可實現(xiàn)磁干擾環(huán)境下的磁梯度張量系統(tǒng)集成校正，對測量環(huán)境的硬、軟磁干擾具有良好的補(bǔ)償效果，同時可有效消除傳感器系統(tǒng)誤差和非對準(zhǔn)誤差。磁傳感器輸出補(bǔ)償公式[12]為

圖9 平面十字形磁梯度張量系統(tǒng)示意圖

圖10 傳感器5的間接測量示意圖

Bc=TF(Br-Iw)

(26)

式中：Bc為理想傳感器輸出；Br為待補(bǔ)償?shù)膶嶋H傳感器輸出。

③以測量平面xoy建立空間直角坐標(biāo)系，x軸向東為正，z軸向上為正。以一小型磁鐵為探測目標(biāo)，張量系統(tǒng)單航線測量方向自西向東，起始坐標(biāo)為(0,0,0)，終點坐標(biāo)為(2.1m,0,0)，目標(biāo)磁鐵坐標(biāo)為(1.10m,0.50m,-0.65m)，采樣間隔為0.05m，則該航線共有43個測量點，傳感器5的讀數(shù)間接等效于第6至第43測點的傳感器4讀數(shù)和第34至38測點的傳感器2讀數(shù)。該單航線測量實驗裝置示意圖見圖11。

由于探測距離遠(yuǎn)大于磁鐵尺度的2.5倍，因此將其視為磁偶極子，并進(jìn)行目標(biāo)單點定位計算。首先，將步驟③中全部原始數(shù)據(jù)直接進(jìn)行磁鐵位置坐標(biāo)演算；然后利用表2中參數(shù)對步驟③中全部測量數(shù)據(jù)進(jìn)行校正處理，再演算磁鐵坐標(biāo)。校正前、后估計的磁鐵位置如圖12所示。將校正前、后估算坐標(biāo)的精度使用均方根誤差(Root mean square error,RMSE)[19]進(jìn)行量化

圖11 單航線磁性目標(biāo)定位實驗裝置示意圖

傳感器1傳感器2傳感器3傳感器4F0.9968 0.0058 0.01230.0058 1.0058 -0.00900.0123 -0.0090 1.0089é?êêêù?úúú 1.0098 -0.0069 0.0176-0.0069 0.9899 0.0084 0.0176 0.0084 1.0108é?êêêù?úúú1.0052 0.0116 0.00700.0116 0.9961 -0.01640.0070 -0.0164 1.0113é?êêêù?úúú 1.0088 -0.0121 0.0016-0.0121 0.9967 0.0008 0.0016 0.0008 1.0072é?êêêù?úúúIw/nT[48.51 35.17 -24.39]T[9.86 -27.77 12.38]T[-24.39 -4.98 19.66]T[-28.87 32.48 10.87]TT- 0.9986 -0.0471 0.0244 0.0462 0.9983 0.0366-0.0261 -0.0355 0.9990é?êêêù?úúú0.9988 -0.0313 -0.03700.0297 0.9986 -0.04470.0384 0.0436 0.9983é?êêêù?úúú 0.9984 0.0436 0.0367-0.0428 0.9988 -0.0227-0.0376 0.0211 0.9991é?êêêù?úúú

圖12 經(jīng)橢球擬合集成校正前、后對磁鐵的定位結(jié)果

(27)

實驗結(jié)果表明，在該實驗工況下，校正后的系統(tǒng)對磁鐵坐標(biāo)的估計精度可控制在均方根誤差10cm以內(nèi)(表3)。當(dāng)背景場更穩(wěn)定、系統(tǒng)校正精度更高時，理論上能實現(xiàn)磁性目標(biāo)厘米級別以上的精確定位。系統(tǒng)誤差校正技術(shù)的應(yīng)用對提升張量數(shù)據(jù)解釋精度尤為關(guān)鍵。

表3 系統(tǒng)校正前、后磁鐵目標(biāo)的坐標(biāo)估計精度對比

6 結(jié)論

(1)推導(dǎo)的張量衍生不變關(guān)系表明：磁偶極子磁矩矢量與測量點的位置矢量夾角是恒定的，且僅與該點處張量矩陣特征值有關(guān)；磁偶極子場源的張量矩陣中間特征值對應(yīng)的特征向量垂直于磁矩矢量與位置矢量，最大、最小特征值對應(yīng)的特征向量與磁矩矢量和位置矢量共面。

(2)提出了平面十字形二階張量系統(tǒng)概念與測量方法，對三維歐拉反褶積公式求導(dǎo)后，利用一階、二階張量與張量衍生不變關(guān)系可求解磁源位置矢量。

(3)相比Nara定位法，本文方法能將磁異常信息與背景地磁場有效分離，坐標(biāo)估計值唯一，在有效探測范圍內(nèi)能實現(xiàn)測區(qū)中任意點對磁性目標(biāo)的單點定位。

(4)仿真和實驗結(jié)果均表明：針對小尺度磁鐵的定位實驗中，在該實驗工況下，經(jīng)誤差校正后的磁梯度張量系統(tǒng)的定位精度控制在均方根誤差10cm以內(nèi)。

基于二階磁張量歐拉反褶積的磁源單點定位方法同樣適用于其他結(jié)構(gòu)類型的二階磁梯度張量系統(tǒng)，且在實測中利用誤差校正技術(shù)極大地提升了磁性目標(biāo)的定位精度，表明誤差校正對于提升張量數(shù)據(jù)解釋精度的重要性。高精度目標(biāo)單點定位技術(shù)與張量系統(tǒng)誤差校正技術(shù)的結(jié)合，為后期針對磁異常體的反演識別、三維重塑等更深層次張量數(shù)據(jù)解釋提供了理論基礎(chǔ)與經(jīng)驗參考。