基于背景噪声和地震面波联合反演华北克拉通中部岩石圈结构＊

黄 翔 丁志峰 宁杰远 徐小明

1）中国北京 100081 中国地震局地球物理研究所

2）中国北京 100871 北京大学地球与空间科学学院

引言

华北克拉通是世界上最古老的太古代克拉通之一，通常认为其主要由三部分组成，包括东部、西部两个太古代地块以及两者在古元古代相互碰撞缝合形成的中部造山带（Zhaoet al，2005）.华北克拉通在早古生代依然保留着典型的太古代克拉通特性（Griffinet al，1998），然而自晚中生代以来经历了剧烈的热构造活化和不同程度的破坏减薄.一般认为，在太平洋板块的西向俯冲作用的影响下，华北克拉通东部地块，即华北盆地，经历了剧烈的构造活化和岩石圈破坏减薄，而西部的鄂尔多斯盆地则未遭受到构造活化作用的破坏，依然保留着克拉通特性，位于鄂尔多斯盆地和华北盆地之间的中部造山带则可能遭受了不同程度的局部改造（Griffinet al，1998；Menzieset al，2007；Zhuet al，2011，2012）.

大量的岩石学、地球化学等研究表明，华北克拉通东部地块在晚中生代表现出伸展作用强烈、热流值大、岩浆活动强烈等特点（Menzieset al，1993；Griffinet al，1998；Fanet al，2000），巨厚的克拉通岩石圈（＞180 km）经受了强烈的破坏改造，导致大部分区域岩石圈厚度显著减薄（＜100 km），失去了典型的克拉通特性.而最近的地震学研究显示，华北克拉通中部造山带也存在局部的改造和减薄（Chen，2010；Tanget al，2013；Aiet al，2019）.从晚中生代至新生代，中部造山带表现出一系列的伸展和裂谷作用以及广泛的岩浆活动.前人在该研究区域进行了大量的研究，揭示了很多重要的信息.如Tang等（2013）基于背景噪声和地震面波成像得到了华北克拉通中部的岩石圈厚度模型，提出华北盆地的岩石圈受到的破坏和减薄比较严重，厚度约为60—65 km，并提出岩石圈破坏和减薄不仅局限于华北盆地，可能对中部造山带和鄂尔多斯盆地也有显著的影响，包括河套、渭河裂谷等.Ai等（2019）基于背景噪声和地震面波成像获得了华北克拉通中部的岩石圈模型，提出中部造山带北部遭受的岩石圈改造相对于南部更加显著且造成两者的原因不同；根据Lei （2012）基于P波成像的结果显示大同火山群下方的低速体穿过太平洋停滞板块中的缺口延伸至下地幔，认为北部的岩石圈减薄原因来源于下地幔，而南部的岩石圈减薄则是由停滞的板块脱水反应引起的.

大同火山群位于山西裂谷北端，是华北克拉通最大的第四纪板内火山之一，最早喷发于早更新世，岩浆活动持续至晚更新世，岩浆规模相对较小（陈文寄等，1992；Xu，Ma，1992；Xuet al，2005）.作为华北克拉通重要的构造单元之一，大同火山群的岩浆活动对华北克拉通破坏以及山西裂谷演化等过程起着至关重要的作用，因此，其形成机制与岩浆深部来源一直是研究人员关注的热点问题，但目前对于相关问题的解释仍存在争议：Huang和Zhao（2006）基于P波走时成像结果提出大同火山群可能与太平洋板块的深部俯冲脱水有关；Lei（2012）基于远震P波到时层析成像结果认为大同火山的岩浆可能来源于下地幔深部上涌的地幔柱；Li等（2018）基于背景噪声和远震面波层析成像结果认为大同火山的岩浆可能来源于青藏高原东北缘下的软流圈.

由于构造背景复杂而且地震台阵不够密集等原因，目前关于华北克拉通的破坏机制和动力学过程认识还不够充分，对于华北克拉通岩石圈精细结构和破坏程度等问题仍存在争议.例如：Huang等（2009）基于面波成像计算的鄂尔多斯盆地岩石圈厚度约为160 km；Tang等（2013）基于S波速度与温度经验公式计算的鄂尔多斯盆地岩石圈厚度约为120 km；汪洋和程素华（2011）基于局部地壳均衡分析得到鄂尔多斯盆地岩石圈厚度为130—140 km，认为下方不存在岩石圈根；这些研究结果之间都存在一定的差异，且与Chen等（2009）基于接收函数得到的结果（＞200 km）相差较大.在地球动力学中，岩石圈通常被定义为对流地幔的热边界，岩石圈底部一般与固定的温度相对应，而在50—250 km深度范围内，地震波的速度主要与温度相关（Goeset al，2000；An，Shi，2006）.因此，利用更加密集的台阵观测资料来获取该区域精细的岩石圈速度结构对于理解华北克拉通岩石圈破坏机制和破坏程度具有重要意义.另外，地幔热物质的上涌通常会产生局部高温，进而造成地震波的低速异常，因此，我们可以利用更加密集的台阵观测资料来获取该区域的高分辨率三维速度模型来研究大同火山群岩浆活动在地壳内的上涌通道及其深部来源.

本研究拟基于ChinArray三期一阶段覆盖华北克拉通中部造山带的密集流动台阵观测数据，利用背景噪声互相关和地震面波层析成像获得的瑞雷面波频散信息，联合反演以期获得研究区域岩石圈尺度精细的三维 S波速度模型.相对于前人的研究（Jianget al，2013；Tanget al，2013；Aiet al，2019），本研究的速度模型有望达到更高的分辨率，可揭示华北克拉通中部岩石圈速度结构的更多细节，以期进一步探讨华北克拉通岩石圈的厚度以及改造减薄程度的局部差异以及破坏机制，并研究山西裂谷北端的大同火山群岩浆活动在地壳内的上涌通道及其在地幔中的深部来源.

1 数据与方法

本文使用中国地震科学台阵探测项目（ChinArray）三期一阶段在华北克拉通中部布设的流动台阵（中国地震科学探测台阵数据中心，2011）的观测数据来进行华北克拉通中部的岩石圈结构成像研究.其中，背景噪声层析成像的数据集采集自2017年1月1日至7月28日台阵记录的连续波形.地震面波程函方程成像的数据集采集自2016年11月1日至2019年1月21日台阵记录的地震事件波形.台阵统一使用了CMG-3T地震计（频带为50 Hz—120 s）和REFTEK数据采集器，平均台间距约为30 km，台站位置分布如图1b所示.该台阵覆盖了北克拉通中部造山带及其西边的鄂尔多斯地块东部和其东边的华北盆地西部，覆盖范围大致为（34°N—44°N，109°E—115°E），是目前为止布设在该研究区域最为密集且均匀的台阵，有助于获取该区域更为精细的岩石圈结构.

图1 研究区地形、构造背景（a）及本研究使用的台阵分布（b）Fig. 1 Topography and tectonic settings of the studied area （a）and the distribution of array used in this study （b）

1.1 背景噪声成像

背景噪声成像的数据预处理过程主要参照Bensen等（2007）提出的方法：包括重采样至1 Hz、去均值、去线性趋势、去仪器响应和带通滤波等处理；然后使用滑动平均法进行波形的时域归一化处理，并进行频域谱白化；最后对每两个台站对每天的波形记录进行互相关计算.

为提高波形信噪比，本文采用两步叠加法进行波形叠加处理：① 将每五天的互相关波形进行线性叠加；② 使用Li等（2018）发展的时频相位加权叠加法（time-frequency domain phase-weighted stacking，缩写为tf-PWS）将每五天的线性叠加结果进一步叠加得到台站对之间最终的互相关波形.为了进一步提高信噪比，将每两个台站对之间的互相关波形的因果信号与果因信号进行对称叠加用于频散曲线的提取，然后使用Levshin和Ritzwoller （2001）发展的自动频时分析方法（frequency-time analysis，缩写为FTAN）进行频散曲线的提取.为控制频散曲线的质量，剔除台站间距小于3倍波长，信噪比小于20的数据，得到的不同周期的频散数量如图2所示.

图2 不同周期的瑞雷面波相速度频散数量Fig. 2 The number of Rayleigh wave phase velocity dispersions at different periods

经过初步质量控制后，我们使用Barmin等（2001）发展的一种基于线性射线理论的阻尼最小二乘反演方法来反演二维瑞雷面波相速度.为了进一步剔除频散中可能存在的异常值，我们进行了两次反演计算：① 对初步得到的频散曲线进行初始的二维反演，得到较为光滑的相速度模型；② 把走时残差大于4 s的测量值剔除，将筛选后的频散再次进行二维反演得到最终的0.25°×0.25°的二维瑞雷面波相速度模型.

1.2 地震面波成像

为获取更长周期的面波速度，本文使用Jin和Gaherty （2015）发展的基于程函方程的面波成像方法进行地震面波层析成像.该方法利用相邻地震台站之间记录的地震波形的互相关函数来计算台站间的相位差，然后基于程函方程将相位差反演成二维相速度模型.本文收集了2016年11月1日至2019年1月21日台阵记录的979个地震事件的垂向波形.数据的前期处理工作主要包括：① 参考前人研究（钟世军等，2017），对地震事件进行了初步筛选，要求地震震级不小于MS5.5，震中距在10°—160°之间，震源深度小于50 km，得到了592个符合条件的地震事件，事件的方位分布见图3；② 对符合条件的地震事件波形进行预处理，处理过程与上述背景噪声数据预处理一致.

图3 地震事件相对于台阵中心的方位分布Fig. 3 The azimuth distribution of seismic events relative to the center of the array

面波相速度计算主要包括以下步骤：① 对于每个事件−台站对，截取波形记录中的基阶瑞雷面波能量；② 对于同一事件，对一定距离内（200 km）的不同台站记录的波形进行互相关计算；③ 对互相关波形进行窄带通滤波，滤波带宽为中心频率的10%；④ 使用包含振幅因子、窄带波形的半带宽、窄带波形的中心频率、群延迟和相延迟五参数的高斯小波来拟合该互相关波形，从而获得两个台站间的相位延迟，图4为波形拟合示意图；⑤ 利用得到的台站间的相位延迟反演慢度矢量；⑥ 依次利用程函方程和亥姆霍兹方程获得研究区域内0.25°×0.25°网格的二维相速度模型（Jin，Gaherty，2015）.相比于使用传统的双台法计算相速度，基于程函方程的面波成像方法不需要进行射线追踪，允许射线传播方向存在弯曲，更符合真实情况，能够得到更可靠的面波相速度模型（Linet al，2009）.但该方法对台阵要求较高，需要台阵分布相对密集且均匀，本研究使用的台阵满足这一条件.

图4 互相关波形拟合示意图Fig. 4 The fitting of the cross-correlation waveform

分别利用背景噪声互相关和地震面波计算得到不同周期的瑞雷面波频散后，将基于两种数据集计算的频散曲线进行拼接.考虑到背景噪声成像可以获得可靠的较短周期的相速度而地震面波成像可以获得可靠的较长周期的相速度，经过对两种数据集获得的频散曲线的观察和测试，我们选取了6—15 s周期的基于背景噪声互相关计算的相速度和30—140 s 周期的基于地震事件程函成像的相速度作为各周期对应的相速度值.对于20 s和25 s的重叠周期，参考Fadel等（2020）的方法采用了3 ∶ 7的比例进行加权平均，即在20 s，使用0.7权重的背景噪声互相关计算结果和0.3权重的地震事件程函成像计算结果的和，而在25 s，使用0.3权重的背景噪声互相关计算结果和0.7权重的地震事件程函成像计算结果的和.最终得到了研究区域6—140 s的二维瑞雷面波相速度模型用于进行S波速度模型的反演.

图5展示了部分周期的二维瑞雷面波相速度分布，图6展示了各周期的瑞雷面波相速度对于 S波速度的敏感核曲线.结果显示，在短周期（6—12 s）相速度分布主要与地表地质构造相关，在河套裂谷、山西裂谷、渭河裂谷以及华北盆地和鄂尔多斯盆地局部表现出低速特征，而阴山造山带、吕梁造山带、太行造山带表现出高速特征；20 s周期的相速度主要对应中下地壳结构，较为明显的速度异常是在大同火山群区域的低速，向西与河套裂谷下方的低速相连，向南扩散至山西裂谷北段；40 s周期的相速度主要与上地幔顶部的速度结构相关，大同火山群区域表现出更加显著的低速特征，低速范围也有所扩大，低速的中心偏移至大同火山群西侧，另外太行山与华北盆地的交界线以及太行山南端也出现了带状低速；60 s以上周期的相速度主要与上地幔结构相关，表现出较为相似的速度分布，而大同火山群区域低速中心随深度增加延伸至华北盆地下方，可能显示了大同火山群岩浆活动的深部来源.

图5 华北克拉通中部不同周期的瑞雷面波相速度分布图（Ave表示每个周期T对应的平均相速度）Fig. 5 Maps of Rayleigh surface wave phase velocity at different periods in the central North China Craton The average phase velocity Ave for each period T is shown in the lower right corner

图6 一维参考模型（a）和不同周期的瑞雷面波相速度对于S波速度的敏感核曲线（b）Fig.6 1-D reference model （a）and sensitive kernel curves of the Rayleigh wave phase velocity at different periods relative to the S wave velocity（b）

1.3 一维S波速度联合反演

本文使用Li等（2019）发展的基于马尔科夫链蒙特卡罗（Markov chain Monte Carlo，缩写为MCMC）随机采样的非线性反演方法对背景噪声和地震面波计算的瑞雷面波相速度进行联合反演，获取每个网格点的一维S波速度模型，并据此整合成研究区域内0.25°×0.25°网格的三维S波速度模型.MCMC反演方法主要通过抽取收敛于贝叶斯后验概率分布的随机样本，对抽取的样本进行统计，从而获得后验分布的一些性质.近些年该方法在地球物理反演中得到了广泛应用，其原理也得到了详细的阐述（Bodinet al，2012；Afonsoet al，2013；Shenet al，2013），在此不作赘述.

对于每个网格点，利用背景噪声和地震事件得到的相速度频散联合反演S波速度.我们将一维初始速度模型设置为三层，包括低速沉积层，地壳和上地幔，每一层的厚度都随S波速度一起反演.其中沉积层由4个参数约束，包括顶部和底部的S波速度、沉积层厚度和沉积层中的波速比；地壳和上地幔的S波速度各由4个三次B样条参数约束，地壳厚度由莫霍面深度约束.因此总共设置了可在固定范围内扰动的13个独立参数用于反演一维S波速度.本文参考Crust1.0全球模型（Laskeet al，2013）来设置沉积层的初始厚度扰动范围.莫霍面的初始深度来源于Wang等（2017）基于接收函数的结果，并允许其在±5 km范围内扰动.本文采用PREM模型（Dziewonski，Anderson，1981）作为地壳和地幔的初始S波速度，其扰动范围设置为±20%.参照Li等（2019）的方法，利用马尔可夫链蒙特卡罗抽样算法在模型空间搜索采样并求解反演问题，对于每个网格点，进行十万次优化采样，利用马尔可夫链中最终被接受的4 000个模型来估计后验概率密度分布，并采用其均值作为每个网格点最终的一维S波速度.图7展示了其中三个网格点的一维S波速度模型反演结果，其位置以蓝色圆点标注在图1中.三个网格点下方的一维S波速度模型表现出不同的特征，其中：网格点1（38°N，109.5°E）（图7a，7d）位于鄂尔多斯盆地，在上地幔整体表现出高速特征；网格点2（40°N，112°E）（图7b，7e）位于大同火山群西侧，表现出相对较深的莫霍面和较低的上地幔速度；网格点3 （37.5°N，115°E）（图7c，7f）位于华北盆地，表现出较薄的地壳和相对低速的岩石圈结构.

图7 一维S波速度模型反演举例，其位置在图1中以蓝色圆点表示（a−c）瑞雷波相速度频散；（d−f）一维S波速度模型.黑色圆点代表观测的频散曲线，灰色线条代表计算的频散曲线，蓝色曲线代表最终的平均S波速度模型Fig.7 Examples of 1-D S-wave velocity model inversion. Their positions are represented by blue dots in Fig. 1（a−c）Rayleigh wave phase velocity dispersions；（d−f）1-D S wave velocity models. The black dots represent the observed dispersion curves，and the gray lines represent the calculated dispersion curves，and the blue curve represents the final average S-wave velocity model

2 结果与分析

计算得到每个网格点的一维S波速度模型后，将其整合成研究区域0.25°×0.25°网格的三维S波速度模型，图8展示了不同深度的S波速度异常结果，各个深度的平均S波速度标注在右下角.可以观察到，在不同深度S波速度都表现出比较明显的横向不均匀性.

图8 华北克拉通中部不同深度的S波速度水平切片（图中Ave表示每个深度h对应的S波平均速度）Fig.8 Horizontal slices of S-wave velocities at different depths in the central North China Craton. The average S wave velocity Ave of each depth is shown in the lower right corner（a）h＝5 km；（b）h＝15 km；（c）h＝30 km；（d）h＝50 km；（e）h＝70 km；（f）h＝100 km；（g）h＝140 km；（h）h＝200 km

在地壳深度范围内，5 km的S波速度切片（图8a）展示了上地壳速度结构，主要与地表地质构造有关.在该深度处，研究区域内的裂谷带和盆地都表现为低速，包括山西裂谷，渭河裂谷，河套裂谷和华北盆地，太原盆地以及鄂尔多斯盆地等，推测是与这些区域浅部覆盖的新生代沉积层有关.而山脉区域的S波都表现为高速特征，包括阴山山脉，吕梁山脉和太行山脉，与这些区域结晶基底的出露相对应.本文模型与地表地质构造表现出很好的相关性，证明了结果的可靠性.

15 km深度处的S波速度切片（图8b）大致对应中地壳的速度结构，其中，一个较为显著的低速区主要集中于大同火山群附近，向西延伸至河套裂谷东部，向北扩散至阴山—燕山造山带，覆盖了整个山西裂谷的北段，主要与位于山西裂谷北端的大同火山群的岩浆活动有关.另外，华北盆地也表现为低速异常，可能与新生代以来地壳活化有关（Tanget al，2013）.渭河裂谷表现为低速异常，反映了裂谷作用导致的较弱的地壳.鄂尔多斯盆地和吕梁、太行造山带整体表现为高速，反映了这些区域较为稳定的地壳结构.

30 km深度处的S波速度切片（图8c）主要反映了下地壳结构，相较于中上地壳，大同火山群附近的低速异常明显增强，且范围增大，低速异常蔓延至整个阴山—燕山造山带，可能代表了在下地壳中岩浆的侵入范围，大同火山群附近区域的下地壳被上涌的地幔热物质侵入加热，表现出大范围的显著低速区域.华北盆地，渭河裂谷，太原盆地转变为高速异常，可能表示这些区域在该深度已经接近或到达上地幔顶部，代表这些区域有着较薄的地壳厚度.鄂尔多斯盆地整体依旧表现为高速，证明其下地壳未遭受破坏，依然保持稳定.

50 km深度处的S波速度切片（图8d）主要反映了上地幔顶部的结构，大同火山群附近的低速异常相较于地壳中有所增强，但范围缩小，低速区的中心向西偏移至鄂尔多斯盆地的东北角以下，与阴山、河套裂谷处的低速异常相连，可能预示了地幔热物质在上地幔顶部的上涌通道.太行造山带和华北盆地的交界线，对应了南北重力梯度线，在沿线附近表现出南北走向的低速带，且与太行山南端的低速异常相连接.鄂尔多斯盆地下方依然表现为稳定的高速，与山西裂谷南段的高速相连.

70，100，140，200 km深度处的S波速度切片（图8e−h）表现出大致相同的分布特征，反映了上地幔中的速度结构.可以观察到，随着深度的增加，低速区的分布位置基本保持不变，即主要集中于大同火山群附近、沿南北重力梯度带以及中部造山带南端，而低速异常的幅度和范围逐渐减小.值得注意的是，在200 km深度，大同火山群东南方向的华北盆地表现出相对较强的低速异常，可能显示了大同火山群下方岩浆活动的深部来源.另外，鄂尔多斯盆地和山西裂谷南段下方则继续表现出较强的高速异常.

3 讨论

3.1 地壳结构

因为本研究利用了背景噪声成像得到的短周期面波速度进行联合反演，所以本文S波速度模型在浅部也有着较高的分辨率.上地壳的速度异常分布主要与地表的构造和地形有关，华北盆地、河套裂谷、渭河裂谷、太原盆地以及鄂尔多斯盆地北缘等区域的低速异常可能与浅部的新生代沉积层有关（Zhanget al，1998）.其中，华北盆地、河套裂谷的低速异常尤为显著，可能代表这些区域的沉积层相对较厚.上地壳中，山西裂谷整体表现为低速特征，NNE−SSW向的低速带与山西裂谷的位置走向相对应，证明我们的模型与地表构造比较吻合.李自红等（2014）基于地震反射剖面结果认为山西裂谷南段的临汾盆地是一个受拉张的纯剪切盆地，而Yin等（2017）基于大地电磁剖面的研究结果则支持简单剪切模型.在本研究中，从东西向的剖面CC′和DD′ （图9c，9d）可以观察到，在中下地壳，山西裂谷南段的低速异常随深度增加逐渐向东倾斜至太行造山带下方，东向倾斜的低速带可能显示了裂谷在地壳内的破裂面.因此，我们的结果为山西裂谷南段属于简单剪切变形模式（Yinet al，2017）提供了证据.经过阴山造山带和鄂尔多斯盆地的南北走向剖面EE′ （图9e）显示，河套裂谷下地壳附近存在显著的低速，并向北倾斜延伸至上地幔中，可能代表了深部上涌的地幔热物质侵入了下地壳，造成河套裂谷下方大范围的低速，并导致了岩石圈的改造和减薄，这与陈凌等（2010a）研究观察到的该区域较薄的岩石圈厚度相符合.

图9 S波速度模型垂直剖面，其位置以黑色虚线标注在图8h中，图中黑色粗实线代表莫霍面深度Fig. 9 Vertical sections of the S-wave velocity model，the positions of which are marked with black dotted lines in Fig. 8h，the thick black lines denote the Moho depth

3.2 岩石圈结构

本文的岩石圈模型中，最为显著的S波速度特征就是在上地幔中华北克拉通中部造山带南北两端的大范围低速异常区域以及连接两者的太行造山带东侧的低速异常带.其中，北段的低速异常一般认为与大同火山群的岩浆活动有关，地幔上涌的热物质形成了第四纪火山的喷发（Tanget al，2013；Liet al，2018）.Xu等（2005）基于第四纪大同火山群碱性玄武岩成分分析表明该处岩石圈发生了减薄，这与本文岩石圈模型在该区域下地壳及上地幔表现为大尺度的低速异常相对应，证明大同火山群区域的岩石圈已经被上涌的岩浆侵蚀发生了改造和减薄.同时本文的岩石圈速度模型显示该低速异常区域扩散至阴山造山带以及河套裂谷下方，可能代表这些区域的岩石圈也遭受到了不同程度的改造和减薄（Huanget al，2009；Tanget al，2013；Donget al，2014）.

中部造山带南端的低速异常代表其南部可能也遭受了岩石圈改造和减薄，但本文的模型显示其低速异常相对北段较弱（图8e−g），可能代表了南部的改造程度相对较低，岩石圈相对较厚，这与Tang等（2013）的研究结果一致.另外，该低速区在前人的地震成像研究中也被观察到（Baoet al，2013；Jianget al，2013；Guo，Chen，2017），但在细节上存在差异.Guo和Chen （2017）基于背景噪声和接收函数得到的模型显示该低速体在约70 km和大于100 km深度侵入了鄂尔多斯盆地下方，提出鄂尔多斯盆地东南部下方的岩石圈地幔存在局部改造；而Bao等（2013）基于背景噪声成像的结果在鄂尔多斯盆地东南部未发现显著的低速区域，认为鄂尔多斯盆地整体保持稳定.在本文的高分辨率模型中，鄂尔多斯盆地东南部下方的岩石圈地幔50—100 km深度处观察到局部的低速体（图8d−f），可能代表了局部的岩石圈改造，因此本文的结果支持鄂尔多斯盆地东南部岩石圈存在局部改造的结论.

另外，沿着太行造山带与华北盆地交界线的大范围低速带与南北重力梯度线的位置相对应，可能代表着华北克拉通中部造山带的岩石圈东侧已经被改造，这与梯度线以东较薄的地壳和岩石圈相对应（陈凌等，2010b）.由西向东的S波速度递减可能代表了中部造山带的岩石圈厚度由西向东逐渐减薄，这与基于接收函数的结果一致（Chenet al，2009），这种由西向东逐渐增强的岩石圈破坏可能与太平洋板块的西向俯冲有关（Huang，Zhao，2006）.值得注意的是，中部造山带的中南段在上地幔中（50—140 km）表现出与鄂尔多斯盆地东南部相连的高速异常（图8d−g），表明该区域的岩石圈可能未遭受严重破坏，整体保持相对稳定.地质年代学研究表明，华北克拉通东、西地块在约1.85 Ga前发生碰撞并最终导致华北克拉通的融合，同时形成了中部造山带（Zhaoet al，2005）.在本文的岩石圈模型中，中部造山带岩石圈不同位置的改造差异可能与华北克拉通复杂的构造演化活动和地幔深部的流体运动以及这些区域早期本身的不均一性有关（朱日祥等，2011）.本文的模型还显示出，在上地幔中，华北克拉通中部造山带南北两端的低速异常通过沿南北重力梯度线分布的低速带相连，且延伸至华北盆地下方，可能代表造成中部造山带南北两端的岩石圈改造减薄的机制相同，即华北盆地下方的地幔上涌的热物质侵入太行造山带下的岩石圈地幔并流向中部造山带南北两端，造成局部的岩石圈的改造减薄作用，这与朱日祥和郑天愉（2009）提出的该区域的岩石圈破坏以岩浆底侵或热侵蚀作用为主相符合.Huang和Zhao（2006）基于体波成像的结果显示，西向俯冲的太平洋板块前缘停滞在南北重力梯度线东侧，与本研究观察到的上述低速带位置相对应.而前人的研究表明，停滞的太平洋俯冲板块的脱水会导致上地幔软流圈热物质上涌，从而造成S波低速异常（Fukaoet al，1992；Zhao，2004；Lei，Zhao，2006）.因此，我们推测中部造山带上地幔中的低速体可能来源于停滞的太平洋俯冲板块脱水作用，也就是说，中部造山带岩石圈的改造减薄主要与太平洋板块的西向俯冲有关，这与朱日祥等（2012）提出的太平洋板块俯冲是华北克拉通破坏的主要动力因素相符合.

鄂尔多斯盆地下方除东南缘局部区域外，整体表现为持续的高速特征直至200 km以下（图8，9），证明其岩石圈厚度可能达到了200 km以上，整体保持了克拉通的属性.这与前人的面波成像结果显示鄂尔多斯盆地下方较薄的岩石圈存在差异（Huanget al，2009；Tanget al，2013），与基于接收函数和体波成像的研究结果相符合（Huang，Zhao，2006；Chenet al，2009）.因此本文的高分辨率S波速度模型倾向于支持位于华北克拉通腹地的鄂尔多斯盆地岩石圈整体未受到破坏减薄，保留了克拉通特性，与陈凌等（2010a）的研究结果相符.

3.3 大同火山群岩浆活动

在本研究得到的岩石圈速度模型中，大同火山群附近区域的显著低速特征由中下地壳一直延伸至上地幔中，其在深度上连续的分布可能显示了上地幔热物质由下而上的上涌通道和侵入范围.图8的水平切片显示，随着深度增加，大同火山群下方的低速区域在下地壳与上地幔顶部呈现向西偏移的特征，前人也观察到类似的特征（Tanget al，2013；Liet al，2018；Aiet al，2019）.在50 km深度，低速异常最强且范围集中于鄂尔多斯盆地东北角下方，可能显示了地幔热物质在上地幔顶部的侵入位置.基于接收函数的研究观察到该区域较深的莫霍面可能与岩浆的板底作用导致的地壳增厚有关（Xuet al，2005；Wanget al，2017）.岩浆从该位置侵入地壳且向周围扩散，将该区域附近的下地壳加热甚至可能导致部分熔融，造成了鄂尔多斯盆地东北角下方大范围的低速区域.而在地壳内继续上涌的过程中受到鄂尔多斯盆地东北角的上地壳阻挡，发生东向流动，随后上升至大同火山群区域正下方，形成了大同火山群的岩浆活动.

在上地幔中，该低速体的范围随深度的增加而逐渐扩大，但中心位置不再向西偏移，而是垂直向下延伸.在200 km深度左右大同火山群下方的低速异常逐渐减弱，而在其东南方向的华北盆地观测到相对较强的低速异常（图8h）.图9a中沿东西向的AA′剖面能够更加直观地观察到该低速体在垂向上的分布变化，剖面显示该低速体位于大同火山群西侧的鄂尔多斯盆地东北角下方，沿垂向一直延伸至约170 km深度，然而在更深处低速异常却逐渐减弱，可能表明该低速体不是来源于垂直向下的更深处，而是来源于其它方向的水平流动.另一条沿南北走向穿过大同火山群的FF′剖面显示，随着深度的增加，大同火山群下方的低速体逐渐向南延伸至约38°N—39°N之间，连续的低速异常到达200 km以下（图9f）.而穿过该低速区域的东西走向的剖面BB′显示，随深度的增加，该低速体向东延伸到华北盆地200 km以下（图9b），可能显示了大同火山群岩浆的深部来源.

综合以上分析，我们推测由华北盆地下方地幔中上升的热物质沿西北方向上涌至中部造山带北段西侧的鄂尔多斯盆地东北角下方，在该处侵入地壳后在继续上升过程中受到该区域上地壳阻挡，发生东向流动，最终到达大同火山群下方，为大同火山群的岩浆活动提供热物质来源.因此，本文的研究结果不支持Li等（2018）提出的大同火山群岩浆来源于青藏高原东北缘的地幔流，而更倾向于支持大同火山群岩浆来源于华北盆地下方的上地幔.考虑到目前没有在大同火山群区域发现明确的地球化学方面的证据支持该区域在新生代存在地幔柱活动（Menzieset al，2007），而深部的俯冲板块脱水作用会造成大规模的软流圈热物质上涌（Fukaoet al，1992；Zhao，2004；Huang，Zhao，2006），从而为大同火山岩浆活动提供热物质来源.因此，基于本研究所得到的S波速度模型，推测大同火山群的岩浆活动深部来源可能与西向俯冲的太平洋停滞板块有关.

4 结论

基于ChinArray三期一阶段密集台阵观测数据，利用背景噪声和地震面波层析成像联合反演了华北克拉通中部造山带及其邻近区域地幔岩石圈的高分辨率三维S波速度模型.本文的速度模型为研究华北克拉通破坏减薄机制和横向的破坏改造程度差异以及大同火山群岩浆活动在地壳内的上涌通道和深部来源提供了新的证据.结果显示，华北克拉通不同地块的岩石圈结构存在很强的横向变化，鄂尔多斯盆地腹地下方除东南缘存在小范围的局部改造外，整体保留了稳定的克拉通属性，北缘的河套裂谷区域存在自下而上的岩石圈改造减薄作用.东部的华北盆地岩石圈被强烈破坏减薄，显示出较薄的低速地壳；中部造山带南北两端以及沿南北重力梯度线的岩石圈也被不同程度地改造减薄，其中北端大同火山群下方的岩石圈改造最为强烈，南端的改造程度相对较弱；而南北两端的低速体通过南北重力梯度线下方的低速带相连并延伸至华北盆地地幔深处，可能代表了相同的改造减薄机制，即华北盆地下方地幔上涌的热物质侵入太行造山带下的岩石圈地幔并流向中部造山带南北两端，造成岩石圈局部的改造减薄作用，地幔热物质的深部来源可能与太平洋板块的西向俯冲有关；大同火山群的岩浆活动可能来源于华北盆地下方的地幔热物质沿西北方向上涌，在鄂尔多斯东北角下方侵入地壳，在地壳内上升过程中受到上地壳阻挡发生东向流动，到达大同火山群下方.因此，本文的结果倾向于支持大同火山群的岩浆活动在深部可能来源于西向俯冲的太平洋停滞板块的脱水作用.

中国地震局地球物理研究所“中国地震科学探测台阵数据中心”为本研究提供了地震波形数据，作者在此表示感谢.