地球物理学报  2012, Vol. 55 Issue (5): 1613-1623   PDF    
鄂尔多斯地区上地幔岩石圈三维速度结构面波反演研究
李多, 周仕勇 , 陈永顺, 冯永革, 李鹏     
北京大学地球与空间科学学院, 北京 100871
摘要: 双平面波拟合法是一种新的面波成像方法, 反演中考虑地震波场中的非平面波成分, 提高反演的分辨率.本文利用双平面波拟合法, 反演获得鄂尔多斯地区上地幔岩石圈的速度结构.所用资料为国家数字地震台网69个宽频带地震仪和北京大学34个流动数字地震台观测到的地震波面波资料.首先从面波记录中提取了研究区域20~125 s瑞利波相速度频散曲线, 进而得到各个周期瑞利波相速度异常分布图.结果显示, 短周期瑞利波相速度异常与地表的构造特征吻合较好, 中长周期的瑞利波相速度可以反映出上地幔岩石圈的速度异常分布以及构造特征.由研究区20~125 s的瑞利波相速度分布图可以反演得到地表到地下200 km范围内的三维剪切波速度结构.结果显示, 鄂尔多斯块体内部稳定均一, 活化或改造的痕迹不明显;鄂尔多斯块体西南缘受到青藏高原的强烈作用, 有大量地幔物质流动的痕迹存在;中央转换带下超过200 km深度存在地幔物质上涌, 可能与太平洋板块的俯冲和青藏高原板块的挤压有关.
关键词: 双平面波拟合法      面波成像      鄂尔多斯      华北克拉通      岩石圈     
3-D lithospheric structure of upper mantle beneath Ordos region from Rayleigh-wave tomography
LI Duo, ZHOU Shi-Yong, CHEN Yong-Shun, FENG Yong-Ge, LI Peng     
Institute of Theoretical and Applied Geophysics, School of Earth and Space Science, Peking University, Beijing 100871, China
Abstract: We apply two-plane-wave tomography, which takes the influences of the non-plane wavefield into consideration, to study the structure of velocity of the upper mantle lithosphere of Ordos region. The source of data consists of two parts, one is from 69 broad-band seismographs of China Earthquake Networks Center, and the other is from 34 mobile digital broad-band seismographs of PKU. At short periods most high and low velocity anomalies correlate well with surface geological features. The difference of the structure of upper mantle lithosphere is revealed by long-period surface waves. We extracted the 3-D structure of shear velocity anomalies of 200 km depth from the Rayleigh wave phase velocities. The results reveal that the Ordos block has the high velocity beyond 200 km depth, and no evidence for reactivation. There is fierce interaction between Ordos block and Tibet block at the southwestern edge of Ordos block, which causes the upper mantle flow there. Upwelling beneath the Central zone may be caused by both the subduction of Pacific plate and India plate.
Key words: Two-plane-wave method      Surface wave tomography      Ordos      North China Craton      Lithosphere     
1 引言

鄂尔多斯块体属于华北克拉通的一部分. 华北克拉通可分为东部块体、西部块体和中间转换带三部分,其中西部块体包括阴山块体和鄂尔多斯块体. 地质学研究表明,华北克拉通的东、西块体在早元古代至中元古代相互分离,各自经历了不同的发育历史. 根据鄂尔多斯块体北缘的孔慈岩证据显示,鄂尔多斯块体在约2.0~1.9Ga与北边的阴山板块碰撞拼合,在约1.85Ga与东部板块碰撞拼合,并形成了现在的中央转换带[1-3]. 晚中生代到新生代时期,伴随着岩石圈减薄、断陷盆地形成、地下热流、火山活动等,发生了华北克拉通的活化[4-7].

最近十几年来,中国东部岩石圈的减薄和华北克拉通的活化成了国内外学者的热门研究问题,相关的研究成果层出不穷. 地震波层析成像是用于研究地球内部结构的最有效手段,在地球深部结构探测研究中具有重要作用. 已有很多地球物理学家利用宽频带地震台网记录,运用各种地震波层析成像技术来探测鄂尔多斯地区的地壳和上地幔结构. 已有的瑞利面波成像和P 波成像得到了华北地区的三维地壳和上地幔结构,结果显示鄂尔多斯块体整体坚固稳定[8-10]. 接收函数的结果显示出华北克拉通西部板块和东部板块的壳-幔过渡带存在明显的结构差异[11-12]. 横波分裂的结果显示山西断陷带上地幔的各向异性方向与地表断陷带的走向一致,这与拉张区地幔物质流动有关[13]. 利用P波和S波层析成像技术得到的华北克拉通速度异常分布结果显示,中央转换带下超过500km 深的低速异常是由于地幔物质上涌引起的,这与华北克拉通东部板块的活化密切相关[14].

利用面波相速度的频散特征反演地下速度结构,是地壳及上地幔结构探测的主要地震学方法之一. 目前,单台法和双台法仍然是国内地震学家提取面波频散特性的主要方法[15]. 单台法测量的是源台间面波的频散特性,在较近的台站上计算得到的面波相速度会存在较大的误差,因此不可能用于区域小尺度速度结构精细反演. 而双台法则需要两个接收台站和震源位于同一大圆路径上,利用满足此条件的双台地震资料可以求得双台间的面波频散特性,然而两个接收台站和震源严格位于同一大圆路径上的记录并不很多,使双台法的应用同样受到局限.

传统的单台法及双台法提取面波频散特性的理论基础也面临挑战. Wielandt和Friederich 自1993 年来[16-19]发表的系列理论研究表明单台法及双台法提取面波频散特性的理论基础所隐含的平面波假设在实际中很难成立. Wielandt[19]研究指出在地震面波的非平面波成分不可忽视的时候,传统面波相速度测量方法所得到的只是波场的动态相速度. 动态相速度除了包含传播介质的物理性质信息之外,还包含了面波波场振幅局部不均匀分布的信息;因此不能直接反映研究区域介质的物理性质. 只有正确估计面波波场局部的结构,才可能得到真正与介质相关的面波结构相速度. 为解决这一问题我们拟应用Forsyth和Li[20]提出的一种通过用双平面面波拟合地震台阵记录的非均匀面波波场,提取台阵覆盖区介质的面波相速度频散曲线的方法,这样因为无平面波假设,减少了“结构性误差”,理论更为严谨. 而且面波波场拟合法同时利用了台站记录的面波振幅和相位信息,且不受大圆路径的局限,利用的资料更丰富.

Jiang等[21]曾经利用双平面面波拟合法开展了藏南地区面波精细成像并发现藏南底下存在约40km 厚的软弱下地壳层,为青藏高原下地壳流模型提供了支持. 本文作者李鹏曾经利用双平面面波拟合法开展了鄂尔多斯块体和中央转换带地区地下结构的面波反演研究[22],探测到稳定的鄂尔多斯块体岩石圈可能延伸至底下120~140km 深度,而山西断陷区及华北盆地的岩石圈可能减薄至70~80km 深度,给出了华北克拉通现代处于活化过程中及其可能动力学图像的面波探测证据,同时也验证了双平面面波拟合法的有效性. 但因为当时没有收集到银川、乌海等地鄂尔多斯西北缘的几个中国国家台网台站的记录(图 1) ,不能对西缘—青藏板块与鄂尔多斯块体交汇区底下结构有清晰成像,并影响研究区成像的整体分辨率. 本研究我们将在李鹏等[22]的工作的基础上,补充银川台等23个鄂尔多斯西北缘的中国国家台网台站的记录,用双平面波拟合法反演开展鄂尔多斯块体和中央转换带地区瑞利波相速度二维分布图像的反演,期望得到研究区更高分辨率的面波成像结果,为认识华北克拉通的形成和演化提供更可靠的地震学证据.

图 1 所用台站分布图及研究区域构造地质构造图.红色三角为国家地震局固定台站,绿色三角为北京大学流动台站.白色方框为研究区域.台站主要分布在山西断陷带和渭河断陷带附近 Fig. 1 Distribution of the stations and tomography and tectonics of the study region. Red triangles denote the stations of CENC; Green triangles denote the stations of PKU. White pane denotes the study region. Stations distribute densely along Shanxi Ritt and WeiheRitt
2 数据处理

本文的数据来源为两部分:中国地震局台网中心在33°N—42°N 和105°E—115°E 的区域范围内的69个固定台站的地震波形数据[23]. 记录时间为2007年8月—2009年11 月;北京大学在山西断陷带布设的34个二维流动数字台网的地震波形数据,记录时间为2007年1—6月. 在数据选取时,只选取了瑞利面波的垂直向记录,同时根据信噪比,筛选了196个震级≥6.0级、震中距为30°~120°、方位角分布较为均匀的远震事件(图 2) .

图 2 所用地震事件后方位角分布图.五角星为 研究区域中心,实心圆点为所用地震事件 Fig. 2 Azimuthal distribution of the teleseismic events used in the study. Black star denotes the center of the study region. Black circles denote the events used in the study

处理数据过程中,将所有台站的波形信号统一校正为CNG-3ESPC的拾震器的仪器响应. 选取带宽为10mHz的窄带滤波器(巴特沃斯型滤波器,级数4,通道数2) . 从0.008Hz到0.05Hz选取了12 个中心频率点,对地震波进行滤波,从而获取了12 个窄频段的波形数据. 对于每个地震事件,反演方程中包含6个波场参数(每个平面波需要用振幅、初相位与传播方位角3个参数描述)和2 个格点的相速度参数,因此,我们保留射线记录数大于9的地震事件,共得到约27000条地震射线,各频段(周期)的记录数条形统计图如图 3. 图 4 显示了50s周期的地震射线分布,可以看出射线覆盖较好的地区为南部和东部,北部和西部由于受到台站分布和地震事件分布的局限,射线分布较为稀疏.

图 3 所用25〜125 s周期瑞利面波记录统计图 Fig. 3 Numbers of the Rayleigh waveforms used at 12 periods of 25〜125 s in the inversion
图 4 周期为50 s的地震射线覆盖分布图. 研究区域东部和南部射线覆盖密集 Fig. 4 Great circle ray path coverage at the period of 50 s. The density on east and south is good
3 方法概述 3.1 动态波数与结构波数的区别

考虑一列沿地球表面传播的地震波(面波),其子波频率为ω 的波场函数可以写为

(1)

(x,y)地表上某观测点坐标,波数矢量Ko = (kx,ky),ko =|Ko|,令

(2)

则有

(3)

由此可得到(x,y)点上频率为ω 的相速度:

(4)

注意到(3) ,(4) 式中的波数和相速度是直接由测量面波波场在空间的分布而得到,我们称之为波场动态波数|相速度(Dynamic wave wavenumber|phase velocity).

由动力学方程出发,我们可以推导出(1) 式表达的波场函数满足Helmholtz方程,即有

(5)

(5) 式中kω/c. 由于(5) 式直接由动力学方程推导,kc由波的传播介质的弹性参数唯一确定,我们将k,c称为结构波数和结构相速度(structural wavenumber|phase velocity). 由(5) 式有

(6)

显然(6) 式与(3) 式所表达的波数在物理上并不完全一样,数学表达上也存在差别.

由于波场函数必须满足Helmholtz方程,我们将(2) 式代入(6) 式有

(7)

将(3) 式代入(7) 则有

(8)

由(8) 式可清楚地看出:当面波不是平面波时,即台阵记录的同一面波因记录台位置不同,其振幅不一样,(1) 式中的A亦即(8) 式中a在空间上不是常量,则观测的波场动态波数k0 与结构波数k是存在系统偏差的. 传统的台源法和双台法测量的面波相速度频散特性只是观测的面波波场动态相速度频散特性,是一种平面波假设下的近似测量.

3.2 双平面波拟合法

为得到精确的面波反演结果,本研究拟采用Forsyth与Li[20]方法原理,其基本原理概述如下:对频率为ω 沿水平方向传播的非平面波可以用两组同频率、不同水平方向传播、具有不同相位和振幅的平面波的干涉效应表达,即非平面波垂直向位移波场可以表达为[20]

(9)

ikU是第k个台站记录的第i个事件垂直方向的位移;iA1iA2 表示两平面波的振幅,为待定参量;i0φ1 和i0φ2 是两平面面波在参考台站的初相位. kiτi0τ是波从研究区域的边界到待求台站和参考台站沿大圆路径的传播时间,iθ1 和iθ2 分别为两平面波入射方向与大圆弧路径的夹角,kiS 为平均慢度. 这样对于频率为ω 的地震面波波场只有6个参数就可以模拟:两列平面波的振幅、入射方向和在参考台站的相位(图 5) .

图 5 研究区域建立的区域直角坐标系.参考台站位于 坐标原点,(r,φ)是与直角坐标系(x,y)相对应的极坐标 系坐标,横轴为射线大圆弧路径方向,平面波以与横轴 夹角为彡的方向入射.虚线方框为研究区域,相对到时零点(τ=0)为沿大圆弧路径传播的面波波前最先接触 到研究区域的时刻[24] Fig. 5 Local coordinate system of two-plane-wave method. Reference station is at (0,0),(r,φ)is the position of polar coordinate as (x,y)in rectangular coordinate.x-coordinate is along the great circle path,while y axis is perpendicular to x axis. The incoming plane wave s along the direction of propagation 6. The dashed square frame denotes the study region. The time that the wavefront along the great circle path arrives at the first point of the study region is 0 (τ=0)[24]

我们将研究区分成N个等尺度网格,一般情形下,介质速度具有各向异性,因此不同地震的面波通过j网格,由于射线穿过的方位不同,速度不同. 因此i地震面波在j网络传播的速度为

(10)

我们将研究区域划分为0.5°×0.5°的网格,忽略每个格点的各向异性成分,即令B1,B2 为0,每个格点的瑞利波相速度为均一值. 对每个格点加入高斯平滑因子,则每个格点的慢度表示为

(11)

其中,i表示地震编号,j表示格点编号,高斯权重函数为

(12)

其中,特征长度Lw 直接影响高斯权重函数的形状,可以控制反演得到的速度结构的光滑程度,Lw 过小会导致空间小尺度内出现突变;过大则会使得合理的速度差异不够显著. 所以,在反演过程中需要调整特征长度来调节结果的分辨率和标准差范围[20]. 经过不同的尝试,以获取研究区最佳分辨率为标准,我们将Lw取为50km.

将6个波场参数和每个格点的慢度参数带入反演,分别利用模拟退火法[25]和广义最小二乘法反演[26],可以得到每个格点的瑞利波相速度.

4 反演结果与分析 4.1 检测板试验

检测板是验证结果的分辨率和可靠性的重要依据[27]. 在本文的研究中,我们将输入的速度模型的异常范围设定为±5% 之间,异常尺度设置为2°×2°,即空间尺度约为200km×200km,这一尺度基本可以分辨出研究区域的构造单元的速度异常.

通过将反演结果与输入模型进行比对(图 6) ,可以看出,短周期时,仅有山西断陷带地区台站密集处可以得到反演结果,这是因为高频信号衰减快、信噪比低. 中长周期时,鄂尔多斯块体、山西断陷带和南部的渭河断陷带处结果较为清晰准确,而北部,尤其是西北部,由于射线分布稀疏,反演可信度不高. 由于本文的目的就是研究鄂尔多斯块体和山西断陷带的上地幔岩石圈结果,西北部和东北部的低分辨检测结果对我们后面的研究结论影响不大.

图 6 二维瑞利波相速度检测板试验.输入速度异常模型为一5%〜5%之间 Fig. 6 Checkerboard test for the 2-D Rayleigh wave isotropic phase velocity. The input velocity model iswith anomalieswithin —5 % 〜5 %
4.2 瑞利波相速度

首先将研究区域内各点瑞利波相速度设为相同,反演得到一维瑞利波相速度频散曲线(图 7) . 在短周期时,该地区相速度基本吻合AK135模型[28],而在中长周期时,该地区相速度低于AK135 模型[28],反映了华北克拉通的活化而引起的岩石圈整体的温度升高.

图 7 —维瑞利波相速度频散曲线.实线为AK135参 模型,长虚线为输入模型,短虚线黑线为反演结果 Fig. 7 Average Rayleigh-wave phase velocities in Ordos region (black circles with error bar) compared with those of global AK135 model (solid black line). The long dashed line is the input velocities for 12 different periods

以一维平均瑞利波速度为初始模型,反演得到二维瑞利波相速度异常分布图(图 8) . 图中红色为低速异常,蓝色为高速异常. 短周期时(20~29s)的瑞利波反映地表到莫霍面的速度异常分布,与地表的构造特征基本吻合,如太原盆地为新生代的断陷盆地,沉积层较厚,表现为瑞利波的低速异常,而太行山、吕梁山等为高速异常.

图 8 二维瑞利波相速度异常分布图 Fig. 8 Maps of Rayleigh wave isotropic phase velocity anomalies

中长周期(33~125s)反映了上地幔岩石圈的速度结构. 这里可以看到鄂尔多斯块体整体为明显的高速异常体,内部速度差异较小,与这里块体内部少地震的事实相吻合. 以山西断陷带为界,东西两侧的高低速异常分割明显. 低速异常沿中央转换带呈南北条带分布.

鄂尔多斯块体是中国大陆最古老的克拉通,自形成以来,没有经历过活化和改造,被认为是最坚固稳定的大陆岩石圈[29]. 晚太古代至元古代时期,华北克拉通的西部板块与东部板块碰撞拼合,拼合部位在中央转换带[30]. 晚中生代到新生代,华北克拉通再次活化,华北地区岩石圈拉张减薄,形成了一系列的断陷盆地[31-32]. 岩石圈的热结构研究[33]显示,山西断陷带存在较高的热值,对应这里的低速异常,很有可能与地幔物质流动有关.

4.3 剪切波速度结构

剪切波速度相对于物质的流变性变化敏感. 不同深度剪切波速度对于不同频率的瑞利面波相速度的敏感程度不同,大体上为面波波长的1/3 处为峰值,向两侧逐渐递减[34]. 并且越深处的剪切波对各频率相速度的敏感程度差别减弱. 本文所用研究方法是Weerarantne等[35]和Li等[27]在Saito[36]的方法基础上而发展起来的. 参考模型选择了AK135模型[28],初始地壳模型为三层(图 9) . 反演结果受到剪切波速度阻尼系数和地壳厚度阻尼系数的共同控制. 阻尼系数过大会导致误差较大,无法得到最优化解;阻尼系数过小则会出现约束不够强,导致解的不稳定. 所以阻尼系数的选择要权衡两方面的因素,反演过程中要根据结果调整参数,以得到较好的结果. 本文中根据多次反演结果的比较,设定剪切波速度阻尼系数和地壳厚度阻尼系数为0.05 和2. 另外,由于莫霍面深度对于反演结果影响较大,所以反演过程不限定莫霍面深度,这样保证了结果的稳定性和可靠性[37].

一维平均剪切波速度曲线图(图 9) 显示,该地区整体上地幔岩石圈速度低于AK135 全球平均模型[28]. 须指出的是,图 7 显示鄂尔多斯的相速度频散曲线在几乎整个地壳(50s以内)内与AK135 非常一致,但是图 9显示的反演得到的一维剪切波速度结构则显示鄂尔多斯内地壳速度明显低于AK135,这尽管可能与我们观测的频散曲线缺乏高频段信息(短周期截止点在20s),对地壳速度的约束较弱有关,但长周期(大于50s)部分低相速度,也可能造成对地壳速度的反演结果系统偏低. 鄂尔多斯内上地幔岩石圈速度明显偏低的事实是清晰的,这可能与晚中生代到新生代的华北克拉通的活化而引起的热物质的流动有关,具体分析将在后面详细说明.

图 9 一维剪切波速度结构. 实线为AK135参考模型,虚线为反演结果 Fig. 9 1-D average shear wave velocities in Ordos region(black dashed line). The solid line is the AK135 global model with three crustal layers. The thickness of the crust s 40 km from the shear wave inversion

二维剪切波速度差异分布图(图 10) 显示出了不同深度处剪切波速度结构. 浅层的剪切波速度与地表的地质构造特征吻合较好. 莫霍面两侧的速度差异反映出了地幔岩石圈的结构变化特征. 可以看出鄂尔多斯块体在莫霍面以下超过200km 深处均为明显的高速异常.

图 10 三维剪切波速度异常分布图.图中数字表示出所在深度范围.黑色实线表示纵向剖面图 Fig. 10 Maps of shear wave velocity anomaly. The depth range of each map has been marked. The thick black lines are vertical profiles in further study

鄂尔多斯块体西南缘出现低速异常,深度超过200km. 这里是青藏板块挤压处,隆起形成了六盘山. 青藏板块受印度板块的挤压作用,印度板块仍以每年6cm 的速度向北运动,由此挤压青藏板块,青藏板块遇到坚硬的鄂尔多斯板块,引起地壳和地幔物质沿边界向东溢出. 横波分裂结果显示这里的上地幔各向异性强烈,证明了地幔存在大量的物质溢出[38-39]. 接收函数的结果显示,青藏板块与鄂尔多斯板块的过渡带地壳变形强烈,说明这里强烈的挤压作用[40]. 我们的观测结果显示,这里的低速异常区深度约为200km,这预示着上地幔较大范围内可能存在热物质的流动.

中央转换带下的低速异常明显,深度超过200km,印证了中央转换带内新生代以来的活化导致了山西断陷带中的一系列断陷盆地的形成[41-42]这一结论. 横波分裂的结果[13-14, 43]显示,中央转换带的各向异性方向与断陷带的走向一致,即为北东-南西向. 而我们的结果给出了上地幔岩石圈的低速异常范围,很好地补充了前人的研究.

值得注意的是,中央转换带下的低速异常并不是连续的,在中段下,存在低速异常的间断. 这个间断在莫霍面以下的各个图像上均存在,超过200km 深. 这个间断面的形成可能与整个华北克拉通的活化与减薄的成因有关.

垂向的剪切波速度分布显示(图 11) ,鄂尔多斯地区整体的平均莫霍面深度为40km,山西断陷带和渭河断陷带的地壳出现较明显的抬升和减薄,这符合拉张区的构造特征. 接收函数结果显示,鄂尔多斯地区的莫霍面深度为40~43km,山西断陷带莫霍面比鄂尔多斯地块薄5km 以上[44-45]. 面波方法对莫霍面的深度约束并不强,只能一定程度上反映莫霍面的变化范围,我们的结果与接收函数的结果较为一致. 六盘山、太行山下的地壳厚度较大,符合造山带的特点. 太行山北段和华北平原的西缘下存在明显的低速异常带,这两处异常带都处在中央转换带上,华北克拉通的活化可能导致原本的板块拼合带出现了拉张减薄,引起地幔物质上涌.

图 11 剪切波速度垂直剖面.蓝线表示地表地势.黑线表示莫霍面位置.剖面位置在图 10中表示出.剖面AA'和BB'为经 向剖面,分别沿山西断陷带的东西两侧.剖面CC'和DD'为纬向剖面,其中CC'剖面经过渭河断陷带,DD'剖面横贯鄂尔多斯块体内部 Fig. 11 Vertical profiles of shear wave velocity. Blue lines show the topography. Black lines show the positions of Moho. Positions of these profiles are shown in Fig. 10. Profile AA'and BB; are along the western and eastern edges of Shanxi Rift. Profile CC' is along the Weihe Ritt,andDD'passes through Ordos block.
5 结论

我们利用双平面波拟合法,得到了鄂尔多斯地区的瑞利波相速度结构,进而反演得到剪切波三维速度结构. 从双平面波波场拟合方法的原理来看,该方法很大程度上还原了波场中的非平面波能量成分,使得反演结果更为准确. 在资料处理方面,抛弃了台源大圆弧路径这一条件,使得每一个地震事件在研究区域内的全部台站资料都可以得到有效利用,大大提高了反演的效率.

短周期瑞利波相速度异常分布图(图 8) 很好地反映了地表的地质构造特征,太原盆地、临汾盆地在新生代接受了大量沉积物,表现为低速异常,太行山、吕梁山为坚固稳定的块体,表现为高速异常. 中长周期瑞利波相速度反映出上地幔岩石圈的构造特征,为分析华北克拉通的地质构造背景提供了重要的依据.

瑞利波相速度结果和剪切波速度结果均显示鄂尔多斯块体具有超过200km 的克拉通根基,坚固稳定,没有活化和改造的迹象. 其西南缘上,青藏高原本身的地壳增厚和向东北方向的推挤,使得这里产生了沿着鄂尔多斯块体边界向东的地幔物质流动.

中部转换带下存在超过200km 的低速异常并不是呈连续条带分布,说明这里的地幔物质上涌并不是连续贯穿整个中央转换带的,这种现象并不是单一因素造成的,很可能要受到太平洋板块俯冲下插和青藏高原板块的俯冲推覆的双重作用. 而这两种作用各自的机制以及两者之间的关联,还需要进一步的观测资料才能解释.

致谢

感谢中国地震局地球物理研究所“国家数字测震台网数据备份中心”为本研究提供地震波形数据. 感谢参与北京大学二维流动台站布设的所有老师和同学. 感谢中国科学院地球物理研究所姜明明博士对本研究提出的意见和帮助.

参考文献
[1] Zhao G, Sun M, Wilde S A, et al. Late Archean to Paleoproterozoic evolution of the North China Craton: key issues revisited. Precambrian Research , 2005, 136(2): 177-202. DOI:10.1016/j.precamres.2004.10.002
[2] Wu C, Zhong C T. The Paleoproterozoic SW-NE collision model for the central North China Craton. Progress in Precambrian Research , 1998, 21(3): 28-50.
[3] Lu L, Xu X, Liu F. Early Precambrian Khondalite Series in North China. Changchun: Changchun Publishing House , 1996: 219-234.
[4] Griffin W L, Andi Z, O'Reilly S Y, et al. Phanerozoic evolution of the lithosphere beneath the Sino-Korean Craton. // Flower M F J, Chung S L, Lo C H, Lee T Y eds. Mantle Dynamics and Plate Interactions in East Asia. Washington: American Geophysical Union, 1998: 107-126.
[5] Menzies M A, Xu Y G. Geodynamics of the North China Craton. // Flower M F J, Chung S L, Lo C H, Lee T Y eds. Mantle Dynamics and Plate Interactions in East Asia. Washington: American Geophysical Union, 1998: 155-165.
[6] 吴福元, 葛文春, 孙德有, 等. 中国东部岩石圈减薄研究中的几个问题. 地学前缘 , 2003, 10(3): 51–60. Wu F Y, Ge W C, Sun D Y, et al. Discussions on the lithospheric thinning in Eastern China. Earth Science Frontier (in Chinese) (in Chinese) , 2003, 10(3): 51-60.
[7] Wu F, Xu Y, Gao S, et al. Lithospheric thinning and destruction of the North China Craton. Acta Petrologica Sinica , 2008, 24(6): 1145-1174.
[8] 陈国英, 宋仲和, 安昌强, 等. 华北地区三维地壳上地幔结构. 地球物理学报 , 1991, 34(2): 172–181. Chen G Y, Song Z H, An C Q, et al. Three-dimensional crust and upper mantle structure of the North China region. Chinese J. Geophys. (in Chinese) (in Chinese) , 1991, 34(2): 172-181.
[9] 易桂喜, 姚华建, 朱介寿, 等. 中国大陆及邻区Rayleigh面波相速度分布特征. 地球物理学报 , 2008, 51(2): 402–411. Yi G X, Yao H J, Zhu J S, et al. Rayleigh-wave phase velocity distribution in China continent and its adjacent regions. Chinese J. Geophys. (in Chinese) (in Chinese) , 2008, 51(2): 402-411.
[10] Li C, van der Hilst R D, Toks?z M N. Constraining P-wave velocity variations in the upper mantle beneath Southeast Asia. Physics of the Earth and Planetary Interiors , 2006, 154(2): 180-195. DOI:10.1016/j.pepi.2005.09.008
[11] Zheng T Y, Zhao L, Zhu R X. Insight into the geodynamics of cratonic reactivation from seismic analysis of the crust-mantle boundary. Geophys. Res. Lett. 2008, 35(8): L08303, doi: 10.1029/2008GL033439.
[12] Zheng T Y, Chen L, Zhao L, et al. Crust-mantle structure difference across the gravity gradient zone in North China Craton: Seismic image of the thinned continental crust. Physics of the Earth and Planetary Interiors , 2006, 159(1-2): 43-58. DOI:10.1016/j.pepi.2006.05.004
[13] Zhao L, Zheng T. Complex upper-mantle deformation beneath the North China Craton: implications for lithospheric thinning. Geophysical Journal International , 2007, 170(3): 1095-1099. DOI:10.1111/gji.2007.170.issue-3
[14] Zhao L, Allen R M, Zheng T, et al. Reactivation of an Archean craton: Constraints from P-and S-wave tomography in North China. Geophys Res Lett , 2009, 36(17): L17306. DOI:10.1029/2009GL039781
[15] 姚华建, 徐果明, 肖翔, 等. 基于图像分析的双台面波相速度频散曲线快速提取方法. 地震地磁观测与研究 , 2004, 25(1): 1–8. Yao H J, Xu G M, Xiao X, et al. A quick tracing method based on image analysis technique for the determination of dual stations phase velocities dispersion curve of surface wave. Seismol. Geomagn. Observ. Res. (in Chinese) (in Chinese) , 2004, 25(1): 1-8.
[16] Friederich W. Wave-theoretical inversion of teleseismic surface waves in a regional network: phase-velocity maps and a three-dimensionalupper-mantle shear-wave-velocity model for southern Germany. Geophysical Journal International , 1998, 132(1): 203-225.
[17] Friederich W, Wielandt E, Stange S. Multiple forward scattering of surface waves: Comparison with an exact solution and Born single-scattering methods. Geophysical Journal International , 1993, 112(2): 264-275. DOI:10.1111/gji.1993.112.issue-2
[18] Friederich W, Wielandt E. Interpretation of seismic surface-waves in regional networks-joint estimation of wave-field geometry and local phase-velocity-method and numerical tests. Geophysical Journal International , 1995, 120(3): 731-744. DOI:10.1111/j.1365-246X.1995.tb01849.x
[19] Wielandt E. Propagation and structural interpretation of non-plane waves. Geophysical Journal International , 1993, 113(1): 45-53. DOI:10.1111/gji.1993.113.issue-1
[20] Forsyth D W, Li A. Array analysis of two-dimensional variations in surface wave phase velocity and azimuthal anisotropy in the presence of multipathing interference. // Levander A, Nolet G eds. Seismic Earth: Array Analysis of Broad-band Seismograms. Washington: American Geophysical Union. 2005: 81-97.
[21] Jiang M M, Zhou S Y, et al. 3-D lithospheric structure beneath southern Tibet from Rayleigh-wave tomography with a 2-D seismic array. Geophysical Journal International , 2010, 185(2): 593-608.
[22] 李鹏, 周仕勇, 陈永顺, 等. 利用双平面波干涉面波层析成像方法研究山西断陷盆地及鄂尔多斯地台三维速度结构. CT理论与应用研究 , 2010, 19(3): 47–60. Li P, Zhou S Y, Chen Y S, et al. 3D velocity structure in Shanxi graben and Ordos from two plane waves method. CT Theory and Applications (in Chinese) (in Chinese) , 2010, 19(3): 47-60.
[23] 郑秀芬, 欧阳飚, 张东宁, 等. "国家数字测震台网数据备份中心"技术系统建设及其对汶川大地震研究的数据支撑. 地球物理学报 , 2009, 52(5): 1412–1417. Zheng X F, Ouyang B, Zhang D N, et al. Technical system construction of Data Backup Centre for China Seismograph Network and the data support to researches on the Wenchuan earthquake. Chinese J. Geophys. (in Chinese) (in Chinese) , 2009, 52(5): 1412-1417.
[24] 姜明明. 藏南岩石圈结构研究及地震活动性统计建模. 北京: 北京大学地球物理系, 2008. Jiang M M. Studies on the lithospheric stucture in souther Tibet and modeling the seismicity of north China by statistics (in Chinese). Beijing: Geophysics Department of Peking University, 2008.
[25] Press W H. Numerical Recipes in FORTRAN: the art of Scientific Computing. Cambridge: Cambridge Univ Pr, 1992.
[26] Tarantola A, Valette B. Generalized nonlinear inverse problems solved using the least squares criterion. Rev. Geophys. Space Phys. , 1982, 20(2): 219-232. DOI:10.1029/RG020i002p00219
[27] Li A, Forsyth D W, Fischer K M. Shear velocity structure and azimuthal anisotropy beneath eastern North America from Rayleigh wave inversion. J. Geophys. Res. , 2003, 108(B8): 2362. DOI:10.1029/2002JB002259
[28] Kennett B L N, Engdahl E R, Buland R. Constraints on seismic velocities in the Earth from traveltimes. Geophysical Journal International , 1995, 122(1): 108-124. DOI:10.1111/gji.1995.122.issue-1
[29] Liu D Y, Nutman A P, Compston W, et al. Remnants of ≥3800 Ma crust in the Chinese part of the Sino-Korean craton. Geology , 1992, 20(4): 339-342. DOI:10.1130/0091-7613(1992)020<0339:ROMCIT>2.3.CO;2
[30] Zhao G, Wilde S A, Cawood P A, et al. Archean blocks and their boundaries in the North China Craton: lithological, geochemical, structural and PT path constraints and tectonic evolution. Precambrian Research , 2001, 107(1-2): 45-73. DOI:10.1016/S0301-9268(00)00154-6
[31] Ren J, Tamaki K, Li S, et al. Late Mesozoic and Cenozoic rifting and its dynamic setting in Eastern China and adjacent areas. Tectonophysics , 2002, 344(3-4): 175-205. DOI:10.1016/S0040-1951(01)00271-2
[32] Zhao L, Zheng T. Seismic structure of the Bohai Bay Basin, northern China: Implications for basin evolution. Earth and Planetary Science Letters , 2005, 231(1-2): 9-22. DOI:10.1016/j.epsl.2004.12.028
[33] 臧绍先, 刘永刚, 宁杰远. 华北地区岩石圈热结构的研究. 地球物理学报 , 2002, 45(1): 56–66. Zang S X, Liu Y G, Ning J Y. Thermal structure of the lithosphere in North China. Chinese J. Geophys. (in Chinese) (in Chinese) , 2002, 45(1): 56-66.
[34] Forsyth D W, Webb S C, Dorman L M, et al. Phase velocities of Rayleigh waves in the MELT experiment on the East Pacific Rise. Science , 1998, 280(5367): 1235-1238. DOI:10.1126/science.280.5367.1235
[35] Weeraratne D S, Forsyth D W, Fischer K M, et al. Evidence for an upper mantle plume beneath the Tanzanian craton from Rayleigh wave tomography. J. Geophys. Res. , 2003, 108(B9): 2427.
[36] Saito M. DISPER80: A subroutine package for the calculation of seismic normal-mode solution. // Doornbos D J eds. Seismological Algorithms: Computational Methods and Computure Programs. New York: Elsevier. 1988: 293-319.
[37] Jiang M M, Zhou S Y, Sandvol E, et al. 3-D lithospheric structure beneath southern Tibet from Rayleigh-wave tomography with a 2-D seismic array. Geophysical Journal International , 2011, 185(2): 593-608. DOI:10.1111/gji.2011.185.issue-2
[38] Yin A, Harrison T M. Geologic evolution of the Himalayan-Tibetan orogen. Annual Review of Earth and Planetary Sciences , 2000, 28(1): 211-280. DOI:10.1146/annurev.earth.28.1.211
[39] Bilham R, Larson K, Freymueller J, et al. GPS measurements of present-day convergence across the Nepal Himalaya. Nature , 1997, 386(6620): 61-64. DOI:10.1038/386061a0
[40] 陈九辉, 刘启元, 李顺成, 等. 青藏高原东北缘—鄂尔多斯地块地壳上地幔S波速度结构. 地球物理学报 , 2005, 48(2): 333–342. Chen J H, Liu Q Y, Li S C, et al. Crust and upper mantle S 2 wave velocity structure across Northeastern Tibetan Plateau and Ordos block. Chinese J. Geophys. (in Chinese) (in Chinese) , 2005, 48(2): 333-342.
[41] Zhang Y Q, Mercier J L, Vergély P. Extension in the graben systems around the Ordos (China), and its contribution to the extrusion tectonics of south China with respect to Gobi-Mongolia. Tectonophysics , 1998, 285(1): 41-75.
[42] Zhang Y, Ma Y, Yang N, et al. Cenozoic extensional stress evolution in North China. Journal of Geodynamics , 2003, 36(5): 591-613. DOI:10.1016/j.jog.2003.08.001
[43] Zhao L, Zheng T. Using shear wave splitting measurements to investigate the upper mantle anisotropy beneath the North China Craton: Distinct variation from east to west. Geophys. Res. Lett. , 2005, 32(10). DOI:10.1029/2005GL022585
[44] 唐有彩, 冯永革, 陈永顺, 等. 山西断陷带地壳结构的接收函数研究. 地球物理学报 , 2010, 53(9): 2102–2109. Tang Y C, Feng Y G, Chen Y S, et al. Receiver function analysis at Shanxi Rift. Chinese J. Geophys. (in Chinese) (in Chinese) , 2010, 53(9): 2102-2109.
[45] 田宝峰, 李娟, 王为民, 等. 华北太行山区地壳各向异性的接收函数证据. 地球物理学报 , 2008, 51(5): 1459–1467. Tian B F, Li J, Wang W M, et al. Crust anisotropy of Taihangshan mountain range in north China inferred from receiver functions. Chinese J. Geophys. (in Chinese) (in Chinese) , 2008, 51(5): 1459-1467.