基于对数时频域双曲平滑的动态反褶积方法

摘要

本发明涉及一种基于对数时频域双曲平滑的动态反褶积方法，时频谱平滑方法的选择直接决定了吸收衰减和震源子波的质量，对动态反褶积结果具有至关重要的影响。常规时频域的双曲平滑方法主要借助除法实现吸收衰减和震源函数的分离，需要在分母中增加白噪因子以确保数值稳定性。这将明显降低迭代效率，致使吸收衰减与震源子波估求不准确，必将对最终结果造成不利影响。鉴于此，本发明将双曲平滑方法推广到对数时频域，运用线性拟合(舍去截距)与减法实现了吸收衰减和震源函数的完美分离，避免了白噪因子的使用，提高了迭代效率。另外，本发明提出了一种分段线性化方法，克服了实际处理中吸收衰减曲线并非呈单一线性的瓶颈。模型测试和实际资料处理结果从诸多方面都证明了本发明相对于传统方法的优越性。

权利要求

1.基于对数时频域双曲平滑的动态反褶积方法，其特征在于采用以下步骤实现：步骤1 ：根据POU 条件，确定一套分析时窗与合成时窗；

步骤2 ：借助分析时窗与快速傅里叶变换实现Gabor 变换，获得输入地震道的时频谱，并对其值谱取对数获得对数频谱；

步骤3 ：根据地震资料的信噪比确定有效振幅水平，确定双曲平滑的时频范围并执行双曲划分；

步骤4 ：在对数时频域执行双曲平滑，优化对吸收衰减的估计，移除吸收衰减的影响，根据纵向平均、横向平滑求震源子波；

步骤5 ：采用最小相位假设构建传播子波谱，形成反褶积因子，在时频域执行反褶积；

步骤6 ：借助合成时窗与逆傅里叶变换实现反Gabor 变换，将反褶积结果转换到时间域，取其实部作为最终输出结果。

2.根据权利要求1 所述的基于对数时频域双曲平滑的动态反褶积方法，其特征在于，步骤4 采用对数时频域双曲平滑方法，解决了常规双曲平滑方法估求吸收衰减和震源子波不准确的缺点，并且减小了白噪因子对最终结果的影响，在对数时频域中，吸收衰减函数转化为关于tf 的线性表达式Lα(t，f) ＝ tf(-π/Q)，地震道的动态褶积模型转换为 ，即动态褶积模型的斜率与表征吸收衰减的Q 值密

切相关，截距则包含子波项与反射系数项，因此通过最佳线性拟合求斜率，同时舍去截距项，可以在估计吸收函数的过程避免子波与反射系数非白噪的影响，实现吸收函数的精确求解。

3.根据权利要求1 所述的基于对数时频域双曲平滑的动态反褶积方法，其特征在于，步骤4 采用迭代思路不断移除震源子波的影响，获得准确的吸收衰减信息，之后，从动态褶积模型中减去该吸收衰减信息，并沿时间方向取均值，频率方向利用矩形窗褶积平滑进而获得震源子波的频率响应，由此实现传播子波的估算。

4.根据权利要求1 所述的基于对数时频域双曲平滑的动态反褶积方法，其特征在于，步骤4 将原时频域中的除法运算转换成对数时频域的减法运算，避免了白噪因子的使用，减小了人为因素对最终结果的影响。

5.根据权利要求1 所述的基于对数时频域双曲平滑的动态反褶积方法，其特征在于，在处理实际地震数据时，简单的线性拟合难以满足要求，需要利用吸收衰减曲线上各点斜率，采用类似于测井曲线分层的思路执行分段线性拟合。

说明书

基于对数时频域双曲平滑的动态反褶积方法

技术领域

[0001] 本发明属于地球物理勘探领域，涉及在对数时频域实现双曲平滑，进而完成动态反褶积，提高地震资料分辨率。

背景技术

[0002] 反褶积是地震资料处理中的一个基础环节，不仅能有效压缩地震记录中的地震子波，并且可以压制鸣震和多次波，能够明显提高地震数据的垂向分辨率。经过半个多世纪的发展，反褶积已演化多种成熟的方法，如最小平方反褶积、预测反褶积、子波整形反褶积、同态反褶积等。然而，目前大多数反褶积方法均建立在常规Robinson褶积模型之上，认为地震子波不随时间变化，忽略了地层的吸收衰减等效应，限制了其适用范围和实际价值，属于静态反褶积的范畴。

[0003] 动态褶积模型考虑了地层的吸收衰减，由Margrave等人推导化简并引入地震勘探领域(Margrave，1998；Margrave等，2002)。在时频域，该模型表示为震源子波、吸收衰减与反射系数相乘的形式，若只考虑幅值关系，其形式为：

[0004]

[0005] 其中 代表震源子波，|α(τ，f)|为吸收衰减， 代表反射系数项。动态反褶积的目的在于从地震道时频谱中移除震源子波与吸收衰减的影响，仅保留反射系数项，从而达到补偿吸收衰减，提高垂向分辨率的目的。

[0006] 毫无疑问，在动态反褶积的整个技术流程中震源子波与吸收衰减的估计是关键，而这两项的估计是通过平滑地震道时频谱获得的，因此平滑方法的选择对最终输出结果具有至关重要的影响。经过近几年的发展(Grossman等，2002；Montana等，2005)，已成熟的主流平滑方法有十字平滑，矩形窗平滑与常规时频域的双曲平滑。在研究中，许多研究对比工作发现十字平滑与矩形窗平滑的AGC效应严重，破坏了地震道垂向的振幅特征，输出结果会严重偏离实际情况(Margrave等，2003)。相比之下，常规双曲平滑利用恒Q与白噪反射系数的假设，在时频域沿tf＝常数的双曲线进行平滑，首先利用迭代思想从地震道时频谱中移除震源子波的影响，不断提纯对吸收衰减的估计。迭代结束后，从地震道时频谱移除最终估计的吸收衰减，并沿时间方向计算均值，沿频率方向平滑进而获得震源子波的频率响应。常规双曲平滑的输出结果稳定，是一种适用于动态反褶积的较合理的平滑方法(Margrave等，2011)。

[0007] 然而，常规双曲平滑方法的不足之处主要在于，利用迭代估计吸收衰减的过程中，需要借助除法实现不断移除震源子波的影响，然而为确保数值稳定性，防止除数很小的现象的出现，需要在分母中人为地增加白噪因子项，这将降低迭代效率，致使求解的吸收衰减与震源子波不准确，这将影响最终的反褶积效果。

发明内容

[0008] 针对常规双曲平滑方法的不足，本发明提供了一种在对数时频域实现双曲平滑的方法。通过本发明的改进，能够准确获得吸收衰减与震源子波，加速迭代收敛效率，减少由白噪因子引起的人为干扰，明显提高地震资料的分辨率，在油气勘探领域具有重要的应用价值。

[0009] 本发明实现上述目的的具体实施方案如下：

[0010] 步骤1：根据POU条件，确定一套分析时窗与合成时窗；

[0011] 步骤2：借助分析时窗与快速傅里叶变换实现Gabor变换，获得输入地震道的时频谱，并对其值谱取对数获得对数频谱；

[0012] 步骤3：根据地震资料的信噪比确定有效振幅水平，确定双曲平滑的时频范围并执行双曲划分；

[0013] 步骤4：用本发明所采用的算法在对数时频域执行双曲平滑，利用迭代思想，通过分段线性拟合并舍去截距的方法，不断优化对吸收衰减的估计。随后移除吸收衰减的影响，根据纵向(时间方向)平均、横向(频率方向)平滑求震源子波；

[0014] 步骤5：采用最小相位假设构建传播子波谱，形成反褶积因子，在时频域执行反褶积；

[0015] 步骤6：借助合成时窗与逆傅里叶变换实现反Gabor变换，将反褶积结果转换到时间域，取其实部作为最终输出结果。

[0016] 本发明提供的基于对数时频域双曲平滑的动态反褶积方法的基本原理如下：在对数时频域，动态褶积模型转换为(仅考虑幅值大小)：

[0017]

[0018] 式中第一项为吸收衰减项，Q代表岩石的品质因子；后两项为分别为子波项与反射系数项，两者共同构成了该线性方程的截距。

[0019] 该方程的斜率直接体现了地层的吸收衰减，子波和反射系数的影响则包含与截距之中。因此，在利用迭代过程求解吸收函数的过程中，对吸收衰减曲线(由各双曲条带内部均值构成)进行线性拟合，舍弃截距项就可以获得精确地吸收衰减。与此同时，在对数频率中减法代替了原时频域的除法，避免了白噪因子的使用，加速了迭代收敛速度，减少了最终结果中的人为干扰。

[0020] 然后利用白噪反射系数的假设，从动态褶积模型中移除经迭代获得的吸收衰减，并按照沿时间方向求均值，频率方向进行平滑得到震源子波的频率响应，进而以此形成震源子波时频谱。至此，吸收衰减时频谱与震源子波时频谱均已获得，使用Hilbert变换求最小相位，即可形成最终反褶积因子。

[0021] 在处理实际地震数据时，由各双曲条带内部均值构成的吸收衰减曲线通常不具有统一的线性规律，在大多数情况下呈现分段线性的特征，在这种情况下需根据该曲线的斜率，借助类似于测井曲线分层方法实现吸收衰减曲线的分段线性化。

[0022] 本发明的特色是：极大地简化了动态褶积模型，能够快速准确地计算吸收衰减与震源子波，避免了白噪因子的使用，减少了人为因素对最终结果的影响，提高了动态反褶积技术在增强地震数据分辨率方面的能力，在储层预测中极具重要价值。

附图说明

[0023] 图1(a)～图1(e)是由三种不同平滑方法所获得时间域的反褶积结果对比。 [0024] 图2(a)～图2(d)是由三种不同平滑方法所获得反褶积结果的时频域响应特征对比。

[0025] 图3是由两种双曲平滑所得震源子波的频率响应与实际响应对比。 [0026] 图4是基于对数时频域双曲平滑的动态反褶积方法的实际处理流程。 [0027] 图5(a)～图5(d)是实际处理中的分段线性化流程。

[0028] 图6(a)～图6(c)是基于两种不同双曲平滑的动态反褶积方法在处理某海上偏移地震剖面的效果对比。

具体实施方式

[0029] 以下结合附图说明本发明具体实施方式。

[0030] 图1、图2、图3来源于同一个反射系数模型。该模型分为三层，中间层的反射系数幅值水平是上、下两层的1/3。模型测试中，使用主频为40hz的最小相位的震源子波与该反射系数序列合成地震道，随后给定Q＝32在并时频域对该地震道进行衰减，由此获得最终的合成地震记录。

[0031] 图1(a)～图1(e)是由三种不同平滑方法所获得时间域的反褶积结果对比。图1(e)中由于矩形窗平滑无视地层的吸收衰减规律，反褶积结果存在严重的AGC效应，破坏了反射系数的纵向幅值关系；而图1(d)所示的常规双曲平滑与图1(c)所示的对数时频域双曲平滑方法遵循了地层的吸收衰减规律，能够很好地恢复原始幅值关系，可信度明显提高，相比之下，图1(c)所示的基于对数时频域双曲平滑的反褶积结果与实际反射系数序列的相似度要高于图1(d)所示的基于常规双曲平滑的反褶积结果。

[0032] 图2(a)～图2(d)是由三种不同平滑方法所获得反褶积结果的时频响应特征对比。图2(a)为实际反射系数序列的时频振幅谱。图2(b)所示的基于对数时频域双曲平滑的反褶积结果的时频响应与图2(a)的能量分布几乎一致，吻合程度最高；图2(c)所示的基于对数时频域双曲平滑的反褶积结果的时频响应次之，部分区域的能量没有得到准确恢复，如图中椭圆圈定部分所示；而图2(d)中所示的矩形窗平滑再次体现了严重的AGC效应，中间层的能量被明显放大。

[0033] 图3是由两种双曲平滑所得震源子波的频率响应与实际响应对比。对数时频域的双曲平滑方法使得在利用迭代过程中求吸收衰减时能够精确地去除震源子波的影响，并且限制了白噪因子的影响，这些优点大大提高了估算子波的精度。如图3中粗点划线所示，对数时频域的双曲平滑方法与粗实线所示的实际震源子波的频率响应具有极高的相似性。而如图3细实线所示，由常规双曲平滑估计的震源子波的频率响应在一定程度上偏离了理论响应。

[0034] 图1、图2、图3分别从不同角度证明了基于对数时频域双曲平滑的动态反褶积方法的优越性。

[0035] 图4是基于对数时频域双曲平滑的动态反褶积方法的实际处理流程： [0036] 步骤1：根据POU条件，确定一套分析时窗与合成时窗；

[0037] 步骤2：借助分析时窗与快速傅里叶变换实现Gabor变换，获得输入地震道的时频谱，并对其值谱取对数获得对数频谱；

[0038] 步骤3：根据地震资料的信噪比确定有效振幅水平，确定双曲平滑的时频范围并执行双曲划分；

[0039] 步骤4：用本发明所采用的算法在对数时频域执行双曲平滑，运用迭代思想，通过分段线性拟合并舍去截距的方法，不断优化对吸收衰减的估计。随后移除吸收衰减的影响，根据纵向(时间方向)平均、横向(频率方向)平滑求震源子波；

[0040] 步骤5：采用最小相位假设构建传播子波谱，形成反褶积因子，在时频域执行反褶积；

[0041] 步骤6：借助合成时窗与逆傅里叶变换实现反Gabor变换，将反褶积结果转换到时间域，取其实部作为最终输出结果。

[0042] 图5(a)～图5(d)是实际处理中的分段线性化流程。实际处理中，吸收衰减曲线通常不具有统一的线性规律，需要执行分段线性化。图5(a)为由各双曲条带内部均值构成的原始吸收衰减曲线；图5(b)为对图5(a)进行掐头去尾，平滑等预处理之后的吸收衰减曲线；图5(c)为图5(b)对应的斜率曲线，这是分段线性化的依据；类似于测井曲线分层方法，在图5(c)中利用峰-谷值之间长度值的某一百分比(本例中为50％)确定分段线性化的各个分界点的位置，如图5(d)所示。

[0043] 图6(a)～图6(c)是基于两种不同双曲平滑的动态反褶积方法在处理某海上偏移地震剖面的效果对比。图6(a)为原始地震剖面，目的层是具有杏仁状构造的生物礁。图6(b)是本发明处理(即对数时频域的双曲平滑)获得反褶积剖面；图6(c)是常规双曲平滑获得反褶积剖面。与传统方法相比，不难发现本发明能够更加清晰地反映礁体内部构造，并且生物礁下方的调谐效应在一定程度上也得到了压制。

[0044] 模型测试和实际处理均表明了本发明相对传统方法的巨大改进，并且在提高地震资料分辨率方面具有得天独厚的优势。

[0045] 以上具体实施方式仅用于说明本发明，而非用于限定本发明。

[0046] 参考文献

[0047] Grossman，J.P.，G.F.Margrave，and M.P.Lamoureux，2002，Constructing adaptive，nonuniform Gabor frames from partitions of unity：CREWES Research Report，14，1-10.

[0048] Margrave，G.F.，1998，Theory of nonstationary linear fltering in Fourier domain with application to time-variant filtering：Geophysics，63，244-259. [0049] Margrave，G.F.，D.C.Henley，M.P.Lamoureux et al，2003，Gabor deconvolution revisited：75th Annual International Meeting，SEG，22，714-717

[0050] Margrave，G.F.，M.P.Lamoureux，D.C.Henley，2011，Gabor deconvolution：Estimating reflectivity by nonstationary deconvolution of seismic data，Geophysics，76，W15-W30.

[0051] Margrave，G.F.，M.P.Lamoureux，J.P.Grossman，and V.Iliescu，2002，Gabor deconvolution of seismic data for source waveform and Q correction：72nd Annual International Meeting，SEG，Expanded Abstracts，2190-2193.

[0052] Montana，C.A.，G.F.Margrave，2005，Phase correction in Gabor deconvolution：75th Annual International Meeting，SEG，Expanded Abstracts，2173-2176。