基于頻率域三維等效源的磁梯度張量轉(zhuǎn)化方法

李金朋， 范紅波， 張英堂， 李志寧， 尹剛， 孫富成

(1.陸軍工程大學(xué)石家莊校區(qū) 車輛與電氣工程系， 河北 石家莊 050003；2.中國(guó)空氣動(dòng)力研究與發(fā)展中心 高速空氣動(dòng)力研究所，四川 綿陽(yáng) 621000； 3.63880部隊(duì)， 河南 洛陽(yáng) 471003)

0 引言

隨著磁異常檢測(cè)理論的不斷進(jìn)步和發(fā)展，磁梯度張量數(shù)據(jù)因?yàn)榫哂懈叻直媛实奶攸c(diǎn)而逐漸成為當(dāng)前的研究熱點(diǎn)[1-2]。地磁背景場(chǎng)的磁梯度較低，張量分量主要反映近地磁源的梯度數(shù)據(jù)，因此利用磁梯度張量分量數(shù)據(jù)能夠?qū)崿F(xiàn)磁源定位、二維邊界識(shí)別、深度反演以及三維姿態(tài)反演等計(jì)算[2-4]，并廣泛應(yīng)用于軍事偵察(未爆彈、潛艇和水雷等)和民用等領(lǐng)域。

等效源計(jì)算方法能夠利用等效源模型對(duì)磁源的觀測(cè)數(shù)據(jù)進(jìn)行重建，是一種應(yīng)用廣泛的數(shù)據(jù)處理方法。在傳統(tǒng)的等效源計(jì)算方法中，常常利用空間域的等效源計(jì)算方法構(gòu)造地下單層或多層等效源模型。Lyrio[5]利用等效源的計(jì)算方法對(duì)不規(guī)則數(shù)據(jù)進(jìn)行插值計(jì)算；Ali等[6]利用三維等效源計(jì)算方法，實(shí)現(xiàn)對(duì)觀測(cè)數(shù)據(jù)進(jìn)行延拓計(jì)算；陳濤等[7]基于三維等效源法對(duì)重力數(shù)據(jù)進(jìn)行等效源模擬，獲得重力梯度張量數(shù)據(jù)。在三維等效源計(jì)算方法中，在利用多層等效源進(jìn)行計(jì)算時(shí)，計(jì)算精度會(huì)提高，模擬場(chǎng)源與真實(shí)場(chǎng)更加接近，但計(jì)算過(guò)程類似于三維反演，提高了計(jì)算難度和計(jì)算成本，在處理大規(guī)模數(shù)據(jù)的過(guò)程中受到一定的限制。李端等[8]提出多層等效源的計(jì)算方法，并應(yīng)用于磁源的磁測(cè)數(shù)據(jù)轉(zhuǎn)換，但是如何確定多層等效源的參數(shù)設(shè)置，需要進(jìn)一步的探索。謝汝寬等[9]提出利用最小擬合差的單層等效源計(jì)算方法，該算法計(jì)算速度快，能夠直接對(duì)磁源的深度進(jìn)行估計(jì)。

頻率域的計(jì)算方法具有簡(jiǎn)單、快速的計(jì)算特點(diǎn)。Mickus等[10]基于快速傅里葉計(jì)算原理將重力垂直分量轉(zhuǎn)化為重力梯度張量；Jiang等[11]提出了利用余弦變換的方法對(duì)重力梯度張量數(shù)據(jù)進(jìn)行計(jì)算。上述頻率域的計(jì)算方法計(jì)算速度快，但是受噪聲影響較為明顯。

本文針對(duì)上述問(wèn)題，提出基于頻率域三維等效源的磁梯度張量轉(zhuǎn)換方法。首先，提出頻率域條件下的磁梯度張量正演計(jì)算方法；然后，構(gòu)建頻率域三維等效源模型，并利用迭代法進(jìn)行計(jì)算，以提高計(jì)算精度。最后，將獲得的地下等效源模型的磁化強(qiáng)度通過(guò)C范數(shù)正則化選擇方法進(jìn)行正演計(jì)算，以提高計(jì)算結(jié)果的穩(wěn)定性。

1 方法原理

1.1 頻率域正演基本理論與方程

假設(shè)觀測(cè)面為平面，將地下劃分為不同深度的水平長(zhǎng)方體層，設(shè)置長(zhǎng)方體層具有相同的尺寸大小。在地理坐標(biāo)系下，x軸正方向指向北向，y軸正方向指向東向，z軸正方向指向垂直向下。觀測(cè)面高度設(shè)置為zs(m). 假設(shè)某一長(zhǎng)方體平面的頂面深度為zh(m)，底面深度為zb(m). 長(zhǎng)方體層的磁異常Bz分量與磁化強(qiáng)度M分布的頻率域計(jì)算公式[12]為

F[Bz]={2πCmΘme|k|zs(e-|k|zh-e-|k|zb)}·
F[M],zs

(1)

式中：F[]為對(duì)數(shù)據(jù)進(jìn)行二維傅里葉變換；Cm=μ0/(4π)為常量，μ0為磁導(dǎo)率，μ0=4π×10-7H/m；Θm=(i·(xkx+yky)+z|k|)/|k|，x、y、z分別為實(shí)際磁化方向與x軸、y軸、z軸方向夾角的方向余弦，kx和ky分別為沿x軸、y軸方向的角頻率，

假設(shè)U是磁勢(shì)異常，則觀測(cè)面上的磁異常場(chǎng)總強(qiáng)度矢量ΔT[13]可以表示為

(2)

(3)

利用二維傅里葉變換以及傅里葉變換的空間域微分定理，可以得到磁總場(chǎng)異常與三分量之間的頻率域轉(zhuǎn)化關(guān)系[13]為

(4)

磁梯度張量是磁場(chǎng)矢量B=[Bx,By,Bz]沿著x軸、y軸、z軸3個(gè)正交方向上的空間變化率。則磁梯張量矩陣G可以寫成分別包含3個(gè)矢量元素的兩個(gè)矩陣的乘積:

(5)

式中：Bab(a,b=x,y,z)為磁梯度張量分量數(shù)據(jù)。對(duì)(5)式進(jìn)行二維傅里葉變換，得到

(6)

根據(jù)(4)式和(6)式，可以得到傅里葉變換后的磁梯度張量全張量數(shù)據(jù)：

F[G]=ξ·F[Bz]，

(7)

1.2 頻率域三維等效源模型

假設(shè)觀測(cè)面為平面，磁場(chǎng)H與磁源之間的關(guān)系[14]可以表示為

(8)

式中：R代表積分函數(shù)在三維空間內(nèi)定義；L=[(xV-xO)2+(yV-yO)2+(zV-zO)2]，(xO,yO,zO)為觀測(cè)點(diǎn)坐標(biāo)，(xV,yV,zV)為三維空間中體積元dv的坐標(biāo)。

由于將待測(cè)區(qū)域劃分為無(wú)數(shù)長(zhǎng)方體層，磁化分布僅在x軸和y軸方向上變化，因此將(8)式轉(zhuǎn)化為兩個(gè)部分：沿東- 北方向的水平部分S，以及沿z軸方向的垂直部分U[15]，二者的關(guān)系為

(9)

(9)式的頻率域計(jì)算公式[15]為

(10)

因此，磁異常Bz分量的頻率域計(jì)算公式[15]可以表示為

(11)

根據(jù)(11)式計(jì)算得到第c(c=1,2,3,…,Q，Q為劃分層數(shù))層地下水平層的磁化強(qiáng)度Mc分布[16]為

(12)

式中：h(kx,ky,zc)=(zce-|k|zc)s，s為常數(shù)。

根據(jù)(2)式，可以獲得不同水平層的磁化強(qiáng)度分布Mc. 根據(jù)獲得磁源的三維磁化強(qiáng)度分布，利用(7)式以及(11)式能夠得到磁源的磁梯度張量數(shù)據(jù)。

1.3 磁梯度張量數(shù)據(jù)正則化轉(zhuǎn)化理論

在對(duì)磁梯度張量數(shù)據(jù)進(jìn)行轉(zhuǎn)化過(guò)程中，觀測(cè)面的磁異常分量Bz數(shù)據(jù)中常常存在噪聲。由于張量轉(zhuǎn)換系數(shù)ξ的噪聲放大特性，導(dǎo)致獲得的張量數(shù)據(jù)計(jì)算結(jié)果不穩(wěn)定。為了抑制高頻噪聲，將正則化參數(shù)引入計(jì)算方程，構(gòu)造一個(gè)極小化正則化泛函[17-18]：

min{‖ξ-1·F[G]-F[Bz]‖2+λ‖F(xiàn)[G]‖2}，

(13)

式中：λ為正則化參數(shù)。上述極小化問(wèn)題的解為

F[G]=Aξ·F[Bz]，

(14)

在上述計(jì)算過(guò)程中，對(duì)正則化參數(shù)的選擇十分重要，數(shù)值過(guò)小將沒(méi)有濾波效果，數(shù)值較大將使計(jì)算的張量數(shù)據(jù)不準(zhǔn)確。本文采用C范數(shù)的計(jì)算方法對(duì)正則化參數(shù)進(jìn)行選擇[19]，計(jì)算不同正則化參數(shù)下相鄰兩個(gè)結(jié)果之間的C范數(shù)值，確定局部極小值對(duì)應(yīng)的正則化參數(shù)為最優(yōu)正則化參數(shù)。因此，通過(guò)選擇最佳的正則化參數(shù)，可以獲得穩(wěn)定且準(zhǔn)確的磁梯度張量數(shù)據(jù)。

1.4 迭代計(jì)算法

為了進(jìn)一步提高三維頻率域等效源計(jì)算方法的精度，本文利用迭代法確定磁源的各項(xiàng)參數(shù)。

步驟1對(duì)觀測(cè)面上的磁異常Bz分量進(jìn)行二維傅里葉變換。

步驟2利用(12)式依次計(jì)算每個(gè)水平層的二維磁化強(qiáng)度分布的頻譜數(shù)據(jù)，得到三維磁化強(qiáng)度分布數(shù)據(jù)Mc。

步驟3利用正演計(jì)算(1)式計(jì)算得到磁異常Bz，n分量，并計(jì)算觀測(cè)數(shù)據(jù)Bz與計(jì)算數(shù)據(jù)Bz，n的差值B′z,n，n表示第n次迭代計(jì)算，n=1，2，…，N，N為最大迭代次數(shù)。

步驟4計(jì)算B′z，n的均方根誤差，當(dāng)均方根誤差大于容差極限且迭代次數(shù)未達(dá)到最大迭代次數(shù)時(shí)，計(jì)算B′z，n數(shù)據(jù)的磁化強(qiáng)度三維反演結(jié)果M′c，對(duì)磁化強(qiáng)度進(jìn)行更新Mc+1=Mc+M′c，重復(fù)計(jì)算，直到均方根誤差小于容差極限，或者迭代次數(shù)達(dá)到最大迭代次數(shù)。

利用頻率域三維等效源方法對(duì)磁梯度張量數(shù)據(jù)進(jìn)行計(jì)算，需要確定地下劃分的網(wǎng)格數(shù)量及尺寸，由于頻率域計(jì)算速度快，因此采用與觀測(cè)點(diǎn)相同的劃分方式對(duì)地下網(wǎng)格的水平方向進(jìn)行劃分，在垂直面的分布深度上，利用迭代法思想設(shè)置不同的深度，對(duì)磁測(cè)數(shù)據(jù)反演精度進(jìn)行計(jì)算，將達(dá)到精度要求對(duì)應(yīng)的深度，作為地下三維網(wǎng)格的垂直深度。

2 數(shù)值仿真

2.1 正演模型及本文參數(shù)選擇

為了對(duì)本文方法進(jìn)行有效分析，構(gòu)建7個(gè)不同模型對(duì)本文方法進(jìn)行驗(yàn)證。將水平觀測(cè)面劃分為43×43=1 849個(gè)網(wǎng)格，網(wǎng)格間距為0.1 m. 模型的地下正演模型分布如圖1所示，各模型具有不同的形狀、位置、磁化方向以及磁化強(qiáng)度，具體物性參數(shù)如表1所示。地磁背景場(chǎng)的磁化方向?yàn)榇艃A角I為60°，磁偏角D為20°. 圖2為組合模型在觀測(cè)面上的磁異常Bz分量數(shù)據(jù)，圖2(a)為模型理論磁異常Bz分量，圖2(b)為加入信噪比為20 dB的高斯噪聲的磁異常Bz分量，其中虛線為觀測(cè)線，觀測(cè)線上不同模型磁異常Bz分量數(shù)據(jù)見(jiàn)圖2(c)所示，模型的理論磁梯度張量數(shù)據(jù)如圖3所示。

圖1 地下正演模型Fig.1 Underground forward model

表1 模型理論物性參數(shù)Tab.1 Physical and geometrical properties of models

圖2 模型產(chǎn)生的磁異常Bz分量數(shù)據(jù)Fig.2 Magnetic anomaly Bz components generated by the models

圖3 磁梯度張量理論數(shù)據(jù)Fig.3 Theoretical data of magnetic gradient tensor

圖4 磁梯度張量與Bz分量在不同迭代次數(shù) 下的均方根誤差Fig.4 RMSE values of magnetic gradient tensor and component Bz under different iteration times

2.2 不同方法的計(jì)算精度對(duì)比

確定最小擬合差單層等效源方法最佳深度時(shí)，設(shè)置地下網(wǎng)格最大深度為1 m，計(jì)算間隔為0.1 m. 利用圖2(b)測(cè)線上的磁異常Bz分量數(shù)據(jù)進(jìn)行計(jì)算，根據(jù)圖像可知，該測(cè)線上的磁測(cè)數(shù)據(jù)主要由模型6和模型7的磁場(chǎng)組成，而模型1、模型2、模型3、模型4、模型5磁異常數(shù)幾乎為0，進(jìn)而可以討論包含兩種不同深度條件下組合模型不同計(jì)算方法的計(jì)算準(zhǔn)確度。在進(jìn)行等效源計(jì)算時(shí)，將地磁背景場(chǎng)的磁化方向作為等效源的磁化方向，不同方法得到的磁梯度張量分量數(shù)據(jù)如圖5所示。根據(jù)圖5可以看出，利用最小擬合差單層等效源方法計(jì)算得到的磁梯度張量數(shù)據(jù)的誤差較為明顯，頻率域三維等效源方法與本文算法獲得的計(jì)算結(jié)果與理論計(jì)算結(jié)果基本保持一致，但頻率域三維等效源計(jì)算方法獲得的計(jì)算結(jié)果受噪聲影響較為明顯，存在一定的波動(dòng)。為了進(jìn)一步對(duì)本文方法的抗噪性能進(jìn)行說(shuō)明，將本文方法以及頻率域三維等效源方法獲得的磁梯度張量計(jì)算輪廓進(jìn)行比較，計(jì)算結(jié)果如圖6所示。由圖6可以看出，本文算法通過(guò)正則化計(jì)算，獲得的計(jì)算結(jié)果與理論磁梯度張量較為匹配。

圖5 測(cè)線上不同計(jì)算方法的磁梯度張量Fig.5 Magnetic gradient tensors of different calculation methods over the observation line

圖6 磁梯度張量數(shù)據(jù)輪廓(單位：nT/m)Fig.6 Contours of magnetic gradient tensor data (unit: nT/m)

圖7 磁梯度張量分量的互相關(guān)系數(shù)與均方根誤差Fig.7 Cross-correlation coefficients and RMSE values of magnetic gradient tensor components

表2 仿真計(jì)算中不同算法得到的參數(shù)及計(jì)算耗時(shí)Tab.2 Parameters and time consumptions obtained bydifferent methods in simulation calculation

3 試驗(yàn)驗(yàn)證

3.1 待測(cè)模型及本文參數(shù)選擇

圖8 測(cè)量裝置及實(shí)測(cè)磁異常Bz分量數(shù)據(jù)Fig.8 Measurement device and measured magnetic anomaly component Bz

對(duì)地下未爆彈進(jìn)行一組實(shí)測(cè)試驗(yàn)。試驗(yàn)地點(diǎn)位于石家莊某地，將水平觀測(cè)面劃分為19×22=418個(gè)網(wǎng)格，網(wǎng)格間距為0.1 m. 地磁背景場(chǎng)的磁化方向?yàn)榇艃A角I為56°，磁偏角D為-18°. 試驗(yàn)所用設(shè)備如圖8(a)所示，利用英國(guó)Bartington公司三軸磁通門傳感器搭建十字型磁梯度張量探頭，測(cè)試系統(tǒng)主要包括十字型探頭和數(shù)字采集模塊及軟件操作終端。試驗(yàn)過(guò)程中探頭固定在無(wú)磁試驗(yàn)臺(tái)架上，利用單點(diǎn)測(cè)量方式，滑動(dòng)張量探頭在無(wú)磁滑桿上對(duì)待測(cè)區(qū)域的每一個(gè)觀測(cè)點(diǎn)進(jìn)行測(cè)量，無(wú)磁試驗(yàn)臺(tái)架的4個(gè)固定點(diǎn)位置為可調(diào)整角度的三腳架結(jié)構(gòu)，利用水平儀將無(wú)磁試驗(yàn)臺(tái)架調(diào)整在同一水平面上，保證探頭在測(cè)試過(guò)程中能夠保持水平。十字型探頭的基線距離為40 cm，在進(jìn)行實(shí)際測(cè)量試驗(yàn)時(shí)，將單傳感器測(cè)量得到的Bz數(shù)據(jù)與背景磁場(chǎng)的Bz數(shù)據(jù)進(jìn)行做差，作為測(cè)區(qū)內(nèi)該點(diǎn)的磁異常Bz數(shù)據(jù)，測(cè)區(qū)內(nèi)實(shí)際測(cè)量數(shù)據(jù)如圖8(b)所示。同時(shí)，該點(diǎn)的張量數(shù)據(jù)利用差分原理可以直接求解，實(shí)際測(cè)量得到的磁梯度張量分量數(shù)據(jù)如圖9所示。為了減少傳感器探頭的噪聲、不一致性、零漂以及非對(duì)準(zhǔn)誤差對(duì)測(cè)試數(shù)據(jù)精度帶來(lái)的影響，該系統(tǒng)利用非線性校正方法進(jìn)行了非線性校正[20]，因此將獲得的觀測(cè)數(shù)據(jù)作為實(shí)測(cè)試驗(yàn)的理論數(shù)據(jù)。

圖9 磁梯度張量實(shí)測(cè)數(shù)據(jù)Fig.9 Measured magnetic gradient tensor data

為了證明本文算法的有效性，向?qū)崪y(cè)磁異常Bz分量數(shù)據(jù)中加入信噪比為40 dB的高斯噪聲，獲得的測(cè)量數(shù)據(jù)如圖8(c)所示。通過(guò)迭代計(jì)算獲得最大劃分深度為8層，在最大劃分層數(shù)下的迭代次數(shù)共16次。在迭代終止前，隨迭代次數(shù)增加獲得的張量和磁異常Bz分量的均方根誤差計(jì)算結(jié)果如圖10所示，測(cè)線上不同方法的磁梯度張量分量數(shù)據(jù)如圖11所示，本文算法以及頻率域三維等效源計(jì)算方法獲得的磁梯度張量計(jì)算輪廓如圖12所示。在試驗(yàn)計(jì)算中將實(shí)測(cè)磁梯度張量數(shù)據(jù)作為實(shí)測(cè)數(shù)據(jù)，根據(jù)計(jì)算結(jié)果可以看出，本文算法獲得的計(jì)算結(jié)果能夠有效降低噪聲的干擾，與實(shí)測(cè)數(shù)據(jù)基本保持一致，計(jì)算效果最好。

圖10 磁梯度張量與Bz分量在不同迭代 次數(shù)下的均方根誤差Fig.10 The e values of magnetic gradient tensor and component Bz under different iteration times

3.2 不同方法的計(jì)算精度對(duì)比

本文算法與頻率域三維等效源算法、最小擬合差單層等效源算法以及FFT計(jì)算方法進(jìn)行比較，磁梯度張量分量的互相關(guān)系數(shù)與均方根誤差如圖13所示。由圖13可見(jiàn)：在噪聲等因素影響下，F(xiàn)FT計(jì)算方法受噪聲影響較為明顯；最小擬合差單層等效源算法計(jì)算得到的磁梯度張量數(shù)據(jù)相關(guān)性較差；本文算法的計(jì)算精度最高，通過(guò)正則化計(jì)算能夠獲得與實(shí)測(cè)磁梯度張量較為匹配的計(jì)算結(jié)果；FFT計(jì)算方法受噪聲影響計(jì)算精度較差，與數(shù)值仿真結(jié)論相一致；在進(jìn)行單層等效源計(jì)算時(shí)，需要已知先驗(yàn)信息，且計(jì)算結(jié)果存在一定的漏磁現(xiàn)象。上述方法在計(jì)算過(guò)程得到的參數(shù)以及計(jì)算耗時(shí)如表3所示。由表3可見(jiàn)，最終耗時(shí)結(jié)果與數(shù)值仿真結(jié)果類似，而基于最優(yōu)擬合差的單層等效源算法用時(shí)最多，本文算法的計(jì)算時(shí)間居中，F(xiàn)FT算法耗時(shí)最短。通過(guò)上述分析可知，本文算法能夠有效對(duì)未爆彈的張量數(shù)據(jù)進(jìn)行計(jì)算，為下一步的定位與識(shí)別[4]提供有效的數(shù)據(jù)支持。

圖11 測(cè)線上不同計(jì)算方法的磁梯度張量Fig.11 Magnetic gradient tensors for different calculation methods over the observation line

圖12 磁梯度張量數(shù)據(jù)輪廓(單位：nT/m)Fig.12 Contours of magnetic gradient tensor data (unit: nT/m)

圖13 磁梯度張量分量的互相關(guān)系數(shù)與均方根誤差Fig.13 Cross correlation and RMSE values of magnetic gradient tensor components

4 結(jié)論

對(duì)磁性異常體進(jìn)行正演計(jì)算時(shí)，當(dāng)噪聲較大時(shí)會(huì)對(duì)張量計(jì)算結(jié)果產(chǎn)生一定的影響。針對(duì)上述問(wèn)題，本文提出了基于頻率域三維等效源的磁梯度張量正演理論，構(gòu)建頻率域三維等效源并利用迭代法自適應(yīng)的計(jì)算網(wǎng)格劃分深度與精度，通過(guò)C范數(shù)正則化選擇方法對(duì)地下等效源模型的磁梯度張量數(shù)據(jù)進(jìn)行正演計(jì)算。通過(guò)仿真和試驗(yàn)證明了本文方法具有如下優(yōu)勢(shì)：

表3 試驗(yàn)計(jì)算中不同算法得到的參數(shù)及計(jì)算耗時(shí)Tab.3 Parameters and time consumptions obtained by differentmethods in simulation calculation

1) 構(gòu)建基于磁異常Bz分量數(shù)據(jù)的頻率域三維等效源計(jì)算模型，提高了計(jì)算結(jié)果的穩(wěn)定性，并減少了計(jì)算成本。

2) 利用迭代法對(duì)地下水平層深度以及磁化強(qiáng)度反演結(jié)果進(jìn)行自適應(yīng)求解計(jì)算，實(shí)現(xiàn)了計(jì)算過(guò)程的自動(dòng)化計(jì)算。

3) 提出了基于C范數(shù)正則化計(jì)算方法的磁梯度張量正演計(jì)算方法，能夠在噪聲干擾下穩(wěn)定的獲得磁梯度張量轉(zhuǎn)換結(jié)果。

磁性異常體磁梯度張量數(shù)據(jù)的正演是對(duì)磁性異常體進(jìn)行數(shù)據(jù)解釋的基礎(chǔ)，通過(guò)仿真和試驗(yàn)證明了本文方法能夠在噪聲條件下實(shí)現(xiàn)對(duì)磁性異常體的張量數(shù)據(jù)正演，為未爆彈藥等隱蔽武器的探測(cè)與數(shù)據(jù)解釋提供了有力的技術(shù)手段。