CN103412324B - 一种估计介质品质因子的epifvo方法 - Google Patents

Abstract

本发明公开了一种估计介质品质因子的EPIFVO方法：当震源子波传播到第j-1层后，根据惠更斯原理，将该子波波前作为子波源，向各方向继续传播，经过时间Δt_j/2后子波到达第j层并反射，在第j-1层接收到反射波；f_p(t_j-1)为第j-1层子波源的EPIF，f_p(t_j)为接收到的反射波的EPIF；利用式求出第j层的层间Q值。该方法采用具有4个待定系数的函数去逼近震源子波，利用粘弹介质中单程波传播理论推导出地震子波包络峰值处的瞬时频率(EPIF)与不同偏移距处旅行时的关系；利用该关系，外推出该同相轴零偏移距处的EPIF，用小波域包络峰值处瞬时频率法(WEPIF)结合层位信息估计Q值。

Description

一种估计介质品质因子的EPIFVO方法

技术领域

本发明属于地震勘探技术领域，涉及估计介质品质因子的方法，尤其是一种估计介质品质因子的EPIFVO方法。

背景技术

地层是粘弹性介质，地震波在地下传播过程中，由于介质的吸收效应，致使产生波的衰减和频散。Q值是度量介质衰减特性的重要参数，也是指示地层含油气性的重要属性之一。由地震资料估计的地层衰减（Q值）不仅可用于岩性及含流体分析、储层识别和烃类检测等，也可用于对地震资料进行吸收补偿，提高地震资料分辨率等。

目前常用的估计Q值的地震资料有VSP资料、井间资料、叠后反射地震资料、叠前反射地震资料等。VSP资料信噪比高但成像范围小；井间地震资料能提供高精度的井间储层信息但其费用昂贵；地面地震资料在横向上可在很大范围内追踪地层厚度和岩性变化，但地面地震资料的纵向分辨率相对较低，且是远离地层界面的间接观测。叠后反射地震资料的信噪比高，但是经过NMO拉伸却损坏了频率信息；叠前反射地震资料虽然信噪比较低，但是比叠后反射地震资料有更加精确的地层信息及丰富的振幅和旅行时信息，频率无损坏，一些细微的地层特征在叠后地震资料上是看不到的，尤其是薄互层存在时，叠后地震资料很难区分散射与吸收。鉴于叠前资料丰富的振幅和频率信息，我们有必要研究基于叠前反射地震资料的介质品质因子估计方法。

已有许多学者对叠前地震资料的衰减估计技术进行了研究。王小杰和吴国忱等提出了基于叠前资料利用时频谱分解技术估算地层吸收参数的方法，文章基于地震子波为一般零相位子波的假设，用质心频率偏移法、斜率法、谱比法实现了衰减估计，然而文章没有讨论射线几何路径不同对吸收参数的影响，即具有不同入射角的射线穿过多层介质时，如何拾取层间旅行时的问题。文章采用频率域方法估计Q值，这些方法大多需加时窗截取地震记录并用Fourier变换计算频谱,而在实际问题中恰当地选择窗函数的类型及长度是比较困难的。Zhangchangjun和UlrychTJ提出了基于叠前CMP资料且用PFS法(峰值频率移动法)估计Q值的方法，即结合先验的层位信息拾取该层位对应的地震道的非零偏移距处的峰值频率，然后用PFVO(峰值频率随偏移距变化关系)方法外推出该同相轴零偏移距处的峰值频率，进而用PFS法计算出Q曲线。该方法把叠前Q值提取问题转换为类似叠后Q值提取问题，解决了直接利用叠前资料时旅行时难提取的问题。zhang提出的PFS法假设震源子波是零相位的Ricker子波，当震源子波满足假设条件时估计的Q值精度较高，然而在实际资料中，震源子波通常不是零相位的，导致该方法估计的Q值误差较大；另外zhang的方法的实现步骤中，用包络峰值处瞬时频率EPIF代替峰值频率PF，二者的定义不同，求取方法不同，稳定性不同，可能会导致误差。高静怀和杨森林等提出用WEPIF(子波包络峰值处瞬时频率)方法估计Q值，该方法利用子波包络峰值处瞬时频率与衰减之间的关系来估计Q值，有更高的纵向分辨率，能避免加时窗问题，抗噪性好而且更易于实现。基于上述讨论，本发明发展WEPIF方法，研究适合一般地震子波的叠前CMP资料衰减估计技术。

发明内容

本发明的目的在于克服上述现有技术的缺点，提供一种估计介质品质因子的EPIFVO方法，该方法采用具有4个待定系数的函数去逼近震源子波，利用粘弹介质中单程波传播理论推导出地震子波包络峰值处的瞬时频率(EPIF)与不同偏移距处旅行时的关系；利用该关系，外推出该同相轴零偏移距处的EPIF，用小波域包络峰值处瞬时频率法(WEPIF)结合层位信息估计Q值。

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

这种估计介质品质因子的EPIFVO方法为：当震源子波传播到第j-1层后，根据惠更斯原理，将该子波波前作为子波源，向各方向继续传播，经过时间Δt_j/2后子波到达第j层并反射，在第j-1层接收到反射波；f_p(t_j-1)为第j-1层子波源的EPIF，f_p(t_j)为接收到的反射波的EPIF；利用(14)式求出第j层的层间Q值：

$\frac{1}{Q_j}=\frac{1}{\mathrm{\Δ}{t}_j}\frac{4\mathrm{\π}(f_p\left(t_{j-1}\right)-f_p\left(t_j\right))}{{\mathrm{\δ}}^{2}k\left(\mathrm{\η}\right)}---\left(14\right);$

其中下标j表示第j层介质。

进一步的，以上在进行计算层间Q值前，首先进行层位信息拾取，包括以下步骤：

(1)选取近偏移距中的N道，N大于1，计算这N道的瞬时振幅IA和瞬时频率IF，由IA和IF拾取EPIF及其位置α，将EPIF的位置等效为同相轴的位置，EPIF的位置为一个列向量；

(2)从近偏移距的N道中选取参考道，根据式(15b)，求参考道的EPIF位置向量的某一元素α_pi与某一非参考道的EPIF位置向量每个元素α_qj的比值以确定w_qj；

$w_{\mathrm{qj}}=\frac{{\mathrm{\α}}_{\mathrm{pi}}}{{\mathrm{\α}}_{\mathrm{qj}}}---\left(15b\right)$

若w_qj达到给定的阈值，则w_qj置1，否则置0；求出所有w_qj后，根据式(15c)确定w_pi；

$w_{\mathrm{pi}}=\frac{\underset{j=1}{\overset{N_0}{\mathrm{\Σ}}}{w}_{\mathrm{qj}}{\mathrm{\α}}_{\mathrm{qj}}}{{\mathrm{\α}}_{\mathrm{pi}}}---\left(15c\right)$

遍历参考道EPIF位置向量的所有元素，最终确定w_p；参考道不变，将上述过程遍历其他非参考道，每次将w_p累加且将w_q清零；

(3)依次从近偏移距的N道中选取其他道做为参考道，重复第（2）步的操作，并将每次求得的权重w_p累加；若某位置w_p的值达到给定的阈值，则认为该位置是同相轴位置；

（4）用同样的方法求取由选取的近偏移距地震道确定的参考同相轴位置与先验层位信息提供的同相轴位置的相关系数，确定最终的同相轴位置。

进一步的，以上步骤（2）中，根据以下式(15a)求参考道的EPIF位置向量的某一元素α_pi与某一非参考道的EPIF位置向量每个元素α_qj的比值以确定w_qj：

$E\left(w_p\right)=\underset{p=1,...,5,q\≠p}{\min }\underset{i=1}{\overset{N_1}{\mathrm{\Σ}}}{(w_{\mathrm{pi}}{\mathrm{\α}}_{\mathrm{pi}}-\underset{j=1}{\overset{N_0}{\mathrm{\Σ}}}{w}_{\mathrm{qj}}{\mathrm{\α}}_{\mathrm{qj}})}^2---\left(15a\right);$

式中w_p及w_q为相邻道层位位置相关性的权重；α_q为除参考道外其余几道提供的同相轴位置，α_p、α_q、w_p、w_q同为列向量，下标i、j分别为列向量的元素，N₀、N₁为向量元素个数。

本发明具有以下有益效果：

本发明采用具有4个待定系数的函数去逼近震源子波，利用粘弹介质中单程波传播理论推导出地震子波包络峰值处的瞬时频率(EPIF)与不同偏移距处旅行时的关系；利用该关系，外推出该同相轴零偏移距处的EPIF，用小波域包络峰值处瞬时频率法(WEPIF)结合层位信息估计Q值。其发展了适合常相位子波的EPIFVO(子波包络峰值处瞬时频率随偏移距变化)法;使估计Q值更精确，更有利于含气性预测。

附图说明

图1(a)最短时间路径;(b)零偏移道的波传播示意图;

图2入射P波、反射P波、透射P波的关系及CMP模型示意图;

图3(a)多层介质;(b)单层介质;

图4单个倾斜界面示意图;

图5多个倾斜界面示意图;

图6EPIFVO示意图(a)零偏EPIF(b)CMP道集的瞬时频率IF剖面;

图7惠更斯原理示意图;

图8(a)模型的参数和观测系统;(b)正演的叠前CMP资料;

图9子波EPIF与偏移距的剖面;

图10各同相轴的子波EPIF与时间的剖面;

图11以相邻两层直接求取方式用WEPIF法估计的衰减曲线;(b)用PFS方法估计的衰减曲线;

图12以逐层递推的方式用WEPIF法估计的衰减;

图13用线性回归得到的斜率估计的衰减;

图14人机交互误差产生示意图;

图15(a)反射系数序列模型;(b)射线追踪法合成的叠前CMP资料;

图16(a)未处理异常的衰减曲线;(b)处理异常后估计的衰减曲线;

图17单个倾斜层状介质合成记录;

图18单个倾斜层叠前CMP合成资料的衰减估计曲线。

具体实施方式

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

1）EPIF与偏移距的关系(EPIFVO)

处理叠前资料的一个难点是波在各层间传播的旅行时难提取。图1是水平层状介质中波传播示意图。图1a中，检波器R₁处接收到的波传播旅行时与入射角θ₁有关，检波器R₂处接收到的波的旅行时与入射角θ₁及透射角θ₂有关，以此类推，多层介质的波传播旅行时与各层入射角及透射角有关，即与介质的速度、密度等有关，但实际中获取精确的角度信息及介质参数比较困难。而零偏移距道的波是垂直入射垂直反射的，如图1b所示，各层的入射角和透射角都已知，故其旅行时可以通过测井曲线或叠后资料等途径得到，因此我们试图寻找不同偏移距处的EPIF与偏移距或旅行时的关系，以期利用该关系外推出零偏移距道的EPIF，进而利用零偏移距道资料求取Q值，这样就可以避免求取不同偏移距处的层间旅行时。

下面推导EPIE与旅行时的关系。

1.1单层介质的情况

假设震源子波可以用如下具有4个参数的常相位子波近似：

$s(0,t)=A{\left(\frac{{\mathrm{\δ}}^2}{\mathrm{\π}}\right)}^{1/4}\exp [i(\mathrm{\σt}+\mathrm{\φ})-{\left(\mathrm{\δt}\right)}^2/2]---\left(1\right)$

其中，σ为调制角频率，δ为能量衰减因子，A和φ分别为振幅和相位。

对(1)式两边做Fourier变换得：

在水平层状粘弹介质中，假设Q值与频率无关，只考虑单程波传播，则由衰减引起的振幅变化的表达式为：

$A\∝A_0{e}^{-\mathrm{\ω}{t}^{*}/2}---\left(3\right)$

其中，A₀及A分别为射线起始点和终止点的振幅谱。

地震子波以入射角θ传播到厚度为Δz的界面又反射回地面(如图2所示)后，将(2)式代入(3)式可得子波的反射波频域表达式，即

$S(\mathrm{\Δz},\mathrm{\θ},\mathrm{\ω})=G\left(\mathrm{\Δz}\right)R\left(\mathrm{\θ}\right)S(0,\mathrm{\ω})\exp [-\frac{i2\mathrm{\ω\Δz}\·\sec \left(\mathrm{\θ}\right)}{c\left(\mathrm{\ω}\right)}-\frac{\mathrm{\ω\Δz}\·\sec \left(\mathrm{\θ}\right)}{c\left(\mathrm{\ω}\right)Q}]---\left(4\right)$

其中ω为角频率，Δz为地层厚度，G(Δz)和R(θ)分别为不依赖于频率和吸收的因子（如几何扩散等）及界面反射系数，c(ω)为相速度，S(0,ω)为初始时刻的震源子波频域表达式，S(Δz,θ,ω)为检波器接收到的反射波的频域表达式。

Barnes指出，常相位子波的包络峰值瞬时频率(EPIF)等于子波振幅谱对Fourier频率加权的平均频率，即

$f_p\left(\mathrm{\τ}\right)=\frac{{\∫}_0^{\∞}\mathrm{fA}(\mathrm{\τ},f)\mathrm{df}}{{\∫}_0^{\∞}A(\mathrm{\τ},f)\mathrm{df}}---\left(5\right)$

其中，A(τ,f)为子波的傅里叶振幅谱，即Fourier谱的模。

将式(2)的模代入式(5)得到零时刻地震子波的EPIF：

$f_p\left(0\right)=\frac{\mathrm{\σ}}{2\mathrm{\π}}+\frac{\frac{{\mathrm{\δ}}^2}{2{\mathrm{\π}}^2}\exp [-\frac{2{\mathrm{\π}}^2}{{\mathrm{\δ}}^2}{\left(\frac{\mathrm{\σ}}{2\mathrm{\π}}\right)}^2]}{{\∫}_0^{\∞}\exp [-\frac{2{\mathrm{\π}}^2}{{\mathrm{\δ}}^2}{(f-\frac{\mathrm{\σ}}{2\mathrm{\π}})}^2]\mathrm{df}}---\left(6\right)$

若忽略速度频散，即c(ω)＝c，则将（2）式代入(4)式可得子波传播时间τ后的振幅谱：

$S(\mathrm{\τ},\mathrm{\θ},f)=\mathrm{AG}\left(\mathrm{\Δz}\right)R\left(\mathrm{\θ}\right){\left(\frac{4\mathrm{\π}}{{\mathrm{\δ}}^2}\right)}^{1/4}\exp (\frac{{\mathrm{\δ}}^2{\mathrm{\τ}}^2}{8Q^2}-\frac{\mathrm{\σ\τ}}{2Q})\exp [-\frac{2{\mathrm{\π}}^2}{{\mathrm{\δ}}^2}{(f-\frac{\mathrm{\σ}}{2\mathrm{\π}}+\frac{{\mathrm{\τ\δ}}^2}{4\mathrm{\πQ}})}^2---\left(7\right)$

其中， $\mathrm{\τ}=\frac{2\mathrm{\Δz}\·\sec \left(\mathrm{\θ}\right)}{c\left(\mathrm{\ω}\right)}$ 为信号传播时间。

将(7)式代入（5）式得到传播时间τ之后子波的EPIF如下：

$f_p\left(\mathrm{\τ}\right)=\frac{{\∫}_0^{\∞}\mathrm{fexp}[-\frac{2{\mathrm{\π}}^2}{{\mathrm{\δ}}^2}{(f-\frac{\mathrm{\σ}}{2\mathrm{\π}}+\frac{{\mathrm{\τ\δ}}^2}{4\mathrm{\πQ}})}^2]\mathrm{df}}{{\∫}_0^{\∞}\exp [-\frac{2{\mathrm{\π}}^2}{{\mathrm{\δ}}^2}{(f-\frac{\mathrm{\σ}}{2\mathrm{\π}}+\frac{{\mathrm{\τ\δ}}^2}{4\mathrm{\πQ}})}^2]\mathrm{df}}$

（8）

$=(\frac{\mathrm{\σ}}{2\mathrm{\π}}-\frac{{\mathrm{\τ\δ}}^2}{4\mathrm{\πQ}})+\frac{\frac{{\mathrm{\δ}}^2}{4{\mathrm{\π}}^2}\exp [-\frac{2{\mathrm{\π}}^2}{{\mathrm{\δ}}^2}{(\frac{\mathrm{\σ}}{2\mathrm{\π}}-\frac{{\mathrm{\τ\δ}}^2}{4\mathrm{\πQ}})}^2]}{{\∫}_0^{\∞}\exp [-\frac{2{\mathrm{\π}}^2}{{\mathrm{\δ}}^2}{(f-\frac{\mathrm{\σ}}{2\mathrm{\π}}+\frac{{\mathrm{\τ\δ}}^2}{4\mathrm{\πQ}})}^2]\mathrm{df}}$

由上式可以看出，影响子波振幅谱的振幅因子A、不依赖于频率的因子G(Δz)及反射系数R(θ)，在求取包络峰值瞬时频率时由于比值关系而相互抵消，因此在层状均匀粘弹介质中，不考虑多次波的情况下，利用EPIF可以消除界面反射系数、振幅等因素的影响。

将式（8）经过一系列的展开、变量替换及舍去高阶项后得到：

$f_p\left(\mathrm{\τ}\right)\stackrel{\·}{=}{f}_p\left(0\right)-(\frac{{\mathrm{\δ}}^2}{4\mathrm{\πQ}}-\frac{\mathrm{D\σ}}{2Q})\mathrm{\τ}\stackrel{\·}{=}{f}_p\left(0\right)-\frac{{\mathrm{\δ}}^{2}k\left(\mathrm{\η}\right)}{4\mathrm{\πQ}}\mathrm{\τ}---\left(9\right)$

其中， $D=\frac{1}{I_1}\frac{{\mathrm{\δ}}^2}{4{\mathrm{\π}}^2}\exp (-2{\mathrm{\π}}^2{\mathrm{\η}}^2),$ $\mathrm{\η}=\frac{\mathrm{\σ}}{2\mathrm{\π\δ}},$ $I_1={\∫}_{-\frac{\mathrm{\σ}}{2\mathrm{\π}}}^{\∞}\exp (-\frac{2{\mathrm{\π}}^2}{{\mathrm{\δ}}^2}{f}^2)\mathrm{df},$

$\mathrm{\κ}\left(\mathrm{\η}\right)=1-\sqrt{2\mathrm{\π}}\mathrm{\η}{\mathrm{\Φ}}^{-1}\left(2\mathrm{\π\η}\right)\exp (-2{\mathrm{\π}}^2{\mathrm{\η}}^2),$ Φ^-1(*)为标准正态分布概率积分函数。

由(9)式可见，传播一定时间后子波的EPIF与旅行时近似为线性关系，其截距是初始时刻子波的EPIF，斜率与Q有关。

1.2多层介质的情况

假设介质有N个反射界面，检波器得到的经第N个界面反射的反射波为：

$S(\mathrm{\Δz},\mathrm{\θ},\mathrm{\ω})=\left[\underset{i=1}{\overset{N}{\mathrm{\Π}}}G\left(\mathrm{\Δ}{z}_i\right)R\left({\mathrm{\θ}}_i\right)\right]S(0,\mathrm{\ω})\exp [-i2\mathrm{\ω}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\left(\frac{\mathrm{\Δ}{z}_i\·\sec \left({\mathrm{\θ}}_i\right)}{c_i\left(\mathrm{\ω}\right)}\right)-\mathrm{\ω}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\left(\frac{\mathrm{\Δ}{z}_i\·\sec \left({\mathrm{\θ}}_i\right)}{c_i\left(\mathrm{\ω}\right)Q_i}\right)]---\left(10\right)$

若忽略速度频散，令 $\mathrm{\Δ}{t}_i=\frac{2\mathrm{\Δ}{z}_i\·\sec \left({\mathrm{\θ}}_i\right)}{c_i\left(\mathrm{\ω}\right)},$ 则反射波的振幅谱为：

$S(\mathrm{\Δt},f)=\left[\underset{i=1}{\overset{N}{\mathrm{\Π}}}G\left(\mathrm{\Δ}{z}_i\right)R\left({\mathrm{\θ}}_i\right)\right]S(0,\mathrm{\ω})\exp [-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\left(\frac{\mathrm{\πf\Δ}{t}_i}{Q_i}\right)]---\left(11\right)$

(11)式代入(5)式得到EPIF的表达式后，将其展开、变量替换及舍去高阶项后得到：

$f_p\left(\mathrm{\τ}\right)=f_p\left(0\right)-\frac{{\mathrm{\δ}}^{2}k\left(\mathrm{\η}\right)}{4\mathrm{\π}}\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\mathrm{\Δ}{t}_i}{Q_i}\right)=f_p\left(0\right)-\frac{{\mathrm{\δ}}^{2}k\left(\mathrm{\η}\right)}{4\mathrm{\π}}(\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}\frac{\mathrm{\Δ}{t}_i}{Q_i}-\frac{\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}\mathrm{\Δ}{t}_i}{Q_N})-\frac{{\mathrm{\δ}}^{2}k\left(\mathrm{\η}\right)}{4\mathrm{\π}{Q}_N}\mathrm{\τ}---\left(12\right)$

其中Δt_i是反射波在第i层的双程旅行时，

从(12)式可以看出，N层介质时子波的EPIF与旅行时仍然满足线性关系，但它的截距除了与零时刻震源子波的EPIF有关外，还与各层介质的旅行时及品质因子有关，而旅行时Δt_i因射线入射角不同而不同，故其截距不是常数，因此不能用此公式直接求衰减。

地震波传播时真正遵循的是“沿最小时间路程传播”原理（费马原理），如图3a所示。在小偏移距小排列的情况下，将多层介质的反射波时距曲线近似地看成双曲线，在一定精度要求下是可以的。讨论层状介质问题的基本思路是：如图3a所示的水平层状介质，我们可以把R₂界面以上的介质设法用等效均匀介质来代替，并令这种假想的均匀介质的波速度及品质因子取某个等效值，使得R₂界面以上的介质简化为均匀介质，即变成单层模型，同样可以把R_N界面以上的N层介质用具有某个等效品质因子的均匀介质来代替，如图3b所示。

对于（12）式中左端的两项，令（12）式简化为：

$f_p\left(\mathrm{\τ}\right)\text{=}{f}_p\left(0\right)-\frac{{\mathrm{\δ}}^{2}k\left(\mathrm{\η}\right)}{4\mathrm{\π}}\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\mathrm{\Δ}{t}_i}{Q_i}\right)=f_p\left(0\right)-\frac{{\mathrm{\δ}}^{2}k\left(\mathrm{\η}\right)}{4\mathrm{\π}}\frac{\mathrm{\τ}}{Q_{\mathrm{eff}}}---\left(13\right)$

上式与（9）式形式相同，不同点是(12)式中的Q值为等效Q值，而(9)式中的Q值为层间Q值。

下面分析(13)式的意义及作用。(13)式表示的是地表上震源与接收点的EPIF与旅行时的关系。沿每个同相轴拾取不同偏移距处的EPIF和旅行时并进行线性回归，计算斜率和截距，进而外推出零时刻子波的EPIF，这就是EPIFVO（包络峰值瞬时频率随偏移距的变化）的定义。由（13）式看出，斜率应该是负的，若拟合出的斜率为正，则认为此处频率的变化不是由衰减引起的，可能是由薄互层导致的同相轴调谐，故需先处理异常值。

前面讨论单层及多层介质中EPIF与旅行时满足的关系时，波是从零时刻开始传播的，为了求取Q值方便，下面将其推广到以任意时刻为起始点时的情况。

设震源子波的振幅谱为S(0,f),那么零偏移距处(M点)接收到的子波振幅谱为：

$S(t_0,f)=G_{1}R\left(\frac{\mathrm{\π}}{2}\right)S(0,f)\exp [-\frac{\mathrm{\π}{\mathrm{ft}}_0}{Q_{\mathrm{eff}}}]---\left(14\right)$

而偏移距x处(R点)接收到的子波振幅谱为：

$S(\mathrm{\τ},f)=G_{2}R\left(\mathrm{\θ}\right)S(0,f)\exp [-\frac{\mathrm{\πf\τ}}{Q_{\mathrm{eff}}}]=G_{2}R\left(\mathrm{\θ}\right)S(0,f)\exp [-\frac{\mathrm{\π}{\mathrm{ft}}_0}{Q_{\mathrm{eff}}}]\exp [-\frac{\mathrm{\πf\Δt}}{Q_{\mathrm{eff}}}]---\left(15\right)$

其中，τ与t₀满足τ＝t₀+Δt。

将(14)式代入(15)式，可得：

$S(\mathrm{\τ},f)=\mathrm{GRS}(t_0,f)\exp (-\frac{\mathrm{\πf\Δt}}{Q_{\mathrm{eff}}})---\left(16\right)$

其中G＝G₂/G₁,R＝R(θ)/R(π/2)。

经过与式(9)的推导过程类似的推导，得到初始时刻为t₀时τ时刻子波的EPIF为：

$f_p\left(\mathrm{\τ}\right)=f_p(t_0+\mathrm{\Δt})=f_p\left(t_0\right)-\frac{{\mathrm{\δ}}^{2}k\left(\mathrm{\η}\right)}{4\mathrm{\π}{Q}_{\mathrm{eff}}}\mathrm{\Δt}---\left(17\right)$

其中，f_p(t₀)为中点处接收到的子波的包络峰值瞬时频率。由式（17）可知，若线性拟合时采用正常时差Δt（即垂直时差），则拟合出的直线的截距为该同相轴零偏移距处的EPIF。

1.3含倾斜反射界面时的情况

界面倾斜时，中心点在界面的投影与真正的反射点有偏移。单个倾斜层的情况下(如图4所示)，假设每个炮点的震源子波相同，从叠前CMP资料的时距曲线上拾取接收信号的旅行时τ，若先验已知中心点M处的自激自收旅行时t₀(射线垂直于界面入射)，可以证明，不管是沿上倾方向还是下倾方向，中心点M处的旅行时t₀是最小的。证明如下：

设H为M点处界面的法线深度，x为炮检距，根据几何关系且应用余弦定理，则偏移距x处射线的路程平方为：

中心点M处的垂直路程平方为I₂＝(2H)²，路径差速度为常数。这样就可以式(17)求出单个倾斜层零偏移距处的EPIF。

由于速度均匀，因此对水平层状介质的EPIFVO分析同样适用于单个倾斜层的情况。根据(9)式及(17)式可以用线性回归的方法求出零偏移距处的EPIF，若线性回归时所用的时间为接收信号的旅行时τ，则拟合得到的截距为震源子波的EPIFf_p(0)，若线性回归时所用的时间为正常时差Δt＝τ-t₀，则拟合得到的截距为中心点M在界面上的投影点R处的EPIFf_p(t₀)，用下节介绍的Q值计算公式可求得Q曲线。

多个倾斜层的情况下，每一层有任意倾斜角，如图5所示。CMP数据的中心点M在目标反射界面上垂直入射到D'，实际路径的反射点为D。Hubral和Krey提出射线沿SDG的旅行时公式为：

其中， $v_{\mathrm{NMO}}^2=\frac{1}{t_0{\cos }^2{\mathrm{\β}}_0}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{v}_i^2\mathrm{\Δ}{t}_i\left(0\right)\underset{k=1}{\overset{i-1}{\mathrm{\Π}}}\left(\frac{{\cos }^2{\mathrm{\α}}_k}{{\cos }^2{\mathrm{\β}}_k}\right),$ t₀为中心点M到D'的垂直双程旅行时，Δt_i(0)为射线从M点到D'点穿过每层时的双程旅行时，当倾斜角很小且小排列（偏移距小于深度）时，旅行时方程可用双曲线方程近似，动校正速度v_NMO可用均方根速度v_rms近似，即式(18)近似为：

$t^2\left(x\right)=t_0^2+\frac{x^2}{v_{\mathrm{rms}}^2}---\left(19\right)$

这种双曲线形式的时距曲线是我们所熟知的，水平层状介质及单个倾斜层模型的时距曲线都是双曲线，因此小倾斜角小排列时，多个倾斜层的介质可以简化为单个倾斜层模型，用EPIFVO方法外推出零偏移距处的EPIF后，由下节介绍的逐层递推法求出Q值。

2)两种求取Q值的方式

2.1逐层递推方式

如图6所示，(a)零偏EPIF(b)CMP道集的瞬时频率IF剖面，在第N个同相轴的零偏移距处，子波的旅行时为τ，拾取的EPIF为f_p(τ)，震源子波的EPIF为f_p(0)，则由式（9）左端前两项可得：

$\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\mathrm{\Δ}{t}_i}{Q_i}=\frac{4\mathrm{\π}}{{\mathrm{\δ}}^{2}k\left(\mathrm{\η}\right)}(f_p\left(0\right)-f_p\left(\mathrm{\τ}\right))---\left(11\right)$

其中，Δf_p＝f_p(0)-f_p(τ)为子波的EPIF之差，Δt_i为零偏移距道第i层的垂直双程旅行时，Q_i为各层间Q值。

将(11)式左边展开移项后得：

$\frac{1}{Q_N}=\frac{1}{\mathrm{\Δ}{t}_N}(\frac{4\mathrm{\π\Δ}{f}_p}{{\mathrm{\δ}}^{2}k\left(\mathrm{\η}\right)}-\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}\frac{\mathrm{\Δ}{t}_i}{Q_i})---\left(12\right)$

(12)式为计算各层衰减的逐层递推公式，其中Δf_p为震源子波的EPIF与检波器接收到的子波的EPIF之差，第N层的衰减估计精度依赖于其上各层的衰减估计精度，因此有误差累积效应。

2.2相邻两层直接求取方式

我们将（9）式做如下变形：

$f_p\left(t_j\right)=f_p\left(0\right)-\frac{{\mathrm{\δ}}^{2}k\left(\mathrm{\η}\right)}{4\mathrm{\π}}(\underset{i=1}{\overset{j-1}{\mathrm{\Σ}}}\frac{\mathrm{\Δ}{t}_i}{Q_i}+\frac{\mathrm{\Δ}{t}_j}{Q_j})=f_p\left(t_{j-1}\right)-\frac{{\mathrm{\δ}}^{2}k\left(\mathrm{\η}\right)}{4\mathrm{\π}}\frac{\mathrm{\Δ}{t}_j}{Q_j}---\left(13\right)$

进而可得出层间Q值的另一个计算公式：

$\frac{1}{Q_j}=\frac{1}{\mathrm{\Δ}{t}_j}\frac{4\mathrm{\π}(f_p\left(t_{j-1}\right)-f_p\left(t_j\right))}{{\mathrm{\δ}}^{2}k\left(\mathrm{\η}\right)}---\left(14\right)$

其中下标j表示第j层介质。式（14）的物理意义是：当震源子波传播到第j-1层后，根据惠更斯原理，将该子波波前作为子波源，向各方向继续传播，经过时间Δt_j/2后子波到达第j层并反射，在第j-1层接收到反射波，如图7所示。f_p(t_j-1)为第j-1层子波源的EPIF，f_p(t_j)为接收到的反射波的EPIF。利用(14)式就可以求出第j层的层间Q值。这样利用零偏移距道相邻两层之间子波的EPIF变化来估计Q值可以避免误差累积效应。

层位信息拾取

处理实际资料时，由叠后数据或测井资料拾取到的层位信息与叠前CMP资料的同相轴位置不是完全吻合，因此本发明不直接用反射系数序列来引导同相轴位置的拾取，而是利用叠前CMP道集中道与道之间的相关性并结合先验的层位信息来拾取。在叠前CMP道集中，每一道的有效信号与邻近道的有效信号具有相关性，而噪声无相关性，故我们用近偏移距处相邻的几道地震记录来确定同相轴位置。目标函数为：

$E\left(w_p\right)=\underset{p=1,...,5,q\≠p}{\min }\underset{i=1}{\overset{N_1}{\mathrm{\Σ}}}{(w_{\mathrm{pi}}{\mathrm{\α}}_{\mathrm{pi}}-\underset{j=1}{\overset{N_0}{\mathrm{\Σ}}}{w}_{\mathrm{qj}}{\mathrm{\α}}_{\mathrm{qj}})}^2---\left(15a\right)$

$w_{\mathrm{qj}}=\frac{{\mathrm{\α}}_{\mathrm{pi}}}{{\mathrm{\α}}_{\mathrm{qj}}}---\left(15b\right)$

$w_{\mathrm{pi}}=\frac{\underset{j=1}{\overset{N_0}{\mathrm{\Σ}}}{w}_{\mathrm{qj}}{\mathrm{\α}}_{\mathrm{qj}}}{{\mathrm{\α}}_{\mathrm{pi}}}---\left(15c\right)$

其中p为选取的参考道的数目，α_p为由近偏移距的几道选出的参考道确定的同相轴位置，w_p及w_q为相邻道层位位置相关性的权重，α_q为除参考道外其余几道提供的同相轴位置，α_p、α_q、w_p、w_q同为列向量，下标i、j分别为列向量的元素，N₀、N₁为向量元素个数。

层位拾取的实现步骤如下：

(1)选取近偏移距中的几道，如选取叠前CMP资料的第1～5道（即p=5），计算这五道的瞬时振幅IA和瞬时频率IF，由IA和IF拾取EPIF及其位置α，我们认为EPIF的位置等效为同相轴的位置，EPIF的位置为一个列向量；

(2)从近偏移距的几道中选取参考道，如选取第一道作为参考，根据式(15b)，求参考道的EPIF位置向量的某一元素α_pi与某一非参考道的EPIF位置向量每个元素α_qj的比值以确定w_qj。若w_qj达到给定的阈值，则w_qj置1，否则置0；求出所有w_qj后，根据式(15c)确定w_pi；遍历参考道EPIF位置向量的所有元素，最终确定w_p；参考道不变，将上述过程遍历其他非参考道，每次将w_p累加且将w_q清零；

(3)依次从近偏移距的几道中选取其他道做为参考道，重复第2步的操作，并将每次求得的权重w_p累加；若某位置w_p的值达到给定的阈值，则认为该位置是同相轴位置。该步骤可以避免参考道的某层反射受干扰的情况；

用同样的方法求取由选取的近偏移距地震道确定的参考同相轴位置与先验层位信息提供的同相轴位置的相关系数，确定最终的同相轴位置。

通过以上公式的推导及水平层状介质和倾斜层状介质中EPIFVO(包络峰值处瞬时频率随偏移距变化)方法的讨论可知，实现该方法首先要确定同相轴的位置，计算各非零偏道的EPIF，外推求出零偏移距处的EPIF，求取子波参数然后利用WEPIF方法求Q值。本发明的具体实现方法如下：

1)层位信息拾取

处理实际资料时，由叠后数据或测井资料拾取到的层位信息与叠前CMP资料的同相轴位置不是完全吻合，因此本发明不直接用反射系数序列来引导同相轴位置的拾取，而是利用叠前CMP道集中道与道之间的相关性并结合先验的层位信息来拾取。在叠前CMP道集中，每一道的有效信号与邻近道的有效信号具有相关性，而噪声无相关性，故我们用近偏移距处相邻的几道地震记录来确定同相轴位置。目标函数为：

$E\left(w_p\right)=\underset{p=1,...,5,q\≠p}{\min }\underset{i=1}{\overset{N_1}{\mathrm{\Σ}}}{(w_{\mathrm{pi}}{\mathrm{\α}}_{\mathrm{pi}}-\underset{j=1}{\overset{N_0}{\mathrm{\Σ}}}{w}_{\mathrm{qj}}{\mathrm{\α}}_{\mathrm{qj}})}^2---\left(25a\right)$

$w_{\mathrm{qj}}=\frac{{\mathrm{\α}}_{\mathrm{pi}}}{{\mathrm{\α}}_{\mathrm{qj}}}---\left(25b\right)$

$w_{\mathrm{pi}}=\frac{\underset{j=1}{\overset{N_0}{\mathrm{\Σ}}}{w}_{\mathrm{qj}}{\mathrm{\α}}_{\mathrm{qj}}}{{\mathrm{\α}}_{\mathrm{pi}}}---\left(25c\right)$

其中p为选取的参考道的数目，α_p为由近偏移距的几道选出的参考道确定的同相轴位置，w_p及w_q为相邻道层位位置相关性的权重，α_q为除参考道外其余几道提供的同相轴位置，α_p、α_q、w_p、w_q同为列向量，下标i、j分别为列向量的元素，N₀、N₁为向量元素个数。

层位拾取的实现步骤如下：

(1)选取近偏移距中的几道，如选取叠前CMP资料的第1～5道（即p=5），计算这五道的瞬时振幅IA和瞬时频率IF，由IA和IF拾取EPIF及其位置α，认为EPIF的位置等效为同相轴的位置，EPIF的位置为一个列向量；

(2)从近偏移距的几道中选取参考道，如选取第一道作为参考，根据式(25b)，求参考道的EPIF位置向量的某一元素α_pi与某一非参考道的EPIF位置向量每个元素α_qj的比值以确定w_qj。若w_qj达到给定的阈值，则w_qj置1，否则置0；求出所有w_qj后，根据式(25c)确定w_pi；遍历参考道EPIF位置向量的所有元素，最终确定w_p；参考道不变，将上述过程遍历其他非参考道，每次将w_p累加且将w_q清零；

(3)依次从近偏移距的几道中选取其他道做为参考道，重复第2步的操作，并将每次求得的权重w_p累加；若某位置w_p的值达到给定的阈值，则认为该位置是同相轴位置。该步骤可以避免参考道的某层反射受干扰的情况；

(4)用同样的方法求取由选取的近偏移距地震道确定的参考同相轴位置与先验层位信息提供的同相轴位置的相关系数，确定最终的同相轴位置。

2)计算瞬时频率及估计子波参数

瞬时频率的计算可以直接采用瞬时频率的定义，即瞬时相位的导数，对于含噪信号，通常采用阻尼瞬时频率以增加计算的稳定性，即

$f\left(t\right)=\left(\frac{1}{2\mathrm{\π}}\right)\frac{s\left(t\right)\frac{d{s}^{*}\left(t\right)}{\mathrm{dt}}-s^{*}\left(t\right)\frac{\mathrm{ds}\left(t\right)}{\mathrm{dt}}}{a^2\left(t\right)+{\mathrm{\ϵ}}^2}---\left(26\right)$

其中s^*(t)是信号s(t)的希尔伯特变换，ds^*(t)是s^*(t)的导数，a(t)是瞬时振幅。

为了进一步稳定瞬时频率，通常还用地震信号幅度的平方对阻尼瞬时频率加权，即

$f_w\left(t\right)=\frac{{\∫}_{t-L}^{t+L}f\left(t^{\′}\right)W\left(t^{\′}\right)d{t}^{\′}}{{\∫}_{t-L}^{t+L}W\left(t^{\′}\right)d{t}^{\′}}---\left(27\right)$

其中，W(t')为信号幅度的平方，L为加权窗的长度。

利用WEPIF法估计Q曲线时，需要知道子波调制频率σ和能量衰减因子δ的值。可用如下公式近似计算：

$\mathrm{\σ}=\frac{{\∫}_0^{\∞}\mathrm{\ω}\left|\hat{s}\left(\mathrm{\ω}\right)\right|\mathrm{d\ω}}{{\∫}_0^{\∞}\left|\hat{s}\left(\mathrm{\ω}\right)\right|\mathrm{d\ω}}---\left(28\right)$

${\mathrm{\δ}}^2=\frac{{\∫}_0^{\∞}{(\mathrm{\ω}-\mathrm{\σ})}^2\left|\hat{s}\left(\mathrm{\ω}\right)\right|\mathrm{d\ω}}{{\∫}_0^{\∞}\left|\hat{s}\left(\mathrm{\ω}\right)\right|\mathrm{d\ω}}---\left(29\right)$

其中为子波s(t)的频谱。

3)EPIFVO计算步骤

用WEPIF方法估计叠前CMP道集的Q值前，先进行预处理，去除不同炮子波的不一致性。EPIFVO的实现步骤如下：

（1）、构造一个CMP超道集；

叠前数据信噪比低，因此在每道数据的邻域内选择极少的几道构造一个CMP超道集，对子波特征破坏不大，又可起到压制部分随机噪声的作用。

（2）、每道数据用三参数小波做小波变换，在小波域计算瞬时振幅IA和瞬时频率IF；

（3）、根据道与道的相关性拾取叠前CMP数据的同相轴，即确定层位信息；

（4）、根据层位信息，在反射界面附近拾取各道IA的包络峰值对应的瞬时频率，即包络峰值处瞬时频率EPIF，由于水平层状介质中各同相轴的时距曲线是双曲线，故同相轴位置满足如下约束条件：

其中，pos|_now(event)为当前道某同相轴的位置，pos|_fore(event)为前一道同一个同相轴的位置，即当前道的同相轴位置必然大于前一道同一同相轴的位置，而有可能小于前一道下一同相轴的位置（小排列时成立，而同相轴有交叉时可能不成立），实际中利用同相轴的梯度信息拾取连续同相轴更准确；

（5）、单个同相轴计算，从叠前剖面上拾取波至时间，用不同偏移距处的EPIF和旅行时差(即正常时差)线性回归计算斜率和截距（截距即零偏移距道的EPIF），并处理异常斜率对应的截距；

（6）、用WEPIF法计算Q曲线。

数值仿真结果

1)合成叠前CMP资料

首先用一个合成的叠前CMP资料来检验文中所提方法的有效性。图8a为模型的参数和观测系统。正演叠前CMP资料时，震源子波采用主频为50Hz的常相位子波，最小炮检距为10m，用49个间距为5m的检波器进行接收，共49道，采样率为1ms。叠前CMP资料的正演采用基于Fermat原理的射线追踪法，正演的数据如图8b所示。合成资料共5个同相轴，各层理论的层间Q值分别为150、200、100、150、250。图9为5个同相轴的子波EPIF与偏移距的剖面，可以看出其为双曲线；图10为5个同相轴的子波EPIF与时间的剖面，周围的5个图分别为5个同相轴的EPIF与时间关系的放大图，可以看出5个同相轴都近似满足线性关系，验证了理论的正确性；

本发明对比了WEPIF法和Zhang文章中所用的PFS法估计Q值的精度。WEPIF法和PFS法对子波有一定的依赖性，WEPIF法对子波的依赖性要小于PFS法。合成记录的子波参数η＝0.467。基于形如公式(1)的震源子波推导出PFS法的Q值计算公式为：

$Q=-\frac{\mathrm{\Δt}{\mathrm{\δ}}^2}{2(2\mathrm{\π\Δ}{f}_p-\mathrm{\σ})}$

其中Δf_p为峰值频率的差，σ,δ为子波参数。图11(a)为根据公式(24)，利用相邻两层的EPIF，用WEPIF法估计的衰减曲线，其中点线是真实的Q值，实线是估计的Q值，估计的误差均在合理范围内，验证了该方法的有效性；图11(b)为PFS法估计的衰减曲线，其中红色虚线是真值，蓝色实线是估计的Q值，从图中可见误差较大。

图12为用逐层递推公式(22)估计的衰减曲线，可以看出深层的误差较大，误差累积效应较明显；

由式(13)可知，线性回归的斜率为与等效衰减Q_eff有关，因此我们可以用斜率直接求取等效Q_eff值，即

$\frac{1}{Q_{\mathrm{eff}}}=-\frac{4\mathrm{\πGf}}{{\mathrm{\δ}}^{2}k\left(\mathrm{\η}\right)}---\left(30\right),$

进而用公式求取各层间Q值。结果如图13所示。

由图13可看出，直接用斜率估计衰减的误差波动较大，估计结果不稳定，原因是线性回归时需要进行人机交互修正部分EPIF值，以图14a所示的序列为例，如图14b所示，只上下移动第一个和第三个点，线性回归得到的斜率和截距差距很大。

式(30)中Q_eff与斜率成正比例关系，故误差较大，且求取层间Q值时有误差累计效应，而式(22)及式(24)中Q与Δf_p(EPIF之差)有关，若人机交互采用的准则相同，两个包络峰值瞬时频率(EPIF)相减可以消去部分人机交互带来的主观影响，且可以直接求取层间Q值。

2)合成含有薄层调谐效应的CMP数据

用射线追踪法合成一个含有薄互层的模型，其中第四层为薄互层，最后一层的Q为50，其上各层的Q为500，其余参数与图8a所示相同。图15a为反射系数序列模型(本图指示反射系数的时间位置，长短不代表反射系数的数值)，图15b为合成的叠前CMP数据，第四层的上下两个同相轴由于调谐使得分辨率降低，波形干扰使得两个同相轴相互耦合成波形展宽的一个同相轴；

5个同相轴线性回归拟合的斜率和截距如表1所示，第3个同相轴的斜率为正，其余同相轴的斜率为负。2.1节指出，拟合的斜率应该是负的，如果为正，说明频率变化不是由衰减引起的，故需处理异常。表1的结果验证了理论的正确性。本发明处理异常的方法是，若斜率为正，则将对应的截距取其相邻两个截距的平均。

图16a为未处理异常时估计的衰减曲线，误差较大，图16b为处理异常后估计的衰减曲线，其误差在合理范围内。

3)合成单个倾斜层状介质的叠前CMP资料

用射线追踪法合成单个倾斜界面的叠前CMP资料，如图17所示:

倾斜角度为3°，理论Q值为75，估计的衰减值为79.7872，误差在合理的范围内，如图18所示。

Claims (2)

1.一种估计介质品质因子的EPIFVO方法，其特征在于：首先进行层位信息拾取，包括以下步骤：

(1)选取近偏移距中的N道，N大于1，计算这N道的瞬时振幅IA和瞬时频率IF，由IA和IF拾取EPIF及其位置α，将EPIF的位置等效为同相轴的位置，EPIF的位置为一个列向量；

(2)从近偏移距的N道中选取参考道，根据式(15b)，求参考道的EPIF位置向量的某一元素α_pi与某一非参考道的EPIF位置向量每个元素α_qj的比值以确定w_qj；

$w_{qj}=\frac{{\mathrm{\α}}_{pi}}{{\mathrm{\α}}_{qj}}---\left(15b\right)$

若w_qj达到给定的阈值，则w_qj置1，否则置0；求出所有w_qj后，根据式(15c)确定w_pi；

$w_{pi}=\frac{\underset{j=1}{\overset{N_0}{\Σ}}{w}_{qj}{\mathrm{\α}}_{qj}}{{\mathrm{\α}}_{pi}}---\left(15c\right)$

遍历参考道EPIF位置向量的所有元素，最终确定w_p；参考道不变，将上述过程遍历其他非参考道，每次将w_p累加且将w_q清零；

(3)依次从近偏移距的N道中选取其他道做为参考道，重复第(2)步的操作，并将每次求得的权重w_p累加；若某位置w_p的值达到给定的阈值，则认为该位置是同相轴位置；

(4)用同样的方法求取由选取的近偏移距地震道确定的参考同相轴位置与先验层位信息提供的同相轴位置的相关系数，确定最终的同相轴位置；

当震源子波传播到第j-1层后，根据惠更斯原理，将该子波波前作为子波源，向各方向继续传播，经过时间Δt_j/2后子波到达第j层并反射，在第j-1层接收到反射波；f_p(t_j-1)为第j-1层子波源的EPIF，f_p(t_j)为接收到的反射波的EPIF；利用(14)式求出第j层的层间Q值：

$\frac{1}{Q_j}=\frac{1}{{\mathrm{\Δt}}_j}\frac{4\mathrm{\π}\left(f_p\right(t_{j-1})-f_p(t_j\left)\right)}{{\mathrm{\δ}}^{2}k\left(\mathrm{\η}\right)}---\left(14\right);$

其中下标j表示第j层介质。

2.根据权利要求1所述的估计介质品质因子的EPIFVO方法，其特征在于：步骤(2)中，根据以下式(15a)求参考道的EPIF位置向量的某一元素α_pi与某一非参考道的EPIF位置向量每个元素α_qj的比值以确定w_qj：

$E\left(w_p\right)=\underset{p=1,...,5,q\≠p}{min}\underset{i=1}{\overset{N_1}{\Σ}}{(w_{pi}{\mathrm{\α}}_{pi}-\underset{j=1}{\overset{N_0}{\Σ}}{w}_{qj}{\mathrm{\α}}_{qj})}^2---\left(15a\right);$

式中w_p及w_q为相邻道层位位置相关性的权重；α_q为除参考道外其余几道提供的同相轴位置，α_p、α_q、w_p、w_q同为列向量，下标i、j分别为列向量的元素，N₀、N₁为向量元素个数。