近年来,在碳酸盐岩裂缝型储层预测、致密油气水平井布设及岩性油气藏勘探等领域,均对地震叠前深度偏移的精度提出了更高的要求[1-2]。由于地下介质各向异性的影响,采用各向同性叠前深度偏移方法得到的偏移结果通常会产生以下2个问题:①偏移剖面中的地震层位与测井分层在深度上相差较大,误差范围达不到地质精度要求;②各向同性叠前深度偏移的共成像点(CIP)道集存在远偏移距没有校平的现象,影响叠加效果,同时给AVO与AVA等叠前反演带来困难。因此,开展各向异性叠前深度偏移研究具有重要的现实意义。
目前,对各向异性叠前深度偏移算法的研究比较活跃,尤其以VTI、TTI和ORT等介质的逆时偏移为代表,将各向异性叠前深度偏移算法推到了一个新的高度[3-4]。研究表明,现阶段影响叠前深度偏移成像效果的瓶颈并不是偏移算法本身,而是速度模型的准确性问题[5]。各向异性介质的速度建模不仅要求取各向异性速度场,还要求取各向异性参数场,常用的各向异性参数包括Thomsen参数[6]和η参数[7]。Alkhalifah[8]和Grechka等[9]先后提出利用纵波非双曲时差分析和旅行时反演获取短排列动校正速度Vnmo和参数η的有效值,并通过剥层法将其换算成层间值,但这2种方法都会产生不稳定现象。随后Dewangan等[10]和Wang等[11]提出了与速度无关的剥层法旅行时反演以解决上述问题,并将该方法由VTI介质扩展到正交各向异性介质。在上述几种方法中,仅利用地面地震资料便可得到Vnmo和η,且这2个参数可满足所有关于P波的时间域处理。
对于P波的深度域处理,通常要求取纵波垂向速度VP0、各向异性参数δ和ε。目前求取这3个参数的方法有实验室岩石物理测量[12]、反射横波(PSV波或SV-SV波)弹性参数反演[13]或基于各向异性假设的全波形反演[14]等,但这些方法或者实现难度大、成本高,或者流程繁琐,在实际生产中难以大范围推广。目前,实际生产中常用的方法是将地震资料和测井数据或VSP数据进行结合,以同时求取这3个参数。Zhang等[15]利用VSP数据进行初至旅行时反演提取垂向速度VP0和各向异性参数δ与ε,在实际应用中均取得了较好的效果。常用的井震联合法[16]先分别计算地震层厚和测井层厚,再在井点位置计算出每一层的δ值,然后进行空间插值生成δ平面图,最后生成δ数据体。该方法计算简单,易于操作,但无法进行迭代,计算结果的精度较低。为提高各向异性速度场和参数场的精度,将网格层析技术应用于VTI介质井震联合各向异性速度建模中,建立并求解由深度层位模型、井震闭合差以及初始各向异性体构成的线性方程组,沿射线路径同时更新每个网格点的各向异性速度和参数,并通过迭代计算进一步提高速度模型的精度。对比常规井震联合法和井震联合网格层析法速度建模的原理及实现过程,并以四川盆地S工区为例,分别应用2种方法进行各向异性速度建模,以期提高各向异性速度场和参数场的精度、减小井震闭合差、改善局部成像质量。
1 常规井震联合法速度建模的实现过程目前,常规井震联合法仍是一种主流的各向异性速度建模方法,它的实现过程有如下几步:
第1步:采用层控网格层析技术建立三维各向同性速度体,使深度偏移的共成像点道集的中、近偏移距被拉平,整体构造形态合理,并得到全区的叠前深度偏移剖面。
第2步:由叠前深度偏移剖面标定的反射层深度可以计算出各反射层的各向同性地层厚度,而根据测井分层数据可以计算出对应地层的各向异性地层厚度,从而计算出井点各个反射层的各向异性参数δ
$ \delta = \frac{1}{2}\left[{{{\left( {\frac{{\Delta {Z_{\text{i}}}}}{{\Delta {Z_{\text{a}}}}}} \right)}^2}-1} \right] $ | (1) |
式中:ΔZi,ΔZa分别为某一地层的各向同性层厚度和各向异性层厚度,m。
结合地质模型,通过空间插值便可以得到工区的各向异性参数δ数据体。
第3步:结合δ数据体,建立垂向速度场,即各向异性速度场
$ {V_{{\text{P0}}}} = \frac{{{V_{{\text{nmo}}}}}}{{\sqrt {1 + 2\delta } }} \approx \frac{{{V_{\text{i}}}}}{{\sqrt {1 + 2\delta } }} $ | (2) |
式中:Vi表示各向同性层速度,m/s;Vnmo为短排列动校正速度,m/s。
地下介质的δ通常为正值,因此,各向同性层速度Vi通常大于垂向速度VP0,各向同性深度偏移的成像深度也大于井上测量深度。
第4步:假设ε = δ,进行叠前深度偏移,在偏移后的共成像点道集上进行剩余ε分析,并通过迭代优化得到最终的ε体。
从上述实现过程可以看出,常规井震联合法各向异性速度建模只是计算出井点位置处每一层的δ值,再通过空间插值得到三维δ体,层内δ垂向上保持不变,是一种简单的垂向拉伸方法,无法进行迭代。同时,从式(2)也可以看出,短排列动校正速度Vnmo只是近似等于各向同性层速度Vi。因此,该方法求取的各向异性速度和参数的精度较低。
2 井震联合网格层析反演原理网格层析是一种反射波层析反演方法,近年来被广泛应用于三维深度域速度建模。该方法利用射线追踪原理,建立并求解由共成像点道集拾取的剩余时差(RMO)及模型参数变化量等组成的线性方程组,通过迭代计算,逐步更新优化速度模型,最终消除共成像点道集的RMO,拉平CIP道集。层析反演的线性方程组如下
$ A\Delta m = \Delta z $ | (3) |
式中:Δm为模型参数的变化量矩阵,模型参数可以为速度,也可以为各向异性参数δ和ε;Δz为数据残差向量,可以是RMO,也可以是每个层位的井震闭合差;A为表示数据残差随模型参数扰动而变化的灵敏度矩阵。
对于各向异性介质,求取δ时,将井震闭合差作为式(3)的右端项;求取ε时,将CIP道集的剩余时差作为式(3)的右端项。
由于层析反演线性方程组具有严重的病态性,为了提高计算的稳定性、减少反演的多解性,通常采用加权的最小二乘算法对其进行求解,此时,层析反演方程组变为
$ \left( {{A^T}C_{\text{d}}^{-1}A + C_{\text{m}}^{-1} + {L^T}L} \right)\Delta m = {A^T}C_{\text{d}}^{-1}\Delta z $ | (4) |
式中:Cd是一个对角协方差矩阵,其对角元素为数据的方差;Cm是模型的协方差矩阵,起到阻尼因子的作用,可以压制过大的模型扰动;L为平滑约束算子,可以沿着构造方向对模型扰动进行平滑,使层析方程组收敛到一个符合地质构造的最优解。
3 井震联合网格层析各向异性速度建模实现过程三维VTI介质各向异性速度建模须要求取各向异性速度、各向异性参数δ和ε,本次研究采用井震联合网格层析方法进行各向异性速度建模,具体实现流程如图 1所示。
第1步:采用层控网格层析技术建立三维各向同性速度体,进行三维各向同性叠前深度偏移。
第2步:在深度偏移剖面上解释出主要的地质层位,与测井数据解释的层位进行匹配对比,计算井点位置处主要地质层位的井震闭合差,并对其进行插值,得到各个层位的井震闭合差平面图。
第3步:通过网格层析迭代,计算各向异性速度和各向异性参数δ。第一次进行网格层析迭代时,输入各向同性速度体、深度层位模型和井震闭合差模型,输出更新后的各向异性速度体、深度层位模型和各向异性参数δ。然后,利用新的深度层位模型,重复第2步,计算各个层位新的井震闭合差,再进行下一轮网格层析迭代,此时,输入上一轮迭代得到的各向异性速度体、深度层位模型和井震闭合差模型,输出更新后的各向异性速度体、深度层位模型和参数δ。一般进行2~3轮网格层析迭代即可使井震闭合差迅速减小,满足精度要求。通过网格层析更新δ时,采用椭圆各向异性假设,保持ε = δ。为保证层析反演的精度和计算效率,须要选取合适的网格大小。网格太大,反演模型的精度和分辨率会降低;网格太小,计算效率又会明显下降。实际应用中,通常将x和y方向上的层析网格取为面元大小的10~20倍,将z方向上的层析网格取为200~400 m。
第4步:网格层析更新各向异性参数ε。利用第3步计算得到的各向异性速度体、各向异性参数δ和ε,进行三维各向异性叠前深度偏移,得到共成像点道集。此时,CIP道集的中、近偏移距均被基本拉平,若实际介质模型不符合椭圆各向异性假设,远偏移距会存在没有拉平的现象。自动拾取CIP道集的剩余时差,进行网格层析。此时,输入第3步得到的各向异性速度、各向异性参数δ和ε,只更新ε,各向异性速度和δ保持不变。网格层析更新各向异性参数ε时,仍然可以进行多轮迭代,直到CIP道集的远偏移距被拉平,剩余时差被消除为止。
第5步:利用前面4步得到的最终各向异性速度、各向异性参数δ和ε,进行各向异性叠前深度偏移,得到最终的叠前深度偏移结果。
4 应用实例四川盆地S区块属山地-丘陵地貌,区内西北部地形起伏较大,目的层褶皱强烈,断裂、断块比较发育,因此,该地区地震资料处理对速度建模和偏移成像的精度要求较高。首先对地震资料进行三维各向同性叠前深度偏移。图 2为各向同性叠前深度偏移联井剖面。从图 2可以看出,工区内3口井的5个主要地质层位与偏移剖面上解释的地震层位在深度上相差较大。这5个地质层位中,J3s为沙溪庙组底界,J1z为自流井组底界,T3x1为须家河组底界,T2l1为雷口坡组底界,P2l为茅口组顶界。
表 1为3口井主要地质层位的深度、井震误差及分别应用前文提到的2种建模方法得到的井震闭合差(B井和C井缺失部分地质层位的钻井深度数据)。从表 1可以看出,各向同性叠前深度偏移的井震误差较大,地震层位与测井分层之间的误差为100~300 m,其中,J3s层的误差达到10%,目的层P2l的误差达到5%,这与地质上小于1%的精度要求存在较大差距。为减小井震闭合差,同时考虑到工区大部分地层的倾角较小,本次研究采用VTI介质井震联合网格层析方法进行各向异性速度建模,并与常规井震联合法进行对比。
采用网格层析方法计算各向异性速度和各向异性参数δ时,共进行2次迭代。图 3为过A井的各向同性速度剖面及2次迭代后的各向异性速度剖面的对比图。图 4为网格层析第2次迭代后的δ剖面和常规井震联合法得到的δ剖面对比图。从图 4看出,采用常规井震联合法得到的δ剖面的结构比较简单,层内的δ在垂向上保持不变,横向上的δ由井点处的δ内插得到,而采用网格层析得到的δ剖面的信息更为丰富,垂向和横向都有变化。
图 5为T3x1层的各向同性井震闭合差以及网格层析第1次迭代、第2次迭代后的井震闭合差。从图 5可以看出,各向同性的井震闭合差为-178~ -170 m;第1次网格层析迭代后,井震闭合差减小到-47~-37 m;第2次网格层析迭代后,减小到-5~-3 m。从各向异性叠前深度偏移联井剖面上也可以看到,第2次网格层析迭代后,T3x1层的地震层位与测井层位基本吻合(图 6)。
对比表 1中井震联合网格层析法与常规井震联合法得到的井震闭合差,可以看到,前者的井震闭合差均在10 m以内,误差小于0.15%,后者的误差最大可达3%,显然,井震联合网格层析方法的精度高于常规井震联合法。
图 7为网格层析更新各向异性参数ε前后的CIP道集。从图 7可以看出,更新各向异性参数ε后,CIP道集远偏移距的拉平程度更高,因此,在进行AVO与AVA等叠前反演时可以利用更多的远偏移距信息。同时,从更新ε前后的叠前深度偏移剖面上也可以看到,工区局部的成像效果得到了明显改善(图 8)。
(1) 将网格层析技术应用于三维VTI介质井震联合各向异性速度建模,可以沿射线路径同时更新模型中每个网格点的各向异性速度和参数,并通过迭代计算进一步提高速度模型的精度。
(2) 实际资料处理表明,相对于常规井震联合法,井震联合网格层析法建立的各向异性速度模型的精度更高,地震层位与测井分层的吻合度更好。
(3) 更新ε参数后,各向异性叠前深度偏移的共成像点道集的远偏移距拉平程度更高,为叠前反演提供了更丰富的远偏移距信息,并有效改善了局部成像效果。
[1] |
李凡异, 狄帮让, 魏建新, 等. 溶洞体宽度对偏移剖面反射振幅影响的定量研究. 岩性油气藏, 2016, 28(5): 113-116. LI F Y, DI B R, WEI J X, et al. Quantitative study of the effect of cavern width on reflection amplitude of migrationsection. Lithologic Reservoirs, 2016, 28(5): 113-116. |
[2] |
魏新善, 胡爱平, 赵会涛, 等. 致密砂岩气地质认识新进展. 岩性油气藏, 2017, 29(1): 11-20. WEI X S, HU A P, ZHAO H T, et al. New geological understanding of tight sandstone gas. Lithologic Reservoirs, 2017, 29(1): 11-20. |
[3] |
ZHANG Y, ZHANG H Z, ZHANG G Q. A stable TTI reverse time migration and its implementation. Geophysics, 2011, 76(3): WA3-WA11. DOI:10.1190/1.3554411 |
[4] |
THOMAS M, MOTHI S, MCGILL P. Improving subsalt images using tilted-Orthorhombic RTM in Green Canyon, Gulf of Mexico. 82th Annual Meeting, SEG Expanded Abstracts, 2012, 693-697. |
[5] |
陈可洋. 逆时成像技术在大庆探区复杂构造成像中的应用. 岩性油气藏, 2017, 29(6): 91-100. CHEN K Y. Application of reverse-time migration technology to complex structural imaging in Daqing exploration area. Lithologic Reservoirs, 2017, 29(6): 91-100. |
[6] |
THOMSEN L. Weak elastic anisotropy. Geophysics, 1986, 51(10): 1954-1966. DOI:10.1190/1.1442051 |
[7] |
ALKHALIFAH T, TSVANKIN I. Velocity analysis for transversely isotropic media. Geophysics, 1995, 60(5): 1550-1566. DOI:10.1190/1.1443888 |
[8] |
ALKHALIFAH T. Velocity analysis using nonhyperbolic moveout in transversely isotropic media. Geophysics, 1997, 62(6): 1839-1854. DOI:10.1190/1.1444285 |
[9] |
GRECHKA V, TSVANKIN I. Feasibility of nonhyperbolic moveout inversion in transversely isotropic media. Geophysics, 1998, 63(3): 957-969. DOI:10.1190/1.1444407 |
[10] |
DEWANGAN P, TSVANKIN I. Velocity-independent layer stripping of PP and PS reflection traveltimes. Geophysics, 2006, 71(4): U59-U65. DOI:10.1190/1.2210975 |
[11] |
WANG X X, TSVANKIN I. Estimation of interval anisotropy parameters using velocity-independent layer stripping. Geophysics, 2009, 74(5): WB117-WB127. DOI:10.1190/1.3157462 |
[12] |
BACHRACH R, OSYPOV K, NICHOLS D, et al. Applications of deterministic and stochastic rock physics modelling to anisotropic velocity model building. Geophysical Prospecting, 2013, 61(2): 404-415. DOI:10.1111/gpr.2013.61.issue-2 |
[13] |
杜丽英, 杜丽娟, 肖乾华, 等. VTI介质中弹性参数地震反演方法. 石油勘探与开发, 2002, 29(3): 59-63. DU L Y, DU L J, XIAO Q H, et al. A new approach to seismic inversion of elastic parameters for VTI media. Petroleum Exploration and Development, 2002, 29(3): 59-63. |
[14] |
PRIEUX V, BROSSIER R, GHOLAMI Y, et al. On the footprint of anisotropy on isotropic full waveform inversion:the Valhallcase study. Geophysical Journal International, 2011, 187(3): 1495-1515. DOI:10.1111/gji.2011.187.issue-3 |
[15] |
ZHANG Z, LIN G, CHEN J, et al. Inversion for elliptically anisotropicvelocity using VSP reflection traveltimes. Geophysical Prospecting, 2003, 51(2): 159-166. DOI:10.1046/j.1365-2478.2003.00362.x |
[16] |
刘茂诚. 一个各向异性速度分析应用实例. 石油地球物理勘探, 2010, 45(4): 525-529. LIU M C. An application case study for anisotropic velocity analysis. Oil Geophysical Prospecting, 2010, 45(4): 525-529. |