CN104500050A - 一种裂缝性储层渗透率张量及各向异性定量预测方法 - Google Patents

Abstract

本发明涉及油气田勘探开发领域，尤其是一种裂缝性储层渗透率张量及各向异性定量预测方法。本发明在应力场数值模拟的基础上，利用裂缝的现今产状将静态坐标系与动态坐标系融合到大地坐标系中，建立了多组裂缝渗透率张量的定量预测模型，给出了渗透率主方向、大小的定量计算公式，通过调整动态坐标系旋转角预测单元体内裂缝在不同方向的渗透率。本发明由严格的数学算法推导组成，可以利用计算机编程语言开发相应的计算程序，实现裂缝性储层渗透率张量及各向异性定量预测，预测成本低廉、可操作性强，预测结果对确定断块渗透率优势方向、合理的部署开发井网、确定注水井与采油井的空间位置关系有一定的指导意义。

Description

一种裂缝性储层渗透率张量及各向异性定量预测方法

技术领域

本发明涉及油气田勘探开发领域，尤其是一种裂缝性储层渗透率张量及各向异性定量预测方法。

背景技术

在低渗透储层开发过程中，裂缝是油气渗流的主要通道，裂缝渗透率的各向异性和各向异性是影响油水流动方向的主要因素，裂缝性油气藏勘探、开发的最大难点在于储层岩体中裂缝发育程度、分布范围的预测以及裂缝渗透率各向异性分析。裂缝形成演化与古今地应力场演化密切相关，其中，有效应力对裂缝性低渗透油藏的渗流能力影响很大，裂缝渗透率随有效应力变化表现出较强的应力敏感性。因此，利用古今应力场数值模拟结果开展裂缝渗透率定量预测工作具有良好的理论基础，是一种比较可靠的裂缝定量预测方法。

在地下储层中，天然裂缝分布极为复杂，裂缝渗透率各向异性主要受裂缝产状、密度和开度等因素的影响。在利用裂缝性油气藏动静态相结合的建模方法中，预测结果大都局限于定性评价裂缝的渗透率，并没有考虑渗透率各向异性的特点；而传统的数值模拟方法，采用向量形式表征裂缝的渗透率，但前提是渗透率主方向与坐标轴方向相同，若渗透率主方向与坐标轴方向不一致或者有较大的差异，采用向量形式可能导致后期量化裂缝渗透率各向异性时出现较大的误差，难以符合油田合理高效开发的要求；刘月田等采用建立物理模型的方法评价裂缝性油藏的非均质性，具有较高的实用性和可行性。在本发明中，借用裂缝的现今产状，将静态坐标系与动态坐标系融合到大地坐标系中，同时利用古今应力场数值模拟结果，建立一套实用完善的裂缝渗透率各向异性定量预测-评价方法。

发明内容

本发明旨在解决上述问题，提供了一种裂缝性储层渗透率张量及各向异性定量预测方法，它解决了裂缝性储层渗透率张量及各向异性无法定量预测的问题。

本发明的技术方案为：一种裂缝性储层渗透率张量及各向异性定量预测方法，具体步骤如下：

第一步对研究区古应力场定时、定向、定量，进行古应力场数值模拟

结合区域构造演化、构造应力场演化，分析断层的活动性、活动强度，评价构造演化各个时期的强度，确定研究区的主要造缝期。

古应力场的定向主要依据岩心裂缝、薄片裂缝的观测结果以及成像测井资料，并结合区域活动断裂的产状，断裂的力学性质，地震资料中的褶皱、提取的地震属性以及薄片显微小构造确定主应力的方向。

古应力场大小确定可以采用声发射法、数学解析等方法，在应力场数值模拟中，可以采用一种等效古应力，即按照岩石破裂准则判断，现今应力场中，岩石不会产生明显的宏观破裂，同时古应力场作用下形成的裂缝还存在应力释放效应，因此现今地应力主要使已有的裂缝进一步延伸、扩张，一般不再形成新的裂缝，但可能会形成少量的微裂纹。等效古应力的确定主要是在确保应力类型与形成裂缝的应力类型一致的前提下，以单井储层裂缝参数为约束条件，结合地应力与储层裂缝参数之间的关系，通过不断地反演、拟合来实现。

在对古应力场定时、定向、定量后，依据恢复的古构造图以及岩石力学参数建立研究区古地质-力学模型，对研究区的古应力场数值模拟。

第二步计算裂缝的现今产状

运用应力场数值模拟方法能够得到古裂缝产状，而裂缝的形成发育主要受岩石力学层控制。所谓的岩石力学层是指一套岩石力学行为相似或者岩石力学性质基本一致的岩层，岩石力学层一般但不总是岩性均一层。因此，借助于古今岩石力学层的产状变化可以得到现今裂缝的产状：设裂缝形成时期所在平面的单位法向矢量为n′，倾角为α，倾向为β；现今裂缝面的单位法向矢量为n；古岩石力学层所在平面的法向矢量为p₁，现今岩石力学层所在平面的法向矢量为p₂，矢量p₁与矢量p₂所在平面的单位法向矢量为p。

矢量n′在大地坐标系中的三个分量为：

$\left\{\begin{array}{c}{n_x}^{\′}=\sin \mathrm{\α}\sin \mathrm{\β}\\ {n_y}^{\′}=\sin \mathrm{\α}\cos \mathrm{\β}\\ {n_z}^{\′}=\cos \mathrm{\α}\end{array}\right.---\left(1\right)$

公式(1)中，(n_x′，n_y′，n_z′)为矢量n′在大地坐标系中三个坐标轴的分量。

矢量n与矢量n′的定量转换关系可以表示为：

$\left[\begin{array}{c}{n}_x\\ n_y\\ n_z\end{array}\right]=T_1\·\left[\begin{array}{c}{n^{\′}}_x\\ {n^{\′}}_y\\ {n^{\′}}_z\end{array}\right]---\left(2\right)$

公式(2)中，(n_x，n_y，n_z)为矢量n在大地坐标系中三个坐标轴的分量。

公式(2)中，旋转矩阵T₁可以表示为：

$T_1=\left[\begin{array}{ccc}(1-{p_x}^2)\cos \mathrm{\ξ}+{p_x}^2& p_x{p}_y(1-\cos \mathrm{\ξ})+p_z\sin \mathrm{\ξ}& p_x{p}_z(1-\cos \mathrm{\ξ})-p_y\sin \mathrm{\ξ}\\ p_x{p}_y(1-\cos \mathrm{\ξ})-p_z\sin \mathrm{\ξ}& (1-{p_y}^2)\cos \mathrm{\ξ}+{p_y}^2& p_y{p}_z(1-\cos \mathrm{\ξ})+p_x\sin \mathrm{\ξ}\\ p_x{p}_z(1-\cos \mathrm{\ξ})+p_y\sin \mathrm{\ξ}& p_y{p}_z(1-\cos \mathrm{\ξ})-p_x\sin \mathrm{\ξ}& (1-{p_z}^2)\cos \mathrm{\ξ}+{p_z}^2\end{array}\right]---\left(3\right)$

公式(3)中，(p_x，p_y，p_z)为矢量p在大地坐标系中三个坐标轴的分量；ξ为矢量p₁与矢量p₂重合的旋转角，(°)。

第三步现今应力场数值模拟

现今地应力场对于油气田的开发来说具有重要意义。研究现今地应力场，可以为油气田的开发方案等提供理论支持，尤其是可以对开发过程中地层物性的变化进行宏观把握，进而制定相应的储层保护措施。新构造期的构造变形相对稳定，构造应力相对较低，在我国大陆板块内部没有发生较为强烈的构造变形，因此对于储层裂缝研究而言，现今地应力的作用主要是对古应力场下所形成的裂缝进行改造，裂缝的开度和延伸长度会随着现今地应力的变化而变化，并进一步引起裂缝孔隙度和裂缝渗透率的变化，但密度和组系特征基本保持不变。水力压裂法、波速法、声发射法、应力解除法、光弹性应力法等。其中水力压裂对环境的要求比较宽松，测量应力的代表地下空间范围较大，受局部因素的影响较小，成功率较高，是一种较为方便直接的进行地应力测量的方法。在钻井过程中，会在井壁附近产生一定数量的钻井诱导缝。钻井诱导缝往往指示了现今主应力的方向，其走向一般平行于现今最大主应力。另外，随着井筒岩心的不断取出，井壁应力发生释放，在井壁的应力集中或者脆弱处往往产生井壁坍塌，井壁坍塌的方位一般平行于现今最小主应力的方向。因此，现今主应力的方向也可以通过钻井诱导缝的走向和井壁坍塌方位来确定。

通常采用多种方法对研究区现今地应力方位进行了判断，并通过测井、压裂资料以及物理实验计算了关键井现今地应力的数值。在此基础上，通过确定岩石力学参数并建立有限元模型，对研究区现今地应力场进行了数值模拟。

第四步现今裂缝开度、平行渗透率预测

利用构造应力场数值模拟方法进行储层构造裂缝预测具有良好的理论基础，是一种比较可靠的裂缝预测方法，国内学者戴俊生、汪必峰、季宗镇和冯建伟等从影响储层构造裂缝形成与发育的地质因素入手，通过对裂缝形成时期的古应力场数值模拟，结合岩石破裂准则和单井裂缝描述成果，利用岩石应变能理论、裂缝表面能理论以及能量守恒定律，建立了古今应力场与储层构造裂缝参数之间的定量关系，实现了裂缝线密度预测。在现今应力场中一般不再产生新的裂缝，裂缝的线密度、产状及组系特征基本保持不变，但在三向挤压应力影响下，已有的裂缝会出现一定程度的闭合，开度降低，渗流性能变差。Willis-Richards等和Jing等同时考虑了正应力及剪应力对裂缝开度的影响，得出在现今地应力场改造后裂缝开度的计算公式为：

$b=\frac{b_0}{1+9{\mathrm{\σ}}_n^{\′}/{\mathrm{\σ}}_{\mathrm{nref}}}+b_{\mathrm{res}}---\left(4\right)$

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

在古今应力场数值模拟的基础上，可以求取单元体内现今裂缝线密度以及开度，因此，单组裂缝的平行渗透率K表示为：

$K=\frac{b^3{D}_{\mathrm{lf}}}{12}---\left(5\right)$

公式(5)中，b为裂缝的现今开度，m；D_lf为单组裂缝的线密度，条/m。

第五步建立裂缝渗透率张量计算模型

地下岩层受地应力作用破裂产生裂缝后，裂缝渗透率张量可以借助于适当的模型进行研究。通过建立微小单元体模型用以表征裂缝渗透率的各向异性，为了满足研究需要，认为单元体模型足够小，所有裂缝均能将其切穿。

如图2所示，以裂缝为参照物建立静态坐标系(O-EFG)，静态坐标系中的三个坐标轴分别对应于裂缝面的法线方向(OE轴)、裂缝走向方向(OF轴)、裂缝面内垂直于裂缝走向的方向(OG轴)。以大地坐标为参考建立动态坐标系(O-XYZ)，定义θ为水平面内OX轴与正东方向的夹角，即θ是动态坐标系的旋转角，通过调整θ的大小，求取裂缝在不同坐标系中的渗透率张量；定义OX轴位于北东方向时θ为负值，南东方向时θ为正值；动态坐标系中OZ轴为铅直方向。

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

$K=(X,Y,Z)\·K_{\mathrm{XYZ}}\·\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]---\left(6\right)$

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

$K=(E,F,G)\·K_{\mathrm{EFG}}\·\left[\begin{array}{c}E\\ F\\ G\end{array}\right]---\left(7\right)$

在静态坐标系中，每个单元体内单组裂缝的渗透率张量可以表示为：

$K_{\mathrm{EFG}}=\left[\begin{array}{ccc}0& 0& 0\\ 0& K& 0\\ 0& 0& K\end{array}\right]---\left(8\right)$

由公式(6)-(8)可以得到：

K_XYZ＝T₂·K_EFG·T₂ ^T      (9)

其中，静态坐标系转换为动态坐标系的旋转矩阵T₂表示为：

$T_2=\left[\begin{array}{ccc}{C}_{\mathrm{\θ}}& \frac{D_{\mathrm{\θ}}}{\sqrt{{n_x}^2+{n_y}^2}}& \frac{-n_z\·C_{\mathrm{\θ}}}{\sqrt{{n_x}^2+{n_y}^2}}\\ D_{\mathrm{\θ}}& \frac{-C_{\mathrm{\θ}}}{\sqrt{{n_x}^2+{n_y}^2}}& \frac{-n_z\·D_{\mathrm{\θ}}}{\sqrt{{n_x}^2+{n_y}^2}}\\ n_z& 0& \sqrt{{n_x}^2+{n_y}^2}\end{array}\right]---\left(10\right)$

在动态坐标系中，每个单元体内单组裂缝的渗透率张量可以表示为：

$K_{\mathrm{XYZ}}=K\left[\begin{array}{ccc}\frac{{C_{\mathrm{\θ}}}^2+{n_z}^2{C_{\mathrm{\θ}}}^2}{\sqrt{{n_x}^2+{n_y}^2}}& -C_{\mathrm{\θ}}{D}_{\mathrm{\θ}}& -n_z{C}_{\mathrm{\θ}}\\ -C_{\mathrm{\θ}}{D}_{\mathrm{\θ}}& \frac{{C_{\mathrm{\θ}}}^2+{n_z}^2{D_{\mathrm{\θ}}}^2}{\sqrt{{n_x}^2+{n_y}^2}}& -n_z{D}_{\mathrm{\θ}}\\ -n_z{C}_{\mathrm{\θ}}& -n_z{D}_{\mathrm{\θ}}& {n_x}^2+{n_y}^2\end{array}\right]---\left(11\right)$

公式(10)-(11)中，C_θ、D_θ满足：

$\left\{\begin{array}{c}{C}_{\mathrm{\θ}}=n_x\cos \mathrm{\θ}-n_y\sin \mathrm{\θ}\\ D_{\mathrm{\θ}}=n_x\sin \mathrm{\θ}+n_y\cos \mathrm{\θ}\end{array}\right.---\left(12\right)$

第六步多组裂缝渗透率主方向、主值预测

在复杂的地质条件下，储层中往往发育多组裂缝，在古今应力场综合作用下，各组裂缝的产状、线密度、开度往往不同，依据公式(5)、公式(11)以及公式(12)，每个单元体内多组裂缝最大渗透率的主方向可以表示为：

${\mathrm{\θ}}_{\max }=\arctan \left(\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{b_i}^3{D}_{\mathrm{lfi}}{n}_{\mathrm{xi}}}{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{b_i}^3{D}_{\mathrm{lfi}}{n}_{\mathrm{yi}}}\right)---\left(13\right)$

每个单元体内裂缝最大渗透率主方向上，渗透率大小可以表示为：

$K_{\max }=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}(\frac{{b_i}^3{D}_{\mathrm{lfi}}}{12}\·\frac{{(n_{\mathrm{xi}}\sin {\mathrm{\θ}}_{\max }+n_{\mathrm{yi}}\cos {\mathrm{\θ}}_{\max })}^2+{n_{\mathrm{zi}}}^2{(n_{\mathrm{xi}}\cos {\mathrm{\θ}}_{\max }-n_{\mathrm{yi}}\sin {\mathrm{\θ}}_{\max })}^2}{\sqrt{{n_{\mathrm{xi}}}^2+{n_{\mathrm{yi}}}^2}})---\left(14\right)$

公式(13)-(14)中，k为裂缝的组数；b_i为第i组裂缝的开度，m；D_lfi为第i组裂缝的线密度，条/m；(n_xi，n_yi，n_zi)分别为第i组裂缝面的单位法向向量在大地坐标系中三个坐标轴的分量，可由公式(1)-(3)得到。

第七步渗透率各向异性定量评价

裂缝的渗透率具有强烈的各向异性，参考表征储层各向异性的参数，定义了评价裂缝各向异性的三个参数-渗透率极差比τ_b、渗透率突进系数τ_k、渗透率变异系数τ_y，由于地质因素的影响具有地区性和随机性，因此针对不同地区裂缝的发育规律，选取具有代表性的非均质性参数，当单元体内不同方向(0°-180°)的渗透率变化为椭圆型时(单峰单谷型)，选用渗透率极差比τb便能较好的评价渗透率非均质性；当单元体内不同方向(0°-180°)的渗透率多次起伏变化时(多峰多谷型)，宜采用突进系数τ_k或变异系数τ_y评价渗透率的非均质性。

渗透率极差比(τ_b)表征的是井点附近渗透率最大值与最小值的比值，渗透率极差比的变化范围为(1-∞)，数值越小越均质，越大越非均质。

${\mathrm{\τ}}_b=\frac{K_{\max }}{K_{\min }}---\left(15\right)$

公式(15)中，K_max为单元体内渗透率最大值，10^-3μm²；K_min为单元体内渗透率最小值，10^-3μm²。

渗透率突进系数(τ_k)表征的是一定井段内渗透率最大值与其平均值的比值，同样，τ_k值变化范围为(1-∞)，数值越小越均质，越大越非均质。

${\mathrm{\τ}}_k=\frac{K_{\max }}{\stackrel{\‾}{K}}---\left(16\right)$

公式(16)中，K_max为单元体内渗透率最大值，10^-3μm²；K为单元体内不同方向的渗透率平均值，10^-3μm²。

渗透率变异系数(τ_y)是一个数理统计的概念，用于度量统计的若干数值相对于平均值的离散程度。其反映的是单元体内裂缝渗透率偏离整体平均值的程度。数值变化范围(0-∞),数值越小越均质，越大越非均质。

${\mathrm{\τ}}_y=\frac{\sqrt{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(K_{\mathrm{\θi}}-\stackrel{\‾}{K})}^2/(n-1)}}{\stackrel{\‾}{K}}---\left(17\right)$

公式(17)中，K_θi为单元体内θi方向的渗透率值，10^-3μm²；n为计算的不同方向裂缝渗透率的数目；K为单元体内不同方向的渗透率平均值，10^-3μm²。

本发明的有益效果是：本发明在应力场数值模拟的基础上，借助于古今岩石力学层，预测裂缝的现今产状；以断裂力学中裂缝表面能以及岩石应变能理论预测现今裂缝的线密度；在现今应力场数值模拟的基础上，计算三向挤压条件下裂缝的开度，进而确定现今裂缝的平行渗透率。利用裂缝的现今产状将静态坐标系与动态坐标系融合到大地坐标系中，建立了多组裂缝渗透率张量的定量预测模型，给出了渗透率主方向、大小的定量计算公式，并且通过调整动态坐标系旋转角预测单元体内裂缝在不同方向的渗透率。定义表征渗透率各向异性的三个参数-裂缝渗透率极差比、渗透率突进系数、渗透率变异系数，定量评价裂缝渗透率各向异性。本发明由严格的数学算法推导组成，对相应的地质信息数字化后，可以利用计算机编程语言开发相应的计算程序，实现裂缝性储层渗透率张量及各向异性定量预测。本发明对于裂缝性储层渗透率张量及各向异性定量预测-评价具有较高的实用价值，并且预测成本低廉、可操作性强，预测结果对确定断块渗透率优势方向、合理的部署开发井网、确定注水井与采油井的空间位置关系有一定的指导意义。

附图说明

图1为一种裂缝性储层渗透率张量及各向异性定量预测方法的流程图。

图2为静态坐标系与动态坐标系关系示意图。

图3为铜城断裂带东翼构造位置图。

图4为铜城断裂带东翼阜宁晚期最小主应力分布图。

图5为铜城断裂带东翼现今水平最小主应力分布图。

图6为铜城断裂带东翼现今水平最大主应力分布图。

图7为铜城断裂带东翼裂缝渗透率主值大小分布图。

图8为铜城断裂带东翼裂缝渗透率主值方向分布图。

图9为铜城断裂带东翼不同井点的不同方向上裂缝渗透率变化图。

图10为铜城断裂带东翼裂缝渗透率极差比分布图。

在图2中，θ为动态坐标系的旋转角，(°)。

在图4中，负值代表挤压应力，正值代表拉张应力(单位：Pa)。

在图5中，负值代表挤压应力，正值代表拉张应力(单位：Pa)。

在图6中，负值代表挤压应力，正值代表拉张应力(单位：Pa)。

在图7中，数值单位：10^-3μm²。

在图8中，数值单位：(°)。

在图10中，数值为无量纲单位。

具体实施方式

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

以苏北盆地金湖凹陷铜城断裂带东翼阜宁组二段(简称阜二段)裂缝性储层渗透率张量及各向异性定量预测为例来说明本发明的具体技术方案：

铜城断裂带地处江苏金湖县境内，研究区主要位于苏北盆地-东台坳陷-金湖凹陷-铜城断裂带的东翼。铜城断层是一条典型的走滑断层，南部交汇于金湖凹陷的边界-杨村断层，向北消失于铜城地区(图3)。研究区主力产油层为阜二段，储层物性差，岩性以低渗透粉砂岩为主。

铜城断裂带自阜二段沉积以来，主要经历吴堡和三垛两期构造运动。研究区阜二段地层主要造缝期为吴堡期，该时期断层强烈活动，主要表现为近南北向拉张，形成了凹陷内部复杂的断裂系统。结合区域地质资料，对铜城断裂带东翼现今构造裂缝产状定量预测的步骤如下：

第1步对铜城断裂带东翼裂缝形成时期的古应力场定时、定向、定量，进行古应力场数值模拟。

对铜城断裂带东翼进行古应力场数值模拟主要依托ANSYS软件。根据铜城断裂带构造演化，分析断层的活动性，评价构造演化各个时期的强度，确定研究区的主要造缝期。由于吴堡运动，已形成的生长断层进一步活动。自阜一期至阜四期断层活动速率逐渐增强，在阜四期达到最高值，即阜四期断层活动性最强；戴南期、三垛期断层活动性逐渐减弱，因此，基本确定研究区的主要造缝期为阜宁晚期。

对铜城断裂带阜宁晚期古应力场的定向主要依据岩心、薄片裂缝的观测结果并结合造缝期活动断裂的产状、断裂的力学性质，确定主应力的方向。根据铜城断裂带阜宁晚期活动断裂走向玫瑰花图以及对铜城断裂带岩心裂缝定向得到的铜城断裂带东翼阜二段裂缝走向玫瑰花图，确定阜宁晚期最小主应力为近南北向的拉张应力，最大主应力为垂直方向，主要由重力应力产生。

古应力场大小确定采用一种等效古应力，即按照岩石破裂准则判断，现今应力场中，岩石不会产生明显的宏观破裂，同时古应力场作用下形成的裂缝还存在应力释放效应，因此现今地应力主要使已有的裂缝进一步延伸、扩张，一般不再形成新的裂缝，但可能会形成少量的微裂纹。等效古应力的确定主要是在确保应力类型与形成裂缝的应力类型一致的前提下，以单井储层裂缝参数为约束条件，结合地应力与储层裂缝参数之间的关系，通过不断地反演、拟合来实现，最终确定在铜城断裂带东翼模型的南北边界施加4.5MPa的拉张应力，东西边界施加1.5MPa拉张应力，垂向应力通过设置重力加速度由ANSYS软件自动生成。

对古应力场定时、定向、定量后，依据恢复的古构造图以及岩石力学参数建立研究区古地质-力学模型，对研究区的古应力场数值模拟(图4，负值代表挤压应力，正值代表拉张应力)。

第2步计算研究区裂缝的现今产状。

通过古应力场数值模拟，结合岩心、薄片观测结果，选择合适的破裂准则，依据岩石力学实验结果求取相应的岩石力学参数，可以得到古裂缝的产状，采用公式(1)-(3)得到裂缝的现今产状。

第3步现今应力场数值模拟。

通过岩心声速实验、古地磁实验以及微地震监测并结合区域构造背景对铜城地区现今地应力方位进行了判断，并利用压裂资料计算了关键井现今地应力数值，确定阜二段现今最小水平主应力为NNW349°方向25MPa的挤压应力，水平最大主应力为NEE79°方向37MPa的挤压应力，通过确定岩石力学参数并建立有限元模型，对铜城地区现今地应力场进行了数值模拟，模拟结果显示：在现今应力场中，水平最小主应力在24MPa～32MPa之间，水平最大主应力在35MPa～67MPa之间；最小主应力梯度为11MPa/km～17MPa/km，最大主应力的梯度为23MPa/km～28MPa/km；水平主应力方向在断块内变化不大，断裂带附近应力方向发生3°～10°的变化(图5-图6)。

第4步现今裂缝开度、平行渗透率预测。

古应力场模拟结果显示，铜城断裂带东翼阜二段储层均已达到岩石破裂条件，断层附近张应力最大，裂缝相对发育，越靠近断层，裂缝发育的概率越大。将应力、应变数据导入裂缝参数计算模型，预测了该地区阜二段地层现今裂缝参数。模拟结果显示，铜城断裂带东翼裂缝线密度为2.85-4.68条/m，断层的附近以及构造低部位裂缝的线密度为高值，主要方位与次要方位的裂缝线密度比值为3:2。由于现今应力场的三向挤压作用，裂缝古今开度差别很大。古应力场中裂缝开度为1.07×10^-4-1.33×10^-4m。现今应力场中裂缝开度为2.0×10^-5-4.4×10^-5m，现今裂缝开度约是古开度的0.2-0.3倍。由于应力场的模拟结果是裂缝参数预测的基础，其中现今应力场又是裂缝现今开度预测关键，因此对比铜城断裂带东翼部分井位实测结果与现今应力场模拟结果，发现两者吻合较好(表1)，由此，认为裂缝参数的预测结果可以较准确地反映了地下裂缝的实际情况。

表1 铜城断裂带东翼实测结果与现今应力场数值模拟对比

第5步建立裂缝渗透率张量计算模型。

利用现今裂缝产状建立静态坐标系，利用公式(6)-(11)建立动态坐标系以静态坐标系的关系，采用裂缝最大渗透率主值、最大渗透率主方向、渗透率极差比三个参数定量评价铜城断裂带东翼阜二段裂缝水平面内渗透率的各向异性。

第6步多组裂缝最大渗透率主方向、主值预测。

利用公式(13)-(14)，求取铜城断裂带东翼裂缝的多组裂缝最大渗透率主方向、主值预测，铜城断裂带东翼阜二段储层最大渗透率主值在断层附近为高值，在构造低部位为低值(图7)；天96断块最大渗透率主方向以北东东向为主，天33断块最大渗透率主方向以南东东向为主(图8)，这与油田的动态开发资料吻合较好。

通过不断地调整动态坐标系中θ值的大小，能够得到井点附近不同方向上的裂缝渗透率数值(图9)。由图得出，在靠近断层的井点，由于裂缝开度大，裂缝的渗透率主值大；井点附近裂缝的渗透率优势方向为近东西向，这也与油田的动态开发资料吻合较好。

第7步定量评价铜城断裂带东翼渗透率各向异性。

如图10所示，利用公式(15)-(17)，在构造高部位，裂缝的渗透率极差比为低值；在构造低部位，渗透率极差比为高值。因此得出，在对构造低部位油气勘探时，同样应该着重考虑裂缝所造成的渗透率非均质性。

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

Claims (4)

1.一种裂缝性储层渗透率张量及各向异性定量预测方法，所述的步骤如下：

1)对研究区裂缝形成时期的古应力场定时、定向、定量，确定研究区岩石力学参数，建立有限元地质模型进行古应力场数值模拟；

2)运用应力场数值模拟方法得到古裂缝产状，利用古裂缝产状计算现今裂缝产状：

设裂缝形成时期所在平面的单位法向矢量为n′，倾角为α，倾向为β；现今裂缝面的单位法向矢量为n；古岩石力学层所在平面的法向矢量为p₁，现今岩石力学层所在平面的法向矢量为p₂，矢量p₁与矢量p₂所在平面的单位法向矢量为p，矢量n′在大地坐标系中的三个分量为：

$\left\{\begin{array}{c}{n_x}^{\′}=\sin \mathrm{\α}\sin \mathrm{\β}\\ {n_y}^{\′}=\sin \mathrm{\α}\cos \mathrm{\β}\\ {n_z}^{\′}=\cos \mathrm{\α}\end{array}\right.---\left(1\right)$

公式(1)中，(n_x′，n_y′，n_z′)为矢量n′在大地坐标系中三个坐标轴的分量；α为裂缝形成时期的倾角，(°)；β为裂缝形成时期的倾向，(°)；

矢量n与矢量n′的定量转换关系可以表示为：

$\left[\begin{array}{c}{n}_x\\ n_y\\ n_z\end{array}\right]=T_1\·\left[\begin{array}{c}{n^{\′}}_x\\ {n^{\′}}_y\\ {n^{\′}}_z\end{array}\right]---\left(2\right)$

公式(2)中，(n_x，n_y，n_z)为矢量n在大地坐标系中三个坐标轴的分量；

公式(2)中，旋转矩阵T₁可以表示为：

$T_1=\left[\begin{array}{ccc}(1-{p_x}^2)\cos \mathrm{\ξ}+{p_x}^2& p_x{p}_y(1-\cos \mathrm{\ξ})+p_z\sin \mathrm{\ξ}& p_x{p}_z(1-\cos \mathrm{\ξ})-p_y\sin \mathrm{\ξ}\\ p_x{p}_y(1-\cos \mathrm{\ξ})-p_z\sin \mathrm{\ξ}& (1-{p_y}^2)\cos \mathrm{\ξ}+{p_y}^2& p_y{p}_z(1-\cos \mathrm{\ξ})+p_x\sin \mathrm{\ξ}\\ p_x{p}_z(1-\cos \mathrm{\ξ})+p_y\sin \mathrm{\ξ}& p_y{p}_z(1-\cos \mathrm{\ξ})-p_x\sin \mathrm{\ξ}& (1-{p_z}^2)\cos \mathrm{\ξ}+{p_z}^2\end{array}\right]---\left(3\right)$

公式(3)中，(p_x，p_y，p_z)为矢量p在大地坐标系中三个坐标轴的分量；ξ为矢量p₁与矢量p₂重合的旋转角，(°)；

3)通常采用多种方法对研究区现今地应力方位进行了判断，并通过测井、压裂资料以及物理实验计算关键井现今地应力的数值，在此基础上，通过确定岩石力学参数并建立有限元模型，对研究区现今地应力场进行了数值模拟；

4)利用岩石应变能理论、裂缝表面能理论以及能量守恒定律，建立了古今应力场与储层构造裂缝参数之间的定量关系，实现了裂缝线密度预测；在现今应力场中改造后裂缝开度的计算公式为：

$b=\frac{b_0}{1+9{\mathrm{\σ}}_n^{\′}/{\mathrm{\σ}}_{\mathrm{nref}}}+b_{\mathrm{res}}---\left(4\right)$

公式(4)中，b为裂缝的现今开度，m；b₀为裂缝的原始开度，m；σ′_n为有效正应力，MPa；b_res代表裂缝面承受最大有效正应力时的裂缝开度，m；σ_nref为使裂缝开度减小90％的有效正应力，MPa；

在古今应力场数值模拟的基础上，可以求取单元体内现今裂缝线密度以及开度，单组裂缝的平行渗透率K表示为：

$K=\frac{b^3{D}_{\mathrm{lf}}}{12}---\left(15\right)$

公式(5)中，b为裂缝的现今开度，m；D_lf为单组裂缝的线密度，条/m；

5)利用裂缝的产状将静态坐标系与动态坐标系统一到大地坐标系中，建立了裂缝渗透率张量的定量预测模型；

6)利用建立缝渗透率张量的定量预测模型，建立多组裂缝渗透率主方向、主值预测计算模型；

7)建立渗透率各向异性定量预测模型，采用的渗透率极差比τ_b、渗透率突进系数τ_k、渗透率变异系数τ_y评价裂缝各向异性。

2.根据权利要求1所述的一种裂缝性储层渗透率张量及各向异性定量预测方法，其特征在于：

所述的裂缝渗透率张量的定量预测模型，原理如下：

以裂缝为参照物建立静态坐标系(O-EFG)，静态坐标系中的三个坐标轴分别对应于裂缝面的法线方向(OE轴)、裂缝走向方向(OF轴)、裂缝面内垂直于裂缝走向的方向(OG轴)；以大地坐标为参考建立动态坐标系(O-XYZ)，定义θ为水平面内OX轴与正东方向的夹角，即θ是动态坐标系的旋转角，通过调整θ的大小，求取裂缝在不同坐标系中的渗透率张量；定义OX轴位于北东方向时θ为负值，南东方向时θ为正值；动态坐标系中OZ轴为铅直方向；每个单元体内单组裂缝的渗透率张量可以表示为：

$K_{\mathrm{XYZ}}=K\left[\begin{array}{ccc}\frac{{D_{\mathrm{\θ}}}^2+{n_z}^2{C_{\mathrm{\θ}}}^2}{\sqrt{{n_x}^2+{n_y}^2}}& -C_{\mathrm{\θ}}{D}_{\mathrm{\θ}}& -n_z{C}_{\mathrm{\θ}}\\ -C_{\mathrm{\θ}}{D}_{\mathrm{\θ}}& \frac{{C_{\mathrm{\θ}}}^2+{n_z}^2{D_{\mathrm{\θ}}}^2}{\sqrt{{n_x}^2+{n_y}^2}}& -n_z{D}_{\mathrm{\θ}}\\ -n_z{C}_{\mathrm{\θ}}& -n_z{D}_{\mathrm{\θ}}& {n_x}^2+{n_y}^2\end{array}\right]---\left(11\right)$

公式(11)中，C_θ、D_θ满足：

$\left\{\begin{array}{c}{C}_{\mathrm{\θ}}=n_x\cos \mathrm{\θ}-n_y\sin \mathrm{\θ}\\ D_{\mathrm{\θ}}=n_x\sin \mathrm{\θ}+n_y\cos \mathrm{\θ}\end{array}\right.---\left(12\right)$

公式(10)-(11)中，(n_x，n_y，n_z)为矢量n在大地坐标系中三个坐标轴的分量；θ为动态坐标系的旋转角，(°)。

3.根据权利要求1所述的一种裂缝性储层渗透率张量及各向异性定量预测方法，其特征在于：

所述的多组裂缝渗透率主方向、主值预测计算模型中具体算法为：

每个单元体内多组裂缝最大渗透率的主方向可以表示为：

${\mathrm{\θ}}_{\max }=\arctan \left(\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{b_i}^3{D}_{\mathrm{lfi}}{n}_{\mathrm{xi}}}{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{b_i}^3{D}_{\mathrm{lif}}{n}_{\mathrm{yi}}}\right)---\left(13\right)$

每个单元体内裂缝最大渗透率主方向上，渗透率大小可以表示为：

$K_{\max }=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}(\frac{{b_i}^3{D}_{\mathrm{lfi}}}{12}\·\frac{{(n_{\mathrm{xi}}\sin {\mathrm{\θ}}_{\max }+n_{\mathrm{yi}}\cos {\mathrm{\θ}}_{\max })}^2+{n_{\mathrm{zi}}}^2{(n_{\mathrm{xi}}\cos {\mathrm{\θ}}_{\max }-n_{\mathrm{yi}}\sin {\mathrm{\θ}}_{\max })}^2}{\sqrt{{n_{\mathrm{xi}}}^2+{n_{\mathrm{yi}}}^2}})---\left(14\right)$

公式(13)-(14)中，k为裂缝的组数；b_i为第i组裂缝的开度，m；D_lfi为第i组裂缝的线密度，条/m；(n_xi，n_yi，n_zi)分别为第i组裂缝面的单位法向向量在大地坐标系中三个坐标轴的分量，可由公式(1)-(3)得到。

4.根据权利要求1所述的一种裂缝性储层渗透率张量及各向异性定量预测方法，其特征在于：

所述的评价裂缝各向异性的三个参数-渗透率极差比τ_b、渗透率突进系数τ_k、渗透率变异系数τ_y具体算法为：

渗透率极差比(τ_b)表征的是井点附近渗透率最大值最小值与的比值，渗透率极差比的变化范围为(1-∞)，数值越小越均质，越大越非均质：

${\mathrm{\τ}}_b=\frac{K_{\max }}{K_{\min }}---\left(15\right)$

公式(15)中，K_max为单元体内渗透率最大值，10^-3μm²；Kmin为单元体内渗透率最小值，10^-3μm²；

渗透率突进系数(τ_k)表征的是一定井段内渗透率最大值与其平均值的比值，同样，τ_k值变化范围为(1-∞)，数值越小越均质，越大越非均质：

${\mathrm{\τ}}_k=\frac{K_{\max }}{\stackrel{\‾}{K}}---\left(16\right)$

公式(16)中，K_max为单元体内渗透率最大值，10^-3μm²；K为单元体内不同方向的渗透率平均值，10^-3μm²；

渗透率变异系数(τ_y)是一个数理统计的概念，用于度量统计的若干数值相对于平均值的离散程度。其反映的是单元体内裂缝渗透率偏离整体平均值的程度；数值变化范围(0-∞),数值越小越均质，越大越非均质：

${\mathrm{\τ}}_y=\frac{\sqrt{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(K_{\mathrm{\θi}}-\stackrel{\‾}{K})}^2/(n-1)}}{\stackrel{\‾}{K}}---\left(17\right)$

公式(17)中，K_θi为单元体内θi方向的渗透率值，10^-3μm²；n为计算的不同方向裂缝渗透率的数目；K为单元体内不同方向的渗透率平均值，10^-3μm²；

针对不同地区裂缝的发育规律，选取具有代表性的非均质性参数，当单元体内不同方向(0°-180°)的渗透率变化为椭圆型时(单峰单谷型)，选用渗透率极差比τ_b便能较好的评价渗透率非均质性；当单元体内不同方向(0°-180°)的渗透率多次起伏变化时(多峰多谷型)，宜采用突进系数τ_k或变异系数τ_y评价渗透率的非均质性。