青藏高原是新生代以来地壳上地幔各圈层的能量和物质流动的集结之地，它的形成和演化是现今地球科学研究的焦点，查明青藏高原地壳上地幔构造是许多地球物理学家的研究目标.应用大地电磁法进行岩石圈的电阻率成像是揭示地壳上地幔三维构造的主要方法之一.自从1980年代以来, 在青藏高原已经开展了大地电磁法调查和二维剖面成像的研究，取得了很好的成果(Sun et al., 2003；谭捍东等，2004；魏文博等，2006；Zhao et al., 2008；金胜等, 2009, 2010；谢成良等，2012；王刚等，2017).但是, 有限的二维剖面成像取得的结果无法揭示青藏高原地壳上地幔的整体结构，难以对青藏高原的岩石圈地质作用给出全面的证据.从2010年起，本文的作者们在青藏高原局部多个地区进行了大地电磁法调查和三维成像研究，经过数据采集、资料处理和三维反演，取得了测点分布比较均匀的大地电磁法调查数据，为研究高原的三维电性结构和刻画高原岩石圈的内部构造，提供了新的可能性.本文报道通过大地电磁法探测和三维反演建立大地电磁法调查三维电阻率模型的初步成果.

青藏高原面积广大，自然条件恶劣，岩石圈电阻率三维结构的成像要求资料的逐渐积累，不可能一跺而就.笔者们采取分期分区采集数据的策略，先在SINOPROB-1项目中由中国地质大学(北京)组队，对西藏地区进行大地电磁法调查(杨文采等，2011；巩秀钢等，2012；金胜等，2012；Ouro-Djobo et al., 2018).“十二·五”项目结束后，本文第一作者又向中国地质调查局申请在青藏高原东部开展大地电磁法调查项目，依然由中国地质大学(北京)组队，用相同的仪器和方法于2017年完成了数据采集.此时中国地质大学(武汉)胡祥云研究团队也完成了在羌塘地区的三条大地电磁法剖面(Zeng et al., 2015).胡祥云教授大力支持本项三维成像研究，提供了羌塘地区的大地电磁法调查数据(图 1a中用Q1、Q2、Q3标明).加上少量其他调查取得的数据，我们的大地电磁法调查的区域覆盖了青藏高原的大部分面积，为三维电阻率成像研究打下了可靠的基础.在2014年“十二·五”项目结束后，中国地质大学(北京和武汉)都对有关局部地区的调查成果做了报道(见上述文献)，本文的报道是上述数据融合后进行全区反演成像的结果.

青藏高原大地电磁测点位置见图 1a，图 1b为测区位置与大地构造简图.全区数据采集观测物理点1065个，分布比较均匀.数据采集按《大地电磁测深法技术规程》等有关技术要求, 测点采用GPS定位(巩秀钢等，2012；金胜等，2012；Zeng et al., 2015).使用的仪器包括MTU-5宽频MT采集系统，和Lemi-417长周期MT采集系统.全部仪器经过标定和一致性实验，同时采集E_x、E_y、H_x、H_y和H_z五个大地电磁场分量数据，电道采用“十”字布极.大多数测点的观测时间在8000 s以上.每天测量不极化电极极差，同一对电极的极差小于1 mV，保持稳定在20 h以上.

图 1

Fig. 1

图 1 (a) 青藏高原大地电磁测点位置平面图，(b)测区位置与构造简图，图中的矩形框表示图(a)的测区范围 Fig. 1 (a) Distribution of magnetotelluric stations and the survey area in the Qinghai-Tibet Plateau; (b) Simplified tectonic map with the survey area shown in the rectangular block

大地电磁的资料融合、处理和反演成像由同济大学负责，他们2014年在塔里木盆地进行了数据融合、处理和反演成像，取得了相当好的成果(杨文采等, 2015d, e).运用的资料处理流程见图 2a(张罗磊等, 2009, 2010；陈晓等，2011).由于不同测区数据格式不完全统一，首先需要对数据进行整理和同化，改写所有三维反演输入数据的格式.同化的数据格式依次为频率、阻抗(Z_xx, Z_xy, Z_yx, Z_yy)的实部、虚部和误差.数据都同化为84个频点，频率范围为(320~0.172)×10^-3Hz.同化时也对数据进行挑选，取出950个误差小的数据用于正反演计算.

图 2

Fig. 2

图 2 (a) 资料处理流程图；(b)大地电磁反演流程框图 Fig. 2 Data processing flow-chart (a) and the inversion computation flow-chart (b)

资料处理主要包括统一数据格式、误差分析与压制、资料预处理、极化模式识别、静校正等.由于复杂的近地表条件和地下地质构造，造成野外观测的大地电磁曲线发生相当程度的畸变，必需对畸变数据进行校正处理，以保证大地电磁反演的质量.单纯使用局部算法直接求解GB分解定义式，初始值选取不当会导致求解过程不稳定.用基于混合优化算法的MT阻抗张量畸变分解方法，对GB张量分解畸变因子和区域阻抗的求解方法进行了改进(张罗磊等, 2009, 2010；陈晓等，2011).改进的单频点GB分解方法提高了求解过程的稳定性, 解决了由于实测数据包含噪音和实际构造无法符合GB分解模型假设所造成的求解结果不稳定的问题.试验表明，张量分解对消除近地表电性不均匀体造成的畸变影响效果显著，可以恢复真实的区域2-D阻抗张量，提供更全面的畸变参数信息，为下一步定量反演和解释打下了可靠的基础.

在完成观测数据的畸变改正后，进行反演二维电阻率模型，并用它作为三维反演的初始模型.传统的光滑模型约束的稳定泛函所得到的反演结果对电性边界刻画不够清晰，而一些考虑尖锐边界的稳定泛函则存在对计算参数选择的依赖性较大、稳定性不好等问题.在改进迭代优化算法后，形成了基于SVD(奇异值分解)迭代优化的二维反演新方法，并更合理考虑了数据误差的加权方式.

对全区九百多个大地电磁测点数据进行三维反演，对算法要求很高，运用的反演流程框图见图 2b.我们应用了交错网格有限差分法进行三维正演模拟计算，可非等间距剖分.由于剖分网格随深度逐渐增加，增加少数网格可减少正演时底部边界的影响.为提高反演空间分辨率，引入了最小梯度支撑泛函代替平滑模型约束，得到尖锐电性分界面，更有利于刻画地下构造.兼顾光滑和尖锐边界的“最小支撑梯度”新稳定泛函应用于正则化反演的目标函数，在使反演更稳定和收敛更快的同时，可以刻画尖锐的电性界面，使三维图像丰满清晰(Zhang et al., 2012；尹曜田等，2012；董浩等, 2012).三维反演采用的网格尺寸为20 km×20 km，垂直方向为26层，不等间距剖分到200 km.x与y方向扩边300 km.反演时主要使用了1~3000 s的数据，用23个周期(T)计算，测点位于网格中心位置.迭代进行了50次，迭代后拟合差减小到4%.由于本文只是青藏高原三维电阻率成像的初步成果，下一步我们准备继续积累数据，并对大于3000 s的数据进行精细处理和分析，以增强反演的深度和可靠性.

在进行青藏高原三维电阻率成像的研究之前，我们已经完成了青藏高原地壳三维密度成像和相关地壳构造的研究(杨文采等, 2015a, b, c；孙艳云，2015)，也在进行青藏高原地壳上地幔三维地震波速成像和相关地壳构造的研究(杨文采等, 2017, 2019a, b, c).笔者希望，通过青藏高原三维电阻率成像的研究，可以把不同方法的成像结果互相验证，取得对大陆碰撞造山作用过程的共同的认识.与密度和地震波速不同，和岩石圈的电阻率关系最密切的因素是流体，因此电阻率的平面和剖面图像与密度和地震波速图像应该有明显的差别，直观的对比和相关分析很难得到可靠的结论.总之，青藏高原三维电阻率图像的地质解释需要做大量的工作，不可能在一篇论文中完成.下面我们只能粗略地介绍一些取得的初步认识，和有兴趣的读者交流讨论.

青藏高原地壳不同深度三维反演电阻率平面图示于图 3，深度分别为2 km，20 km，40 km和60 km.先看2 km的表层电阻率平面图，电阻率的分区和地面的构造单元有很好的对应关系.低电阻率区分布在青藏高原东部、羌塘地块的中北部、喜马拉雅地块和西部阿里高原，与地表含水的沉积地层覆盖有关.拉萨地块和察隅地块电阻率都很高，说明它们的浅层结晶程度高，含水分很少.在深度20 km的上地壳中，低阻区还是分布在青藏高原东部、羌塘地块的中部和北部、喜马拉雅地块和西部阿里高原；但是，拉萨地块南部的电阻率变低了，因为这里也是新生代火山和岩浆活动比较活跃的地区.在深度40 km的中下地壳，低阻区依然分布在青藏高原东部、羌塘地块的中部和北部、喜马拉雅地块和西部阿里高原；但是，拉萨地块的电阻率除东部外全变低了，可能是因为下地壳流活动造成岩石变形引起的(杨文采等, 2017, 2019a).在60 km的下地壳，电阻率扰动的平面模式大体上和深度40 km平面相同，但是，喜马拉雅地块的电阻率大部分变高了.喜马拉雅地块的电阻率变高，可能反映了高电阻率的印度克拉通板块的俯冲的影响.

图 3

Fig. 3

图 3 青藏高原地壳不同深度三维反演电阻率平面图 (a)深度2 km; (b)深度20 km; (c)深度40 km; (d)深度60 km. Fig. 3 The plane resistivity images obtained by the 3-D inversion at depth 2 km (a), 20 km (b), 40 km (c), and (d) 60 km respectively, corresponding to the crust of Qinghai-Tibet Plateau

现在来看地幔电阻率扰动的平面展示，图 4b—d分别为深度100 km、140 km和200 km的三维反演电阻率平面图.和深度60 km平面上的电阻率扰动(图 4a)比较, 深度100 km地幔的电阻率扰动图中, 喜马拉雅完全变为高电阻率地块，代表的是印度克拉通板块的俯冲前沿，俯冲前沿到达了拉萨地块东部和南部，这里的电阻率也变高了.同时青藏高原东部的低阻区面积沿100°E有所扩大，说明这里岩石圈的变形在加剧.深度140 km已经到岩石圈地幔的底部，察隅地块、喜马拉雅地块和拉萨地块东部联成统一的高电阻率地块，它们反映了向北东俯冲的印度克拉通地幔.青藏高原东部的低阻区沿100°E继续扩大，反映在金沙江断裂带有切穿岩石圈的趋势.深度200 km已经到软流圈地幔了，可以看到青藏高原分布了大面积的低电阻率异常，分布在羌塘盆地北部、拉萨地块中东部，阿里—帕米尔和100°E两侧，它们很可能反映了软流圈流体的上涌.从平面图上看到的这些构造特征，和以前对密度和地震波速成像结果是一致的.

图 4

Fig. 4

图 4 青藏高原地幔不同深度三维反演电阻率平面图 (a)深度60 km; (b)深度100 km; (c)深度140 km; (d)深度200 km. Fig. 4 The plane resistivity images obtained by the 3-D inversion at depth 60 km (a), 100 km (b), 140 km (c), and (d) 200 km respectively, corresponding to the mantle of Qinghai-Tibet Plateau

地幔电阻率平面扰动的模式明确显示，青藏高原东、西部的地体碰撞和拼合形式是不同的.在青藏高原西部，羌塘、拉萨和喜马拉雅等地体从北到南碰撞拼合.在青藏高原东部，羌塘—拉萨、察隅、印支、雅安和扬子等地体拼合的方向变化大，有一部分从西到东拆离拼贴.西部的北南向拼合发生在先，已经在地壳形成了广泛的变形构造，包括正交的拆离带和裂谷系.东部的拆离拼合发生在后，正在地壳造成广泛的变形，包括不正交的拆离带和压扭构造系.从高阻-低阻区的分布看，东部的地体拼合根源应该在地幔，今后还会进一步发展.

大地电磁法反演得到的三维电阻率图像上提供了许多岩石圈局部构造信息.为了快速考察青藏高原岩石圈的局部构造，要对比不同方向的电阻率剖面图.图 5为青藏高原不同纬度的东西向电阻率剖面图.由图 5显示，青藏高原从西向东有4条低电阻率带，它们分布在阿里、定日—色林错、罗布莎—唐古拉山口和金沙江.从位置上看，西部的3条和西藏的南北向裂谷对应，表明裂谷带在深处可能有地幔热流体侵入(杨文采等, 2019a, b, c).金沙江低电阻率带的位置和印支地体边界对应，反映了2亿年前的古板块拼合的边界，可能至今没有完全封闭，依然有热流体活动，造成现今西藏地面较高的大地热流(汪集旸, 2015).从图 5看来，在33°N的羌塘电阻率扰动的模式发生了突然的变化，和南部的4条岩石圈剖面明显不同，电阻率异常带不和南部相连，可能表明羌塘地体现今受印度板块的陆-陆俯冲影响很小，不同于南部的喜马拉雅、阿里和拉萨地体.

图 5

Fig. 5

图 5 青藏高原不同纬度的东西向电阻率剖面图，从上到下纬度分别为33°N, 32°N, 31°N，30°N和29°N Fig. 5 The resistivity profiles obtained by the 3-D inversion along latitude of 33°N, 32°N, 31°N, 30°N and 29°N respectively in the Qinghai-Tibet Plateau

在南部的4条剖面中，经度80°E—81°E、87°E—88°E、92°E—93°E和99°E—101°E都有明显的低电阻率异常，其中92°E—93°E的低电阻率异常幅度比较弱.这些异常在软流圈都有低电阻率的根，向上面有的发生分岔.联系到青藏高原的大地热流值很高(汪集旸, 2015)，它们可能反映了来自软流圈的热流体上涌.从低电阻率异常的形状看，热流体上涌到下地壳的深度后发生了水平扩展，然后沿上中地壳的裂隙带向上升腾，最后可以在地表形成高温热泉.不过，如羊八井那样的高温热泉，虽然在上地壳有低电阻率异常显示，但是位置并不在软流圈的热流体上涌的正上方，而是位移了百十千米.

再来对比南北方向的电阻率剖面图.图 6为青藏高原西部地区的3条南北向电阻率剖面图.由图 6可见，在西部阿里地区中下地壳的电阻率都比较低，而上地壳和岩石圈地幔电阻率都高.图 7为青藏高原中部地区的5条南北向电阻率剖面图.由图可见，在喜马拉雅和拉萨地体之间的雅鲁藏布江缝合带下方岩石圈的电阻率比较低，低电阻率层倾向于从喜马拉雅地体向拉萨地体下方俯冲，可能是印度大陆板块俯冲顶面发生摩擦生热，造成了岩石圈部分熔融.羌塘地体下方的岩石圈也有宽大的低电阻率带发育，地壳的低电阻率异常分布比较发散，不集中.结合青藏高原三维地震波速成像结果可以认为，羌塘地体下方的低电阻率异常带是由于特提斯洋在软流圈拆沉，引起的浮力激发热流体上涌和火山爆发引起的(杨文采等, 2019a, b, c).

图 6

Fig. 6

图 6 青藏高原西部南北向电阻率剖面图, 对应青藏高原阿里地区地壳和最上地幔 Fig. 6 The resistivity profiles obtained by the 3-D inversion along longitude of 80°E—83°E, corresponding to the crust and uppermost mantle of the Ali area in the Qinghai-Tibet Plateau

图 7

Fig. 7

图 7 青藏高原西部南北向电阻率剖面图，剖面W6—W13位置见图 1a Fig. 7 The resistivity profiles obtained by the 3-D inversion along longitude of 83°E—92°E, corresponding to the crust and uppermost mantle of the west Qinghai-Tibet Plateau

转向青藏高原东部，在图 8显示了8条东部的南北向电阻率剖面图.从剖面E2—E6可见，94°E到东边的金沙江断裂带，岩石圈的电阻率比较高.不过，羌塘地块的上地壳和杂多以北的岩石圈地幔电阻率还是比较低的，和西边的低电阻率异常带是相连的.松潘—甘孜复理石盆地下方的上地壳和岩石圈地幔电阻率比较低；巴塘和雅江地块下方的中下地壳和岩石圈地幔电阻率比较高，而且岩石圈也变薄了.结合三维地震波速和密度扰动成像结果分析可以认为(杨文采等, 2015a, b, c，2019a, b, c)，位于青藏高原东构造结的察隅地块岩石圈是高密度、高波速和高电阻率的，它对青藏高原东部的楔入，使其北部和东部地块的岩石圈发生拆离撕裂，也可能造成热流体上涌的低电阻率异常带.

图 8

Fig. 8

图 8 青藏高原东部南北向电阻率剖面图，剖面E2—E12位置见图 1a Fig. 8 The resistivity profiles obtained by the 3-D inversion along longitude of 93°E—103°E, corresponding to the crust and uppermost mantle of the east Qinghai-Tibet Plateau

青藏高原的大地电磁法调查区域覆盖了高原大部分面积，为全区三维电阻率成像研究打下了可靠的基础.对多个测区大地电磁数据进行精细的同化处理和反演成像，取得了青藏高原可靠的岩石圈三维电阻率结构图像.成像的区域为28°N—35°N，80°E—104°E.三维反演计算时采用的网格尺寸为20 km×20 km，垂直方向不等间距剖分为26层.结果表明，青藏高原现今岩石圈电阻率扰动主要反映印度克拉通对亚欧大陆板块俯冲引起的热流体运动和大陆碰撞和拆离产生的构造.在岩石圈地幔，察隅地块、喜马拉雅地块和拉萨地块东部联成统一的高电阻率地块，它们反映了向北东俯冲的印度克拉通.雅鲁藏布江、班公—怒江和金沙江缝合带都有明显的低电阻率异常，表明岩石圈深处有热流体活动.雅鲁藏布江、班公—怒江和金沙江缝合带都有明显的低电阻率异常，也表明它们的岩石圈还有流体活动.青藏高原东部的低阻区沿100°E向地幔下方扩大，反映了金沙江断裂带有切穿岩石圈的趋势.地幔电阻率平面扰动的模式显示，青藏高原东西部的地体碰撞拼合形式和方向是不同的.在青藏高原西部，羌塘、拉萨和喜马拉雅等地体从北到南碰撞拼合.在青藏高原东部，羌塘—拉萨、察隅、印支、雅安和扬子等地体多方向拆离拼合，在地壳造成不正交的拆离带和压扭构造系.从高阻-低阻区的分布看，东部的地体拼合有地幔的根源，今后还会进一步发展.察隅地块岩石圈对青藏高原东部的楔入，使其北部和东部地块的岩石圈发生拆离撕裂，也造成热流体上涌的低电阻率异常.