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

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

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

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

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

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

1 頻率-波數域的數字濾波

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

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

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

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

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

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

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

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

2 Butterworth 濾波器

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

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

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

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

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

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

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

3 線性噪聲壓制方法的實現

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

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

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

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

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

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

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

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

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

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

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

5 結論

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

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