重力异常二维正演中的无网格方法_李俊杰.pdf

第 46 卷 第 6 期煤田地质与勘探Vol. 46 No.6 2018 年 12 月COALGEOLOGY 2. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China Abstract Meshfree as a kind of new numerical has the advantages of high precision, high order function structure and convenient loading of physical properties, is widely applied in the field of computational mechanics. The meshfree including PIM, RPIM and EFGM was used in the two-dimensional forward calculation of gravity anomaly field. Firstly, from the 2-D variational problem of gravity anomaly, the corresponding meshfree discrete system matrix expression was deduced by combining Galerkin and Gauss integral ula. Then suggested values of the shape parameters of RPIM-MQ, RPIM-exp and EFGM-exp were obtained by numerical experiments, the calculation results of different meshfree s under optimal shape parameters were compared and analyzed. The conclusions are verified as follow Meshfree is suitable for 2-D forward modeling of gravity anomaly with big variation of the distribution of medium physical properties, the optimal interval of exponential function shape parameter αcis in range of 1.5 to 1.7, the proposed value of β is 0.6 and the proposed interval of q in MQ function is in range of – 4.1 to 1.9, EFGM has higher computational accuracy than PIM and RPIM. Keywords meshfree ; point interpolation ; radial point interpolation ; element-free Galerkin ; gravity anomaly 重力异常二维正演一般采用解析积分法， 然而 对于非均匀、 几何形状复杂的地质体需要利用数值 积分进行近似的正演计算， 常用的方法是基于三角 剖分的有限元法[1-4]，但高质量的三角剖分程序实 现过程较复杂。无网格法[5]作为一类新型节点数值 算法，虽然其发展历史仅 20 余年，但其具有形函 数构造无需网格，精度高，计算一般复杂模型及连 续介质模型较网格方法便利等优势， 引起了国内外 学者浓厚的兴趣， 现已成为数值计算领域一大研究 热点。 ChaoXing 182煤田地质与勘探第 46 卷 无网格法在弹性波场[6-8]、雷达波场[9]、大地电 磁场[10-17]、直流电场[18]及磁场[19]的正演模拟中取得 了一定的成效， 研究结果表明当地下介质为一维时， 数值方法能取得很高的计算精度[10-11]，但对应的高 维问题不存在解析解，计算复杂模型时只能通过断 面图异常的形态来近似反映算法的优劣。由于规则 二度体的重力异常存在解析解，将无网格法用于此 类模型的数值计算有利于探究高维问题下无网格法 的计算效果。 本文以重力异常二维变分问题为例，简述了无 网格法的近似原理，推导了变分问题的无网格总体 矩阵表达式，通过圆柱二度体的计算与数值试验， 得出了无网格法形状参数的最优取值区间，验证了 算法的有效性。 1重力异常二维变分问题 取走向为 y 轴，z 轴垂直向下，x 轴水平向右， o 为坐标原点，求解域为矩形区域，4 个顶点依次以 ABCD 逆时针编号图 1。 图 1重力异常二维模型 Fig.12-D model of gravity anomaly 当地下电性结构为二维时，重力异常满足如下 的变分问题[2]     2 2 sin dd 2 11 4 2sin 0 π gi FggG r F g g z                  1 式中  是二维哈密顿算子， xy    ＝ xy ee ；为 密度差；为求解域，为无穷远边界，H为求 解域深度，h为异常体质心的深度，r为质心到无 穷远边界上的距离；为r与水平线的夹角，i为 边界外法线n与矢径 r 的夹角图1。 2重力变分问题的无网格解法 2.1无网格法近似方式 无网格法利用位于积分点支持域内的场节点构 造形函数[10-14]，用Galerkin法构造系统方程，近似 手段多采用点插值法Point Interpolation ，简 称PIM、径向基点插值法Radial Point Interpolation ，简称RPIM及移动最小二乘法Moving Least squares Mothed，简称MLSM，系统方程表达 式中包含的积分项可采用含背景网格的高斯积分求 解图2。 图 2无网格法背景网格、支持域、积分点与 场节点示意图 Fig.2Background grid，support domain，integral points and field nodes in meshfree 2.1.1点插值近似 求解域Ω中任意一点X处场变量uX的PIM 近似表达式为  1 T m j h jj upcu  XXXpX c2 式中c为系数向量； T 1 xyxypX为二 维线性基函数；m为单项式的个数。将式2写成矩 阵的形式有 UPc3 式中 T 1234 {}uuuuU， T 1234 {}c c c cc ，P的展 开式为 111 1 2222 3333 4444 1 1 1 1 xyx y xyx y xyx y xyx y         P4 支持域内仅含4个节点的PIM方程组未知数与 节点个数相等mn， 只需进行如式5所示的矩阵 逆运算便可求得系数向量c的值。 1 c P U5 ChaoXing 第6期李俊杰等重力异常二维正演中的无网格方法183 将式5代入式2得  T1T n h i i i uu    XpXUX UP 6 式中 T X 即为PIM形函数，其表达式为  TT 1234 1   XpXXXXPX  7 2.1.2径向基点插值近似 RPIM与PIM的差别在于前者引入了径向基函 数Radial Basis Function，简称RBF，改善了形函 数构造过程的奇异性问题。求解域Ω中任意一点X 处场变量u X的RPIM近似表达式为   1 TT 1 nm i ij h ijj uRbpcu     XXXX RX bpX c 8 式中 点插值部分各参数含义与式2相同， T RX 为径向基函数，本文RBF取MQ函数及exp函数， 其表达式如式9。 22 2 MQ exp[ / ]exp q iicc icic R x,yrd R x,yrd         9 式中 c d为与节点间距有关的特征长度，当节点均匀分 布时可取 22 ccxcy ddd；ri 22 ii xxyy， x与y为高斯点的坐标， i x与 i y为支持域内节点的 坐标； c 与q均为径向基函数的形状参数，不同学 科领域形状参数的最优值一般不同，本文的数值计算 中MQ函 数 取[ 200, 200] c  且0 c ， [ 4.1,1.9]q  且0q ，exp函数取[1.5 1.7] c ，。为 方便描述，本文将基于MQ函数的RPIM记为 RPIM-MQ，基于exp函数的RPIM记为RPIM-exp。 将式8写成矩阵的形式 URb Pc10 式中 3 T 12 {} n u uuuU ， 3 T 12 {} n b bbbb ，c T 312 {} m c ccc。R与P的展开式分别为式11、 式12。 11121131111 12222232222 13323333333 123 ,,,, ,,,, ,,,, ,,,, n n n nnnnnnnnn R x yRx yRx yRx y R xyRxyRxyRxy R xyRxyRxyRxy R xyRxyRxyRxy             R 11 111 11 22222 33333 1 1 1 1 m m m nnnnmn xyx yp xyx yp xyx yp xyx yp             x x Px x 12 式10中有nm个变量，需添加如式13所示的m 个约束条件使其变为方阵。  T 1 0 ,1, 2 , , iji n i jpbm    XP b 13 联立式10与式13可得到如式14的矩阵方程 T g           RP b Ga 0cP0 U U14 式中 T 132 {0}00 ng uu uuU； 121 { n b bbca T 2 } m cc。求解式14得式15 1 g         aG c U b 15 将式15代入式8得。  1TTT h ggg u  XRXpGXUX U 16 式中  T 121 gnn   XXXXX   n m   X，取 T g X 前n项即为RPIM形函数。 2.1.3移动最小二乘近似 无单元Galerkin法采用移动最小二乘法MLS 构造形函数，MLS近似与PIM近似区别在于前者的 系数向量为一变量，MLS近似中的节点数通常大于 基函数个数 n m ， 系数向量c X需通过如式17 的加权残差法求得  T2 1 [] 0 n i i fWu f            ii XXpXc X c   17 式中W i XX 为权函数， 简写为WiX， 本文选 择三次样条函数spline及指数函数exp，表达式如 式18、式19。  23 2 3 2 440.5 3 44 440.51 3 3 01 iii ii i i LLL WLLL L L             ≤ ≤ ii X18  2 / e1 01 i L i i L W L       ≤ i X19 式中/ i LDr＝ i ；Di为高斯点X与节点 i X之间的 距离，|| i D＝ i XX，r为权函数支持域尺寸，为 形状参数，本文数值计算中取0.6。与RPIM一 样，将基于exp函数的EFGM记为EFGM-exp， 将 基于spline函数的RPIM记为RPIM-spline。 通过式17的运算将产生如式20的线性关系。  A X c XB X U20 式中 T 12 {} n u uuU，A X与B X表达式如式 21、式22。  1 T n i i W   ii A XX p XpX21 ChaoXing 184煤田地质与勘探第46卷    12n WWW  12n B XX p XX p XX p X 22 求解式20得  1 c XAX B X U23 将式23代入式17得  1 T n h i i i uu    XX U 24 T X 即为EFGM形函数，节点i的形函数 i X表达式为   1 1 T1 m ij ji j i p        XXAX B X pXA B 25 2.2无网格离散系统方程的构造 将变分符号代入式1，代入无网格法构造的形 函数[10-14]，得   T T T 4d sin d0 sin π g F ggg z i g r G g                  26    T 1212 n 11 I nn IInn gNNN N ggg g      gN 27 11 n IInn I gNg    g 28 式中 为形函数矩阵；n为支持域内的节点数；g 为支持域内n个节点的重力场向量。将式27、式28 代入式8，式26变为如下形式 T T d sin dd sin 4π nn JJII I IJ n I IJJI I G NNNN g xxzz Ni N Ngg rz                            29 式29采用支持域内n个场节点编号， 当节点I 和节点J位于不同支持域时， 相应的被积函数为零， 故可得到如式30所示按求解域全部节点1M编号 的总体刚度矩阵。 T T T T T d sin d si 0 n 4 d M M JJII I IJ IJJ M I I I NNNN g xxzz i N Ng r N g z G                                                K G KG PG P KPGG 30 由于 T G是任意的，故式30成立的条件为 KGP31 式31即为无网格法构造的系统方程。由于求 解域为矩形，在左右边界有i，此情况下左右边 界的积分为零，如图1所示，有 11 22 1 11 22 2 22 22 2 22 22 2 sincos |||| cossin sincos |||| cossin HhHh i r Hhx xx i r Hhx hh i r hx xx i r hx                         32 式中H为求解域深度，h为异常体质心的深度， 1 r为 质心到底边界的距离， 2 r为质心到地面的距离。将K包 含的无穷边界积分项展开成底边界与表面界积分之和 的形式，利用式32，式30中的边界积项可表示为 1 2 22 32 22 22 sinsin coscos sin dd sinsin sin coscos sin d sin d d BC AD BC AD IJIJ IJ IJ IJ iii N NN N rr ii N N r H hx N N H hx H h hx N N hx h                         33 2.3求解域背景网格积分 K的表达式中包含对求解域与求解域边界 的积分，可采用如图2所示含背影网格的高斯积 分求解，如式34 222 1 22 1 2 1 || ]|| |4| [ g c c g c n n DJJII IJiIJik k nn B iIJilI l n n DI iik k i i Qi i NNNN KN NJ xxyy Hhxhx N NJP rHhr h N GJ z                                X g 34 式中 c n与 c n  分别为背景网格及边界网格的个数； g n与 g n  分别为背景网格及边界网格内的高斯点数 量，本文中EFGM每个背景网格采用4个高斯点， 每个边界网格采用2个高斯点；PIM与RPIM采用 的高斯数目为EFGM的2倍； i 为权系数； D ik J与 B il J 分别为对应的雅可比Jacobian矩阵。 3算法验证与精度分析 设均匀介质中存在一个圆柱二度异常体，异常 ChaoXing 第6期李俊杰等重力异常二维正演中的无网格方法185 体密度差 3 2 g/cm，截面半径25 mR ，埋深 50 mh 。场节点等间距分布于求解域，PIM与 RPIM采用6 5618181个场节点，6 4008080个 背影网格，横纵向节点间距均为12.5 m。EFGM采 用1 6814141个场节点，1 6004040个背影网格， 横纵向节点间距均为25 m。 如图3所示，exp函数形状参数的取值对RPIM计 算精度影响较大， 不同形参计算结果在200 mX 的区 域附近存在较大差异， c 最优取值区间为[1.5，1.7]。 图 3不同形状参数 RPIM-exp 圆柱二度体重力异常 计算结果 Fig.3Calculation results of gravity anomaly of RPIM-exp for cylindrical two-dimensional abnormal body in different shape parameters 图4为不同形状参数的EFGM-exp计算结果， exp函数形状参数的取值对RPIM计算精度影响较 大， 随着值的增大， 曲线向上偏离解析解，0.8 时异常体中心横坐标0 mX 处的异常值明显小于 解析解，0.6计算效果最佳。 图4不同形状参数EFGM-exp圆柱二度体重力异常场计 算结果 Fig.4Calculation results of gravity anomaly of EFGM-exp for cylindrical two-dimensional abnormal body in different shape parameters 图5为最优形状参数下不同无网格法的圆柱二度体 重力异常场计算结果，插值型无网格法采用四点高斯积 分公式，EFGM采用二点高斯积分公式。EFGM能达到 插值型无网格法高斯点与计算点增加一倍的计算精度， 原因在于EFGM能够通过权函数的适当选取获得高阶 次的形函数，从而提高了近似精度与计算效率。 数值模拟表明无网格法计算结果与解析解吻合， 均较好地反映出了重力异常场的分布特性，PIM计算 结果与RPIM-MQ相同， 验证了无网格法在重力异常二 维正演中的有效性。数值试验过程中还发现MQ函数 形状参数具有较大的选取区间， c 可在[ 200,200]区 间取值，q取值区间[ 4.1,1.9]，不同形状参数的 RPIM-MQ计算出的重力异常曲线重合，超过此区间， RPIM计算会产生奇异。为保证RPIM-MQ在各种模型 下的普适性， c 与q不建议在区间两端附近取值。 图 5最优形状参数值的无网格法圆柱二度体重力 异常场计算结果 Fig.5Calculation results of gravity anomaly of meshfree for cylindrical two-dimensional abnormal body in optimal shape parameters 4结 论 a. 将PIM、RPIM及EFGM 3种无网格法应用 于重力异常场二维正演，简述了无网格法形函数的 构造过程，推导了重力异常二维变分问题的无网格 系统矩阵表达式。 b.RPIM-MQ形状参数具有较大的选择区间， 本 文 的 数 值 试 验 表 明 当0.98q 时 ， c 在 [ 200,200]区间取不为零的值均能得到一致的计算 结果，q的取值区间远小于 c ，当1.5 c 时，q值 可在[ 4.1,1.9]区间取非零值，超出此区间计算将产 生奇异。RPIM-exp形状参数 c 最优取值区间为 [1.5,1.7]，EFGM-exp形状参数建议值为0.6。 c.PIM与RPIM-MQ计算曲线重合，相同节点 分布下，EFGM模拟二维重力异常场精度高于插值 型无网格法。基于二点高斯积分公式的EFGM计算 效果可与基于四点高斯积分公式及网格加密一倍的 插值型无网格法相当。无网格法计算结果与解析解 较吻合，验证了算法的正确性。 ChaoXing 186煤田地质与勘探第46卷 参考文献 [1] 金钢燮，胡祥云，超敬来，等. 复杂形体重磁异常的等参数有 限元积分算法研究[J]. 石油地球物理勘探，2009，442 231–239. KIM Kangsop，HU Xiangyun，CHO Gyonglae，et al. Study on isoparametric finite-element integral algorithm of gravity and magneticanomalyforbodywithcomplexshape[J].Oil Geophysical Prospecting，2009，442231–239. [2] 朱自强，曾思红，鲁光银，等. 二度体的重力张量有限元正演 模拟[J]. 物探与化探，2010，345668–672. ZHU Ziqiang，ZENG Sihong，LU Guangyin，et al. Finite element forward simulation of the two-dimensional gravity gradient tensor[J]. GeophysicalandGeochemicalExploration，2010，345668–672. [3] 曾思红，朱自强，鲁光银. 基于 DCT 有限元二维重力张量的 正演[J]. 物探化探计算技术，2010，324392–396. ZENG Sihong ， ZHU Ziqiang ， LU Guangyin. The forward modeling of two-dimensional gravity gradient tensor by finite element based on discrete cosine trans[J]. Computing Techniques for Geophysical and Geochemical Exploration ， 2010，324392–396. [4] 李焓， 邱之云， 王万银. 复杂形体重、 磁异常正演问题综述[J]. 物探与化探，2008，32136–43. LI Han，QIU Zhiyun，WANG Wanyin. A review of the forward calculation of gravity and magnetic anomalies caused by irregular models[J]. Geophysical and Geochemical Exploration， 2008，32136–43. [5] 李俊杰，严家斌. 无网格法进展及其在地球物理学中的应用[J]. 地球物理学进展，2014，291452–461. LI Junjie，YAN Jiabin. Developments of meshless and applications in geophysics[J]. Progress in Geophysics，2014， 291452–461. [6] 贾晓峰， 胡天跃，王润秋. 复杂介质中地震波模拟的无单元法[J]. 石油地球物理勘探，2006，41143–48. JIA Xiaofeng ， HU Tianyue ， WANG Runqiu. A mesh-free algorithm of seismic wave simulation in complex medium[J]. Oil Geophysical Prospecting，2006，41143–48. [7] 周震，贾晓峰. 基于 OpenMP 加速的无单元逆时偏移成像[J]. 物探化探计算技术，2016，386821–831. ZHOU Zhen，JIA Xiaofeng. OpenMP-accelerated element-free reverse-timemigration[J].ComputingTechniquesfor Geophysical and Geochemical Exploration ， 2016 ， 386 821–831. [8] 蒋婵君，王有学，李川，等. 求解声波方程的径向基函数无网 格法[J]. 煤田地质与勘探，2016，445136–141. JIANG Chanjun，WANG Youxue，Li Chuan，et al. Numerical solution of the acoustic wave equation using meshless with radial basis functions[J]. Coal Geology Exploration ， 2016，445136–141. [9] 冯德山， 王洪华， 戴前伟. 基于无单元 Galerkin 法探地雷达正 演模拟[J]. 地球物理学报，2013，561298–308. FENG Deshan，WANG Honghua，DAI Qianwei. Forward simulation of ground penetrating radar based on the element-free Galerkin [J]. Chinese Journal of Geophysics ，2013 ， 561298–308. [10] 严家斌，李俊杰. 无网格法在大地电磁正演计算中的应用[J]. 中南大学学报自然科学版，2014，45103513–3520. YAN Jiabin，LI Junjie. Magnetotelluric forward calculation by meshless [J]. Journal of Central South UniversityScience and Technology，2014，45103513–3520. [11] 李俊杰， 严家斌. 无网格点插值法大地电磁二维正演数值模拟[J]. 石油物探，2014，535617–626. LI Junjie， YAN Jiabin. Magnetotelluric two-dimensional forward numerical modeling by meshfree point interpolation [J]. Geophysical Prospecting for Petroleum，2014，535617–626. [12] LI Junjie， YAN Jiabin， HUANG Xiangyu. Precision of meshfree sandapplicationtoforwardmodelingof two-dimensionalelectromagneticsources[J].Applied Geophysics，2015，124503–515. [13] WITTKE J， TEZKAN B. Meshfree magnetotelluric modeling[J]. Geophysical Journal International，2014，19821255–1268. [14] 李俊杰，严家斌. 大地电磁二维正演中的有限元–径向基点插 值法[J]. 中国有色金属学报，2015，2551314–1324. LI Junjie， YAN Jiabin. Magnetotelluric two-dimensional forward by finite element-radial point interpolation [J]. The Chinese Journal of Nonferrous Metals， 2015， 255 1314–1324. [15] 李俊杰， 严家斌， 皇祥宇. 无单元 Galerkin 法大地电磁三维正 演模拟[J]. 地质与勘探，2015，515946–952. LI Junjie，YAN Jiabin，HUANG Xiangyu. Three-dimensional forward modeling of magnetotellurics using the element-free Galerkin [J]. Geology and Exploration，2015，515 946–952. [16] LI Junjie ， YAN Jiabin. Two dimensional magnetotelluric forward by meshfree local radial point interpolation [C]// Proceedingsofthe6thInternationalConferenceon Environmental and Engineering Geophysics ， Science Press USA，2014623–628. [17] 卢杰， 李予国. 无网格局部 Petrov-Galerkin 法大地电磁场二维 正演模拟[J]. 地球物理学报，2017，6031189–1200. LU Jie，LI Yuguo. Two-dimensional magnetotelluric forward modeling using the meshfree local Petrov-Galerkin [J]. Chinese Journal of Geophysics，2017，6031189–1200. [18] 麻昌英，柳建新，刘海飞，等. 基于全局弱式无单元法直流电 阻率正演模拟[J]. 地球物理学报，2017，602801–809. MA Changying，LIU Jianxin，LIU Haifei，et al. A global weak element free for direct current resistivity forward simulation[J]. Chinese Journal of Geophysics，2017，602 801–809. [19] 孔倩，李鹏，李晶. 二维位场延拓的无网格方法研究[J]. 物探 化探计算技术，2017，392155–160. KONG Qian ， LI Peng ， LI Jing. Research on element free Galerkin for 2-D potential field extension[J]. Computing Techniques for Geophysical and Geochemical Exploration ， 2017，392155–160. 责任编辑 聂爱兰 ChaoXing