引言

区域构造情况和地壳应力状态是决定地震活动水平的重要因素，地壳应力变化可以改变地壳介质的属性，而地震波速比可以反映地壳介质性质的变化，因此利用波速比的变化预测未来较大地震的发生是重要的监测预报手段。早在20世纪50年代就已经有日本学者对大地震前的波速比变化情况进行了研究（冯德益等，1974，1981），地震波速比的变化有望提供有关震源和地震孕育过程的某些信息（Scholz等，1973；冯德益等，1975，1976）。近年来的研究显示在较大地震发生前波速比会有明显变化，陈丽娟等（2022）利用纵横波震源谱参数对2013年岷县-漳县6.6级地震前波速比进行了研究，发现震前一年波速比有较明显上升，到最高值后开始下降，临震前下降到平均值以下的较低水平；龙海英等（2011a，2011b）运用多台和达法对新疆乌苏5.1级与和静5.6级地震震前波速比进行了研究，发现震前波速比有较长时间的低值异常；张琳琳等（2016）运用多台和达法研究了2009—2014年天山地区中强地震震前的波速比变化，发现多数中强地震发生前波速比出现低值状态；王林瑛等（2014）和李艳娥等（2014，2016）运用单台、多台和直达波视速度等多种方法研究了汶川地震和芦山地震前后的波速比变化，发现强震前波速比均存在不同程度的异常变化；刘文邦等（2016）运用多台和达法和单台法对2016年青海门源6.4级地震震前波速比进行了研究，2种方法均显示波速比在震前出现了不显著的低值现象，这可能与多个台站将异常平均以及地震射线路径对震源区的覆盖有关。

目前常采用多台、多震方法（黎明晓等，2004；王林瑛等，2014；李艳娥等，2014，2016）或波速比与其他参数联合求解方法（Smith等，1983；Koch，1992）对波速比进行研究，由此求出的波速比是记录到地震震相的台站到地震震源之间的平均波速比。由于该段地震波传播距离一般较大，因此震源区波速比的变化量被稀释得很小，不易分析出异常情况。震源区的距离尺度较小，其地壳介质可认为具有相对一致性，地震纵波和横波速度在此距离上的变化较小，因此在震源区附近的波速比变化更清晰和准确。震群是一定范围内相对较短的时间里发生的密集地震，与大地震孕育过程相比，震群持续时间较短（数天到数年不等），震源区小震的密集发生，使得震源区介质不断摩擦、破碎，因此在震群持续时间内其震源区波速比的变化必然是显著的，如何消除地震波传播路径效应的影响而得到震源区波速比的变化是本文研究的重点。借鉴双差定位的思想，Lin等（2007）首次应用2次差分的方法分析了震源区的波速比变化；Dahm等（2014）采用2步反演技术，首先计算区域的平均波速比，然后使用最小中值平方法（LMS）基于双差后的到时数据计算了震源区的波速比；贾漯昭等（2017）应用双差波速比法分析了2014、2015年安徽金寨震群；郑建常等（2019）运用基于误差分布的双差波速比法分析了乳山震群的波速比变化。

贾漯昭等（2017）应用的双差波速比法未去除发震时刻的影响，客观上扩大了波速比计算的误差来源；郑建常等（2019）运用基于误差分布的双差波速比法在乳山震群中有较好的应用，但据作者反馈，有读者将其运用在其他震群中，计算结果偏差较大。本文为了提高计算波速比方法的普适性，应用郑建常等（2019）提出的双差波速比方法，运用去除误差分布后的纵横波双差到时二维正态分布拟合了震源区波速比，并认为反演波速比的关键在于二维正态函数相关系数的确定，最后本文运用以上方法对2017、2018年长岛震群震源区波速比进行了分析。

1.   理论与方法

1.1   双差波速比法简介

参考Lin等（2007）、Dahm等（2014）及郑建常等（2019）的研究，假设震群共发生了n次地震并被周边m个地震台站监测到，其中第i（i=1,2,3，…,n）个地震的发震时刻为$ {t}_{i} $，被第j（j=1,2,3，…,m）个台站记录到的地震纵横波到时分别为$ {T}_{ij}^{{\rm{P}}} $、$ {T}_{ij}^{{\rm{S}}} $，则第i个地震震源到第j个台站的纵横波走时${t}_{ij}^{{\rm{P}}}$、${t}_{ij}^{{\rm{S}}}$分别为：

第i个和第k个地震（记为“事件对d”）到共同记录到它们的第j个台站的纵横波走时差$\Delta T_{dj}^{\rm{P}}$、$\Delta T_{dj}^{\rm{S}}$分别为：

同理，第i个和第k个地震到共同记录到它们的第（j+1）个台站的纵横波走时差$\Delta T_{d,j + 1}^{\rm{P}}$、$\Delta T_{d,j + 1}^{\rm{S}}$分别为：

则上述2个走时差之差$ {\delta }_{d}^{{\rm{P}}} $、$ {\delta }_{d}^{{\rm{S}}} $分别为：

即：

最后得出波速比r计算公式为：

一般情况下，震群范围远小于震群到地震台站的距离，所以震群中的两次地震在震源区之外到台站的路径几乎相同（图1），由式（3）~式（6）计算的走时差约等于震群中2次地震震源之间某段距离的走时。式（3）与式（5）之差和式（4）与式（6）之差（即式（7）和式（8））可认为是震群中更小距离的走时，以消除式（7）、式（8）中第i、k个地震的发震时刻，公式中只剩下2个地震到2个台站的到时。式（9）、式（10）计算结果的误差仅与台站记录的地震到时有关，大大减少了误差的来源。

图 1  双差波速比法示意

Figure 1.  Schematic diagram of double differential wave velocity ratio method

下载: 全尺寸图片 幻灯片

1.2   函数随机误差计算

随机误差是用表征其取值分散程度的标准差评定的，函数的随机误差是用函数的标准差评定的。因此，函数随机误差计算是研究函数y标准差与各观测值$ \mathop x\nolimits_1 ,\mathop x\nolimits_2 ,...,\mathop x\nolimits_n $标准差之间的关系。

若函数的一般形式为$ {{y}} = f(\mathop x\nolimits_1 ,\mathop x\nolimits_2 ,...,\mathop x\nolimits_n ) $，则函数的标准差$ {\delta }_{y} $为：

当观测值的随机误差为正态分布时，式（12）中的标准差$ {\delta }_{y} $用极限误差$ {\delta }_{{\rm{lim}}y} $代替，可得函数的极限误差公式为：

上式中的$\dfrac{\partial f}{\partial {x}_{i}}$为各观测值的误差传递系数；$ {\delta }_{{{\rm{lim}}}xi} $为各观测值的极限误差；$ {\rho }_{ij} $为第i个观测值与第j个观测值之间的误差相关系数（费业泰，2017）。

2.   数据与计算

2.1   基本情况

长岛震群位于山东胶东半岛北部海域内（图2），自2017年2月14日再次形成震群活动，此次震群一直持续到2018年3月20日，之后该区地震活动逐渐减弱。此次震群余震丰富，很适合应用波速比变化揭示震源区介质的物性变化。通过震群的震中图（图3）可知，长岛震群的震源区最大长度在20 km左右，这相对于震群到台站的距离是较小的，因此可运用双差波速比法对长岛震群进行波速比变化分析。

图 2  长岛震群位置

Figure 2.  Location map of Changdao earthquake swarm

下载: 全尺寸图片 幻灯片

图 3  长岛震群M_L≥1.0地震震中分布（双差定位）

Figure 3.  Distributio of earthquake with M_L≥1.0 in Changdao earthquake swarm(Double-Difference Hypocenter Locations)

下载: 全尺寸图片 幻灯片

本文选取此次震群中M_L≥2.0地震296个，选取震中距40～150 km的14个台站，分别为LOK、YTA、LZH、DL2、LAY、WEH、WED、RSH、HAY、PENL、ZHY、MUP、LAIZ、BHC，可用Pg、Sg震相到时4 438个，经过2次差分后，一共得到纵横波双差数据481 855对。该方法充分利用不同地震事件不同台站的到时差信息，丰富了数据量。

分别使用0.1 s×0.1 s直方图和0.1 s×0.1 s网格对的481 855对双差数据进行统计（图4），结果显示，双差数据对主要集中在到时数据覆盖范围的中部区域，数据对最多的单元格包含10 717个数据对，占总双差数据对的2.22%。双差数据分布呈现二维正态分布置信椭圆的形状如图4（b）所示，长岛震群震源深度随时间的变化如图5所示。

图 4  长岛震群P、S震相双差数据对统计结果

Figure 4.  Statistical plots for data pairs of double differential P and S phases of the Changdao earthquake swarm

下载: 全尺寸图片 幻灯片

图 5  长岛震群震源深度随时间变化曲线（双差定位）

Figure 5.  Time variation curve of focal depth of Changdao earthquake swarm（Double-difference hypocenter locations）

下载: 全尺寸图片 幻灯片

2.2   双差波速比方法的应用前提及误差影响

本文使用的双差波速比法有几个前提条件，一是震群中地震之间的位置足够近；二是地震波速度在震群所在区域内相对稳定；三是P波、S波传播路径一致；四是地震震中到台站的距离不能太近。按照双差波速比计算原则，在经过两次差分消除传播路径和地震发震时刻的影响后，对波速比反演影响较大的只剩下P波、S波到时误差。双差数据对被约束在以波速比为斜率的直线附近，但有试验表明（贾漯昭等，2017），在震相到时数据误差较大时，拟合得到的波速比与真实值之间存在较大偏差。实际情况中，在震相噪声和误差的影响下，双差数据对在平面内的分布可能严重离散（Palo等，2016）。

波速比的计算受震相读取误差的影响较大，有研究人员认为这样求得的波速比不能反映真实波速比的变化。虽然人为读取震相的误差和观测系统的误差对波速比计算的影响是显著的，但它们又是不可避免的，因此很多研究人员多采用波形互相关技术改善震相到时的准确度，但波形互相关技术在震群地震的应用中也有不足，即震群震源区介质破碎、裂隙发育，不是所有的地震都能找到合适的模板地震进行相关（郑建常等，2019）。

此前的研究求解双差波速比，多使用最小二乘法对双差数据对进行直线拟合，但最小二乘法易受明显偏离异常点的影响，使拟合的直线偏离数据分布趋势，为此本文应用去除误差分布、基于相关系数的二维正态分布对双差波速比进行拟合。通过以上方法，可以较好地消除异常点的影响，达到稳健反演的目的。

2.3   正态分布的组合

根据现行区域台网编目震相评比标准，一般初至震相误差最大为0.2 s，续至震相误差不能大于0.5 s（贾漯昭等，2017）。本文为了突出数据规律，将纵横波到时的误差极限值分别设置为0.2、0.3 s。由式（13）可以计算出式（9）、式（10）的误差极限值分别为0.4、0.6 s。纵横波到时的误差可以认为是符合正态分布的随机变量。

因为真值为测量值与误差之差，所以可以得到$ {r}_{真}=\dfrac{{\delta }_{d}^{{\rm{S}}}-{E}_{{\rm{S}}}}{{\delta }_{d}^{{\rm{P}}}-{E}_{{\rm{P}}}} $，$ \delta _d^{\rm{P}} $ 和 $ \delta _d^{\rm{S}} $ 分别为P波、S波的双差到时数据，$ {E_{\rm{P}}} $和$ {E_{\rm{S}}} $分别为P波、S波的双差到时误差。

因为 $ \delta _d^{\rm{S}} $ 和 $ {E_{\rm{S}}} $ 相互独立，且$ \delta _d^{\rm{S}}{\text{～N}}({\mu _{\text{1}}}{\text{,}}\sigma _1^2) $、$ {E_{\rm{S}}}{\text{～N}}({\mu _{\text{2}}}{\text{,}}\sigma _2^2) $，则根据有限个相互独立的正态随机变量线性组合仍然服从正态分布，（$ \delta _d^{\rm{S}} - {E_{\rm{S}}} $）符合$(\delta _d^{\rm{S}} - {E_{\rm{S}}}){{\text{～N}}}({\mu _{\text{1}}} - {\mu _{\text{2}}}{\text{,}}\sigma _1^2 + \sigma _2^2)$，（$ \delta _d^{\rm{P}} - {E_{\rm{P}}} $）的正态分布同理可得。

由双差到时数据的数值拟合可知，S波双差数据 $ \delta _d^{\rm{S}} $ 符合${\delta }_{d}^{{\rm{S}}}\text{～N}（0.083\;1,{0.758}^{2}）$，S波双差数据误差$ {E_{\rm{S}}} $符合$ {E}_{s}\text{～N}（0,{0.6}^{2}） $，则$ (\delta _d^{\rm{S}} - {E_{\rm{S}}}){{\text{～N}}}(0.083\;1{\text{,}}0.758_{}^2 + 0.6_{}^2) $，即$(\delta _d^{\rm{S}} - {E_{\rm{S}}}){{\text{～N}}}(0.083\;1{\text{,}}0.966\;7_{}^2)$。同理可知P波双差数据 $ \delta _d^{\rm{P}} $ 与P波双差数据误差$ {E_{\rm{P}}} $之差符合$(\delta _d^{\rm{P}} - {E_{\rm{P}}}){{\text{～N}}}(0.039\;3{\text{,}}0.652\;3_{}^2)$。综上，（$ \delta _d^{\rm{S}} - {E_{\rm{S}}} $）和（$ \delta _d^{\rm{P}} - {E_{\rm{P}}} $）符合相关系数为 $ \rho $ 的二维正态分布。设（$ \delta _d^{\rm{S}} - {E_{\rm{S}}} $）的数据为$ x $，其均值为$ {\mu _x} $，标准差为$ {\sigma _x} $，（$ \delta _d^{\rm{P}} - {E_{\rm{P}}} $）的数据为$ y $，其均值为$ {\mu _y} $，标准差为$ {\sigma _y} $，则$ x,y $符合的二维正态分布函数为$f(x,y)=\dfrac{1}{2{\text{π}} {\sigma }_{x}{\sigma }_{y}\sqrt{1-{\rho }^{2}}}\cdot \mathrm{exp}\left(-\dfrac{1}{2(1-{\rho }^{2})}\left[\dfrac{(x-{\mu }_{x}{)}^{2}}{{\sigma }_{x}^{2}}-\right.\right. \left.\left.\dfrac{2\rho (x-{\mu }_{x})(y-{\mu }_{y})}{{\sigma }_{x}{\sigma }_{y}}+ \dfrac{(y-{\mu }_{y}{)}^{2}}{{\sigma }_{y}^{2}}\right]\right)$。

2.4   相关系数$ \rho $的确定

由前文可知双差波速比的理想真值$ {r}_{真}=\dfrac{{\delta }_{d}^{{\rm{S}}}-{E}_{{\rm{S}}}}{{\delta }_{d}^{{\rm{P}}}-{E}_{{\rm{P}}}} $，又根据双差到时数据 $ \delta _d^{\rm{S}} $ 和 $ \delta _d^{\rm{P}} $ 可知两者的相关系数为0.703，$ {E_{\rm{S}}} $和$ {E_{\rm{P}}} $是分别符合${{\rm{N}}}(0,{0.6}^{2})$和${{\rm{N}}}(0,{0.4}^{2})$的随机变量，那么（$ \delta _d^{\rm{S}} - {E_{\rm{S}}} $）和（$ \delta _d^{\rm{P}} - {E_{\rm{P}}} $）的相关系数最大值求解如下。

已知相关系数$ \rho = \dfrac{{\displaystyle\sum\limits_{i = 1}^n {({x_i} - \overline x )({y_i} - \overline y )} }}{{\sqrt {\displaystyle\sum\limits_{i = 1}^n {{{({x_i} - \overline x )}^2}} } \sqrt {\displaystyle\sum\limits_{i = 1}^n {{{({y_i} - \overline y )}^2}} } }} $，对于随机生成的1组变量，不论其排序如何，$ \rho $的分母是不变的，所以要求相关系数的最大值，只要求出$ \rho $计算式分子最大值即可。已知$ ab \leqslant \dfrac{{{a^2} + {b^2}}}{2} $，则$\displaystyle\sum\limits_{i = 1}^n {({x_i} - \overline x )({y_i} - \overline y )} \leqslant \displaystyle \sum\limits_{i = 1}^n {\dfrac{{{{({x_i} - \overline x )}^2}}}{2}} + \displaystyle \sum\limits_{i = 1}^n {\dfrac{{{{({y_i} - \overline y )}^2}}}{2}}$，代入本文数据得到$\displaystyle\sum\limits_{i = 1}^n {((\delta {{_d^{\rm{S}}}_i} - {E_{\rm{S}}}_i) - \overline x )((\delta {{_d^{\rm{P}}}_i} - {E_{\rm{P}}}_i) - \overline y )} \leqslant \displaystyle\sum\limits_{i = 1}^n {\dfrac{{{{((\delta {{_d^{\rm{S}}}_i} - \overline x ) - {E_{\rm{S}}}_i)}^2}}}{2}} + \sum\limits_{i = 1}^n {\dfrac{{{{((\delta {{_d^{\rm{P}}}_i} - \overline y ) - {E_{\rm{P}}}_i)}^2}}}{2}}$，随机生成1组符合正态分布的$ {E_{\rm{S}}} $，显然将$ (\delta_{di}^{\rm{S}} - \overline x ) $按照降序排列、$ {E_{\rm{S}}} $按照升序排列，可以得到$\displaystyle \sum\limits_{i = 1}^n {\dfrac{{{{((\delta {{_d^{\rm{S}}}_i} - \overline x ) - {E_{\rm{S}}}_i)}^2}}}{2}} $的最大值，同理可以得到$\displaystyle \sum\limits_{i = 1}^n {\dfrac{{{{((\delta {{_d^{\rm{P}}}_i} - \overline y ) - {E_{\rm{P}}}_i)}^2}}}{2}} $的最大值，此时得到的相关系数也是最大的。将此过程随机模拟10 000次，得到相关系数区间为[0.684 8,0.722 9]，取其平均值为0.703，与双差数据$ \delta _d^{\rm{S}} $和$ \delta _d^{\rm{P}} $的相关系数相同。由此可以得到（$ \delta _d^{\rm{S}} - {E_{\rm{S}}} $）和（$ \delta _d^{\rm{P}} - {E_{\rm{P}}} $）的相关系数最大值为双差到时数据$ \delta _d^{\rm{S}} $和$ \delta _d^{\rm{P}} $的相关系数。

2.5   双差波速比的计算

相关系数为$ \rho $的二维正态分布函数协方差矩阵为$\displaystyle \sum =\left[\begin{array}{cc}{\sigma }_{x}^{2}& \rho {\sigma }_{x}{\sigma }_{y}\\ \rho {\sigma }_{x}{\sigma }_{y}& {\sigma }_{y}^{2}\end{array}\right] $，代入本文数据得到协方差矩阵为$\displaystyle \sum =\left[\begin{array}{cc}0.4254& 0.4433\\ 0.4433& 0.9346\end{array}\right] $。

对上述协方差矩阵求特征向量，得到2个特征向量：$ {\left[-0.865 5 0.501\right]}^{\mathrm{T}} $和$ {\left[0.501 0.865 5\right]}^{\mathrm{T}} $。前者为以置信椭圆中心为原点的短轴向量（图6绿色箭头线所示），后者为以置信椭圆中心为原点的长轴向量（图6红色箭头线所示）。长轴与$ x $轴夹角的正切值为$ \dfrac{0.865\;5}{0.501}=1.727\; 5 $，这是长轴斜率。因为置信椭圆以长轴为对称轴，所以长轴两侧的对应数据对到长轴的距离是相等的，长轴上的数据对是整个椭圆内数据对的均值；因为（$ \delta _d^{\rm{S}} - {E_{\rm{S}}} $）为坐标竖轴，（$ \delta _d^{\rm{P}} - {E_{\rm{P}}} $）为坐标横轴，且$ {r}_{真}=\dfrac{{\delta }_{d}^{{\rm{S}}}-{E}_{{\rm{S}}}}{{\delta }_{d}^{{\rm{P}}}-{E}_{{\rm{P}}}} $，所以长轴斜率可以认为是双差波速比的拟合值。

图 6  二维正态分布置信椭圆及长短轴向量

Figure 6.  Two dimensional normal distribution confidence ellipses and major axis and minor axis vectors

下载: 全尺寸图片 幻灯片

2.6   反演数据统计分析

运用正态分布进行统计分析，其结果随着样本量的增加而趋于稳定，因此为保证有足够多的震相数据参与计算，选择等地震数滑动的方法计算长岛震群波速比随时间的变化。最终通过比较，选取60个地震作等地震数滑动计算，双差数据对个数变化曲线如图7所示，其中最小样本量为12 534个，最多数据量为29 927个，满足了正态分布拟合的要求。

图 7  长岛震群波速比分析结果

Figure 7.  Results of wave velocity ratio of Changdao earthquake swarm

下载: 全尺寸图片 幻灯片

相关系数 $ \rho $ 对二维正态分布形态有很大影响，由协方差矩阵可知相关系数直接参与波速比计算，因此相关系数是反演双差波速比过程中重要的一环。一般情况下，$ \rho \geqslant 0.6 $为强相关，0.4＜$ \rho $＜0.6为相关，0.2＜$ \rho $＜0.4为弱相关。相关系数越大，数据对越集中于置信椭圆长轴附近；相关系数越小，数据对越离散。由图7可知，相关系数最大值为0.884 5，最小值为0.302，51.3%的相关系数在0.4以上。

通过等地震数滑动方法计算得到此次震群波速比变化范围为[1.488 7，3.054 1]，前期计算的长岛震群整体波速比为1.727 5。

长岛震群双差波速比在95%置信区间内的误差范围如图8所示，在平均值以下的波速比误差中，最大误差为0.060 7，最小误差为0.023 1，误差＜0.04的占68.4%；在平均值以上的波速比误差中，最大误差为0.061 4，最小误差为0.023 3，误差＜0.04的占68.4%。

图 8  长岛震群双差波速比在95%置信区间内的误差范围

Figure 8.  Error range within 95% confidence intervals for double difference wave velocity ratios in Changdao earthquake swarm

下载: 全尺寸图片 幻灯片

3.   结果与分析

3.1   传统方法计算波速比分析

运用多台和达法计算了长岛震群301个M_L≥2.0地震的波速比，并给出了30点逐1点滑动的波速比曲线（图9红色曲线），结果显示该次震群的平均波速比为1.708 6。与本文方法相比，多台和达法计算的波速比未呈现出明显的趋势性变化，且整体值及均值均较小，除个别点外，整体处于均值附近2倍标准差的范围。在震群强烈地震活动的背景下，波速比显得波动较小且无趋势性，似乎与实际经验不符，多台和达法结果表明来自震源介质物性的变化被区域路径上的信息掩盖了。

图 9  多台和达法计算的长岛震群波速比曲线

Figure 9.  Plot of wave velocity ratio calculated by the multiple summation method

下载: 全尺寸图片 幻灯片

3.2   长岛震群波速比变化分析

本文应用基于相关系数的二维正态分布拟合得到的波速比曲线显示，长岛震群波速比最小值为1.488 7，最大值为3.054 1，除震群开始和结束阶段外，整体高于理论值1.732。2017年8月之前，震群主要分布在38°N以北的区域（称为“北部震群”），波速比除起始阶段外整体较高；2017年8月之后，震群主要分布在38°N以南区域（称为“南部震群”），由于本文采用60个地震滑动，且南部震群较北部震群地震数量较少，所以在南部震群开始阶段波速比仍呈现较高的值，但南部震群波速比整体上较北部震群低。

4.   讨论与结论

4.1   长岛震群震源区波速比变化的物理机制

岩石波速比是孔隙状态、流体饱和条件、有效压力和温度的函数（Domenico等，1984），且岩石物理性质，如岩性、矿物成分、胶结和固结状态也是影响波速比的重要因素（Wilkens等，1984）。本文认为在北部震群和南部震群（以38°N为界）存在不同程度的流体作用（北部震群较强，南部震群较弱），震群波速比变化与介质流体饱和程度、裂隙密度和状态有较大相关性。北部震群开始阶段地震密集发生，且震级较大，发生了3次4.0级以上地震，地震导致断层面附近介质破碎化，地下岩体裂隙增多，出现岩体孔隙度增大和流体不饱和状态，波速比呈现低值变化；随着地下流体沿岩体裂隙的渗透，岩体流体饱和度逐渐增大，波速比呈现升高趋势，从2017年3月下旬到4月上旬波速比总体增大；2017年4月8日发生了3.8级地震，导致地下岩体孔隙度继续增大，波速比有下降趋势，之后地下流体在前期裂隙贯通良好的情况下以较快速度渗透到刚形成的裂隙中，波速比又增大，并达到最大值3.054 1；之后震群地震稀疏，波速比有下降趋势，这与流体的再扩散有关系；2017年6月，地震移至北部震群西南侧，且整个6月小震密集发生，因此波速比小幅上升，直至6月18日震群发生4.0级地震后波速比下降；2017年8月，震群迁移至38°N以南区域，形成南部震群，从9月开始，波速比出现较大下降，这与南部震群地震较少且流体作用不强有关。

本文计算的长岛震群波速比为1.727 5，较理论值1.732稍低，等地震数滑动得到的波速比最大值超过3.0，除个别时段外，整体波速比相对较高，这与郑建常等（2019）得到的乳山震群波速比有较好的一致性，对此有几种可能的解释。Wang等（2012）试验表明，在流体饱和状态下，孔隙压力增大和孔隙密度增加均会使波速比超过2.0，且孔隙排列导致的介质各向异性对前两者引起的波速比变化有放大作用。另外，Nicholson等（1985）研究表明，波速比随深度的增加而减小；Catchings等（2014）研究发现，由于近地表的复杂性，浅部断层带往往具有高波速比。这意味着长岛震群的高波速比对应着震源的较浅分布，长岛震群震源深度集中在6～15 km（申金超等，2019），图5可较好地体现该结论。

4.2   结论

参考双差波速比方法，应用基于相关系数的二维正态分布拟合，得到了长岛震群整体波速比及等地震数计算的分步波速比。该方法充分利用了不同台站的纵横波到时数据，有效消除了地震波从震源区到台站传播路径效应的影响，同时该法在去除误差分布之后，拟合得到的波速比更趋近实际情况。

采用上述方法对2017年2月至2018年3月长岛震群波速比进行分析，选择震中距40～150 km范围内14个台站记录到的长岛震群296个M_L≥2.0地震震相到时数据，运用误差传播理论计算了经过2次差分后的纵横波到时差误差极限值，应用基于相关系数的二维正态分布拟合得到了长岛震群波速比为1.727 5，该值较理论值1.732稍小。使用等地震数（60个）进行滑动计算，得到了长岛震群2017年2月至2018年3月的波速比变化区间为[1.488 7, 3.054 1]。分析认为，波速比的波动变化与震群活动密切相关，且反映了介质的流体饱和度、裂隙密度和状态变化。相比于传统的多台和达法，本文方法得到的震源区波速比由于消除了区域地壳介质传播路径效应的影响，更能真实地反映震源区介质性质。