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

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

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

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

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

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

(1)

(2)

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

圖1 實(shí)驗(yàn)測(cè)量設(shè)備

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

(3)

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

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

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

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

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

(4)

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

(5)

二維表示磁標(biāo)量的逆傅里葉變換:

在這里,把這兩種變換縮寫(xiě)為:

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

(6)

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

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

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

(7)

這是在波數(shù)域中,由磁場(chǎng)矢量推導(dǎo)垂直方向的梯度張量Gkz的公式。其傅立葉變換展開(kāi)如式(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ù)域上建立起磁場(chǎng)矢量三分量和所有磁梯度張量之間的解析推導(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和磁場(chǎng)矢量三分量在波數(shù)域上的關(guān)系:

(17)

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

3 插值處理

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

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

圖3 空間域梯度張量

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

4 實(shí)驗(yàn)分析

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

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

圖2 實(shí)驗(yàn)測(cè)量臺(tái)

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

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

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

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

5 結(jié)論

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

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