2013-10-31~11-23吉林省松原市前郭尔罗斯蒙古族自治县连续发生了一系列地震，其中M_S≥5.0地震就有5次，最大震级为M_S5.8(表 1)。根据历史地震目录记载，松原地区地壳稳定，地震发生的频率远小于其他地区，强度也不大，该地震群发生之前震区附近仅发生过两次中强地震^[1]，震中位置见图 1。松原地区地震的成因，引起了研究者们的关注，高金哲等^[2]推测可能与油气开采有关。对于2013年吉林松原地震群的发震机制，不少学者从地震精定位、震源机制解、发震构造等多个方面进行研究，得到不同的结论^[3]。

前人利用蕴含丰富地震信息的数字化地震波形资料，在不同类型地震的识别方面做了一定工作^[4-5]，本文在参考前人研究成果的基础上，从地震波形的时域、频域及震源参数等地震学特征入手，综合震区周围环境特征，探究2013年松原地区几次M_S≥5地震的成因。

1 研究方法及数据 1.1 研究方法

由于震源体、传播路径、场地等的不同，台站接收到的地震波形信息十分丰富，通过对地震波进行分析，可对构造地震和因爆破、抽水、蓄水、采矿等人为因素引起的非天然地震加以区分。在台站j的地震仪记录的第i个地震在时间域里可表达为^[6]：

$ \begin{array}{l} U_{ij}^{obs}(t) = ({S_{ij}}(t) * P{}_{ij}{\rm{(}}t{\rm{)}} * \\ {L_j}{\rm{(}}t{\rm{)}} * {I_j}{\rm{(}}t{\rm{)}}) \cdot {K_j} \cdot Rij{\kern 1pt} \end{array} $

(1)

式中，$U_{ij}^{obs}(t)$为位移，S_ij(t)为第i个地震的震源项，P_ij(t)为第i个地震的震源与第j个台站之间的传播路径效应，L_j(t)为第j个台站的局部场地效应，I_j(t)为第j个台站的仪器响应，K_j为自由表面放大因子，R_ij为地震波的辐射花样(取决于震源机制和台站相对于震源的位置)，表示褶积。可通过傅里叶变换将式(1)在频率域中表达为：

$ \begin{array}{l} U_{ij}^{obs}(f) = {S_{ij}}(f) \cdot P{}_{ij}{\rm{(}}f{\rm{)}} \cdot {L_j}{\rm{(}}f{\rm{)}} \cdot \\ {I_j}{\rm{(}}f{\rm{)}} \cdot {K_j} \cdot \left| {{R_{ij}}} \right| \end{array} $

(2)

去除仪器响应后为：

$ U_{ij}^{obs}(f) = {S_{ij}}(f) \cdot P{}_{ij}{\rm{(}}f{\rm{)}} \cdot {L_j}{\rm{(}}f{\rm{)}} \cdot {K_j} \cdot \left| {Rij} \right| $

(3)

通过式(3)可以看出，若x、y地震的震源足够接近，即两次地震震源的距离比震源距R_ij小得多，则可认为台站j记录到的x、y地震的观测位移谱U_xj、U_yj只与各自的理论震源谱S_ij(f)和辐射花样R_ij中的震源机制部分有关，传播路径效应P_ij(f)和台站的局部场地效应L_j(f)及自由表面放大系数K_j可以消除，也就是说，此时两次地震的观测位移谱可以反映各自的震源特征。

根据Brune模型^[7]，地震i的理论震源谱可表示为：

$ {S_i}{\rm{(}}f{\rm{)}} = \frac{{{\Omega _{{\rm{0}}i}}}}{{{\rm{1}} + {{{\rm{(}}f/{f_{Ci}}{\rm{)}}}^{\rm{2}}}}} $

(4)

式中，Ω_0i为地震i的震源谱低频渐近线值，即零频极限值；f_Ci为低频渐近线和高频渐近线交点处的频率，即拐角频率。

对于位移震源谱的计算，本文首先截取地震波Sg到达至基本衰减结束为S波段，经带宽为0.2~22 Hz的4阶Butterworth带通滤波后，将S波段分段，对每小段扣除仪器响应，进行快速傅里叶变换，加汉宁窗以弥补能量损失，再按照邻段50%重叠的规则进行移窗。同样操作后得到每小段的傅里叶谱，叠加起来即得到S波段傅里叶谱。利用同样方法，对P波初动前256个采样点进行处理，得到噪声段的傅里叶谱，扣除噪声后将速度型数据除以2πf换算成位移型数据，即得到地震波的观测位移震源谱。根据Brune的ω^－2震源模型^[7]，使拟合的理论位移震源谱和观测位移震源谱残差最小，即可得到零频极限值Ω_0i和拐角频率f_Ci。

结合地震的理论震源谱，通过计算2013年松原5次M_S≥5地震的零频极限值Ω_0i之间的相关系数分析这5次地震震源特性的相似性，若两次地震的震源特性相似，则相关系数r_xy应较大，接近1；若震源特性不相似，则相关系数r_xy较小^[8]。地震x、y的零频极限值相关系数r_xy可表示为^[9]：

$ {r_x}y{\rm{ = }}\frac{{\sum\limits_{j = 1}^n {\sum\limits_{l = 1}^2 {({x_{jl}} - \bar x)({y_{jl}} - \bar y)} } }}{{\sqrt {\sum\limits_{j = 1}^n {\sum\limits_{l = 1}^2 {{{({x_{jl}} - \bar x)}^{\rm{2}}}} } } \sqrt {\sum\limits_{j = 1}^n {\sum\limits_{l = 1}^2 {{{({y_{jl}} - \bar y)}^{\rm{2}}}} } } }} $

(5)

式中，x_jl、y_jl分别为台站j记录的地震x、y的l分量波形资料反演得到的零频极限值Ω_0xj和Ω_0yj，x、y分别为x_jl、y_jl的平均值，n为台站数目，l取1~2(分别为垂向和水平向S波)。此外，理论震源谱的拐角频率和高频衰减系数常用来区分不同成因的地震，本文计算了2013年5次松原地震和一些天然构造地震及塌陷地震的拐角频率和高频衰减系数，试图对2013年松原5次M_S≥5地震的成因机制作出判断。

不同成因地震的震源特性不同，在地震波形时域分布特征及频率域的分布特征上都有所体现^[10]。一般来说，在时域特征方面，塌陷地震较构造地震长周期面波发育、衰减较慢，波形比较简单；在频域上的差异体现在塌陷地震高频成分不发育，能量主要集中在低频段。在分析地震波形时域和频域的特征时，计算P波和S波的振幅比A_P/A_S及优势频率是一种常用的方法，但这些参数仅能反映某个时间点或某个频率点的幅值情况，具有较大的偶然性，分析结果不够准确。本文通过计算波形复杂度C及谱比值SR来反映地震波时域和频域的分布情况。波形复杂度C和谱比值SR的定义为^[11]：

$ {\rm{SR}} = \frac{{\int_{f = {L_1}}^{f = {H_1}} {\left| {X(f)} \right|} df}}{{\int_{t = {L_2}}^{t = {H_2}} {\left| {X(f)} \right|df} }} $

(6)

$ C = \frac{{\int_{t = L}^{t = H} {{x^2}} (t)dt}}{{\int_{t = {t_0}}^{t = L} {{x^2}} (t)dt}} $

(7)

式中，x(t)为地震信号t时刻的幅值，t₀为P波到达台站的时刻，L为P波结束时刻，H为S波基本衰减完时刻，X(f)为地震信号在f频率处的幅值，L₁、H₁、L₂、H₂分别为高频段和低频段的上、下限(本文采用的带限赫兹数为：L₁=0.2 Hz，H₁=1.2 Hz，L₂=1.2 Hz，H₂=6 Hz)。

地震的震源深度是地震学研究的核心问题之一，准确确定地震震源深度对于确定地震成因有非常重要的意义^[12]。因此，本文根据震中距较小的流动台资料，利用Pg波、Sg波到时差的方法检验松原地震的震源深度。

1.2 数据选取

本文以2013年松原5次5级以上的地震事件作为研究对象，再选取两次构造地震事件及两次塌陷地震事件作为对比，具体信息见表 2。由于2013年松原地区5次5级以上地震震中位置接近，为尽量减少传播路径和场地对计算震源谱的影响，选用震中距较小的固定台站乾安台接收到的地震波数据。

2 研究结果 2.1 震源谱特征分析

在位移震源谱中，拐角频率是对位移谱宽度的一种估计，反映低频与高频成分的分布特征，高频衰减系数则反映地震波中高频部分衰减的快慢程度，这两个参数基本上反映了震源体的特征^[13]。震源谱参数受方位角的影响不大，因此选择震中距最小的固定台站波形数据进行震源参数的计算(表 3)。可以看出，2013年松原地区M_S≥5.0地震的拐角频率和高频衰减系数较构造地震的小，与塌陷地震的较为接近。根据前人研究资料^[5]，构造地震拐角频率和高频衰减斜率比塌陷地震的大，因此从震源谱特征的角度来分析比较，松原地震的震源参数与典型构造地震的震源参数有一定差异，同典型构造地震相比，松原地震的震源谱具有衰减早、衰减快的特征，这种特征与塌陷地震更加类似。

此外，选择2013-10-31~11-23松原地区5次5级以上地震，根据乾安台(QAT)、长岭台(CLT)和松原台(SYT)等3个台站记录到的这5次地震的波形记录，反演各地震S波段垂向和水平向的零频极限值，并计算5次地震的相关系数，结果见表 4。

由表 4可知，在震源特性上，2013年松原5次M_S≥5地震事件中，第1次和第5次地震事件比较相似，第2、3、4次地震事件比较相似。

2.2 时域特征分析

首先将松原地震波形与构造地震及塌陷地震的波形进行宏观上的比较(图 2)，可以看出，松原地震的波形与构造地震的波形存在一定差异。一方面，松原地震衰减速度要比典型的构造地震衰减速度慢；另一方面，松原地震的波形周期较大，“毛刺”较少，主要是因为面波发育，频率成分比较单一，无高频波叠加，这些特征都与塌陷地震较为相似。通过计算波形复杂度(表 5)可以看出，塌陷地震的波形复杂度明显大于构造地震，而松原地震一些分量的波形复杂度又比塌陷地震大得多，说明虽然松原地震的波形在时域特征上与塌陷地震较为相似，但比典型的塌陷地震更为复杂。

2.3 傅氏谱特征分析

一般塌陷地震较构造地震的震源深度浅，而地震波在松软的沉积层中传播时高频部分会被吸收，因此塌陷地震在频域中的特征也与构造地震有着明显的区别^[4]。同样，将松原地震与构造地震及塌陷地震的傅氏谱进行宏观上的比较(图 3)可以发现，松原地震和塌陷地震的傅氏谱主要在低频段发育，优势频段集中在2 Hz以内；构造地震傅氏谱的优势频段较宽，大约为4 Hz以内，没有只集中在低频段的特点。对比地震的谱比值(表 5)发现，松原地震的谱比值明显小于构造地震的谱比值，与塌陷地震更为接近。由此可知，2013年松原地震与典型的构造地震明显不同，更加符合塌陷地震的特征。

2.4 震源深度检验

本文研究的5次松原地震中，只有2013-11-22发生的M_S5.3地震有近震流动台站记录，因此选取此次地震事件的一个流动台数据对震源深度进行检验。当震中距较小时(小于1~2倍震源深度)，可假设台站就位于震源的正上方，这样根据直达波法测得的震源距就可以看作震源深度的下限(真正的震源深度要浅于假定的震源深度)^[14]。利用这种方法，选择震中距为2 km的流动台站L2202台记录到的地震波形，计算得到Pg波和Sg波的到时差Δt为2 s，结合当地速度模型^[15]，选择P波传播速度为v_P=4 km/s，S波传播速度v_S=1.32 km/s，根据式(8)计算得到震源深度的下限d为6.9 km，即震源深度不超过6.9 km，而不是8 km(表 1)。

$ {\rm{d = }}\frac{{{\rm{\Delta t}} \cdot {\rm{v_p}} \cdot {\rm{v_s}}}}{{{\rm{v_p - v_s}}}} $

(8)

地震勘探资料显示，松原市乾安地区的采油井深度多在3.8 km左右^[2]，松辽盆地中生代沉积底部、基岩顶层深度约为6 km^[3]。此次地震的震源深度较浅，仍在油气圈闭的沉积盖层内。因此，在震源深度特征上，松原地区2013-11-22发生的M_S5.3地震可能与该地区油气开采有关。由于近震台站数量较少，不能准确地对5次地震事件的震源深度进行精确检验，但综合地震的时域和频谱特性来看，5次地震呈现出一致性，且根据地震的零频极限频率的相关系数来看，第2、3、4次地震的震源特性具有相似性，大致推测松原地区地震群与该地区油气开采有关。若日后震区台站加密，对震源深度的计算可以更加准确，对地震动的时频特性和震源特性的认识会更加清晰^[16]。

3 结语

本文结合2013年松原地区M_S≥5地震的台站记录信息，从不同方面对地震群的成因进行分析，得出以下结论：

1) 时空特征方面，震区附近历史上仅发生过两次5级以上地震，而2013年仅一个月内就发生5次5级以上地震，通过检验其中一次地震的震源深度发现，这种短暂突发、震群密集、震源浅的特点不符合该地区构造地震的特征。

2) 地震学参数方面，松原地震拐角频率为2.4 Hz左右，高频衰减斜率在－3.0左右，与构造地震相比呈现衰减早、衰减快的特征，而这种特征与塌陷地震类似。通过计算松原M_S≥5地震零频极限频率的相关系数发现，5次地震中第1次和第5次地震的震源特性相似性较高，第2、3、4次地震的震源特性相似性较高。此外，松原地震的波形周期大、面波发育、频率成分单一，优势频率主要集中在2 Hz以内，通过对波形复杂度和谱比值的计算发现，松原地震明显不同于构造地震，而是更加符合塌陷型地震特征。

通过对2013年松原地区M_S≥5地震的数字化波形时域、频域、震源参数进行分析，松原地震与典型构造地震有明显差别，其特征更加符合塌陷型地震。综合震区周边环境及震源深度推测，松原地震的成因与长期油气开采活动有关，但因近震台站数据及采油作业的具体资料有限，为进一步对地震成因机理进行分析，日后将加密震区观测，收集更多采油作业详细资料。

致谢: 感谢吉林省地震局和江西省地震局地震监测中心提供相关地震波形资料。