孫赫軒裴東興*祗會(huì)強(qiáng)雒茁君

(1.中北大學(xué)儀器科學(xué)與動(dòng)態(tài)測(cè)試教育部重點(diǎn)實(shí)驗(yàn)室，山西 太原030051；2.國網(wǎng)山西省電力公司電力科學(xué)研究院，山西 太原030051)

磁異常探測(cè)技術(shù)被廣泛應(yīng)用于地質(zhì)勘測(cè)、變形監(jiān)測(cè)、未爆炸物檢測(cè)、近岸目標(biāo)檢測(cè)、醫(yī)療內(nèi)鏡定位、水下目標(biāo)體定位等場(chǎng)景[1-5]，具有極高的民用價(jià)值與軍事意義。 磁性目標(biāo)產(chǎn)生明顯的磁場(chǎng)異常是磁異常探測(cè)的前提，而描述目標(biāo)磁性特征最直觀重要的物理量就是磁矩[6]。 在排雷或醫(yī)療等相關(guān)領(lǐng)域中應(yīng)用磁異常探測(cè)技術(shù)時(shí)，磁矩計(jì)算誤差引發(fā)的一系列偏差很可能會(huì)造成無法挽回的嚴(yán)重事故。 磁梯度張量相比于以往的磁場(chǎng)強(qiáng)度包含了更豐富的磁場(chǎng)信息，同時(shí)有效地克服了地磁場(chǎng)的影響，成為了國內(nèi)外研究的熱點(diǎn)[7-8]。

周建軍[9]從磁性產(chǎn)生機(jī)理出發(fā)，推導(dǎo)了鐵磁體磁矩計(jì)算公式，該方法使用起來較為繁瑣，且僅適用于粗略估算。 Shutyǐ[10]對(duì)具有偶極磁矩的三球體和四球體系統(tǒng)進(jìn)行了數(shù)值分析，研究了控制由交變磁場(chǎng)的幅值或頻率變化引起的感應(yīng)磁矩的可能性。Zhou[11]驗(yàn)證了感應(yīng)磁矩與地磁場(chǎng)成正比，并設(shè)計(jì)了一種鐵磁力矩的實(shí)時(shí)補(bǔ)償方法，該方法可用于凈化航磁測(cè)量中的磁環(huán)境。 Inamori T[12]使用地磁場(chǎng)信息計(jì)算磁矩，進(jìn)而對(duì)衛(wèi)星姿態(tài)進(jìn)行磁補(bǔ)償，使姿態(tài)穩(wěn)定精度提高了5 倍。 Nara[13]研究了磁矩與位置矢量的夾角對(duì)磁矩的影響。 楊慶宇[6]、姜浩[14]等人優(yōu)化了磁矩測(cè)量計(jì)算方法，原理上都基于磁場(chǎng)強(qiáng)度的測(cè)量，沒有考慮到測(cè)量點(diǎn)處的磁場(chǎng)強(qiáng)度中疊加了地磁場(chǎng)值，這給磁矩計(jì)算帶來很大誤差。 周[15]提出了磁場(chǎng)積分法和最小二乘正則化方法反演鐵磁物體的固有磁矩，初步解決了固有磁矩引起的感應(yīng)磁矩問題。 Chadebec O[16]將磁矩量法與有限元法相比，認(rèn)為該方法具有較高的分辨率和高精度的復(fù)雜場(chǎng)計(jì)算。 但使用這種方法要對(duì)物理現(xiàn)象和數(shù)值方法本身有良好的了解，需要較高的操作技巧。

本文提出了一種基于磁梯度張量與Levenberg-Marquardt 優(yōu)化的磁矩計(jì)算方法。 該方法利用磁偶極子模型，采用磁梯度張量概念與Levenberg-Marquardt 優(yōu)化算法(下文簡(jiǎn)稱L-M 算法)結(jié)合的方式，有效消除地磁場(chǎng)的干擾與測(cè)量計(jì)算奇異點(diǎn)，提高磁矩計(jì)算精度與抗噪聲能力，具有更高的魯棒性。 為后期針對(duì)磁性目標(biāo)的反演識(shí)別、三維重建等更深層次磁異常數(shù)據(jù)的應(yīng)用提供了理論基礎(chǔ)與經(jīng)驗(yàn)參考。

1 磁偶極子模型與磁梯度張量理論

為簡(jiǎn)化磁場(chǎng)問題，可以將磁場(chǎng)目標(biāo)抽象為一個(gè)僅具有磁性特性相關(guān)的理想化模型——磁偶極子。理論上講，任意復(fù)雜源都可視為由一定數(shù)量偶極子組成，這為復(fù)雜源的加入及復(fù)雜環(huán)境下的瞬態(tài)場(chǎng)計(jì)算提供條件[17]。 測(cè)量距離約為被測(cè)物尺寸2.5 倍以上時(shí)，磁偶極子模型可用于替代形狀較為規(guī)則的磁性物體[18]。 假設(shè)磁偶極子磁矩為m，則在距離偶極子r處的磁感應(yīng)強(qiáng)度可以表示為:

式中:μ0為真空磁導(dǎo)率，ro為位移r的方向向量，r＝｜r｜，mo為磁矩m的方向向量，m＝｜m｜。

磁梯度張量是磁場(chǎng)矢量正交三分量在空間各個(gè)方向上的變化率構(gòu)成的二階張量[19]，共九個(gè)分量，可表示如下:

地磁場(chǎng)強(qiáng)度約為50 000 nT，且會(huì)隨空間時(shí)間發(fā)生變化，但短時(shí)間小區(qū)域內(nèi)可認(rèn)為是一個(gè)恒定磁場(chǎng)，其梯度值很小[20-21]，約小于0.02 nT/m。 在無源靜磁場(chǎng)下，靜磁場(chǎng)磁感應(yīng)強(qiáng)度的散度和旋度都等于零[22]，故磁梯度張量G可視為僅含有Bxx、Bxy、Byy、Bzx、Bzy5 個(gè)獨(dú)立分量。 即:

2 磁矩測(cè)量模型與優(yōu)化算法

2.1 基于十字形磁傳感器陣列的磁矩測(cè)量模型

磁性目標(biāo)的磁矩是一個(gè)既有大小又有方向的矢量，是描述物體磁性特征最直觀的物理量。 磁矩矢量可以寫做磁矩模值m與磁矩方向向量mo的乘積，也可依照坐標(biāo)系方向分解為正交三分量mx、my、mz。

磁偶極子磁矩測(cè)量模型示意圖見圖1，設(shè)磁矩為m＝(mx，my，mz)的磁偶極子位于坐標(biāo)系原點(diǎn)，則根據(jù)式(1)和式(2)可得在距離磁偶極子r＝(x，y，z)處的測(cè)量點(diǎn)P的磁梯度張量5 個(gè)獨(dú)立元素可表示為:

圖1 磁偶極子磁矩測(cè)量模型

磁梯度張量測(cè)量系統(tǒng)采用十字架結(jié)構(gòu)，如圖2所示，測(cè)量中心位于十字架中心，XY軸分別與十字架兩個(gè)軸平行，Z軸向上垂直于十字架所在平面，四個(gè)磁傳感器分別位于十字架四個(gè)頂點(diǎn)，傳感器的坐標(biāo)軸與測(cè)量系統(tǒng)坐標(biāo)軸保持一致。

圖2 十字形磁傳感器陣列

利用差分代替微分的思想，可得磁場(chǎng)中測(cè)量中心處的磁梯度張量獨(dú)立分量值為:

式中:Sij代表傳感器i在j軸方向上的磁場(chǎng)分量，d為十字架基線長(zhǎng)度。 將以上獨(dú)立分量代入式(3)可得測(cè)量中心點(diǎn)處的磁梯度張量。

在同一磁場(chǎng)環(huán)境下聯(lián)立式(4)和式(5)中任意三個(gè)分量，可以使用解析方法計(jì)算出磁矩矢量m＝(mx，my，mz)的值。 之后可以解出磁矩模值m和方向矢量mo:

然而求解過程中發(fā)現(xiàn)，不管如何選取組合其中的三個(gè)分量，在某些位置都會(huì)出現(xiàn)方程組自由度不足的現(xiàn)象，導(dǎo)致出現(xiàn)奇異點(diǎn)。

2.2 基于Levenberg-Marquardt 算法的磁矩解算優(yōu)化

針對(duì)解析方法求解會(huì)出現(xiàn)奇異點(diǎn)的問題，考慮使用L-M 優(yōu)化算法對(duì)磁矩相關(guān)數(shù)據(jù)進(jìn)行求解。 LM 算法是一種迭代優(yōu)化算法，綜合了高斯牛頓法局部收斂和梯度法全局收斂的特點(diǎn)，解決了梯度法收斂速度慢的問題，也優(yōu)化了高斯牛頓法的Hessian矩陣不正定的問題[23-24]。

以磁梯度張量測(cè)試系統(tǒng)測(cè)的測(cè)量點(diǎn)張量數(shù)據(jù)作為范例，磁偶極子模型計(jì)算的包含磁源磁矩信息的磁梯度張量數(shù)據(jù)作為目標(biāo)數(shù)據(jù)，運(yùn)用L-M 優(yōu)化算法進(jìn)行優(yōu)化處理。

令F(mx，my，mz)為目標(biāo)函數(shù)，形式見下式:

利用L-M 算法對(duì)該非線性目標(biāo)函數(shù)進(jìn)行參數(shù)優(yōu)化，求(mx，my，mz)使得輸出的如下最小二乘表達(dá)式成立:

設(shè)定初值對(duì)目標(biāo)函數(shù)做非線性最小乘二擬合，所得(mx，my，mz)即為磁性目標(biāo)的磁矩優(yōu)化計(jì)算結(jié)果。

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

為驗(yàn)證所提方法的反演性能，本文使用MATLAB軟件進(jìn)行多組仿真實(shí)驗(yàn)，仿真對(duì)象針對(duì)可能影響反演結(jié)果的地磁背景場(chǎng)、測(cè)量系統(tǒng)所在位置、環(huán)境噪聲等因素。 仿真中預(yù)設(shè)磁性目標(biāo)磁矩模值為1 000 Am2，磁傾角為π/4，磁偏角為π/6。 地磁場(chǎng)取經(jīng)驗(yàn)?zāi)Ｖ?0 000 nT，地磁傾取-0.925，地磁偏角-0.30進(jìn)行模擬。 十字架磁梯度張量測(cè)量系統(tǒng)基線為0.4 m，傳感器精度為10 nT。 為避免偶然性，測(cè)量系統(tǒng)中心在Z＝3 m 的平面內(nèi)半徑為5 m 的圓形軌跡上等間隔取64個(gè)測(cè)量點(diǎn)，對(duì)位于坐標(biāo)系原點(diǎn)處的磁源點(diǎn)磁場(chǎng)信息進(jìn)行測(cè)量，測(cè)量系統(tǒng)位置示意圖見圖3。

圖3 測(cè)量系統(tǒng)位置示意圖

利用測(cè)量點(diǎn)的磁場(chǎng)信息，分別使用姜浩等人提出的基于磁場(chǎng)強(qiáng)度B 的磁矩計(jì)算方法(簡(jiǎn)稱MOB)、基于磁梯度張量的磁矩計(jì)算解析方法(簡(jiǎn)稱MOG)、和本文所提優(yōu)化方法對(duì)磁矩進(jìn)行計(jì)算，比較分析其磁矩反演效果。 仿真實(shí)驗(yàn)中，采用計(jì)算方向向量與理論方向向量之間的夾角度數(shù)來作為衡量磁矩方向反演效果的評(píng)價(jià)指標(biāo)，采用計(jì)算模值與理論模值之間的相對(duì)誤差作為衡量磁矩大小反演效果的評(píng)價(jià)指標(biāo)。

3.1 地磁背景場(chǎng)影響分析

在地球磁場(chǎng)背景下，磁傳感器測(cè)得的磁場(chǎng)強(qiáng)度為磁異常信息與地球磁場(chǎng)的混疊。 由于地磁場(chǎng)強(qiáng)度的數(shù)量級(jí)遠(yuǎn)大于磁異常強(qiáng)度，如果不能很好的分離地磁場(chǎng)與目標(biāo)信號(hào)，將會(huì)對(duì)磁矩反演造成困難與誤差。 本組仿真對(duì)同一組傳感器在同一時(shí)間測(cè)量到的磁場(chǎng)數(shù)據(jù)分別采用方法B 與本文所提方法進(jìn)行磁矩求解。 仿真結(jié)果見圖4 與圖5。

圖4 測(cè)量點(diǎn)的磁矩方向向量角度誤差

由圖4 和圖5 可知，基于磁感應(yīng)強(qiáng)度的磁矩求解結(jié)果與地磁場(chǎng)方向有明顯的相關(guān)關(guān)系，地磁場(chǎng)方向上誤差值最大，說明此方法受地磁場(chǎng)影響極大，無法將磁源信息從背景場(chǎng)中提取出來。 而利用磁梯度張量方法的結(jié)果與地磁場(chǎng)無明顯相關(guān)關(guān)系且角度誤差與相對(duì)誤差都較小，有效解決了磁測(cè)環(huán)境背景地磁場(chǎng)干擾較大且難以分離的問題。

圖5 測(cè)量點(diǎn)的磁矩模值相對(duì)誤差

3.2 測(cè)量系統(tǒng)位置影響分析

方法應(yīng)用中發(fā)現(xiàn)，測(cè)量系統(tǒng)與磁源點(diǎn)相對(duì)位置不同，磁矩反演效果也可能不同。 由于上一節(jié)證實(shí)方法B 存在明顯缺陷，故在下文主要對(duì)基于磁梯度張量的磁矩計(jì)算方法進(jìn)行探究。 本組仿真針對(duì)不同的測(cè)量系統(tǒng)位置，分別采用解析方法與L-M 優(yōu)化算法對(duì)磁梯度張量測(cè)量系統(tǒng)的數(shù)據(jù)進(jìn)行后續(xù)處理，分析兩種方法受測(cè)量系統(tǒng)位置的影響情況，仿真結(jié)果見圖6 與圖7。

圖6 測(cè)量點(diǎn)的磁矩方向向量角度誤差

圖7 測(cè)量點(diǎn)的磁矩模值相對(duì)誤差

由圖6 和圖7 可知，使用解析方法處理梯度張量數(shù)據(jù)時(shí)，在特定坐標(biāo)軸附近存在誤差明顯異于常值的現(xiàn)象，而采用了L-M 優(yōu)化算法的方法則很好的解決了這一問題，數(shù)據(jù)波動(dòng)平穩(wěn)，誤差分布均勻，處理結(jié)果非常理想，使磁矩在磁源點(diǎn)附近有效距離內(nèi)全方位精準(zhǔn)可測(cè)。

3.3 環(huán)境噪聲影響分析

實(shí)際測(cè)量中存在的環(huán)境噪聲會(huì)導(dǎo)致磁矩反演困難，因此抗噪聲能力是衡量一個(gè)方法實(shí)用性的重要指標(biāo)。 在傳感器測(cè)量值中加入不同信噪比的高斯白噪聲。 高斯白噪聲存在隨機(jī)性，為觀測(cè)普遍規(guī)律，保持其他條件不變，每種信噪比情況下選取64 個(gè)測(cè)量點(diǎn)的數(shù)據(jù)，并對(duì)64 組數(shù)據(jù)求取均值以代表此信噪比下解析方法與L-M 優(yōu)化方法的反演效果。 仿真結(jié)果見圖8 和圖9。

圖8 不同信噪比下的磁矩方向角度誤差

圖9 不同信噪比下的磁矩模值相對(duì)誤差

由圖8 和圖9 可知，隨著高斯噪聲信噪比的增大，磁矩計(jì)算的誤差迅速減小。 綜合來看，當(dāng)信噪比大于等于60 dB 時(shí)，反演誤差達(dá)到最小并趨于平穩(wěn)。對(duì)比兩種方法，所提方法的數(shù)據(jù)總是先于解析方法平穩(wěn)，且總體誤差更小，抗干擾能力更強(qiáng)。 這是由于所提方法在解算過程中使用數(shù)值法計(jì)算完成，對(duì)噪聲有更強(qiáng)的魯棒性。

3.4 實(shí)驗(yàn)驗(yàn)證分析

選擇開闊無擾動(dòng)、地磁場(chǎng)穩(wěn)定的室外環(huán)境，劃定空間坐標(biāo)系，X軸方向?yàn)榈乩頄|向，Y軸方向?yàn)榈乩砟舷?，Z軸垂直向下。 已知尺寸為100 mm×100 m×50 mm，出廠磁矩參數(shù)范圍為510 ～525Am2，方向沿A 面垂直向上的鐵磁體安裝在原點(diǎn)處，由于本文算法基于磁偶極子模型，校準(zhǔn)后的磁梯度張量測(cè)量系統(tǒng)需位于磁體0.25 m 之外的位置。 磁性目標(biāo)姿態(tài)變換三次，本次實(shí)驗(yàn)測(cè)量系統(tǒng)的位置在距磁鐵直線距離為1.5 m～3.4 m 的范圍內(nèi)，每個(gè)磁鐵姿態(tài)下隨機(jī)放置40 個(gè)測(cè)試系統(tǒng)位置，采用上位機(jī)系統(tǒng)實(shí)時(shí)采集測(cè)量點(diǎn)處的多組磁場(chǎng)數(shù)據(jù)。 分別使用MOB、MOG和本文所提方法計(jì)算磁矩，并與出廠磁矩參數(shù)進(jìn)行對(duì)比，驗(yàn)證本文方法的可行性；在每一種磁性目標(biāo)姿態(tài)下多次變換測(cè)量系統(tǒng)位置，查看磁矩計(jì)算結(jié)果的變化情況，驗(yàn)證所提方法的可重復(fù)性。

由表1 可知，實(shí)驗(yàn)結(jié)果基本與仿真情況相符。MOB 方法混疊了地磁場(chǎng)，誤差極大。 MOG 方法測(cè)量結(jié)果處于可接受范圍內(nèi)，但穩(wěn)定性較差，個(gè)別位置偏差較大，在少量樣本情況下需慎重使用。 使用本文方法對(duì)三組姿態(tài)磁體反演得到的磁矩幅值數(shù)據(jù)穩(wěn)定，且處于廠家給出的磁矩范圍內(nèi)，認(rèn)為實(shí)測(cè)結(jié)果有效，且數(shù)據(jù)波動(dòng)范圍可接受，能夠滿足工程需要。 偏角誤差的實(shí)測(cè)值比仿真值要大一些，下一步實(shí)驗(yàn)擬使用精度更高的磁傳感器，設(shè)計(jì)更優(yōu)良的測(cè)量系統(tǒng)，選擇更理想的測(cè)試環(huán)境，預(yù)計(jì)可進(jìn)一步提高所提方法的應(yīng)用精度。

圖10 基于磁梯度張量測(cè)量系統(tǒng)的磁矩反演實(shí)驗(yàn)

表1 磁矩計(jì)算平均相對(duì)誤差與平均標(biāo)準(zhǔn)差表

4 結(jié)論

本文提出了一種基于磁梯度張量與Levenberg-Marquardt 優(yōu)化的磁矩計(jì)算方法，優(yōu)化了現(xiàn)有的磁矩計(jì)算方法使用不靈活、精度較低、受地球背景磁場(chǎng)影響大的問題。 該方法采用磁梯度張量概念消除了地磁背景場(chǎng)的干擾，采用L-M 優(yōu)化算法解決了解析方法求解存在奇異點(diǎn)，誤差相對(duì)較大的弊端。 仿真實(shí)驗(yàn)與現(xiàn)場(chǎng)實(shí)驗(yàn)的結(jié)果表明，本文所提磁矩計(jì)算方法相比于傳統(tǒng)方法消除了地磁背景場(chǎng)干擾，解決了奇異點(diǎn)問題，可以實(shí)現(xiàn)地磁環(huán)境和一定信噪比的高斯環(huán)境噪聲情況下的磁矩快速反演。 整體測(cè)量方案方便簡(jiǎn)單，精準(zhǔn)易用，具有較強(qiáng)的可行性與可靠性。 所提方法同樣適用于其他結(jié)構(gòu)類型的磁梯度張量測(cè)量系統(tǒng)。 高精度抗噪聲的磁矩計(jì)算技術(shù)為后期針對(duì)磁性目標(biāo)的反演識(shí)別、三維重建等更深層次磁異常數(shù)據(jù)的應(yīng)用提供了理論基礎(chǔ)與經(jīng)驗(yàn)參考。