改進(jìn)的磁源體邊界識(shí)別方法

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

(陸軍工程大學(xué)石家莊校區(qū)車(chē)輛與電氣工程系， 河北 石家莊 050003)

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

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

1 邊界識(shí)別理論

1.1 磁梯度張量

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

(1)

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

1.2 曲率磁梯度張量

曲率磁梯度張量矩陣C為

(2)

求解其特征值矩陣，可得

(3)

(4)

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

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

(5)

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

1.3 傾斜角法

傾斜角是位場(chǎng)總水平導(dǎo)數(shù)和垂向?qū)?shù)比值的反正切函數(shù)，其零值等值線可突出場(chǎng)源體的邊界特征。張雙喜等[6]對(duì)傾斜角進(jìn)行了改進(jìn)，將改進(jìn)形式定義為T(mén)1，其表達(dá)式為

(6)

式中：U為磁標(biāo)量勢(shì)。

將式(6)拓展成張量形式，并定義為T(mén)2，其表達(dá)式為

(7)

將增強(qiáng)型傾斜角定義為T(mén)3，其表達(dá)式為

(8)

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

1.4 Theta圖法

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

(9)

(10)

(11)

(12)

1.5 最小二乘法

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

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

2 仿真分析

針對(duì)磁性組合體模型進(jìn)行一組試驗(yàn)。由長(zhǎng)方體(目標(biāo)1)、正方體(目標(biāo)2)和圓柱體(目標(biāo)3)構(gòu)成的磁性組合體模型及其真實(shí)邊界位置如圖1所示，其幾何參數(shù)如表1所示。磁性組合體的磁化強(qiáng)度為30 A/m，磁傾角為56°，磁偏角為-16°。網(wǎng)格大小為121×121，水平方向上的采樣間隔為1 m。

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

目標(biāo)中心坐標(biāo)/m軸向長(zhǎng)度/mxcyczcxayaza130-101620402021030142020203-25-518205020

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

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

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

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

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

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

圖4 磁性組合體的邊界輪廓擬合結(jié)果

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

3 實(shí)例驗(yàn)證

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

對(duì)實(shí)測(cè)TMI數(shù)據(jù)進(jìn)行化極處理，化極TMI和特征邊界等值線的計(jì)算結(jié)果如圖6所示?？梢钥闯觯?/p>

圖5 圓筒和圓盤(pán)的探測(cè)試驗(yàn)

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

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

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

圖7 圓筒和圓盤(pán)的邊界擬合結(jié)果

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

根據(jù)圖7和表3可知：擬合所得橢圓的中心坐標(biāo)、長(zhǎng)軸和短軸等參數(shù)與實(shí)際模型大體吻合，橢圓長(zhǎng)軸的旋轉(zhuǎn)角為判斷實(shí)際模型的走向提供了大致參考，這在一定程度上提高了邊界識(shí)別結(jié)果的準(zhǔn)確度。

4 結(jié)論

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

[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] 王萬(wàn)銀.位場(chǎng)總水平導(dǎo)數(shù)極值位置空間變化規(guī)律研究[J].地球物理學(xué)報(bào),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] 馬國(guó)慶,黃大年,于平,等.改進(jìn)的均衡濾波器在位場(chǎng)數(shù)據(jù)邊界識(shí)別中的應(yīng)用[J].地球物理學(xué)報(bào),2012,55(12):4288-4295.

[6] 張雙喜,張青杉,陳超,等.利用磁異常模量進(jìn)行磁性體邊界檢測(cè)[J].地質(zhì)與勘探,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] 邰振華.位場(chǎng)數(shù)據(jù)高精度處理方法的研究與應(yīng)用[D].長(zhǎng)春:吉林大學(xué),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.