大地电磁测深(MT)是一种利用天然交变电磁场研究地球电性结构的地球物理勘探方法,其探测范围可从近地表十几米向下至壳幔过渡带.大地电磁数据解释通常假设平面电磁波场源垂直入射水平方向无限延伸的地表,研究区域可利用笛卡尔坐标系近似地表示(Berdichevsky and Dmitriev, 2010; Chave and Jones, 2012).以此为前提,大地电磁正反演经历了由二维到三维的快速发展,从数十个站点组成的单一测线扩展到包含数百个测点的MT阵列,越来越多的大地电磁数据使得该方法不再局限于小范围的自然资源调查和地质结构解释.
MT方法在不同尺度的地球内部结构和动力学研究中得到了广泛应用(Zhang et al., 2016; Ye et al., 2018).近十年来,一些横跨大陆规模的MT阵列调查项目在全世界范围内开展.美国的USArray、澳大利亚的Aus-Lamp及中国的SinoProbe(魏文博等, 2010; 杨文采等, 2011)等综合地球深部探测计划都专门部署有大规模的MT阵列,旨在揭示大陆岩石圈的组成、结构与演化,获取地球深部的物理性质(Dong et al., 2016; Meqbel et al., 2014; Robertson et al., 2016; Sedikou et al., 2018; Yang et al., 2015).一旦这些项目全部完成,整体研究区域将横跨数千公里,并且考虑到区域外延以满足边界条件这一实际情况,模型空间的水平距离将远大于地球半径.因此这些大陆尺度的MT调查带来了一个新的挑战:当地球的曲率不可忽视时,传统笛卡尔坐标是否仍然适用?
Srivastava(1966)分别计算了地球平面和球面状态下的阻抗和相位,结果显示地球弧度只对周期大于一天的响应值有影响.Weidelt(1972)给出了不同坐标系下MT响应函数相互转换的关系式,根据此转换公式,可以计算球形弧度的影响(Berdichevsky and Dmitriev, 2002).罗威等(2012)推导了基于球体层状介质模型的大地电磁正演公式,并计算了若干理论模型.然而,上述公式推导及计算都基于一维层状模型,地球真实的三维电性结构必定会导致更加复杂的情况.Grayver(2019)利用高阶有限元数值方法计算了三维球体模型的MT阻抗张量,并提出了适用于球体模型的响应函数.
此外,大地电磁理论建立于平面波垂直入射这一假设之上,其场源为两个均匀平面波的线性组合,可解耦成两种极化模式.该假设成立的前提条件是模型空间足够小,场源的空间变化可以忽略.为获得大尺度下精确的大地电磁三维正演结果,本文发展了一种基于球坐标系的大地电磁交错网格有限差分三维正演算法,并且利用极向-环向分解方法构建了场源模型,采用P10球谐函数使模拟的场源既满足响应函数定义又与电离层场源结构类似.首先利用经纬度信息建立三维地电模型,将构建的场源位于模型空间的正上方,然后通过直接求解球坐标系下麦克斯韦方程来获得大地电磁响应.通过与一维模型的大地电磁响应对比验证了本文算法的正确性.在此基础上,利用多个三维地电模型的对比结果分析了不同坐标系下大地电磁响应差异的影响因素.
1 球坐标大地电磁三维正演理论 1.1 数值模拟本研究使用的正演方法为基于交错网格剖分的有限差分算法(Egbert and Kelbert, 2012; Kelbert et al., 2014; 董浩等, 2014),该方法具有计算精度和效率高等优点,是电磁法领域采用较多的数值模拟手段.在大地电磁测深频率范围内,通常忽略位移电流,麦克斯韦方程组的积分形式为(取时谐因子eiωt):
(1) |
(2) |
(3) |
其中,E, H, J分别为电场强度(V·m-1)、磁场强度(A·m-1)和电流密度(A·m-2),μ为真空磁导率(H·m-1),σ为介质电导率,ω为角频率,dS为由线元dl包围的面元,与电流密度J方向垂直.为了计算上述电场或磁场,需将麦克斯韦方程组离散化.如图 1a所示,按照经度ϕ、纬度θ、深度-r(与地球半径r方向相反)将整体模型空间包括上覆空气层离散为M=Nθ×Nϕ×Nr个曲面六面体.根据交错网格设置,定义电场分布于每个单元的边缘中点,磁场分布于面元中点,以该方式互相匹配的电磁场能最大限度满足能量守恒定律(Mitsuhata and Uchida, 2004).
与笛卡尔坐标系下剖分得到的规则六面体单元不同,曲面单元的边长、面积和体积会随着深度和纬度而变化,因此每个单元都需要单独计算它们的几何参数(Uyeshima and Schultz, 2000).如图 1b所示,以中心电导率为σ(i, j, k)的曲面单元为例,其在纬度、经度和深度方向的边长计算公式为
(4) |
(5) |
(6) |
其中i=1,2,3,…,Nθ,j=1,2,3,…,Nϕ,k=1,2,3,…,Nr,r(k)为地球中心到网格第k层的距离,定义Sθ(i, j, k)为垂直于θ方向的面元面积:
(7) |
同理,可得另外两个面元面积:
(8) |
(9) |
将式(4)—(9)代入麦克斯韦方程组,取顺时针为正方向,则式(1)(2)数值离散形式为
(10) |
(11) |
其中指数i+表示i与i+1的中点,i-为i-1与i的中点,以此类推.l_σ(i, j, k)为体电导率在单元边缘lθ(i, j, k)的加权值,加权因子为环绕该边缘的四个单元体积,以图 1b所示单元为例,其体积计算公式为
(12) |
其中r1=[r(k)+r(k+1)]/2, θ1=[θ(i+1)+θ(i)]/2.将式(11)代入式(10)消去磁场Hϕ(i+, j, k+).同理根据式(2)可得到Hϕ(i+, j, k-),H-r(i+, j+, k),H-r(i+, j-, k)关于电场E的表达式:
(13) |
(14) |
(15) |
将式(11),(13)—(15)代入式(10)得到仅关于电场的二阶差分方程式.对于某一指定角频率ω,上述问题转换为求解线性方程组:
(16) |
展开矩阵形式为:
(17) |
A为3M × 3M大型稀疏矩阵,E为需要求解的电场3M阶未知量,B则为边界条件的3M阶矢量.边界条件由一维模型计算得到,初始条件由二维模拟计算获得.
因空气与地下介质电导率存在巨大差异,对式(16)迭代求解时(Smith, 1996),可能导致矩阵A产生严重的病态,本文算法采用一阶不完备LU分解进行矩阵预优以避免此问题.另外,大地电磁测深理论显式地要求电场散度为0,即求解的电场需满足▽·E=0,假设第n次迭代电场为En,其散度ψ=▽·En,ψ与电势φe存在以下关系(Mackie et al., 1994; Smith, 1996):
(18) |
则电场散度校正公式为:
(19) |
电场校正并不会影响磁场,关系式▽×(En-▽ψ)=▽×E=iωμH始终满足.
在测点处对求解的电磁场进行线性插值,水平电场Eh和水平磁场Hh有如下关系:
(20) |
阻抗张量Z是大地电磁方法的基本测量参数.进一步计算可得视电阻率ρ和相位Φ:
(21) |
(22) |
将地球及其上覆空气层划分为球体层状电导率模型,每一层电场的散度为零,则其极向-环向分解形式为
(23) |
其中P是极向电势,T为环向电势,分别对应不同模式的电流源,
(24) |
(25) |
其中Ψ代表P或者T,
(26) |
τl(λr)和ξl(λr)分别为第一类和第三类Riccati贝塞尔函数,与第一类和第三类球形贝塞尔函数有关.
对于上覆空气层介质,其满足准静态条件λr << 1,贝塞尔函数τl(λr)和ξl(λr)可分别退化为常数因子rl+1和r-l,则Ψl(r)可简化为
(27) |
Sun和Egbert (2012)通过电磁场球谐分解证明外部场源主要与环向电势T相关,极向电势P的影响可以忽略.将式(27)代入式(24),消去与极向电势P相关的因式项,得到准静态条件下电场表达式:
(28) |
空气层中环向系数αlT,βlT可以通过递归方法求得.理论上式(28)可采用任意阶(l≥1)球谐函数,但是P10参数设置(l=1, m=0)构建的模型最接近静日状态下电离层场源结构(Berdichevsky and Dmitriev, 2002; Grayver et al., 2019).图 2显示了模拟的场源分布情况,图 2a、2b分别为式(28)计算结果的平面和球面示意图,对结果进行归一化处理,使得其值在0~1之间变化,作为东西方向极化模式的场源.将球谐函数系数进行旋转得到的结果如图 2c、2d所示,作为南北方向极化模式的场源.
考虑到程序的通用性及便捷性,本文在ModEM框架(Kelbert et al., 2014)基础上实现了上述所发展的球坐标系大地电磁三维正演算法.为验证该算法的有效性与准确性,本文分别计算了10 Ωm,100 Ωm及1000 Ωm均匀半空间模型在球坐标和笛卡尔坐标的响应值.前人研究表明对于简单的一维模型在周期小于104 s范围内地球本身的弧度可以忽略(Berdichevsky and Dmitriev, 2002; Grayver et al., 2019),因此不同坐标系下响应值的差异十分微小.计算结果如图 3a所示, 除了1000 Ωm模型的响应曲线在104~105s呈现明显的不同,其他模型的视电阻率和相位角差异不超过1%.图 3b为全球一维层状电导率模型的响应曲线,该层状模型来源于磁层卫星和潮汐磁信号联合反演(Grayver et al., 2017),结果显示二者拟合情况较好,大部分误差在0.5%以内.同时,本文响应曲线与利用球坐标高阶有限元计算的结果吻合(Grayver et al., 2019),进一步证明了本文算法的正确性.三维模型的情况相对于一维会更加复杂,除了电导率的三维特性还存在地理坐标与笛卡尔转换带来的畸变问题,下一节将对三维球体电导模型进行模拟及对比研究.
本节利用所开发的球坐标大地电磁正演算法对大尺度的三维电导模型进行模拟计算.如图 4所示的三维模型包含两个位于中心,东西排列的异常体,每个异常体尺寸为30°×20°×20 km,电阻率分别为1000 Ωm和10 Ωm,周围为100 Ωm均匀半空间,测点以2°为间隔均匀分布于异常体上方.对该模型分别以1°,0.5°和0.25°为间隔进行剖分得到三种不同分辨率的模型.
大地电磁数据处理流程首先是按照某一投影方法将研究区域转换到笛卡尔坐标系,然后在此坐标系下进行后续的反演解释工作.在电磁勘探领域,研究人员或采用成熟的商业软件或自己编写算法进行投影,并没有统一的规则.已公开发表的地图投影方法大概有上百种,按照其投影方式可大致分为圆柱投影、圆锥投影和方位投影三大类,其中每一类又可细分为等距离、等面积和等角投影.一般来说采用何种投影方式需要考虑研究区域位置、形状、投影用途等各种因素,恰当的投影方法可以很好地维持球体模型的某一特性,但不论何种都不可避免会产生畸变(Snyder, 1987).
大地电磁数据处理中常用的投影方法为横轴墨卡托投影(Transverse Mercator projection),又称高斯-克吕格投影(Gauss-Kruger projection),属于等角横切椭圆柱类型.该投影将地球从0°子午线起自西向东划分为多个6度或3度带,这种投影方法适用于小尺度勘探区域,可以最大程度保证投影的真实性.本节三维模型位于北半球中纬度区域,远离赤道和北极点,主要呈现东西向延伸形态(图 4).为了保持异常体的形态特征,本文采用兰勃特等角圆锥投影(Lambert Conformal Conic),其基本思想为以正圆锥切于或割于球面,应用等角条件将地球面投影到圆锥面上,然后沿某一经线展开成平面.投影后纬线为同心圆弧,经线为同心圆半径,没有角度变形,兰勃特投影常用于我国1:100万地形图绘制.
从球体模型转换到笛卡尔模型分为两个步骤:第一步经纬度网投影到公里网,第二步反向投影插值获得公里网的电导率值,二者组成笛卡尔三维电导模型.本文采用北纬30°和60°为标准纬线(即地球与投影圆锥相割的两条纬线,标准纬线没有任何变形,离开标准纬线愈远,长度和面积变形就越大),模型中心坐标(45°N,120°E)为原点投影得到笛卡尔三维电导率模型(如图 5所示,该模型对应分辨率为0.5°×0.5°球体模型).可以看到投影后的笛卡尔模型保持了异常体的形态,体现了等角类投影方法的特质,但牺牲了异常体面积和距离的准确性.表 1给出了异常体真实面积和不同分辨率投影后的面积.位于地表的测点按照相同的投影方法转换到笛卡尔坐标,使异常体与测点相对位置保持不变.
分别对球坐标和笛卡尔模型进行三维大地电磁正演模拟,计算平台为曙光高性能刀片服务器,计算方式为按频点并行,表 1列出不同分辨率模型的正演时间,球体模型耗时大约是笛卡尔模型的两倍.图 6显示了不同极化模式下计算时间随周期变化曲线,整体表现为球坐标系的正演计算所耗费的时间要长于笛卡尔坐标系.分析其原因主要有以下两点:(1)球坐标下每个单元的几何参数都需要单独计算;(2)在球坐标体系下网格单元的边长、面积随纬度和深度变化,由它们组成的系数矩阵要比笛卡尔坐标系的系数矩阵复杂,因此在求解线性方程组和散度改正时需要更多迭代次数.
为了对比的一致性,我们首先将笛卡尔模型正演的电场反向投影回球坐标,然后在同一坐标下与球体模型电场进行比较.图 7给出了不同分辨率的模型在周期1000 s时电场差异百分比E_diff,其在电性横向变化的边界处有明显的波动.Eθ_diff从低纬度到高纬度逐渐变大,反映了越往极点靠近,因坐标体系不同导致的模拟差异愈明显.Eϕ_diff背景值在20%左右,在进行投影时我们设置原点在研究区中点,因此Eϕ_diff在中央经线处呈现两端大中间小的分布.
图 8显示了地表测点处视电阻率的差异(取以10为底的对数),其分布与E_diff相关,位于边界周围的测点误差较大.值得注意的是,随着模型分辨率增大即网格剖分愈细致,视电阻率差异会增大.视电阻率差异与网格和电导率这两个因素相关,当模型的分辨率增大,即网格剖分越细致,不同坐标系下网格单元的几何差异会减小,电导率差异会增大,视电阻率计算的是单点的值,因此在这种情况下,电导率的影响占主导地位.另外,随着周期的增大,视电阻率的差异呈现先减小后增大的趋势.长周期通常反映深部信息,如表 1所示,兰勃特投影使得笛卡尔异常体表面积缩小,随着深度的增加,球坐标异常体的横向面积会持续减少,在这过程中与笛卡尔模型的差异会逐渐减小然后又增大,是导致上述变化趋势的原因.
一般来说,大陆尺度的模型通常会包含部分海洋信息,我们利用ETOPE2海洋测深数据建立了中国大陆及邻海电导模型(图 9).陆地电导率设置为100 Ωm,海洋部分取海水平均值为0.2 Ωm,模型中心区域按照0.5°×0.5°精度进行网格划分,测点以2°×2°间隔均匀分布于陆地上.设定兰勃特等角圆锥投影标准纬线分别为26°N和42°N,原点为(34°N, 106°E)对模型和测点分别投影得到匹配的笛卡尔电导模型.如图 10所示,投影失真导致笛卡尔模型显示出更多的海域,例如南海、孟加拉湾和日本海面积都有所增加.
我们取5个频点进行并行计算,球坐标模型耗时94 min全部达到收敛,笛卡尔模型耗时69 min收敛.同样地,我们对电场和视电阻率进行比较,如图 11、图 12所示,整体来说电场差异比视电阻率差异更加明显,图 11中Eθ_diff的分布具有从低纬度到高纬度增大的趋势,而Eϕ_diff大部分区域保持固定的20%误差,沿海地带模型电导率变化致使电场差异值波动,视电阻率在沿海和东北地区显示出明显的不同.本文设置中国大陆为均匀半空间,虽然被海洋包围但是内陆三维性不甚显著,如果考虑真实地下电导模型,不论是电场还是视电阻率差异都会明显增大.
本文基于交错网格有限差分方法,推导并实现了球坐标系的大地电磁三维正演计算,利用极向-环向分解构建了球坐标场源模型,该模型符合电离层场源空间结构同时也满足响应函数定义,可取代传统平面波场源.与笛卡尔计算结果的对比表明,对于一维模型,地球弧度对MT常用频段的影响很小,与前人理论吻合从而证明了算法的正确性.对于三维模型,不同坐标体系计算的电场和视电阻率存在明显的差异,主要由电导率的三维特性和投影畸变及二者的耦合效应导致的.因此,在进行大尺度MT数据处理和解释时,应优先考虑采用球坐标系.
根据本文三维球体模型的位置、形状等特点,我们采用兰勃特等角圆锥投影进行模型转换,其他投影方法可能会得到不同的对比结果,所以当忽略地球弧度进行模拟时,需要依据实际情况谨慎地选择合适的投影方法以及投影参数.但无论何种投影,笛卡尔模型畸变始终存在,本文提出的球坐标算法的优势即在于可以跳过这一步骤,直接根据研究区经纬度信息构建模型并进行后续的数据处理和计算工作.
球体模型因其网格的复杂性,其正演时间约是笛卡尔的两倍,下一步的研究工作主要为提高球坐标三维正演算法的计算效率,在此基础上,最终实现球坐标大地电磁三维反演解释.
致谢 感谢评审专家对本文提出的修改意见,感谢美国地质调查局Anna Kelbert研究员对模拟工作的指导.
Berdichevsky M N, Dmitriev V I. 2002. Magnetotellurics in the context of the theory of Ill-posed problems. Society of Exploration Geophysicists. |
Berdichevsky M N, Dmitriev V I. 2010. Models and Methods of Magnetotellurics. Berlin: Springer Science & Business.
|
Chave A D, Jones A G. 2012. The Magnetotelluric Method:Theory and Practice. Cambridge: Cambridge University Press.
|
Dong H, Wei W B, Ye G F, et al. 2014. Study of three-dimensional magnetotelluric inversion including surface topography based on finite-difference method. Chinese Journal of Geophysics (in Chinese), 57(3): 939-952. DOI:10.6038/cjg20140323 |
Dong H, Wei W B, Jin S, et al. 2016. Extensional extrusion:Insights into south-eastward expansion of Tibetan Plateau from magnetotelluric array data. Earth and Planetary Science Letters, 454: 78-85. DOI:10.1016/j.epsl.2016.07.043 |
Egbert G D, Kelbert A. 2012. Computational recipes for electromagnetic inverse problems. Geophysical Journal International, 189(1): 251-267. DOI:10.1111/j.1365-246X.2011.05347.x |
Grayver A V, Munch F D, Kuvshinov A V, et al. 2017. Joint inversion of satellite-detected tidal and magnetospheric signals constrains electrical conductivity and water content of the upper mantle and transition zone. Geophysical Research Letters, 44(12): 6074-6081. DOI:10.1002/2017GL073446 |
Grayver A V, Van Driel M, Kuvshinov A V. 2019. Three-dimensional magnetotelluric modelling in spherical Earth. Geophysical Journal International, 217(1): 532-557. DOI:10.1093/gji/ggz030 |
Kelbert A, Meqbel N, Egbert G D, et al. 2014. ModEM:A modular system for inversion of electromagnetic geophysical data. Computers & Geosciences, 66: 40-53. |
Luo W, Wang X B, Qin Q Y. 2012. MT 1D forward theory based on sphere layered model. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 34(4): 384-389. |
Mackie R L, Smith J T, Madden T R. 1994. Three-dimensional electromagnetic modeling using finite difference equations:The magnetotelluric example. Radio Science, 29(4): 923-935. DOI:10.1029/94RS00326 |
Meqbel N M, Egbert G D, Wannamaker P E, et al. 2014. Deep electrical resistivity structure of the northwestern U. S. derived from 3-D inversion of USArray magnetotelluric data. Earth and Planetary Science Letters, 402: 290-304. |
Mitsuhata Y, Uchida T. 2004. 3D magnetotelluric modeling using the T-Ω finite-element method. Geophysics, 69(1): 108-119. |
Robertson K, Heinson G, Thiel S. 2016. Lithospheric reworking at the Proterozoic-Phanerozoic transition of Australia imaged using AusLAMP magnetotelluric data. Earth and Planetary Science Letters, 452: 27-35. DOI:10.1016/j.epsl.2016.07.036 |
Sedikou B O D, Wei W, Ye G, et al. 2018. Experiments of magnetotelluric observation network on North China and lithospheric conductivity structure from fast imaging method. Chinese Journal of Geophysics (in Chinese), 61(8): 2508-2524. DOI:10.6038/cjg2018L0695 |
Smith J T. 1996. Conservative modeling of 3-D electromagnetic fields, Part Ⅱ:Biconjugate gradient solution and an accelerator. Geophysics, 61(5): 1319-1324. DOI:10.1190/1.1444055 |
Snyder J P. 1987. Map Projections: a working manual. United States Government Printing Office, Washington.
|
Srivastava S P. 1966. Theory of the magnetotelluric method for a spherical conductor. Geophysical Journal International, 11(4): 373-387. DOI:10.1111/j.1365-246X.1966.tb03090.x |
Sun J, Egbert G D. 2012. Spherical decomposition of electromagnetic fields generated by quasi-static currents. GEM-International Journal on Geomathematics, 3(2): 279-295. DOI:10.1007/s13137-012-0039-0 |
Uyeshima M, Schultz A. 2000. Geoelectromagnetic induction in a heterogeneous sphere:a new three-dimensional forward solver using a conservative staggered-grid finite difference method. Geophysical Journal International, 140(3): 636-650. DOI:10.1046/j.1365-246X.2000.00051.x |
Wei W B, Jin S, Ye G F, et al. 2010. On the Conductive Structure of Chinese Continental Lithosphere-Experiment on "Standard Monitoring Network" of Continental EM Parameters (SinoProbe-01). Acta Geologica Sinica (in Chinese), 84(6): 788-800. |
Weidelt P. 1972. The inverse problem of geomagnetic induction. Zeitschriftfiir Geophysik, 38: 257-289. |
Yang B, Egbert G D, Kelbert A, et al. 2015. Three-dimensional electrical resistivity of the north-central USA from EarthScope long period magnetotelluric data. Earth and Planetary Science Letters, 422: 87-93. DOI:10.1016/j.epsl.2015.04.006 |
Yang W C, Wei W B, Jin S, et al. 2011. Experimental study of the continental standard grid of electromagnetic parameters:an introduction to project sinoprobe-01. Acta Geoscientica Sinica (in Chinese), 32(S1): 24-33. |
Ye T, Huang Q H, Chen X B, et al. 2018. Magma chamber and crustal channel flow structures in the Tengchong volcano area from 3-D MT inversion at the intracontinental block boundary southeast of the Tibetan Plateau. Journal of Geophysical Research:Solid Earth, 123(12): 11112-11126. DOI:10.1029/2018JB015936 |
Zhang H Q, Huang Q H, Zhao G Z, et al. 2016. Three-dimensional conductivity model of crust and uppermost mantle at the northern Trans North China Orogen:Evidence for a mantle source of Datong volcanoes. Earth and Planetary Science Letters, 453: 182-192. DOI:10.1016/j.epsl.2016.08.025 |
董浩, 魏文博, 叶高峰, 等. 2014. 基于有限差分正演的带地形三维大地电磁反演方法. 地球物理学报, 57(3): 939-952. DOI:10.6038/cjg20140323 |
罗威, 王绪本, 覃庆炎. 2012. 基于球体层状介质模型的大地电磁正演. 物探化探计算技术, 34(4): 384-389. DOI:10.3969/j.issn.1001-1749.2012.04.03 |
Sedikou B O D, 魏文博, 叶高峰, 等. 2018. 华北大地电磁测深阵列观测实验与岩石圈导电性快速成像模型. 地球物理学报, 61(06): 2508-2524. DOI:10.6038/cjg2018L0695 |
魏文博, 金胜, 叶高峰, 等. 2010. 中国大陆岩石圈导电性结构研究——大陆电磁参数"标准网"实验(SinoProbe-01). 地质学报, 84(6): 788-800. |
杨文采, 魏文博, 金胜, 等. 2011. 大陆电磁参数标准网实验研究——SinoProbe-01项目介绍. 地球学报, 32(S1): 24-33. |