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

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

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

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

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

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

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

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

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

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

圖2 時(shí)間-空間域?yàn)V波過(guò)程示意圖Fig. 2 Schematic map of the time-space domain filtering process

以上時(shí)間-空間域?yàn)V波過(guò)程可以用二維褶積公式（1）表示:

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

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

2 Butterworth 濾波器

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

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

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

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

筆者對(duì)公式（3）進(jìn)行簡(jiǎn)化后可用公式（4）表示:

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

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

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

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

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

在時(shí)空域和f-k 域之間的數(shù)據(jù)轉(zhuǎn)換依賴常規(guī)的傅里葉變換算法, 此處不再贅述。 本算法的核心為Butterworth 濾波器在頻率-波數(shù)域的應(yīng)用, 此處簡(jiǎn)要介紹調(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-程序運(yùn)行結(jié)果指示信息； IPAM濾波參數(shù)信息: fN-壓制截頻, vN-壓制視速度, dx-道間距。

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

4 線性噪聲壓制方法的應(yīng)用效果對(duì)比

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

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

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

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

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

5 結(jié)論

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

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