利用接收函数反演东北地区地壳结构

李秀丽，朱叶琳，高业欣，燕 云，赵龙梅

（辽宁省地震局，辽宁 沈阳 110034）

0 引言

东北地区是指黑龙江省、吉林省、辽宁省和内蒙古自治区东部地区。东北做为全球鲜有的深震孕育区，以及火山活动剧烈地区，地貌交错、构造情况较为复杂备受地学家关注，能够获取该区的速度结构，对于研究深层地震孕育机制、板块运动、火山活动等都有重要的意义[2]。长期以来，前人利用不同的观测资料及方法在东北地区开展了大量的研究，并取得了一系列成果[4]。

以东北区域内不同的地块作为研究范围，分别采用P 波接收函数CCP 法、体波层析成像法、H-K 叠加法、近震及远震走时反演法等方法获得了东北地区1100 km 深度的地幔结构模型、东北地区地壳上地幔三维P 波速度模型、东北下方达到800 km 处的P 波和S 波速度结构、中朝地壳上地幔三维P 波速度结构、地壳厚度与平均速度比、地壳厚度图像等研究成果。发现了东北地区壳幔结构上体现的一些共性问题：中国东北地区地壳厚度东西布局差异明显。异常方向北东向与间断面方向一致。壳幔介质存在明显的横向不均匀性：低速突出异常区在长白山、五大连池火山区域下方，高速异常区在松辽盆地地块下方。

1 数据

受区域国界及海岸线的限制，我们选取的地震台站都位于东北域内以及西部邻省。我们选取了黑龙江、吉林、辽宁、内蒙东部共111个测震台站的宽频带数字化地震计观测资料，观测资料由中国CDSN 波形数据中心提供。查阅相关资料及用2012 年数据进行初算，发现中等规模的浅源地震得到的结果比较好，最后向国家数据备份中心选取上述台站自2012 年1 月1 日至2015 年10 月1 日，震中距在2000～18000 km 范围内，震级在6.0 以上，共542 个地震事件（事件震中分布图见图1）。

首先对本文研究所用数据进行初步处理，依次读取地震事件观测波形，删除伪地震事件记录，调整各台站接收波形的到时，调整台站方位角翻转波形，删除初至不明显事件，最后得到287 个波形发育良好，初至明显的地震事件用于接收函数研究。经过计算和筛选共获得4683 条接收函数，较好地反演了研究区域台站下方地壳厚度及速度比。

图1 接收函数所用事件震中分布图Fig.1 The epicenter distribution of the events used by the receiver function

2 接收函数提取

远震P 波在传播时会受到震源函数、测震仪器响应函数、传播路径、传播介质等因素共同影响，而当垂直向可以近似为脉冲响应时，就可以将远震P 波接收函数三分量地震记录其他两个分量进行反褶积，进而消除了传播路径和震源时间函数的影响，分离出了台站下方介质所形成的间断面对入射P 波的响应，得到台站下方的介质信息。接收函数的提取，实际上就是得到包含Ps 转换波及其他多次反射震相函数的过程，过程中要分离出路径及震源因素序列。

接收函数的计算方法因研究域不同分为频域计算和时间域计算，时间域计算方法多由频率域经傅里叶变换得到。由于实际观测仪器采用频率不同导致现有观测资料带宽有限，而且地震波形记录中包含了噪声，体波垂直分量频谱存在零点和很小值，针对以上情况，学者提出了各种改进方法以获得较可靠稳定的接收函数。Wiggin 提出了“水准量反褶积”方法，在计算中引入高斯滤波因子和水准量，使得接收函数结果更加精确和可靠。然而该方法虽然获得较稳定的接收函数，但由于水准量的引入降低了弱转换震相的分辨率，因此学者们发展了多种时间域计算方法，如同步时间反褶积法、最大熵反褶积法、脉冲反褶积法、迭代反褶积法、Wiener 滤波反褶积法等。

本文接收函数提取采用最大熵反褶积法。首先对三分量的观测记录进行预处理，然后截取-20 s 至120 s 的波形，并对两个水平分量进行坐标旋转，以及反褶积计算，得到径向和切向接收函数。本部分计算内容均通过脚本调用sac 命令实现。

首先：截取波形数据。利用IASP’91 模型得到计算P 波理论到时，读取SAC 文件，去除理论到时，并修改头段变量中时间参数，取初至前10s，P 波后100s 为波形长度截取波形。

第二步：对所选波形去除仪器响应、去均值处理，进行带通滤波。

第三步：旋转坐标。检查并修改两个分量的头段变量cmpinc 和cmpaz，使其正交。使用rotate 命令旋转坐标。

第四步：利用反褶积计算接收函数。查阅文献及该区域相关研究资料，设置高斯参数为3.0，迭代次数设置为150 次。设置拟合度大于90%，初步对接收函数进行挑选。

最后，对程序计算的接收函数进行二次挑选，删除Ps 震相不清楚、幅值偏离总体幅值的接收函数，同时检查程序处理相反的情况，挑选完成后重新绘制台站接收函数及结果图。

经过计算和二次挑选（部分台站经过多次挑选）共得到111 个台站4683 条接收函数，其中5 条以下台站仅2 个，10 条以下台站7 个，10～50 条台站58 个，44 个台站接收函数条数大于50，获得了较好的结果，部分台站接收函数结果如图2 所示。

3 接收函数HK 扫描

地壳厚度H 与地壳Vp/Vs 速度比存在类似反比关系，当地壳比较薄时，波速比较大，当波速比较小时，地壳比较厚，震相对直达P 波和Ps 波区别不大，然而从单条接收函数上获取较为准确的Ps 多次转换震相是非常不科学的。因此，我们通过叠加取平均的方法来获得较为可靠的Ps 震相到时。有人提出了H-Kappa 方法来求取地壳厚度。该方法利用Ps、PpPs、PpSs+PpS 震相到时同时叠加来求取地壳厚度，叠加过程中，有效地消除了介质的横向不均匀引起的散射效应，所以该方法可以获得更加可靠的速度界面及深度信息。本文采取该方法进行研究。

本部分计算内容均通过脚本调用sac 命令实现，设置权系数w1为0.5，w2为0.3，w3为0.2 进行叠加运算及绘制波速比图。

最后得到东北地区MOHO 面分布如图3 所示，各台站波速比信息如表1 所示。

图4 部分台站H-K 扫描结果Fig.4 H-K scan results of some seismic stations

4 结论

从计算结果可知，研究区最大地壳厚度55 km，最小地壳厚度25 km。但从图4 可以看出，东北地区地壳厚度多分布在25～40 km 中间，MOHO 面与断层方向一致呈北东方向展布，并呈现中部隆起，东西两侧下凹的特点。东部地壳厚度介于35～40 km 之间，中部介于25～32 km 之间，西部佳木斯地块介于35～37 km之间，三江盆地附近在25～32 km 之间，该总体趋势与张广成等[6]所述一致。

重力梯度带以西地壳较厚，介于35～40 km，重力梯度带以东，地壳厚度较薄，介于26～33 km之间，反映了研究区地壳厚度横向差异较大。该结果与江为为等[7]重力异常测量结果一致，杨宝俊等[9]人工地震测深研究表明，重力梯度带是地壳的陡变带，地壳厚度变化是引起重力改变的一个原因。

在重力梯度带以东到兴蒙槽地，地壳厚度逐渐减小，在小兴安岭北端达到最小25 km，小兴安岭与松辽平原过渡地带，地壳厚度也较周边地块薄，推测与复杂构造拉伸关系密切。

佳木斯地块地壳厚度较厚，约为36～40 km，较其周边地区呈现出偏厚的特点，推测可能与小兴安岭及长白山褶皱带形成相关，褶皱带挤压变形过程使得该地区地壳厚度加厚。

而三江盆地区域地壳又变薄，约为28～30km，反映地壳横向差异较大。

松辽盆地内部也表现出不同区域的莫霍面差异，除几个台站下方显示地壳较厚外，其北、东、地壳较薄，莫霍面的形态与地表形态比较一致，厚度约为25～30 km，可能是板块运动过程中的拉张力作用结果。

下辽河盆地除一个台站外，地壳厚度较薄，在28～32 km 深度上变化，反映了该区域复杂的地壳变形过程。该结果与张广成等[6]相一致。

大兴安岭褶皱带地壳厚度变化幅度较小，厚度值在33～36km 之间。小兴安岭、张广才岭地壳厚度总体较薄，厚度值在25～32 km 之间。

长白山地区地壳厚度较厚，对比泊松信息，发现该地区泊松比值也很高，熊小松等[4]曾经在研究中推测该现象可能与火山岩浆活动有关，使地壳处于热状态。

总体来说，重力带以西、佳木斯地块及长白山地块地壳较厚、松辽盆地北、东、南地区地壳厚度较薄，与获得的地壳厚度结果基本一致。

从莫霍面分布图可以看出，研究区四个大型盆地（海拉尔盆地、松辽盆地、下辽河盆地、三江盆地）及主要山脉（大兴安岭、小兴安岭、张广才岭）的地壳厚度与地形起伏有非常好的“镜像”关系，与江为为等[7]重力磁场资料分析结果相映证。