不同地球物理方法的反演对地下结构的岩石物理性质表达不同，每种方法只是从一种物性参数的角度来看待地质体，从而只能反映地质体的一个侧面.联合反演不仅能克服单一方法的局限性，减少地球物理反演的多解性，甚至还能降低环境噪声的影响.近十年来，利用地震和电法数据联合反演的优秀研究成果层出不穷.国内外的研究可归为两类:一类是联合反演大地电磁(或直流电法)与折射地震法的初至走时数据^[1-11].由于大地电磁与折射地震单独的二维甚至三维正反演算法都已成熟，因而其联合反演算法也趋向于高维.另一类则是大地电磁数据与远震接收函数或者结合地震面波资料的联合反演^[12]，这类研究同时拟合视电阻率和相位曲线，接收函数和(或)面波群速度频散曲线，主要应用于地球深部壳幔结构的研究.Moorkamp等^[13]首先利用遗传算法联合反演接收函数和大地电磁数据，并采用加拿大Slave克拉通的实测资料，取得了满意的效果，而随后通过加入地震面波资料加强了对地震绝对速度的约束^[14].Zevallos等^[15]也运用远震接收函数和大地电磁数据联合反演研究巴西巴拉那沉积盆地的地下构造.最近，Roux将宽频带大地电磁测深和地震面波资料的联合反演应用于德国中部地区的岩石圈结构研究^[16].

在深部探测中，大地电磁提供的非常低的电阻率一般是地下连通的流体的反映，这通常与宽频地震所获得的低速体相对应，在综合解释中相互验证，互为补充.尽管两者物性存在差异(如地下介质的波速大小取决于岩石组分和矿物类型，而电阻率对地下导电体更为敏感)，但岩石圈内大的地层界限(如Moho 面)总会导致波速和电阻率同时变化. Kozlovskaya 对岩石圈内岩石的数值模拟结果表明:当岩石孔隙通道被水或部分熔融物充填时，该岩层总是表现低阻伴随低速的物性特征^[17].而这一兼容性成为本文运用接收函数和大地电磁资料联合反演研究的前提.

南迦巴瓦地处喜马拉雅山脉东端，这里地质界限、构造格架和地貌水系发生急剧转折，在构造地质学上被称之为喜马拉雅东构造结，是研究的热点地区.区内划分为三个构造单元:冈底斯(拉萨)地体、 雅鲁藏布江缝合带和印度地体(图 1).构造结外侧为欧亚大陆岩系，其核部变质地体为印度大陆高喜马拉雅岩系.构造结变质体被认为发育叠置的岩片构造，呈现十分复杂的复式背形构造样式.自NW 到SE 由比鲁构造岩片、直白构造岩片、南派乡构造岩片和多雄拉变质穹隆组成，其中由高压麻粒岩相组成的直白构造岩片被韧性拆离断裂和逆冲断裂所夹持，为挤出构造岩片.组构运动学的研究认为:南迦巴瓦变质体相对拉萨地体的运动转为水平走滑运动，并向北插入拉萨地体之下^[18].南迦巴瓦变质体中大片高压麻粒岩和石榴辉石岩的发现揭示出古近纪的高压(超高压)变质作用，认为印度板块向欧亚板块之下的俯冲深度达到80~100km^[19].尽管南迦巴瓦的构造格架已初具雏形，但其形成和演化过程及其动力学机制还有待进一步研究.

图 1

Fig. 1

图 1 南迦巴瓦宽频地震流动台站和大地电磁测点位置示意图图中构造划分依据南迦巴瓦构造图[18],左上角图构造划分依据青藏高原结构图[20]. Fig. 1 Sketch map showing locations of broadband seismic stations and magnetotelluric sites in Namche Barwa, Tibet The figure of tectonic division is based on tectonic map of the Namche Barwa[18] and the upper-left figure of tectonic division is based on the map of the structural architecture of the Qinghai-Tibet Plateau[20].

在深入分析Moorkamp算法的基础上，本文拟通过对南迦巴瓦变质体中部四个宽频地震流动台站和同位置大地电磁测深点数据的联合反演研究，获取深部壳幔的电性和横波波速结构特征，以探讨南迦巴瓦构造结中部的岩石圈构造及其地球动力学机制.

2010年4 月至10 月，项目组在西藏林芝地区开展了南迦巴瓦深部构造的探测工作，围绕南迦巴瓦峰分阶段实施了宽频带地震流动台站的布设和大地电磁测深数据的野外采集工作.大地电磁数据观测使用加拿大凤凰公司生产的V5-2000 大地电磁系统，期间共采集83 个点的320~0.0005Hz频段内长周期数据，测点间距约5km；宽频带地震观测采用美国Ref-Tek 130 地震仪，配备STS-2 和Trillium-40检波器，期间共布设三分量地震台站24 台次，台站间距约20km,其中部分地震台至今仍在持续观测中(图 1).考虑到部分台站观测时间有限，挑选了数据量较足、质量较好的郎嘎、崩嘎、拉格和直白四个宽频地震台分别与同位置高质量大地电磁测深点数据做一维联合反演研究.

大地电磁法(Magnetotellurics)是一种天然源的频率域电磁法，它以天然平面电磁波为场源，通过在地表观测相互正交的电磁场分量来获取地下地电构造信息.近年来大地电磁测深的正反演技术及资料处理方法已经发展相当成熟，成为研究地壳和上地幔结构的主要方法之一.

大地电磁均匀各向同性介质的构造假设通常在实际资料处理中给反演结果带来畸变.本文则通过引入别尔季切夫斯基不变量(Berdichevskiy invariant)作为联合反演的MT 参数来减小一维构造假设对反演结果的影响.它是一种对波阻抗的非对角元素求算术平均的方法:

(1)

数值模拟结果已证明，只要地下介质的横向变化不太剧烈，其反演结果是可靠的^[21].当然，该不变量同样可以分别通过波阻抗与电阻率(相位)的关系式(2)计算相应视电阻率和相位值:

(2)

其中，Z是波阻抗；ω 是电磁场的圆频率(ω=2π/T，T为周期，单位为s)；μ 是磁导率(为常数，单位为H/m).

实际资料处理时，对电磁场各个分量的时间序列进行频谱分析，并采用Robust估计方法估算阻抗张量元素.为消除背景噪声干扰等因素引发的曲线畸变，采用Rhoplus理论^[22]检查视电阻率和相位资料的一致性.以拉格测点为例，图 2为采用Rhoplus 分析获得的视电阻率和相位曲线.可见与原始数据相比，校正结果很好地压制了不相关噪声的影响，为联合反演提供了高质量的大地电磁测深数据.

图 2

Fig. 2

图 2 拉格台Berdichevskiy不变量的视电阻率和相位及其Rhoplus分析 Fig. 2 Berdichevskiy apparent resistivity and phase of Lage station and its Rhoplus analysis

接收函数(Receiver Function)是研究地壳、上地幔速度精细结构的一种天然地震处理方法.它是去除了震源、地震波传播路径以及仪器响应等因素后的时间序列，主要包含地震台站下方地壳和上地幔速度间断面所产生的转换波及其多次反射波的信息.由于对横波波速的垂向变化敏感，因而能被用来反演台站下方垂向横波速度结构.本文通过提取P 波接收函数做联合反演.

P波接收函数需从三分量远震P波波形中提取关于接收介质的水平响应，用远震记录的垂直分量对水平分量作反褶积得到.其计算方法众多，主要分为两类:一类是以水准量反褶积(Water Level Deconvolution)^[23]为代表的频率域类方法，这类方法比较简单，但算法不够稳定，效果不佳.另一类是在时间域进行反褶积，近年来这一类方法在接收函数研究中较为流行.本文引入一种时间域迭代反褶积法(Time Domain Iterative Deconvolution)^[24].

在接收函数偏移成像研究中，一般对台站不同方位和震中距的接收函数进行偏移叠加以增强信噪比并消除散射波和背景噪声等影响.但对于反演来说，进行接收函数的直接叠加处理反而不合适.因为后续多次波包含地下结构信息并有助于对速度间断面的约束，而叠加会减小甚至消除多次波的振幅.一种较好的处理则是通过选择一定震中距和反方位角范围内的接收函数进行叠加平均.

处理时，首先对三分量的原始记录作去均值和去线性化趋势处理，并进行0.05~0.8Hz的带通滤波；然后截取直达P 波前30s为起点、后90s为止的波形数据；对水平分量作坐标旋转得到径向和切向分量；最后，通过选取3.0的高斯滤波因子对垂直分量和径向分量作时间域迭代反褶积，得到接收函数.以直白台的数据处理为例，共处理了2010 年5 月至8月震级大于5.2级的远震事件84个，并从中筛选出反方位角在62°至67°之间，震中距在39.5° 至40.6°之间的接收函数作叠加平均.图 3为筛选出的21条高质量接收函数，其中最下方为这些接收函数的叠加平均结果.联合反演时，将采用这一叠加平均接收函数作为实测数据.

图 3

Fig. 3

图 3 直白台宽频地震接收函数及叠加平均 Fig. 3 Receiver functions of Zhibai station and stacking average result

本文对远震接收函数和大地电磁数据做一维联合反演选用求解多目标优化问题的非优超排序遗传算法(Non-dominated Sorting Genetic Algorithm II，缩写为NSGA-II)^[25].虽然这一算法相比各类线性化反演方法更为耗时，但是可以避免陷入局部极小，能够获得更为可靠的反演结果.而且，NSGA-II 的优势还体现在:一次运行最终能得到一组优化模型集，并通过拟合差权衡图来衡量不同数据间的拟合程度，从而能直观地寻找到最佳模型.

对于线性化反演而言，总目标函数一般写成各分量加权相加的形式，是单一目标的优化算法，如各类地球物理正则化反演就是最典型例子.一般来说，单一目标的优化算法在处理多种数据时往往表现欠佳，而采用多目标优化解决此类问题成为当今发展趋势.这种算法不需要事先指定各目标函数的权重，本文中，多目标优化的总目标函数写成如下分量形式:

(3)

NSGA-II采用Pareto排序策略对每一代模型群体搜索排序，而且保留不同目标函数分量的数值信息.Pareto排序策略可以描述为:对于任意模型a和b，若a的偏序小于b(a$ \prec $b)，当且仅当模型a所有目标函数分量均小于或等于模型b各目标函数分量值，而且模型a至少有一个目标函数分量小于模型b相应目标函数分量值时才成立.

(4)

在方程(4)中，可见该算法通过比较模型间各目标函数分量值的大小，将目标函数最小化反映到优化模型的选取上.一般而言，满足上述方程的模型是一组模型的集合，这组模型集合称之为Pareto最优前沿.在遗传算法中，每一代的模型群体都通过这种选择方式获得Pareto最优前沿，同时选取排序级别靠前的模型作为新一代种群，并反复迭代，朝更优的Pareto最优前沿演化.最终，只需从最后迭代的最优前沿中评价不同数据间的拟合权重以寻找最佳模型.理论模型试验表明，当两种反演模型兼容时，结果模型很少，拟合差权衡图的分布较为集中，具有点、线的形态特征；而当模型不兼容时，结果模型较多，且分布散乱，无法找到统一模型同时使两套数据都较好拟合^[14].应该指出，NSGA-II算法所获得的反演结果给我们提供了评价地下导电性和横波速度结构模型之间兼容性的有利信息.我们不仅可以通过权衡图的分布形态了解数据是否兼容，而且还能通过改变不同拟合差权重来分析相应模型的变化情况.

以上着重阐述了反演算法的模型选择策略，而在遗传算法的实现上首先要产生初始种群，这就需要对模型参数进行二进制编码.二进制编码的长度取决于模型参数的范围和分辨率，因而在反演前需设定每一模型参数的最小值m_i^min ，最大值m_i^max和分辨率Δm_i.然后，对初始种群进行模型繁殖，依次经过交换、变异和更新三步遗传操作使模型进行有效演化.经过多次试验，本文遗传算法采用的交换概率和变异概率分别设为0.6 和0.2，以保证算法搜索效率和质量达到最优.

还应指出，本文中电阻率和波速模型的耦合建立在具有相同的地层界面上，但又不严格限制电阻率值(波速值)相对波速(电阻率)变化在相邻层间必须改变.这种结构化的耦合关系类似于一种交叉梯度(cross-gradient)的耦合方法^[1-3]，它要求不同方法的模型值只能朝着相同方向发生改变，或者其中一种模型值改变而另一种保持不变.该方法是一种松约束的耦合方式，已更多地运用到大地电磁与地震折射层析的二维甚至三维联合反演^{[1-3, 6]}.

联合反演采用的基因算法种群数(每次迭代的模型个数)设为1000，迭代200次.联合反演时需设置每一层的层厚、地震波横波速度以及电阻率参数，并且每种参数都由最小值、最大值和步长组成(表 1).本文通过比较多次反演结果的权衡图及其最佳模型，选择其中最稳定可靠的一次作为联合反演的结果.一般来说，增大联合反演模型层数能同时减小大地电磁和地震接收函数反演的拟合差，但同时也会使电阻率模型和波速模型的耦合变弱.考虑到这一特性，本文通过比较不同层数的联合反演结果来分析层数对联合反演的拟合差和模型耦合的影响. 每一个台站都将采用7层、9层和11层分别做联合反演.图 4为各台站联合反演结果模型的大地电磁与接收函数拟合差(MT RMS-RFRMS)权衡图.下面依次对郎嘎台、崩嘎台、拉格台和直白台的权衡图及其最佳联合反演模型作分析和讨论.

图 4

Fig. 4

图 4 联合反演结果模型的大地电磁与接收函数拟合差权衡图图中分别以三角形、星号和五角星标识所选取的7层、9层和11层联合反演模型:(a)郎嘎台；(b)崩嘎台；(c)拉格台；(d)直白台. Fig. 4 Trade-off plots of joint inversion of magnetotellurics and receiver function In the figure, the 7-layer joint inversion model is represented by a triangle，the 9-layer model by a asterisk and the 11-layer model by a pentagram: (a) Langga station; (b) Bengga station; (c) Lage station; (d) Zhibai station.

表 1(Table 1)

表 1 遗传算法联合反演参数列表 Table 1 Parameters of genetic algorithm for joint inversion

表 1 遗传算法联合反演参数列表 Table 1 Parameters of genetic algorithm for joint inversion

图 4a中从左至右分别代表郎嘎台7层，9层和11层的拟合差权衡图，其中9 层和11 层模型的权衡图优化模型分布形态较为相似.大量模型的MT 拟合差(RMS)小于2，部分模型甚至低于1.5，同时接收函数RMS 也基本控制在2 以下，表明数据的拟合情况良好.对比9 层与11 层权衡图看到，层数增加没有同时明显降低MT 和接收函数的RMS.而与此同时由于层数增加势必会使联合反演模型的耦合变弱，因此适合采用9层模型表达郎嘎台地下结构.以拟合差同时最小为选取原则，从7层，9层和11层的拟合差权衡图中分别选取一个模型作为该层的联合反演模型.从模型结构图看到:三个不同层模型差异不大，只在深部地幔出现差异；都显示在地下5km 至15km 之间含有高速层，而且接收函数(见图 6)在1s至2.5s之间也显示正振幅，证实了高速层的存在，同时电阻率模型指示该层电阻率值也相对偏高.此外，中地壳和下地壳(25~60km)的电阻率相对偏高；波速结构指示Moho面位于地下60 km 深度处，波速从下地壳的3.5km/s增大到上地幔4.5km/s左右，其中11层模型显示下地壳与上地幔之间存在一个过渡层；上地幔顶部岩石圈的电阻率总体偏低，接近10Ωm；上地幔约130km 深处横波速度急剧下降，但11层速度模型指示另一可能的结构:分别在115~140km 以及约180km 以下存在两个低速层.

图 5

Fig. 5

图 5 各台站联合反演模型 各图中左侧为大地电磁的电阻率模型，右侧为地震横波速度模型:(a)郎嘎台；(b)崩嘎台；(c)拉格台；(d)直白台. Fig. 5 Joint inversion models of all stations In each figure，resistivity model show on the left side and shear wave velocity model on the right side: (a) Langga station； (b) Bengga station； (c) Lage station； (d) Zhibai station.

图 6

Fig. 6

图 6 各台站9层模型的正演响应与观测数据对比 图中左侧为大地电磁视电阻率和相位，右侧为接收函数:(a)郎嘎台；(b)崩嘎台；(c)拉格台；(d)直白台. Fig. 6 Comparison of predicted and measured data of the 9-layer model for each station In the figure, magnetotellurics show on the lett side and receiver functions show on the right side: (a) Langga station;b) Bengga station;c) Lage station;d) Zhibai station.

崩嘎台各拟合差权衡图(图 4b)都显示接收函数最低RMS均明显高于大地电磁最低RMS(接近1)，而且与其他台的接收函数最低RMS 相比也相对较高.其中7层和9层模型的分布较为相似:优化模型在权衡图中具有分化的特征(MT 拟合较好模型与接收函数拟合较好模型在权衡图上出现分离)，这表明该测点电阻率与地震波速模型可能不兼容. 三个模型(图 5b)都显示地下14km 以上存在高阻高速层，而中地壳(20~40km)电阻率和波速值则保持相对稳定.从各横波速度模型均可以推断Moho面在接近60km 处，而位于Moho面以上包含约15km 厚的过渡层.9层模型和11层模型均指示在上地幔约150km 以下存在低速区，该深度以下的电阻率也有下降趋势，但在何处下降各层模型却没有给出较统一深度，这可能与大地电磁的探测深度不足有关.

拉格台的7层、9层和11层的拟合差权衡图的分布都较为相似，基本都能找到使大地电磁RMS 和接收函数RMS同时最小的最优模型(图 4c).部分较优模型的MT 和接收函数拟合差都控制在2 以内，说明数据的拟合情况良好.层数的增加并未显著减小拟合差，表明7 层模型已足以同时拟合好MT 和地震数据.分别选取权衡图中“L"拐角处模型作为拉格台7 层、9 层和11 层的联合反演模型. 联合反演模型(图 5c)具有一定相似性:地下9km 以上有一高阻高速层(地震横波波速达到3.6km/s,而电阻率超过10000Ωm)；从波速结构推断Moho 面位于60km 附近，下地壳含有约20km 厚度的过渡层；电阻率值在中下地壳总体偏低，维持在100至200Ωm 之间；进入上地幔顶部电阻率值略微上升，直到约96km才开始急剧升高，达到2000Ωm 以上； 9层和11层模型均显示上地幔岩石圈中100~120km 深度内具有低速结构特征，同时在大约150km 以下，横波速度显著降低.

直白台的权衡图中大量优化模型集中在直线(7层)或“L"形折线(9层和11层)上(图 4d)，这样的分布特征说明电阻率模型和速度模型比较兼容. 直白台接收函数的RMS 最低能到1.5 以下，相比其他台拟合得更好.层数增加对MT 和接收函数RMS的改善不大，说明7层模型已能较好地表达电阻率和速度结构.三个模型(图 5d)都显示:在地表 4k m以下、16km 以上具有低阻低速层，电阻率最低到10Ωm 以下，而波速最低至3km/s；中地壳电阻率值相对偏高，达200Ωm(或以上)；下地壳同样存在厚度达20km 的壳幔过渡层，至60km 附近到达其底界(Moho面)；上地幔顶部(60~90km)的电阻率和速度结构都比较稳定，但在约140km 以下速度突然下降，这一特点在上述朗嘎、崩嘎和拉格的速度结构都有所表现，因此推测可能为岩石圈和软流圈的分界面.

对9层联合反演模型中的电阻率模型和速度模型分别做一维正演计算可得到相应模型的正演响应，并将其与各台站观测数据一一对比.对比结果表明，图 6所显示的各地震台站接收函数的主要波形特征都得到很好拟合，而大地电磁反演模型的正演响应也与实测数据基本吻合.各台站联合反演模型都同时较好拟合了大地电磁和接收函数实测数据，说明结果稳定可信.

将四个台站的结果对比来看，其相似结构特征主要有:(1)壳幔过渡层底界的埋深非常一致，约为60km.(2)下地壳都存在15km 至20km 厚的壳幔过渡层，该层地震横波波速在4km/s左右，而直白和拉格的电阻率模型都显示为高导的电性特征.(3) 除直白台外，郎嘎台、崩嘎台和拉格台的上地壳都存在约9km 至14km 厚的高阻高速层；在中地壳的15km 至40km 范围内都含有低阻低速层，但中心埋深和层厚不一，平均中心埋深在28km 左右，平均厚度约为20km.(4)各台上地幔约130~150km 之间存在速度急剧下降界面，推断为岩石圈与软流圈分界面.

位于南迦巴瓦邻近地区的错那至墨竹工卡剖面的大地电磁已有研究成果表明:印度斯-雅鲁藏布江缝合带以南存在延深近20km 的巨厚高阻层，其下发育较大规模的良导体，中心埋深在30km 左右^[26].本文的联合反演结果不仅符合这一电性结构特征，而且其地震波速也同样表现为类似的变化，说明藏南雅江缝合带以南地区的中上地壳可能普遍存在高速高阻层之下的低速高导体，而其中的低速高导体可能与深部流体和局部熔融共同作用有关.

徐果明等^[27]采用Rayleigh 面波研究青藏高原的三维横波速度结构，其中的a-a′纵剖面横穿南迦巴瓦，层析成像结果显示南迦巴瓦的下地壳的横波速度明显偏高；上地幔存在较为明显的低速软流层，其顶部埋深大致在120km 附近，中心深度位于160km 附近.Chen等^[28]给出的面波层析成像结果中B-B′ 南北向纵剖面也显示普遍存在约20km 厚的壳幔过渡带，而且29°N 附近地区的上地幔150km 以下也发现横波波速约4.1km/s的低速软流圈.本文的联合反演结果与上述两者的这两个波速特征十分吻合，进一步增强了该区存在壳幔过渡带和低速软流圈的可靠性.此外，本文显示的壳幔过渡层偏低的电阻率这一结果可能还对过渡层的成因研究具有很好启示.

经分析，郎嘎、崩嘎和拉格的上地壳所具有的高阻高速层很可能对应于多雄拉组混合岩化角闪岩相变质岩，而联合反演结果显示该岩层是有根的，其底界不超过20km.值得注意的是，直白台所显示上地壳4km 至16km 之间为低阻低速层的结果与其他三个台站完全不同.根据直白构造岩片为高压麻粒岩及其典型的“挤出岩片"的构造特征^[18]，我们认为低阻低速层与直白组高压麻粒岩的少量部分熔融有关，其来源可能与下地壳壳幔过渡带中的镁铁质岩石及其高压(超高压)变质作用产生的麻粒岩和榴辉岩的向上运移有关，甚至能更深地追溯到幔源玄武质岩浆的底侵作用.

本文利用南迦巴瓦中部地区郎嘎、崩嘎、直白和拉格4个台站的宽频地震和大地电磁数据联合反演得到了各台站下方的一维速度结构和电阻率结构. 联合反演中通过采用同层的松约束耦合方式使速度和电阻率在地层界面间同时发生变化，其反演结果不仅同时较好地拟合了观测数据，而且各台速度和电阻率的变化在部分地层也相当一致.从4 个台的横向对比来看，其壳幔结构具有一定相似性:

(1) 郎嘎台、崩嘎台和拉格台的上地壳都存在9km至14km 厚的高阻高速层，在中地壳的15km 至40km 范围内都含有低阻低速层，这种高阻层覆盖于低阻层之上的电性结构在南迦巴瓦以西的邻近地区已有报道，表明藏南雅江缝合带以南地区的中上地壳可能普遍存在这种结构，其低速高导体可能与深部流体和局部熔融共同作用有关.

(2) 各台的下地壳都存在15km 至20km 厚的壳幔过渡层，而且该层直白台和拉格台电阻率偏低； 过渡层与典型上地幔的分界面比较一致，均在60km左右.上地幔约130~150km 存在速度急剧下降界面，解释为低速软流圈的存在.这一结果与前人面波层析的邻近剖面的结构基本一致.

(3) 郎嘎台、崩嘎台和拉格台的上地壳所指示的高阻高速层应对应于多雄拉组混合岩化角闪岩相变质岩，而直白台所显示4km 至16km 之间的低阻低速层可能与高压麻粒岩的少量部分熔融有关.我们认为这可能源于壳幔过渡带镁铁质岩石的相变或更深处幔源玄武质岩浆的底侵作用.

需指出，本文的研究是基于地壳水平分层的情形下获得的，但由于南迦巴瓦地区中上地壳构造的复杂性和横向不均匀性，倾斜岩层必定存在，这势必给反演结果带来一定误差.因此，进一步的研究任务将克服一维构造假设的局限性，开展二维甚至三维天然地震和大地电磁联合反演的研究.

本文研究使用了德国学者Dr.Moorkamp 提供的联合反演程序，并得到他的热心帮助和有益的建议，史大年研究员在数据处理上给予了指导，匿名评审专家为本文修改提出了细致而宝贵的意见，在此一并表示感谢.