一种基于多种地震属性的地层孔隙压力预测方法技术领域
本发明属于油气地球物理勘探技术领域，特别涉及一种基于多种地震属性的砂泥岩地层孔隙压力预测方法。
背景技术
地层孔隙压力是石油勘探开发工作中基础数据之一，地层孔隙异常压力(特别是异常高压)的预测研究，愈来愈引起人们注意。在油气地球物理勘探中，地层孔隙压力为油气的分布、运移、储集提供了重要信息；石油钻井工程中，地层孔隙压力不仅是确定钻井液密度和井深结构的依据，而且关系到能否安全、快速、经济的钻井，甚至会影响钻井的成败。
总的来说，地层孔隙压力预测方法可分为两类，一类是利用测井资料进行地层孔隙压力预测，另一类是运用地震层速度进行地层孔隙压力预测。利用测井资料的地层孔隙压力预测方法的预测精度较高，且测井资料受人为影响因素较少，是公认理想的地层孔隙压力预测方法。目前，利用测井资料的地层孔隙压力预测方法主要有等效深度法、Eaton法，Stone法、综合参数法等。
基于综合参数法的地层孔隙压力预测方法的主要原理如下：
Eberhart-Phillips,D.,Han,D.-H.andZoback,M.D.(1989.Empiricalrelationshipsamongseismicvelocity,effectivepressure,porosity,andclaycontentinsandstone.Geophysics,54(1):82-89)根据大量试验数据分析，影响地层砂泥岩中声波传播速度的因素主要有孔隙度、泥质含量和有效应力，从而构建纵波速度的经验模型。
该模型描述了孔隙度、有效应力和泥质含量对地层岩石中声波速度的综合影响规律:声波速度随孔隙度和泥质含量的增加而减小,随垂直有效应力的增加而增大。这与地层岩石对声波速度测井的响应规律一致。纵波速度经验模型涉及多种测井数据，因此基于上述纵波速度经验模型的地层孔隙压力预测方法命名为综合参数法。重新拟合上述纵波速度经验模型系数并进行推广，获得如下模型：
式(1)中，A0、A1、A2、A3、D为模型系数。
综合参数法首先利用相关的测井数据和实测压力数据确定模型系数A0、A1、A2、A3和D；然后，利用孔隙度测井数据确定目标段的孔隙度利用自然伽玛或自然电位测井数据确定泥质含量Vsh，利用声波时差测井数据确定目标段纵波速度Vp；接着，基于模型(1)计算有效应力Pe；最后，根据地层孔隙压力定义计算目标段的地层孔隙压力Pf。
基于地震资料的地层孔隙压力预测方法主要包括：(1)云美厚(地震地层压力预测[J].石油地球物理勘探,1996,31(4)：575～586)利用单一属性纵波速度计算地层孔隙压力的Fillippone公式法及其改进；(2)石万忠,何生,陈红汉等(多地震属性联合反演在地层压力预测中的应用[J].石油物探,2006年11月,第45卷第6期)利用井约束反演的波阻抗和瞬时频率联合对库车坳陷内的超压带进行预测。
云美厚提出的地层孔隙压力预测方法基于二维或三维地震数据实现，具体为：反演获得层速度，依据Gardner公式将反演层速度转换为密度，并计算上覆地层平均密度ρ，从而获得上覆地层压力。通过给定或搜索目标范围内的最大速度(即基质速度)和最小速度(即孔隙流体速度)，最终根据公式(1)所示的Fillippone公式计算地层孔隙压力Pf。
式(2)中，Pf为地层孔隙压力，单位：MPa；h表示上覆地层深度，单位：m；ρ为上覆地层平均密度，单位：g/cm3；g为重力加速度；Vmax为最大层速度，即岩层有效孔隙度近于零时的纵波速度，单位：m/s；Vmin为最小层速度，即岩层刚性近于零时的纵波速度，单位：m/s；Vi为第i层的层速度，单位：m/s。
石万忠等提出的地层孔隙压力预测方法也是基于二维或三维地震数据实现，具体为：首先，根据瞬时频率的低频特性圈定超压带的范围。然后，进行超压预测，构建各超压带内波阻抗与过剩压力的关系，这样一方面可剔除岩性对波阻抗的影响，另一方面可实现不同成因超压的预测。接着，根据各超压带过剩压力与波阻抗的统计关系，求取各超压带内的过剩压力。最后，根据各超压带内过剩压力绘制压力分布图。
上述云美厚提出的地层孔隙压力预测方法仅利用了地震层速度，紧密依靠纵波速度，通过建立纵波速度与地层孔隙压力间的关系式(即Fillippone公式)进行地层孔隙压力预测，预测效果受反演纵波速度的影响。上述石万忠等提出的地层孔隙压力预测方法仅利用了泊松比，低瞬时频率带仅起到对超压带的横向约束作用。本发明是通过建立泊松比与地层孔隙压力间的统计关系进行地层孔隙压力预测，反演泊松比对其预测效果影响较大，而泊松比为速度与密度的乘积。综上，上述两种方法均是基于建立单一地震属性与的地层孔隙压力预测，且预测效果均与地震层速度有极大相关性，一旦反演地震层速度不准确，将会导致地层孔隙压力预测的不准确。
发明内容
针对现有技术存在的不足，本发明提供了一种基于多种地震属性的砂泥岩地层孔隙压力预测方法，其间接构建了地层孔隙压力与多种地震属性间的非线性关系，从而实现地层孔隙压力预测精度的提高。
为解决上述技术问题，本发明采用如下技术方案：
一种基于多种地震属性的砂泥岩地层孔隙压力预测方法，包括步骤：
步骤1，根据测井所得数据，选取其中的砂泥岩层段的纵波速度V、孔隙度φ、伽玛值GR、密度ρ、深度h和地层孔隙压力Pf。
根据上覆地层压力的计算公式计算得到上覆地层压力。其中，h表示上覆地层深度，ρ(h)为深度为h时所对应的密度，g为重力加速度。
根据压实平衡方程计算有效应力Pe，Pe＝Pov-βPf，Pf为测井所得地层孔隙压力，β为Biot系数。其计算公式我们采用基于临界孔隙度的线性Biot系数模型：
其中φc为临界孔隙度，此处砂泥岩地层我们取值0.4。
步骤2，利用砂泥岩层的纵波速度V、密度ρ，拟合Gardner公式。如(2)式所示：
ρ＝αVβ(2)
其中ρ为密度，V为纵波速度，α、β分别为拟合的系数和指数(单位v＝m/s，ρ＝g/cm3)。
步骤3，利用Vc++6.0进行编程，构建了自适应逐步向前多元回归的算法，首先将步骤1中的变量纵波速度V、孔隙度φ、伽玛值GR代入(3)式；接着将公式(3)中D值范围设置为从0到40，搜索间隔为0.01；再用最小二乘法将超定方程组化成正定方程组；然后用乔里斯基分解法解此方程组；最后依据(4)式所示目标函数求取不同D值反演所得的复相关系数r；当目标函数取得最大值时，A0、A1、A2、A3和D值即为反演所得系数。
其中V为测井所得纵波速度，GR为测井所得伽玛值，GRmax、GRmin分别为所取井数据伽玛值中的最大、最小值。
其中n为步骤1所取测井数据个数，V实为测井实测纵波速度，V拟为公式(3)反演纵波速度，V平为测井所取纵波速度的平均值。
步骤4，根据二维或三维地震数据，反演孔隙度φ'、伽马GR'、纵波速度V1'和V2'的二维或三维数据。其中V1'为使用高精度波阻抗反演所得目的砂泥岩层纵波速度，V2'为使用DIX公式反演所得整个数据体的背景纵波速度。
步骤5，在二维或三维伽马数据体中，对所有道循环寻找最大伽马值GR'max和最小伽马值GR'min，然后再对每一个采样点循环，读取该点的孔隙度φ'、纵波速度V1'和伽玛值GR'。并带入(5)式中，根据牛顿迭代法求解非线性方程得到有效应力P'e。
步骤6，利用(6)式所示的Gardner公式，计算每一个采样点所对应的密度ρ’，再利用(7)式计算P'ov
ρ'＝αV2'β(6)
步骤7，根据压实平衡方程求取地层孔隙压力Pf',其中βPf'＝P'ov-Pe'，β取值参考步骤1中的公式(1)；
步骤8，最后对反演所得二维或三维地层孔隙压力Pf进行二维或三维中值滤波。
字母加’为地震数据，没加是测井数据。
和现有技术相比，本发明具有如下优点和有益效果：
1、综合多地震属性，利用反演的纵波速度进行地层压力预测时，扣除了其他影响因素，例如孔隙度、泥质含量等，从而提高了地层孔隙压力的预测精度。
2、创造性的使用两种计算方法反演的纵波速度，高精度波阻抗反演所得目的砂泥岩层纵波速度和使用DIX公式反演的背景纵波速度。前者能较为准确反演目的层内纵波速度，从而较为准确计算目的层内的有效应力Pe，后者能从宏观上反映上纵波速度与密度的关系，从而较为准确的计算上覆地层压力Pov。
具体实施方式
本发明可以针对二维地震数据和三维地震数据进行处理，具体步骤包括：
步骤1，根据测井所得数据，选取其中的砂泥岩层段的纵波速度V、孔隙度φ、伽玛值GR、密度ρ、深度h和地层孔隙压力Pf。
根据上覆地层压力的计算公式计算得到上覆地层压力。其中，h表示上覆地层深度，ρ(h)为深度为h时所对应的密度，g为重力加速度。
根据压实平衡方程计算有效应力Pe，Pe＝Pov-βPf，Pf为测井所得地层孔隙压力，β为Biot系数。其计算公式我们采用基于临界孔隙度的线性Biot系数模型：
其中φc为临界孔隙度，此处砂泥岩地层我们取值0.4。
步骤2，利用砂泥岩层的纵波速度V、密度ρ，拟合Gardner公式。如(2)式所示：
ρ＝αVβ(2)
其中ρ为密度，V为纵波速度，α、β分别为拟合的系数和指数(单位v＝m/s，ρ＝g/cm3)。
步骤3，利用Vc++6.0进行编程，构建了自适应逐步向前多元回归的算法，首先将步骤1中的变量纵波速度V、孔隙度φ、伽玛值GR代入(3)式；接着将公式(3)中D值范围设置为从0到40，搜索间隔为0.01；再用最小二乘法将超定方程组化成正定方程组；然后用乔里斯基分解法解此方程组；最后依据(4)式所示目标函数求取不同D值反演所得的复相关系数r；当目标函数取得最大值时，A0、A1、A2、A3和D值即为反演所得系数。
其中V为测井所得纵波速度，GR为测井所得伽玛值，GRmax、GRmin分别为所取井数据伽玛值中的最大、最小值。
其中n为步骤1所取测井数据个数，V实为测井实测纵波速度，V拟为公式(3)反演纵波速度，V平为测井所取纵波速度的平均值。
步骤4，根据二维或三维地震数据，反演孔隙度φ'、伽马GR'、纵波速度V1'和V2'的二维或三维数据。其中V1'为使用高精度波阻抗反演所得目的砂泥岩层纵波速度，V2'为使用DIX公式反演所得整个数据体的背景纵波速度。
步骤5，在二维或三维伽马数据体中，对所有道循环寻找最大伽马值GR'max和最小伽马值GR'min，然后再对每一个采样点循环，读取该点的孔隙度φ'、纵波速度V1'和伽玛值GR'。并带入(5)式中，根据牛顿迭代法求解非线性方程得到有效应力P'e。
步骤6，利用(6)式所示的Gardner公式，计算每一个采样点所对应的密度ρ’，再利用(7)式计算P'ov
ρ'＝αV2'β(6)
步骤7，根据压实平衡方程求取地层孔隙压力Pf',其中βPf'＝P'ov-Pe'，β取值参考步骤1中的公式(1)；
步骤8，最后对反演所得二维或三维地层孔隙压力Pf进行二维或三维中值滤波。
字母加’为地震数据，没加是测井数据。
其中，目标函数：把判断若干可行设计方案优劣的评价标准用涉及变量的函数形式来表示，所述的涉及变量的函数形式即目标函数。衡量可行设计方案优劣的标准一般根据实际情况确定，可以是在满足一定约束条件下求目标函数的极值。
自适应：在新的环境或运行条件下，依据某种判别结果适当改变反演公式的运行参数，从而保持对新事物的良好预测。
地层孔隙压力：地层岩石孔隙内流体所承受的压力。
上覆地层压力：覆盖在地层以上的岩石及其岩石孔隙中流体的总重量造成的压力。
有效应力：即骨架应力，固体矿物之间的接触应力(或支撑应力)。一般有效应力等于上覆地层压力减去等效孔隙压力。
孔隙度：岩样中所有孔隙空间体积之和与岩样体积的比值。
伽玛测井：沿井身测量岩层的天然伽马射线强度的方法。可以划分出钻孔的地质剖面、确定砂泥岩剖面中砂岩泥质含量和定性地判断岩层的渗透性。
BISQ系数：相同孔隙压力下，孔隙体积变化与总体积变化之比。
上述具体实施例用来解释本发明，而不是对本发明进行限制，在本发明的精神和权利要求的保护范围内，对本发明做出的任何修改和改变，都落入本发明的保护范围。