王本鋒,陳小宏,李景葉,陳增保

(1.中國石油大學(北京)油氣資源與探測國家重點實驗室,北京102249;2.中國石油大學(北京)海洋石油勘探國家工程實驗室,北京102249)

由于地下介質的非均勻性以及黏滯性對地震波產生了吸收衰減作用(即正Q濾波效應),使其振幅減弱、相位畸變,地震記錄的分辨率降低[1]。Varela等[2]利用Futterman衰減頻散關系分析了衰減與頻散效應,驗證了地下介質的吸收衰減效應。衰減補償是提高地震記錄分辨率,進行精細油藏描述的有效途徑之一[3]。反Q濾波是衰減補償?shù)闹匾椒?其相位校正無條件穩(wěn)定且具有較強的抗噪性,但振幅補償因子隨著傳播時間的增大以及頻率的增加呈指數(shù)放大趨勢,穩(wěn)定性及抗噪性較差[4-7]。Hargreaves等[8]利用類似Stolt頻率-波數(shù)域偏移的方法,高效地實現(xiàn)了反Q濾波算法,有效地校正了由速度頻散引起的相位畸變,但為了避免振幅補償?shù)牟环€(wěn)定性,僅對相位進行了校正。趙建勛等[9]通過反Q濾波與串聯(lián)偏移相結合,將常Q算法推廣至深變Q值模型；通過串聯(lián)頻散補償算法與限幅補償算法相結合,更好地提高了地震記錄的分辨率。

針對反Q濾波振幅補償?shù)牟环€(wěn)定性,Wang[4,6]提出了振幅增益控制方法以及穩(wěn)定的反Q濾波方法(可同時補償振幅以及校正相位),設計了Gabor域算法以提高計算效率。對含噪地震記錄,振幅增益控制方法在提高分辨率的同時,較大地放大了噪聲,導致振幅補償不穩(wěn)定;穩(wěn)定的反Q濾波方法可較好地處理含噪地震記錄,但一定程度上導致振幅欠補償。Yan等[10]將穩(wěn)定的反Q濾波公式推廣至多分量地震勘探中,利用疊前轉換波道集提取S波的Q值,利用疊前縱波道集提取P波的Q值,實現(xiàn)了疊前PP波以及PS波的衰減補償,提高了疊前地震記錄的分辨率。Zhao等[11]為克服振幅補償過程中噪聲被放大的缺陷,在時頻域定義信噪比,對信噪比大于1的區(qū)域進行振幅補償,其它區(qū)域不進行補償,從而在保證信噪比的前提下,提高了地震記錄的分辨率。為減弱反Q濾波對Q值的依賴性,Braga等[12]在小波域實現(xiàn)反Q濾波方法,在Q值誤差較大時,補償后的地震記錄能在提高分辨率的同時依然保持原有的AVO特性,從而有利于疊前反演。Zhang等[13]為克服反Q濾波振幅補償不穩(wěn)定或欠補償?shù)娜毕?利用最小二乘策略以及貝葉斯理論,給定子波或提取最小相位子波來迭代反演稀疏的反射系數(shù)序列,提取Q值作相位校正,得到高分辨率的地震記錄。Wang[14]利用Futterman吸收衰減模型,基于爆炸反射面模型推導了正Q濾波公式,利用正則化策略以及反演理論反演得到補償后的地震記錄,克服了常規(guī)反Q濾波方法振幅補償不穩(wěn)定或欠補償?shù)娜毕?且具有一定的抗噪性。Margrave等[15]將維納反褶積進行推廣,提出Gabor反褶積策略;Margrave等[16]在Q值未知以及反射系數(shù)白噪的前提下,對衰減地震記錄的Gabor譜進行雙曲平滑,統(tǒng)計衰減函數(shù)以及最小相位子波,最終得到高分辨率的反射系數(shù)序列。Reine等[17]基于窗口變化與否分析了時頻分析方法度量地震波衰減的穩(wěn)健性,提出變窗口的時頻變換方法可降低衰減估計的不確定性以及估計偏差,使其具有較高的精度,但該時頻分析方法均假設子波是最小相位的,一定程度上偏離實際情況。

我們借鑒反Q濾波思想僅進行吸收衰減補償,在Gabor域基于能量均分雙曲平滑方法估計衰減函數(shù)譜,基于非穩(wěn)態(tài)褶積理論[18],利用正則化方法反演得到衰減補償后的高分辨率地震記錄。

2 方法原理

2.1 非穩(wěn)態(tài)褶積理論

假設震源子波為w(t),反射系數(shù)序列為r(t),根據(jù)褶積原理,理想的無噪聲地震記錄s(t)的表述為

(1)

由于地下介質的非均勻性以及黏滯性,地震子波在傳播過程中被吸收、衰減產生頻散效應,實際上為非穩(wěn)態(tài)子波,因此基于穩(wěn)態(tài)褶積的處理結果往往偏離實際情況[16]。Margrave[18]將穩(wěn)態(tài)褶積模型推廣到非穩(wěn)態(tài)褶積情況,公式為

(2)

根據(jù)Gabor變換的性質知

(3)

其中,S(f,τ),R(f,τ),α(f,τ),W(f)分別為衰減地震記錄以及反射系數(shù)的Gabor譜、衰減函數(shù)、穩(wěn)態(tài)子波頻譜。在時頻域,可以較容易地估計出非穩(wěn)態(tài)子波的振幅譜,并在最小相位的假設下,借助希爾伯特變換求得相位譜,從而實現(xiàn)非穩(wěn)態(tài)反褶積。因此將非穩(wěn)態(tài)褶積公式(2)轉換到時頻域:

(4)

其中,S(f)為衰減記錄的頻譜,αw(f,τ)為衰減子波的時頻譜。

由于實際子波大多為混合相位子波,基于最小相位子波假設的非穩(wěn)態(tài)反褶積方法一定程度上偏離實際情況,因此我們忽略子波的影響,僅消除吸收衰減的效應,提高地震記錄的分辨率。非穩(wěn)態(tài)褶積公式(4)演變?yōu)?/p>

(5)

其中,α(f,τ)為復值衰減函數(shù),s0(τ)是原始未衰減的地震記錄。

公式(5)可離散成矩陣方程組的形式,即

(6)

其中,S為包含地震數(shù)據(jù)帶限頻率成分的復值向量,Φ為含有衰減以及傅里葉變換信息的復值矩陣,s0為原始未衰減記錄(補償后記錄),為一實值向量。

將公式(6)轉換為實值方程,即

(7)

(8)

其中,μ為正則化因子,D為微分算子,可選為單位算子、一階微分算子或二階微分算子。

將(8)式分別對s0求導并令之為0,即可獲得關于s0的正則化方程:

(9)

可利用共軛梯度算法、牛頓法等對正則化方程進行求解,最終得到衰減補償后的地震記錄。但是復值衰減函數(shù)α(f,τ)一般是未知的,如何得到精確的衰減函數(shù),成為衰減補償?shù)年P鍵問題。

2.2 衰減函數(shù)估計

衰減函數(shù)的理論表達式為

(10)

其中,Q為品質因子,H(·)為沿頻率方向的希爾伯特變換,其假設衰減函數(shù)為最小相位。如果已知Q值,則衰減函數(shù)具有解析表達式。但是精確的Q值估計相對困難且其誤差會影響衰減補償?shù)木?導致錯誤的解釋結果。因此,避開Q值估計問題,直接估計衰減函數(shù)可提高補償精度。令c=τf,根據(jù)衰減函數(shù)的理論表達式知,在c等值線上衰減值相同,因此在反射系數(shù)白噪的假設前提下,對衰減地震記錄的Gabor譜進行雙曲平滑可估計衰減函數(shù)的振幅譜;在最小相位假設下,利用希爾伯特變換可求得相應的相位譜,利用公式(9)反演得到衰減補償后的地震記錄。

不同的雙曲平滑方法對衰減函數(shù)的估計精度不同。傳統(tǒng)的雙曲平滑方法[16]將時頻域分成N個雙曲型條帶,條帶邊界為ci=(i-1)dc,i=2,…,N+1,c1=0,dc=cmax/N,如圖1所示。

圖1 傳統(tǒng)雙曲平滑方法條帶分割a N=10; b N=30

由圖1a可知,傳統(tǒng)的雙曲平滑分割方法對c值較小的區(qū)域欠采樣，對c值較大的區(qū)域過采樣,而地震信號的主要能量集中在c值較小的區(qū)域,使得估計的衰減函數(shù)譜精度較低。增加劃分條帶的數(shù)目,如圖1b所示,一定程度上可提高估計精度,但是光滑方法屬于統(tǒng)計方法,條帶數(shù)目增多,則每一條帶內的樣點數(shù)減少,統(tǒng)計效應減弱。

另外,以ci=exp[(i-1)dc],i=2,…,N+1,c1=0為條帶邊界,也可將時頻域分成N個雙曲型條帶,如圖2所示,其中dc=(lncmax)/N。得到的雙曲型條帶寬度從低c值到高c值逐漸增加,克服了傳統(tǒng)雙曲光滑方法在c值較小區(qū)域的欠采樣問題。c值較大區(qū)域的過采樣問題也得到了一定程度上的緩解，但是位于低c值條帶內的樣點數(shù)減少,統(tǒng)計效應減弱,且其只是運用數(shù)學方法機械地將條帶寬度進行調整,缺乏物理涵義。

圖2 對數(shù)域雙曲平滑方法條帶分割a N=10; b N=30

針對雙曲平滑方法遇到的問題,Li等[19]基于能量均分思想,提出變步長的雙曲平滑方法。由理論衰減函數(shù)振幅譜A=exp(-πτf/Q)=exp(-πc/Q)知,c在[0,cmax]變化時,A的取值范圍為(0,1],將能量區(qū)間進行均分，即dA=1/N,得到條帶邊界ci=-{ln[(i-1)dA]}Q/π,c1=0。品質因子Q未知,但由于平滑得到的估計值對Q值不敏感,所以Q取經驗常數(shù)值即可。能量均分雙曲平滑分割方法具有一定的物理涵義,其按照能量進行區(qū)域分割與疊前反演抽取近、中、遠偏移距道集類似。圖3展示了分割份數(shù)N=10,30時分別對應的c條帶分割結果。分割份數(shù)較少時,如圖3a所示,在保持統(tǒng)計效應的同時,克服了c值較小的區(qū)域欠采樣以及c值較大的區(qū)域過采樣問題。

圖3 能量平滑雙曲平滑方法條帶分割a N=10; b N=30

將時頻域分割成N個雙曲型條帶之后,在每一條帶內估計均值Ek,并將該平均值放置在雙曲條帶的中心線mk上:Ek=mean({|S(τ,f)|,(τ,f)∈Ωk}),其中Ωk為第k個雙曲型條帶,mk=0.5×(ck+ck+1)。

假設時頻域(τ,f)點位于[mk,mk+1],其中mk為第k個條帶的中心,由雙曲平滑方法得其時頻譜值Ek,則(τ,f)點處的時頻譜值可由插值方法獲得，即

(11)

3 數(shù)值模型算例與實際資料處理

本文方法為單道處理方法,為驗證方法的可行性，首先設計反射系數(shù)序列,分別與最小相位以及零相位子波褶積合成無衰減的地震記錄,利用非穩(wěn)態(tài)褶積理論合成衰減的地震記錄；然后基于傳統(tǒng)雙曲平滑以及能量均分雙曲平滑方法進行衰減函數(shù)估計,利用反演策略進行衰減補償。通過對實際資料進行處理,驗證了方法的有效性。

3.1 理論模型分析

設計的反射系數(shù)序列如圖4a所示;零相位子波為30Hz的雷克子波,如圖4b所示;最小相位子波為與零相位子波振幅譜相對應的最小相位子波,如圖4c所示。

基于褶積理論以及非穩(wěn)態(tài)褶積理論(品質因子Q=50)的最小相位子波合成記錄、零相位子波合成記錄及其理論衰減函數(shù)譜分別如圖5，圖6和圖7 所示。由圖5b，圖6b可知,隨著傳播時間的增加,子波的能量逐漸減弱,且高頻損失相對嚴重;由圖7可知,隨著時間的增加和頻率的提高,衰減量越來越大。

圖4 反射系數(shù)序列(a)與零相位子波(b)以及最小相位子波(c)

對最小相位子波合成的未衰減以及衰減地震記錄(圖5)進行Gabor譜分析,如圖8所示。由圖8a 可知,未衰減地震記錄Gabor譜的深、淺層能量相對均衡,由于帶限子波的影響,能量主要集中在中、低頻區(qū)域;由圖8b可知,衰減地震記錄的Gabor譜能量隨著傳播時間的增加逐漸減小,且高頻成分衰減相對較快,與理論衰減規(guī)律一致。

圖5 衰減前(a)、后(b)的最小相位子波合成記錄

圖6 衰減前(a)、后(b)的零相位子波合成記錄

圖7 理論衰減函數(shù)譜

圖8 最小相位子波合成的未衰減記錄Gabor譜(a)以及衰減記錄Gabor譜(b)

基于傳統(tǒng)的雙曲平滑方法,將時頻區(qū)域分割成N=10,30,50,100份得到的衰減函數(shù)譜如圖9a至圖9d所示?？梢钥闯?隨著分割份數(shù)的增加,估計得到的衰減譜與理論衰減譜(圖7)越來越接近;但是條帶數(shù)越多,每個條帶內的采樣點越少,統(tǒng)計效應越弱?；谀芰烤值碾p曲平滑方法(N=10)估計得到的衰減譜如圖10所示,該方法在分割份數(shù)較少的情況下與理論衰減模型(圖7)具有較好的一致性,在保持統(tǒng)計性的同時,具有較高的精度。由估計的衰減函數(shù)振幅譜在最小相位假設下利用希爾伯特變換得到相位譜,基于公式(9)反演得到補償后的地震記錄,如圖11所示。

圖9 傳統(tǒng)雙曲平滑方法不同分割份數(shù)估計的衰減譜a N=10; b N=30; c N=50; d N=100

圖10 能量均分雙曲平滑方法估計的衰減譜(N=10)

分析圖9至圖11可見,傳統(tǒng)的雙曲平滑方法估計的衰減譜的精度隨著分割份數(shù)的增加而提高,基于該衰減譜反演得到的補償記錄的精度也相應提高;基于能量均分的雙曲平滑方法在分割份數(shù)較少、保證統(tǒng)計性的前提下,估計得到的衰減譜(圖10)與理論衰減譜(圖7)具有較高的一致性,利用該衰減譜反演得到的補償后的地震記錄(圖11g)與未衰減地震記錄(圖11a)具有較高的一致性。綜合分析可知,基于能量均分的雙曲平滑方法相對傳統(tǒng)的雙曲平滑方法具有較大的優(yōu)勢,估計得到的衰減譜精度較高,反演得到的補償后地震記錄與原始未衰減記錄具有較好的一致性。

圖11 理論模型數(shù)據(jù)衰減補償前后的地震記錄a 原始未衰減地震道; b 品質因子Q=50時的衰減記錄; c—f 基于傳統(tǒng)雙曲平滑方法在分割份數(shù)分別為10，30，50，100時估計的衰減譜反演得到的補償后地震記錄; g 基于能量均分雙曲平滑方法估計的衰減譜反演得到的補償后地震記錄

將能量均分雙曲平滑方法以及反演策略應用于零相位子波合成的地震記錄(圖6)中,處理過程與處理最小相位子波合成記錄類似,補償后的地震記錄如圖12所示。由圖12可見,基于能量均分雙曲平滑方法以及反演策略得到的補償?shù)卣鹩涗浥c原始未衰減地震記錄具有較好的一致性,補償?shù)木容^高。理論模型數(shù)據(jù)試算的結果表明，我們研究并提出的方法對地震子波的相位沒有特殊要求,且不需要精確Q值的信息,衰減函數(shù)從地震記錄的Gabor譜中統(tǒng)計得到,利用反演的策略得到衰減補償后的地震數(shù)據(jù)避免了振幅補償?shù)牟环€(wěn)定性,實現(xiàn)了提高地震記錄分辨率的目的。由于忽略子波的影響,該方法具有數(shù)據(jù)自適應性。

3.2 實際地震資料處理

圖13a為某油田實際疊后地震資料,共有251道,每道651個采樣點,時間采樣率為2ms?？梢钥闯?由于吸收衰減作用,深層地震反射同相軸能量相對較弱,波形畸變且連續(xù)性較差(黑色橢圓以及箭頭所示位置)。利用本文方法,首先從地震記錄的Gabor譜中估計衰減函數(shù)譜,在衰減函數(shù)最小相位的假設下,利用希爾伯特變換求得相位譜,利用公式(9)實現(xiàn)地震記錄的衰減補償。補償后的地震剖面如圖13b所示,深層反射能量相對增強,波形得到一定程度的校正,且連續(xù)性變好,剖面分辨率得到提高。實際地震資料處理驗證了本文方法的有效性。

圖12 零相位子波合成記錄以及補償結果a 原始地震記錄; b Q=50對應的衰減地震記錄; c 補償后的地震記錄

圖13 實際地震資料衰減補償前(a)、后(b)的剖面對比

4 結論與認識

我們基于非穩(wěn)態(tài)地震道模型,利用能量均分雙曲平滑方法從衰減地震記錄的Gabor譜中估計衰減函數(shù),基于反演策略得到補償后的地震數(shù)據(jù)。研究得出以下幾點結論與認識:

1) 該方法忽略子波的影響,對子波的相位無要求,適用于任意相位子波合成的地震記錄,克服了Gabor反褶積中子波最小相位假設的缺陷;

2) 該方法不需要精確的Q值信息,衰減信息從觀測地震記錄的Gabor譜中提取,克服了常規(guī)反Q濾波需要精確Q值的弊端;

3) 該方法基于反演思想得到補償后的地震數(shù)據(jù),避免了常規(guī)反Q濾波振幅補償?shù)牟环€(wěn)定性;

4) 相對傳統(tǒng)雙曲平滑方法，基于能量均分雙曲平滑方法估計衰減譜的精度高,據(jù)此反演得到的衰減補償記錄分辨率更高;

5) 該方法沒有消除子波的影響,建議對補償后的地震記錄進行穩(wěn)態(tài)盲反褶積處理,同時得到子波信息以及高分辨率的反射系數(shù)序列;

6) 該方法為單道處理方法,沒有考慮橫向約束,有待拓展到多道處理,提高橫向的連續(xù)性;

7) 基于衰減函數(shù)最小相位假設和利用希爾伯特變換求取相位譜時會帶入計算誤差,下一步需要研究如何減弱最小相位假設和提高精度。

參 考 文 獻

[1] Futterman W I.Dispersive body waves[J].Journal of Geophysical Research,1962,67(13):5279-5291

[2] Varela C L,Rosa A L R,Ulrych T J.Modeling of attenuation and dispersion[J].Geophysics,1993,58(8):1167-1173

[3] Van der Baan M.Bandwidth enhancement:inverse Q filtering or time-varying Wiener deconvolution?[J].Geophysics,2012,77(4):V133-V142

[4] Wang Y.A stable and efficient approach of inverse Q filtering[J].Geophysics,2002,67(2):657-663

[5] Wang Y.Quantifying the effectiveness of stabilized inverse Q filtering[J].Geophysics,2003,68(1):337-345

[6] Wang Y.Inverse Q-filter for seismic resolution enhancement[J].Geophysics,2006,71(3):V51-V60

[7] Zhang X,Han L,Zhang F,et al.An inverse Q-filter algorithm based on stable wavefield continuation[J].Applied Geophysics,2007,4(4):263-270

[8] Hargreaves N,Calvert A.Inverse Q filtering by Fourier transform[J].Geophysics,1991,56(4):519-527

[9] 趙建勛,倪克森.串聯(lián)反Q濾波及其應用[J].石油地球物理勘探,1992,27(6):722-730

Zhao J X,Ni K S.Cascaded inverse Q filtering and the application[J].Oil Geophysical Prospecting,1992,27(6):722-730

[10] Yan H,Liu Y.Estimation of Q and inverse Q filtering for prestack reflected PP- and converted PS-waves[J].Applied Geophysics,2009,6(1):59-69

[11] Zhao Y,Liu Y,Li X,et al.An approach of inverse Q filtering considering time-frequency domain signal-to-noise ratio[J].74th EAGE Conference & Exhibition,2012,P103

[12] Braga I L S,Moraes F S.High-resolution gathers by inverse Q filtering in the wavelet domain[J].Geophysics,2013,78(2):V53-V61

[13] Zhang C,Ulrych T J.Seismic absorption compensation:a least squares inverse scheme[J].Geophysics,2007,72(6):R109-R114

[14] Wang S.Attenuation compensation method based on inversion[J].Applied Geophysics,2011,8(2):150-157

[15] Margrave G F,Dong L,Gibson P,et al.Gabor deconvolution:extending Wiener’s method to nonstationarity[R].Calgary:CREWES Research Report,2003

[16] Margrave G F,Lamoureux M P,Henley D C.Gabor deconvolution:estimating reflectivity by nonstationary deconvolution of seismic data[J].Geophysics,2011,76(3):W15-W30

[17] Reine C,van der Baan M,Clark R.The robustness of seismic attenuation measurements using fixed-and variable-window time-frequency transforms[J].Geophysics,2009,74(2):WA123-WA135

[18] Margrave G F.Theory of nonstationary linear filtering in the Fourier domain with application to time-variant filtering[J].Geophysics,1998,63(1):244-259

[19] Li F,Wang S,Chen X.Gabor deconvolution-hyperbolic smoothing with variable -step sampling[J].75thEAGE Conference & Exhibition,2013,Th-06-04