无网格法精度分析及在电磁法二维正演中的应用

2015 年 12 月
APPLIED GEOPHYSICS（应用地球物理）
第 12 卷，第 4 期
无网格法精度分析及在电磁法二维正演中的应用//Precision of meshfree methods and application to forward modeling of two-dimensional electromagnetic sources, APPLIED GEOPHYSICS, 2015, 12(4), P. 503-515. DOI:10.1007/s11770-015-0511-3
全域弱式无网格法近似原理
1. 无网格参数
全域弱式无网格法求解变分问题主要包括
504
(1)
式 中 ： a = [a1 a2 a3 an ]T 为 系 数 向 量 ；
U [u1 u2 u3 un ]T 为节点场值向量； p X 为二
第 12 卷
李俊杰等，无网格法精度分析及在电磁法二维正演中的应用
(5)
3. 径向基点插值近似
RPIM 与 PIM 的差别在于前者引入了径向基 函数(radial basis function，RBF)，改善了形函数 构造过程的奇异性问题。 求解域 Ω 中任意一点 X 处 场 变 量 u ( X) 的 RPIM 近 似 表 达 式 如 式 (3) (Wang and Liu, 2002)
a， 将其代回式(1)可求得对应的形函数 ( X) 表
T
c 不同， 泊松方程数值试验表明当 q = 1.02 时，
取值区间约 [858 858] ；当 c =0.3 时， q 取值区 间约 [1.4 5.2] 且 q 0 ， q 1 ，在此区间取值结 果无明显差异，超出此区间计算将出现奇异， q 对精度的影响较 c 大， 为避免奇异性， 形状参数 取值不宜过大，推荐 c 50 ， q 取 1 或 1 附近 的值。 式(3)中有 n m 个变量， 需添加 m 个约束条 件将其转化为方阵如式(5)
u ( X) Ri X ai p j X b j i 1 j 1 T T R X a p X b U Ra + Pb
n m
Байду номын сангаас
式 中 ： U g {u1 u2 u3 un 0 0 0}
T
T
；
a g = {a1 a2 an b1 b2 bm } 。 采用与 PIM 相
元法类似， 相比其它类型无网格法其稳定性好且 计算精度高， 缺点在于系统矩阵的形成过程涉及 对求解域及边界的积分，计算效率低。 经历 20 余年的发展，无单元 Galerkin 法 (Element-free Galerkin Method, EFGM) (Belytschko et al., 1994b)作为最具影响力的全域 弱式无网格方法， 已被广泛应用于网格方法如有 限差分法(Finite deference method, FDM) (Liu et al., 2013; Dong et al., 2014) 、有限元法 (Finite element method, FEM) (Hu et al., 2015; Ren and Yan, 2010)所触及的领域，但其在地球物理学领 域的研究文献仍较少。 (Li and Yan, 2014) 综述了
无网格法精度分析及在电磁法二维正演中的应用*
李俊杰 1，严家斌 2，皇祥宇 2 （1. 浙江省水利水电勘测设计院，浙江 杭州 310002；2. 中南大学地球科学与信 息物理学院有色资源与地质灾害探查湖南省重点实验室，湖南 长沙 410083） 摘要：无网格法形函数构造不依赖预定义的单元，具有计算精度高、处理复杂模型 便利等优点。本文介绍了无单元 Galerkin 法(EFGM)、点插值法(PIM)与径向基点插 值法(RPIM)三种全域弱式无网格法的近似原理及特点；以二维泊松方程为例研究了 支持域无量纲尺寸、场节点与背景网格设置对无网格法计算精度的影响。将 RPIM 与 EFGM 应用于频率域线源二维正演,给出了 RPIM 形状参数的推荐值； 分析了均匀 介质模型大地电磁(MT)二维正演无网格法边界条件直接加载与罚函数法加载的精 度差异，结合 PIM 与 RPIM 边界条件加载便利及 EFGM 计算复杂模型精度高的优 势，提出了 EFG-PIM 及 EFG-RPIM 耦合算法，数值计算结果验证了耦合算法的有 效性。研究发现：无网格法及其耦合方法适用于电磁法数值模拟；支持域无量纲尺 寸取 1.0 时无网格法精度与效率高，场节点与背景网格重合时计算效果佳；泊松方 程求解 PIM 及 RPIM 精度较 EFGM 低，计算均匀介质 MT 响应精度较 EFGM 高； RPIM 改善了 PIM 计算涉及的奇异性问题，对应支持域无量纲尺寸选择空间大。 关键词：无单元 Galerkin 法，点插值法，径向基点插值法，泊松方程，线源二维 正演，无网格耦合法，大地电磁
第4期
维直角坐标 XT [ x, y ] 的基函数；m 为单项式个 数。 p X 可通过帕斯卡(Pascal)三角形确定，其 多项式基函数 p T X [1 x y xy ] 。 PIM 方 程 组 未 知 数 与 节 点 个 数 相 等 ( m n )，即 P 为一方阵，通过矩阵逆运算求取
式中： d c 为与节点间距有关的特征长度，节点
2 2 d cy ， d cx 与 d cy 为横 均匀分布时可取 d c d cx
纵向节点间距；ri ( x xi ) 2 ( y yi ) 2 ，x 与 y 为高斯点的坐标， xi 与 yi 为支持域内节点的坐 标； c 与 q 为形状参数，不同学科其最优值一般
引言
无网格法是计算力学领域的研究热点， 其采 用离散于求解域内的节点构造形函数， 在网格方 法求解存在网格畸变或单元分裂的特殊力学问 题如高速冲击与大爆炸(Feng et al., 2014; Ma et al., 2009)、金属成型(Liu et al., 2009; Lu et al., 2008) 、 裂 纹 扩 展 (Zhang and Chen, 2008; Belytschko et al., 1994a)上都取得了良好的应用 效果。 根据选用加权残量法的不同无网格法主要 分为全域弱式、局部弱式和配点三类。全域弱式 无网格法求解变分问题系统矩阵的加载与有限
(3)
取 ΦT g ( X) 前 n 项即得计算点 X 在支持域内的
式中点插值部分各参数含义与式(1)相同，a 与 b 为系数向量， R T X 为径向基函数，常用 RBF 有 MQ(multi-quadrics)函数、薄板样条(thin plate spline,TPS) 函 数 、 高 斯 (Gaussian) 函 数 及 对 数 (Logarithmic)函数(Liu and Gu, 2005)， (Li and Yan, 2015)通过数值试验表明基于 MQ 函数的 RPIM 在大地电磁二维正演计算中具有很好的普适性， 本文选择 MQ 函数， 其表达式如式(4) (Franke and Schaback, 1998)
收稿日期：2015-07-17；修改稿收到日期：2015-09-09 *基金项目：本研究由国家自然基金项目（编号：40874055）和湖南省自然基金项目（编号：14JJ2012）资助。 © 2015 应用地球物理编辑部，保留所有版权
第 12 卷
李俊杰等，无网格法精度分析及在电磁法二维正演中的应用
图 1 全域弱式无网格法高斯点、场节点、支持域 与背景单元示意图
2.点插值近似
求解域 Ω 中任意一点 X 处场变量 u ( X) 的 PIM 近似表达式及其矩阵形式如式(1) (Liu and Gu, 2001)
m X u p j X a j pT X a j 1 U Pa
U R P a U g T Ga g 0 P 0 b n T p j Xi ai P a = 0 , j = 1 , 2 , , m i 1
达式 u X p T X P 1U n i ui T ( X)U (2) i T ( X) p T X P 1 [1 X 2 X 3 X n X ] 点插值近似过程简单高效， 边界条件可直接 加载， 缺点在于节点分布不当时矩阵 P 易呈现奇 异性或支持域内的节点数大于基函数单项式个 数( n m )，P 不再为方阵，无法直接逆运算。
节点生成、求解域背景单元离散、形函数构造、 总体矩阵组装、边界条件加载及场值后处理过 程。图 1 为全域弱式无网格法相关参数示意图， 背景单元节点与场节点可以不重合，与 FEM 单 元概念不同，它能够独立于场节点存在，用于求 解系统方程中包含的区域及边界积分。 电磁场计 算时物性参数加载于背景单元内的高斯点中， 同 一背景单元内可具有几种不同的电阻率值。 常用支持域形状主要有矩形与圆形两种，对 于本文节点规则分布的情况多采用矩阵支持域， 其包含横纵两方向的尺寸，一般定义为支持域无 量纲尺寸 与节点间距的乘积， 对无网格法精 度与效率影响很大， 不同学科 最优值存在差异， 需通过数值试验获得。支持域内场节点被用于形 函数的构造，常用无网格近似方式主要有点插值 近似、径向基点插值近似及移动最小二乘近似。
T u X [R T X p T X ]G 1U g g ( X)U g (6)
同的矩阵逆运算求解 a g ，并代回式(3)得
式中
1 X , 2 X , , n X , gT ( X) ， n 1 X , n 2 X ,, n m X
第4期
无网格法研究进展，Jia 等(2005, 2006)研究了基 于不同吸收边界的 EFGM 弹性波场二维正演， 结果表明 EFGM 较 FDM 稳定且精度更高； (Feng et al., 2013; Dai and Wang, 2013)将 EFGM 用于地 质雷达二维正演，并与 FEM 计算结果对比，得 出了类似的结论；Yan 和 Li(2014)将 EFGM 应 用于 MT 二维正演， 得出了适用于 MT 正演支持 域无量纲尺寸与高斯点数量的建议值。 Wittke 和 Tezkan (2014) 采用结合了径向基函数的无网格 局部 Petrov–Galerkin 法(Meshless Local Petrov– Galerkin method,MLPGM) 模拟了若干电导率模 型的平面电磁波场响应,分析了 MLPGM 的特点， 显示了 MLPGM 处理复杂模型的便利性。 EFGM主要缺陷在于形函数采用拟合方法构 造，边界条件无法直接加载。(Liu and Gu, 2001; Wang and Liu, 2002)提出了形函数采用插值法构 造的无网格点插值法(point interpolation method, PIM)及径向基点插值法(Radial point interpolation method, RPIM)。PIM插值形式简单，在节点均匀 分 布 下 能 取 得 较 高 的 计 算 精 度 (Liu and Gu, 2001)，然而其计算过程涉及矩阵的逆运算，当节 点分布不合理时易产生奇异性问题。RPIM与PIM 区别在于前者引入了径向基函数， 改善了奇异性， 由于径向基函数表达式包含形状参数，不同学科 其最优值需通过数值试验获得。Li 和Yan (2014b, 2015)将PIM与RPIM用于MT二维正演， 并将RPIM 与有限元法耦合，提高了算法的计算效率。 本文讨论了影响无网格法计算精度与效率 的因素，将PIM、RPIM及EFGM用于频率域线源 及大地电磁二维正演， 分析了各种无网格法的优 缺点，结合EFGM与插值型无网格法的优势形成 了EFG-PIM及EFG-RPIM，数值计算结果验证了 算法的正确性。