球坐标系下的大地电磁正演模拟研究现状_李金良.pdf

2020年第12期西部探矿工程 * 收稿日期 2020-03-17修回日期 2020-03-18 第一作者简介 李金良 （1995-） ， 男 （汉族） ， 河南漯河人， 中南大学地球科学与信息物理学院在读硕士研究生， 研究方向 勘查地球物理。 球坐标系下的大地电磁正演模拟研究现状 李金良*1， 2， 郭荣文 1， 2， 柳建新1， 2 （1.中南大学地球科学与信息物理学院， 湖南 长沙 410083； 2.有色资源与地质灾害探查湖南省重点实验室， 湖南 长沙 410083） 摘要 为满足地球动力学研究的需要， 以及深部资源勘探的需要， 大尺度的电磁法作为一种深部勘 探技术受到了越来越多的关注。在此背景下， 地球的曲率影响增大， 传统的笛卡尔坐标系不再适用 实际情况， 因此研究人员开展了球坐标系的大尺度大地电磁法方法正、 反演理论研究。介绍三种常 用的地球物理数值模拟方法的基本原理， 以及在球坐标系下大地电磁法正演研究的现状和进展， 并 在此基础上对未来球坐标系下大地电磁法正演研究提出展望。 关键词 球坐标； 大地电磁； 正演； 数值模拟； 积分方程法； 有限元法； 有限差分法 中图分类号 P631 文献标识码 A 文章编号 1004-5716202012-0099-04 地球物理学主要通过研究地球的弹性、 磁性、 导电 性及密度分布等物理特性来判断地下物性结构特征， 它的服务对象， 大体上分几类 科学探索、 工程建设、 资 源开发等。经历了几十年的发展， 地球浅层的地球物 理勘探技术已经相对成熟， 浅层的资源开采程度也较 高， 出于科学探索和资源开发的需求， 我们必须进行更 深部的勘探技术研究。 从科学探索的角度考虑， 深部勘探主要是为了满 足地球动力学研究的需要。地球动力学的研究目的是 为了深刻认识地球本体[1]。这一领域的研究， 需要深部 勘探技术为其提供数据支持。若以20世纪70年代美 国开展的大陆反射地震探测计划COCORP计划为起 点计算, 各国对地壳和岩石圈结构的探测活动已持续 近50年 （如表1所示） 。 由于电磁场的趋肤效应， 为了获得较深部的电性 信息， 主要使用低频电磁波， 深部勘探实质上就是低频 勘探， 即长周期勘探。国内外已有许多长周期勘探的 研究和应用， 然而， 这些理论和实践中大多沿用了传统 的笛卡尔坐标系， 忽略了地球曲率对大尺度地球物理 勘探的影响。 早在大地电磁法诞生后不久， 就有研究者注意到 地球曲率对勘探效果可能造成的影响， 20世纪60年 代， Srivastava研究了平面波场源入射层状球体介质的 大地电磁正演理论， 指出对于地壳尺度范围内的大地 电磁测深， Cagniard的假设是可行的， 但随着探测深度 的加大， 地球曲率的影响也逐渐增大[2]。1972年， Wei- delt建立了一维地球的球形和笛卡尔阻抗之间的数学 关系[3]。Schmucker总结了若干地球物理模型的构建和 使用过程中的一些处理原则[4]。后续的研究者也不断 提及地球曲率对地球物理数值模拟的影响， 其中， 在球 坐标系下进行数值模拟也是一个探索的方向[5]。 国家 美国 欧洲 德国 英国 意大利 瑞士 俄罗斯 加拿大 澳大利亚 中国 项目 大陆反射地震探测计划COCORP 地球透镜计划EarthScope 地球探测计划EUROPROBE 大陆反射地震计划DEKORP 反射地震计划BIRPS 地壳探测计划CROP 大陆反射地震计划DEKORP 深部探测计划 岩石圈探测计划LITHOPROBE 四维地球动力学计划AGCRC 玻璃地球计划Glass-Earth 地球探测计划AuScope 深部探测技术与实验研究专项SinoProbe 表1世界各国 （地区） 的深部地球勘探计划 99 2020年第12期西部探矿工程 1简介 大地电磁法正演， 常使用数值模拟方法进行， 主流 的数值模拟方法是积分方程法、 有限元法和有限差分 法， 这几种方法各有所长， 也各有不足[6]。对这三种方 法在球坐标系下的应用情况， 2014年Kelbert等人做了 详细的对比[7]。 积分方程法是从基本电磁理论开始推导， 得到对 应的积分方程并进行求解， 其优点是剖分对象不包括 背景场， 未知变量少， 具有较高的计算精度； 但仅适合 求解稠密矩阵， 且异常体增大时， 求解难度也随之增 大， 此外还有占用内存过大的缺点。 有限元法则是通过变分方法， 将计算区域离散为 若干个内部均匀的有限元， 并生成对应的线性方程组， 进而求解得到场值， 其优点在于剖分灵活， 对异常体的 适应性较好， 能够更准确地模拟复杂地电模型； 但在处 理大型问题时存储量大、 计算速度慢。 有限差分法是结合边界条件将连续微分方程离散 为若干个差分方程， 进而求解得到场值， 其优点在于实 现方式简单、 迭代稳定、 计算效率高等， 但对异常体的 适应性较差， 不适用于求解复杂地电模型。 为求解全球电磁感应或地球大尺度的局部电磁响 应， 有不少研究者使用以上三种方法进行了数值模拟。 2球坐标系下的积分方程法 积分方程法正演模拟的基本原理如下 根据麦克斯韦方程组和积分方程理论， 结合电磁 张量格林函数， 得到三维大地中的电磁响应函数； E r Ep r Δσ∫ VG E r,r′ E r′ dV′（1） H r Hp r Δσ∫ VG H r,r′ E r′ dV′（2） 式中E r、H rr处的总电场和总磁场； Ep r、Hp r以天然平面电磁波为场源， 在 大地中产生的一次电场和一次磁场； Δσ三维异常体与大地电导率差值； G E r,r′、G H r,r′电并矢格林函数和磁并 矢格林函数； E r′异常体内电场。 在数值求解阶段， 对异常体进行网格剖分， 分成若 干个单元， 并假设每个单元内部的电阻率分布是均匀 的， 在此假设下， 求解各单元网格的电场， 则各网格单 元内部的电场可近似为 ErmEprm∑ n1 N σn-σ1∫ VG E r,r′ dV′En（3） 转换为矩阵形式为 [ ]L [ ]E -[]Ep（4） 求解该方程， 即可得到异常体内的电场分布。 根据电场分布， 结合电磁响应函数 （1） 、（2） ， 即可 得到区域内的电磁场分布[8]。推导过程可参考Ting和 Hohmann的文章[9-10]。 为了不同的目的， 研究者们相继使用积分方程法 进行了球坐标系下的大地电磁正演模拟研究。Koya- ma为了研究D″层 （横向不均匀性在地幔中最明显的地 层） 的横向电导率异质性对核心成因电磁场波动的可 能影响Koyama et al., 2002； Kuvshinov则是为了研究 海洋全球电磁感应模拟过程中的影响AV Kuvshinov, 2008； 为了研究有限电导的不均匀薄片覆盖的层状球 形地球上的频域电磁 （EM） 感应问题， Sun与 Egbert也 进行了球坐标系下的积分方程法大地电磁法正演模拟 J Sun et al., 2012。他们都采用了基于带有收缩核的 收缩积分方程 （contracting integral equation， CIE） 进行 求解。 Fainberg 和Zinger于1980年首次提出收缩积分方 程法， 也称为迭代耗散法 （iterative dissipative ， IDM）， 并给出了相应的推导过程Fainberg et al., 1980。在此基础上， Fainberg 于1990年使用IDM计算 了球形地球模型中的电磁场Fainberg et al., 1990a, 1990b。随后研究人员对收缩积分方程法进行了改 进， 使得它的计算效率和内存消耗大幅减少A Kuvshi- nov et al., 2012。其中应用较广的有将改进的Neu- mann 级数 （modified Neumann series， MNS） 与收缩积 分方程法相结合以提高正演效率AV Kuvshinov et al., 1999。Avdeev采用Krylov子空间迭代法来代替 收缩积分方程法求解过程中的 Neumann 级数求和 Avdeev, 2000。Singer 的工作显示采用 Krylov 子空 间迭代法可以有效减少迭代次数B. Singer, 2008， 提 高正演效率。 Sun与 Egbert利用球坐标系下的对称性， 提出采 用FFT变换与积分方程相结合， 并采用Krylov子空间 迭代法进行求解J Sun et al., 2012。为了进一步提高 效率， Sun与Kuvshinov提出采用SVD分解的收缩积分 方程法， 使得积分方程法的效率大大提高J Sun et al., 2015。Chen等人采用嵌套式细网格的方式， 在大区域 用粗网格， 而在局部感兴趣区域采用细网格进行积分 方程法正演， 这样就大大提高了正演效率， 同时又保证 了局部异常信息Chen et al., 2020。该方法大大提升 100 2020年第12期西部探矿工程 了积分方程法的适用范围， 进一步提高了积分方程法 的实用性。 3球坐标系下的有限元法 有限元法正演模拟的基本原理如下 根据麦克斯韦方程组， 结合辅助方程， 选取时谐因 子为e-iωt， 得到背景场所满足的表达式 ∇∇Epk2Ep-iωμJf（5） 式中Ep背景场的电场强度。 进而得到二次场的表达式 ∇∇Es--iωμσEs -iωΔσEp（6） 式中Es二次场的电场强度； Δσ异常体电导率与背景电导率差值。 结合第一类边界条件， 利用变分原理， 得到对应的 泛函 F[Es] 1 2∫V[∇Es∇Es-k 2E sEs]dV （7） 对求解区域进行网格剖分， 将其剖分成若干个长 方体， 则在每个单元内， 电场满足 Ee∑ j1 12 N e jE e j∑ i1 4 N e xiE e xix ∑ i1 4 N e yiE e yiy ∑ i1 4 N e ziE e ziz （8） 将 （7） 式离散化， 得到每一个单元内部的系数矩 阵， 汇总为如下形式 KEB（9） 求解该方程即可得到区域内每个单元的电场。利 用求得的电场值， 代入麦克斯韦方程组即可计算磁场 Domnikov, 2018; 童孝忠 et al., 2009。 有限元法在工程领域的应用十分广泛， 由于其网 格剖分要求比较高， 通常需要专业的网格剖分软件进 行剖分。因此相关的软件也很丰富， 目前已有许多软 件可以用来进行电磁场有限元数值模拟， 如ALBER- TA、 ANSYS、 COMSOL MULTIPHYSICS、 FlexPDE 等陈小燕, 2014。这些软件大大方便了有限元的地 球物理正演， 目前有研究结合以上软件开展球坐标系 下的电磁正演模拟工作。Ribaudo基于FlexPDE的基 础上， 开展全球地磁感应问题的二维和三维模拟， 认为 使用该技术可以模拟几乎任意的电导率和背景场， 并 且很容易适应地球自转， 可以通过地空界面的边界 条件施加横向变化的表面电导， 进而总结出FlexPDE 的 优 势 Ribaudo et al., 2012。 Haviland 等 人 利 用 COMSOL等工具开展月球的有限元电磁正演模拟研 究， 实现了月球上的大尺度电磁正演模拟Haviland et al., 2019。Grayver 等人利用基于正交网格的有限元 法实现了基于球坐标系下的大地电磁正演， 在这基础 上研究了不同投影对大地电磁正演结果的影响Gray- ver et al., 2019。 4球坐标系下的有限差分法 大地电磁有限差分正演主要使用交错网格有限差 分法， 它是一种结构化网格数值方法， 通常被认为是大 地电磁有限差分正演模拟方法中最简单的方法之一 李建平 et al., 2018， 其基本原理如下Mackie et al., 1994 首先对求解区域进行交错采样网格化， 将其分割 为交错的电场网格和磁场网格， 在区域内部， 二者的网 格节点互为彼此的网格中心； 然后结合交错网格， 对麦 克斯韦方程组的积分形式进行离散化， 得到电场与磁 场各分量间的离散关系， 形如 |Ex0-Ex0|Δx ˉ-|Ez1-Ez0|Δz ˉiωμ0Hy0Δx ˉΔz ˉ 10 |Hz1-Hz0|Δz-|Ey1-Ey0|Δy 2 ρ0ρ1 Ex0Δx ˉΔz ˉ 11 进一步合并方程， 即可得到与单个网格上某一磁 场分量与其周围相邻12个磁场分量的关系式， 形如 C0H0C1Hx1C2Hx2C3Hx3C4Hx4 C5Hy1C6Hy2C7Hy3C8Hy4 C9Hz1C10Hz2C11Hz3C12Hz40 12 汇总全部网格的磁场离散方程， 则得到大地电磁 正演模拟需要求解的关于磁场的离散方程组， 求解即 可得到区域内的磁场分布； 将其代入麦克斯韦方程组， 即可求得区域内的电场分布陈辉 et al., 2011。 有限差分法由于其方法原理简单， 算法实现容易， 计算效率高， 被广泛应用于全球电磁正演模拟， 目前已 经成为实用化大尺度电磁勘探数据解释的主要数值方 法之一。交错网格有限差分法由于它的网格定义很好 地表达了电磁感应空间关系， 在Mackie等人将该方法 应用到大地电磁后， 该方法就成为大地电磁有限差分 的标准方法。Uyeshima与Schultz把Mackie等的笛卡 尔坐标系下的交错网格有限差分大地电磁数值模拟方 法扩展到求解球坐标系下全球尺度的电磁感应问题 Uyeshima et al., 2000。为了使大地电磁法正演模拟 更具通用性， 研究人员提出了基于模块化的基于球坐 标系的大地电磁交错网格有限差分大地电磁正演， 目 前该部分工作已经包含在 ModEM 软件中Egbert et al., 2012。由于该算法实现简单， 计算效率高， 已经广 泛应用于实际的深部探测数据的正反演解释中罗威 101 2020年第12期西部探矿工程 et al., 2019。 5结论与展望 球坐标系下的大地电磁正演模拟研究， 在三大数 值模拟方法上已经各有进展。其中有限元法已经催生 出一大批较为成熟的软件， 但是由于全球尺度的大地 电磁正反演解释涉及计算区域大， 计算量巨大， 针对球 坐标系下 （全球尺度） 的大地电磁有限元正演模拟的研 究较少； 而积分方程法和有限差分法由于具有相对高 效， 在球坐标系下 （全球尺度） 的大地电磁正演模拟中 得到了广泛的应用。 随着资源勘探开发的深入， 全球尺度的地球物理 勘探研究得到了越来越多的重视， 相信在同行们的努 力下， 这一技术在精度、 速度、 适应性等方面将有所突 破。人工智能在各个科学研究领域都得到了广泛的应 用， 相信不久的将来， 机器学习也将在球坐标系下的大 地电磁正、 反演中得到广泛的应用。 参考文献 [1]滕吉文,等.中国地球动力学研究的方向和任务[J].岩石学报, 2010,2611. [2]Srivastava S.Theory of the Magnetotelluric for a Spheri- cal Conductor[J]. Geophysical Journal International,1966,114 373-387. [3]Weidelt P.The Inverse Problem of Geomagnetic Induction[J]. Geophysical Journal International,1973,351-3379. [4]Schmucker U.Substitute Conductors for Electromagnetic Re- sponse Estimates[J].Pure and Applied Geophysics,1987,1252- 3341-367. [5]Grayver AV, et al.Three-dimensional Magnetotelluric Model- ling in Spherical Earth[J]. Geophysical Journal International, 2019,2171532-557. [6]徐凌华,等.基于有限单元法的二维∕三维大地电磁正演模拟 策略[J].物探化探计算技术， 200955,19-23. [7]Kelbert A, et al.Global 3-D Electromagnetic Forward Model- lingABenchmark Study[J]. Geophysical Journal International, 201422. [8]任政勇,等.一种新的三维大地电磁积分方程正演方法[J].地 球物理学报,2017,60114506-4515. [9]Hohmann GW.Three- dimensional Induced Polarization and Electromagnetic Modeling[J]. Geophysics,1975,402309-324. [10]Ting SC,Hohmann GW.Integral Equation Modeling of Three- dimensional Magnetotelluric Response[J].Geophysics,1981,46 2182-197. （上接第98页） （1） 构造等深度图的结果中可以发现， 该勘探区在 上、 下两部分区域存在地势突变较大的区域， 说明该区 域经历过大规模的构造运动。 （2） 从深度剖面图中的结果中识别出地堑及地垒 地堑复合型断层， 可以看出该勘探区发育大量反转构 造,背斜发育,面积大。在B、 C两区域存在圈闭构造， 在穿过C区域的剖面上有明显的强反射层的出现， 可 能是油气或者油水的分界面， 且在剖面图上可以看出 发育有良好的盖层， 所以推断在B、 C两块区域发育有 油气。通过两幅剖面图的对比， 可以看出剖面2的质量 更好。这主要是以下两个原因造成的 第一， 测井资料 的质量存在差异； 第二， 在进行剖面的时深转换时， 测 井的位置距离剖面的位置不一致， 有近有远， 从而导致 测井资料对剖面时深关系的约束不同。 （3） 由得出的构造等深度图和深度剖面图可以看 出， 在构造等深度图上出现的圈闭以及深度等值线变 化剧烈的地方可以较好的与深度剖面图上出现的断层 等构造的位置相吻合， 所以根据提取出的勘探区的测 井数据， 拟合的时深公式是准确可靠的。因此， 本文中 采用的时深转换的方法是切实可行并且行之有效的。 参考文献 [1]张银国.东海西湖凹陷花港组油气地质条件与油气分布规 律[J].石油实验地质,2010,323223-226. [2]姜文斌.东海西湖凹陷中央反转构造带成藏特征研究[D].中 国地质大学北京,2012. [3]徐国盛,等.西湖凹陷中央反转构造带花港组致密砂岩储层 成岩环境演变与孔隙演化[J].成都理工大学学报自科版, 2016,434385-395. [4]宋小勇.东海盆地西湖凹陷构造样式及其对油气聚集的控 制[D].中国地质大学北京,2007. [5]贾凌云,林年添,赵蕾.地震波时深转换方法初探[J].图书情报 导刊,2011,219191-193. [6]陈传仁,周熙襄.一种精确层速度反演方法[J].物探化探计算 技术,1999,213212-215. [7]胡英,姚逢昌.用于叠前深度偏移速度建模的一种新方法[J]. 石油勘探与开发,2000,27262-64. [8]刘启凤.地震资料成像和构造成图中速度的分析与应用[J]. 江汉石油职工大学学报,2004,174 52-53. [9]彭嫦姿,周惠群.湘中地区高陡度地层速度研究方法技术探 讨[J].国土资源导刊,2006S188-92. 102