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

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

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

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

1 邊界識別理論

1.1 磁梯度張量

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

(1)

式中：G為磁梯度張量矩陣；Bi(i=x,y,z)為磁場矢量；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 傾斜角法

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

(6)

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

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

(7)

將增強(qiáng)型傾斜角定義為T3，其表達(dá)式為

(8)

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

1.4 Theta圖法

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

(9)

(10)

(11)

(12)

1.5 最小二乘法

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

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

2 仿真分析

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

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

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

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

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

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

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

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

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

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

目標(biāo)中心坐標(biāo)/mxcyc長軸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 實例驗證

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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