基于射线路径的叠前高精度Q值估计方法

王德利，戴建芳

（吉林大学地球探测科学与技术学院，吉林长春130026）

地震波传播过程中存在着能量衰减。引起地震波衰减的因素很多，总的来说可分为两大类［1］：一类与地震波传播特性有关，如球面扩散、界面透射损失、层间多次波、振幅调谐引起的衰减、波的相位转换导致的衰减等；另一类衰减与储层岩石物理参数的变化相关，液体的性质（如粘度、可压缩性）及其饱和度，孔隙的性质（包括孔隙度、渗透率）以及裂隙的几何参数都会对这类衰减产生影响。对于后一类衰减我们通常用品质因子Q来衡量。正确估计Q值对于提高地震记录的分辨率及对储层的预测与描述具有重要意义。

提取Q值的方法有很多种［2］，包括振幅衰减法、子波模拟法、谱模拟法、上升时间法、解析信号法和质心频移法等。Tonn［3］认为这几种方法都有一定的适用性，至今没有任何一种方法是普遍适用的。现在实际应用中最常用的方法是谱比法，这是一种频率域方法，其基本思路是分别截取目标地层顶、底面一个时窗内的地震记录，通过线性拟合频率和振幅之间的关系得出Q值［4］。

许多因素使得衰减测定的准确度降低，这主要源于地震谱中不同波至走时不同，计算出的衰减受相邻层位影响很大。Dasgupta等［5］提出了一种利用振幅值随炮检距的变化从CMP道集中提取Q值的方法；Hackert等［6］对该方法进行了改进，提出了对振幅谱进行校正后提取Q值的方法（LSR）；王小杰等［7］使用两步拟合算法（QVO 方法）通过偏移距的归零处理进一步提高Q值提取的精度。这些方法的有效性已经被证明，包括岩性的识别［8］和裂缝的方位识别［9］。但是，在 QVO方法的实现过程中，时窗内的地震记录会受调谐效应和多次波反射的影响，使得Q值提取精度不高；另外，计算过程中所用到的时差并非仅仅是目标层的旅行时。由于介质的各向异性以及震源与检波器的方向性的影响［10－11］，使得目标层的衰减受相邻层位的影响，所以QVO方法第2步拟合中的归零处理求取的零偏移距上的频谱比斜率并不精确。

考虑到叠前地震资料比叠后地震资料的振幅和旅行时信息更加丰富，更能反映地层的一些细微特征，我们在QVO方法的基础上提出一种基于射线路径的叠前地震数据高精度Q值估计方法，其核心是通过τ－p变换对同一水平慢度的道集进行处理，采用广义S变换计算地层顶、底面反射波的瞬时频谱，从而克服传统Q值估计方法的两步拟合所带来的统计性误差，以提高Q值计算的精度。

1 基本理论

1.1 常规Q值计算方法

振幅谱为S（f）的地震波在层状介质中传播时，不仅有地层的粘弹性所引起的e指数衰减，还有透射和反射的能量损失P和球面扩散能量损失G，而且这两项与频率无关。对于一个初始时刻振幅谱为S0（f）的地震子波，传播t时间后，其振幅谱S（f）表示为

式中：Q是传播过程中的平均衰减效应。

同一道来自同一地层顶、底面反射波的瞬时振幅谱可分别表示为

式中：S1（f）为接收到的地层顶面反射的频谱；S2（f）为接收到的地层底面反射的频谱；t21为接收到的地层底面反射在地层1中的旅行时；t22为其在地层2中的旅行时（图1）。

图1 同一道接收到的地层2顶、底面的反射

（3）式除以（2）式且两边取对数得

通常的处理方法是假设接收到的同一地层顶、底面反射在地层1中的旅行时是一样的，即t21＝t1，从而得到

1.2 高精度Q值计算方法

为了更正确地计算走时，解决地层各向异性和震源、检波器的方向性问题［12］，叠前高精度Q值计算方法是把数据变换到τ－p域，对同一水平慢度的道集进行处理，以避免传统方法的两步拟合所带来的统计性不精确；并且引入走时各向异性参数，使计算出的衰减仅仅由目标层位引起。

在以上常规Q值计算方法中做了一个假设，即t21＝t1，然而由图1可知事实并不是这样的，故计算出的Q值有一定的误差。而在τ－p域中，同一水平慢度的道集接收到的地层2的反射波路径在地层1中的走时是相同的，如图2所示，由此可知t21＝t1。

图2 同一水平慢度道集接收到的地层2的反射波路径

为了把计算过程中相邻层位的干涉效应减到最小，叠前高精度Q值计算方法对振幅谱的计算是利用时频分布图［13］恰当地计算地震波到达地层顶、底面的瞬时频谱。S变换是一种介于短时傅里叶变换和小波变换之间的对非平稳信号处理的时频分析技术［14］，广义S变换是在S变换的基础上发展而来［15］，其窗函数加入了调节系数，是一种分辨率比较高的时频分析方法。利用广义S变换的变化时窗对地层上、下界面对应的瞬时频谱进行分析，能够克服传统Q值计算方法中时窗长度和宽度选择对计算结果的影响。

1.2.1 τ－p 变换

对于同一水平慢度的道集数据，Snell定律保证了在水平各向同性速度和衰减的假设下其在地层1中的反射沿着同一路径。Behura等［16］指出，具有同一水平慢度的反射在t－x域中表示为具有相同瞬时斜率的不同反射。同时，在把数据变换到τ－p域的过程由于综合效应能压制部分随机噪声。

二维空间的τ－p变换的数学表达式为［17－18］

对于离散数据有

Van der Baan等［19］指出，对于τ－p 域 一个给定的层，其校正量τi为

式中：τ0是垂直入射波的旅行时；υi是层速度；ηi是Alkhalifah等［20］定义的各向异性参数。一个水平层的总校正量为

τ－p校正曲线表示的仅是平面波的垂直传播时间部分。然而，为了提取Q值，需知总的传播时间。由τ－p变换公式t＝τ＋px和（9）式可知

将（8）式代入（10）式有

由此，对于τ－p域中任意给定的一个点，给定一个p值，运用（11）式和（12）式就可以找到t－x域中对应的点［21］，这样就可以在t－x域追踪到图2中相同水平斜率的反射波。

（11）式和（12）式是针对各向异性走时方程建立的，其中含有各向异性参数ηi，所以，我们提出的方法对各向异性介质具有适用性。如果仅对各向同性介质数据进行处理，可默认ηi为0。

1.2.2 广义S变换

广义S变换是由S变换发展而来。地震波在地层中传播时，高频成分要比低频成分衰减严重，主频变低；S变换的窗函数固定不变，不能随着信号本身的特点进行调整。我们采用了一种改进的广义S变换方法，其变化的时窗可以满足吸收衰减分析的需求，即

当λa和p变化时，广义S变换的窗函数形状随频率而变化。具体来说，窗函数的宽度随频率呈反比变化，幅值与频率呈正比变化。低频段时窗较宽，幅值较低，具有较高的频率分辨率；高频段时窗较窄，幅值较高，可获得较高的时间分辨率［22］。

1.2.3 算法实现步骤

图3给出了叠前高精度Q值估计方法的计算流程。如图3所示，对一个叠前CMP道集数据，在给定的水平慢度范围内依照设置的间隔选取p1，p2，…，pn，依照（11）式和（12）式，把一个恒定水平慢度pi产生的反射（pi，τi）映射到t－x 域（ti，xi），公式中的τ0i，υi，ηi由高精度速度分析得到［23］，也就是选取目标层位顶、底面在相同的表层路径（即恒定pi）下所对应的道集xi1，xi2。这里仅对各向同性介质进行处理，即

图3 叠前高精度Q值估计方法计算流程

2 模型数据测试

对图4a所示的4层地质模型（各层的Q值依次为100，50，100，150）采用射线追踪法进行正演模拟。模拟采用的道间距为20m，炮检距为0～1 180m；时间采样间隔为1ms；所用子波为50Hz的雷克子波。模拟出的含衰减零偏移距CMP地震道集记录见图4b。

将图4b的CMP道集记录变换到τ－p域（图5），地震记录中接收到的第1层的反射波衰减程度仅由第1层的走时决定，故可用t－x域中不同地震道之间的频谱比斜率值反演求取该层Q值。对第2层、第3层和第4层分别选取的水平慢度范围为0～0.3，0～0.2，0～0.1，间隔都为0.001；采用本文提出的高精度方法计算这3层的Q值。

图4 地质模型（a）及其正演模拟地震记录（b）

图5 转换到τ－p域中的地震记录

表1给出了采用本文方法与用QVO（两步拟合算法）方法计算的结果及其误差比较，可以明显看出本文方法的精确性。

为了验证本文方法的抗噪性，给上述模型的理论数据加入30dB的随机噪声，结果如图6所示。图7是水平慢度为0.09的频谱。图8给出了第2层、第3层的振幅谱比值与频率的拟合关系。

表2给出了加入噪声后用本文方法与用QVO方法计算的结果及其误差，可见加入噪声后本文方法仍保持较高的精度。比较加入噪声前（表1）、后（表2）的误差可以看出，本文方法抗噪性较强，加入噪声后计算的Q值波动不大。

表1 两种方法计算出的Q值及其误差

表2 加入30dB噪声后两种方法计算出的Q值及其误差

图6 加入30dB噪声后的t－x 域（a）和τ－p 域（b）地震记录

图7 p＝0.09时模型4个界面的反射波频谱

图8 模型第2层（a）和第3层（b）的振幅谱比值与频率的拟合关系

3 实际资料处理

为了验证本文方法在实际应用中的可行性和有效性，选取了一个实际CMP道集数据进行Q值试算。图9为经过去噪预处理之后的CMP道集，经高精度速度分析后扫描出t0和vNMO，解释出反射层位，图9中A，B，C，D，E为各层反射界面。

为了验证本文方法的精确性，用计算出的Q值对CMP道集数据作反Q滤波处理。图10显示了截取的一段反Q滤波前、后的数据，通过对比可以看出，处理后的数据纵向分辨率有了较明显的提高，很多层间小薄层经过反Q滤波后能明显分辨。可见本文方法在实际数据应用中的效果较好，可用来做精确的反Q滤波处理，以提高地震记录的分辨率。

图9 经去噪预处理后的实际CMP道集

图10 反Q滤波前（a）、后（b）的一段道集记录

4 结束语

我们提出的叠前高精度Q值估计方法在理论上具有一定的先进性，相对于传统谱比法的两步拟合计算，叠前高精度Q值估计方法作出了两个重要的改进：①为使参考信号和测量信号具有相同的表层传播路径，对同一水平慢度的道集进行处理，避免了传统方法的两步拟合所带来的统计性不精确；②采用变时窗的时频变换方法，最大限度地减少相邻层位的干扰效应。

通过理论模型数据测试证明了本文方法的精确性，实际数据的应用结果也证明了该方法的有效性。

［1］李生杰，施行觉，王宝善，等.地层衰减在地震记录上的特征分析［J］.石油地球物理勘探，2002，37（3）：248－253 Li S J，Shi X J，Wang B S，et al.Characteristics analysis of formation attenuation in seismic record［J］.Oil Geophysical Prospecting，2002，37（3）：248－253

［2］魏文，王小杰，李红梅.基于叠前道集小波域Q值求取方法研究［J］.石油物探，2011，50（4）：355－360 Wei W，Wang X J，Li H M.Analysis of Qcalculating method based on the prestack gather in wavelet domain［J］.Geophysical Prospecting for Petroleum，2011，50（4）：355－360

［3］Tonn R.The determination of the seismic quality factor Qfrom VSP data：a comparison of different computational methods［J］.Geophysical Prospecting，1991，39（1）：1－27

［4］付勋勋，徐峰，秦启荣，等.基于改进的广义S变换求取地层品质因子Q值［J］.石油地球物理勘探，2012，47（3）：451－467 Fu X X，Xu F，Qin Q R，et al.Study on quality factors Qby modified generalized Stransform［J］.Oil Geophysical Prospecting，2012，47（3）：451－467

［5］Dasgupta R，Clark R A.Estimation of Qfrom surface seismic reflection data［J］.Geophysics，1998，63（6）：2120－2128

［6］Hackert C L，Parra J O.Improving Qestimation from seismic reflection data using well－log－based localized spectral correct ion［J］.Geophysics，2004，69（6）：1521－1529

［7］王小杰，印兴耀，吴国忱.基于叠前地震数据的地层Q 值估计［J］.石油地球物理勘探，2011，46（3）：423－428 Wang X J，Yin X Y，Wu G C.Estimation of quality factors Qfrom pre－stack data［J］.Oil Geophysical Prospecting，2011，46（3）：423－428

［8］Clark R A，Carter A J，Nevill P C，et al.Attenuation measurements from surface seismic data—azimuthal variation and timelapse case studies［J］.Expanded Abstracts of 63rdEAGE Annual Conference，2001，L28

［9］Clark R，Benson P M，Carter A J，et al.Anisotropic P－wave attenuation measured from a multi－azimuth surface seismic reflection survey［J］.Geophysical Prospecting，2009，57（5）：835－845

［10］Reine C，Clark R，van der Baan M.Robust prestack Q－determination using surface seismic data：part 1—method and synthetic examples［J］.Geophysics，2012，77（1）：R45－R56

［11］Reine C，Clark R，van der Baan M.Robust prestack Q－determination using surface seismic data：part 2—3Dcase study［J］.Geophysics，2012，77（1）：B1－B10

［12］Reine C.A robust prestack Q－inversion in theτ－p domain using variable window spectral estimates［D］.England：University of Leeds，2009

［13］Reine C，van der Baan M，Clark R.The robustness of seismic attenuation measurements using fixed－and variable－window timefrequency transforms［J］.Geophysics，2009，74（2）：WA123－WA135

［14］刘喜武.基于广义S变换研究地震地层特征［J］.地球物理学进展，2006，21（2）：215－218 Liu X W.Study on characteristics of seismic stratigrephy by generalized S transform［J］.Progress in Geophysics，2006，21（2）：215－218

［15］Pinnegar C R，Mansinha L.The S transform with windows of arbitrary and varing shape［J］.Geophysics，2003，68（1）：381－385

［16］Behura J，Tsvankin I.Estimation of interval anisotropic attenuation from reflection data［J］.Geophysics，2009，74（6）：A69－A74

［17］Chapman C H.Generalized radon transforms and slant stacks［J］.Geophysical Journal of the Royal Astronomical Society，1981，66（2）：445－453

［18］Schultz P S，Claerbout J F.Velocity estimation and downward continuation by wavefront synthesis［J］.Geophysics，1978，43（4）：691－714

［19］Van der Baan M，Kendall J M.Estimating anisotropy parameters and traveltimes in theτ－p domain［J］.Geophysics，2002，67（4）：1076－1086

［20］Alkhalifah T，Tsvankin I.Velocity analysis for transversely isotropic media［J］.Geophysics，1995，60（5）：1550－1566

［21］Van der Baan M.Processing of anisotropic data in theτ－p domain：I—geometric spreading and moveout corrections［J］.Geophysics，2004，69（3）：719－730

［22］Stockwell R G，Mansinha L，Lowe R P.Localization of the complex spectrum：the S transform ［J］.IEEE Transactions on Signal Processing，1996，44（4）：998－1001

［23］Tsvankin H D，Thomsen I.Nonhyperbolic reflection moveout in anisotropic media［J］.Geophysics，1994，59（8）：1290－1304