王清振 姜秀娣 翁 斌 刘志鹏 印海燕 陈剑军

（中海油研究总院 北京 100028）

王清振，姜秀娣，翁斌，等.基于变时窗Gabor变换的高分辨率处理技术研究与应用［J］.中国海上油气，2015，27（6）：19-26.

目前常用的提高地震资料纵向分辨率的方法，如最小平方反褶积、预测反褶积、变模反褶积、同态反褶积、最小熵反褶积、基于互信息率准则的盲反褶积以及谱白化等，其理论基础是传统的褶积模型［1-3］，该模型的基本假设之一就是子波在地下传播过程中不随时间变化。然而，实际地震子波常常是时变的，这就使得以该模型为理论基础的提高分辨率的方法在许多情况下难以取得好的效果。为此，一些学者提出了反射地震记录的另一种模型［4-5］，认为子波在地下传播过程中随时间发生变化，不同时刻到达检波器的子波的波形是不同的，反射地震记录是这些具有不同到达时间的子波的叠加，因此这种模型被称为反射地震记录变子波模型。

地震记录的子波时变性主要是由波前扩散和地层吸收衰减效应引起的［6］，其中波前扩散可以用几何扩散函数来校正，而地层吸收衰减引起的子波时变效应（即所谓的Q效应）常采用反Q滤波和时变谱白化方法进行补偿［7-8］。但是，地下介质构造复杂导致Q值往往难以求准，从而影响反Q滤波的准确程度，而时变谱白化方法的参数需要人为选择，对补偿效果影响较大，并且处理后的数据的振幅相对关系难以描述。2001年，G F Margrave等［9］提出了一种Gabor反褶积方法，将地层视为单Q值的均匀黏弹性介质，且仅考虑吸收衰减条件下直接将Wiener反褶积算法扩展到待分析信号为时变的情况。这种方法不需要预先估计Q值就可完成反Q滤波，但对平滑窗依赖性较强，而且经这种方法处理后的地震记录不能很好地刻画反射系数的局部能量相对关系（直观地讲会产生等幅效应，很像经自动增益控制（AGC）处理后的结果），在某些情况下不利于地震资料解释。2002年，G F Margrave等［10］对上述方法进行了改进，用双曲型平滑代替矩形窗平滑，压制了AGC效应，并于2005年系统地阐述了Gabor反褶积方法的原理［11］。由于Gabor反褶积方法采用均匀Gabor标架分析信号，分析窗的长度难以适应待分析信号各段的情况，并且按照单Q值模型同时求解不同时刻的反褶积算子，使得该方法不能用于地层Q值随深度变化的情况。2009年，高静怀 等［12］提出了基于地震记录变子波模型提高地震资料分辨率的方法，但该方法是在单时窗内采用维纳反褶积方法进行拓频，没有考虑不同深度时窗由吸收引起的振幅衰减影响。

笔者在高静怀 等［12］方法的基础上，考虑地层吸收的振幅衰减影响，基于反射地震记录变子波模型与Gabor变换，在地震数据地层结构的约束下构造一组自适应于地震记录的分析时窗，对地震记录进行自适应非均匀划分，在时频域利用非线性压缩映射提取时变子波，对非平稳地震记录做频谱拓宽和振幅校正，估计补偿因子，对不同时窗之间由吸收引起的能量衰减进行补偿，反变换回时间域得到相对保幅的高分辨率处理结果。本文方法在渤海某油田进行了实际应用，显著提高了该油田地震资料的分辨率，并在波阻抗反演、层序检测等方面发挥了重要作用。

1 基于变时窗Gabor变换的高分辨率处理技术

本文提出的基于变时窗Gabor变换的高分辨率处理技术包括自适应时窗分解、变时窗Gabor正变换、时变子波提取、高分辨率处理和变时窗Gabor反变换等5个步骤，其中最关键的3个步骤是自适应时窗分解、时变子波提取和高分辨率处理。

1.1 自适应时窗分解

1）Gabor变换。

如果一个函数集｛Ψj∶j∈Z｝的和为1，即

则此函数集构成对单位1的分割［11］。令φ（t）满足

则集合｛φm（t）∶m∈Z｝构成一个单位1的分割，取分析时窗函数g（t）=φ（t）u，综合时窗函数γ（t）=φ（t）1-u，其中0≤u≤1，可以证明由｛gm∶m∈Z｝和｛γm∶m∈Z｝能够构造一对Gabor标架。取u=1，则g（t）=φ（t），γ（t）=1，对于待分析信号s（t），相应于上述Gabor标架的Gabor变换为

Gabor逆变换为

可以看出，这种Gabor变换方法是通过设置u=1，使得反变换所用的综合时窗为单位1，不需要额外计算综合时窗，并且式（3）、（4）可以通过快速Fourier变换实现，计算效率较高。为了表述直观，这里用连续变量表示时间和频率，用离散变量表示各时窗的中心在时间轴上的位置。对于满足一定条件的和不为1的函数集，也可以通过归一化的方法将其变为对单位1的分割。

2）构造分析时窗。

由于Gabor变换在每个时间采样点都有一个分析时窗，冗余度很高，因此如果在Gabor域逐窗处理数据，则计算量会很大。而且如果选择的时窗太窄，则频率分辨率低，窗内的时变子波无法提取；如果为了提高频率分辨率，将窗口变宽，则窗内平稳假设的近似程度可能会变差。针对上述问题，提出了构造自适应于非平稳地震记录的分析时窗的方法，具体步骤为：

①提取包络峰值。地震道的包络峰值与反射界面有一定的对应关系，在一定程度上可以大致反映地层的层序结构，因此，如果使用地震道的包络峰值约束分子窗的构造，那么构成的分子窗在横向上将与地层结构有关，有利于保持处理后资料的横向连续性（图1）。根据复道分析原理，设s＊（t）为s（t）的Hilbert变换，则s（t）的包络为

图1 实际地震剖面及其包络峰值剖面Fig.1 Field seismic section and peak envelope section

②生成满足单位分解的原子窗族。在每个采样点建立的分析时窗被称为原子窗，所有采样点对应的分析时窗的集合被称为原子窗族。选择基本原子窗函数G（t）为Lamoureux函数，N为地震道采样点个数，用G（t）=G（t-jΔt）表示中心位于第j个采样点上的原子窗，按式（6）对原子窗族｛Gj（t）∶1≤j≤N｝进行归一化［13］处理，由式（1）可得到一组满足单位分解的原子窗族｛gj（t）∶1≤j≤N｝，即

③构造初始分子窗。在包络的控制下将地震记录自适应地分成多个时间段，通过对包络求两级差分可以得到这些时间段的分割点，然后将相邻分割点间的小原子窗叠加起来就得到了与各个分段对应的分析时窗（形象的被称为分子窗）。设第k个分子窗的起点和终点分别为Mk-1和Mk，则第k个分子窗为

④分子窗能量归一化。第k个分子窗的能量为

则能量归一化后第k个分子窗为Ψk（t）／Ek。

图2给出了自适应时窗的分解过程。

1.2 时变子波提取

在很多情况下（如爆炸震源等），地震记录震源子波的振幅谱是频率的单峰函数。对于子波振幅谱非单峰函数的地震记录，也可以将其转化为单峰函数的情况。假定地震记录中地震子波的振幅谱是频率的单峰、光滑函数，可以构造一个压缩算子，把提取子波振幅谱的问题转化为求解该算子的不动点问题，从而根据不动点迭代算法从地震记录中提取子波振幅谱。

设x0∈（a，b），取δ0＞0足够小，对于任意给定的0＜δ≪1，若

则Ωx0，δ0，δ是L2［a，b］空间中的子集，由于区间［a，b］是有限的，可以推导出L2［a，b］⊂Lp［a，b］，0＜p＜1，

图2 自适应时窗的分解过程Fig.2 Adaptive time-window decomposition process

显然，0≤F（ŝ；f）≤1。

设c＞0，α＞1，β＞1，任 取ŝ∈Ωx0，δ0，δ，定 义Ωx0，δ0，δ上的算子P为

显然，δ≤P（ŝ；f）≤c＋δ，可以推导出P（ŝ；f）∈L2［a，b］，因此，P是Ωx0，δ0，δ⊂L2［a，b］到L2［a，b］的非线性算子。由Banach空间算子不动点定理可以证明，P是Ωx0，δ0，δ到Ωx0，δ0，δ的压缩算子，并且在Ωx0，δ0，δ中存在唯一不动点。

设ŝ（f）为地震记录的振幅谱，f为频率，fcut为振幅谱的截止频率，0＜p＜1，设迭代初值为

因此，Ωx0，δ0，δ⊂Lp［a，b］。

任取ŝ∈Ωx0，δ0，δ，令

其中压缩算子P的参数由最小二乘得到。因此，通过建立如下迭代

最终得到不动点Faim，则估计的子波振幅谱｜ŝw（f）｜为F1／paim（下文称为F函数）。

图3给出了图2中合成记录时变子波振幅谱的估计效果，其中反射系数采样点数为500，采样间隔2 ms，子波为50 Hz雷克子波。从图3可以看出，估计的子波振幅谱和实际子波振幅谱吻合很好。

图4展示了在图2中各分子窗内提取的时变子波的振幅谱，可以看出，受地层吸收的影响，随着深度增加，子波振幅谱主频逐渐降低。因此，下一步的高分辨率处理就是要把这些传播过程中损失的信息补偿回来。

1.3 高分辨率处理

各分子窗对应的振幅谱既含有等效子波的信息，也含有该段反射系数序列的信息。对于待分析的非平稳地震记录来说，分子窗都是自适应于地震记录构造的，因此可认为各分子窗内片段的F函数近似等于该段的等效子波的振幅谱乘以一个常数。等效子波的振幅谱引起该段地震记录频带变窄，使分辨率降低，因此从该段地震记录中剔除等效子波的效应就可提高其分辨率。

图3 图2中合成记录时变子波振幅谱估计效果Fig.3 Time-varying wavelet amplitude spectrum estimation of the synthetic signal in Fig.2

图4 图2中合成记录各分子窗内提取的时变子波的振幅谱Fig.4 Amplitude spectrum extracted in each molecules window of the synthetic signal in Fig.2

1）不带衰减补偿的拓频。

以第k段为例，为了提高该段的分辨率，采用零相位Wiener反褶积算法拓宽该段的频带，即

式（14）中：ŝk（f）为原始数据的频谱；ŝBk（f）为数据拓频后的频谱；｜ŝwk（f）｜为第k段对应的F函数，它与第k段的等效子波振幅谱＾wk（f）相似；Amax为｜ŝwk（f）｜的最大值；ε是一个小于1的无量纲常数。

由于各片段对应的F函数和等效子波的振幅谱之间存在一个比例系数，这个比例系数和该段内反射系数的能量有关。也就是说，｜ŝwk（f）｜中既包含第k段等效子波的能量，也包含该段内反射系数的能量，因此，式（14）的结果会产生类似于AGC的等幅效应［10］。为了让得到的高分辨结果的振幅有意义，采用式（15）校正ŝBk（f）的能量，即

2）带衰减补偿的拓频。

对于第k个分子窗截出的地震记录片段，把从参考子波（震源子波）到第k个分子窗之间的介质视为均匀黏弹性介质，介质的等效品质因子记为Qk。令参考子波从震源传播到中心为k的分子窗处所用的时间为Tk，则平面波在频率域满足因果律的传播算子可表示为［14-16］

式（16）中：H为在某个时刻t对频率f的Hilbert变换。

用（f）表示参考子波的Fourier频谱，用（f）表示第k段上等效子波的Fourier频谱，根据波的传播理论有

将式（16）代入式（17）中，且仅考虑振幅部分，可得

由于本文构造的分子窗所覆盖的信号段近似平稳，所以每个窗截取的信号段有一个近似不变的等效子波，该段的F函数可表示为该等效子波的振幅谱乘以一个常数。记（f）的F函数为｜（f）｜，Ck为第k段的常数比例因子，则

将式（18）代入式（19）中，可得

对式（20）两边取对数，采用对数谱比值法可计算出Qk，从而得到｜h（f，Tk）｜。令

用βk乘以ŝk（f），得到补偿后的片段为

记（f）的F函数为｜（f）｜，将式（14）和式（15）分别作用于（f），可得拓频结果为

第k段对应的带衰减补偿的拓频结果为

图5 图2中合成记录高分辨率处理效果Fig.5 High resolution processing result of the synthetic signal in Fig.2

对比图5可以看出，反射系数①②③④（图5a）在合成地震记录（图5b）上不清晰。反射系数①②在谱白化处理结果（图5c）中比在合成地震记录上清晰，而反射系数③④在谱白化处理结果中仍不清晰。在图5d、e中，反射系数①②③④都得到了很好的刻画，图5d中的高分辨地震记录保持了图5b原始地震记录的相对能量关系；图5e中的高分辨地震记录恢复了衰减的能量，与图5a中的反射系数序列有较好的对应关系。这表明，该技术不仅能够提高地震资料分辨率，而且能够补偿由地层吸收引起的能量衰减。

2 应用效果

图6、7是渤海A油田提高分辨率前后地震剖面及频谱对比图，可以看出，提高分辨率后剖面波组特征保持得较好，频谱的主频和带宽都得到了提升。

图8是渤海A油田提高分辨率处理前后波阻抗反演结果对比图，可以看出，提高分辨率后反演结果和井的吻合度更高，波阻抗资料的分辨率也有了较大提高，反演精度更高。

图6 渤海A油田提高分辨率前后地震剖面对比Fig.6 Comparison between field data section and high resolution section in A oilfield，Bohai sea

图7 渤海A油田提高分辨率前后地震资料频谱对比Fig.7 Comparison of amplitude spectrum between field data and high resolution data in A oilfield，Bohai sea

图8 渤海A油田提高分辨率前后波阻抗反演结果对比Fig.8 Comparison of impedance inversion result between field data and high resolution data in A oilfield，Bohai sea

图9是渤海A油田提高分辨率前后层序检测（瞬时频率属性）对比图，可以看出，提高分辨率后红圈范围内的几套三角洲砂体的叠置关系更加明显，砂体边界也更加清晰。

图9 渤海A油田提高分辨率前后层序检测（瞬时频率属性）对比图Fig.9 Comparison of sequence detection（instantaneous frequency）between field data and high resolution data in A oilfield，Bohai sea

3 结束语

本文提出的基于变时窗Gabor变换的高分辨处理技术克服了时变谱白化等方法中窗长难以确定和等幅效应的问题，以及反Q滤波等方法中Q值难以求准的问题，理论分析表明其适用条件更接近地下实际情况，是一种相对保幅的提高地震资料分辨率的新方法。实际应用表明，该方法显著提高了地震资料的分辨率，能够在波阻抗反演、层序检测等领域发挥重要作用。

［1］YILMAZ O.Seismic data processing［M］.Tulsa：Society of Exploration Geophysicists，1987：49-61，85-94.

［2］孙雷鸣，万欢，陈辉，等.基于广义S变换地震高分辨率处理方法的改进及在流花11-1油田的应用［J］.中国海上油气，2011，23（4）：234-237.Sun Leiming，Wan Huan，Chen Hui，et al.An improved method of seismic high-resolution processing based on generalized S transform and its application in LH11-1 oilfield［J］.China Offshore Oil and Gas，2011，23（4）：234-237.

［3］陈宝书，汪小将，李松康，等.海上地震数据高分辨率相对保幅处理关键技术研究与应用［J］.中国海上油气，2008，20（3）：162-166.Chen Baoshu，Wang Xiaojiang，Li Songkang，et al.Key technique research and application in relative amplitude preservation processing with high resolution for offshore seismic data［J］.China Offshore Oil and Gas，2008，20（3）：162-166.

［4］RICKER N.The form and nature of seismic wavelets and the structure of seismograms［J］.Geophysics，1940，5（4）：348-366.

［5］DOBRIN MB，SAVIT C H.Introduction to geophysical prospecting［M］.New York：McGraw Hill Book Co.，1988：200-220.

［6］ZIOLKOWSKI A.Why don′t we measure seismic signatures？［J］.Geophysics，1991，56（2）：190-201.

［7］HALE I D.Q adaptive deconvolution［R］.Palo Alto：Stanford Exploration Project Report Number 30，1982：133-158.

［8］GIBSON B，LARNER K L.Comparison of spectral flattening techniques［R］.Western Geophysical Company，1982：10-30.

［9］MARGRAVE G F，LAMOUREUX MP.Gabor deconvolution［R］.Calgary：CREWES Research Report，2001：241-276.

［10］MARGRAVE G F，HENLEY D C，LAMOUREUX MP，et al.An update on Gabor deconvolution［R］.Calgary：CREWES Research Report，2002：11-27.

［11］MARGRAVE G F，GIBSON P C，GROSSMAN J P，et al.The Gabor transform，pseudo differential operators，and seismic deconvolution［J］.Integrated Computer-Aided Engineering，2005，12（3）：43-45.

［12］高静怀，汪玲玲，赵伟.基于反射地震记录变子波模型提高地震记录分辨率［J］.地球物理学报，2009，52（5）：1289-1300.Gao Jinghuai，Wang Lingling，Zhao Wei.Enhancing resolution of seismic traces based on the changing wavelet model of the seismogram［J］.Chinese Journal of Geophysics，2009，52（5）：1289-1300.

［13］GROSSMAN J P，MARGRAVE G F，LAMOUREUX MP.Constructing adaptive，nonuniform Gabor frames from partitions of unity［R］.Calgary：CREWES Research Report，2002：1-10.

［14］WANG Yanghua.A stable and efficient approach of inverse Q filtering［J］.Geophysics，2002，67（2）：657-664.

［15］WANG Yanghua.Inverse Q-filter for seismic resolution enhancement［J］.Geophysics，2006，71（3）：51-61.

［16］余连勇，胡光义，赵岩，等.稳定的反Q滤波统一算法及其在地震资料高分辨率处理中的应用［J］.中国海上油气，2014，26（4）：29-33.Yu Lianyong，Hu Guangyi，Zhao Yan，et al.A unified algorithm of stable inverse Q filtering and its application to highresolution processing of seismic data［J］.China Offshore Oil and Gas，2014，26（4）：29-33.