孫赫軒裴東興*祗會強雒茁君

(1.中北大學儀器科學與動態(tài)測試教育部重點實驗室，山西 太原030051；2.國網山西省電力公司電力科學研究院，山西 太原030051)

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

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

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

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

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

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

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

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

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

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

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

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

圖1 磁偶極子磁矩測量模型

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

圖2 十字形磁傳感器陣列

利用差分代替微分的思想，可得磁場中測量中心處的磁梯度張量獨立分量值為:

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

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

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

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

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

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

令F(mx，my，mz)為目標函數，形式見下式:

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

設定初值對目標函數做非線性最小乘二擬合，所得(mx，my，mz)即為磁性目標的磁矩優(yōu)化計算結果。

3 仿真分析與實驗驗證

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

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

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

3.1 地磁背景場影響分析

在地球磁場背景下，磁傳感器測得的磁場強度為磁異常信息與地球磁場的混疊。 由于地磁場強度的數量級遠大于磁異常強度，如果不能很好的分離地磁場與目標信號，將會對磁矩反演造成困難與誤差。 本組仿真對同一組傳感器在同一時間測量到的磁場數據分別采用方法B 與本文所提方法進行磁矩求解。 仿真結果見圖4 與圖5。

圖4 測量點的磁矩方向向量角度誤差

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

圖5 測量點的磁矩模值相對誤差

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

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

圖6 測量點的磁矩方向向量角度誤差

圖7 測量點的磁矩模值相對誤差

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

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

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

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

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

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

3.4 實驗驗證分析

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

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

圖10 基于磁梯度張量測量系統(tǒng)的磁矩反演實驗

表1 磁矩計算平均相對誤差與平均標準差表

4 結論

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