李 婧， 王修田,2,3❋❋， 宋 鹏,2,3， 姜秀萍,2,3， 赵 波,2,3

(1.中国海洋大学海洋地球科学学院，山东 青岛 266100；2.青岛海洋科学与技术试点国家实验室 海洋矿产资源评价与探测技术功能实验室，山东 青岛 266237；3.海底科学与探测技术教育部重点实验室，山东 青岛 266100)

反射系数反演方法可以消除干涉效应等影响，拓宽地震记录的频带范围，有效地将反射系数序列恢复出来，从而提高地震地质解释精度。近些年来，为了提高地震反射系数反演的可靠性和准确度，许多学者相继提出和发展了多种时间域和频率域的反射系数反演方法[1-5]。但大多数方法都基于平稳褶积模型[6]，所用的地震子波是时不变子波。在实际地震勘探中，由于受到大地滤波作用等影响，地震子波具有时变和空变特征，因此，地震记录具有非平稳性，采用平稳褶积模型无法准确表达非平稳地震信号。为解决此问题，Clarke[7]提出了非平稳褶积模型，并发展了一种基于最优维纳滤波的时域非平稳反褶积方法。之后，学者们基于非平稳褶积模型研究了不同的反射系数反演方法，这些方法在实际地震资料的反演处理中取得了较好的应用效果[8-12]。Margrave等[13]提出了Gabor域的反褶积方法，其通过在Gabor域估计震源子波和衰减函数进而实现时变反褶积处理，但该方法要求反射系数序列需要满足高斯白噪假设。一般而言，若仅考虑吸收衰减情况，时变子波在某时刻的振幅谱值可视为震源子波的振幅谱值与该时刻衰减函数值之积。据此，在进行非平稳地震资料反射系数反演时，可以由地震记录提取时变子波，进而可实现更高精度的时变反射系数反演。

尽管依据高分辨率测井资料和井旁地震记录可估计出较高精确度的确定性子波，但该子波并不适用于全工区的反演处理。根据地震记录统计特征实现的分段时变子波提取，需要假设各时窗段内的地震记录近似为平稳信号，在这种情况下，地震子波在各时窗段内被认为是时不变的。虽然分段提取地震子波最终可以达到时变子波估计的目的[14-15]，但受时窗段的划分方式(例如时窗长度)与子波平均效应的影响，分段提取的时变子波无法精细地反映出相邻地层反射波的子波变化情况。鉴于时频分析方法可以较为准确地描述地震数据时间与频率之间的关系，以此为基础，学者们提出了适应地震资料时变特征的逐点时变子波提取方法。该类方法首先对地震记录进行时频变换，再在时频域内利用平滑方法实现逐点时变子波的提取。邬世英等[16]在Gabor域基于动态褶积模型研究了反射系数(满足白噪和非白噪情况)、地震记录和子波之间的关系，并指出地震记录对数时频谱为子波和反射系数二者对数时频谱之和，因此可利用平滑地震记录的对数时频谱求取子波谱。此外，由Rosa和Ulrych[17]提出的谱模拟法因其无需预先假设反射系数振幅谱的形态，对于白噪和非白噪反射系数序列也都有很好的适用性，因而可将其拓展到基于时频变换的时变子波提取方法[18-20]中。在标准S变换基础上，李振春等[21]采用谱模拟法提取时变子波，实现了S变换域的时变反褶积。由于S变换其采用的高斯窗函数的尺度因子仅与频率有关，变换趋势单一，王元君等[22]利用广义S变换具有更高灵活性和时频分辨率的优势，提出了一种基于广义S变换的动态反褶积方法，其无需直接求取品质因子(Q值)，且适用于变Q值情况，但广义S变换经常会出现信号高频和低频段精度变差的情况，会影响反演结果的准确性。

本文首先根据地震反射系数的稀疏性构建L1范数约束的反射系数反演模型，利用稀疏约束正则策略提高反演精度；然后，以影响反演精度和稳定性的地震子波为主要考虑因素，研究发展了基于改进广义S变换时变子波提取与时变反射系数反演方法。该方法通过提取适应地震数据时变特征的地震子波，重建地震子波褶积矩阵，并求解相应的稀疏约束正则化目标函数，从而可获取与地下真实情况更为匹配的反射系数。

1 L1范数约束下的统一子波反射系数反演

地震波反射发生在地下介质声波阻抗具有明显变化的位置，因此，野外采集的地震数据中蕴含丰富的地层特征信息。Robinson提出的地震记录褶积模型，阐明了自激自收的地震波场近似于地震子波同地层反射系数序列的褶积[6]。在考虑噪声的情况下，平稳地震记录褶积模型的矩阵形式可以写为：

s=Wr+n。

(1)

式中：s表示单道地震记录观测值s(t)组成的列向量；W表示地震子波w(t)组成的褶积矩阵；r表示反射系数序列r(t)构成的列向量；n表示该道噪声值n(t)组成的列向量。用符号i表示时间采样点(i=1,…,N，N为总时间采样点数)，则式(1)按时间展开的具体形式可表示为：

(2)

根据地层的构造特点，大反射系数通常被认为对应的是主要岩性界面或不整合面，其个数远小于整个地震记录采样点的数目，即反射系数序列具有稀疏性，因而可以采用稀疏约束策略解决反演问题的不确定性和多解性。为了获得稀疏的反射系数解，将L1范数用于正则化目标函数，其数学公式如下：

(3)

2 基于时频变换的时变反射系数反演

2.1 基于改进广义S变换的时变子波提取

Rosa和Ulrych认为地震子波的振幅谱是接近平滑的，而反射系数的振幅谱是相对震荡的，据此提出了在最小平方意义下的谱模拟法[17]，其主要思想是通过多项式拟合的平滑方式由地震记录的振幅谱获得子波振幅谱。

现假设在时刻τ的地震子波振幅谱|W(τ,f)|是单峰光滑曲线，将谱模拟法扩展到基于时频变换的时变子波提取，则其在频率域的数学模型可构建为：

(4)

其中：f表示频率；k表示常数；n表示阶数；在一般情况下，0≤k≤3，4≤N≤7[21]。固定N和k的值，经过多项式拟合求出an(n=0,1,…,N)后，即可根据式(4)得到一条光滑的拟合曲线|W(τ,f)|，改变时刻τ的值，即可按时间顺序估算地震子波，这个过程就是时变子波的提取过程。

((a) 反射系数序列；(b) 非平稳地震记录。(a) Reflection coefficient sequence；(b) Non-stationary seismogram.)

((a) 主频为50 Hz；(b) 主频为30 Hz。红色线为反射系数反演结果，黑色线为原设计的反射系数。(a) 50 Hz dominant frequency；(b) 30 Hz dominant frequency. The red lines show the inversion results and the black lines represent the designed reflection coefficients.)

图2 采用固定主频子波的常规反射系数反演效果对比图

Fig.2 Comparisons of traditional reflectivity inversion results using wavelets with fixed dominant frequencies

Stockwell结合短时傅里叶变换的相位特性与小波变换的多分辨率特性，提出了无损可逆的S变换(S Transform)[24]。对于信号s(t)，其时间域S变换数学表达式为：

(5)

MGST(S(τ,f))=

(6)

在应用(6)式时可保持参数q=1不变，让p随频率线性减小(可根据不同类型信号和用途通过a和b调节其变化范围)，这样可使得改进广义S变换能够对所有频率成分的时频谱保持较好的聚焦性，由此可提取不同时刻的时变子波，其主要提取过程如下：

①对时间域的非平稳地震记录s(t)沿时间t方向进行改进广义S变换，得到相应的时频谱|MGST(S(τ,f))|；

②固定时间T，对该时刻的地震记录振幅谱|MGST(S(T,f))|进行多项式拟合，通过合理设置参数k和N，得到对应时刻的子波振幅谱|MGST(W(T,f))|；

③改变τ的取值，逐点计算子波振幅谱|MGST(W(τ,f))|，实现“时变子波”的时频谱提取；

④拟合出“时变子波”时频谱后，若子波相位恒为零，则各时刻的子波振幅谱经过傅里叶反变换，即是所求零相位时变子波；当子波为混合相位时，可以利用希尔伯特变换(Hilbert Transform)获得子波最小相位谱，通过Z变换进一步求得相同振幅谱的全部相位谱系列[27]，然后采用子波相位谱扫描技术进行地震子波的混合相位提取[28]。在以下的实验中假定地震数据已通过零相位化处理，因此仅考虑零相位时变地震子波的提取即可。

2.2 时变反射系数反演

式(2)中所示的子波褶积矩阵是由非时变子波构成的，现改用时变子波重新构建地震子波褶积矩阵，其数学表达式可写为：

(7)

在式(7)表示的矩阵里，每一个列向量wi均可不同，即每一个时刻的地震子波均可不同。由此，可以将标准褶积模型改写为时变褶积模型：

(8)

同样，相应时变反射系数反演的目标函数可表示为：

(9)

通过矩阵向量化，时变反射系数反演方法可以由单道地震记录拓展到多道情况，多道地震数据的褶积模型可表示成如下矩阵-向量乘积形式：

(10)

(11)

3 数值模型实验

对图1(b)所示的非平稳地震记录(为方便对比分析，将其以时间轴为纵向坐标重绘于图3(a))。首先分别采用标准S变换(ST)和改进广义S变换(MGST)求取地震记录的时频谱图(分别见图3(b)和(d))，再由其求取相应的不同“时变子波”的时频谱图(分别见图3(c)和(e))。通过对比四幅时频谱图可以看出，基于改进广义S变换计算的地震记录时频谱及提取的“时变子波”时频谱的聚焦效果均好于基于标准S变换的计算结果，这将有利于时变地震子波提取的精确性。

图4(a)和(b)分别给出了在0.10和0.91 s两个不同时刻所提取的时变子波示例。由其可见，基于改进广义S变换提取的时变子波与原始子波吻合度更好。图5给出了应用不同时变子波进行时变反射系数反演的结果，比对分析表明基于改进广义S变换的时变反射系数反演解的精度较高，其与原设计的反射系数吻合程度明显优于基于标准S变换的计算结果。

为进一步检验基于改进广义S变换的时变反射系数反演方法的抗噪能力，对图3(a)所示的地震记录加入不同程度的随机噪声(图6(a)和(d)分别为含噪10%和含噪20%的地震记录)，再进行基于改进广义S变换的时变子波提取与时变反射系数反演。图6(b)和(e)分别为对应含噪10%和含噪20%的地震记录的时频谱图，图6(c)和(f)分别为对应(b)和(e)的“时变子波”时频谱图。可见，即使在含噪10%和20%的情况下，基于改进广义S变换计算的时频谱及“时变子波”的时频谱仍然具有相对较好的聚焦性。

图7分别给出了在0.10和0.91 s两个不同时刻估算的地震子波示例，将其与原设计子波(见图7中的黑色虚线)比较可知，所提取的时变子波并没有因为噪声的加强而剧烈变化。由图8所示的不同噪声水平下的反演结果可知，除了噪声水平较大时产生的一定局部误差外，时变反射系数反演的结果与原设计的反射系数模型仍具有较好的一致性，由此表明基于改进广义S变换的时变反射系数反演方法具有较好的容噪能力和稳定性。

((a) 时变子波合成地震记录；(b) 标准S变换时频谱；(c) 由标准S变换时频谱提取的“时变子波”时频谱；(d) 改进广义S变换时频谱；(e) 由改进广义S变换时频谱提取的“时变子波”时频谱。(a) Synthetic seismogram by the time-varying wavelets；(b) Time-frequency spectrum of ST；(c) Time-frequency spectrum of the ‘time-varying wavelets’ based on ST；(d) Time-frequency spectrum of MGST；(e) Time-frequency spectrum of the ‘time-varying wavelets’ based on MGST.)

图3 时变子波合成地震记录及相应时频谱图

Fig.3 Synthetic seismogram by the time-varying wavelets and the corresponding time-frequency spectrums

((a) 0.10 s；(b) 0.91 s。红色线为基于改进广义S变换提取的时变子波；蓝色线为基于标准S变换提取的时变子波；黑色线为原设计子波。(a) 0.10 s；(b) 0.91 s. Red curves represent the wavelets extracted based on MGST；Blue curves show the wavelets extracted based on ST；Black curves are the designed wavelets.)

图4 不同时刻时变子波与原设计子波对比图

Fig.4 Comparisons of time-varying wavelets at different time points with the designed wavelets

((a) 基于标准S变换；(b) 基于改进广义S变换。蓝色线为基于标准S变换的时变反射系数反演结果；红色线为基于改进广义S变换的时变反射系数反演结果；黑色线为原设计的反射系数。(a) ST-based；(b) MGST-based. Blue curve shows the ST-based time-varying inversion result；Red curve represents the MGST-based time-varying inversion result；Black curves are the designed reflection coefficients.)

图5 基于不同时频变换的时变反射系数反演效果对比图

Fig.5 Comparisons of time-varying reflectivity inversion results based on different time-frequency transforms

((a) 含噪10%的合成地震记录；(b) 含噪10%时改进广义S变换时频谱图；(c) 含噪10%的“时变子波”时频谱图；(d) 含噪20%的合成地震记录；(e) 含噪20%时改进广义S变换时频谱图；(f) 含噪20%的“时变子波”时频谱图。(a) Synthetic seismogram with 10% noise added；(b) MGST-based time-frequency spectrum at 10% noise level；(c) Time-frequency spectrum of the ‘time-varying wavelets’ based on MGST at 10% noise level；(d) Synthetic seismogram with 20% noise added；(e) MGST-based time-frequency spectrum at 20% noise level；(f) Time-frequency spectrum of the ‘time-varying wavelets’ based on MGST at 20% noise level.)

图6 不同噪声水平下的合成地震记录及相应时频谱图

Fig.6 Synthetic seismograms at different noise levels and the corresponding time-frequency spectrums

((a) 0.10 s；(b) 0.91 s。红色线为含噪10%情况下提取的时变子波；蓝色线为含噪20%情况下提取的时变子波；黑色线为原设计子波。(a) 0.10 s；(b) 0.91 s. Red curves represent the wavelets extracted at 10% noise level；Blue curves show the wavelets extracted at 20% noise level；Black curves are the designed wavelets.)

图7 不同含噪水平下基于改进广义S变换时变子波提取效果对比图

Fig.7 Comparisons of time-varying wavelets extracted based on MGST with the designed wavelets at different noise levels

((a) 含噪10%；(b) 含噪20%。红色线为含噪10%时基于改进广义S变换时变反射系数反演结果；蓝色线为含噪20%时基于改进广义S变换时变反射系数反演结果；黑色线为原设计的反射系数序列。(a) At 10% noise level；(b) At 20% noise level. Red curve represents the inversion result based on MGST at 10% noise level；Blue curve shows the inversion result based on MGST at 20% noise level；Black curves are the designed reflection coefficients.)

图8 不同含噪水平的基于改进广义S变换时变反射系数反演效果对比图

Fig.8 Comparisons of time-varying reflectivity inversion results based on MGST with the designed reflection coefficients at different noise levels

4 实际资料处理

图9所示为一实际地震叠前时间偏移剖面，其CDP道号范围为500～1 600，截取的时窗范围为2.40～4.40 s，图中红线标注位置为第800道地震记录。

基于改进广义S变换由地震剖面逐道提取时变地震子波，图10给出了第800道地震记录及基于改进广义S变换的时频谱图，图11中的红色、绿色与蓝色线分别表示在2.80、3.65和4.00 s三个时刻提取的时变子波。可见，即使在同一位置，随着时间的增加，地震记录的高频成分不断损失，地震子波波形逐渐变宽，这符合其在地下传播时的频率衰减特性，也反映了地震记录的时变特征；与传统的谱模拟法提取的时不变子波(由图11中的黑色虚线表示)比较可知，不同时刻的地震子波波形均有明显差异。

图9 地震叠前时间偏移剖面

((a)第800道地震记录；(b)改进广义S变换时频谱；(c)“时变子波”时频谱。(a) Seismic record at CDP 800；(b) Time-frequency spectrum of seismic record based on MGST；(c) Time-frequency spectrum of ‘time-varying wavelets’.)

图10 第800道地震记录及时频谱图

Fig.10 Seismic record at CDP 800 and the corresponding time-frequency spectrums

(红色线为2.80 s时刻提取的时变子波；绿色线为3.65 s时刻提取的时变子波；蓝色线为4.00 s时刻提取的时变子波；黑色线为常规提取的时不变子波。Red curve represents time-varying wavelet extracted at 2.80 s；Green curve shows time-varying wavelet extracted at 3.65 s；Blue curve shows time-varying wavelet extracted at 4.00 s；Black curve is the conventional time-invariant wavelet.)

图11 由第800道地震记录提取的时变子波及常规提取的时不变子波示例

Fig.11 Time-varying wavelets of CDP 800 at different time points and the conventional time-invariant wavelet

考虑到地震子波的空间连续性，以提取的时变子波为基础，对所估计的地震子波采取多道加权空变处理[19]，再利用本文方法反演求取反射系数剖面如图12(a)所示；同时采用时不变子波进行常规反射系数反演，所得反射系数剖面如图12(b)所示。

为比较本文时变子波反演方法与常规时不变子波反演方法的精确度，将以上两种方法得到的反射系数分别与对应的时变/时不变地震子波褶积合成地震记录，以合成的地震记录与叠前时间偏移剖面中的原始地震记录之间的平均绝对误差(见图13)来评价两种方法反演结果的精确度。由图13可见，基于改进广义S变换时变反射系数反演结果合成的地震记录与原始地震记录的平均绝对误差，明显小于基于常规时不变子波反演结果的计算误差。由此表明，由于利用了地震子波的时/空变特性，基于改进广义S变换的整体反演精度得以提升。

5 结论

将基于改进广义S变换提取的时变地震子波用于重建褶积模型的核矩阵，可通过求解相应的L1范数稀疏约束问题实现时变反射系数反演。数值模拟和实际资料处理结果均表明：

((a) 基于改进广义S变换的时变反射系数反演结果； (b) 常规反射系数反演结果。(a) Time-varying reflectivity inversion result based on MGST；(b) Traditional reflectivity inversion result.)

图12 不同反射系数反演方法结果对比图

Fig.12 Comparison of different reflectivity inversion results

(红色线为基于改进广义S变换的时变反射系数反演结果合成地震记录与原始地震剖面的平均绝对误差；黑色线为传统反射系数反演结果合成地震记录与原始地震剖面的平均绝对误差。Red curve represents the mean absolute error between seismogram synthesized by time-varying reflectivity inversion result based on MGST and the seismic section；Black curve shows the mean absolute error between seismogram synthesized by traditional reflectivity inversion result and the seismic section.)

图13 反射系数反演解的误差评价

Fig.13 Error evaluation of inverted reflection coefficients

(1) 基于改进广义S变换相应时频谱提取的时变子波与原始子波吻合度较高，符合地震数据的时变/空变特征，且具有较好的容噪能力。

(2) 相比采用统一子波的常规反演方法，基于改进广义S变换的时变反射系数反演方法利用了地震子波的时/空变特性，可获得更高精度的反射系数剖面。