随着观测精度不断提高，利用区域数字地震资料分析中强地震前后波速比变化特征的研究和震例逐渐增多^[1]。单台多震和达法采用固定台站记录的震相数据，参与计算的地震方位和范围可以根据需要进行选择，提高了空间分辨率。孕震区附近台站单台波速比异常范围与孕震区域有一定的相关性^[2]，如2010年河津4.8级地震前临汾台波速比出现持续2~3 a的异常，地震台站的方位与孕震区范围明显相关^[3]。唐山地区自2000年以来地震活动水平在4级左右，近年来连续发生多次有感地震，王亚茹等^[4]利用多台多震和达法研究认为，唐山地区平均波速比是具有一定预报效能的中短期指标。

本文搜集河北数字地震台网2008~2018年的地震震相数据，利用单台多震和达法和多台多震和达法对唐山地区不同台站的波速比进行研究，结合唐山地区的地震活动特征讨论不同台站波速比变化与地震的关系，探讨唐山地区显著性地震发生前后波速比的变化特征，为唐山地区中等地震的预报提供一定的参考。

1 计算原理

假定震源到地表介质为理想的均匀弹性体，利用一组地震多个台站的P波走时和P、S波到时差数据(多台多震和达法)或者对一个固定的台站记录到的不同地震的P波走时和P、S波到时差数据(单台多震和达法)拟合直线求得斜率，再加1即为震源区到台站间介质的平均波速比，主要用来研究上地壳的介质物性。单台多震和达法固定了台站的位置，其结果为多次地震震源至某个固定台站的一定范围内的平均波速比，该方法可以根据需要选择不同方位和范围的地震参与计算，有效避免地震分布不均的影响，增加空间分辨能力，在大量中小地震频发和台站分布相对密集的地区应用广泛。

波速比根据天然地震的P波与S波到时差Δt和P波到时的和达图的线性关系确定：

$ \frac{{{V_{\rm{P}}}}}{{{V_{\rm{S}}}}} = 1 + \frac{{n\sum\limits_{i = 2}^n \Delta t_i^2 - {{\left( {\sum\limits_{i = 1}^n \Delta {t_i}} \right)}^2}}}{{n\sum\limits_{i = 1}^n \Delta {t_i}{t_{{{\rm{P}}_i}}} - \sum\limits_{i = 1}^n {{t_{{{\rm{P}}_i}}}} \sum\limits_{i = 1}^n \Delta {t_i}}} $

(1)

$ R = \frac{{\sum\limits_{i = 1}^n {\left( {{t_{{{\rm{P}}_i}}} - {{\bar t}_{\rm{P}}}} \right)} - \left( {\Delta {t_i} - \Delta {{\bar t}_i}} \right)}}{{\sqrt {\sum\limits_{i = 1}^n {{{\left( {{t_{{{\rm{P}}_i}}} - {{\bar t}_{\rm{P}}}} \right)}^2}} \sum\limits_{i = 1}^n {{{\left( {\Delta {t_i} - \Delta {{\bar t}_i}} \right)}^2}} } }} $

(2)

式中，t_Pi为P波走时，t_Si为S波走时，Δt_i=t_Si-t_Pi，i=1, 2，…，n，其中，在多台单震和达法中n为每个台站到时的数据个数，而在单台多震和达法中n为每个地震到时的数据个数，R为相关系数。波速比计算误差γ为：

$ \begin{array}{c} \gamma = {\left( {\frac{{{V_{\rm{P}}}}}{{{V_{\rm{S}}}}} - 1} \right)^2} \cdot \\ \sqrt {\left| {\frac{{\sum\limits_{i = 1}^n {{{\left( {\delta {t^\prime }_{{{\rm{P}}_i}}} \right)}^2}} }}{{(n - 2)[n\sum\limits_{i = 2}^n \Delta t_i^2 - {{(\sum\limits_{i = 2}^n \Delta {t_i})}^2}]}}} \right|} \end{array} $

(3)

2 资料选取

唐山块陷形态上为NEE走向的菱形块体，四周被活动断裂包围：南界是宁河-昌黎深断裂，北界是丰台-野鸡坨断裂，二者走向均为NEE；东界是滦县-乐亭断裂，西界是蓟运河深断裂，二者走向均为NW。菱形块体的基底为前古生代的变质岩系，盖层为古生代的沉积岩系，中生代地层缺失，新生代沉积了大量的第四系松散层，其沉积厚度向南逐渐增加到800 m以上，显示明显的S倾特征。在块体的中央存在发育于古生代盖层中的唐山断裂带，走向NW，由陡河断裂、唐山-巍山-长山南坡断裂和唐山-古冶断裂3条平行的断裂组成^[5]。唐山菱形块体是一个强烈活动的断块，1976-07-28唐山7.8级地震发生在唐山断裂带上，2000年之后该区地震活动水平维持在4级左右。相关研究^[6-7]表明，唐山地区地壳厚度在28~36 km，平均厚度约33 km，区内横向波速比无显著变化。

研究区域(39°~41°N，117.5°~119.5°E)位于首都圈东部(图 1)，共设有16个测震台站，台站分布密集且较为合理，地震监测水平较高。河北地震统一速报目录显示，2008年以来，研究区内共发生M_L4.0以上地震9次(表 1)，中小地震活动频繁，积累了大量高精度的震相数据，为波速比的计算提供了数据基础。选择研究区内连续性较好的16个台站2008~2018年的震相数据，分别采用单台多震和达法和多台多震和达法计算波速比，结合唐山地区的几次显著地震活动对波速比的变化特征进行分析。

3 波速比计算结果分析 3.1 单台多震波速比特征分析

选择研究区内台站周边100 km范围内至少4个地震的震相数据，利用单台多震和达法进行计算，筛选相关系数R≥0.98、计算误差γ≤0.05且震级M_L≥1.0的地震的波速比结果进行分析。由于地震分布的时空不均匀性，导致某些台站周边可以计算波速比的地震较少，数据连续性较差，因此挑选其中连续性较好的台站分析其时序变化特征(图 2)。图 2中，左边为地震射线图，红色圆表示台站记录到的符合上述筛选条件的地震，蓝色三角形表示台站，黄色五角星表示研究区6次M_L≥4.5地震；右边为所选台站的波速比时序曲线，曲线上部标注了台站名称和滑动步长，如“DOH-50”(图 2(a))代表DOH台波速比50点滑动均值。除LUT和ZUH台因数据量较少分别采用20点和10点滑动外，其余台站单台波速比结果均采用50点的滑动窗口进行平滑平均。黑色虚线为标准差，蓝色曲线为滑动均值曲线，灰色竖线为每个地震的计算误差。下面结合唐山地区的中等地震活动分析单台波速比的变化特征。

2004~2014年唐山地区M_L≥4.0地震显示出2 a左右的发震规律，同时兼有双震特征^[8]。选取2008~2018年的震相数据进行分析，2008-03-11卢龙M_L4.4地震因震前波速比数据较少，不作研究。2010-03-06和2010-04-09研究区分别发生滦县M_L4.8和丰南M_L4.6地震，几乎所有台站的地震射线都较好地覆盖了滦县M_L4.8地震的震中，DOH、LUT、QIX、FTZ和XAZ台的地震射线基本覆盖丰南M_L4.6地震。震前FTZ和DOH台均出现明显的低值下降，低值持续1 a左右后缓慢回升。LUT和XAZ台周边的地震相对稀疏，导致波速比数据连续性不佳，但震前也出现低值回升的趋势。出现异常的这4个台站中，DOH台位于唐山断裂带附近菱形块体的中间部位，距离滦县M_L4.8地震最近；FTZ、LUT和XAZ台位于蓟运河断裂附近，大体位于菱形块体的南部区域，其地震射线均较好地覆盖了滦县M_L4.8和丰南M_L4.6地震，这两次地震也分别发生在菱形块体的中部和南部，这表明异常台站与地震具有一定的相关性。震后DOH、FTZ和LUT台的波速比从高值快速回落。研究表明，大震后由于受震源区应力释放、介质产生破裂及液体流入等因素的影响，震源区地震波速会出现同震或者震后大幅度快速下降的现象^[9]，而岩石实验研究也表明^[10]，地下岩石的波速和波速比随时间的变化跟其周围应力变化存在一定的相关性，这可能表明，地震发生后台站附近的应力累积有所释放，应力水平降低。2012-05-28唐山M_L5.1地震为研究区近年来最为显著的地震事件，距2010-03-06滦县M_L4.8地震仅5 km，距最近的DOH台约16 km。该次地震位于唐山-滦县断裂带附近，后天津宝坻又发生M_L4.5地震，距离其最近的台站为FTZ。这两次地震发生前，DOH、LUT、QIL、QIX、FTZ和XAZ台在震前1.5 a左右出现同步下降，以DOH台异常最为显著，持续时间1 a左右，随后慢慢回升到均值线附近并小幅度波动直至地震发生，表现为下降-低值-恢复-发震的典型波速比异常形态。

2014-10-14滦县M_L4.0地震发生于该区2 a准周期的时间节点上，相较于前两次准周期活动，本次地震序列的强度和频度都较低。随后，唐山地区的地震活动出现新的变化。2015-09-14河北卢龙M_L4.7地震打破该区2 a准周期的发震规律，1 a后在2016-09-10唐山又发生M_L4.6地震。其中，滦县M_L4.0和卢龙M_L4.7地震均位于滦县-乐亭断裂带附近，靠近菱形块体的北部；唐山M_L4.6地震发生于唐山断裂带附近，位于菱形块体的中间。DOH、CLI、BDH、TLK、QIX和QIL台的地震射线较好地覆盖了这3次地震，但是卢龙M_L4.7地震前并无明显的低值异常，这可能与2014-10-14滦县M_L4.0地震缓解唐山地区的应力累积有一定的关系。在唐山M_L4.6地震发生前，CLI、BDH、TLK、QIL台波速比出现低值异常，其中BDH台低值持续近2 a，CLI、TLK、QIL台则在2016年初出现同步下降，持续0.5 a左右后又逐步回升，随后发震。这4个台站主要位于菱形块体的东北部，也较好地对应了这3次地震的位置。目前唐山地区M_L≥4.0地震已经平静了1.5 a，具有发生4级地震的背景。DOH台自2017年出现低值回升，CLI台在2017年底也出现小幅度的低值，目前已回升到均值线附近，其他台站无明显异常。综上所述，唐山地区单台波速比在中等地震发生前存在一定的异常变化，异常持续时间多为1 a左右，异常台站的方位与地震具有一定的相关性。

3.2 多台多震波速比特征分析

为了进行对比分析，选择至少4个P波、S波到时差Δt≤12 s的台站的地震震相数据，利用多台多震和达法计算波速比，筛选相关系数R≥0.98、计算误差γ≤0.05、震级M_L≥1.0的地震的波速比结果进行分析。图 3(a)为地震射线图，其中黄色五角星代表 2008年以来6次震级M_L≥4.5的地震，蓝色三角形代表台站，红色圆代表波速比符合相关系数R≥0.98、计算误差γ≤0.05的地震，黑色射线为地震射线，基本覆盖了研究区的6次震级M_L≥4.5的地震。图 3(b)为多台波速比时序曲线，为了对比分析，采用与单台波速比一致的限定范围(1.55~1.85)，其中黑色实线为波速比均值(1.71)，黑色虚线为标准差，蓝色曲线为50点滑动均值曲线，灰色竖线为每个地震的计算误差。由图可见，多台波速比结果较为平稳，6次M_L≥4.5的地震发生前均未发现明显的低值异常。为了显示更多的细节，将多台波速比滑动曲线限定范围缩小为1.68~1.75(图 3(c))，结果也呈现出明显波动。王亚茹等^[4]研究认为，唐山地区M_L≥1.5地震多台波速比低值异常持续3个月后，该区0.5 a内存在发生M_L≥4.0地震的可能。

将地壳视为理想弹性介质时，其泊松比为0.25，对应的波速比为1.73。唐山地区10个台站的单台波速比均值最大为1.72，最小为1.68，集中分布在1.71~1.72(7个台)，多台波速比均值为1.71，总体来说波速比波动范围不大，侧面反映该区的地壳各项异常并不显著。

3.3 讨论

王林瑛等^[10]采用多台法研究汶川地震前后波速比变化特征认为，地震前波速比低值异常持续时间仅为2 a多，且幅度较小；李艳娥等^[2]利用单台法对汶川地震波速比进行研究认为，震前4个台站出现长达7 a左右的中长期波速比异常，且与汶川地震孕震范围大体一致，表明汶川地震前，附近台站单台波速比异常持续的时间要远长于多台波速比，且变化幅度更大。同样，文中多台波速比相对稳定，表明该区地下介质的物性变化较小，而单台波速比的波动幅度更大，异常持续时间也更长，且异常台站的方位与地震具有一定的相关性。这可能是因为多台多震和达法利用的震相数据不仅包含选定区域的台站震相，还包含更大范围内能够记录到这些地震的台站震相数据，计算结果是较选定区域范围更大的波速比平均值，反映的是研究区域波速比变化的总体趋势，忽略了较小范围的波速比变化细节；而单台多震和达法采用固定台站记录的震相数据，参与计算的地震的方位和范围可以根据需要进行选择和限定，计算结果是台站至地震的射线区域的波速比平均值，极大地提高了空间分辨率。但是该方法也存在不可忽视的缺陷，单台波速比计算过程受主观因素的影响很大。若选定特定范围或者特定方位的地震来计算波速比，通过一些条件的限定能够与地震很好地对应，改变限定条件则可能得到完全不同的结果，这为判断地震危险性带来很大困难。在日常工作中将这两种方法相结合，可达到更好的效果。

4 结语

本文通过唐山地区单台和多台波速比结果的对比研究认为，多台波速比结果较为稳定，而单台波速比波动范围大；异常台站的方位与地震具有一定的相关性。由于单台波速比的地震射线是台站附近地震震源到台站下方集中穿透，提高了波速比的空间分辨能力，所携带震源至台站的地壳介质信息更加丰富，可以为地震位置的判断提供一定的参考。