基于波數(shù)域上Fourier變換的磁異常信號處理

鄭建擁,范紅波*,尹 剛,李志寧

(1.陸軍工程大學(xué)石家莊校區(qū),石家莊 050003;2.中國空氣動力研究與發(fā)展中心高速所,四川 綿陽 621000)

在地球物理學(xué)領(lǐng)域里,Fourier變換是一種常用的數(shù)據(jù)處理方式,在不同的研究方向上,都有很深的研究潛力。如Nabighian(1972,1984),Nabighian and Hansen(2001),Nelson(1986)[1-3]利用Fourier變換在磁場的水平和垂直分量之間進行轉(zhuǎn)換;Nelson(1988)[4-5]計算了總磁場異常的磁梯度張量;Lourenco和Morrison(1973)[6]描述了磁總場模量(TMI)和磁場矢量分量之間轉(zhuǎn)換的數(shù)學(xué)推導(dǎo);Mohan,Sundararajan和Seshagiri Rao(1982)[7]利用Fourier變換建立地質(zhì)框架;Sundararajan和Srinivas(1996)[8-9]提出內(nèi)稟磁標勢方法以及在地震上的應(yīng)用。(Hood,1980;Overton,1980)[10]也做過關(guān)于地球物理應(yīng)用的磁梯度張量數(shù)據(jù)的分析工作。目前在磁異常正演領(lǐng)域,在數(shù)據(jù)的時間域和頻率域分析上都有成熟的推導(dǎo)理論,但這些理論也僅限于磁異常數(shù)據(jù)正演,實際應(yīng)用上并無巨大突破,也尚未涉及波數(shù)域的推演。而且在磁異常信號預(yù)處理方面,并沒有成熟的方法或者研究的思路、方向。

傳統(tǒng)的地磁勘測,是通過標量磁力儀采集地磁總場強度的值,即磁總場模量(TMI),所含信息量比較少,且信號質(zhì)量也不高。而近些年磁通門傳感器和超導(dǎo)量子干涉器件得到一系列新發(fā)展,可以直接測量磁場矢量及其空間變化率;全張量磁力梯度儀系統(tǒng)可以測量所有線性獨立梯度張量分量,也即磁場分量的一階導(dǎo)數(shù)。因此所得到的磁異常信息的質(zhì)量和數(shù)量得到很大改善,精確度明顯提高。而且磁性目標產(chǎn)生的總場、矢量場、總場梯度和梯度張量場在空間域中的計算方法已經(jīng)非常成熟。但以上途徑的數(shù)據(jù)處理過程都是根據(jù)不同儀器的測量方式,利用簡單的空間上的差分或求和來完成TMI、磁場矢量和梯度張量之間的變換,實際應(yīng)用上存在很大的誤差。因此,有必要運用數(shù)學(xué)方法推導(dǎo)TMI、磁梯度張量和磁場矢量等之間的解析關(guān)系,系統(tǒng)地進行各量之間的轉(zhuǎn)換?；诖?在Nabighian等人所做的數(shù)學(xué)框架的基礎(chǔ)上,提出了利用波數(shù)域上的Fourier變換,推導(dǎo)TMI、磁場矢量和磁梯度張量相互之間關(guān)系的方法。并在實驗數(shù)據(jù)上填充反距離加權(quán)插值算法,得到了噪聲更小、更清晰、分辨率更高的磁異常數(shù)據(jù),有利于之后的磁異常信號的解釋和應(yīng)用。

1 空間域的推導(dǎo)

在空間域數(shù)據(jù)中做頻譜分析,即可得到波數(shù)域數(shù)據(jù)。因此需要先得原始空間域數(shù)據(jù)。式(1)為傳統(tǒng)空間域上磁標量φ、磁場矢量三分量(Bx、By、Bz)和磁梯度張量矩陣G三者的關(guān)系公式,而式(2)為從磁場矢量推導(dǎo)磁總場模量(TMI)的求解公式(公式中表示為Bt)[11-14]

(1)

(2)

圖1所示的4磁通門十字形陣列探測儀是較為普遍的地磁測量設(shè)備,除此外還有長方形陣列、四面體陣列等根據(jù)不同需要和算法設(shè)計的儀器。除磁通門外,還有超導(dǎo)量子干涉儀等精度更高的地磁測量儀器[15-20]。

圖1 實驗測量設(shè)備

磁通門直接測量的僅為磁場矢量三分量(Bx、By、Bz),可以根據(jù)式(2)求解出磁總場模量(TMI)。而磁梯度張量的算法,是根據(jù)磁通門擺放位置的不同而差分計算得來的。以圖1的陣列為例,磁梯度張量矩陣G的計算公式如下[11]:

(3)

Bij(i=1,2,3,4;j=x,y,z)表示第i個傳感器在j方向上的磁場矢量;d為基線距離。

可知,空間域上實際磁異常數(shù)據(jù)的測量和磁梯度張量數(shù)據(jù)的計算受儀器設(shè)備的影響比較大,推導(dǎo)計算過程本身存在一定的誤差,如基線的選取、4個傳感器的對稱對中,以及數(shù)據(jù)計算方法等因素不可避免的產(chǎn)生了誤差。

2 波數(shù)域的推導(dǎo)

波數(shù)域方法首先對探測得到空間域數(shù)據(jù)進行頻譜分析(一般為傅里葉變換)后,在波數(shù)域上推導(dǎo)磁異常數(shù)據(jù)各種量的解析關(guān)系式,然后利用逆傅里葉變換得到磁異常數(shù)據(jù)。相對于直接從空間域得來的數(shù)據(jù),其優(yōu)點是異常頻譜表達式較為簡潔緊湊,計算效率更高。

磁場矢量是磁標量φ=φ(x,y,z)的負梯度,而磁場梯度張量元素是磁標量φ的二次導(dǎo)數(shù)[6]:

(4)

式(5)是空間域磁標量的傅里葉變換公式:

(5)

二維表示磁標量的逆傅里葉變換:

在這里,把這兩種變換縮寫為:

φ(p,q)=FT{φ(x,y)}
φ(x,y)=FT-1{φ(p,q)}

(6)

式中:p和q為水平波數(shù)。

二維傅立葉變換的使用要求將測量的數(shù)據(jù)采樣到規(guī)則的空間平面。因此,積分平面被限制為垂直坐標為常數(shù)的某一平面。

對于平面的垂直梯度張量Gkz[3],計算方法見式(7):

(7)

這是在波數(shù)域中,由磁場矢量推導(dǎo)垂直方向的梯度張量Gkz的公式。其傅立葉變換展開如式(8):

(8)

(9)

聯(lián)立式(7)～式(9)可以得到波數(shù)域的垂直梯度張量公式。以上這些關(guān)系可用于磁梯度垂直張量Gkz和磁矢量三分量Bk之間的變換。

為了獲得水平張量公式,有必要考慮波數(shù)域中的水平導(dǎo)數(shù):

FT{Gkx}=ipFT{Bk}

(10)

FT{Gky}=iqFT{Bk}

(11)

(12)

把式(7)、式(10)、式(11)代入式(12),可得波數(shù)域上的磁梯度全張量公式[3]:

(13)

聯(lián)立式(8)、式(13)可以在波數(shù)域上建立起磁場矢量三分量和所有磁梯度張量之間的解析推導(dǎo)。

還可以把式(13)重新定義為下式[4-5]:

FT{Bx}=FT{-?φ/?x}=-ipFT{φ}
FT{By}=FT{-?φ/?y}=-iqFT{φ}
FT{Bz}=FT{-?φ/?z}=-rFT{φ}

(14)

(15)

(16)

將此式整理為FT{φ}的形式,然后代入到式(14),可以得到TMI和磁場矢量三分量在波數(shù)域上的關(guān)系:

(17)

至此,完成在波數(shù)域上TMI、磁場矢量和磁梯度張量相互關(guān)系的解析公式推導(dǎo)。

3 插值處理

通過測量平面上的垂直或水平的磁場矢量,可以計算出梯度張量。利用張量梯度代替總場或者矢量場的優(yōu)勢在于張量分量不是地球磁場方向的函數(shù),所以數(shù)據(jù)的等高線不會受到總場或垂直梯度所造成的傾斜問題。因此,張量梯度可以更容易解釋或增強從全域中獲得的信息。計算張量分量的信噪比本質(zhì)上是原始磁場矢量測量的信噪比。但是在信號轉(zhuǎn)換計算過程中,高頻噪聲會被加強,且TMI轉(zhuǎn)換為梯度張量時對噪聲的放大作用,遠大于矢量場轉(zhuǎn)換為梯度張量,因此在轉(zhuǎn)換過程中,需要進行消噪處理。同時還應(yīng)該要注意數(shù)據(jù)網(wǎng)格的間距,數(shù)據(jù)點的總個數(shù)除以網(wǎng)格間距決定了二維數(shù)字傅里葉變換中的分辨率。

因此,使用反距離加權(quán)算法(IDWA)對數(shù)據(jù)圓滑網(wǎng)格化處理。IDWA計算一個格網(wǎng)結(jié)點時給予一個特定數(shù)據(jù)點的權(quán)值,權(quán)值與指定方向的結(jié)點到觀測點的距離的倒數(shù)成比例,較近的數(shù)據(jù)點被給定一個較高的權(quán)值[4-5],同時在此過程中,數(shù)據(jù)的高頻噪聲得到抑制,數(shù)據(jù)等高線的邊界更加柔和,消除了數(shù)據(jù)點的畸變。

圖3 空間域梯度張量

由于基于FFT的算法在波數(shù)域中是離散化的,采樣點之間數(shù)據(jù)連續(xù)性可以通過以下過程自動實現(xiàn)。一般來說,在應(yīng)用波數(shù)域中的任何轉(zhuǎn)換之前,輸入數(shù)據(jù)網(wǎng)格必須以特殊的方式進行處理[15]。采用的方法是將減去了所有線性二維趨勢的數(shù)據(jù)擴展到具有2N個點的方陣,并采用反距離加權(quán)算法處理,使其與常規(guī)數(shù)據(jù)保持一定距離,以確保滿足FFT算法要求的連續(xù)性和周期性邊界條件。這些預(yù)處理措施可以快速準確地計算網(wǎng)格。在對變換后的數(shù)據(jù)進行逆傅里葉變換后,將被移除的線性趨勢解析后添加到輸出里。為了將TMI轉(zhuǎn)換成磁場矢量,將減去的趨勢與地磁背景場方向矢量的相應(yīng)分量相乘并添加到結(jié)果中。最后一步,對于磁梯度張量的垂直和整體積分,將水平梯度的傅里葉變換除以算子核p和q后進行高通濾波。這可以消除由于合并而在相鄰網(wǎng)格邊緣的磁異常上產(chǎn)生的人為影響,最后簡單地加上地磁背景場的絕對值。

4 實驗分析

基于圖1和圖2所示的磁異常信號探測設(shè)備,設(shè)計了相對應(yīng)的實驗。以驗證上述方法對磁異常信號的預(yù)處理作用。

實驗測量了直徑15 cm,長1.2 m的鐵圓柱在2 m的埋深下,其正上方20×20個點的水平測量面的磁異常數(shù)據(jù)。如圖2所示的正方形實驗臺,邊長兩米,每隔0.1 m設(shè)置一個測量點,通過橫縱移動磁通門傳感器測量每個點的矢量三分量。再分別利用空間域和本文提出的波數(shù)域兩種方法,由磁通門測得的磁矢量場信號計算出TMI和梯度張量信號,分別畫出相應(yīng)的磁異常信號的等高線圖,比較兩種方法所計算得到的信號的質(zhì)量和效用。

圖2 實驗測量臺

圖3為空間域方法計算所得的5個獨立張量元素[11-13](Bxx、Bxy、Bxz、Byy、Byz)的等高線圖。

圖4為利用本文提出的波數(shù)域方法得到的相應(yīng)5個元素。圖5中,(a)為空間域上計算所得的TMI,(b)則是在波數(shù)域上計算所得的TMI。很明顯可見,磁異常數(shù)據(jù)中高頻噪聲被抑制,圖形塊邊緣的畸變點得到柔和處理。數(shù)據(jù)的整體層次分布更加簡潔明顯。更加方便從中提取目標形狀和位置特征信息。相比于圖3和圖5(a),圖4和圖5(b)磁場的特征更加明顯,并且結(jié)合了磁梯度張量分量的兩個特征:低噪聲水平和與磁通門讀出的磁場矢量分量相比的更清晰的磁場描述。由于疊加的水平梯度張量分量增加了它們對降低磁場矢量分量產(chǎn)生的噪聲水平的貢獻,所以水平波數(shù)p和q的分割具有更低的噪聲,因此圖中顯示的邊界更加柔和,特征更加突出,條紋條理更加清晰。

而在波數(shù)域上積分的另一個優(yōu)點是它突出了更短的波長細節(jié),并且可以通過快速傅立葉變換算法(例如,在MATLAB中)快速執(zhí)行。因此,在波數(shù)域上積分比在空間域中更加準確和方便。

圖5 空間域和波數(shù)域TMI

5 結(jié)論

本文提出了在波數(shù)域上利用Fourier變換,計算推導(dǎo)TMI、磁場矢量和梯度張量三者關(guān)系的方法,并結(jié)合反距離加權(quán)插值算法,對磁通門所測量的的磁場矢量數(shù)據(jù)進行了預(yù)處理,較為準確地計算出了磁性目標的張量信號,并在一定程度上抑制了噪聲,證明了波數(shù)域上計算的可行性。相比于傳統(tǒng)的空間域上得計算方法,波數(shù)域所計算的梯度張量信號質(zhì)量更高、邊界更清晰、噪音更小。

證明了可以通過對某類測量系統(tǒng)采集的信號進行適當?shù)臄?shù)據(jù)處理,然后通過波數(shù)域的推導(dǎo)計算得到高質(zhì)量和空間分辨率的另一形式的磁信號。并且這種轉(zhuǎn)換是雙向的,我們可以從磁場矢量中得到磁梯度張量,反之亦然。