2019, Vol. 62
2. 东华理工大学放射性地质与勘探技术国防重点学科实验室, 南昌 330013
, DENG JuZhi1,2
, CHEN Hui1,2, CHEN Xiao1,2, WANG XianXiang1,2, ZHANG ZhiYong1,2, YE YiXin1,2, CHEN ShuNi1,2
2. Key Laboratory of Radioactive Geology and Exploration Technology Fundamental Science for National Defense, East China University of Technology, Nanchang 330013, China
大地电磁法(Magnetotelluric，简称MT)是一种天然电磁场勘探方法(Tikhonov, 1950；Cagniard, 1953)，它具有穿透深度大和不受高阻屏蔽影响等优点，是构造研究(魏文博等，2003；赵国泽等，2004；Bedrosian and Feucht, 2014；Meqbel et al., 2014；Dong et al., 2016；Yin et al., 2017)、石油勘探(Mansoori et al., 2015；Sarvandani et al., 2017)、矿产及地热勘察(Zhang et al., 2015；邓居智等，2015；Rosenkjaer et al., 2015；Mohan et al., 2017)的有效方法，同时还可以结合地震数据(Roux et al., 2011；陈晓等，2011；Beka et al., 2016)进行综合地质解释.近年来，计算机技术的迅速发展促使大地电磁场二维正演算法的研究趋于成熟；且当前大地电磁测深数据反演仍然以二维反演为主.得益于众多学者在大地电磁场三维正反演研究中的贡献，逐步实现了大地电磁场三维反演(Avdeev, 2005；Siripunvaraporn, 2012)，但在完全实现大地电磁场三维反演实用化的进程中仍然存在诸多挑战(Börner, 2010)，其中，起伏地形条件下大地电磁场三维正反演算法的研究尤为关键.
在地球电磁法勘探中，近地表起伏地形会引起地表电流聚集而产生电场畸变，给地形起伏剧烈的山区大地电磁数据反演解释造成极大困难(Jiracek, 1990；Vallianatos, 2002).为了解决该难题，通常采用畸变校正和带地形反演两种方式，其中畸变校正虽然可以在一定程度上消除地形引起的静位移现象，但也可能造成反演模型中异常体几何形态的畸变，难以完全解决电流畸变带来的问题(Gürer and Ìlkşık, 1997).而带地形三维反演能够真实模拟实际地形起伏情况，从根本上解决地形影响问题，成为近年来主要的发展方向.在三维正演中，对起伏地形的模拟和逼近极为关键，求解区域的剖分方法包括规则网格和非规则网格两种.规则网格剖分将空气和地下介质当成一个整体求解区域，利用正六面体进行剖分并在地形以上的网格单元中填充一个极大电阻率值，通常应用于有限差分或者有限体积算法中(Chen et al., 1998；Haber and Ruthotto, 2014)；非规则网格利用四面体网格对起伏地形进行逼近，通常应用于有限单元法(Wannamaker et al., 1986；徐世浙等，1997；Nam et al., 2008；曹晓月等，2018)中，该方法更能适应各种复杂地形且模拟准确度更强，尤其是自适应网格剖分技术(Ren et al., 2013；Grayver, 2015；Jahandari and Farquharson, 2017；Yoshimura et al., 2018).但由于非规则网格剖分的计算量巨大，反演稳定性不足，仍难应用于三维反演中(Usui, 2015).在大地电磁三维反演中，各种不同的常用最优化算法都得到应用，如OCCAM法(Siripunvaraporn et al., 2005)、非线性共轭梯度法NLCG(Newman and Alumbaugh, 2002；董浩等，2014；Zhang et al., 2017)、L-BFGS法(Avdeev and Avdeeva, 2009；阮帅，2015；秦策等，2017)、高斯牛顿法(GN)(Haber et al., 2004；Kordy et al., 2016；Yang and Oldenburg, 2016)等.其中L-BFGS法和NLCG法因内存需求最少成为大地电磁法三维反演的常用方法，但L-BFGS算法的下降方向更接近牛顿法，又无需求解和存储Hessian矩阵，从而大量缩短了搜索时间，因此更适合求解大规模的三维反演问题(刘云鹤和殷长春，2013).但是，受正则化因子选择、带地形的协方差矩阵的处理、步长搜索策略等因素影响，造成L-BFGS算法稳定性不足，在大地电磁三维反演中并未得到广泛的应用(Newman and Boggs, 2004).
为此，本文开展基于L-BFGS法的带地形大地电磁三维反演研究.设计一种既能保证空气电阻率固定不变又能保证模型平滑约束的协方差矩阵统一表达式以解决带地形反演的问题.同时在反演算法中采用正则化因子冷却法以及基于Wolf条件的步长搜索策略，提升反演的稳定性.通过对一系列带地形理论模型正演合成数据进行带地形三维反演试算，并将其结果与开源的非线性共轭梯度反演结果进行对比，验证了本文算法的正确性和可靠性.最后，利用本文开发的算法对华南某山区大地电磁实测数据进行带地形三维反演，检验了该算法的实用性.
1 基于L-BFGS的大地电磁三维反演众所周知，反演的基础和前提是快速高效的正演模拟，本文选用常用的规则网格剖分策略，在地表以上空气层网格中填充一个极大电阻率的方式实现起伏地形逼近，并采用有限体积法(Haber and Ascher, 2001；Haber and Ruthotto, 2014)对双旋度电场方程进行离散，结合第一类边界条件得到庞大的稀疏复数线性方程组；通过采用代数多重网格算法——聚集多重网格算法(AGMG)(陈辉等，2018)对该线性方程组进行迭代求解，并且在迭代求解过程中施加散度校正技术(Henson and Yang, 2002；Liu and Yin, 2013)来减小迭代次数，实现大地电磁三维正演模拟.在此基础之上，本文开展了基于L-BFGS法的带地形大地电磁三维反演的研究.
1.1 大地电磁三维反演目标函数构建一般而言，地球物理反演是根据观测到的地球物理数据推测地球内部介质的物性结构或空间变化.根据Backus-Gilbert反演理论可知，地球物理反演问题是非唯一的病态问题.为了减少反演的非唯一性和提高稳定性，通常采用Tikhonov正则化反演思想，在目标函数中施加模型约束项.本文采用平滑约束模型，则大地电磁反演问题可以归结为求解如下目标函数的极小值：
|
(1) |
其中d=(d1, d2, …, dN)T，为数据空间，N为观测数据大小；模型空间表示为m=(m1, m2, …, mM)T, M为模型参数大小；f(m)为正演算子，表示模型空间m到数据空间d的映射，采用有限体积算法实现；Cd是数据协方差矩阵，为一个与数据误差相关的对角矩阵；m0为初始模型参数；Cm是模型协方差矩阵，约束模型光滑度；λ为正则化因子，是平衡数据拟合项和模型约束项的一个参数.为了简化目标函数的表达形式，消除数据协方差矩阵Cd和模型协方差矩阵Cm.根据Egbert和Kelbert(2012)提出的变换方法，将模型空间m及数据空间d和正演算子f(m)分别变换为Cm-1/2(m-m0)、Cd-1/2d和Cd-1/2f(m)，并记为
|
(2) |
由(2)式可知，目标函数变成标准的数据拟合项和模型约束项的
根据最优化理论思想，反演问题可以归纳为给定一个初始模型值
|
(3) |
其中ak为每次最优化迭代的搜索步长，通常以模型更新量
|
(4) |
其中，
|
(5) |
|
(6) |
|
(7) |
由上式可以看出，迭代过程中无需存储近似Hessian矩阵逆矩阵
|
(8) |
其中
|
(9) |
γk是估计搜索方向的比例因子，能够确保适当的搜索方向pk.为了得到模型更新量
对于模型更新量
|
(10) |
从(10)式可以发现，关键问题就是计算数据灵敏度矩阵J.对于大规模问题的三维反演，一方面通过计算数据灵敏度矩阵J或JT与某个向量乘积的方法避免数据灵敏度矩阵J存储；另一方面采用伴随正演实现数据灵敏度矩阵J或JT的快速计算，该方法只需进行两次正演就能够得到目标函数梯度G.本文通过计算灵敏度矩阵J与数据空间向量的乘积JTr来实现灵敏度矩阵的求解.
1.3 模型协方差在目标函数(1)式的模型约束项中包含一个模型协方差矩阵Cm.该矩阵通常既可以用来衡量反演迭代模型与初始模型变化程度，也可以当作一个模型平滑约束算子.本文考虑模型光滑约束，首先引入一个扩散系数γ，描述各个方向模型光滑程度；然后采用递归自回归方法定义(Egbert and Kelber, 2012)，其Cm-1/2表达式为：
|
(11) |
其中n为光滑次数；Cx、Cy和Cz为三个不同方向扩散系数γ构成的矩阵；Σ为一个对角矩阵，是一个模型空间变化对角矩阵，用来控制每个模型单元的变化率.在起伏地形反演中，空气中每个模型单元电阻率保存为一个极大值固定不变，Σ矩阵主对角元素取值为0；地下每个模型单元电阻率变化率由三个方向扩散系数的乘积决定.Σ矩阵主对角元素为
|
(12) |
从式(3)和式(12)可以看出，本文设计的协方差矩阵既解决了起伏地形处理问题，又保证了模型的平滑约束，建立了起伏地形下三维反演中的模型协方差矩阵统一表达式.
1.4 算法流程为了直观展示L-BFGS法反演的实现过程，本文给出L-BFGS法反演的具体流程(图 1).从图中可以看出，每次迭代只需计算一次梯度和若干次目标函数，但需存储m个向量对(sk, yk).对于正则化因子λ的选取采用冷却法.首先给定较大正则化因子λ0，冷却因子ρ>1，截断拟合差η，最小正则化因子λmin，然后根据前后两次数据拟合差的变化是否小于截断拟合差η，若成立则正则化因子进行冷却λk=λk-1/ρ.重复上述过程直到满足反演条件或达到正则化因子下界λmin结束反演.另外采用基于Wolf条件下的步长搜索策略，它由充分下降条件(Armijo条件)和曲率条件两个构成，能够在每次迭代中选择一个合适的步长ak.
|
图 1 L-BFGS法反演流程示意图 Fig. 1 Flow chart of L-BFGS inversion |
为了验证本文所实现的带地形大地电磁三维反演算法的适用性并分析L-BFGS算法的优缺点，本文设计几个带地形的地电模型(单个异常体模型，“棋盘”模型)，并与开源大地电磁三维反演程序ModEM反演结果对比.为了更加真实的模拟野外实测情况，所有理论数据加入3%高斯白噪声.本文设计的所有模型反演工作均在小型服务器上完成，该服务器为两个Intel(R) Xeon(R) CPU E5-2620 v2 @ 2.10 GHz处理器，64 GB内存，64位Red Hat Enterprise Linux Server Release 6.4操作系统.算法开发和运行平台为Intel(R) Visual Fortran 13.1.1，MPI并行库采用Intel(R) MPI Library 4.1，所有理论数据反演都采用12个进程的MPI并行程序进行.
2.1 带山峰地形的单个异常体模型图 2为一个山峰地形下单个异常体模型.山峰地形及异常体平面视图和测点分布情况如图 2a所示，图中小圆圈表示测点，AB表示沿X-Z方向(Y=2 km)横穿山峰地形及异常体的测线；图 2b所示为该理论模型X-Z方向(Y=0 km)剖面图；图 2c所示为异常体三维等值面视图，异常体为产状水平的板状体.本文将山峰地形设计为一个上底边长为5 km，下底边长为20 km，山峰高度距地面1 km的正四棱台状模型.地下围岩电阻率被设置为100 Ωm，而空气中的电阻率固定为1010Ωm.为了考虑不同导电性异常体的响应情况，分别设计了低阻异常体理论模型和高阻异常体理论模型.因此，板状异常体电阻率为10 Ωm或4000 Ωm，该异常体顶部距离地面为3 km，其几何尺度在x、y、z方向分别是11 km、11 km和3 km.大地电磁测量区域位于异常体上方，如图 2a所示，布设等距分布的6条主测线，线距为4 km，此外还增设了2条旁测线，共计8条测线；每条测线上15个测点，共计120个测点.测量数据集为全阻抗(Zxx, Zxy, Zyy, Zyx)数据，大地电磁测量观测信号频率范围为100~0.01 Hz，并按对数等间隔取23个测量频率，共计11040个观测数据.
|
图 2 三维山峰地形及异常体位置示意图 (a)山峰地形及异常体平面视图和测点分布情况；(b)山峰地形下高/低阻异常模型X-Z方向剖面图；(c)山峰地形下异常体三维等值面视图. Fig. 2 Schematic illustrations of single anomaly model with topographic relief (a) X-Y plane view of topography and single anomaly model with measurement points; (b) X-Z section view of the conductive or resistive anomaly model; (c) 3D view of anomaly body with topographic relief. |
反演网格采用不等距剖分，网格大小为37×37×52，共计71188个单元.反演模型初始电阻率设置为100 Ωm，初始正则化因子λ0设为10，冷却因子ρ设为10，最小正则化因子λmin设为10-3，模型协方差三个方向的平滑因子设置为0.2，目标拟合差RMS设为0.5，L-BFGS法存储向量对个数m设为6.低阻理论模型反演计算经过32次迭代后，目标拟合差RMS从初始的1.94下降到0.69(图 3)，耗时211 min；高阻理论模型反演计算经过21次迭代后，目标拟合差RMS从初始的1.01下降到0.69(图 3)，历时79 min.图 4所示为沿AB剖面上周期为0.01 s、0.12 s、1 s、3.52 s、8.11 s、100 s时各测点反演前后视电阻率和相位响应拟合情况.图中散点表示反演前的数据响应，曲线表示反演后的数据响应，散点和曲线的填充色由冷色(蓝色)转为暖色(黑色)表示由短周期过度到长周期.对比结果显示无论是XY模式还是YX模式，L-BFGS反演前后视电阻率和相位响应在整个周期内(0.01 s to 100 s)均吻合良好，表明计算结果是可靠的.
|
图 3 三维山峰地形下单个异常体理论模型反演迭代拟合差曲线 Fig. 3 Misfit curves of iterative inversion for single anomaly model with topographic relief |
|
图 4 三维山峰地形下单个异常体理论模型反演前后数据拟合对比曲线 Fig. 4 MT original and predicted data from section A-B in the single anomaly model with topographic relief |
图 5为三维山峰地形下单个异常体理论模型反演结果，其中，图 5a是低阻异常理论模型的平地形三维反演结果X-Z方向(Y=0 km)剖面图，图 5b是带地形三维反演结果X-Z方向(Y=0 km)剖面图，图 5c是三维等值面视图.从中不难发现，采用带地形三维反演计算方法所获得的反演模型，其导电性结构特征与真实模型结构基本吻合，而且异常体周围没有出现冗余的虚假异常(图 5b).而不带地形反演结果虽然也能大致反映低阻异常体的轮廓和相对准确的埋深，但其对围岩电阻率的复原程度较差，在异常体周围存在虚假异常，尤为显著的是山峰地形边坡下出现的高阻假异常.对于山峰地形下存在高阻异常体的理论模型，其反演结果如图 5d—f所示.从中可以看出，不带地形的三维反演结果难以反映出真实模型中异常体的结构特征，异常体响应与山峰地形下虚假异常连成一片，表现为一个长条形的高阻异常带，其轮廓与模型存在较大差异；而带地形三维反演结果基本与模型一致，但高阻异常体边界范围和数值相对偏小，这是由于高阻异常分辨率差引起的.试算结果表明本文开发的起伏地形下三维大地电磁L-BFGS反演代码是正确可行的.
|
图 5 三维山峰地形下单个异常体理论模型反演结果 (a)低阻异常模型平地形三维反演X-Z方向剖面图；(b)低阻异常模型带地形三维反演X-Z方向剖面图；(c)低阻异常模型带地形三维反演结果等值面视图；(d)高阻异常模型平地形三维反演X-Z方向剖面图；(e)高阻异常模型带地形三维反演X-Z方向剖面图；(f)高阻异常模型带地形三维反演结果等值面视图. Fig. 5 Inversion results of data on single anomaly model with topographic relief (a) X-Z section of conductive anomaly model by inversion without topography; (b) X-Z section of conductive anomaly model by inversion with topography; (c) 3D view of the conductive anomaly model with iso-surfaces of 30 Ωm; (d) X-Z section of resistive anomaly model by inversion without topography; (e) X-Z section of resistive anomaly model by inversion with topography; (f) 3D view of resistive anomaly model with iso-surfaces of 250 Ωm. |
为了进一步评价本文反演算法的有效性和稳定性，设计了一个三维峰-谷组合地形下高-低阻交错分布的“棋盘”模型，其中图 6a所示为三维峰-谷组合地形及异常体平面视图和测点分布情况，图中测点用小圆圈表示，MN表示沿X-Z方向(Y=-10 km)横穿峰-谷地形和高-低阻异常体的测线；图 6b和6c所示分别为该理论模型X-Z方向(Y=-10 km)剖面示意图和三维等值面视图.其中山峰地形同上文所述山峰地形模型结构一致，山谷地形则设计为一个下底边长为5 km，上底边长为20 km，山谷深度距地面1 km的倒四棱台状模型；地下围岩和空气中的电阻率分别为100 Ωm和1010 Ωm.图中共分布4个尺寸完全相同的水平产状异常体.其中，高阻和低阻异常体电阻率分别为2000 Ωm和10 Ωm，各个异常体之间相距5 km，顶部距地面均为3 km，其几何尺度在x、y、z方向分别是15 km、15 km和3 km.大地电磁测量区域位于异常体上方，如图 6a所示，布设等距分布的12条主测线，线距为4 km；每条测线上25个测点，共计300个测点.观测数据集、观测信号频率范围、测量频点数与前文所述原则一致，本次观测数据共计27600个.反演模型网格剖分大小为61×61×58，共计215818个单元.在对该理论模型进行反演时，除最小正则化因子λmin设为10-6以外，其他参数与前文所述山峰地形下单个异常体理论模型反演时一致.对该理论模型采用本文研发的L-BFGS算法进行反演计算，经过79次迭代后，目标拟合差RMS从初始的2.67下降到0.72(图 7)，耗时671 min.相同条件下，利用开源的NLCG算法进行试算，经100次迭代后目标拟合差RMS从2.67衰减到0.69(图 7)，全程耗时994 min.图 8所示为沿MN剖面上周期为0.01 s、0.12 s、1 s、3.52 s、8.11 s、100 s时各测点反演前后视电阻率和相位响应拟合情况.图中散点表示反演前的数据响应，实线和虚线分别表示L-BGFS法和NLCG法反演后的数据响应，，颜色由冷色(蓝色)转为暖色(黑色)表示由短周期过度到长周期.对比结果显示，对于短周期(0.01 s to 0.12 s)数据，两种算法反演前后的视电阻率响应均吻合良好，相位响应在山峰和山谷处有细微的差别，对于长周期(1 s to 100 s)数据，两种反演方法的视电阻率和相位响应均吻合良好，表明计算结果是可靠的.同时可以发现，L-BGFS法和NLCG法反演前后拟合效果相当.
|
图 6 三维峰-谷组合地形下“棋盘”模型示意图 (a)峰-谷组合地形及异常体平面视图和测点分布情况；(b)峰-谷组合地形下“棋盘”模型X-Z方向剖面图；(c)峰-谷组合地形下“棋盘”模型三维等值面视图. Fig. 6 Schematic illustrations of resistive-conductive anomalies model with peak-valley topography (a) X-Y plane view of peak-valley topography and anomaly model with measuring points; (b) X-Z section view of the resistive-conductive anomaly model; (c) 3D view of the resistive-conductive anomaly with iso-surfaces(10 Ωm and 2000 Ωm) |
|
图 7 三维峰-谷组合地形下“棋盘”模型反演迭代拟合差曲线 Fig. 7 Misfit curves of. Iterative inversion on resistive- conductive anomaly model with peak-valley topography |
|
图 8 三维峰-谷组合地形下“棋盘”模型反演前后数据拟合对比曲线 Fig. 8 MT original and predicted data from section M-N in the resistive-conductive anomaly model with peak-valley topography |
图 9为三维峰-谷地形下“棋盘”模型反演结果，图 9A是X-Y方向(Z=-5 km)剖面图，图 9B是X-Z方向(Y=-10 km)剖面图，图 9C是Y-Z方向(X=10 km)剖面图，图 9D是三维等值面视图；从左到右分别是真实模型、开源三维非线性共轭梯度法反演结果和采用本文研发的起伏地形下三维反演算法获得的结果.从中可以发现，两种反演计算方法所获得的三维反演模型，其导电性结构特征都能够与真实模型吻合良好，而且异常体周围没有出现冗余的虚假异常.此外，相较于前文所述山峰地形下高阻异常体模型的反演结果，本次反演模型对山谷地形下高阻异常体边界的刻画更为清晰，高阻异常体电阻率和围岩电阻率之间的差异更接近真实模型，高阻体电阻率复原程度更佳.对比这两种反演方法获得的结果，可以发现，本文提出的起伏地形下三维大地电磁L-BFGS反演方法能够达到开源三维NLCG反演方法的同等效果，一定程度上改善了电流畸变现象对反演结果造成的偏差，使得反演结果更接近真实的地下导电性结构特征.
|
图 9 三维峰-谷组合地形下“棋盘”模型反演结果 Fig. 9 Inversion results of data on resistive-conductive anomaly model with peak-valley topography |
本文在总进程数为24的MPI并行程序上，分别采用NLCG法和本文开发的L-BFGS算法对该组理论模型进行反演试算，反演结果表明，L-BFGS方法79次迭代过程中进行了80次正演计算，而NLCG方法100次迭代过程中需要191次正演计算.将23个频点分配到不同的节点上进行独立地正演计算，每次计算仅包含300个位置的子正演，统计结果表明，L-BFGS方法每次计算大约耗时83.42 s，整个计算耗时为111.22 min，而NLCG方法每次计算则需要耗时87.2 s，计算总耗时为277.62 min.虽然二者反演过程中所需内存相当(约1 G)，但在同等目标拟合差的设定下，L-BFGS法相较NLCG法能够相对节省近一半的计算时间.
3 华南某山区实测数据带地形三维反演我国幅员辽阔，自西向东地形起伏变化大，特别是在我国的华南山区，地形复杂程度更高，而在这些地区进行的大地电磁剖面观测难免跨过多个起伏地形单元.因此，为了提高类似地区复杂地形下大地电磁数据的解释精度，实现带地形大地电磁三维反演显得尤为迫切.将本文开发的L-BFGS大地电磁带地形三维反演算法应用于我国华南某山区的大地电磁数据反演，并对反演结果进行了地质解释，检验了本文算法的实用性和适用性.
华南某山区实测大地电磁数据测点位置如图 10b所示，图中共两条大地电磁测深剖面(GG1和GG2)，合计81个测点.反演数据集选用非对角阻抗(Zxy、Zyx)数据，为了提高反演效率，对观测数据集在320 Hz~0.1 Hz频率范围内进行对数等间距抽稀为23个频点，网格剖分大小为82×60×90.根据观测数据的视电阻率和一维反演结果，将初始模型电阻率设置为100 Ωm.初始正则化因子λ0设置为10，冷却因子ρ为10，最小正则化因子λmin为10-7.本次反演计算经215次迭代后，RMS从10.27衰减为1.82，如图 10d中RMS衰减曲线所示.图 10c所示为部分测点大地电磁观测数据与L-BFGS法三维反演模型响应阻抗和相位曲线的拟合情况，图中显示了测区7个测点的观测数据和反演模型响应结果对比情况，对比结果表明，两种极化模式下的阻抗和相位响应曲线都与观测数据吻合较好.图 10a以GG1剖面为例展示了该剖面实测数据与反演模型响应的拟断面对比情况，包括阻抗和相位的拟断面.从图中可以发现，不论是XY极化模式还是YX极化模式，反演模型响应的阻抗拟断面反映的高阻和低阻异常位置与观测数据揭示的情况基本一致，相位响应除了低频部分存在数值上的偏差外，整体上与观测数据趋同，异常的区域位置反映良好.总体上拟合情况较为理想，表明L-BFGS大地电磁带地形三维反演结果是正确可靠的.
|
图 10 实测数据视电阻率和相位与三维反演模型响应对比及误差衰减曲线和测点布置图 (a) GG1剖面大地电磁观测数据和反演模型响应对比图；obs表示实测数据，pred表示反演模型响应数据；(b)测点布置图；(c)观测数据与反演模型响应视电阻率和相位曲线拟合情况；(d)反演拟合差结果. Fig. 10 Comparison of apparent resistivity and phase between observed data and the results of 3D inversion, measurement sites and fitting error curves (a) Pseudo-sections of observed and modelled data of 3D inversion for MT sites along profile GG1, obs: observed data, pred: modelled data; (b) Arrangement of measurement sites; (c) Comparison of apparent resistivity and phase curves between observed and the responses of 3D inversion model for 7 sites; (d) Fitting errors of inversion. |
图 11a所示为测区GG1剖面大地电磁三维反演断面图.图上左侧为地理西方向，从图中可以看出，该反演剖面深度约6 km，自西向东大致上可以分出5个高阻异常和4个低阻异常电性区块，分别用R1—R5和C1—C4表示，根据上述电性变化，在该剖面上推测出5条断层，图中用F1a—F1e表示.从图 11a中可以发现，剖面西侧的R1和山峰范围深部的R3、R4以及剖面东侧的R5高阻异常区间的电阻率值都在5000 Ωm以上，结合其电性特征和测区地表地质资料，推测为青白口系库里组Qbk变质岩基底；而山峰范围浅部的高阻异常区间R2与其下方低阻异常区间C3形成两套不同电性分布的层位，鉴于R2的电性特征和地表出露地质单元，推测解释为K1e火山岩，其主体为碎斑熔岩，而C3则可以解释为不整合接触面和以流纹英安岩为主的K1d火山岩；位于剖面西侧的低阻异常区间C1和C2，其电阻率值均在100 Ωm以下，其中浅部的C2低阻异常区间所揭示的空间位置刚好与K2红盆范围重合，深部的C1低阻异常空间上表现为从上往下逐渐扩大的趋势，推测深部可能存在一条隐伏构造断裂带；该剖面东侧的低阻异常区间C4与地表地质资料揭示的Qbk上段变质岩电性特征吻合，其主体为泥质千枚岩类.图 11b所示为测区GG2剖面大地电磁三维反演的断面图.图中左侧为地理南方向，根据电性差异划分出5个高阻异常和2个低阻异常电性区块，分别用R1—R5和C1—C2表示，F2a—F2d为推测断层.其中标高-2 km以下的高阻异常区间推测为Qbk变质岩基底，山峰范围浅部的R2则推测为K1e火山岩；而低阻异常区间C1和C2与GG1剖面反演结果断面图上C3和C4低阻异常区间情况类似，综合其电性特征和地表地质资料，C1解释为不整合接触面和K1d火山岩，C2推测为Qbk上段变质岩.
|
图 11 华南某山区实测大地电磁测深数据L-BFGS三维反演结果 (a) GG1剖面大地电磁观测数据三维反演结果；(b) GG2剖面大地电磁观测数据三维反演结果；(c)三维视图. Fig. 11 3D L-BFGS inversion results with topography of MT data in a South China mountainous area (a) 3D inversion results for MT sites along profile GG1; (b) 3D inversion results for MT sites along profile GG2; (c) 3D view of slices of inversion results. |
综上所述，华南某山区大地电磁数据采用L-BFGS带地形三维反演获得的结果可以较好地反映出该区的电性结构.揭示了研究区以高阻介质为基底，中间以低阻不整合面和相对低阻介质连续分布，浅部覆盖高阻介质的电性结构特征.反演模型能够反映出研究区地下地质体的电阻率变化和空间展布，与该区地表地质资料和图切地质剖面相对应，表明本文提出的带地形三维反演方法可以适用于复杂地形下实测大地电磁数据的带地形三维反演.本次反演获得的电性结构成功鉴别出该地区的基底构成，初步查明已知断裂在深部的接触关系，一定程度上有助于深化认识研究区成矿作用，为深部找矿突破提供了新的理论依据.
4 结论本文实现了一种基于L-BFGS法的带地形大地电磁三维反演算法，通过对一系列理论数据和实测数据进行反演试算，并与已有大地电磁三维反演程序(ModEM)对比，验证了本文反演算法的可靠性和实用性.通过本文研究得到以下几点认识：
(1) 对于地形起伏较大时，带地形三维反演基本能够恢复理论模型的基本形态和电阻率值，因此，在实际资料处理中尽量选择带地形反演，以免造成解释和推断的偏差.
(2) 将L-BFGS法成功应用在华南某山区大地电磁实测数据的三维带地形反演中，查明了研究区的电性结构特征，表明该方法可以满足复杂地形条件下大地电磁数据带地形三维反演的需求.
(3) 本文实现的基于L-BFGS法的大地电磁带地形三维反演算法，其反演适用性和效果基本与广泛应用于实际中的非线性共轭梯度算法(NLCG)相当.但对于地形起伏剧烈的情况下，本文采用的有限体积正演算法需要对地形起伏层进行精细剖分，造成反演效率急剧下降，有待进一步改进.
致谢 作者在此感谢Anna Kelbert和Gary D. Egbert教授提供其三维反演代码ModEM.感谢审稿专家提供的宝贵意见.
Avdeev D, Avdeeva A. 2009. 3D magnetotelluric inversion using a limited-memory quasi-Newton optimization. Geophysics, 74(3): F45-F57. DOI:10.1190/1.3114023 |
Avdeev D B. 2005. Three-dimensional electromagnetic modelling and inversion from theory to application. Surveys in Geophysics, 26(6): 767-799. DOI:10.1007/s10712-005-1836-x |
Bedrosian P A, Feucht D W. 2014. Structure and tectonics of the northwestern United States from EarthScope USArray magnetotelluric data. Earth and Planetary Science Letters, 402: 275-289. DOI:10.1016/j.epsl.2013.07.035 |
Beka T I, Smirnov M, Birkelund Y, et al. 2016. Analysis and 3D inversion of magnetotelluric crooked profile data from central Svalbard for geothermal application. Tectonophysics, 686: 98-115. DOI:10.1016/j.tecto.2016.07.024 |
Börner R U. 2010. Numerical modelling in geo-electromagnetics:advances and challenges. Surveys in Geophysics, 31(2): 225-245. |
Cagniard L. 1953. Basic theory of the magneto-telluric method of geophysical prospecting. Geophysics, 18(3): 605-635. DOI:10.1190/1.1437915 |
Cao X Y, Yin C C, Zhang B, et al. 2018. A goal-oriented adaptive finite-element method for 3D MT anisotropic modeling with topography. Chinese Journal of Geophysics (in Chinese), 61(6): 2618-2628. DOI:10.6038/cjg2018L0068 |
Chen H, Yin M, Yin C C, et al. 2018. Three-dimensional magnetotelluric modelling using aggregation-based algebraic multigrid method. Journal of Jilin University(Earth Science Edition) (in Chinese), 48(1): 261-270. |
Chen P F, Hou Z Z, Fan G H. 1998. Three-dimensional topographic responses in MT using finite difference method. Acta Seismologica Sinica (English Edition), 11(5): 631-635. DOI:10.1007/s11589-998-0079-6 |
Chen X, Yu P, Zhang L L, et al. 2011. Adaptive regularized synchronous joint inversion of MT and seismic data. Chinese Journal of Geophysics (in Chinese), 54(10): 2673-2681. DOI:10.3969/j.issn.0001-5733.2011.10.024 |
Deng J Z, Chen H, Yin C C, et al. 2015. Three-dimensional electrical structures and significance for mineral exploration in the Jiujiang-Ruichang District. Chinese Journal of Geophysics (in Chinese), 58(12): 4465-4477. DOI:10.6038/cjg20151211 |
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 |
Gürer A, Ìlkışık O M. 1997. The importance of topographic corrections on magnetotelluric response data from rugged regions of anatolia. Geophysical Prospecting, 45(1): 111-125. DOI:10.1046/j.1365-2478.1997.3160236.x |
Grayver A V, Kolev T V. 2015. Large-scale 3D geoelectromagnetic modeling using parallel adaptive high-order finite element method. Geophysics, 80(6): E277-E291. DOI:10.1190/geo2015-0013.1 |
Haber E, Ascher U M. 2001. Fast finite volume simulation of 3D electromagnetic problems with highly discontinuous coefficients. SIAM Journal on Scientific Computing, 22(6): 1943-1961. DOI:10.1137/S1064827599360741 |
Haber E, Ascher U M, Oldenburg D W. 2004. Inversion of 3D electromagnetic data in frequency and time domain using an inexact all-at-once approach. Geophysics, 69(5): 1216-1228. DOI:10.1190/1.1801938 |
Haber E, Ruthotto L. 2014. A multiscale finite volume method for Maxwell's equations at low frequencies. Geophysical Journal International, 199(2): 1268-1277. DOI:10.1093/gji/ggu268 |
Henson V E, Yang U M. 2002. BoomerAMG:A parallel algebraic multigrid solver and preconditioner. Applied Numerical Mathematics, 41(1): 155-177. |
Jahandari H, Farquharson C G. 2017. 3-D minimum-structure inversion of magnetotelluric data using the finite-element method and tetrahedral grids. Geophysical Journal International, 211(2): 1189-1205. DOI:10.1093/gji/ggx358 |
Jiracek G R. 1990. Near-surface and topographic distortions in electromagnetic induction. Surveys in Geophysics, 11(2-3): 163-203. DOI:10.1007/BF01901659 |
Kordy M, Wannamaker P, Maris V, et al. 2016. 3-dimensional magnetotelluric inversion including topography using deformed hexahedral edge finite elements and direct solvers parallelized on symmetric multiprocessor computers-Part Ⅱ:direct data-space inverse solution. Geophysical Journal International, 204(1): 94-110. DOI:10.1093/gji/ggv411 |
Liu Y H, Yin C C. 2013. 3D inversion for frequency-domain HEM data. Chinese Journal of Geophysics (in Chinese), 56(12): 4278-4287. DOI:10.6038/cjg20131230 |
Liu Y H, Yin C C. 2013. Electromagnetic divergence correction for 3D anisotropic EM modeling. Journal of Applied Geophysics, 96: 19-27. DOI:10.1016/j.jappgeo.2013.06.014 |
Mansoori I, Oskooi B, Pedersen L B. 2015. Magnetotelluric signature of anticlines in Iran's Sehqanat oil field. Tectonophysics, 654: 101-112. DOI:10.1016/j.tecto.2015.05.004 |
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. DOI:10.1016/j.epsl.2013.12.026 |
Mohan K, Kumar G P, Chaudhary P, et al. 2017. Magnetotelluric investigations to identify geothermal source zone near Chabsar hotwater spring site, Ahmedabad, Gujarat, Northwest India. Geothermics, 65: 198-209. DOI:10.1016/j.geothermics.2016.10.001 |
Nam M J, Kim H J, Song Y, et al. 2008. Three-dimensional topography corrections of magnetotelluric data. Geophysical Journal International, 174(2): 464-474. DOI:10.1111/j.1365-246X.2008.03817.x |
Newman G A, Alumbaugh D L. 2002. Three-dimensional magnetotelluric inversion using non-linear conjugate gradients. Geophysical Journal International, 140(2): 410-424. |
Newman G A, Boggs P T. 2004. Solution accelerators for large-scale three-dimensional electromagnetic inverse problems. Inverse Problems, 20(6): S151-S170. DOI:10.1088/0266-5611/20/6/S10 |
Qin C, Wang X B, Zhao N. 2017. Parallel three-dimensional forward modeling and inversion of magnetotelluric based on a secondary field approach. Chinese Journal of Geophysics (in Chinese), 60(6): 2456-2468. DOI:10.6038/cjg20170633 |
Ren Z Y, Kalscheuer T, Greenhalgh S, et al. 2013. A goal-oriented adaptive finite-element approach for plane wave 3-D electromagnetic modelling. Geophysical Journal International, 194(2): 700-718. DOI:10.1093/gji/ggt154 |
Rosenkjaer G K, Gasperikova E, Newman G A, et al. 2015. Comparison of 3D MT inversions for geothermal exploration:Case studies for Krafla and Hengill geothermal systems in Iceland. Geothermics, 57: 258-274. DOI:10.1016/j.geothermics.2015.06.001 |
Roux E, Moorkamp M, Jones A G, et al. 2011. Joint inversion of long-period magnetotelluric data and surface-wave dispersion curves for anisotropic structure:Application to data from Central Germany. Geophysical Research Letters, 38(5): L05304. DOI:10.1029/2010GL046358 |
Ruan S. 2015. Three-dimensional Magnetotelluric Inversion Based on Limited Memory Quasi-Newton Method (Doctorate thesis). Chengdu University of Technology.
|
Sarvandani M M, Kalateh A N, Unsworth M, et al. 2017. Interpretation of magnetotelluric data from the Gachsaran oil field using sharp boundary inversion. Journal of Petroleum Science and Engineering, 149: 25-39. DOI:10.1016/j.petrol.2016.10.019 |
Siripunvaraporn W, Egbert G, Lenbury Y, et al. 2005. Three-dimensional magnetotelluric inversion:data-space method. Physics of the Earth and Planetary Interiors, 150(1-3): 3-14. DOI:10.1016/j.pepi.2004.08.023 |
Siripunvaraporn W. 2012. Three-dimensional magnetotelluric inversion:an introductory guide for developers and users. Surveys in Geophysics, 33(1): 5-27. |
Tikhonov A N. 1950. On determining electrical characteristics of the deep layers of earth's crust. Doklady Akademii Nauk, 73(2): 295-297. |
Usui Y. 2015. 3-D inversion of magnetotelluric data using unstructured tetrahedral elements:applicability to data affected by topography. Geophysical Journal International, 202(2): 828-849. DOI:10.1093/gji/ggv186 |
Vallianatos F. 2002. A note on the topographic distortion of magnetotelluric impedance. Annals of Geophysics, 45(2): 313-320. |
Wannamaker P E, Stodt J A, Rijo L. 1986. Two-dimensional topographic responses in magnetotellurics modeled using finite elements. Geophysics, 51(11): 2131-2144. DOI:10.1190/1.1442065 |
Wei W B, Jin S, Ye G F, et al. 2003. Methods to study electrical conductivity of continental lithosphere. Earth Science Frontiers (in Chinese), 10(1): 15-23. |
Xu S Z, Ruan B Y, Zhou H, et al. 1997. Numerical modeling of 3-D terrain effect on MT field. Science in China Series D:Earth Sciences, 40(3): 269-275. DOI:10.1007/BF02877535 |
Yang D K, Oldenburg D W. 2016. Survey decomposition:A scalable framework for 3D controlled-source electromagnetic inversion. Geophysics, 81(2): E69-E87. DOI:10.1190/geo2015-0217.1 |
Yin Y T, Jin S, Wei W B, et al. 2017. Lithospheric rheological heterogeneity across an intraplate rift basin (Linfen Basin, North China) constrained from magnetotelluric data:Implications for seismicity and rift evolution. Tectonophysics, 717: 1-15. DOI:10.1016/j.tecto.2017.07.014 |
Yoshimura R, Ogawa Y, Yukutake Y, et al. 2018. Resistivity characterisation of Hakone volcano, Central Japan, by three-dimensional magnetotelluric inversion. Earth, Planets and Space, 70: 66. DOI:10.1186/s40623-018-0848-y |
Zhang K, Yan J Y, Lü Q T, et al. 2017. Three-dimensional nonlinear conjugate gradient parallel inversion with full information of marine magnetotellurics. Journal of Applied Geophysics, 139: 144-157. DOI:10.1016/j.jappgeo.2017.02.008 |
Zhang K, Lü Q T, Yan J Y. 2015. The lower crust conductor from Nanjing (Ning)-Wuhu (Wu) area in the Middle and Lower Reaches of Yangtze River:Preliminary results from 3D inversion of mag netotelluric data. Journal of Asian Earth Sciences, 101: 20-29. DOI:10.1016/j.jseaes.2015.01.007 |
Zhao G Z, Tang J, Zan Y, et al. 2005. Relation between electricity structure of the crust and deformation of crustal blocks on the northeastern margin of Qinghai-Tibet Plateau. Science in China Ser. D Earth Sciences, 48(10): 1613-1626. DOI:10.1360/02YD0047 |
曹晓月, 殷长春, 张博, 等. 2018. 面向目标自适应有限元法的带地形三维大地电磁各向异性正演模拟. 地球物理学报, 61(6): 2618-2628. DOI:10.6038/cjg2018L0068 |
陈辉, 尹敏, 殷长春, 等. 2018. 大地电磁三维正演聚集多重网格算法. 吉林大学学报(地球科学版), 48(1): 261-270. |
陈晓, 于鹏, 张罗磊, 等. 2011. 地震与大地电磁测深数据的自适应正则化同步联合反演. 地球物理学报, 54(10): 2673-2681. DOI:10.3969/j.issn.0001-5733.2011.10.024 |
邓居智, 陈辉, 殷长春, 等. 2015. 九瑞矿集区三维电性结构研究及找矿意义. 地球物理学报, 58(12): 4465-4477. DOI:10.6038/cjg20151211 |
董浩, 魏文博, 叶高峰, 等. 2014. 基于有限差分正演的带地形三维大地电磁反演方法. 地球物理学报, 57(3): 939-952. DOI:10.6038/cjg20140323 |
刘云鹤, 殷长春. 2013. 三维频率域航空电磁反演研究. 地球物理学报, 56(12): 4278-4287. DOI:10.6038/cjg20131230 |
秦策, 王绪本, 赵宁. 2017. 基于二次场方法的并行三维大地电磁正反演研究. 地球物理学报, 60(6): 2456-2468. DOI:10.6038/cjg20170633 |
阮帅. 2015.三维大地电磁有限内存拟牛顿反演. (博士).成都理工大学.
|
魏文博, 金胜, 叶高峰, 等. 2003. 大陆岩石圈导电性的研究方法. 地学前缘, 10(1): 15-23. DOI:10.3321/j.issn:1005-2321.2003.01.003 |
徐世浙, 阮百尧, 周辉, 等. 1997. 大地电磁场三维地形影响的数值模拟. 中国科学(D辑), 27(1): 15-20. |
赵国泽, 汤吉, 詹艳, 等. 2004. 青藏高原东北缘地壳电性结构和地块变形关系的研究. 中国科学D辑:地球科学, 34(10): 908-918. |