CN103760602A - 一种三维水力压裂及微地震事件正演模拟的方法及系统 - Google Patents

Abstract

本发明提供一种三维水力压裂及微地震事件正演模拟的方法及系统，方法包括：采集当前地层的地层数据；根据当前地层裂缝的初始破裂位置、当前地层的区域范围设置网格；根据基于压裂缝的流量守恒关系构建基于含压裂缝网格的流量守恒关系；根据基于含压裂缝网格的流量守恒关系确定三维孔隙压力分布；从地层数据中获取网格所在位置的临界压力值；根据三维孔隙压力分布与网格所在位置的临界压力值生成当前地层的微地震分布；根据微地震分布开发当前地层。拓展了常规水力压裂模拟技术，在获得压裂缝信息之外，还能得到整个地层的孔隙压力分布；在水力压裂模拟基础上，融入孔隙压力控制型微地震的正演模拟，更好的表征压裂过程中微地震的变化特征。

Description

一种三维水力压裂及微地震事件正演模拟的方法及系统

技术领域

本发明关于地球物理勘探技术，特别是关于水力压裂及微地震技术，具体的讲是一种三维水力压裂及微地震事件正演模拟的方法及系统。

背景技术

在非常规油气开采过程中，由于地层的渗透率过低，常规油气开采技术无法获得客观产量，往往采用水力压裂的方式来提高采收率。水力压裂技术利用流体传压特性，将流体高压注入地层，使岩石破裂并延伸一定长度，改善地层渗透率以达到提高油气产量的目的。在水力压裂过程中，微地震监测技术是通过观测、分析微地震信息来确定裂缝走向、破裂类型和地下信息，为压裂过程监测、分析的核心技术。因此，构建水力压裂与微地震综合模拟技术对非常规油气勘探开发来说是一项必不可少的关键技术。

目前水力压裂模拟中，往往忽略压裂缝内流体漏失或仅通过漏失率近似压裂缝向基质的渗流过程，模拟结果只能得到裂缝形态、缝内孔隙压力信息，无法得到压裂过程中周边地层的裂隙、孔隙压力变化。

在压裂过程中，压裂液从压裂缝中渗流到周围地层中，改变地层的孔隙压力，降低有效应力，会引起岩石破裂和改变地层弹性性质，而传统的水力压裂模拟方法无法获得这些变化信息。微地震模拟研究多集中于常规油气注采、地热开发和应力挤压等过程，往往通过PKN等公式近似压裂过程中的微地震信号变化。公式拟合方法只能够给出整体分布趋势，无法得到精确和微观的微地震变化过程。

因此，在水力压裂过程中，如何综合模拟压裂缝破裂以及地层孔隙压力变化对微地震事件的分布进行模拟成为水力压裂及微地震研究中急需解决的问题。

发明内容

为了克服现有技术存在的上述问题，本发明提供了一种三维水力压裂及微地震事件正演模拟的方法及系统，拓展了常规水力压裂模拟技术，在获得压裂缝信息之外，还能够得到整个地层的孔隙压力分布；同时在水力压裂模拟基础上，融入孔隙压力控制型微地震的正演模拟，能够更好的表征压裂过程中微地震的变化特征。

本发明的目的之一是，提供一种三维水力压裂及微地震事件正演模拟的方法，包括：采集当前地层的地层数据；根据当前地层裂缝的初始破裂位置、当前地层的区域范围设置网格；根据基于压裂缝的流量守恒关系构建基于含压裂缝网格的流量守恒关系；根据基于含压裂缝网格的流量守恒关系确定当前地层的三维孔隙压力分布；从所述的地层数据中获取网格所在位置的临界压力值；根据所述三维孔隙压力分布与网格所在位置的临界压力值生成当前地层的微地震分布；根据所述的微地震分布开发当前地层。

本发明的目的之一是，提供了一种三维水力压裂及微地震事件正演模拟的系统，包括：地层数据采集装置，用于采集当前地层的地层数据；网格设置装置，用于根据当前地层裂缝的初始破裂位置、当前地层的区域范围设置网格；流量守恒构建装置，用于根据基于压裂缝的流量守恒关系构建基于含压裂缝网格的流量守恒关系；孔隙压力分布确定装置，用于根据基于含压裂缝网格的流量守恒关系确定当前地层的三维孔隙压力分布；临界压力值获取装置，用于从所述的地层数据中获取网格所在位置的临界压力值；微地震分布确定装置，用于根据所述三维孔隙压力分布与网格所在位置的临界压力值生成当前地层的微地震分布；地层开发装置，用于根据所述的微地震分布开发当前地层。

本发明的有益效果在于，提供了一种三维水力压裂及微地震事件正演模拟的方法及系统，将常规基于压裂缝的流量守恒关系推广至含压裂缝网格，以保证模拟过程中压裂缝与基质网格在尺度上的匹配。同时综合流量守恒、渗流方程，同步模拟压裂缝形态、缝内孔隙压力和基质孔隙压力变化。在孔隙压力模拟的基础上，采用断裂力学模拟压裂缝的破裂，同时基于临界压力理论模拟孔隙压力控制型的微地震事件。该方案与常规水力压裂方法和微地震相比，能够同步模拟压裂缝破裂、缝内孔隙压力变化、基质孔隙压力变化和孔隙压力控制型的微地震事件，获得整个地层的微地震事件和流体压力分布，为进一步的微地震研究和地层参数研究提供基础。

为让本发明的上述和其他目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附图式，作详细说明如下。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的一种三维水力压裂及微地震事件正演模拟的方法的实施方式一的流程图；

图2为图1中的步骤S104的具体流程图；

图3为本发明实施例提供的一种三维水力压裂及微地震事件正演模拟的方法的实施方式二的流程图；

图4为图3中的步骤S313的具体流程图；

图5为本发明实施例提供的一种三维水力压裂及微地震事件正演模拟的系统的实施方式一的结构框图；

图6为图5中的孔隙压力分布确定装置400的具体结构框图；

图7为本发明实施例提供的一种三维水力压裂及微地震事件正演模拟的系统的实施方式二的结构框图；

图8为本发明提供的一种三维水力压裂及微地震事件正演模拟的系统中微地震分布确定装置600的结构框图；

图9为100s时刻水力压裂缝三维空间分布示意图；

图10为500s时刻水力压裂缝三维空间分布示意图；

图11为水力压裂过程中三维孔隙压力分布示意图；

图12为孔隙压力分布的x-y剖面示意图；

图13为水力压裂过程中的微地震分布示意图；

图14为微地震事件分布的x-y剖面示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明提出的一种三维水力压裂及微地震事件正演模拟的方法，图1为该方法的具体流程图，由图1可知，所述的方法包括：

S101：采集当前地层的地层数据。

S102：根据当前地层裂缝的初始破裂位置、当前地层的区域范围设置网格。在具体的实施方式中，将裂缝的初始破裂位置设在射孔点上，射孔点可以有多个。根据所述的射孔点、当前地层的区域范围设置网格。

S103：根据基于压裂缝的流量守恒关系构建基于含压裂缝网格的流量守恒关系。

由于常规水力压裂过程基于压裂缝构建流量守恒关系，而裂缝宽度较小，在渗流模拟过程中无法与基质网格匹配，不能够进一步用于地层渗流模拟，因此本发明需要基于含压裂缝网格（包含压裂缝的网格单元）构建新的流量守恒关系。基于含压裂缝网格的流量守恒关系为：

dmass_f+dmass_l=△t·M_in

$\begin{array}{c}{\mathrm{dmass}}_f={\∫}_{A_f}(w(p_0^{*}+\mathrm{\Δ}{p}^{*},H\left(x\right))\·{\mathrm{\ρ}}_f(p_0^{*}+\mathrm{\Δ}{p}^{*})-w(p_0^{*},H\left(x\right))\·{\mathrm{\ρ}}_f\left(p_0^{*}\right))d{A}_f\\ +{\∫}_{V_f}(({\mathrm{\ρ}}_f(p_0^{*}+\mathrm{\Δ}{p}^{*})-{\mathrm{\ρ}}_f\left(p_0^{*}\right))\·\mathrm{\φ})d{V}_f\\ {\mathrm{dmass}}_l=\mathrm{\Δt}\·{\∫}_{A_f}\stackrel{----}{{\mathrm{\ρ}}_f{v}_D}d{A}_f\end{array}---\left(1\right)$

其中，dmass_f为含压裂缝网格中的新增流体质量，dmass_l为渗流到基质中的流体质量，p₀ ^*为n△t时刻的孔隙压力，p₀ ^*+△p^*为(n+1)△t时刻的孔隙压力，V_f为含压裂缝网格的体积，w为裂缝宽度，H(x)为网格对应的裂缝高度，M_in为压裂液的注入速率，A_f为裂缝网格与基质的接触面积，

为平均流体密度，为平均渗流速度，△t为模拟时间间隔，ρ_f为流体密度。

S104：根据基于含压裂缝网格的流量守恒关系确定当前地层的三维孔隙压力分布。图2为步骤S104的具体流程图，由图2可知，该步骤具体包括：

S201：根据基于含压裂缝网格的流量守恒关系、流体粘度、注入速率构建含压裂缝网格孔隙压力变化式。

由于裂缝宽度可以表示为含压裂缝网格孔隙压力的函数，结合注入速率、流体粘度，能够构建关于含压裂缝网格孔隙压力变化式，可用于直接计算该变化值。即上述的公式（1）最后可转化为含裂缝网格孔隙压力变化值的函数，如下所示：

$f(A_f,p,\mathrm{\Δp},H,\mathrm{\μ},w(p+\mathrm{\Δp}),\stackrel{\‾}{E})=M_{\mathrm{in}}\·\mathrm{\Δt}---\left(2\right)$

其中，A_f为裂缝网格与基质的接触面积，p为孔隙压力，△p为含压裂缝网格的孔隙压力变化量，H为网格对应的裂缝高度，μ为流体粘度，w为裂缝宽度，

为弹性参数，M_in为注入速率，△t为模拟时间间隔。

S202：根据所述的含压裂缝网格孔隙压力变化式确定压裂缝里的孔隙压力变化量；

S203：根据压裂缝里的孔隙压力变化量，结合渗流方程模拟当前地层的孔隙压力变化。渗流方程如下所示：

$-\frac{\∂}{\∂x}\left({\mathrm{\ρ}}_f{v}_{\mathrm{Dx}}\right)-\frac{\∂}{\∂y}\left({\mathrm{\ρ}}_f{v}_{\mathrm{Dy}}\right)-\frac{\∂}{\∂z}\left({\mathrm{\ρ}}_f{v}_{\mathrm{Dz}}\right)=\frac{\∂\left(\mathrm{\φ}{\mathrm{\ρ}}_f\right)}{\∂t}---\left(3\right)$

其中，v_Dx，v_Dy，v_Dy为x，y，z方向上的渗流速度。

为了保证整个模拟过程中流量守恒，取前一时刻和更新后孔隙压力的平均值赋予含裂缝网格，并结合渗流方程采用有限差分法模拟整个地层的孔隙压力变化。

S204：根据当前地层的孔隙压力变化构建三维孔隙压力分布。

得到裂缝网格的孔隙压力变化量后，结合渗流方程（公式3），采用有限差分算法模拟整个模型的孔隙压力变化，最后得到如图11所示的三维孔隙压力分布示意图。

由图1可知，该方法还包括：

S105：从所述的地层数据中获取网格所在位置的临界压力值；

S106：根据所述网格的孔隙压力值与网格所在位置的临界压力值生成当前地层的微地震分布。

S107：根据所述的微地震分布开发当前地层。

图3为本发明实施例提供的一种三维水力压裂及微地震事件正演模拟的方法的实施方式二的流程图，由图3可知，在实施方式二中，步骤S301至步骤S304与实施方式一中的步骤S101至步骤S104相同，此处不再赘述。该方法还包括：

S305：从所述的地层数据中获取最大地应力、最小地应力以及岩石临界应力强度因子。在具体的实施方式中，最大地应力可用σ₁来表示，最小地应力可用σ₃来表示，岩石临界应力强度因子可用K_IC来表示。

S306：从所述的网格上获取裂缝长度，获取裂缝传播方向与最大地应力的夹角，在具体的实施方式中，裂缝长度可通过L来表示，β为裂缝传播方向与最大地应力的夹角，该参数根据裂缝分布与应力方向可确定。

S307：根据所述的含压裂缝网格孔隙压力变化式确定孔隙流体压力，即根据公式（2）确定孔隙流体压力p。

S308：根据所述的最大地应力、最小地应力、孔隙流体压力、裂缝长度以及夹角确定裂缝尖端的应力强度因子。在具体的实施方式中，裂缝尖端的应力强度因子通过K_I来表示，则K_I(σ₁,σ₃,L,p)，具体的， $K_I=(\frac{{\mathrm{\σ}}_1-{\mathrm{\σ}}_3}{2}\cos \left(2\mathrm{\β}\right)+p-\frac{{\mathrm{\σ}}_1+{\mathrm{\σ}}_3}{2})\sqrt{\mathrm{\πL}}.$

S309：根据断裂力学理论，判断裂缝尖端的应力强度因子是否小于岩石临界应力强度因子，即判断K_I(σ₁,σ₃,L,p)<K_IC。

S310：当判断为否时，裂缝发生破裂，计算破裂后压裂缝内的平衡流体压力。在每一步模拟过程中，得到孔隙压力分布后，结合地应力和岩石物理参数，采用断裂力学作为判别条件判断压裂缝破裂。根据断裂力学理论，当裂缝尖端的应力强度因子达到岩石临界应力强度因子时，裂缝发生破裂。计算破裂后压裂缝内的平衡流体压力可通过如下公式进行：

$\begin{array}{c}{\&Integral;}_{A_1}(w(p_0^{*}+\mathrm{\&Delta;}{p}^{*},H_1\left(x\right))\&CenterDot;{\mathrm{\&rho;}}_f(p_0^{*}+\mathrm{\&Delta;}{p}^{*}))\mathrm{dA}+{\&Integral;}_{V_{f1}}({\mathrm{\&rho;}}_f(p_0^{*}+\mathrm{\&Delta;}{p}^{*})\&CenterDot;\mathrm{\&phi;})\mathrm{dV}+{\&Integral;}_{V_{f2}}({\mathrm{\&rho;}}_f\&CenterDot;\mathrm{\&phi;})\mathrm{dV}\\ ={\&Integral;}_{A_2}(w(p_b^{\#},H_2\left(x\right))\&CenterDot;{\mathrm{\&rho;}}_f\left(p_b^{\#}\right))\mathrm{dA}+{\&Integral;}_{V_{f1}+V_{f2}}({\mathrm{\&rho;}}_f\left(p_b^{\#}\right)\&CenterDot;\mathrm{\&phi;})\mathrm{dV}\end{array}---\left(4\right)$

其中，A₁为压裂缝破裂前的裂缝面积、V_f1为裂缝网格体积，H₁为裂缝高度，V_f2为新增的发生破裂的网格体积,A₂为破裂后的裂缝面积，H₂为破裂后的高度，p_b ^#为裂缝网格内平衡后的孔隙压力。图9和图10中的黑色网格展示压裂过程中的裂缝分布。

S311：根据破裂后压裂缝内的平衡流体压力更新所述的三维孔隙压力分布。

也即，模拟过程中，如果压裂缝不满足破裂条件，则前述步骤得到的孔隙压力、裂缝分布均保持不变；如果裂缝满足破裂条件，需要计算破裂后缝内的平衡流体压力，并用计算得到的平衡流体压力更新裂缝网格内的孔隙压力。

S312：从所述的地层数据中获取网格所在位置的临界压力值；

S313：根据所述网格的孔隙压力值与网格所在位置的临界压力值生成当前地层的微地震分布。图4为步骤S303的具体流程图。由图4可知，该步骤具体包括：

S401：根据所述的三维孔隙压力分布确定网格对应的孔隙压力值；

S402：依次判断各个网格的孔隙压力值是否大于该位置的临界压力值。即判断p(x,t)>C(x)。

水力压裂过程中，一部分压裂液在裂缝内流动，还有一部分会漏失到基质中。漏失的流体会提升地层孔隙压力和降低有效应力，也会引起微地震事件。本发明中采用临界压力理论模拟此类微地震事件，一旦孔隙流体压力p达到临界压力C，地层会发生微小的破裂或滑移并产生微地震信号。

如果网格对应的孔隙压力值为压裂缝里的孔隙压力值，由于该点位于压裂缝内，已经发生破裂，认为该位置不再发生破裂；如果网格对应的孔隙压力值为压裂缝外的孔隙压力值，则判断压裂缝外的孔隙压力值是否大于该位置的临界压力值。

S403：当判断为是时，所述网格所在位置产生孔隙压力控制型的微地震事件。

S404：所述的微地震事件组成微地震分布。多个网格对应的孔隙压力控制型的微地震事件组成微地震分布。如图13和图14所示，图中的菱形、圆形和星号网格展示孔隙压力控制型的微地震分布。图11、图12为孔隙压力分布。

由图3可知，该方法还包括：

S314：根据所述的微地震分布开发当前地层。

在本发明中，射孔点、流体粘度以及注入速率是可调的，因此，设定一组射孔点、流体粘度以及注入速率组成的压裂参数，即可根据本发明的方法得到微地震分布，改变压裂参数即可得到另一个微地震分布。如此，在多个微地震分布中可以选出优选微地震分布，对应的压裂参数即为优选压裂参数。该组优选压裂参数（射孔点、流体粘度、注入速率）即可用于指导当前地层的开发。

如上所述，即为本发明提供的一种三维水力压裂及微地震事件正演模拟的方法，能够同步模拟压裂缝形态、裂缝孔隙压力、基质孔隙压力及基质中微地震事件。

本发明还提出的一种三维水力压裂及微地震事件正演模拟的系统，图5为该系统的实施方式一的结构框图，由图5可知，所述的系统包括：

地层数据采集装置100，用于采集当前地层的地层数据。

网格设置装置200，用于根据当前地层裂缝的初始破裂位置、当前地层的区域范围设置网格，在具体的实施方式中，将裂缝的初始破裂位置设在射孔点上，射孔点可以有多个。根据所述的射孔点、当前地层的区域范围设置网格。

流量守恒构建装置300，用于根据基于压裂缝的流量守恒关系构建基于含压裂缝网格的流量守恒关系。

由于常规水力压裂过程基于压裂缝构建流量守恒关系，而裂缝宽度较小，在渗流模拟过程中无法与基质网格匹配，不能够进一步用于地层渗流模拟，因此本发明需要基于含压裂缝网格（包含压裂缝的网格单元）构建新的流量守恒关系。基于含压裂缝网格的流量守恒关系为：

dmass_f+dmass_l=△t·M_in

$\begin{array}{c}{\mathrm{dmass}}_f={\&Integral;}_{A_f}(w(p_0^{*}+\mathrm{\&Delta;}{p}^{*},H\left(x\right))\&CenterDot;{\mathrm{\&rho;}}_f(p_0^{*}+\mathrm{\&Delta;}{p}^{*})-w(p_0^{*},H\left(x\right))\&CenterDot;{\mathrm{\&rho;}}_f\left(p_0^{*}\right))d{A}_f\\ +{\&Integral;}_{V_f}(({\mathrm{\&rho;}}_f(p_0^{*}+\mathrm{\&Delta;}{p}^{*})-{\mathrm{\&rho;}}_f\left(p_0^{*}\right))\&CenterDot;\mathrm{\&phi;})d{V}_f\\ {\mathrm{dmass}}_l=\mathrm{\&Delta;t}\&CenterDot;{\&Integral;}_{A_f}\stackrel{----}{{\mathrm{\&rho;}}_f{v}_D}d{A}_f\end{array}---\left(1\right)$

其中，dmass_f为含压裂缝网格中的新增流体质量，dmass_l为渗流到基质中的流体质量，p₀ ^*为n△t时刻的孔隙压力，p₀ ^*+△p^*为(n+1)△t时刻的孔隙压力，V_f为含压裂缝网格的体积，w为裂缝宽度，H(x)为网格对应的裂缝高度，M_in为压裂液的注入速率，A_f为裂缝网格与基质的接触面积，

为平均流体密度，

为平均渗流速度，△t为模拟时间间，ρ_f为流体密度。

孔隙压力分布确定装置400，用于根据基于含压裂缝网格的流量守恒关系确定当前地层的三维孔隙压力分布。图6为孔隙压力分布确定装置400的具体结构框图，由图6可知，孔隙压力分布确定装置具体包括：

压力变化式构建模块401，用于根据基于含压裂缝网格的流量守恒关系、流体粘度、注入速率构建含压裂缝网格孔隙压力变化式；

由于裂缝宽度可以表示为含压裂缝网格孔隙压力的函数，结合注入速率、流体粘度，能够构建关于含压裂缝网格孔隙压力变化式，可用于直接计算该变化值。即上述的公式（1）最后可转化为含裂缝网格孔隙压力变化值的函数，如下所示：

$f(A_f,p,\mathrm{\&Delta;p},H,\mathrm{\&mu;},w(p+\mathrm{\&Delta;p}),\stackrel{\&OverBar;}{E})=M_{\mathrm{in}}\&CenterDot;\mathrm{\&Delta;t}---\left(2\right)$

其中，A_f为裂缝网格与基质的接触面积，p为孔隙压力，△p为含压裂缝网格的孔隙压力变化量，H为网格对应的裂缝高度，μ为流体粘度，w为裂缝宽度，E为弹性参数，M_in为注入速率，△t为模拟时间间隔。

孔隙压力变化量确定模块402，用于根据所述的含压裂缝网格孔隙压力变化式确定压裂缝里的孔隙压力变化量；

模拟模块403，用于根据压裂缝里的孔隙压力变化量，结合渗流方程模拟当前地层的孔隙压力变化。渗流方程如下所示：

$-\frac{\&PartialD;}{\&PartialD;x}\left({\mathrm{\&rho;}}_f{v}_{\mathrm{Dx}}\right)-\frac{\&PartialD;}{\&PartialD;y}\left({\mathrm{\&rho;}}_f{v}_{\mathrm{Dy}}\right)-\frac{\&PartialD;}{\&PartialD;z}\left({\mathrm{\&rho;}}_f{v}_{\mathrm{Dz}}\right)=\frac{\&PartialD;\left(\mathrm{\&phi;}{\mathrm{\&rho;}}_f\right)}{\&PartialD;t}---\left(3\right)$

其中，v_Dx，v_Dy，v_Dy为x，y，z方向上的渗流速度。

为了保证整个模拟过程中流量守恒，取前一时刻和更新后孔隙压力的平均值赋予含裂缝网格，并结合渗流方程采用有限差分法模拟整个地层的孔隙压力变化。

三维孔隙压力分布构建模块404，用于根据当前地层的孔隙压力变化构建三维孔隙压力分布。

得到裂缝网格的孔隙压力变化量后，结合渗流方程（公式3），采用有限差分算法模拟整个模型的孔隙压力变化，最后得到如图9所示的三维孔隙压力分布。

由图5可知，该系统还包括：

临界压力值获取装置500，用于从所述的地层数据中获取网格所在位置的临界压力值；

微地震分布确定装置600，用于根据所述三维孔隙压力分布与网格所在位置的临界压力值生成当前地层的微地震分布；

地层开发装置700，用于根据所述的微地震分布开发当前地层。

图7为本发明实施例提供的一种三维水力压裂及微地震事件正演模拟的系统的实施方式二的结构框图，由图7可知，在实施方式二中，该系统还包括：

岩石临界应力强度因子获取装置800，用于从所述的地层数据中获取最大地应力、最小地应力以及岩石临界应力强度因子。在具体的实施方式中，最大地应力可用σ₁来表示，最小地应力可用σ₃来表示，岩石临界应力强度因子可用K_IC来表示。

裂缝长度获取装置900，用于从所述的网格上获取裂缝长度，获取裂缝传播方向与最大地应力的夹角，在具体的实施方式中，裂缝长度可通过L来表示。β为裂缝传播方向与最大地应力的夹角，该参数根据裂缝分布与应力方向可确定。

孔隙流体压力确定装置1000，用于根据所述的含压裂缝网格孔隙压力变化式确定孔隙流体压力，即根据公式（2）确定孔隙流体压力p。

应力强度因子确定装置1100，用于根据所述的最大地应力、最小地应力、孔隙流体压力、裂缝长度以及夹角确定裂缝尖端的应力强度因子。在具体的实施方式中，裂缝尖端的应力强度因子通过K_I来表示，则K_I(σ₁,σ₃,L,p)。

具体的， $K_I=(\frac{{\mathrm{\&sigma;}}_1-{\mathrm{\&sigma;}}_3}{2}\cos \left(2\mathrm{\&beta;}\right)+p-\frac{{\mathrm{\&sigma;}}_1+{\mathrm{\&sigma;}}_3}{2})\sqrt{\mathrm{\&pi;L}}.$

判断装置1200，用于根据断裂力学理论，判断裂缝尖端的应力强度因子是否小于岩石临界应力强度因子，即判断K_I(σ₁,σ₃,L,_p)<K_IC。

平衡流体压力计算装置1300，用于当所述的判断装置判断为否时，裂缝发生破裂，计算破裂后压裂缝内的平衡流体压力。在每一步模拟过程中，得到孔隙压力分布后，结合地应力和岩石物理参数，采用断裂力学作为判别条件判断压裂缝破裂。根据断裂力学理论，当裂缝尖端的应力强度因子达到岩石临界应力强度因子时，裂缝发生破裂。计算破裂后压裂缝内的平衡流体压力可通过如下公式进行：

$\begin{array}{c}{\&Integral;}_{A_1}(w(p_0^{*}+\mathrm{\&Delta;}{p}^{*},H_1\left(x\right))\&CenterDot;{\mathrm{\&rho;}}_f(p_0^{*}+\mathrm{\&Delta;}{p}^{*}))\mathrm{dA}+{\&Integral;}_{V_{f1}}({\mathrm{\&rho;}}_f(p_0^{*}+\mathrm{\&Delta;}{p}^{*})\&CenterDot;\mathrm{\&phi;})\mathrm{dV}+{\&Integral;}_{V_{f2}}({\mathrm{\&rho;}}_f\&CenterDot;\mathrm{\&phi;})\mathrm{dV}\\ ={\&Integral;}_{A_2}(w(p_b^{\#},H_2\left(x\right))\&CenterDot;{\mathrm{\&rho;}}_f\left(p_b^{\#}\right))\mathrm{dA}+{\&Integral;}_{V_{f1}+V_{f2}}({\mathrm{\&rho;}}_f\left(p_b^{\#}\right)\&CenterDot;\mathrm{\&phi;})\mathrm{dV}\end{array}---\left(4\right)$

其中，A₁为压裂缝破裂前的裂缝面积、V_f1为裂缝网格体积，H₁为裂缝高度，V_f2为新增的发生破裂的网格体积,A₂为破裂后的裂缝面积，H₂为破裂后的高度，p_b ^#为裂缝网格内平衡后的孔隙压力。图9和图10中黑色-网格展示压裂过程中的裂缝分布。

三维孔隙压力分布更新装置1400，用于根据破裂后压裂缝内的平衡流体压力更新所述的三维孔隙压力分布。

也即，模拟过程中，如果压裂缝不满足破裂条件，则前述步骤得到的孔隙压力、裂缝分布均保持不变；如果裂缝满足破裂条件，需要计算破裂后缝内的平衡流体压力，并用计算得到的平衡流体压力更新裂缝网格内的孔隙压力。

图8为本发明提供的一种三维水力压裂及微地震事件正演模拟的系统中微地震分布确定装置600的结构框图。由图8可知，微地震分布确定装置具体包括：

孔隙压力值确定模块601，用于根据所述的三维孔隙压力分布确定网格对应的孔隙压力值；

判断模块602，用于依次判断各个网格的孔隙压力值是否大于该位置的临界压力值。即判断p(x,t)>C(x)。

水力压裂过程中，一部分压裂液在裂缝内流动，还有一部分会漏失到基质中。漏失的流体会提升地层孔隙压力和降低有效应力，也会引起微地震事件。本发明中采用临界压力理论模拟此类微地震事件，一旦孔隙流体压力p达到临界压力C，地层会发生微小的破裂或滑移并产生微地震信号。

如果网格对应的孔隙压力值为压裂缝里的孔隙压力值，由于该点位于压裂缝内，已经发生破裂，认为该位置不再发生破裂；如果网格对应的孔隙压力值为压裂缝外的孔隙压力值，则判断压裂缝外的孔隙压力值是否大于该位置的临界压力值。

微地震事件产生模块603，用于当所述的判断模块判断为是时，所述网格所在位置产生孔隙压力控制型的微地震事件。

微地震分布组成模块604，用于所述的微地震事件组成微地震分布。多个网格对应的孔隙压力控制型的微地震事件组成微地震分布。如图13、图14所示，图中菱形、圆形和星号网格展示孔隙压力控制型的微地震分布。图11、图12为孔隙压力分布。

在本发明中，射孔点、流体粘度以及注入速率是可调的，因此，设定一组射孔点、流体粘度以及注入速率组成的压裂参数，即可根据本发明的方法得到微地震分布，改变压裂参数即可得到另一个微地震分布。如此，在多个微地震分布中可以选出优选微地震分布，对应的压裂参数即为优选压裂参数。该组优选压裂参数（射孔点、流体粘度、注入速率）即可用于指导当前地层的开发。

如上所述，即为本发明提供的一种三维水力压裂及微地震事件正演模拟的系统，能够同步模拟压裂缝形态、裂缝孔隙压力、基质孔隙压力及基质中微地震事件。

下面结合具体的实施例，详细介绍本发明的技术方案。采集国内某地层的地层数据，对其进行三维水力压裂及微地震事件正演模拟分析。

图9为100s时刻水力压裂缝三维空间分布示意图，图10为500s时刻水力压裂缝三维空间分布示意图。由图9、图10可知，由于上下存在隔层，压裂缝沿水平方向延伸。

图11为水力压裂过程中三维孔隙压力分布示意图，图12为孔隙压力分布的x-y剖面示意图。

图13为水力压裂过程中的微地震分布示意图，图14为微地震事件分布的x-y剖面示意图，由图13、图14可知，图中黑色点状网格为压裂缝，菱形、圆形和星号网格分别为不同时刻孔隙压力控制型的微地震事件。

由此可知，本发明的技术方案能够同步模拟压裂缝破裂、缝内孔隙压力变化、基质孔隙压力变化和孔隙压力控制型的微地震事件，获得整个地层的微地震事件和流体压力分布，为进一步的微地震研究和地层参数研究提供基础。

综上所述，本发明提出的一种三维水力压裂及微地震事件正演模拟的方法及系统，将常规水力压裂基于压裂缝的流量守恒关系推导为基于含压裂缝网格的形式；根据质量守恒定律，构建关于含压裂缝网格孔隙压力变化的等式，利用该公式可直接计算含压裂缝网格的孔隙压力变化。得到该孔隙压力变化值后，结合有限差分法求解渗流方程模拟基质中孔隙压力变化；根据更新后的孔隙压力、地应力和岩石参数，采用断裂力学理论模拟压裂缝破裂，破裂位置对应微地震事件；同时根据临界压力理论模拟压裂缝周边地层由于孔隙压力升高产生的微地震事件（孔隙压力控制型微地震）。本发明拓展了常规水力压裂模拟技术，在获得压裂缝信息之外，还能够得到整个地层的孔隙压力分布；同时本发明在水力压裂模拟基础上，融入孔隙压力控制型微地震的正演模拟，构建了一种水力压裂及微地震综合模拟技术，能够更好的表征压裂过程中微地震的变化特征。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一般计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体（Read-Only Memory，ROM）或随机存储记忆体（Random AccessMemory，RAM）等。

本领域技术人员还可以了解到本发明实施例列出的各种功能是通过硬件还是软件来实现取决于特定的应用和整个系统的设计要求。本领域技术人员可以对于每种特定的应用，可以使用各种方法实现所述的功能，但这种实现不应被理解为超出本发明实施例保护的范围。

本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (14)

1.一种三维水力压裂及微地震事件正演模拟的方法，其特征是，所述的方法具体包括：

采集当前地层的地层数据；

根据当前地层裂缝的初始破裂位置、当前地层的区域范围设置网格；

根据基于压裂缝的流量守恒关系构建基于含压裂缝网格的流量守恒关系；

根据基于含压裂缝网格的流量守恒关系确定当前地层的三维孔隙压力分布；

从所述的地层数据中获取网格所在位置的临界压力值；

根据所述三维孔隙压力分布与网格所在位置的临界压力值生成当前地层的微地震分布；

根据所述的微地震分布开发当前地层。

2.根据权利要求1所述的方法，其特征是，基于含压裂缝网格的流量守恒关系为：

dmass_f+dmass_l=△t·M_in

$\begin{array}{c}{\mathrm{dmass}}_f={\&Integral;}_{A_f}(w(p_0^{*}+\mathrm{\&Delta;}{p}^{*},H\left(x\right))\&CenterDot;{\mathrm{\&rho;}}_f(p_0^{*}+\mathrm{\&Delta;}{p}^{*})-w(p_0^{*},H\left(x\right))\&CenterDot;{\mathrm{\&rho;}}_f\left(p_0^{*}\right))d{A}_f\\ +{\&Integral;}_{V_f}(({\mathrm{\&rho;}}_f(p_0^{*}+\mathrm{\&Delta;}{p}^{*})-{\mathrm{\&rho;}}_f\left(p_0^{*}\right))\&CenterDot;\mathrm{\&phi;})d{V}_f\\ {\mathrm{dmass}}_l=\mathrm{\&Delta;t}\&CenterDot;{\&Integral;}_{A_f}\stackrel{----}{{\mathrm{\&rho;}}_f{v}_D}d{A}_f\end{array}$

其中，dmass_f为含压裂缝网格中的新增流体质量，dmass_l为渗流到基质中的流体质量，p₀ ^*为n△t时刻的孔隙压力，p₀ ^*+△p^*为(n+1)△t时刻的孔隙压力，V_f为含压裂缝网格的体积，w为裂缝宽度，H(x)为网格对应的裂缝高度，M_in为压裂液的注入速率，A_f为裂缝网格与基质的接触面积，

为平均流体密度，

为平均渗流速度，△t为模拟时间间隔，ρ_f为流体密度。

3.根据权利要求1所述的方法，其特征是，根据基于含压裂缝网格的流量守恒关系确定当前地层的三维孔隙压力分布具体包括：

根据基于含压裂缝网格的流量守恒关系、流体粘度、注入速率构建含压裂缝网格孔隙压力变化式；

根据所述的含压裂缝网格孔隙压力变化式确定压裂缝里的孔隙压力变化量；

根据压裂缝里的孔隙压力变化量，结合渗流方程模拟当前地层的孔隙压力变化；

根据当前地层的孔隙压力变化构建三维孔隙压力分布。

4.根据权利要求3所述的方法，其特征是，所述的含压裂缝网格孔隙压力变化式为：

$f(A_f,p,\mathrm{\&Delta;p},H,\mathrm{\&mu;},w(p+\mathrm{\&Delta;p}),\stackrel{\&OverBar;}{E})=M_{\mathrm{in}}\&CenterDot;\mathrm{\&Delta;t}$

其中，A_f为裂缝网格与基质的接触面积，p为孔隙流体压力，△p为含压裂缝网格的孔隙压力变化量，H为网格对应的裂缝高度，μ为流体粘度，w为裂缝宽度，

为弹性参数，M_in为注入速率，△t为模拟时间间隔；

所述的渗流方程为：

$-\frac{\&PartialD;}{\&PartialD;x}\left({\mathrm{\&rho;}}_f{v}_{\mathrm{Dx}}\right)-\frac{\&PartialD;}{\&PartialD;y}\left({\mathrm{\&rho;}}_f{v}_{\mathrm{Dy}}\right)-\frac{\&PartialD;}{\&PartialD;z}\left({\mathrm{\&rho;}}_f{v}_{\mathrm{Dz}}\right)=\frac{\&PartialD;\left(\mathrm{\&phi;}{\mathrm{\&rho;}}_f\right)}{\&PartialD;t}$

其中，v_Dx，v_Dy，v_Dy为x，y，z方向上的渗流速度。

5.根据权利要求1或3所述的方法，其特征是，所述的方法还包括：

从所述的地层数据中获取最大地应力、最小地应力以及岩石临界应力强度因子；

从所述的网格上获取裂缝长度，获取裂缝传播方向与最大地应力的夹角；

根据所述的含压裂缝网格孔隙压力变化式确定孔隙流体压力；

根据所述的最大地应力、最小地应力、孔隙流体压力、裂缝长度以及夹角确定裂缝尖端的应力强度因子；

根据断裂力学理论，判断裂缝尖端的应力强度因子是否小于岩石临界应力强度因子；

当判断为否时，裂缝发生破裂，计算破裂后压裂缝内的平衡流体压力；

根据破裂后压裂缝内的平衡流体压力更新所述的三维孔隙压力分布。

6.根据权利要求5所述的方法，其特征是，计算破裂后压裂缝内的平衡流体压力通过如下公式进行：

$\begin{array}{c}{\&Integral;}_{A_1}(w(p_0^{*}+\mathrm{\&Delta;}{p}^{*},H_1\left(x\right))\&CenterDot;{\mathrm{\&rho;}}_f(p_0^{*}+\mathrm{\&Delta;}{p}^{*}))\mathrm{dA}+{\&Integral;}_{V_{f1}}({\mathrm{\&rho;}}_f(p_0^{*}+\mathrm{\&Delta;}{p}^{*})\&CenterDot;\mathrm{\&phi;})\mathrm{dV}+{\&Integral;}_{V_{f2}}({\mathrm{\&rho;}}_f\&CenterDot;\mathrm{\&phi;})\mathrm{dV}\\ ={\&Integral;}_{A_2}(w(p_b^{\#},H_2\left(x\right))\&CenterDot;{\mathrm{\&rho;}}_f\left(p_b^{\#}\right))\mathrm{dA}+{\&Integral;}_{V_{f1}+V_{f2}}({\mathrm{\&rho;}}_f\left(p_b^{\#}\right)\&CenterDot;\mathrm{\&phi;})\mathrm{dV}\end{array}$

其中，A₁为压裂缝破裂前的裂缝面积、V_f1为裂缝网格体积，H₁为裂缝高度，V_f2为新增的发生破裂的网格体积,A₂为破裂后的裂缝面积，H₂为破裂后的高度，p_b ^#为裂缝网格内平衡后的孔隙压力。

7.根据权利要求5所述的方法，其特征是，根据所述三维孔隙压力分布与网格所在位置的临界压力值产生孔隙压力控制型的微地震事件具体包括：

根据所述的三维孔隙压力分布确定网格对应的孔隙压力值；

依次判断各个网格的孔隙压力值是否大于该位置的临界压力值；

当判断为是时，所述网格所在位置产生孔隙压力控制型的微地震事件；

所述的微地震事件组成微地震分布。

8.一种三维水力压裂及微地震事件正演模拟的系统，其特征是，所述的系统具体包括：

地层数据采集装置，用于采集当前地层的地层数据；

网格设置装置，用于根据当前地层裂缝的初始破裂位置、当前地层的区域范围设置网格；

流量守恒构建装置，用于根据基于压裂缝的流量守恒关系构建基于含压裂缝网格的流量守恒关系；

孔隙压力分布确定装置，用于根据基于含压裂缝网格的流量守恒关系确定当前地层的三维孔隙压力分布；

临界压力值获取装置，用于从所述的地层数据中获取网格所在位置的临界压力值；

微地震分布确定装置，用于根据所述三维孔隙压力分布与网格所在位置的临界压力值生成当前地层的微地震分布；

地层开发装置，用于根据所述的微地震分布开发当前地层。

9.根据权利要求8所述的系统，其特征是，基于含压裂缝网格的流量守恒关系为：

dmass_f+dmass_l=△t·M_in

$\begin{array}{c}{\mathrm{dmass}}_f={\&Integral;}_{A_f}(w(p_0^{*}+\mathrm{\&Delta;}{p}^{*},H\left(x\right))\&CenterDot;{\mathrm{\&rho;}}_f(p_0^{*}+\mathrm{\&Delta;}{p}^{*})-w(p_0^{*},H\left(x\right))\&CenterDot;{\mathrm{\&rho;}}_f\left(p_0^{*}\right))d{A}_f\\ +{\&Integral;}_{V_f}(({\mathrm{\&rho;}}_f(p_0^{*}+\mathrm{\&Delta;}{p}^{*})-{\mathrm{\&rho;}}_f\left(p_0^{*}\right))\&CenterDot;\mathrm{\&phi;})d{V}_f\\ {\mathrm{dmass}}_l=\mathrm{\&Delta;t}\&CenterDot;{\&Integral;}_{A_f}\stackrel{----}{{\mathrm{\&rho;}}_f{v}_D}d{A}_f\end{array}$

其中，dmass_f为含压裂缝网格中的新增流体质量，dmass_l为渗流到基质中的流体质量，p₀ ^*为n△t时刻的孔隙压力，p₀ ^*+△p^*为(n+1)△t时刻的孔隙压力，V_f为含压裂缝网格的体积，w为裂缝宽度，H(x)为网格对应的裂缝高度，M_in为压裂液的注入速率，A_f为裂缝网格与基质的接触面积，

为平均流体密度，为平均渗流速度，△t为模拟时间间隔，ρ_f为流体密度。

10.根据权利要求8所述的系统，其特征是，所述的孔隙压力分布确定装置具体包括：

压力变化式构建模块，用于根据基于含压裂缝网格的流量守恒关系、流体粘度、注入速率构建含压裂缝网格孔隙压力变化式；

孔隙压力变化量确定模块，用于根据所述的含压裂缝网格孔隙压力变化式确定压裂缝里的孔隙压力变化量；

模拟模块，用于根据压裂缝里的孔隙压力变化量，结合渗流方程模拟当前地层的孔隙压力变化；

三维孔隙压力分布构建模块，用于根据当前地层的孔隙压力变化构建三维孔隙压力分布。

11.根据权利要求10所述的系统，其特征是，所述的含压裂缝网格孔隙压力变化式为：

$f(A_f,p,\mathrm{\&Delta;p},H,\mathrm{\&mu;},w(p+\mathrm{\&Delta;p}),\stackrel{\&OverBar;}{E})=M_{\mathrm{in}}\&CenterDot;\mathrm{\&Delta;t}$

其中，A_f为裂缝网格与基质的接触面积，p为孔隙流体压力，△p为含压裂缝网格的孔隙压力变化量，H为网格对应的裂缝高度，μ为流体粘度，w为裂缝宽度，

为弹性参数，M_in为注入速率，△t为模拟时间间隔；

所述的渗流方程为：

$-\frac{\&PartialD;}{\&PartialD;x}\left({\mathrm{\&rho;}}_f{v}_{\mathrm{Dx}}\right)-\frac{\&PartialD;}{\&PartialD;y}\left({\mathrm{\&rho;}}_f{v}_{\mathrm{Dy}}\right)-\frac{\&PartialD;}{\&PartialD;z}\left({\mathrm{\&rho;}}_f{v}_{\mathrm{Dz}}\right)=\frac{\&PartialD;\left(\mathrm{\&phi;}{\mathrm{\&rho;}}_f\right)}{\&PartialD;t}$

其中，v_Dx，v_Dy，v_Dy为x，y，z方向上的渗流速度。

12.根据权利要求8或11所述的系统，其特征是，所述的系统还包括：

岩石临界应力强度因子获取装置，用于从所述的地层数据中获取最大地应力、最小地应力以及岩石临界应力强度因子；

裂缝长度获取装置，用于从所述的网格上获取裂缝长度，获取裂缝传播方向与最大地应力的夹角；

孔隙流体压力确定装置，用于根据所述的含压裂缝网格孔隙压力变化式确定孔隙流体压力；

应力强度因子确定装置，用于根据所述的最大地应力、最小地应力、孔隙流体压力、裂缝长度以及夹角确定裂缝尖端的应力强度因子；

判断装置，用于根据断裂力学理论，判断裂缝尖端的应力强度因子是否小于岩石临界应力强度因子；

平衡流体压力计算装置，用于当所述的判断装置判断为否时，裂缝发生破裂，计算破裂后压裂缝内的平衡流体压力；

三维孔隙压力分布更新装置，用于根据破裂后压裂缝内的平衡流体压力更新所述的三维孔隙压力分布。

13.根据权利要求12所述的系统，其特征是，所述的平衡流体压力计算装置通过如下公式计算破裂后压裂缝内的平衡流体压力：

$\begin{array}{c}{\&Integral;}_{A_1}(w(p_0^{*}+\mathrm{\&Delta;}{p}^{*},H_1\left(x\right))\&CenterDot;{\mathrm{\&rho;}}_f(p_0^{*}+\mathrm{\&Delta;}{p}^{*}))\mathrm{dA}+{\&Integral;}_{V_{f1}}({\mathrm{\&rho;}}_f(p_0^{*}+\mathrm{\&Delta;}{p}^{*})\&CenterDot;\mathrm{\&phi;})\mathrm{dV}+{\&Integral;}_{V_{f2}}({\mathrm{\&rho;}}_f\&CenterDot;\mathrm{\&phi;})\mathrm{dV}\\ ={\&Integral;}_{A_2}(w(p_b^{\#},H_2\left(x\right))\&CenterDot;{\mathrm{\&rho;}}_f\left(p_b^{\#}\right))\mathrm{dA}+{\&Integral;}_{V_{f1}+V_{f2}}({\mathrm{\&rho;}}_f\left(p_b^{\#}\right)\&CenterDot;\mathrm{\&phi;})\mathrm{dV}\end{array}$

其中，A₁为压裂缝破裂前的裂缝面积、V_f1为裂缝网格体积，H₁为裂缝高度，V_f2为新增的发生破裂的网格体积,A₂为破裂后的裂缝面积，H₂为破裂后的高度，_pb ^#为裂缝网格内平衡后的孔隙压力。

14.根据权利要求12所述的系统，其特征是，所述的微地震分布确定装置具体包括：

孔隙压力值确定模块，用于根据所述的三维孔隙压力分布确定网格对应的孔隙压力值；

判断模块，用于依次判断各个网格的孔隙压力值是否大于该位置的临界压力值；

微地震事件产生模块，用于当所述的判断模块判断为是时，所述网格所在位置产生孔隙压力控制型的微地震事件；

微地震分布组成模块，用于将所述的微地震事件组成微地震分布。

Images