任朗寧，張鳳旭，郝夢(mèng)成，李銀飛

吉林大學(xué) 地球探測科學(xué)與技術(shù)學(xué)院,長春 130026

0 引言

垂向?qū)?shù)作為常規(guī)的位場數(shù)據(jù)處理方法，不僅在分離疊加異常、突出局部異常和識(shí)別地質(zhì)體邊界等方面具有重要作用[1--3]，還是其他一些位場數(shù)據(jù)處理方法的重要環(huán)節(jié)。例如歐拉反褶積、向下延拓[4]和重力歸一化總梯度等方法都需要先進(jìn)行垂向?qū)?shù)的換算，并且導(dǎo)數(shù)換算的結(jié)果精度將直接影響到上述位場處理方法的最終計(jì)算精度。由于導(dǎo)數(shù)換算對(duì)噪聲的敏感性，常常得不到較好的求導(dǎo)結(jié)果，影響后續(xù)異常的處理與解釋。為此，許多改進(jìn)位場垂向?qū)?shù)換算的方法被提出，它們大致可分為兩類：一是空間域法；二是波數(shù)域法?？臻g域法包括量板法、邊界單元法和樣條函數(shù)法等[5--7]。波數(shù)域法包括常規(guī)FFT法、維納濾波法[8]、補(bǔ)償圓滑濾波法[9]和正則化方法[10]等。

為改善導(dǎo)數(shù)換算的不穩(wěn)定，常規(guī)做法是通過引入低通濾波器，將其與垂向?qū)?shù)算子相結(jié)合，構(gòu)建新的垂向?qū)?shù)算子，來壓制垂向?qū)?shù)計(jì)算過程中噪聲的放大作用。Butterworth低通濾波、Hanning窗濾波以及向上延拓濾波是3種常見的低通濾波，在導(dǎo)數(shù)計(jì)算中取得了較好的效果，但各方法的濾波參數(shù)常常需要反復(fù)嘗試，降低了計(jì)算效率，筆者利用相關(guān)系數(shù)法來選擇合適的濾波參數(shù)，實(shí)現(xiàn)位場數(shù)據(jù)3種方法快速垂向?qū)?shù)換算，降低了人為因素對(duì)導(dǎo)數(shù)結(jié)果的影響。

1 方法原理

1.1 垂向?qū)?shù)中常用的低通濾波因子

位場G(x,y)與其各階垂向?qū)?shù)VDm(x,y)(m代表階數(shù))在波數(shù)域的關(guān)系式為[11]：

(1)

VDm(u,v)=φm·G(u,v)

(2)

由前人研究可知[12]，φm是一個(gè)放大型因子，求導(dǎo)過程會(huì)對(duì)高頻噪聲放大，導(dǎo)致垂向?qū)?shù)結(jié)果精度變差。通常的做法是在常規(guī)導(dǎo)數(shù)因子中附加一個(gè)低通濾波器來減小φm在高頻部分的放大作用，保留低頻有效信息。下面給出三種常規(guī)的低通濾波方法：

①Butterworth低通濾波器是一個(gè)在通帶內(nèi)擁有最大限度平坦度的濾波器，其頻率響應(yīng)為：

(3)

式中：ω為頻率；ωb為Butterworth濾波器的截止頻率;N為濾波器階次。

②Hanning窗低通濾波是位場信號(hào)處理中常用的一類窗函數(shù)濾波因子,其表達(dá)式[13]為：

(4)

式中：ωh為Hanning窗濾波器的截止頻率。

③向上延拓是將觀測平面上的異常換算到觀測平面之上，相當(dāng)于一個(gè)低通濾波[14]，在波數(shù)域中，向上延拓因子是:

(5)

式中：h是上延高度。

由上述的式(2～5)可各自構(gòu)造出新的求導(dǎo)公式：

VDB(u,v)=B(ω)·VDm(u,v)

(6)

VDH(u,v)=H(ω)·VDm(u,v)

(7)

VDY(u,v)=Y(h)·VDm(u,v)

(8)

將結(jié)果進(jìn)行Fourier反變換，即得到垂向?qū)?shù)結(jié)果。

從上式及前人研究結(jié)果可知，這些方法的濾波效果主要與低通濾波器濾波參數(shù)(截止頻率或上延高度)的選擇有關(guān)。以截止頻率的選擇為例，當(dāng)截止頻率過小時(shí)，求導(dǎo)導(dǎo)致低頻部分的有效異常損失；當(dāng)截止頻率過大時(shí)，求導(dǎo)過程中對(duì)噪聲的壓制作用較弱。這兩種情況都會(huì)導(dǎo)致導(dǎo)數(shù)換算結(jié)果的精度不高，因此需要選擇一個(gè)合適的截止頻率來平衡這兩種狀況，即在保留有用信號(hào)的同時(shí)，最大限度的減少噪聲干擾。

1.2 基于相關(guān)系數(shù)法的濾波參數(shù)選取方法

為了快速確定濾波參數(shù)，這里借鑒曾華霖確定最佳上延高度時(shí)采用的相關(guān)系數(shù)法[15]。即利用取不同濾波參數(shù)時(shí)，求導(dǎo)結(jié)果之間的相關(guān)性來確定最佳濾波參數(shù)。具體做法是：

①求原始異常的Fourier變換譜；

②在一定范圍內(nèi)，設(shè)置n個(gè)離散的濾波參數(shù)，并計(jì)算相應(yīng)的垂向?qū)?shù)；

③利用下列公式計(jì)算相鄰兩個(gè)濾波參數(shù)求導(dǎo)結(jié)果的相關(guān)系數(shù)；

(9)

④當(dāng)相關(guān)系數(shù)曲線開始趨于穩(wěn)定的時(shí)候所對(duì)應(yīng)的參數(shù)確定為最佳濾波參數(shù)。

2 模型試驗(yàn)

為檢驗(yàn)濾波參數(shù)選取的相關(guān)系數(shù)法在垂向?qū)?shù)計(jì)算中的效果， 本文設(shè)計(jì)包含4個(gè)大小埋深各不相同的長方體的組合模型進(jìn)行試驗(yàn)。 模型參數(shù)見表1。

表1 模型參數(shù)

在實(shí)際中，由于高頻噪聲干擾的始終存在，故對(duì)該組合模型產(chǎn)生的重力添加均值為0，標(biāo)準(zhǔn)差為0.1 g.u.的高斯白噪聲，其含噪聲重力異常Δg見圖1a。從圖1a中無法準(zhǔn)確識(shí)別地質(zhì)體的真實(shí)邊界位置，需要對(duì)模型重力異常進(jìn)行處理。圖1b是該模型的理論垂向二階導(dǎo)數(shù)。對(duì)比圖1a和圖1b，可以看出異常的垂向二階導(dǎo)數(shù)換算，對(duì)地質(zhì)體的異常信息進(jìn)行突出。

對(duì)含噪聲重力異常(圖1a)采用式(6～8)進(jìn)行垂向二階導(dǎo)數(shù)計(jì)算。其結(jié)果Δgzz見圖1d、圖1f和圖1h。圖1c、圖1e和圖1g分別是Butterworth低通濾波、Hanning窗濾波以及向上延拓濾波進(jìn)行低垂向二階導(dǎo)數(shù)計(jì)算的濾波參數(shù)選取圖，為驗(yàn)證相關(guān)系數(shù)法選取的濾波參數(shù)是最佳的，本文還求取了各方法導(dǎo)數(shù)計(jì)算結(jié)果與理論垂向二階導(dǎo)數(shù)之間的均方根誤差。從圖1c、圖1e和圖1g中可以看出，Butterworth低通濾波、Hanning窗濾波以及向上延拓濾波的濾波參數(shù)分別是ωb=3.413 8，ωh=5.862 0和h=0.206 9，其所在位置與均方根誤差最低的位置相對(duì)應(yīng)，說明利用相關(guān)系數(shù)法確定濾波參數(shù)是合理的。對(duì)比圖1d、1f和1h可以看出，經(jīng)過本文選取合適濾波參數(shù)，各方法垂向?qū)?shù)計(jì)算中的高頻噪聲均得到了有效壓制，圖1d和1f中的結(jié)果相較于圖1h中的處理結(jié)果更好，能夠準(zhǔn)確地看出異常體的位置和形態(tài)，與圖1b中的結(jié)果相似。圖1d和1f中的導(dǎo)數(shù)結(jié)果雖然都能夠較好顯出異常體的位置和形態(tài)，但是圖1d異常的幅值與理論值(圖1b)的差距較小，結(jié)果精度較高。

為了進(jìn)一步檢驗(yàn)各方法的垂向二階導(dǎo)數(shù)效果，對(duì)導(dǎo)數(shù)結(jié)果進(jìn)行二次積分反算，并將積分結(jié)果與原異常求殘差，其結(jié)果分別見圖2a、2c、2e和圖2b、2d、2f。對(duì)比圖2a、2c、2e和圖1a，可以看出圖2a中Butterworth低通濾波的二次積分結(jié)果與原異常最相近。再對(duì)比圖2b、2d、2f，可以看出圖2b中的殘差結(jié)果里包含的地質(zhì)體異常信息最少，而Hanning窗濾波和向上延拓處理的結(jié)果中異常信息較多，說明這兩種方法在垂向二階導(dǎo)數(shù)計(jì)算中對(duì)異常的轉(zhuǎn)換不夠徹底。綜上所述，使用相關(guān)系數(shù)法選取最佳的濾波參數(shù)的做法是合適的，利用確定的濾波參數(shù)獲得垂向二階導(dǎo)數(shù)結(jié)果中，Butterworth濾波的導(dǎo)數(shù)精度最高。

3 實(shí)際應(yīng)用

為了檢驗(yàn)上述方法在實(shí)際數(shù)據(jù)處理中的效果，對(duì)位于大興安嶺外圍的松遼盆地鎮(zhèn)賚地區(qū)的布格重力異常數(shù)據(jù)(圖4a)進(jìn)行了垂向二階導(dǎo)數(shù)計(jì)算。數(shù)據(jù)網(wǎng)格點(diǎn)距0.5 km，網(wǎng)格點(diǎn)數(shù)124×231。首先，利用相關(guān)系數(shù)法選取本文各方法最佳的濾波參數(shù)分別為ωb=0.9，ωh=0.9，h=1.8(圖3a、3b、3c)。再將選取的濾波參數(shù)代入到各個(gè)濾波器中進(jìn)行垂向二階導(dǎo)數(shù)計(jì)算，結(jié)果見圖4b、4c、4d。對(duì)比上述結(jié)果可以看出，Butterworth低通濾波的垂向二階導(dǎo)數(shù)(圖4b)處理效果較好， 其異常曲線圓滑，對(duì)噪聲的壓制作用較好，同時(shí)對(duì)淺部異常的突出作用也較強(qiáng)；Hanning窗濾波的垂向二階導(dǎo)數(shù)結(jié)果(圖4c)精度較差，淺部異常較寬緩，說明相鄰異常的分離不夠徹底。而向上延拓濾波的垂向二階導(dǎo)數(shù)(圖4d)雖然突出了大部分的淺部異常特征，但是異常曲線明顯存在畸變現(xiàn)象，說明噪聲干擾被壓制的不夠徹底。為了進(jìn)一步檢驗(yàn)求導(dǎo)結(jié)果的準(zhǔn)確性，對(duì)各方法垂向二階導(dǎo)數(shù)結(jié)果求取二次積分并與原異常做差，結(jié)果分別見圖5a、5c、5e和圖5b、5d、5f。二次積分的結(jié)果與該區(qū)原異常進(jìn)行對(duì)比，可以看出Butterworth低通濾波的結(jié)果(圖5a)與原異常形態(tài)相似，幅值相近，而其他兩種方法的異常形態(tài)與原異常存在較大差異，幅值也有不同。從圖5b、5d、5f的殘差結(jié)果中可看出，圖5b中的異常特征最少，幅值最低，說明Butterworth低通濾波求取垂向二階導(dǎo)數(shù)過程中異常轉(zhuǎn)換地較為徹底，間接證明了其垂向二階導(dǎo)數(shù)的計(jì)算結(jié)果比其他兩種方法的結(jié)果精度高，這與模型試驗(yàn)的效果相似。

(a)含噪聲重力異常；(b)理論垂向二階導(dǎo)數(shù)；(c)Butterworth濾波參數(shù)選擇；(d)垂向二階導(dǎo)數(shù)(Butterworth濾波器)；(e)Hanning窗濾波參數(shù)選擇；(f)垂向二階導(dǎo)數(shù)(Hanning窗濾波)；(g)向上延拓濾波參數(shù)選擇；(h)垂向二階導(dǎo)數(shù)(向上延拓)。圖1 模型重力異常及各方法截止頻率處理效果對(duì)比分析圖Fig.1 Comparative analysis of model gravity anomaly and processing effect of cutoff frequency of each method

(a)Butterworth二次積分結(jié)果；(b)二次積分與原異常殘差(Butterworth濾波器)；(c)Hanning窗濾波二次積分結(jié)果；(d)二次積分與原異常殘差(Hanning窗濾波)；(e)向上延拓二次積分結(jié)果；(f)二次積分與原異常殘差(向上延拓)。圖2 各方法求導(dǎo)結(jié)果二次積分及殘差圖Fig.2 Secondary integration and residual plot of each method

(a)Butterworth截止頻率選取；(b)Hanning窗截止頻率選??；(c)向上延拓上延高度選取。圖3 松遼盆地鎮(zhèn)賚地區(qū)各方法截止頻率選取圖Fig.3 Selection chart of cutoff frequency for each method of Zhenlai area, Songliao Basin

(a)西部斜坡布格重力異常；(b)Butterworth濾波器垂向二階導(dǎo)數(shù)；(c)Hanning窗濾波垂向二階導(dǎo)數(shù);(d)向上延拓濾波器垂向二階導(dǎo)數(shù)。圖4 松遼盆地鎮(zhèn)賚地區(qū)重力異常及垂向?qū)?shù)處理圖Fig.4 Gravity anomaly and vertical derivatives processing map in Zhenlai area, Songliao Basin

(a)Butterworth二次積分結(jié)果；(b)二次積分與原異常殘差(Butterworth濾波器)；(c)Hanning窗濾波二次積分結(jié)果；(d)二次積分與原異常殘差(Hanning窗濾波)；(e)向上延拓二次積分結(jié)果；(f)二次積分與原異常殘差(向上延拓)。圖5 松遼盆地鎮(zhèn)賚地區(qū)各方法求導(dǎo)結(jié)果二次積分及殘差圖Fig.5 Secondary integration and residual map of results of Zhenlai area, Songliao Basin

4 結(jié)論

(1)使用相關(guān)系數(shù)法確定了Butterworth低通濾波、Hanning窗濾波和向上延拓濾波的最佳濾波參數(shù)，實(shí)現(xiàn)了3種方法位場數(shù)據(jù)垂向?qū)?shù)的快速計(jì)算，降低了人為因素對(duì)計(jì)算結(jié)果的影響。

(2)模型試驗(yàn)結(jié)果表明，使用相關(guān)系數(shù)法選取最佳的濾波參數(shù)的做法是合理的，利用確定的濾波參數(shù)計(jì)算得到的垂向二階導(dǎo)數(shù)結(jié)果中，Butterworth低通濾波法的導(dǎo)數(shù)結(jié)果與其他兩種方法的結(jié)果相比，精度較高。

(3)松遼盆地鎮(zhèn)賚地區(qū)的實(shí)測布格重力異常數(shù)據(jù)處理檢驗(yàn)了3種垂向?qū)?shù)換算方法的實(shí)際效果，獲得了與模型試驗(yàn)相似的效果，為進(jìn)一步數(shù)據(jù)處理與解釋提供了支持。