CN106019376A - 一种频率驱动空变q值模型构建的地震波补偿方法 - Google Patents

Abstract

本发明公开了一种频率驱动空变Q值模型构建的地震波补偿方法，通过加伯变换生成对应于时间域地震数据每个样点的时间频率谱，逐点计算中心频率，由中心频率估算衰减值Q^‑1。品质因子Q值即衰减值Q^‑1的倒数，由此构建的Q值模型连续并且空变，计算过程稳定、高效。在此基础上实现高精度的地震波衰减补偿，解决了地震数据的振幅衰减和相位畸变问题，并提高地震数据的分辨率和信噪比。

Description

一种频率驱动空变Q值模型构建的地震波补偿方法

技术领域

本发明涉及地震信号处理领域，尤其涉及一种频率驱动空变Q值模型构建的地震波补偿方法。

背景技术

地震波由浅至深传播时，大地介质的滤波效应使地震波发生吸收衰减，引起地震波的振幅减弱、波形畸变、频带变窄以及主频降低，地下介质的这种岩石物理特性通常与品质因子Q值密切相关。研究估算品质因子Q值并对地震波补偿，可以消除地震波在传播过程中的衰减，恢复地震信号中地下岩层的反射系数特征和提高地震数据的分辨率。

地震波补偿，需要合理地估算Q值，计算Q值的常规方法有以下几种：

1)常Q值扫描法

常规地震波补偿处理中最常用到的就是常Q值扫描法，该方法简单易操作，但是其时变和空变性较差，主观性强，准确度低。

2)VSP求值法

VSP数据可以直接观测地震波的变化，准确度较高，但通常VSP数据缺乏，而且VSP法求得的Q值反应的仅是井点附近的地震波衰减的情况，难以通过少量的VSP数据计算全工区的Q值。

3)经验公式法

品质因子Q可以采用李庆忠院士的经验公式估算，该公式描述的是Q值与地震波传播速度V的函数，如果测井数据丰富，可以利用测得的速度V估算Q值，但缺乏测井资料的情况下，通过此方法估算Q值失效。

4)谱比法

谱比法使用较广泛，它假设Q值与频率无关，根据地震波振幅谱比值的对数随频率的变化关系，利用最小二乘法拟合斜率实现Q值的估算。谱比法原理简单，但要求地震波具有较高的信噪比，并且层顶、底界面两个时间处的地震子波不受附近其他子波的干扰，这限制了该方法计算Q值的适用性。

发明内容

本发明的发明目的是针对现有技术的缺陷，提供了一种频率驱动空变Q值模型构建的地震波补偿方法，用于高精度地提取空变Q值并实现地震波补偿。

根据本发明的一个方面，提供了一种频率驱动空变Q值模型构建的地震波补偿方法，包括如下步骤：

(1)通过加伯变换生成对应于时间域地震数据每个样点的时间频率谱；

(2)计算每个样点的中心频率，生成所有样点的中心频率数据体或剖面；

(3)根据所述每个样点的中心频率估算每个样点的衰减值Q^-1，品质因子Q值即衰减值Q^-1的倒数，由此构建Q值模型；

(4)在频率域实现地震波衰减补偿；

(5)利用加伯反变换将补偿后的频率域数据转换重构到时间域。

可选地，步骤(1)中加伯变换可逆，通过时间频率谱反变换能恢复原始时间域地震数据。

可选地，步骤(1)中使用空间上相邻的地震道组合，用于计算各个地震道的时间频率谱。

可选地，步骤(2)中对地震数据中相邻空间位置的时间频率谱进行组合，生成所有样点的中心频率数据体或剖面。

可选地。步骤(3)还包括由中心频率估算衰减值Q-1时，结合相邻地震道，使用中值滤波方法获取所述时窗起点处的中心频率和所述时窗起点的中心频率与终点的中心频率的比值。

根据本发明提供的频率驱动空变Q值模型构建的地震波补偿方法，通过加伯变换从叠后地震数据转换得到每个样点的频率信息进而逐点构建空变Q值模型，计算量小，提取Q值方法稳定，Q值模型纵向和横向上连续变化，符合实际地质情况。在高精度的Q值模型基础上，通过地震波补偿不仅尽可能地恢复了振幅信息和纠正了相位畸变，还提高了地震数据的分辨率和信噪比。

上述说明仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，而可依照说明书的内容予以实施，并且为了让本发明的上述和其它目的、特征和优点能够更明显易懂，以下特举本发明的具体实施方式。

附图说明

图1为本发明提供的频率驱动空变Q值模型构建的地震波补偿方法的流程图；

图2为实际三维地震数据体沿线方向抽取的地震剖面；

图3为三维中心频率数据体沿线方向抽取的剖面；

图4为三维Q数据体沿线方向抽取的剖面；

图5为三维Q数据体在近地表的水平切片；

图6为地震数据进行地震波衰减补偿前的地震剖面；

图7为地震数据进行地震波衰减补偿后的地震剖面；

图8为地震数据进行地震波衰减补偿前的频谱对比；

图9为地震数据进行地震波衰减补偿后的频谱对比。

具体实施方式

下面将参照附图更详细地描述本公开的示例性实施例。虽然附图中显示了本公开的示例性实施例，然而应当理解，可以以各种形式实现本公开而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更透彻地理解本公开，并且能够将本公开的范围完整的传达给本领域的技术人员。

从地震数据信息中提取Q值是一个重要的思路。叠后地震数据经过处理，无法保证每个步骤都是真振幅处理，基于地震数据的振幅或能量的变化来估算地层衰减时精度受限。相对而言，处理对频率成份的损害很小，比如，动校正之后的子波拉伸切除旨在保护地震记录的真实频率，时间域中的几何扩散补偿及自动增益均衡会改变其振幅大小，但是并不影响地震频率成份。在叠后地震数据中，频率成分比振幅变化更为可靠，因此，可以通过频率成分的变化进行Q分析。

下面为利用Q值实现地震波衰减补偿的基本原理：

根据波动方程理论，当地震波传播了Δx的距离之后，其平面波分量U(x,ω)可以表述为：

U(x+Δx,ω)＝U(x,ω)exp[-ik(ω)Δx] (1)

式中i为虚部单位，ω为角频率，k(ω)为平面波的波数。将波数变为复数，则可以表达地层的Q效应：

$k\left(\mathrm{\ω}\right)=\left(1-\frac{i}{2Q_r}\right)\frac{\mathrm{\ω}}{v_r}{\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_h}\right)}^{-\mathrm{\γ}}---\left(2\right)$

其中，为经验公式，Q_r和v_r分别为某一参考频率下的Q值和相速度，ω_h是与勘探地震频带中最高可能频率相对应的调谐参数。将复波数k(ω)代入(1)式并且将距离Δx增量替换为时间增量Δτ，即可得到地震波补偿(反Q滤波)的表达式如下：

$U(\mathrm{\τ}+\mathrm{\Δ}\mathrm{\τ},\mathrm{\ω})=U(\mathrm{\τ},\mathrm{\ω})\exp \[{\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_h}\right)}^{-\mathrm{\γ}}\frac{\mathrm{\ω}\mathrm{\Δ}\mathrm{\τ}}{2Q_r}\]\exp \[i{\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_h}\right)}^{-\mathrm{\γ}}\mathrm{\ω}\mathrm{\Δ}\mathrm{\τ}\]---\left(3\right)$

式中两个指数项分别为振幅效应(即能量吸收)的补偿算子和相位效应(即速度频散)的校正算子。(3)式是地震波补偿(反Q滤波)算法的基本表达式。

考虑地层Q模型为随旅行时变化的一维函数Q(τ)，利用(3)式对波场自地表τ₀＝0向下延拓到深度时间τ，此时的波场可以表示为：

$\begin{array}{c}U\left(\mathrm{\τ},\mathrm{\ω}\right)=U\left(0,\mathrm{\ω}\right)\exp \[{\∫}_0^{\mathrm{\τ}}{\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_h}\right)}^{-\mathrm{\γ}\left({\mathrm{\τ}}^{\′}\right)}\frac{\mathrm{\ω}}{2Q_r\left({\mathrm{\τ}}^{\′}\right)}{\mathrm{d\τ}}^{\′}\]\\ \×\exp \[i{\∫}_0^{\mathrm{\τ}}{\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_h}\right)}^{-\mathrm{\γ}\left({\mathrm{\τ}}^{\′}\right)}{\mathrm{\ωd\τ}}^{\′}\],\end{array}---\left(4\right)$

其中γ(τ)＝(1/π)Q^-1(τ)，dτ′为τ的积分参数。从(4)式可看出，振幅补偿算子是频率和旅行时的指数函数，包括振幅补偿的地震波补偿会引起数值不稳定，会在地震数据中产生噪音。因此，在任何地震波补偿算法中数值稳定性都是值得关注的主要问题。

为了稳定整个计算过程，我们将(4)式写成：

$\mathrm{\β}(\mathrm{\τ},\mathrm{\ω})U(\mathrm{\τ},\mathrm{\ω})=U(0,\mathrm{\ω})\exp \[i{\∫}_0^{\mathrm{\τ}}{\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_h}\right)}^{-\mathrm{\γ}\left({\mathrm{\τ}}^{\′}\right)}{\mathrm{\ωd\τ}}^{\′}\]---\left(5\right)$

其中，

$\mathrm{\β}(\mathrm{\τ},\mathrm{\ω})=\exp \[-{\∫}_0^{\mathrm{\τ}}{\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_h}\right)}^{-\mathrm{\γ}\left({\mathrm{\τ}}^{\′}\right)}\frac{\mathrm{\ω}}{2Q_r\left({\mathrm{\τ}}^{\′}\right)}{\mathrm{d\τ}}^{\′}\]---\left(6\right)$

然后将方程(5)作为一个反问题来进行稳定化求解，则生成了如下的稳定化的表达式：

$U(\mathrm{\τ},\mathrm{\ω})=U(0,\mathrm{\ω})\mathrm{\Λ}(\mathrm{\τ},\mathrm{\ω})\exp \[i{\∫}_0^{\mathrm{\τ}}{\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_h}\right)}^{-\mathrm{\γ}\left({\mathrm{\τ}}^{\′}\right)}{\mathrm{\ωd\τ}}^{\′}\]---\left(7\right)$

其中，

$\mathrm{\Λ}(\mathrm{\τ},\mathrm{\ω})=\frac{\mathrm{\β}(\mathrm{\τ},\mathrm{\ω})+{\mathrm{\σ}}^2}{{\mathrm{\β}}^2(\mathrm{\τ},\mathrm{\ω})+{\mathrm{\σ}}^2}---\left(8\right)$

σ²为稳定因子。

将(7)式应用于所有频率，并且将这些平面波叠加(即成像条件)，我们得到一个时间域的地震样点：

$U\left(\mathrm{\τ}\right)=\frac{1}{\mathrm{\π}}{\∫}_0^{\mathrm{\∞}}U(0,\mathrm{\ω})\mathrm{\Λ}(\mathrm{\τ},\mathrm{\ω})\exp \[i{\∫}_0^{\mathrm{\τ}}{\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_h}\right)}^{-\mathrm{\γ}\left({\mathrm{\τ}}^{\′}\right)}{\mathrm{\ωd\τ}}^{\′}\]d\mathrm{\ω}---\left(9\right)$

这是2006年王仰华提出的稳定完备地震波补偿的表达式，通过此表达式实现地震波衰减补偿。

图1为本发明提供的频率驱动空变Q值模型构建的地震波补偿方法的流程图，如图1所示，频率驱动空变Q值模型构建的地震波补偿方法包括：

步骤S101，通过加伯变换生成对应于时间域地震数据每个样点的时间频率谱。

原始单炮地震数据经过精细预处理和偏移成像后生成叠后地震数据u(t)。叠后地震数据体是经过叠加、偏移等若干精细处理流程之后得到的，相对于叠前数据而言数据量小，计算量小。加伯变换可以从单个的时间域的叠后地震数据u(t)来计算其时间频率域

加伯变换公式为：

$\tilde{U}(\mathrm{\τ},\mathrm{\ω})={\∫}_{-\mathrm{\∞}}^{\mathrm{\∞}}u\left(t\right)w(t-\mathrm{\τ})\exp \[-i\mathrm{\ω}t\]dt---\left(10\right)$

上式中w(t)是加伯分析窗，其中τ为该窗的中心位置。可定义加伯切片为：

$\tilde{u}(\mathrm{\τ},t)=u\left(t\right)w(t-\mathrm{\τ})---\left(11\right)$

叠后地震数据u(t)的加伯变换可看成是加伯切片对所有可能τ位置关于时间t的Fourier变换：

$\tilde{U}(\mathrm{\τ},\mathrm{\ω})=F\left\{\tilde{u}(\mathrm{\τ},t)\right\}---\left(12\right)$

加伯分析窗为高斯窗：

其中T为半窗宽度。

步骤S102，计算每个样点的中心频率，生成所有样点的中心频率数据体或剖面。

中心频率ω_c指的是频带中央所在频率。在实际计算中，我们用中值滤波法从相邻5×5(或3×3)道的样值中选取中值，可实现中心频率求取的稳定化。图3显示了对应图2地震剖面的中心频率剖面。

步骤S103，根据所述每个样点的中心频率估算每个样点的衰减值Q^-1，品质因子Q值即衰减值Q^-1的倒数，由此构建Q值模型。

在参照时间点τ₀，假设子波可以用雷克子波表示：

$R\left(\mathrm{\τ}\right)=(1-\frac{1}{2}{\mathrm{\ω}}_0^2{\mathrm{\τ}}^2)\exp (-\frac{1}{4}{\mathrm{\ω}}_0^2{\mathrm{\τ}}^2)---\left(14\right)$

其中ω₀是子波在初始时间位置τ＝τ₀时的峰值频率，该子波的均值为零,雷克子波的频谱为：

$R\left(\mathrm{\ω}\right)=\frac{2{\mathrm{\ω}}^2}{\sqrt{\mathrm{\π}}{\mathrm{\ω}}_0^2}\exp \left(-\frac{{\mathrm{\ω}}^2}{{\mathrm{\ω}}_0^2}\right)---\left(15\right)$

将对应于最大振幅值的频率称之为峰值频率ω_p，我们可设定频谱对频率的导数为零:

$\frac{\∂A(\mathrm{\ω},\mathrm{\τ})}{\∂\mathrm{\ω}}=0---\left(16\right)$

在任意时间τ位置振幅谱可以表示为：

$A(\mathrm{\ω},\mathrm{\τ})=\frac{2}{\sqrt{\mathrm{\π}}}\frac{{\mathrm{\ω}}^2}{{\mathrm{\ω}}_0^3}\exp \left(-\frac{{\mathrm{\ω}}^2}{{\mathrm{\ω}}_0^2}\right)\exp \left(-\frac{\mathrm{\ω}(\mathrm{\τ}-{\mathrm{\τ}}_0)}{2Q}\right)---\left(17\right)$

导数为零的公式(15)将会生成下列二次方程：

$\frac{{\mathrm{\ω}}^2}{{\mathrm{\ω}}_0^2}+\frac{(\mathrm{\τ}-{\mathrm{\τ}}_0)\mathrm{\ω}}{4Q}-1=0---\left(18\right)$

将其解析为：

${\mathrm{\ω}}_p\left(\mathrm{\τ}\right)={\mathrm{\ω}}_0^2\left(\sqrt{\frac{1}{{\mathrm{\ω}}_0^2}+{\left(\frac{\mathrm{\τ}-{\mathrm{\τ}}_0}{8Q}\right)}^2}-\frac{\mathrm{\τ}-{\mathrm{\τ}}_0}{8Q}\right)---\left(19\right)$

子波在初始时间位置τ＝τ₀时，峰值频率ω_p(τ₀)＝ω₀。

对于雷克子波而言，根据峰值频率与中心频率的理论关系ω_P＝0.9442ω_c，得到峰值频率ω_p之后，可估算地层衰减值Q^-1。对上述子波由(17)式推导出：

$Q^{-1}(\mathrm{\τ}-{\mathrm{\τ}}_0)=\frac{4}{(\mathrm{\τ}-{\mathrm{\τ}}_0){\mathrm{\ω}}_0}\left(\frac{{\mathrm{\ω}}_0}{{\mathrm{\ω}}_p\left(\mathrm{\τ}\right)}-\frac{{\mathrm{\ω}}_p\left(\mathrm{\τ}\right)}{{\mathrm{\ω}}_0}\right)---\left(20\right)$

式中ω₀＝ω_p(τ₀),子波在初始时间点τ₀处的峰值频率。

从上式可以看出，层Q值在[τ₀,τ]时段取决于两个基本元素，一个是ω₀，即时窗起点处的频率，二是ω₀/ω_p,即时窗两端点处峰值频率之比值。图4从剖面方向显示了Q值的空间变化，图5显示了Q值在近地表的横向变化。

步骤S104，利用加伯变换对频率域实现地震波衰减补偿。

从(7)式开始进行推导，假设地震波传播的介质满足Q^-1＝0，在深度时间τ记录到的波场可以由(7)得到：

$\tilde{U}(\mathrm{\τ},\mathrm{\ω})=U(0,\mathrm{\ω})\exp \[i\mathrm{\ω}\mathrm{\τ}\]---\left(21\right)$

利用(7)式变为：

$U(\mathrm{\τ},\mathrm{\ω})=\tilde{U}(\mathrm{\τ},\mathrm{\ω})\mathrm{\Λ}(\mathrm{\τ},\mathrm{\ω})\exp \[i{\∫}_0^{\mathrm{\τ}}\left({\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_h}\right)}^{-\mathrm{\γ}\left({\mathrm{\τ}}^{\′}\right)}-1\right){\mathrm{\ωd\τ}}^{\′}\]---\left(22\right)$

上式是基于加伯变换的地震波补偿(反Q滤波)算法的核心公式。U(τ,ω)和都是在深度时间τ记录到的波场，区别在于是没有进行过地震波，而如(10)式是通过加伯变换由叠后地震数据u(t)转换得到。

步骤S105,利用加伯反变换将补偿后的频率域数据转换重构到时间域。

在加伯变换方法中，当对波场进行地震波补偿之后，我们可用反加伯变换来重构时间域的地震道，即经过地震波补偿的地震道。

从加伯变换域U(τ,ω)中重构叠后地震数据u(t)的反加伯变换定义为：

$u\left(t\right)=h\left(t\right){\∫}_{-\mathrm{\∞}}^{\mathrm{\∞}}{\∫}_{-\mathrm{\∞}}^{\mathrm{\∞}}U(\mathrm{\τ},\mathrm{\ω})\exp \[i\mathrm{\ω}t\]d\mathrm{\ω}d\mathrm{\τ}---\left(23\right)$

其中h(t)是加伯合成窗，合成窗是为了减少当在地震道上滑动分析时窗时，由加伯分析窗的数字化和边缘效应导致的指数部分的数字误差。经推导：

$h\left(t\right)={\[{\∫}_{-\mathrm{\∞}}^{\mathrm{\∞}}w\left(t-\mathrm{\τ}\right)d\mathrm{\τ}\]}^{-1}---\left(24\right)$

图6为地震数据进行地震波衰减补偿前的地震剖面，图7为地震数据进行地震波衰减补偿后的地震剖面。图6和图7的对比说明了利用本发明提取的Q值的地震波补偿效果，经过补偿后，地震数据的信噪比和分辨率都有一定程度的提高，断裂构造的成像也更清晰。

图8为地震数据进行地震波衰减补偿前的频谱对比，图9为地震数据进行地震波衰减补偿后的频谱对比。图8和图9的对比说明了经过本发明提取的Q值的地震波补偿效果，经过补偿后，地震数据的主频提高了约10Hz,频带也拓宽了大概10Hz。

根据本发明提供的频率驱动空变Q值模型构建的地震波补偿方法，通过加伯变换从叠后地震数据转换得到每个样点的频率信息进而逐点构建空变Q值模型，计算量小，提取Q值方法稳定，Q值模型纵向和横向上连续变化，符合实际地质情况。在高精度的Q值模型基础上，通过地震波补偿不仅尽可能地恢复了振幅信息和纠正了相位畸变，还提高了地震数据的分辨率和信噪比。

上述说明示出并描述了本申请的若干优选实施例，但如前所述，应当理解本申请并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文所述发明构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离本申请的精神和范围，则都应在本申请所附权利要求的保护范围内。

Claims (5)

1.一种频率驱动空变Q值模型构建的地震波补偿方法，其特征在于，包括如下步骤：

(1)通过加伯变换生成对应于时间域地震数据每个样点的时间频率谱；

(2)计算每个样点的中心频率，生成所有样点的中心频率数据体或剖面；

(3)根据所述每个样点的中心频率估算每个样点的衰减值Q^-1，品质因子Q值即衰减值Q^-1的倒数，由此构建Q值模型；

(4)在频率域实现地震波衰减补偿；

(5)利用加伯反变换将补偿后的频率域数据转换重构到时间域。

2.根据权利要求1所述的方法，其特征在于，所述步骤(1)中加伯变换可逆，通过时间频率谱反变换能恢复原始时间域地震数据。

3.根据权利要求1所述的方法，其特征在于，所述步骤(1)中使用空间上相邻的地震道组合，用于计算各个地震道的时间频率谱。

4.根据权利要求1所述的方法，其特征在于，所述步骤(2)中对地震数据中相邻空间位置的时间频率谱进行组合，生成所有样点的中心频率数据体或剖面。

5.根据权利要求1所述的方法，其特征在于，所述步骤(3)还包括：由中心频率估算衰减值Q^-1时，结合相邻地震道，使用中值滤波方法获取所述时窗起点处的中心频率和所述时窗起点的中心频率与终点的中心频率的比值。