基于Gmsh的起伏地形下井—地直流电法正演模拟

2022-02-26张宇哲孟麟王智

物探与化探 2022年1期

张宇哲，孟麟，王智

(长江大学电子信息学院，湖北荆州 434023)

通讯作者：王智(1985-)，男，湖北省武汉市人，博士，硕士生导师，主要研究方向为电磁法数值模拟和机器学习。Email:1324385898@qq.com

0 引言

Gmsh是一个具有CAD内核和后处理器的开源三维网格生成软件，提供了一种带有参数输入和高级可视化功能的快速、轻便和用户友好的网格剖分工具[1]，它由几何建模、网格剖分、求解器和后处理4个部分组成。自1971年，Coggon[2]首次将有限单元法应用到二维线电流场的正演数值模拟过程后，很多学者对利用有限单元法解决直流电法正演问题进行了深入研究。如：Rijo[3]引入通用性网格，提高了正演的计算速度和精度；Dey等[4]对井—地电法和地表电法的各自装置类型的特点进行了对比分析；Wu X P[5]将不完全Cholesky共轭梯度方法应用在有限单元法的正演模拟中，提升了运算速率；Cardarelli等[6]利用艾米特插值进行傅里叶反变换，提高了正演结果的精度；Tang J T等[7]使用非结构化网格来拟合复杂地形的电阻率模型；Ren Z Y等[8-10]提出自适应有限元算法自动进行网格加密，提高正演的准确性；Pan K J等[11]将外推瀑布式多重网格法应用在2.5D和3D直流电法正演中；Zhang Q J等[12]研究了加密—收缩的网格在有限元模拟中的应用；Kuang X T等[13]对起伏地形条件下长方体磁场无解析解奇点表达式进行了研究；王智等[14]结合井—地观测方法和改进的异常电位法，确定了适用于起伏地表的自然边界条件；武建平等[15]将有限单元法应用在电磁法三维正演模拟中，取得了不错的结果；严波等[16]将对偶加权后验误差估计方法应用在自适应有限元正演模拟中，提高了自适应有限元网格加密的有效性。本文采用不规则网格剖分拟合起伏地形，使用井—地联合观测方法来进行起伏地形下的异常响应的观测，研究使用不同观测装置时山谷地形对下方异常响应的影响。

1 2.5D有限元正演模拟

在三维无限半空间Ω中，点电源A的坐标为(xA,yA,zA)，若要求解直流电法勘探的点源2.5D边值问题，需要通过傅里叶变换将空间域电位u(x,y,z)转化为波数域电位U(x,k,z)。波数域电位U满足的偏微分方程[17]如下：

(1)

式中：σ是地下介质电导率；I是供电电流强度；δ是Delta脉冲函数；ΓS和Γ∞分别为三维无线半空间的地表边界和无穷远边界；n代表边界的外法向向量；r是Γ∞上任意一点到点电源的向量；K1和K0分别为一阶和零阶修正贝塞尔函数；k为波数。该偏微分方程的边界条件分别为：①地表边界采用Newman边界条件；②无穷远边界采用混合边界条件。

用有限单元法解决上述边值问题时，需推导边值问题对应的变分问题，式(1)对应的变分问题为：

(2)

得到对应变分问题后，用三角形单元将求解区域离散化；然后在每个离散单元e上用多项式函数和对应的形函数来近似待求量；最后将每个单元合并，求解大型稀疏线性方程组，得到各个节点处的波数域电位值U，进行傅里叶反变换便可得到空间中的电位值u。本文选用的波数ki与傅里叶反变换系数gi选自文献[17]中采用最优化方法得到的5个波数和对应的傅里叶反变化系数。

2 建模及观测装置

在Gmsh建立模型和网格剖分过程中，有2种建模方式：①采用Gmsh软件左侧的GUI进行交互建模；②在geo格式文件中采用Gmsh自己的脚本语句进行建模。本文选择后一种方式，Gmsh脚本建模原理请参考本文附录A。

井—地电阻率观测装置由井下供电电极和地表测量电极共同组成，由于井中的观测电极距离地下的异常体更近，所以得到的观测数据量更大且更精确，还可以分别移动供电电极和测量电极的相对位置，来观测和研究不同水平位置和不同深度的地下介质情况。根据实际使用的观测装置的不同，井—地电阻率观测装置可以细分为井—地二极观测装置、井—地三极观测装置、井—地偶极观测装置等。在记录和作图时，井—地二极装置(图1a)以井中供电电极A的纵坐标为记录点纵坐标，地表观测电极M的横坐标为记录点横坐标；井—地三极装置(图1b)以井中供电电极A的纵坐标为记录点纵坐标，地表观测电极MN中点的横坐标为记录点横坐标。

图1 观测装置示意[18]Fig.1 Example of observation device[18]

3 数值模拟及分析

地质模型的建立和网格剖分是使用有限单元法解决地球物理正演问题的重要环节。用Gmsh软件对地质模型进行不规则的网格剖分时进行局部加密，能更精确地拟合各种地形情况，提高正演结果的准确性。2.5D有限单元正演数值模拟程序采用CSR格式来压缩储存大型稀疏矩阵；选用Eigen库中的BiCGSTAB(biconjugate gradient stabilized method)求解器解正演过程中的线性方程组。首先，为了验证Gmsh在正演问题中的适用性和有限单元正演数值模拟算法的正确性，选用具有解析解的层状模型进行正演模拟；然后，使用井—地二极和井—地三极电阻率观测装置进行观测，对两种山谷地形下异常体的异常响应进行正演模拟，讨论山谷地形对下方异常体的异常响应的影响。

3.1 算法验证

水平层状模型如图2所示，采用二极装置进行观测，得到的电位解析解和数值解的对比如图3所示。从结果可以看出数值解和解析解的误差主要集中在靠近供电电极的位置，其余位置的误差较小，证明了Gmsh在正演问题中的适用性和有限单元正演数值模拟算法的正确性。

图2 水平层状模型Fig.2 Two-layers geo-electric cross section

图3 电位数值解和解析解对比Fig.3 The curve for analytical solution and calculation solution of potential

3.2 计算实例

所设3个模型如图4所示，分别为平坦地形、底部平缓的山谷地形以及底部尖锐的山谷地形下的异常体模型，井中供电电极A的变换范围为0～20 m，间距1 m；地表观测电极的变换范围为0～50 m，间距1 m;山谷地形底部深度为5 m;异常体顶部埋深为8.5 m，长4 m，高2 m。为了比较各种地形和装置的观测视电阻率异常，围岩电阻率设定为100 Ω·m，低阻体电阻率为5 Ω·m，高阻体电阻率为5 000 Ω·m。在作图时，超出测量点范围的数据由Surfer软件根据测量范围内的数据推导画出。

图4 模型示意Fig.4 Sketch map of model

采用井—地二极装置对模型一进行正演模拟所得视电阻率断面如图5所示。在图5a中，低阻异常体的视电阻率断面中只存在独立的椭圆形低阻异常，是由于低阻异常体对电流的吸引作用而形成；垂直方向上低阻异常中心深度z=8.5 m，与异常体的顶部埋深一致，最低幅值为76 Ω·m，水平方向上低阻异常中心的范围x=23～27 m，与异常体的宽度相同。图5b显示,高阻异常体的视电阻率断面垂直方向上存在对称的高阻异常和低阻异常，当源点在异常体上方时，由于高阻体对电流的排斥作用，异常体上方地表的电流密度大于正常电流密度，因此呈现出大于背景值的高阻异常，当源点在异常体下方时，由于高阻体的屏蔽电流作用，异常体上方地表的电流密度小于正常电流密度，因此呈现小于背景值的低阻异常；上部高阻异常的最高幅值为120 Ω·m，下部低阻异常的最低幅值为84 Ω·m，异常分界面为z=9.5 m，和异常体中心深度相同;在水平方向上，高阻、低阻异常的分布范围都为x=23～27 m，与异常体的宽度相同。两种情况下均可通过视电阻率断面图中的异常响应准确推断异常体位置。

图5 模型一的视电阻率断面(井—地二极装置)Fig.5 Apparent resistivity cross-section view of model 1(borehole-ground pole-pole device)

图6为模型二对应的纯山谷地形的视电阻率断面。采用井—地二极装置观测时，在水平方向上视电阻率等值线最为密集的地方在山谷地形和平坦地形的分界处(图6a)，山谷地形最低处下方存在一个相对左右两侧的高阻异常，异常幅值为95 Ω·m；在垂直方向上，相对左右两侧的高阻异常的顶部和山谷地形的底部的深度平齐。采用井—地三极装置，在水平方向从左到右依次为低阻—高阻—低阻趋势(图6b)，两侧的低阻脉冲异常位于山谷地形和平坦地形的分界处，最低异常幅值为46 Ω·m；高阻异常位于山谷地形最低处下方，最高异常幅值为106 Ω·m；在垂直方向上，两侧低阻脉冲异常最低幅值等值线顶部高度和山谷地形的底部的深度相等。

图6 模型二的视电阻率断面Fig.6 Apparent resistivity cross-section view of model 2

模型二正演结果的视电阻率断面如图7所示。采用井—地二极观测装置，低阻异常体的视电阻率断面图(图7a)上存在一个独立的椭圆形低阻异常，和平坦地形下的低阻异常相比，该异常的范围明显缩小，其中心深度位于z=8.5 m，与异常体的顶部埋深一致，最低幅值为54 Ω·m。山谷地形底部存在相对左右两侧的高阻异常，异常中心深度和山谷地形底部深度一致，幅值为83 Ω·m；低阻异常中心的范围为x=23～27 m，与异常体的宽度相同。高阻异常体的视电阻率断面图(图7b)中高阻、低阻并存，和平坦地形下的异常响应相比，上部分高阻异常的范围有所缩小;高阻异常和低阻异常不完全对称,上部高阻异常的最高幅值为128 Ω·m，下部低阻异常的最低幅值为66 Ω·m,异常分界线为z=9.5 m，和异常体中心深度相同;在水平方向上，高阻、低阻异常的分布范围都在x=23～27 m中，与异常体的宽度相同。

山谷地形下异常体模型断面图中，低阻异常体所产生的低阻异常区域基本和低阻异常体模型重合，高阻异常区域中心位于山谷底部最深处；高阻异常体产生的高阻异常区域和低阻异常区域基本对称，分界线和异常体中心埋深一致，据此可以判断山谷地形下方为低阻异常体还是高阻异常体。

图7c、d为采用井—地三极观测装置所得。图7c中可见，低阻异常体的视电阻率断面存在山谷地形下方的高阻异常和独立的椭圆形低阻异常，高阻异常中心位于山谷地形的最低处，最高幅值为105 Ω·m；低阻异常中心深度z=9.5 m，与异常体中心的埋深平齐，高、低阻异常分界线为z=8.5 m，与异常体的顶部埋深一致；在水平方向上，低阻异常中心的范围与异常体的宽度相同，位于山谷地形和平坦地形的分界处的低阻脉冲异常依然明显存在。图7d显示高阻异常体的视电阻率断面中同样存在高阻和低阻异常，高阻异常中心深度z与异常体的顶部埋深一致，其电阻率最高幅值为180 Ω·m，低阻异常最低为10 Ω·m，高、低阻异常的分界线为z=10.5 m，与异常体的底部深度相同；在水平方向上，两侧的低阻脉冲异常位于山谷地形和平坦地形的分界处，高、低阻异常中心的范围均为x=23～27 m，与异常体的宽度相同。

图7 模型二的正演结果Fig.7 Apparent resistivity cross-section view of the forward result of model 2

图8为模型三对应的纯山谷地形的视电阻率断面。采用井—地二极装置观测时，视电阻率断面存在方向向下的低阻脉冲异常(图8a),从左到右依次为高阻—低阻—高阻趋势，等值线最密集的位置在山谷地形和平坦地形的分界处;低阻异常位于山谷地形下方，最低幅值为76 Ω·m，向下的低阻脉冲异常的底部和山谷地形的底部的深度平齐。采用井—地三极装置观测时，视电阻率断面与井—地二极装置观测得到的相反，存在方向向上的脉冲异常(图8b)，异常范围x=10～40 m，与山谷地形范围一致;低阻异常位于山谷地形正下方，其顶部和山谷地形的底部的深度平齐，低阻异常最低幅值为25 Ω·m。

模型三正演结果的视电阻率断面如图9所示。图9a、b为采用井—地二极观测装置所得。低阻异常体的视电阻率断面(图9a)中存在一个独立的椭圆形低阻异常,异常中心深度与异常体的顶部埋深一致，位置与异常体中心水平位置相同，最低幅值为48 Ω·m。高阻异常体的视电阻率断面(图9b)中高阻、低阻异常并存，高阻异常的最高幅值为108 Ω·m，低阻异常的最低幅值为60 Ω·m，异常分界线和异常体中心深度相同；高阻异常的中心由于山谷地形的影响被阻断，低阻异常的中心位于x=25 m处，与异常体的水平中心位置平齐。

图8 模型三的视电阻率断面Fig.8 Apparent resistivity cross-section view of model 3

图9 模型三的正演结果Fig.9 Apparent resistivity cross-section view of the forward result of model 3

图9c、d 为采用井—地三极观测装置所得。在低阻异常体的视电阻率断面(图9c)中，垂直方向上有存在一个椭圆形低阻异常和一个脉冲异常。在垂直方向上，椭圆形低阻异常中心位于z=9.5 m，与异常体中心埋深相同；纯山谷地形显示的低阻脉冲异常顶部应在z=5 m处，由于存在低阻异常体的影响，低阻脉冲异常的顶部位置下降到z=12 m处；在水平方向上，低阻异常和低阻脉冲异常中心均为x=25 m，最低幅值为15 Ω·m，低阻脉冲异常的分布范围仍与山谷地形范围一致。高阻异常体的视电阻率断面(图9d)中高阻、低阻异常并存，在垂直方向上，上面的高阻异常中心深度z=8.5 m，与异常体的顶部埋深一致，其最高幅值为115 Ω·m，下面是低阻异常最低幅值为10 Ω·m，高、低阻异常的分界线z=9.5 m，与异常体的中心深度相同；在水平方向上，高阻异常在x=25 m被阻断，下面的低阻脉冲异常形态相对于纯山谷地形的的异常响应情况来说有所变化，但影响范围仍然在x=10～40 m，与山谷地形范围一致。