复杂地形下高密度激电法2.5维有限单元法数值模拟

麻昌英，柳建新*1,，刘海飞1,，郭荣文1,

(1. 中南大学 有色金属成矿预测教育部重点实验室，长沙 410083；2. 中南大学 地球科学与信息物理学院，长沙 410083)

0 引言

激发极化法是勘探金属、非金属矿产等具有激发极化效应矿产的主要物探方法之一，在工程物探中也有着重要的应用。1971年Coggon[1]首先将有限元法应用到电磁场模拟中，自此有限单元法在地球物理模拟计算中得到广泛的应用。在我国，徐世哲[2]对有限单元法进行了深入研究。近年来，对于激发极化法正演模拟研究取得了较大的成果，主要是对激发极化法应有优先发进行三维模拟。A.Dey等[3]对三维电阻率任意形状模型展开了研究；黄俊革等[4-8]对有限元法三维地电断面电阻率正演模拟展开了深入研究，实现了齐次条件下及介质参数连续条件下三维电阻率有限元模拟；熊彬等[9-11]对复杂地形下电阻率/极化率有限单元法三维数值模拟进行了研究，实现了复杂地形三维电阻率/极化率有限元模拟。

虽然电阻率/极化率三维有限元模拟取得了很好的成果，但用有限元法解决三维问题是比较麻烦的，网格剖分量大，导致计算时间较慢，对计算机内存要求较高。在2.5维问题中，介质为二维构造电性特征，场源为三维场，通过傅里叶变换，将三维电场问题转化成带参数的二维问题求解，从而降低网格剖分量，节省计算时间。目前，2.5维有限元模拟在直流电阻率模拟中应用已较成熟，徐世哲[2,12]对2.5维稳定电流场问题进行了深入研究；阮百尧等[13]实现了三角元部分电导率分块连续变化2.5维电场有限元模拟；汤井田等[14-15]实现了基于非结构化网格的2.5维直流电阻率法自适应有限元模拟。对于高密度激电法模拟，一条剖面上一般模拟点数较大，因此对模拟计算效率要求较高，2.5维模拟相对于三维模拟计算量要小很多，因而计算效率相对于三维模拟较明显。本次研究在参考上述的基础上，利用等效电阻率法实现高密度激电法2.5维有限元数值模拟。

1 点源场的变分问题

点源二维稳定电流场的边值问题

(1)

采用有限元法求解点源2.5D稳定电流场的边值问题，首先需将方程(1)转化为等价的变分问题

(2)

然后再利用有限元法解变分问题式(2)。

2 网格剖分

对传统的二维矩形剖分，其不适合于对复杂地形的模拟，对于非结构化网格，其虽然可以很好地模拟复杂地形，正演计算精度也较高，但其网格不规则，剖分难度较大，计算量较大，不易应用于反演计算。采用结构化三角形网格，可以较好地解决网格剖分复杂度及复杂地形模拟问题，在高密度激电法中，采用规则的网格剖分有利于模拟高密度激电法的滚动测量。在纵向方向上，从上到下网格剖分间距由小到大，在保证计算精度的情况下，减少网格剖分数，可节省计算时间。

在具体实现过程中，根据同一条测线上所有测点的相对位置构成横向上的网格剖分，并根据测深点所能达到的勘探深度，确定研究区域的纵向网格剖分，再根据每个测深点的供电和测量极距关系进行网格剖分，构成正演模拟网格(如图1中的粗实线网格)。为保证正演模拟的精度，在网格的左、右和下边缘做一定程度的网格外延，如图1所示。

图1 网格剖分示意图Fig．1 Mesh sketch map

3 有限单元法

将积分区域剖分成许多三角单元，总节点数为N(图1)，则式(2)中对区域Ω和边界Γ∞的积分可分界为对各三角单元e和Γe的积分之和则有：

K0(λr)dΓ

(3)

假设三角单元内转换电位V和电导率σ线性变化，即在每个三角单元有：

(4)

式(4)中，下标1、2和3是三角单元的三角节点号；σ=(σ1,σ2,σ3)T为节点上的电导率；V=(V1,V2,V3)T为三角单元节点上的波数域电位；N=(N1,N2,N3)T，Nl=(alx+blz+cl)/2Δ为x和z的线性函数(l=1,2,3)，其中Δ=(a1b2-a2b1)为三角单元的面积，其中a1=z2-z3、a2=z3-z1、a3=z1-z2；b1=x3-x2、b2=x1-x3、b3=x2-x1；c1=x2z3-x3z2、c2=x3z1-x1z3、c3=x1z2-x2z1。

根据上面所述的网格剖分，根据参考文献[2]和[13]，得到变分问题式(2)的泛函F(V)的数值表达式：

(5)

式中K是由全部三角单元和边界单元的(Ke1+Ke2+Ke3)相加组成的N×N阶对称带宽矩阵，V是由所有N个三角网格节点上的转换电位组成的列矢量，S是与源有关的列矢量。式(5)泛函F(V)对V求变分，并令其为零，就得到求波数域中各节点上电位的线性方程组:

KV=S

(6)

通过变带宽存储乔里斯基法解线性方程组(6)，便得各节点的波数域电位V。

对波数域电位V进行付氏逆变换，付氏逆变换公式为

(7)

4 模型计算

为了验证程序的正确性，计算均匀半空间ρ0=100 Ω·m的电阻率测深曲线。如图2所示，实点为本文程序模拟结果，实线为解析解，从图2中可看出，除了前三个点由于电源点影响外，其余各点拟合很好。

图2 均匀半空间测深曲线Fig．2 Sounding curve of homogeneous half space

1)模型一。水平均匀半空间下，围岩电阻率ρ0=100 Ω·m，极化率η0=1%，存在一个电阻率ρ1=100 Ω·m，极化率η2=10%，沿Y方向无限延伸的4 m×4 m的正方形极化异常体，其与围岩极化率在接触带处逐渐变化(线性)。异常体顶部埋深6 m，用对称四极装置进行断面测量，点距2 m，极化率正演结果如图3所示。

图3 模型一：正演视极化率拟断面图Fig．3 Model one: polarization forward profile

2)模型二。山脊地形起伏高度4 m，左右各延伸6 m，围岩电阻率ρ0=100 Ω·m，极化率η0=1%，在山脊正下方存在一个电阻率ρ1=100 Ω·m，极化率η2=10%，沿Y方向无限延伸的4 m×4 m的正方形极化异常体，其与围岩极化率在接触带处逐渐变化(线性)。异常体顶部距离水平地面6 m，用对称四极装置进行断面测量，点距2 m，极化率正演结果如图4所示。

图4 模型二正演视极化率拟断面图Fig．4 Model two: polarization forward profile

3)模型三。山谷地形起伏高度4 m，左右各延伸6 m，围岩电阻率ρ0=100 Ω·m，极化率η0=1%，在山谷正下方存在一个电阻率ρ1=100 Ω·m，极化率η2=10%，沿Y方向无限延伸的4 m×4 m的正方形极化异常体，其与围岩极化率在接触带处逐渐变化(线性)。异常体顶部距离水平地面6 m，用对称四极装置进行断面测量，点距2 m，极化率正演结果如图5所示。

图5 模型三正演视极化率拟断面图Fig．5 Model three: polarization forward profile

4)模型四。山脊-山谷地形，起伏高度均为4 m，左右各延伸均为6 m，山脊与山谷之间为4 m水平地形连接，在此段水平地形正下方存在一个电阻率ρ1=100 Ω·m，极化率η2=10%，沿Y方向无限延伸的4 m×4 m的正方形极化异常体，其与围岩极化率在接触带处逐渐变化(线性)，围岩电阻率ρ0=100 Ω·m，极化率η0=1%。异常体顶部埋深6 m，用对称四极装置进行断面测量，点距2 m，极化率正演结果如图6所示。

图6 模型四正演视极化率拟断面图Fig．6 Model four: polarization forward profile

5 结论

在实际地质环境中，岩石、矿物的物性参数多为变化的，本研究基于矩形-三角形网格，设计单元内极化率及电阻率为双线性变化，实现了连续介质的激发极化法2.5维有限元数值模拟。计算了激发极化模型，通过三角元模拟了地形的起伏对激发极化法的影响，从上述模型计算结果可以分析，当地形起伏及起伏地形下存在极化异常体时，山脊地形与山谷山谷都影响视极化率的观测，从模型二与模型三模拟结果分析比较可知，山谷地形影响较大，山脊地形影响较小。在非结构网格系统中，网格的的不规则性质导致其较难应用于反演计算，而在矩形-三角形网格系统中，利用三角元进行正演，可以模拟复杂的地形起伏，当应用于反演计算时，可以将其矩形网格用于反演，使得反演问题相对简化。

参考文献：

[1] COGGON J H. Electromagnetic and eletrical modeling by the finite element method[J].Geophysics, 1971, 36(1): 132-155.

[2] 徐世浙. 地球物理中的有限单元法[M]. 北京：科学出版社, 1994.

[3] DEY A, MORRISION H F. Resistivity modeling for arbitrarily shaped 3-D structures[J].Geophysics, 1979, 44( 4): 753-780.

[4] 阮百尧, 黄俊革, 鲍光淑. 齐次边界条件下三维地电断面电阻率有限元数值模拟法[J].桂林工学院学报, 2002,22: 11- 14.

[5] 黄俊革, 阮百尧,鲍光淑. 基于有限单元法的三维地电断面电阻率反演[J]. 中南大学学报：自然科学版, 2004, 35(2): 295-299.

[6] 阮百尧, 熊彬, 徐世浙. 三维地电断面电阻率测深有限元数值模拟[J]. 地球科学-中国地质大学学报, 2001, 26(1): 73-77.

[7] 阮百尧, 熊彬. 电导率连续变化的三维电阻率测深有限元模拟[J]. 地球物理学报,2002, 45(1): 131-138.

[8] 黄俊革. 三维电阻率/极化率有限元正演模拟与反演成像[D]. 长沙：中南大学,2003.

[9] 熊彬, 阮百尧, 罗延钟. 复杂地形条件下直流电阻率异常三维数值模拟研究[J]. 地质与勘探, 2003, 30(4): 60-64.

[10] 林家勇, 汤井田, 丁茂斌. 复杂地形条件下激发极化有限单元法三维数值模拟[J]. 吉林大学学报：地球科学版, 2010, 40(5): 1183-1187.

[11] 林家勇, 汤井田, 李开壁，等. 双频激电法有限元数值模拟研究[J]. 物探化探计算技术, 2009, 31(3): 189-192.

[12] 徐世浙. 点电源二维电场问题中付氏变换的波数k的选择[J]. 物探化探计算技术, 1988, 10(3): 235-239.

[13] 阮百尧. 三角单元部分电导率分块连续变化点源二维电场有限元数值模拟[J]. 广西科学, 2001, 30(5): 1-3.

[14] 汤井田, 王飞燕. 基于非结构化网格的2.5D直流电阻率模拟[J]. 物探化探计算技术, 2008, 8(1): 413-418.

[15] 王飞燕. 基于非结构化网格的2.5-D直流电阻率法自适应有限元数值模拟[D]. 长沙：中南大学, 2009.