基于磁梯度張量的地下小目標相關成像方法

李金朋，張英堂，范紅波，李志寧，張　帥

(1.解放軍軍械工程學院七系， 河北 石家莊 050003;2.解放軍72465部隊，山東 濟南 250022)

?

基于磁梯度張量的地下小目標相關成像方法

李金朋1，張英堂1，范紅波1，李志寧1，張帥2

(1.解放軍軍械工程學院七系， 河北 石家莊 050003;2.解放軍72465部隊，山東 濟南 250022)

摘要:針對小目標鐵磁物質反演方法研究中，磁總場探測包含信息不完整并且傳統反演方法無法通過先驗信息與目標建立聯系的問題，提出了基于磁梯度張量的地下小目標相關成像方法。該方法首先利用有限元方法對鐵磁目標在測量面的磁梯度張量場進行正演計算以作為實測數據；然后將地下待成像空間劃分為三維規則網格并將每個網格節點等效為磁偶極子，計算每個磁偶極子在測量面產生的磁梯度張量數據與實測數據的互相關系數，作為磁偶極子在待成像空間的密度分布圖；最后利用迭代計算將密度分布圖轉化為實際的物性參數分布圖。仿真驗證表明：該方法不僅能夠獲得地下小目標的分布情況，還能得到目標位于不同位置的實際物性參數值。

關鍵詞:磁梯度張量；地下小目標；相關成像;迭代計算

0引言

地下小目標磁性體反演技術是通過地下磁性目標(未爆彈藥、掩埋水雷等)所引起的磁異常值解算出其在地下的分布信息，其特點是定位速度快、經濟性比較高，是應用比較廣泛的一種目標探測方法。三維相關成像技術是一種通過計算實測磁異常與地下待成像磁偶極子掃描函數之間的歸一化互相關系數來獲得地下目標分布的方法。1997年由Patella首次提出將概率成像方法應用于自然電位異常的解釋[1-2]。Mauriello等將概率成像法推廣到大地電磁和重力領域[3-4]。國內許令周、王緒本、毛立峰等利用概率成像方法對自然電場和大地電磁數據進行了研究[5-7]，取得了一定進展。中國地質大學郭良輝等根據概率成像的概念提出重力數據的三維相關成像法并將其推廣到磁場數據的相關成像，提出了磁總場異常相關成像法[8]；石磊等利用磁總場異常垂直梯度數據來解釋異常源的位置，提高了成像分辨率，一定程度上改善了復雜異常的成像效果。

然而，在以上的目標反演研究中，由于磁總場探測包含信息不完整并且相關成像法所得到的結果為互相關系數分布圖，無法通過先驗信息與之建立聯系[9]，最終影響了相關成像法在實際中的應用。針對此問題，本文提出基于磁梯度張量的地下小目標相關成像方法。

1相關成像法基本原理

笛卡爾坐標的(x,y)平面位于基準面上，z取向下為正，假設測區地下任意一點坐標為(x,y,z),體積為v的均勻磁化球體p的磁化強度為Jq，磁矩M=Jv，由于均勻磁化的球體外部磁場與位于球心的偶極子磁場等效[10]，因此我們把均勻磁化的球體視為磁偶極子。磁偶極子在測區上任意點(x,y,z)的磁總場異常B(x,y,z)可以表示為：

(1)

式(1)中，u0為真空磁導率。(xi,yi,zi)為測區第i個點的坐標，

借鑒Abdelrahman等的相鄰最小二乘剩余異常的歸一化互相關公式，我們定義實測磁異常與磁偶極子p在測區上的理論磁異常歸一化互相關公式為：

式(2)中，(xi,yi,zi)為測區第i個觀測點坐標標，B(xi,yi,zi)為該觀測點實測磁異常，Bq(xi,yi,zi)為地下第q個點在該觀測點的磁異常，N為測區觀測點總數。

根據Schwarz不等式可知:

(3)

因此可知，公式(3)的相關系數C的的取值范圍為-1≤C≤+1。其值的意義是實測磁異常由偶極子p產生的可能性。C絕對值越高，說明q點存在偶極子的可能性越大。

2基于磁梯度張量的地下小目標相關成像方法

2.1磁梯度張量要素

與傳統的磁總場異常相比，磁梯度張量能提供更加豐富的異常信息，能更好地反映地下磁性體的細節特征。磁梯度張量可以由一個二階矩陣表示為：

(4)

式(4)中,G為磁梯度張量；Bαβ(α,β=x,y,z)為磁梯度張量的分量，共有9個；根據場論可知:Bxy=Byx,Bxz=Bzx,Byz=Bzy,Bxx+Byy+Bzz=0。因此，磁梯度張量的9個分量只有5個是獨立的。

根據磁總場異常三維相關成像原理，同理可推導出各磁梯度張量分量Cαβ的歸一化互相關公式：

αβ=x,y,z

(5)

地下小目標三維相關成像的步驟為：

1)將地下待成像空間劃分為三維均勻網格;

2)根據公式(5)求網格節點的異常與實測異常的歸一化互相關；

3)將所得結果置于網格節點；

4)從淺層到深層逐層計算各網格節點的相關系數，得到相關系數數據庫；

5)對數據進行可視化，獲得待測小目標鐵磁物質在地下的分布信息。

2.2由相關系數向實際物性參數轉化

通過相關成像法計算得到的是等效的物性參數。然而，在進行反演計算時，我們更加希望獲得真實的物性參數的分布。針對這一問題，根據相關成像結果與實際物性參數之間的對應關系，在相關成像結果的基礎上反演物性參數。具體做法是：假設待成像區剖分后的某個節點磁偶極子磁性參數為M，通過正演計算獲得位于觀測面上的磁梯度張量數據Data1，并Data1與實測數據進行相關成像得到相關系數值C1；然后，用一個小的磁性參數乘C1，此時待測鐵磁物質的磁性參數更新為M+C1·ΔM1；再次對該節點磁偶極子進行正演計算，獲得正演數據Data2，對Data1-Data2的差值與實測數據進行相關成像獲得相關系數C2。這個過程便可以理解為是反演過程中磁性參數的擾動量，而M2=M1+C2·ΔM為更新后的磁性參數值。具體流程見圖1。

圖1　相關系數到實際物性參數轉化流程Fig.1　The conversion process of correlationto actual physical parameters

3仿真驗證

3.1鐵磁目標正演計算

本文采用COMSOL軟件模擬地下未爆彈模型。COMSOL是一款通用的有限元建模軟件[11-12]，利用COMSOL軟件可以對未爆彈進行幾何建模和網格劃分。所建未爆彈模型頭部為尖拱弧形，彈翼為梯形板狀翼，十字布局，沿南北走向，其中心埋深距離地面3.5 m。對未爆彈模型自動網格剖分后的四面體單元數為5 333，平均單元質量為0.685 6。根據場論關系設定未爆彈的求解域屬性和邊界條件。在外界無電流環境下，▽×H=0，本文中設定未爆彈模型的相對磁導率為3.5，背景區域的相對磁導率為1。圖2展示了未爆彈在垂直磁化作用下的磁感應強度切片圖、等勢面圖以及磁化的方向。H為背景磁場，設定H=50 000nT。磁通密度和磁場強度的關系為B=u0(H+M)。

利用COMSOL對未爆彈模型的正演計算，以獲得未爆彈位于觀測面上的磁梯度張量數據。圖3為對未爆彈進行正演計算后得到的位于觀測面上的磁梯度張量數據分布圖。

圖2　未爆彈正演模型Fig.2　Forward model of UXO

圖3　未爆彈模型的磁梯度張量數據分布圖Fig.3　Magnetic gradient tensor model of UXO

3.2模型測試

為驗證本文反演方法的有效性，利用仿真所得磁梯度張量數據對未爆彈模型進行測試。觀測面位于距地面0.1m處，反演網格間距為0.5m×0.5m×0.2m。分別對未爆彈的磁梯度張量異常Bxz，Byz，Bzz進行相關成像,圖4為Z=3.5m時磁梯度張量相關成像結果圖，圖4(a)、(c)、(e)為相關成像結果圖。圖4(b)、(d)、(f)為迭代8次之后的實際物性參數成像的結果。

圖4　未爆彈磁梯度數據相關成像結果Fig.4　The correlation imaging results of magnetic gradient tensor

(a)磁梯度張量加權后相關成像結果圖

(b)加入先驗信息后相關成像結果圖

圖5加入約束后反演切片圖

Fig.5inversionresultssliceswithconstraint

通過圖4測試的結果可以看出：未進行迭代時，在深度為3.5m處相關成像的結果出現了正的成像值，基本能夠表征未爆彈的存在，但是，異常范圍偏大，成像的聚焦效果較差；采用迭代計算的方法，將相關系數轉換為實際的物性參數值，同時改善了聚焦效果。通過對比兩個結果可以看出：本文方法的反演結果首先將相關系數轉換為真正的物性參數值并獲得實際物性參數的分布圖；其次，改善了目標的聚焦效果，能更好的實現未爆彈的定位。

圖5(a)為對不同張量信息進行加權所得到的成像結果，(b)為加入已知先驗信息所得的成像結果。由于本文采用迭代計算的方法能夠引入已知的地質信息作為約束。因此，可以對每一個反演數值設置解區間[Mmin,Mmax]。對于本例而言,對不同張量數據進行加權并對相對磁導率設置統一解區間為[3.2,3.8]。對比圖5(a)、(b)可以看出：在反演計算中加入已知信息可以縮小解的范圍，減小解的非唯一性，提高反演的分辨率；此外，通過迭代法計算所獲得的物性參數值也更加接近真值。

4結論

本文提出了基于磁梯度張量地下小目標相關成像方法。該方法首先利用有限元方法對待測未爆彈模型在測量面的磁梯度張量場進行正演計算以作為實測數據；然后將地下待成像空間劃分為三維規則網格并將每個網格節點等效為磁偶極子，計算每個磁偶極子在測量面產生的磁梯度張量數據與實測數據的互相關系數，作為磁偶極子在待成像空間的密度分布圖；最后利用迭代計算將密度分布圖轉化為實際的物性參數分布圖。仿真驗證表明：該方法能更好的提高未爆彈的反演分辨率、改善聚焦效果，而且獲得的物性參數值更加接近真實值。不足之處在于只是通過實驗驗證了該方法的有效性，在工程實踐中的應用還有待進一步的補充和改進。

參考文獻:

[1]ReneRM.Gravityinversionusingopen,reject,and“shapetotheanomaly”fillcriteria[J].Geophysics, 1986,51(4):988-994.

[2]PatellaD.Introductiontogroundsurfaceself-potentialtomography[J].GeophysicalProspecting, 1997,45(4):653-681.

[3]MaurielloP，PatellaD.PrincipleofProbabilitytomgraphyfornatural-sourceelectromagneticinductionfields[J].Geophysics,1999,64(5):1404-1417.

[4]MaurielloP，PatellaD.Localizationofmaximum-depthgravityanomalysourcebyadistributionofe-quivalentpointmasses[J].Geophysics.2001,66(5):1431-1437.

[5]許令周，關繼騰，房文靜.高次導數的概率成像原理[J].青島大學學報，2003,16(4):32-36.

[6]王緒本，毛立峰，高永才.電磁導數場概率成像方法研究[J].成都理工大學學報(自然科學版)2004,31(6):679-684.

[7]王緒本，毛立峰，高永才.大地電磁概率成像效果評價[J].地球物理學報，2005,48(2):429-433.

[8]郭良輝，孟小紅，石磊.磁異常ΔT三維相關成像[J].地球物理學報，2010，53(2):435-441.

[9]孟小紅，劉國峰，陳召曦，等.基于剩余異常相關成像的重磁物性反演方法[J].地球物理學報， 2012，35(1):304-309.

[10]吳招才.磁力梯度張量技術及其應用研究[D].武漢：中國地質大學地球物理與空間信息學院，2008.

[11]ButlerSL.ForwardmodelingofappliedgeophysicsmethodsusingComsolandcomparisonwithanalyt-icalandlaboratoryanalogmodels[J].Computer&Geosciences,2012,42:168-176.

[12]孟慧.磁梯度張量正演、延拓、數據解釋方法研究[D].吉林：吉林大學儀器科學與電器工程學院,2012.

*收稿日期:2015-11-28

作者簡介:李金朋(1991—)，男，吉林長春人，碩士研究生，研究方向：目標探測與識別。E-mail:18626648671@163.com。

中圖分類號:P631

文獻標志碼:A

文章編號：1008-1194(2016)03-0075-04

CorrelationImagingMethodofSmallSubsurfaceTargetBasedonMagneticGradientTensor

LIJinpeng1,ZHANGYingtang1,FANHongbo1,LIZhining1,ZHANGShuai2

(OrdnanceEngineeringCollege,Shijiazhuang050003,China)

Abstract:The total magnetic field detection anomaly include a single exception message and can not establish contact with the target through prior information in the inversion study of small ferromagnetic material target. A correlation imaging method based on magnetic gradient tensor data was put forward in this paper. Firstly, forward modeling data was obtained by using finite element method as measured magnetic gradient tensor data of ferromagnetic target. Then the subsurface space was divided into a 3D regular grid and each grid node was equivalent to a magnetic dipole. The cross correlation was calculated between the magnetic gradient tensor at the measurement surface and the observed data. The resultant correlation coefficients was used to describe equivalent magnetic dipole distribution in a probability sense. Finally, an iterative method was proposed to finish the transformation from probability value to actual physical parameter distribution. Theoretical analysis showed that this method could not only obtain the distribution of the small subsurface target, but also get the actual material parameter at different location.

Key words:magnetic gradient tensor; small subsurface target; correlation imaging; iterative calculation