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

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

?

基于磁梯度張量的地下小目標相關(guān)成像方法

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

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

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

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

0引言

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

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

1相關(guān)成像法基本原理

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

(1)

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

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

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

根據(jù)Schwarz不等式可知:

(3)

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

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

2.1磁梯度張量要素

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

(4)

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

根據(jù)磁總場異常三維相關(guān)成像原理，同理可推導(dǎo)出各磁梯度張量分量Cαβ的歸一化互相關(guān)公式：

αβ=x,y,z

(5)

地下小目標三維相關(guān)成像的步驟為：

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

2)根據(jù)公式(5)求網(wǎng)格節(jié)點的異常與實測異常的歸一化互相關(guān)；

3)將所得結(jié)果置于網(wǎng)格節(jié)點；

4)從淺層到深層逐層計算各網(wǎng)格節(jié)點的相關(guān)系數(shù)，得到相關(guān)系數(shù)數(shù)據(jù)庫；

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

2.2由相關(guān)系數(shù)向?qū)嶋H物性參數(shù)轉(zhuǎn)化

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

圖1　相關(guān)系數(shù)到實際物性參數(shù)轉(zhuǎn)化流程Fig.1　The conversion process of correlationto actual physical parameters

3仿真驗證

3.1鐵磁目標正演計算

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

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

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

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

3.2模型測試

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

圖4　未爆彈磁梯度數(shù)據(jù)相關(guān)成像結(jié)果Fig.4　The correlation imaging results of magnetic gradient tensor

(a)磁梯度張量加權(quán)后相關(guān)成像結(jié)果圖

(b)加入先驗信息后相關(guān)成像結(jié)果圖

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

Fig.5inversionresultssliceswithconstraint

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

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

4結(jié)論

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

參考文獻:

[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]許令周，關(guān)繼騰，房文靜.高次導(dǎo)數(shù)的概率成像原理[J].青島大學學報，2003,16(4):32-36.

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

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

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

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

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

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

[12]孟慧.磁梯度張量正演、延拓、數(shù)據(jù)解釋方法研究[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