2020, Vol. 63
地下介质的非弹性和非均质性会导致地震波传播过程中发生衰减和频散(Ricker，1953；Futterman，1962；White，1983；Kneib and Shapiro, 1995).不同频率的地震波衰减程度不同，这种与频率有关的地震波衰减特征可以通过品质因子Q来描述.地震波衰减与介质的物性参数有关，如岩性、饱和度、渗透率和孔隙压力等(Müller et al., 2010).研究表明地震波衰减对含气饱和度参数比较敏感(Winkler and Nur, 1982)，因此Q值可以用于气藏的探测与识别，同时根据地震衰减随方位角变化特征还可以预测裂缝方向等(Carter et al., 2001；Maultzsch et al., 2007).在地震数据处理领域，由于地震波衰减引起的频谱变化会降低地震数据的分辨率(Wang，2004)，Q值一直是地震衰减补偿提高分辨率处理的重要参数.地震衰减因子既可以用于地震数据的衰减补偿，还可以为油藏描述和流体预测提供参考信息，因此准确估计Q值对地震数据处理和解释都非常重要.
近几十年来，许多地球物理学家对地震品质因子Q值估计进行了研究并提出了一些方法.最常见的方法之一是对数谱比方法(Bth，1974)，该方法通过利用传播过程中的两个地震波波形的频谱比来获得Q，但在实际应用中谱比法对噪声比较敏感，具有不稳定性(Tonn，1991；Wang et al., 2015).Tonn(1991)利用谱模拟和谱比法研究了垂直地震剖面(Vertical Seismic Profiling，VSP)数据的Q值估计方法.Quan和Harris(1997)提出了频移法估计VSP数据的Q值.Zhang和Ulrych(2002)提出了峰值频率法估计地震数据的层Q值.Blias(2012)针对VSP数据通过建立待估计的Q值和理论计算的Q值的误差目标函数提出了VSP数据层Q值估计方法.Wang等(2015)基于地震波的对数谱面积差提出了一种Q值估计的方法.Liu等(2018a)利用整形正则化技术在谱比法中引入频率相关的线性拟合，提出了一种改进的稳定谱比法估计Q值的方法.
在谱比法估计Q值中一般采用时频分析获取地震数据不同传播时刻的瞬时振幅谱.Reine等(2009)比较了谱比法估计Q值中的四种时频变换，即短时傅里叶变换(STFT)、Gabor变换、S变换和连续小波变换，指出不同的变换具有不同的时频表征能力.Hao等(2016)利用广义S变换研究了地层Q值估计方法.Sun等(2015)利用改进的S变换研究了叠前道集的Q值估计方法.在时频分析领域，Liu等(2011)通过最小平方反演结合正则化约束提出了一种新的时频表征方法，并在地震数据砂体探测及低频阴影识别中得到了应用.Liu和Fomel(2013)将其扩展到更稳健的形式，并在地震数据噪声压制和振幅均衡中得到了应用.基于最小平方方法的时频分析方法已经被应用到叠后地震数据Q值估计(刘喜恒等，2016)和衰减补偿中(刘洋等，2016)，本文在多射线联合反演叠前数据Q值估计中也采用了这种基于最小平方反演的局部时频分析方法，该方法可以更直观的调整正则化参数控制时频谱的平滑特征.
叠前或叠后的地震数据均可用于估计Q值，叠后地震数据信噪比较高，在采用谱比法估计Q值时稳定性较好，但叠后地震数据没有考虑真实的射线传播路径，动校正叠加处理也会引起波形发生变化，这些因素均会影响Q值估计的精度(Dasgupta and Clark, 1998；Reine，2012).叠前地震数据可以保持原始的地震波形，如果不考虑信噪比因素，理论上叠前地震数据估计的Q值更精确，但叠前地震数据信噪比较低，因此针对叠前地震数据需要研究稳定可靠的Q值估计方法.另外，由于地震波在介质中的射线路径不一定是直的，叠前Q值估计还需要考虑射线路径问题.针对叠前地震数据Q值估计问题，Dasgupta和Clark(1998)在动校正后道集上提出了Q随炮检距变化的估计方法(QVO).Zhang和Ulrych(2002)提出了频移法在叠前共中心点(CMP)道集上估计地层Q值，该方法建立了峰值频率和品质因子之间的关系，但该方法将波场的传播路径假设为直射线，降低了最终Q值估计的精度.Behura和Tsvankin(2009)提出了一种逐层剥离的方法来估计反射地震数据的层Q值.Reine等(2012)考虑到τ-p域地震数据和射线的关系提出了τ-p域叠前地震数据Q值估计方法.Oliveira等(2017)提出了基于波场延拓的叠前Q值估计方法，该方法首先利用速度信息将波场延拓至目标层处，然后利用常规的方法(峰值频率法)对其进行估计.
当前存在的叠前Q值估计的方法一般需要速度信息将非零炮检距的Q值映射到零炮检距处进行反演，而速度信息的准确性会影响到Q值估计的精度.同相轴局部斜率隐含了地下介质的构造和速度等信息，已经在地震数据处理领域得到了广泛的应用，如速度无关时间域地震成像(Fomel，2007)、各向异性动校正(Burnett and Fomel, 2009)、走时参数估计(Liuet al., 2018b)及VSP波场分离(杜婧等，2009)等.局部斜率与射线参数具有明确的关系，具有相同局部斜率的地震反射波具有相同的传播射线参数，因此局部斜率信息可以被用来找到具有相同射线参数(即相同射线路径)的反射波形，这些波形可以用于谱比法估计Q值.基于此，本文在局部斜率属性基础上提出了多射线联合反演的速度无关Q值估计方法，主要利用局部斜率信息确定多射线路径的反射波形，通过预测映射技术(Fomel，2010)将非零偏移距走时映射到零炮检距走时，然后利用局部时频变换计算时频谱，将不同射线参数相同层位的反射联合反演求解Q值.本文首先给出了衰减层状介质中射线传播与同相轴局部斜率的关系，然后回顾了同相轴斜率估计和预测映射技术基本理论，在此基础上建立了多射线参数反射联合反演Q值方法，最后利用模拟和实际算例验证了本文提出方法的有效性和可靠性.
1 方法原理 1.1 衰减介质中的射线传播考虑地震波在水平层状介质中的传播情况，图 1显示了CMP道集的走时曲线和相应的射线路径.假设第二层为目的层，射线路径CO1C′和ABO2B′A′具有相同的射线参数(对应CMP中M和N点处的反射)，射线ABO2B′A′分别在第一层和第二层中传播，其中BO2B′是在第二层中传播的部分，如果计算第二层的Q值，应该考虑波场在BO2B′这一段传播引起的波形振幅谱的变化.假设震源子波一致，传播路径AB段和CO1段波场传播影响一样，而B′A′段和O1C′段影响一样，因此N点和M点的频谱衰减差异应该是由第二层的传播BO2B′引起的，这里说的频率衰减差异仅考虑频率有关的频谱衰减，不包括反射系数等引起的频谱差异.因此需要在CMP道集中找到具有相同射线参数的这些反射波形(如M点和N点)，注意到走时曲线中M和N点具有相同的斜率(视速度相同)，如果我们选用经典的谱比法估计Q值，在某一射线参数p下有如下的对数谱比公式(Tonn，1991)：
|
图 1 波场射线传播和CMP道集走时示意图 Fig. 1 A sketch of ray propagation and traveltime curves in CMP gather |
|
(1) |
其中St2(f, p)和St1(f, p)分别是位于N和M处的反射波形的振幅谱，Δt(p)是该射线参数p的传播走时差，Δt(p)=TABO2B′A′－TCO1C′，TABO2B′A′和TCO1C′分别对应走时曲线中的TN和TM.式(1)中G是频率无关的，代表着反射、透射及球面扩散引起的振幅变化.
在叠前CMP道集上估计Q值，走时差Δt(p)是不容易获得的，因为它与射线参数、层速度和层厚度有关.常规的叠前估计Q值的方法一般在某一个偏移距道上选取波形做谱比，这样走时差不是真实的传播时差，并且波形也是不沿着同一条射线传播的，不是真实在第二层中传播的路径，因此估计的第二层的Q值精度不高.从CMP道集上获得Δt(p)需要确定CMP道集中M和N的位置，而注意到M和N点处的地震反射同相轴具有相同的局部斜率，因此可以利用局部斜率相同来确定具有相同射线参数反射的位置，找到这些反射后根据它们的走时可以计算走时差Δt(p)=TN－TM，而该走时差是目的层中射线路径BO2B′的传播时间，如果已知同相轴局部斜率就可以不用速度信息获得相同射线参数下来自不同反射界面的反射波.
1.2 地震同相轴局部斜率估计在本方法中，局部斜率不仅被用来计算相同射线参数的走时差Δt(p)，而且在预测映射中也用到了局部斜率参数，因此局部斜率的估计是本文方法关键的一步.目前地震数据处理领域有许多估计走时局部斜率的方法，主要包括结构张量法(Bakker et al., 1999；Bakker，2002；Hale，2009；Wu and Janson, 2017)、平面波预测法(plane-wave destruction，PWD)(Fomel，2002)和相似扫描法(Marfurt，2006)等.由于平面波预测法(Fomel，2002)比其他数值微分的方法在估计局部斜率时更具有抗噪性(Schleicher et al., 2009；Khoshnavaz et al., 2016)，因此本文选用平面波预测法来估计同相轴局部斜率.
Fomel(2002)在Z变换域中构造了一个二维预测误差滤波器
|
(2) |
其中
|
(3) |
该滤波器C(Zt, Zx)是局部斜率p(t, x)的函数，可以通过求解如下的最小二乘方程来实现局部斜率估计：
|
(4) |
其中d是地震道，C[p(t, x)]代表预测误差滤波器C(Zt, Zx)卷积算子，滤波器C(Zt, Zx)是局部斜率的非线性函数.式(4)表明地震数据是可以预测的，通过算子C(Zt, Zx)可以建立预测误差方程，在该预测误差方程中p(t, x)是待求的未知参数.由于滤波器C(Zt, Zx)是局部斜率p(t, x)的非线性函数，因此利用式(2)求解局部斜率p(t, x)本质上是非线性优化问题.Fomel(2002)利用线性化方法(Gauss-Newton迭代)将其转化为线性最小二乘问题
|
(5) |
其中Δp(t, x)是迭代求解局部斜率过程中的增量.如果给定初始局部斜率p0(t, x)，则可以根据式(5)通过线性共梯度迭代方法求解Δp(t, x)，考虑到Δp(t, x)的平滑特征，可以对其进行正则化约束以保证估计的局部斜率是平滑的.
1.3 零炮检距走时预测映射尽管局部斜率解决了寻找相同射线参数的反射波形问题，但是还没有将其映射到准确的零偏移距走时处，而期望获得的层Q值是零偏移距的函数，即Qt0，因此需要采用预测映射技术将非零炮检距的反射映射到零偏移距的走时刻，才可以获得零炮检距Q值，即时间域的层Q值Q(t0).
预测映射技术(Fomel，2010)通过走时局部斜率信息将具有相同零炮检距走时的同相轴追踪出来.考虑CMP道集的某条反射走时曲线t=φ(t0, x)，它是零炮检距走时t0和炮检距x的函数，与走时局部斜率有如下的关系(Wu and Fomel, 2018)：
|
(6) |
可以看出，走时局部斜率p(t, x)是走时曲线t=φ(t0, x)对炮检距x的偏导数，同时走时曲线满足在零炮检距处为t0.如果已知走时局部斜率p(t, x)通过数值求解微分方程(6)可以获得零炮检距t0对应的走时曲线t=φ(t0, x).预测映射技术是通过给出的种子点(如零偏移距走时，图 1中的t0, M或t0, N)来预测非零炮检距处的走时曲线(Fomel，2010；Wu and Fomel, 2018)，基本算法是通过递归的方式从第1道预测第k道
|
(7) |
其中Pi, j描述了由第i道数据预测第j道数据的预测算子.地震道之间的预测是通过沿着同相轴的方向利用局部斜率对地震道进行时移完成的，预测算子可以看作微分方程(6)的直接求解.利用预测映射获得非零炮检距处反射对应的零炮检距时间，这样可以将非零炮检距估计的Q值映射到零炮检距时间上.
1.4 多射线Q值联合反演通过局部斜率属性和预测映射，可以获得目的层顶底界面的反射波频谱St2(f, p)和St1(f, p)，然后就可以利用谱比法求取对应的零炮检距走时的Q值.从CMP道集中获取反射波的频谱由于受到波形干涉的影响，一般采用时频分析的方法减弱波形互相干涉的影响.本文采用局部属性时频分析方法(Liu et al., 2011；Liu and Fomel, 2013)获得时频谱.考虑[0, L]范围内的信号f(t)，利用傅里叶级数可以表示为
|
(8) |
其中øk(t)=ejkω0t表示傅里叶基函数，ck(t)表示时频系数.Liu等(2011)考虑到傅里叶基函数的正交性建立了如下的最小平方优化问题求解时变傅里叶系数ĉk(t)：
|
(9) |
其中R是正则化算子，目的是约束时变傅里叶系数的平滑性.
假设地震CMP道集为d(t, x)，通过局部属性时频分析可以获得时频谱s(t, x, f)，在CMP道集时频谱中利用预测映射技术将目的层顶底反射自动追踪出来，然后将相同射线参数的反射频谱进行谱比反演求取Q.为了利用多射线信息同时约束反演，我们将公式(1)写成如下最小平方优化问题：
|
(10) |
其中ȃ和
|
(11) |
其中M是采用的射线总个数，N是选用的频率采样总个数.在实际计算过程中，一般选用一定频带范围内的频率采样和一定范围内的射线，频带选取对Q值估计非常重要，在讨论部分会对其进一步测试分析.为了获得更稳定和鲁棒的反演结果，本文使用正则化技术约束ȃ和
|
(12) |
其中λ1和λ2是分别对应ȃ和
|
(13) |
简记为
|
(14) |
其中FT代表矩阵共轭转置，N有如下的形式：
|
(15) |
在公式(10)—(15)所示的叠前地震数据Q值估计中，考虑了所有的传播方向和所有频率分量的地震波的衰减情况，而叠后地震数据Q值估计没有利用所有的传播方向的信息，动校正叠加处理实际上将衰减特征不同的波形进行了叠加，影响了Q值估计效果.采用叠前地震数据可以充分考虑波传播方向上的衰减特征来估计Q值，因此更符合实际波场传播特征，估计的Q值也更准确.
综上，基于多射线联合反演的速度无关叠前Q值估计方法计算过程如下：
(1) 通过平面波预测法(Plane-Wave Destruction，PWD)算法(Fomel，2002)获得CMP道集的局部斜率p(t, x)，利用局部属性时频分析方法得到CMP道集的局部时频振幅谱s(t, x, f).
(2) 从零偏移距剖面中选取目的层顶底走时(实际应用中利用解释层位或自动选取比较强的反射)作为种子点用于预测映射技术，获得该零偏移距的走时曲线φ(t0, x).
(3) 在叠前走时曲线φ(t0, x)上找到具有相同射线参数的反射频谱对，如图 1中所示的M和N点对，然后从时频谱s(t, x, f)中提取它们的振幅谱St2(f, p)和St1(f, p).
(4) 将步骤(3)获得的多射线振幅谱St2(f, p)和St1(f, p)应用于公式(12)—(15)联合反演估计ȃ=1/Q.
(5) 对所有的层位循环进行上述(1)—(4)步，获得不同层位的Q值，如果想要得到全部地震数据的Q值剖面，可以通过层间插值和平滑建立整个剖面对应的Q值模型.
2 模拟CMP道集算例我们首先用一个模拟的CMP道集验证本文方法，图 2显示了模型参数(速度、密度和Q值)及合成的CMP记录，随深度增加层速度从1300 m·s-1增大到2200 m·s-1，品质因子Q从10增大到100，震源子波采用40 Hz的Ricker子波，在合成的CMP道集中添加了少量的随机噪声.图 3a和3b分别是通过PWD算法估计的局部斜率及预测映射算法获得的零偏移距走时剖面(非零炮检距处反射对应的零炮检距时间)，可以发现局部斜率和零偏移距走时剖面是平滑的.为了说明预测映射技术对零炮检距走时的映射作用，我们选取了几个比较强的振幅作为种子点，采用预测映射技术对其进行拾取(图 3c)，可以发现预测映射技术将具有相同零炮检距走时的同相轴很好的追踪出来.图 4显示了同相轴局部斜率相等的线(蓝线)和零炮检距走时相等的线(黑线)，蓝线和黑线的交点对应着图 1中的M和N点，通过选取目的层顶底界面的零炮检距走时，然后通过预测映射获得等零炮检距线(黑线)，再通过与等局部斜率线(蓝色)交点获得反射点对，最后提取它们的时频谱参与谱比法联合反演Q值.注意对同一目的层(如图 4箭头所示)有很多类似于图 1中M，N点的这样的点对，这些点对对应的Q值都是这一目的层的Q值，因此利用多射线联合反演这一层的Q值，这一计算过程是不需要人工干预的全自动化过程.
|
图 2 模型参数及CMP合成记录 (a)速度；(b)密度；(c) Q值；(d)合成CMP道集. Fig. 2 The parameters of the model and synthetic CMP record (a) Velocity; (b) Density; (c) Q value; (d) Synthetic CMP gather. |
|
图 3 (a) PWD算法估算的局部斜率; (b)通过预测映射技术估计的零偏移距走时剖面; (c)通过预测绘画标记的数据 Fig. 3 (a) Local slopes estimated by PWD; (b) Zero-offset traveltimes by predictive mapping; (c) The marked events by predictive painting |
|
图 4 相同射线参数和相同零炮检距走时示意图 蓝线表示相同的局部斜率(射线参数)，黑线表示相同的零偏移距走时. Fig. 4 An illustration for the same ray parameter and the same zero-offset traveltime The blue lines represent the same local slopes (ray parameter) and the black lines represent the same zero-offset traveltimes. |
图 5是利用局部属性时频分析获得某一道的时频谱，可以看出局部属性方法能够反映地震数据的衰减特征，高频分量被严重衰减，主频从高频向低频移动.带宽的选取对谱比法很重要，因为带宽严重影响谱比法的稳定性，图 5中的黑线是选用的带宽范围，后文讨论部分将对带宽做进一步的探讨.为了分析本文方法的有效性和精度，我们对比了叠后谱比法(Tonn，1991)和Q随偏移距变化法(Q Variation with Offset，QVO)(Dasgupta and Clark, 1998)两种Q值估计的方法.叠后谱比法通过叠加提高了信噪比，但平均效应改变了波形信息(振幅和频率)，从而会影响Q值估计精度，尤其对于浅层数据由于动校正影响较大，Q值估计的误差也较大.QVO方法尽管利用了叠前的地震数据，但没有考虑射线的实际传播路径，且从等效Q转化为层Q过程中存在等效误差的影响.从结果和误差图上看，与其他方法相比，本文方法得到了更为准确的Q值估计结果(图 6)，表明本文方法可以有效稳定的从叠前地震数据中估计Q值.
|
图 5 模拟CMP道集(图 2d)中零炮检距记录道的时频谱 黑线表示所选带宽的范围. Fig. 5 The time-frequency spectra of zero-offset trace from synthetic CMP gather (Fig. 2d) The black lines indicate the range of the selected bandwidth. |
|
图 6 Q值估计的结果及误差 (a)本文方法；(c)叠后谱比法；(e)叠前QVO方法；(b)、(d)和(f)分别是(a)、(c)和(e)的误差.虚线是估计的Q值，实线是真实Q值. Fig. 6 The results of (a) the proposed method, (c) poststack spectral-ratio method, and (e) prestack QVO method. The subplots (b), (d) and (f) are the errors of (a), (c) and (e), respectively The dashed line indicates the estimated Q values and the solid line indicates the Q values used in modeling. |
首先我们将本文方法应用于信噪比较低的二维陆地地震数据Q值估计中，对实际数据经过简单的预处理，如抽道集、面波压制、随机噪声衰减等(Liu et al., 2015)，得到去噪后的CMP道集(图 7a).图 7b和图 7c分别显示了局部斜率和零偏移距走时剖面，图 8a显示了共偏移距道集的时频谱，可以看出振幅谱从浅层到深层高频分量迅速衰减(如图 8a箭头所示).图 8b显示了20 Hz的频率切片，发现深层20 Hz能量变弱，同样在远炮检距能量也变弱，表明随着地震波传播20 Hz的地震波能量被逐渐衰减.选取主要层位(反射振幅较强层)进行Q值估计，然后利用插值和平滑得到最终的Q剖面(图 9)，可以看出估计的Q位于50~180之间，不同的地层Q值有差异，注意到尽管在叠后地震剖面上(图 9a)有些同相轴信号信噪比低，但在估计Q过程中由于采用了叠前数据多射线联合反演，因此得到Q值还是比较稳定可靠的.
|
图 7 (a) 陆上地震数据CMP道集；(b)估计的局部斜率；(c)零偏移距走时 Fig. 7 (a) The CMP gathers of land seismic data; (b) The estimated local slopes; (c) Zero-offset traveltimes |
|
图 8 (a) 共偏移距道集的时频谱；(b) 20 Hz时频切片 Fig. 8 (a) Time-frequency spectra of one common offset gather; (b) The 20 Hz slice |
|
图 9 (a) 叠后地震数据；(b)本文方法估计Q值剖面 Fig. 9 (a) The seismic image after stacking; (b) The final Q section estimated by the proposed method |
然后，将其应用于信噪比较高的海洋地震数据Q值估计中，图 10显示了海洋地震勘探某测线的CMP道集，虚线是两口VSP井的位置CMP 808和CMP 1572，这两口井提供了VSP数据估计Q值的结果(Keys and Foster, 1998).图 11显示CMP 1572处叠前零炮检距道和叠后道的时频谱，通过对比发现叠后时频谱信噪比高，但与叠前原始数据的频谱存在一些差异.我们每隔15个CMP进行一次Q值估计，然后把CMP估计结果利用插值得到整体剖面的Q值.图 12显示了其中四个CMP(808，1195，1045和1572)的Q的结果，其中图 12a和图 12d是VSP井的位置.蓝线表示我们选取的几个层位用本文方法估计的Q值，经过CMP间的横向插值和平滑后的结果如绿线所示，红色点是VSP井估计的Q值(VSP数据只有部分结果)，粉色线是叠后数据估计的Q值.通过比较可以看出，本文方法结果和VSP井方法基本接近，叠后谱比法结果在浅层(1 s处)Q值变化较大，主要是由于动校正和叠加引起的.最终整体剖面的Q值估计结果如图 13所示，叠前的估计的方法整体上连续性较好，叠后方法横向差异较大，因为叠前采用了多射线约束提高了Q值估计的可靠性.
|
图 10 海上地震数据CMP道集 虚线表示CMP 808和CMP 1572处VSP井的位置. Fig. 10 Marine field seismic CMP gathers The dashed lines indicate the locations of two VSP wells at CMP 808 and CMP 1572. |
|
图 11 (a) CMP 1572处叠前地震道的时频谱；(b) CMP 1572处叠后地震道的时频谱 Fig. 11 (a) Time-frequency spectrum of prestack seismic trace at CMP 1572; (b) Time-frequency spectrum of CMP 1572 from poststack seismic trace |
|
图 12 不同方法估计的Q值 (a) CMP 808；(b) CMP 1195；(c) CMP 1405；(d) CMP 1572.红点：从井资料(VSP)中得到Q值；蓝线：本文方法没有插值和平滑得到的Q值；绿线：本文方法经过插值和平滑后得到的Q值；粉线：从叠后数据中用谱比法得到Q值. Fig. 12 The estimated Q values using different methods (a) CMP 808; (b) CMP 1195; (c) CMP 1405 and (d) CMP 1572. The red points: Q obtained from well data (VSP data); the blue line: Q from the proposed method without interpolation and smooth; the green line: Q from the proposed method after interpolation and smooth; the pink line: Q obtained by spectra-ratio method from the poststack data. |
|
图 13 (a) 叠后地震数据；(b)本文方法估计的Q值剖面；(c)叠后谱比法估计的Q值剖面 Fig. 13 (a) The poststack seismic data; (b) The Q profile estimated by the proposed method; (c) The Q profile estimated by poststack spectra-ratio method |
在谱比法Q值估计中，地震噪声和频谱带宽对结果精度有很大的影响，下面针对噪声和带宽对本方法的影响进行讨论分析.为了分析噪声水平对所提出算法的影响，在模拟CMP道集(图 2)中加了不同方差水平的噪声进行测试分析(图 14).尽管PWD算法与其他数值微分的方法在抗噪性上更具有优越性，但随着噪声水平的增大，PWD估计的局部斜率精度也会降低.在PWD算法中采用高阶算子能够提高局部斜率估计的精度，但同时也降低了计算效率(Fomel，2002).与图 3a相比，较强的噪声导致局部斜率准确性降低，这意味着当有强噪声时，从CMP道集中找到具有相同射线参数的反射波形更困难，从而影响Q值估计的结果.
|
图 14 (a) 加方差为5×10－7的噪声后的CMP道集；(c)加方差为10－6的噪声后的CMP道集；(b)，(d)分别为由(a)，(c)估计的局部斜率 Fig. 14 (a) Noised CMP gather with noise variance 5×10－7; (b) The estimated local slope of (a); (c) Noised CMP gather with noise variance 10－6; (d) The estimated local slope of (c) |
图 15显示了不同噪声水平下的Q值估计结果.可以看出噪声越强，Q值估计的误差越大，当高斯分布噪声方差为5×10－7时(地震有效信号振幅的最大值是2×10－2)，估计的Q值结果与理论解还比较接近，但噪声增大后(方差为10－6)，一些层位的Q值估计结果误差就很大，因为有效信号被淹没在噪声中，频谱比结果变差，如0.6 s处层位Q结果误差很大就是因为顶界面反射太弱，被淹没在噪声中造成的.因此，在应用本方法时建议首先做好噪声压制处理，提高叠前数据的信噪比，然后再采用本文方法来估计Q值，提高Q值估计的精度.
|
图 15 (a) 由图 14a CMP道集估计的Q值(虚线)和真实Q值(实线)；(b)图(a)估计Q的误差；(c)由图 14c CMP道集估计的Q值(虚线)和真实Q值(实线)；(d)图(c)估计Q的误差 Fig. 15 (a) The estimated Q (dashed line) and real Q (solid line) for the CMP gather shown in Fig. 14a; (b) The error of (a); (c) The estimated Q (dashed line) and real Q (solid line) for the CMP gather shown in Fig. 14c; (d) The error of (c) |
为了测试带宽对Q估计的影响，我们以图 2模型中Q=30的目的层为例，分析了选用不同带宽下Q估计的结果.设置带宽最低限从5 Hz到20 Hz，最高限分别从21 Hz到60 Hz，分别进行多射线联合反演Q值估计.图 16a显示了不同带宽估计的1/Q结果，横坐标代表带宽的低限值，纵坐标代表带宽高限值，颜色代表 1/Q值.可以看出，完全用低频的带宽时(如14 Hz到22 Hz)结果很差，甚至出现了负数，原因是分布在这一频段内的频谱能量很弱.我们统计了所有不同的带宽下估计的1/Q结果，采用直方图形式显示在图 16b中，可以看出所有带宽选择中，1/Q=0.033的概率是最高的，说明本方法对不同带宽的选择还是比较稳定的.在实际应用中，为了降低带宽的影响可以对所有的带宽选择进行处理，最终保留直方图中概率最大的1/Q值，这样可以降低带宽对算法的影响，但同时也增加计算量，因为对每一个带宽就需要进行一次反演.
|
图 16 (a) 不同带宽估计的1/Q；(b)采用不同带宽估计的1/Q值概率分布柱状图 Fig. 16 (a) The result of 1/Q with different bandwidths; (b) Histogram of 1/Q estimated with different bandwidths |
另外，本方法假设Q值是与频率无关的，如果考虑Q值频变的情况，可以采用类似于Beckwith等(2017)的方法，将本文方法扩展到频变的情况，注意Beckwith等(2017)采用的频变Q模型是非线性的，因此联合Q反演变为非线性问题，数值求解过程更复杂.本方法是在CMP道集上做的Q值估计，假设界面是局部水平的，如果界面倾角较大，需要进一步的研究.
5 结论本文针对叠前地震数据，将局部斜率和预测映射技术引入到谱比法Q值估计中，提出了一种更稳定鲁棒的叠前地震资料Q值估计方法.利用局部斜率信息确定多射线路径的反射波形，通过预测映射技术将非零偏移距走时映射到零炮检距走时，通过局部时频变换计算时频谱，将不同射线参数相同层位的反射频谱联合反演求解Q值.该方法是一种速度无关的Q值自动估计方法，避免了波场传播直射线假设，减小了速度分析误差和人工拾取叠前走时曲线对Q值结果估计的影响.在联合反演中采用了多射线地震反射信息，对不同物理参数引入不同正则化参数，提高了反演结果的精度和稳定性.模拟算例和实际算例表明，与常规方法相比，本方法可以获得更精确、更稳定的Q值估计结果.本文假设衰减是各向同性和频率无关的，将本方法扩展到各向异性和频变的情况是我们后续研究的重要内容之一.
致谢 感谢Sergey Fomel教授关于平面波预测滤波器和预测映射算法方面的讨论.
Bakker P, van Vliet L J, Verbeek P W. 1999. Edge preserving orientation adaptive filtering.//Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition. Fort Collins, CO, USA: IEEE, 535-540. https://www.researchgate.net/publication/3813449_Edge_preserving_orientation_adaptive_filtering
|
Bakker P. 2002. Image structure analysis for seismic interpretation[Ph. D. thesis]. Delft, Netherlands: Delft University of Technology.
|
Båth M. 1974. Developments in Solid Earth Geophysics. Amsterdam:Elsevier Science Publishing Co. |
Beckwith J, Clark R, Hodgson L. 2017. Estimating frequency-dependent attenuation quality factor values from prestack surface seismic data. Geophysics, 82(1): O11-O22. |
Behura J, Tsvankin I. 2009. Estimation of interval anisotropic attenuation from reflection data. Geophysics, 74(6): A69-A74. |
Blias E. 2012. Accurate interval Q-factor estimation from VSP data. Geophysics, 77(3): WA149-WA156. |
Burnett W, Fomel S. 2009. 3D velocity-independent elliptically anisotropic moveout correction. Geophysics, 74(5): WB129-WB136. |
Carter A J, Clark R A, Neivill P C, et al. 2001. Attenuation measurements from surface seismic data-Azimuthal variation and time-lapse case studies.//63rd Conference and Technical Exhibition. EAGE.
|
Dasgupta R, Clark R A. 1998. Estimation of Q from surface seismic reflection data. Geophysics, 63(6): 2120-2128. |
Du J, Wang S X, Liu G C, et al. 2009. VSP wavefield separation using local slopes attribute. Chinese Journal of Geophysics (in Chinese), 52(7): 1867-1872. |
Fomel S. 2002. Applications of plane-wave destruction filters. Geophysics, 67(6): 1946-1960. |
Fomel S. 2007. Velocity-independent time-domain seismic imaging using local event slopes. Geophysics, 72(3): S139-S147. |
Fomel S. 2010. Predictive painting of 3D seismic volumes. Geophysics, 75(4): A25-A30. |
Futterman W I. 1962. Dispersive body waves. Journal of Geophysical Research, 67(13): 5279-5291. |
Hale D. 2009. Structure-oriented smoothing and semblance. CWP Report 635, Colorado School of Mines. |
Hao Y J, Wen X T, Zhang B, et al. 2016. Q estimation of seismic data using the generalized S-transform. Journal of Applied Geophysics, 135: 122-134. |
Keys R G, Foster D J. 1998. Comparison of seismic inversion methods on a single real data set. Society of Exploration Geophysicists. |
Khoshnavaz M J, Bóna A, Urosevic M. 2016. Velocity-independent estimation of kinematic attributes in vertical transverse isotropy media using local slopes and predictive painting. Geophysics, 81(5): U73-U85. |
Kneib G, Shapiro S A. 1995. Viscoacoustic wave propagation in 2-D random media and separation of absorption and scattering attenuation. Geophysics, 60(2): 459-467. |
Liu G C, Fomel S, Chen X H. 2011. Time-frequency analysis of seismic data using local attributes. Geophysics, 76(6): P23-P34. |
Liu G C, Chen X H, Rao Y. 2018a. Seismic quality factor estimation using frequency-dependent linear fitting. Journal of Applied Geophysics, 156: 1-8. |
Liu G C, Li C, Liu X Y, et al. 2018b. Automatic stacking-velocity estimation using similarity-weighted clustering. Geophysical Prospecting, 66(4): 649-663. |
Liu X H, Liu G C, Cui Y Q, et al. 2016. Q estimation with lateral constraint based on local time-frequency transform. Oil Geophysical Prospecting (in Chinese), 51(5): 868-874. |
Liu Y, Fomel S. 2013. Seismic data analysis using local time-frequency decomposition. Geophysical Prospecting, 61(3): 516-525. |
Liu Y, Fomel S, Liu C. 2015. Signal and noise separation in prestack seismic data using velocity-dependent seislet transform. Geophysics, 80(6): WD117-WD128. |
Liu Y, Li B X, Liu C, et al. 2016. Attenuation compensation method of seismic wave in the local time-frequency transform domain. Journal of Jilin University (Earth Science Edition) (in Chinese), 46(2): 594-602. |
Marfurt K J. 2006. Robust estimates of 3D reflector dip and azimuth. Geophysics, 71(4): P29-P40. |
Maultzsch S, Chapman M, Liu E R, et al. 2007. Modelling and analysis of attenuation anisotropy in multi-azimuth VSP data from the Clair field. Geophysical Prospecting, 55(5): 627-642. |
Müller T M, Gurevich B, Lebedev M. 2010. Seismic wave attenuation and dispersion resulting from wave-induced flow in porous rocks-A review. Geophysics, 75(5): 75A147-75A164. |
Oliveira F de S, de Figueiredo J J S, Oliveira A G, et al. 2017. Estimation of quality factor based on peak frequency-shift method and redatuming operator:Application in real data set. Geophysics, 82(1): N1-N12. |
Quan Y L, Harris J M. 1997. Seismic attenuation tomography using the frequency shift method. Geophysics, 62(3): 895-905. |
Reine C, van der Baan M, Clark R. 2009. The robustness of seismic attenuation measurements using fixed-and variable-window time-frequency transforms. Geophysics, 74(2): WA123-WA135. |
Reine C, Clark R, van der Baan M. 2012. Robust prestack Q-determination using surface seismic data:Part 1-Method and synthetic examples. Geophysics, 77(1): R45-R56. |
Ricker N. 1953. The form and laws of propagation of seismic wavelets. Geophysics, 18(1): 10-40. |
Schleicher J, Costa J C, Santos L T, et al. 2009. On the estimation of local slopes. Geophysics, 74(4): P25-P33. |
Sun S Z, Sun X K, Wang Y G, et al. 2015. Q estimation using modified S-transform based on pre-stack gathers and its applications on carbonate reservoir. Journal of Geophysics and Engineering, 12(5): 725-733. |
Tonn R. 1991. The determination of the seismic quality factor Q from VSP data:A comparison of different computational methods. Geophysical Prospecting, 39(1): 1-27. |
Wang S D, Yang D F, Li J N, et al. 2015. Q factor estimation based on the method of logarithmic spectral area difference. Geophysics, 80(6): V157-V171. |
Wang Y H. 2004. Q analysis on reflection seismic data. Geophysical Research Letters, 31(17): 159-180. DOI:10.1029/2004GL020572 |
White J E. 1983. Underground Sound:Application of Seismic Waves. Amsterdam:Elsevier Science Publications, 83-137. |
Winkler K W, Nur A. 1982. Seismic attenuation:Effects of pore fluids and frictional-sliding. Geophysics, 47(1): 1-15. |
Wu X M, Janson X. 2017. Directional structure tensors in estimating seismic structural and stratigraphic orientations. Geophysical Journal International, 210(1): 534-548. |
Wu X M, Fomel S. 2018. Least-squares horizons with local slopes and multigrid correlations. Geophysics, 83(4): IM29-IM40. |
Zhang Z J, Ulrych T J. 2002. Estimation of quality factors from CMP records. Geophysics, 67(5): 1542-1547. |
杜婧, 王尚旭, 刘国昌, 等. 2009. 基于局部斜率属性的VSP波场分离研究. 地球物理学报, 52(7): 1867-1872. |
刘喜恒, 刘国昌, 崔永谦, 等. 2016. 局部时频变换横向约束Q值估计方法. 石油地球物理勘探, 51(5): 868-874. |
刘洋, 李炳秀, 刘财, 等. 2016. 局部时频变换域地震波吸收衰减补偿方法. 吉林大学学报(地球科学版), 46(2): 594-602. |