改进的衰减旅行时层析方法估计Q值

金子奇,孙赞东

(中国石油大学(北京),北京102249)

根据应用空间的不同,求取层间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]:

(1)

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

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

(2)

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

(3)

(4)

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

(5)

2　方向梯度射线追踪计算衰减旅行时

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

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

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

(7)

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

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

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

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

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

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

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

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

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

(12)

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

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

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

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

(13)

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

(14)

(15)

得到的衰减旅行时。求解(15)式可得模型更新量,代入公式(5)求得最终Q值分布。

4　模型数据测试

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

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

表1　水平层状介质模型参数

图5　水平层状模型合成炮集记录

表2　两层介质模型中改进前后射线追踪算法得到折射角对比

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

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

图7　衰减旅行时计算结果

图8　水平层状模型中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所示。

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

图10　共炮点道集记录(拾取反射轴较明显位置作为Q值层析反演输入)

图11　改进前(a)、后(b)衰减旅行时Q值层析法计算得到的Q值剖面

图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 ofQfrom surface seismic reflection data[J].Geophysics,1998,63(6):2120-2128

[2]陈文爽,管路平,李振春,等.基于广义S变换的叠前Q值反演方法研究[J].石油物探,2014,53(6):706-712

CHEN W S,GUAN L P,LI Z C,et al.PrestackQ-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 accuracyQestimation 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 ofQestimation 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.StabilityQestimation 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

[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 estimatingQPlogs from full-waveform sonic logs[J].Geophysics,2007,72(4):E107-E117

[13]LIU Y,WEI X.IntervalQinversion from CMP gathers:Part I absorption equation[J].Expanded Abstracts of 75thAnnual Internat SEG Mtg,2005:1717-1720

[14]LIU Y,WEI X.IntervalQinversion from CMP gathers:Part II absorption equation[J].Expanded Abstracts of 75thAnnual 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

[16]LIAO Q,MCMECHAN G A.Tomographic imaging of velocity andQ,with application to crosswell seismic data from the Gypsy Pilot Site,Oklahoma[J].Geophysics,1997,62(6):1804-1811

[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 andQvalue[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

[19]樊计昌,刘明军,赵成彬,等.岫岩陨石坑三维Q值层析成像[J].地球物理学报,2010,53(10):2367-2375

FAN J C,LIU M J,ZHAO C B,et al.Three-dimensionalQtomography for Xiuyuan meteorite impact center[J].Chinese Journal of Geophysics,2010,53(10):2367-2375

[20]CAVALCA M,MOORE I,ZHANG L.Ray-based tomography forQestimation andQcompensation in complex media[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:3989-3992

[21]XIN K F,HE Y,XIE Y.RobustQtomographic inversion through adaptive extraction of spectral features[J].Expanded Abstracts of 84thAnnual Internat SEG Mtg,2014:3726-3730

[22]JIN Z Q.An improvedQtomography inversion and its application[J].Expanded Abstracts of 85thAnnual Internat SEG Mtg,2015:1721-1724

[23]魏文,王小杰,李红梅.基于叠前道集小波域Q值求取方法研究[J].石油物探,2011,50(4):356-360

WEI W,WANG X J,LI H M.Qestimation 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 factorQfrom VSP data:a comparison of different computational techniques[J].Geophysical Prospecting,1991,45(1):1-27