利用相關(guān)系數(shù)法確定三種位場垂向?qū)?shù)換算的濾波參數(shù)

任朗寧，張鳳旭，郝夢成，李銀飛

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

0 引言

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

為改善導數(shù)換算的不穩(wěn)定，常規(guī)做法是通過引入低通濾波器，將其與垂向?qū)?shù)算子相結(jié)合，構(gòu)建新的垂向?qū)?shù)算子，來壓制垂向?qū)?shù)計算過程中噪聲的放大作用。Butterworth低通濾波、Hanning窗濾波以及向上延拓濾波是3種常見的低通濾波，在導數(shù)計算中取得了較好的效果，但各方法的濾波參數(shù)常常需要反復嘗試，降低了計算效率，筆者利用相關(guān)系數(shù)法來選擇合適的濾波參數(shù)，實現(xiàn)位場數(shù)據(jù)3種方法快速垂向?qū)?shù)換算，降低了人為因素對導數(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是一個放大型因子，求導過程會對高頻噪聲放大，導致垂向?qū)?shù)結(jié)果精度變差。通常的做法是在常規(guī)導數(shù)因子中附加一個低通濾波器來減小φm在高頻部分的放大作用，保留低頻有效信息。下面給出三種常規(guī)的低通濾波方法：

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

(3)

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

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

(4)

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

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

(5)

式中：h是上延高度。

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

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é)果進行Fourier反變換，即得到垂向?qū)?shù)結(jié)果。

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

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

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

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

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

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

(9)

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

2 模型試驗

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

表1 模型參數(shù)

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

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

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

3 實際應用

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

(a)含噪聲重力異常；(b)理論垂向二階導數(shù)；(c)Butterworth濾波參數(shù)選擇；(d)垂向二階導數(shù)(Butterworth濾波器)；(e)Hanning窗濾波參數(shù)選擇；(f)垂向二階導數(shù)(Hanning窗濾波)；(g)向上延拓濾波參數(shù)選擇；(h)垂向二階導數(shù)(向上延拓)。圖1 模型重力異常及各方法截止頻率處理效果對比分析圖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 各方法求導結(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濾波器垂向二階導數(shù)；(c)Hanning窗濾波垂向二階導數(shù);(d)向上延拓濾波器垂向二階導數(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ū)各方法求導結(jié)果二次積分及殘差圖Fig.5 Secondary integration and residual map of results of Zhenlai area, Songliao Basin

4 結(jié)論

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

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

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