地震波在传播的过程中，大地对信号吸收衰减.实际介质中引起地震波衰减的因素很多，总的来说主要分为固有衰减和非固有衰减(熊翥，2008).固有衰减是反映介质内在属性的地层本征衰减，即地层的粘弹性衰减；非固有衰减与地震波传播特性有关，包括球面扩散引起的衰减、介质非均匀性引起的衰减、界面反射和透射系数及薄层调谐效应引起的衰减等.通常用品质因子Q值来描述介质本身所固有的衰减特性，与其内部的结构特征，以及孔隙度，流体饱和度和渗透率等介质参数密切相关.Q值有大量的应用空间，能用来研究岩性及含流体属性分析、储层识别及烃类检测，以及应用反Q滤波来提高地震资料的分辨率(马昭军和刘洋，2005；Raji and Rietbrock，2013).

关于Q值的估计，前人进行了大量的研究，根据垂直地震剖面(VSP)数据、叠后地震记录和叠前CMP数据等，研究了多种方法来估计Q值，分为时间域方法、频率域方法和反演类Q值估计方法，并根据数据类型，对一些方法进行了改进和创新.相比叠后的资料来说，叠前的资料含有更多、更精确的层位信息以及振幅、走时和频率信息，目前，国内外已有许多学者对叠前数据的Q值估计进行了研究(Hauge，1981；Tonn，1991；Reine et al.，2012).

Dasgupta和Clark(1998)利用谱比法从叠前CMP反射道集估计了品质因子Q值，并对计算Q值的流程进行了详细说明；Hackert和Parra(2004)对该方法进行了改进，提出了对振幅谱进行校正后提取Q值的方法(LSR).Zhang和Ulrych(2002)基于直射线路径假设提出了一种新的使用峰值频率从叠前CMP道集求Q值的方法，陈爱萍等(2004)将此方法应用于叠前CMP资料估计Q值，验证了其可行性，但是应用于实际数据时，受噪音影响较大.李雪英等(2009)提出了叠前等效Q值扫描的方法，该方法利用广义S变换进行时频分析，具有一定的鲁棒性，是宏观介质的等效Q值分析方法，但是扫描范围难以确定，估计精度和效率会受到影响.Yan和Liu(2009)提出了利用谱比法估算叠前反射纵波和转换横波的方法，然后利用估算的Q值进行反Q滤波，在实际数据中估计转换横波的Q值难度较高.袁恩辉和顾汉明(2010)提出了基于小波域的叠前Q值估计方法，该方法结合小波变换，利用谱比法提取等效地层品质因子，然后通过剥层法估算层品质因子，实际计算时，应注意频带范围选择的问题.Zhao和Gao(2011)，赵静等(2013)利用子波包络峰值瞬时频率法估计叠前CMP道集的Q值，该方法精度较高，抗噪性也较强，但是对时窗和子波都具有一定的依赖性.魏文等(2011)提出了一种从叠前道集中估计Q值的QVO(品质因子Q随偏移距变化)方法，他改进了谱比法，在小波域中提取振幅谱，计算非零偏移距频谱比斜率值，然后根据求得的频谱比斜率值随偏移距的平方的近似线性关系，求得零偏移距处的Q值.王小杰等(2011)提出了一种与S变换结合的估计叠前地震资料的地层Q值的方法，该方法先借助于S变换对地震道做时频分析，求取频谱比斜率随非零偏移距处的变化关系，经炮检距归零处理后得到零炮检距处的地层Q值，也属于QVO方法.但是，在QVO方法的实现过程中，Q值估计的精度会受到时窗和频谱比斜率的影响，精度不高.Reine等(2012a)提出了一种抗噪性较强的叠前Q值估计方法，首先利用变化时窗的时频变换减少干扰谱，然后利用频率和走时各向异性同时反演减少剩余干扰的影响后，进行变换，逐道的计算其水平慢度，该方法已经用于3D地震数据.Raji和Rietbrock(2013)年提出了一种计算质心频率法估计Q值的新方法，该方法借助变化的时窗长度和傅里叶变换计算单个信号的质心频率，然后根据目的层上下界面的反射波质心频率和旅行时间的差异来估计Q值，该方法需要估计准确的质心频率.王德利和戴建芳(2013)提出了一种基于射线路径的叠前高精度Q值估计方法，该方法在谱比法的基础上，首先在域中对叠前数据进行处理，并引入走时各向异性参数和利用了广义S变换，避免了传统方法中时窗长度和宽度选择对计算结果的影响.

对数谱比法、质心频移法等方法估计Q值都受到加时窗问题的影响，峰值频率法和解析信号法对噪声很敏感，多数改进的方法都是在峰值频率法、谱比法和质心频移法的基础上，与不同的时频分析方法结合，但仍然不能突破每种方法的局限性.本文在总结上述方法优缺点的基础上，选取小波变换，利用小波理论与地震波传播方程相结合形成的尺度能量方程，把小波振幅谱与峰值频率结合起来估计叠前的Q值，其核心就是根据峰值频率法估计的Q值给定一个Q值的估计范围，利用补偿前后的小波能量谱的相似性来搜索这个范围内的Q值(Yang and Gao，2008)，得到精确的Q值，该方法不受几何扩散的影响，能避免频带范围选择的问题，而且具有一定的抗噪性，估计的Q值精度较高.

假设地下介质为水平层状均匀粘弹性介质，且Q值与频率无关，地震波传播一段时间后其振幅谱为(张显文，2007)

其中，B(t)是与频率和吸收无关的影响.对(1)式求导数，推导出品质因子Q与峰值频率f_p之间的关系为

式中f_m为源子波的峰值频率，f_p为接收信号的峰值频率.

式(2)是利用沿着不同炮检距峰值频率来估计各道的Q值，因此可以消除地表起伏效应和随机噪声，提高Q值的精度.

实信号f(t)对于解析小波ψ(t)的连续小波变换为(李宏兵等，2004；陈文超等，2009)

式中，参数a，b分别为尺度因子和平移因子.母小波ψ(t)的傅氏变换为(ω)，满足“容许”条件，地震资料处理时，通常要研究信号的振幅谱和相位谱，这里采用修正后的复数型Morlet小波(高静怀等，1996)：

式中，Ω₀为调制频率，Ω₀≥5，c为调幅因子，它控制了小波函数的长度.

由于地震波波形的变化与地层吸收有关，地层吸收的强弱可用品质因子Q来定量表示.所以可用目的层上下界面的反射波波形的相似程度来估计该地层的品质因子.假设震源子波为理想的脉冲源，在粘弹性介质中传播一段时间t后，且忽略由频散效应引起的衰减(张显文等，2010；张固澜等，2014)，可以推导出小波的尺度能量谱为(赵伟和葛艳，2008；魏文等，2011)

式中E_a为尺度能量.由上式知，在小波尺度域信号的能量与品质因子Q、尺度因子a和传播时间t有关.

根据地震波在粘弹性介质中的吸收衰减传播规律，t₁时刻的尺度能量谱为E_a1，对t₂时刻的小波尺度能量谱E_a2进行补偿.补偿后的t₂时刻的尺度能量谱为

若由公式(7)估计的经过补偿后的t₂时刻的尺度能量谱的波形、峰值频率与t₁时刻的实际尺度能量谱E_a1的波形、峰值频率相似，则认定此时的Q值便为该目的层的Q值.

确定两条曲线的相似程度时，用相关系数来定量表示，当相关系数为0时，相似度最差，即不相关.当相关系数为1时，则两个信号相似度很好，这里给出一个定量的衡量两条曲线的相似程度的公式为

因此，只要搜索R达到最大值时的Q值，即两条曲线最相似时的Q值，即为估计的品质因子Q.R是与Q有关的非线性函数，理论上，当两条曲线完全相似时，R取得最大值为1，即经该Q值补偿后的t₂时刻的尺度能量谱的波形、峰值频率与t₁时刻的实际波形、峰值频率完全相似.文中，把补偿后的t₂时刻的能量谱称为补偿谱，t₁时刻的实际能量谱称为参考谱.

在无衰减介质中，对于Ricker子波，其反射信号的能量在每个尺度都相同，在衰减介质中，其反射信号的能量随尺度的增大呈指数升高.图 1a为在无衰减介质以及Q=50的衰减介质中的尺度能量E_a的分布图.取Ω₀为100，c=30时，取Q从1到300，根据公式(7)，计算衰减补偿后的尺度能量谱，与无衰减介质中的实际的参考尺度能量谱E_a1比较，利用式(8)计算相关系数，每一个Q值的相关系数如图 1b所示，在Q=50时，相关系数为1.此时，补偿的尺度能量谱与参考能量谱最相似.

图 1

Fig. 1

图 1(a)Ricker子波尺度能量;(b)相关系数曲线图Fig. 1(a)Scales energy spectrum of Ricker wavelet;(b)The correlation coefficient curve

搜索R最大时的Q值，需要确定一个Q值的搜索范围，在这个范围内，计算每一个经过Q值补偿后的t₂时刻的估计尺度能量谱，比较和E_a1之间的相似性，当他们的振幅谱、峰值频率最相似时，即R最大时，补偿的Q值即为此目的层的品质因子.

这里，引入峰值频率法，利用峰值频率法估计出的Q值作为初始Q₀值，由这个初始Q₀值确定一个Q值的搜索范围，在这个范围内，利用公式(8)，搜索出精确的Q值.

设计五层水平层状介质，得到含衰减的叠前CMP道集(图 2)，震源子波为Ricker子波，主频为30 Hz，炮检距为50 m，层间Q值为30、50、70、90和60.通过小波变换得到小波尺度能量谱，并且估计接收信号的振幅谱的峰值频率.

图 2

Fig. 2

图 2 合成含5层反射同相轴的叠前CMP衰减道集Fig. 2 Synthetic 5 layers pre-stack CMP gather

这里以第一道第二层为例，应用本文方法估计Q值.首先通过小波变换，得到第二层上下界面的振幅谱和峰值频率，如图 4所示，图 3为零偏移距道的小波变换的频谱.如图 5为第一道用PFS法的估计的Q值与理论Q值.可以看出，估计的Q值与理论的Q值相差较大，精度不高.这是由于峰值频率法受到峰值频率的精度的影响很严重，所以估计的Q值精度较低，因此，引入小波能量谱.

图 3

Fig. 3

图 3 第一道的小波变换频谱Fig. 3 The wavelet transform spectrum of the first trace

图 4

Fig. 4

图 4 无噪情况下，第二层上下界面的振幅谱Fig. 4 The second layer amplitude spectrum of the front and bottom with no noise

图 5

Fig. 5

图 5 峰值频率法估计的Q值与理论Q值的关系曲线Fig. 5 The curve of theoretical Q value and estimated Q value using PFS method

把由峰值频率法估计的Q值，作为零偏移距道的每层Q₀值，根据初始Q₀值给定Q值的大致范围，这里给定10~400.在给定的Q值范围内，利用公式(7)对下界面的振幅谱进行补偿，比较补偿谱和上界面的参考谱的振幅谱和峰值频率，利用公式(8)，搜索出振幅谱最相似时的Q值.根据小波能量谱和峰值频率估计的Q值，即满足了峰值频率的关系，又满足了小波能量谱的相似性，确保了Q值估计的准确可靠.第一道第二层的补偿谱和参考谱的相关系数曲线图如图 6所示，当估计Q值为50时，相关系数最大，R为0.9983，接近于1，估计的Q值精度很高.

图 6

Fig. 6

图 6 相关系数曲线图Fig. 6 The correlation coefficient curve

同时，测试初值Q₀对该方法的计算效率和精度的影响.计算初值Q₀与估计Q值之间的偏差，与搜索次数之间的关系.如图 7a所示，当初值Q₀为58时，搜索次数为3时，即可得到相关系数为最大时的Q值，此时偏差为8，估计的Q值为初值Q₀减去最终偏差为50.当初始Q₀值为5时，需要搜索13次，如图 7b所示，最终偏差为-45，估计的Q值为初值Q₀减去最终偏差为50.因此采用本文方法估计Q值，在保证估计的精度的情况下，可大大减少搜索次数，提高计算效率.

图 7

Fig. 7

图 7 搜索过程中的偏差与搜索次数的关系(a)初值Q0为58(b)初值Q0为5Fig. 7 The relation of differ and searched times(a)original Q0is 58(b)original Q0 is 5

为测试本文方法的抗噪性，在合成地震记录中加入35 dB的噪声.表 1和表 2分别为无噪情况下和加噪情况下，本文方法与峰值频率法估计的Q值.可以看出与峰值频率法相比，本文方法估计的Q值精度有明显的提高，说明本文方法有更强的抗噪性.

表 1(Table 1)

表 1 无噪情况下，峰值频率法与本文方法估计的Q值结果及误差 Table 1 The estimated Q values and error by two methods with no noise 层数 峰值频率法 误差(%) 本文方法估计Q值 误差(%) 1 32.80 10.67 30.50 1.67 2 58.22 16.44 50.00 0 3 80.99 15.70 72.10 3.00 4 87.31 2.99 91.88 2.09 5 58.10 3.17 59.52 0.80

表 1 无噪情况下，峰值频率法与本文方法估计的Q值结果及误差 Table 1 The estimated Q values and error by two methods with no noise

表 2(Table 2)

表 2 加入噪声后峰值频率法与本文方法估计的Q值结果及误差 Table 2 The estimated Q values and error by two methods when synthetic record contains noise 层数 峰值频 率法 误差 (%) 本文方法 估计Q值 误差 (%) 1 32.82 10.67 31.12 3.73 2 58.22 16.44 51.92 3.84 3 50.42 27.97 73.20 4.57 4 80.30 10.78 88.76 1.38 5 81.63 36.05 62.38 3.97

表 2 加入噪声后峰值频率法与本文方法估计的Q值结果及误差 Table 2 The estimated Q values and error by two methods when synthetic record contains noise

为验证本文方法在实际应用中的可行性和有效性，把该方法应用与实际数据，选取某工区保幅的叠前CMP道集，如图 8a所示，共85道，其目的层大约在900 ms至1130 ms之间，储层岩性为砂岩，可能含有油气.提取该CMP道集的第6道，如图 8b所示.

图 8

Fig. 8

图 8 叠前CMP道集(a)及其第6道(b)Fig. 8 Prestack CMP gather(a) and the sixth trace(b)

用本文方法求取该道的Q值，选取977 ms和1152 ms两个反射时间，通过小波变换，得到经过归一化处理后的振幅谱，如图 9a所示，用峰值频率估计的Q值为26，作为本文方法估计的初始Q₀值.选取Q值的搜索范围为5至300，本文方法估计的该层的Q值为60，相关系数曲线如图 9b所示.提取该道724 ms到2280 ms之间的各层的Q值，通过峰值频率法估计的Q值曲线，如图 10a所示，本文方法估计的Q值曲线，如图 10b所示.从图 10b可以看出在目标层附近有较大的衰减，可能有储层.但是通过峰值频率法估计的多层的Q值都较小，不能分辨出储层.这里需要注意的是，当峰值频率法估计的Q值较大时，如图 10a中Q值为350时，不能作为初始Q值，此时给定的初始Q值为5.

图 9

Fig. 9

图 9(a)目的层上下界面的振幅谱;(b)相关系数曲线图Fig. 9(a)The second layer amplitude spectrum of the front and bottom;(b)The correlation coefficient curve

图 10

Fig. 10

图 10(a)峰值频率法估计的叠前CMP道集的衰减曲线;(b)本文方法估计的衰减曲线Fig. 10 The estimated Q-curve from prestack CMP data(a)Peak frequency method(b)the new method

为进一步验证本文方法求取的Q值的精确性，用计算出的Q值对CMP道集进行反Q滤波.截取一段反Q滤波前后的数据，如图 11所示，可以看出反Q滤波后的剖面纵向分辨率有了明显提高，可以分辨出很多层间的小薄层.因此，本文方法在实际数据中求取的Q值，可用来做精确的反Q滤波处理，能有效的提高地震资料的纵向分辨率.

图 11

Fig. 11

图 11 CMP道集反Q滤波前后的效果对比图(a)反Q滤波前(b)反Q滤波后Fig. 11 The results contrast figure of inversing Q filtering(a)before inversing Q filtering(b)after inversing Q filtering

5.1    利用本文方法估计叠前CMP道集的品质因子Q，充分的把峰值频率和小波振幅谱结合起来.理论模型测试得出，与峰值频率法相比，有效的提高了叠前Q值估计的精度及效率，并且该方法有较好的抗噪性和较高的计算效率.实际数据结果表明，该方法能很好估计实际叠前CMP道集的Q值，通过反Q滤波处理很好的提高了地震资料的纵向分辨率.

5.2    本文方法估计实际数据的Q值时，选取Q值的范围最大值为350，若Q值大于350时，相关系数还不能取到最小值，即把最大的Q值作为该层的Q值，选取Q值的最大范围时，350~500左右即可，大于这个范围的即可认为无衰减.

5.3    选取的叠前CMP道集应为做过精细去噪，静校正等保幅处理的数据，没有进行过反褶积处理.

致 谢    感谢张生强博士、凌云博士给予第一作者的指导和帮助.