龙门山断裂带南段地壳电性特征：来自速度结构约束下大地电磁反演的证据

许乐红, 谭捍东*, 吴萍萍, 彭淼，王帅军，姜枚

1 中国地质大学(北京)地球物理与信息技术学院，北京 100083 2 中国地震局地球物理研究所，北京 100081 3 中国地震局地球物理勘探中心，郑州 450002 4 中国地质科学院矿产资源所，北京 100037

0 引言

龙门山断裂带形成于印支期松潘—甘孜地块与扬子地体的碰撞以及逆冲推覆作用(林茂炳，1992)，三条主断裂呈近南北向分布，其南段位于华南地块、松潘—甘孜地块和川滇地块交界处，紧邻东西向的鲜水河断裂(图1b).区域内地形复杂、演化历史及板块动力学机制等问题仍存在很大争议，因此该区域成为研究青藏高原物质向东运动及青藏高原东缘隆升机制的重要场所(许志琴等, 1991；Worley and Wilson,1996; Hubbard et al.,2009; Xu et al.,2009).目前青藏高原剧烈隆升的主要机制有：中下地壳管道流模型(Royden et al.,1997,2008; Clark and Royden,2000)、侧向逃逸模式(Tapponnier et al.,2001)与上地壳缩短模型(Hubbard et al.,2009).近年来随着地球物理方法的不断改进和多资料综合分析研究的深入，对传统经典模型的质疑也逐渐增加.Clark和Royden(2000)认为青藏高原东缘的地壳增厚无法用断层与褶皱作用缩短上地壳的理论来解释，大部分抬升区上地壳没有沿着大型逆冲断层方向的显著缩短(Burchfiel et al.，1995)，高原东缘的地势升高可能是深部地壳增厚与中央高原下方下地壳疏散的直接结果.朱介寿等(2017)认为由青藏高原腹地羌塘高原流出的黏滞性地壳流，带动上地壳块体水平移动，遇四川盆地阻挡后的分层形成了上地壳的隆升与下地壳的增厚.He等(2014)认为印度板块持续向北和印度—亚洲碰撞引起的东西向的相对压缩是龙门山逆冲断裂带和松潘—甘孜地块地壳增厚的主导因素.

大地电磁法(Magnetotelluric sounding，MT)，探测深度大、不受高阻层屏蔽，广泛应用于壳幔结构探测领域.横跨龙门山断裂带已有多条大地电磁剖面的成果.Zhao等(2012)利用MT数据与其他地质和地球物理观测相结合，提出松潘—甘孜地块隆起是由四川盆地和龙门山下方稳定岩石圈阻挡的东南向地壳流引起的.张乐天等(2012)分析了青藏高原、龙门山断裂带及四川盆地下方电性结构的差异及分布特征.詹艳等(2013)在龙门山中段、南段布设两条大地电磁剖面，分析了中段与南段的电性差异.王绪本等(2018)结合长周期大地电磁数据分析了龙门山南段电性结构与中、北段的区别，推断芦山地震形成于深部隐伏壳幔韧性剪切带向上扩展的影响.不同电阻率剖面均显示松潘—甘孜地块下方存在高导层，龙门山断裂带下方扬子地体呈现整块高阻特征，但大地电磁常规带地形反演，垂向分辨率不高，受体积效应影响边界不够清晰，对深部构造分析可能造成误差.

近年来大地电磁数据与地震数据之间的约束、联合反演算法已日趋成熟并在实际资料中验证了其有效性与实用性.吴萍萍等(2021)利用VP/VS波速比地震学模型约束二维大地电磁反演，并应用于张渤地震带的大地电磁资料处理中，分析了区域内地质构造及孕震环境.陈高等(2010)利用基于不同物性参数随机分布网格模型的大地电磁与地震数据同步联合反演，对黔中隆起地区古生界分布及构造特征进行了分析.同时，多种地球物理资料综合分析在龙门山断裂带构造解释中也有进展.李大虎等(2015)利用台网及流动台站天然地震数据，P波到时方法结合重力、航磁数据反演，得到了鲜水河断裂区三维P波速度结构，认为松潘—甘孜地块地壳物质较软弱，盆地西缘密度较高、物质坚硬.徐啸等(2016)利用宽角反射、折射地震数据与重力数据联合反演，结合地壳结构与密度分布，分析了龙门山及邻近区域上中下地壳厚度、速度及密度的分布规律，推断四川盆地刚性基底西缘因挤压产生的弯曲应力是该地区抬升的重要条件.目前，大地电磁方法与速度结构综合分析解释还未在龙门山南段开展.对于地球深部结构及其地球动力学研究而言，多资料联合解释、互相验证很可能是未来的发展趋势.

本文实现了一种速度模型约束下的带地形大地电磁二维反演算法，并应用于龙门山南段大地电磁实测数据中.基于交叉梯度结构约束方法，利用人工地震高分辨率宽角反射/折射P波的速度结构约束大地电磁的电阻率模型，获得了相较于常规反演方法更可靠的电性结果，为分析龙门山断裂带南段的深部结构提供了更准确的约束.

1 区域地质构造和大地电磁数据

1.1 区域地质构造

龙门山断裂带沿走向大致以卧龙—怀远、北川—安县一线为界将其划分为南、中、北三段(李智武等, 2008; 刘树根等, 2009).此次研究区域为龙门山南段，区域包括整个造山带及其前陆盆地和松潘—甘孜地块，横跨区域内主要构造单元，且相对于中北段，其形变历史相对较晚，现今活动性强(李智武等，2008).区域包含了三个主构造单元，松潘—甘孜褶皱带、龙门山断裂带以及四川盆地，是龙门山南段构造发育的代表性地区.

龙门山三大主断裂包括：北川—映秀断裂(YBF)、安县—灌县断裂(GAF)和茂县—汶川断裂(WMF).北川—映秀断裂在龙门山南段为耿达—陇东断裂(LDF)；安县—灌县断裂在南段为盐井—五龙断裂(WLF)；茂县—汶川断裂南段是双石—大川断裂(SSF).其中北川—映秀断裂形成于三叠世早期松潘—甘孜造山带与扬子地体碰撞期，三叠世晚期松潘—甘孜地块与上扬子地体发生陆内俯冲，三叠世末期松潘—甘孜地块左旋走滑逆冲于扬子地体之上形成了松潘—甘孜山与前陆盆地(乔秀夫等, 2012).安县—灌县断裂与茂县—汶川断裂的逆冲活动发生于侏罗纪-早白垩世，古龙门山形成一个由三条逆冲断裂带形成的山体.现今的龙门山经历了中生代地震造山与平静期山脉剥蚀降低的过程，以及新生代多次古地震瞬时造山活动(乔秀夫等, 2016).龙门山南段的前缘及其纵深地区——宝兴地区的宝兴杂岩，存在着多次推覆和多层滑覆现象，以及推覆滑覆叠加样式(林茂炳，1992).印支期和喜马拉雅期的造山运动，使得这一地区地震频发、构造极不稳定.先后发生了2013年4月20日芦山地震以及2014年11月22日康定地震.芦山地震震源精定位的震中为30.28°N，103.98°E，震源深度16.38 km(王小娜等，2015).初步调查认为，主震源位于双石—大川断裂附近，多数余震分布在主断裂西南侧区域(刘杰等，2013).龙门山断裂南段的速度结构、地壳电性结构分布、断裂带产状和深部延伸特征均与大地震活动密切相关.

1.2 大地电磁数据来源

雅安地震后由汶川科钻项目的“横穿龙门山地震反射剖面探测及深部结构的地震解释”课题支持，于2013年至2015年间分三次穿越龙门山断裂带进行大地电磁测深野外数据采集工作.测线南起雅安市严桥镇，沿北西向经过雅安市芦山县、宝兴县，泸定县，穿过小金县，向北延伸至阿坝县大寨村，横穿四川盆地以及龙门山三大主断裂带并深入松潘—甘孜地块，直线距离近190 km.从四川盆地到青藏高原西南缘，高程变化约4500 m，测点设计间距约5 km,由于地形复杂断裂带多破碎严重，山地居多伴随河流且植被丰茂，沿途有水电站、高压塔等强干扰因素的影响，部分区域测点较稀疏(测点位置如图1中黑色三角形所示).采集仪器为加拿大凤凰公司的MTU-5宽频大地电磁仪，采集频率范围为320～0.0005 Hz，电极采用标准“十”字形布设，平均采集时间为24 h.

图1 龙门山断裂带南段大地电磁测深剖面位置图及区域构造情况简图(a) 青藏高原及周缘块体分布简图.图中HLB—喜马拉雅地块；LSB—拉萨地块；QTB—羌塘地块；S-GB—松潘—甘孜地块；QDB—柴达木盆地；QL—祁连地块; ODB—鄂尔多斯盆地; SCB—四川盆地；YTT—扬子地体；SCFS—华南褶皱系.黑色方框为(b)图范围，红色方框为(c)图范围；(b)龙门山断裂带南段区域板块构造简图.图中块体边界数据来自张培震等(2003)、张国民等(2005)和王辉等(2003).青藏高原及周缘板块主要运动特征据Tapponnier等(2001)和李英强(2018)修改；(c)大地电磁和地震测点位置图.图中红色蓝色圆点为该区域历史地震位置(2000-11-01—2020-01-01)，蓝色代表MS3.0～4.9之间的地震，共计744次；红色代表MS5.0～6.9之间的地震，共计52次.地震目录来自USGS(United States Geological Survey，美国地质勘探局).黑色三角形为大地电磁测点位置；黑色圆点为宽角折/反射地震接收台站位置；红色五角星为宽角折/反射地震炮点位置.紫色直线为图12中所参考人工反射地震资料(李英强，2018)测线位置.汶川MS8.0地震、芦山MS7.0地震、康定MS6.3地震的震源机制解，数据来自GCMT(Global Centroid-Moment-Tensor，全球质心矩张量).断裂带位置(红色实线)来自于邓起东等(2003)，张量分解玫瑰图显示电性主轴方位为近北东向.Fig.1 Topography map showing MT stations layout and tectonic structures in the south part of the Longmenshan fault zone(a)Distribution diagram of Qinghai-Tibet Plateau and surrounding plates. HLB—Himalaya Block; LSB—Lhasa Block；QTB—Qiangtang Block；S-GB—Songpan-Garzê Block； QDB—Qaidam Basin；QL—Qilian Block; ODB—Ordos Basin; SCB—Sichuan Basin；YTT—Yangtze Terrane；SCFS—South China Fold System. The black box shows the scope of figure (b), the red box shows the area of figure (c); (b) Regional plate tectonics in the south part of the Longmenshan fault zone. The block boundary data according to Zhang P Z et al.,2003; Zhang G M et al.,2005; Wang et al.,2003. The main movement characteristics of the Qinghai-Tibet Plateau and its surrounding plates were modified from Tapponnier et al., 2001 and Li et al.,2018;(c) Location map of MT and seismic profile. Red and blue circles are the position of the historical earthquake in this region (2000-11-01—2020-01-01). Blue circles: 744 earthquakes between MS3.0～4.9; red circles: 52 earthquakes between MS5.0～6.9. The earthquake catalog comes from USGS. Black triangles represent the MT sites; Black circles represent receivers of the wide-angle seismic reflection/refraction profile. The red stars represent the shot points of the wide-angle seismic reflection/refraction profile; The straight purple line is the artificial seismic reflection profile shown in Fig.12, according to Li et al., 2018. The mechanism solutions of Wenchuan MS8.0, Lushan MS7.0, and Kangding MS6.3 earthquakes are obtained from GCMT. The location of faults (red solid lines) refers to Deng et al.(2003). The rose diagrams of MT tensor decomposition show that the electrical strike direction orients nearly northeast.

1.3 大地电磁数据预处理

野外采集到的大地电磁数据，通过预处理过程剔除信噪比低、误差较大的频点及数据并将时间序列转换至阻抗张量数据.首先通过快速傅里叶变换将时间域数据转换到频率域.利用ROBUST估计等技术对数据进行筛选，计算得到初始阻抗张量，手工剔除个别畸变的频点数据.对于个别频点或连续频点值分布不正常，或者仪器不稳的情况，运用RHOPLUS理论对其进行视电阻率和相位曲线的一致性校正(谭捍东等, 2004).校正后的数据采用MT-pioneer软件分析电性主轴方向(蔡军涛等, 2010)，利用相位张量方法计算最佳主轴方向角为北东向16°，部分测点的单台全频段玫瑰图(见图1)显示电性主轴方向与龙门山断裂构造走向相符，四川盆地区域测点及小金地区个别测点三维性较强，其余测点二维性较好，测线整体二维性较强.利用MT-editor软件将所有测点旋转角度至电性主轴方向上.典型测点视电阻率与相位曲线如图2所示，其中龙门山至松潘—甘孜地块区域，地形起伏变化剧烈，数据静位移现象明显.

2 速度结构约束下的带地形大地电磁二维反演方法

自从Gallardo和Meju(2003, 2004)提出基于两种物性参数梯度的叉乘获取两种物性的结构相似性以来，该方法在各种地球物理反演约束或联合反演中广泛应用(彭淼等，2013；Peng et al., 2019; Wu et al., 2020；闫政文等, 2020).但已有的研究成果主要针对不带地形的交叉梯度结构约束，而四川盆地至青藏高原测区抬升剧烈，地形起伏高差达5 km，地形是影响反演效果和精度的主要因素之一(王绪本等，1999).因此本文根据龙门山区域实际地形反演的需求开发了带地形的交叉梯度约束算法，并将已有的速度模型作为固定约束加到算法中，进行约束反演.以下介绍算法的实现.

图2 部分典型大地电磁台站校正后的实测视电阻率和阻抗相位曲线红色曲线为TE模式；蓝色曲线为TM模式.Fig.2 Apparent resistivity and impedance phase curves after correction at several typical MT stationsThe red curve represents the TE mode and the blue curve represents the TM mode.

2.1 带地形的交叉梯度项离散方法

根据交叉梯度函数的定义，可以将两种物性的交叉梯度项写成：

(1)

mr表示电阻率模型，ms表示速度模型.在x-z坐标(二维坐标系)中，交叉梯度项分量中只有ty有值，tx和tz为零(吴萍萍，2019).根据矢量叉积求解公式，可将ty写成：

(2)

采用中心差分法将公式(2)离散化成：

(3)

其中i表示z轴第i个网格点，j表示x轴方向第j个网格点.Δxj表示x轴方向第j个网格的网格间距，Δzi表示z轴方向第i个网格的网格间距.

在考虑地形的交叉梯度函数的离散化过程中，主要考虑非边界网格、地形边界、两个侧边界和底界面，跟不带地形的交叉梯度函数离散化主要区别在于地形界面的处理方式.非边界网格、两个侧边界和底界面的计算公式可以参考吴萍萍(2019).以下主要介绍地形边界，公式(3)的计算方式.

在二维地形界面中，以网格(C)为中心，该网格四周有4个网格，分别为上(U)、下(D)、左(L)、右(R).本文主要将地形边界网格分成3种情况：

(1)当C网格周边只有一个空气网格时，则对应的情况有图3(a)、(b)和(c)三种情况.以(a)为例，简要说明该网格的ty的计算方式，图中(i-1,j)为空气网格，则将公式(3)中有关(i-1,j)的物性和网格参数用中心网格(i,j)代替，则公式可以改写成：

(4)

同理，(b)和(c)的公式可以写成如下：

(5)

(6)

(2)当C网格周边有两个空气网格时，则对应的情况有图3(d)、(e)和(f)三种情况.以(d)为例，当(i,j-1)和(i-1,j)为空气时，公式(3)可以改写成：

(7)

同理，(e)和(f)的公式可以写成如下：

(8)

图3 地形边界网格示意图Fig.3 Sketch map of terrain boundary grid

ty(i,j)=0.

(9)

(3)当C网格周边有三个空气网格时，只有一种情况(图3(g))，即(i,j-1)和(i-1,j)和(i+1,j)为空气网格，则该网格的ty可以写成：

ty(i,j)=0.

(10)

在约束反演过程中，ms是已知的速度模型，为常数.将交叉梯度函数的所有网格组合成一个线性矩阵，

ty=Wcg,smr,

(11)

其中mr中只包含非空气的介质电阻率，Wcg,s为常系数矩阵.

取ty的平方和作为结构耦合约束项，则约束项可表示为：

(12)

2.2 基于速度结构约束的带地形大地电磁反演算法

带地形大地电磁反演目标函数由数据项和模型项组成，公式如下：

f(·)是二维大地电磁正演响应函数，λ为拉格朗日因子，Cd为数据协方差矩阵，Cm为模型光滑矩阵，dobs为观测数据矩阵.

将交叉梯度项通过权重因子η加到无约束反演目标函数中，形成速度模型约束下的带地形大地电磁二维反演目标函数，具体表示如下：

+η‖Wcg,smr‖2.

(14)

对目标函数(14)求导数，可得目标函数的梯度：

(15)

其中J为雅可比矩阵,它是正演响应函数f(mr)对模型mr的一阶偏导数矩阵.采用非线性共轭梯度法对加有交叉梯度约束的目标函数(14)进行最小化.在求取公式(15)所示的梯度时，该方法避免了直接求解和存储雅可比矩阵J，先通过解一次拟正演问题求雅可比矩阵的转置和一个向量的乘积，即JT(dobs-f(mr))，进而再按照公式(15)去计算梯度，避免了每次迭代都计算完整的大型矩阵，减少了计算时间和内存需求.图4为本文采用的速度结构约束下大地电磁二维带地形反演算法流程图.

本文在传统的大地电磁二维非线性共轭梯度(NLCG)反演(Rodi and Mackie, 2001)基础上，设计了带地形的网格，将地震资料插值到网格得到约束矩阵，引入了带地形的交叉梯度项，形成了速度结构约束的带地形大地电磁二维反演算法.

图4 速度结构约束下大地电磁二维带地形反演算法流程图Fig.4 Flow chart of the 2D MT inversion algorithm in consideration of both the terrain and the velocity structure constraints

3 理论模型合成数据的反演算例

为了验证本文算法的可靠性和有效性，设计了一个长度为280 km的地堑-地垒推覆构造模型(图5)，地形起伏8 km，沿地表均匀布设46个测点，测点间距6 km，图中以黑色三角形表示.电阻率模型共分为五部分，A部分为高阻高速异常体，A 块体ρA=1000 Ωm，vA=6.0 km·s-1；B部分为低阻低速异常体，B块体ρB=1 Ωm，vB=4.0 km·s-1；C部分为高阻低速异常体，C块体ρC=400 Ωm，vC=4.0 km·s-1；D部分为低阻高速异常体，D块体ρD=10 Ωm，vD=5.5 km·s-1.速度模型中单独设计了一块高速异常体E，vE=6.5 km·s-1，此区域无电阻率异常.A块体上涌推覆至B块体之上，背景场值ρ=100 Ωm，v=5.0 km·s-1.设计带地形的网格，正演计算出频段0.00057～320 Hz之间40个频点的TE、TM数据，加入5%的高斯误差作为观测数据.

反演网格横向采用等间隔剖分，网格间距2 km，纵向为不等间隔剖分，地形部分采用200 m间隔剖分，地形以下按1.5倍间距逐渐增大.初始模型选用电阻率为100 Ωm的均匀半空间.反演结果如图6所示，其中图6a为大地电磁单方法结果，可见带地形的反演方法效果较好，无明显由地形起伏引入的假异常现象，推覆体构造部分形态、位置与模型基本一致.但深部体积效应明显，高低阻分界面及轮廓圆滑不清晰，高阻体相邻两块体反演结果有连通的假象.图6b为速度结构约束下的大地电磁带地形交叉梯度反演结果，相较于单方法，高阻体轮廓恢复较好，与模型区域基本一致，阶梯状推覆构造清晰，高低阻体分界面分辨率有所提高.其中对于块体C与A，常规带地形反演结果有部分连通趋势，通过约束反演后，分离为单个独立块体.对于低阻体部分，块体D轮廓有所改善，块体B顶部轮廓恢复较好.对于速度异常体E(图5b)，此处大地电磁结果无异常，约束后的结果在该区域内仍然未出现异常，说明约束反演只作用于大地电磁有异常区域，不会引入虚假异常.电阻率异常体对应的速度异常体无论是高速还是低速，均有较好的约束效果.

图5 电阻率及速度结构的理论模型图(a) 电阻率模型图； (b) 速度模型图.Fig.5 The theoretical model plots of resistivity structure and velocity structure(a) The resistivity model; (b) The velocity model.

图6 理论模型合成数据的大地电磁二维反演结果对比(a) 常规带地形反演； (b) 本文算法进行的约束反演结果.Fig.6 The comparison of 2D MT inversion results of synthetic data from theoretical models(a) The regular inversion result with terrain; (b) The constrained inversion result of the algorithm in this study.

图7a为大地电磁常规带地形反演的L-curve曲线图，采用了多个正则化因子τ进行多次试算，计算所得RMS与粗糙度分布规律如图所示，最终选定处于曲线拐点位置τ=10时的反演结果.图7b为约束反演过程中，权重系数η=1×1016时，约束反演与无约束反演的RMS下降曲线对比，约束反演迭代次数增加，最终RMS值略高于常规带地形反演.图7c为无约束反演电阻率模型与速度模型的交叉梯度值空间分布图，可见模型边界与两种方法没有共同异常体的区域E，这两部分交叉梯度值较大，即常规带地形反演结果与速度模型在边界区域有较大的非零值.图7d为约束反演与速度模型的交叉梯度值空间分布图，图中交叉梯度值在异常体边界部分明显减小.约束反演有效的提升了大地电磁反演在轮廓边界位置的精确度.同时，参考Moorkamp 等(2011)论文，对于反演模型分别做了模型还原度求解，其公式如下：

(16)

图7 反演RMS曲线及交叉梯度值空间分布图(a) 大地电磁反演L-curve曲线； (b) 无约束和有约束反演RMS迭代曲线对比； (c) 无速度模型约束的反演结果交叉梯度分布图; (d) 速度结构约束的反演结果交叉梯度分布图.Fig.7 RMS curves in the inversion and the spatial distribution of the cross-gradient values(a) The L-curve subplot for MT inversion; (b) The RMS iteration curves for no-constrained and constrained inversion; (c) The spatial distribution map of the cross-gradient values from the inversion results without the velocity structure constraints; (d) The spatial distribution map of the cross-gradient values from the inversion results with the velocity structure constraints.

4 小金—雅安剖面实测数据的反演

4.1 P波速度结构

利用P波速度模型约束实测数据的大地电磁反演中充分考虑了两种方法在分辨率以及勘探深度、测线区域的匹配度.对比分析后认为宽角折/反射地震资料与大地电磁方法更匹配.我们选择了中国地震局地球物理勘探中心在芦山震后布设的一条长约410 km的人工地震宽角折/反射探测剖面的走时反演结果(王帅军等, 2015)，作为约束大地电磁数据的速度结构.地震观测系统，由沿剖面的8个炮点激发和268台三分量PDS-2型地震仪器同步接收组成，接收台间距在0.8～2 km之间.对观测资料进行了二维非均匀介质动力学射线追踪、走时拟合计算，最终获得二维壳幔速度结构模型.其分辨率约为1 km×1 km，勘探深度达到70 km，与大地电磁深度大致相当.由于地震测线较长，故截取其中长度约240 km的P波速度数据(图8)，保证全部覆盖了大地电磁测点区域且向测线两侧各有几十公里延伸.依据大地电磁方法的反演网格，采用不等间隔插值，将速度剖面数据插值到大地电磁反演网格中，约束大地电磁反演.

如图8所示，速度模型呈层状结构，Moho界面清晰.松潘—甘孜地块与扬子地体速度结构迥异，并以龙门山三大主断裂为界.松潘—甘孜地块上覆盖层已部分出露地表，20 km以上至地表附近有一个低速异常区，且中下地壳速度以阶梯状垂向增大，其Moho界面深度在60 km左右.扬子地体地表下沉积盖层较厚，约7～8 km，呈低速特征；Moho界面深度浅，约50 km左右.波速异常区域与大地电磁电阻率异常区域对应较好.

图8 沿大地电磁测线的P波速度模型截面Moho面深度及各区块分布位置参考王帅军等(2015).Fig.8 The section of VP model along the MT profileThe Moho depth and the distribution of each block refers to Wang et al.(2015).

4.2 速度结构约束下的大地电磁二维反演

反演网格部分，纵向地形区域以50 m等间隔剖分，不参与计算的地形网格，采取将实测高程数据插值到网格节点的方法来剔除(图9中红色区域)，地形以下区域网格间距以1.5倍逐渐增大.横向网格剖分采用不等间隔，测点附近加密，测点两侧以1.5倍逐渐加大.在保证了模型计算边界纵向区域足够深、横向跨度足够大的前提下，大大减少了网格块的数量，从而提高了正反演计算效率，减少了内存需求.反演的初始模型参考同区域内大地电磁常规带地形反演结果(张乐天等，2012; 王绪本等，2018)，并通过反复试算后确定选用电阻率为100 Ωm的均匀半空间模型.

加入地形网格后的单独大地电磁反演结果如图10b所示，可见区域内三大主体构造走向、位置清晰，电阻率值过渡光滑.松潘—甘孜地块下方40 km以上普遍存在低阻异常，低阻体以逆冲形态延伸至龙门山断裂带下方大块高阻体之上，龙门山深部高阻异常分布于70 km以上，四川盆地浅层5 km左右为低阻异常带，深部20～60 km左右存在低阻异常区.三大块体具有明显“低-高-低”的宏观电性特征，常规带地形反演结果与王绪本等(2018)在该区域内长周期大地电磁反演结果比对，60 km以上部分主体构造分布基本一致，说明单方法反演结果可靠.加入P波速度模型，进行交叉梯度约束反演后的结果如图10c所示，约束反演后的结果有效克服了大地电磁单方法边界收敛不好的情况，高低阻异常体分布区域及轮廓、走向均有变化.松潘—甘孜块体中下地壳低阻体从50 km至地表均有显示，低阻区域深度有所增加，两块异常体有连通通道.低阻体向龙门山下方高阻体顶部逆冲并深入地表盖层之下，部分已出露至地表.龙门山断裂带高阻体范围较单方法有了明显的收敛，底部边界收敛至50 km左右深度，Moho界面上部，等值线倾斜方向显示与人工地震资料中龙门山三条主断裂的深部延伸产状一致(图12)，更加符合真实地质情况.四川盆地浅地表低阻层约束后薄层特征明显，浅部低阻体与深部低阻体中间有一条形高阻带，高低阻过渡清晰，断裂带构造明显.中下地壳大范围低阻体分布，其中西侧紧邻龙门山高阻体的低阻异常区域向深部延伸至Moho界面以下，东侧低阻体下部为高阻区，低阻分布于40 km以上.

图9 大地电磁反演网格剖分蓝色三角形为测点位置，红色填充部分为地形网格.Fig.9 Mesh discretization for the MT inversionBlue triangles represent the MT sites，and the red fills are terrain grids.

图10 实测数据大地电磁二维反演结果(a) 实际地形剖面； (b) 常规带地形反演结果； (c) 本文算法进行的约束反演结果.Fig.10 2D MT inversion results of the measured data(a) The actual topographic section; (b) The regular inversion results with terrain; (c) The constrained inversion results of the algorithm in this study.

图11 反演RMS曲线及交叉梯度值空间分布图(a) 大地电磁反演L-curve曲线，红色空心圆表示不同τ的RMS值，蓝色五角星表示最终选择值； (b) 无约束和有约束反演RMS迭代曲线； (c) 无速度模型约束的反演结果交叉梯度分布图; (d) 速度结构约束的反演结果交叉梯度分布图.Fig.11 RMS curves in the inversion and the spatial distribution of the cross-gradient values(a) The L-curve subplot for MT inversion; The red hollow circle curve represents the RMS value for different τ, the blue pentagram represents the optimal value; (b) The RMS iteration curves for no-constrained and constrained inversion; (c) The spatial distribution map of the cross-gradient values from the inversion results without the velocity structure constraints; (d) The spatial distribution map of the cross-gradient values from the inversion results with the velocity structure constraints.

图12 综合解释图地形数据来自GMT(Generic Mapping Tools)软件，全球地形数据包earth_relief；断裂带深部产状及延伸位置依据人工反射地震结果(李英强，2018)修改；地表断裂带位置(红线)来自于邓起东等(2003)；沙滩球为芦山MS7.0地震震源机制解，数据来自于GCMT.Fig.12 Integrated interpretation mapTerrain data comes from the global terrain data package (earth_relief) of GMT software. The fault geometry and the extended position underground are modified based on the results of artificial seismic reflection profiles (Li,2018); The location of faults (red lines) refers to Deng et al.(2003); The beachball is the focal mechanism solutions of Lushan MS7.0 earthquake obtained from GCMT.

反演选定的参数、迭代次数与对应的RMS值如图11a、b所示，常规带地形反演试算多个正则化因子，最终选取τ=10的结果，因为曲线没有明显拐点，τ=10与τ=5相比，RMS差别不大但粗糙度相差较大，故选择τ=10的结果更光滑.约束前后RMS曲线对比如图11b所示，RMS值略有增加但差别不大，说明拟合情况较好.图11c与图11d分别为约束反演前、后电阻率模型与速度模型的交叉梯度值的空间分布图，可见约束后整体梯度值都有下降，其中上地壳浅层结构改善显著.

5 构造解释分析讨论

龙门山南段区域断裂带展布较北部更宽，分支断裂多.本次研究测线自西向东穿过的主要断裂带有：硗碛断裂(QQF)、陇东断裂(LDF)、五龙断裂(WLF)、小关子组断裂(XGZF)、双石—大川断裂(SSF)、新开店断裂(XKDF)、大邑隐伏断裂(DYF).以双石—大川断裂为界，西侧为韧性较强的基底卷入构造，东侧为基底卷入变形隐伏于沉积盖层之下(李英强，2018)，硗碛乡至小金段分布多条弧形断裂.

5.1 松潘—甘孜褶皱带逆冲推覆构造形成机制

松潘—甘孜褶皱带位于青藏高原东部，由北部的劳亚板块、西部的昌都—羌塘微板块向上仰冲，东部的扬子地体向下俯冲而成(许志琴等，1991).众多学者研究表明，松潘—甘孜地块岩石圈普遍存在低阻异常带，但低阻层分布不均.对比龙门山中、北段大地电磁结果(詹艳等，2013；彭淼等，2015)，低阻层位置及规模与南段相比差异较大，北段低阻层更接近于地表且深度浅，分布在15 km以上，地表出露较多.中段低阻层规模大，深度更深，位置大约在20～50 km左右，地表出露几乎没有.南段低阻层位置较浅，规模更大，浅地表至50 km范围内均有分布，西侧小金附近有部分出露地表.中部、南部低阻体逆冲趋势明显，覆盖于扬子地体表层并向东侧有一定延伸.

深部地球物理探测中，脱水和部分熔融可用于解释壳内低速-高导层的存在(Meissner and Wever,1986).岩石实验结果证实，少量的熔融体便可以大量吸收弹性波能量，造成地震波速的明显下降.Sato和Sacks(1989)发现地幔橄榄岩中只要有5%的熔融体存在就可以导致地震波速下降约5%.结合该区域较高的纵横波速度比以及存在的壳内负极性震相，说明青藏高原东缘地壳内部存在低黏滞性通道(Zhang et al., 2009)，支持“下地壳流”或“构造逃逸”模型.郭晓玉等(2014)推断北段扬子板块边界位于龙日坝断裂带附近，王绪本等(2018)推测南段扬子板块西缘位于大地电磁结果中的高低阻过渡区.从大地电磁反演结果来看，南段边界符合位于高低阻分界面的判断，但并未有中上地壳物质深入扬子地体下方的趋势，另外一支热流应转向南流入三江块体(熊熊等，2001).综合分析后，推断龙门山南段逆冲推覆构造，形成于下地壳流向东运移，造成了青藏高原隆升，在侧向挤压扬子地体过程中，向东侧运移的热流分支逆冲至扬子地体之上；另一分支向南移动形成弧形断裂带.龙门山断裂带的产生源于印度板块与欧亚大陆发生陆陆碰撞造成青藏高原深部物质运移，运移方向如图1b所示.

5.2 龙门山断裂带南段构造特征

龙门山三条主断裂延伸至南段，宽度逐渐加大，次级断裂增多，破碎严重.大地电磁反演结果显示，薄皮盖层下方高阻块体厚约50 km，三大断裂均分布于高阻体内部，没有高、低阻过渡区.宽角折/反射资料显示，转换带内部震相紊乱、不清晰，速度结构异常呈现低速异常特征(王帅军，2015).人工地震资料显示，断裂带呈叠瓦状排列，交汇于地下20 km处逆冲构造滑脱面，符合逆冲推覆特征.

地震精定位结果显示，芦山地震震中位于芦山向斜下方，双石断裂与新开店断裂附近，但震源深度更深，电阻率结果中，该区域为靠近高阻块体边缘的位置，附近没有汶川地震类似的低阻区，人工地震资料显示大邑断裂带东侧的隐伏断裂从震源附近穿过.稳定的扬子地体，阻挡了松潘—甘孜地块内的物质运移.下地壳流沿大地电磁低阻异常带运移，推测西侧浅层滑脱面沿高低阻过渡带纵向深入至滑脱层.在高原物质运移过程中，挤压扬子地体，驱动了芦山地震的发生.

5.3 四川盆地西南缘构造特征

测线位置穿过四川盆地南端边缘，此处的沉积规律与北部基本相符却更复杂.该区域分布有多处地热资源，推测深部有热源上涌通道.反演结果显示浅地表低阻体分布范围约5～10 km.低阻异常区域内，P波、S波速度结构均为明显的低速异常，Liu等(2014)推断该区域为四川盆地上三叠纪-第四纪沉积物的反映.深部低阻带分为两个区域，靠近龙门山的低阻带有向深部延伸的趋势，推断为川西壳幔低阻带；测线东侧接近边缘部分的低阻体下方为高阻体，推测为川中高阻体边缘.地表低阻体与中上地壳低阻带之间有一层高阻过渡区域，依据人工地震结果划分的隐伏断裂位置，可见多条隐伏断裂存在于此，断裂带沿电阻率变化界面分布，推断芦山地震可能由隐伏断裂驱动产生.

6 结论

本文利用人工地震资料反演获得的P波速度约束大地电磁反演，有效提高了大地电磁方法的纵向分辨率.新算法应用于龙门山南段地区大地电磁实测资料中，得到了该区域更为可靠的电性结构分布.本次研究表明：

(1)大地电磁方法体积效应明显，纵向分辨率低，利用速度模型层状结构及边界清晰的特点，可以有效改善大地电磁方法的体积效应，更清晰刻画异常体轮廓；

(2)带地形的约束反演方法，尤其适用于地形起伏变化剧烈的研究区域开展多资料反演和解释；

(3)推断龙门山南段逆冲推覆构造，形成于下地壳流向东运移，造成了青藏高原隆升并侧向挤压扬子地体，向东侧逃逸的热流分支逆冲至扬子地体之上，另一分支向南移动形成附近弧形断裂带；

(4)推断松潘—甘孜地块低阻异常区下方为深部滑脱层，高原物质运移过程中，挤压扬子地体，驱动了芦山地震的发生.

本文提出的基于地震速度约束的大地电磁反演方法尤其适用于类似龙门山等地形起伏剧烈研究区域.反演结果结合已有物探资料分析了龙门山南段与北段电性结构的差异，推断出扬子板块前缘位于高低阻过渡带附近，分析了芦山地震由隐伏断裂驱动，为龙门山南段的地质构造解释提供了新资料.