技术领域
[0001] 本
发明涉及一种随钻预测地层孔隙压力的方法，更具体地说，涉及一种利用
地震资料随钻预测地层孔隙压力的方法。
背景技术
[0002] 在石油钻井过程中，不确定性的异常地层孔隙压力经常造成井下作业复杂，且容易发生事故，给钻井作业的人力、物力带来不同程度的损失。 为此如何在钻井前事先预测地层孔隙压力是非常必要的，因若在钻井前能预测到地层孔隙压力就可以科学设计
套管程序，合理选择
钻井液密度，以实现快速、安全、高效、低成本钻井和完井。 [0003] 为此本发明中的创作人凭借其多年从事相关行业的经验与实践，并经潜心研究与开发，终于创造出一种利用地震资料来随钻预测地层孔隙压力的方法。 发明内容
[0004] 本发明的主要目的在于提供一种利用地震资料随钻预测钻头底下地层孔隙压力的方法。
[0005] 本发明中利用地震资料随钻预测邻井钻头底下地层孔隙压力的方法，包括下列步骤：
[0006] 步骤1，针对同处一区
块的待钻井和已钻井，分别提取各自井旁若干道地震记录，作加权处理获得待钻井和已钻井的地震记录；
[0007] 步骤2，利用已钻井不同层系地层的
测井数据和地震记录，建立地震记录预测地层
声波时差测井曲线的分层模型；
[0008] 步骤3，预测待钻井钻头底下地层的声波时差测井曲线；
[0009] 步骤4，利用已钻井的地层密度测井数据、声波时差测井数据和自然伽
马测井数据建立该地质区块的上覆地
层压力计算模型、孔隙度与泥质含量计算模型和分层段的地层孔隙压力计算模型；
[0010] 步骤5，结合已建立的本区块的上覆地层压力计算模型、该层系的孔隙度与泥质含量计算模型和地层孔隙压力计算模型预测当前钻遇地层的孔隙压力。 [0011] 所述步骤2包括：
[0012] 步骤21，提取已钻井的声波时差测井数据、地层密度测井数据和自然伽马测井数据；
[0013] 步骤22，提取已钻井相邻若干道地震记录，作加权处理获得已钻井地震记录； [0014] 步骤23，查阅已钻井的地质录井资料，获取准确的地质分层数据(地层地质分层的顶部深度和底部深度)；
[0015] 步骤24，取处理完毕的已钻井每一地质分层层段的地震记录； [0016] 步骤25，分层段求取已钻井的关联维数、R/S分数维、振幅谱维、最大李雅普洛夫指数、突变参数和综合反射地震参数序列；
[0017] 步骤26，结合已钻井地质资料和测井资料，得到已钻井分层段的测井曲线； [0018] 步骤27，以关联维数、R/S分数维、振幅谱维、最大李雅普洛夫指数、突变参数和综合反射地震参数为
输入层的神经元，以声波时差值为
输出层的神经元，分层段创建地震参数与声波时差关系的反向传播神经网络模型；
[0019] 步骤28，分层段训练神经网络模型，获得分层段预测测井曲线的反向传播神经网络模型。
[0020] 所述步骤3包括：
[0021] 步骤31，根据地质录井数据，确定当前钻井钻头所处的地质层段； [0022] 步骤32，选取该层段的声波时差测井曲线预测的神经网络模型； [0023] 步骤33，将该层段的地震参数代入已训练的神经网络模型，预测出该层系井段的声波时差测井曲线。
[0024] 所述步骤4包括：
[0025] 步骤41，建立区域上覆地层压力梯度计算模型；
[0026] 步骤42，建立区域孔隙度和泥质含量计算模型；
[0027] 步骤43，建立区域分层孔隙压力计算模型。
[0028] 所述步骤41包括：
[0029] 1)利用密度测井资料计算不同深度的上覆地层压力梯度
[0030]
[0031] 其中G0i为某深度的上覆地层压力梯度(单位：g/ml)，ρi，Δhi分别为该深度上部各层的密度平均值(单位：g/ml)和层间深度间隔(单位：km)，ρ0，h0分别为上部无密度测井数据层段的平均密度(单位：g/ml)和厚度(单位：km)，
[0032] 2)利用下式建立上覆地层压力梯度和深度的非线性关系模型
[0033]
[0034] 其中，G0为上覆地层压力梯度值(单位：g/ml)，h为测井深度值(单位： km)，a1，a2，a3，a4为模型参数。
[0035] 所述步骤42包括：
[0036] 利用下式建立孔隙度、泥质含量与纵横
波速度的非线性关系模型 [0037] φ＝a1+a2vp+a3vs
[0038] Vsh＝b1+b2vp+b3vs
[0039] 其中φ为孔隙度，由声波测井资料计算得到；Vsh为泥质含量，由自然伽马测井资料计算得到；a1，a2，a3，b1，b2，b3为模型参数，vs为横波速度(单位km/s)，可以利用纵波速度vp由下面的经验关系确定
[0040]
[0041] 纵波速度可以直接由声波时差测井数据计算出
[0042]
[0043] 式中，vp为纵波速度(单位：km/s)，Δt为声波时差测井值(单位：μs/ft)。 [0044] 所述步骤43包括：
[0045] 1)建立纵波速度、泥质含量、孔隙度以及垂直有效
应力之间非线性关系的分层模型
[0046]
[0047] 上式中vp为纵波速度(单位：km/s)，Vsh为泥质含量，φ为孔隙度，Pe为垂直有效应力(单位MPa)，a1，a2，a3，a4，a5为模型参数。
[0048] 2)利用有效应力定律预测地层孔隙压力
[0049] Pp＝Po-Pe
[0050] 其中，Pe为垂直有效应力(单位MPa)，Po为上覆地层压力(单位MPa)，Pp为地层孔隙压力(单位MPa)。
[0051] 本发明中利用地震资料随钻预测邻井地层孔隙压力的方法主要是找出同 一井眼中不同层系地层的已钻井段测井结果与地震特征参数的关联，而根据地质分层随钻预测地层孔隙压力，以有效防止井下复杂情况的发生。
附图说明
[0052] 图1是已钻井通过相邻若干道地震记录加权处理获得的地震记录； [0053] 图2是利用已钻井资料建立的该地质区块地震双程旅行时和测井深度之间的对应关系；
[0054] 图3是已钻井第三系安居安组井段关联维的计算结果图；
[0055] 图4为已钻井第三系安居安组井段R/S分数维的计算结果图；
[0056] 图5为已钻井第三系安居安组井段振幅谱维的计算结果图；
[0057] 图6为已钻井第三系安居安组井段最大李雅普洛夫指数的计算结果图； [0058] 图7为已钻井第三系安居安组井段突变参数的计算结果图；
[0059] 图8为已钻井第三系安居安组井段综合反射地震参数的计算结果图； [0060] 图9为根据已钻井测井资料和地质资料提取的已钻井第三系安居安组井段的声波时差测井曲线；
[0061] 图10是待钻井通过相邻若干道地震记录加权处理获得的地震记录； [0062] 图11是待钻井第三系安居安组井段关联维的计算结果图；
[0063] 图12为待钻井第三系安居安组井段R/S分数维的计算结果图；
[0064] 图13为待钻井第三系安居安组井段振幅谱维的计算结果图；
[0065] 图14为待钻井第三系安居安组井段最大李雅普洛夫指数的计算结果图； [0066] 图15为待钻井第三系安居安组井段突变参数的计算结果图；
[0067] 图16为已钻井第三系安居安组井段综合反射地震参数的计算结果图； [0068] 图17为待钻井所预测的第三系安居安组井段的声波时差测井曲线； [0069] 图18是根据已钻井密度测井资料建立的该地质区块的上覆地层压力计算模型； [0070] 图19为待钻井第三系安居安组井段的地层孔隙压力；
[0071] 图20为本发明预测钻头底下地层孔隙压力的方法
流程图。
具体实施方式
[0072] 以下参照附图20对本发明中利用地震记录随钻预测井眼待钻井段孔隙压力的方法进行详细说明，主要包括下列步骤：
[0073] 步骤1，对已钻井和待钻井分别提取各自井旁若干道地震记录，作加权处理获得已钻井和待钻井地震记录(图1)。
[0074] 步骤2，对已钻井进行声波时差、
中子密度测井和自然伽马测井，获得声波时差测井曲线、地层密度测井曲线和自然伽马测井曲线。
[0075] 步骤3，建立地震特征参数与测井数据间的关系，即建立利用地震记录构建测井曲线的模型：
[0076]
地震波的反射特征与地下介质的波阻抗有关，波阻抗的自然对数对时间增量的导数定义为反射率函数，地震记录的反射率函数深刻反映了地震的线性与非线性特征；而地层的声波速度直接决定着反射率函数。 因此利用地震的线性与非线性特征可以映射地震的反射率函数，也就是说地震的线性与非线性特征通过适当的映射模型可以反映地层的声波特征。
[0077] 其中包括：
[0078] 1)提取已钻井的测井数据：声波时差测井数据和地层密度测井数据； [0079] 2)查阅已钻井的地质录井资料，获取准确的地质分层数据：地层地质分层的顶部深度和底部深度；
[0080] 3)取已处理完毕的已钻井的每一地质分层层段的地震记录；
[0081] 4)利用所述每一地质分层层段的地震记录，分层段求取地震特征参数；其中，所述地震特征参数包括：关联维数、R/S分数维、振幅谱维、最大李雅普洛夫指数、突变参数和综合反射地震参数；其中，
[0082] a.用声波测井时差进行深时转换，其转换公式为
[0083]
[0084] 式中为H0起始深度，H为终止深度，T(H)为深度H所对应的时间，Δt(h)为深度h处的声波测井的时差值。 将测井的深度
采样点变成了时间序列，然后根据地震记录时间采样间隔对声波、密度测井值进行重采样，将深度采样值转换为时间采样值，以此时深关系(图2)为标准来实现地震和测井数据在时间和深度上的对应，将地震记录的时间刻度转换为深度刻度，得到分层段的地震记录，用于提取地震特征参数。 [0085] b.关联维数
[0086] 从时间序列直接计算关联维数的
算法，即G-P算法。
[0087] 第一步是将地震记录的时间序列xi＝x(ti)i＝1，2，......，N重新排列创立一个m维的向量相空间
[0088] Xn(m，τ)＝(xn，xn+τ，......，xn+(m-1)τ) n＝1，2，......，Nm (2) [0089] 式中τ＝kΔt是固定时间间隔，即时间延迟，Δt是两次相邻采样的间隔，k是整数
[0090] Nm＝N-(m-1)τ (3)
[0091] 从这Nm个点中任意
选定一个参考点Xi，计算其余Nm-1个点到Xi的距离 [0092]
[0093] 对所有点重复这一过程，得到关联积分函数
[0094]
[0095] 式中H是Heaviside函数
[0096]
[0097] 在lgCm(r)-lgr的坐标图上，如点分布在一条直线附近，则分形存在，采用直线拟合的方法，拟合直线的斜率就是关联维D2。 图3为关联维的计算结果图。 [0098] c.R/S分数维
[0099] 重标度极差分析方法(R/S分析法)是Hurst于1965年提出的一种新的统计方法，它在分形理论中有着广泛的应用。
[0100] 设已知地震记录的时间序列为xi＝x(ti)，i＝1，2，......，N，则τ个时间数据的均值为
[0101]
[0102] 由此可以求得累积离差
[0103]
[0104] 极差
[0105]
[0106] 标准差
[0107]
[0108] 根据Hurst分析的统计规律的关系式为H
[0109] R/S∝(τ/2) (10)
[0110] 具体计算时，计算出R/S的一组值。在lg(R/S)-lg(τ/2)的坐标图上拟合 出一条直线，其斜率就是R/S分数维。
[0111] 图4为R/S分数维的计算结果图。
[0112] d.振幅谱维
[0113] 对地震记录的时间序列xi＝x(ti)，i＝1，2，......，N作
功率谱分析，功率谱S(f)随着f的变化可以表示为S(f)∝fβ。因为振幅谱 所以有A(f)∝fβ/2。令Df＝β/2，则有A(f)∝fDf。 绘制一张lgA(f)-lgf的图，然后用直线拟合，其斜率即为振幅谱维Df。
[0114] 图5为振幅谱维的计算结果图。
[0115] e.李雅普洛夫指数
[0116] 李雅普洛夫指数是指相空间中邻近轨道发散或收敛的平均指数率，反映了系统的性态对于初值的敏感程度。
[0117] 一般计算最大李雅普洛夫指数采用A.Wolf重构法，其步骤如下： [0118] A.重建相空间
[0119] 设已知地震记录的时间序列为xi＝x(ti)，i＝1，2，......，N，i为采样点的序号，xi为在该时间下的振幅值，输入嵌入空间维数m及延迟时间τ得相空间： [0120] Xn(m，τ)＝(xn，xn+τ，......，xn+(m-1)τ) n＝1，2，......，Nm (11) [0121] Nm＝N-(m-1)τ (12) [0122] B.在相空间中，以初始点A(t0)为参考点，选取A(t0)的最近邻点B(t0)，设在t1＝t0+kΔt时，A(t0)和B(t0)分别演化到A(t1)和B(t1)，计算从t0到t1时的指数增长率 [0123]
[0124] 式中， 分别为相空间A(t1)与B(t1)，A(t0)与B(t0)两点之间的距离。
[0125] C.在A(t1)的若干邻近点中，找出一个与A(t1)的夹
角θ1很小的邻近点C(t1)，如果找不到，仍然选取B(t1)，设在t2＝t1+kΔt时，A(t1)和C(t1)分别演化到A(t2)和C(t2)，则
[0126]
[0127] 将这一过程一直进行到点集的终点，然后取的平均值作为最大李雅普洛夫指数的估计值LE(m)。
[0128] D.增加嵌入空间维数m，重复A-C步，直到LE(m)保持平稳为止，此时的LE(m)即为所求的最大李雅普洛夫指数。
[0129] 图6为最大李雅普洛夫指数的计算结果图。
[0130] f.突变参数
[0131] 将地震
信号看成对时间变量的连续函数x(t)，x(t)可展成级数形式如下 t 2 n
[0132] y＝x(t)＝a0+a1+a2t+…+ant+… (14)
[0133] 式中，t为时间，y为对应t的位移，a0，a1，a2， …an为待定的系数。 实际分析发现，对具有一定趋势规律的时间序列，截取到4次项时，
精度已足够高。 则对上式可近似表示为t 2 3 4
[0134] y＝x(t)＝a0+a1+a2t+a3t+a4t (15)
[0135] 对上式作变量代换，化为尖点突变的标准形式，先令
[0136] t＝zt-q (16)
[0137] 将其代入上式，可得
[0138]4 3 2
[0139] b0＝a4q-a3q+a2q-a1q+a0
[0140] 式中， b1＝-4a4q3+3a3q2-2a2q+a1
[0141] b2＝6a4q2-3a3q+a2
[0142] b4＝a4
[0143] 上面仍不是尖点突变的标准形式，作进一步变量代换，令
[0144]
[0145] 则
[0146]
[0147] 式中：c＝b0
[0148]
[0149]
[0150] z为状态变量，a，b为控制变量，由突变理论可知，平衡曲面方程为3
[0151] z+az+b＝0 (20)
[0152] 分叉集方程为3 2
[0153] 4a+27b ＝0 (21)
[0154] 只有在控制变量满足分叉集方程时，系统才是不稳定的，才有可能从一个平衡态突变到另一个平衡态。
[0155] 图7为突变参数的计算结果图。
[0156] g.综合反射地震参数
[0157] 将常规的反射地震参数，即瞬时振幅、瞬时
频率、瞬时
相位、自相关函数极小值、自相关函数极大值以及五阶
自回归模型系数共十种参数进行综合评判，优选出一种综合反射参数参与建模。
[0158] 图8为综合反射地震参数的计算结果图。
[0159] 5)根据已钻井的地质资料和测井资料，提取已钻井分层段的声波时差测井曲线(图9)。
[0160] 6)将以上提取的各种地震特征参数以及对应的声波时差测井数据输入神经网络，分层创建反向传播神经网络模型。
[0161] 步骤4，利用待钻井段的地震特征参数预测待钻井段的声波时差测并曲线： [0162] 1)根据地质录井数据，确定待钻井钻头所处的地质层段；
[0163] 2)选取该层段的声波时差测井曲线神经网络模型；
[0164] 3)将根据待钻井的地震记录(图10)计算出的待钻井段的关联维数(图11)、R/S分数维(图12)、振幅谱维(图13)、最大李雅普洛夫指数(图14)、突变参数(图15)和综合反射地震参数(图16)代入已训练好的对应层段神经网络模型，预测出待钻井在不同深度地层的声波时差测井曲线(图17)。
[0165] 步骤5，利用已钻井的地层密度测井数据建立该地质区块的上覆地层压力计算模型：
[0166] 上覆地层压力计算模型首先要求确定上覆地层压力梯度，上覆地层压力梯度可以由密度测井的散点数据计算得到
[0167]
[0168] 上式中，G0i为某深度的上覆地层压力梯度(单位：g/ml)，ρi，Δhi分别为该深度上部各层的密度平均值(单位：g/ml)和层间深度间隔(单位：km)，ρ0，h0分别为上部无密度测井数据层段的平均密度(单位：g/ml)和厚度(单位：km)。 [0169] 由上覆地层压力梯度可计算各深度的上覆地层压力
[0170] P0i＝0.01G0ihi (23)
[0171] 上式中，P0i为地下某点上覆地层压力(单位：MPa)，G0i为该点上覆地层压力梯度(单位：g/ml)，hi为该点深度(单位：m)。
[0172] 通常利用已钻井密度测井的散点数据将上覆地层压力梯度回归为深度的函数，研究认为将上覆地层压力梯度回归为下式效果最好
[0173]
[0174] 为了克服通常密度测井数据不足的问题，在只有n个密度测井散点数据的情况下，可以建立变量W和G0的关系
[0175] Wi＝0 (i＝1)
[0176]
[0177]
[0178] 根据以上关系，将W回归为深度h的函数
[0179]
[0180] 利用上式进行非线性回归，得到a1，a2，a3，a4代入(24)式中进行计算即可以得到随深度变化的上覆地层压力梯度剖面(图18)，利用该剖面可以方便地计算该地质区块任一深度的上覆地层压力。
[0181] 步骤6，利用已钻井的声波时差测井曲线和自然伽马测井曲线建立该地质区块的孔隙度和泥质含量计算模型：
[0182] 大量实验证明，声波速度、孔隙度和泥质含量之间存在良好的线性关系，这种关系随地质区块的差异而有所不同，所以可以将声波速度、孔隙度和泥质含量代入下式，建立适合于本区块的孔隙度和泥质含量计算模型
[0183] φ＝a1+a2vp+a3vs (27)
[0184] Vsh＝b1+b2vp+b3vs
[0185] 其中φ为孔隙度，Vsh为泥质含量，vs为横波速度(单位km/s)，可以利用纵波速度vp由下面的经验关系确定
[0186]
[0187] 纵波速度可以直接由声波时差测井数据计算出
[0188]
[0189] 式中，vp为纵波速度(单位：km/s)，Δt为声波时差测井值(单位：μs/ft)， [0190] 在获取已钻井的声波时差测井数据和自然伽马测井数据的
基础上，就可以利用测井解释的结论来求取孔隙度和泥质含量。 通常选取声波时差测井数据由下式可以计算孔隙度
[0191]
[0192] 式中，Δt，Δtma，Δtf分别为目的层、
岩石骨架和地层
流体的声波时差测井值(单位：μs/ft)。
[0193] 在获取自然伽马测井数据的基础上利用相对值法计算泥质含量 [0194]
[0195] 式中，Vsh为泥质含量，GCUR为希尔奇指数，与地质时代有关，可根据取芯分析资料和自然伽马测井值统计确定，一般对于第三系地层取3.7，老地层取2。IGR为泥质含量指数， GR，GRmin，GRmax分别表示目的层、纯泥岩层和纯砂岩层的自然伽马测井值(API)。
[0196] 将计算出的声波速度、孔隙度和泥质含量按照(27)式通过多元线性回归的方法确定a1，a2，a3，b1，b2，b3，从而得到该区块的孔隙度和泥质含量计算模型。 [0197] 步骤7，利用已钻井资料建立该区块分层段的地层孔隙压力计算模型： [0198] 地层孔隙压力预测的理论基础是饱和多孔介质的有效应力定律 [0199] Pe＝Po-Pp (32)
[0200] 其中，Pe为垂直有效应力(单位MPa)，Po为上覆地层压力(单位MPa)，Pp为地层孔隙压力(单位MPa)。从上式可以知道，在已知上覆地层压力的情况下，只要求得垂直有效应力就可以计算出地层孔隙压力，所以建立地层孔隙压力计算模型实际上就是建立垂直有效应力计算模型。
[0201] 大量的
岩心测试实验已经证明，影响声波在岩石中传播速度的主要因素 [0202] 是岩性、孔隙度和垂直有效应力，所以可以将纵波速度表示为岩性、孔隙度和垂直有效应力的函数
[0203] vp＝f(Vsh，φ，Pe) (33)
[0204] 上式中vp为纵波速度(单位：km/s)，Vsh为泥质含量，φ为孔隙度，Pe为垂直有效应力(单位MPa)。
[0205] 实验工作表明建立以下模型来表征上面的函数关系最为理想
[0206]
[0207] 在建立分层垂直有效应力计算模型时，vp的值由声波时差测井资料直接计算，孔隙度和泥质含量的计算利用已建立的孔隙度和泥质含量分层计算模型。 [0208] 利用已建立的该区块的上覆地层压力模型和静液压力或实测压力就确定用于建模的分层垂直有效应力
[0209] Pe＝Po-gwh (35)
[0210] 式中，gw表示静液压力梯度，一般取0.0105MPa/m。
[0211] 将vp，φ，Vsh，Pe分层数据代入式(34)中，利用多元非线性回归的方法计算出a1，a2，a3，a4，a5，从而建立起该地质区块分层的地层孔隙压力计算模型。 [0212] 在实际应用中还经常使用以下简化的模型，即忽略孔隙度和泥质含量对垂直有效应力的影响，而只考虑声波速度对垂直有效应力的影响，在这种情况下，就无需建立孔隙度和泥质含量计算模型，这就大大简化了计算程序。 这种简化计算模型如下式所示 [0213]
[0214] 步骤8，利用预测得到待钻井段的声波时差测井数据结合已建立的区域上覆地层压力计算模型、孔隙度与泥质含量计算模型和分层段的地层孔隙压力 计算模型计算当前钻遇层段的地层孔隙压力：
[0215] 实际钻进中，通过岩屑录井资料和分层神经网络
预测模型随钻预测出钻头下地层的声波时差测井数据，将其转换为声波速度，结合分层孔隙度与泥质含量计算模型计算待钻地层的孔隙度与泥质含量，再利用分层地层孔隙压力计算模型反算待钻地层的垂直有效应力，最后利用区域上覆地层压力计算模型计算待钻地层的孔隙压力(图19)。 [0216] 上述
实施例仅用于说明本发明，而并非用于限定本发明。
