0 引言

煤层气是世界公认的三大非常规天然气之一，储层的压裂改造是提高采收率的主要手段^[1]。其中，微地震监测技术可以通过储层压裂改造引发的震源点进行成像，从而为储层压裂改造提供依据。然而，煤层气水力压裂区围岩或煤层破碎严重、煤层较软，致使微地震信号微弱，资料解释困难^[2]。由于电磁信号对导电流体和裂缝中的压裂剂较敏感，Beskardes等^[3]采用三维直流电法对煤层气压裂效果进行了深入研究，其中高精度三维电法数值模拟是关键。

目前，三维直流电正、反演计算中应用最广的数值模拟方法主要有积分方程法^[4-6]、有限差分法^[7-8]和有限元法^[9-13]。这三种方法各有优缺点：积分方程法计算速度快，但只能应用结构化网格处理异常体，对复杂地质体的模拟具有局限性；有限差分法精度不高，且只能应用结构化网格；有限元法精度较高，且可以结合非结构化网格^[14-17]处理复杂模型，但相较于前两种方法，所需的计算内存较高。近些年来，随着计算机硬件计算能力的提高和数值计算方法的进步，有限元法得到广泛应用^[18]。

在三维直流电数值模拟中，有限元法通常可通过两种策略提高解的精度，即加密网格和应用高阶形函数。为快速得到最优的网格分布，目前正演通常采用h型自适应加密算法^[19-22]对网格进行加密。该算法首先固定网格单元上形函数的阶数(通常较低)，然后利用后验误差估计估算每个单元的误差，最后对误差较大的单元进行加密。但是，若单元上解的光滑度较高，对单元进行加密仍无法快速提高解的精度。为解决这一问题，Grayver等^[23]在求解三维地电模型时，对网格中的所有单元应用高阶形函数，提高了解的精度，并证明高阶自适应有限元法可以用最少的自由度得到最精确的有限元解。

本文基于高阶自适应有限元实现了直流电阻率模型的正演。首先阐述了有限元算法的实现步骤和关键技术，包括高阶形函数生成方法、后验误差估计技术和悬挂点上的约束条件；然后通过数值算例验证了算法的精确性，并且对比了形函数阶数不同时(p=1、2、3)的模拟结果；最后，将此方法应用于沁水盆地南部某煤气层压裂监测区的三维模型，通过三维直流电阻率模拟，验证了方法的有效性。

1 正演理论 1.1 控制方程

假设在地面A点放置一个电流强度为I的点电源(图 1)，空间电位分布满足

$ \nabla \cdot (\sigma {\kern 1pt} \nabla U) = - I{\delta _A} $

(1)

式中：σ表示电导率；δ_A表示A点的Dirac delta函数；U表示电位。在地表Γ_s处，电流沿地表流动，所以电流密度的法向分量为0，表示为

$ {\left. {\frac{{\partial U}}{{\partial \mathit{\boldsymbol{n}}}}} \right|_{{\varGamma _{\rm{s}}}}} = 0 $

(2)

式中n是边界上外法向向量。

假定区域Ω内，电性不均匀性对无穷远边界Γ_∞(图 1)处的电位分布不产生影响，即Γ_∞处的电位为点电源所产生的电位，表示为

$ {\left. U \right|_{{\varGamma _\infty }}} = \frac{c}{r} $

(3)

式中：c是比例系数；r是A点到边界Γ_∞的距离。将式(3)对n求导，消去常数c，可得

$ \frac{{\partial U}}{{\partial \mathit{\boldsymbol{n}}}} + \frac{{\cos (\mathit{\boldsymbol{r}},\mathit{\boldsymbol{n}})}}{r}U = 0 $

(4)

式中r是由点A指向边界Γ_∞的方向向量。

综上所述，三维直流电阻率模型边界值问题可由以下偏微分方程描述

$ \left\{ {\begin{array}{*{20}{l}} {\nabla \cdot (\sigma \nabla U) = - I{\delta _A}}&{ \in \varOmega }\\ {\frac{{\partial U}}{{\partial \mathit{\boldsymbol{n}}}} = 0}&{ \in {\varGamma _{\rm{s}}}}\\ {\frac{{\partial U}}{{\partial \mathit{\boldsymbol{n}}}} + \frac{{\cos (\mathit{\boldsymbol{r}},\mathit{\boldsymbol{n}})}}{r}U = 0}&{ \in {\varGamma _\infty }} \end{array}} \right. $

(5)

1.2 有限元离散

将区域Ω离散为一系列非结构化六面体单元，每个单元e上电位的近似解为

$ {U_{\rm{e}}} = \sum\limits_{i = 1}^m {{u_{\rm{i}}}} {\phi _{\rm{i}}} = \mathit{\boldsymbol{u}}_{\rm{e}}^{\rm{T}}\mathit{\boldsymbol{\phi}} = {\mathit{\boldsymbol{\phi}} ^{\rm{T}}}{\mathit{\boldsymbol{u}}_{\rm{e}}} $

(6)

式中：m是单元上自由度个数；u_i是单元中第i个自由度上的电位；u_e是由m个自由度上未知电位组成的向量；ϕ_i是第i个自由度上的形函数；ϕ是由m个形函数组成的向量。

对式(5)应用伽辽金有限元法^[24]，可得单元e上的弱解形式为

$ \begin{array}{*{20}{c}} {\int_{{\varOmega ^{\rm{e}}}} \sigma \nabla {\phi _j} \cdot \nabla {U_{\rm{e}}}{\rm{d}}{\varOmega ^{\rm{e}}} + \int_{\varGamma _\infty ^{\rm{e}}} {\frac{{\cos (\mathit{\boldsymbol{r}},\mathit{\boldsymbol{n}})}}{r}} \cdot {\phi _j}{U_{\rm{e}}}{\rm{d}}\varGamma _\infty ^{\rm{e}} - }\\ {\int_{{\varOmega ^{\rm{e}}}} {{\phi _j}} \cdot I{\delta _A}{\rm{d}}{\varOmega ^{\rm{e}}} = 0\quad j = 1,2 \cdots ,m} \end{array} $

(7)

式中：Ω^e是单元e所围区域；Γ_∞^e是无穷远边界上单元e的面。据式(7)对每个单元积分，可得

$ {\mathit{\boldsymbol{k}}_{\rm{e}}}{\mathit{\boldsymbol{u}}_{\rm{e}}} = {\mathit{\boldsymbol{f}}_{\rm{e}}} $

(8)

式中：k_e是m×m阶矩阵；f_e是m维向量。矩阵k_e中第i行、第j列元素k_ij^e和向量f_e的第i行元素fi^e分别为

$ k_{ij}^{\rm{e}} = \sigma \int_{{\varOmega ^{\rm{e}}}} {\nabla {\phi _i} \cdot } \nabla {\phi _j}{\rm{d}}{\varOmega ^{\rm{e}}} + \frac{{\sigma \cos (\mathit{\boldsymbol{r}},\mathit{\boldsymbol{n}})}}{r}\int_{\varGamma _\infty ^{\rm{e}}} {{\phi _i}{\phi _j}{\rm{d}}\varGamma _\infty ^{\rm{e}}} $

(9)

$ f_i^{\rm{e}} = \int_{{\varOmega ^{\rm{e}}}} {{\phi _i}I{\delta _A}{\rm{d}}{\varOmega ^{\rm{e}}}} $

(10)

对所有单元分别组装矩阵，可以得到系统线性方程组

$ \mathit{\boldsymbol{KU = F}} $

(11)

式中：K是一个n×n阶的稀疏对称正定矩阵，n是网格中自由度总数；U是网格中所有自由度上未知电位组成的n维向量；F是一个n维向量，描述源的分布。

1.3 形函数的生成

任意空间中、任意阶数的形函数可由一维多项式的张量积表示^[25]

$ \phi _i^{(p)}(x) = \prod\limits_{0 \le d < {\rm{dim}}} {\theta _{{j_d}(i)}^{(p)}({x_d})} $

(12)

式中：θ_j^(p)(·)是一维多项式，p表示多项式的阶数；x_d表示空间中心点的坐标，下标“d”表示空间维度序号；dim表示维数；j_d(i)表示将d维中按字典排序的形函数的下标映射到一维形函数中的下标。本文取dim=3，有

$ \left\{ {\begin{array}{*{20}{l}} {{j_0}(i) = \left\lfloor {i/(p + 1)} \right\rfloor }\\ {{j_1}(i) = i \times {\rm{mod}}(p + 1)}\\ {{j_2}(i) = \left\lfloor {i/{{(p + 1)}^2}} \right\rfloor } \end{array}} \right. $

(13)

式中：$\left\lfloor \cdot \right\rfloor $为向下取整；mod(·)表示取模。多项式θ_j^(p)(·)可以根据多项式阶数p动态计算得到

$ \theta _j^{(p)}\left( {\frac{n}{p}} \right) = {\delta _{nj}};\quad 0 \le n \le p $

(14)

单元上形函数的个数m可由空间维度d和形函数的阶数p确定

$ m = {(p + 1)^d} $

(15)

以二维形函数(p=2)在单元上的分布为例(图 2)，这些自由度上形函数的数学表达式为

$ \left\{ \begin{array}{l} \begin{array}{*{20}{l}} {{\phi _0} = (2x - 1)(x - 1)(2y - 1)(y - 1)}\\ {{\phi _1} = x(2x - 1)(2y - 1)(y - 1)} \end{array}\\ \begin{array}{*{20}{l}} {{\phi _2} = (2x - 1)(x - 1)y(2y - 1)}\\ {{\phi _3} = x(2x - 1)y(2y - 1)}\\ {{\phi _4} = (2x - 1)(x - 1) \times 4(1 - y)y} \end{array}\\ \begin{array}{*{20}{l}} {{\phi _5} = 4x(2x - 1)(1 - y)y}\\ {{\phi _6} = 4(1 - x)x \times (2y - 1)(y - 1)}\\ {{\phi _7} = 4(1 - x)x \times y(2y - 1)}\\ {{\phi _8} = 4(1 - x)x \times 4(1 - y)y} \end{array} \end{array} \right. $

(16)

根据式(16)，分别绘制自由度0和自由度7上形函数在单元上的分布形态(图 3)。图中黑色区域表示z=0平面，所以形函数为负值的区域在图 3中显示不完整。可以看出，形函数所在自由度上函数值为1，在其他自由度上函数值均为0。

1.4 后验误差估计方法

后验误差估计可以无需人工干预，以最快的速度得到最优的网格分布。目前主要有两类后验误差估计方法，即基于梯度恢复^[15]和基于残差^[26]的后验误差估计方法。本文采用前者将单元中每个面上有限元解的梯度跳跃量^[27]作为单元上的基本误差

$ {E_{{F_w}}} = {c_{{F_w}}}\int_{\partial {K_w}} {{{\left[ {\frac{{\partial {u_{hp}}}}{{\partial \mathit{\boldsymbol{n}}}}} \right]}^2}} $

(17)

式中：w表示单元面的编号；E_F表示一个面上的误差估计量；c_F=h_e/24，h_e是单元e中最长的对角线长度；符号[·]代表在单元面上的跳跃值。单元e上的误差估计量η_e²可由公式

$ \eta _{\rm{e}}^2 = \sum\limits_{w = 1}^l {{E_{{F_i}}}} $

(18)

计算，其中l是单元中面数目。

得到所有单元的误差估计后，对误差前10%的单元进行加密得到新的网格，然后在新的网格上重新计算有限元解，直到迭代次数或自由度个数达到设定值，迭代结束。h型自适应有限元算法流程如图 4所示。

1.5 悬挂点上的约束条件

本文采用八叉树加密方法^[28](图 5)对网格进行自适应加密，所以当一个单元一个面被加密时，在相邻两个单元的公共面上会出现悬挂点。为保证有限元的连续性，需要在这些悬挂点上添加一些约束条件。为方便起见，以二维网格上应用2阶形函数为例说明悬挂点上所施加的约束条件(图 6)。三维空间及其不同阶数形函数的情况与此类似。

图 6中，自由度0、1和2属于单元M₁，自由度3、4和5属于单元M₂。为保证单元M₁与单元M₂的公共面上有限元的连续性，需满足

$ {u_0}{\phi _0} + {u_1}{\phi _1} + {u_2}{\phi _2} = {u_3}{\phi _3} + {u_4}{\phi _4} + {u_5}{\phi _5} $

(19)

其中自由度2与自由度3、自由度1与自由度5表示同一自由度，即

$ \left\{ {\begin{array}{*{20}{l}} {{u_2} = {u_3}}\\ {{u_1} = {u_5}} \end{array}} \right. $

(20)

当通过1.3节中的步骤求得自由度0、1和2上的形函数表达式后，可以得到

$ {u_4} = - \frac{1}{8}{u_0} + \frac{3}{4}{u_1} + \frac{3}{8}{u_2} $

(21)

上式即为悬挂点4上的约束条件。

2 数值算例

本文在程序实现中使用了开源有限元框架deal.II^[29-30]。所使用的计算平台配备了Intel Xeon E5 2680 V4 CPU，包含14个处理器核心，主频为2.4GHz，内存为128GB。

2.1 正确性验证

采用均匀半空间模型验证形函数阶数为3时算法的正确性。模型的计算区域为200m×200m×200m，点源位于原点O(0, 0, 0)，沿x轴正方向以2m的间隔布置20个测点。

由于此模型较为简单，电位解析解^[31]计算公式为

$ u = \frac{I}{{2{\rm{ \mathit{ π} }}\sigma R}} $

(22)

式中R为测点与点源的距离。

通过有限元模拟各测点的电位，再计算得到视电阻率曲线(图 7a)；通过与解析解对比，得到所有测点的相对误差(图 7b)。从图 7a可以看出：距离点源越远，有限元的解越精确；随着迭代次数的增加，有限元解快速收敛于解析解。从图 7b可以看到，随着自由度个数的增加，所有测点的相对误差大幅度降低。表 1给出了第1次、第3次以及第5次网格剖分方案下的自由度个数、计算时间和测点的最大相对误差，最大相对误差从124.8%降至0.3%，证明本文算法具有很高的精度。

2.2 收敛速度

应用均匀半空间模型分析形函数阶数p=1、2、3时自适应有限元解的收敛速度。计算区域为500m×500m×500m，背景电阻率为100Ω·m，初始剖分网格为8000个单元。从源点S(0, 0, 0)处向地下注入电流强度为1A的电流，沿x轴正方向100~400m范围内等间隔布置31个测点，间距为10m。对不同的网格剖分方案分别计算各测点的有限元解，然后参照解析解计算所有测点的相对误差，最终用所有测点的视电阻率平均相对误差展示有限元解的收敛速度，结果如图 8所示。表 2列出了形函数阶数p=1、2、3时最终网格上的自由度个数和所有测点的电阻率有限元解的平均相对误差。

从图 8和表 2可以看出，对于粗糙网格，形函数的阶数越高，有限元解的精度越高；有限元解的相对误差随着迭代次数的增加逐渐下降，p=1时有限元解曲线下降最慢，且最终网格上自由度个数最多为8256929，而所有测点的有限元解的平均相对误差仅降至1.5×10^-4；p=2时有限元解的收敛速度有所提升，最终网格上的自由度个数为3277344，大约是p=1时最终网格上自由度个数的5/8，所有测点有限元解的平均相对误差降至1.4×10^-5，即有限元解的精确度相较于p=1时提高了约11倍；p=3时有限元解的收敛性能最佳，最终网格上自由度个数为1968695，所有测点的有限元解的平均相对误差为8.8×10^-6，相较于p=2时，自由度个数减少了1308649，精度提高了约1.5倍。综上所述，p=3时有限元程序可以用最少的自由度个数得到最精确的解。

图 9展示了形函数阶数p=1、2、3时最终网格剖分数目和误差分布。众所周知，在点源附近电势变化剧烈，导致解的光滑度较低，在距离点源较远的区域中解的光滑度较高。从图中可以看到，在距离点源较远、解比较光滑的区域中，形函数的阶数越高，最终网格上有限元解的误差越小；形函数的阶数p越高，最终整体误差越小。在距点源非常近的区域，p=1、2时网格得到了充分加密，而p=3时网格未得到充分加密。由于点源附近网格加密程度不同，从图 9c(右)可以清楚看到，p=3时在距离点源非常近的区域，有限元解的误差较高；而p=1、2时点源附近的误差较p=3时明显降低。

2.3 异常体模型

对低阻异常体模型(图 10)进行模拟，以验证算法的有效性。异常体是一个电阻率为10Ω·m的低阻长方体，位于坐标原点正下方。背景电阻率为100Ω·m。

采用对称四极装置模拟计算沿x轴的视电阻率拟断面。x方向的计算范围为-400~400m，测点间距为20m；z方向为100~600m。供电电极的极距由200m逐渐增加到1200m，测量电极的极距为相应供电电极距的1/3。

首先，计算得到沿x轴的视电阻率拟断面(图 11a)，可以清楚地看到，在x轴-100~100m范围内有一个低阻异常区，其宽度与模型异常体基本一致。

然后，在地表平行于x轴布置41条测线，测线范围为y=-200~200m，测线间距为10m，沿x方向测点的布设同图 11a。供电电极极距仍为400m。截取视深度z=200m的数据绘制视电阻率平面图(图 11b)。可以看到，大致由x轴-100~100m、y轴-50~50m所围区域是一个低阻区，其中心位置与模型异常体的中心位置吻合，覆盖范围与模型也基本一致。

2.4 合成数据模拟

为验证高阶自适应有限元算法对较复杂模型的模拟效果，本文基于沁水盆地南部某煤层气压裂监测资料，建立三维电性模型验证本文算法的有效性。

沁水盆地南部煤层埋藏浅、厚度大、煤储层中煤层气含量高，且具有相对较高的渗透率，压裂后富水地层较上覆和下伏地层通常表现为低阻特征^[32-33]。模型的计算区域为2000m×2000m×1000m。煤层气异常体在x、y和z轴方向的范围分别为-76~92m、-50~50m、124~190m，异常体示意图如图 12所示。模型的背景电阻率为100Ω·m，压裂前异常体电阻率设为10000Ω·m，压裂后为1Ω·m。

采用对称四极电阻率测量装置。测线沿x方向布设，测线范围y为-80~80m，测线间隔5m，测点范围x为-100~100m，间隔为5m。供电电极的极距为300m。分别计算煤层压裂前、后的视电阻率分布，并从计算结果中提取x方向-100~100m、y方向-80~80m、z=150m的视电阻率数据。通过定量计算得到压裂后的视电阻率相对异常

$ R = \frac{{C - B}}{B} $

(23)

式中B、C分别为压裂前、后的视电阻率。图 13a、图 13b分别为第一次和压裂结束后的视电阻率相对异常分布。

第一次压裂后，x方向-74~-40m范围内充满压裂剂，该区域电阻率约1Ω·m(即压裂剂的电阻率)，其他范围内异常体的电阻率值仍很高，约为10000Ω·m。计算得到对称四级装置条件下z=150m平面的视电阻率后，通过定量计算得到了视电阻率相对异常图(图 13左)。因为压裂范围较小，所以图中左侧区域显示出小范围的低阻异常(图 13中蓝色的低阻区域)。

压裂结束后，异常区域内全部充满压裂剂，这些区域视电阻率主要体现为压裂液的低阻(1Ω·m)特征(图 13右)。图中的蓝色低异常区域与模型(图 12)中的煤层气分布区域基本一致。

分析图 13可知，采用本文算法并结合对称四级电阻率测量装置，可以定性地反演异常体的水平分布，计算结果与模型相吻合。因此，该方法在煤层气压裂监测领域具有重要应用价值。

3 结论

本文应用高阶自适应有限元算法实现了直流电阻率模型正演，算例表明该算法具有较高的精度。

(1) 形函数阶数p=3时，数值模拟有限元解收敛最快，可以用少量的自由度个数得到精确的有限元解。

(2) 在点源附近电势剧烈变化的区域，可通过加密网格提高解的精度；在离点源较远处，解析解比较光滑，可通过提高单元上形函数的阶数提高解的精度，效果显著。

(3) 合成数据模型的模拟结果证明三维直流电阻率法是煤层气储层压裂监测的一项重要补充技术。

感谢河南理工大学对本项目提供了高性能计算资源！