一种基于井筒裂缝观测的三维裂缝渗透率张量计算模型的制作方法

本发明涉及油气田勘探开发领域，尤其是一种基于井筒裂缝观测的三维裂缝渗透率张量计算模型。

背景技术：

在裂缝性储层勘探开发过程中，天然裂缝的渗流方向与裂缝的产状、充填性以及现今地应力大小、方向密切相关，裂缝的主渗流方向影响井网部署、注水开发等措施实施。在油藏中，天然裂缝分布极为复杂；在一个单元内，不同期次、不同展布以及不同开度的裂缝在油气导流中所起的作用不一，如何考虑三维的地质因素，输出三维的裂缝渗透率张量是裂缝描述的一大难点。刘月田等(刘月田丁祖鹏屈亚光赵辰军.油藏裂缝方向表征及渗透率各向异性参数计算[J].石油学报,2011,32(5):842-846.)所提的裂缝渗透率计算方法仅考虑了特殊情况下(剪切裂缝，对称，等开度)裂缝渗透张量的计算，并且该方法也无法定量输出主渗透率方向；刘敬寿等(刘敬寿,戴俊生,邹娟,杨海盟,汪必峰,周巨标.裂缝性储层渗透率张量定量计算方法[J].石油与天然气地质,2015,36(6):1022-1029.)所提的裂缝渗透率张量计算方法考虑的三维的因素，只能输出二维的结果，无法准确的计算三维等渗透率面、最大差渗透率面、不同方向的渗透率剖面等。本发明专利基于岩心统计、薄片观测以及测井解释，在应力场模拟的基础上，确定裂缝地下开度，利用古地磁定向、岩层产状法或地层倾角法确定裂缝的三维展布，进而建立真三维裂缝渗透率计算模型，在推导相应算法的基础上，编制相应的程序，实现准确的计算三维等渗透率面、最大差渗透率面、不同方向的渗透率剖面等。

本发明旨在解决上述问题，提供了一种基于井筒裂缝观测的三维裂缝渗透率张量计算模型，它解决了无法准确输出井筒裂缝真三维空间等渗透率面、最大差渗透率面以及不同方向的渗透率剖面的问题。

本发明的技术方案为：一种基于井筒裂缝观测的三维裂缝渗透率张量计算模型，具体步骤如下：

第一步 井筒裂缝发育程度的观测

利用岩心观测、测井资料解释对裂缝的开度测量，确定裂缝的组系特征；依据成像资料、测井或者岩心观察，求取裂缝的线密度，为后面计算裂缝平行渗透率准备。

第二步 裂缝定向

在岩心裂缝观测中，对非定向岩心采用古地磁定向、岩层产状法或地层倾角法确定每条裂缝的产状；在成像资料、声波测井解释裂缝时，利用地层倾角测井或者成像测井资料对裂缝定向。

第三步 应力场模拟确定裂缝地下的开度

利用水力压裂法、波速法、声发射法、应力解除法、光弹性应力法等可以对工区现今地应力方位进行了判断，并通过测井、压裂资料以及物理实验计算了关键井现今地应力的数值，利用地震反演、测井解释、岩石力学实验确定岩石的力学属性。在此基础上，利用软件建立有限元模型，采用多井约束的方法对研究区现今地应力场进行了数值模拟。利用模拟结果，确定井筒裂缝发育处的现今地应力大小、方向；国外专家Willis-Richards等和Jing等同时考虑了正应力及剪应力对裂缝开度的影响，得出在现今地应力场改造后的裂缝开度计算公式：

公式(1)中，b为裂缝的现今地下实际开度，单位：m；b₀为裂缝的原始开度，单位：m；σ′_n为有效正应力，单位：MPa；b_res代表裂缝面承受最大有效正应力时的裂缝开度，单位：m；σ_nref为使裂缝开度减小90％的有效正应力，单位：MPa，是一个与岩性相关的系数，已有学者[Willis-Richards J，1996；Jing Z，1998；Durham WB，1994；Chen Z，2000；秦积舜；2002；]给出了不同条件下的测试数值。

第四步 建立真三维裂缝渗透率张量计算几何模型

如图2所示，设裂缝面Φ的倾角为δ，倾向为ω，以裂缝为参照物建立O-ABC坐标系(O-ABC)，O-ABC坐标系中的三个坐标轴分别对应于裂缝面的法线方向(OA轴)、裂缝走向方向(OB轴)、裂缝面内垂直于裂缝走向线的方向(OC轴)。设任意渗流面ψ的倾角为α，倾向为β，其单位法向矢量为m；以渗流面ψ为参照物建立坐标系(O-XYZ)，其中，OX轴为渗流面ψ的法线方向，OY轴、OZ轴位于渗流面ψ内，OP为渗流面ψ的走向线，HQ为渗流面ψ的倾斜线，HR为渗流面ψ倾斜线的水平投影；线OS、ON、OE、OW、HR、OH、OP、OB在水平面内，定义θ为渗流面ψ内OY轴与OP轴的夹角，通过调整θ的大小，可以求取裂缝在θ中不同方向的渗透率。

第五步 真三维裂缝渗透率张量数学模型

在已知裂缝线密度、地下开度的基础上，单组(单条)裂缝的平行渗透率K可以表示为：

公式(2)中，b为裂缝的现今地下实际开度，单位：m；D_lf为单组裂缝的线密度，条/m。

设第i条裂缝的渗透率张量为K，则其在坐标系O-XYZ中的表达式可表示为：

同理，渗透率张量K在坐标系O-ABC中的表达式可表示为：

OA轴在大地坐标系中的三个分量为：

OX轴在大地坐标系中的三个分量为：

设裂缝的渗透率张量为K，其在O-ABC坐标系中的表达式可表示为：

同理，渗透率张量K则其在O-XYZ坐标系中的表达式可表示为：

在O-XYZ坐标系中，每个单元体内单组裂缝的渗透率张量K_ABC可表示为：

由公式(7)-(9)可以得到：

K_XYZi＝T·K_ABC·T^T (10)

其中，O-ABC坐标系转换为O-XYZ坐标系的旋转矩阵T可表示为：

公式(11)的旋转矩阵T中各参数分别表示为：

在O-XYZ坐标系中，单元体内每组(每条)裂缝的渗透率张量K_XYZi可表示为：

在复杂的地质条件下，储层裂缝经历多期构造应力场演化，每组裂缝的产状、线密度、开度往往不同，依据公式(10)、公式(11)及公式(13)，单元体内发育多组裂缝时，渗透率张量K_XYZ可表示为：

公式(14)中，k为单元体内裂缝的组数；b_i为第i组裂缝的开度，单位：m；T_i为第i组裂缝在O-ABC坐标轴分量转换为O-XYZ坐标轴分量的旋转矩阵，单位：m；D_lfi为第i组裂缝的线密度，单位：条/m。

第六步 编制程序，输出渗透率剖面

利用公式(1)-(14)推导的算法编制相应的程序，通过循环迭代变量α(渗流面ψ的倾角)、β(渗流面ψ的倾向)、θ(OY轴与OP轴的夹角)，其中的α范围为[0-90°]，β的范围为[0-360°)，θ的范围为[0-90°)，根据计算精度需求，设置循环迭代步长(⊿sp，例如：5°，2°，1°，0.5°，0.2°，0.1°等)将α、β、θ逐次累加⊿sp，利用公式(12)-(14)，计算坐标系(O-XYZ)不同的α、β、θ对应的K_XXαβθ(X方向渗透率主值)、K_YYαβθ(Y方向渗透率主值)、K_ZZαβθ(Z方向渗透率主值)。

利用公式(1)-(14)，设置α＝0°，β为[0-360°)的任意数值；调节θ的大小，步长⊿sp，范围为[0-90°)，输出K_YYαβθ(Y方向渗透率主值)、K_ZZαβθ(Z方向渗透率主值)，得到水平面的渗透率变化。

利用公式(1)-(14)，设置α＝90°，β＝0°；调节θ的大小，步长⊿sp，范围为[0-90°)，输出K_YYαβθ(Y方向渗透率主值)、K_ZZαβθ(Z方向渗透率主值)，得到东西向剖面的渗透率变化。

利用公式(1)-(14)，设置α＝90°，β＝90°；调节θ的大小，步长⊿sp，范围为[0-90°)，输出K_YYαβθ(Y方向渗透率主值)、K_ZZαβθ(Z方向渗透率主值)，得到南北向剖面的渗透率变化。

利用公式(1)-(14)，设置α范围为[0-90°]，β的范围为[0-360°)，步长，⊿sp；计算不同的α、β(渗流面ψ确定)的平面内——即不同的渗流面ψ内，K_YYαβθ(Y方向渗透率主值)、K_ZZαβθ(Z方向渗透率主值)，每个渗流面ψ内计算得到的K_YYαβθ、K_ZZαβθ数目为2n，2n个K_YYαβθ、K_ZZαβθ的平均值为K_aver；计算不同的α、β对应的方差函数W_αβ：

公式(15)中，K_YYαβθi、K_ZZαβθi分别为特定的α、β对应的第i个Y方向渗透率主值、Z方向渗透率主值，单位：10^-3μm²；n为每个渗流面ψ计算的Y或Z方向渗透率主值的数目，单位：个；K_aver为2n个K_YYαβθ、K_ZZαβθ的平均值，单位：10^-3μm²。

通过公式(15)，计算不同的α、β对应的方差函数W_αβ，比较不同的α、β的W_αβ。其中，W_αβ最大值对应的渗流面ψ为最大差渗透率面，记录对应的倾角α，倾向β；W_αβ最小值对应的渗流面ψ为等渗透率面，记录对应的倾角α，倾向β。

本发明的有益效果是：本发明专利基于岩心统计、薄片观测以及测井解释，通过现今地应力分析，确定裂缝地下开度，利用古地磁定向、岩层产状法或地层倾角法确定裂缝的三维展布，进而建立真三维裂缝渗透率计算模型，推导可编程算法，编制相应的程序，实现准确的计算三维等渗透率面、最大差渗透率面、不同方向的渗透率剖面等。本发明由严格的数学算法推导组成，对相应的地质信息数字化后，可以利用计算机编程语言开发相应的计算程序，实现井筒裂缝真三维空间渗透率张量定量计算。本发明对于真三维空间裂缝渗透率张量及井网设计、确定注水井与采油井的空间位置关系、寻找最优钻井轨迹等油田开发方案的部署具有较高的实用价值，并且预测成本低廉、可操作性强，预测结果可靠。

附图说明

图1为一种基于井筒裂缝观测的三维裂缝渗透率张量计算模型的流程图。

图2为一种基于井筒裂缝观测的三维裂缝渗透率张量数学模型示意图。

图3为苏北盆地金湖凹陷铜城断裂带天33断块构造位置图。

图4为天33断块阜二段岩心裂缝开度统计分布图。

图5为天33断块阜二段岩心裂缝走向玫瑰花图。

图6为天33断块水平剖面裂缝渗透率变化图。

图7为天33断块东西剖面裂缝渗透率变化图。

图8为天33断块南北剖面裂缝渗透率变化图。

图9为天33断块最大差渗透率面(252°∠46°)变化图。

图10为天33断块剖面等渗透率面(350°∠81°)变化图。

具体实施方式

下面结合附图说明本发明的具体实施方式：

以苏北盆地金湖凹陷铜城断裂带东翼阜宁组二段(简称阜二段)岩心裂缝三维渗透率张量定量计算为例来说明本发明的具体技术方案：

天33断块地处江苏金湖县境内，研究区主要位于苏北盆地-东台坳陷-金湖凹陷-铜城断裂带的东翼。铜城断层是一条典型的走滑断层，南部交汇于金湖凹陷的边界-杨村断层，向北消失于铜城地区(图3)。研究区主力产油层为阜二段，储层物性差，岩性以低渗透粉砂岩为主。裂缝对油气的渗流起到至关重要的作用。结合区域地质资料，对天33断块阜二段岩心裂缝渗透率张量定量计算的步骤如下：

第1步 对天33断块井筒裂缝发育程度的观测，利用岩心观测、测井资料解释对裂缝的开度测量(图4)，确定裂缝的组系特征；依据成像资料、测井或者岩心观察，求取裂缝的线密度，为后面计算裂缝平行渗透率准备。

第2步 采用斜井岩心裂缝产状校正方法对工区的裂缝进行定向(图5)。

第3步 通过确定岩石力学参数、井点地应力状态，建立地质力学模型，利用公式(1)计算裂缝在地下的开度。

第4步 如图2所示，建立真三维裂缝渗透率张量计算几何模型。

第5步 建立真三维裂缝渗透率张量数学模型，如公式(2)-公式(14)所示。

第6步 编制程序，输出渗透率剖面，利用公式(1)-(14)推导的算法编制相应的程序，通过循环迭代变量α(渗流面ψ的倾角)、β(渗流面ψ的倾向)、θ(OY轴与OP轴的夹角)，其中的α范围为[0-90°]，β的范围为[0-360°)，θ的范围为[0-90°)，根据计算精度需求，设置循环迭代步长⊿sp＝1°，将α、β、θ逐次累加⊿sp，利用公式(12)-(14)，计算坐标系(O-XYZ)不同的α、β、θ对应的K_XXαβθ(X方向渗透率主值)、K_YYαβθ(Y方向渗透率主值)、K_ZZαβθ(Z方向渗透率主值)。

利用公式(1)-(14)，设置α＝0°，β＝0°；调节θ的大小，步长⊿sp，范围为[0-90°)，输出K_YYαβθ(Y方向渗透率主值)、K_ZZαβθ(Z方向渗透率主值)，得到水平面的渗透率变化(图6)。

利用公式(1)-(14)，设置α＝90°，β＝0°；调节θ的大小，步长⊿sp，范围为[0-90°)，输出K_YYαβθ(Y方向渗透率主值)、K_ZZαβθ(Z方向渗透率主值)，得到东西向剖面的渗透率变化(图7)。

利用公式(1)-(14)，设置α＝90°，β＝90°；调节θ的大小，步长⊿sp，范围为[0-90°)，输出K_YYαβθ(Y方向渗透率主值)、K_ZZαβθ(Z方向渗透率主值)，得到南北向剖面的渗透率变化(图8)。

利用公式(1)-(14)，设置α范围为[0-90°]，β的范围为[0-360°)，计算不同的渗流面ψ内计算得到的K_YYαβθ、K_ZZαβθ数目为180，180个K_YYαβθ、K_ZZαβθ的平均值为K_aver；计算不同的α、β对应的方差函数W_αβ；通过公式(15)，计算不同的α、β对应的方差函数W_αβ，比较不同的α、β的W_αβ。其中，W_αβ最大值对应的渗流面ψ为最大差渗透率面(图9)；W_αβ最小值对应的渗流面ψ为等渗透率面(图10)。

上面以举例方式对本发明进行了说明，但本发明不限于上述具体实施例，凡基于本发明所做的任何改动或变型均属于本发明要求保护的范围。