基于Butterworth 濾波器的線性噪聲壓制方法研究及實現(xiàn)

李子偉，鄒明亮，曹成寅，吳曲波

（1. 核工業(yè)北京地質(zhì)研究院, 北京100029； 2. 核工業(yè)二三〇研究所, 長沙 410007）

作為發(fā)展核電和鞏固國防的戰(zhàn)略資源,鈾礦有著不同于其他礦產(chǎn)的特殊地位。 21 世紀初我國鈾礦戰(zhàn)略轉(zhuǎn)移至北方, 此后的十余年落實和擴大了一批砂巖型鈾礦床, 資源增長量很快, 在新增探明的資源量中約占4/5,也使得砂巖型鈾礦在我國鈾資源總量中所占的比例超過了40%, 因此我國砂巖型鈾礦存在巨大的發(fā)展?jié)摿1-2]。

砂巖型鈾礦勘查實踐表明利用地球物理方法查明地層結(jié)構(gòu)、 構(gòu)造和砂體分布規(guī)律對于砂巖型鈾礦的研究具有重要意義[2-6]。 在眾多地球物理方法中, 地震勘探方法具有較大的探測深度和較高的探測精度[5-6], 是未來砂巖型鈾礦地球物理勘探中不可或缺的主要方法之一。 隨著勘探深度的增加, 砂巖型鈾礦的地質(zhì)條件越來越復雜, 這對地震資料的處理也提出了更高的要求, 其中地震資料的信噪比是保證地震資料品質(zhì)的基礎, 因此提高信噪比是地震數(shù)據(jù)處理的首要任務。

渥·伊爾馬茲等人將地震資料中的噪聲分為兩類: 隨機噪聲和相干噪聲, 其中相干噪聲又包括線性噪聲、 交混回響和多次波[7-9]。相干線性噪聲在陸上資料中以頻散瑞雷波的形式存在, 通常把它稱為面波, 這些噪聲的存在直接影響地震資料疊加剖面的成像效果,從而增加了地震資料解釋的難度。 為滿足當前砂巖型鈾礦勘查的要求, 研究與開發(fā)壓制線性噪聲的技術方法顯得十分重要。 筆者基于Butterworth 濾波器, 開展了頻率-波數(shù)域（f-k 域） 的線性噪聲壓制方法研究, 并利用Fortran 語言編制了計算機程序, 實現(xiàn)了地震數(shù)據(jù)線性噪聲的壓制, 新疆準噶爾盆地的實際地震數(shù)據(jù)處理結(jié)果表明, 該方法能夠較好地壓制線性噪聲, 有效提高地震資料的信噪比。

1 頻率-波數(shù)域的數(shù)字濾波

在地震勘探中, 用數(shù)字儀記錄地震波時,為了保持更多的地震波的特征, 通常利用寬頻帶進行數(shù)據(jù)記錄, 因此, 在寬頻帶范圍記錄了各種有效反射波的同時, 也記錄了各種干擾波。 有效波和干擾波的差異表現(xiàn)在多個方面, 例如頻譜、 傳播方向、 能量等, 在地震數(shù)字處理中利用頻譜特征的不同來壓制干擾波, 以突出有效波的方法就是數(shù)字濾波,是地震數(shù)據(jù)處理方法的基礎[9]。

圖1 頻率-波數(shù)域（f-k 域）示意圖Fig. 1 Schematic map of f-k domain

對于線性干擾而言, 它與有效信號的差異主要體現(xiàn)在頻率（f）、 視速度（v）方面, 該差異在頻率-波數(shù)域（f-k 域）易于區(qū)分, 如圖1所示, 通過原點的兩條直線的斜率就是視速度, I 區(qū)為高速干擾區(qū), II 區(qū)是有效信號區(qū),III 區(qū)是低速干擾區(qū), 常見的面波干擾即位于該區(qū)域。 可以看出, 視速度存在差異的有效信號和干擾信號在f-k 域能有效區(qū)分, 因此利用f-k 域濾波可以壓制各種線性干擾。

用X（t, x）、 Y（t, x）表示濾波前后的地震記錄, 用H（t, x）表示濾波因子, 對地震記錄X（t, x）施加濾波因子H（t, x）, 經(jīng)處理后即可得到濾波后的地震記錄Y（t, x）, 該時空域的濾波過程可以用圖2 簡要表示。

圖2 時間-空間域濾波過程示意圖Fig. 2 Schematic map of the time-space domain filtering process

以上時間-空間域濾波過程可以用二維褶積公式（1）表示:

根據(jù)褶積原理, 兩個時間函數(shù)在時間域褶積的頻譜等于兩個函數(shù)的頻譜在頻率域的乘積, 因此式（1）的時間-空間域的二維褶積運算可以在頻率-波數(shù)域完成, 即輸出信號的頻-波譜是輸入信號的頻-波譜與濾波器的頻-波譜的乘積。

在進行濾波處理時, 通過二維傅里葉變換將時空域的地震數(shù)據(jù)變換到f-k 域,得到地震數(shù)據(jù)在f-k 域的頻-波譜, 與構(gòu)建的濾波器的頻-波譜進行乘積運算, 即可得到輸出信號的頻-波譜, 再利用二維反傅里葉變換, 獲得濾波后的時空域地震數(shù)據(jù)。

2 Butterworth 濾波器

在f-k 域數(shù)字濾波過程中, 濾波器選取的好壞, 直接決定了濾波的效果。 常規(guī)的扇形濾波器, 由于通帶和阻帶間尖銳邊界的存在,會產(chǎn)生嚴重的吉布斯效應[8-9], 從而影響濾波后的地震信號的保真度。 為最大程度地消除吉布斯效應, Hale & Claerbout 提出了一種在時空域?qū)崿F(xiàn)的傾角濾波方法---Butterworth 類傾角濾波方法[10], 該方法采用的Butterworth類傾角濾波器在通帶與阻帶間變化緩慢, 其算子帶內(nèi)、 帶外幅頻特性最佳平穩(wěn), 可以避免濾波后因吉布斯效應產(chǎn)生的斷斷續(xù)續(xù)的假同相軸現(xiàn)象。

考慮到砂巖型鈾礦勘查中的地震資料的干擾主要以面波和聲波為主, 在地震剖面上表現(xiàn)為較低視速度、 較高傾角的同相軸, 因此本文在Hale & Claerbout 提出的Butterworth濾波器基礎上進行了簡化, 僅采用高通傾角濾波器濾壓制低視速度的噪聲[10-11]。

Hale & Claerbout 提出的Butterworth 高通傾角濾波器在f-k 域的頻譜響應可以用公式（3）表示:

式中: ω-弧度, ω＝2πf。

筆者對公式（3）進行簡化后可用公式（4）表示:

根據(jù)LI Jianchao 等提出的濾波參數(shù)選擇結(jié)果[12], 令α＝fNνN, 代入公式（4）, 可得到簡化后的Butterworth 高通傾角濾波器在f-k域的頻譜響應, 用公式（5）表示:

基于公式（5）表示的Butterworth 濾波器,利用Fortran 語言進行計算機實現(xiàn), 完成f-k域的線性噪音壓制。

3 線性噪聲壓制方法的實現(xiàn)

在前述數(shù)字濾波和Butterworth 濾波器的數(shù)學模型基礎上, 構(gòu)建了f-k 域的Butterworth濾波器, 在f-k 域通過施加該濾波器進行線性噪聲的濾波處理, 并采用Fortran 語言編寫了以上濾波程序的核心算法, 實現(xiàn)了線性噪聲的壓制, 達到了提高地震數(shù)據(jù)信噪比的目的,具體的計算機實現(xiàn)流程見圖3。

圖3 線性噪聲壓制方法的計算機實現(xiàn)流程及函數(shù)調(diào)用Fig. 3 Computer implementation process and function call of linear noise suppression method

在時空域和f-k 域之間的數(shù)據(jù)轉(zhuǎn)換依賴常規(guī)的傅里葉變換算法, 此處不再贅述。 本算法的核心為Butterworth 濾波器在頻率-波數(shù)域的應用, 此處簡要介紹調(diào)用函數(shù)的基本情況。

函 數(shù) 定 義 -SUBROUTINE BUTDIP（IHEAD, GATHER, IPAM, TABLE, ISTAT）；函數(shù)功能-核心算法主子程序, 執(zhí)行在頻率-波數(shù)域的傾角濾波； 返回值-IHEAD-地震數(shù)據(jù)道集道頭, GATHER-地震數(shù)據(jù)道集數(shù)據(jù)體,IPAM-濾波參數(shù)信息, TABLE-FFT 用余鉉函數(shù)表, ISTAT-程序運行結(jié)果指示信息； IPAM濾波參數(shù)信息: fN-壓制截頻, vN-壓制視速度, dx-道間距。

利用以上核心算法程序, 對讀取的地震數(shù)據(jù)開展傅里葉變換, 變換至f-k 域后施加濾波器, 然后再利用二維反傅里葉變換, 獲取濾波后的時空域地震數(shù)據(jù)。 根據(jù)線性噪聲與有效信號之間的視速度差異, 確定輸入壓制截頻參數(shù)（fN）和壓制視速度（vN）, 調(diào)整Butterworth 濾波器參數(shù), 壓制指定速度的線性噪聲, 進而實現(xiàn)不同類型線性噪聲的壓制,達到提高地震數(shù)據(jù)信噪比的目的。

4 線性噪聲壓制方法的應用效果對比

采用本文開發(fā)的地震數(shù)據(jù)處理程序, 對新疆準噶爾盆地的實際地震數(shù)據(jù)進行了處理測試, 經(jīng)處理后的單炮記錄和時間剖面的線性干擾得到了較好的壓制, 突出了有效信號,提高了地震資料的信噪比, 這對后續(xù)的地震解疑具有重要的意義。

圖4 展示了開發(fā)的處理程序在準噶爾盆地的單炮地震記錄上的應用效果。 圖4a 為原始的單炮地震記錄, 可以看到, 在原始單炮記錄上線性干擾嚴重, 雙曲線特征的有效信號被淹沒在線性干擾中。 采用同一組濾波參數(shù)（帶通濾波15～150 Hz, 壓制視速度2 000 m/s）,利用某商業(yè)軟件的f-k 濾波模塊和研發(fā)的濾波程序分別進行濾波處理, 結(jié)果分別如圖4b 和圖4c 所示。 圖4b 為某商業(yè)軟件的f-k 濾波模塊的處理效果, 可以看到線性干擾受到一定程度的壓制, 有效信號凸顯, 但線性干擾仍然有所殘留； 圖4c 為開發(fā)的濾波程序的處理效果, 可以看到, 相對某商業(yè)軟件的f-k 濾波效果, 開發(fā)的濾波程序?qū)€性噪聲的壓制效果更佳, 經(jīng)處理后的單炮記錄有效信號更加明顯。

圖5 展示開發(fā)的處理程序在新疆準噶爾盆地地震疊加剖面上的應用效果。 圖5a 為原始的地震疊加剖面, 可以看到, 在原始疊加剖面上線性干擾嚴重, 影響了有效信號的成像。 利用某商業(yè)軟件的f-k 濾波模塊和研發(fā)的濾波程序分別進行濾波處理, 結(jié)果如圖5b 和圖5c 所示。 圖5b 為某商業(yè)軟件的f-k 濾波模塊的處理效果, 可以看到線性干擾受到一定程度的壓制, 但仍然存在部分殘留, 尤其是在剖面的右側(cè), 線性干擾仍然比較嚴重, 有效信號的成像受到了較大的影響； 圖5c 為開發(fā)的濾波程序的處理效果, 可以看到, 相對某商業(yè)軟件f-k 濾波效果, 開發(fā)的濾波程序?qū)€性噪聲的壓制效果更佳, 經(jīng)處理后的疊加剖面成像效果更好, 尤其是在剖面的右側(cè),有效信號的成像效果明顯增強。

通過對比分析開發(fā)的濾波方法和某商業(yè)軟件的f-k 濾波模塊的去噪效果, 結(jié)果表明開發(fā)的基于Butterworth 濾波器的去噪方法能夠有效地壓制線性噪聲, 效果優(yōu)于常規(guī)商業(yè)軟件的f-k 濾波模塊。

圖4 準噶爾盆地去噪前后的地震單炮記錄對比Fig. 4 Comparison of seismic single-shot record before and after filtering in Junggar Basin

5 結(jié)論

開展了Butterworth 濾波器的研究, 并基于該濾波器開展了f-k 域的數(shù)字濾波研究, 利用Fortran 語言編制了計算機程序, 實現(xiàn)了線性噪聲的壓制。 利用新疆準噶爾盆地的地震數(shù)據(jù)進行了實際測試, 單炮記錄和疊加剖面的測試結(jié)果表明, 地震數(shù)據(jù)的線性噪聲得到了較好的壓制, 突出了有效信號, 成像效果得到了較大改善, 有效地提高了地震資料的信噪比。

研究成果可推廣應用到砂巖型鈾礦地震資料的去噪處理中, 可以有效提高地震剖面的信噪比和成像效果, 有利于降低地震解釋的難度, 保障地質(zhì)解釋的可靠性。