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

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

0 引言

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

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

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

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

1 方法原理

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

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

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

(1)

式中：F[]為對數(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是磁勢異常，則觀測面上的磁異常場總強(qiáng)度矢量ΔT[13]可以表示為

(2)

(3)

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

(4)

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

(5)

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

(6)

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

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

(7)

1.2 頻率域三維等效源模型

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

(8)

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

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

(9)

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

(10)

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

(11)

根據(jù)(11)式計算得到第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)化理論

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

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

(13)

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

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

(14)

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

1.4 迭代計算法

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

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

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

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

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

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

2 數(shù)值仿真

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

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

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

圖5 測線上不同計算方法的磁梯度張量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 仿真計算中不同算法得到的參數(shù)及計算耗時Tab.2 Parameters and time consumptions obtained bydifferent methods in simulation calculation

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

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

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

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

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

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

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

3.2 不同方法的計算精度對比

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

圖11 測線上不同計算方法的磁梯度張量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é)論

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

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

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

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

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

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