CN102183787A - 一种基于地震记录变子波模型提高地震资料分辨率的方法 - Google Patents

Abstract

本发明涉及一种基于地震记录变子波模型提高地震资料分辨率的方法，其特征在于，包括以下步骤：1)构造分子-Gabor窗，将非平稳地震记录自适应地划分为若干平稳的片段，每个片段中拥有一个近似不变的等效子波，且所述等效子波易于从该片段中提取出来；2)用第1)步构造得到的分子窗生成分子-Gabor标架，将非平稳地震记录变换到分子-Gabor域；3)在分子-Gabor域，对每个分子-Gabor窗内的地震记录片段所对应的分子-Gabor系数，进行拓频和能量补偿处理；4)将处理后的分子-Gabor系数反变换到时间域，得到提高分辨率后的地震记录。该方法以现代拟微分算子理论为基础，以自适应时-频分析方法为工具，处理后的地震资料具备高分辨率和相对保持振幅特性。

Description

一种基于地震记录变子波模型提高地震资料分辨率的方法

技术领域

本发明属于油气勘探中地震资料分析与处理领域，涉及该领域的一种基于自适应时-频变换的提高地震资料分辨率处理的技术，具体地说，是关于一种基于地震记录变子波模型提高非平稳地震资料分辨率的方法。

背景技术

常用的提高地震资料纵向分辨率的方法(例如谱白化方法、各种反褶积方法)，其理论基础是传统的褶积模型。这一模型基于若干基本假设，其中之一就是假设子波是平稳的，即假设子波在地下传播过程中不随时间变化。然而，实际中子波常常是非平稳的，这使得以该模型为理论基础的提高分辨率的方法，在许多情况下难以取得好的效果。为此，一些学者提出了反射地震记录的另一种模型，这种模型认为子波在地下传播过程中随着传播时间发生变化，不同时刻到达检波器的子波的波形是不同的，反射地震记录是这些具有不同到达时的子波的叠加。这种模型被称为非平稳地震记录模型。地震记录的非平稳性主要是由波前扩散和频率衰减效应引起的，波前扩散可以用几何扩散函数来校正。由频率衰减引起的非平稳效应，即所谓的Q效应。反Q滤波和时变谱白化方法是常用的频率衰减补偿方法。但是，由于Q值较难求准，而限制了反Q滤波方法的应用；而时变谱白化方法处理的结果难以确保地震记录局部能量相对关系。将地层视为单Q值的黏弹性介质，Margrave和Lamoureux(Margrave，G.F.，M.P.Lamoureux，2001，“Gabor deconvolution，”CREWES Research Report，vol 13，pp.241-276，pp.252-253)给出了一个非平稳地震道模型的显式表达式。以这一表达式为基础，在假定反射系数序列满足白噪假设的前提下，他们又提出了一种Gabor反褶积方法，直接将Wiener反褶积算法扩展到待分析信号为非平稳的情况。然而，由于按照单Q值模型同时求解不同时刻的反褶积算子，使得Gabor反褶积方法对地层Q值随深度变化的情况效果不佳。

发明内容

基于现有技术中所存在的上述问题，本发明目的是提出一种基于地震记录变子波模型提高地震资料分辨率的方法。

为实现上述发明目的，本发明采用了以下技术方案：一种基于地震记录变子波模型提高地震资料分辨率的方法，其特征在于，包括以下步骤：

1)构造分子-Gabor窗，将非平稳地震记录自适应地划分为若干平稳的片段，每个片段中拥有一个近似不变的等效子波，且所述等效子波易于从该片段中提取出来；

2)用第1)步构造得到的分子-Gabor窗生成分子-Gabor标架，将非平稳地震记录变换到分子-Gabor域；

3)在分子-Gabor域，对每个分子-Gabor窗内的地震记录片段所对应的分子-Gabor系数进行拓频和能量补偿处理；

4)将处理后的分子-Gabor系数反变换到时间域，得到提高分辨率后的地震记录。

所述步骤1)中，分子-Gabor窗通过单位分解法构造，具体分为三个步骤：

①生成满足单位分解的原子窗族

选择基本原子窗函数G(t)为Lamoureux函数，用G_j(t)＝G(t-jΔt)表示中心位于第j个采样点上的原子窗，对原子窗族{G_j(t)：1≤j≤N}按下式归一化：

$g_j\left(t\right)=G_j\left(t\right)/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{G}_i\left(t\right)---\left(1\right)$

N为地震道采样点的个数，由(1)式可得一组满足单位分解的原子窗族{g_j(t)：1≤j≤N}；

②构造初始分子-Gabor窗

首先定义地震信号s(t)的瞬时频率如下：

f(t)＝1/2π{[s(t)ds^*(t)/dt-s^*(t)ds(t)/dt]/[a(t)²+ξ²]}(2)

这里s^*(t)为s(t)的Hilbert变换，a(t)＝[s(t)²+s^*(t)²]^1/2为s(t)的包络，ξ为无量刚的小常数；

其次用加权瞬时频率f′(t)代替局部均值频率来表征地层的吸收效应，加权瞬时频率定义如下：

$f^{\′}\left(t\right)={\∫}_{t-T}^{t+T}f\left(t\right)W\left(t\right)\mathrm{dt}/{\∫}_{t-T}^{t+T}W\left(t\right)\mathrm{dt}---\left(3\right)$

这里W(t)取为s(t)包络幅值的平方；

然后，采用保边缘平滑拟合的方法消除子波干涉对加权瞬时频率的影响：选取信号包络峰值处的加权瞬时频率，对它们作保边缘平滑拟合，检验出偏离拟合曲线较远的数据点，将这些较远的数据点去掉，然后用线性插值去填充这些去掉的点，如此反复迭代直到拟合曲线上的值不再有大的变化，从而得到最终拟合的曲线；

最后根据拟合曲线上相邻两点之间的差值和距离计算各分子窗的分割点，然后将相邻分割点间的小原子窗叠加起来就得到了自适应分子窗，设第k个分子窗的起点和终点分别为M^k-1和M^k，则第k个分子窗可以表示为：

${\mathrm{\ψ}}_k\left(t\right)=\underset{j=M^{k-1}}{\overset{M^k-1}{\mathrm{\Σ}}}{g}_j\left(t\right)---\left(4\right)$

③对初始分子-Gabor窗进行能量归一化，得到分子-Gabor窗

令E_k表示第k个分子窗的能量，则有

$E_k=\left|\right|{\mathrm{\ψ}}_k\left|\right|={\left[{\∫}_{-\∞}^{+\∞}{\left|{\mathrm{\ψ}}_k\left(t\right)\right|}^2\mathrm{dt}\right]}^{1/2}---\left(5\right)$

能量归一化以后第k个分子窗ψ_k(·)/E_k即为我们所构造的分子-Gabor窗。

所述步骤2)中，由步骤1)中得到的分子-Gabor窗构造分子-Gabor标架，相应的分子-Gabor变换定义如下：

${\hat{s}}_k\left(f\right)=\frac{1}{E_k}\underset{-\∞}{\overset{\∞}{\∫}}s\left(t\right){\mathrm{\ψ}}_k\left(t\right)e^{-2\mathrm{i\πft}}\mathrm{dt}---\left(6\right)$

将非平稳地震记录s(t)按照(6)式变换到分子-Gabor域。

所述步骤3)中，在分子-Gabor域对每个地震记录片段进行频带拓宽和振幅校正，方法之一是保持原始地震记录相对能量关系的拓频：

对于第k个分子窗截出的地震记录片段，采用下式拓宽该片段的频带：

${\hat{s}}_k^h\left(f\right)={\hat{s}}_k\left(f\right)/\left(\right|{\hat{s}}_k^w\left(f\right)|+\mathrm{\ϵ}{A}_{\max })---\left(7\right)$

式中ε是一个小于1的无量纲的常数， 

是第k个片段的F-function函数，它与第k个片段的等效子波的振幅谱相似，A_max为 

的最大值，εA_max用来增强计算的稳定性，采用

${\hat{s}}_k^{\′}\left(f\right)=\left|\right|{\hat{s}}_k\left|\right|{\hat{s}}_k^h\left(f\right)/\left|\right|{\hat{s}}_k^h\left|\right|---\left(8\right)$

来校正 

的能量， 

即为第k个地震记录片段对应的保持原始地震记录相对能量关系的拓频结果。

所述步骤3)中，在分子-Gabor域对每个地震记录片段进行频带拓宽和振幅校正，方法之二是带衰减补偿的拓频：

对于第j个分子窗截出的地震记录片段，把由参考子波到第j个分子之间的介质视为均匀粘弹性介质，介质的等效品质因子记为Q_j，令参考子波从震源传播到中心为j的分子窗处所用的时间为T_j，则平面波在频率域满足因果律的传播算子表示为

$h(f,T_j)=e^{-\mathrm{\πf}{T}_j/Q_j}{e}^{\mathrm{iH}[\mathrm{\πf}{T}_j/Q_j]}---\left(9\right)$

H为在某个时刻t对频率f的Hilbert变换， 用 表示第j个片段上的等效子波的Fourier频谱，根据波传播理论有：

${\hat{w}}_j\left(f\right)=h(f,T_j)\hat{w}\left(f\right)---\left(10\right)$

为参考子波的Fourier频谱，把(9)式代入(10)式，且仅考虑振幅部分得到：

$\left|{\hat{w}}_j\left(f\right)\right|=e^{-\mathrm{\πf}{t}_j/Q_j}\left|\hat{w}\left(f\right)\right|---\left(11\right)$

由于构造的分子窗所覆盖的信号段近似平稳，所以每个窗截取的信号段有一个近似不变的等效子波，该等效子波可用该段的F-function(拟合函数)乘以一个常数来表示，记 的F-function为 

则有：

$\left|{\hat{s}}_j^w\left(f\right)\right|\≈C_j\left|{\hat{w}}_j\left(f\right)\right|---\left(12\right)$

这里C_j为第j个片段的常数比例因子，由(11)和(12)式可得：

$\left|{\hat{s}}_j^w\left(f\right)\right|\≈C_j{e}^{-\mathrm{\πf}{T}_j/Q_j}\left|\hat{w}\left(f\right)\right|---\left(13\right)$

对(13)式两边取对数，采用对数谱比值法可计算出Q_j，从而可得到|h(f，T_j)|，

令

β_j(f)＝1/|h(f，T_j)|(14)

用β_j乘以 

得到补偿后的片段：

${\hat{s}}_j^c\left(f\right)={\hat{s}}_j\left(f\right){\mathrm{\β}}_j\left(f\right)---\left(15\right)$

令 

表示 

的F-function，将(7)式和(8)式作用于 

可得拓频结果

${\hat{s}}_j^{c,h}\left(f\right)={\hat{s}}_j^c\left(f\right)/\left(\right|{\hat{s}}_j^{c,w}\left(f\right)|+\mathrm{\ϵ}{A}_{\max })---\left(16\right)$

和能量校正后的结果

${\hat{s}}_j^{\′\′}\left(f\right)=\left|\right|{\hat{s}}_j^c\left|\right|{\hat{s}}_j^{c,h}\left(f\right)/\left|\right|{\hat{s}}_j^{c,h}\left|\right|---\left(17\right)$

即为第j个地震记录片段对应的带衰减补偿的拓频结果。

所述步骤4)中，分子-Gabor逆变换定义如下：

$s\left(t\right)=\underset{k=1}{\overset{L}{\mathrm{\Σ}}}\left[E_k\underset{-\∞}{\overset{\∞}{\∫}}{\hat{s}}_k\left(f\right)e^{2\mathrm{i\πft}}\mathrm{df}\right],---\left(18\right)$

对(8)式处理的结果做分子-Gabor逆变换，得到保持原始地震记录相对能量关系的高分 辨地震道。对(17)式处理的结果做分子-Gabor逆变换，得到带衰减补偿的高分辨地震道。

本发明由于采取了以上技术方案，其具有的有益效果是：本发明给出了一种基于地震记录变子波模型提高非平稳地震资料分辨率的方法，该方法以现代拟微分算子理论为基础，以自适应时-频分析方法为工具，处理后的地震资料具备高分辨率和相对保持振幅特性。将该方法用于处理非平稳合成地震记录，结果表明该技术较之普通的谱白化技术不仅能够更好地保持原始记录的相对能量关系，而且能够恢复衰减的高频能量，且该技术能够更好地适用于地层Q值随深度变化的情况。将该方法用于处理叠加后三维实际地震数据，结果也表明该技术能够有效拓宽非平稳地震记录的频带，并且恢复被吸收的能量。此外，在高分辨率处理后数据上得到的波阻抗反演结果比在处理前数据上得到的波阻抗反演结果的分辨率高，进一步证明了该技术的有效性。

附图说明

图1是构造分子窗的自适应划分原理图；

图2是合成地震记录算例图；

图3是实际地震资料高分辨处理前后与合成记录的对比图；

图4(a)是实际地震资料高分辨处理前振幅谱图；

图4(b)是实际地震资料高分辨处理后振幅谱图。

具体实施方式

下面结合附图和实施例对本发明进行详细的说明。

本发明提出的基于地震记录变子波模型，在自适应时-频域提高非平稳地震资料分辨率的方法，包括以下步骤：

1)构造分子-Gabor窗，将非平稳地震记录自适应地划分为若干平稳的片段，每个片段中拥有一个近似不变的等效子波，且该等效子波易于从该片段中提取出来；

2)用第1步构造得到的分子窗生成分子-Gabor标架，将非平稳地震记录变换到分子-Gabor域；

3)在分子-Gabor域，对每个分子-Gabor窗内的地震记录片段所对应的分子-Gabor系数进行拓频和能量补偿处理；

4)将处理后的分子-Gabor系数反变换到时间域，得到提高分辨率后的地震记录。

所述步骤1)中，分子-Gabor窗可以通过单位分解法构造，具体可分为三个步骤：

①生成满足单位分解的原子窗族

恰当地选择基本原子窗函数G(t)，本发明G(t)选为Lamoureux函数(Margrave，G.F.，Peter C.Gibson，Jeff P.Grossman，David C.Henley，Victor Iliescu，and  Michael P.Lamoureux，2005，The Gabor transform，pseudodifferential operators，and seismic deconvolution：Integrated Computer-Aided Engineering，12：43-45)，用G_j(t)＝G(t-jΔt)表示中心位于第j个采样点上的原子窗。对原子窗族{G_j(t)：1≤j≤N}按下式归一化：

$g_j\left(t\right)=G_j\left(t\right)/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{G}_i\left(t\right)---\left(1\right)$

这里N为地震道采样点的个数。由(1)式可得一组满足单位分解的原子窗族

{g_j(t)：1≤j≤N}。

②构造初始分子-Gabor窗

地层吸收使得地震子波主频随着传播时间慢慢变小，这一特性可由信号的局部均值频率来描述。为了节省计算量，本发明采用加权瞬时频率f′(t)代替局部均值频率来表征地层的吸收效应。地震信号s(t)的瞬时频率定义如下：

f(t)＝1/2π{[s(t)ds^*(t)/dt-s^*(t)ds(t)/dt]/[a(t)²+ξ²]}(2)

其中s^*(t)为s(t)的Hilbert(希尔波特)变换，a(t)＝[s(t)²+s^*(t)²]^1/2为s(t)的包络，ξ为无量刚的小常数。则定义加权瞬时频率如下：

$f^{\′}\left(t\right)={\∫}_{t-T}^{t+T}f\left(t\right)W\left(t\right)\mathrm{dt}/{\∫}_{t-T}^{t+T}W\left(t\right)\mathrm{dt}---\left(3\right)$

当W(t)取为s(t)包络幅值的平方时，可以证明当T逐渐变大时，f′(t)趋于信号s(t)的局部均值频率。

实际上，信号加权瞬时频率随时间的变化除了与地层吸收有关外，还与窗内子波干涉有关，为此，本发明进一步采用保边缘平滑拟合的方法(AlBinHassan，N.M.，Luo，Y.，and Al-Faraj，M.N.，July-August 2006，3D edge-preserving smoothing and applications：Geophysics，71(4)，5-11)，消除子波干涉对加权瞬时频率的影响。保边缘平滑拟合的方法是用一个矩形窗在信号的目标点周围滑动形成一组矩形窗模板，然后选择与目标点最同类(模板内的点的值与目标点的值的均方差最小)的模板内那些点的均值代替目标点，从而实现对信号的保边缘平滑。可见，该方法可以将目标点和其附近最同类的点归为一类，据此我们可以用该方法对包络峰值处的加权瞬时频率进行平滑分类。

图1为构造分子窗的自适应划分原理图。其中图1(a)是信号s(t)(实线)及其包络(点线)和包络峰值(*)图；(b)是信号s(t)的瞬时频率(实线)、包络峰值处 的瞬时频率(*)和保边缘平滑拟合后的拟合点(○)关系图，以及依据拟合点随时间的衰减量和它们之间的间距确定的分界点(+)；(c)是将分界点(+)间的原子窗叠加起来形成的初始分子-Gabor窗，这些初始分子-Gabor窗的叠加和为1；(d)是对初始分子-Gabor窗进行能量归一化，得到的分子-Gabor窗；(e)是用(d)所示的分子-Gabor窗构造分子-Gabor标架对s(t)做分子-Gabor变换得到的时频振幅谱图，图右侧为振幅谱值对应的色标，深色代表大幅值；(f)是由分子-Gabor变换得到的时频系数重构得到的信号；(g)是重构误差，误差数量级为10^-15，可见，分子-Gabor变换方法的重构精度很高。

如图1(b)所示，选取信号包络峰值处的加权瞬时频率(*处所示)，对它们作保边缘平滑拟合，检验出偏离拟合曲线较远的数据点，将它们去掉以便消除子波干涉对加权瞬时频率的影响，然后用线性插值去填充去掉的点，如此反复迭代直到拟合曲线上的值不再有大的变化，从而得到图1(b)中的“○”点。根据两个“○”点间的差值和“○”点之间的距离计算各窗的分割点，如图1(b)中“+”所示，将两个“+”间的相邻小原子窗叠加起来就得到了图1(c)所示的自适应分子窗。通过保边缘平滑拟合的方法对提取出的瞬时频率进行筛选和分类，使得每个分子窗所覆盖的信号段是近似平稳的。设第k个分子窗的起点和终点分别为M^k-1和M^k，则第k个分子窗可以表示为：

${\mathrm{\ψ}}_k\left(t\right)=\underset{j=M^{k-1}}{\overset{M^k-1}{\mathrm{\Σ}}}{g}_j\left(t\right)---\left(4\right)$

③对初始分子窗进行能量归一化，得到分子-Gabor窗

令E_k表示第k个分子窗的能量，则有

$E_k=\left|\right|{\mathrm{\ψ}}_k\left|\right|={\left[{\∫}_{-\∞}^{+\∞}{\left|{\mathrm{\ψ}}_k\left(t\right)\right|}^2\mathrm{dt}\right]}^{1/2}---\left(5\right)$

能量归一化以后第k个分子-Gabor窗为ψ_k(·)/E_k。

所述步骤2)中，由步骤1中得到的分子-Gabor窗ψ_k(·)/E_k构造分子标架ψ_k(t)e^2iπft/E_k，相应的分子-Gabor变换定义如下：

${\hat{s}}_k\left(f\right)=\frac{1}{E_k}\underset{-\∞}{\overset{\∞}{\∫}}s\left(t\right){\mathrm{\ψ}}_k\left(t\right)e^{-2\mathrm{i\πft}}\mathrm{dt}---\left(6\right)$

将非平稳地震记录s(t)按照(6)式变换到分子-Gabor域。

所述步骤3)中，在分子-Gabor域对每个地震记录片段进行频带拓宽和振幅校正。本发明针对不同情况给出两种方法，可根据实际需要选其中一种。

方法1：保持原始地震记录相对能量关系的拓频

以第k个分子窗截出的地震记录片段为例，为了提高第k个地震记录片段的分辨率，我们采用下式拓宽该片段的频带：

${\hat{s}}_k^h\left(f\right)={\hat{s}}_k\left(f\right)/\left(\right|{\hat{s}}_k^w\left(f\right)|+\mathrm{\ϵ}{A}_{\max })---\left(7\right)$

式中ε是一个小于1的无量纲的常数， 

是第k个片段的F-function(拟合函数)，它与第k个片段的等效子波的振幅谱相似，A_max为 

的最大值，εA_max用来增强计算的稳定性。由于F-function和等效子波的振幅谱之间存在一个比例系数，因此，(7)式的结果会产生类似于AGC(幅度均衡)的等幅效应。我们采用

${\hat{s}}_k^{\′}\left(f\right)=\left|\right|{\hat{s}}_k\left|\right|{\hat{s}}_k^h\left(f\right)/\left|\right|{\hat{s}}_k^h\left|\right|---\left(8\right)$

来校正 

的能量。 

即为第k个地震记录片段对应的保持原始地震记录相对能量关系的拓频结果。

方法2：带衰减补偿的拓频

这里，我们首先对每个地震记录片段进行衰减补偿，然后对补偿后的结果进行频带拓宽和振幅校正。以第j个分子窗截出的地震记录片段为例。把由参考子波(震源子波)到第j个分子之间的介质视为均匀粘弹性介质，介质的等效品质因子记为Q_j。令参考子波从震源传播到中心为j的分子窗处所用的时间为T_j，则平面波在频率域满足因果律的传播算子可表示为(Wang，Y.，“Seismic Inverse Q Filtering”.Singapore：Blackwell Publishing Ltd.，2008)：

$h(f,T_j)=e^{-\mathrm{\πf}{T}_j/Q_j}{e}^{\mathrm{iH}[\mathrm{\πf}{T}_j/Q_j]}---\left(9\right)$

此处H为在某个时刻t对频率f的Hilbert变换， 

如前所述，用 

表示第j个片段上的等效子波的Fourier(傅立叶)频谱，根据波传播理论有：

${\hat{w}}_j\left(f\right)=h(f,T_j)\hat{w}\left(f\right)---\left(10\right)$

这里 

为参考子波的Fourier频谱。把(9)式代入(10)式，且仅考虑振幅部分得到：

$\left|{\hat{w}}_j\left(f\right)\right|=e^{-\mathrm{\πf}{t}_j/Q_j}\left|\hat{w}\left(f\right)\right|---\left(11\right)$

由于本发明构造的分子窗所覆盖的信号段近似平稳，所以每个窗截取的信号段有一个近似不变的等效子波，该等效子波可用该段的F-function乘以一个常数来表示。记 

的F-function为 

则有：

$\left|{\hat{s}}_j^w\left(f\right)\right|\≈C_j\left|{\hat{w}}_j\left(f\right)\right|---\left(12\right)$

这里C_j为第j个片段的常数比例因子，由(11)和(12)式可得：

$\left|{\hat{s}}_j^w\left(f\right)\right|\≈C_j{e}^{-\mathrm{\πf}{T}_j/Q_j}\left|\hat{w}\left(f\right)\right|---\left(13\right)$

对(13)式两边取对数，采用对数谱比值法可计算出Q_j，从而可得到|h(f，T_j)|。令

β_j(f)＝1/h(f，T_j)|(14)

用β_j乘以 

我们可得到补偿后的片段：

${\hat{s}}_j^c\left(f\right)={\hat{s}}_j\left(f\right){\mathrm{\β}}_j\left(f\right)---\left(15\right)$

令 

表示 

的F-function，将(7)式和(8)式作用于 

可得拓频结果

${\hat{s}}_j^{c,h}\left(f\right)={\hat{s}}_j^c\left(f\right)/\left(\right|{\hat{s}}_j^{c,w}\left(f\right)|+\mathrm{\ϵ}{A}_{\max })---\left(16\right)$

和能量校正后结果

${\hat{s}}_j^{\′\′}\left(f\right)=\left|\right|{\hat{s}}_j^c\left|\right|{\hat{s}}_j^{c,h}\left(f\right)/\left|\right|{\hat{s}}_j^{c,h}\left|\right|---\left(17\right)$

即为第j个地震记录片段对应的带衰减补偿的拓频结果。

所述步骤4)中，分子-Gabor逆变换定义如下：

$s\left(t\right)=\underset{k=1}{\overset{L}{\mathrm{\Σ}}}\left[E_k\underset{-\∞}{\overset{\∞}{\∫}}{\hat{s}}_k\left(f\right)e^{2\mathrm{i\πft}}\mathrm{df}\right]---\left(18\right)$

对(8)式处理的结果做分子-Gabor逆变换，可得到保持原始地震记录相对能量关系的高分辨地震道。对(17)式处理的结果做分子-Gabor逆变换，可得到带衰减补偿的高分辨地震道。

图2给出了合成地震记录算例图，其中图2(a)是反射系数序列；(b)是非平稳合成地震记录；(c)是(b)中信号对应的能量归一化后的分子-Gabor窗；(d)是谱白化方法处理的结果；(e)和(f)分别为(8)式和(17)式对应的自适应拓频处理后的高分辨地震记录。对比图2(a)和(b)可见，反射系数①，②，③和④在(b)中不清晰。①和②在(d)中比在(b)中清晰，然而③和④在(d)中仍然不清晰。在图2(e)和(f)中，反射系数①，②，③和④都得到很好的刻画。并且图2(e)中的高分辨地震记录保持了图2(b)所示原始地震记录的相对能量关系；图2(f)中的高分辨地震记录恢复了衰减的能量，与图2(a)中的反射系数序列有较好的对应关系。

图3(a)、3(b)是实际地震记录片段图及其对应的高分辨处理后的地震记录片段图，由测井资料得到的合成地震记录已经插入到图(b)中井标记 

的位置。对比图3(a)和(b)可见，高分辨处理后的地震记录纵向分辨率明显提高，对地层的刻画更加精细， 且和测井资料得到的合成地震记录有很好地对应关系，如图中箭头标识的位置所示。

图4(a)和4(b)分别给出了图3中实际地震资料高分辨处理前后振幅谱图。由振幅谱对比可见，高分辨处理后，地震资料主频从原来的20Hz提高到40Hz，3dB带宽从原来的18Hz提高到40Hz，地震资料频带得到有效拓宽。

Claims (6)

1.一种基于地震记录变子波模型提高地震资料分辨率的方法，其特征在于，包括以下步骤：

1)构造分子-Gabor窗，将非平稳地震记录自适应地划分为若干平稳的片段，每个片段中拥有一个近似不变的等效子波，且所述等效子波易于从该片段中提取出来；

2)用第1)步构造得到的分子窗生成分子-Gabor标架，将非平稳地震记录变换到分子-Gabor域；

3)在分子-Gabor域，对每个分子-Gabor窗内的地震记录片段所对应的分子-Gabor系数，进行拓频和能量补偿处理；

4)将处理后的分子-Gabor系数反变换到时间域，得到提高分辨率后的地震记录。

2.如权利要求1所述的一种基于地震记录变子波模型提高地震资料分辨率的方法，其特征在于，所述步骤1)中，分子-Gabor窗通过单位分解法构造，具体分为三个步骤：

①生成满足单位分解的原子窗族

选择基本原子窗函数G(t)为Lamoureux函数，用G_j(t)＝G(t-jΔt)表示中心位于第j个采样点上的原子窗，对原子窗族{G_j(t)：1≤j≤N}按下式归一化：

$g_j\left(t\right)=G_j\left(t\right)/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{G}_i\left(t\right)---\left(1\right)$

N为地震道采样点的个数，由(1)式可得一组满足单位分解的原子窗族{g_j(t)：1≤j≤N}；

②构造初始分子-Gabor窗

首先定义地震信号s(t)的瞬时频率如下：

f(t)＝1/2π{[s(t)ds^*(t)/dt-s^*(t)ds(t)/dt]/[a(t)²+ξ²]}(2)

这里s^*(t)为s(t)的Hilbert变换，a(t)＝[s(t)²+s^*(t)²]^1/2为s(t)的包络，ξ为无量刚的小常数；

其次用加权瞬时频率f′(t)代替局部均值频率来表征地层的吸收效应，加权瞬时频率定义如下：

$f^{\′}\left(t\right)={\∫}_{t-T}^{t+T}f\left(t\right)W\left(t\right)\mathrm{dt}/{\∫}_{t-T}^{t+T}W\left(t\right)\mathrm{dt}---\left(3\right)$

这里W(t)取为s(t)包络幅值的平方；

然后，采用保边缘平滑拟合的方法消除子波干涉对加权瞬时频率的影响：选取信号包络峰值处的加权瞬时频率，对它们作保边缘平滑拟合，检验出偏离拟合曲线较远的数据点，将这些较远的数据点去掉，然后用线性插值去填充这些去掉的点，如此反复迭代直到拟合曲线上的值不再有大的变化，从而得到最终拟合的曲线；

最后根据拟合曲线上相邻两点之间的差值和距离计算各分子窗的分割点，然后将相邻分割点间的小原子窗叠加起来就得到了自适应分子窗，设第k个分子窗的起点和终点分别为M^k-1和M^k，则第k个分子窗可以表示为：

${\mathrm{\ψ}}_k\left(t\right)=\underset{j=M^{k-1}}{\overset{M^k-1}{\mathrm{\Σ}}}{g}_j\left(t\right)---\left(4\right)$

③对初始分子-Gabor窗进行能量归一化，得到分子-Gabor窗

令E_k表示第k个分子窗的能量，则有

$E_k=\left|\right|{\mathrm{\ψ}}_k\left|\right|={\left[{\∫}_{-\∞}^{+\∞}{\left|{\mathrm{\ψ}}_k\left(t\right)\right|}^2\mathrm{dt}\right]}^{1/2}---\left(5\right)$

能量归一化以后第k个分子窗ψ_k(·)/E_k即为我们所构造的分子-Gabor窗。

3.如权利要求1或2所述的一种基于地震记录变子波模型提高地震资料分辨率的方法，其特征在于，所述步骤2)中，由步骤1)中得到的分子-Gabor窗构造分子标架，相应的分子-Gabor变换定义如下：

${\hat{s}}_k\left(f\right)=\frac{1}{E_k}\underset{-\∞}{\overset{\∞}{\∫}}s\left(t\right){\mathrm{\ψ}}_k\left(t\right)e^{-2\mathrm{i\πft}}\mathrm{dt}---\left(6\right)$

将非平稳地震记录s(t)按照(6)式变换到分子-Gabor域。

4.如权利要求1所述的一种基于地震记录变子波模型提高地震资料分辨率的方法，其特征在于，所述步骤3)中，在分子-Gabor域对每个地震记录片段进行频带拓宽和振幅校正，方法之一是保持原始地震记录相对能量关系的拓频：

对于第k个分子窗截出的地震记录片段，采用下式拓宽该片段的频带：

${\hat{s}}_k^h\left(f\right)={\hat{s}}_k\left(f\right)/\left(\right|{\hat{s}}_k^w\left(f\right)|+\mathrm{\ϵ}{A}_{\max })---\left(7\right)$

式中ε是一个小于1的无量纲的常数，

是第k个片段的F-function函数，它与第k个片段的等效子波的振幅谱相似，A_max为

的最大值，εA_max用来增强计算的稳定性，采用

${\hat{s}}_k^{\′}\left(f\right)=\left|\right|{\hat{s}}_k\left|\right|{\hat{s}}_k^h\left(f\right)/\left|\right|{\hat{s}}_k^h\left|\right|---\left(8\right)$

来校正

的能量，

即为第k个地震记录片段对应的保持原始地震记录相对能量关系的拓频结果。

5.如权利要求4所述的一种基于地震记录变子波模型提高地震资料分辨率的方法，其特征在于，所述步骤3)中，在分子-Gabor域对每个地震记录片段进行频带拓宽和振幅校正，方法之二是带衰减补偿的拓频：

对于第j个分子窗截出的地震记录片段，把由参考子波到第j个分子之间的介质视为均匀粘弹性介质，介质的等效品质因子记为Q_j，令参考子波从震源传播到中心为j的分子窗处所用的时间为T_j，则平面波在频率域满足因果律的传播算子表示为

$h(f,T_j)=e^{-\mathrm{\πf}{T}_j/Q_j}{e}^{\mathrm{iH}[\mathrm{\πf}{T}_j/Q_j]}---\left(9\right)$

H为在某个时刻t对频率f的Hilbert变换，

用

表示第j个片段上的等效子波的Fourier频谱，根据波传播理论有：

${\hat{w}}_j\left(f\right)=h(f,T_j)\hat{w}\left(f\right)---\left(10\right)$

为参考子波的Fourier频谱，把(9)式代入(10)式，且仅考虑振幅部分得到：

$\left|{\hat{w}}_j\left(f\right)\right|=e^{-\mathrm{\πf}{t}_j/Q_j}\left|\hat{w}\left(f\right)\right|---\left(11\right)$

由于构造的分子窗所覆盖的信号段近似平稳，所以每个窗截取的信号段有一个近似不变的等效子波，该等效子波可用该段的F-function乘以一个常数来表示，记的F-function为

则有：

$\left|{\hat{s}}_j^w\left(f\right)\right|\≈C_j\left|{\hat{w}}_j\left(f\right)\right|---\left(12\right)$

这里C_j为第j个片段的常数比例因子，由(11)和(12)式可得：

$\left|{\hat{s}}_j^w\left(f\right)\right|\≈C_j{e}^{-\mathrm{\πf}{T}_j/Q_j}\left|\hat{w}\left(f\right)\right|---\left(13\right)$

对(13)式两边取对数，采用对数谱比值法可计算出Q_j，从而可得到|h(f，T_j)|，

令

β_j(f)＝q/|h(f，T_j)|(14)

用β_j乘以

得到补偿后的片段：

${\hat{s}}_j^c\left(f\right)={\hat{s}}_j\left(f\right){\mathrm{\β}}_j\left(f\right)---\left(15\right)$

令

表示的F-function，将(7)式和(8)式作用于

可得拓频结果

${\hat{s}}_j^{c,h}\left(f\right)={\hat{s}}_j^c\left(f\right)/\left(\right|{\hat{s}}_j^{c,w}\left(f\right)|+\mathrm{\ϵ}{A}_{\max })---\left(16\right)$

和能量校正后的结果

${\hat{s}}_j^{\′\′}\left(f\right)=\left|\right|{\hat{s}}_j^c\left|\right|{\hat{s}}_j^{c,h}\left(f\right)/\left|\right|{\hat{s}}_j^{c,h}\left|\right|---\left(17\right)$

即为第j个地震记录片段对应的带衰减补偿的拓频结果。

6.如权利要求1或5所述的一种基于地震记录变子波模型提高地震资料分辨率的方法，其特征在于，所述步骤4)中，分子-Gabor逆变换定义如下：

$s\left(t\right)=\underset{k=1}{\overset{L}{\mathrm{\Σ}}}\left[E_k\underset{-\∞}{\overset{\∞}{\∫}}{\hat{s}}_k\left(f\right)e^{2\mathrm{i\πft}}\mathrm{df}\right],---\left(18\right)$

对(8)式处理的结果做分子-Gabor逆变换，得到保持原始地震记录相对能量关系的高分辨地震道。对(17)式处理的结果做分子-Gabor逆变换，得到带衰减补偿的高分辨地震道。

Images