0 引言

Maverhofer等2008年第1次正式地提出了改造油藏体积的概念，并在以往微地震研究的基础上对改造体积(SRV)进行了计算^[1]。SRV即Stimulated reservoir volume，常称为压裂体积或改造体积，是通过压裂的方式对储层实施改造，形成天然裂缝和人工裂缝相互交错的裂缝网络，极大地改造储集层有效泄油体积，从而提高储层动用率和最终采收率^[2]。大量研究表明，压裂油藏的油气生产力不仅取决于被识别的水力压裂，还与改造体积有关，储层改造体积越大，增产效果越明显，储层的改造体积与增产效果具有显著的正相关性^[3]。Fisher等用通道长度和通道宽度来表征裂缝扩展的长度和宽度^[4]。Warpinski等用缝网长度、缝网宽度及缝高来表征SRV^[5-6]。此外，应用裂缝半长和裂缝间距表征SRV的宽度，并用水平井长度表征SRV长度来计算SRV也较为常用，以上传统的SRV计算方法均是用立方体来拟合SRV，故拟合结果较为粗略^[7]。目前，主要有离散格子法^[2]、凸包^[8]等算法拟合微地震事件点集获得微地震SRV。

水力压裂利用高压向致密低渗透储层注入流体，造成裂缝周围储层的张裂或错动，各种张裂或错动会形成相互交错的裂缝网络，从而增加改造体积，同时，各种张裂或错动会向外辐射弹性地震波能量^[9]。微地震监测技术是通过观测、分析生产活动中所产生的微小地震事件来监测生产活动的影响、效果及地下状态的地球物理技术，故微地震监测技术可以通过监测和收集地震波信号来定位微地震事件点^[10]。传统PKN、GDK等水力压裂模型，仅考虑了裂缝体积，却忽略了压裂液在基质中的漏失，导致对改造区域及改造体积的认识均不准确^[11-12]。

本文首先基于物质平衡、压裂缝网格(包含压裂缝和基质的网格)、渗流方程和断裂力学来拟合孔隙压力的变化和裂缝的扩展，同时拟合裂缝扩展，以及裂缝和基质中的孔隙压力变化^[13]；然后，根据水力压裂模型得到的三维孔隙压力分布，依据断裂准则和临界孔隙压力分别判断主裂缝破裂和流体漏失引起的微地震事件(P型微地震事件)，从而确定微地震事件所占网格体积即为水力压裂SRV^[14]。这些微地震事件位于产油层，以微地震事件点为基础可构建SRV。

通常，微地震定位方法构建的速度模型只是压裂前的地层速度，未考虑水力压裂过程引起的地层速度变化和速度分析中的误差^[15]，因此，微地震正演过程中检波器接收到的初至旅行时会有偏差。初至旅行时的偏差会影响微地震事件点集的确定，使得计算的储层改造体积与现场监测结果相差较大。

本文在微地震定位的反演过程中通过对初至旅行时进行不同程度地加噪处理，来反映水力压裂对微地震定位的影响。由于加噪后得到的微地震事件点必定不准确，基于小概率事件概念，通过去除微地震事件点集中的异常点，尽量减少噪声对反演定位以及拟合SRV区域的影响。最后，相较于长方体拟合微地震事件点的传统方法，采用更为精确的面元(BIN)、三维狄洛尼三角剖分和最小体积覆盖椭球(Minimum volume enclosing ellipsoid，MVEE)算法拟合逼近SRV区域，并通过对比去异常点前后各拟合算法的拟合结果，对上述几种基于微地震事件点拟合SRV的算法进行了比较分析。

1 水力压裂模拟及压裂过程中产生的微地震信号模拟

建立三维油藏模型模拟水力压裂过程，包括主裂缝扩展和压裂液由裂缝向基质扩散过程。渗流模型假设微可压缩流体，压裂模拟过程中，压裂液除了渗流到裂缝中，还有一部分会漏失到基质中，根据物质平衡方程和渗流方程可以模拟压裂层的三维孔隙压力分布。

渗流控制方程(物质平衡方程)为

$ \dfrac{\partial }{{\partial t}}\left( {\phi \rho } \right) + \nabla\left( {\rho u} \right) = {M_{\text{in}}} $

(1)

式中：

$\phi$—孔隙度，无因次；

$\rho$—流体密度，g/cm$^{3}$；

$u$—渗流速度，m/s；

$M_{\text{in}}$—压裂液注入速率^[14-16]，m$^{3}$/s。

利用达西方程描述渗流速度(渗流方程)

$ u = - \dfrac{K}{\mu }\nabla p $

(2)

式中：

$K$—渗透率，mD；

$p$—孔隙流体压力，MPa；

$\mu$—流体黏度，mPa$\cdot$s。

对上述模型划分网格，从而进行数值离散。采用有限体积法、全隐式格式进行求解，求解变量为各个网格的压力($p$)。

本文在压裂过程中假设地应力及岩石性质不变，依据应力强度因子断裂准则(K准则)判断主裂缝扩展。通过研究裂纹问题所采用的线弹性力学方法，对裂纹尖端附近区域的应力状态进行研究，可以得到裂纹尖端附近各点的应力分量，并引出了“应力强度因子$K_{\text{I}}$”概念。

在一个双轴应力场中(图 1)，假设$x$方向应力为$\sigma_{\alpha}$，$y$方向应力为$\sigma$，存在一个长为2$a$倾斜的裂缝，岩石裂缝尖端的应力强度因子可表示为^[17]

$ K_{\text{I}} = \sigma \sqrt{{\rm{ \mathsf{ π} }} a} \left (\cos^{2} \beta + \alpha \sin ^{2} \beta \right ) $

(3)

式中：

$\beta$—裂缝与$x$轴之间的夹角，°。

在岩石力学中要考虑孔隙压力的影响，通过式(3)，可推导出岩石裂缝尖端的应力强度因子^[18]

$ K_{\text{I}} = \left (\dfrac{\sigma_{1}-\sigma_{3}}{2} \cos{2 \alpha} + p-\dfrac{\sigma_{1} + \sigma_{3}}{2} \right)\sqrt{{\rm{ \mathsf{ π} }} L} $

(4)

式中：

$L$—裂缝长度，m；

$\sigma_{1}$，$\sigma_{3}$—最大主应力和最小主应力，MPa；

$\alpha$—裂缝面与最大主应力之间的夹角，°。

依据应力强度因子断裂准则(K准则)，裂缝扩展的临界条件为

$ K_{\text{I}} = K_{\text{IC}} $

(5)

式中：

$K_{\text{IC}}$—岩石的断裂韧性，无因次。

压裂过程中的微地震信号可分为两类：P类(Pore pressure diffusion controlled type)，该类微地震信号受流体压力控制；H类(Hydraulic fracturing controlled type)，该类微地震信号由压裂缝破裂引起。在水力压裂过程中，压裂液渗流到基质中引起的微地震事件均属于P类微地震信号^[19-21]。Rothert and Humme等采用临界孔隙压力做为P类微地震的判断准则^[22-23]，如果孔隙流体压力($p$)超过临界值($p_{\text{C}}$)，储层会产生微小的破裂或滑移会发生并发出微地震信号^[24-25]。由于实际地层中岩石属性为非均质的，为了体现这种非均质性，故将临界孔隙压力设置成一个范围。因此，基质中微地震信号的判定条件为^[14]

$ p \left (z_{i}, t \right)>p_{\text{C}} \left (z_{i} \right) $

(6)

式中：

$z_{i}$—第$i$个网格，$i = 1$，2，$\cdots$，$n$。

综上，可得到水力压裂过程中产生的所有微地震信号，考虑到发出微地震信号的点处发生了破裂或滑移，水力压裂模拟的SRV为所有产生微地震信号的网格体积之和。

2 拟合微地震事件点计算SRV 2.1 三维狄洛尼三角剖分拟合微地震事件点集计算SRV

对于三维空间给定的点集$\boldsymbol{F}$，它的三角剖分是对三维空间中给定的某一包含$\boldsymbol{F}$的一个区域$\Omega$，它将$\Omega$划分为多个四面体，这些四面体除了在边界上(包括顶点、边、面)可以共享外不能有其他交点；这些四面体的顶点都是$\boldsymbol{F}$中的点，且$\boldsymbol{F}$中的每个点都是某一四面体的顶点；每个四面体内的每一个点一定属于区域$\Omega$，区域$\Omega$中每一个点也一定属于某一个四面体。如果没有指定区域，一般取点集$\boldsymbol{F}$的凸包作为这个区域$\Omega$。

三维凸包的基本原理为：首先构造一个初始凸包，然后逐步添加点，形成新的凸包。每次添加点时需要判断该点是否位于已知凸包内，如果是，则放弃该点；如果在凸包外，则计算点到凸包的相切锥面，构造新的凸包。循环上述过程，直到凸包外没有点为止。

目前，已有诸多学者开始采用隐式曲面法构造三角剖分^[26-27]。尽管如此，作为三角剖分的基础，Delaunay法仍然集成于主流算法中^[28]。

三维点集的Delaunay三角剖分是其Voronoi图的对偶图^[29]。对于给定三维空间中任意$n$个位置互异的离散点集合$\boldsymbol{F} = \left \{ f_{1}, f_{2}, \cdots, f_{n} \right \}$，每个站点$f_{i} \left (1 \leqslant i \leqslant n \right)$的Voronoi区域为

$ R_{\text{V}} \left (f_{i} \right) = \left \{q \in \boldsymbol{R}^{3}, \left \vert q f_{i} \right \vert \leqslant \left \vert q f_{j} \right \vert, \forall j \not = i \right \} $

(7)

式中：

$\left \vert q f_{i} \right \vert$—点$q$和$f_i$之间的欧氏距离；

$\left \vert q f_{j} \right \vert$—点$q$和$f_j$之间的欧氏距离；

$R_{\text{V}} \left (f_{i} \right)$—站点$f_{i}$的Voronoi区域，是一个凸多面体。

$\boldsymbol{F}$的Voronoi图即为所有站点的Voronoi区域的并集。与BIN方法相比，三维Delaunay三角剖分提供了更详细、更定量化的SRV几何结构，并且可以通过求其内部每个四面体的体积得到SRV区域的总体积^[11]。

Delaunay三角剖分算法主要分为3类：分治算法、逐点插入算法和三角网生长法^[30]。本文采用逐点插入算法，该方法的基本步骤如下。

(1) 建立包含微地震事件点集的初始四面体；

(2) 对点集按横坐标值大小进行排序，依次插入点集；

(3) 寻找外接球包含当前插入点的所有四面体，并将这些四面体删除形成空腔；

(4) 连接空腔边界面上三角形与当前插入点，形成新的四面体；

(5) 重复上述步骤，依次插入所有的微地震事件点。

四面体集合中的每一个四面体，若有顶点跟超四面体重合则删除该四面体。将基于上述算法所求得的所有四面体的体积相加求和，即可求得改造体积SRV。

2.2 最小体积覆盖椭球拟合微地震事件点计算SRV

通过将微地震事件点拟合成椭球体，可以获得具有更详细几何特征的SRV结构。对于椭球拟合算法，常用的是建立椭球曲面的一般方程。

椭球曲面的一般方程为

$ a'x^{2} + b'y^{2} + c'z^{2} + 2d'xy + 2e'xz + 2f'yz + 2p'x + 2q'y + 2r'z = 1 $

(8)

式中：

$a'$，$b'$，$c'$，$d'$，$e'$，$f'$，$p'$，$q'$，$r'$—控制椭圆形态的9个参数。

故而求解上述方程的充分条件为至少需要9个点，可应用最小二乘法、高斯消除等进行参数求解。

该方法稳定性较差，无法控制椭球的形状，甚至生成不为椭球的其他形状，且无法控制包含的待拟合点数。针对以上问题，本文基于最小体积覆盖椭球(Minimum volume enclosing ellipsoid，MVEE)概念，即采用具有最小体积的椭球覆盖微地震点集中的所有点^[31]，通过构建椭球中心表达式，应用上升法优化椭球拟合。该方法不仅稳定性好，而且可以得到覆盖微地震点集的最小椭球，得到更为合理精确的SRV。

求解最小体积覆盖椭球(MVEE)原理如下^[32]。

假设微地震点集为$\boldsymbol{S} = \left \{x_{1}, x_{2}, \cdots, x_{m} \right\} \in \boldsymbol{R}^{n}$，定义MVEE($\boldsymbol{S}$)为包含集合$\boldsymbol{S}$的最小椭球体积。为了确保获得的包含微地震点集的椭球体积为正，需要假设$m$个点的线性包张成$\boldsymbol{R}^{n}$空间。

因此，椭球的中心形式的表达式可表示为

$ \varepsilon = \left \{\boldsymbol{x} \in \boldsymbol{R}^{n} \vert \left (\boldsymbol{x}-\boldsymbol{c}\right)^{\text{T}}\boldsymbol{E} \left (\boldsymbol{x}-\boldsymbol{c}\right)\leqslant 1\right \} $

(9)

式中：

$\boldsymbol{c} \in \boldsymbol{R}^{n}$—椭球的中心；

$\boldsymbol{E} \in \boldsymbol{S}_{++}^{n}$—对称正定矩阵，确定椭球的形状。

椭球的体积表达式为

$ V\left (\epsilon \right) = \dfrac{v_{0}}{\sqrt{\det \boldsymbol{E} }} = v_{0}\sqrt{\det \boldsymbol{E}^{-1}} $

(10)

式中：

$v_{0}$—$n$维空间中的单位超球体的体积。

从式(10)可以看出，求取包含微地震点集$\boldsymbol{S}$的最小椭球体积问题等价于求取元素$\boldsymbol{c} \in \boldsymbol{R}^{n}$和$n$阶对称正定矩阵，使$\det \boldsymbol{E}^{-1}$最小。

故最小体积覆盖椭球问题可以表示为如下对偶问题

$ \left \{ \begin{array}{l} \max\left [\log \det \boldsymbol{V}\left( u \right)\right ]\\ \text{s.t.} \ \ \ \ \boldsymbol{I}^{\text{T}}\boldsymbol{u} = 1 \end{array} \right. $

(11)

式中：

$\boldsymbol{V} \left (u \right ) = \boldsymbol{Q}\text{diag}(\boldsymbol{u})\boldsymbol{Q}^{\text{T}}$；

$\boldsymbol{u} = \left [u_1, u_2, \cdots, u_m \right]^{\text{T}}$，且$u_i \geqslant 0$，$i = 1$，2，$\cdots$，$m$；

${I^{\text{T}}} = \left[{1, 1, \cdots, 1} \right] \in {\boldsymbol{R}^m}$；

$\boldsymbol{Q} = \left[q_1, q_2, \cdots, q_m \right] = \left[\begin{array}{ cc} \boldsymbol{P} \\ \boldsymbol{I}^{\text{T}} \end{array} \right]\in \boldsymbol{R}^{\left (n + 1 \right)\times m} $；

$\boldsymbol{P} = \left [p_1, p_2, \cdots, p_m \right] \in {\boldsymbol{R}^{n \times m}}$。

上式为凸优化问题，本文采用条件梯度上升法进行求解。

梯度方向$\Delta \boldsymbol{u}$($\Delta \boldsymbol{u} = \boldsymbol{e}_{j}-\bar{\boldsymbol{u}}$，$\boldsymbol{e}_{j}$为第$j$个单位向量，$j = \arg \max\limits_{i}g_{i} \left (\bar{\boldsymbol{u}} \right)$，$g_{i} \left (\boldsymbol{u} \right) = \boldsymbol{q} _{i}^{\text{T}}V \left (\boldsymbol{u} \right)^{-1}\boldsymbol{q} _{i}$，$\bar{\boldsymbol{u}}$是一个满足对偶问题的$\boldsymbol{u}$值)为函数变化最陡的方向，要求函数的最小值，最好的方法就是沿梯度方向按步长$\omega = \dfrac{g_{i} \left (\bar{\boldsymbol{u}} \right)-\left (n + 1\right)}{\left (n + 1\right) \left [g_{i} \left (\bar{\boldsymbol{u}} \right)-1\right]} $进行线性搜索，直至满足收敛精度，该精度可控制椭球覆盖点的范围。

最后，令

$ \boldsymbol{c} ^{*} = \boldsymbol{P} \boldsymbol{u} ^{*} $

由所求对偶问题的最优解$\boldsymbol{u}^{*}$，可得

$ f_{\text{MVEE}} \left (\boldsymbol{S} \right) = \left \{x \in \boldsymbol{R}^{n} \vert \left (\boldsymbol{x}-\boldsymbol{c}^{*} \right)^{\text{T}}\boldsymbol{E}^{*} \left (\boldsymbol{x}-\boldsymbol{c}^{*} \right)\leqslant 1 \right \} $

(12)

式中：

$\boldsymbol{E}^{*} = \dfrac{1}{d} \left [\boldsymbol{P}\boldsymbol{U}^{*}\boldsymbol{P}^{\text{T}}-\boldsymbol{P}\boldsymbol{u}^{*} \left (\boldsymbol{P}\boldsymbol{u}^{*} \right)^{\text{T}} \right]^{-1}$；$d$为该对偶问题的维度，本文取3。

2.3 面元拟合微地震事件点集并计算SRV

油藏储层具有非均质性，故微地震事件点具有一定的离散性。因此，三维狄洛尼三角剖分和最小椭球拟合微地震事件点得到的SRV可能较大。

基于上述原因，本文采用两种联系较紧密的面元算法计算压裂油藏面积(Stimulated reservoir area，SRA)。

第1种面元方法(BIN1)从距离原点最近的事件点开始，通过简单地连接相邻的事件点来形成面元(图 2a)。因此，通常面元会形成梯形。此情况下约束了面元的高度为事件点到$x$轴的实际距离，最小化了面元，是最保守的SRV计算方法^[33]。

与第1种面元方法一样，第2种面元方法(BIN2)同样从距离原点最近的事件点开始，但不是简单地链接相邻面元，而是选择相邻事件点到$x$轴的最大高度创建矩形面元(图 2b)。通过上述方法得到的SRA与$z$方向的平均高度($h$)相乘，即可得到SRV。正如预期的那样，这种方法得到的SRA值比梯形面元法的计算值更大。

具体步骤为：(1)将微地震事件点集投射到$x-y$平面，计算其方位角，并将方位线旋转使其与$x$轴重合，应用上述两种面元算法计算SRA；(2)确定$z$方向的平均高度，将微地震事件点集投影到穿过方位线且与$x-y$平面正交的方位面上，并拟合出一条回归线；(3)在此方位面上求回归线上方和下方的事件点关于回归线的平均高度，上、下方事件点平均高度之和为总厚度($h$)，从而确定SRV。

3 模拟试算 3.1 水力压裂和微地震模拟

建立了均质且各向异性的三维水力压裂模型，模型大小为100 m$\times$20 m$\times$30 m，网格尺寸为1 m$\times$1 m$\times$1 m，射孔点为(50 m，10 m，15 m)，地层的拉梅系数($\lambda$)为9.312$\times{10}^{9} $kg$\cdot$m$^{-1}\cdot$s$^{-2}$，黏度为6.144$\times{10}^{9}$ mPa$\cdot$s，储层渗透率$K_x$ = $K_y$ = 5$K_z$ = 5 mD，孔隙度为0.1，密度2.4 g/cm$^{3}$，压裂液注入量为3.125 m$^{3}$/min。地层初始孔隙压力19 MPa，最大水平地应力($x$方向)25 MPa，最小水平地应力($y$方向)20 MPa，临界应力强度因子$2M$ N$\cdot$m$^{-3/2}$，临界孔隙压力设置为3 000~4 500 psi^[22](1 psi=6.895 kPa)。

模拟得到了水力压裂0.1 d时刻的主裂缝、储层孔隙压力及微地震事件点集分布(图 3~图 5)。

图 3展示了$x-z$方向主裂缝的扩展情况，由于模型设置的储层渗透率为均值，主裂缝也呈较为规则的椭圆。由于裂缝渗透率远大于基质渗透率，流体主要沿裂缝传播，裂缝形态控制着流体压力分布。

由于缝内流体压力传导速率远大于基质，故从图 4和图 5中可以看到，在压力扩散过程中出现$y$方向较窄、$x-z$方向呈高压区且裂缝主要在$x-z$剖面内二维扩展。

图 6展示了0.1 d时刻在$x-y$剖面上的微地震事件点分布。根据临界孔隙压裂理论，从三维孔隙压力分布(图 5)得到微地震事件分布(图 7)。根据传统方法，应用囊括所有破裂点的立方体求解SRV，可得SRV的体积为11 200 m$^{3}$。上述结果包含了过多的未破裂区域，使得SRV计算结果过于乐观，故本文采取了更为精确保守的方法如下：由于模拟得到的破裂点即微地震事件点占据的网格为1 m$^{3}$，可近似认为网格中心即为微地震事件点的坐标，且认为网格单元即为泄流区域，故可将微地震事件点个数做为水力压裂模型得到的SRV，即$V_{\text{SRV}}^{\text{HF}} = 5 484$ m$^{3}$。

3.2 基于微地震事件点拟合SRV

本文基于Eikonal(程函)方程构建目标函数来定位微地震事件。

模型设置为水平层状模型，分为6层，层界面位于2 000，2 040，2 100，2 130，2 180 m处，其中，第2层、第4层、第6层为致密砂岩，第1层、第3层、第5层为泥岩。纵波速度分别设为3 740，4 630，3 750，4 650，3 750，4 650 m/s，横波速度分别设为2 030，2 060，2 040，2 600，2 040，2 600 m/s。

将水力压裂模型得到的微地震事件点坐标采取相同的平移操作，如射孔点(50 m，10 m，15 m)平移为(550 m，510 m，2 115 m)，故2 100~2 130 m层段为压裂层。共有12个检波器，所在的网格位置从(1 000 m，1 000 m，2 050 m)到(1 000 m，1 000 m，2 160 m)，纵向间隔为10 m。

(1) 旅行时计算

首先，建立各向异性介质的二维程函方程^[34]

$ \dfrac{1}{v \left (\theta\right)^{2}} = \left (\dfrac{\partial\tau}{\partial x} \right)^{2} + \left (\dfrac{\partial\tau}{\partial z} \right)^{2} $

(13)

式中：$v \left (\theta \right)$—地震波的相速度，m/s；

$\theta$—地震波传播的相角，(°)；

$\tau$—地震波旅行时，s。

然后，应用有限差分算法求解程函方程，计算初至波旅行时^[35]

$ t_{4} = t_{1} + \sqrt{2 \left (l/v \left (\theta \right)\right)^{2}-\left (t_{3}-t_{2} \right)^{2}} $

(14)

式中：

$t_{1}$，$t_{2}$，$t_{3}$—同一个网格3个顶点的旅行时，s；

$t_{4}$—同一网格未知顶点旅行时，s；

$l$—网格大小，m。

上述推导过程为求解一个网格的旅行时，针对三维介质，建立包含炮点M和一个检波器R的垂直于水平面的切面，对所有检波器做如上操作。根据$t = \sum il/v_{i}$(其中，$i$为网格点，$v_{i}$为网格点的速度)，确定沿炮点垂直与水平方向各网格点的旅行时。应用式(13)进行扫描，可得所有切面网格点的旅行时。

(2) 反演定位微地震事件

基于最小误差原理，构建目标函数如下

$ \begin{array}{ll} O \left (\boldsymbol{X} \right) = &\dfrac{1}{2} \left (\boldsymbol{X} -\boldsymbol{X} ^{\text{p}} \right)^{\text{T}}\boldsymbol{C} _{\text{X}}^{-1} \left (\boldsymbol{X} -\boldsymbol{X} ^{\text{p}} \right) \\ &+ \dfrac{1}{2} \left (\boldsymbol{d} _{\text{obs}}-\boldsymbol{d} _{\text{calc}} \right)^{\text{T}}\boldsymbol{C}_{\text{D}}^{-1} \left (\boldsymbol{d}_{\text{obs}}-\boldsymbol{d}_{\text{calc}} \right) \end{array} $

(15)

式中：

$\boldsymbol{X}$—参数向量，且$\boldsymbol{X} = \left [x, y, z, t_{0}\right]^{\text{T}}$，($x$，$y$，$z$)为事件位置；

$t_{0}$—先验的事件发生时间，s；

$\boldsymbol{X}^{\text{p}}$—先验数据的均值；

$\boldsymbol{C}_{\text{X}}$—先验数据的协方差矩阵；

$\boldsymbol{d}_{\text{obs}}$—观测数据向量；

$\boldsymbol{d}_{\text{calc}}$—计算数据；

$\boldsymbol{C}_{\text{D}}$—目标函数中各数据提供权重系数的测量误差协方差矩阵。

目标函数右式中的的第1项是基于先验信息的修正项。先验事件位置($x$，$y$，$z$)的均值为相应压裂点处射孔位置，先验事件发生时间$t_{0}$是产生裂缝所需要的总时间的一半。协方差矩阵是一个对角矩阵，其对角线上的值与事件位置参数的范围有关。因此，假设水平井沿$x$轴方向，$x$，$y$，$z$方向的允许偏差分别为裂缝宽度、长度、高度的三分之一，事件发生时间的允许偏差为完成该级水力压裂总时间的三分之一。

目标函数右式中的第2项是观测数据$\boldsymbol{d}_{\text{obs}}$与计算数据$\boldsymbol{d}_{\text{calc}}$的匹配项。应用高斯—牛顿法求解上述目标函数的最小值，可求得微地震事件的位置坐标，即定位微地震事件。

为反映初至旅行时的偏差对微地震定位的影响，对正演得到的初至旅行时添加了最大不超过1%，2%，3%，4%，5%的高斯白噪声。通过对加噪后的初至旅行时进行反演，得到在初至旅行时偏差影响下的微地震事件点云。然而，受噪声影响的微地震点集具有较强的离散性，因此，本文以小概率事件为理论基础，鉴定并排除由噪声引起的微地震点云中的异常点，此时，微地震点云中的异常点与平均位置的距离大于计算平均值与3倍的标准偏差之和^[36]。

采用两种面元方法(BIN1、BIN2)、三维狄洛尼三角剖分和最小体积覆盖椭球算法逼近拟合微地震事件点集，得到改造体积$V_{\text{SRV}}^{\text{fitted}}$。图 8至图 12展示了MVEE算法拟合加噪1%至5%和去除异常点情况下得到的微地震事件点云，可以看出去异常点效果明显。这是由于异常点明显孤立于其他微地震事件点集，导致最后形成的改造体积$V_{\text{SRV}}^{\text{fitted}}$包含大片未破裂区域，使得$V_{\text{SRV}}^{\text{fitted}}$偏大，相对误差也相应增大。

表 1展示了基于微地震事件点拟合计算的改造体积($V_{\text{SRV}}^{\text{fitted}}$)与水力压裂SRV($V_{\text{SRV}}^{\text{HF}}$)的相对误差，即${\left (V_{\text{SRV}}^{\text{fitted}}-V_{\text{SRV}}^{\text{HF}} \right)}/{V_{\text{SRV}}^{\text{HF}}} $。从噪声对微地震事件定位的影响来看，相对误差随所加高斯白噪声的增加而增大。其原因是：随着所加随机高斯白噪声的增加，微地震事件点集的离散性随之增强，拟合所得的$V_{\text{SRV}}^{\text{fitted}}$随之增大，与水力压裂的SRV的相对误差也随之加大。

如表 1所示，两种面元算法得到的$V_{\text{SRV}}^{\text{fitted}}$均随着噪声的增大而增大，而在相同加噪情况下，面元算法拟合的$V_{\text{SRV}}^{\text{fitted}}$明显小于MVEE和三维狄洛尼三角剖分的拟合结果，更接近$V_{\text{SRV}}^{\text{HF}}$，相对误差更小。但是，即便是面元算法(BIN1和BIN2)得到的SRV相对误差仍较大，而去异常点后拟合的$V_{\text{SRV}}^{\text{fitted}}$与$V_{\text{SRV}}^{\text{HF}}$的相对误差较小。可以看出这两种面元算法虽然保守，却较稳定，可以为其他拟合算法提供参考。表 1还展示了MVEE和三维狄洛尼三角剖分拟合算法的相对误差随加噪的变化趋势，两种算法得到的$V_{\text{SRV}}^{\text{fitted}}$均随着噪声的增大而增大，加噪后拟合得到的$V_{\text{SRV}}^{\text{fitted}}$与$V_{\text{SRV}}^{\text{HF}}$的相对误差很大，去异常点后拟合的$V_{\text{SRV}}^{\text{fitted}}$与$V_{\text{SRV}}^{\text{HF}}$的相对误差虽明显减小，但仍然较大。可见这两种算法抗噪性不强、稳定性较差，但两种方法不仅从数学领域清楚地划分了微地震活动区域，还提供了更详细、更定量化的SRV几何结构，有助于对压裂效果的评价。

4 结论

(1) 利用水力压裂模型来模拟储层由于主裂缝扩展和压裂液在基质中的漏失，并基于应力强度因子断裂准则(K准则)和临界孔隙压力概念得到实际微地震事件点集，求得了水力压裂模拟的SRV。

(2) 应用面元(BIN1、BIN2)、三维狄洛尼三角剖分和最小体积覆盖椭球算法拟合微地震事件点集。高斯白噪声的加入使得反演定位的微地震事件点集离散性较强，故各算法拟合得到的SRV均与水力压裂模拟得到的SRV相差较大，且随着噪声的增强，相对误差均加大。

(3) 基于小概率事件概念，通过去除微地震事件点集中的异常点，合理降低了微地震点集的离散程度，减小了高斯白噪声对反演定位及拟合计算SRV的影响。

(4) 应用面元、三维狄洛尼三角剖分和最小体积覆盖椭球算法拟合去异常点后的微地震点集，去异常点后拟合得到的SRV值明显变小，与水力压裂模拟得到的SRV相对误差也明显减小。基于BIN1、BIN2两种面元算法得到的SRV与水力压裂的SRV相差不大，表明基于该类方法求解SRV虽然较为保守，但具有较好的稳定性和较强的抗噪性。三维狄洛尼三角剖分和最小体积覆盖椭球两种算法拟合的SRV虽然与水力压裂计算得到的SRV误差相对较大，但这两种方法不仅从数学领域清楚地划分了微地震活动区域，还提供了更详细、更定量化的SRV几何结构。故在实际应用中，可将上述拟合算法中的多种拟合算法结合，有助于压裂效果评价及油气产量预测。