0 引言

唐山地区位于阴山-燕山隆起带东段的燕山隆起和华北平原拗陷带北部的冀渤凹陷的结合部位，在燕山隆起一侧，区内浅层活动性断裂发育，地质构造复杂，地震活动频繁(徐杰等，1996；罗艳等，2008)。1976年唐山7.8级大地震发生在唐山断裂附近，震后由于老震区地壳介质相对破碎，不易于积累更大能量，唐山地区地震序列强度和频次总体随时间呈现阶段性衰减的趋势，近10年来仅发生一次M_L≥5地震，即2012年5月28日唐山M_L5.2地震。此次地震是在老震区5级地震平静17年背景下发生的，且与地震序列强度衰减趋势相比震级明显偏高。研究表明地震孕育与发生是应力积累与释放的过程，在这个过程中通常会引起地壳介质的物性改变，这种改变可通过地震波速的变化体现(Birch，1961)。因此，获得该地震前后不同时间段的唐山老震区的三维速度结构，总结速度结构的变化特征，对分析和探讨该地震的形成、产生以及诱发因素具有意义。

诸多学者对唐山地区的速度结构进行了研究，并取得大量成果。陈立华等(1990)根据构造特征将华北地区划分为4个区域，研究该区P波速度结构，结果表明横向不均匀性和低速的存在可能是构造活动区的最主要特征；朱露培等(1990)利用空间频率域离散参数化反演方法，计算京津唐张地区的地壳和上地幔深部结构，发现研究区内大多数强震发生在低速区向高速区过渡的梯度带；孙若昧等(1999)采用分类表格数据的AIC分析方法，对华北地区强震发生位置与速度结构之间的关系进行定量化评估；卢造勋等(2002)采用正交投影法重建了研究区(110°~126°E，36°~44°N)的地壳和上地幔的三维速度结构，结果显示了“易震层”的速度分布特征，强震多发生在低速层上方的高速层内；王志铄等(2008)利用LSQR方法反演了华北及邻区地壳上地幔三维速度结构，发现唐山地区在25km、42km和60km深处存在连续的低速异常；赖晓玲等(2013)利用唐山临时地震台网数据5000多个P波到时进行速度成像，结果表明唐山断裂带为速度陡变带，其两侧的速度存在显著差异，西北侧为高速区，东南侧为低速区；杨峰(2019)采用双差层析成像的方法反演得到华北北部地区地壳三维P波速度结构，发现唐山震源区下方存在明显的低速异常区。

虽然前人在唐山地区速度结构方面开展了大量研究，但是上述研究基本是静态成像，未分析地震前后波速的动态演化，然而正是这种动态演化过程，可能有助于深入探索地震的发震机理(周龙泉等，2007；雷建设等，2010)。本研究采用结合波形互相关的双差层析成像方法，以2年为时间间隔对唐山地区2010年5月~2016年5月的地震数据进行反演，获得2012年唐山M_L5.2地震前后介质的三维地震P波、S波速度的空间分布，结合此次地震的精定位震源参数，探讨唐山地区中强震发生前该区介质的波速特征。

1 理论方法与数据 1.1 双差层析成像方法

双差层析成像方法(Zhang et al，2003、2007)结合绝对走时和差分走时资料对三维速度结构和震源参数的联合反演。假设有2个相邻地震的距离远小于其到同一台站的距离，则可忽略速度差异，认为它们到同1个台站的路径相似。对地震对的走时残差做差即可得到双差。

根据射线理论，地震事件i到地震台站k的体波到时可以表示为积分形式

$ T_{k}^{i}=\tau^{i}+\int_{i}^{k} u \mathrm{d} s $

(1)

式中，τⁱ为地震事件i的发震时刻；u代表慢度；其中震源坐标(x₁，x₂，x₃)、发震时刻、射线路径及慢度未知。将式(1)在震源处展开为一阶Taylor级数得到线性方程

$ r_k^i = \sum\limits_{l = 1}^3 {\frac{{\partial T_k^i}}{{\partial x_l^i}}} \Delta x_l^i + \Delta {\tau ^i} + \sum\limits_{m = 1}^{{M_{ik}}} {\sum\limits_{n = 1}^N {{w_{mn}}} } \Delta {u_n}\Delta {s_m} $

(2)

式中，r_kⁱ为地震事件i到台站k观测到时与理论到时的残差；M_ik为地震事件i到台站k的射线路径的分段数；w_mn为以第m个节点为中心的第n个网格节点插值权系数。同理地震事件j到台站k的走时方程也可表示为式(2)的相同形式，对式(1)、(2)做差即可得到双差方程

$ \begin{array}{l} r_k^i - r_k^j = \sum\limits_{l = 1}^3 {\frac{{\partial T_k^i}}{{\partial x_l^i}}} \Delta x_l^i + \Delta {\tau ^i} + \sum\limits_{m = 1}^{{M_{ik}}} {\sum\limits_{n = 1}^N {{w_{mn}}} } \Delta {u_n}\Delta {s_m}\\ - \sum\limits_{l = 1}^3 {\frac{{\partial T_k^j}}{{\partial x_l^j}}} \Delta x_l^j - \Delta {\tau ^j} - \sum\limits_{m = 1}^{{M_{jk}}} {\sum\limits_{n = 1}^N {{w_{mn}}} } \Delta {u_n}\Delta {s_m} \end{array} $

(3)

式(3)也可表示为

$ d r_{k}^{i j}=r_{k}^{i}-r_{k}^{j}=\left(t^{\mathrm{obs}}-t^{\mathrm{cal}}\right)_{k}^{i}-\left(t^{\mathrm{obs}}-t^{\mathrm{cal}}\right)_{k}^{j}=\left(t_{k}^{i}-t_{k}^{j}\right)^{\mathrm{obs}}-\left(t_{k}^{i}-t_{k}^{j}\right)^{\mathrm{cal}} $

(4)

其中，$\left(t_{k}^{i}-t_{k}^{j}\right)^{\mathrm{obs}} $表示观测到时差，可通过地震观测报告或波形互相关计算求得；$ \left(t_{k}^{i}-t_{k}^{j}\right)^{\mathrm{obs}}$表示理论到时差，可通过正演求得。将所有地震事件对联合，便可以得到线性方程组，应用双差算子Q_DD便可求解方程组(Wolfe，2002)。

双差层析成像利用绝对走时数据确定震源区外的速度结构，利用到时差数据确定震源区的精细速度结构。其通过速度结构的适当调整提高地震定位的精度，而地震定位精度的提高又可以降低由于射线路径差异造成的速度结构的误差，进而提高速度结构反演的精度(吕子强等，2016)。因此，双差层析成像方法可以得到比标准层析成像更精确的定位和速度结构结果。

1.2 地震数据与模型

唐山老震区小震序列丰富，且台站分布较为合理，能有效地监测该区的地震活动，为双差层析成像提供了良好的数据基础。本文研究区域为39°00′~40°30′N、117°30′~119°30′E，收集河北地震台网2010年5月~2016年5月唐山地区的地震事件波形及震相报告，根据观测数据，2010年5月~2012年5月共1398个地震事件，P波走时16822条，S波走时17616条；2012年5月~2014年5月共1628个地震事件，P波走时18958条，S波走时20546条；2014年5月~2016年5月共3219个地震事件，P波走时22937条，S波走时23697条。2014年5月~2016年5月地震事件的数量与前2个时段有明显差别，这可能是因为近年来地震目录精度提高，存在一定数量的单台地震，为保证对比效果，需对数据进行筛选。

结合走时曲线(图 1)剔除误差较大的震相数据，且挑选出满足每个地震需要的最少记录震相数(8个)，去掉单台地震，筛选后，2010年5月~2012年5月可用的地震事件1271个，P波走时16280条，S波走时16598条；2012年5月~2014年5月可用的地震事件1423个，P波走时18274条，S波走时19396条；2014年5月~2016年5月可用的地震事件1520个，P波走时19246条，S波走时20001条。3个时段台站分布、地震数量及到时资料数基本相当。

在双差层析成像中，事件对之间的距离对反演结果至关重要。该阈值越小则反演结果精度越高，但取值越小则地震对的数量越少，不利于反演运算，因此确定合适的事件对距离是双差层析成像计算过程中重要的一环。通过统计研究区所有台站记录到的事件对的震源距和相应的波形互相关系数(图 2)，确定事件对距离为30km。经过筛选后的地震震中主要集中在唐山断裂附近，震源深度多数小于30km，为确保每个网格有足够的射线通过，将网格间距设置为0.2°(图 3)。根据于湘伟等(2010)的一维速度模型研究成果(表 1)，垂直向下的深度分别取0、5km、10km、15km、20km、25km和30km，与速度模型的分层一致。

2 结果与讨论

结合初始一维速度模型，采用带阻尼的最小二乘正交分解法，通过8组23次迭代获得唐山地区的三维速度结构，采用DWS(Derivative Weighted Sum)对反演结果的可靠性进行检验。

2.1 结果可靠性分析

DWS能很好地衡量节点的射线密度(Thurber et al，1999)，前人的研究认为在DWS大于100的区域具有很高的模型分辨率(Scarfì et al，2007)。本文采用DWS检测3个时段反演速度模型的分辨率。根据反演结果得到了P波DWS分布图(图 4~6)，结果表明，反演网格大部分区域的DWS大于100，具有较好的模型分辨率，3个时段的震源深度均在5~20km，尤其在唐山断裂附近，具有很高的射线密度分布，能保证得到较好的分辨率。

2.2 2012年唐山M_L5.2地震前后三维速度结构对比分析

精定位后2012年唐山M_L5.2地震的基本参数如表 2所示。精定位后该地震的震源深度为11km，且震中位于唐山断裂附近。DWS结果表明在5~20km的震源深度，尤其在唐山断裂附近具有较高的模型分辨率，且前后3个时间段的地震分布基本一致，获得的射线密度分布差异不大，反演过程从数据处理到反演阈值的选取均保持一致，因此认为2012年唐山M_L5.2地震前后三维速度结构的对比结果具有一定的参考意义。选取前后3个时期10km的水平速度结构切片进行观察比较。

由图 7可以看出，在10km深度，研究区P波平均速度结构横向不均性显著，存在2个NNW向的高低速交界区，分别与滦县-乐亭断裂和蓟运河断裂位置对应。沿NE向的唐山断裂带具有串珠状高低速的线性特征，反映了唐山断裂带是研究区内的主要构造带，这与赖晓玲等(2013)的研究成果基本一致。2012唐山M_L5.2地震前，震中位于在P波高低速过渡带偏向高速的位置(图 7)，出现这种现象的一种解释是，低速体难以积累能量，而高速体有可能积累足以引发大地震的应变能(孙若昧等，1995)。震后震中附近的波速有一定程度的降低，西北方向高速块体范围减小；而震中东南方向的低速区变为高速区，形成了1个高速块体(图 8)；之后，震中附近高低速过渡带逐渐消失。

2012唐山M_L5.2地震发生在唐山断裂附近，根据震源机制计算结果发现，有1个节面与唐山断裂走向一致，因此沿唐山断裂方向(AA′)和垂直唐山断裂方向(BB′)做剖面(图 9)，观察该地震前后在深度上的速度结构变化。由图 10、11可以看出，5~15km深度为P波低速区与高速区对比强烈的交界地带，唐山地区中下地壳存在低速异常，这与王椿镛等(2017)的研究成果一致。

王椿镛等(2017)的地震层析成像结果表明，华北大震大多发生在高速和低速区的过渡带上。水平速度结构切片与垂直剖面速度结构切片均显示，此次地震发生前，震中位置处于高低速过渡带，震后震中附近的波速出现一定程度的降低，与王椿镛等(2017)的研究成果保持一致。2012年唐山M_L5.2地震前，震中位置处在高低速过渡带附近，震后震中附近的波速出现的小范围的降低，震中下方的低速区范围逐渐增大，震中附近高低速过渡带逐渐消失。

波速比能够综合反映P波和S波变化的相对大小，单纯地反映介质的泊松比，所受影响比P波、S波速度小很多。图 12给出了10km深度的研究区波速比分布，唐山M_L5.2地震前震源区附近存在明显的波速比低值，震后波速比值回升。王亚茹等(2016)研究结果表明，2012年唐山M_L5.2地震前，震源区附近波速比出现了下降—低值的过程。由于本研究在震前仅计算1个时段的数据，无法判断震前是否存在下降—低值的过程，但是低值—恢复的过程与已有研究结果(王亚茹等，2016)一致。

3 结论

采用波形互相关的双差层析成像方法对唐山地区2010年5月~2016年5月的地震数据以2年为时间间隔进行反演，获得2012年唐山M_L5.2地震前后介质的三维地震P波、S波速度的空间分布。结果表明：

(1) 在10km深度，研究区P波平均速度结构横向不均性明显，沿NE向的唐山断裂带有串珠状高低速的线性特征，反映了唐山断裂带为研究区内的主要构造带。

(2) P波和S波水平速度结构切片与垂直剖面速度结构切片均表明，2012唐山M_L5.2地震前，震中位置处于高低速过渡带，震后震中附近的波速出现一定程度的降低，原有的高低速过渡带逐渐消失。

(3) 唐山M_L5.2地震前，震源区波速比存在低值异常，震后波速比值回升。

地震的形成是一个复杂的过程，用层析成像的方法探讨2012唐山M_L5.2地震前后波速演化过程时，虽然尽量保持不同时段地震数量、台站分布及反演参数的一致性，但是地震射线覆盖的差别，仍会对成像结果造成一定影响。因而，本研究中的速度变化是否仅仅受此次地震的影响而产生，仍需要进一步研究。