基于时变子波的基追踪地震反演

周玉毅, 周怀来,2, 廖璐瑶, 王元君,2, 张 舜

(1.成都理工大学 地球物理学院，成都 610059；2.“油气藏地质及开发工程”国家重点实验室，成都 610059)

0 引言

地震子波提取的研究由来已久，许多研究者相继提出和发展了多种提取地震子波方法[1-3]。地震子波估计准确程度直接影响着反演的精度与可靠性，传统的地震子波估计方法应用前提是,地震波在传播的过程中其波形以及能量保持不变，然而由于地下介质对地震波能量的衰减吸收，导致地震子波在地层传播过程中将产生高频成分缺失、相位畸变以及频散现象，即实际地震记录中的地震子波具有时变性，因此基于Robinson平稳褶积模型估计子波的理论是不完善的。

针对地震子波的非平稳特性，Clarke[4]首次提出非平稳褶积模型;Margrave等[5]在非平稳褶积模型基础上研究了地震子波具有时变性的原因;彭才等[6]建立符合地震子波衰减特性的动态褶积模型，在随机白噪反射系数和地震子波最小相位的假设前提下，通过对非平稳地震道进行加高斯窗傅里叶变换后，可以获得地震道振幅谱，再经方形窗平滑处理得到动态子波；Margrave等[7]在Gabor变换的基础上，模拟各个时窗内地震记录振幅谱提取动态子波。上述方法应用的前提是反射系数序列为白噪声的假设，后续研究者采用谱模拟技术对其进行规避。Rosa等[8]结合Ricker的研究结果，认为子波在频率域内的曲线趋势是光滑连续的，与之相对应的是反射系数在频率域内是振荡的，地震记录的振幅谱主要受地震子波影响，而反射系数只是对细节部分有所影响，因此可以通过模拟地震记录振幅谱来提取子波。由于谱模拟方法对于反射系数序列有很好的适应性，因此，诸多学者将该方法拓展到高精度时频分析方法中来提高时频聚焦性。李振春等[9]利用S变换将地震记录转换到频率域，在频率域内用数学方法的多项式拟合出时变子波的振幅谱，在地震子波相位为零的假设下，得到频率域中的时变子波;张漫漫等[10]利用广义S变换结合时频滤波求取时频谱，对时频谱进行多项式谱模拟来提取时变子波;Zhou等[11]对谱模拟反褶积方法进行推广，提出基于广义S变换的动态反褶积算法，估算出的子波准确性有所提高。在考虑实际地层吸收衰减基础上，时变子波谱可近似表示为震源子波谱与衰减函数的乘积形式，因此可以从实际地震记录中提取时变子波进一步提高反射系数反演精度。姚振岸等[12]对非平稳地震道进行基追踪谱分解获得各时刻的局部功率谱，采用广义地震子波构建方法提取时变子波，最终实现非平稳地震道的基追踪反演;李婧等[13]基于改进广义S变换提取符合地震数据非平稳特性的时变子波，建立时变子波矩阵，通过求解对应的稀疏约束正则化目标函数获得与实际接近的反射系数。

笔者以前人的研究为基础，利用广义S变换谱模拟方法从地震资料中提取时变子波，用估算的时变子波逐点替换时不变子波构建时变子波库，对地震记录采用时变子波进行基追踪反演。理论模型和实际工区地震资料测试结果表明，基于本文方法提取的时变子波基追踪地震反演的结果,与传统时不变子波基追踪反演的结果相比，井震标定吻合度更高、横向连续性更好以及细节展示更为丰富，说明对非平稳地震信号进行基追踪反演时采用时变子波库的有效性。

1 方法原理

1.1 基于广义S变换的时变子波提取

Rosa结合Ricker的研究结果分析认为：子波振幅谱的曲线趋势是光滑连续的，与之相对应，反射系数振幅谱的曲线趋势是波动的。因此，只要设置适当的拟合参数，就能采用数学手段对地震记录的振幅谱利用多项式拟合进行处理，最后得到子波振幅谱[8]。现在假设一个地震子波的振幅谱W(t,f)曲线趋势与雷克子波单峰光滑连续的振幅谱相似，则可以将其扩展为基于时频变换的时变子波，从而在频率范围内构成一种可用的子波数学分析模式：

(1)

其中:f为频率；k、n为常数；an为关于f的多项式的系数。通常情况下，0≤k≤3，4≤n≤7[9]。

描述信号时间和频率联合分布的函数有很多种，如STFT、CWT、Hilbert-Huang Transform等，为了实现局部时频分析，Stockwell等[14]继承和扩展了STFT和CWT的优势方面，提出了一种与其他时频分析方法相同的无损可逆S变换，该方法对信号进行处理后，可以观察到信号在频域和时域的能量分布。设x(t)为平方可积的函数，则x(t)的S变换可以定义为：

(2)

式中：τ表示时间;f表示频率;ST(x(τ,f))表示x(t)的S变换。

其中S变换的基本小波函数定义为：

gf(t)exp(-i2πft)

(3)

式中gf(t)为高斯窗函数，可定义为式(4)。

(4)

由于基本小波窗函数形态固定，S变换不能够依据实际情况来调整窗口的大小，在高频部分分辨率低，导致在实际应用中不能得到较好的效果。针对S变换的问题，陈学华等[15]在Stockwell的S变换基础上提出了广义S变换(Generalized S Transform)，通过加入两个参数rgs和ρ来调节高斯窗函数的大小，其数学表达式为式(5)。

exp(-i2πft)dt

(5)

通过调整控制rgs和ρ的大小，该方法能够对所有的时频谱具有较高的时频聚焦性，可以用来提取不同时间点的时变子波，利用广义S变换时频谱提取子波过程(图1)如下：

图1 时变子波提取流程图

1)对沿t方向的非平稳地震记录进行广义S变换，得到时频谱|GST(x(τ,f))|。

2)固定某一个时间T，通过设置合理的参数k和n，采用数学手段对该时间点地震记录的振幅谱|GST(x(T,f))|利用多项式拟合进行处理，拟合得到相应的子波振幅谱|GST(W(T,f))|。

3)改变τ的取值，逐点计算子波振幅谱，从而实现“时变子波”的时频谱提取。

4)在子波相位为零的假设前提下[9]，对各时刻的子波振幅谱进行傅里叶逆变换后，可求得各时刻零相位时变子波时域表达式。

1.2 时变子波基追踪反演

地震记录通常是由一个地震子波和地层内部反射界面的反射系数进行褶积而得到的，用公式表示如下：

s(t)=w(t)*r(t)

(6)

其中：w(t)为一个地震子波；r(t)为一个反射系数序列；“*”为褶积运算的符号。

将式(6)写为矩阵形式:

s=Wr

(7)

其中：假设地震数据的采样点为M；s是由地震数据所构成的M×1矩阵;W是由地震子波w(t)经过逐点移动构成的M×M对角矩阵；r是由地下反射系数构成的M×1矩阵。虽然式(6)和式(7)一个是褶积运算，一个是矩阵运算，但是它们的运算结果是完全相等的。

基追踪是一种基于范数约束的最优化准则，基追踪反演结果与实际地层相接近，分辨率较高，近年来逐渐被应用于地震反演，其中目标函数被定义为式(8)[16]。

J=argmin||s-Wr||2+λ||r||1

(8)

其中：J为目标函数；下标1和2为目标向量，分别用于表示L1和L2范数；λ是影响调整值的因子。满足J→min即为一个解函数向量r的最终函数值。

为了表征地下介质的稀疏性，而不是反射系数的稀疏，根据Helholm定理，将原始反射系数对分解成如图2所示的偶脉冲对与奇脉冲对的线性组合，则任意反射系数序列可表示为式(9)。

bm,nro(t,m,n,Δt)}

(9)

其中：ro(t,m,n, Δt)为奇楔反射系数对分量；re(t,m,n, Δt)为偶楔反射系数对分量；am,n和bm,n为对应奇、偶反射系数对对应的系数；Δt为采样间隔；m、n分别为地下介质层顶、底界面所对应的采样点位置；由于不同反射对出现的位置不固定，顶底反射系数以平移量mΔt沿着图2中t的方向进行逐点移动，其中m的取值范围为“0”到M，M、N由实际地震资料给定。

根据褶积定理，将反射系数对与地震子波进行褶积可得到地震信号:

re(t,m,n,Δt)}+bm,n[w(t,m,n)*

ro(t,m,n,Δt)]}

(10)

用矩阵的形式表示为式(11)。

s=Dwr

(11)

也就是说，地震响应是通过奇、偶楔形反射系数模型和地震子波褶积获得的奇、偶地震响应的线性组合。式(11)中，Dw是奇、偶楔形模型地震响应。传统的基追踪反演是在褶积模型假设的前提下，用统计方法或井震标定方法从实际地震数据提取出地震子波w(t)。由于随机噪声干扰、地下介质吸收衰减作用的影响，从实际地层采集的地震数据呈现非平稳特性，运用本文方法逐道提取各时刻的子波，逐列替代时不变子波构建时变子波核矩阵，可得：

re(t,m,n,Δt)]+bm,n[w(t,m,n)*

ro(t,m,n,Δt)]}

(12)

反演过程中，目标函数变为：

(13)

采用基追踪最优化准则从式(12)中求得参数am,n和bm,n，因为式(9)和式(12)采用相同的系数，因此将am,n、bm,n代入式(9)通过求和重构可得到有效的反射系数。在准确提取地震时变子波的前提下，通过基追踪最优化准则可求得最优的反射系数值。

2 模型测试

2.1 不含噪声模型

为了检验时变子波库的建立对于基追踪反演的重要性，设计如图3(a)所示以1 ms为采样间隔且时间长度为1 000 ms的反射系数序列，其中各反射系数脉冲的值分别设定为(120 ms，-0.15)、(165 ms，0.09)、(380 ms，-0.06)、(540 ms，-0.06)、(760 ms，0.13)、(900 ms，0.18)，非平稳地震记录由主频从60 Hz单调下降到40 Hz的Ricker子波与设计反射系数褶积合成，如图3(b)所示。分别采用60 Hz子波字典、50 Hz子波字典、40 Hz子波字典以及由本文方法提取的随时间变化子波构建的子波字典，对合成的非平稳地震记录进行基追踪反演，其反演结果如图4所示。由图4(a)～图4(c)可知，反演所用子波的主频越接近合成地震记录所用的子波主频时，反演出的反射系数才更接近实际，否则反演出的反射系数的值与实际相差较远，甚至出现虚假的反射系数。由此可知，针对非平稳地震记录而言，若采用固定主导频率的时不变子波字典用于基追踪反演，则难以确保反演结果的准确度和可靠性，而采用时变子波字典进行基追踪反演时能获得准确的反射系数系列(图4(d))，说明时变子波的估计准确程度对于基追踪反演精度有重要影响。

图3 反射系数序列及合成地震记录

图4 基追踪反演结果

2.2 含噪声模型

为检验本文方法的抗噪性，对非平稳地震道分别添加SNR=15 dB和SNR=25 dB的随机噪声后,得到如图5(a)与图5(d)所示的含噪非平稳地震记录，对加噪后的非平稳地震道采用上述方法处理,可得到含噪非平稳地震道的时频谱和“时变子波”时频谱。对图5(b)和图5(e)进行分析可知，广义S变换在含噪的情况下依旧具有较好的时频聚焦性。利用本文方法分别对含噪地震数据提取时变子波，进而基追踪反演，反演结果如图6所示，将反演出的反射系数与原始反射系数作对比可知，时变子波基追踪反演结果与原始反射系数结果近似，由此可以得出本文子波提取及反演方法具有较好的抗噪性。

图5 含噪地震记录及时频谱

图6 含噪地震记录基追踪反演结果

3 实际应用

这里所用数据为南海某实际工区的叠后过井地震数据剖面如图7所示，该地震剖面共有201道，每道采样点数为256，截取的时窗范围为2 401 ms～2 912 ms，其中测井曲线为密度曲线。

图7 实际地震数据剖面

选择图7中紫色虚线标注的第86道地震记录(图8(a))进行广义S变换处理后，可以获得如图8(b)所示的时频谱图，对时频谱进行谱模拟后得到如图8(c)所示的 “时变子波”时频谱图。对比图8(b)和图8(c)可知，时变子波时频谱消除了地层反射系数的影响，与地震波形吻合度较高。

图8 第86道地震记录及时频谱图

利用广义S变换谱模拟方法对该地震剖面逐道提取时变子波，其不仅沿时间轴方向变化，而且还沿着联络测线的方向变化。在整个地震剖面中，对每道地震记录的所有时间点提取随时间变化的地震子波，而对于地震剖面而言，提取的地震子波随时间变化的同时，还随空间变化。考虑到地震子波在不同道之间的差异，以提取的时变子波为基础，对所估计的地震子波采用多地震道加权空间变化处理[17]。为了研究时不变子波和时变子波对反演结果精度的影响和差异，在此基础上采用本文方法提取的时变子波和时不变子波,用于实际工区地震资料的基追踪反演，反演结果如图9所示。对比图9(a)与图9(b)可以看出，与时变子波基追踪结果图中的绿色方框部分相比，时变子波基追踪地震反演结果的分辨率要高于时不变子波基追踪反演结果的分辨率；时不变子波基追踪反演出的反射系数出现虚假反射，时变子波基追踪反演出的反射系数与原始剖面更为接近；同时，与时不变子波基追踪结果图中的绿色方框部分相比，时变子波基追踪反演出的反射系数横向连续性更好，究其原因，是因为本文方法在求取时变子波时还考虑了地震数据纵横向的变化特性。

图9 不同子波字典基追踪反演剖面对比图

4 结论

1)本文方法提取的时变子波符合地震子波在地下传播时的频率衰减特性，利用时变子波建立的子波库具有较高的精度，为后续的基追踪反演提供了保证。

2)笔者提取时变子波过程中，没有考虑子波的相位影响，后续可以结合子波相位的估计来提取时变子波。

3)实际地震剖面的反演结果表明，基于本文方法提取的时变子波基追踪地震反演的结果，与传统时不变子波基追踪反演得到的反演结果相比，能进一步提高反演精度，可用于后期储层预测与精细刻画。