基于二阶自适应同步挤压S变换的时变子波提取方法

江雨濛,曹思远,陈思远,马敏瑶,郑 铎,黄 芳,曹国明

(1.中国石油大学(北京)油气资源与探测国家重点实验室,北京102249;2.中国石油大学(北京)地球物理学院,北京102249;3.中国电子科技集团第二十九研究所,四川成都610036;4.中国石油大港油田公司,天津300280)

地震波在地下介质中传播时,由于地层吸收衰减和噪声干扰等因素的影响,子波高频成分能量衰减更快,因此实际采集的地震数据往往是非稳态的,且地震子波的频谱形态会随传播时间变化。由于地震子波是反射地震勘探的重要参数[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]为:

(1)

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

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

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

(2)

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

(3)

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

(4)

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

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

s=Wr+n

(5)

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

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

s=W*r+n

(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,进一步验证了本文方法所提子波的精确性。

为了测试该方法的抗噪声能力,在上述非稳态地震记录模型中加入信噪比为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。定性和定量分析发现,即使是在含有中等噪声的情况下,本文方法也可以提供准确可靠的时变子波估计结果,由此表明该方法具有较好的抗噪能力和稳定性。

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

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

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