婁 兵,陳建友,黃 堅(jiān),程志國,李曉峰,張 蓉

(1.中國石油新疆油田勘探開發(fā)研究院,新疆克拉瑪依834000;2.酒泉衛(wèi)星發(fā)射中心,甘肅酒泉732750;3.西安交通大學(xué)電子與信息學(xué)部,陜西西安710049;4中國石油新疆油田勘探開發(fā)研究院地球物理研究所,新疆烏魯木齊830013)

傳統(tǒng)的炸藥震源對(duì)環(huán)境破壞大且可控性差?？煽卣鹪词遣煌谡ㄋ幷鹪吹牧硪环N形式的激發(fā)源?？煽卣鹪醇ぐl(fā)方式具有安全環(huán)保、組織靈活、激發(fā)信號(hào)人為可控等優(yōu)點(diǎn),同時(shí)能夠滿足寬方位、高密度、點(diǎn)激發(fā)等高精度地震勘探的需求,在國內(nèi)外得到越來越廣泛的關(guān)注和應(yīng)用[1-3]。目前,滑動(dòng)掃描是可控震源用于地震勘探的一種常見的高效采集方式,該方法壓縮了相鄰兩次信號(hào)掃描的間隔時(shí)間,極大地提高了野外采集效率[4-7]。

可控震源滑動(dòng)掃描采集技術(shù)日趨成熟,采用該方法能夠得到高密度炮點(diǎn)、高覆蓋次數(shù)的地震數(shù)據(jù),但其最大問題是諧波畸變[8]。李振春等[9]研究了可控震源產(chǎn)生的兩種諧波形成機(jī)理與特征,認(rèn)為采用可控震源滑動(dòng)掃描采集的地震數(shù)據(jù)產(chǎn)生諧波噪聲的主要原因有兩個(gè):一是機(jī)械系統(tǒng)的畸變導(dǎo)致產(chǎn)生諧波,二是表層的介質(zhì)畸變導(dǎo)致產(chǎn)生諧波。根據(jù)諧波干擾產(chǎn)生原因,李振春等[9]將這兩類諧波分別稱為機(jī)械系統(tǒng)諧波與表層響應(yīng)諧波,同時(shí),還歸納總結(jié)了兩種諧波噪聲的幾個(gè)重要特征:①機(jī)械系統(tǒng)諧波頻率為基波頻率整數(shù)倍,表層響應(yīng)諧波頻率集中在一段狹窄的頻率范圍內(nèi);②當(dāng)使用線性升頻信號(hào)作為掃描信號(hào)時(shí),機(jī)械系統(tǒng)諧波出現(xiàn)在負(fù)時(shí)間軸,表層響應(yīng)諧波主要出現(xiàn)在正時(shí)間軸;③機(jī)械系統(tǒng)諧波噪聲在整個(gè)剖面出現(xiàn),表層響應(yīng)諧波主要出現(xiàn)在本炮近偏移距附近,兩種諧波噪聲出現(xiàn)位置都具有一定的隨機(jī)性,難以預(yù)測(cè)。

針對(duì)滑動(dòng)掃描采集方式產(chǎn)生的諧波干擾,國內(nèi)外學(xué)者提出了不同的解決思路與方案,主要有改善可控震源構(gòu)造、相關(guān)前諧波噪聲壓制和相關(guān)后諧波噪聲壓制3類方法。SERIFF[10]分析了可控震源的基板和地面耦合導(dǎo)致高階諧波這一現(xiàn)象。張宏樂[11]提出了通過減小可控震源的驅(qū)動(dòng)幅度、增大系統(tǒng)阻尼系數(shù)、采用變相位的掃描信號(hào)3種改善可控震源的方式來減小可控震源諧波畸變。針對(duì)相關(guān)前的地震記錄諧波噪聲壓制,近年來許多學(xué)者提出了多種思路與方案。LI等[12]提出了一種利用純相移方法壓制滑動(dòng)掃描地震數(shù)據(jù)中諧波干擾的方法,該方法比較簡(jiǎn)單而且容易實(shí)現(xiàn),可以有效壓制近場(chǎng)和遠(yuǎn)場(chǎng)諧波噪聲。黃建平等[13]、周學(xué)鋒[14]、李鳳磊[15]設(shè)計(jì)了純相移濾波器來濾除諧波,該方法比較穩(wěn)定且容易實(shí)現(xiàn)。ZHANG等[16]提出反相關(guān)法,依次壓制二次諧波、三次諧波以及更高階次的諧波噪聲,該方法只需先對(duì)相關(guān)后的數(shù)據(jù)做反相關(guān)處理且省時(shí)、高效,可以處理相關(guān)前的地震數(shù)據(jù),也可以處理相關(guān)后的地震數(shù)據(jù)。針對(duì)可控震源基板與地面的諧波畸變,XIA等[17]提出了利用時(shí)頻域奇異值分解(Singular Value Decomposition,SVD)壓制諧波噪聲的方法,可以達(dá)到壓制諧波噪聲的目的,但該方法針對(duì)的是可控震源基板與地面耦合產(chǎn)生的諧波畸變,無法壓制可控震源自身產(chǎn)生的諧波干擾。相關(guān)前地震數(shù)據(jù)諧波噪聲的壓制方法可以有效壓制諧波噪聲,但相關(guān)前的地震數(shù)據(jù)量大,計(jì)算效率相對(duì)低,需要的數(shù)據(jù)存儲(chǔ)空間更大,因此,為滿足需求,研究人員將大量精力集中在壓制相關(guān)后地震數(shù)據(jù)諧波噪聲方面。YU等[18]提出了分頻異常振幅壓制方法,該方法相當(dāng)于在時(shí)頻域進(jìn)行濾波處理,能夠壓制一部分諧波與面波,但會(huì)明顯損傷有效信號(hào)。SHARMA等[19]提出使用模擬退火法壓制諧波噪聲,該方法能夠壓制一定的諧波噪聲,但算法效率較低。SICKING等[20]提出地面力信號(hào)濾波法,此方法需要相關(guān)后數(shù)據(jù)、地面力信號(hào)和每一炮數(shù)據(jù)初至?xí)r間。WANG等[21]基于震源力信號(hào)設(shè)計(jì)諧波預(yù)測(cè)算子,求取相關(guān)后記錄中的諧波干擾,然后根據(jù)求得的諧波干擾,在被干擾區(qū)域中剔除諧波噪聲,以達(dá)到壓制地震數(shù)據(jù)中諧波干擾的目的,但該方法需要震源力信號(hào),且需要在上一炮數(shù)據(jù)中減掉本炮對(duì)上一炮的諧波干擾,因此必須從沒有受到諧波干擾的最后一炮數(shù)據(jù)開始處理,所以該算法無法進(jìn)行單炮數(shù)據(jù)處理。LI等[22]和羅勇等[23]利用時(shí)頻域稀疏優(yōu)化方法壓制相關(guān)后地震記錄數(shù)據(jù)中的諧波噪聲,該方法在確保有效信號(hào)保真性的同時(shí),不需要震源力信號(hào)就能夠直接有效壓制諧波噪聲,且可以進(jìn)行單炮處理,降低了采集成本。

由于可控震源與檢波器距離遠(yuǎn)近及地表介質(zhì)結(jié)構(gòu)的差異,導(dǎo)致地震數(shù)據(jù)中諧波噪聲位置分布與各道數(shù)據(jù)中諧波噪聲強(qiáng)弱難以預(yù)測(cè)。在利用時(shí)頻域稀疏優(yōu)化方法對(duì)諧波干擾噪聲進(jìn)行壓制時(shí),如果不考慮諧波噪聲強(qiáng)弱,往往會(huì)導(dǎo)致對(duì)某些道的諧波噪聲壓制不徹底或損傷有效信號(hào)的問題。因此需要考慮地震數(shù)據(jù)所含諧波噪聲的強(qiáng)弱,針對(duì)性地選擇算法參數(shù)。如果只憑直覺和經(jīng)驗(yàn)判斷諧波噪聲強(qiáng)弱,不僅效率很低,而且也不夠科學(xué)和精確。因此,迫切需要一種能夠自動(dòng)衡量可控震源滑動(dòng)掃描采集地震數(shù)據(jù)中諧波噪聲強(qiáng)弱的方法。為此,本文提出了一種衡量自動(dòng)相關(guān)后地震資料中諧波噪聲能量的振幅譜比值方法,并應(yīng)用于諧波噪聲的自適應(yīng)壓制。

1 諧波噪聲振幅譜比值衡量方法

1.1 諧波噪聲分析

在可控震源滑動(dòng)掃描施工中,采集到的地震數(shù)據(jù)往往受諧波噪聲的干擾。圖1為西部某油田可控震源滑動(dòng)掃描采集的相關(guān)后炮集記錄,每條測(cè)線480道,記錄長度為7.0s,采樣間隔為2ms。從圖1中可以看到,諧波噪聲的影響較為嚴(yán)重,尤其在黑框區(qū)域,諧波噪聲掩蓋了部分同相軸,降低了數(shù)據(jù)信噪比,嚴(yán)重影響地震數(shù)據(jù)后續(xù)的分析與處理。

圖1 可控震源滑動(dòng)掃描實(shí)際地震數(shù)據(jù)

可控震源地震勘探使用的線性掃描信號(hào)頻率是隨時(shí)間線性變化的,如下式所示:

(1)

式中:f1為掃描信號(hào)起始頻率;q為掃描信號(hào)頻率變化率,其中,q>0表示s(t)為線性升頻掃描信號(hào),q<0表示s(t)為線性降頻掃描信號(hào),本文選擇q>0;T為掃描信號(hào)持續(xù)時(shí)間,即掃描長度;a(t)為掃描信號(hào)的窗函數(shù)。a(t)可表示為:

a(t)=

(2)

式中:T1為該窗函數(shù)在掃描信號(hào)兩端增加的一個(gè)逐漸變化的過渡帶的長度。

諧波畸變掃描信號(hào)的波形與時(shí)頻譜如圖2所示。諧波畸變掃描信號(hào)可以定義為:

圖2 諧波畸變信號(hào)波形(a)及其時(shí)頻譜(b)

(3)

式中:an(t)表示基波以及各次諧波的振幅包絡(luò)(n=1表示參考掃描信號(hào),n=2表示二階諧波,n=3表示三階諧波)。一般從地面力信號(hào)能夠分析得到諧波畸變階數(shù),由于更高階的諧波噪聲經(jīng)過地層衰減后殘留能量很低,因此本文只考慮二階、三階的諧波噪聲。

可控震源勘探施工包括采集和相關(guān)兩個(gè)過程。所謂相關(guān),是指將檢波器上原始采集的地震記錄數(shù)據(jù)與參考掃描信號(hào)進(jìn)行互相關(guān),使得相關(guān)后的諧波噪聲分布在負(fù)的時(shí)間軸上(對(duì)于線性升頻掃描信號(hào))。假設(shè)x(t)為檢波器記錄的信號(hào),s(t)為參考掃描信號(hào),兩者做互相關(guān)處理后,獲得相關(guān)后的數(shù)據(jù)s1(t):

(4)

式中:τ為信號(hào)的時(shí)間。

根據(jù)滑動(dòng)掃描采集理論,為了提高采集效率,本炮尚未接收完畢就可開始進(jìn)行下一炮的激發(fā)。圖3 為相關(guān)后的時(shí)域波形及其時(shí)頻譜(以升頻掃描信號(hào)為例),可見,相關(guān)后的諧波噪聲出現(xiàn)在負(fù)的時(shí)間軸上。圖4為相關(guān)后的模擬地震記錄時(shí)頻譜,圖中紅圈處參考掃描信號(hào)與下一炮諧波噪聲有交叉重疊部分,也就是本炮數(shù)據(jù)主要受到下一炮諧波噪聲干擾。

圖4 相關(guān)后模擬地震信號(hào)(兩炮數(shù)據(jù))時(shí)頻譜

圖3 相關(guān)后模擬地震信號(hào)(單炮數(shù)據(jù))波形(a)及其時(shí)頻譜(b)

1.2 振幅譜比值計(jì)算

圖5展示了起始頻率為5Hz,終止頻率為90Hz的模擬數(shù)據(jù)滑動(dòng)掃描相關(guān)后的時(shí)頻譜,可以清晰地看到有效信號(hào)(基波)及噪聲干擾(二次諧波),也可清晰觀察到二次諧波起始頻率為10Hz,三次諧波起始頻率為20Hz,與理論推導(dǎo)結(jié)果一致。在實(shí)際勘探中,為了節(jié)省存儲(chǔ)及計(jì)算時(shí)間,一般根據(jù)最深目標(biāo)層確定相關(guān)后數(shù)據(jù)長度,會(huì)將完整的相關(guān)后數(shù)據(jù)截短,如圖5 中黃色虛線所示。此時(shí),截取后的數(shù)據(jù),主要保留了諧波噪聲的高頻能量部分,即諧波干擾的起始頻率將抬高,不再是原來的10Hz及20Hz。

圖5 單道信號(hào)的時(shí)頻譜

自動(dòng)衡量可控震源滑動(dòng)掃描采集地震記錄中諧波噪聲強(qiáng)弱的振幅譜比值方法的計(jì)算流程如圖6所示,主要步驟如下。

圖6 振幅譜比值方法的計(jì)算流程

1) 假設(shè)x[n]表示單道地震數(shù)據(jù),X[k]為x[n]的離散傅里葉變換(DFT):

X[k]=DFT(x[n])

(5)

假設(shè)單道地震數(shù)據(jù)x[n]的采樣點(diǎn)個(gè)數(shù)為N,則通過DFT得到的頻譜離散采樣點(diǎn)個(gè)數(shù)也為N;本文記采樣頻率為fs,由于DFT算法得到的頻譜以尼奎斯特(Nyquist)頻率對(duì)稱,因此最多只需考慮前N/2個(gè)離散點(diǎn)的頻譜。

2) 假設(shè)基于震源力信號(hào)起始低頻值及相關(guān)后數(shù)據(jù)長度確定的高、低頻分界點(diǎn)的頻率為fp,可以確定地震記錄高頻能量計(jì)算起始點(diǎn)p的頻率為Nfp/fs,即諧波噪聲在DFT域主要能量分布起始位置,計(jì)算總頻帶能量及分界點(diǎn)以上的高頻帶能量。

(6)

(7)

式中:Efh為地震記錄高頻能量;Ef為地震數(shù)據(jù)總頻帶能量;假設(shè)諧波噪聲頻帶范圍上限為fmax,則M可由Nfmax/fs確定。

3) 計(jì)算分界點(diǎn)高頻帶能量與總頻帶能量的比值,即為振幅譜比值α:

α=Efh/Ef

(8)

2 地震數(shù)據(jù)實(shí)驗(yàn)

2.1 模型數(shù)據(jù)

用簡(jiǎn)單的層狀介質(zhì)模型數(shù)據(jù)進(jìn)行測(cè)試。圖7a、圖7b分別為合成的有效信號(hào)和諧波噪聲模型數(shù)據(jù),共300道,記錄長度為6.0s,采樣間隔為2ms。從圖7b 中可以看到,近偏移距的諧波噪聲強(qiáng)于遠(yuǎn)偏移距。將圖7a的有效信號(hào)與圖7b的諧波噪聲相加,得到圖7c所示的含噪模型數(shù)據(jù),信噪比為23.86。將圖7a 的有效信號(hào)加上10倍的圖7b的諧波噪聲,得到圖7d的含噪模型數(shù)據(jù),信噪比為-22.19。很明顯,圖7d合成數(shù)據(jù)的諧波噪聲要強(qiáng)于圖7c合成數(shù)據(jù)的諧波噪聲。

分別計(jì)算這兩種含噪合成數(shù)據(jù)各道的振幅譜比值。基于震源力信號(hào)起始低頻值及相關(guān)后數(shù)據(jù)長度確定高、低頻的分界點(diǎn)的頻率fp為40Hz,即40～100Hz作為高頻帶,0～100Hz作為全頻帶。計(jì)算的振幅譜比值如圖8所示。從圖8 可以看到,含強(qiáng)諧波噪聲合成數(shù)據(jù)每一道的振幅譜比值高于同一道的含弱諧波噪聲數(shù)據(jù),說明單道地震數(shù)據(jù)諧波噪聲越強(qiáng),其振幅譜比值越大;此外,無論是圖7c的弱噪模型數(shù)據(jù)還是圖7d的強(qiáng)噪模型數(shù)據(jù),其振幅譜比值曲線都是先增大后減小,與含噪模型數(shù)據(jù)的諧波噪聲在近偏移距要強(qiáng)于遠(yuǎn)偏移距的特點(diǎn)一致。分別抽取圖7c與圖7d所示數(shù)據(jù)的第50道、第150道、第200道,計(jì)算其振幅譜比值,對(duì)比結(jié)果如表1 所示,可以看到,諧波噪聲越強(qiáng),振幅譜比值越大。

圖7 模型數(shù)據(jù)

圖8 不同信噪比合成數(shù)據(jù)的振幅譜比值

表1 單道合成數(shù)據(jù)的振幅譜比值

2.2 實(shí)際數(shù)據(jù)

計(jì)算圖1所示的實(shí)際地震數(shù)據(jù)的振幅譜比值?；谡鹪戳π盘?hào)起始低頻值及相關(guān)后數(shù)據(jù)長度,我們選擇50Hz作為高、低頻帶的分界頻率,即選擇50～100Hz作為高頻帶,0～100Hz作為全頻帶,計(jì)算振幅譜比值,結(jié)果如圖9所示。從圖9可以看到,在黑框內(nèi)對(duì)應(yīng)的區(qū)域振幅譜比值較大,其它區(qū)域振幅譜比值較小,也就是說,諧波噪聲越強(qiáng)的區(qū)域,振幅譜比值越大。

圖9 實(shí)際地震數(shù)據(jù)的振幅譜比值

為了進(jìn)一步驗(yàn)證振幅譜比值方法,我們從圖1所示的實(shí)際地震資料中抽取第50道、第100道和第210道數(shù)據(jù)進(jìn)行分析,其波形和時(shí)頻譜分別如圖10和圖11所示。可以看出,第50道含有較弱諧波噪聲,第100道含有中等強(qiáng)度諧波噪聲,第210道含有較強(qiáng)諧波噪聲。采用本文方法計(jì)算振幅譜比值,第50道數(shù)據(jù)振幅譜比值為0.6095,第100道數(shù)據(jù)振幅譜比值為0.6947,第210道數(shù)據(jù)振幅譜比值為0.8912,反映出單道數(shù)據(jù)的諧波噪聲越強(qiáng),其振幅譜比值越大。因此,采用本文方法可以自動(dòng)衡量滑動(dòng)掃描地震數(shù)據(jù)的諧波噪聲強(qiáng)弱。

圖11 第50道(a)、第100道(b)及第210道(c)數(shù)據(jù)時(shí)頻譜

圖10 第50道(a)、第100道(b)及第210道(c)數(shù)據(jù)波形

自動(dòng)衡量諧波噪聲強(qiáng)弱的振幅譜比值方法是一種高效算法,采用該方法處理10×104道(每道3500個(gè)采樣點(diǎn))地震數(shù)據(jù)僅耗時(shí)62.7s(單核處理)。

從模型數(shù)據(jù)與實(shí)際數(shù)據(jù)的實(shí)驗(yàn)結(jié)果可以看出,諧波噪聲越弱,振幅譜比值越小,反之,諧波噪聲越強(qiáng),振幅譜比值越強(qiáng)。因此,使用振幅譜比值可以對(duì)三維工區(qū)地震數(shù)據(jù)的每一道數(shù)據(jù)中諧波噪聲能量的強(qiáng)弱進(jìn)行自動(dòng)估計(jì),并且可以根據(jù)諧波噪聲強(qiáng)弱自動(dòng)調(diào)節(jié)諧波噪聲壓制算法中的一些關(guān)鍵參數(shù)。例如,在使用基于稀疏優(yōu)化壓制諧波噪聲算法中,當(dāng)諧波噪聲能量較強(qiáng)時(shí),則可調(diào)小稀疏表示諧波噪聲字典的閾值系數(shù)參數(shù),反之則調(diào)大閾值系數(shù),以實(shí)現(xiàn)諧波噪聲的自適應(yīng)壓制。

圖12為可控震源滑動(dòng)掃描采集的實(shí)際疊后地震數(shù)據(jù),共有420道,采樣長度為6.0s,采樣間隔為2ms,可見,該地震數(shù)據(jù)受到較強(qiáng)諧波噪聲的干擾。圖13a為該數(shù)據(jù)中局部區(qū)域(201～300道,4.0～6.0s)的放大顯示。首先采用固定參數(shù)的稀疏優(yōu)化方法(即處理各道數(shù)據(jù)時(shí)參數(shù)相同,通過多次實(shí)驗(yàn)得到最佳的固定參數(shù))對(duì)該地震數(shù)據(jù)進(jìn)行去噪處理,并選擇201～300道,4.0～6.0s的時(shí)窗區(qū)域進(jìn)行放大觀察,結(jié)果如圖13b所示。采用本文的振幅譜比值方法衡量每道數(shù)據(jù)的諧波噪聲強(qiáng)弱,然后根據(jù)噪聲強(qiáng)弱自適應(yīng)地選擇壓噪算法的參數(shù),并同樣選擇201～300道,4.0～6.0s的時(shí)窗區(qū)域進(jìn)行觀察,去噪結(jié)果如圖13c 所示?？梢钥吹?兩種方法都能夠壓制5s以內(nèi)的淺層區(qū)域諧波噪聲,但是觀察圖13中黑框區(qū)域可知,采用固定參數(shù)方法進(jìn)行去噪處理后具有明顯的諧波噪聲殘留,而采用自適應(yīng)壓噪方法對(duì)諧波噪聲的壓制效果則更加徹底,幾乎沒有諧波噪聲殘留。為了進(jìn)一步驗(yàn)證此結(jié)論,抽取第240道數(shù)據(jù)分析其時(shí)頻譜,結(jié)果如圖14 所示,同樣可以看到,采用固定參數(shù)的稀疏優(yōu)化方法殘留諧波噪聲比較明顯,采用自適應(yīng)壓噪方法能夠更加徹底壓制諧波噪聲。因此,基于振幅譜比值的諧波噪聲自適應(yīng)壓噪方法有效解決了固定參數(shù)的稀疏優(yōu)化方法可能帶來的諧波噪聲壓制不徹底的問題,具有很強(qiáng)的自適應(yīng)性。

圖12 實(shí)際地震數(shù)據(jù)

圖13 采用不同方法壓制諧波噪聲的結(jié)果

圖14 采用不同方法壓制諧波噪聲的時(shí)頻譜

3 結(jié)論

單道地震數(shù)據(jù)從時(shí)域變換到頻率域,基于震源力信號(hào)起始低頻值及相關(guān)后數(shù)據(jù)長度確定高、低頻分界點(diǎn),計(jì)算高頻能量占總頻帶能量的比值,作為衡量諧波噪聲強(qiáng)弱的指標(biāo)。模型數(shù)據(jù)與實(shí)際數(shù)據(jù)實(shí)驗(yàn)結(jié)果表明,諧波噪聲能量越強(qiáng),振幅譜比值越大。振幅譜比值能夠自動(dòng)估計(jì)可控震源滑動(dòng)掃描地震記錄中諧波噪聲能量的強(qiáng)弱。根據(jù)振幅譜比值自動(dòng)選擇諧波噪聲壓制算法中的參數(shù),可以實(shí)現(xiàn)諧波噪聲的自適應(yīng)壓制。