0 引言

辽西M工区目的层位于潜山上方，有效地层圈闭发育在不整合面附近。潜山不整合面整体上地震反射较强、连续性较好，下伏地层削截明显，同向轴能量减弱，上覆地层同相轴具有上超的特征，地层圈闭识别困难。如何有效提高该地区的地震剖面分辨率、准确识别地层尖灭点，是描述和刻画该区地层圈闭的关键。为此，本文通过时频(S)域相位校正反褶积处理提高工区地震资料分辨率，以达到准确识别尖灭点的目的。

由于介质的吸收效应，随着传播距离增加，地震记录纵向分辨率降低。Wang^[1-3]认为非弹性介质的衰减可以通过反Q滤波补偿，因此很多学者对反Q模型进行了研究，在叠加偏移之前或者在偏移的过程中对地震数据进行反Q滤波，以恢复地震波能量和校正地震记录相位。该处理方法是通过补偿信号高频段的能量，提高原始信号频带范围，从而提高地震资料的高分辨率。反Q滤波最大的困难在于如何准确估计地层品质因子Q，以及保证反Q过程中算法的稳定性。反Q滤波的不稳定性主要由振幅补偿时指数增益导致，但其相位校正过程是稳定的。Margrave等^[4]基于Gabor变换提出了非平稳反褶积，通过在每个时间点的反褶积处理，避免了子波衰减问题。金明霞等^[5]对Gabor反褶积进行了研究，并应用于深水数据处理。在Gabor反褶积算法中，由于Gabor变换时窗固定，影响了时频域反褶积的结果。Zhou等^[6-7]、郭廷超等^[8]在S域实现反褶积，避免了Gabor变换时窗固定的缺陷，但是没有解决非平稳地震记录的相位畸变导致信号从时频域转换到时间域中产生的误差。

M工区地震数据整体频带较宽，高频信号能量弱。常规S变换基函数固定，高频信号识别能力低，达不到精细时频分析的要求。为此，对S变换窗函数进行扩展，提出并推导了改进的广义S变换。调节改进广义S变换参数，可获取非平稳地震记录的精细时频谱，有效地提高时频聚焦度，达到对信号的高频成分进行精确分析的目的。

本文首先利用相位校正将地震记录零相位化；再使用改进的广义S变换对没有相位畸变的地震记录进行精细时频谱分析，从而在S域中精确提取子波，实现S域反褶积。与传统时间域反褶积相比，相位校正S域反褶积能得到更高分辨率的地震资料，可有效识别出M工区的地层尖灭点位置和薄互层，为地层圈闭刻画和储层综合研究提供了优质的资料基础。

1 方法原理 1.1 改进广义S变换

Stockwell在短时傅里叶变换和小波变换的基础上提出了S变换，基于高斯窗能量归一以及其在时域和频域形态相同的优点^[9]，实现了无损的S正、反变换，提高了信号时频分析精度。信号x(t)的S变换为

$ \begin{array}{*{20}{l}} {S(\tau ,f) = \int_{ - \infty }^\infty x (t)\left\{ {\frac{{|f|}}{{\sqrt {2\pi } }}\exp \left[ {\frac{{ - {f^2}{{(\tau - t)}^2}}}{2}} \right] \times } \right.}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left. {{\rm{exp}}( - {\rm{i}}2\pi ft)} \right\}{\rm{d}}t} \end{array} $

(1)

式中τ为变换域时间。S变换的基本小波函数定义为

$ \begin{array}{*{20}{l}} {h(t) = \frac{{|f|}}{{\sqrt {2\pi } }}{\rm{exp}}\left( {\frac{{ - {t^2}{f^2}}}{2} - {\rm{i}}2\pi ft} \right)}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = g(t){\rm{exp}}( - {\rm{i}}2\pi ft)} \end{array} $

(2)

其中g(t)为高斯窗函数

$ g(t) = \frac{{|f|}}{{\sqrt {2\pi } }}{\rm{exp}}\left( { - \frac{{{t^2}{f^2}}}{2}} \right) $

(3)

由上式可知，S变换将窗函数与信号的频率建立起直接联系，因此可根据信号频率变化调节窗函数的宽窄，提高时频分辨率。但S变换窗函数对信号频率变化不够敏感，在实际处理的过程中缺乏灵活性。为了适应海上宽频地震数据，对S变换的窗函数进行改进，以提高高斯窗函数对频率的敏感程度。

引入参数λ、r调节高斯函数的形状，即

$ g(t) = \frac{{|f{|^r}}}{{\lambda \sqrt {2\pi } }}{\rm{exp}}\left( { - \frac{{{t^2}|f{|^{2r}}}}{{2{\lambda ^2}}}} \right) $

(4)

则信号x(t)的改进广义S变换为

$ \begin{array}{*{20}{l}} {G(\tau ,f) = \int_{ - \infty }^{ + \infty } x (t)\frac{{|f{|^r}}}{{\lambda \sqrt {2\pi } }}{\rm{exp}}\left[ {\frac{{ - {{(\tau - t)}^2}|f{|^{2r}}}}{{2{\lambda ^2}}}} \right] \times }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{exp}}( - {\rm{i}}2\pi ft){\rm{d}}t} \end{array} $

(5)

1.2 相位校正S域反褶积

Futterman^[10]基于品质因子Q与频率无关的假设和一维波动方程，提出了描述地震波振幅衰减和速度频散的方程。基于Futterman模型的地震波吸收衰减方程为

$ S(\tau + \Delta t,f) = S(\tau ,f){\rm{exp}}\left[ {{\rm{i}}\varphi (f) + \frac{{\varphi (f)}}{{2Q}}} \right] $

(6)

式中：S(τ，f)为S域τ时刻的地震波波场；φ(f)为频率域的衰减函数，即

$ \varphi (f) = - {\left( {\frac{f}{{{f_{\rm{h}}}}}} \right)^{ - \gamma }}f\Delta t $

(7)

其中$\gamma=\frac{1}{\pi Q}$；f_h为有效频带的最高频率，Δt是时窗长度。

由式(6)可以看出，Futterman的衰减模型表现为，随着传播时间的增加, 信号的振幅衰减和相位产生畸变。为了消除这种影响，一般对信号进行反Q滤波处理^[11]实现相位校正和振幅补偿。针对反Q滤波特点，本文将稳定的相位校正算法与时频域反褶积算法进行结合。时频域反褶积的精度主要是取决于精细的时频谱分析以及时频域子波谱的准确提取。S域子波谱一般利用多项式平滑的方式进行提取^[6-7]，即

$ |W(t,f)| = |f{|^k}{\rm{exp}}\left[ {\sum\limits_{n = 0}^N {{a_n}} (t){f^n}} \right] $

(8)

式中k、N为正整数，取值范围为1 < k < 3、2 < N < 7。利用式(7)提取子波谱均基于子波是零相位的假设，当子波不满足该假设时，时频域反褶积结果从S域转换到时间域会产生误差。

根据式(5)，信号相位校正可表示为

$ S(\tau + \Delta t,f) = S(\tau ,f){\rm{exp}}[ - {\rm{i}}\varphi (f)] $

(9)

对S(τ+Δt，f)进行逆Fourier变换，可得相位校正后的信号

$ {x_{\rm{p}}} = {\rm{IFT}} [S(\tau + \Delta t)] $

(10)

对校正相位后的地震信号x_p进行改进的广义S变换，在S域实现反褶积^[4]。在S域，地震信号x_p的时频振幅谱等于子波的时频振幅谱|G_w(f, τ)|与反射系数时频振幅谱|G_R(f, τ)|的乘积，即

$ |{G_{{x_{\rm{p}}}}}(f,\tau )| \approx |{G_w}(f,\tau )| \cdot |{G_R}(f,\tau )| $

(11)

假设反射系数满足白谱特性，地震记录振幅谱的锯齿波动是反射系数引起的，利用式(7)平滑每个时间点τ处的非平稳地震记录时频振幅谱值|G_xp(f, τ)|，可估计|G_w(f, τ)|，则在S域估计的反射系数时频振幅谱^[12-13]为

$ |{G_R}(f,\tau )| = \frac{{|{G_{{x_{\rm{p}}}}}(f,\tau )|}}{{|{G_w}(f,\tau )| + u(\tau )}} $

(12)

式中u(τ)为引入的调谐参数，以免分母出现零值。则最终的时间域反射系数为

$ R(t) = \sum\limits_\tau {{\rm{ IS }}} [{G_R}(f,\tau )] $

(13)

式中IS(·)为逆S变换。

2 模型数据验证

Chirp信号^[14-15]如图 1a所示，在30~75Hz频率范围呈二次函数递增变化；图 1b~图 1d分别为Gabor变换、S变换^[16-19]和改进的广义S变换(λ=1.5、r=0.8)的时频分析结果。可以看出：短时傅里叶变换存在整体频率分辨率低的问题；S变换改善了短时傅里叶变换的时频分辨率，但是随着频率的升高，其高频分辨能力明显降低；改进的广义S变换在信号处于低频和高频时，都体现出较高的时频分辨能力，且能量更聚焦，可为后续时频域反褶积奠定良好的基础。

模拟单道衰减地震记录验证相位校正对地震波形的影响以及其抗噪性能。图 2a为Q=30用零相位Ricker子波合成的衰减地震记录，随着传播时间的增大，子波波形发生畸变。对图 2a记录应用式(8)进行相位校正，结果如图 2b所示。可见相位校正后，衰减地震记录子波已零相位化，基本消除波形畸变。图 2c和图 2d分别为图 2a、图 2b记录的瞬时相位谱，校正后相位谱得以归位，且更精确描述了地震信号随时间的突然相位变化关系。对图 2a记录加入信噪比为35dB的随机噪声，如图 2e所示，图 2f是其相位校正结果。图 2g和图 2h分别为图 2e和图 2f记录的瞬时相位谱。对比图 2e和图 2f可以看出，噪声对相位校正的过程影响不大，基本不影响相位校正的稳定性。对比图 2g和图 2h可以看出，该相位校正算法抗噪能力强，但噪声无法通过相位校正去除。

建立简单模型说明本文相位校正方法可以有效改善反褶积效果。图 3a为反射系数序列，用主频为35Hz的Ricker子波与其合成的地震记录如图 3b所示^[20]，波形和相位不随时间变化。Q=40的衰减记录如图 3c所示，随着传播时间的增长，能量衰减、波形畸变^[21-23]。图 3c的衰减地震记录的相位校正结果如图 3d所示。对比图 3c和图 3d可以看出，相位校正成功地消除了地震记录传播过程中的相位畸变。分别对相位校正前、后的地震记录进行S域反褶积，结果如图 3e、图 3f所示。在图 3e红框内，由于没有校正相位，反褶积结果出现两个不正确的尖脉冲，这是由于地震子波在传播过程中相位发生畸变，而S域反褶积算法无法准确提取子波的相位谱，导致信号从S域转换到时间域出现了误差。对相位校正后的地震记录进行S域反褶积，能够很好地压缩子波旁瓣、获取准确的反射系数，与真实的反射系数具有很好地对应关系(图 3f红框所示)。单道模型试算结果表明，将相位校正与S域反褶积相结合能够得到更加准确的反褶积结果，避免了子波相位谱提取不准的缺陷。

建立二维地震模型(图 4a)验证本文的相位校正S域反褶积的正确性。该模型在浅部存在一个较小的超覆体，深部存在一个薄层。对图 4a的速度模型进行正演获得衰减地震剖面(图 4b)，其中Q值从上往下依次为30、35、50、60、80，随着传播时间的增大，地震记录的能量发生衰减，振幅减弱，薄层不清晰^[24-25]。对图 4b的正演记录分别用传统反褶积方法和本文方法处理，结果如图 4c和图 4d所示。传统的反褶积有效地补偿了深层地震波能量，同时一定程度上提高了剖面的分辨率，但整体效果不明显。

本文的相位校正S域反褶积算法可以有效地压缩地震记录旁瓣，获得分辨率更高的地震剖面。从正演剖面识别的超覆体尖灭点位于第14道；传统反褶积处理后，可在第11道识别出尖灭点位置(图 4c)，但不能识别出深部的薄层；而经过相位校正S域反褶积处理后，在第4道就可以识别出尖灭点(箭头所示)，与真实的尖灭点更接近，同时能识别出深部薄层的顶、底界(图 4d)。

3 实际资料处理

选取渤海辽西M工区地震资料(图 5a)验证相位校正S域反褶积方法在实际资料处理中的适用性。该工区目的层位于潜山顶界上，砂岩超覆现象明显，多处砂岩受潜山隆起影响，出现了同相轴能量减弱、层位消失的现象，导致后期井位设计出现困难。

如何有效地识别该地区地层实际尖灭位置，刻画地层尖灭线是亟需解决的问题。因此，将本文的反褶积方法应用于该工区，旨在恢复该地区地震波能量，提高地震剖面分辨率，达到识别该地区地层尖灭点和薄互层的目的。分别对原始数据使用传统反褶积和相位校正S域反褶积处理，结果如图 5b和图 5c所示。传统反褶积处理后，整个剖面的能量更强，改善了同相轴连续性，潜山上覆地层的弱反射能量得到了一定的恢复(图 5b箭头所示)，提高了原始地震剖面的分辨率。与传统反褶积方法处理结果(图 5b)相比，本文方法更好地恢复了潜山上覆地层反射波组的能量，且层间连续性要优于传统反褶积处理结果(图 5c箭头所示)。

为了凸显处理结果的细节部分，将图 5中的箭头处放大显示，如图 6所示。在图中红色虚线框处，原始剖面上为一套厚层，经传统反褶积处理后，压缩了该处的同相轴，但并没有识别出薄层，而本文方法处理后，有效地识别出该处的薄互层，地震剖面分辨率明显更高。在原始地震剖面上，储层弱反射在第130道消失，出现地层尖灭(红色虚线处箭头所示)。经过传统的反褶积处理后，有效将该地层尖灭点前推至第118道处，而经过本文方法处理后，在108道可见该地层反射，能更准确地识别出地层尖灭点。可见，本文的相位校正S域反褶积方法能有效提高地震资料纵、横向分辨率，为精确刻画岩性圈闭提供有力支撑。

为了展示高分辨处理前、后时频能量变化，分别从图 6中抽取的第120道地震记录进行时频分析，结果如图 7所示。对比原始信号(图 7a)与两种反褶积方法处理结果时频谱(图 7b和图 7c)可见，本文方法处理效果好，在保持强信号能量的基础上恢复了弱信号的能量(红色箭头所示)；主频(图中白线为主频随时间的变化)明显提升，大幅提高了地震记录分辨率。图 8为图 5剖面的振幅谱，原始数据、传统反褶积方法和本文方法处理结果的有效频带范围分别为4~50Hz(图 8a)、10~100Hz(图 8b)和1~140Hz(图 8c)，可见本文方法处理后，在保留原始数据的低频信息基础上，有效拓宽了地震资料的频带宽度，结果更优。

4 结束语

本文在相位校正基础上进行时频域反褶积处理，形成了一套适用性广的高分辨处理方法，并通过模型和实际资料处理取得了良好的应用效果，得到以下几点认识：

(1) 针对数据不同频率成分能量不同的特点，对S变换的窗函数进行扩展，可以提高其对信号不同频率成分的识别能力，达到灵活时频分析的目的。

(2) 在反褶积的处理过程中，本文结合了传统的反Q滤波中的相位校正方法，对地震数据进行相位校正后，可以更准确地将时频域反褶积结果转换到时间域中，使反褶积的效果更好，薄层识别能力更强。

(3) 经本文的相位校正S域反褶积处理后，能有效拓宽原始地震资料的频带、提高分辨率。