石油物探  2018, Vol. 57 Issue (2): 222-230  DOI: 10.3969/j.issn.1000-1441.2018.02.007
0
文章快速检索     高级检索

引用本文 

金子奇, 孙赞东. 改进的衰减旅行时层析方法估计Q值[J]. 石油物探, 2018, 57(2): 222-230. DOI: 10.3969/j.issn.1000-1441.2018.02.007.
JIN Ziqi, SUN Zandong. Improved attenuated traveltime tomography for Q estimation[J]. Geophysical Prospecting for Petroleum, 2018, 57(2): 222-230. DOI: 10.3969/j.issn.1000-1441.2018.02.007.

作者简介

金子奇(1988—), 男, 博士在读, 主要从事地震数据处理研究。Email:13945920057@126.com

文章历史

收稿日期:2017-08-22
改回日期:2017-12-12
改进的衰减旅行时层析方法估计Q
金子奇, 孙赞东     
中国石油大学(北京), 北京 102249
摘要:衰减旅行时层析作为一种估计Q值的方法, 其衰减旅行时计算的准确程度影响最终Q值反演结果的精度。为了提高衰减旅行时层析求解精度而得到更真实的地下Q值分布, 提出了一种改进的衰减旅行时层析方法, 提高了模型中的衰减旅行时计算精度和实际资料中的衰减旅行时拾取精度。在射线追踪计算模型中的衰减旅行时的过程中, 用方向速度梯度代替传统的速度梯度, 考虑不同射线在不同传播方向上的速度差异, 提高了射线追踪计算旅行时和反射角的精度, 从而提高衰减旅行时计算精度。从实际资料中拾取衰减旅行时的时候, 采用对数谱比反演方法, 利用多道反射波信息同时反演多道衰减旅行时, 从而避免了传统方法中人为选择最优频带带来的误差, 使反演结果更稳定。在反演过程中考虑反射波在不同界面内的旅行时之差, 可以有效解决上覆地层的影响。模型数据测试和实际资料应用结果表明, 该方法Q值估算结果更加准确, Q补偿偏移后的剖面更加符合实际地下介质情况, 验证了方法的可行性和有效性。
关键词品质因子    Q值估计    衰减旅行时层析    射线追踪    Q补偿偏移    
Improved attenuated traveltime tomography for Q estimation
JIN Ziqi, SUN Zandong     
China University of Petroleum, Beijing 102249, China
Abstract: Attenuated traveltime tomography is an efficient method for Q estimation, and the calculation accuracy of attenuated traveltimes is essential.An improved attenuation traveltime tomography method is proposed.The proposed method improves the calculation accuracy of attenuated traveltimes for model data, and the picking accuracy of attenuated traveltimes for field data.Considering various properties of seismic rays in different directions, the method defines the direction gradient in ray-tracing instead of the conventional velocity gradient.The calculated reflection angles and ray paths are more accurate, thus more accurate Q results could be obtained.Log spectral ratio inversion is designed for field data attenuation traveltime calculation.Reflections of multi offset traces are used for attenuated traveltime estimations of multi offset traces simultaneously.This method could avoid artificial error by manually selecting optimistic frequency band, and so can derive more robust results.In addition, the method can avoid the effect from overburden by considering the traveltime difference between different layers.Synthetic and field data tests demonstrated that Q estimation using the proposed method is more accurate, and that the migration profile after Q compensation is more reliable, compared with the conventional method.
Key words: quality factor    Q estimation    attenuated traveltime tomography    ray-tracing    Q compensation    

根据应用空间的不同, 求取层间Q值的方法可以分成两类, 即时间域Q值估计方法和频率域Q值估计方法。时间域Q值估计方法很难区分吸收衰减和其它类型的衰减, 不建议使用。依据固有衰减作用与频率相关的特点, 频率域Q值估计方法能够很好地解决这个问题。通过拾取不同时间位置的两个子波并进行傅里叶变换, 频率域Q值估计方法能够根据子波振幅谱的变化来估计Q值。DASGUPTA等[1]首次利用谱比法估计Q值, 在叠前反射波道集中通过建立子波对数谱和Q值间的关系估计Q值。随后, 许多学者对其做了改进, 陈文爽等[2]、王小杰等[3]、王德利等[4]以S变换的频谱作为谱比法的输入, 避免了通常计算地层吸收参数时的平均效应。张繁昌等[5]将自适应分解得到的子波应用谱比法解决薄层和反射波干涉问题。李伟娜等[6]在微测井记录中应用谱比法并借鉴QVO(Q versus offset)的思想以适应Q值横向变化的情况。QUAN等[7]在VSP资料中, 最先利用质心频率位移量估算Q值, 建立了衰减后的子波峰值频移量与Q值之间的关系。与质心频移法类似, ZHANG等[8]在叠前CMP道集中采用峰值频率法估计Q值。随后, 许多学者对上述两类频移方法做了改进, 王宗俊[9]利用谱模拟的子波对质心频移法进行约束; 张立彬等[10]引入积分中值参变量避免了质心频移法中对震源子波的先验假设; 高静怀等[11]利用特征结构法提高地震子波峰值频率的估计精度。

然而, 许多因素影响频域Q值估计方法的精度。震源强度和带宽会影响Q值估计结果, 在主频范围之外估计得到的Q值会变得不稳定[12]Q值估计结果随着信噪比的降低而显著畸变。对于一些频率成分, 远偏移距检波器记录的振幅有时会比近偏移距检波器记录的振幅强, 这种效应会引起Q值估计结果的误差[12], 最严重的缺点是这些方法计算得到的Q值与路径有关[13-14], Q值估计结果随观测系统的变化而变化。而真实地下介质Q值分布应仅与构造、岩性和流体有关, 与观测系统无关。另外, 频域Q值估计方法得到的是大套地层层间的等效Q值, 无论是纵向上还是横向上, 分辨率都较低。BRZOSTOWSKI等[15]建立了不同传播距离地震波振幅的变化与Q值之间的关系, 采用联合迭代重建技术SIRT求解方程实现振幅衰减Q值层析。振幅衰减层析方法在时间域实现, 由于地震波振幅受多种因素(几何扩散、仪器响应、震源/检波器耦合特性、反射/透射等)的影响, 因此其稳定性较差。LIAO等[16]和严又生等[17]利用质心频移法结合层析反演方法计算Q值。该方法利用波传播路径下Q值与质心位移量的对应关系来迭代反演速度和Q值信息, 对绝对振幅信息不敏感, 只适应较窄频带范围内的数据。

天然地震学中, NOWACK等[18]将衰减旅行时应用于地震折射波数据, 反演浅层地壳内部衰减及速度结构。樊计昌等[19]利用三维Q值层析成像方法对岫岩坑的地震波衰减结构进行了研究。随后, 勘探地震学中也开始利用衰减旅行时层析方法反演地层Q值, CAVALCA等[20]、XIN等[21]和JIN[22]将衰减层析方法应用于地震反射波数据, 利用对数谱比信息拟合衰减旅行时, 并利用非线性反演方法计算地下真实Q值分布。衰减旅行时层析方法反演地下真实Q值分布, 能适应Q值剖面纵横向变化, 且对绝对振幅信息不敏感。衰减旅行时层析法, 分别在实际资料和初始Q值模型中, 求取实际衰减旅行时和模型衰减旅行时, 建立二者之差与Q值模型更新量之间的关系, 通过非线性反演的方法, 计算Q模型更新量, 更新初始模型直到接近真实地下Q值分布[22]。因此, 衰减旅行时的计算是该方法的核心。传统方法在计算实际衰减旅行时的时候, 需要人工拾取最优频带, 人为误差对Q值反演结果造成很大影响, 并且, 模型衰减旅行时的计算精度受射线追踪精度影响, 有偏差的旅行时和不合理的传播路径会导致反演结果产生误差。魏文等[23]和王小杰等[24]采用时频分析结合谱比法解决人为拾取最优频带的问题, 但该方法仍属于频率域Q值估计方法, 得到的是大套地层的等效Q值。

因此, 为了得到更为可靠和准确的Q值估计结果, 本文采用改进后的衰减旅行时层析方法反演地下Q值分布。该方法在实际资料中拾取衰减旅行时时, 采用对数谱比反演方法, 利用多道反射波信息同时反演多道衰减旅行时, 避免了传统方法中人为选择最优频带可能引入的误差, 反演结果更稳定。在计算模型衰减旅行时过程中, 为提高衰减旅行时计算精度, 采用方向梯度射线追踪方法计算旅行时。改进的算法用方向速度梯度代替传统的速度梯度, 考虑了不同射线在不同传播方向上的速度差异, 得到更加准确的旅行时和传播路径。而且, 通过在反演过程中考虑反射波在不同界面内的旅行时之差, 可以有效解决上覆地层的影响。改进后的衰减旅行时层析方法, 经模型数据和实际资料测试, Q值反演结果更为准确, 为后续吸收补偿和成像提供可靠的衰减信息。

1 衰减旅行时层析反演算法

在地震勘探中, 震源激发的地震波经过地下介质传播, 经反射层反射回地面, 再由地面检波器接收。地震波本身携带了振幅和旅行时信息, 即地震波的动力学信息和运动学信息。我们利用旅行时信息进行层析反演求取地下Q值分布。衰减旅行时层析反演求取Q值流程如图 1所示, 具体算法可以分为以下几步:①对数频谱法反演方法计算炮集记录上的反射波衰减旅行时; ②建立初始Q值模型, 利用方向梯度射线追踪模拟地震波传播, 计算模型中的反射波衰减旅行时; ③建立实际衰减旅行时与模型衰减旅行时之差与Q值模型更新量的关系, 采用联合迭代重建方法求解大型稀疏矩阵, 更新Q值模型。

图 1 衰减旅行时层析反演求取Q值流程

图 2为单层二维介质模型, 由震源激发的射线经反射后到达地表被接收。衰减旅行时层析方法, 首先将地下介质划分成网格(图 2), 并在每一网格内定义独立的速度和Q值, 速度和Q值均为常量。一条由震源激发经模型底界面反射并由地表检波器接收的射线, 在网格中被分为多段。

图 2 衰减旅行时层析反演中的网格划分示意

k条射线的衰减旅行时定义为[22]:

$ t_k^ * = \sum\limits_k {\frac{{{l_{kl}}}}{{{v_{kl}}}}Q_l^{ - 1}} = \sum\limits_k {{t_{kl}}Q_l^{ - 1}} $ (1)

式中:lkl, tkl, vkl分别为第k条射线在第l个网格中的传播距离、传播时间和传播速度。

将方程(1)改写为矩阵形式:

$ \mathit{\boldsymbol{TQ}} = {\mathit{\boldsymbol{t}}^ * } $ (2)

式中:矩阵T包含由射线追踪得到的每一网格内的传播时间tkl; Q为包含Q值的向量; t*向量为衰减旅行时。

选择试验点利用谱比法计算Q值, 插值后得到工区内的初始Q值模型Qm。通常, 初始Q值模型不准确, 初始模型与真实值之差记做δQl-1, 即模型更新量。建立模型更新量与衰减旅行时之间的关系如下:

$ \delta \mathit{\boldsymbol{t}}_k^ * = t_{\rm{r}}^ * - t_{\rm{m}}^ * = \sum {{t_{kl}}{\rm{ \mathit{ δ} }}Q_l^{ - 1}} $ (3)

式中:tr*为实际资料中计算的衰减旅行时; tm*为模型中射线追踪得到的衰减旅行时。将公式(3)表示为完整的矩阵形式:

$ \left[ {\begin{array}{*{20}{c}} {{t_{11}}}& \cdots &{{t_{1n}}}\\ \vdots &{}& \vdots \\ {{t_{m1}}}& \cdots &{{t_{mn}}} \end{array}} \right] \times \left[ {\begin{array}{*{20}{c}} {\delta Q_1^{ - 1}}\\ \vdots \\ {\delta Q_n^{ - 1}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\rm{ \mathit{ δ} }}t_1^ * }\\ \vdots \\ {{\rm{ \mathit{ δ} }}t_m^ * } \end{array}} \right] $ (4)

公式(4)为m个方程n个未知数组成的大型稀疏矩阵, 应用正交分解最小二乘法(LSQR)求解该方程组。最终Q值结果Qf由初始模型和模型更新量共同表示:

$ {Q_{\rm{f}}} = \frac{{{Q_{\rm{m}}}}}{{1 + {Q_{\rm{m}}} \cdot {\rm{ \mathit{ δ} }}{Q^{ - 1}}}} $ (5)
2 方向梯度射线追踪计算衰减旅行时

由公式(3)可知, 衰减旅行时层析方法涉及到两种类型的衰减旅行时, 一种是从实际资料中拾取的衰减旅行时tr*, 另一种是由射线追踪得到的模型中的衰减旅行时tm*

计算模型中的衰减旅行时tm*, 首先根据野外观测系统定义炮点和检波点位置, 射线追踪计算每条射线的传播时间。采用传统衰减旅行时层析方法的射线追踪算法计算射线在分界面处的折射时有一定误差, 由于不准确的旅行时信息会导致最终反演结果的误差, 因此我们提出一种方向梯度射线追踪方法。在射线的传播过程中考虑不同射线在不同传播方向上的速度差异, 定义方向速度梯度。采用改进后的算法求取旅行时精度得到提高。

传统算法中, 利用中心差分公式计算网格内的速度梯度:

$ {\mathit{\boldsymbol{\lambda }}_x} = \frac{{{v_{i + 1}} - {v_{i - 1}}}}{{2{\rm{d}}x}} $ (6a)
$ {\mathit{\boldsymbol{\lambda }}_y} = \frac{{{v_{j + 1}} - {v_{j - 1}}}}{{2{\rm{d}}y}} $ (6b)

式中:下标i, j分别表示在x方向和y方向的网格索引; v为定义的网格内速度参数。

当射线穿过分界面时, 由于不准确的速度梯度定义, 传统算法计算得到的射线方向会偏离真实射线路径, 尤其是当网格间的速度变化超过10%时, 偏离误差更加明显, 使得传播时间的计算也随之偏离真实值。射线传播距离s后的传播时间定义为:

$ \begin{array}{*{20}{c}} {t\left( s \right) = \frac{s}{{{v_0}}}\left\{ {1 + \frac{{{s^2}}}{{6v_0^2}}\left[ {{\mathit{\boldsymbol{\lambda }}^2} + {{\left( {\mathit{\boldsymbol{\lambda }} \cdot {\mathit{\boldsymbol{n}}_0}} \right)}^2}} \right] - } \right.}\\ {\left. {\left( {\mathit{\boldsymbol{\lambda }} \cdot {\mathit{\boldsymbol{n}}_0}} \right)\frac{s}{{2{v_0}}}} \right\} + o\left( {{\mathit{\boldsymbol{\lambda }}^3}} \right)} \end{array} $ (7)

式中:向量n0为射线出射位置s=0处单位向量, 指示射线传播方向; 向量λ=(λx, λy)为沿射线方向的速度梯度; v0为原点处的速度; o(λ3)为包含了λ的三阶及更高阶项。将旅行时结果代入公式(1)可得衰减旅行时, 旅行时的误差最终导致衰减旅行时的计算误差。

改变传统算法中对速度梯度的定义, 考虑速度随传播方向的变化, 用方向速度梯度代替公式(7)中的速度梯度。图 3展示了4种不同的射线传播方向, 重新定义速度梯度如下。

图 3 不同射线传播方向示意 a 射线沿西北方向传播; b 射线沿东北方向传播; c 射线沿东南方向传播; d 射线沿西南方向传播

1) 射线沿西北方向传播时(图 3a):

$ {\mathit{\boldsymbol{\lambda }}_x} = \frac{{{v_{i - 1}} - {v_i}}}{{{\rm{d}}x}} $ (8a)
$ {\mathit{\boldsymbol{\lambda }}_y} = \frac{{{v_{j - 1}} - {v_j}}}{{{\rm{d}}y}} $ (8b)

2) 射线沿东北方向传播时(图 3b):

$ {\mathit{\boldsymbol{\lambda }}_x} = \frac{{{v_{i + 1}} - {v_i}}}{{{\rm{d}}x}} $ (9a)
$ {\mathit{\boldsymbol{\lambda }}_y} = \frac{{{v_{j - 1}} - {v_j}}}{{{\rm{d}}y}} $ (9b)

3) 射线沿东南方向传播时(图 3c):

$ {\mathit{\boldsymbol{\lambda }}_x} = \frac{{{v_{i + 1}} - {v_i}}}{{{\rm{d}}x}} $ (10a)
$ {\mathit{\boldsymbol{\lambda }}_y} = \frac{{{v_{j + 1}} - {v_j}}}{{{\rm{d}}y}} $ (10b)

4) 射线沿西南方向传播时(图 3d):

$ {\mathit{\boldsymbol{\lambda }}_x} = \frac{{{v_{i - 1}} - {v_i}}}{{{\rm{d}}x}} $ (11a)
$ {\mathit{\boldsymbol{\lambda }}_y} = \frac{{{v_{j + 1}} - {v_j}}}{{{\rm{d}}y}} $ (11b)

根据对不同射线出射方向的方向梯度定义, 利用公式(7)计算网格内射线旅行时。

3 利用对数谱信息反演衰减旅行时

传统层析方法在实际资料中拾取衰减旅行时tr*, 采用的是谱比法[25]。考虑通常震源子波未知, 拾取不同传播时间位置的反射波, 谱比法计算公式可以表示为:

$ \begin{array}{l} b = \ln \left[ {\frac{{A\left( {{t_2},f} \right)}}{{A\left( {{t_1},f} \right)}}} \right] = \ln G - \frac{{{\rm{ \mathit{ π} }}f\left( {{t_2} - {t_1}} \right)}}{Q}\\ \;\;\; = B - {\rm{ \mathit{ π} }}f\Delta {t^ * } \end{array} $ (12)

式中:f为地震信号频率; A(t1, f), A(t2, f)分别为地震波在传播时间t1t2处的振幅谱; G为频率独立因子, 代表与频率无关的振幅衰减; B为独立因子项; Δt*为两个反射波衰减旅行时之差。传统层析方法通过拟合对数谱比与频率得到斜率, 建立斜率与衰减旅行时的关系式, 求得衰减旅行时。拟合斜率时, 由于无法准确选取最优拟合频带, 且无法消除噪声等影响, 通常会引入人为误差。为解决上述问题, 在共炮点道集中采用多道反射波同时计算多道衰减旅行时。同时反演频率独立因子GQ值, 可以压制非吸收衰减作用(如球面扩散等)的干扰。而且, 传统方法得到的是衰减旅行时之差。由于拾取不同传播时间位置的反射波, 若拾取的子波在上覆地层具有不同的传播路径, 上覆地层内的衰减影响则不能相互抵消, 从而影响Q值估计结果。在反演过程中考虑反射波在不同界面内的旅行时之差, 可以有效解决上覆地层的影响。

以两层介质模型为例, 沿偏移距方向, 拾取t1时刻多个反射波A1, A3, …, A2n-1t2时刻多个反射波A2, A4, …, A2n, 如图 4所示。

图 4 两层介质模型内不同界面反射波射线路径示意(这些反射波均来自同一炮集记录)

拾取的反射波经傅里叶变换得到其振幅谱, 整理并分别沿纵向和横向取对数可得:

$ \left\{ \begin{array}{l} {b_1} = \ln \left( {\frac{{{A_2}}}{{{A_1}}}} \right) = {B_1} - {\rm{ \mathit{ π} }}f\left( {t_2^ * - t_1^ * } \right)\\ {b_2} = \ln \left( {\frac{{{A_3}}}{{{A_1}}}} \right) = {B_2} - {\rm{ \mathit{ π} }}f\left( {t_3^ * - t_1^ * } \right)\\ {b_3} = \ln \left( {\frac{{{A_4}}}{{{A_3}}}} \right) = {B_3} - {\rm{ \mathit{ π} }}f\left( {t_4^ * - t_3^ * } \right)\\ {b_4} = \ln \left( {\frac{{{A_5}}}{{{A_3}}}} \right) = {B_4} - {\rm{ \mathit{ π} }}f\left( {t_5^ * - t_3^ * } \right)\\ {b_{2n - 1}} = \ln \left( {\frac{{{A_{2n}}}}{{{A_{2n - 1}}}}} \right) = {B_{2n - 1}} - {\rm{ \mathit{ π} }}f\left( {t_{2n}^ * - t_{2n - 1}^ * } \right)\\ {b_{2n}} = \ln \left( {\frac{{{A_{2n + 1}}}}{{{A_{2n - 1}}}}} \right) = {B_{2n}} - {\rm{ \mathit{ π} }}f\left( {t_{2n + 1}^ * - t_{2n - 1}^ * } \right) \end{array} \right. $ (13)

其中, b1, b2, …, b2n中包含对数谱比信息。将公式(13)表示为d=Fm的矩阵形式:

$ \left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{d}}_1}}\\ {{\mathit{\boldsymbol{d}}_2}}\\ \vdots \\ {{\mathit{\boldsymbol{d}}_n}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} \mathit{\boldsymbol{v}}&\mathit{\boldsymbol{z}}& \cdots &\mathit{\boldsymbol{z}}&\mathit{\boldsymbol{z}}&{ - {\rm{ \mathsf{ π} }}\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}}&{ - {\rm{ \mathsf{ π} }}\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}}&\mathit{\boldsymbol{z}}& \cdots &\mathit{\boldsymbol{z}}&\mathit{\boldsymbol{z}}&\mathit{\boldsymbol{z}}\\ \mathit{\boldsymbol{z}}&\mathit{\boldsymbol{v}}& \cdots &\mathit{\boldsymbol{z}}&\mathit{\boldsymbol{z}}&{ - {\rm{ \mathsf{ π} }}\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}}&\mathit{\boldsymbol{z}}&{ - {\rm{ \mathsf{ π} }}\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}}& \cdots &\mathit{\boldsymbol{z}}&\mathit{\boldsymbol{z}}&\mathit{\boldsymbol{z}}\\ \vdots&\vdots &{}& \vdots&\vdots&\vdots&\vdots&\vdots &{}& \vdots&\vdots&\vdots \\ \mathit{\boldsymbol{z}}&\mathit{\boldsymbol{z}}& \cdots &\mathit{\boldsymbol{v}}&\mathit{\boldsymbol{z}}&\mathit{\boldsymbol{z}}&\mathit{\boldsymbol{z}}&\mathit{\boldsymbol{z}}& \cdots &{ - {\rm{ \mathsf{ π} }}\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}}&{ - {\rm{ \mathsf{ π} }}\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}}&\mathit{\boldsymbol{z}}\\ \mathit{\boldsymbol{z}}&\mathit{\boldsymbol{z}}& \cdots &\mathit{\boldsymbol{z}}&\mathit{\boldsymbol{v}}&\mathit{\boldsymbol{z}}&\mathit{\boldsymbol{z}}&\mathit{\boldsymbol{z}}& \cdots &{ - {\rm{ \mathsf{ π} }}\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}}&\mathit{\boldsymbol{z}}&{ - {\rm{ \mathsf{ π} }}\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{D}}_1}}\\ {{\mathit{\boldsymbol{D}}_2}}\\ \vdots \\ {{\mathit{\boldsymbol{D}}_{2n}}}\\ {\mathit{\boldsymbol{\tau }}_1^ * }\\ {\mathit{\boldsymbol{\tau }}_2^ * }\\ \vdots \\ {\mathit{\boldsymbol{\tau }}_{2n}^ * } \end{array}} \right) $ (14)

式中:d1, d2, …, d2n分别包含不同频率对应的对数谱比信息b1, b2, …, b2n; v=(1 1 … 1)T, Ω=(f1 f2fl)T, z=(0 0 … 0)T, 向量v, Ω, z的长度均为l, l为满足信噪比要求频带范围内的频率个数; 求解的向量m中, D1, D2, …, D2n分别包含l个独立因子项B1, B2, …, B2n; τ1, τ2, …, τ2n分别包含l个衰减旅行时t1*, t2*, …, t2n*。可以采用LSQR法求解公式(14)。最终提取的n个实际衰减旅行作为输入求解Q值模型更新量, 考虑上覆地层影响, 改变公式(3)可得:

$ \begin{array}{l} {\rm{ \mathit{ δ} }}t_{12}^ * = \left( {t_{{\rm{r1}}}^ * - t_{{\rm{r2}}}^ * } \right) - \left( {t_{{\rm{m1}}}^ * - t_{{\rm{m2}}}^ * } \right)\\ \;\;\;\;\;\; = \sum {\left( {{t_1} - {t_2}} \right){\rm{ \mathit{ δ} }}{Q^{ - 1}}} \end{array} $ (15)

式中:tr1*, tr2*分别为实际资料中不同分界面处反射波的衰减旅行时, tm1*, tm2*分别为对应模型中射线追踪得到的衰减旅行时。求解(15)式可得模型更新量, 代入公式(5)求得最终Q值分布。

4 模型数据测试

建立4层介质模型, 测试改进后的层析反演方法估计Q值的准确性, 模型参数如表 1所示。震源子波选择主频为60Hz的雷克子波, 射线追踪计算传播路径和旅行时, 并计算振幅衰减项和相位畸变项, 得到合成共炮点道集(为避免动校正拉伸作用对子波拾取的影响, 正演过程中未考虑正常时差的变化), 如图 5所示。

表 1 水平层状介质模型参数
图 5 水平层状模型合成炮集记录

以第1层和第2层为例, 说明改进后的射线追踪算法在射线路径计算精度上的改善。将模型划分网格, 每一网格大小为10m×10m, 模型尺度为80×40。对比改进算法与传统算法得到的折射角, 结果如表 2所示。图 6为两层介质模型下采用改进前、后的射线追踪算法计算得到的射线路径对比结果。从图 6中可以看出, 两种传播路径有明显差异, 射线传播至200m深处界面时发生折射, 且改进后的算法由于考虑方向梯度的变化, 得到的结果更接近真实射线路径, 提高了射线路径的计算精度, 尤其对于大角度入射的射线改进更加明显。与传统算法误差相比, 在小入射角时(小于20°), 改进算法得到的折射角误差小于传统算法误差的2%, 在大入射角时(大于50°), 误差之比也不会超过10%。

表 2 两层介质模型中改进前后射线追踪算法得到折射角对比
图 6 两层介质模型下改进前、后射线追踪算法计算得到的射线路径对比

随后, 采用公式(1)和公式(7)计算得到相应的旅行时及模型中射线追踪计算tm*。由于射线路径更加准确, 衰减旅行时也更准确。改进后射线追踪算法更适合建模和旅行时的计算。

选取一炮地震记录, 拾取来自不同界面的反射波同相轴, 分别采用多道同时反演法和传统方法计算衰减旅行时, 以第1层反射为例, 计算结果如图 7所示, 可以看出, 与传统方法相比, 改进后的方法更加稳定且计算结果更加精确。

图 7 衰减旅行时计算结果

得到衰减旅行时tr*后, 利用其与模型衰减旅行时tm*之差反演Q值, 图 8为采用改进的层析反演方法与传统层析反演方法得到的Q值结果(初始Q值设为100)。由于传统方法得到的衰减旅行时有较大误差, 使得反演计算的Q值偏离真实Q值, 尤其在深层, 即使是很小的衰减旅行时误差也导致最终Q值计算结果出现很大误差。从图 8可以看出, 改进的层析反演方法得到的Q值结果与真实值非常相近。

图 8 水平层状模型中Q值反演结果

改变模型第1层参数, 允许Q值横向变化, 目的是测试改进后的方法是否适应Q值变化的介质, 为了模型和计算上的简单性, 这里只考虑横向上8个有效单元(表 3)。同样利用对数谱比反演方法和方向梯度射线追踪方法计算实际衰减旅行时tr*和模型衰减旅行时tm*, 随后对Q值进行层析反演, 得到的结果如图 9所示。从图 9可以看出, 对于Q值横向变化的情况, 传统方法得到的Q值由于不精确的衰减旅行时, 存在不可忽略的误差; 而改进后的方法可以得到更加准确的Q值, 反演结果非常理想。

表 3 Q值横向变化介质模型参数
图 9 Q值横向变化模型中Q值反演结果
5 实际资料应用

将衰减旅行时层析法应用于实际资料求取Q值分布。我们将某一地区深度超过8000m的地层作为目标进行层析计算, 地震剖面中对应同向轴位置位于时间4.0s。抽取共炮点道集, 选择反射波同向轴明显的位置3.2s和4.0s计算衰减旅行时(图 10), 同时在初始Q值剖面上射线追踪得到相应的正演衰减旅行时信息。初始Q值模型中Q设置为常数100, 目标区域以80m×80m进行网格划分, 划分网格后的模型横向上包含80网格、纵向上包含80网格。最终反演Q值分布如图 11所示。

图 10 共炮点道集记录(拾取反射轴较明显位置作为Q值层析反演输入)
图 11 改进前(a)、后(b)衰减旅行时Q值层析法计算得到的Q值剖面

利用层析反演法求得的Q值分布, 对地震数据的衰减做Q偏移补偿。图 12a图 12b分别以传统的和改进后的衰减旅行时层析求得的Q值作为输入进行Q补偿偏移得到的CRP道集。对比可以看出, 图 12b的道集质量得到进一步提高, 振幅能量得到更好恢复, 相位畸变得到更合理校正, 尤其3.5s位置改善十分明显; 同相轴更清晰, 连续性增强, 原本无法识别的薄层得以识别。

图 12 以改进前(a)、后(b)衰减旅行时层析法计算的Q值为输入进行Q补偿偏移后的道集对比

对CRP道集做叠加得到偏移叠加剖面如图 13所示。从频谱上看, 改进后的衰减旅行时层析方法的频带更宽, 主频由30Hz提高到40Hz, 分辨率得到提高。从剖面上看, 对于一些由于吸收衰减作用能量较弱、模糊不清的层位, 经吸收衰减补偿后得以正确显示, 地下介质结构刻画更加清晰。

图 13 改进前(a)、后(b)衰减旅行时层析法计算的Q值进行Q偏移补偿的叠加剖面
6 结论

本文采用改进的衰减旅行时Q层析方法计算地下纵、横向Q值分布。衰减旅行时层析法的关键在于准确求取衰减旅行时, 衰减旅行时上的较小误差都会导致最终Q值反演的较大偏差。计算模型衰减旅行时, 射线追踪的精度是关键, 在追踪过程中考虑速度梯度随传播方向的变化, 能够更准确地追踪射线路径, 获得更准确的传播时间。从实际资料中拾取衰减旅行时, 最优频带的选择通常会引入人为误差, 对数谱比反演方法不仅无需人为选择参与反演频率成分, 反演结果也更稳定。而且, 通过在反演过程中考虑反射波在不同界面内的旅行时之差, 可以有效解决上覆地层的影响。因此得到Q值反演结果更加符合地下介质的吸收衰减分布。

采用改进后的方法在数值模型和实际资料应用中取得了较好的Q值估计结果。经Q偏移补偿后, 提高了地震资料分辨率, 恢复了衰减造成的能量衰减和相位畸变, 为后续AVO分析和储层预测提供更可靠资料。

参考文献
[1] DASGUPTA R, CLARK R A. Estimation of Q from surface seismic reflection data[J]. Geophysics, 1998, 63(6): 2120-2128 DOI:10.1190/1.1444505
[2] 陈文爽, 管路平, 李振春, 等. 基于广义S变换的叠前Q值反演方法研究[J]. 石油物探, 2014, 53(6): 706-712
CHEN W S, GUAN L P, LI Z C, et al. Prestack Q-inversion based on generalized S transform[J]. Geophysical Prospecting for Petroleum, 2014, 53(6): 706-712
[3] 王小杰, 印兴耀. 基于零相位子波地层Q值估计[J]. 地球物理学进展, 2011, 26(6): 2090-2098
WANG X J, YIN X Y. Estimation of layer quality factors based on zero-phase wavelet[J]. Progress in Geophysics, 2011, 26(6): 2090-2098
[4] 王德利, 戴建芳. 基于射线路径的叠前高精度Q值估计方法[J]. 石油物探, 2013, 52(5): 475-481
WANG D L, DAI J F. Prestack high accuracy Q estimation based on ray path[J]. Geophysical Prospecting for Petroleum, 2013, 52(5): 475-481
[5] 张繁昌, 张汛汛, 张立强, 等. 基于自适应子波分解的品质因子Q提取方法[J]. 石油物探, 2016, 55(1): 41-48
ZHANG F C, ZHANG X X, ZHANG L Q, et al. Extraction method of Q estimation based on adaptive wavelet decompesation[J]. Geophysical Prospecting for Petroleum, 2016, 55(1): 41-48
[6] 李伟娜, 云美厚, 党鹏飞, 等. 基于微测井资料的双线性回归稳定Q估计[J]. 石油物探, 2017, 56(4): 483-490
LI W N, YUN M H, DANG P F, et al. Stability Q estimation by dual linear regression based on uphole survey data[J]. Geophysical Prospecting for Petroleum, 2017, 56(4): 483-490
[7] QUAN Y, HARRIST J M. Seismic attenuation tomography using the frequency shift method[J]. Geophysics, 1997, 62(3): 895-905 DOI:10.1190/1.1444197
[8] ZHANG C, ULRYCH T J. Estimation of quality factors from CMP records[J]. Geophysics, 2002, 67(10): 1542-1547
[9] 王宗俊. 基于谱模拟的质心法品质因子估算[J]. 石油物探, 2015, 54(3): 267-273
WANG Z J. Quality factor estimation by centroid frequency shift of spectrum fitting[J]. Geophysical Prospecting for Petroleum, 2015, 54(3): 267-273
[10] 张立彬, 王华忠, 马在田. 基于积分中值参变量法的质心频移Q值估算[J]. 石油物探, 2010, 49(3): 214-223
ZHANG L B, WANG H Z, MA Z T. Quality factor estimation by centroid frequency shift of integration median[J]. Geophysical Prospecting for Petroleum, 2010, 49(3): 214-223
[11] 高静怀, 杨森林. 利用零偏移VSP资料估计介质品质因子方法研究[J]. 地球物理学报, 2007, 50(4): 1198-1209
GAO J H, YANG S L. On the method of guaiity factors estimation from zero-offset VSP data[J]. Chinese Journal of Geophysics, 2007, 50(4): 1198-1209
[12] PARRA J O, XU P C, HACKERT C L. A borehole-model-derived algorithm for estimating QP logs from full-waveform sonic logs[J]. Geophysics, 2007, 72(4): E107-E117 DOI:10.1190/1.2734109
[13] LIU Y, WEI X. Interval Q inversion from CMP gathers:Part Ⅰ absorption equation[J]. Expanded Abstracts of 75th Annual Internat SEG Mtg, 2005: 1717-1720
[14] LIU Y, WEI X. Interval Q inversion from CMP gathers:Part Ⅱ absorption equation[J]. Expanded Abstracts of 75th Annual Internat SEG Mtg, 2005: 1721-1724
[15] BRZOSTOWSKI M A, MCMECHAN G A. 3-D tomographic imaging of near-surface seismic velocity and attenuation[J]. Geophysics, 1992, 57(3): 396-403 DOI:10.1190/1.1443254
[16] LIAO Q, MCMECHAN G A. Tomographic imaging of velocity and Q, with application to crosswell seismic data from the Gypsy Pilot Site, Oklahoma[J]. Geophysics, 1997, 62(6): 1804-1811 DOI:10.1190/1.1444281
[17] 严又生, 宜明理, 魏新, 等. 井间地震速度和Q值联合层析成像及应用[J]. 石油地球物理勘探, 2001, 36(1): 9-17
YAN Y S, YI M L, WEI X, et al. Joint tomographic imaging for cross-hole seismic velocity and Q value[J]. Oil Geophysical Prospecting, 2001, 36(1): 9-17
[18] NOWACK R L, MICHAEL P M. Inversion of seismic attributes for velocity and attenuation structure[J]. Geophysics of Journal International, 1997, 128(3): 689-700 DOI:10.1111/gji.1997.128.issue-3
[19] 樊计昌, 刘明军, 赵成彬, 等. 岫岩陨石坑三维Q值层析成像[J]. 地球物理学报, 2010, 53(10): 2367-2375
FAN J C, LIU M J, ZHAO C B, et al. Three-dimensional Q tomography for Xiuyuan meteorite impact center[J]. Chinese Journal of Geophysics, 2010, 53(10): 2367-2375 DOI:10.3969/j.issn.0001-5733.2010.10.010
[20] CAVALCA M, MOORE I, ZHANG L. Ray-based tomography for Q estimation and Q compensation in complex media[J]. Expanded Abstracts of 81st Annual Internat SEG Mtg, 2011: 3989-3992
[21] XIN K F, HE Y, XIE Y. Robust Q tomographic inversion through adaptive extraction of spectral features[J]. Expanded Abstracts of 84th Annual Internat SEG Mtg, 2014: 3726-3730
[22] JIN Z Q. An improved Q tomography inversion and its application[J]. Expanded Abstracts of 85th Annual Internat SEG Mtg, 2015: 1721-1724
[23] 魏文, 王小杰, 李红梅. 基于叠前道集小波域Q值求取方法研究[J]. 石油物探, 2011, 50(4): 356-360
WEI W, WANG X J, LI H M. Q estimation in wavelet domain using prestack data[J]. Geophysical Prospecting for Petroleum, 2011, 50(4): 356-360
[24] 王小杰, 印兴耀, 吴国忱. 基于变换的吸收衰减技术在含气储层预测中的应用研究[J]. 石油物探, 2012, 51(1): 38-42
WANG X J, YIN X Y, WU G S. Application on reservoir prediction using attenuation technique based on transform[J]. Geophysical Prospecting for Petroleum, 2012, 51(1): 38-42
[25] TONN R. The determination of seismic quality factor Q from VSP data:a comparison of different computational techniques[J]. Geophysical Prospecting, 1991, 45(1): 1-27