2. 中国石油大学(北京)海洋石油勘探国家工程实验室, 北京 102249
2. National Engineering Laboratory for Offshore Oil Exploration, China University of Petroleum, Beijing 102249, China
地震波在地下传播过程中,由于介质的非均质性导致地震信号的频率、振幅以及相位随时间变化而变化,即地震信号是一种非稳态信号(Wang et al., 2014; Wang, 2016).谱分解技术可将非稳态信号由单一的时间域转换到包含多种信息的时-频域,描述信号的多方面特征.传统时频分析方法通常被分为两大类,即线性变换方法和双线性变换方法.
谱分解的研究始于19世纪40年代,Gabor(2005)把解析信号(复信号)的重要概念引入信号分析中,提出了短时Fourier变换(STFT),第一次将时域信号转换到时-频域.Morlet等(1982)在1982年提出小波变换,使得非稳态地震信号分析有了飞速的发展,目前小波分析仍是地震勘探的重要分析手段之一(Sinha et al., 2005; Bi et al., 2014).Stockwell等(1996)在STFT和小波变换的基础上,提出了具有更高时频分辨率的S变换.随后,很多学者(Gao et al., 2003; Pinnegar and Mansinha, 2003; Huang et al., 2014; Wang, 2016)通过改变S变换中的窗函数特征,提出了广义S变换(GST).以上这些方法,均属于第一类时频分析方法,即线性变换方法.
双线性变换方法则主要是Wigner-Ville分布及其改进的算法.Ville在Wigner分布(Wigner, 1932)的基础上,提出了Wigner-Ville分布(WVD)(Ville, 1948),有效解决了STFT不能调节时窗所造成的时频聚焦性差的问题,但WVD交叉项问题也严重地限制了其应用.为了解决交叉项的干扰,学者们(Sattar and Salomonsson, 1999; Li and Zheng, 2001; de Souza Pereira Neto et al., 2001; Zhang, 2008; Wu and Liu, 2009; Li et al., 2010; Wang et al., 2012)提出了许多改进方法,如PWVD、SPWVD等,但这些方法在减小交叉项干扰的同时,也降低了谱分解的分辨率.Mallat和Zhang(1993)以及Chakraborty和Okaya(1995)将匹配追踪(MP)算法结合谱分解技术,得到了时频分辨率更高的MP-WVD方法.Wang(2007)则将匹配追踪时频分析技术应用于地震信号分析当中,取得了很好的效果.
压缩感知技术的提出使得L1范数求解问题的不断改进,基于稀疏约束的反演谱分解方法应运而生.Portniaguine和Castgna(2004)基于小波变换提出了反演谱分解方法(ISD-IWT).随后,很多学者将其应用于地震勘探领域中,取得了很好的应用效果.
在应用ISD-IWT分析信号时,通常采用Ricker子波经过Hilbert变换得到的解析子波作为复数小波.而利用小波变换进行信号分析时,应根据不同信号调整小波基的类型.因此,基于小波变换的ISD方法在实际信号分析中存在信号分析不准确的情况.
在前人研究的基础上,本文首先详细推导了基于离散小波变换(CWT)的反演谱分解(ISD)方法,并创新性地将其拓展到基于STFT和GST的ISD方法中.高时-频分辨率的ISD不再仅限于以往的时-频振幅谱研究,还可同时获得高分辨率的时-频相位谱,能为非稳态信号分析提供更多信息.分别采用常规谱分解方法和ISD方法对1-D合成信号进行分析并对比不同方法的振幅谱和相位谱结果.在此基础上,提出将ISD方法与FAVO反演相结合的储层流体识别的方法,并将其应用于实际数据,验证该方法的有效性.
2 反演谱分解基本原理反演谱分解是建立在常规谱分解方法可逆性的基础上,借鉴反演的思想和稀疏约束正则化策略,提出的一种具有高时频分辨率的谱分解方法.随着压缩感知理论的提出,L1范数求解算法飞速发展,这也使得稀疏约束的谱分解反演成为可能.
2.1 基于CWT的反演谱分解算法小波变换始于19世纪80年代,也是目前应用较为广泛的时频分析方法之一.信号f(t)的小波变换表达式为
(1) |
其中,a为尺度因子,b为平移因子,函数φ(t)为基本小波(母小波)应满足平方可积条件Cφ=
通过式(1), 我们就得到了小波变换的小波谱C(a, b),然后通过换算尺度因子与频率f的关系,即可得到时-频谱.
由于母小波函数的平方可积条件,所以小波变换存在逆变换,其表达式为
(2) |
通过观察可以发现,式(2) 可以表达为褶积形式和矩阵形式,即:
(3) |
(4) |
由此,小波反变换公式可以简化为
(5) |
其中,X为地震数据,W为母小波矩阵,C为待求的小波时-频谱.此时,式(5) 与常规地震反演所基于的线性正演公式一致,可采用相同的线性反演方法进行求解.为提高时-频谱的分辨率,以L1范数作为约束条件,采用SPGL1方法进行求解得到基于小波变换的反演谱分解结果.
2.2 基于其他变换的反演谱分解方法前面介绍了小波变换的反演谱分解方法,其实短时傅里叶变换(STFT)、Gabor变换、S变换以及广义S变换等方法均可以推导其相应的反演谱分解方法.如STFT的表达式为
(6) |
其中,f(t)为原始信号,g(t)表示窗函数.若g(t)高斯窗函数,则式(6) 可表示Gabor变换; 若g(t)=
而短时傅里叶变换是基于窗口的傅里叶变换,因此存在逆变换,其表达式为
(7) |
与小波变换类似,式(7) 同样可以表达成褶积形式和矩阵形式为:
(8) |
(9) |
同理,可以得到Gabor变换、S变换以及广义S变换的矩阵形式,均可结合反演的思想,以L1范数作为约束,采用SPGL1进行求解,最终得到基于STFT(或基于Gabor、ST、GST)的时频谱结果.
3 频变AVO基本原理根据Aki和Richard (1980)提出的AVO理论,平面波入射到界面处,反射波强度取决于入射角度和界面两侧介质的弹性参数.基于此推导出经典的Zoeppritz方程.
随着岩石物理的不断发展,许多学者逐渐发现地震反射系数不仅与入射角度和界面两侧介质弹性参数有关,也受到频率的影响.随后,频变AVO (FAVO)被提出.建立一个双层模型,上层介质为泥岩,下层介质为含气砂岩,采用斑块饱和模型对含气砂岩进行模拟.图 1为主频为25Hz的地震波到界面处的反射系数,颜色代表不同纵波反射系数值.如图 1所示,随着入射角度的增大,地震反射波反射系数逐渐发生频散现象,且入射角越大频散现象越严重.因此频变AVO效应(FAVO)在中、大角度更为明显.
由于储层流体的存在通常会产生更明显的FAVO效应,因此FAVO反演被提出并应用于流体的识别(见附录A).Wilson等(2009)将频率引入到Smith & Gidlow近似式中,基于谱分解结果提出了频变AVO(FAVO)反演方法,成功从地震资料中提取出对流体有很好指示作用的地震波频散属性信息.Wu等(2010)将平滑伪Wigner-Ville(SP-WVD)与FAVO反演相结合,并在实际资料中取得了良好的应用效果.Zhang等(2013)考虑流体因子f的频散因素,基于Russell近似提出了多尺度频变AVO反演方法.Guo等(2013)和Liu等(2014)分别利用频变AVO正演对非弹性层状介质和多尺度裂缝介质进行了流体的响应分析和预测.李坤等(2016)将匹配追踪的谱分解算法与FAVO反演相结合,提高了FAVO反演的分辨率.
以往的FAVO反演通常仅利用子波的振幅谱和地震数据的时频分析振幅谱来替代的子波矩阵W和地震记录S,表达式为
(10) |
而ISD方法不但可以获得高时频分辨率的振幅谱信息,还能获得高时频分辨率的相位谱信息.因此,基于ISD的FAVO反演能够利用振幅谱和相位谱的完整信息,保证反演输入数据的完整性,增强反演结果的准确性.由于流体会引起纵波较明显的FAVO现象而对横波并不敏感,因此通常利用频变AVO反演得到的纵波频散项
对常用的谱分解方法进行了简单的分类(见图 2):
(1) 短时傅里叶变换及近似变换,当STFT变换中的窗函数为高斯函数时,即为Gabor变换,因此Gabor变换与STFT同属一类.
(2) S变换及改进方法,GST是基于S变换的改进方法,因此GST与S变换同属一类.
(3) WVD及改进方法,MP-WVD通过匹配追踪方法(MP)将信号分解为多个子波原子,再利用WVD对子波原子进行分析,从而达到消除交叉项的目的,是基于WVD的一种改进算法,因此与WVD同属一类.
(4) 小波变换(CWT).
(5) ISD方法.
基于以上分类,选取STFT、CWT、GST、MP-WVD以及ISD方法分别对复杂合成信号进行频谱分解,并对比分析振幅谱结果.测试信号一个由11个不同主频子波叠加合成的复杂信号,如图 3所示.
图 4是利用五种谱分解方法得到的振幅谱,通过对比分析可知:
(1) 在时-频分辨率方面,MP-WVD和ISD要明显优于其他三种线性谱分解方法.
(2) 在结果准确性方面,MP-WVD的分析结果与原始合成信号的时频关系存在一定误差,分析错误的位置已在图 4e中用红框标出.这是由于MP算法在匹配合成信号时对于复合波的重构出现错误所引起的.此外,MP方法的准确性受限于迭代次数,迭代次数不足会导致细节被忽略,增大分析误差.
(3) 在计算效率方面,STFT、CWT和GST具有较快的计算速率,而MP-WVD需要对多个参数(时移、频率、相位以及尺度因子)进行全局最优搜索,因此计算速率最低.ISD方法的计算速率虽大于STFT等方法,但却远小于MP-WVD算法.
综上所述,ISD方法具有更高的时频分辨率,能使非平稳信号能量在时-频谱上更好聚焦.此外,ISD方法计算速率适中,分析结果具有较高准确度和可信度.合成一个带有相位信息的信号,利用STFT、GST与ISD方法得到信号相位谱,如图 5所示.
图 5a为测试合成信号,由5个不同主频和相位的子波所组成.子波中心时刻分别对应0.05、0.15、0.27、0.375 s和0.425 s,主频自上而下分别为55、50、30、40 Hz和40 Hz,相位依次为30°、-90°、-180°、0°和60°.图 5b为瞬时相位属性,可以看出瞬时相位属性能较为准确地指示子波中心位置相位值,但不能反映波形的信息,即没有波形时仍存在相位信息,因此利用瞬时相位属性指示信号相位信息存在局限性.
图 5c—e分别为由STFT、GST和ISD方法得到的时频相位谱.由于STFT和GST分辨率不足,能量不聚焦,导致在不应该存在能量的位置,仍存在相位的变化,这显然会给相位的分析带来影响.而高时-频分辨率的ISD方法,使得能量聚焦,从而得到更加准确的相位信息.
FAVO反演是基于地震数据的时频分析结果完成的,因此时频分析结果的分辨率、准确性是影响FAVO反演结果的关键因素.ISD方法不但能提供准确的高时-频分辨率的振幅谱信息,还可提供高时-频分辨率的相位谱信息,保证了时频分析结果的分辨率、准确性和完整度,也为更好地利用FAVO反演进行油气预测奠定了基础.
5 实际数据应用基于前面的分析认识,本文提出基于ISD的FAVO反演(ISD-FAVO)方法,并将该方法分别应用于中国东部某油田的实际井旁地震道和叠前地震数据中,以验证ISD方法的稳定性和ISD-FAVO反演方法的可行性.研究区位于中国东部某油田,目的层段为砂泥岩互层,共发育三套连续的含油储层.
5.1 井旁地震数据测试选取井旁叠前地震道数据进行测试分析,首先利用ISD方法对叠前部分叠加数据进行时频分析,得到不同角度地震数据的谱分解结果.图 6a—c分别为井旁5°、15°和25°叠前部分叠加数据的反演谱分解结果.
基于反演谱分解结果,开展FAVO反演,得到纵波频散数据,见图 7.图 7a为FAVO反演得到的纵波频散属性,图 7b—f分别为自然伽马GR、纵波速度VP、横波速度VS、泥质体积VCL和含水饱和度Sw曲线.结合测井曲线有助于验证反演结果的准确性.测井中通常利用伽马曲线指示泥质含量,低伽马值指示低泥质含量的砂岩储层;砂岩储层中烃类的存在会引起纵波频散现象,因此砂岩储层应表现为高纵波频散特征;图 7中低伽马与高纵波频散部分相对应,均指示了砂岩储层的位置,与理论认识相吻合,验证了ISD-FAVO方法的有效性.
选用同一研究区的叠前地震数据对ISD-FAVO方法进行测试,图 8为近、中、远偏移距的叠前部分叠加数据,入射角度分别为5°、15°及25°.分别采用STFT、GST以及ISD方法进行谱分解分析,并在此基础上进行FAVO反演.
由前面的分析可知,ISD方法相较于文中提到的其他常用谱分解方法,不仅能获得更高时-频分辨率的振幅谱,还可以同时获得可靠的时-频相位谱信息.与传统FAVO反演不同之处在于,本文不再仅使用时-频振幅谱,而是采用时频振幅谱与相位谱相结合的完整信息开展反演工作,提高反演结果的可靠性.
图 9a—c分别为利用STFT、GST以及ISD的结果进行FAVO反演得到的纵波频散属性Ivp,图中黑色曲线为测井解释得到的含水饱和度曲线,饱和度曲线的右端对应含水饱和度的高值.目的层发育的三套砂岩储层应表现为低含水饱和度的特征.对比图 9,结果表明:
(1) 基于STFT得到的纵波频散结果时间分辨率最低,三套砂岩储层混叠无法区分,不能提供储层的边界信息,储层预测结果(高纵波频散区域)与测井解释结果(低含水饱和度值)对应不佳.
(2) GST结果的分辨率高于STFT,但明显低于ISD,能简单区分三套砂岩储层,能够粗略提供储层边界信息,储层预测结果与测井解释结果的对应性好于STFT.
(3) ISD属性结果具有最高的时间分辨率,能清晰区分三套储层,缩小储层预测范围,降低预测的不确定性,提供更为真实可信的储层边界信息,储层预测结果能更好地与测井解释结果相对应.
由此可以得出,ISD-FAVO反演效果明显好于基于常规时频分析方法的FAVO反演,具有更高的时间分辨率,降低储层预测的不确定性,提供更为可信的边界信息,这也说明了ISD-FAVO反演方法在烃类识别具有更大的优势和潜能.
6 结论本文提出了基于ISD的FAVO反演技术,得到了较好的应用效果,同时也得出以下结论:
(1) 通过模型测试与实际数据应用发现ISD方法的时-频分辨率明显高于常规线性转换方法(以STFT、GST、CWT为例),相比于双线性转换(以MP-WVD为例)具有更高的结果准确度和较低的运算时间.
(2) FAVO技术能将由于储层烃类存在所引起的纵横波速度频散变化率信息从地震数据中提取出来,用以实现储层烃类的识别.FAVO反演需基于时频分析结果进行,时频分析结果的分辨率和准确性对反演结果影响明显.ISD-FAVO反演的结果具有更高的时间分辨率,能降低储层预测的不确定性,提供更为可信的边界信息.
(3) 基于稀疏约束的反演谱分解方法,不仅能提供高时频分辨率的振幅谱,还能为非稳态信号分析提供可靠的时-频相位谱信息,这能为利用相位信息进行地下构造、岩性分析提供更可靠的信息.
附录A
Smith和Gidlow(1987)将Gardner G H F和Gardner L W (1974)给出的含水岩石的密度与速度关系
(A1) |
Wilson等(2009)在此基础上,将式(A1) 中的
fdom表示地震信号某一时段的主频,
(A2) |
其中:
(A3) |
取一阶近似,可将式(A1) 表示为
(A4) |
在不考虑频散情况下,纵波反射系数为
(A5) |
令:
(A6) |
则频散情况和无频散情况的纵波反射系数的关系为
(A7) |
通过褶积公式可知:
(A8) |
(A9) |
对比式(A8) 与式(A9) 发现,通过约去相同项即可得到频散项
(A10) |
其矩阵形式表示为
(A11) |
感谢外审专家对本文的审阅.
Aki K, Richards P G. 1980. Quantitative Seismology:Theory and Methods (New York:Freeman).
|
|
Bi J W, Wu Z J, Wang Z J, et al. 2014. Wavelet transform and its application in earthquake engineering.//2014 5th International Conference on Intelligent Systems Design & Engineering Applications. Hu'nan, China:IEEE, 1126-1128.
|
|
Chakraborty A, Okaya D.
1995. Frequency-time decomposition of seismic data using wavelet-based methods. Geophysics, 60(6): 1906-1916.
DOI:10.1190/1.1443922 |
|
Gabor D. 2005. Theory of communication. New York:Dover Publications, 58.
|
|
Gao J H, Chen Y C, Li Y M, et al.
2003. Generalized S transform and seismic response analysis of thin interbedss surrounding regions by GPS. Chinese Journal of Geophysics, 46(4): 759-768.
DOI:10.1002/cjg2.v46.4 |
|
Gardner G H F, Gardner L W.
1974. Gregory A R. Formation velocity and density; the diagnostic basics for stratigraphic traps. Geophysics, 39(6): 770-780.
|
|
Guo Z Q, Li X Y, Liu C, et al.
2013. A shale rock physics model for analysis of brittleness index, mineralogy and porosity in the Barnett Shale. Journal of Geophysics & Engineering, 10(2): 025006.
|
|
Huang H, Feng N, Wang Y, et al.
2014. High-resolution seismic processing based on generalized S transform. Oil Geophysical Prospecting, 49(1): 82-88.
|
|
Li K, Yin X Y, Zong Z Y.
2016. Time-frequency-domain FAVO fluid discrimination method based on matching pursuit spectrum decomposition. Acta Petrolei Sinica , 37(6): 777-786.
|
|
Li Y D, Li J S, Zheng X D. 2010. Channel system characterization using Wigner-Ville distribution-based spectral decomposition.//SEG Technical Program Expanded. SEG, 1418-1422.
|
|
Li Y D, Zheng X D.
2008. Spectral decomposition using Wigner-Ville distribution with applications to carbonate reservoir characterization. The Leading Edge, 27(8): 1050-1057.
DOI:10.1190/1.2967559 |
|
Liu C, Li B N, Zhao X, et al.
2014. Fluid identification based on frequency-dependent AVO attribute inversion in multi-scale fracture media. Applied Geophysics, 11(4): 384-394.
DOI:10.1007/s11770-014-0454-0 |
|
Mallat S G, Zhang Z.
1993. Matching pursuits with time-frequency dictionaries. IEEE Transactions on Signal Processing, 41(12): 3397-3415.
DOI:10.1109/78.258082 |
|
Morlet J, Arens G, Fourgeau E, et al.
1982. Wave propagation and sampling theory-Part Ⅰ:Complex signal and scattering in multilayered media. Geophysics, 47(2): 203-221.
DOI:10.1190/1.1441328 |
|
Pereira de Souza Neto E, Custaud M A, Frutoso J, et al.
2001. Smoothed pseudo Wigner-Ville distribution as an alternative to Fourier transform in rats. Autonomic Neuroscience:Basic & Clinical, 87(2-3): 258-267.
|
|
Pinnegar C R, Mansinha L.
2003. The S-transform with windows of arbitrary and varying shape. Geophysics, 68(1): 381-385.
DOI:10.1190/1.1543223 |
|
Portniaguine O, Castagna J. 2004. Inverse spectral decomposition.//SEG Technical Program Expanded Abstracts. SEG, 1786-1789.
|
|
Sattar F, Salomonsson G.
1999. The use of a filter bank and the Wigner-Ville distribution for time-frequency representation. IEEE Transactions on Signal Processing, 47(6): 1776-1783.
DOI:10.1109/78.765169 |
|
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.
DOI:10.1190/1.2127113 |
|
Smith G C, Gidlow P M.
1987. Weighted stacking for rock property estimation and detection of gas. Geophysical Prospecting, 35(9): 993-1014.
DOI:10.1111/gpr.1987.35.issue-9 |
|
Stockwell R G, Mansinha L, Lowe R P.
1996. Localization of the complex spectrum:the S transform. IEEE Transactions on Signal Processing, 44(4): 998-1001.
DOI:10.1109/78.492555 |
|
Ville J.
1948. Theorie et applications de la notion de signal analytique. Cables et Transmissions, 2(1): 61-74.
|
|
Wang B F, Chen X H, Li J Y. 2014. Inversion based data-driven attenuation compensation method.//SEG Technical Program Expanded Abstracts 2014. SEG, 3267-3271.
|
|
Wang B F.
2016. An amplitude preserving S-transform for seismic data attenuation compensation. IEEE Signal Processing Letters, 23(9): 1155-1159.
DOI:10.1109/LSP.2016.2586445 |
|
Wang Y H.
2007. Seismic time-frequency spectral decomposition by matching pursuit. Geophysics, 72(1): V13-V20.
DOI:10.1190/1.2387109 |
|
Wang Z W, Wang X L, Xiang M, et al.
2012. Reservoir information extraction using a fractional fourier transform and a smooth pseudo Wigner-Ville distribution. Applied Geophysics, 9(4): 391-400.
DOI:10.1007/s11770-012-0352-2 |
|
Wigner E P.
1932. On the quantum correction for thermodynamic equilibrium. Physical Review, 40(5): 749-759.
DOI:10.1103/PhysRev.40.749 |
|
Wilson A, Chapman M, Li X Y. 2009. Frequency-dependent AVO inversion.//SEG Technical Program Expanded Abstracts. SEG, 341-345.
|
|
Wu X Y, Chapman M, Wilson A, et al. 2010. Estimating seismic dispersion from pre-stack data using frequency-dependent AVO inversion.//SEG Technical Program Expanded. SEG, 425-429.
|
|
Wu X Y, Liu T Y.
2009. Spectral decomposition of seismic data with reassigned smoothed pseudo Wigner-Ville distribution. Journal of Applied Geophysics, 68(3): 386-393.
DOI:10.1016/j.jappgeo.2009.03.004 |
|
Zhang R. 2008. Spectral decomposition of seismic data via wigner-ville distribution and chirplet transform.//70th EAGE Conference & Exhibition. SPE, EAGE.
|
|
Zhang Z, Yin X Y, Zong Z Y. 2013. A new frequency-dependent AVO attribute and its application in fluid identification.//Proceedings of the 75th EAGE Conference & Exhibition. London:EAGE.
|
|
李坤, 印兴耀, 宗兆云.
2016. 基于匹配追踪谱分解的时频域FAVO流体识别方法. 石油学报, 37(6): 777–786.
|
|