CN103728661B - 一种高精度反q滤波地震资料处理方法 - Google Patents

Abstract

本发明提供了一种高精度反Q滤波地震资料处理方法，属于地震信号处理领域。本发明方法包括：(1)，从原始地震道信号u₀(τ)的频率域中提取品质因子Q值；(2)，通过Gabor正变换将地震信号u₀(τ)转换成时频谱其中τ为旅行时，ω为角频率；(3)，对所述时频谱进行反Q滤波，得到新的时频谱U(τ，ω)；(4)，通过Gabor反变换，将U(τ，ω)转换成反Q滤波补偿后的时间域地震道信号u(τ)。利用本发明不仅能够实现全频Q值稳定补偿，而且能够使地震振幅信息尽可能地恢复，在振幅补偿方面，是一种有效的地震信号处理方法。

Description

一种高精度反Q滤波地震资料处理方法

技术领域

本发明属于地震信号处理领域，具体涉及一种高精度反Q滤波地震资料处理方法。

背景技术

当今石油勘探的主要地质目标已经从简单构造油气藏变成了隐蔽性油气藏和复杂构造油气藏，这对勘探的精度提出了更高的要求。然而，地震波在地下介质中传播会产生吸收衰减现象，其振幅随着传播时间路径增大而快速减弱，从而降低了地震资料的信噪比和分辨率。为了解决这种现象，提出了反Q滤波，它能够恢复高频信息，提高弱反射波的能量和地震资料的信噪比、分辨率。

针对反Q滤波算法，国内外专家、学者进行了大量研究，如Hale提出的基于Futterman数学模型反Q滤波、McCarley提出的基于常数Q值滤波的自适应模型算法、Hargreaves与Calvert提出的基于常数Q值模型的相位反Q滤波、裴江云提出的基于Kjartansson模型的反Q滤波、王仰华提出的一个全局稳定有效反Q滤波等等。其算法大致可以分成三大类：用级数展开作近似高频补偿的反Q滤波方法、基于波场延拓的反Q滤波方法和其他反Q滤波方法。相比较而言，基于波场延拓的反Q滤波方法由于使用了快速傅里叶变换，计算速度较快，生产效率高，能够有效地校正由频散引起的相位畸变，而且无条件稳定，但是在振幅补偿方面却存在一定的不稳定性，还有一些尚待解决的问题，需要进一步研究探讨，以满足复杂地震资料高分辨率处理的要求。本发明属于基于波场延拓的反Q滤波方法范畴，下面介绍该类型反Q滤波基本原理。

基于波动理论，平面波U(x，ω)沿波前方向传播距离Δx后，当前能量为：

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

其中，i为单位虚数，ω为角频率，k(ω)为复数波数。

当引入品质因子Q，k(ω)变形为：

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

其中，v_γ为参考频率对应的相速度，为经验公式，ω_h为调谐主频。

由于反Q滤波是波前传播逆过程，将(2)式代入(1)式，同时用Δτ替代Δx，得出常规反Q滤波基本公式(3)：

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

定义QP、QA分别为相位项与振幅项，有：

$\mathrm{QP}\left(\mathrm{\ω}\right)=\exp \left[i{\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_h}\right)}^{-\mathrm{\γ}}\mathrm{\ω\Δ\τ}\right]---\left(4\right)$

$\mathrm{QA}\left(\mathrm{\ω}\right)=\exp \left[{\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_h}\right)}^{-\mathrm{\γ}}\frac{\mathrm{\ω\Δ\τ}}{2Q}\right]---\left(5\right)$

如是，(3)式可以表示为：

U(τ+Δτ，ω)＝U(τ，ω)×QA(ω)×QP(ω)(6)

根据前人研究结果，在反Q滤波计算过程中，相位项QP恒稳定，但是振幅项QA出现不稳定现象。

因此，为了克服常规反滤波不稳定性，许多学者做了大量研究和试验，得出一个恒稳定的条件公式：

$\mathrm{QA}\left(\mathrm{\ω}\right)=\exp \left({\left|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_0}\right|}^{-\mathrm{\γ}}\frac{\mathrm{\ω}}{2Q}\mathrm{\Σ\ΔT}\right)\≤e---\left(7\right)$

对应一个恒稳定频率范围：

即， $\mathrm{QA}\left(\mathrm{\ω}\right)=\left\{\begin{array}{c}\exp \left({\left|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_0}\right|}^{-\mathrm{\γ}}\frac{\mathrm{\ω\ΔT}}{2Q}\right),\mathrm{\ω}\≤{\mathrm{\ω}}_q\\ \mathrm{QA}\left({\mathrm{\ω}}_q\right),\mathrm{\ω}>{\mathrm{\ω}}_q\end{array}\right.---\left(8\right)$

该方法被称为增益限制反Q滤波算法。可以看出，该公式类似于带通滤波，即截取高频段振幅项为常数不变。

为了既能克服常规反Q滤波算法不稳定性，又能尽可能地恢复地震波振幅，学者王仰华(2006年)提出了一个数学近似的稳定全局算法公式：

$\mathrm{QA}\left(\mathrm{\ω}\right)=\frac{1}{\mathrm{\β}(\mathrm{\τ},\mathrm{\ω})}\≈\frac{\mathrm{\β}(\mathrm{\τ},\mathrm{\ω})+{\mathrm{\σ}}^2}{{\mathrm{\β}}^2(\mathrm{\τ},\mathrm{\ω})+{\mathrm{\σ}}^2}---\left(9\right)$

其中，且σ²为正常数，取值范围为10^-2～10^-10。

该方法与2009年甘其刚等人发明专利“一种对地震波信号进行反Q滤波的方法”相同，后者的表达式为：

$U(\mathrm{\τ},\mathrm{\ω})=\tilde{U}(\mathrm{\τ},\mathrm{\ω})\×(\mathrm{\τ},\mathrm{\ω})\×\exp \left[i{\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_h}\right)}^{-\mathrm{\γ}}\mathrm{\ω\Δ\τ}\right]---\left(12\right)$

其稳定化处理方--振幅项为：

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

其中，且σ²为正常数，取值范围为10^-2～10^-10。

相比较增益限制算法，从理论上得出稳定全局算法能够获取更多高频信息，不仅克服了常规算法的不稳定性，而且振幅恢复效果有一定程度提高。

但是，可以看出，当β(τ，ω)值小于0.03时，公式(9)式数学近似计算误差较大，例如当品质因子Q值较小或地震波传播时间较大时，导致稳定全局算法反Q滤波振幅恢复效果不理想。因而，需要一种更高精度反Q滤波算法，尽可能地恢复更多高频地震振幅信息。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种高精度反Q滤波地震资料处理方法，通过合理选择参数，实现全频Q值稳定补偿，同时在高信噪比地震资料处理中能够使地震振幅信息尽可能地恢复。

本发明是通过以下技术方案实现的：

一种高精度反Q滤波地震资料处理方法，所述方法包括以下步骤：

(1)，从原始地震道信号u₀(τ)的频率域中提取品质因子Q值；

(2)，通过Gabor正变换将地震信号u₀(τ)转换成时频谱其中τ为旅行时，ω为角频率；

(3)，利用下式对所述时频谱进行反Q滤波，得到新的时频谱U(τ，ω)：

$U(\mathrm{\τ},\mathrm{\ω})=\tilde{U}(\mathrm{\τ},\mathrm{\ω})\×\mathrm{QA}\left(\mathrm{\ω}\right)\×\exp \left[i{\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_h}\right)}^{-\mathrm{\γ}}\mathrm{\ω\Δ\τ}\right],$

且振幅项 $\mathrm{QA}\left(\mathrm{\ω}\right)=\left\{\begin{array}{c}\exp \left(k\right),k\≤k0\\ \frac{\mathrm{\α}+l}{{\mathrm{\α}}^2+l},k>k0\end{array}\right.,$ 其中 $k={\left|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_h}\right|}^{-\mathrm{\γ}}\frac{\mathrm{\ω\ΔT}}{2Q}$

α＝exp(-(k-k0))，γ＝1/(πQ)，ω_h为主频，i为虚部单位，指数参数k0与稳定系数l均为正常数，ΔT表示地震波传播时间深度；

(4)，通过Gabor反变换，将U(τ，ω)转换成反Q滤波补偿后的时间域地震道信号u(τ)。

所述步骤(3)中的指数参数k0的取值范围为6～18。

所述步骤(3)中的稳定系数l的取值范围为10^-2～10^-10。

所述步骤(1)中的提取品质因子Q值采用的是谱比法。

所述步骤(2)中的Gabor变换采用的是短时窗Gabor变换。

与现有技术相比，本发明的有益效果是：本发明通过合理选择参数，不仅能够实现全频Q值稳定补偿，而且使地震振幅信息尽可能地恢复，在振幅补偿方面，是一种有效的地震信号处理方法。

附图说明

图1是本发明高精度反Q滤波地震资料处理方法的步骤框图；

图2(a)是单道常数Q值模型正演地震道；

图2(b)是对图2(a)数据做基础反Q滤波处理；

图2(c)是对图2(a)数据做增益限制反Q滤波处理；

图2(d)是对图2(a)数据做稳定全局反Q滤波处理；

图2(e)是对图2(a)数据做本发明反Q滤波处理；

图3是图2(a)数据0.95秒处反Q滤波前后3种方法能量误差对比；

图4(a)是单道多层Q值模型示意图；

图4(b)是图4(a)模型正演Q吸收地震道；

图4(c)是图4(b)地震道Gabor时频谱；

图5(a)是对图4(b)地震道做增益限制反Q滤波处理；

图5(b)是对图4(b)地震道做稳定全局反Q滤波处理；

图5(c)是对图4(b)地震道做本发明反Q滤波处理；

图6(a)是图5(a)地震道Gabor时频谱；

图6(b)是图5(b)地震道Gabor时频谱；

图6(c)是图5(c)地震道Gabor时频谱；

图7(a)是反Q滤波处理前实际资料地震叠加剖面；

图7(b)是稳定全局算法反Q滤波实际资料处理后地震叠加剖面；

图7(c)是本发明算法反Q滤波实际资料处理后地震叠加剖面；

图8(a)是第20道反Q滤波处理前后地震道波形显示对比；

图8(b)-1是第20道反Q滤波处理前原始地震第20道的时频谱后相应地震道显示对比；

图8(b)-2是采用现有技术对图8(b)-1的原始地震第20道进行处理后得到的Gabor时频谱；

图8(b)-3是采用本发明方法对图8(b)-1的原始地震第20道进行处理后得到的Gabor时频谱；

图9(a)是第40道反Q滤波处理前后地震道波形显示对比；

图9(b)-1是第40道反Q滤波处理前原始地震第20道的时频谱后相应地震道显示对比；

图9(b)-2是采用现有技术对图8(b)-1的原始地震第40道进行处理后得到的Gabor时频谱；

图9(b)-3是采用本发明方法对图8(b)-1的原始地震第40道进行处理后得到的Gabor时频谱；

具体实施方式

下面结合附图对本发明作进一步详细描述：

针对目前常用的稳定全局算法中的振幅项在反Q滤波提高振幅存在局限性条件下，本文发明提供了一种高精度反Q滤波地震资料处理方法，对地震信号进行高精度振幅项补偿。

如图1所示，本发明提供的高精度反Q滤波地震资料处理方法包含：

首先，从原始地震道信号频率域中提取品质因子Q值；反Q滤波处理效果好坏与Q值的准确拾取有关。而品质因子Q的提取方法有很多，如谱比法、质心偏移法、解析法、基于衰减函数法、基于振幅补偿Q值提取法(2009年甘其刚等人发明专利“一种对地震波信号进行反Q滤波的方法”)等等，本发明采用谱比法，具体如下：基于相邻2个接收器之间相邻反射时间对应的地层Q值是一个常量为假设，该方法要求选取的反射波同相轴完整“干净”，特点是简单实用，具体步骤如下：

选取相邻道相邻反射时间t₁、t₂对应的同相轴振幅A(t₁)、A(t₂)；

分别对A(t₁)、A(t₂)做傅立叶变换，得出相应的频谱A(f，t₁)、A(f，t₂)；

利用谱比法原理，建立关系式：

$A(f,t_2)=A(f,t_1)e^{-\mathrm{\πf}(t_2-t_1)/Q},$

将方程两边求取自然对数，得到线性关系式：

$\ln \left[\frac{A(f,t_2)}{A(f,t_1)}\right]=-\mathrm{\πf}(t_2-t_1)/Q,$

对多个频率f，进行线性拟合，计算出斜率m，最后换算成相应Q值：

Q＝-π(t₂-t₁)/m，即为所提取的品质因子Q值。

然后，通过Gabor变换，将地震信号转换成时频谱其中τ为旅行时，ω为角频率；具体实施时，可采用短时窗Gabor变换，具体公式如下：

$\tilde{U}(\mathrm{\τ},\mathrm{\ω})=\underset{-\∞}{\overset{+\∞}{\∫}}{u}_0\left(t\right)w(t-\mathrm{\τ})\exp (-\mathrm{i\ωt})\mathrm{dt}$

其中w(t)为Gabor时窗，τ为该窗中心位置。

其次，考虑到大地介质吸收，在波场延拓波动方程算法上引入品质因子Q值，建立高精度振幅项指数因子，对地震信号进行反Q滤波，得到新的时频谱U(τ，ω)：

$U(\mathrm{\τ},\mathrm{\ω})=\tilde{U}(\mathrm{\τ},\mathrm{\ω})\×\mathrm{QA}\left(\mathrm{\ω}\right)\×\exp \left[i{\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_h}\right)}^{-\mathrm{\γ}}\mathrm{\ω\Δ\τ}\right]---\left(10\right)$

且振幅项 $\mathrm{QA}\left(\mathrm{\ω}\right)=\left\{\begin{array}{c}\exp \left(k\right),k\≤k0\\ \frac{\mathrm{\α}+l}{{\mathrm{\α}}^2+l},k>k0\end{array}\right.,$ 其中 $k={\left|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_h}\right|}^{-\mathrm{\γ}}\frac{\mathrm{\ω\ΔT}}{2Q}$ (两者的乘积)，α＝exp(-(k-k0))，γ＝1/(πQ)，ω_h为主频，i为虚部单位，k0与l均为正常数；

最后，通过Gabor反变换，将U(τ，ω)转换成反Q滤波补偿后的时间域地震道信号u(t)，具体公式如下：

$u\left(t\right)=h\left(t\right)\underset{-\∞}{\overset{+\∞}{\∫}}U(\mathrm{\τ},\mathrm{\ω})\exp \left(\mathrm{i\ωt}\right)\mathrm{d\ω}$

$h\left(t\right)=1/\underset{-\∞}{\overset{+\∞}{\∫}}w(t-\mathrm{\τ})\mathrm{d\τ}$

其中h(t)为Gabor合成窗。

下面分别从理论模型数据和实际地震资料应用中比较本发明反Q滤波处理效果。

图2(a)至图2(e)和图3说明了本发明的反Q滤波方法进行常数Q值模型数据处理的效果图与其他方法对比。其中，图2(a)为常数Q值模型，子波为零相位雷克子波，而图2(b)-图2(e)分别对应4种不同反Q滤波算法的处理效果。从图中可以看出，相比较其他3种方法，本发明振幅恢复效果最佳。图3为0.95s时间处的能量误差对比，同样显示出本发明振幅恢复误差最小。

图4(a)至图4(c)、图5(a)至图5(c)以及图6(a)至图6(c)说明了本发明的反Q滤波方法进行多层Q值模型数据处理的效果图与其他方法对比。其中，图4(a)为多层Q值模型示意图，选择主频25Hz雷克子波做正演，而图4(b)为正演Q值吸收地震道，相应图4(c)为其时频谱图，可以看出，随着波传播时间增大，地震波振幅受Q值吸收影响逐步减弱。分别对该模型做3种反Q滤波方法处理，所采取的Q值为模型真实值，通过对比地震振幅恢复效果(图5(a)至图5(c))与相应的时频谱(图6(a)至图6(c))，本发明的反Q滤波尽可能地将地震振幅信息恢复到最佳效果。

图7(a)至图7(c)、图8(a)至图8(b)-3和图9(a)至图9(b)-3说明了本发明的反Q滤波方法在实际地震资料处理的效果图以及不同反Q滤波方法对比。其中，图7(a)为输入的原始地震叠加剖面，根据图1所示处理流程对其进行反Q滤波信号处理。其中图7(b)为稳定全局方法处理效果，图7(c)为本发明的反Q滤波处理效果，相比较，本发明处理的剖面地震振幅恢复更好。另外，抽取单个地震道进行方法对比，图8(a)至图8(b)-3、图9(a)至图9(b)-3分别为第20道、第40道单道数据反Q滤波处理前后对比图以及相应的时频谱，再次说明了本发明中的反Q滤波是一种更有效的高精度振幅补偿方法。

由以上描述的模型数据和实际资料数据反Q滤波的处理效果分析，可以看出，相比较目前的稳定全局方法(2009年甘其刚等人发明的“一种对地震波信号进行反Q滤波的方法”，专利申请号200910236232.6)，本发明实现了更高精度振幅项补偿的地震信号处理。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (5)

1.一种高精度反Q滤波地震资料处理方法，其特征在于：所述方法包括以下步骤：

(1)，从原始地震道信号u₀(τ)的频率域中提取品质因子Q值；

(2)，通过Gabor正变换将地震道信号u₀(τ)转换成时频谱其中τ为旅行时，ω为角频率；

(3)，利用下式对所述时频谱进行反Q滤波，得到新的时频谱U(τ,ω)：

$U(\mathrm{\τ},\mathrm{\ω})=\tilde{U}(\mathrm{\τ},\mathrm{\ω})\×QA\left(\mathrm{\ω}\right)\×\exp \[i{\left(\frac{\mathrm{\ω}}{{\mathrm{\ω}}_h}\right)}^{-\mathrm{\γ}}\mathrm{\ω}\mathrm{\Δ}\mathrm{\τ}\],$

且振幅项 $QA\left(\mathrm{\ω}\right)=\left\{\begin{array}{c}\exp \left(k\right),k\≤k0\\ \frac{\mathrm{\α}+l}{{\mathrm{\α}}^2+l},k>k0\end{array}\right.,$ 其中 $k=|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_h}{|}^{-\mathrm{\γ}}\frac{\mathrm{\ω}\mathrm{\Δ}T}{2Q}$

α＝exp(-(k-k0))，γ＝l/(πQ)，ω_h为主频，i为虚部单位，指数参数k0与稳定系数l均为正常数，ΔT表示地震波传播时间深度；

(4)，通过Gabor反变换，将U(τ,ω)转换成反Q滤波补偿后的时间域地震道信号u(τ)。

2.根据权利要求1所述的高精度反Q滤波地震资料处理方法，其特征在于：所述步骤(3)中的指数参数k0的取值范围为6～18。

3.根据权利要求2所述的高精度反Q滤波地震资料处理方法，其特征在于：所述步骤(3)中的稳定系数l的取值范围为10^-2～10^-10。

4.根据权利要求1至3任一所述的高精度反Q滤波地震资料处理方法，其特征在于：所述步骤(1)中的提取品质因子Q值采用的是谱比法。

5.根据权利要求1至3任一所述的高精度反Q滤波地震资料处理方法，其特征在于：所述步骤(2)中的Gabor变换采用的是短时窗Gabor变换。