2019, Vol. 62
2. 中南大学有色金属成矿预测与地质环境监测教育部重点实验室, 长沙 410083;
3. 中南大学地球科学与信息物理学院, 长沙 410083;
4. 江苏省有色金属华东地质勘查局八一四队, 江苏镇江 212005
2. Key Laboratory of Metallogenic Prediction of Nonferrous Metals and Geological Environment Monitoring, Ministry of Education, Changsha 410083, China;
3. Institute of Applied Geophysics, School of Geosciences and Info-Physics, Changsha 410083, China;
4. Jiangsu Nonferrous Metal East China Geological Survey Bureau 814 Team, Zhenjiang Jiangsu 212005, China
射频大地电磁法(Radio-magnetotelluric, RMT)是通过无线电发射机发射10~300 kHz的场源信号，然后在远区采集电磁资料进而研究百米内地下构造的一种物探方法.近年来，射频大地电磁法被广泛地应用于工程勘探及环境监测(Pedersen et al., 2005, 2006; Tezkan et al., 2000, 2005; Tezkan and Saraev, 2008; Ismail et al., 2011; Bastani et al., 2013)等，已成为浅地表勘探的重要手段之一.
射频大地电磁法最早是由Müller及其课题组于1994年提出(Turberg et al., 1994)，当时，为了探测一地下孔隙层的含水情况，他们研制了探测频率为12~240 kHz的标量RMT-R仪器.随后，Tezkan等(1996, 2000)采用该仪器进行了废弃物的排查；Persson(2001)利用该仪器研究了某区域的断裂构造.Bastani(2001)研发了一个新的RMT探测系统(称为EnviroMT)，该系统实现了RMT法的张量测量，进一步推广了RMT法的实际应用.Bastani(2001)详细阐述了张量RMT法的数据处理方法，并采用该仪器探测了地下20 m处的砂岩中赋存的地下水；Bastani等(2009)将RMT法与可控源方法相结合探测某区域的热液型铜矿；Tezkan等(2005)采用RMT法进行了地表石油泄漏污染的排查，相比于传统的钻孔方法，既节省了勘探成本又提高了勘探效率；Bastani等(2013)采用RMT法判断了某地区石灰岩脉的走向，为进一步开采提供帮助.受仪器的限制，RMT法在国内的实际应用尚未开展.
由于RMT法是MT法的频带拓展，它通过无线电发射机发射10~300 kHz的场源信号，然后在远区探测电磁场实现百米内的电性勘探，因此RMT法的数据采集与处理手段与MT法类似(Bastani et al., 2013; Li et al., 2017, 2018, 2019).很多学者在进行RMT数据处理与解释时直接采用MT的数值模拟方法及反演程序进行(Tezkan et al., 2000;Newman et al., 2003; Linde and Pedersen, 2004; Candansayar and Tezkan, 2006, 2008; Ren et al., 2013a, 2013b, 2014a, 2014b).然而，直接套用MT反演软件进行RMT数据解释会导致反演得到的地下电阻率值出现不切实际的极值(Persson and Pedersen, 2002; Kalscheuer et al., 2008)，从而得不到合理的解释结果.Kalscheuer等(2008)基于有限差分正演算法，开发了考虑位移电流的2D RMT反演程序，并分析对比了准静态反演和全电流反演结果的差异，但是文中反演过程中并未考虑介电常数的变化.
在进行地球物理电磁数据反演时，通常都是根据观测资料(视电阻率、相位、倾子等)来反演得到地下介质的电阻率分布情况.这对于传统准静态假设下(中低频)的电磁反演是合理的，因为此时介电常数对电磁响应的影响微乎其微.对于射频大地电磁数据而言，是否考虑介电常数会使视电阻率产生高达20%以上的相对误差，而相位也会有十几度的误差(原源等，2015)，尤其是当地下岩体为结晶类高阻岩石(如石英岩等)(Chave and Jones, 2012).这时，不考虑地下介电常数分布，仅通过反演电阻率信息来拟合观测数据可能得到错误的资料解释结果.
部分学者研究了探地雷达(GPR)数据的多参数同步反演.Busch等(2012)研究了两层模型下GPR数据的全波形同步反演，对三种同步策略进行了对比；Lavoué等(2014)采用拟牛顿法实现了2D频率域GPR数据的同步反演；Da Gomes等(2014)采用两种反演算法(蚁群优化算法和拟牛顿法)，对比了二者的同步反演效果；Operto等(2013)研究了多参数全波形反演中合理的参数化选择策略、不同参数权重对反演的影响；Meles等(2010)基于时域有限差分正演实现了GPR数据的电导率-介电常数同步反演.然而，由于GPR数据受介电常数参数的影响更为明显，因此已发表的文献中部分研究结论并不适用于RMT数据的同步反演.
本文拟研究RMT数据同步反演中参数的变换、灵敏度的统一等问题，实现RMT数据的电导率-介电常数同步反演.
1 考虑位移电流的2D RMT正演问题反演迭代中的正演算子F(m)采用基于非结构网格的有限单元法进行求解(原源等，2015).假设x方向为构造走向方向，z方向垂直向下，根据Maxwell方程组推导出全电流条件下的2D电场磁场满足的Helmholtz方程(取时间因子eiwt)：
|
(1a) |
|
(1b) |
其中，ω为角频率，σ为电导率，ε为介电常数，μ0为磁导率.式(1a)为TE模式，式(1b)为TM模式.文中采用层状解析解作为边界条件，经过有限元离散后可得到如下线性方程组：
|
(2) |
其中x为TE模式的电场Ex或TM模式的磁场Hx，A为刚度矩阵，B为与边界条件有关的源项.
文中TE模式的磁场和TM模式的电场分量通过四点数值微分得到：
|
(3a) |
|
(3b) |
视电阻率和相位计算公式如下：
|
(4) |
|
(5) |
|
(6) |
在反演算法中，存在两种网格：正演网格和反演网格(Günther et al., 2006;韩琦等, 2015;吴小平等, 2015).正演网格用于计算反演迭代后的地电模型的电磁响应或阻抗、视电阻率等参数.通常，为了得到地表高精度的数值解，正演网格在地表附近剖分较密；对于反演网格，由于我们关心的是地下电阻率的分布规律，因而只需较疏的网格，既降低了反演自由度同时减小了反演求解过程中的不适定性.所以，在合理的反演策略中通常采用相互独立的正反演网格，这样既保证了正演的计算精度，同时避免了由网格引起的不必要的反演计算量.
2.1 双网格生成策略本文首先采用开源代码Triangle(Shewchuk, 1996)对带地形模型进行网格离散，即生成反演粗网格.反演网格中，通过设定三角形最下角不小于30°来保证每个单元网格质量，同时在反演目标区域的网格相对较密，越靠近边界条件网格越疏.由于在反演过程中涉及到正演计算和灵敏度的求取，每次反演迭代后得到的修正模型参数需要传递给正演单元.因此在生成正演网格时我们希望反演网格是正演网格的子集，即不破坏原来的反演网格，从而实现参数的无偏传递.生成正演网格的具体方法如下：首先根据反演网格的边信息构建模型输入文件；然后通过限定每个反演单元的最大面积建立单元面积局部加密文件；最后判断每个反演单元面积，若超过设定的最大面积则进行单元剖分，最终保证每个反演单元中所包含的子单元面积均小于设定的最大面积，从而生成局部加密的正演网格，如图 1所示.需要特别指出的是，由于RMT数值计算中考虑了位移电流的影响，因而其网格剖分与准静态近似条件下的MT模拟略有不同.
|
图 1 非结构双网格示意图 (a)起伏地形模型; (b)反演粗网格; (c)正演密网格. Fig. 1 Illustration of the unstructured duple mesh (a) A 2-D terrain model; (b) Coarse mesh for inversion; (c) Refined mesh for forward modeling. |
(1) 通常，电磁法数值计算时网格剖分区域的大小一般为5~10倍的趋肤深度，对于考虑了位移电流的RMT计算，其趋肤深度δ′=1/β′= 
(2) 对于网格单元节点间距，由于考虑位移电流情况下1/α′ < 1/β′，因此应保证最小网格单元节点间距小于最大观测频率下所计算的1/α′，其中
(3) 与MT法不同，RMT法计算TM模式的磁场时，空气中的磁场不能看做常数，因此TM模式也需要考虑空气层.在高频情况下，电磁场在空气中以电磁波的形式向外传播，空气中电磁波的传播速度为光速c=3×108 m·s-1，根据c=λ·f，其中λ为波长，f为频率，据此可计算出不同频率下的波长，本文中网格剖分时最短波长内保证有20~30个节点.此外，在网格剖分时应注意缩小测点及电阻率突变区域附近的网格单元大小，以保证该区域的计算精度.
2.2 正反演网格映射在得到修正后的反演网格上的电性参数后，需要将反演网格上的参数映射到正演网格中以进行下一次反演迭代.对于正反演网格参数映射，最简单直接的办法就是循环搜索，假设正演网格数为MFWD，反演网格数为MINV，将每一个正演网格对所有反演网格进行搜索，当第i个正演网格的三角形中心点与某一个反演网格的三个节点相连得到的三个角度之和为360°，那么这个正演网格就位于该反演网格内，反之则位于该反演网格之外.这种循环搜索的时间复杂度为O(MFWDMINV)，当模型规模较大时，这样遍历搜索就会比较耗时(例如3D情况).对于本文而言，在作者2.3 GHz CPU，2GB RAM的笔记本上，假设反演迭代10次，网格映射总耗时仅占反演总耗时的2%.因此作者采用遍历法实现网格参数映射.
3 双参数同步反演算法 3.1 目标函数本文双参数反演迭代采用光滑约束的Gauss-Newton法，构建目标函数时同时考虑地下电阻率和介电常数两个参数，双参数目标函数如下所示:
|
(7) |
其中目标函数Φ(m1, m2)是与地下介电常数m1和电阻率m2均相关的函数；d为N维观测资料向量，d=(d1, d2, …, dN)T；F(m1, m2)是由m1和m2共同作用得到的正演响应.Cd－1为观测资料d和模型正演响应F(m1, m2)的数据拟合差权重，Cd－1=diag(1/εi)；Cm1-1和Cm2-1分别为与电阻率和介电常数相关的光滑度矩阵，其取值为相邻单元i，j的模型标准差(Günther et al., 2006).m10和m20分别为地下电阻率和介电常数的先验信息.λ和γ为调节两个模型拟合差与数据拟合差之间的权重因子.
对正演算子F(m1, m2)进行二元Taylor展开：
|
(8) |
将二元目标函数Φ(m1, m2)分别对m1和m2求导数：
|
(9a) |
|
(9b) |
令式(9a)和(9b)分别为零，得到如下等式：
|
(10a) |
|
(10b) |
联立(10a)和(10b)后可得如下反演迭代方程组：
|
(11) |
令(11)式中，
|
(12) |
其中，A、B、C、D均为M×M阶矩阵，x1、x2、y1、y2均为M×1阶向量.通过求解方程组(12)即可得到第k次迭代的模型修正量Δm1k和Δm2k，然后根据(13)式即可得到新的模型参数：
|
(13a) |
|
(13b) |
在单参数反演中，为保证反演过程的稳定性，反演模型参数通常取电阻率对数值.而本文双参数反演中，如何通过对反演参数的合理转换来保证反演过程的稳定是一大难点.
3.2 反演参数尺度变换在进行双参数同步反演前，首先要研究不同地电参数(σ, ε)对正演响应的灵敏度，然后根据参数灵敏度的不同对σ和ε进行合理的尺度变换，以保证反演的稳定性.我们首先以均匀半空间模型为例来讨论电导率和介电常数的模型改变对电场场值的影响.假设均匀半空间电导率为σb，介电常数为εb，那么根据地表电场表达式Ex=E0xeiωte－kz，将Ex分别对σ和ε求导可得
|
(14) |
其中，
|
(15) |
根据差分公式，电导率改变δσb和介电常数改变量δεb引起的电场场值改变量δExσb和δExεb分别为
|
(16a) |
|
(16b) |
将(14)、(15)式代入(16)式，并将(16a)除以(16b)得
|
(17) |
通常，令模型改变量为原来的Δ倍，即
|
那么(17)式可转换为
|
(18) |
从(18)式可看出，在MT的勘探频率下σb≫ωεb，即|δExσb|≫|δExεb|，也就是说此时观测数据对电导率更为灵敏；类似地，在GPR的勘探频率下σb≪ωεb，即|δExσb|≪|δExεb|，这就意味着GPR观测数据对介电常数更为灵敏(尤其是在高阻地区)；而在RMT勘探频段，σb≈ωεb，此时，观测数据对二者灵敏度相当.然而，考虑到实际中电导率的变化范围更广(10－4~0.1 S·m-1)，而介电常数的变化范围为1~81，这使得实际反演中，电导率更容易引起观测数据的改变，而削弱了介电常数的影响.此外，从(17)式中我们可看出，电导率的灵敏度与频率的一次方成正比，而介电常数的灵敏度与频率的二次方成正比，因此，随着频率的升高，介电常数对观测数据的影响更大.为此，我们期望寻求一种电导率变换式，使得其数据尺度缩小，同时保证频率对电导率灵敏度的影响与介电常数相当.
我们定义复介电常数ε′为
|
(19) |
(14) 式的实部代表位移电流的贡献，它在电磁波的传播过程中不会产生能量损耗；而虚部代表传导电流的贡献，由于传导电流的传播过程中会产生焦耳热，从而使能量有所损耗.根据物理学中的定义，介质损耗角正切tanδ0是表征电磁波在介质中传播时将电磁能量转化为热能所消耗的能量，其数学表达式为
|
(20) |
令(20)式中的tanδ0=1，此时，位移电流与传导电流对场的贡献相等，我们定义tanδ0=1时的地电参数为参考场.σ0为参考场的电导率，ω0为参考频率，ε0为真空中的介电常数.由此可得参考场的电导率σ0=ω0ε0.与相对介电常数εr=ε/ε0的定义类似，我们可定义一个相对电导率参数σr=σ/σ0，然后在双参数反演中反演相对电导率和相对介电常数.经过这一变换后，观测数据对σr和εr的灵敏度的频率依赖性相当，并且缩小了电导率参数的变化范围.
为测试观测数据对经过尺度变换后的反演参数的灵敏度，我们以图 2中的均匀半空间中赋存一矩形异常体为例，研究视电阻率和相位对电导率参数变换前后的灵敏度.矩形异常体埋深为20 m，长宽分别为100 m和40 m.图 2中黑线为计算灵敏度时所用的正演网格，为保证计算的准确性这里采用较密的网格剖分，网格剖分范围为x∈[－5 km, 5 km]，z∈[－5 km, 5 km]，总节点数为19143，单元数为38032.
|
图 2 理论模型及正演网格剖分示意图 半空间电阻率为10000 Ωm，相对介电常数为5，矩形电阻率为1000 Ωm，相对介电常数为1. Fig. 2 A synthetic model and the unstructured mesh for forward calculation The resistivity and relative permittivity of half-space medium are 10000 Ωm and 5. The resistivity and relative permittivity of anomalous block are 1000 Ωm and 1. |
图 3为TE模式下计算得到的视电阻率及相位对不同反演参数的灵敏度.观测数据点位于(0, 0)处，计算频率为250kHz，计算耗时12.76s.图 3(a—c)为视电阻率对以10为底的对数电阻率log10(ρ)、相对电导率σr及相对介电常数εr的灵敏度；图 3(d—f)为相位的灵敏度.从3a和3d两图中可看出，观测数据对log10(ρ)的灵敏度在10-2到10-3数量级上，而3c和3f图中，观测数据对εr的灵敏度较小，为10-4到10-5数量级上，在这种情况下，由于观测数据对log10(ρ)的灵敏度远大于εr，也就是说εr对观测数据的影响微乎其微，导致反演迭代中对相对介电常数的修正几乎为零(反演结果参看图 8).图 3b和3d为通过本节中所定义的相对电导率对电阻率参数进行变换后的灵敏度，参数变换中参考频率ω0取100 kHz.从图中可看出，经过尺度变换后，视电阻率和相位对σr的灵敏度也处于10-4到10-5数量级上，从而保证σr和εr对观测数据有相当的贡献，为双参数的反演奠定基础.
|
图 3 图 2所示模型的TE模式下视电阻率及相位对对数电阻率、相对电导率、相对介电常数的灵敏度.其中观测数据点位于(0, 0)处，计算频率为250 kHz (a)—(c)分别为视电阻率对对数电阻率、相对电导率及相对介电常数的灵敏度，(d)—(f)为相位对三个参数的灵敏度. Fig. 3 The sensitivity of the TE-mode apparent resistivity and phase to logarithmic resistivity, relative conductivity and relative permittivity for the model shown in Fig. 2 The measurement point is located at (0, 0) and the frequency is 250 kHz. Subplots (a—c) shows sensitivity of the apparent resistivity and subplots (d—f) shows sensitivity of the phase. |
从3.2节中的相对电导率的定义可看出，相对电导率的大小随着参考频率的变化而改变，进一步而言，参考频率直接影响观测数据对相对电导率的灵敏度的大小(如附录B所示).为此，我们有必要研究参考频率的选取对反演过程的影响.为简单起见，在研究参考频率时我们令λ和γ均为零，即忽略模型拟合项，同时令Cd－1等于单位矩阵I.经过上述简化后，(12)式左端的系数矩阵可表示为
|
(21) |
由于参考频率的大小只影响J2k，因此(21)式中J1kTJ1k与参考频率无关.
图 4为三个不同参考频率下计算得到的系数矩阵.系数矩阵的维数为2MINV×2MINV，各个子图的横纵坐标取系数矩阵的角标，图中符号S1、S2、S3、S4分别表示系数矩阵的左上部分J1kTJ1k、右上部分J1kTJ2k、左下部分J2kTJ1k和右下部分J2kTJ2k.图 4a中参考频率ω0=10 kHz，由于参考频率较小，使得系数矩阵中S4的数值小于S1，也就是说εr对观测数据的贡献更大；图 4b中参考频率ω0=50.1 kHz，在这个参考频率下，系数矩阵中S4的数值和S1相当，也就是说此时εr和σr对观测数据的贡献差不多；而图 4c中，当参考频率ω0=100 kHz时，系数矩阵中S4的数值要大于S1，这意味着此时σr对观测数据的贡献更大.
|
图 4 不同参考频率的选取对反演系数矩阵的影响 Fig. 4 The effect of reference frequency on the inverse coefficient matrix |
从反演的角度来看，我们通过引入相对电导率的概念就是希望使双参数εr和σr对观测数据的影响相当，避免某一参数的影响被掩盖，从而保证双参数反演迭代能够进行.根据图 4的结果，我们进一步从更直观的角度研究不同参考频率ω0对反演结果的影响.为了避免其他反演参数对结果的影响，在反演中我们同样令λ和γ均为零，即忽略模型拟合项，同时令Cd－1等于单位矩阵I，经过这一假设后，(12)式的反演方程组可简化为
|
(22) |
图 6为在上述简化后的情况下不同参考频率的同步反演结果.三个参考频率分别取10、50.1和100 kHz，反演迭代次数均为15次，最终三个参考频率下的数据拟合差分别为0.113、0.10和0.097，反演迭代总耗时分别为6266、6322、6223 s.反演网格仍采用非结构的双网格，经网格剖分后的反演网格数为3605，正演网格数为11905.网格剖分如图 5所示，其中黑线为反演网格，红线为正演网格，网格映射耗时11.45 s.反演的初始模型取背景电阻率10000 Ωm和相对介电常数5.图 6(a—c)是参考频率分别取10、50.1和100 kHz时反演得到的电导率分布，从图中可看出，当参考频率ω0为10 kHz时，反演结果无法得到异常电导率分布信息(图 6a)，这是因为参考频率较小时会降低相对电导率的灵敏度，使得每次反演迭代中电导率的改变非常微弱，进而无法得到异常信息；当ω0增加到50.1 kHz时，电导率异常有所显示，但是异常仍不太明显(图 6b)；而当参考频率增至100 kHz时，电导率异常分布位置和大小均能准确地反演出来(图 6c)，这是因为较大的参考频率增加了相对电导率的灵敏度.图 6d为参考频率为100 kHz时反演得到的相对介电常数分布，这里仅给出100 kHz下的相对介电常数反演结果是因为介电常数灵敏度大小不受参考频率的影响，不同参考频率反演得到的介电常数分布类似.从反演结果来看，相对介电常数的反演能够准确地区分出异常体的上下边界，但是对左右边界区分度不够.
|
图 5 反演网格剖分示意图，其中黑线为反演网格，灰线为正演网格 Fig. 5 An illustration of mesh for inversion calculation. The inversion mesh is shown in black line and the forward mesh is shown in gray line |
|
图 6 不同参考频率下σr和εr同步反演结果 三个参考频率分别为10、50.1和100 kHz；(a—c)为三个参考频率下反演得到的电导率图像，(d)为100 kHz下同步反演得到的相对介电常数分布图.图中下三角为测点位置，虚线框为异常体所在位置. Fig. 6 The σr and εr simultaneous inversion with different reference frequencies Three reference frequencies are 10, 50.1 and 100 kHz. Subplots (a—c) shows the inversed conductivity imaging and subplot (d) shows the relative permittivity imaging which calculated with reference frequency of 100 kHz. The lower triangle shows the location of measurement points and the dotted line indicated the position of anomalous body. |
综上分析，我们可得到以下结论：在双参数反演中，平衡两个参数对观测数据的贡献对反演结果至关重要.若反演参数选取为εr和σ，由于σ的灵敏度远远大于εr(相差几个数量级)，在反演过程中会更多地进行σ的修正，而εr几乎不变，因而无法得到正确的反演结果；若反演参数选取为εr和σr，由于σr的灵敏度与ω0成正比，因此ω0的选取对反演结果影响很大，从图 6的算例中可看出，ω0的选择至少要保证σr对观测数据的贡献大于或等于εr才能得到正确的反演结果.同时，考虑到RMT数据对电导率参数的依赖程度更高，因此在实际中取略大的ω0，使σr的贡献略大于εr能够得到更为合理的反演结果.
3.4 双参数反演流程图 7为双参数同步反演流程图.
|
图 7 基于GN算法的双参数同步反演流程图 Fig. 7 The flow chart of dual-parameters simultaneous inversion based on GN algorithm |
本节对第3节所述的双参数同步反演算法进行理论模型测试，分析反演算法的正确性，同时探讨不同参考频率ω0及权重因子λ、γ的选择对反演结果的影响.
4.1 算例一 4.1.1 理论模型及反演数据的生成理论模型如图 2所示，均匀半空间中赋存一矩形异常体，背景电阻率为10000 Ωm，相对介电常数εr=5；矩形异常体电阻率为1000 Ωm，相对介电常数εr=1.计算时，空气层的电阻率为1.0×1016Ωm.图 2中黑线为非结构正演网格剖分示意图，为保证计算精度，对正演响应的计算采用较密的网格剖分.网格剖分范围为x∈[－5 km, 5 km]，z∈[－5 km, 5 km]，总节点数为17945，单元数为35626.网格剖分在测点附近进行了加密，以保证计算精度.正演计算频点数为25，从1 kHz到251 kHz，每个数量级上对数等间距分布25个频点.测点从-500 m到500 m，间距为50 m，总测点数为21个.对所有频点均进行TE和TM模式计算，总耗时为28.43 s.得到所有视电阻率和相位数据后，我们对理论响应数据添加2%的高斯噪声作为反演数据，反演数据量N=25×21×4.
4.1.2 反演结果分析本节讨论正则化因子λ和γ的选取对反演结果的影响，并对双参数反演效果进行探讨.
我们首先测试不进行参数变换的同步反演效果.反演数据为4.1.1节中所给数据，反演网格如图 5所示，反演初始模型为电阻率10000 Ωm和相对介电常数5的均匀半空间.图 8为不进行参数变换情况下得到的反演结果，即反演双参数直接取为log10ρ和εr，反演迭代15次总耗时6530 s，最终数据拟合差RMS为4.2%，反演迭代时不考虑模型约束项，同时令数据权重Cd－1=I.图 5a为电阻率反演结果，从图中可看出同步反演能够准确地反演出电阻率的异常位置和大小；图 5b为同步反演得到的相对介电常数模型，反演结果无法反映出介电常数异常，这是因为观测数据对电阻率的灵敏度远远大于介电常数，从而在反演迭代步中侧重于电阻率参数的修正，而相对介电常数几乎不变，因而最终相对介电常数分布仍与初始模型相当.图 9为反演方程组的系数矩阵，图中显示S1区域的数值远大于S4区域，也就是说相对介电常数对观测数据的贡献远不及电阻率，因此介电常数在反演迭代中几乎无变化.
|
图 8 对数电阻率和相对介电常数同步反演结果 Fig. 8 The result of simultaneous inversion for logarithmic resistivity and relative permittivity parameters |
|
图 9 log10ρ和εr同步反演时的系数矩阵 其中S1区域为J1TJ1，J1是正演响应对log10ρ的灵敏度; S4区域为J2TJ2，J2是正演响应对εr的灵敏度; S2和S3区域分别为J1TJ2和J2TJ1. Fig. 9 The coefficient matrix of simultaneous inversion for log10ρ and εr parameters The region S1, S2, S3 and S4 shows the element values of J1TJ1, J1TJ2, J2TJ1 and J2TJ2. Where J1 and J2 are sensitivities of forward responses to log10ρ and εr. |
根据上述分析，在进行RMT双参数反演时，首先要对待反演模型参数(电导率和介电常数)进行尺度变换，如本文所采用的相对电导率σr和相对介电常数参数εr，使得两个反演参数对观测数据的贡献相当，从而保证双参数反演结果是有效的(如图 6c、6d所示).
从图 6中可看出，当令λ和γ均为零时，得到的相对介电常数模型左右边界分辨率较差，为此，我们对目标函数增加介电常数约束.图 10为考虑介电常数模型约束的反演结果，反演迭代次数为15次.图 10(a—d)分别为λ取5和0.5时得到的电导率和相对介电常数模型.从图中可看出，当λ取5时，能够反演出电导率和相对介电常数异常，但是异常数值非常小，电导率范围为0.95~1.05×10－4，相对介电常数范围为4~6，这是因为当λ取值较大时，反演迭代步长很小，收敛速度很慢，因而在15次迭代后无法显著地突出异常；当λ取0.5时，反演结果能够准确地反映出电导率和相对介电常数异常，从图 10c中可看出反演得到的电导率模型与图 6c类似，而图 10d中经过对相对介电常数进行模型约束后，发现反演得到的相对介电常数模型的边界分辨率较图 6d有所提高；此外，作者计算了当λ取0.05和0.005时的反演(限于篇幅，未给出反演结果图)，结果显示虽然电导率异常有所反映，但是相对介电常数的反演结果完全是失败的，尽管在这两种情况下模型正演结果与观测数据具有更小的拟合差(如图 11所示).
|
图 10 只考虑介电常数模型约束的反演结果 (a—b)是λ=5时反演得到的σ和εr模型; (c—d)是λ=0.5时反演得到的σ和εr模型. Fig. 10 The simultaneous inversion result which calculated without consideration of permittivity model constraint The subplots (a—b) show the inversed conductivity σ and relative permitivity εr imaging which computed with λ=5 and subplots (c—d) show the inversion results computed with λ=0.5. |
|
图 11 λ取不同值时得到的反演迭代收敛曲线 Fig. 11 The inversion convergence curves for different λ |
综上所述，相对介电常数模型约束的权重因子λ的取值既不能过大，使得收敛速度非常缓慢(图 11中λ=5的曲线)，在合理的迭代步中无法突显出异常体，同时，λ的取值也不能太小，虽然当λ取值很小时会得到更小的数据拟合差，但是，过小的λ使得反演结果不可靠，甚至是错误的.
图 12和图 13分别为λ取0.5，γ分别取0.1和0.01时εr和σr同时反演的结果.反演迭代次数均为15次.从图中可看出，不同γ的取值影响电导率反演模型分布，而相对介电常数模型无明显差异，这是因为介电常数的模型目标函数权重λ相同；当γ取值较大时，电导率模型在异常体上方出现了高阻异常，且数据拟合差无法稳定收敛；而当γ相对较小时，电导率模型能够准确反演出异常体，同时数据拟合差也逐步收敛.此外，我们也测试了γ取0.001和0.0001的反演结果，得到的模型分布与图 13类似.
|
图 12 λ=0.5, γ=0.1时εr和σr同时反演结果 Fig. 12 The simultaneous inversion for εr and σr parameters with regularization factors of λ=0.5, γ=0.1 |
|
图 13 λ=0.5, γ=0.01时εr和σr同时反演结果 Fig. 13 The simultaneous inversion for εr and σr parameters with regularization factors of λ=0.5, γ=0.01 |
根据前面的对比分析及诸多测试，在进行εr和σr同时反演时，一般情况下选取λ为JT2J2的最大特征值，γ为J1TJ1的最大特征值能够得到较为合理的结果.
4.2 算例二：低阻高极化模型本节通过对更多的理论模型进行反演计算来验证所开发的程序的正确性.所有理论模型反演计算中，均对理论模型的正演响应添加2%的高斯噪声生成反演数据.反演初始模型均为均匀半空间，测点在图中采用下三角表示，所有算例的计算频点均为25个，从1 kHz到251 kHz对数等间距分布.
本算例为理论模型反演算例.理论模型为山谷下方存在一低阻高极化异常体，异常体位于x∈[－100, 100], z∈[－130, －180]处，山脊背景电阻率为10000 Ωm，相对介电常数为5，异常体电阻率为5000 Ωm，相对介电常数为30.对该模型进行TE/TM模式的25个频点进行正演计算，计算得到的视电阻率和相对添加2%的高斯噪声作为反演数据，总反演数据量为2500.
图 14为相对电导率-相对介电常数同步反演结果.初始模型取背景电阻率和介电常数值.图 14a为反演中双网格剖分示意图，反演网格单元数为5322，正演网格单元数为20126.图 14b为迭代12次得到的电阻率单参数反演结果，反演中相对介电常数固定为1，初始正则化因子λ为50，然后在每个迭代步中以0.7的倍数衰减.图 14c和14d为迭代15次后的双参数反演得到的电导率分布和相对介电常数分布，参考频率f0=100 kHz，λ和γ分别取0.5和0.1，总计算耗时16788 s.图 14b中，由于理论模型的电阻率差异不大，单参数的电阻率反演在地表浅部出现了小范围的低阻异常，低阻体下方出现了高阻异常，这与理论模型出现了较大的出入；从图 14c中可看出，双参数反演得到的电导率分布与理论模型吻合的很好，准确地反演出低阻体的位置；图 14d中双参数反演得到的相对介电常数分布能够反映出地下的高介电异常体，但是对异常体的边界区分不够明显.该算例表明在地下电阻率差别不太明显的区域，双参数的反演能够在一定程度上提高反演准确性.
|
图 14 低阻高极化模型的双参数同步反演结果与单参数反演结果对比 Fig. 14 A comparison of dual-parameters simultaneous inversion result and single-parameter inversion result for model with anomalous body of low resistivity and high polarizability |
本文实现了基于非结构双网格的双参数同步反演程序.由于位移电流在高阻背景下或起伏明显的地形中会对RMT正演响应带来不可忽视的影响，也就是说，实际中观测数据是由地下电阻率和相对介电常数共同决定的，因此采用传统的单参数电阻率反演来拟合观测数据存在一定的偏差.为此，本文研究了双参数同步反演算法.
首先，为了能够处理复杂地形，本文中网格采用非结构的双网格.然后，通过定义参考频率，提出了相对电导率的概念，确保双参数对正演响应的灵敏度相当，进而保证反演迭代能够稳定收敛；构建了双参数目标函数，包括数据拟合项、相对电导率模型拟合项、相对介电常数模型拟合项，讨论了模型正则化因子的选取.最后通过一理论模型算例研究了各种参数对反演结果的影响.得出如下结论：(1)参考频率ω0的选择应以统一相对电导率和相对介电常数的灵敏度为目的，考虑到RMT数据对电导率参数的依赖程度更高，因此在实际中取略大的ω0，使σr的贡献略大于εr能够得到更为合理的反演结果.(2)相对介电常数模型约束的权重因子λ的取值既不能过大，使得收敛速度非常缓慢，在合理的迭代步中无法突显出异常体，同时，λ的取值也不能太小，因为过小的λ使得反演结果不可靠，甚至是错误的.(3)相对电导率模型约束的权重因子γ取为0或者较小的值都可得到合理的反演结果，但是γ不能取值过大，这会导致反演无法稳定收敛.
由于实际观测数据是由地下介质的多种电性参数共同作用产生，因此从理论上来说，双参数反演相比单参数反演更为合理.然而，在相同的观测数据下，双参数反演待求解的未知模型参数个数是单参数的2倍，这就增加了反演的非唯一性.因此，在实际应用中，建议先进行单参数的电阻率反演，然后再进行双参数反演，综合对比分析来降低反演的多解性.
附录 灵敏度计算A相对介电常数灵敏度
计算观测数据对电阻率的灵敏度时，视电阻率取以10为底的对数值.对TE模式而言，根据正演方程组Ax=B求得所有节点上的电场Ex后，地表测点i处的电场和磁场可由下式求取：
|
(A1) |
其中ai是在测点i处元素为1、其余元素为0的向量；bi是在四点差分节点上有值、其余位置为0的向量，
|
(A2a) |
|
(A2b) |
双参数反演时，模型的介电常数参数取相对值εr，εr=ε/ε0, 其中ε为实际介电常数值，ε0为真空中的介电常数.
|
(A3) |
|
(A4) |
对TE模式而言，根据(A1)和(A2)式的定义，有
|
(A5) |
其中，
|
(A6) |
|
(A7) |
对TM模式而言，根据正演方程组Ax=B求得所有节点上的电场Hx后，地表测点i处的电场和磁场可由下式求取：
|
(A8) |
其中ai是在测点i处元素为1、其余元素为0的向量；bi是在四点差分节点上有值、其余位置为0的向量，
|
(A9a) |
|
(A9b) |
根据(A8)和(A9)式的定义，
|
(A10) |
其中，
|
(A11) |
|
(A12) |
B相对电导率灵敏度
双参数反演时，模型的电导率参数取相对电导率σr，σr=σ/σ0，其中σ为地下实际电导率，σ0=ε0ω0为等效电导率.
|
(B1) |
|
(B2) |
对TE模式而言，根据(A1)和(A2)式的定义，
|
(B3) |
其中，
|
(B4) |
|
(B5) |
对TM模式而言，根据(A8)和(A9)式的定义，
|
(B6) |
其中，
|
(B7) |
|
(B8) |
|
(B9) |
Bastani M. 2001. EnviroMT-a new controlled source/radio magnetotelluric system[Ph.D. thesis]. Uppsala: Acta Universitatis Upsaliensis.
|
Bastani M, Malehmir A, Ismail N, et al. 2009. Delineating hydrothermal stockwork copper deposits using controlled-source and radio-magnetotelluric methods:A case study from northeast Iran. Geophysics, 74(5): B167-B181. DOI:10.1190/1.3174394 |
Bastani M, Persson L, Beiki M, et al. 2013. A radio magnetotelluric study to evaluate the extents of a limestone quarry in Estonia. Geophysical Prospecting, 61(3): 678-687. DOI:10.1111/gpr.2013.61.issue-3 |
Busch S, Van Der Kruk J, Bikowski J, et al. 2012. Quantitative permittivity and conductivity estimation using full-waveform inversion of on-ground GPR data. Geophysics, 77(6): H79-H91. DOI:10.1190/geo2012-0045.1 |
Candansayar M E, Tezkan B A. 2006. A comparison of different radiomagnetotelluric data inversion methods for buried waste sites. Journal of Applied Geophysics, 58(3): 218-231. DOI:10.1016/j.jappgeo.2005.07.001 |
Candansayar M E, Tezkan B. 2008. Two-dimensional joint inversion of radiomagnetotelluric and direct current resistivity data. Geophysical Prospecting, 56(5): 737-749. DOI:10.1111/gpr.2008.56.issue-5 |
Chave A D, Jones A G. 2012. The Magnetotelluric Method:Theory and Practice. Cambridge: Cambridge University Press.
|
Da Gomes M G, Souto R P, De Athayde A S, et al. 2014. Estimating dielectric permittivity and electric conductivity from simulated multichannel GPR pulses using aco and quasi-newton inversion techniques. Revista Brasileira de Geofísica, 32(4): 595-614. DOI:10.22564/rbgf.v32i4.535 |
Günther T, Rucker C, Spitzer K. 2006. Three-dimensional modeling and inversion of dc resistivity data incorporating topography-Ⅱ. Inversion. Geophysical Journal International, 166(2): 506-517. DOI:10.1111/gji.2006.166.issue-2 |
Han Q, Hu X Y, Cheng Z P, et al. 2015. A study of two-dimensional MT inversion with steep topography using the adaptive unstructured finite element method. Chinese Journal of Geophysics (in Chinese), 58(12): 4675-4684. DOI:10.6038/cjg20151228 |
Ismail N, Schwarz G, Pedersen L B. 2011. Investigation of groundwater resources using controlled-source radio magnetotellurics (CSRMT) in glacial deposits in Heby, Sweden. Journal of Applied Geophysics, 73(1): 74-83. |
Kalscheuer T, Pedersen L B, Siripunvaraporn W. 2008. Radiomagnetotelluric two-dimensional forward and inverse modelling accounting for displacement currents. Geophysical Journal International, 175(2): 486-514. DOI:10.1111/gji.2008.175.issue-2 |
Lavoué F, Brossier R, Métivier L, et al. 2014. Two-dimensional permittivity and conductivity imaging by full waveform inversion of multioffset GPR data:A frequency-domain quasi-Newton approach. Geophysical Journal International, 197(1): 248-268. DOI:10.1093/gji/ggt528 |
Li G, Xiao X, Tang J T, et al. 2017. Near-source noise suppression of AMT by compressive sensing and mathematical morphology filtering. Applied Geophysics, 14(4): 581-589. DOI:10.1007/s11770-017-0645-6 |
Li J, Zhang X, Gong J Z, et al. 2018. Signal-noise identification of magnetotelluric signals using fractal-entropy and clustering algorithm for targeted de-noising. Fractals, 26(2): 1840011. DOI:10.1142/S0218348X1840011X |
Li J, Zhang X, Tang J T, et al. 2019. Audio magnetotelluric based on multifractal spectrum and matching pursuit. Fractals, 27(1): 1940007. DOI:10.1142/S0218348X19400073 |
Linde N, Pedersen L B. 2004. Characterization of a fractured granite using radio magnetotelluric (RMT) data. Geophysics, 69(5): 1155-1165. DOI:10.1190/1.1801933 |
Meles G A, Van Der Kruk J, Greenhalgh S A, et al. 2010. A new vector waveform inversion algorithm for simultaneous updating of conductivity and permittivity parameters from combination crosshole/borehole-to-surface GPR data. IEEE Transactions on Geoscience and Remote Sensing, 48(9): 3391-3407. DOI:10.1109/TGRS.2010.2046670 |
Newman G A, Recher S, Tezkan B, et al. 2003. 3D inversion of a scalar radio magnetotelluric field data set. Geophysics, 68(3): 791-802. DOI:10.1190/1.1581032 |
Operto S, Gholami Y, Prieux V, et al. 2013. A guided tour of multiparameter full-waveform inversion for multicomponent data:From theory to practice. The Leading Edge, 32(9): 1040-1054. DOI:10.1190/tle32091040.1 |
Pedersen L, Bastani M, Dynesius L. 2005. Groundwater exploration using combined controlled-source and radiomagnetotelluric techniques. Geophysics, 70(10): G8-G15. |
Pedersen L, Bastani M, Dynesius L. 2006. Some characteristics of the electromagnetic field from radio transmitters in Europe. Geophysics, 71(6): G279-G284. DOI:10.1190/1.2349222 |
Persson L. 2001. Plane wave electromagnetic measurements for imaging fracture zones[Ph.D. thesis]. Uppsala: Acta Universitatis Upsaliensis.
|
Persson L, Pedersen B. 2002. The importance of displacement currents in RMT measurements in high resistivity environments. Journal of Applied Geophysics, 51(1): 11-20. |
Ren Z Y, Kalscheuer T, Greenhalgh S, et al. 2013a. A goal-oriented adaptive finite-element approach for plane wave 3-D electromagnetic modeling. Geophysical Journal International, 194(2): 700-718. DOI:10.1093/gji/ggt154 |
Ren Z Y, Kalscheuer T, Greenhalgh S, et al. 2013b. Boundary element solutions for broad-band 3-D geo-electromagnetic problems accelerated by an adaptive multilevel fast multipole method. Geophysical Journal International, 192(2): 473-499. DOI:10.1093/gji/ggs043 |
Ren Z Y, Kalscheuer T, Greenhalgh S, et al. 2014a. A finite-element-based domain-decomposition approach for plane wave 3D electromagnetic modeling. Geophysics, 79(6): E255-E268. DOI:10.1190/geo2013-0376.1 |
Ren Z Y, Kalscheuer T, Greenhalgh S, et al. 2014b. A hybrid boundary element-finite element approach to modeling plane wave 3D electromagnetic induction responses in the Earth. Journal of Computational Physics, 258: 705-717. DOI:10.1016/j.jcp.2013.11.004 |
Shewchuk J. 1996. Triangle:Engineering a 2D quality mesh generator and Delaunay triangulator. Lecture Notes in Computer Science, 1148: 203-222. DOI:10.1007/BFb0014474 |
Tezkan B, Goldman M, Greinwald S, et al. 1996. Joint application of radiomagnetotellurics and transient electromagnetics to the investigation of a waste deposit in Cologne (Germany). Journal of Applied Geophysics, 34(3): 199-212. DOI:10.1016/0926-9851(95)00016-X |
Tezkan B, Hördt A, Gobashy M. 2000. Two-dimensional radiomagnetotelluric investigation of industrial and domestic waste sites in Germany. Journal of Applied Geophysics, 44(2-3): 237-256. DOI:10.1016/S0926-9851(99)00014-2 |
Tezkan B, Georgescu P, Fauzi U. 2005. A radiomagnetotelluric survey on an oil-contaminated area near the Brazi Refinery, Romania. Geophysical Prospecting, 53(3): 311-323. DOI:10.1111/gpr.2005.53.issue-3 |
Tezkan B, Saraev A. 2008. A new broadband radiomagnetotelluric instrument:Applications to near surface investigations. Near Surface Geophysics, 6(4): 245-252. DOI:10.3997/1873-0604.2008019 |
Turberg P, Müller I, Flury F. 1994. Hydrogeological investigation of porous environments by radio magnetotelluric-resistivity (RMT-R 12-240 kHz). Journal of Applied Geophysics, 31(1-4): 133-143. DOI:10.1016/0926-9851(94)90052-3 |
Wu X P, Liu Y, Wang W. 2015. 3D resistivity inversion incorporatingtopography based on unstructured meshes. Chinese Journal of Geophysics (in Chinese), 58(8): 2706-2717. DOI:10.6038/cjg20150808 |
Yuan Y, Tang J T, Ren Z Y, et al. 2015. Two-dimensional complicated radio-magnetotelluric finite-element modeling using unstructured grids. Chinese Journal of Geophysics (in Chinese), 58(12): 4685-4695. DOI:10.6038/cjg20151229 |
韩骑, 胡祥云, 程正璞, 等. 2015. 自适应非结构有限元MT二维起伏地形正反演研究. 地球物理学报, 58(12): 4675-4684. DOI:10.6038/cjg20151228 |
吴小平, 刘洋, 王威. 2015. 基于非结构网格的电阻率三维带地形反演. 地球物理学报, 58(8): 2706-2717. DOI:10.6038/cjg20150808 |
原源, 汤井田, 任政勇, 等. 2015. 基于非结构化网格的任意复杂2DRMT有限元模拟. 地球物理学报, 58(12): 4685-4695. DOI:10.6038/cjg20151229 |