改進的磁源體邊界識別方法

張 琪， 張英堂， 李志寧， 范紅波

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

邊界識別是磁異常數據解釋的重要組成部分,其中單個磁源體的邊界描繪方法有解析信號模的極大值[1]、水平導數的極大值[2]和垂向導數的零值[3]等。但當存在多個不同埋深的磁源體時，深部弱異常體的磁異常往往會被淺部強異常體所掩蓋，難以同時描繪各磁源體的邊界。傾斜角法[4-6]和Theta圖法[7-8]雖然能夠描繪不同埋深的磁源體的邊界，但識別結果趨于模糊、發散，限制了其進一步發展。場源體的邊界可通過張量數據的特征值描繪，如：ORUC等[9]利用曲率重力梯度張量的最大特征值來解釋Erzurum盆地的地質結構；ZHOU等[10]對其進行了改進，將重力異常數據結合到最大特征值公式中，以描繪具有正負重力異常場源體的邊界。然而，與重力異常簡單的對應關系不同，磁異常對磁化方向十分敏感，這給準確識別磁源體的邊界增加了難度。

鑒于此，筆者首先利用化極算法消除斜磁化帶來的不利影響，然后采用基于曲率磁梯度張量矩陣的最大特征值的歸一化改進算法，提取得到多個磁源體的特征邊界等值線，最后通過對離散的邊界等值點進行最小二乘擬合，得到最接近磁源體真實邊界的連續曲線。

1 邊界識別理論

1.1 磁梯度張量

磁梯度張量表示磁場三分量在三個互相正交方向上的空間變化率，其矩陣形式為

(1)

式中：G為磁梯度張量矩陣；Bi(i=x,y,z)為磁場矢量；Bij(i,j=x,y,z)為磁梯度張量分量。

1.2 曲率磁梯度張量

曲率磁梯度張量矩陣C為

(2)

求解其特征值矩陣，可得

(3)

(4)

式中：λ1和λ2分別表示曲率磁梯度張量矩陣的最大特征值和最小特征值。

在此基礎上，筆者提出了基于曲率磁梯度張量矩陣的最大特征值的歸一化改進算法，并將其定義為λ，其表達式為

(5)

式中：λ1_max為λ1的最大值;k1的取值為0.001～0.1。

1.3 傾斜角法

傾斜角是位場總水平導數和垂向導數比值的反正切函數，其零值等值線可突出場源體的邊界特征。張雙喜等[6]對傾斜角進行了改進，將改進形式定義為T1，其表達式為

(6)

式中：U為磁標量勢。

將式(6)拓展成張量形式，并定義為T2，其表達式為

(7)

將增強型傾斜角定義為T3，其表達式為

(8)

式中：k=min(|Bz|)/max(Bzz)

1.4 Theta圖法

Theta圖法的最大值等值線(趨近于1)可識別場源體的邊界。對Theta圖法進行改進[5]，得到4種改進形式，將其依次定義為θ1、θ2、θ3和θ4，其表達式分別為

(9)

(10)

(11)

(12)

1.5 最小二乘法

受實際探測方式的影響，利用實測磁異常數據計算得到的邊界等值線通常由一系列離散的邊界等值點組成，其離散性特點給磁異常數據的解釋增加了困難。最小二乘法[11]是一種常見的擬合方法，最小二乘解是在假設測量點產生的隨機誤差為正態分布的前提下，采用最大似然估計法得到的一個最優估計解，該解使得測量誤差的平方和最小,因此可采用最小二乘法對離散等值點進行擬合，以求解得到最接近磁源體邊界輪廓曲線，進而利用擬合曲線及其參數對磁源體的水平分布范圍做出定量推斷。

由于水平放置的規則長方體和圓柱體等磁源體的特征邊界等值線的形狀非常接近橢圓，因此可對其進行最小二乘橢圓擬合。其基本過程為：首先假設橢圓參數，然后計算每個待擬合點到該橢圓的誤差距離的平方和，最后求出使該平方和最小的橢圓參數。

2 仿真分析

針對磁性組合體模型進行一組試驗。由長方體(目標1)、正方體(目標2)和圓柱體(目標3)構成的磁性組合體模型及其真實邊界位置如圖1所示，其幾何參數如表1所示。磁性組合體的磁化強度為30 A/m，磁傾角為56°，磁偏角為-16°。網格大小為121×121，水平方向上的采樣間隔為1 m。

圖1 磁性組合體模型及其真實邊界位置

目標中心坐標/m軸向長度/mxcyczcxayaza130-101620402021030142020203-25-518205020

為了消除斜磁化給邊界識別結果帶來的不利影響，首先對磁總場強度(Total Magnetic Intensity, TMI)數據進行化極處理。磁性組合體的化極TMI和特征邊界等值線的計算結果如圖2所示，可以看出化極TMI和特征邊界等值線與磁性組合體真實邊界位置的對應關系良好：λ、T1、T2、T3的零值等值線識別效果較好，基本能夠突出磁性組合體的邊界位置；θ1、θ2、θ3、θ4的邊界等值線識別效果較差，其中θ1、θ2、θ3的識別結果出現了假邊界，而θ4的識別結果由一系列封閉圓環組成，辨識度較差。

為了使仿真更加真實，向TMI數據中加入幅值為TMI最大值1%的高斯白噪聲，磁性組合體的化極TMI和特征邊界等值線的計算結果如圖3所示?？梢钥闯觯菏茉肼暩蓴_的影響，T2、T3、θ1、θ2、θ3、θ4識別結果的信噪比較低，辨識度較差；λ和T1受噪聲干擾的影響較小，其中λ的識別結果較準確，T1的識別結果較發散。

圖2 磁性組合體的化極TMI和特征邊界等值線

結合圖2、3可知：基于曲率磁梯度張量矩陣的最大特征值的歸一化改進方法受噪聲干擾的影響較小，其零值等值線更接近磁源體的真實邊界輪廓，識別精度更高。

為了進一步提高邊界識別精度，將圖3(b)中λ的零值等值線作為待擬合對象，對其進行濾波處理，然后對離散的邊界等值點進行最小二乘橢圓擬合，最終得到的磁性組合體邊界輪廓擬合結果如圖4所示，其擬合所得橢圓的幾何參數如表2所示。

由圖4可以看出：最小二乘擬合結果能夠直觀地描繪磁性組合體的邊界，且與其真實邊界輪廓的吻合度較高。對比表1、2可知：擬合所得橢圓的中心坐標、長軸和短軸與磁性組合體的實際邊界尺寸基本一致，且橢圓長軸的旋轉角為判斷磁源體的走向提供了近似參考。由此可以得出結論：最小二乘擬合結果能夠通過定量描繪磁源體的邊界特征提高邊界識別結果的精度。

圖4 磁性組合體的邊界輪廓擬合結果

目標中心坐標/mxcyc長軸2a/m短軸2b/m旋轉角θ/(°)130．21-10．3542．821．281．21210．0230．1120．9420．6151．683-25．74-5．2751．4718．251．3

3 實例驗證

為驗證上述方法的實際應用效果，在某2.1 m×2.1 m的測區內分別放置一個圓筒(目標1)和一個圓盤(目標2)進行探測試驗，如圖5所示。設觀測面的z坐標為0，則圓筒的中心坐標為(1.28,0.63,0.4)，母線和直徑長度分別為0.46、0.11 m；圓盤的中心坐標為(0.43,1.45,0.6)，母線和直徑長度分別為0.1、0.33 m；網格大小為22×22，水平方向上的采樣間隔為0.1 m；測區的地磁場傾角為56°，地磁場偏角為-16°。

對實測TMI數據進行化極處理，化極TMI和特征邊界等值線的計算結果如圖6所示?？梢钥闯觯?/p>

圖5 圓筒和圓盤的探測試驗

T2、T3、θ2、θ3的邊界等值線能夠大致描繪出磁源體的輪廓，但識別結果較發散；T1、θ1、θ4的識別結果更發散，已嚴重干擾真實邊界的辨別：相比之下，λ的識別效果最好，其零值等值線與實際模型真實邊界輪廓的吻合度較高。因此，可將λ的零值等值線作為邊界擬合對象。

為了提高識別精度，首先對圖6(b)中λ的零值等值線進行濾波處理，然后利用最小二乘法對其進行擬合，最終得到圓筒和圓盤的邊界擬合結果如圖7所示，擬合所得橢圓的幾何參數如表3所示。

圖6 實際模型的化極TMI和特征邊界等值線

圖7 圓筒和圓盤的邊界擬合結果

目標中心坐標/mxcyc長軸2a/m短軸2b/m旋轉角θ(°)11．250．60．550．44178．0820．431．410．610．4560．85

根據圖7和表3可知：擬合所得橢圓的中心坐標、長軸和短軸等參數與實際模型大體吻合，橢圓長軸的旋轉角為判斷實際模型的走向提供了大致參考，這在一定程度上提高了邊界識別結果的準確度。

4 結論

最小二乘擬合結果為判斷磁源體的中心位置、水平分布范圍和大致走向等提供了定量參考，提高了邊界識別結果的準確度，可用于靶場中地雷和未爆彈的型號判斷，海洋環境中潛艇、水雷和地下管道的輪廓識別以及油氣礦產的勘探范圍估計等。但受測量誤差等因素的影響，最小二乘擬合結果與實際磁源體模型的真實邊界仍存在一定誤差，下一步將圍繞磁源體的三維反演展開研究。

[1] NABIGHIAN M N.The analytic signal of two-dimensional magnetic bodies with polygonal cross-section:its properties and use for automated anomaly interpretation[J].Geophysics,1972,37(3):507-517.

[2] 王萬銀.位場總水平導數極值位置空間變化規律研究[J].地球物理學報,2010,53(9):2257-2270.

[3] WANG W,PAN Y,QIU Z.A new edge recognition technology based on the normalized vertical derivative of the total horizontal derivative for potential field data[J].Applied geophysics,2009,6(3):226-233.

[4] MILLER H G,SINGH V.Potential field T：a new concept for location of potential field sources[J].Journal of applied geophy-sics,1994,32(2/3):213-217.

[5] 馬國慶,黃大年,于平,等.改進的均衡濾波器在位場數據邊界識別中的應用[J].地球物理學報,2012,55(12):4288-4295.

[6] 張雙喜,張青杉,陳超,等.利用磁異常模量進行磁性體邊界檢測[J].地質與勘探,2015,51(6):1025-1032.

[7] WIJNS C,PEREZ C,KOWALCZYK P.Theta map:edge detection in magnetic data[J].Geophysics,2005,70(4):L39-L43.

[8] 邰振華.位場數據高精度處理方法的研究與應用[D].長春:吉林大學,2016.

[9] ORUC B,SERTCELIK I,KAFADAR?,et al.Structural interpretation of the Erzurum basin,eastern Turkey,using curvature gravity gradient tensor and gravity inversion of basement relief[J].Journal of applied geophysics,2013(88):105-113.

[10] ZHOU W,DU X,LI J.The limitation of curvature gravity gradient tensor for edge detection and a method for overcoming it[J].Journal of applied geophysics,2013(98):237-242.

[11] FITZGIBBON A,PILU M,FISHER R B.Direct least square fitting of ellipses[J].IEEE transactions on pattern analysis and machine intelligence,1999,21(5):476-480.