作为重要的震源参数之一，震源深度对于孕震深部构造、断层展布、震源过程、余震发展、核爆监测等有着重要意义，同时震源深度也是地震灾害评估的重要参数^[1]。在同等条件下，震源越浅产生的灾害就会越大，所以震源深度是地震成灾的关键指标^[2]。

北京时间2015年9月14日18时10分河北昌黎地区发生M_L4.5地震，秦皇岛、唐山地区震感强烈，北京、天津、河北承德及辽宁部分地区有震感，速报结果显示这次地震深度为14 km，而河北台网编目结果显示震源深度为5 km，二者差异较大，为后续发震构造、断层展布、灾害评估等相关研究造成了一定的困难，因此需要重新测定此次地震的震源深度。由于缺乏有效的地下约束条件，且缺少近台时，仅利用一种定位方法难以得到准确的震源深度，因此大部分学者会考虑综合多种定位方法以得到准确度较高的深度值^[2-4]。

目前，深度定位方法主要分为走时反演与波形反演两类。其中，基于走时反演的核心思想是：首先假定一个初始震源位置与发震时刻，基于速度模型和台站位置计算台站的理论到时，并与实际到时拟合比对；之后对初值不断修改迭代，当理论到时与实际到时的差异满足一定判据时，修改后的初值即为最终的定位结果。昌黎地震速报与编目过程中分别采用的盖革法、单纯型法等均属于此类^[5]。这类方法对于台站分布、速度模型要求较高，只有当震中距小于1.4倍的震源深度时，才能计算得到较高精度的震源深度。但是实际情况中，地震发生是随机的，台站附近20 km范围内地震发生几率并不高，所以震源深度小于15 km的浅源地震其深度定位误差一般较大^[6]。

相较于震相走时，地震波形则包含了更加丰富的信息，如震相到时、振幅、周期、衰减时长等。当波形质量高时，可利用波形反演方法得到精度更高的震源深度，比如全波形反演CAP方法以及矩张量反演TDMT方法。其主要思想是计算并搜索反演得到的理论波形与观测波形之间拟合误差，当误差最小时，可同时得到最佳震源机制解与最佳拟合深度^[7]。地震波形中的震相多种多样，部分震相因其对深度较为敏感，被研究者用于研究震源深度。比如近震深度震相有sPg、sPmP、sPn、sPL以及各自对应震相Pg、PmP、Pn、Pg等，这两类震相的到时差与震中距关系不大，而与震源深度相关^[8]，所以被很多人用于震源深度计算，且取得了很好的成果。因此，为得到昌黎M_L4.5地震较为准确的深度结果，将采用上述的CAP、TDMT以及深度震相sPL等3种方法计算震源深度，并对这3种方法在该地区的适用情况进行分析。

CAP方法原型是由和赫尔姆伯格于1994年提出的，后经朱露培、和赫尔姆伯格进一步改进发展为CAP方法，即“Cut and Paste”（剪切与粘贴）。该方法是将近震记录波形分成Pnl段和面波段，在反演的过程中分别对两段赋予不同的权重并拟合。反演过程中，以理论波形与观测波形f(t)一致性作为判断标准，当拟合结果残差最小时，对应的走向、倾角、滑动角作为震源机制解，此时拟合所得的残差最小的深度值作为最终震源深度。

该方法为了降低速度模型带来的不确定性，在拟合过程中将Pnl波和面波分别作时间平移，同时相对提高Pnl波的权重，对震源深度测定有较好的约束^[10]。

采用的矩张量反演方法（time domain moment tensor inversion，TDMT_INV）是由Dreger等^[12]提出的，利用区域范围长周期体波三分量波形在时间域反演地震矩张量解。反演基本原理如下：利用最小二乘法计算不同深度的理论波形，并与实测波形拟合，拟合程度越高，得到震源机制解越可靠，对应的即为最佳震源深度。

sPL震相是由崇加军等人提出的一种可以表征震源深度的震相。该震相是由SV波入射到地表转换而来（图1）。当震中距大于$\sqrt {10} $倍的震源深度时，t_sPL-P与震中距关系不大，而与震源深度线性相关^[14]。因此，可利用sPL震相确定震源深度。该震相的优势震中距范围为30～50 km，传播速度小于P波、大于S波，表现为低频，径向能量大于垂向，切向能量最小且不易观察。

图  1  sPL传播路径

2015年9月14日18时10分，河北昌黎地区发生M_L4.5地震，震中为118.79°E、39.73°N。当地及附近震感强烈，河北地震台网大部分台站均清晰记录到该地震事件，为本次深度研究提供了高质量的波形数据（图2）。

图  2  昌黎M_L4.5地震震中分布图

采用的3种测深方法均是基于理论波形与观测波形拟合得到最终结果，在计算理论波形的过程中，速度模型是基础资料，其一致性与准确度均会影响最终计算结果。为了保证研究结果的一致性，本文参考了前人研究成果，使用表1中的速度模型^[15]，该模型是参考孙若昧等^[16]首都圈地区地壳一维P波速度模型所得。

选择震中距250 km范围内的宽频带台站，对其波形去除仪器响应并积分至位移记录，再将位移记录的东西向（EW）、南北向（NS）、垂直向（UD）分别旋转至径向（R）、切向（T）、垂向（Z），手动标识P波到时。之后将波形分成体波与面波两部分，分别赋予不同的权重，并对体波、面波分别采取0.05～0.20 Hz、0.05～0.1 Hz滤波，以减小地壳精细结构和噪声带来的影响^[17]。基于表1中的速度模型，利用F-K方法计算得到不同台站的理论波形，并与观测波形相拟合，通过比较理论与观测波形拟合结果的相关系数，最终得到昌黎地震震源机制解（图3a）以及对应拟合深度（图3b）。

图  3  昌黎M_L4.5地震CAP震源机制解相关结果图

由图3可看出，理论与观测波形拟合相关系数大部分均大于50%，且大部分在80%以上，拟合度较好，另外台站的初动方向与震源机制解结果吻合较好，最终得到震源机制解结果为节面I走向286.1°，倾角82.0°，滑动角37.4°；节面II走向190°，倾角53°，滑动角170°。对比其他研究此次地震的相关文献可知，本文计算的震源机制解与郭蕾等^[18]、李慧等^[19]研究结果基本一致，说明此次震源机制解结果较为可靠。图3b为对应的拟合深度，结果显示为4.5 km，与河北台网编目结果接近。

与CAP方法类似，TDMT方法也是利用观测波形与理论波形的拟合程度判定最终的震源机制解。观测波形是指经过一系列处理之后得到的数据波形，处理过程包括去均值、去倾斜分量、反褶积仪器传递函数、积分并将三分向分别旋转到切向、径向与垂向，以及滤波等一系列操作。理论波形是利用F-K方法计算格林函数得到的，由于计算格林函数时，震源深度对理论波形影响较大，会影响最终的反演结果^[20]，所以在计算过程中会进行深度迭代，反演计算不同深度对应的理论波形^[21]，再与处理后的观测波形拟合，对比拟合残差与方差缩减等参数寻求最佳震源机制解以及对应的最佳深度。

在TDMT计算过程中，每一个震源机制解会对应几个参数，包括方差缩减值（VR）、纯双力偶源（DC）、补偿线性矢量偶极（CLVD）及残差（RES/Pdc）等等。其中，VR较大、RES/Pdc较小时，对应的结果为最佳结果；DC与CLVD二者互补，DC值越大，则说明地震符合纯双力偶机制。

针对此次地震给定的深度范围为4～13 km，昌黎地震矩张量解中VR、DC、CLVD、RES/Pdc随深度变化见图4。图中可清楚地看到，当震源深度为5～6 km时，VR最大为91.3%，RES/Pdc最小为3.85$ \times $10⁻¹³，DC均为100%。以深度5 km结果为例（图5），最佳震源机制解节面I的走向、倾角、滑动角分别为287°、75°、17°，与CAP计算结果很是接近。上述说明此次地震为纯双力偶机制，震源深度为5～6 km。

图  4  昌黎M_L4.5地震矩张量解中VR、DC、CLVD、RES/Pdc随深度变化

图  5  昌黎M_L4.5地震TDMT震源机制解（震源深度为5 km）

据之前研究发现，不同地区sPL震相发育的优势震中距范围不同^[13]。为寻找到较多的sPL震相，将震中距范围扩大为20～80 km，此范围内共有8个台站：CHL（震中距26 km）、DOH（震中距43 km）、TLK（震中距50 km）、JTG（震中距58 km）、BDH（震中距59 km）、TAH（震中距62 km）、QIX（震中距66 km）、QIL（震中距77 km）。对上述8个台站波形进行转置、积分、去均值以及1～4 Hz低频滤波处理，最终在CLI（昌黎台）、BDH（北戴河台）两个台站的P波之后，S波之前观测到一个低频震相，对比R、Z、T三分向发现，该震相在R方向振幅最大，Z方向次之，T方向振幅最小，符合sPL震相特征，所以认为此震相即为sPL震相（图6）。

图  6  两台站处理后波形

对比CAP与TDMT计算的震源机制解结果，二者相差不大，所以本文以CAP计算得到的节面I结果作为计算理论波形的震源机制解。利用F-K方法分别计算河北昌黎地震的CLI、BDH等2个台站深度范围1～16 km的理论波形，按照起始震相到时对齐，与2个台站处理后的观测波形比对。对比结果显示，当震源深度为5 km时，理论波形与观测波形拟合程度最高，故认为利用sPL震相计算此次地震的震源深度为5 km（图7）。

图  7  两台站sPL震相观测波形与理论波形拟合图

利用CAP、TDMT、sPL三种方法测得河北昌黎M_L4.5地震震源深度分别为4.5 km、5～6 km、5 km，与速报结果有所差异，与河北台网编目结果基本一致。根据河北台网目录可知，昌黎M_L4.5地震发生13 min后，同一震中位置再次发生M_L3.9地震，速报深度为13 km，编目结果为5 km，与本次地震情况类似。由地震波形可知，两次地震记录台站一致，最近3个台站分别为LUX、CHL、DOH台。其中，LUX台所记录到两次地震的Sg与Pg到时差分别为1.42 s、1.43 s，而CHL、DOH两台所记录到两次地震的Sg与Pg到时差相同，分别为3.33 s、5.50 s。以上数据均表明，两次地震发震位置基本一致，震源深度应相差不大。因此，针对M_L3.9地震也采用上述3种方法、同一速度模型计算了震源深度。结果显示，利用CAP方法计算得到M_L3.9地震震源机制解为节面I走向282°、倾角52°、滑动角61°，节面II走向144°、倾角46.4°、滑动角121.8°，拟合深度为4.2 km（图8a）；在利用TDMT计算过程中由于台站数据信噪比低，未能成功计算出此次地震的震源机制解及震源深度；利用sPL震相计算过程中，仅在CHL单台识别到sPL震相（图8b），利用F-K方法计算CHL台理论波形，与观测波形拟合（图9），显示当震源深度为5 km时，理论与观测拟合效果最好。故M_L3.9地震震源深度应在4～5 km范围之间，与M_L4.5地震计算结果基本一致，符合这两次地震的波形特征。

图  8  河北昌黎M_L3.9地震

图  9  昌黎M_L3.9地震CHL台观测与理论sPL震相拟合图

2008年1月至2019年9月期间，河北昌黎地区发生M_L1.0以上地震251次。统计发现，所有地震震源深度基本呈正态分布，主要集中在4～11 km范围内，与本文计算结果基本一致。

综上认为，河北昌黎M_L4.5地震震源深度范围为4～6 km，深度较浅，位于上地壳。

1）本文采用CAP、TDMT、sPL等3种方法测定昌黎M_L4.5地震震源深度，其中CAP方法中，当震源深度为4.5 km时，理论与观测波形拟合最好，TDMT方法测定的深度为5～6 km，sPL震相拟合的深度为5 km。3种结果中，CAP方法和TDMT方法计算得到震源深度为地震的矩心深度，通过计算地震矩、应力降等震源参数得到昌黎M_L4.5地震的震源尺度约600 m，所以结合两种方法的计算结果认为，CAP、TDMT两种方法计算的深度结果为4～6 km。sPL震相来源于震源，表征的也是矩心深度，而且据李志伟等人^[22]研究结果认为，sPL震相主要受限于速度模型，模型参数变化10%，震源深度为10 km的地震测定结果偏差1 km，所以假定本文所用的速度模型与真实结构相差10%，深度结果也相应增减1 km为4～6 km，3种结果一致性较好。为验证所得结果，采用同样方法对发震时间与该地震相差13 min的昌黎M_L3.9地震进行了测定，结果一致，故综合认为河北昌黎M_L4.5地震的震源深度范围为4～6 km。

2）本文计算结果表明，昌黎M_L4.5地震为浅源地震，根据昌黎地震震中位置可知，本次地震位于唐山震源区内卢龙断裂与滦县-乐亭断裂交汇处^[18]。刘启元等人^[23]研究认为唐山大震区中上地壳具有明显的非均匀壳内低速体，会导致震源深度变浅；王晓山^[24]也曾指出，唐山震区自西南至北东方向震源深度表现为深一变浅一深一更浅的分布，且滦县-乐亭断裂之后的卢龙断裂处震源深度越来越浅，本文结果与以上研究结论相符。

3）河北台网速报与编目结果的不一致性，应该是因为二者使用的台站数量、速度模型、定位方法不同导致结果的差异性，所以在实际工作或者研究中可采用不同的定位方法，尤其是对争议较大的地震可采用多种方法综合判定震源参数。

4）本文所采用的CAP、TDMT、sPL3种方法中CAP与TDMT方法对台站数量、布局有较高要求，台站数量越多且包围性越好，结果越可靠。但是这两种方法计算过程较为繁琐，需要不断地优选台站，不确定因素较多。从文中可看出，当地震震级较小时，波形数据信噪比降低，TDMT无法完成计算；而sPL震相仅靠单台识别即可确定震源深度，相对来说应用更为方便，计算过程也较简单，尤其是有多台可识别sPL震相时，可互相印证以提高结果的可靠性与准确性。不过该方法也具有局限性，要求使用宽频带台站记录波形，且只有在结构相对简单、一维地壳速度模型能够较好描述的研究区域，sPL震相才较为发育^[9]。因此，在实际应用中可针对不同情况，选择多种方法计算，以得到可靠性较高的深度结果。