0 引言

波速比是地震学中最重要的运动学参数之一，它与地震波传播路径上介质的物理性状密切相关，反映了岩石受应力作用以及其他物理条件改变的影响.岩石物理和实验研究表明，波速比与岩石的岩性之间存在密切关系，岩石微观结构如裂隙、孔隙的形状、密度以及岩石含水的饱和程度等都对波速比起着重要影响(Domenico, 1984; 耿乃光等，1992；Christensen, 1996; Bhakta and Landrϕ, 2014; Catchings et al., 2014; Lin et al., 2015; Saito et al., 2016).波速比也是地震学领域重要的研究方向，地震学家很早就开始这方面的研究，并认为研究波速(比)的变化有望提供有关震源和地震孕育过程的某些信息(Scholz et al., 1973；冯德益，1975；冯德益等，1976).近些年的研究显示在大地震发生前后观测到明显的速度变化(~0.06%，Brenguier et al., 2008)，认为是可能的地震前兆(Xu and Song, 2009).国内的研究如：嘉世旭和刘昌铨(1996)研究发现邢台地震震源区存在高波速比异常；王林瑛和李艳娥等分别利用单台、多台及直达波视速度等多种方法研究了汶川地震、芦山地震前后波速比的变化特征，发现强震前均存在不同程度的异常变化(王林瑛等，2014；李艳娥等, 2014, 2016).

常用计算波速比的方法是利用多台或多震的方法计算平均波速比(黎明晓和张晓东，2004；王林瑛等，2014；李艳娥等, 2014, 2016)，或者进行波速比和其他参数的联合求解(Smith et al., 1983；Koch, 1992).平均波速比是包括震源区和射线传播路径介质效应的综合反映，通常被解释为地壳结构和构造性质的差异(Walck, 1988; Hauksson and Haase, 1997; 李志伟等，2009；危自根等，2011；危自根和陈凌，2012；张洪艳等，2015).这里我们考虑面向震群活动的特殊情况：震群是局部地区在相对较短时间范围内发生的密集破裂，相比于大震的孕育过程，震群活动时间持续较短(数天、数月，个别情况能持续数年)，由于震源区频繁的小震活动，震源区介质持续不断地摩擦、破裂、破碎，因此其波速比的变化也必然是显著的.

在一般的观测条件下，相对震源到台站的传播路径，震群的震源区尺度很小；因此即使震源区出现较大的波速比变化，在利用台站震相到时计算波速比时，最终计算结果中很难反映出源区的波速比变化(例如，震源区尺度10 km，震源到台站的传播路径100 km，即使震源区波速比变化0.1，而最终在台站观测到的波速比变化可能只有0.01，因此在实际情况下使用平均波速比很难识别出震源区内波速比的变化)；此时，区域地壳的速度结构中介质的横向不均匀性，是影响台站震相到时和波速比值的重要因素.在震群活动的情况下，震源区介质在持续的破碎和错动，从而介质性质包括波速和波速比等出现快速变化，而传播路径上的介质性质在较短的时间尺度上可以认为是相对稳定的，因此如何消除传播路径效应的影响而获取来自震源区的介质性质变化是本文研究的重点所在.

借鉴双差定位的思想，Lin和Shearer(2007)首次讨论和检验了应用两次差分计算震源区波速比的可行性，随后Lin与她的合作者应用该方法研究并讨论了震源区附近的波速比异常，认为可能波速比的异常变化与地下流体的作用有关(Lin and Shearer, 2009; Lin, 2013; Zhang and Lin, 2014; Lin et al., 2015)；Palo等(2016)严格检查了Lin等提出的方法，通过不同模型下的理论分析和数值检验，发现Lin的方法存在一定的限制，在某些情况下最终的估计值可能出现系统偏差，并讨论了消除偏差的方法.Dahm和Fisher(2014)采用两步反演技术，先网格搜索计算区域平均波速比，然后使用LMS(Least Median of Squares)方法基于双差后的到时数据估计震源区的波速比和泊松比；基于他们的方法，Bachura和Fischer(2016)发现震源区波速比值与震群活动的时间演化特征相一致，可能表明了震群活动期间震源区介质性质的变化；基于同样的方法，Novotny等(2016)利用震中距小于25 km的23个台站，选择了158个M_L1.5-3.8的事件，研究了西波希米亚2008年的一次震群活动的波速比变化.此外，Kato等(2010)、Akuhara等(2013)、Doi等(2013)研究表明波速比异常在研究地震机理、识别地震活动模式以及分析断层或地下介质性质等方面都有重要的作用.

本文在前人研究的基础上，基于误差分布的思想和概率统计分析的方法，借鉴双差定位程序中两次差分的技术，实现了一种新的求解震群活动震源区的双差波速比方法，并将该方法应用于2013年活动至今的乳山震群；基于震源区波速比以及震群活动随时间的变化，对乳山震群的源区破裂性质进行了讨论.

1 理论与方法 1.1 震源波速比的双差求解方法

参考Lin和Shearer(2007)及Dahm和Fischer(2014)的研究，设震群中共发生了m次事件被区域范围内n个台站记录到，第i(i=1, …, m)次事件在台站j(j=1, …, n)的纵横波观测到时分别为T_ij^P，T_ij^S，设发震时刻为t_i⁰，从震源到台站纵横走时分别记为t_ij^P和t_ij^S，则

(1)

任选两个台站组成台站对，则对每次事件，共有C_n²个台站对的组合，台站对组合记为k(k=1, …, C_n²).不失一般性，考虑k=(j, j+1)的台站对，定义

(2)

注意该公式中消除了事件的发震时刻，即台站之间到时差仅和走时差有关.

任选震群中的两次事件组成事件对，则共有C_m²个组合，事件对记为h(h=1, …, C_m²).考虑如图 1所示震群中的事件对(i, i+1)，设震源区的纵横波速度为α₀, β₀，波速比为γ₀；而在源区外围射线路径上的纵横波速度和波速比分别为α, β, γ，则公式(1)中的台站到时可以写为

(3)

定义

(4)

在一般情况下，震群破裂区范围要远小于台站的震中距，例如一次6级地震的破裂半径可能约10 km，而区域范围内大多数台站的震中距可能大于100 km.因此震群中两次事件在震源区之外到台站的路径非常接近，近乎相同，并且在震群活动期间，震源区外路径上的介质性质可视为稳定，α, β, γ均保持不变，因此可以认为公式中t_ij^P外－t_{i+1, j}^P外≈0、t_ij^S外－t_{i+1, j}^S外≈0.则某一台站不同事件的到时差与发震时刻以及震群事件在震源区内的相对位置有关.

将公式(4)中的台站到时替换为公式(2)给出的台站到时差，即对初始观测到时进行两次差分.同样不失一般性，此处选择h=(i, i+1);k=(j, j+1)，定义

(5)

设震源破碎区内的波速比为γ₀，则显然

(6)

如此，我们经过观测到时的两次差分，可以得到震源破碎区的波速比.与传统的和达线等方法相比，两次差分的方法更充分地利用了台站到时信息：考虑(T^P, T^S)数据集，传统方法中可用的数据量为2·m·n条，而经过两次差分后可用的数据量则为2·C_m²·C_n²条，在独立测量的情况下，30个地震8个台站，可用数据就超过2万个(20880).

1.2 基于正态分布的误差加权

1.1节介绍了源区波速比的双差求解方法(Lin, 2007, Dahm和Fischer，2014).该方法有效依赖于几个基本前提：(a)震群中的事件之间的距离足够近；(b)地震波传播速度在区域范围内是相对的常量；(c)P波、S波路径一致；(d)事件到台站的距离不能太近.在理想的情况下，经过两次差分消除路径效应和发震时刻的差异后，对反演波速比存在显著影响的仅仅是观测到时的误差，数据对被约束在斜率为波速比r的直线附近，但有数值实验显示(贾漯昭等，2017)，震相误差较大时，拟合得到的波速比可能与真实值严重偏离.在实际情况下，震相读取的噪声与误差可能使得数据对在平面内严重离散(Palo et al., 2016).

此前研究者在求解双差波速比时，都是通过对平面上的数据对进行直线拟合来得到结果，而拟合过程通常使用的是最小二乘法.最小二乘法的缺点是容易受明显偏离的异常点影响，导致拟合直线偏离真正的数据变化趋势，为此基于误差的正态分布函数对观测点到拟合解的距离进行加权.首先设观测误差满足正态分布，对于离散分布的观测点x(i)，定义权重：

(7)

对比正态分布的概率密度函数

(8)

注意该权重函数与之在形式上基本相同，因此该权重是基于正态分布的误差加权.在我们的反演中，定义公式(7)中的τ为台站到时的误差范围，在100 Hz采样率的情况下，τ取0.01 s、0.1 s、0.3 s时的权重函数形态分别如图 2所示.

通过不同的误差权重，可以明显消除偏离异常点的影响，从而达到稳健反演的目的.

1.3 误差分析与二维高斯分布

波速比的计算容易受到震相读取误差的影响，因此有研究人员对波速(比)的变化持有怀疑态度(Geller, 1997)，认为难以区分真实的异常变化和误差效应.观测系统误差和人为读数偏差对波速比的影响显而易见(顾瑾平和盛国英，1991；蔡静观，2000)，但这些误差又是不可避免的.为了最大可能消除误差的影响，多数研究者采用的方法是使用波形互相关来改善震相到时的准确度(Lin and Shearer, 2009; Lin, 2013; Zhang and Lin, 2014; Lin et al., 2015; Bachura and Fischer, 2016).但该方法仍然存在很大的局限性，尤其是应用到震群活动时：震群活动期间，震源区介质高度破碎，细小的断层、裂纹四处发展，不是所有的小震都能找到合适的模板地震来进行相关.

从数理统计的角度，大量次数的测量误差服从一定的统计规律.对于波速比计算的误差，蔡静观(2000)曾讨论过其统计分布，也有一些研究学者通过不同的统计手段来估计或消除波速比计算中震相误差的影响(如，Bootstrap方法，Koch, 1992; Jack knife方法，Wu and Lees ,1999, ).本文基于误差分布的思想和概率统计分析，提出一种新的思路实现双差波速比的稳健反演，介绍如下：

理论上，随机误差满足正态分布(高斯分布).对单个的随机变量x，其概率密度函数如公式(8).那么，对于彼此独立的变量x和y，其测量误差的分布是典型的椭圆形二维正态分布(图 3a)，但当变量x和y存在正相关关系时误差分布则呈现类似图 3b的形态.

设变量x和y的联合分布为二维正态分布，假设其相关系数为ρ，则两个变量的联合概率密度分布函数为

(9)

其中，σ为标准差，并且

(10)

μ为数学期望，Σ为协方差矩阵.

对于满足二维正态分布的数据，通过数值拟合，可以得到联合概率密度分布函数f(x, y)，并给出概率密度等值线图，即不同置信水平的误差椭圆(如图 3c).特别地，对于标准化(μ_x=μ_y=0)的二维正态分布，如图 3c所示，其误差椭圆有三个重要的参数：倾角θ、长轴u和短轴v，且有以下关系：

(11)

(12)

(13)

考虑T^P, T^S的到时误差(ε(T^P), ε(T^S))的分布，则呈现出如图 3b的椭圆形态的分布；两个到时的误差过程可以认为相互独立，则(ε(T^P), ε(T^S))满足二维高斯分布，且误差的概率密度分布函数呈现如图 3c的以原点为中心的误差椭圆，更进一步，由于纵横坐标轴分别代表t^P, t^S，由于t_P=γt_S，显然误差椭圆的长轴的斜率即等于波速比，即γ=tan(θ).如前所述，可视为等于ε(T^P)，，所以我们在对观测到时进行两次差分的基础上可以利用二维高斯分布的概率密度函数拟合来估计震源区波速比.

2 数据与资料 2.1 基本情况

乳山震群自2013年10月1日开始活动，截至2016年中，共记录到地震事件超过12000次，其中M_L≥3.0有31次，最大为2015年5月22日M4.6级地震(图 4).这次震群的余震事件丰富，持续时间长，因此很适于通过波速比来研究震源区介质物性的变化.

乳山震群震中距200 km范围内有山东测震台网的固定台19个(图 4a)，考虑到观测数据的连续性以及信噪比等情况，我们选择了震中距150 km以内的HAY、WED、LAY、LZH、YTA、WEH、RCH、ZHY等8个台的资料进行分析，RSH台震中过近(12 km)，因此没有选用.

2.2 数据初步分析

为保证Pg、Sg震相的数量和可靠性，选择M_L≥1.0级的地震，共计350次(见图 4b).使用双差方法重新定位后，可以发现乳山震群集中在很小的破裂范围内(图 4c)，破裂区的长度不超过5 km(郑建常等，2015)，相对于胶东半岛地区台站震中距是很小的，因此可以使用双差波速比的方法来求取震源区的介质波速比变化.

在本文选用的8个台站中，可用Pg、Sg震相到时共4622个；经过两次差分以后，可用数据348637个，共338917个，极大地丰富了数据量，并且充分利用了不同事件不同台站之间的到时差信息.

使用0.01 s×0.01 s网格进行统计，可用数据对322980个，空间分布见图 5.方便起见，在后文中，我们使用ddtp和ddts分别代替和.根据统计结果，数据对最多的单元格2430个，约占总数据对数量的0.752%(图 5a给出了0.1 s网格的统计结果)；图 5b给出了0.01 s网格的统计密度，可以看出震相数据对的空间分布呈现出很好的误差分布图像.

使用本文提出的方法进行二维高斯分布拟合，结果见图 6：图中等值线为拟合函数概率密度分布，误差椭圆的长轴，即图中红色实线，其斜率即为计算得到的震源区波速比值(r=1.8935).

3 结果与分析 3.1 传统平均波速比变化分析

为了便于与传统的波速比计算结果进行比较，我们给出了使用单震多台法给出的乳山震群350次事件的平均波速比(见图 7).结果显示：在乳山震群活动期间，平均波速比未出现明显的趋势性变化，除了个别点外，整体都处于均值1.708附近两倍标准差范围内.如果从震源的角度分析，在总计一万多次、日频次超过100次的强烈活动下，这一结果似乎与物理经验不符.因此正如我们前面所分析的，来自震源的变化被区域路径上的信息掩盖掉了.

3.2 反演质量分析

基于误差的统计分析，其结果随着样本量的增加而趋于稳定.因此，为保证足够的震相数据量，我们选择等地震数滑动的方法来计算乳山震群震源区波速比随时间的变化.虽然在理论上，30个地震8个台站情况下可用数据就超过2万个，足以满足稳定性要求，但实际情况中，往往出现部分台站的报告中震相缺失的现象，并且也有个别台站由于种种原因导致没有记录.最终通过比较，我们将用于滑动计算的地震数选为60个.

图 8给出了60个地震初次滑动的计算结果.图 8c是计算过程中使用的震相双差数据对{ddtp, ddts}的个数变化曲线，震群活动开始后的大约1年时间，可用数据对在14000个左右；随后有些下降，但整体仍在11000个以上，最少为10813个，满足了用于概率分布拟合的数据量.

由公式(9)可见，相关系数ρ对二维高斯分布的形态具有重要影响，因而ρ也同时影响了最终双差波速比的反演质量，如果ρ接近0或者小于0，即ddtp和ddts不相关甚至反相关，则说明原始数据中可能存在问题，反演得到的震源区波速比值不可靠.图 8b给出了乳山震群计算过程中数据ddtp和ddts的相关系数变化曲线，根据一般认识：相关系数ρ≥0.6为强相关，ρ∈[0.4, 0.6]为相关，ρ∈[0.2, 0.4]为弱相关，计算结果显示：所有数据点的ρ均大于0.3；63.1%的相关系数大于0.4，近1/5的计算结果达到了强相关的程度，说明计算结果的可靠性较高.相关系数的最低值ρ=0.305，仍然是可以接受的结果.

3.3 波速比变化分析

反演得到的波速比变化曲线(图 8a)显示：乳山震群震源区波速比值在[1.647, 2.921]之间变化，除了个别时段外，整体高于1.732；从震群活动开始到其后一年多的时间内，波速比值呈现上升趋势，而后回落，到2015年底再次出现上升.震源区波速比的变化与震群活动过程密切相关.震群活动的初期，波速比的持续上升，可能反映了震源区持续的破碎和地下流体逐渐趋于饱和，小震重新定位的结果也显示，在此期间，小震分布逐渐向外围扩散(图 4c).而从震群活动第400至600天，也即到震群最大事件，2015年5月22日M4.6地震前，震源区波速比存在一个下降的趋势，或许和此次最大地震的孕育有关.另外，Zheng等(2017)研究发现，从2014年5月到2015年4月，在震源区存在一个M_L≥2.0级地震的“空区”；这一时段内，M_L≥2.0的余震几乎都发生在空区外围；这一现象恰好与我们计算的震源区波速比下降变化过程相一致.

4 讨论与结论 4.1 乳山震群源区波速比变化的物理机制

波速比是理解岩石微结构及其岩性的重要工具.岩石的波速比是孔隙和流体条件、有效压力和温度的函数(Domenico, 1984)；此外，岩石的物性特征如岩性、矿物组份、胶结和固结程度等也是决定波速比的重要因素(Wilkens, 1984).Nicholson和Simpson(1985)研究发现，波速比随深度的增加而减小.按照传统的观点，波速比下降是介质微裂隙增多导致体积膨胀、岩石中流体变得不饱和所致.Christensen(1996)研究发现不同温度和压力条件下，V_P/V_S变化约在1.478~2.119范围内；一般情况下，压力增大可能使波速比升高，但变化幅度并不大(0.0077每100 MPa).因此，Lin等(2016)认为波速比值的变化对温度和压力的变化不敏感，而对存在流体与否更敏感.其他的研究如：Wang等(2012)通过岩石实验和数值模拟发现高孔隙压力和破裂的各向异性的累积效应是出现高波速比的主要原因；Bhakta和Landro(2014)研究认为，在岩性相同和外在条件变化不大的情况下，通常波速比强烈依赖于孔隙度和流体饱和条件，并且孔隙形态的变化对波速比也有显著影响.

如果我们假设乳山震群震源区存在流体作用，则我们计算得到的震源区波速比变化可以很好地匹配乳山震群的活动过程.例如：在乳山震群活动初期，断层错动造成断层面附近介质破碎，岩石破碎化和水不饱和，导致了震源区波速比的低值(早期余震)；随着震群活动的持续(晚期余震)，地下流体向震源区裂隙渗透，震源区介质含水趋于饱和，波速比显示高值变化(周惠兰等，1982；Bachura and Fischer, 2016).对应到乳山震群，活动开始一个月的时间内，波速比在理论值1.732附近，之后由于介质的破裂，岩石孔隙度的发展使得波速比降低(Saito et al., 2016)，属于所谓的“早期余震”，而100天以后至2014年9月之前，震源区波速比整体呈现上升趋势，与理论分析的高压流体逐渐渗透到饱和过程相吻合；而2014年10月至2015年5月的波速比下降恢复，应当与震源区流体在集中饱和后，随着流体的再次扩散，由饱和到与周边达到平衡过程有关，另外这一变化可能还与应力作用有关，因为随后发生了震群的最大地震，并且最大地震后，震源区波速比又再次上升到高值.

精定位的结果显示，乳山震群的余震集中在一个狭小的空间范围内(郑建常等，2015)，可以认为乳山震群活动过程中，震源区的岩性、压力、温度等条件变化不大，因此我们计算得到的震源区波速比的变化，更主要地应该是岩石孔隙度和饱和程度的反映.Gritto和Jarpe(2014)研究地热库区的震群波速比与注水水压的关系时发现，波速比与饱和程度正相关；钟羽云等(2015b)研究珊溪水库震群的波速比变化发现，珊溪水库震中区岩石始终处于接近水饱和的饱水状态，波速比的变化实质上反映了震中区岩石破碎后孔隙度增大(饱和度减小)到饱和度增大的变化.我们对乳山震群的分析表明，波速比的变化与前面讨论的震群活动过程有很好的一致性.

乳山震群震源区波速比最高值超过2.5(虽然相关系数不高)，整体也相对较高；对此，有多种可能的解释.Wang等(2012)的岩石实验显示，在液体饱和情况下，孔隙压力增大和孔隙密度的增加都会使得岩样的波速比值超过2.0，并且数值计算显示，孔隙排列导致的各向异性对前二者引起的波速比变化有放大作用.Audet等(2009)对俯冲洋壳的研究以及Bezacier等(2010)对蛇纹岩的试验显示，由裂纹或者矿物结晶体的优势方向导致的各向异性会造成很高的波速比；而乳山震群震源区恰好是残存俯冲洋壳蛇纹岩分布的区域.另外，Catchings等(2014)研究发现，浅部断层带，由于近地表复杂性，往往具有高波速比的特征；或许相对较高的波速比也意味着乳山震群的震源区较浅，与我们使用现场台阵数据和深度震相得到的结果相吻合(郑建常等，2015).

4.2 结论

参考Lin和Shearer(2007)和Dahm和Fischer(2014)的双差波速比方法，基于误差分布和概率统计分析，通过对两次差分后的震相残差数据对进行二维高斯分布拟合，发展了一种稳健的估计震群震源区波速比的方法.该方法充分利用了不同台站Pg、Sg到时差的所有信息，有效消除了震源区到台站传播路径效应的影响.在误差不可避免的情况下，鉴于对双差后的数据直接线性拟合计算波速比的方法可能存在系统性偏差(Palo et al., 2016; 贾漯昭等，2017)，在大样本量的前提下，我们充分考虑了误差分布和误差加权，使得该方法显得更加稳定.

将该方法应用到2013年起活动至今的乳山震群.选择震中距150 km以内的8个台记录到的乳山震群350次M_L≥1.0级的震相记录，使用0.01 s×0.01 s网格对两次差分后的震相对数据进行统计，利用二维高斯分布拟合，得到震源区平均波速比为1.8935；相对较高的波速比值或许意味着乳山震群活动过程中震源区存在流体作用.使用等地震个数(60)进行滑动计算，结果显示：乳山震群震源区波速比值在[1.647, 2.921]之间变化；震源区波速比的变化与震群活动过程密切相关，震源活动初期，随着破碎区的扩散，震源区波速比升高；波速比的波动变化反映了序列活动的阶段性特征.乳山震群的应用结果显示，相对于传统的平均波速比，本文方法得到的震源波速比，由于消除了区域地壳介质的传播路径效应，更能真实地反映震源区介质的性质.

致谢  感谢审稿专家的意见和建议.本文使用的震相数据来自山东省地震台网中心提供的观测报告，对此表示感谢.