基于改进广义S变换二次时频谱反褶积在南海工区的应用

周 慰， 黄 饶， 周怀来,2， 牛 聪， 巫南克

(1.成都理工大学 地球物理学院，成都 610059；2.“油气地质及开发工程”国家重点实验室，成都 610059；3.中海油研究总院，北京 100028)

0 引言

传统的反褶积基于地震波在地下传播的过程中能量没有发生变化的前提,但地震波在实际的传播过程中主频会随着传播时间的增长向低频偏移，传统的褶积模型不能有效模拟这个动态过程[1-2]。

考虑到地下介质的衰减作用，彭才[3]假设子波是最小相位，利用非平稳地震记录提取时变子波；Margrave[4]将衰减地震记录转换到Gabor域，结合Rosa提出的平滑函数，对每个时窗振幅谱平滑，直接估计出衰减子波振幅谱，提出了一种时频域褶积模型；王蓉蓉[5]对比各种时频分析方法，对时变子波提取方法进行研究并总结。反褶积的关键是能否更好的分离子波和反射系数，在频率域中分析关键在于能否更准确地提取子波振幅谱。Rosa[6]提出了直接对地震记录的振幅谱进行谱模拟，提取子波振幅谱，实现了谱模拟反褶积；周怀来[7-9]在S域反褶积，结合了Rosa谱模拟提取子波是思想，对子波时频谱进行谱模拟，提取动态子波时频谱，实现S域反褶积。

地震信号的时频谱分析方法有Gabor变换、小波变换、S变换等。其中Gabor变换受时窗的影响，不能根据实际信号的频率变化进行灵活处理；小波变换需要对小波基进行合理选择，高频区分辨率差；S变换基本函数固定，高频信号识别能力低。笔者对S变换窗函数进行扩展，提高其对高频信号的分辨率，推导出一种改进广义S变换，并在此基础上使用傅里叶变换得到非平稳地震记录二次时频谱，通过低通滤波提取子波时频谱，实现二次时频谱反褶积，弥补了传统平滑函数拟合阶数不确定性的缺陷，摆脱了反射系数满足白谱特性的假设，提高了算法的实用性，从而达到在不影响地震资料信噪比的前提下，恢复地震波被地层吸收所引起的衰减，并提高对薄储层地识别能力。

1 方法原理

1.1 改进广义S变换

由于S变换基本窗函数基本形态固定，不能根据实际情况调节窗口的大小，近年来广大专家学者提出很多改进广义S变换，以达到提高时频聚焦性的目的。为了弥补S变换在信号高频部分分辨率低的缺陷，这里将其高斯窗函数进行扩展，令σ=λ/|f|r，使得当信号处于高频时，加快窗口的收缩速度[10]。

(1)

式中：t为时间;f为频率;λ、r用来调节窗函数的宽度和高度。

因此得到对给定信号x(t)的一种两参数改进广义S变换为：

exp(-i2πt)dt

(2)

1.2 二次谱反褶积模型

由于地震波在地下传播的过程中会受到地下介质的影响，主频向低频偏移。为了更加符合这一实际情况，笔者引入一个时间延迟因子τ和衰减函数来描述这个衰减的动态过程，推导出一个频率域的衰减褶积模型[11]：

(3)

其中：ω(f)是震源子波的振幅谱；r(t)是时间域的反射系数；α(τ，f)用来描述地震子波在频率域中的衰减，表示为式(4)。

a(τ，f)=exp(πfτ/Q)

(4)

将式(4)扩展到S域中，可以将衰减地震道的时频谱G(τ，f)近似等同子波频谱ω(f)、衰减函数a(τ，f)和反射系数改进广义S变换后时频谱R(τ，f)之间的乘积，如果只从振幅谱出发则可以表示成式(5)。

|G(τ，f)|≈|ω(f)‖a(τ，f)‖R(τ，f)|

(5)

非平稳地震记录的动态子波振幅谱W(τ，f)可以用其静态子波振幅谱ω(f)与衰减函数振幅谱a(τ，f)的乘积表示，即：

|G(τ，f)|≈|W(τ，f)‖R(τ，f)|

(6)

基于Rosa提出的传统的子波提取函数为式(7)[12]。

(7)

由于式(7)谱模拟时受数据体干扰大，平滑函数受拟合阶数k、n的不确定性影响，当多项式拟合阶数低时，拟合误差大，阶数高时计算效率低且容易溢出。为了弥补这一缺陷，我们将非平稳地震记录时频谱转换到二次时频谱谱中，通过滤波的方式提取子波时频谱[13]。提取非平稳地震记录时频谱|G(τ，f)|每个时间点振幅谱，对其傅里叶变换得到二次时频谱，其结果记为|G(τ，f)|(2)，则式(6)可以表示成为式(8)。

(|R(τ，f)|)(2)

(8)

其中:|W(τ，f)|(2)表示动态子波的二次时频谱；|R(τ，f)|(2)表示反射系数的二次时频谱。

基于Rosa理论，地震子波振幅谱是光滑的，子波和反射系数是可以分开的，在非平稳地震记录二次时频谱中，低频稳定的部分是子波的能量，而高频扰动的部分是由于反射系数引起的。因此可以通过滤波的方式来保留低频的子波能量，再通过相位恢复，反变换得到地震子波时频谱。这里记低通滤波后的地震记录二次时频谱为|G(τ，f)|fit(2)，通过对其傅里叶反变换即可以提取动态子波振幅谱：

|W(τ，f)|fit≈ifft(|G(τ，f)|(2))

(9)

则反射系数时频谱值R(τ，f)可以表示成式(10)。

(10)

为了防止分母出现零值引入参数μ，Amax表示动态子波时频谱|ωa(τ，f)|fit的最大值。再对R(τ，f)进行S反变换就可以获得时间域的反射系数r(t)为式(11)。

(11)

2 模型论证

使用30 Hz到120 Hz线性递增变化的chirp信号与60 Hz到250 Hz二次函数递增的chirp信号组合成一个复合信号如图1(a)所示。利用图1(a)来检测Gabor变换、S变换和改进广义S变换(GST)的时频分析效果。我们可以通过调节GST参数，调节时频分析结果,图1(b)、图1(c)、图1(d)相应的是Gabor变换、S变换(λ=1、r=1)、改进广义S变换(λ=2，r=0.8)对图1(a)的时频分析结果。与图1(b)中的Gabor变换相比，图1c中的S变换分析结果在低频时分辨率高，随着信号频率的升高，时频分辨率降低，图1(d)中的GST方法能精确的检测到两个chirp信号的时间-频率变化趋势，从低频到高频，时频分析结果优于Gabor变换和S变换。

图1 复合信号时频分析对比图Fig.1 The comparison of time-frequency analysis of composite signal(a)Chirp信号;(b)对a短时傅里叶变换结果;(c)对aS变换结果;(d)对a改进广义S变换结果

图2 模拟地震记录二次时频谱反褶积结果图Fig.2 The quadratic time-spectrum deconvolution results of simulation seismic records(a)模拟反射系数；(b)平稳地震记录；(c)非平稳地震记录；(d)反褶积之后的地震记录

为了检验二次时频谱反褶积实用性，建立一维模型对其分析。如图2(a)是一个加微弱噪声的随机反射系数，将其与主频为40 Hz的雷克子波进行卷积得到平稳地震记录(图2(b))；为了满足实际衰减地震模型，引入Q取值50的衰减函数，结合式(3)得到衰减地震道如图2(c)所示。将衰减地震记录进行改进广义S变换分析得到其时频谱如图4(b)，取其中一个时刻的时频谱值进行单独子波提取(图3)，图3(a)中蓝色线表示原始的非平稳地震道振幅谱，其明显受反射系数影响很大，将其傅里叶变换得到二次谱如图3(a)红色线表示，将图3(a)中的二次谱进行低通滤波并反傅里叶变换得到提取后的子波振幅谱(图3(b)中红色线表示)，与原始子波振幅谱(图3(b)绿色线表示)进行对比能很好地得到吻合，这说明二次谱提取子波振幅谱的准确性相比能量有明显的衰减；再用二次谱模拟(式(8)～式(9))依次对每个时间点进行二次时频谱提取，得到衰减子波时频谱如图4(c)所示，再利用(式(10)～式(11))并结合地震道时频谱实现二次时频谱反褶积得到结果如图2(d)所示。可以发现，经过本文的反褶积方法处理后，原始衰减剖面的能量不仅得到了恢复，同时与原始的平稳地震记录相比结果更收敛；对比图2中两个红色框位置可以看出，原始平稳地震记录的两个叠加在一起的地震波被有效的分开,且能和图2(a)相对应，验证了本文的方法能在恢复地震波的能量的基础上，提高对薄层地识别能力。

图3 子波振幅谱提取结果图Fig.3 The result of wavelet amplitude spectrum extraction(a)子波二次谱；(b)提取子波振幅谱对比

图4 不同信号时频谱Fig.4 The time spectrum of different signal(a)平稳信号时频谱；(b)非平稳信号时频谱；(c)对(b)平滑后信号时频谱；(d)反褶积后时频谱

图5 二维地震模型反褶积结果对比图Fig.5 The comparison of two-dimensional seismic model deconvolution results(a)模拟反射系数;(b)平稳地震记录;(c)衰减地震记录;(d)反褶积地震记录

设计一个二维楔形模型验证算法的鲁棒性，如图5(a)设定一个二维反射系数，在200 ms和300 ms的位置分别是是一负一正的两个薄层，400 ms位置建立一个楔形模型；用一个主频为35 Hz雷克子波和图5(a)中的反射系数进行褶积得到地震记录如图5(b)；利用式(3)得到一个衰减二维地震剖面(如图5(c))，其中衰减函数Q值取50；图5(d)是用本文方法对图5(c)中的衰减地震记录进行反褶积得到的结果。对比分析可以得出结论：经过反褶积处理后整个剖面被衰减的能量得到了恢复，子波被有效压缩，提高了地震记录的纵向分辨率。在楔形模型位置，观察到原始地震剖面(图5(b))在第19道之后两个子波才能被分离(红色实线位置)，经过反褶积处理后的剖面(图5(d))在11道(红色实线位置)之后就可以完全的将两个子波识别出来，这验证了本文的方法可以有效地提高地震资料的时间分辨率。综上模型论证说明，基于改进广义S变换二次时频谱反褶积能够恢复地震波能量，同时压缩子波，提高地震记录分辨率。

3 实际资料处理

图6 反褶积前后剖面对比图Fig.6 The comparison of the deconvolution result(a)原始地震剖面；(b)Gabor反褶积结果；(c)二次时频谱反褶积结果

图7 井资料解释对比图Fig.7 The comparison of well data to profile(a)原始地震剖面；(b)Gabor反褶积结果；(c)二次时频谱反褶积结果

图8 反褶积前后剖面振幅谱对比Fig.8 The comparison of the amplitude spectra before and after deconvolution

选取中国南海流花工区数据进行二次时频谱反褶积处理。流花构造位于珠江口盆地白云凹陷东部，北部为白云东凹，处于白云东区深水扇砂体成藏带与东沙隆起复式成藏带交界处，为了达到更好地识别储层的目的，使用本文的方法进行高分辨处理。选取目的层区域如图6(a)所示，对其进行使用Gabor反褶积和二次时频谱反褶积处理得到结果图6(b)、图6(c)，可以发现，二次时频谱反褶积处理之后整个剖面的同相轴被压缩，能量聚焦性更好，原始不明显的地层经过处理后更加地清楚，出现了更多的层位。为了证实处理结果的正确性，将图6中红色长方形标记的过井目的层，结合密度测井曲线进行井震标定，并放大如图7(a)、图7(b)和图7(c)所示。对比原始剖面和Gabor反褶积处理的剖面(图7中红色椭圆位置)，二次时频谱反褶积处理后层位与井曲线对照更加准确且清晰，在2.5 s位置，原始剖面两气层之间的薄盖层也被准确分辨出来。结合测井资料我们对地震剖面进行解释，从二次时频谱反褶积处理后的剖面可以准确的区分出四个气层和两个油层，其中Gas2和Gas4是薄气层，而原始剖面只能区分出两个气层和两个油层，这为后期储层预测提供借鉴。

通过对剖面进行综合振幅谱分析，定量评估了两种反卷积方法的应用效果，对图7的三个剖面求平均振幅谱如图8所示。其中二次时频谱反褶积后的频谱(红色线)比Gabor反褶积数据(绿曲线)的频谱更宽，且都优于原始数据的频谱(蓝色线)。另外，对比三个频谱的主频，二次时频谱反褶积的主频明显高于Gabor反褶积高于原始剖面，说明本文的方法应用效果更优。

4 结论与认识

1)引入改进广义S变换，将时间域反褶积扩展到时频域，更好适应地震波传播时的时变性。

2)基于改进广义S变换的二次时频域反褶积后，能有效地恢复地层对地震波的衰减。同时能有效地压缩子波，提高对薄储层的识别能力。

3)利用非平稳地震记录二次谱提取子波摆脱了反射系数要满足白谱特性的假设，同时没有平滑函数阶数不确定性的问题。

4)由于不同单道记录的时频谱最大值不同，可以求出整个剖面时频谱最大值的平均值，以此来代替Amax进行反褶积效果更好。