張力起,張 猛,2,王華忠,秦廣勝

(1.波現(xiàn)象與智能反演成像研究組(WPI),同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院,上海200092;2.中國(guó)石油化工股份有限公司勝利油田分公司物探研究院,山東東營(yíng)257000;3.中國(guó)石油化工股份有限公司中原油田分公司物探研究院,河南濮陽(yáng)457001)

近十余年,油氣地震勘探技術(shù)的顯著進(jìn)展主要體現(xiàn)在“兩寬一高”地震數(shù)據(jù)采集技術(shù)的普遍使用和地震波反演成像方法技術(shù)的逐漸實(shí)用化。噪聲在地震波反演成像中影響巨大,它決定了反演成像方法的收斂性及成像質(zhì)量?！皟蓪捯桓摺钡卣饠?shù)據(jù)采集得到了(盡可能)無(wú)假頻且高維的疊前地震數(shù)據(jù),為在高維數(shù)據(jù)(信號(hào))空間中進(jìn)行徹底的噪聲壓制提供了良好的數(shù)據(jù)基礎(chǔ)。

勘探地震中的地震波場(chǎng),可以表達(dá)成信號(hào)、相干噪聲和隨機(jī)噪聲3部分的疊加,可記為D=S+Nc+Nr,其中,D為觀測(cè)數(shù)據(jù),S為地震信號(hào),Nc為相干噪聲,Nr為隨機(jī)噪聲。噪聲壓制的思想框架通常是將含噪信號(hào)視為Hilbert空間中的一個(gè)平方可積函數(shù),選取合適的基函數(shù)對(duì)信號(hào)進(jìn)行稀疏、近似的表達(dá)(信號(hào)預(yù)測(cè)的正問(wèn)題);對(duì)上述信號(hào)預(yù)測(cè)正問(wèn)題,在一定的估計(jì)理論基礎(chǔ)上,求解最佳的信號(hào)預(yù)測(cè)器,從而實(shí)現(xiàn)信號(hào)的最佳估計(jì)。在這樣的函數(shù)逼近思想下,地震數(shù)據(jù)的去噪、插值、壓縮和特征提取等正問(wèn)題是統(tǒng)一的。

在Bayes估計(jì)理論框架下,假定信號(hào)可以用既定的預(yù)測(cè)器進(jìn)行預(yù)測(cè),噪聲滿足某種統(tǒng)計(jì)分布,信號(hào)自身滿足一些先驗(yàn)信息(例如具有稀疏性),可以通過(guò)后驗(yàn)概率最大化估計(jì)出最佳預(yù)測(cè)信號(hào)。若假定信號(hào)是線性可預(yù)測(cè)的且噪聲滿足高斯分布,那么可以在均方誤差或L2范數(shù)誤差最小準(zhǔn)則下對(duì)信號(hào)實(shí)現(xiàn)最佳預(yù)測(cè)。常見(jiàn)的方法包括f-x濾波[1]、奇異譜分析(依據(jù)觀測(cè)數(shù)據(jù)的自相關(guān)函數(shù)分解)[2]。信號(hào)在一些變換空間如小波基、曲波基、物理小波基和字典基[3-5]中,具有更好的稀疏性,將其與L0/L1范數(shù)下稀疏反演理論[6]結(jié)合可以更好地實(shí)現(xiàn)信號(hào)預(yù)測(cè),達(dá)到信號(hào)分析的目的。這是當(dāng)前現(xiàn)代信號(hào)分析的重點(diǎn)發(fā)展方向。此外,按照地震數(shù)據(jù)(隨機(jī)信號(hào))的統(tǒng)計(jì)特征設(shè)計(jì)濾波器是另一種重要的濾波方法,例如高斯濾波器、中值濾波器。將此類統(tǒng)計(jì)濾波方法擴(kuò)展到高維數(shù)據(jù)的情況是必然的,例如非局部均值濾波器[7]、三維塊匹配濾波器[8-9],此時(shí)數(shù)據(jù)或圖像中結(jié)構(gòu)信息在構(gòu)建統(tǒng)計(jì)濾波器時(shí)是非常關(guān)鍵的。據(jù)此,基于信號(hào)逼近理論的濾波方法和基于信號(hào)結(jié)構(gòu)信息的統(tǒng)計(jì)濾波方法逐漸有融合的趨勢(shì)。

基于AR模型的預(yù)測(cè)濾波方法通常是在頻率空間域(f-x)采用前向或后向預(yù)測(cè)濾波方法,或利用前向(后向)預(yù)測(cè)算子的共軛對(duì)稱性同時(shí)進(jìn)行后向(前向)預(yù)測(cè)濾波,從而達(dá)到壓制噪聲的目的。GüLüNAY[10]將非因果濾波(類似高維Wiener中心濾波)應(yīng)用于三維疊后數(shù)據(jù)的去噪;WANG[11]將f-x域地震數(shù)據(jù)插值推廣至三維的f-x-y域進(jìn)行數(shù)據(jù)插值;LIU等[12]利用ARMA模型和非因果空間預(yù)測(cè)算子進(jìn)行f-x-y濾波。

“兩寬一高”地震勘探中采集的巨量地震數(shù)據(jù)(一塊上百平方千米的工區(qū)疊前數(shù)據(jù)達(dá)到十幾到幾十TB),對(duì)高維空間、高效且高質(zhì)量的噪聲壓制技術(shù)的研發(fā)有著迫切的需求。本文未采用更復(fù)雜的信號(hào)預(yù)測(cè)器(譬如高維小波變換),而是將AR模型的f-x濾波器(前向或后向預(yù)測(cè)濾波器)修改為Wiener中心預(yù)測(cè)濾波器,同時(shí)將高維數(shù)據(jù)排成Hankel矩陣的形式使其與Wiener中心預(yù)測(cè)濾波器對(duì)應(yīng),在預(yù)測(cè)誤差的L2范數(shù)定義下求解對(duì)應(yīng)的法方程獲得最佳Wiener中心濾波器系數(shù)來(lái)構(gòu)建高維線性預(yù)測(cè)濾波器,實(shí)現(xiàn)最佳濾波。對(duì)于含有復(fù)雜波場(chǎng)的地震數(shù)據(jù),為滿足局部線性假設(shè),采用對(duì)數(shù)據(jù)局部取窗的方式進(jìn)行去噪。數(shù)據(jù)窗長(zhǎng)度以及濾波算子長(zhǎng)度的選取會(huì)明顯地影響去噪效果和執(zhí)行效率。為了提高本文方法的計(jì)算效率,利用MPI進(jìn)程級(jí)并行和OpenMP線程級(jí)并行來(lái)加速計(jì)算,將本文方法應(yīng)用于合成數(shù)據(jù)和實(shí)際數(shù)據(jù),結(jié)果表明其有良好的去噪效果。

1 Bayes框架下的地震數(shù)據(jù)分析理論

地震數(shù)據(jù)預(yù)處理(地震信號(hào)處理)中包含的大多數(shù)問(wèn)題,譬如去噪、插值、多次波壓制等正問(wèn)題的構(gòu)建都可以在Bayes反演理論框架下統(tǒng)一。實(shí)際采集的地震數(shù)據(jù)被認(rèn)為是隨機(jī)信號(hào),滿足一定的概率分布。地震數(shù)據(jù)去噪的正問(wèn)題是建立對(duì)數(shù)據(jù)中所蘊(yùn)含信號(hào)的預(yù)測(cè)理論,通常利用函數(shù)逼近理論建立包含參數(shù)的信號(hào)預(yù)測(cè)模型。預(yù)測(cè)結(jié)果與實(shí)測(cè)數(shù)據(jù)的殘差是隨機(jī)的,滿足一定的先驗(yàn)概率分布。因此,隨機(jī)地震觀測(cè)數(shù)據(jù)可表示為:

(1)

(2)

式中:P(m|dobs)=P(dobs|m)P(m)/P(dobs),P(m|dobs)為后驗(yàn)概率密度函數(shù),其中,P(dobs|m)為實(shí)測(cè)數(shù)據(jù)dobs對(duì)模型參數(shù)m的先驗(yàn)概率密度函數(shù),P(m)為模型參數(shù)m的先驗(yàn)概率密度函數(shù),P(dobs)為實(shí)測(cè)數(shù)據(jù)dobs的先驗(yàn)概率密度函數(shù)。P(m|dobs)可以理解為當(dāng)觀測(cè)數(shù)據(jù)為實(shí)測(cè)數(shù)據(jù)dobs時(shí),對(duì)應(yīng)不同的參數(shù)m時(shí)的概率。在高維情況下,后驗(yàn)概率密度函數(shù)P(m|dobs)幾乎無(wú)法由實(shí)際計(jì)算得到。因此,常取后驗(yàn)概率密度函數(shù)P(m|dobs)最大值時(shí)對(duì)應(yīng)的參數(shù)m為最終的估計(jì)結(jié)果:

(3)

一般地,在假設(shè)P(dobs|m)和P(m)這兩個(gè)概率密度函數(shù)都是高斯函數(shù)且正問(wèn)題是線性的情況下,有:

(4)

式中:C為常數(shù);S(m)為代價(jià)函數(shù)。

其中,代價(jià)函數(shù)S(m)為:

(5)

式中:mprior是模型參數(shù)m的先驗(yàn)值;CD是數(shù)據(jù)自相關(guān)矩陣;CM是模型參數(shù)自相關(guān)矩陣。采用各種數(shù)值優(yōu)化算法求解公式(5),均可得到估計(jì)的模型參數(shù)m。

信號(hào)預(yù)測(cè)算子G(·)可記為G,它由各種選定的基函數(shù)族中的基函數(shù)線性疊加產(chǎn)生?；瘮?shù)族可分為兩類,第1類是預(yù)先選定的基函數(shù)簇,包括Fourier基函數(shù)、余弦基函數(shù)、Radon基函數(shù)、小波基函數(shù)和曲波基函數(shù)等;第2類是由數(shù)據(jù)驅(qū)動(dòng)獲得的基函數(shù)簇,包括K-L變換、奇異譜分析和字典學(xué)習(xí)獲得的基函數(shù)等。模型參數(shù)m為線性信號(hào)預(yù)測(cè)器的系數(shù)?；瘮?shù)實(shí)質(zhì)上是信號(hào)的特征成分?；瘮?shù)選擇恰當(dāng),則可以用很少的特征成分組合很好地表達(dá)信號(hào),這是信號(hào)(與圖像)分析中最根本的基礎(chǔ)。在上述Bayes框架下,對(duì)模型參數(shù)進(jìn)行最佳估計(jì)后,即可實(shí)現(xiàn)對(duì)數(shù)據(jù)中所蘊(yùn)含信號(hào)的最佳表達(dá)。

上述抽象的理論框架建立了數(shù)據(jù)分析的理論基礎(chǔ),原則上適用于任何信號(hào)和圖像分析。本文中提出的高維Wiener中心濾波方法則是基于信號(hào)是線性可預(yù)測(cè)且噪聲滿足高斯分布的假設(shè),在L2范數(shù)誤差最小準(zhǔn)則下對(duì)信號(hào)進(jìn)行建模和最佳預(yù)測(cè)。

2 f-x域中線性信號(hào)模型的建立與信號(hào)預(yù)測(cè)

Wiener濾波可實(shí)現(xiàn)對(duì)線性信號(hào)或相干噪聲的最佳預(yù)測(cè),達(dá)到壓制隨機(jī)噪聲(也包括相干噪聲)和提高信噪比的目的。在頻率空間域,這些線性信號(hào)可以表示成具有線性相位移的諧波疊加。時(shí)間空間域中的二維信號(hào)s(t,x)在f-x域表示為:

(6)

式中:p為水平射線參數(shù);Δx為道間距;m為道數(shù);Δt為相鄰兩道間的時(shí)移量;A(f)為子波的頻譜;φ=2πfpΔx,φ代表由p確定的線性同相軸的線性相位移。

對(duì)于含有一個(gè)線性同相軸的任意一個(gè)單頻,相鄰道之間存在如下關(guān)系:

(7)

式中:a=ei2πfpΔx;Si-1為頻率域第i-1道的值;Si為頻率域第i道的值。(7)式說(shuō)明各道信號(hào)之間存在可預(yù)測(cè)性。類比AR模型,當(dāng)含有p個(gè)線性同相軸時(shí),可以得到前向線性預(yù)測(cè)濾波器{gi},i∈[1,p]:

(8)

(9)

(10)

圖1 算子長(zhǎng)度為3時(shí)前向預(yù)測(cè)濾波(圓圈)和后向預(yù)測(cè)濾波(方框)示意

(11)

(12)

為求解濾波器系數(shù)g,需要建立如下誤差泛函J(g):

(13)

(13)式的法方程為:

(14)

(15)

式中:Df為前向預(yù)測(cè)濾波器構(gòu)建的矩陣;gf為前向預(yù)測(cè)濾波器系數(shù);Db為后向預(yù)測(cè)濾波器構(gòu)建的矩陣;gb為后向預(yù)測(cè)濾波器系數(shù)。

3 高維Wiener中心濾波

根據(jù)實(shí)際數(shù)據(jù)各自構(gòu)建法方程,然后求解得到的前向預(yù)測(cè)濾波器和后向預(yù)測(cè)濾波器(基于AR模型構(gòu)建),只適用于噪聲滿足高斯分布、振幅不隨空間變化的局部線性信號(hào)預(yù)測(cè),若不滿足該條件,則在壓制噪聲的同時(shí)會(huì)破壞信號(hào)的振幅和相位信息。Wiener中心預(yù)測(cè)方法適用于振幅緩變的局部線性信號(hào)。利用數(shù)據(jù)的多級(jí)Hankel矩陣排布來(lái)構(gòu)建Wiener中心預(yù)測(cè)濾波器的法方程,進(jìn)而實(shí)現(xiàn)高維數(shù)據(jù)的預(yù)測(cè)濾波,該方法稱為高維Wiener中心濾波方法。

3.1 Hankel矩陣的構(gòu)建

對(duì)于一維數(shù)據(jù)(向量),其構(gòu)建的Hankel矩陣是指每一條交叉對(duì)角線上的元素都相等的矩陣,即Hankel矩陣的項(xiàng)hi,j=hi-1,j+1。對(duì)于高維數(shù)據(jù),也同樣能夠構(gòu)建出Hankel矩陣,即塊狀Hankel矩陣(block Hankel matrix)。多維信號(hào)以二維數(shù)據(jù)為例,D∈Rm×n所代表的矩陣為:

(16)

(17)

(18)

式中:l2+k2-1=m。圖2為將3×3的二維數(shù)據(jù)排成塊狀Hankel矩陣結(jié)構(gòu)。

圖2 由二維數(shù)據(jù)構(gòu)建的塊狀Hankel矩陣結(jié)構(gòu)

3.2 Wiener濾波中心預(yù)測(cè)濾波

將二維數(shù)據(jù)變換到f-x域后,單個(gè)頻率片數(shù)據(jù)為一維數(shù)據(jù)。采用圖3中的一維Wiener中心預(yù)測(cè)濾波方法進(jìn)行預(yù)測(cè)濾波。假定數(shù)據(jù)長(zhǎng)度為6,Wiener中心濾波算子長(zhǎng)度為2,一維Wiener中心濾波器的構(gòu)建過(guò)程如圖4所示,可以看出,只需要對(duì)一維數(shù)據(jù)進(jìn)行Hankel排布就可以得到公式(14)中的D和S。

構(gòu)建多級(jí)Hankel矩陣可以將一維Wiener中心濾波方法推廣到高維。以二維Wiener中心預(yù)測(cè)濾波

s

s

k,l

(19)

圖5 二維Wiener中心預(yù)測(cè)濾波示意(算子大小為5×5)

圖6 根據(jù)單頻空間二維數(shù)據(jù)矩陣構(gòu)建D和向量d

若是按行優(yōu)先的方式進(jìn)行數(shù)據(jù)的Hankel矩陣排布,首先對(duì)窗內(nèi)(紅色框)數(shù)據(jù)按行優(yōu)先方式排成一行;然后將窗沿行方向滑動(dòng)一個(gè)單位長(zhǎng)度;再將窗內(nèi)數(shù)據(jù)按行優(yōu)先方式排成一行,依照上述方式滑動(dòng)窗并取數(shù)據(jù)進(jìn)行排布;直到窗滑動(dòng)到數(shù)據(jù)末端;再將窗按列滑動(dòng)一個(gè)單位長(zhǎng)度,重復(fù)上述按行方向滑動(dòng)窗取數(shù)據(jù)并排布的操作;重復(fù)以上過(guò)程直到窗滑動(dòng)到數(shù)據(jù)末端。

高維Wiener中心濾波的計(jì)算步驟與f-x濾波計(jì)算步驟類似,具體如下:

1) 將數(shù)據(jù)從時(shí)間空間域變換到頻率空間域;

2) 對(duì)每個(gè)頻率切片數(shù)據(jù)進(jìn)行Hankel矩陣排布,構(gòu)造法方程,估計(jì)濾波器系數(shù);

3) 利用濾波器系數(shù)對(duì)變換后的復(fù)數(shù)域數(shù)據(jù)進(jìn)行褶積運(yùn)算;

4) 將濾波后的結(jié)果變換回時(shí)間空間域。

為了滿足信號(hào)的空間局部線性假設(shè)條件,對(duì)數(shù)據(jù)進(jìn)行取窗處理。因三維地震數(shù)據(jù)的單個(gè)頻率片數(shù)據(jù)為二維復(fù)數(shù)域(如圖7藍(lán)色框所示),故對(duì)二維復(fù)數(shù)域數(shù)據(jù)取窗數(shù)據(jù)(圖7中黑色框所示),對(duì)缺少的數(shù)據(jù)進(jìn)行補(bǔ)零后,將其排成Hankel矩陣,最終濾波后輸出數(shù)據(jù)的區(qū)域如圖7紅色框所示,具體算法流程如圖8所示。

圖7 二維單個(gè)頻率片數(shù)據(jù)的I/O示意

圖8 實(shí)際數(shù)據(jù)處理流程

為了滿足局部線性假設(shè)條件,需對(duì)地震數(shù)據(jù)進(jìn)行取窗處理。不同的數(shù)據(jù)窗長(zhǎng)度和不同的算子長(zhǎng)度會(huì)造成去噪結(jié)果和效率的明顯差異。長(zhǎng)數(shù)據(jù)窗會(huì)破壞局部線性假設(shè),無(wú)益于提高去噪能力,同時(shí)還會(huì)增加計(jì)算量;而短數(shù)據(jù)窗則導(dǎo)致該算法去噪能力下降,無(wú)法適應(yīng)復(fù)雜波場(chǎng)的情況。因此,需要根據(jù)數(shù)據(jù)波場(chǎng)的復(fù)雜程度選取合適的數(shù)據(jù)窗長(zhǎng)度和算子長(zhǎng)度。一般情況下,算子長(zhǎng)度選取為3或5,空間窗長(zhǎng)度選取范圍為15～30,時(shí)間窗長(zhǎng)度選取范圍為100～200個(gè)采樣點(diǎn)。由于采用了分塊處理數(shù)據(jù)的方式并且預(yù)測(cè)濾波是在單個(gè)頻率片上進(jìn)行的,故可以采用并行策略加速運(yùn)算。并行主要利用MPI和OpenMP進(jìn)行加速。利用MPI多進(jìn)程分塊讀取高維窗數(shù)據(jù)體,各個(gè)進(jìn)程對(duì)各自的數(shù)據(jù)體獨(dú)立進(jìn)行高維去噪;在每個(gè)進(jìn)程均利用OpenMP多線程分別處理單個(gè)頻率片數(shù)據(jù)。

4 測(cè)試應(yīng)用及分析

對(duì)三維數(shù)據(jù)體取窗后獲得的數(shù)據(jù)體同樣是三維的,而頻率空間域中濾波算子的維度則是二維的。公式(20)為信噪比(SNR)的定義。

(20)

式中:{xi},i∈[1,N]為無(wú)噪聲數(shù)據(jù);{ni},i∈[1,N]為含噪聲數(shù)據(jù);N為數(shù)據(jù)長(zhǎng)度。

4.1 三維合成地震數(shù)據(jù)

圖9中合成的三維數(shù)據(jù)大小為600×60×60,時(shí)間采樣間隔為4ms,空間采樣間隔為1m,子波主頻為30Hz。圖9a為三維合成數(shù)據(jù)中某條測(cè)線的地震記錄,圖9b為三維合成數(shù)據(jù)加入-6dB高斯白噪聲后該測(cè)線的地震記錄(對(duì)應(yīng)圖9a同一條測(cè)線)。采用不同長(zhǎng)度的一維前向和后向預(yù)測(cè)濾波器沿兩個(gè)方向各做一次一維預(yù)測(cè)濾波,其濾波結(jié)果分別如圖9c、圖9e 所示,數(shù)據(jù)殘差剖面中均出現(xiàn)了強(qiáng)同相軸(圖9d、圖9f),說(shuō)明這種濾波方法破壞了信號(hào)結(jié)構(gòu),壓制噪聲的能力差。分別采用不同大小的數(shù)據(jù)窗和二維Wiener中心濾波算子進(jìn)行預(yù)測(cè)濾波,結(jié)果如圖10 所示?？梢钥闯?當(dāng)數(shù)據(jù)窗(100×10×10)和Wiener濾波算子(3×3)較小時(shí),濾波效果較差;當(dāng)數(shù)據(jù)窗(100×30×30)和Wiener濾波算子(5×5)較大時(shí),其數(shù)據(jù)殘差剖面包含的有效信號(hào)較小且信噪比較大,濾波效果較好。

圖11中合成的三維合成單炮炮集數(shù)據(jù)大小為1500×120×120,時(shí)間采樣間隔為4ms,空間采樣間隔為10m,子波主頻為30Hz。圖11a為三維合成炮集數(shù)據(jù)中的某條測(cè)線的地震記錄,圖11b為加入-10dB高斯白噪聲后三維合成炮集數(shù)據(jù)的地震記錄(對(duì)應(yīng)圖11a同一條測(cè)線)。采用二維Wiener中心濾波方法進(jìn)行預(yù)測(cè)濾波,圖11c顯示了數(shù)據(jù)時(shí)窗大小為300×30×30,濾波算子大小為3×3時(shí)的二維Wiener中心濾波后的地震記錄,圖11d為采用二維Wiener

圖9 三維合成數(shù)據(jù)及一維前向和后向預(yù)測(cè)濾波結(jié)果a 三維合成數(shù)據(jù)中某條測(cè)線的地震記錄; b 加入-6dB高斯白噪聲后三維合成數(shù)據(jù)中某條測(cè)線的地震記錄(對(duì)應(yīng)圖9a同一條測(cè)線); c 數(shù)據(jù)窗大小為100×10×10,算子長(zhǎng)度均為3時(shí)前向預(yù)測(cè)和后向預(yù)測(cè)濾波后的結(jié)果; d 圖9c與圖9a的數(shù)據(jù)殘差剖面(SNR為15.69dB); e 數(shù)據(jù)窗大小為100×10×10,算子長(zhǎng)度均為5時(shí)前向預(yù)測(cè)和后向預(yù)測(cè)濾波后的地震記錄; f 圖9e與圖9a的數(shù)據(jù)殘差剖面(SNR為12.98dB)

中心濾波方法得到的圖11c與圖11a的數(shù)據(jù)殘差剖面,圖中無(wú)明顯的連續(xù)的同相軸,并且噪聲幅值小?？梢钥闯鯳iener中心濾波方法在壓制噪聲的同時(shí)較好地保留了有效信號(hào)。

4.2 三維實(shí)際地震資料

圖12a為某地區(qū)實(shí)際采集的四維數(shù)據(jù)(疊前CDP道集)中某個(gè)CDP道集(三維數(shù)據(jù)體),其數(shù)據(jù)大小為1501×120×34。當(dāng)數(shù)據(jù)窗大小為200×30×20,濾波算子為3×3時(shí),采用二維Wiener中心濾波方法對(duì)數(shù)據(jù)進(jìn)行預(yù)測(cè)濾波,得到的CDP道集和殘差剖面分別如圖12b和圖12c所示?？梢钥闯龆SWiener中心濾波方法能夠壓制隨機(jī)噪聲和部分相干噪聲,在振幅存在局部變化時(shí)(圖12a和圖12b中的紅圈)能較好地保留信號(hào)。

圖11 三維合成炮集數(shù)據(jù)二維Wiener中心濾波結(jié)果a 三維合成炮集數(shù)據(jù)中某條測(cè)線的地震記錄; b 加入-10dB高斯白噪聲后三維合成炮集數(shù)據(jù)中某條測(cè)線的地震記錄(對(duì)應(yīng)圖11a同一條測(cè)線); c 數(shù)據(jù)窗大小為300×30×30,濾波算子大小為3×3時(shí)二維Wiener中心濾波后的地震記錄; d 圖11c與圖11a的數(shù)據(jù)殘差剖面(SNR為12.24dB)

圖12 實(shí)際數(shù)據(jù)二維Wiener中心濾波結(jié)果a 野外采集得到的某個(gè)CDP道集; b 數(shù)據(jù)窗大小為200×30×20,濾波算子大小為3×3時(shí)二維Wiener中心濾波后的CDP道集; c 圖12b與圖12a的數(shù)據(jù)殘差剖面

某地區(qū)的三維實(shí)際數(shù)據(jù)疊后剖面如圖13a所示,數(shù)據(jù)大小為1501×567×10。當(dāng)數(shù)據(jù)窗大小為300×20×12,濾波算子為3×3時(shí)的采用二維Wiener中心濾波方法對(duì)數(shù)據(jù)進(jìn)行預(yù)測(cè)濾波,得到的疊后剖面和殘差剖面分別如圖13b和圖13c所示。可以看出,二維Wiener中心濾波方法能夠壓制中、深層雜亂的噪聲,并且很好地保留淺、中層連續(xù)的同相軸以及斷層信息。

圖13 實(shí)際數(shù)據(jù)二維Wiener中心濾波結(jié)果a 三維野外實(shí)際數(shù)據(jù)疊后剖面; b 數(shù)據(jù)窗大小為300×20×12,濾波算子大小為3×3時(shí)二維Wiener中心濾波后的疊后剖面; c 圖13b與圖13a的數(shù)據(jù)殘差剖面

5 討論與結(jié)論

“兩寬一高”地震勘探針對(duì)的是小尺度的精細(xì)油藏和復(fù)雜油藏,相應(yīng)的成像方法也逐漸進(jìn)入以FWI+LS_RTM為標(biāo)志的反演成像階段。噪聲對(duì)成像的影響需要被重新審視。偏移成像中的隨機(jī)噪聲可通過(guò)高覆蓋的成像疊加被壓制,對(duì)成像結(jié)果影響較小;而反演成像過(guò)程中噪聲則會(huì)嚴(yán)重影響反演成像方法的收斂性及反演結(jié)果的精度(分辨率)。另外,針對(duì)當(dāng)前“兩寬一高”地震數(shù)據(jù)采集得到的巨量數(shù)據(jù)體,在實(shí)際生產(chǎn)中實(shí)用的、高維空間中可實(shí)現(xiàn)的且效率和效果達(dá)到均衡的去噪方法還是相當(dāng)缺乏的。

本文所討論的高維Wiener中心濾波方法,是由基于AR模型中的前向預(yù)測(cè)濾波器和后向預(yù)測(cè)濾波器發(fā)展得到的,其要求信號(hào)是線性可預(yù)測(cè)的并且噪聲滿足高斯分布的。對(duì)于非線性的信號(hào),可采用對(duì)數(shù)據(jù)局部取窗的方式滿足局部線性可預(yù)測(cè)的假設(shè)條件;對(duì)于噪聲不滿足高斯分布假設(shè)的信號(hào),需要修正Bayes理論框架下的信號(hào)預(yù)測(cè)理論和方法以達(dá)到更好的去噪效果。地震數(shù)據(jù)基本都是高維數(shù)據(jù)體,采用低維空間中的預(yù)測(cè)算子很難利用其高維空間信息,也無(wú)法達(dá)到滿意的去噪效果。本文通過(guò)對(duì)數(shù)據(jù)體進(jìn)行Hankel矩陣排布解決了這一問(wèn)題,結(jié)合Wiener中心預(yù)測(cè)濾波方法,實(shí)現(xiàn)了三維或者更高維數(shù)據(jù)空間中的噪聲壓制。實(shí)驗(yàn)結(jié)果表明,高維Wiener中心濾波方法能夠在有效壓制噪聲的同時(shí)較好地保留信號(hào),但存在計(jì)算效率低的問(wèn)題。本文采用MPI和OpenMP并行計(jì)算可以部分地提高計(jì)算效率。但是本文方法也存在一些問(wèn)題,如對(duì)于非規(guī)則數(shù)據(jù)、復(fù)雜波場(chǎng)(振幅變化較為強(qiáng)烈)以及在信噪比相對(duì)較低的情況下則很難最佳地預(yù)測(cè)信號(hào)并壓制噪聲。

到目前為止,地震數(shù)據(jù)去噪(插值)真正的問(wèn)題依然是所建立的地震信號(hào)預(yù)測(cè)算子不能很好地預(yù)測(cè)實(shí)測(cè)數(shù)據(jù)中的信號(hào),以及噪聲的概率分布是未知的。當(dāng)前,信號(hào)建模的合理性、參數(shù)估計(jì)的精度、算法計(jì)算效率等方面遠(yuǎn)未達(dá)到理想的程度,這使得去噪結(jié)果達(dá)不到反演成像的要求?！皟蓪捯桓摺钡臄?shù)據(jù)采集方式提供了更好的數(shù)據(jù)基礎(chǔ),但(尤其在復(fù)雜近地表探區(qū))地震數(shù)據(jù)在空間時(shí)間上的復(fù)雜變化使得當(dāng)前已有的去噪技術(shù)并沒(méi)有達(dá)到令人滿意的、滿足反演成像要求的程度,未來(lái)仍需要繼續(xù)深入研究。

致謝:感謝中石油勘探開(kāi)發(fā)研究院及西北分院、中海油研究院和湛江分公司、中國(guó)石油化工股份有限公司石油物探技術(shù)研究院和勝利油田分公司對(duì)波現(xiàn)象與智能反演成像研究組(WPI)研究工作的資助與支持。