射频大地电磁法(Radio-magnetotelluric)是近几年才逐步发展的一种浅地表地球物理勘探方法.RMT法的测量频段为10~250 kHz，勘探深度为数米到数十米，目前多应用于工程勘探及环境监测(Pedersen et al.，2005，2006; Tezkan et al.，2000，2005; Tezkan and Saraev，2008; Ismail et al.，2011a，2011b; Bastani et al.，2013).

由于RMT法与MT法都是在相应源的远区观测水平电磁场，因而目前在进行RMT资料处理与解释时多是直接套用MT的正反演程序(Newman et al.，2003; Linde and Pedersen，2004; Cand ansayar and Tezkan，2006，2008)，这会导致反演得到的地下电阻率值出现不切实际的极小极大值(Persson and Pedersen，2002; Kalscheuer et al.，2008)，从而得到不合理的解释结果.由于MT法测量频率在数百赫兹内，在该频段，电磁场以感应扩散为主，因此在数值模拟时通常进行准静态假设，忽略位移电流.RMT法的测量频段为10~250 kHz，在该频段内，位移电流不是远远小于传导电流，因而不可忽略位移电流，当地下岩体为结晶类高阻岩石(如石英岩等)时尤其如此(Chave and Jones，2012).因此，实现考虑位移电流的RMT正演程序具有重要意义.

Persson和Pedersen(2002)讨论了一维RMT响应中，位移电流对视电阻率和相位的影响，表明考虑位移电流得到的视电阻率和相位值均低于准静态计算结果，同时文中对层状模型进行反演，得出结论：不考虑位移电流的反演电阻率与实际偏差较大.继均匀半空间和层状模型后，Kalscheuer等(2008)首次研究了二维介质中位移电流对RMT响应的影响，说明了在RMT频段内，当地下电阻率大于10000 Ωm后，位移电流会对计算精度产生较大影响，视电阻率和相位值均低于准静态结果，同时准静态条件下的反演结果会带来虚假构造，导致错误的解释.然而，Kalscheuer等(2008)的工作均是基于平地形的有限差分正演程序，在野外实际工作中，通常面临的是复杂起伏地形，因而开发带地形的RMT程序更具实际意义.

为了解决复杂起伏地形情况下RMT的数据解释问题，本文开发了基于非结构三角形网格的2D带地形RMT有限元正演程序.首先通过一算例对比验证了程序的正确性.由于RMT模拟中存在位移电流，那么TM模式的计算也必须考虑空气层的影响，为此，文中讨论了高精度RMT模拟中空气层厚度的最佳选择.接着，文中讨论了起伏地形下位移电流对数值解的影响.最后通过一实际算例进一步 说明在RMT数值模拟中必须考虑位移电流的影响.

假设x方向为构造走向方向(图 1)，z方向垂直向下，电磁场分量只沿y、z方向有变化，因此二维构造下，电磁场方程可写为(取时间因子e^iωt)：

其中ω为角频率，σ为电导率，ε为介电常数，μ₀为磁导率.(1)—(3)式定义了TE模式，(4)—(6)式定义了TM模式.

图 1

Fig. 1

图 1 求解域Ω=Ω1∪Ω2，上边界Γ1，外边界ΓD为无穷远边界，Γint为内边界，x为二维构造走向方向 Fig. 1 Computational domain Ω=Ω1∪Ω2，Γ1 is the upper boundary，ΓD is the infinite boundary，Γint is the inner boundary，the geo-electrical strike direction is along with the x-axis

将(6)式两边求/y，(7)式两边求/z，并代入(5)式即得TE模式下电场E_x分量满足的Helmholtz方程：

同理可得TM模式下磁场H_x分量满足的Helmholtz方程：

(7)式(8)式可统一表示为如下所示的偏微分方程：

对于TE模式，(9)式中的各参数代表：τ= 1 /iωμ，λ=-(σ+iωε)，u=E_x.

对于TM模式，(9)式中的各参数代表：τ= 1 /σ+iωε，λ=-iωμ，u=H_x.

在低频情况下，由(4)、(5)式可知，磁场在空气中几乎不变化，因此对于低频问题，如MT问题，TM模式一般不考虑空气.对于RMT来说，由于位移电流在空气中占主导，因此由(4)、(5)式可知，磁场在空气中存在变化，因此RMT问题的TE和TM模式，均要考虑空气层.

本文中计算区域如图 1所示，外边界(Γ_D)上电位通过层状介质解析解给出，上边界(Γ₁)取u=1，内边界(Γ_int)根据场切向分量连续给出.综上，边值问题中的边界条件如下：

网格剖分可分为两种类型，结构化网格与非结构化网格.结构化的网格在区域加密时会导致不必要的外部节点增加，从而增加额外的计算量，且结构化网格在对复杂地形及地下构造模拟中不够灵活(Li and Spitzer，2002; Mitsuhata and Uchida，2004； 谭捍东等，2003；胡祖志等，2006；胡祥云等，2012；董浩等，2014).因而，非结构网格以其灵活性近些年广泛应用于地球物理数值模拟领域(Franke et al.，2007; Ren and Tang，2010；Ren et al.，2013；刘长生等，2010；吴小平等，2015).

本文中，求解域Ω进行网格离散时采用非结构的三角形网格(Shewchuk，1996)，以保证灵活模拟复杂地质构造及地形.在离散过程中，通过设定每个三角形单元的最小角不低于30°来保证单元质量，同时通过给定局部区域单元片面积最大值来实现网格的局部加密.进行网格剖分时，首先根据最低频率确定最大趋肤深度δ_max，剖分区域边界距离不低于10δ_max；在电阻率突变界面附近控制单元大小进行局部加密；与MT法不同，RMT法计算TM模式的磁场时，空气中的磁场不能看做常数，因此TM模式也需要考虑空气层，而空气中的网格剖分也要仔细考虑.在高频情况下，电磁场在空气中以电磁波的形式向外传播，空气中电磁波的传播速度为光速c=3×10⁸ m·s^-1，根据c=λ·f，其中λ为波长，f为频率，据此可计算出不同频率下的波长，本文中网格剖分时最短波长内保证有20~30个节点.

基于线性三角形有限单元，通过求离散边值问题，从而获得关于方程(9)中未知数u的线性方程组：

式中，U 为方程(9)中u的一维向量，刚度矩阵 A 为稀疏对称复数矩阵.本文通过三个一维数组仅存储矩阵的非零元素，从而节省计算内存.采用Intel Math Kernel Library(Intel，2010)中的直接求解器PARDISO(Schenk and Gärtner，2004)来求解线性方程组.求出TE及TM模式的电场E_x和磁场H_x后，本文采用四点差分求得相应的磁场H_y和电场E_y.最后，阻抗、视电阻率ρ_ij和相位φ_ij通过如下计算公式获得：

其中，Re、Im分别表示复数的实部与虚部.

本节首先计算了Kalscheuer等(2008)文中的模型，并进行了结果对比，验证了程序正确性；然后以Dike模型为例研究了空气层厚度对数值解的影响；接着讨论了起伏地形下，位移电流的影响规律；最后通过舒家店的实际地质模型说明了位移电流对 实际数据的影响.文中所有计算均在CPU为2.3 GHz，RAM为2 GB的个人计算机上完成.

图 2中，矩形异常体的电阻率ρ₁为1000 Ωm，背景电阻率ρ₂为10000 Ωm，相对介电常数ε_r为5.该模型的三个计算频点分别为10 kHz、100 kHz、250 kHz，计算中均考虑位移电流.区域大小为x∈ -5000，5000，z∈ -5000，5000 .图 3为网格剖分示意图，经非结构网格剖分后，总节点数为19693，单元数为39122，三个频点的TE/TM计算总耗时为6.68 s.

图 2

Fig. 2

图 2 均匀半空间中赋存一矩形异常体，半空间电阻率为1 0000 Ωm，相对介电常数为5，矩形电阻率为1000 Ωm Fig. 2 A 2D model with a conductive block of ρ1=1000 Ωm in a half-space with a resistivity of ρ2=10000 Ωm and the constant relative dielectric permittivity εr=5

图 3

Fig. 3

图 3 局部网格剖分示意图 Fig. 3 Illustration of the local grid for the 2D block model

图 4、图 5分别为TE、TM模式下计算得到的视电阻率及相位曲线，其中线条为本文有限元计算结果，符号为Kalscheuer等(2008)的有限差分程序计算结果，二者均考虑位移电流.从图中可看出本文计算结果与Kalscheuer等(2008)的结果吻合，表明程序正确.

图 4

Fig. 4

图 4 图 2所示模型考虑位移电流的TE模式(a)视电阻率及(b)相位曲线，频率分别为10 kHz、100 kHz、250 kHz线条为本文程序计算结果，符号为Kalscheuer等(2008)计算结果. Fig. 4 The apparent resistivity(a) and phase(b)curves of TE-mode for the 2D block model showed in Fig. 2. Three frequencies are calculated，which are 10 kHz，100 kHz，250 kHzCurves with line are solutions calculated with displacement current and curves with symbol are solutions calculated by Kalscheuer et al.(2008).

图 5

Fig. 5

图 5 图 2所示模型考虑位移电流的TM模式(a)视电阻率及(b)相位曲线，图例与图 4相同 Fig. 5 The apparent resistivity(a) and phase(b)curves of TM-mode for the 2D block model showed in Fig. 2. The legends are similar to Fig. 4

在RMT有限元模拟中，由于位移电流项的存在，即使TM模式也不可忽略空气层，为此，本节讨论空气层厚度对数值解的影响.图 6为Dike模型示 意图，背景电阻率ρ₂=10000 Ωm，相对介电常数为5，中间断层电阻率ρ₁=1000 Ωm.

图 6

Fig. 6

图 6 Dike模型示意图 Fig. 6 Geometry of the Dike model

表 1为不同空气层厚度下网格剖分后的单元节 点信息、计算误差及计算耗时.计算频率为250 kHz，空气层厚度从0.01λ~10λ，其中λ为波长，根据c=λ·f 计算而来，c为光速，因此λ=1200 m.从表 1中可看出，当空气层厚度大于0.2λ，TE/TM模式视电阻率及相位的最大误差均在1%以内；当空气层厚度小于0.2λ时，TE/TM模式视电阻率及相位的最大相对误差均随空气层厚度的减小而增大.因此，在数值计算中，保证空气层厚度大于0.2λ即可满足精度要求.

表 1 (Table 1)

表 1 对不同空气层厚度剖分后的网格信息、计算误差及计算量 Table 1 Mesh information，maximum relative errors and computational time of numerical solutions for different thickness of the air layer 空气层厚 (m) 节点数 单元数 误差(%) 计算 时间 (s) ρTE ρTM φTE φTM 0.01λ 72578 14384 10.34 4.54 8.40 2.13 41.49 0.05λ 700064 13958 2.72 2.23 3.40 1.85 39.07 0.1λ 71038 14155 1.66 1.65 2.08 1.50 40.32 0.2λ 71581 142725 0.84 0.73 0.89 0.92 41.29 0.25λ 71753 143076 0.77 0.66 0.82 0.49 41.48 0.5λ 72528 144616 0.60 0.60 0.52 0.59 42.23 λ 74126 147796 0.63 0.59 0.79 0.61 44.31 5λ 79710 158876 0.36 0.47 0.25 0.58 51.27 10λ 91967 182725 - - - - 66.02

表 1 对不同空气层厚度剖分后的网格信息、计算误差及计算量 Table 1 Mesh information，maximum relative errors and computational time of numerical solutions for different thickness of the air layer

图 7、图 8分别为不同空气层厚度下计算的所有测点的TE/TM模式视电阻率、相位及其相应误差.图中展示的是几个典型空气层厚度的结果.所有的误差计算均与10λ的结果进行对比，无论TE还是TM模式，当空气层厚度大于0.2λ时，TE/TM 模式的视电阻率及相位所有测点的误差均在1%以 内，而当空气层厚度为0.1λ时，TE模式视电阻率最大误差为1.66%，相位最大误差为2.08%；TM模式视电阻率最大误差为1.65%，相位最大误差为1.50%.对比TE模式所有测点的误差曲线可发现，除了局部个别测点外，几乎所有测点的误差均随着空气层厚度的减小而增大，而TM模式的误差曲线对这一规律的反映则不够明显，这是因为TM模式受横向异常体影响较大所致.

图 7

Fig. 7

图 7 不同空气层厚度下得到的TE模式(a)视电阻率、(c)相位及其相应误差(b)(d).计算频率f=250 kHz，空气层厚度分别取10λ、1λ、0.25λ、0.1λ；计算误差均与10λ的结果进行对比 Fig. 7 Plots(a) and (c)are the apparent resistivity and phase for TE-mode which calculated with different thicknesses of air layer(10λ，1λ，0.25λ，0.1λ).(b) and (d)are the relative error of apparent resistivity and phase. The frequency is 250 kHz

图 8

Fig. 8

图 8 同图 7，但为Dike模型的TM模式结果对比 Fig. 8 Similar to Fig. 7，but for the TM mode

本节采用图 9所示的模型(Wannamaker et al.，1986)讨论不同高程下位移电流的影响.图 9中共给出了3个测试模型，余弦型山脊的高程分别为100 m、300 m、600 m，电阻率均为10000 Ωm，相对 介电常数为5，计算的三个频率为10 kHz、100 kHz、250 kHz.

图 9

Fig. 9

图 9 三个模型均为余弦型山脊，山脊高程分别为100 m、300 m、600 m，三个模型的电阻率均为1 0000 Ωm，相对介电常数为5 Fig. 9 Geometry of the triangular hill models. The heights of three hill model are 100 m，300 m，600 m，respectively. The resistivity is 10000 Ωm and the relative dielectric permittivity is 5

为保证结果可靠，本节三个模型均采用较密的网格剖分，以降低计算误差.网格剖分参数和计算耗时如表 2所示.

表 2 (Table 2)

表 2 起伏地形模型网格剖分参数表 Table 2 Mesh informations for hill models 模型 节点数 单元数 计算时间(s) model-1 44116 87967 22.9 model-2 43433 86607 22.4 model-3 59222 118183 37.2

表 2 起伏地形模型网格剖分参数表 Table 2 Mesh informations for hill models

图 10为考虑位移电流和准静态情况下计算的 视电阻率及相位曲线，计算频率为10 kHz、 100 kHz、250 kHz.从不同频点的结果对比来看，10 kHz下考虑位移电流和不考虑位移电流的曲线都基本重合，说明本算例中位移电流在10 kHz的频段处影响不大，而计算频率为100 kHz、250 kHz时，考虑位移电流的视电阻率曲线和相位曲线均出现下掉现象，其中相位曲线更为严重，表明对于本算例中的模型参数而言，当频率大于100 kHz后，位移电流对视电阻率的影响不可忽略，而对相位曲线的影响从频率为20 kHz后就必须考虑.

图 10

Fig. 10

图 10 (a)(c)为model-1的TE/TM模式的视电阻率曲线，(b)(d)为model-1的TE/TM模式的相位曲线，计算频率为10 kHz、100 kHz、250 kHz各图中线条为考虑位移电流的结果，符号为准静态条件下的计算结果. Fig. 10 Plots(a) and (c)are the apparent resistivity curves of TE and TM mode for model-1，(b) and (d)are phase curves of TE and TM mode for model-1.Three frequencies are calculated，which are 10 kHz，1 00 kHz，250 kHzCurves with line are solutions calculated with displacement currents and curves with symbol are solutions calculated without displacement currents.

图 11是频率为250 kHz下考虑位移电流和准静态条件下计算得到的视电阻率及相位绝对误差，将三个模型计算结果对比，发现当山脊高程越大，坡度越陡，位移电流的影响就越大.从图 11a，11c可看出，TE模式的视电阻率及相位的最大误差发生在山脊两侧，而山脊顶部则出现局部极大值；而图 11b，11d显示TM模式视电阻率及相位的最大误差也发生在山脊两侧，但山脊顶部无较明显的极值.

图 11

Fig. 11

图 11 考虑位移电流和不考虑位移电流的视电阻率及相位误差对比，频率为250 kHz Fig. 11 The absolute error of apparent resistivity and phase which calculated with and without displacement currents. The frequency is 250 kHz

实际模型来自安徽铜陵矿集区舒家店地区某剖面(图 12)，该剖面长3.8 km，具备以下两个特点：该区域存在高阻闪长岩，位移电流可能会造成显著影响；地形起伏明显，有两个山脊，地形对位移电流的影响不可忽略.

图 12

Fig. 12

图 12 舒家店地质综合剖面图(安徽省地质调查院提供) Fig. 12 Geological section map of the Shujiadian ore deposit (Provided by Geology Investigation Institute of Anhui Province)

基于野外露头及室内岩性测量，建立了电阻率模型(图 13).计算频率从10~250 kHz，共12个频点.网格剖分区域为x∈ -5000，5000，z∈[-5000，5000]， 经网格剖分后的总节点数为26055，单元数为51801，所有频点TE/TM计算总耗时为23.48 s.

图 13

Fig. 13

图 13 舒家店地球物理模型及网格剖分示意图红色代表空气，其中花岗闪长岩(粉红色)电阻率ρ1取1800 Ωm，砂岩、粉砂岩(蓝色)电阻率ρ2取420 Ωm，石英闪长岩(绿色)电阻率ρ3取12000 Ωm. Fig. 13 Illustration of geophysical model and mesh densities Red represent air，pink represent grano-diorite with resistivity of 1800 Ωm，blue represent siltstone with resistivity of 420 Ωm，green represent quartz diorite with resistivity of 12000 Ωm.

图 14为TE/TM模式下考虑位移电流和准静态条件下的视电阻率及相位误差.从图中可看出，误差整体随着频率的增加而增大；在x=-1000 m和x=500 m的高阻石英闪长岩体区域，考虑位移电流的结果和准静态结果存在明显的误差，而在x=-1000 m的山脊处这种误差更大.值得一提的是，相位受位移电流影响比视电阻率更为显著，即使x=1300 m处，粉砂岩电阻率为420 Ωm的情况下，在高频和山脊的共同作用下，相位也出现了明显的误差. ρ_TE的最大相对误差为14.5%，位于x=-950 m处，而在该处附近的误差也都在10%左右，此外，在x=-1600 m、x=500 m以及x=1500 m附近的误 差也都在5%以上；ρ_TM的最大相对误差高达78%，这是因为TM模式横向分辨率高，在电阻率分界面处二者计算结果有较大误差，ρ_TM也主要分布在两个山脊处和x=400 m的分界面处.

图 14

Fig. 14

图 14 TE/TM模式下考虑位移电流和准静态条件下的视电阻率及相位误差曲线横坐标x为测点，纵坐标y为频率的对数值，误差均为绝对误差. Fig. 14 The absolute error of apparent resistivity and phase which calculated in the quasi-static case and in the general case with displacement currents

图 15为考虑位移电流时计算的位移电流密度和传导电流密度，计算频率为250 kHz.图 15a，15b为电流密度的模.图 15c为位移电流密度占总电流密度的百分比.从图 15c中可看出，空气中位移电流密度占100%，也就是说空气中仅存在位移电流，无传导电流.在地下砂岩、粉砂岩对应的低阻区，位移电流密度在总电流密度中所占的比例为1%.花岗闪长岩区域的位移电流密度占3%左右.对于电阻率很高的石英闪长岩，位移电流密度所占比例达到10%，因此，在RMT数值模拟中不可忽略位移电流.

图 15

Fig. 15

图 15 考虑位移电流条件下计算的位移电流密度和传导电流密度，计算频率f=250 kHz(a)为传导电流密度的模；(b)为位移电流密度的模；(c)为位移电流密度占总电流密度的百分比. Fig. 15 The displacement currents density and conduction currents density which calculated with the consideration of displacement currents at f=250 kHz of the TE-mode (a)Module of displacement currents density;(b)Module of conduction currents density; (c)The percentage of displacement currents density.

大多数RMT的资料解释多是直接应用MT程序，忽略了位移电流的影响，这可能导致错误的解释结果；并且，在起伏地形条件下，位移电流的影响更为显著.针对这一问题，本文编写了起伏地形下的全电流2D RMT有限元程序.通过非结构的三角形单元实现了任意复杂模型的离散化，同时，程序采用局部加密策略获得了高精度的数值解.

由于位移电流的存在，空气层中磁场的偏导数不能近似为0，因而在对TM模式进行有限元离散时也必须考虑空气层.通过Dike模型讨论了空气层 厚度对数值解的影响，当空气层厚度大于1/4波长 即可保证TE/TM模式的有限元解的误差均在1%以内.

根据电流公式，位移电流随频率的升高而增大，传导电流随介质电阻率的增大而减小，因此在高频高阻情况下，位移电流在总电流中所占的比重就越大.文中结果表明，考虑位移电流情况下得到的视电阻率和相位曲线相比于准静态条件下的计算结果会出现下掉现象，且频率越高下掉越明显，这是因为位移电流项的引入增大了总电流密度，而总电流密度的增大会使得地下介质的电导率特性放大，从而导致视电阻率曲线下掉.数值算例表明，当背景电阻率为10000 Ωm时，频率大于100 kHz后，位移电流对视电阻率的影响不可忽略，而对相位曲线的影响从频率为20 kHz后就必须考虑.在背景电阻率与频率不变的情况下，地形起伏也会影响位移电流，地形高程越大，考虑位移电流计算的视电阻率与准静态条件下得到的视电阻率差异也越大，同时，数值结果表明差异最大的地方位于山脊的两侧.

通过铜陵舒家店矿床的实际地质模型，分别计算了传导电流和位移电流，结果表明在高阻花岗闪长岩地区位移电流不可忽略，会对模拟、观测结果造成显著的影响.如忽略位移电流，采用现有的MT程序进行RMT反演解释，势必带来结果的偏差，进 而影响资料解释的准确度.而采用本文提供的全电流起伏地形模拟程序，可以获得更为准确的模拟结果.

基于本文的研究成果，作者将进一步开发针对RMT的反演程序.

致谢 感谢Shewchuk提供的非结构三角形网格剖分开源代码Triangle，感谢Schenk提供的方程组求解器PARDISO，感谢网格剖分可视化软件PARAVIEW的作者，感谢安徽省地质调查院提供舒家店矿床地质模型.同时感谢舒家店野外物性测量人员王显莹、张弛等.