大地电磁测深法(MT)是通过天然交变电磁场研究地球模型电性结构的一种地球物理方法.大地电磁测深方法具有较大的穿透深度和分辨能力，它不仅是研究地壳与上地幔深部构造的强有力工具(李冉等，2014；吕庆田等，2011)，而且可以结合地震数据(刘彦等，2012；詹艳等，2014)和重磁数据(张交东等，2012)进行综合解释分析，已成为矿产勘查(Hu et al.，2013；陆桂福和吴新刚，2014)、地热资源(Wu et al.，2012；高景宏等，2010)和油气勘探(夏训银等，2012)的有效方法.

地形对大地电磁测深的影响是地球物理工作者难以避免的问题.特别是我国山地区域地形起伏剧烈，大地电磁数据采集受各种因素的限制，当测点布置在山沟时，垂直山体方向受高电阻体对电场的排斥作用的影响；若山体不规则，地形影响更为复杂(张翔等，1999).有限元方法是模拟MT正演响应的有效手段，许多研究者利用矩形网格(Rodi，1976；陈乐寿，1981；童孝忠等，2009；张昆等，2008)和矩形对分三角网格(Wannamaker et al.，1987；Zyserman et al.，1999；刘小军等，2007；柳建新等，2009)模拟地形正演响应，结果表明地形对TM模式的影响远远大于TE模式(王绪本等，1999)，并且视电阻率曲线表现为畸变、静态平移或者二者的混合，没有一定规律可循.但矩形及其对分三角网格对复杂地质结构(起伏地形和交错地层等)模拟能力有限，且上述规则网格在对目标区域进行加密时，其背景区域也会受到不必要细化，增加了计算时间.

近年来许多学者(Key and Weiss，2006；Franke et al.，2007；Li and Key，2007；Li and Pek，2008)采用自适应非结构三角网格有限元方法计算电磁场.非结构三角网格能够真实地模拟各种复杂地质结 构，该方法利用Delaunay三角划分原理(Shewchuk，2002)对目标区域进行细密网格剖分，而背景区域则粗糙网格剖分.自适应有限元对计算结果进行误差评估(Ovall，2006)，误差较大的单元网格自动细分直至有限元计算结果满足给定的误差范围.自适应非结构有限元方法因其强大的适应性，越来越受到地球物理研究者的青睐，其方法技术也日臻完善.现已成功实现了海洋二维MT正演模拟(Key and Weiss，2006；Franke et al.，2007；赵慧等，2014)，海洋2.5D可控源正演模拟(Key and Ovall，2011).在此基础上，本文采用非结构三角网格剖分、自适应有限元与Occam算法结合的策略来研究陆地MT二维带地形反演.首先讨论了不同地形对MT正演的影响；其次建立3个典型地电模型进行带地形反演，并与平地反演结果进行对比；最后对实测MT数据进行反演并解释.

假设地下电性结构是二维的，x沿构造走向方向，电导率σ只在(y，z)平面变化，对于平面波源，可以将电磁场(E_x，H_x)分解成一次场(E⁰_x，H⁰_x)和二次场(E′ _x，H′ _x)(Coggon，1971)：

事先给出一次场对应的电导率为σ₀，一次场为均匀半空间模型其电磁场值由解析解算出，二次场则是由σ-σ₀的差异引起，采用自适应有限元解出场值，最后将一次场和二次场相加算出总场值.该方法可通过提高二次场解的精度来提高整体解的精度(Rodi，1976)，与总场法求解的结果对比验证了该方法的正确性(赵广茂等，2008；刘小军等，2007).这种方法最大优势在于有限元网格只需要精确地模拟二次场的变化，因此可采用包含更少单元的网格(Key and Ovall，2011)，提高计算速度.利用Helmholtz 方程得到二次场TE和TM偏微分方程表达式如下：

式中ω为角频率，μ=4π×10^-7 H·m^-1为磁导率常数.对研究区域进行三角剖分后生成包含N个顶点的N维向量空间V，式(3)和(4)对应的变分问题可以统一写成：B(v，u)=F(v)，(5)其中 v 、 u ∈ V，u 是场值，v 是向量空间任意值，即对于任意 v，场值 u 都满足式(5).利用有限元方法解上述变分问题即可得到E′ _x与H′ _x.根据偏导关系，求出二次场其他分量的值

二次场与一次场相加求得总场值，则MT阻抗为

MT视电阻率和相位角为

相对于有限元的解u_n，解的误差G(u-u_n)更被研究者重视，其中u是(5)式的真实解，G为误差评估的对偶函数.自适应网格划分的实质是计算出每个三角单元有限元解的误差，对误差较大的三角单元再次细分，直到误差达到允许范围.本文采用对偶加权残差(DWR)法(Ovall，2006)计算单元误差，属于一种后验误差处理技术，通过分级函数对残差进行加权.选择一个合适G函数十分重要，DWR给出了一个比较合适的G函数，单元误差计算过程中多次利用有限元数值方法进行求解，所求误差是逼近真实解的近似值.后验误差方法在工程学领域是一项非常成熟的技术，近几十年被引入到电磁场模拟中，已应用于海洋可控源电磁场 (Li and Key，2007； Key and Ovall，2011)、海洋二维MT电磁场(Franke et al，2007； Key and Weiss，2006；赵慧等，2014)和二维电阻率各向异性介质大地电磁场(Li and Pet，2008)模拟.

反演问题采用Occam算法(Constable et al.，1987)，该算法自发表以来由于其稳定性在电磁法反演领域得到了广泛的应用(Ghaedrahmati et al.，2014；Siripunvaraporn and Sarakorn，2011；Sudha et al.，2014；何梅兴等，2011).Occam反演实质为带光滑约束的最小二乘法，其泛函U表达式为：U=‖ m ‖²+λ { ‖ Wd - W F〖 m 〗‖²-X^2* }，(10)(10)式右边第一项和第二项分别为模型粗糙度因子和数据拟合误差.其中 m 为模型参数表示电阻率或电导率，d 为实测数据， 为粗糙度系数矩阵，W 为数据误差矩阵，F〖 m 〗为正演函数，X^2*为期望的数据拟合误差，λ为拉格朗日乘子，用以平衡粗糙度和数据拟合误差，当λ取较小值时，反演以搜索光滑模型为主，反之则以搜索最小拟合误差为主.反演过程中自动改变λ值，最后得到拟合误差最小情况下的最光滑模型.

自适应有限元方法进行电磁场数值模拟的主要优势在于可根据每次有限元计算的结果自动进行网格精炼，确保了计算结果的精度，避免了以往正反演过程中网格绘制的繁琐工作，图 1给出了自适应有限元与Occam算法结合的计算流程.首先给出反演初始模型作为反演细化的初始网格和电阻率模型参数，一般来说，目标区域网格划分细致，选择小三角单元；背景区域则采用大三角单元.随后进行内存空间分配，即采用多少线程参与计算，每个线程负责多少测点的自适应有限元运算.本文设置单元误差阈值0.02，采用有限元方法计算每个单元的场值，由DWR方法估计单元局部误差，将局部误差大于误差阈值的单元进行细化，得到新的三角网格.重复以上过程直至所有单元误差达到误差阈值或者已达到最大网格细化次数，则网格精炼停止，输出正演结果.设置反演收敛精度1.5，计算数据均方根误差RMS，若未达到反演收敛精度，则更新模型重新开始自适应有限元计算过程直至达到反演收敛精度或者达到最大迭代次数，输出反演结果.

图 1

Fig. 1

图 1 自适应有限元算法反演流程 Fig. 1 Inversion process of adaptive finite element algorithm

为了研究地形对MT勘探的影响，在100 Ωm的均匀半空间中设置一10 Ωm的矩形低阻体，在图 2a中设置包含山谷、山峰起伏地形，最大高差300 m，图 2b采用水平地形.从左往右依次设置11个测点，间距300 m，频率范围为0.1~10000 Hz，共22个频点.

图 2

Fig. 2

图 2 地形影响研究模型 (a)起伏地形模型；(b)平地模型. Fig. 2 Model for research of terrain effect (a)Model with terrain;(b)Flat model.

采用自适应有限元方法计算各测点的视电阻率和相位响应值，选取代表典型地形的三个测点 MT04、MT05和MT06正演结果进行比较.如图 3所示，红色曲线为TM模式，蓝色曲线为TE模式.图 3a中TE模式视电阻率和相位曲线相对于平地模型(图 3b)波动不大；但是对于TM模式，在山峰处(MT04)视电阻率曲线低频陡然下降，山坡处(MT05)视电阻率曲线整体向上平移，且在高频处上扬，山谷处(MT06)视电阻率曲线高频处上扬.不同地形对正演响应产生的干扰不同，视电阻率曲线会产生畸变、平移或二者的结合，没有十分确切的规律可循，整体来说地形对视电阻率的影响大于相位角，对TM模式的影响远远大于TE模式，推测其原因是两种极化模式中，磁场一般变化平缓，TE模式电场平行于构造走向，TM模式电场垂直走向，垂直于走向的电场会产生感应电场，该感应电场在空间上的分布受地形起伏的影响严重，而在TE模式中只有与走向平行的电场.因此，实现带地形反演消除因地形产生的数据异常十分重要.

图 3

Fig. 3

图 3 MT04，MT05，MT06测点视电阻率和相位响应 (a)图 2a模型响应；(b)图 2b模型响应. Fig. 3 Apparent resistivity and phase response at MT04，MT05，MT06 site (a)Response for Fig. 2a model;(b)Response for Fig. 2b model.

在山区进行MT工作主要以寻找金属矿为目 的，金属矿一般呈现低阻性质.如图 4所示，在100 Ωm均匀半空间埋藏两个形态不一，埋深不同的10 Ωm的低阻体，1号低阻体中心埋深400 m，呈块状；2号低阻体中心埋深750 m，呈薄板状.测线垂直走向方向分布，全长2000 m，其中包含40个测点，频率范围0.1~10000 Hz，共60个频点，地形采用一组实测高程数据，最大高差达250 m.

图 4

Fig. 4

图 4 金属矿模型(其中1、2为低阻体编号) Fig. 4 Metal mine model(1，2 indicates number of low resistivity body)

以图 4模型进行网格剖分，初始网格划分为2035个三角单元(图 5a)，对单元误差取对数值显示，其分布如图 5b所示.采用四线程对初始网格细化，以其中一个线程为例，该线程负责最后10个测点的自适应正演运算，经过8次细化，生成18854个三角单元(图 5c)，单元误差达到误差阈值以内(图 5d)则停止计算.网格主要在电性分界处即低阻体边界和测点附近进行细化，图 5e显示测点附近网格细化信息.其余三个线程同时进行自适应细化过程，共耗时123 s.

图 5

Fig. 5

图 5 网格细化 Fig. 5 Procedure of mesh refinement

利用上述模型的正演数据添加4%的随机噪声进行Occam反演，对于目标区域采用小三角单元细致剖分，周围背景区域则粗糙剖分.反演初始模型如图 6，其中包含5958个电阻率为100 Ωm三角单元.

图 6

Fig. 6

图 6 反演初始模型 Fig. 6 Triangulation for inversion

带地形反演结果(图 7a)显示1号和2号低阻体中心位置，可圈定出1号低阻体的轮廓，2号低阻体边界较模糊.图 7b为相同数据平地反演结果，相对于带地形反演，低阻体出现中心位置偏移及范围扩大的情况，且在近地表处存在虚假异常也称为静态效应.地形是引起静态效应的原因之一，常规处理中，反演前会进行静态效应校正，用以消除浅层假异常现象，常用方法为整体平移数据曲线，这需要丰富的经验且平移结果不一定准确，带地形反演不能完全消除但可以有效地抑制静态效应.

图 7

Fig. 7

图 7 金属矿模型反演结果(其中1、2是低阻体编号) (a)带地形反演；(b)平地反演. Fig. 7 Inversion result of metal mine model (1，2 indicates number of low resistivity body)(a)Inversion with terrain；(b)Inversion without terrain.

为了验证该算法对浅部薄层低阻矿床的反演效果，我们建立图 8a所示的地电模型，模型分为两层，上层电阻率为100 Ωm，层厚约870 m；下层为1000 Ωm，上层存在10 Ωm的低阻体，中心埋深450 m.测点、地形、频点设置与图 4模型一致.反演结果(图 8b)较好地显示薄层低阻体的中心位置，可识别下层高阻带，且地层分界明显.

图 8

Fig. 8

图 8 薄层低阻体模型 (a)真实模型；(b)反演结果. Fig. 8 Thin layer model with low resistivity(a)True model;(b)Inversion result.

在我国地形起伏的山区中，构造发育复杂，常见褶皱、断层分布.建立相应的综合褶皱模型(图 9a)进行反演，其中以向斜为主，连接2个小型背斜，10和100 Ωm地层互层分布.测点、地形、频点设置与图 4模型一致.反演结果(图 9b)显示该褶皱的大致轮廓，地层分层清晰，可看出中部向斜构造.

图 9

Fig. 9

图 9 综合褶皱模型 (a)真实模型；(b)反演结果. Fig. 9 Complex fold model (a)True model;(b)Inversion result.

表 1给出了上述三个模型反演耗时，RMS误差及粗糙度.由于本文模型运算采用四线程，模型每更新一次，便重新开始自适应细化，导致反演时间大大增加.在实际应用中，处理大量的数据需要更多的时间，如何优化算法是亟待解决的问题.

表 1 (Table 1)

表 1 反演耗时、误差及粗糙度列表 Table 1 List of inversion runtime，RMS and roughness with four processors 模型 耗时(s) 误差(RMS) 粗糙度 金属矿体 11370 1.48 1.337 薄层低阻体 15100 1.84 10.66 综合褶皱 18520 13.93 316.7

表 1 反演耗时、误差及粗糙度列表 Table 1 List of inversion runtime，RMS and roughness with four processors

工区位于中国西南山地矿区，区内地形起伏剧烈.该区地表水系发育，矿区出露地层比较单一，主要为中三叠统许满组及第四系，但构造发育复杂，主 要为褶皱、断裂和节理.数据采集使用加拿大phoenix 公司V8多功能电法仪，通过不极化电极测量E_x、E_y两道正交的电场水平分量数据，磁传感器测量H_x、H_y、H_z三道正交的磁场分量数据，从而实现五个电磁场分量的电磁测深数据采集.采集频率范围设置为10400~0.35 Hz，在低频处截屏到1 Hz，探测深度约地下2 km以内电性结构.数据噪声干扰不强，在高频有一天然死频带，采用D+方法进行处理.受地形影响，TM和TE曲线首端较多分离，拟断面图“挂面条”现象严重.

选择工区中A3测线中部进行反演.测线由西南向东北方向布设，最大高差242 m，地形起伏明显.工区布置如图 10所示(中间测线为A3)，图左下角红色点表示反演起始测点，右上角红色点为终点. A3测线所在位置出露地层为中三叠统许满组第三段(T₂xm³)，其包含的主要岩矿石电性见表 2.图 10显示A3测线中部穿越一大型向斜构造和一背斜构造，测线末端发育一系列小型褶皱，与两大褶皱主轴近似平行，测线中部和末端分别穿越小型河流.

图 10

Fig. 10

图 10 工区地形地质及工程布置图 Fig. 10 Geological and arrangement map in work area

表 2 (Table 2)

表 2 测区主要岩矿石电阻率统计表 Table 2 Resistivity of mineral and rock in work area 地层 岩矿石 名称 平均厚度 (m) 测定 块数 平均电阻率 (Ωm) 变化范围 (Ωm) T2xm3 粘土岩 400 25 95.38 9.48~300.3 砂岩 500 25 592.7 88~2401 泥灰岩 800 25 901.1 110.3~2688

表 2 测区主要岩矿石电阻率统计表 Table 2 Resistivity of mineral and rock in work area

A3测线中部共40个测点，测点间距约40 m，全长1600 m，频率范围1~10400 Hz，共53个频点，反演共耗时31230 s，RMS为3.74.结果如图 11所示，大致可划分为三层电性层，黑色虚线表示地层线，浅地表处为低阻层，随深度增加电阻率逐渐上升，变化范围10~1000 Ωm，与表 2所列T₂xm³地 层电阻率变化情况一致；反演结果显示背斜和向斜构造，背斜上部被剥蚀，整体电性表现为高阻，在其内部包含低阻体.工区地质资料显示该处断层、裂隙、水系发育丰富，推测该低阻体成因可能是地表水沿裂隙渗漏使该处呈现低阻性质.

图 11

Fig. 11

图 11 A3测线反演结果 Fig. 11 Inversion result of line A3

另一组数据为寻找锰矿实测电磁数据，锰矿成矿地质体主要由炭质、粉砂质粘土岩夹条带状、块状菱锰矿组成，其综合电性特征表现为低电阻率特征，其上覆地层为含砾砂岩、细晶白云岩等高电阻率岩性.

测线共25个测点，全长1200 m，0.35~10400 Hz 共60个频点.其反演结果如图 12所示，地层分层明显，浅部地层为高阻，厚度约200~300 m，随深度增加电阻率逐渐减小，在200 m深度处电阻率有增大趋势，之后又逐渐减小.该测线处分布有三口钻孔，如图中黑色圆柱所示，分别位于测线前段(ZK01)、 中段(ZK02)和后段(ZK03).ZK01见矿深度214 m； ZK02见矿深度261 m；ZK03见矿深度501 m.反演结果中钻孔处视电阻率随深度变化曲线如图 13，括号内的数字代表钻孔所在水平位置，箭头指向见矿深度.三条测深曲线变化趋势一致，电阻率随深度迅速减小后又缓慢回升，见矿深度并不在视电阻率最小处，而是在视电阻率缓慢回升带，结合反演结果和钻孔资料用红色虚线在图 12中圈定出成矿构造带大致范围.

图 12

Fig. 12

图 12 锰矿实测MT数据反演结果 Fig. 12 Inversion result of real MT data from manganese ore

图 13

Fig. 13

图 13 电阻率测深曲线 Fig. 13 Resistivity sounding curves

本文利用非结构自适应有限元算法研究地形对MT正演数据的影响，该方法能够真实地模拟起伏地形、不规则异常体及复杂形态的地质构造.研究结果表明地形对TE模式影响很小，对TM模式影响远大于TE模式，与平地视电阻率正演曲线相比，TM模式在山峰处低频部分下降，斜坡处曲线整体向上平移，在山谷处高频部分上扬.

将自适应有限元与Occam算法结合进行带地形反演运算，通过理论模型试算和分析，结果证明了该算法能清晰识别出异常体的中心位置且分层性良好，与平地反演结果相比，不会在浅部产生虚假异常，能有效抑制静态效应.应用于电磁实测数据反演结果显示了地层分层信息、背斜构造、见矿层位，探测的地质结构和圈定的成矿构造带均与已有地质资料相符合，表明该算法具有较强的实用性，可用于实测资料解释分析.

在此基础上，将继续研究带地形反演条件下如何正确进行静态校正，优化算法缩短反演计算时间，并根据其他已知信息进行模型约束反演和联合反演，最终将非结构网格引入起伏地形下大地电磁三维正反演研究中.

致谢 感谢中国地质科学院吕庆田研究员和中南大学汤井田教授对文章的悉心指导，对手稿提出的意见丰富了文章的内容.