2018-02-12河北省永清县发生M4.3地震(简称永清地震)，是继2006-07-04文安M_S5.1地震之后，首都圈近10 a来发生的规模最大、震感最强烈的一次地震事件^[1]。地震发生在河北平原地震构造带中南段，震区发育的断裂主要有河西务断裂、桐柏断裂、牛驼镇凸起东缘断裂、固安断裂和大兴凸起东缘断裂等^[2]。地震发生后，中国地震台网中心(CENC)、美国地质调查局(USGS)、河北省数字地震台网中心定位的震源深度结果分别为20 km、10 km、17 km，结果差异较大。

震源深度是地震学研究的关键参数之一，精确的震源深度对于地壳流变学性质、孕震机理、发震构造、灾害评估等研究具有重要意义^[3-13]。然而，由于台网分布相对稀疏，且在深度上台站位置均处于震源一侧，传统的定位方法在测定震源深度上经常会出现较大误差。研究表明，基于波形测深的深度震相法可以有效提高测量精度。

为精确测定永清地震的震源深度，本文利用CENC提供的近震Sp沉积层基底界面转换波与USGS提供的震中距在30°~90°的远震深度震相sP、pP联合分析，确定地震的震源深度，并在计算结果的基础上，结合首都圈地区的地质结构，研究讨论区域地震的发震机理。

1 地震数据选取 1.1 近震波形数据

本文收集了河北省和天津市区域数字地震台网中心数字地震波形资料，根据Sp沉积层基底界面转换波的优势震中距为100 km以内，选取信噪比较高、沉积层基底界面转换波震相较为清晰的天津武清台(WUQ)的波形数据作为研究对象，台站距离震中约为26 km。

1.2 远震波形数据

为利用远震深度震相sP、pP测定震源深度，在USGS提供的数据中挑选震中距30°~90°范围内高信噪比的远震台站波形数据，将数据去仪器响应、分量旋转，再进行频率为0.7~2 Hz的带通滤波，以便在垂直(Z)分量上识别到较为清晰的远震深度震相sP、pP。通过逐台数据处理，挑选了BRVK、KURK、SIMI、KBL、TIXI、COLA、KMBO等7个台站数据，具体信息见图 1和表 1。

1.3 速度模型

针对近震波形数据，使用F-K波数与频率方法^[14]拟合波形，研究区一维速度结构模型参考董一兵^[15]的modWUQ模型及首都圈地区地壳速度模型^[16]构建而成(图 2)，震源区沉积盖层内P波波速、S波波速梯度增加，覆盖层下方P波、S波的平均波速分别约为6.1 km/s和3.5 km/s。针对远震波形数据，采用Teleseis远震波形拟合程序^[17]拟合波形。Teleseis速度结构模型(表 2)在震源处选用近震研究时构建的区域模型，中间点处和台站处则选用PREM标准地球分层模型。

2 原理和方法 2.1 基于近震沉积层基底转换波测定震源深度的原理和方法

Sp沉积层基底界面转换波震相是指从震源出发的Sv波入射到沉积覆盖层下方时，将一部分能量转换为P波，经在沉积层内传播后被台站接收到的波形，射线路径见图 3(a)，其优势震中距为100 km范围内。Sp基底转换波到时位于P波与S波之间，优势频率为2~5 Hz，能量主要集中在垂向和径向，其中垂向能量要强于径向，通常情况下切向能量很弱。研究表明，沉积层基底界面转换波Sp与P波到时差对震源深度敏感，随深度呈线性变化，可用来约束震源深度。

F-K波数与频率方法可以合成近震理论地震波，该方法通过对频率和波数分别进行积分，并采用传播矩阵计算地震的全波场位移分布，能够计算各种频率下的体波和面波波形，对于研究震源参数非常有效^[18]。在双力偶点源条件下，需要计算倾角为45°、滑动角为90°的倾滑断层DD，倾角为90°、滑动角为90°的倾滑断层DS，倾角为90°、滑动角为0°的走滑断层SS等3个基本断层的8个格林函数。根据地震的震源倾角、滑动角和台站方位角，震源时间函数及地震矩张量可以合成台站处R、T、Z三分量的理论地震波形，其实质是震源辐射花样系数和格林函数的卷积，具体公式为^[19-20]：

$u_{j}(t)=M_{0} S\left(t_{d}\right) \sum\limits_{i=1}^{3} A_{i j}((\varphi-\theta), \delta, \lambda) G_{i j}(h, \Delta, t) $

式中，下标i表示垂直走滑断层、垂直倾滑断层、45°倾滑断层3种地震断层模型，下标j表示自由界面上一点的R、T、Z三分量，M₀为标量地震矩，S(t_d)为震源的时间函数，A_ij为辐射花样系数，G_ij为格林函数。

2.2 基于远震深度震相sP、pP测定震源深度的原理和方法

远震深度震相pP、sP是指地震震源辐射出的上行P波(S波)在地表自由界面发生反射(偏振)而形成反射P波(转换波S波)，其传播路径与垂直入射的直达P波几乎相同，此类震相称之为远震深度震相pP、sP^[21]，射线路径见图 3(b)。远震深度震相pP、sP的优势震中距为30°~90°，P波与远震深度震相到时差受震中距影响不大，而对震源深度敏感。

Teleseis程序是合成远震理论地震波形的有力工具。在一维平行分层地壳速度结构模型中，地震震源释放的能量在震源处与台站处的脉冲响应可以根据Haskell传播矩阵方法计算，在地幔中(地震波传播发生转折的中间点处)的脉冲响应可以根据射线追踪理论中的几何扩散因子和衰减因子计算，台站处远震体波频率域的格林函数为^[22]：

$G_{i j}^{k}(h, \Delta, \omega)=R_{s}^{i j}(h, \Delta, \omega) *\\ R_{m}^{k}(h, \Delta, \omega) * R_{r}^{i j}(h, \Delta, \omega) $

式中，R_s^ij(h, Δ, ω)、R_m^k(h, Δ, ω)、R_r^ij(h, Δ, ω)分别为震源处、中间点处、台站处地震波的脉冲响应，h为地震震源深度，Δ为地震台站震中距，ω为角频率，i=1、2、3为断层类型，j=1、2、3为R、T、Z三分量。

中间点处地震波脉冲响应可表示为：

$R_{m}^{k}(h, \mathit{\Delta }, \omega)=G(h, \Delta) Q(h, \mathit{\Delta }, \omega)[\operatorname{Ref}(h, \mathit{\Delta })]\\ Q(h, \mathit{\Delta }, \omega)=\mathrm{e}^{-\frac{\omega}{2} t^{*}(h, \mathit{\Delta }, t)} $

式中，G(h, Δ)为几何扩散，Q(h, Δ, ω)为衰减因子，Ref(h, Δ)为反射系数因子。

根据以上公式可以计算出3个基本断层的格林函数，根据弹性波理论，将震源的辐射花样系数同3个基本断层的格林函数卷积，获得台站处R、T、Z三分量的位移。

3 计算结果 3.1 近震数据确定震源深度

对武清台观测数据进行去仪器响应、分量转换等预处理，选用带通滤波器进行滤波，滤波范围为2~5 Hz，通道数选择2以减小相移，极数选择1以压制旁瓣效应。采用F-K波数与频率方法计算武清台不同震源深度的格林函数，选用走向302°、倾角70°、滑动角-30°^[2]作为震源机制解，观测台站位于地表下246 m，震源持续时间为0.3 s，计算合成理论波形，并通过改变震源深度的方式调节理论波形中P波与Sp基底转换波的到时差。图 4为震源深度测试，观察不同深度的理论波形可以清晰地发现，震源深度对P波与Sp基底转换波的到时差敏感。随着震源深度的增加，近震沉积层基底转换波逐渐后移，震源深度为17 km，实际观测波形和理论波形Sp震相拟合度较好。利用武清台Sp沉积层基底界面转换波震相测定的永清地震的震源深度应该为17 km，由于近震沉积层基底转换波振幅比较小，为能有效识别该震相，需继续对武清台波形数据进行分析研究。

为验证Sp沉积层基底界面转换波震相，绘制武清台三分量理论地震波形(图 5)，其中Sp沉积层基底界面转换波震相的到时在P波之后约3.4 s，Z向和R向上有明显振幅，而T向上没有振幅，符合其震相的特征。以10°为步长，计算并绘制震中距26.1 km、方位角0°~360°的垂向理论地震波形，结果见图 6。由图可知，方位角在50°~100°范围内，Sp沉积层基底界面转换波震相能量较弱，0°~40°和140°~170°时能量较强，而武清台处于震中70°方位角方向，因此Sp沉积层基底界面转换波震相振幅较小，反映了永清地震的辐射花样导致在不同方位角上转换波能量不同。

3.2 远震数据确定震源深度

BRVK台震中距为34.40°，将三分量原始数据进行去仪器响应、变换为位移、分量转换和0.7~2 Hz带通滤波，去仪器响应时4个频率参数f₁、f₂、f₃、f₄分别选择为0.1 Hz、0.5 Hz、10 Hz、20 Hz。BRVK台R、T、Z三分量实际数据见图 7，P波在Z分量上振幅最强，R分量上振幅稍弱，而在T分量上几乎没有振幅，在Z分量上P波之后7 s可较清晰地看到深度震相pP。

利用Teleseis程序计算BRVK台不同震源深度的格林函数，选择总点数n_t为2 048、采样间隔dt为0.1 s，选用走向302°、倾角70°、滑动角-30°作为震源机制解，震源持续时间为0.3 s，0.7~2 Hz带通滤波，计算合成理论地震波形，并对Z分量进行震源深度测试，结果见图 8。通过理论波形可以看出，随着震源深度的增加，P波与深度震相pP到时差呈线性增加，在16 km处实际观测波形和理论波形pP震相的拟合度最好。从测定原理来看，利用BRVK台的远震深度震相pP测定的永清地震震源深度应该为16 km。

图 9为震源深度16 km时KURK、KBL、SIMI、COLA、KMBO、TIXI等台站实际观测波形与理论地震波形的对比，其中KURK台震中距近似30°。通过对比分析发现，KURK、KBL、SIMI、COLA台的理论及实际P波和深度震相pP的到时、振幅基本吻合，说明利用远震深度震相计算得到的永清地震的震源深度结果较为可靠；KMBO、TIXI台的信噪比相对较差，波形拟合效果不理想。COLA台位于美国阿拉斯加，震中距为58.42°，图 10为COLA台Z分量的震源深度测试，为得到更清晰的远震深度震相，将COLA台Z分量的原始数据去仪器响应及0.7~2 Hz带通滤波。由理论地震波形可知，随着震源深度逐渐增加，深度震相pP的到时也逐渐增大，但COLA台P波初动较弱，震相到时拾取不准，利用COLA台远震深度震相难以有效确定永清地震震源深度。

4 讨论

利用近震沉积层基底转换波Sp测定震源深度会受到速度模型不确定度的影响，将近震速度模型沉积层厚度分别增减1 km，产生2个扰动模型，测试震源深度结果对模型偏差的依赖性。通过合成理论地震波形，拟合近震沉积层基底转换波Sp与P波到时差，确定震源深度结果见图 11，沉积层厚度增减1 km后，震源深度分别为16 km和19 km，偏差分别为1 km和2 km。测试结果说明，利用近震深度震相进行深度定位受速度结构影响较大，尤其是沉积层的速度结构，因此精确的速度结构模型可以提高震源深度的测定精度。

将王晓山等^[2]的永清地震震源机制解的走向、倾角、滑动角分别增减10°，测试近震台站武清台和远震台站BRVK台的波形拟合情况，结果见图 12和13。结果表明，其对震相到时的影响不明显。

近震沉积层基底转换波Sp震相与其参考震相走时差的精度直接影响震源深度测定的准确性，有效识别和确认近震沉积层基底转换波Sp震相的方法一般有3种：1)通过分析地震波到时和偏振特征，近震沉积层基底转换波Sp震相在Z分量非常明显，T分量几乎没有，震相到时在直达P波和直达S波之间；2)通过TauP计算地震的震相走时，确定近震沉积层基底转换波Sp震相；3)当转换波振幅较小时，以震中为圆点、震中距为半径，构建许多不同方位角的虚拟台站，利用F-K波数与频率方法拟合各个台站的波形，由于地震震源机制解的不同，每个地震会有各自的地震辐射花样，不同方位角辐射的转换波能量会不同，将不同方位角的理论地震波形按P波到时对齐，进行震相识别。

5 结语

渤海湾盆地具有较厚的沉积覆盖层，使得永清地震的震源深度不易精确测定。本文基于近震沉积层基底界面转换波Sp震相和远震深度震相sP、pP对震源深度敏感的特征，运用F-K波数与频率方法拟合近震理论地震波形，采用Teleseis方法拟合远震理论地震波形，约束此次地震的震源深度。综合分析认为，永清地震的震源深度为16~17 km。该方法对于沉积层较厚地区地震震源深度的定位有较好的效果。

受区域张性断陷作用影响，首都圈及邻近地区活动断裂控制的地堑和半地堑断陷盆地是破坏性地震的主要活动场所，永清地震发生在华北克拉通东部破坏区，是一个发生在中上地壳的板内地震。结合已有地质资料表明，震源下方存在明显的上地幔隆起，上地幔软流物质上涌及该深度下温度和压力作用出现的物质部分熔融和蠕变，可能是引发地震活动的原因之一。同时，渤海湾盆地构造演化历史复杂，地震的孕育也可能有多种因素，未来还需对永清地震的成因及区域内的地震活动性进一步研究。