引言

地震是在地下深处介质发生的破裂错动或破裂扩展，依据目前科学技术条件，还无法实现直接测量震源区应力场大小，因此需根据震源机制、应力降、视应力、破裂半径等参数研究震源区应力场大小（陈学忠，2005）。随着数字化地震台站的建设，大量的中小地震甚至微震资料被完整记录，利用中小地震资料求解震源参数，已逐渐成为数字地震学中的日常工作。从原理上来讲，震源参数的求解较简单且较成熟，Brune（1970，1971）给出了通过拟合观测谱得到震源参数的理论基础，在此基础上利用数字化地震资料对震源参数结果的应用研究也是近年来国内外学者感兴趣的焦点（Jin等，2000；Moya等，2000；Bindi等，2001；Choy等，2004；华卫等，2009；赵翠萍等，2011；李艳娥等，2015；郑建常等，2016；王鹏等，2020；张正帅等，2020；臧阳等，2021）。

郯庐断裂带是我国东部地区规模最大的NNE走向的深大活动断裂带（图1（a）），在我国境内长达2 400多km（高维明等，1980；晁洪太等，1994；朱成林，2020），规模宏大，结构复杂，是我国东部地区重要的地震构造带，位于我国东部经济发达人口稠密地区，其地震危险性分析有着强烈的社会需求，一直为地学界和地震学界关注的重点。郯庐断裂带活动性具有明显的分段性，各段的地震活动及构造运动存在较大的差异（郑朗荪，1988）。郯庐断裂带在山东境内的部分又称为沂沭断裂带（以下简称沂沭带，YSFZ），其为郯庐断裂带中切割最深、规模最大、新构造活动最强烈的段落（满洪敏，2005），介于鲁西断块（LXUP）与鲁东隆起（LDUP）和苏鲁造山带（SLOB）之间。同时，沂沭带是强活动地震带，自古至今沂沭带及其附近两侧大大小小的地震从未间断，公元前70年安丘M_S7.0地震和1668 年7月25日郯城M_S8½地震更是波及半个中国，1829年11月19日益都M_S6.0地震、1995年苍山M_S5.2地震及2022年5月1日青州M_L4.1地震群，其震中无一例外地落在沂沭带及其分支断裂上，这显示了沂沭带是现今仍在活动的断裂，也是山东地区最大的地震源地。

图 1  数据资料

Figure 1.  Data sheet

下载: 全尺寸图片 幻灯片

针对沂沭带具有的强震背景，本文搜集整理了2010—2020年沂沭带范围内$ {M}_{\mathrm{L}}\geqslant2.0 $共93次地震事件波形和震相观测报告，采用高频截止（High-cut）震源模型，稳健拟合$ {M}_{\mathrm{L}}\geqslant 2.0 $地震事件的震源谱参数（拐角频率）和震源参数（应力降、地震矩、矩震级、破裂半径及视应力等），得到了谱参数和震源参数之间的定标关系，并对沂沭带应力的时空演化特征进行分析，结合速度结构和地震事件分布对沂沭带强震孕育和构造机理进行研究，对进一步了解沂沭带孕震环境具有现实意义与科学意义。

1.   方法原理

地震波传播过程中不仅随着传播距离出现几何扩散，还会受到传播路径上介质的吸收和散射等影响，接收台站下方的浅层介质也会对地震波产生影响，所以地震计记录的地震波数据是上述各种过程的综合反映（赵翠萍等，2011；张正帅等，2020）。在频率域内，地面运动的位移谱可表示为：

式中，D_ij为第i个地震在第j个台站的观测谱；D_i0为第i个地震的震源谱；P_ij为地震波的传播路径响应项，包括地震波的几何扩散和非弹性衰减；地震波的几何扩散中G_ij为球面几何衰减模型，${G}_{i j}={R}_{i j}^{-1}$，描述了地震波在传播过程中随着波阵面的扩大，波阵面单位面积的能量随着震源距R的增大逐渐减小的参数；R_ij为震源距。地震波的非弹性衰减中Q_ij为介质的品质因子，描述介质的不均匀性和非完全弹性，代表了地震波在传播中能量的损耗；$ \nu $为地震波速，本文使用S波进行计算，取$ \nu =3.2\;\text{}\text{km/s} $；S_j为第j个台站的场地响应；I_ij为第j个台站的仪器响应函数。

消除震源至台站间的传播路径响应、场地响应及消除仪器响应后，可通过地震观测谱D_ij获得震源谱D_i0。与震源参数有关的很多物理量可在震源谱中得到，而这些物理量的测定结果与选取的震源模型有关，经典Brune模型中，理论震源位移谱表示为式（3），其中γ为高频衰减系数。

使用式（4）对地震观测谱进行拟合，可获得4个特征参数，分别为位移谱中的零频极值$ {\varOmega }_{0} $、谱衰减系数为2时的拐角频率$ {f}_{{\rm{c}}} $、截止频率$ {f}_{\mathrm{m}\mathrm{a}\mathrm{x}} $及高于截止频率$ {f}_{\mathrm{m}\mathrm{a}\mathrm{x}} $部分的频谱衰减系数$ N $，High-cut震源模型充分考虑理论震源谱的4个特征参数，理论上可对震源谱的拟合有更好的约束。

2.   数据与计算

本文选择山东测震台网记录到的沂沭带及邻区范围内2010—2020年112个${M}_{\mathrm{L}}\geqslant 2.0$的地震事件波形和地震观测报告用于计算震源参数。为保证计算结果的可靠性，对数据进行预处理时，遵循以下原则：首先，选择震中距＜200 km的台站数据，以确保波形有较高的信噪比，且至少有4个台站记录到波形，经筛选共得到93个地震事件；其次，选择无短时畸变噪声的地震时间波形，避免噪声对观测谱造成污染；最后，由于S波段具有较P波段更高的信噪比，本文截取S波段数据并进行Kaiser加窗处理，以保证数据幅值和频率的精度（王鹏，2019）。

从地震波形得到地震观测谱的步骤如下：首先，对台站地震记录波形进行去趋势、去除直流分量、消除仪器响应处理；然后，根据台网观测报告的S波到时信息自动截取S波段，通过震中与台站位置计算方位角，从而将原始NEZ坐标系三分量地震记录旋转到RTZ坐标系，使用T分量波形，即SH波段进行分析；最后，对观测数据进行余弦边瓣加窗处理，使用快速傅里叶变换得到观测谱，获得观测谱后依次消除几何衰减、非弹性衰减和场地响应后得到地震震源谱。苗庆杰等（2016）使用Moya方法计算得到了山东地区56个台站场地响应（图2），本研究取平均值1.2用于消除场地响应。另外，根据G(R)=R⁻¹消除几何衰减效应后，使用$ Q\left(f\right)=457.1{f}^{0.4317} $消除非弹性衰减效应（苗庆杰等，2016），最终由至少4个台站观测谱得到平均观测震源谱。2012年5月10日莒县M_L2.3地震计算结果如图3所示，由图3可知，记录到本次地震的台站具有较好的包围性，说明得到的平均震源谱较可靠。

图 2  山东地区56个台站的场地响应（苗庆杰等，2016）

Figure 2.  The site responses of 56 stations in Shandong province（Miao et al., 2016）

下载: 全尺寸图片 幻灯片

图 3  2012年5月10日莒县M_L2.3地震计算实例

Figure 3.  Calculation example of the Juxian M_L2.3 earthquake on May 10, 2012

下载: 全尺寸图片 幻灯片

3.   结果分析

利用得到的震源谱参数，结合经验公式，可得到该地震对应的震源参数，主要包括地震矩$ {M}_{0} $、破裂半径$ R $、应力降$ \mathrm{\Delta }\delta $和视应力$ {\delta }_{\mathrm{a}\mathrm{p}\mathrm{p}} $等：

式中，ρ为介质密度，本文取ρ=2.67 g/cm³；β为S波速度，β本文取=3.2 km/s；$ {R}_{\theta \phi } $为震源辐射图型因子，在全震源球上的平均震源辐射因子$ {R}_{\theta \phi } $=0.48；μ为剪切模量（对于地壳介质，μ取30,000 MPa）；E_S为地震能量，可由对速度谱的平方积分求得，即

式中，V(f )为速度谱，利用$V(f)=\sqrt{{V\left(f\right)}_{{\rm{Z}}}^{2}+{V\left(f\right)}_{{\rm{NS}}}^{2}+{V\left(f\right)}_{{\rm{EW}}}^{2}}$求得。

本研究共计算了2010—2020年沂沭带范围内$ {M}_{\mathrm{L}}\geqslant 2.0 $共93次地震的震源参数。随机选取的部分地震震源谱如图4所示，图中红色线为所有台站的平均震源谱，蓝色线为理论震源谱。由图4可知，各次地震平均震源谱和理论震源谱较吻合，说明此次计算得到的震源谱和震源参数具有可靠性。93次中小地震震源参数之间的相关关系如图5所示。

图 4  部分地震震源谱拟合结果

Figure 4.  Partial seismic source spectrum fitting results

下载: 全尺寸图片 幻灯片

图 5  震源参数之间的定标关系

Figure 5.  Scaling relation among source parameters

下载: 全尺寸图片 幻灯片

3.1   拐角频率

拐角频率是震源谱中高频和低频的交点，主要反映了地震波高低频能量的分布特征，本研究拟合得到的沂沭带地区中小地震的拐角频率为2～15 Hz，与矩震级有较好的相关关系（图5（a）），震级越大，拐角频率越低，这与很多研究结果一致（Droue等，2005；Izutani，2008）。拐角频率、应力降和地震矩之间存在明显的依赖关系，如图5所示。另外，拐角频率不仅与震源区应力状态有关，也与发震过程有一定关系。

拐角频率随时间的变化曲线如图6所示，图中给出了拐角频率95%的置信区间。由图6可知，总体上来说，震级越小，受噪声干扰影响越大，拐角频率不确定性越大。2016年1月8日M_L3.6地震发生后，至2020年12月31日研究区内未发生M_L3.0以上地震的时间段内，震级偏小，拐角频率不确定性较大。

图 6  拐角频率随时间的变化曲线

Figure 6.  Variation of corner frequency with time

下载: 全尺寸图片 幻灯片

3.2   地震矩与破裂半径

地震矩用于表示地震强度，与地震释放的能量有关。由图5可知，本次计算的沂沭带地区地震矩M₀为10¹¹～10¹⁴ N·m，与赵翠萍等（2011）给出的中国大陆地区地震矩为l0¹²～10¹⁶ N·m（M_W2.0～5.0）结论相当。$ {M}_{0} $与$ {M}_{\mathrm{L}} $在单对数坐标系下存在较好的线性关系，利用稳健函数线性拟合，二者关系为$ \mathrm{log}{M}_{0}=1.096{M}_{\mathrm{L}}+ 9.78 $。由震源谱得到的破裂尺度虽存在争议（Beresnev，2001，2002），但其仍为认识地震破裂尺度提供了新的途径，图5（d）给出了本文得到的结果，中小地震的破裂半径为80～600 m，破裂半径与矩震级之间存在明显的相关关系，震级越大，释放的能量和产生的破裂尺度越大，其相关关系形态与理论关系形态一致（陈运泰等，1976）。

3.3   视应力

视应力的概念是由Wyss等（1968）于1968年提出来的，是与震源动力学过程有关的物理量，属于震源的动力学参数。对某个区域中引起地震滑动的视应力水平进行区域平均，可作为当地应力水平的间接估计（吴忠良，2001），即视应力越高，震源区的应力水平越高。

视应力是否随地震强度的增大而增加是震源参数研究中广泛讨论的重要问题之一，不同的研究学者得出了不同的结论，目前仍没有确切的定论。沂沭带地区得到的视应力与震级之间存在一定的正相关关系（图5（e）），在低震级的部分视应力较离散，考虑到震级小时波形信噪比低，计算误差可能较大，这种趋势关系在许多研究中是较常见的。本文计算的沂沭带地震视应力为0.04～1.4 MPa，均值为0.38 MPa，与王鹏等（2015）计算的沂沭带范围内43个M_L≥2.0地震视应力结果相当，沂沭带目前属于低应力释放地区。

视应力随时间的演化如图7所示，由图7可知，视应力相对较小，为0.04～1.4 MPa，大部分集中在1 MPa以下。2010年以来，3.0级以上地震共发生9次，其中M_L≥3.5地震5次，分别为2011年6月20日3.7级地震、2012年7月13日3.8级地震、2012年8月8日3.9级地震、2013年12月11日3.6级地震和2016年1月8日3.6级地震，上述5次地震视应力均高于研究区视应力平均值，且上述地震中的3次地震视应力＞1 MPa，由于视应力反映了当地应力水平的下限估计，2012—2014年视应力较高，表明沂沭带内应力水平有一定增强现象，随后发生的一系列中小地震中，2015年小震活动较弱，2016年1月3.6级地震发生后，应力积累水平得到一定程度的释放，此后视应力均＜1 MPa，至2020年12月31日，沂沭带地区地震强度有所降低，未发生M_L3.0以上的地震。

图 7  视应力随时间演化曲线

Figure 7.  Variation of apparent stress with time

下载: 全尺寸图片 幻灯片

本文将得到的93个地震差视应力在空间范围内以0.1°的间隔进行线性插值，得到了差视应力空间分布，如图8所示。由图8可知，沂沭带中南部的莒南-临沂和沂沭带北段的安丘段存在明显高视应力集中区，与李艳娥等（2015）研究结论相似。视应力是对研究区近期应力状态的估计，而长期的应力状态由b值体现，低b值反映高应力（易桂喜等，2013）。由王鹏等（2015），视应力及b值有较好的相关性，2011年5月20日安丘M_L3.7地震和2016年1月8日临沭M_L3.6地震的视应力高值反映了安丘段和莒南-临沂附近的高应力水平，同时，该区域处于低b值状态，与本研究结果一致。

图 8  视应力空间分布特征

Figure 8.  Spatial distribution of differential apparent stress

下载: 全尺寸图片 幻灯片

将1970年以来研究区内M_L≥3.5地震投影到图8中，可发现这些地震多发于高视应力内部或边缘地带，尤其是2010年以来研究区内共发生3.5级以上地震5次，均位于高视应力的边缘地带。

由图9（a）可知，地震基本分布在深度10、23 km左右，较深的地震视应力整体普遍偏大。由图9（b）、（c）可知，1970—2016年32次M_L4.0以上地震有70%发生在高低速过渡带内。本文得到的93次地震高视应力空间分布也主要集中于高低速过渡带或低速体附近，根据已有研究（齐诚等，2006；Chen等，2014），高低速过渡带中强地震的频发主要源于流变边界内物质不均匀易于应力集中，而视应力可作为区域应力水平的间接估计，从沂沭带内高视应力来看，也能进一步体现出应力在高低速过渡带内的集中特征。

图 9  视应力深度分布与地震波速度结构关系

Figure 9.  Relationship between apparent stress depth distribution and seismic wave velocity structure

下载: 全尺寸图片 幻灯片

4.   讨论与结论

本文利用山东台网2010—2020年的中小地震波形资料，采用基于Brune模型的高频截止震源模型，使用稳健的最小二乘法拟合研究区内的震源谱，获得了93个沂沭带及邻区的震源参数，经过分析可得到以下结论：

（1）地震矩$ {M}_{0} $与震级$ {M}_{\mathrm{L}} $在单对数坐标系下存在较好的线性关系，利用稳健函数线性拟合，可表示为$ \mathrm{log}{M}_{0}=1.096{M}_{\mathrm{L}}+9.78 $。

（2）沂沭带地区中小地震的拐角频率为2～15 Hz，与矩震级有较好的相关关系，震级越大，拐角频率越低；震级越小，受噪声干扰影响越大，拐角频率的不确定性越大。

（3）沂沭带地震视应力为0.04～1.4 MPa，均值为0.38 MPa。

（4）地震视应力的时间分布特征表明，2012—2014年沂沭带地震视应力偏高，地震发生频次明显增加，2016年后地震视应力有所降低，地震强度和频次有所降低。

（5）地震视应力的空间分布特征表明，地震多发生在高视应力内部和边缘地区，沂沭带中南部的莒南—临沂和沂沭带北段的安丘段存在明显高视应力集中区，这些区域为沂沭带未来可能发生中强地震的潜在震源区。