2. 吉林大学, 长春 130026
2. Jilin University, Changchun 130026, China
地震谱分解技术可以将一个时域信号分解为时间和频率的二维函数,在储层刻画(Partyka et al.,1999)、烃类检测(Castagna et al.,2003; Han et al.,2012,2014)及频率依赖的AVO分析(Yoo and Gibson,2005)等地震解释领域得到了广泛的应用.
常用的时频谱分解方法有短时Fourier变换、Gabor变换、连续小波变换(CWT)(Mallat,1989; Sinha et al.,2005)、S变换(Pinnegar and Mansinha,2003)等. 这些方法在储层预测和烃类检测领域一直起着非常重要的作用,但同时也存在两方面的不足:一是时间分辨率低,难以对薄层进行精细刻画.当地层比较厚时,地层顶底反射可以通过常规谱分解方法识别,而当地层较薄时,常规谱分解方法很难识别顶底反射.事实上,这些基于变换的方法,不论使用固定时窗还是变尺度时窗,都受测不准原理限制,时间分辨率和频率分辨率此消彼长,即当确定频率分辨率的同时,就决定了最大时间分辨率难以突破(韩利,2013).二是在烃类检测中多解性强.含水储层和含烃类储层常常产生相似的频率响应,仅通过频率信息难以准确识别流体类型.
波在含流体介质中传播时,受到介质衰减和频散的影响(张固澜等,2014),不但会改变子波的频率,而且也会使子波相位发生改变.子波频率的变化特征已被用于储层预测,但子波相位的变化特征一直没有得到应用,其主要原因是现阶段对储层引起的相位响应特征还不够清楚,且缺乏从地震剖面中提取时变子波相位信息的技术手段.近些年,学者们不断提出高分辨率的谱分解方法(Wang,2007;朱振宇等2009;武国宁等;2012;Han and van derBaan,2013; Yuan and Wang,2013; Han et al.,2014),但仍然没有提取和利用相位信息的有效手段.针对这 个问题,本文将基于稀疏反演技术(Bonar and Sacchi,2010)的复谱分解方法用于烃类检测.本文方法具有两方面的优势,一方面提高了谱分解刻画薄储层的时间分辨率,另一方面可以求取用于气藏检测的时变子波相位信息.文章结构如下:首先阐述了复谱分解方法的原理,并用拟合例子展示该方法分辨薄层的能力和子波相位信息求取的准确性;然后通过数 值模拟方法对衰减引起相位改变的原因进行了分析;最后将本文提出的方法应用于实际资料的烃类检测.
2 复谱分解方法原理传统的地震褶积模型将地震记录 s(t)描述为一个固定频率、相位的地震子波 w(t)和反射系数 r(t)的褶积与随机噪声 n(t)之和:
其中,*表示褶积运算.但实际上地震子波受地下介质衰减和频散的影响,频率和相位都是随时间变化的.为了从地震记录中提取这些与岩性和流体有关的频率和相位信息,我们建立了一个比公式(1)更为复杂的复数褶积模型,
其中,w 不再是一个子波,而是一个由不同频率的零相位Ricker子波构成的复数子波库;K表示参与计算的子波的个数; w k 表示复子波库中第k个复子波,其频率和相位分别为f和φ w=0; r k 为对应的携 带有地震记录中子波频率和相位信息的复反射系数.
子波库的构成如图 1所示.图 1a展示了由不同频率的零相位Ricker子波构成的实子波库;通过对实子波库进行Hilbert变换可以得到复子波库,如图 1b所示,实线为复子波库实部,虚线为复子波库虚部.
将式(2)的褶积形式写成如下线性相乘形式:
其中,W 表示复子波库的褶积矩阵.式(3)中 r 的元素个数多于地震记录的元素个数,是一个欠定问题.文中通过如式(4)所示的L1范数进行正则化(Yuan et al.,2015),并使用快速迭代软阈值算法(FISTA)(Beck and Teboulle,2009)对式(4)进行求解,
式中,‖ r ‖1 为L1范数约束项,λ为稀疏度控制参数.
通过对 r 取模的平方和反正切得到地震记录的时频主频谱 F sdom和时频相位谱 φ s,
为了保持复谱分解的高时间分辨率优势,同时避免过高的频率分辨率带来的负面影响,可将时频主频谱乘上Ricker子波库的频谱转置得到时频能量谱,
其中 F Ricker为Ricker子波库频谱,T为转置运算; F s时频能量谱.
3 拟合算例图 2为一个展示复谱分解方法分辨薄层能力的 算例.图 2a为由不同时间间隔的30 Hz主频的Ricker 子波对组成的拟合信号,用来表征不同厚度地层的顶底反射.图 2b为由CWT方法得到的时频能量谱.当地层较厚时(60 ms间隔),CWT方法可以区分地层顶底反射;随着地层厚度变薄(<30 ms间隔),CWT时频谱不能很好地区分开顶底反射.图 2c和2d分别为复谱分解方法得到的时频主频谱和时频能量谱,可以看出复谱分解方法的时间分辨率很高.由于复谱分解的求解过程是一个迭代收敛的过程,稀疏度参数确定后,区分越薄的地层反演需要的迭代次数越多.在牺牲计算速度的情况下,复谱分解有能力将超薄层(<10 ms间隔)顶底反射识别出来.
图 3为一个由不同频率和相位的Ricker子波组成的拟合地震信号例子.拟合信号的子波频率从上到下依次为60,40,20,30 Hz和30 Hz,相位依次为0°,-90°,45°,-180°和-180°.图 3b为CWT方法得到的时频能量谱,代表着常规谱分解的分辨率水平.图 3c为CWT方法得到的时频相位谱,从中难以提取与衰减有关的子波相位信息.图 3d—3f分别为地震复谱分解方法得到的时频主频谱、时频能量谱和时频相位谱.地震复谱分解得到的三种时 频分布都具有很高的时间分辨率,并且所求结果与 拟合信号真实情况一致.拟合算例表明,相对于传统CWT方法,复谱分解方法的优势不仅在于具有高时频分辨率特征,而且还在于能够准确地计算出时变子波的相位信息.
地震波在能耗介质中传播会产生衰减已经成为一个共识,并且这种衰减会引起子波频率和相位的改变.介质衰减程度越大,引起的频率和相位变化量越大.在勘探中还发现另一个常被忽略的现象,即气藏的顶界面也会引起相位的改变.为此,本文利用基于Kelvin-Voigt模型(Morozov,2011)的黏弹波动方程对波在界面产生相位改变的原因进行了研究.在正演中对界面上下介质阻抗差异的影响和衰减差异的影响进行了独立分析,如图 4所示.层状模型中,界面的左侧介质较右侧介质波阻抗低、衰减程度大.图 4a和4b中蓝线表示仅考虑阻抗差异存在时产生的透射波和反射波,可以看到,透射波与入射波一致,仍然是零相位,而反射波相位发生180°改变,说明阻抗差异不会改变透射波的相位,仅引起反射波的极性改变.图 4a中的红线表示仅考虑衰减差异存在时产生的透射波和反射波,可以看出波在衰减介质中传播时相位发生了改变,且衰减程度的差异是引起反射波相位旋转的主要原因.图 4b中黑线表示同时考虑阻抗差异和衰减差异存在时产生的透射波和反射波.从模拟结果中我们可以得出,衰减引起的相位改变不仅发生在波传播过程中,还发生在存在衰减差异的界面处,储层界面反射是阻抗差异和衰减差异共同作用的结果.通常气藏与盖层有明显的阻抗差异和衰减差异,因此,通过相位信息检测气藏存在可能性.
图 5展示了一个实际资料烃类检测的算例.图 5a为目标区过探井的地震剖面,剖面上所投的测井曲线为含水饱和度曲线.根据探井揭示和地质解释的气层和水层位置已在图中标示.图 5b为CWT方法得到的用于烃类检测的频率异常剖面,时间分辨率低于地震分辨率,对薄层层位无法准确刻画;图 5c为使用本文提出的复谱分解方法得到的频率异常剖面,具有很高的时间分辨率,且与测井上揭示的地层位置一致.含气储层在频率异常剖面中有响应,但水层在频率异常剖面上也有响应,所以这会干扰气藏预测,造成解释的多解性.图 5d为复谱分解方法求取的时频相位谱异常剖面,从图中可以看到,气层和水层响应在相位异常剖面上存在一定差异,可用于降低解释的多解性.
本文提出了一种基于地震复谱分解技术的烃类检测方法.借助稀疏反演技术,复谱分解能够得到高分辨率的时频能量谱和时频相位谱.拟合算例证明了复谱分解不仅具有刻画薄层的能力,而且还能够准确求取时变子波相位信息.通过基于Kelvin-Voigt模型的黏弹波动方程数值模拟,证明了与衰减有关的相位变化不仅发生在波传播的过程中,也发生在存在阻抗和衰减差异的界面处.阻抗差异只引起反射子波的极性变化,而衰减差异是引起子波相位旋转的的关键因素.因此,可将子波相位变化用于具有高衰减特征的气藏检测.在实际资料烃类检测应用中,通过复分解方法求取的频率异常剖面具有很高的时间分辨率,并且进一步结合相位谱异常有效地区分了该区含气和含水储层的响应,降低了仅利用频率异常进行烃类检测的多解性.
致谢 感谢萨斯喀彻温省大学博士生邓武兵在黏弹波动方程数值模拟上做的贡献.[1] | Beck A, Teboulle M. 2009. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences, 2(1):183-202. |
[2] | Bonar D C, Sacchi M D. 2010. Complex spectral decomposition via inversion strategies.//SEG Technical Program Expanded Abstracts. SEG:1408-1412. |
[3] | Castagna J P, Sun S J, Siegfried R W. 2003. Instantaneous spectral analysis:Detection of low-frequency shadows associated with hydrocarbons. The Leading Edge, 22(2):120-127. |
[4] | Han J J, van der Baan M. 2013. Empirical mode decomposition for seismic time-frequency analysis. Geophysics, 78(2):O9-O19. |
[5] | Han L, Han L G, Li Z. 2012. Inverse spectral decomposition with the SPGL1 algorithm. Journal of Geophysics and Engineering, 9(4):423-427. |
[6] | Han L. 2013. Research on the methods of high-resolution full spectrum decomposition[Ph. D. thesis](in Chinese). Changchun:Jilin University. |
[7] | Han L, Sacchi M D, Han L G. 2014. Spectral decomposition and de-noising via time-frequency and space-wavenumber reassignment. Geophysical Prospecting, 62(2):244-257. |
[8] | Mallat S G. 1989. A theory for multiresolution signal decomposition:the wavelet representation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 11(7):674-693. |
[9] | Morozov I B. 2011. Anelastic acoustic impedance and the correspondence principle. Geophysical Prospecting, 59(1):24-34. |
[10] | Partyka G, Gridley J, Lopez J. 1999. Interpretational applications of spectral decomposition in reservoir characterization. The Leading Edge, 18(3):353-360. |
[11] | Pinnegar C R, Mansinha L. 2003. The S-transform with windows of arbitrary and varying shape. Geophysics, 68:381-385. |
[12] | Sinha S, Routh P S, Anno P D, et al. 2005. Spectral decomposition of seismic data with continuous-wavelet transform. Geophysics, 70(6):P19-P25. |
[13] | Wang Y H. 2007. Seismic time-frequency spectral decomposition by matching pursuit. Geophysics, 72(1):V13-V20. |
[14] | Wu G J, Cao S Y, Sun N. 2012. Matching pursuit method based on complex seismic traces and its application of hydrocarbon exploration. Chinese J. Geophys.(in Chinese), 55(6):2027-2034, doi:10.6038/j.issn.0001-5733.2012.06.024. |
[15] | Yoo S, Gibson R L Jr. 2005. Frequency dependent AVO analysis after target oriented stretch correction.// SEG:293-296SEG Technical Program Expanded Abstracts.. |
[16] | Yuan S Y, Wang S X. 2013. Spectral sparse Bayesian learning reflectivity inversion. Geophysical Prospecting, 61(4):735-746. |
[17] | Yuan S Y, Wang S X, Luo C M, et al. 2015. Simultaneous multitrace impedance inversion with transform-domain sparsity promotion. Geophysics, 80(2):R71-R80. |
[18] | Zhang G L, He Z H, Wang X M, et al. 2014. Seismic wave dispersion effects and inverse Q-filter phase compensation. Chinese J. Geophys.(in Chinese), 57(5):1655-1663, doi:10.6038/cjg20140528. |
[19] | Zhu Z Y, Lü D Y, Sang S Y, et al. 2009. Research of spectrum decomposition method based on physical wavelet transform and its application. Chinese J. Geophys.(in Chinese), 52(8):2152-2157, doi:10.3969/j.issn.0001-5733.2009.08.025. |
[20] | 韩利. 2013. 高分辨率全谱分解方法研究[博士学位论文]. 长春:吉林大学. |
[21] | 武国宁, 曹思远, 孙娜. 2012. 基于复数道地震记录的匹配追踪算法及其在储层预测中的应用.地球物理学报, 55(6):2027-2034, doi:10.6038/j.issn.0001-5733.2012.06.024. |
[22] | 张固澜, 贺振华, 王熙明等. 2014. 地震波频散效应与反Q滤波相位补偿. 地球物理学报, 57(5):1655-1663, doi:10.6038/cjg20140528. |
[23] | 朱振宇, 吕丁友, 桑淑云等. 2009. 基于物理小波的频谱分解方法及应用研究. 地球物理学报, 52(8):2152-2157, doi:10.3969/j.issn.0001-5733.2009.08.025. |