一种保地层结构的地震波衰减补偿方法与流程

本发明属于地球物理勘探领域，涉及一种地震波衰减补偿方法，特别涉及一种保地层结构的地震波衰减补偿方法。

背景技术：

实际地下介质是粘弹性的，当地震子波在地下介质中传播时，波前扩散、地层吸收衰减及散射等效应会使子波能量衰减、频带变窄，即引起地震子波时变。地层吸收衰减效应是引起地震子波时变的主要因素。在地震资料处理及解释中，常常需要对地震记录进行吸收补偿。例如，为了提高地震数据深层弱反射信号的分辨率，常常需要对地层吸收进行补偿；常用的提高地震记录纵向分辨率的方法，大都假定地震记录是平稳的(即地震子波在地下传播过程中是不变的)，为了使用这些方法，首先必须对地震记录进行吸收补偿。再如，我们知道，avo分析的前提假设之一是，介质是完全弹性的，而实际介质是粘弹性的，为了进行avo分析，也需要对地震记录进行吸收补偿。关于地层吸收补偿方法，一直是国内外学者研究的热点之一。

在已知地下介质中q值分布的情况下，可通过反q滤波完成吸收补偿。地震子波随着传播时间变化的地震记录称之为非平稳地震记录。时频分析是研究非平稳信号的有力工具，因此，在时频域估计介质的吸收因子(或q值)对反射地震记录进行地层吸收补偿是一条重要途径。

1999年，白桦和李鲲鹏提出了基于短时fourier变换(stft)的地层吸收补偿方法，该方法不需要预先知道地层的q值，对于q值随深度变化的情况也适用。为了提高吸收补偿的精度，李鲲鹏(2000)和刘喜武(2006)，分别探讨了基于小波包分解(wt)和基于广义s变换(gst)的地层吸收补偿方法。上述时频域地层吸收补偿方法的假设条件是：“地层足够厚，相邻界面的反射波在地震记录上互不重叠”(以下简称厚层)。在此假设下，这类方法的原理基于这样的事实：若没有地层吸收，反射波各频率分量的能量对时间的分布具有相似性，地层吸收破坏了这一相似性，若能在时频平面上恢复这一相似性，就可补偿地层吸收效应。然而，在我国陆相含油气盆地中，地层较薄(相对于地震波长)，相邻波阻抗界面的反射波在地震记录上常常不能完全分开，这使得这类方法的假设在许多情况下不能成立。

margrave等人以gabor变换为工具，在假定地层为均匀粘弹性介质，且仅考虑地层吸收衰减的情况下，研究了非平稳地震记录的模型。他们指出：地震记录的gabor谱(即地震记录的gabor变换)可近似的等于震源子波的频谱、地层q滤波器和反射系数的gabor谱三者的乘积。据此，他们讨论了均匀粘弹性介质中的q值估计问题，并用于提高地震资料的分辨率。然而，实际地层常常不能视为均匀粘弹性介质，不同深度的地层q值不同，使得该方法在实际应用中受到限制。2010年，汪玲玲等人根据地震信号振幅谱的时变特性，选择地震信号包络峰值处的加权瞬时频率随时间的变化量及它们相互之间的距离的大小为参考量，构造自适应分子标架，然后基于分子分解，利用相关系数法求取补偿滤波器，实现对地震记录的吸收补偿。该方法不需要预先知道q值，对于地层q值随深度变化的情况也适用，然而，由于噪声干扰的影响可能造成分窗点在横向上不连续，从而影响补偿后资料的横向连续性。因此，对于地层结构变化复杂的大规模实际三维叠后地震资料，有必要进一步调整该方法。2011年，margrave等人用层间q计算得到的等效q重新表述了他们在2001年给出的非平稳地震道模型，从而将他们的方法推广到地层q值随深度变化的情况，然而，这一方法使用固定长度的分子分析时窗，窗长需要人为设置，不能适应于地震记录各段不同的特性，若窗长选择不合适，在处理具有复杂结构的地震资料时会存在问题。

技术实现要素：

本发明的目的在于克服上述现有技术的缺点，提供一种基于自适应分子分解的可保地层结构的地震波衰减补偿方法，该方法采用地震包络局部峰值作为地层结构约束来构造自适应于地层结构的分子时窗，然后采用非线性压缩映射子波振幅谱估计方法从地震道分子分解时频谱中估计正比于时变子波振幅谱的分量，接着计算质心频率，并对质心频率做筛选，减少子波干涉和噪声的影响，然后用质心频率偏移法估计得到稳健的地下介质q场，并计算得到各时窗对应的吸收补偿滤波器，最后用这些滤波器对相应的时频谱进行加权补偿，再反变换回时间域，就可得到衰减补偿后的地震记录。

本发明的目的是通过以下技术方案来解决的：

一种保地层结构的地震波衰减补偿方法，采用地震包络局部峰值作为地层结构约束来构造自适应于地层结构的分子时窗，然后采用非线性压缩映射子波振幅谱估计方法从地震道分子分解时频谱中估计正比于时变子波振幅谱的分量，接着计算质心频率，并对质心频率做筛选，减少子波干涉和噪声的影响，然后用质心频率偏移法估计得到稳健的地下介质q场，并计算得到各时窗对应的吸收补偿滤波器，最后用这些滤波器对相应的时频谱进行加权补偿，再反变换回时间域，得到衰减补偿后的地震记录。

本发明进一步的改进在于，具体包括以下步骤：

1)提取地震道包络局部峰值：

2)生成满足单位分割的原子窗集：

选择基本原子窗函数g(t)，令gj(t)＝g(t-jδt)表示中心位于第j个采样点的原子窗，对原子窗族{gj:1≤j≤n}按下式归一化

得一组单位分割原子窗集{gj:1≤j≤n}，这里n为地震道采样点的个数；

3)构造自适应分子时窗：

选择每相邻两个包络局部峰值点的中点作为各分子窗的边界点，将边界点间的满足单位分割的原子窗叠加起来形成了分子窗；设第k个分子窗对应的包络局部峰值点位于pk，前一个包络局部峰值点位于pk-1(p0＝-p1)，后一个包络局部峰值点位于pk+1，则此分子窗ψk(t)，对应的第一个原子窗的中心位于mk-1+1＝(pk-1+pk)/2+1(m0＝0)，对应的最后一个原子窗的中心位于mk＝(pk+pk+1)/2，分子窗由这中间的mk-mk-1个原子窗叠加得到，即分子窗ψk(t)由下式表示

令l为分子窗的个数，分子窗族{ψk(t):1≤k≤l}也构成单位分割；

4)对步骤3)得到的自适应分子时窗进行能量归一化：

令ek表示第k个分子窗的能量，即

能量归一化以后的分子分析时窗为{ψk(t)/ek:1≤k≤l}；

对分子分析时窗进行平移和调制后，得到一组分子标架；此时，分子分解时频变换定义为

其中f为频率，分子分解逆变换定义为

5)利用非线性压缩映射提取时变子波振幅谱：

设x0∈(a,b)，取δ0＞0，使得对于任意给定的δ＞0，定义l²[a,b]空间中的子集

在[a,x0-δ0]上单增,且在[x0+δ0,b]上单减}如果且q≤2，则有任取定义函数

此处0＜q≤1；显然，0≤fq(u；x)≤1；设cq＞0，α,β＞1，对于定义上的算子p为

p是到的非线性压缩映射算子；

对于第k个分子窗中的地震记录片段，其振幅谱为设定迭代初值为：

其中0＜q≤1，

其中fm为振幅谱的截止频率；建立迭代uk＝p(uk-1)，其中算子参数由最小二乘得到；由于算子是压缩映射，得到不动点u^*，则估计的子波振幅谱为记为lk(f)；

6)计算质心频率：

对于第k个分子窗中的地震记录片段，质心频率fc,k为

式中fc为子波振幅谱的截止频率；对计算得到的质心频率处理后得到最终的质心频率

7)估算地下介质q场：

用质心频率偏移法估算q值，相应的估算公式如下

式中：分别为t1时刻频谱的质心频率和方差；是t2时刻频谱的质心频率；时差δt＝t2-t1；为t1和t2之间地层介质的品质因子；进一步推导，得到t时刻地下介质的等效q值

式中fc,0和为初始时刻地震子波振幅谱的质心频率和方差，即

用步骤6)中计算得到的质心频率代替fc(t)，得到各个窗中心点处的q值

tk表示第k个分子窗的中心点所对应的时间位置；对窗中心点处估得的q值插值后即得到稳定的地下介质q场；

8)计算各个时窗对应的补偿滤波器：

将由参考子波到第k个分子gabor时窗之间的介质视为均匀粘弹性介质，介质的等效品质因子即为步骤7)计算得到的qe,k；参考子波从震源传播到第k个分子gabor时窗的中心所用的时间即为tk，则在第k个时窗内平面波在频率域满足因果律的传播算子hk(f)表示为

hk(f)＝exp{-πftk/qe,k+ih[πftk/qe,k]}

此处h(·)表示在某个时刻t对频率f的hilbert变换，相应的补偿滤波器为

βk(f)＝1/hk(f)＝1/exp{-πftk/qe,k+ih[πftk/qe,k]}

9)修正补偿滤波器：

记参考子波的截止频率为fc1，对于第k个片段，设补偿滤波器的值等于某个常数时对应的频率为fc,k，取fc＝min{fc1,fc,k}，设计高斯边缘低通滤波器如下

其中，σ为高斯函数的标准差；从而将补偿滤波器修订为

10)用修正后的补偿滤波器按下式校正时频系数

11)将步骤10)得到的时频系数反变换回时间域，得到衰减补偿后的地震记录；即

其中sc(t)即为衰减补偿后的地震记录。

本发明进一步的改进在于，步骤1)中提取地震道包络局部峰值具体过程为：设s^*(t)为地震道s(t)的hilbert变换，则

a(t)＝[s(t)²+s^*(t)²]^1/2

为地震道s(t)的包络；由上式计算地震道的包络，并提取包络局部峰值点。

本发明进一步的改进在于，步骤6)中对计算得到的质心频率做多次光滑滤波，并计算滤波前后的绝对误差，去掉绝对误差较大的一些质心频率后，再次拟合滤波，得到最终的质心频率

本发明进一步的改进在于，步骤9)中某个常数为50-100之间的数。

与现有技术相比，本发明具有以下有益效果：

本发明利用地震包络峰值作为地层结构约束来构造自适应分子时窗，从而使得构造的时窗在横向上与地层结构有关，有利于保持衰减补偿后地震资料的横向连续性；采用非线性压缩映射子波振幅谱估计方法可以从较短数据中估计子波振幅谱，这为本发明从地震道分子分解时频谱中估计正比于时变子波振幅谱的分量，从而得到稳定的质心频率提供了保证；在计算衰减因子之前对质心频率做了滤波筛选，可以有效减少子波干涉和随机噪声的影响，进一步保障了处理后资料的横向连续性。本发明不仅能够有效地补偿地层吸收衰减，提高地震资料纵向分辨率和深层弱信号识别能力，而且能够很好地保持地震资料的横向连续性和复杂地层结构特征。

附图说明

图1是合成地震道时反射系数示意图；

图2是合成地震道的包络局部峰值示意图；其中，实线表示非平稳合成地震道，虚线表示地震道包络，星号表示包络局部峰值；

图3是构成单位分割的原子gabor时窗示例图；

图4是地震道包络局部峰值约束下寻找分子gabor时窗边界点示例图；

图5是将边界点间的原子gabor时窗叠加形成分子gabor时窗示例图；

图6是能量归一化后得到的分子gabor时窗示例图；

图7为原始地震资料剖面图；

图8为从原始地震资料剖面中提取的质心频率图；

图9为对质心频率做滤波筛选后得到的质心频率图；

图10为基于滤波筛选后的质心频率计算得到的地下介质q场图；

图11为地层衰减补偿后的地震资料剖面图；

图12为图7中矩形框中区域放大图；

图13为图11中矩形框中区域放大图。

具体实施方式

本发明针对复杂的地质结构和噪声干扰常常使得从地震数据中估计出的地下介质的q值不稳定，导致衰减补偿后的地震数据不能保持横向连续性及复杂地层结构特征的问题，提出了一种基于自适应分子分解的可保地层结构的地震波衰减补偿方法。该方法用地震道包络局部峰值作为地层结构约束来构造自适应分子时窗，并用这些时窗生成分子标架对地震记录做分子分解；接着在自适应分子分解时频域计算各时窗对应的吸收补偿滤波器；最后用这些滤波器对相应的时频谱进行加权补偿，再反变换回时间域，从而得到地震波衰减补偿后的地震记录。

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

本发明的具体步骤如下：

1)提取地震道包络局部峰值；

根据复道分析原理，设s^*(t)为地震道s(t)的hilbert变换，则

a(t)＝[s(t)²+s^*(t)²]^1/2

为s(t)的包络。对比图1中合成地震道所用的反射系数和图2中合成地震道包络局部峰值点可见，地震道的包络局部峰值与反射界面有一定的对应关系，在一定程度上可以大致反映地层的层序结构。如果使用地震道的包络局部峰值约束分子窗的构造，那么构成的分子窗在横向上将与地层结构有关，有利于保持处理后资料的横向连续性。

2)生成满足单位分割的原子窗集；

恰当地选择基本原子窗函数g(t)，令gj(t)＝g(t-jδt)表示中心位于第j个采样点的原子窗，对原子窗族{gj:1≤j≤n}按下式归一化

可得一组单位分割原子窗集{gj:1≤j≤n}，这里n为地震道采样点的个数。如图3中的一组窗就是满足单位分割的原子窗集。

3)构造自适应分子时窗；

选择每相邻两个包络局部峰值点的中点作为各分子窗的边界点，将边界点间的满足单位分割的原子窗叠加起来就形成了分子窗。如图4所示，设第k个分子窗对应的包络局部峰值点位于pk，前一个包络局部峰值点位于pk-1(p0＝-p1)，后一个包络局部峰值点位于pk+1，则此分子窗，ψk(t)，对应的第一个原子窗的中心位于mk-1+1＝(pk-1+pk)/2+1(m0＝0)，对应的最后一个原子窗的中心位于mk＝(pk+pk+1)/2，分子窗由这中间的mk-mk-1个原子窗叠加得到，即ψk(t)可由下式表示

令l为分子窗的个数，易知分子窗族{ψk(t):1≤k≤l}也构成单位分割，如图5所示。这样，可以保证每个分子窗内至少有一个反射子波，在减少时窗数量的同时，可以有效减少窗端点对子波的截断效应。并且由于包络峰值可以大致反映地层的层序结构，构成的分子窗在横向上与地层结构有关，有利于保持处理后资料的横向连续性。

4)对步骤3)得到的自适应分子时窗进行能量归一化；

分子分解时频变换是保能量的变换，由图5可以看出，在步骤三得到的分子窗族中，不同的窗能量不同。如果用这些窗直接构造分子gabor标架，则一个地震道的时频能量表示不仅与该地震道有关，还与分子窗有关。希望这一时频能量表示只与地震道有关，为此，需要对每个分子窗进行能量归一化。令ek表示第k个分子窗的能量，即

能量归一化以后的分子分析时窗为{ψk(t)/ek:1≤k≤l}，如图6所示。

对分子分析时窗进行平移和调制后，可得到一组分子标架。此时，分子分解时频变换可定义为

其中f为频率。分子分解逆变换可定义为

5)利用非线性压缩映射提取时变子波振幅谱；

设x0∈(a,b)，取δ0＞0足够小，使得对于任意给定的δ＞0，定义l²[a,b]空间中的子集

在[a,x0-δ0]上单增,且在[x0+δ0,b]上单减}如果且q≤2，则有任取定义函数

此处0＜q≤1。显然，0≤fq(u；x)≤1。设cq＞0，α,β＞1，对于定义上的算子p为

显然，p[u；x]在区间[a,b]上是连续的，并且δ≤p[u；x]≤cq+δ，由此可得p[u；x]∈l²[a,b]。因此，p是到l²[a,b]的非线性算子。根据banach空间算子不动点定理，可以证明：p是到的压缩算子，并且在中存在唯一不动点。

对于第k个分子窗中的地震记录片段，其振幅谱为设定迭代初值为：

其中0＜q≤1，

其中fm为振幅谱的截止频率。建立迭代uk＝p(uk-1)，其中算子参数由最小二乘得到。由于算子是压缩映射，可以得到不动点u^*，则估计的子波振幅谱为记为lk(f)。

6)计算质心频率；

对于第k个分子窗中的地震记录片段，质心频率为

式中fc为子波振幅谱的截止频率。理论上，地震子波在传播的过程中质心频率是随着传播时间衰减的，然而实际中(如图7所示)，当地层较薄时，地震波之间的互相干涉会使得质心频率突然增大或者突然减小(如图8所示)。为此，本发明在计算q值之前对计算得到的质心频率做多次光滑滤波，并计算滤波前后的绝对误差，去掉绝对误差较大的一些质心频率后，再次拟合滤波，得到最终的质心频率以减少地震波干涉的影响(如图9所示)。

7)估算地下介质q场；

用质心频率偏移法估算q值，相应的估算公式如下

式中：分别为t1时刻频谱的质心频率和方差；是t2时刻频谱的质心频率；时差δt＝t2-t1；为t1和t2之间地层介质的品质因子。进一步推导，可得t时刻地下介质的等效q值

式中fc,0和为初始时刻地震子波振幅谱的质心频率和方差，即

用步骤6)中计算得到的质心频率代替fc(t)，可得到各个窗中心点处的等效q值

tk表示第k个分子窗的中心点所在的时间位置。对窗中心点处的等效q值插值后，就可以得到稳定的地下介质q场(如图10所示)。

8)计算各个时窗对应的补偿滤波器；

把由参考子波(震源子波)到第k个分子gabor时窗之间的介质视为均匀粘弹性介质，介质的等效品质因子即为步骤7)计算得到的qe,k；参考子波从震源传播到第k个分子gabor时窗的中心所用的时间即为tk，则在第k个时窗内平面波在频率域满足因果律的传播算子hk(f)可表示为

hk(f)＝exp{-πftk/qe,k+ih[πftk/qe,k]}

此处h(·)表示在某个时刻t对频率f的hilbert变换，相应的补偿滤波器为

βk(f)＝1/hk(f)＝1/exp{-πftk/qe,k+ih[πftk/qe,k]}

9)修正补偿滤波器；

为了防止过补偿，本发明对补偿滤波器作了修订，只对有效频带范围内的时频系数作补偿。记参考子波的截止频率为fc1，对于第k个片段，设补偿滤波器的值等于某个常数(视情况而定，通常取50-100之间的数)时对应的频率为fc,k，取fc＝min{fc1,fc,k}，设计高斯边缘低通滤波器如下

其中，σ为高斯函数的标准差；从而可将补偿滤波器修订为

10)用修正后的补偿滤波器按下式校正时频系数

11)将步骤10)得到的时频系数反变换回时间域，即

其中sc(t)即为衰减补偿后的地震记录。

对比衰减补偿前后的地震资料剖面图7和图11，以及图7和图11中矩形框区域放大图即图12和图13可见，补偿后的地震资料纵向分辨率显著提高，深层弱信号的能量明显增强，对地层结构的刻画更加精细，横向连续性保持的很好，并且补偿后结果能够很好地保持原剖面中的复杂地质结构特征。

本发明用地震道包络局部峰值作为地层结构约束来构造自适应分子时窗，并用这些时窗生成分子标架对地震记录做分子分解；接着在自适应分子分解时频域计算各时窗对应的吸收补偿滤波器；最后用这些滤波器对相应的时频谱进行加权补偿，再反变换回时间域，就得到地震波衰减补偿后的地震记录。所述分子时窗的构造方法，可以保证每个时窗内至少有一个反射子波，在减少时窗数量、提高计算效率的同时，使得这些时窗在横向上与地层结构有关，有利于保持衰减补偿后地震资料的横向连续性。本发明能够在有效补偿地震波衰减、提高地震资料纵向分辨率和深层弱信号识别能力的同时，很好地保持地震资料的横向连续性和复杂结构特征。