2. 中石化胜利油田分公司勘探开发研究院, 山东东营 257015;
3. 长江大学地球科学学院, 武汉 430100
2. Exploration and Development Research Institute of Shengli Oilfield, SINOPEC, Shandong Dongying 257015, China;
3. School of Geosciences, Yangtze University, Wuhan 430100, China
速度是勘探地震学中最核心的问题,介质速度包括纵波速度和横波速度,横波速度是进行叠前AVO反演(Buland and Omre,2003; 杨培杰和印兴耀,2008)、弹性阻抗反演(Connolly,1999; 王保丽等,2005)以及流体识别(Russell et al.,2003; 印兴耀等,2013;杨培杰等,2015)非常重要的信息.然而,实际测井数据中大多缺少横波速度信息,特别是老井更是没有横波速度曲线(李维新等,2009).以胜利油田济阳坳陷为例,截止2013年,济阳坳陷有偶极横波测井的探井仅70多口,远远无法满足精细勘探开发的需要,因此有必要开展横波速度预测方面的研究.
目前,在没有横波信息的情况下,估计横波速度的方法主要有两种,一是统计拟合法(Castagna et al.,1985; Han et al.,1986),通过对纵波速度和横波速度进行统计分析,拟合得到横波速度-纵波速度关系式;二是理论公式法(Kuster and Toksoz,1974; Xu and White,1995),给出横波速度和其他参数之间的一种理论关系,进行横波速度的估计.这些方法的局限性在于它们往往是大量数据统计的结果,只反映了一般性的规律,应用于具体区域研究时往往存在较大的误差.
因此,针对统计拟合法和理论公式法横波估计存在的局限性,很多学者基于岩石物理模型,进行了横波速度预测方法的综合研究(Greenberg and Castagna,1992; Xu and White,1996; Jørstad et al.,1999).Greenberg等假设纵横波速度间有稳健的关系,基于Biot-Gassmann理论预测横波速度;Xu等使用Kuster-Tøksoz理论和微分等效介质理论相结合预测横波速度,并运用孔隙纵横比的概念来表征干岩石颗粒的接触关系;Jørstad等运用有效介质理论预测横波速度,并且认为基于岩石物理方法的横波预测的准确度要高于统计拟合法得到的横波速度.国内也有很多学者(孙福利等,2008; 白俊雨等,2012; 李晓明,2012; 吴志华,2012; 熊晓军等,2012;张广智等,2012)开展了基于岩石物理模型的横波预测研究,并取得了较好的研究成果.
然而,大多数基于岩石物理的横波速度预测方法需要对孔隙形态进行假设.Brown和Korringa(1975)等的实验室数据分析表明,与砂岩有关的孔隙纵横比并不是定值,通过实际电镜观察也会发现,很难用确定的纵横比来描述孔隙的变化,这在一定程度上增加了横波预测过程复杂性和不确定性,也在一定程度上降低了结果的可靠性.同时,对于横波预测目标函数的求解问题,大多数学者选择遗传算法或模拟退火等非线性优化算法(吴志华,2012),其不足之处在于计算效率低,并且横波预测结果不太稳定;或是采用传统的迭代优化方法(Lee,2006; 张广智等,2012),其不足之处在于计算效率较低,算法易陷入局部最优解.以Lee(2006)提出的P-L模型横波估计方法为例,该方法通过P-L模型来描述岩石基质弹性模量、干岩石弹性模量以及孔隙度之间的关系,增加了问题的复杂性,同时,该方法采用牛顿-拉普森迭代方法直接求解目标函数,因此在求解效率和结果的稳定性方面并不高.
针对上述问题,本文在前人研究的基础上(Pride,2004; Lee,2006; 白俊雨等,2012),开展了精细横波预测方法的研究,开发了基于变形P-L模型矩阵方程迭代的横波预测新方法.首先,对Pride模型和Lee模型(P-L模型)进行改进,提出拟固结指数的概念,拟固结指数可以认为是岩石的胶结程度和孔隙形状的一种综合的响应,拟固结指数将干岩石模量和岩石基质模量相联系.基于拟固结指数的变形P-L模型的优势在于,在没有降低P-L模型准确度的情况下简化了岩石物理模型的复杂度;其次,采用V-R-H平均估算岩石基质模量,利用Wood公式计算混合流体体积模量;然后,基于Gassmann方程建立起了饱和流体岩石弹性模量与干岩石模量、岩石基质模量、混合流体模量之间的关系,并计算理论上的纵波速度;最后,借鉴地震反演的思路,通过比较实测纵波速度与计算的理论纵波速度大小,建立起了横波预测的目标函数,将该目标函数的最优化问题转化为线性矩阵方程组迭代求解问题,经过几次迭代就可以求解出合适的拟固结指数,进而得到横波速度预测结果.
2 变形P-L模型2.1 公式推导Pride等(2004)引入固结指数的概念,由岩石基质弹性模量计算干岩石弹性模量,如下式所示:
Lee(2006)对式(1)中计算剪切模量公式做了如下修正:
Lee通过实际应用证明,在应用γ代替常数1.5后,会有更好的横波预测效果.将式(1)和式(2)总称为P-L模型.Lee利用反演方法由实测纵波速度与计算的纵波速度误差得到固结指数,进而计算横波速度.
可以看出,在P-L模型中,固结指数α处于分母的位置,因此,这在一定程度上增加了求解α的复杂度.
因此,笔者对P-L模型进行如下的修改:
将式(3)称为变形P-L模型,其中的φ和ξ被称为拟固结指数,0<φ,ξ≤1.
变形P-L模型用拟固结指数代替固结指数,是对P-L模型的一种改进,可以看出,变形P-L模型与P-L模型在本质上是一样的,但是却将基质模量和干岩石骨架模量由复杂的关系变得简单,在没有降低公式准确度的同时简化了问题的复杂度,易于下面问题的求解.
2.2 拟固结指数与Biot系数目前常用的岩石物理模型中,如果知道了岩石的矿物成分和孔隙流体的组成,就可以计算出岩石基质的弹性模量(体积模量Km和剪切模量μm)和混合流体的体积模量(Kf),但是,干岩石的弹性模量(体积模量Kd和剪切模量μd)却难以获得,很多学者用Biot系数来表征干岩石的弹性模量和岩石基质的弹性模量之间的关系:
对比式(4)的拟固结指数,可以得到下面的关系式:
可以看出,拟固结指数和Biot系数之间是一种线性的关系,并且具有相同的物理意义,都可以看作是岩石胶结程度和孔隙形状的一种综合的响应.
本算法就是通过反演得到拟固结指数φ和ξ,然后带入横波计算公式,从而得到了最终的横波预测结果.
3 基于变形P-L模型的纵波速度计算公式横波预测的目标函数是通过比较实测纵波速度与计算的纵波速度的大小而建立的,即,当计算的纵波速度和实测纵波速度的误差达到最小时,饱和流体岩石纵波速度公式中的体积模量和剪切模量就是正确的,此时通过拾取纵波速度中的饱和岩石的剪切模量信息,结合密度测井数据,就可以得到横波速度的预测结果.
根据纵波速度的计算公式,计算饱和岩石的纵波速度需要用到饱和岩石体积模量、饱和岩石剪切模量和密度信息(Russell et al.,2003).
Ks由低频Biot-Gassmann理论获得:
下面介绍岩石骨架的弹性模量、混合流体体积模量所采用的岩石物理模型.
3.1 岩石骨架的体积模量Km和剪切模量μm计算Km,首先需要知道岩石矿物的组成部分.这个可以通过实验室岩芯测试获得,也可以通过统计数据获得(Mavko et al.,1998),如砂泥岩的基质可以认为是石英和黏土组成的,黏土的百分含量可以由Gr测井曲线获得.对于由N种成分组成的岩石:
一旦岩石矿物的组成部分确定了,Km和μm就可以由Voigt-Reuss-Hill(VRH)平均来获得,Hill(1952)提出对上下界限求取算术平均的方法来近似岩石有效弹性模量值,如下所示:
VRH平均是求取岩石有效弹性模量的一种非常简单的方法,在等效模量计算方面得到了广泛的使用与推广.
3.2 混合流体体积模量Kf利用Wood(1955)公式计算得到混合流体体积模量Kf:
储层中流体的体积模量会随地层压力的增加而增大,随温度的升高而降低,对于油和水,这种影响可以忽略不计,但是对于气体来说,压力和温度会对其体积模量产生很大的影响,因此,不能忽略,Batzle和Wang(1992)给出了气体体积模量计算公式:
由于流体不传播横波,所以饱和岩石的剪切模量与岩石骨架的剪切模量相同,即:
该式即为基于变形P-L模型的纵波速度计算公式.
4 目标函数的建立与求解式(15)基于岩石物理模型,得到了纵波速度的计算公式,下面,借鉴地震反演的思路(杨培杰,2008),通过比较实测纵波速度与计算纵波速度大小,建立非线性的横波预测目标函数,并将目标函数的非线性最优化问题转化为迭代求解一个线性矩阵方程组的问题.
首先给出(16)式:
公式的含义为,实测纵波速度与计算纵波速度是不完全一样的,它们之间有个误差,定义为ε.加平方的目的是为了便于后面公式的推导.
对公式(15)进行变形,并与式(16)相结合:
进一步,将式(17)写成下面的形式:
将式(19)扩展为矩阵的形式:
则,定义最终的目标函数为
对式(22)进行如下的处理:
将该矩阵对拟固结指数求导(杨培杰,2008),最终得到了问题求解的矩阵方程:
通过求解该方程,在得到向量m后,拾取向量中的剪切模量拟固结指数ξ(详见式(20)),通过公式(25)最终可实现横波曲线的预测.
该方法通过若干测井曲线,以及压力和温度等参数实现了精细的横波速度预测,具体步骤如下:
步骤1:输入5条测井曲线以及岩石基质参数
测井曲线包括纵波速度、密度、泥质含量、孔隙度、含水饱和度,岩石基质参数包括基质的模量,并输入压力和温度参数.
步骤2:计算岩石基质模量、混合流体体积模量、岩石骨架模量
采用V-R-H平均估算岩石基质模量,利用Wood公式计算混合流体体积模量.
步骤3:计算饱和岩石体积模量
基于Gassmann方程建立起了饱和流体岩石弹性模量与干岩石模量、岩石基质模量、混合流体模量之间的关系,并计算理论上的纵波速度.
步骤4:如下迭代求解矩阵方程组
(1)给定岩石骨架模量的初始值;
(2)将岩石骨架模量和步骤2中计算的岩石基质模量、混合流体体积模量一起代入Gassmann方程;
(3)组合得到矩阵方程公式(24),求解得到拟固结指数;
(4)将求解得到拟固结指数代入变形P-L模型公式(3),得到迭代后的岩石骨架模量;
(5)循环(2)—(4),直到达到设定的迭代次数,输出最终的拟固结指数.
步骤5:输出横波预测结果
拾取拟固结指数中的剪切模量拟固结指数,代入公式(25),最终得到得到横波速度预测结果.
基于变形P-L模型的矩阵方程迭代横波预测流程如图 1所示.
方法验证的井来自于国内某油田的C井,深度段为1400~1500 m,此处的压力大约在15 MPa,温度约65 ℃,所在地层为馆陶组,属于河流相沉积.
将纵波速度、密度、泥质含量、孔隙度、含水饱和度等5条曲线以及压力和温度参数作为输入,如图 2所示.从泥岩含量曲线可以看出,该层段共有四套砂岩,其余都是泥岩,孔隙度曲线与泥岩含量成反比,含水饱和度曲线在深度1470 m左右含水饱和度不为1.0,为一套厚度约为3 m的含油储层,其他的三套砂岩为含水砂岩.
首先分析一下横波预测迭代过程及其收敛性,如图 3所示,图(A—D)分别为1~4次迭代结果,每次迭代的图中,(a1,b1,c1,d1)黑色实线为实际测井的纵波速度,红色虚线为计算的纵波速度,(a2,b2,c2,d2)黑色实线为预测出的横波速度.给定的岩石骨架模量初始值为0,从图中看出,该方法仅仅通过4次迭代,计算的理论纵波速度实际测井的纵波速度已经非常接近,说明该法有很好的收敛性,同时,预测出的横波速度也非常的稳定.
将本文提出的方法与某商业软件的横波预测结果进行对比分析,说明和验证本文所提出的新方法的收敛性和有效性.该软件是一款用于储层预测的综合软件,主要包括地震反演、属性分析、时频分析以及横波预测等模块,在实际勘探生产中得到了较为广泛的应用.如图 4所示,红色曲线为应用本文方法预测的横波速度以及纵横波速度比(VP/VS),蓝色曲线为应用该软件预测的横波速度以及VP/VS.
仅从横波速度预测结果来看,两种方法预测的横波速度基本重合,无法分辨哪种方法效果更好.下面重点分析图中的黑框所示的两套储层(即水层和油层)的VP/VS,进一步的分析对比.首先,对于新方法预测的VP/VS,水层和油层之间存在着一个差距,即红色箭头所指处的两条绿色虚线间的距离Do.而对于某软件预测的VP/VS,这两套储层之间并没有差距,也就是说,水层和油层具有了一样的VP/VS.
根据流体替代理论,砂岩在含气、油、水时的VP/VS是不一样的,一般来说,含油砂岩的VP/VS要低于含水砂岩的VP/VS,而含气砂岩的VP/VS要低于含油砂岩的VP/VS.从上面的横波预测结果分析可以看出,对于含油储层和含水储层,本文提出的方法所预测的VP/VS是有差别的(Do),即,含油储层的VP/VS比含水储层的VP/VS低了Do,这和理论的结果是一致的,而通过某软件预测的VP/VS在含油储层和含水储层却是一样的,这是不合理的,因此,新方法横波曲线预测准确性要高于某软件的横波曲线预测结果.
进一步,图 5为假设储层含气时,用新方法预测的横波曲线和VP/VS,图中的Dg表示含气储层和含水储层VP/VS的差别,对比Do可以看出,Dg>Do,即,新方法预测的横波曲线和VP/VS能够很好的区分含气、含油和含水储层,这也进一步说明了新方法横波曲线预测的有效性和准确性.
实际应用来自于国内某油田的Q井,录井图如图 6所示.该井在馆陶组钻遇2 m气层,深度1023~1025 m,图 6中黄色矩形1所示;油层7 m,深度1027~1033 m,图中红色矩形2所示;水层3 m,深度1045~1048 m,图中蓝色矩形3所示.输入5条测井数据以及岩石各基质参数,并输入压力和温度参数,设置迭代求解矩阵方程的次数为4次.
横波预测结果如图 7所示,从图中可以看出,横波曲线预测结果稳定,油、气、水层的VP/VS在数值的大小关系上也和理论值是一致的.进一步用预测的横波速度进行了该地区叠前流体因子直接提取方面的研究,并通过反演结果进行了该地区的流体识别,取得了较好的应用效果.
本文提出的横波预测方法具有以下两个创新点:
(1)对P-L模型进行改进,提出了拟固结指数的概念,变形P-L模型在没有降低P-L模型准确性的同时简化了岩石物理模型的复杂度,使得下一步的横波估计过程更加简洁,提高了横波预测的实用性;
(2)借鉴地震反演的思路,通过对目标函数求导,将非线性的目标函数最优化问题转化为迭代求解线性矩阵方程组问题,通过几次迭代求解矩阵方程组,就可以得到横波预测结果,提高了横波预测过程的效率和结果的稳定性.
[1] | Bai J Y, Song Z X, Su L, et al. 2012. Error analysis of shear-velocity prediction by the Xu-White model. Chinese J. Geophys. (in Chinese), 55(2): 589-595, doi: 10.6038/j. issn.0001- 5733. 2012.02.021. |
[2] | Batzle M L, Wang Z J. 1992. Seismic properties of pore fluids. Geophysics, 57(11): 1396-1408. |
[3] | Brown R J S, Korringa J. 1975. On the dependence of the elastic properties of a porous rock on the compressibility of the pore fluid. Geophysics, 40(4): 608-616. |
[4] | Buland A, Omre H. 2003. Bayesian linearized AVO inversion. Geophysics, 68(1): 185-198. |
[5] | Castagna J P, Batzle M L, Eastwood R L. 1985. Relationships between compressional-wave and shear-wave velocities in clastic silicate rocks. Geophysics, 50(4): 571-581. |
[6] | Connolly P. 1999. Elastic impedance. The Leading Edge, 18(4): 438-452. |
[7] | Gassmann F. 1951. Elastic waves through a packing of spheres. Geophysics, 16(4): 673-685. |
[8] | Greenberg M L, Castagna J P. 1992. Shear-wave velocity estimation in porous rocks: Theoretical formulation, preliminary verification and applications. Geophysical Prospecting, 40(2): 195-209. |
[9] | Han D H, Nur A, Morgan D. 1986. Effects of porosity and clay content on wave velocities in sandstones. Geophysics, 51(11): 2093-2107. |
[10] | Hill R. 1952. The elastic behaviour of a crystalline aggregate. Proceedings of the Physical Society. Series A, 65(5): 349-354. |
[11] | Jørstad A, Mukerji T, Mavko G. 1999. Model-based shear-wave velocity estimation versus empirical regressions. Geophysical Prospecting, 47(5): 785-797. |
[12] | Kuster G T, Toksöz M N. 1974. Velocity and attenuation of seismic waves in two-phase media, Part 1. Theoretical formulations. Geophysics, 39(5): 587-606. |
[13] | Lee M W. 2006. A simple method of predicting S-wave velocity. Geophysics, 71(6): F161-F164. |
[14] | Li W X, Wang H, Yao Z X, et al. 2009. Shear-wave velocity estimation and fluid substitution by constraint method. Chinese J. Geophys. (in Chinese), 52(3): 785-791. |
[15] | Li X M, Chen S Q, Li X Y. 2012. An inversion method for near-surface shear-wave velocity using multi-component seismic data. OGP (in Chinese), 47(4): 532-536. |
[16] | Mavko G, Mukerji T, Dvorkin J. 1998. The Rock Physics Handbook: Tools for Seismic Analysis in Porous Media. Cambridge: Cambridge University Press. |
[17] | Pride S R, Berryman J G, Harris J M. 2004. Seismic attenuation due to wave-induced flow. Journal of Geophysical Research, 109: B01201. |
[18] | Russell B H, Hedlin K, Hilterman F J, et al. 2003. Fluid-property discrimination with AVO: A Biot-Gassmann perspective. Geophysics, 68(1): 29-39. |
[19] | Sun F L, Yang C C, Ma S H, et al. 2008. An S-wave velocity predicted method. Progress in Geophysics (in Chinese), 23(2): 470-474. |
[20] | Wang B L, Yin X Y, Zhang F C. 2005. Elastic impedance inversion and its application. Process in Geophysics (in Chinese), 20(1): 89-92. |
[21] | Wood A W. 1955. A Textbook of Sound. New York: The MacMillan Co. |
[22] | Wu Z H. 2012. The study on rock physics analysis for porous media with fluid[Master thesis](in Chinese). Qingdao: China University of Petroleum. |
[23] | Xiong X J, Lin K, He Z H. 2012. A method for S-wave velocity estimation based on equivalent elastic modulus inversion. OGP (in Chinese), 47(5): 723-727. |
[24] | Xu S Y, White R E. 1995. A new velocity model for clay-sand mixtu res. Geophysical Prospecting, 43(1): 91-118. |
[25] | Xu S Y, White R E. 1996. A physical model for shear-wave velocity prediction. Geophysical Prospecting, 44(4): 687-717. |
[26] | Yang P J. 2008. Seismic wavelet blind extraction and non-linear inversion[Ph. D. thesis](in Chinese). Dongying: China University of Petroleum. |
[27] | Yang P J, Yin X Y. 2008. Non-linear quadratic programming Bayesian Prestack inversion. Chinese J. Geophys. (in Chinese), 51(6): 1876-1882. |
[28] | Yang P J, Wang C J, Bi J F, et al. 2015. Direct extraction of the fluid factor based on variable point-constraint. Chinese J. Geophys. (in Chinese), 58(6): 2188-2200, doi: 10.6038/cjg20150631. |
[29] | Yin X Y, Zhang S X, Zhang F. 2013. Two-term elastic impedance inversion and Russell fluid factor direct estimation method for deep reservoir fluid identification. Chinese J. Geophys. (in Chinese), 56(7): 2378-2390, doi: 10.6038/cjg20130724. |
[30] | Zhang G Z, Li C C, Yin X Y, et al. 2012. A shear velocity estimation method for carbonate rocks based on the improved Xu-White model. OGP (in Chinese), 47(5): 717-722. |
[31] | 白俊雨, 宋志翔, 苏凌等. 2012. 基于Xu-White模型横波速度预测的误差分析. 地球物理学报, 55(2): 589-595, doi: 10.6038/j.issn.0001-5733.2012.02.021 . |
[32] | 李维新, 王红, 姚振兴等. 2009. 基于约束条件横波速度反演和流体替代. 地球物理学报, 52(3): 785-791. |
[33] | 李晓明, 陈双全, 李向阳. 2012. 利用多分量地震数据反演近地表横波速度. 石油地球物理勘探, 47(4): 532-536. |
[34] | 孙福利, 杨长春, 麻三怀等. 2008. 横波速度预测方法. 地球物理学进展, 23(2): 470-474. |
[35] | 王保丽, 印兴耀, 张繁昌. 2005. 弹性阻抗反演及应用研究. 地球物理学进展, 20(1): 89-92. |
[36] | 吴志华. 2012. 含流体孔隙介质岩石物理分析方法研究[硕士论文]. 青岛: 中国石油大学. |
[37] | 熊晓军, 林凯, 贺振华. 2012. 基于等效弹性模量反演的横波速度预测方法. 石油地球物理勘探, 47(5): 723-727. |
[38] | 杨培杰. 2008. 地震子波盲提取与非线性反演[博士论文]. 东营: 中国石油大学. |
[39] | 杨培杰, 印兴耀. 2008. 非线性二次规划贝叶斯叠前反演. 地球物理学报, 51(6): 1876-1882. |
[40] | 杨培杰, 王长江, 毕俊凤等. 2015. 可变点约束叠前流体因子直接提取方法. 地球物理学报, 58(6): 2188-2200, doi: 10.6038/cjg20150631. |
[41] | 印兴耀, 张世鑫, 张峰. 2013. 针对深层流体识别的两项弹性阻抗反演与Russell流体因子直接估算方法研究. 地球物理学报, 56(7): 2348-2390, doi: 10.6038/cjg20130724. |
[42] | 张广智, 李呈呈, 印兴耀等. 2012. 基于修正Xu-white模型的碳酸盐岩横波速度估算方法. 石油地球物理勘探, 47(5): 717-722. |