石油物探  2021, Vol. 60 Issue (5): 721-731  DOI: 10.3969/j.issn.1000-1441.2021.05.003
0
文章快速检索     高级检索

引用本文 

江雨濛, 曹思远, 陈思远, 等. 基于二阶自适应同步挤压S变换的时变子波提取方法[J]. 石油物探, 2021, 60(5): 721-731. DOI: 10.3969/j.issn.1000-1441.2021.05.003.
JIANG Yumeng, CAO Siyuan, CHEN Siyuan, et al. Time-varying wavelet estimation method based on second-order adaptive synchro-squeezing S transform[J]. Geophysical Prospecting for Petroleum, 2021, 60(5): 721-731. DOI: 10.3969/j.issn.1000-1441.2021.05.003.

基金项目

国家重点研发计划专项(2017YFB0202900)、国家自然科学基金(41674128)和国家重点研发计划(SQ2017YFGX030021)共同资助

第一作者简介

江雨濛(1994—), 女, 博士在读, 主要从事地震信号处理方法研究工作。Email: jymc16-2@163.com

文章历史

收稿日期:2020-09-08
改回日期:2020-12-10
基于二阶自适应同步挤压S变换的时变子波提取方法
江雨濛1,2, 曹思远1,2, 陈思远1,2, 马敏瑶1,2, 郑铎3, 黄芳4, 曹国明4    
1. 中国石油大学(北京)油气资源与探测国家重点实验室, 北京 102249;
2. 中国石油大学(北京)地球物理学院, 北京 102249;
3. 中国电子科技集团第二十九研究所, 四川成都 610036;
4. 中国石油大港油田公司, 天津 300280
摘要:鉴于时不变子波的假设条件不能反映实际地震数据的非稳态特性, 且传统分段提取时变子波的方法不能满足勘探所需的精度要求, 提出了一种基于二阶自适应同步挤压S变换(ASST2)的时变子波提取方法, 可以从非稳态地震数据中直接提取时变子波。首先对非稳态地震数据进行ASST2时频分析得到高质量的时频谱图, 再通过建立局部谱的统计参数与加权指数型子波的解析关系实现逐点提取子波, 最后基于时空域的时变子波矩阵建立非稳态地震记录褶积模型。与常规提取子波的方法相比, 该方法不需要地震数据稳态或分段稳态的假设条件, 根据数据自身特点自适应估计子波参数, 克服了对子波形态的限制条件, 同时避免了误差累计。模型数据测试和实际地震数据应用结果表明, 该方法能够有效地提取时变子波, 即使在相邻地层也能保证子波估计的准确性, 从而提高地震资料处理的精度。
关键词时变子波    非稳态数据    二阶自适应同步挤压S变换    加权指数型子波    统计参数    能量谱    高分辨率    
Time-varying wavelet estimation method based on second-order adaptive synchro-squeezing S transform
JIANG Yumeng1,2, CAO Siyuan1,2, CHEN Siyuan1,2, MA Minyao1,2, ZHENG Duo3, HUANG Fang4, CAO Guoming4    
1. State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum, Beijing 102249, China;
2. Institute of Geophysics, China University of Petroleum, Beijing 102249, China;
3. Southwest China Research Institute of Electronic Equipment, Chengdu 610036, China;
4. Dagang Oilfield Company, PetroChina, Tianjin 300280, China
Abstract: Conventional wavelet estimation methods involve the assumption that seismic data are stationary when a constant wavelet is considered; thus, they ignore the temporal variability of seismic wavelets.In reality, seismic data are nonstationary, and their frequency spectra change from shallow to deep formations.In this study, a time-varying wavelet extraction method was developed on the basis of the second-order adaptive synchro-squeezing S transform (ASST2).The ASST2 is a highly energy-concentrated time-frequency representation technique that allows for greater accuracy and stability of the proposed method.First, the ASST2 transform was carried out to obtain a high-quality time-spectrum image of the nonstationary seismic data; next, the wavelet was extracted point-by-point by establishing the analytic relationship between the statistical parameters of the local spectrum and the weighted exponential wavelet.Finally, the convolution model of the nonstationary seismic records was established based on the time-varying wavelet matrix in the time-space domain.The proposed method is suitable for handling the nonstationary nature of actual seismic data because it is a fully data-driven method; that is, the parameters of the time-varying wavelet are generated according to the nonstationary frequency spectrum.Moreover, synthetic tests demonstrate that the proposed method is reliable and robust, even under noise contamination.Finally, the application to field data confirms the superiority of the proposed method to conventional methods.
Keywords: time-varying wavelet    nonstationary seismic data    ASST2    FEW    statistical parameters    power spectrum    high resolution    

地震波在地下介质中传播时, 由于地层吸收衰减和噪声干扰等因素的影响, 子波高频成分能量衰减更快, 因此实际采集的地震数据往往是非稳态的, 且地震子波的频谱形态会随传播时间变化。由于地震子波是反射地震勘探的重要参数[1-2], 也是地震资料处理和解释的基础, 在很多环节中都起着关键作用, 因此地震子波提取方法的精度直接影响后续地震资料处理及解释的精度和准确性。传统的地震子波提取方法大多基于时不变子波的假设条件, 忽略了地震数据的非稳态特征, 使得子波估计结果存在误差, 进而限制了后续反演结果的分辨率。为了满足高精度地震勘探的需求, 如何从非稳态地震数据中准确提取时变子波成为近年来研究的热点之一[3-4]

针对上述问题, 前人提出了基于分段处理的时变子波估计方法[5-10], 即首先将地震记录划分为若干段, 并假设每一段内地震子波是稳态的, 进而结合不同方法提取每一段的时不变子波。此类方法依赖于分段方法和分段长度的选择, 且每段提取的平均意义上的子波与实际子波存在一定的误差, 使得后续反演结果不精确[11]。近年来, 为了消除时窗长度对提取子波的影响, 提高子波估计精度, 王蓉蓉等[12-13]和张漫漫等[14]相继提出基于时频谱模拟估计时变混合相位子波, 实现逐点提取子波。在此基础上, 李婧等[15]利用广义S变换改进了该方法并将提取的子波应用于地震反演。此类方法利用经验公式拟合时变子波振幅, 要求子波振幅谱满足类似雷克子波的假设条件, 限制了子波振幅谱的形态, 因此实际应用存在局限。姚振岸等[16]通过点谱模拟提取广义地震子波继而实现非稳态基追踪反演。ZHANG等[17]利用局部谱提取技术求取时变子波并应用于地震反演, 该方法基于自相关技术, 要求反射系数是白噪的, 但大多数反射系数并不符合这一假设条件。

时频分析技术能刻画时间与瞬时频率之间的非平稳关系, 是分析非稳态数据的重要工具之一[18-20]。S变换[21-23]将短时傅里叶变换和小波变换相结合, 既满足了多分辨率的特点, 又避免了小波变换的容许条件, 其实现过程具有较好的便捷性与适应性。同步挤压变换[24-26]是一种较新的具有较高分辨率的时频分析方法, 在原有时频分析基础上按照某种关系重排压缩, 获得更高精度更好聚焦性的时频谱图。TAO等[27]提出了一种基于二阶自适应同步挤压S变换(ASST2)的时频分析方法, 通过改进瞬时频率的计算公式实现高精度时频定位, 其时频分析结果能量聚焦性更强, 有利于精细刻画地震记录每一时刻的频率分量。该方法已经成功地应用于地震处理中的谱分解和面波压制环节。

本文提出了一种基于二阶自适应同步挤压S变换(ASST2)的时变子波提取方法, 该方法从非稳态地震数据中直接逐点提取子波, 避免了误差累计, 同时, 对于实际地震资料具有更好的适应性。首先利用ASST2技术对地震记录进行谱分解, 获取每一时刻的局部谱; 然后通过建立频谱统计参数与加权指数型子波(FWE)[28]之间的数学关系, 根据每一时刻的频谱信息得到对应的时变子波及其解析式, 而不是每一条地震道提取一个时不变子波; 最终得到基于时变子波矩阵的非稳态地震记录模型, 为后续反演提供理论基础。本文方法规避了传统子波估计方法要求地震记录稳态或分段稳态的假设条件, 从时频域逐点提取子波, 充分考虑了地震资料的非稳态特性。此外, 该方法是数据驱动的, 克服了常规方法对地震子波形态的限制条件, 同时精确的数学推导为本文方法提供了理论依据, 提高了估计子波的精度, 为后续实现高精度地震资料处理及解释提供了保障。

1 方法理论 1.1 ASST2时频分析方法

时频分析方法是研究非稳态信号的有效工具, 可以准确定位每一时刻的频率成分, 获取地震记录每一采样点的子波振幅谱, 进而实现逐点提取子波, 更精细地刻画地震记录的非稳态特性。

ASST2的数学表达式[27]为:

$ \operatorname{ASST} 2(t, f)=\int \operatorname{AST}(t, f) \mathrm{e}^{\mathrm{j} 2 \pi t f} \delta(f-\tilde{f}(t, f)) \mathrm{d} f $ (1)

式中: tf分别为时间和频率变量; AST(t, f)=$\mathrm{e}^{-\mathrm{j} 2 \pi t f} \int S(\alpha) \widetilde{W}(\alpha, f) \mathrm{e}^{\mathrm{j} 2 \pi a t} \mathrm{~d} \alpha, \alpha $为相对于t的傅里叶域变量, S(α)为信号频谱, $\widetilde{W} $为自适应时窗; δ为狄拉克函数; $\widetilde{f} $为瞬时频率。

同步挤压变换在实现过程中可以加入阈值来控制谱系数的利用, 在一定程度上可以降低原有时频表示对噪声的敏感程度。ASST2将自适应S变换与同步挤压变换相结合, 通过改进瞬时频率和自适应时窗, 实现时频分析结果能量聚焦并保证定位准确[27]。相较于基于短时傅里叶变换和小波变换的同步挤压变换时频分析方法, ASST2具有更高的时频分析精度, 即使在有噪声的情况下, 也能够有效地分辨出地震数据中频率成分的变化。此外, 该方法是数据驱动的, 可根据信号自身特点自适应调节窗口, 有利于分析具有动态衰减特性的非稳态地震数据。ASST2已经成功地应用于面波去除, 说明了该方法的可行性和有效性[27]

1.2 基于时频谱提取时变子波

对于一个任意振幅谱w(f), 统计参数(中心频率fm和方差σ2)是用来定量描述振幅谱形状及带宽的重要参数, 其定义为:

$\left\{\begin{array}{l} f_{m}=\frac{\int_{0}^{\infty} f w^{2}(f) \mathrm{d} f}{\int_{0}^{\infty} w^{2}(f) \mathrm{d} f} \\ \sigma^{2}=\frac{\int_{0}^{\infty}\left(f-f_{m}\right) w^{2}(f) \mathrm{d} f}{\int_{0}^{\infty} w^{2}(f) \mathrm{d} f} \end{array}\right. $ (2)

为了更好地表示实际子波特征, HU等[28]构建了加权指数型子波(FWE)表达式, 该式对参数要求相对简单, 实际应用相对灵活, 波形丰富可广泛用于多种形态的振幅谱[29]。加权指数型子波频率域解析式如下:

$ w(f)=a f^{n} \exp \left(\frac{-f}{f_{0}}\right) $ (3)

式中: a是振幅归一化系数; n是控制对称性的对称指数; f0是控制带宽的特征频率。将公式(3)代入公式(2), 建立FWE公式的理论参数(nf0)与振幅谱统计参数(fmσ2)之间的解析关系:

$ \left\{ {\begin{array}{*{20}{l}} {{f_m} = \frac{{2n + 1}}{2}{f_0}}\\ {{\sigma ^2} = \frac{{2n + 1}}{4}f_0^2} \end{array}\Rightarrow } \right.\left\{ \begin{array}{l} n = \frac{1}{2}\left( {\frac{{f_m^2}}{{{\sigma ^2}}} - 1} \right)\\ {f_0} = \frac{{2{\sigma ^2}}}{{{f_m}}} \end{array} \right. $ (4)

由公式(4)可知, FWE中的理论参数(nf0)可以由频谱的中心频率fm和方差σ2唯一确定, 也即, 对于任意频谱, 可以利用公式(4)和公式(3)进行定量表征, 再做反傅里叶变换得到最终时变子波。由于能量谱(即振幅谱的平方)对噪声具有一定压制作用, 因此, 公式(4)中通过能量谱得到的统计参数具有较好的抗噪性, 进而提高了子波估计的稳定性。同时, 公式(4)的推导过程没有进行近似, 避免了误差累计, 最大程度地保证子波提取的精度。

根据传统的褶积模型[30], 地震记录的矩阵形式可表示为地震子波矩阵W与反射系数向量的乘积, 即:

$ \boldsymbol{s}=\boldsymbol{W r}+\boldsymbol{n} $ (5)

式中: s, rn分别代表了地震记录、反射系数和噪声的向量; W是由时不变地震子波w(t)组成的托普利兹矩阵。这一模型建立在子波稳定的假设条件上, 即假设地震子波在传播过程中不会随传播距离的增加而变化。实际传播过程中, 由于地层吸收衰减效应, 地震子波会随着波的传播不断变化, 因此, 公式(5)不能准确描述实际地震数据。

考虑实际地震数据的非稳态特性, 利用ASST2可以得到地震记录每一时刻的频率成分, 即每一时刻的频谱, 根据公式(4)箭头左侧公式得到每一时刻频谱的统计参数(fmσ2), 再利用公式(4)箭头右侧统计参数(fmσ2)和加权指数型子波理论参数(f0n)之间的解析关系计算f0n, 最终由公式(3)唯一确定该时刻子波表达式。将提取的时变子波沿对角线排列, 并逐列取代W中的时不变子波w(t), 建立时变子波矩阵W*。由此, 可以将公式(5)中稳态褶积模型改写为基于时变矩阵的非稳态地震记录模型:

$ \begin{array}{*{20}{c}} \boldsymbol{s}=\boldsymbol{W}^{*} \boldsymbol{r}+\boldsymbol{n}\\ \boldsymbol{s}=\left[\begin{array}{c} s_{1} \\ s_{2} \\ \vdots \\ s_{m} \end{array}\right] \quad \boldsymbol{W}^{*}=\left[\begin{array}{llll} {\boldsymbol w_1^*}&{}&{}&{}\\ {}&{\boldsymbol w_j^*}&{}&{}\\ {}&{}& \ddots &{}\\ {}&{}&{}&{\boldsymbol w_m^*} \end{array}\right] \\ \boldsymbol{r}=\left[\begin{array}{c} r_{1} \\ r_{2} \\ \vdots \\ r_{m} \end{array}\right] \quad \boldsymbol{n}=\left[\begin{array}{c} n_{1} \\ n_{2} \\ \vdots \\ n_{m} \end{array}\right] \end{array} $ (6)

式中: m为采样点数。W*中的每一列代表了通过该列传播时间的子波wj*(j=1, 2, …, m), 而不是每一地震道提取一个时不变子波。虽然W*不再是托普利兹矩阵形式, 但它具有了代表地层非稳态特性的物理意义, 允许子波随时间变化而变化。因此, 相较于公式(5)中基于时不变子波矩阵的地震记录模型, 公式(6)中基于时变子波矩阵的非稳态地震记录模型更接近真实地震数据, 有利于构建后续反射系数反演的目标函数。

考虑到地震数据的非稳态特性, 本文提出的基于高精度时频分析的时变子波提取方法的处理流程如图 1所示, 根据公式(2)至公式(4)得到每一时刻子波的解析表达式, 精确的推导过程为该方法的可行性提供了理论依据, 同时提高了时变子波提取的准确性。值得注意的是, 本文方法实现了从非稳态地震记录中直接提取时变子波, 最大程度地避免了误差累计, 与常规方法相比, 本文方法摒弃分段平稳和子波振幅的假设条件, 提高了子波提取方法的适用性, 更有利于精细刻画实际地震资料的时变特征。

图 1 基于ASST2时频分析提取时变子波方法流程
2 模型数据测试

考虑地层吸收衰减因素的影响, 建立如图 2a所示的非稳态地震记录模型, 地震子波初始主频为40Hz, 随着时间增加子波主频不断下降。图 2a中蓝色实线为不含噪声情况下的合成地震记录(理论地震记录), 图 2b为对图 2a原始合成记录利用ASST2得到的时频谱图。由图可见, 地震记录的频谱从浅至深逐渐变化, 体现了地震记录的非稳态特性。同时, 基于ASST2得到的高精度时频分析结果能量聚焦性较好, 有利于精确提取地震记录的局部谱和实现时变子波的准确估计。从图 2b中分别提取每一时刻的频谱, 并通过公式(2)至公式(4)得到每一时刻对应的地震子波。图 2c图 2f分别为沿着图 2b中白色虚线提取的189, 192, 450, 800ms时刻时变子波振幅谱和时变子波(图中红色虚线所示)与相应理论子波振幅谱和理论子波(图中蓝色实线所示)的对比结果。如图中所示, 利用本文方法提取的子波不仅在频率域与理论振幅十分吻合, 反变换到时域亦是如此, 证明了方法的可行性和准确性。对比图 2c图 2d, 上下相邻地层之间虽然存在相互影响, 但本文方法仍然可以准确地提取每一时刻的子波, 为后续高精度反演提供保障。随着时间的增加, 提取的子波主频逐渐降低, 频带变窄, 相应时域波形变宽, 说明该方法真实地反映了地震子波在传播过程中的动态衰减特性。图 2a中的红色虚线为根据公式(6)利用所提取的时变子波构建时变子波矩阵W*与已知反射系数乘积重构的地震记录, 与理论地震记录吻合度较好。重构地震记录与理论地震记录的相关系数c=0.99, 进一步验证了本文方法所提子波的精确性。

图 2 基于ASST2时频分析提取时变子波 a 理论地震记录(蓝色实线)和重构地震记录(红色虚线); b ASST2时频谱; c 189ms时刻的子波振幅谱(上)和子波(下); d 192ms时刻的子波振幅谱(上)和子波(下); e 450ms时刻的子波振幅谱(上)和子波(下); f 800ms时刻的子波振幅谱(上)和子波(下)

为了测试该方法的抗噪声能力, 在上述非稳态地震记录模型中加入信噪比为8dB的随机噪声, 测试结果如图 3所示。图 3a中蓝色实线为含有噪声的非稳态地震记录(理论地震记录), 图 3b为对图 3a理论地震记录利用ASST2得到的时频谱图。对比图 2b图 3b的时频分析结果可以发现, 由于ASST2算法具有一定的抗噪声能力, 因此即使在有噪声的情况下, 仍然可以得到高分辨率、高精度的时频分析结果, 进而保证逐点提取子波的有效性。为了保证振幅谱统计参数(fmσ2)计算的准确性, 在信噪比较低时, 我们选取频谱的有效频带范围(0~100Hz)进行计算。同样, 从图 3b中分别提取每一时刻的子波, 图 3c图 3f分别为189, 192, 450, 800ms时刻提取的时变子波振幅谱和时变子波(图中红色虚线所示)与理论子波振幅谱和理论子波(图中蓝色实线所示) 对比结果。如图中所示, 在含有噪声的情况下, 利用本文方法所得时变子波和振幅都与理论子波基本吻合。对比图 2图 3可知, 所提取的时变子波并没有因为噪声的加入而发生剧烈变化, 波形保持了很好的稳定性, 这是因为本文方法利用能量谱(即频谱的平方)来估计时变子波的理论参数, 对噪声有一定的压制作用。此外, 利用所提取的时变子波构建时变子波矩阵与已知反射系数乘积重构的地震记录, 如图 3a中的红色虚线所示, 与理论地震记录仍具有较好的一致性, 其中, 重构地震记录与理论地震记录相关系数c=0.97。定性和定量分析发现, 即使是在含有中等噪声的情况下, 本文方法也可以提供准确可靠的时变子波估计结果, 由此表明该方法具有较好的抗噪能力和稳定性。

图 3 含噪声情况下基于ASST2时频分析提取时变子波 a 含噪声的理论地震记录(蓝色实线)和重构地震记录(红色虚线); b ASST2时频谱; c 189ms时刻的子波振幅谱(上)和子波(下); d 192ms时刻的子波振幅谱(上)和子波(下); e 450ms时刻的子波振幅谱(上)和子波(下); f 800ms时刻的子波振幅谱(上)和子波(下)

为验证本文方法的优越性, 利用图 2a所示地震记录对基于自相关逐点提取子波的方法[17]进行了处理测试, 结果如图 4所示。图 4a中蓝色实线为不含噪声情况下的合成地震记录(理论地震记录), 图 4b为利用短时傅里叶变换得到的时频谱图。对比图 2b图 4b可以很直观地感受到ASST2在对地震信号时频分析中所具有的优势, 其时频谱图品质更好、能量聚焦性更好, 由于图 2b中能量团集中在较小的时间范围内, 有助于时变子波的精准提取。图 4c图 4f分别为189, 192, 450, 800ms时刻提取的时变子波振幅谱和时变子波(图中红色虚线所示)与理论子波振幅谱和理论子波(图中蓝色实线所示)的对比结果。如图中所示, 基于自相关法逐点提取的子波与理论子波存在较大的差异, 且易受到相邻反射系数的影响(图 4c图 4d)。利用所提时变子波与反射系数重构地震记录的结果如图 4a红色虚线所示, 由于子波估计的不准确导致合成地震记录与理论值存在明显的误差, 幅值也不能得到很好地恢复, 其中, 重构地震记录与理论地震记录相关系数c=0.87。因此利用该子波进行后续反演处理, 其结果必然不准确。

图 4 基于自相关法提取时变子波 a 理论地震记录(蓝色实线)和重构地震记录(红色虚线); b 短时傅里叶时频谱; c 189ms时刻的子波振幅谱(上)和子波(下); d 192ms时刻的子波振幅谱(上)和子波(下); e 450ms时刻的子波振幅谱(上)和子波(下); f 800ms时刻的子波振幅谱(上)和子波(下)
3 实际资料应用

图 5显示了实际地震资料单道处理结果(第11道), 其中, 采样间隔为2ms, 截取的时窗范围为0.5~1.5s。分别应用自相关法提取时不变子波、基于自相关法提取时变子波[17]以及本文方法提取时变子波。对于实际地震资料而言, 无法直接从波形判断子波提取的准确性, 为此根据所得子波构建子波矩阵再利用L1范数稀疏约束正则化问题反演反射系数[31], 进一步对所提取的子波进行评价。图 5a为第11道地震记录, 图 5b图 5d中红色虚线分别为基于自相关方法提取时不变子波、基于自相关方法提取的时变子波和本文方法所提子波的单道反演结果。对比基于3种方法提取子波的反演结果和利用测井数据计算得到的反射系数(图 5b图 5d中蓝色实线所示)可见, 采用本文方法提取的时变子波的反演结果与基于测井数据计算的结果匹配性最好。分别计算反演结果与测井数据得到的反射系数的相关系数, 其中, 采用时不变子波时, 相关系数为0.39, 采用自相关法提取的时变子波时, 相关系数为0.42;而采用本文方法提取的时变子波时, 相关系数为0.56。对实际资料的处理结果与理论结果一致, 时变子波更符合地震信号的非稳态特性。对比图 5中反射系数反演结果可知, 子波提取的准确性直接影响反射系数结果的精度, 前两种方法的反褶积结果存在由于子波的误差引起的假反射系数, 利用本文方法处理后的噪声较弱, 同时相邻反射系数更加清晰、层位置更容易分辨(图 5d箭头所示)。

图 5 第11道地震记录(a)和基于自相关方法提取时不变子波(b)、自相关方法提取时变子波(c)、本文方法提取子波(d)的反演结果(红色虚线)以及基于测井资料计算的反射系数(蓝色实线)

图 6为分别利用上述3种方法提取的子波矩阵。基于自相关法提取的时不变子波(图 6a)其波形不随时间变化, 且存在一定的旁瓣, 基于自相关法提取的时变子波(图 6b)虽然其波形随传播时间逐渐变宽, 符合地下传播的动态衰减特性, 但存在较多的旁瓣(图中箭头所示), 会引起调谐效应。与前两种方法提取子波结果相比, 本文方法得到的时变子波(图 6c)不仅可以反映地震资料的时变特征, 同时其旁瓣能量较弱, 避免了由于旁瓣能量过强而引起的假象。

图 6 自相关法所提取的时不变子波矩阵(a)、自相关法所提取的时变子波矩阵(b)和本文方法所提取的时变子波矩阵(c)

图 7为上述3种方法处理后的振幅谱, 其中蓝色实线为处理前地震记录的振幅谱, 红色虚线为处理后的振幅谱。由图 7可知, 应用3种方法对实际地震资料处理后, 振幅谱的有效频带均得到拓宽, 而本文方法在保持地震资料信息的同时, 实现了高频、低频信息同时恢复。根据已有研究[32], 相对于高频信息, 低频信息的恢复难度更大, 而低频信息对相对频宽的贡献较大, 能较好地降低地震子波旁瓣效应, 且在地震反演中起到非常重要的作用。图 8为分别利用上述3种方法提取子波后反演的二维地震剖面, 道数为50, 可见, 与基于时不变子波处理后的剖面(图 8b)相比, 基于时变子波的反褶积后的剖面质量得到改善(图 8c图 8d), 分辨率较高。基于时变子波的反演结果可以分辨原始数据中不可分辨的薄层, 比如图 8a中单个同相轴的反射(图中箭头所指位置), 在基于本文方法的反演结果中呈现为两个反射同相轴。对比图 8c图 8d发现, 本文方法效果最佳, 不仅增强了同相轴能量(图中矩形框区域), 而且呈现了更多的薄层细节, 同时提高了反演结果的分辨率, 进一步验证了低频信息在地震反演中所起的重要作用。定性和定量分析可知, 相较于常规子波提取方法, 本文方法具有较好的可行性和优越性, 能够进一步提高地震资料反映薄层真实细节的能力, 有利于后续地震资料的精细解释。

图 7 处理前、后地震记录振幅谱 a 基于时不变子波处理结果; b 基于自相关法时变子波处理结果; c 基于本文方法处理结果
图 8 实际地震剖面(a)、基于时不变子波的反演剖面(b)、基于自相关法提取时变子波的反演剖面(c)以及基于本文方法的反演剖面(d)
4 结论

针对地震资料的非稳态特性, 本文提出了基于二阶自适应同步挤压S变换的时变子波提取方法, 实现了从地震数据中直接提取时变子波, 并建立了反映子波动态衰减特性的非稳态地震数据褶积模型。相较于常规子波提取方法, 本文方法通过推导局部谱的统计参数与加权指数型子波之间的数学关系, 准确得到时变子波的解析式, 最大程度避免了误差累计。模型数据测试结果证明了本文方法的可行性和有效性, 即使在信噪比较低的情况下, 也能够保证子波提取的精度。实际资料应用进一步验证了该方法能够适应地震资料的时变特性, 准确估计地震子波, 有利于提高地震资料分辨薄层细节的能力和后续地震资料解释的准确性。

本文方法是基于子波相位时不变提出的, 而实际地震资料中子波相位也会随传播时间变化, 如何将时变子波相位提取与本文方法相结合是我们下一步的研究方向。

参考文献
[1]
李福元, 韦成龙, 邓桂林, 等. 从海洋地震资料直达波提取震源子波[J]. 石油地球物理勘探, 2019, 54(3): 512-521.
LI F Y, WEI C L, DENG G L, et al. Source signature extraction from marine seismic direct waves[J]. Oil Geophysical Prospecting, 2019, 54(3): 512-521.
[2]
张鹏, 戴永寿, 谭永成, 等. 利用EMD和子波振幅谱与相位谱关系的时变子波提取方法[J]. 地球物理学报, 2019, 62(2): 680-696.
ZHANG P, DAI Y S, TAN Y C, et al. A time-varying wavelet extraction method using EMD and the relationship between wavelet amplitude and phase spectra[J]. Chinese Journal of Geophysics, 2019, 62(2): 680-696.
[3]
戴永寿, 张彧豪, 张鹏, 等. 时变地震子波提取研究方法综述[J]. 石油物探, 2020, 59(2): 169-176.
DAI Y S, ZHANG Y H, ZHANG P, et al. A review on time-varying seismic wavelet extraction[J]. Geophysical Prospecting for Petroleum, 2020, 59(2): 169-176. DOI:10.3969/j.issn.1000-1441.2020.02.002
[4]
王蓉蓉, 戴永寿, 张亚南, 等. 非平稳地震记录中时变子波提取方法研究[J]. 地球物理学进展, 2015, 30(2): 700-708.
WANG R R, DAI Y S, ZHANG Y N, et al. Time-varying wavelet extraction methods in non-stationary seismogram[J]. Process in Geophysics, 2015, 30(2): 700-708.
[5]
胡启宇. 同态反褶积的一种可能途径[J]. 石油物探, 1984, 23(2): 109-111.
HU Q Y. A possible way for hommorphic deconvolution[J]. Geophysical Prospecting for Petroleum, 1984, 23(2): 109-111.
[6]
冯晅, 刘财, 杨宝俊, 等. 分时窗提取地震子波及在合成地震记录中的应用[J]. 地球物理学进展, 2002, 17(1): 71-77.
FENG X, LIU C, YANG B J, et al. The extractive method of seismic wavelet in different time window and the application in synthetic seismogram[J]. Progress in Geophysics, 2002, 17(1): 71-77. DOI:10.3969/j.issn.1004-2903.2002.01.008
[7]
彭才, 朱仕军, 黄中玉, 等. 基于动态反褶及模型的动态子波估计[J]. 石油物探, 2007, 46(4): 224-328.
PENG C, ZHU S J, HUANG Z Y, et al. Dynamic wavelet estimation based on the dynamic convolution model[J]. Geophysical Prospecting for Petroleum, 2007, 46(4): 224-328.
[8]
VAN DER MIRKO B. Time-varying wavelet estimation and deconvolution by kurtosis maximization[J]. Geophysics, 2008, 73(2): V11-V18. DOI:10.1190/1.2831936
[9]
高静怀, 汪玲玲, 赵伟. 基于反射地震记录变子波模型提高地震记录分辨率[J]. 地球物理学报, 2009, 52(5): 1289-1300.
GAO J H, WANG L L, ZHAO W. Enhancing resolution of seismic traces based on the changing wavelet model of the seismogram[J]. Chinese Journal of Geophysics, 2009, 52(5): 1289-1300. DOI:10.3969/j.issn.0001-5733.2009.05.018
[10]
黄林军, 郭欣, 刘力辉, 等. 基于对时变子波进行分频段处理的反褶积方法在薄层预测中的应用[J]. 物探化探计算技术, 2019, 41(5): 586-592.
HUANG L J, GUO X, LIU L H, et al. A new deconvolution based on sub-band processing with the time-varying wavelet applied in thin reservoir prediction[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2019, 41(5): 586-592. DOI:10.3969/j.issn.1001-1749.2019.05.04
[11]
戴永寿, 张漫漫, 张亚南, 等. 基于时频谱模拟的时变混合相位子波提取[J]. 石油地球物理勘探, 2015, 50(5): 830-838.
DAI Y S, ZHANG M M, ZHANG Y N, et al. Time-variant mixed-phase seismic wavelet estimation based on spectral modeling in the time-frequency domain[J]. Oil Geophysical Prospecting, 2015, 50(5): 830-838.
[12]
王蓉蓉, 戴永寿, 李闯, 等. 时频分析与自适应分段相结合的时变子波提取方法[J]. 石油地球物理勘探, 2016, 51(5): 850-862.
WANG R R, DAI Y S, LI C, et al. Time-varying wavelet extraction based on time-frequency analysis and adaptive segmentation[J]. Oil Geophysical Prospecting, 2016, 51(5): 850-862.
[13]
王蓉蓉, 戴永寿, 李闯, 等. 基于奇异值分解的时变子波提取准确性评价方法[J]. 石油物探, 2015, 54(5): 531-540.
WANG R R, DAI Y S, LI C, et al. An evaluation criterion on the accuracy of time-varying wavelet extraction based on singular value decomposition[J]. Geophysical Prospecting for Petroleum, 2015, 54(5): 531-540. DOI:10.3969/j.issn.1000-1441.2015.05.006
[14]
张漫漫, 戴永寿, 张亚南, 等. 基于时频域谱模拟的时变子波估计方法[J]. 石油物探, 2014, 53(6): 675-682.
ZHANG M M, DAI Y S, ZHANG Y N, et al. Time-variant seismic wavelet estimation method based on spectral modeling in time-frequency domain[J]. Geophysical Prospecting for Petroleum, 2014, 53(6): 675-682. DOI:10.3969/j.issn.1000-1441.2014.06.007
[15]
李婧, 王修田, 宋鹏, 等. 基于改进广义S变换的时变反射系数反演[J]. 中国海洋大学学报(自然科学版), 2020, 50(6): 109-118.
LI J, WANG X T, SONG P, et al. Time-varying reflectivity inversion based on modified generalized S transform[J]. Periodical of Ocean University of China, 2020, 50(6): 109-118.
[16]
姚振岸, 孙成禹, 李红星, 等. 基于基追踪的时变子波提取与地震反射率反演[J]. 石油地球物理勘探, 2019, 54(1): 137-144.
YAO Z A, SUN C Y, LI H X, et al. Time-variant wavelet extraction and seismic reflectivity inversion based on basis pursuit[J]. Oil Geophysical Prospecting, 2019, 54(1): 137-144.
[17]
ZHANG R, FOMEL S. Time-variant wavelet extraction with a local-attribute-based time-frequency decomposition for seismic inversion[J]. Interpretation, 2017, 5(1): SC9-SC16. DOI:10.1190/INT-2016-0060.1
[18]
ALLEN J. Short term spectral analysis, synthesis, and modification by discrete Fourier transform[J]. IEEE Transactions on Acoustics, Speech, and Signal Processing, 1977, 25(3): 235-238. DOI:10.1109/TASSP.1977.1162950
[19]
SINHA S, ROUTH P S, ANNO P D, et al. Spectral decomposition of seismic data with continuous-wavelet transform[J]. Geophysics, 2005, 70(6): P19-P25. DOI:10.1190/1.2127113
[20]
张繁昌, 兰南英, 李传辉, 等. 地震匹配追踪技术与应用研究进展[J]. 石油物探, 2020, 59(4): 491-504.
ZHANG F C, LAN N Y, LI C H, et al. Geophysical Prospecting for Petroleum[J]. Geophysical Prospecting for Petroleum, 2020, 59(4): 491-504. DOI:10.3969/j.issn.1000-1441.2020.04.001
[21]
STOCKWELL R G, MANSINHA L, LOWE R P. Localization of the complex spectrum: The S transform[J]. IEEE Transactions on Signal Processing, 1996, 44(4): 998-1001. DOI:10.1109/78.492555
[22]
杨子鹏, 宋维琪, 刘军, 等. 联合广义S变换和压缩感知提高地震资料分辨率[EB/OL]. [2020-10-02]. http://kns.cnki.net/kcms/detail/11.2982.P.20200608.1427.130.html
YANG Z P, SONG W Q, LIU J, et al. Combine generalized s transform with compressed sensing to improve the resolution of seismic data[EB/OL]. [2020-10-02]. http://kns.cnki.net/kcms/detail/11.2982.P.20200608.1427.130.html
[23]
朱晓刚, 郑鹏飞, 王海萍, 等. 广义S变换在地震信号处理中的应用研究[J]. 西安文理学院学报(自然科学版), 2019, 22(3): 109-112.
ZHU X G, ZHENG P F, WANG H P, et al. Research on application of generalized S-transform in seismic signal processing[J]. Journal of Xi'an University(Natural Science Edition), 2019, 22(3): 109-112. DOI:10.3969/j.issn.1008-5564.2019.03.025
[24]
LIU W, CAO S Y, WANG Z M, et al. A novel approach for seismic time-frequency analysis based on high-order synchrosqueezing transform[J]. IEEE Geoscience and Remote Sensing Letters, 2018, 15(8): 1159-1163. DOI:10.1109/LGRS.2018.2829340
[25]
张琪, 刘彦萍. 基于同步挤压变换的地震勘探信号时频特性分析[J]. 石化技术, 2019, 26(8): 116-118.
ZHANG Q, LIU Y P. The time-frequency characteristic analysis of seismic prospecting signals based on synchrosqueezing transform[J]. Petrochemical Industry Technology, 2019, 26(8): 116-118. DOI:10.3969/j.issn.1006-0235.2019.08.072
[26]
刘晗, 张建中, 黄忠来. 基于同步挤压S变换的地震信号时频分析[J]. 石油地球物理勘探, 2017, 52(4): 689-695.
LIU H, ZHANG J Z, HUANG Z L. Time-frequency analysis of seismic data using synchro-squeezing S transform[J]. Oil Geophysical Prospecting, 2017, 52(4): 689-695.
[27]
TAO Y, CAO S Y, MA Y Y, et al. Second-Order adaptive synchrosqueezing S transform and its application in seismic ground roll attenuation[J]. IEEE Geoscience and Remote Sensing Letters, 2020, 17(8): 1308-1312. DOI:10.1109/LGRS.2019.2946368
[28]
HU W Y, JONATHAN L, BEAR L, et al. A robust and accurate seismic attenuation tomography algorithm[J]. Expand Abstracts of 81st SEG Annual Internat Mtg, 2011, 2727-2731.
[29]
崔琴, 于新霞, 潘龙, 等. 基于加权指数型子波的质心频移法估计品质因子[J]. 物探与化探, 2016, 40(6): 1203-1210.
CUI Q, YU X X, PAN L, et al. Q estimation by centroid frequency shift method based on frequency-weighted-exponential function[J]. Geophysical and Geochemical Exploration, 2016, 40(6): 1203-1210.
[30]
MARGRAVE G F, LAMOUREUX M P, HENLEY D C. Gabor deconvolution: Estimation reflectivity by nonstationary deconvolution of seismic data[J]. Geophysics, 2011, 76(3): W15-W30. DOI:10.1190/1.3560167
[31]
LI F Y, XIE R, SONG W Z, et al. Optimal seismic reflectivity inversion: Data-driven lp-loss-lq -regularization sparse regression[J]. IEEE Geoscience and Remote Sensing Letters, 2018, 16(5): 806-810.
[32]
曹思远, 袁殿. 高分辨地震资料处理技术综述[J]. 新疆石油地质, 2016, 37(1): 112-119.
CAO S Y, YUAN D. A review of high-resolution seismic data processing approaches[J]. Xinjiang Petroleum Geology, 2016, 37(1): 112-119.