三维电磁扩散场数值模拟及磁化效应的影响*

李志旋 岳明鑫† 周官群2)

1) (中国科学技术大学地球和空间科学学院, 合肥 230026)

2) (安徽惠洲地质安全研究院股份有限公司, 合肥 231202)

(2018 年 8 月 21 日收到; 2018 年 12 月 3 日收到修改稿)

采用矢量有限元法实现了三维电磁扩散场数值模拟, 并成功将其应用在大地电磁的正演研究中. 为灵活精确地拟合起伏地形和地下不规则构造, 采用由不规则四面体单元组成的非结构化网格, 可根据模型设计的需要调整网格的大小. 引入了基于二次场理论, 将解析的一次场从总场中扣除, 直接计算二次场, 使得误差仅局限于相对较小的二次场, 以提高总场计算精度. 常规的节点有限元法不满足电性分界面上法向电场不连续和无源区单元内电流密度无散, 违反麦克斯韦方程组. 为克服节点有限元法的弊端, 使用矢量有限元法求解基于二次电场的偏微分方程. 另外, 在算法设计中, 考虑了磁导率参数的变化, 可以模拟磁导率不均匀的模型.通过与COMMEMI模型已发表的结果对比, 证明了本文算法的正确性和精确性. 为突显非结构网格优势, 计算了椭球异常体模型和任意地形模型的MT响应, 并详细讨论了地形和磁化效应对三维数值模拟结果的影响.

1 引 言

电磁扩散场数值模拟方法中, 积分方程法、有限差分法和有限元法是三种经典并广泛使用的方法. 积分方程法[1−3]受并矢格林函数的复杂性限制,仅适用于简单规则模型的模拟. 有限差分法[4−6]的优点是易于实现, 其方程的离散形式决定了该方法只能模拟可四边形或六面体剖分的计算域, 对于复杂模型, 如起伏地形和不规则地下结构, 采取台阶状近似误差大, 且剖分难实现. 相对而言, 采用非结构化网格时, 有限元法则能够灵活精确地处理各种复杂模型[7].

1971年, 有限元法最早由Coggon[8]引入地球物理电磁场数值模拟. 早期的有限元法使用规则四边形和六面体剖分计算域. 为处理倾斜的地形,Fox等[9]和Wannamaker等[10]将四边形单元分成四个三角单元, 该方法本质上仍是基于结构化网格剖分, 对于复杂的模型仍难以剖分. 非结构化网格的采用, 使得这一问题迎刃而解, 拓展了有限元模拟的应用范围[11−13]. 另外, 直接计算二次场以提高计算精度, 基于二次场的算法其一次场可解析地求得, 不存在误差, 虽然二次场由数值法解得, 存在误差, 但二次场相对于一次场较小, 误差比重不大,从而整体上提高了总场的精度.

可从两个角度来求解三维电磁场问题: 其一,基于势(场的矢势和标势)求解, 然后对势微分求场, 微分过程会影响解的精度[14]; 其二, 直接基于场求解. 本文采用后者, 从二次电场的偏微分方程来求解. 节点有限元法在处理三维矢量电磁场时,存在一些问题: 第一, 不满足电性分界面上法向电场不连续; 第二, 不满足无源区单元内电流密度散度为零. 上述问题直接违反麦克斯韦方程组, 造成伪解. Jin[15]和柳建新等[16]引入散度校正来压制伪解, 但不能完全消除. 本文采用矢量有限元法, 使用Nédélec型矢量形函数, 很好地克服了节点有限元法的弊端. 另外, 在算法设计中, 考虑了磁导率参数的变化, 可以模拟磁导率不均匀地对地球电磁扩散产生的影响.

作为研究地球内部电性构造的一种重要的地球物理手段, 大地电磁(MT) 法广泛应用于多个领域. 深部低频段, 地壳和上地幔的MT电阻率结构成像, 对岩石圈深部和地球动力学的演化过程的研究起着重要作用[17,18]; 浅部高频段, 广泛应用于金属矿产资源勘探、地热勘探、工程调查和地下水检测等方面[19−21]. 本文首先在前人的基础上, 采用更为灵活的非结构化网格剖分方式, 导出了三维MT基于二次电场的矢量有限元公式及阻抗求取公式, 对复杂不规则模型进行了模拟, 取得了理想的结果; 然后, 讨论了三维地形对模拟结果的影响,通过与二维的模型正演结果的对比, 表明了对三维数据进行二维处理的不合理; 最后, 详尽分析了不均匀磁导率对模拟结果的影响.

2 正演理论分析

2.1 方程推导

式 中E和H为 分 别 电 场 和 磁 场 ;σ为 电 导 率 ;µ=µrµ0为磁导率,µ0为真空中磁导率,µr为相对磁导率;ε为真空中介电常数;ω为角频率. 由方程 (1)可导出

将总场分成一次场和二次场, 即 E =Ep+Es , 其中下标p和s分别表示一次场(primary field)和二次场(secondary field). 在背景模型中(见图1),背景电导率和磁导率分别为 σ0 和 µ0 :

由 (2)和(3)式联合, 考虑磁导率的差异, 可得

图1 背景模型示意图Fig.1. Schematic diagram for background model.

2.2 边界条件

得到微分方程后, 还要有边界条件才能组成边值问题. 当背景模型和实际一维模型不同时, 比如实际模型为三层模型(含异常体), 将背景模型设为两层模型, 假设外边界距离异常体足够远, 异常体对边界的影响可忽略, 边界上二次场为实际一维模型和背景模型一次场的差:

式中下标1和2分别表示背景模型和实际一维模型. 当背景模型和实际一维模型相同时, 外边界上仅有一次场, 二次场为零, (5)式变成为齐次第一类边界条件; 当地下模型为二维模型, 其中包含一个三维异常体时, 仍将一维模型设为背景模型, 用二维有限元法计算出边界上二维模型的数值解, 代入(5)式中的本文利用走向方向延伸的长度远大于剖面(与走向垂直)上的尺寸的三维模型来近似二维模型, 故背景模型均为一维层状介质模型.

2.3 一次场的计算

背景模型(见图1), 其一次场可解析地求出[16].为避免由电场和磁场中指数增加项引起的计算机数据类型的越界和计算机对较大数据的截断误差,假设源为x方向, 将一次电场和磁场的形式解设为

令

空气层中, 在模型顶面z=z1上:

可解得

由(7)式边界条件, 可得

将a1和b1代入(9)式依次求出求各层的系数ai和bi, 将ai和bi代入 (6) 式, 可计算模型中任意点的一次场.

2.4 加权余量法

通过加权余量法将微分方程过渡到有限元方程, 进而采取有限元法来求解. 本文使用加权余量法,取权重为 δ Es , dv表示单元体积分, 由 (4)式可得

利用矢量运算公式∇·(A×B)=∇×A·B-∇×B·A, 边界上满足第一类边界条件, 且单元内 物性参数为常数, (10)式可化为

2.5 矢量有限元法求解

节点有限元法能有效地模拟直流问题, 但处理交变的矢量电磁场时, 存在一些问题: 第一, 不满足无源的区域单元内电流密度的散度为零; 第二,在电性分界面上, 不满足法向电流密度连续, 即法向电场不连续. 这些问题明显违反麦克斯韦方程组, 是节点有限元法理论上的缺陷. 虽然引入一些校正手段, 如散度校正, 可以在一定程度上压制伪解, 但并不能从根本上消除. 矢量有限元法由于将场赋到棱边上, 自然满足切向场的连续, 对法向场没有强制性的约束, 允许法向场的不连续, 为法向电流密度连续提供了条件. 本文采用非结构四面体单元网格进行空间离散(Gmsh开源软件产生), 选用Nédélec型矢量形函数[15]:

式中 li 为第 i条边的长度; L i1 和 L i2 分别为第 i 条边的 i1 和 i2 节点的局部坐标(体积坐标).

由(13)式可知, 矢量有限元法自然满足无源区单元内电流密度的无散.

三维四面体单元中棱边的编号规则如图2所示, 单元内任一点的场可由形函数(12)式和插值公式(14)得到:

图2 四面体单元及其棱边编号示意图Fig.2. Schematic diagram for tetrahedral element and edge numbers.

将(14)式代入有限元公式(11)可得

式 中Es和Ep为 6×1 的 标 量 列 向 量 ;K1 和K2 为6×6的矩阵,

方程(15)的右端项为源项, 其中Ep为已知的一次场棱边值. 在计算一次场时, 只能计算具体某节点的值, 而不能直接得到棱边值, 需要将节点的一次场转化到棱边上. 对于四面体单元, 已知四个节点的一次场, 如何求六个棱边的值, 有两种处理方法: 方法一, 解方程法, 四个节点, 每个节点一次场分为x,y和z三个分量, 即总共 12 个方程, 求6个未知量(棱边值), 显然是一个矛盾方程, 可在方程的两边同时左乘系数矩阵的转置, 得到一个近似解; 方法二, 投影法, 棱边中点的一次场向棱边方向的投影值近似为棱边值. 节点往棱边上转化时, 不是严格的一一对应, 上述两种方法都要近似,本文采用易于操作实现的方法二.

将局部单元刚度矩阵扩展成总体刚度矩阵, 其

有限元方程:

边界条件(5)式的加载, 对方程(16)式系数矩阵和右端项做部分修改[22]. (16)式方程条件数较大, 采用迭代法求解收敛慢[23], 特别是低频部分.本文利用PARDISO求解器(MKL库函数)进行求解, 系数矩阵采用CSR格式压缩存储, 算法速度快, 稳定性强.

3 MT 响应的计算

由(16)式得到所有棱边值后, 通过(14)式的插值公式可计算计算域内任一点二次电场, 将二次电场代入(1)式麦克斯韦方程, 可计算出二次磁场(辅助场). 将二次场和解析的一次场相加得到总场,通过(17)式计算阻抗Z:

目前国内学者仍多沿用处理二维的方式处理三维问题, 默认Zxx=0 和Zyy=0 , 则Zxy=Ex/Hy和Zyx=Ey/Hx. 在三维情况下Zxx和Zyy都存在,不等于零的, 上述方法必然会引入误差. 本文采用两个正交的源, 即x和y方向, 对模型两次模拟得到 两组数据, 由(18)式计算阻抗张量:

式中电场和磁场的下标1和2分别表示x和y方向不同的源模拟所得数据. 得到阻抗张量后, 可通过(19)式计算视电阻率和相位:

4 算法验证

20世纪80年代以来, 随着计算机计算能力的提升, 三维电磁场数值模拟方法随之快速发展. 对于MT问题, 一维层状介质和个别简单二维模型具有解析解, 三维模型则不存在解析解. 对于三维模型, 当不同的数值模拟方法得到不同的结果时,孰优孰劣, 没有一个客观的比对标准. COMMEMI(Comparison of Modelling Methods for Electromagnetic Induction)是一个众多学者参与的国际合作项目[23], 不同的学者采用不同数值模拟方法计算相同模型的MT响应, 然后计算均值和标准差,为新的算法提供一个比对标准, 当新的算法的模拟结果落在参考结果的标准差以内, 认为该算法模拟结果正确.

4.1 一维层状模型

一维模型如图3所示. 第一层电阻率为100 Ω·m ,层厚度为2 km; 第二层为基底层延伸到到无穷远,电阻率为400 Ω·m; 该模型属于G型地电模型.该模型剖分为31154个节点182806个单元, 单频的计算时间为8.84 s, 测点位于模型中心, 离边界距离为20 km. 图 4为基于总场和二次场的二次插值的模拟结果与解析解的对比. 由图4可看出, 总

图3 一维层状介质模型Fig.3. One-dimensional layered medium model.

场和二次场算法的模拟结果都与解析解吻合, 进一步研究发现, 相对于总场算法, 二次场算法的模拟结果更为精确, 其中相位曲线的差异较为明显. 这是由于基于二次场的算法, 其一次场可解析地求得, 不存在误差, 虽然二次场由数值法解得存在误差, 但二次场相对于一次场较小, 误差比重不大,从而提高了总场的精度. 因此相比于总场算法, 二次场算法的误差较小, 结果也较为理想.

图4 基于总场和二次场的二次插值的模拟结果与解析解的对比T为周期) (a) 视电阻率 ρs; (b) 相位 ϕFig.4. Comparisons of modelling results of quadratic interpolation and analytical solutions based on total fields and secondary fields: (a) Apparent resistivity ρs; (b) phase ϕ .

图5 COMMEMI-3D1 模型, 其中 (a)图中虚线为 x 和 y 测线Fig.5. COMMEMI-3D1 model. The dashed lines in (a) are the x and y survey lines.

4.2 三维模型COMMEMI-3D1

对于三维地电模型, 本文采用COMMEMI-3D1 (见图 5, 其中ρ0,ρ分别代表背景和异常体的电阻率)模型验证算法的正确性, 生成的四面体网格82992个, 由图6可知(三角符号为平均值, 竖线为标准差), 对于 0.1 Hz的剖面, XY 模式和YX 模式,x测线和y测线, 模拟的结果均在参考结果的标准差以内, 且较为靠近均值; 对于10 Hz的剖面(见图7), YX模式勉强落在标准差以内, 对于XY模式,x测线上−500—500 m和y测线上−1000—1000 m(即异常所在的区域), 模拟的结果明显比参考结果偏低. Mitsuhata 和 Uchida[25]基于势 (电流密度的矢势和标势)的算法, 采用结构化六面体网格, 混合使用矢量有限元和节点有限元法, 模拟得到的XY模式10 Hz剖面的结果相对于参考结果偏低, 为了验证其结果, Mitsuhata 和 Uchida[25]还用交错网格有限差分算法进行了模拟, 结果也偏低, 本文结果与其结论相吻合. 综合以上讨论, 验证了本文算法的正确性.

4.3 三维模型COMMEMI-3D2

为了进一步验证本文算法的精度, 本文计算了COMMEMI-3D2 (见图8)并将结果和基于电场标量势和磁场矢量势的T-Ω算法进行对比[25]. 该模型生成的四面体网格384992个, 计算频率0.001 Hz,耗时 23.86 s, 由图 9 可知, 对于x方向测线, XY模式下本文结果和T-Ω算法吻合良好.

图6 (a), (b) x 和 (c), (d) y 测线模拟结果和 Zhdanov 等[24]的结果对比 (剖面为 0.1 Hz)Fig.6. Comparisons of modelling results and Zhdanov et al[24] of (a), (b) x survey line and and (c), (d) y survey line (profile is 0.1 Hz).

图7 (a), (b) x 和 (c), (d) y 测线模拟结果和 Zhdanov 等[24]的结果对比 (剖面为 10 Hz)Fig.7. Comparisons of modelling results and Zhdanov et al[24] of (a), (b) x survey line and and (c), (d) y survey line (profile is 10 Hz).

图8 COMMEMI-3D2 模型Fig.8. COMMEMI-3D2 model.

图9 模拟结果和 T -Ω 算法的结果[25]对比 (剖面为 0.001 Hz)(a) 视电阻率; (b)相位Fig.9. Comparisons of modelling results and T -Ω algorithm results[25] (profile is 10 Hz): (a) Apparent resistivity;(b) phase.

5 算 例

5.1 椭球体模型

相比于规则六面体网格, 非结构化网格的长处在于模拟不规则地下结构, 现设计六面体网格难以实现的复杂模型以突显本文算法的优点. 椭球体模型如图10所示, 椭球体模型模拟结果切片图如图 11 所示, 频率代表似深度. 由图 11 可看出, 在异常体中心处及其附近的切片, XY和YX模式的视电阻率和相位剖面对异常体均有明显的反映, 且随着切片远离异常体, 异常逐渐减小.

图10 椭球体模型 (a) xy 剖面; (b) xz 剖面Fig.10. Ellipsoidal model: (a) xy profile; (b) xz profile.

图11 椭球体模型模拟结果切片图 (a) XY 模式视电阻率; (b) YX 模式视电阻率; (c) XY 模式相位; (d) YX 模式相位Fig.11. Slices of modelling results for ellipse model: (a) XY-mode apparent resistivities; (b) YX-mode apparent resistivities;(c) XY-mode phase; (d) YX-mode phase.

5.2 起伏地形模型

讨论三维地形模型之前, 先计算二维纯地形模型(见图12), 通过对比三维和二维算法的模拟结果, 来验证本文算法对带地形模型的正确性. 三维数值模拟时, 用一个xoz剖面上同图12和y方向(走向方向)延伸20 km的三维模型来近似二维模型, 二维数值模拟采用有限元法. 由图13可知, 三维和二维算法模拟结果吻合得很好, 且二维地形所引起的异常特征和Franke等[14]的结论一致, 从而证实了本文算法对带地形模型的正确性.

为对比三维地形对MT响应的影响, 仍沿用图10模型, 仅将其水平地表改成起伏地表(见图14),山峰底部长宽均为 3 km, 顶部长宽均为 1 km, 高

图12 二维纯地形模型Fig.12. Two-dimensional pure topographical model.

图13 三维和二维模拟结果对比Fig.13. Comparisions of three-dimensional and two-dimensional modelling results.

图14 山峰模型及其非结构化网格剖分(左图点为测线)Fig.14. Peak model and discretization using unstructured grids (the points in left diagram are survey line).

为 0.5 km, 以原点为中心对称分布, 类同 Ren 等[26]的带地形模型, 测线沿着x轴分布, 生成有限元网格20355个. 对于图14带地形地电模型, 其背景模型为两层模型, 水平地表为空气和地下介质的分界面, 山峰看作空气层中的异常体来处理. 由图14可看出, 在异常体边界及地表附近场变化剧烈的区域对网格进行了加密, 另外将黄色计算域进行了适当的向外延伸, 使边界条件更为贴近实际.

在实际生产中, 考虑到数据三维处理的复杂性和效率, 往往将三维数据逐个剖面二维处理, 然后拼接在一起. 现将三维地形模型简化成测线所在剖面 (y= 0 剖面)二维地形模型 (见图 12), 对比二维算法与三维算法正演模拟的结果差异, 为全面讨论, 另加一个同图14一样尺寸的山谷模型. 为突显地形影响, 背景电阻率和椭球电阻率均设为100 Ω·m , 仅含地形因素导致视电阻率异常. 由图 15(a)和图 15(b)可知, 对于山峰模型, 三维的XY模式和二维的TM模式的模拟结果趋势大致相同, 但具体量值有区别; 三维的YX模式和二维的TE模式几乎完全不同, 随着频率的降低差别越来越大, 电阻率异常完全相反, 剖面为 0.1 Hz时,二维显示轻微的正异常, 三维则显示为负异常, 这是由于二维模型认为走向方向上无限延伸, 不存在地形起伏, 而三维模型包含地形在走向上的变化,深部低频段电阻率负异常恰恰是由于走向方向的地形起伏造成的. 由图13(c)和图13(d)可看出,对于山谷模型, 结论和山峰模型一致, 但三维的XY模式和二维的TM模式差异相比山峰模型更大. 相比于二维地形, TE模式高频段和TM模式全频段受地形影响, 显然, 三维地形对模拟结果的影响更为严重和复杂.

图15 三维和二维模拟结果对比 (a), (b)山峰模型; (c), (d)山谷模型Fig.15. Comparisions of three-dimensional and two-dimensional modelling results: (a), (b) Peak model; (c), (d) valley model.

5.3 磁导率异常模型

以往的三维数值模拟将磁导率假设为真空中的磁导率, 且认为在整个计算域均匀分布, 在磁导率异常较大的矿区该假设明显不合理. Mukherjee和Everett[27]和Ren等[26]指出磁导率对数值模拟的结果有重要的影响. 模型仍沿用图10中的椭球体模型, 将背景电阻率设为 100 Ω·m, 背景磁导率设为真空磁导率, 测点分布在x轴上, 该模型生成的四面体单元个数为15499个. 对于图16(a)和图16(b), 为突出磁导率因素的影响, 将椭球体电阻率也设为100 Ω·m, 该模型地层即为电性均匀半空间, 仅含磁导率不均匀, 模拟结果表现为视电阻率正异常, 且随着相对磁导率的增加, XY和YX模式的视电阻率异常也随着增加; 对于图16(c)和图16(d), 讨论电性参数和磁性参数均存在异常的模型响应, 将椭球体电阻率设为 0.5 Ω·m, 椭球体电阻率引起的视电阻率负异常被磁导率引起的正异常抵消, 且负异常随着相对磁导率的增加逐渐减小. 观察有限元方程(15)的右端源项, 可分为两部分: 一部分为电阻率异常引起, 另一部分为磁导率异常引起, 且电阻率异常引起的源项随频率的降低而减小, 而磁导率异常引起的源项不随频率变化. 因此随着频率的降低, 磁导率异常对模拟结果影响越来越大, 在磁导率异常较大的区域有可能占主导地位. 以上讨论表明, 在磁导率存在异常的区域, 忽略磁导率因素, 仅考虑电阻率参数的MT数值模拟会引入很大的误差, 磁导率必须作为一个独立的参数来对待.

图16 不同磁导率模拟结果对比 (a), (b) 无电性异常模型; (c), (d)含电性异常模型Fig.16. Comparisons of modelling results for different magnetic permeability: (a), (b) Model without electrical anomaly;(c), (d) model with electrical anomaly.

6 结 论

本文采用非结构化网格来离散计算域, 能更灵活精确地拟合起伏地形和地下不规则异常体, 网格的大小可根据模型设计的需求调整, 如高频段加密网格、低频段增大网格和电性分界面加密网格等,提高精度的同时尽可能地减小计算量. 同时, 将解析的一次场从总场中分离, 直接计算二次场, 使得数值误差仅局限于相对一次场较小的二次场, 从而提高了总场的精度. 矢量有限元法满足电性分界面上法向电场不连续和无源区单元内电流密度散度为零, 克服了节点有限元法的弊端. 另外, 在推导基于二次电场的双旋度矢量方程过程中考虑了介质电导率和磁导率的不均匀性, 可以模拟磁导率不均匀模型的电磁场. 通过和COMMEMI模型的已发表结果对比, 证明了本文算法的正确性. 另外,对不规则异常体和任意地形的复杂模型的模拟, 均取得了理想的结果. 三维地形影响比二维更为严重和复杂, 用二维算法处理三维模型或数据会引入很大的误差. 在磁导率存在异常的区域, 磁导率对数值模拟的结果有着重要的影响, 磁导率必须作为一个独立的参数来对待.