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,9,10]提出自适应有限元算法自动进行网格加密,提高正演的准确性;Pan K J等^[11]将外推瀑布式多重网格法应用在2.5D和3D直流电法正演中;Zhang Q J等^[12]研究了加密—收缩的网格在有限元模拟中的应用;Kuang X T等^[13]对起伏地形条件下长方体磁场无解析解奇点表达式进行了研究;王智等^[14]结合井—地观测方法和改进的异常电位法,确定了适用于起伏地表的自然边界条件;武建平等^[15]将有限单元法应用在电磁法三维正演模拟中,取得了不错的结果;严波等^[16]将对偶加权后验误差估计方法应用在自适应有限元正演模拟中,提高了自适应有限元网格加密的有效性。本文采用不规则网格剖分拟合起伏地形,使用井—地联合观测方法来进行起伏地形下的异常响应的观测,研究使用不同观测装置时山谷地形对下方异常响应的影响。

1 2.5D有限元正演模拟

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

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

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

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

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

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

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

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

所设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软件根据测量范围内的数据推导画出。

采用井—地二极装置对模型一进行正演模拟所得视电阻率断面如图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,与异常体的宽度相同。两种情况下均可通过视电阻率断面图中的异常响应准确推断异常体位置。

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

模型二正演结果的视电阻率断面如图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,与异常体的宽度相同。

图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处,与异常体的水平中心位置平齐。

图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,与山谷地形范围一致。

采用Gmsh软件进行建模和网格剖分,使用井—地二极和井—地三极电阻率装置进行观测,用有限单元法对2种山谷地形下异常体的异常响应进行正演模拟后,对正演模拟结果进行分析,得到了以下结论:

1)底部尖锐的纯山谷地形比底部平缓的纯山谷地形产生的异常响应更加明显和突出,对下方异常体的异常响应的影响也更大。尖锐山谷地形产生的影响主要集中在尖锐处水平2 m内,且异常幅度大,平缓山谷地形的影响范围较大而异常幅度小。

2)井—地三级装置观测结果比井—地二级装置观测结果的异常响应幅值更大,异常形态更复杂。就水平分辨率来说,井—地三极装置的分辨率比井—地二极更强;2种装置的垂向分辨率相差不大,但是异常形态各有特点,在实际工程中应注意区分。

3)采用不规则网格剖分拟合起伏地形和使用井—地联合观测方法来进行起伏地形下的地质情况勘探能得到较好的结果,同时也证明有限元软件Gmsh在地球物理有限单元法正演建模和网格剖分方面有良好的应用价值。

点、线、围线和面在geo格式文件中的设置语句依次为:

Point(P_number)={x,y,z,expression};

Line(L_number)={P_number, P_number};

Curve Loop(CV_number)={L_number,…,L_number};

Plane Surface(PS_number)={CV_number,…, CV_number};

其中:点设置语句中P_number是点的编号, x、y、z代表点在三维坐标系中的位置,expression用来设置点周围剖分网格的大小;线设置语句中L_number是线的编号,大括号内是构成线段的起点和终点的编号;围线由多条线段连接闭合而成,围线设置语句中CV_number是围线的编号,大括号内是组成围线的线段的编号;平面由单条(或多条)围线组成,面设置语句中PS_number内代表面的编号,大括号内代表组成面的围线的编号。

在建立二维模型的过程中,依次对点、线、围线和面进行设置。建立一个底边为2,高为0.5的内凹示例模型,其geo脚本文件如下,运行后产生的几何模型如图10a所示:

Point(1) = {0, 0, 0, 0.5};

Point(2) = {0, 2, 0, 0.5};

Point(3) = {0.5, 2, 0, 0.5};

Point(4) = {1.5, 1.5, 0, 0.5};

Point(5) = {2.5, 2, 0, 0.5};

Point(6) = {3, 2, 0, 0.5};

Point(7) = {3, 0, 0, 0.5};

Line(1) = {1,2} ;

Line(2) = {2,3} ;

Line(3) = {3,4} ;

Line(4) = {4,5} ;

Line(5) = {5,6} ;

Line(6) = {6,7} ;

Line(7) = {7,1} ;

Curve Loop(1) = {1,2,3,4,5,6,7} ;

Plane Surface(1) = {1} ;

选用非结构化剖分算法对所建立模型作网格剖分。Gmsh提供了几种非结构化剖分算法:自适应算法、Delaunay算法和波前法,本文采用Delaunay算法对示例模型进行剖分,结果如图10b所示。在地球物理正演建模过程中,需要在建立地质模型时加上测点和源点,同时在测点和源点周围进行网格加密以减少计算误差。Gmsh中能进行网格加密的工具有很多,最基础的就是对点周围的剖分网格大小进行控制,即改变点设置语句中的大括号内expression的值。由于地球物理数值模拟的计算程序里还需要模型的属性信息,用户可以通过Gmsh 提供的物理组来定义,将基本的几何实体组合成不同的物理组。在示例模型中加入2个孤立的点作为源点,并对其周围的网格进行加密,如图10c、d所示。

最后,需要将网格数据输出。Gmsh提供多种类型的网格数据输出,如msh、med、mesh等。本文采用mesh数据格式,将示例模型进行加密和剖分后得到的mesh文件如下:

MeshVersionFormatted 2

Dimension

3

Vertices

209

0 0 0 1

…

1.21 1.21 0 1

Edges

22

1 10 0

…

24 1 0

Triangles

394

121 53 142 0

…

131 69 202 0

End

其中,第一行代表mesh格式网格剖分数据,第二行和第三行代表为三维数据。从Vertices开始输出节点数据,先输出节点总数,然后依次输出x、y、z坐标和点的编号,若是点处在线上(面上)则输出线(面)的编号。从Edges开始输出边界单元信息,先输出边界单元总数,然后依次输出边界单元位于边界上边的2个节点的编号和所处边设置的物理组标号。从Triangles开始输出单元信息,先输出边单元总数,然后依次输出单元的3个节点的编号和所处面设置的物理组标号。以End结束数据。

mesh格式网格剖分数据满足地球物理正演的需要,其操作的便捷性也大大简化了地质建模和网格剖分的流程。Gmsh建立的二维地质模型和网格剖分示例如正文中图4a所示。