奥地利学者康拉德(Conrad, 1925)首次发现了康拉德界面后，该界面在大陆上其他许多地方相继被发现且其相关性质也得到了研究(如Berry and Fuchs, 1973; Kerr, 1989；Pavlenkova, 1993).康拉德界面是否普遍存在一直是一个具有争议性的问题, 一些研究(如Litak and Brown, 1989; Kerr, 1989; Pavlenkova, 1993)认为该界面与研究区的地壳结构及其形成和演化的历史密切相关, 在海洋地壳内和有些构造活动比较强烈的大陆地壳内康拉德界面可能是缺失的.然而近期已有研究者在地质构造活动比较强烈的地区发现来自康拉德界面的震相(如He and Hong, 2010；郭杰等，2013；焦煜媛等，2017).

Pb震相是近震波形记录中能够记录到的重要近震震相之一，其成因也是一直存在争议.Bormann(2012)认为Pb震相是康拉德(Conrad)不连续界面上产生的能识别的折射波或反射波.张少泉(1984)认为Pb可能是地壳内某种震相，其生成机制不是单一的，它并不一定相应于一个稳定的速度间断面或一个稳定的速度高梯度带.IASPEI(International Association of Seismology and Physics of the Earth′s Interior，国际地震学与地球内部物理学会)(Storchak et al., 2003)和地震编目规范(中国地震局，2017)将Pb震相定义为来自下地壳内震源的上行P波，或其底部到达下地壳的P波.由于识别困难，国内外关于该震相的专门研究结果不是很多，也较少有相关的定量分析研究.国内编目工作中有极少数区域台网分析了较少量的Pb震相.在有康拉德界面的区域，台网分析人员常常由于对Pb震相时频特征的认识不足，存在将该震相误识别为Pg或者PmP等其他震相的情况，给后续的地球物理研究工作(如地震定位，利用走时的层析成像)带来误差.对于Pb震相的典型时频特征进行综合分析与总结具有重要的实践和理论意义，将为台网编目人员提供分析该震相的技术标准，丰富台网观测报告震相产出种类并提高震相准确性，还为分析和研究康拉德界面的空间特征等提供了数据信息.

傅淑芳和刘宝诚(1991)指出，在震中距大于200 km时，震相出现的顺序为Pn, Pb, Pg.郭永霞等(2010)对首都圈台网记录的M_L≥2.0地震事件重新分析, 得到313个Pb震相，绘制时距曲线分析显示该震相到时读数离散度较大，很难确认Pb震相识别的准确性；他们还经理论走时计算及人工震相分析结果验证，推测Pb震相在华北地区几乎不作为初至震相到达，但是没有给出实测例子.张蕙(2012)在分析太原盆地波形特征时指出，太原盆地地区在震中距187 km以后，出现较清晰的Pb震相，震相到达的顺序为Pn，Pb，Pg；220 km左右，各震相的界限清晰，很好识别.焦煜媛等(2017)利用发生在甘东南地区岷县、漳县交界处的一个近震事件的波形数据进行了细致研究和分析, 并且借助波形数值模拟和走时曲线拟合准确确定了观测资料中来自康拉德界面的相关震相.

在Pn震相作为初至波出现后，有时Pb和Pg会紧随其后按顺序出现.这时读取的Pb具有相对可靠的信度.但是在Pn作为初至波出现的距离之前，Pb和Pg都有可能作为初至波出现.Pb震相携带了康拉德界面的信息，通过对Pb震相的分析研究，我们不仅为震相识别分析提供参考准则，还可以得出康拉德界面的详细信息.首都圈地区高密度测震台网长期连续的观测，积累了丰富的地震事件，然而，较少有相关首都圈Pb震相的详细分析与研究的报道.首都圈不同台站记录的Pb震相时间域的特征是什么？频率域的特征又是什么？与其他P波类震相相比，不同震中距内的到时顺序是什么？这些都是具有重要意义的科学问题.本研究将利用首都圈测震台网记录到的波形资料，识别Pb震相，分析其时频特征，为首都圈的Pb震相分析及康拉德界面的分布提供技术方案和科学支持.

我们搜集了首都圈及邻区M_L≥2.5的非天然地震10个(包括爆破和塌陷), 天然地震359个(为了使地震事件在空间上尽量有所分布，还额外增加了包括山西、内蒙古、河南等地区的部分地震事件，图 1)，共计369个地震事件, 用于Pb震相的识别和分析工作.

图 1

Fig. 1

图 1 用于Pb震相分析的首都圈及邻区ML≥2.5的369个地震事件的空间分布 Fig. 1 Spatial distribution of 369 seismic events with ML≥2.5 in the Capital Circle region around Beijing and adjacent areas used in this study

Pb震相识别是比较困难的.图 2为震源在4层水平均匀地壳速度模型中震源深度为10 km时，Pg，Pb，PmP，Pn的传播路径示意图(国内Pg定义包括上行的Pg和下行的Pg，这里国内的定义和国外的定义有些差别).我们的模型是一维水平分层均匀速度模型, Pb和Pn震相的射线包含了沿界面和稍微折射到界面之下由于垂向速度梯度作用而上折再次穿过康拉德界面达到地表的射线.当界面的速度对比不强，且速度模型为各向非均匀的3-D模型时，Pb和Pn震相可能存在从震源出发下潜穿过康拉德界面的深度更深，然后因垂向速度梯度作用而上折再次穿过康拉德界面达到地表(Sun et al., 2016；孙伟家等，2018).为了保证震相识别的准确性，我们在基于现有区域台网模型的理论走时作为参考的基础上，选择只有P波段至少出现3种震相时的数据，才进行Pb震相的挑选和识别工作.在河北省地震局与北京市地震局编目结果的基础上，根据震中距范围，重新校正Pg，Pn和PmP的到时，去掉到时相对不清晰的震相，并进行重新定位.

图 2

Fig. 2

图 2 基于本研究联合反演所得4层地壳速度模型，震源深度为10 km时Pg、Pb、PmP、Pn的传播路径示意图 Fig. 2 Schematic diagram showing propagation paths of Pg, Pb, PmP, and Pn phases computed using a 10-km depth source and the 4-layer crustal velocity model from joint inversion in this study

在确定Pb等相应震相的到时顺序时，我们采用多项式拟合方法(如边刚等，2015)，拟合相应的震相走时曲线，并结合震源深度和理论到时来判定在不同震中距范围内Pb所在的顺序.在分析Pb等相应震相频域的特征时，我们主要采用基于短时傅里叶变换的方法(如李稳等，2016)，进行有关频谱和能量的分析；在分析主频时，我们采用归一化功率谱密度方法pmtm(Percival and Walden, 1993)进行分析.

在分析康拉德界面分布特征时，我们利用Pb的射线路径勾勒康拉德界面空间分布.在去除因钟差引起的Pb震相到时异常的数据之后，采用基于遗传算法的联合反演方法(万永革和李鸿吉，1996)和最小二乘法分别确定康拉德界面的深度和相对速度.从而用实测地震事件数据确定首都圈及邻区康拉德界面的空间分布与相对深度等信息.

在利用已知速度模型计算震相理论走时和绘制传播路径时，我们采用Pyrocko方法(Heimann et al., 2017).

我们分析了186个台站记录的369个地震事件的数据.为了保证Pb震相识别的准确性，即使P波有3个震相的情况下，有争议的数据我们也不选择.如图 3a, 为2014年9月6日19时河北涿鹿M_L2.6地震北京市地震局LBP台的记录，由于第二个震相相对于Pg间隔一个半周期，即使有频率的明显改变，但是也有可能是观测中由于仪器滤波器的原因导致的频率不一致.图 3b是2015年11月3日00时河北唐山M_L2.8地震北京市地震局LLM台的记录，也是由于第二个震相距离Pn太近，不能确定是Pn正常衰减还是Pb的到来.图 3c是2015年12月18日20时河北怀安M_L2.5地震河北省地震局XIT台的记录，虽然第二个震相的振幅和频率都有稍微改变，但是不能确定是Pn衰减还是Pb的到来.如以上这些典型的具有争议的事件, 我们没有进行Pb的读取工作.在严格执行以上标准的基础上，我们在157个台站所观测的数据中识别到了Pb的记录情况，其余29个台站没有识别到Pb的观测记录.

图 3

Fig. 3

图 3 Pb识别过程中有争议不选的典型例子 Fig. 3 Typical examples of disregarded waveforms showing problems in identifying the Pb phase

图 4是2013年10月27日03时河北唐山M_L3.9天然地震不同台站记录到的Pb震相典型记录，从图 4a—f，震中距依次增加.图 5是2015年9月14日18时河北昌黎M_L3.9地震Pb震相典型记录，从图 5a—n震中距依次增加.图 6是2010年6月11日13时河北涞源M_L2.6爆破观测到的Pb震相典型例子.通过对本研究所识别出292个拥有Pb记录的地震事件(如图 4，图 5，图 6典型的天然和非天然地震记录)的震相进行分析，我们大概可以得出Pb震相相对Pg或Pn的高频部分频率通常相对高一些(如图 4—6)，振幅增大(如图 4a, 图 5c等)或减小、变化不明显(如图 4d，图 5e，图 5f等)的情况都有.对于Pb震相并未表现出十分明显的强振幅(如图 4c)，而是有较为明显的频率特征变化情况，根据孙伟家等(2018)对莫霍面以及岩石圈内部界面(MLD)的研究结果，认为康拉德界面在这些地区的速度对比并不强，而可能是一个物质和速度渐变的界面.对于可以看到明显的Pb震相(如图 4f)，可能意味着这条射线路径经过的康拉德界面速度对比较强.但相对来说，振幅相对增大时我们更容易识别，而振幅相对减小时，就为我们识别带来了更多的困难，有时候很难分清是否是Pb的到来(如图 3c).从形态上看通常Pb震相相对于Pg和Pn “毛刺”稍微多一些，但是也有特殊的情况(如图 4f，图 5b, c)Pb震相的确已经到来，但是频率却比Pn小，周期相对不变或较长，这可能与台站相对于震中方位、震源机制、场地响应等情况有关.

图 4

Fig. 4

图 4 2013年10月27日03时河北唐山ML3.9天然地震观测到的Pb震相典型记录 Fig. 4 Examples of selected Pb seismic phases for an ML3.9 earthquake at 03h 27 October 2013 in Tangshan, Hebei Province

图 5

Fig. 5

图 5 2015年9月14日18时河北昌黎ML3.9地震观测到的Pb震相典型记录 Fig. 5 Same as Fig. 4 but for an ML 3.9 earthquake at 18h 14 September 2015 in Changli, Hebei Province

图 6

Fig. 6

图 6 2010年6月11日河北涞源ML2.6爆破观测到的Pb震相典型记录 Fig. 6 Same as Fig. 4 but for an ML2.6 explosion event on 11 June 2010 in Laiyuan, Hebei Province

为了更好地了解Pb在频域内的特征，我们利用短时傅里叶变换的方法(如李稳等，2016)对同一地震不同台站观测的Pb记录，以及同一台站观测到的不同事件的Pb记录进行了时频分析.图 7与图 8分别是2012年5月28日10时22分河北唐山M_L5.1和2015年11月3日00时30分河北唐山M_L2.8地震不同台站记录到的Pb时频分析结果.从图可以看出，在频域内，Pb震相相对于Pn和Pg震相频率展布范围稍微高一些，能量相对强一些，但是不算十分显著，也有一些台站没有明显区别(图 7—8).图 9和图 10分别为位于平原和山地交界处的河北局XIT台与山西山区的SHZ台记录到的不同地震事件的时频分析结果.从图 9和图 10可以看出在Pb震相出现的时段内，其相对于Pg和Pn高频部分的能量有轻微的加强，高频成分相对较多.通过对其他事件和地震台站记录的Pb情况的分析，Pb存在相似的特征.

图 7

Fig. 7

图 7 2012年5月28日10时河北唐山ML5.1地震不同台站记录到的Pb时频分析结果(红色竖线为Pg到时，黄绿色竖线为Pb到时，黄色竖线为PmP到时，紫色竖线为Pn到时) Fig. 7 Time-frequency analysis of the Pb phase recorded at different stations for the ML5.1 earthquake at 10h 28 May 2012 in Tangshan, Hebei Province. Vertical lines are arrival times of Pg (red), Pb (yellow-green), PmP (yellow), and Pn (purple), respectively

图 8

Fig. 8

图 8 2015年11月3日00时河北唐山ML2.8地震不同台站记录到的Pb时频分析结果(竖线代表到时同图 7) Fig. 8 Same as Fig. 7 but for an ML 2.8 earthquake at 00h 3 November 2015 in Tangshan, Hebei Province

图 9

Fig. 9

图 9 河北省地震局XIT台记录到Pb震相的时频分析结果(竖线代表到时同图 7) Fig. 9 Same as Fig. 7 but for Pb phases from different events recorded at XIT station of Hebei Earthquake Agency

图 10

Fig. 10

图 10 山西省地震局SHZ台记录到Pb震相的时频分析结果(竖线代表到时同图 7) Fig. 10 Same as Fig. 7 but for Pb phases from different events recorded at SHZ station of Shanxi Earthquake Agency

为了更好地得到各种相应震相的定量主频特征，我们首先对台站随机噪声进行了分析.采用965个垂直向的5 s背景噪声，在去均值、去倾斜处理后，为了减小端点效应，对数据的两边加了COSINE 0.05宽度的Taper窗，然后对所有数据进行归一化处理，利用pmtm方法进行噪声主频分析, 并取平均值(图 11).考虑到首都圈部分地震计观测最小频率为1 Hz(周期1 s), 大部分仪器观测为0.0167 Hz (周期60 s), 还有少部分是0.0083 Hz(周期120 s), 我们采用1 Hz作为参考.由图 11可以看出，噪声的主频主要分布在0~1.5 Hz左右，后续高频部分均有分布，功率谱密度1.6~6 Hz左右平缓降低，6~48 Hz降低极其微小.

图 11

Fig. 11

图 11 首都圈及邻区台站随机噪声平均归一化功率密度图(黑色虚线为1.0 Hz参考线) Fig. 11 Normalized and averaged power density spectrum of station noises in the Capital Circle region around Beijing and adjacent areas. The vertical black dashed line is 1.0 Hz for reference

由于首都圈及邻区平原地区和西部山区的莫霍面深度不一致，在东部平原区深度为33 km左右，而在西部太行山区以及北部燕山地区的莫霍面深度约为40 km(如李松林等, 2001).因此，在平原地区Pn大概在120 km之后开始出现，而在太行山及燕山地区莫霍面相对较深的地方，Pn大概在140 km之后开始出现.在校核震相时，我们考虑了发生在不同区域的地震Pn出现的相对距离，根据地震发生的区域，将数据分为平原区和山区两类.

为了得到在相应震中距范围内的Pg，Pb，PmP，Pn的主频定量特征，在考虑震中距范围的情况下，我们采取了将相应清楚的震相波形各自单独切割出来分析(Pg震相作为后续波时截取长度为5 s)，并对切割出来的震相波形进行去均值，去倾斜处理.为了减小端点效应，对数据的两边加了COSINE 0.05宽度的Taper窗，并将所有数据按该数据集合中最长的一条波形为参考长度，将其余波形以波形中点为参考点两边补零到同等长度，对所有数据进行归一化处理，然后求功率谱密度，并画到同一张图中进行对比分析.由于波形中低频部分的很大一部分频率与随机噪声重叠，为了区别出随机噪声和端点效应的影响，我们又采用人工测量最大周期的办法来确定本研究中各个震相的肉眼能分辨的最小频率.选取了6个地震事件(2012年5月28日10时22分52秒发生的河北唐山M_L5.1地震，2013年4月22日17时11分52秒发生的内蒙古科尔沁左翼后旗M_L5.6地震，2015年12月6日8时7分48秒发生的河北丰南M_L3.4地震，2015年11月25日0时49分19秒发生的天津宝坻M_L2.5地震，2015年11月28日2时10分38秒发生的河北丰南M_L4.0地震，2015年9月14日18时10分10秒发生的河北昌黎M_L4.7地震)进行分析.通过实际测量表明，随着地震震级的增大，地震波形的低频成分相对越来越丰富.在震相到时顺序为Pg，Pb，PmP时，我们得到14组数据，得到Pg最大实测周期为0.6 s(最低实测频率为1.666 Hz)，Pb的最大实测周期为0.52 s(最低实测频率为1.923 Hz)，PmP最大实测周期为0.56 s(最低实测频率为1.785 Hz).在震相到时顺序为Pn，Pb，Pg时，我们得到34组数据，得到Pn最大实测周期为0.66 s(最低实测频率为1.515 Hz)，Pb最大实测周期为0.82 s(最低实测频率为1.219 Hz)，Pg最大实测周期为0.54 s(最低实测频率为0.925 Hz).鉴于实测数据几乎囊括了研究数据集合中的较大地震，因此得到的最低实测频率基本就是相应震中距范围内的肉眼能分辨的最低频率.因此在进行功率谱密度分析时，考虑相应误差，在震相到时顺序为Pg，Pb，PmP时，我们采用肉眼能分辨的最低参考频率Pg的为1.66 Hz，Pb的为1.92 Hz，PmP的为1.78 Hz；在震相到时顺序为Pn，Pb，Pg时，我们采用肉眼能分辨的最低参考频率Pn的为1.51 Hz，Pb的为1.21 Hz，Pg的为0.92 Hz.图 12为平原地区不同震中距范围内的Pb，Pg，Pn，PmP震相不同震源深度平均归一化功率密度分布图.同时对震源深度进行了分类，得出我们的事件的震源深度绝大部分都在20 km以内.从图 12中可以看出，P波的主要能量分布在低频部分，绝大部分在肉眼能分辨的最低频率之下，并随震中距的增加，主频最高频率逐渐降低.P波在肉眼能分辨的最低频率之上的高频部分，总体特征为归一化功率谱密度随频率的增加降低缓慢，同样随震中距的增加高频部分逐渐减少.在Pg作为初至波时，Pb震相的低频主频部分与Pg震相的低频主频部分差不多(受到包含Pg震相的影响)，但是高频主频部分频率更高(图 12a, f).Pn作为初至波时，Pn震相低频主频部分宽度比Pb宽，但是Pb高频主频部分频率相对更高(图 12b, c, d, e, k, l, m, n); 作为续至波的Pg(包含了Pb，Pn和噪声)，主频低频部分随震中距的增加与高频主频部分优势变得不明显，主频高频部分也随震中距的增加高频成分减少；PmP(包含了Pg，Pb和噪声的影响)，与作为初至波的Pg相比高频成分丰富.图 13为山区不同震中距范围内的Pb，Pg，Pn，PmP震相不同震源深度平均归一化功率密度分布图.与平原区相比，由于地壳厚度的区别，Pn出现距离不同，数据量少，因此选择了不同的震中距分析.从图 13依然可以得到与图 12一样的结果.虽然Pb，Pn，Pg，PmP的能量集中的低频部分主频与噪声的主频有很大部分的重合，但是随着震级的增加，事件波形的低频成分也增加，震中距越远能记录到事件的震级也越大.在M_L4.0左右台基好的台站宽频带地震计记录到的地震P波低频成分0.02~0.05 Hz频段信噪比都很强并广泛用于震源机制反演(林向东等, 2013, 2017; Lin et al., 2018)，本研究的地震震级在M_L2.5~5.6之间，因此即使所得的功率谱密度在低频的主频部分包含了噪声，但是它也叠加了地震信号，比主频峰值大的主频部分代表震相的主频可靠性还是强的，小于峰值的部分，可能受噪声的影响会强一些，可靠性差.

图 12

Fig. 12

图 12 平原地区不同震中距范围内的Pb，Pg，Pn，PmP震相平均归一化功率密度分布图(黑色虚线为1.0 Hz参考线，红色虚线为观测到的肉眼可识别最低频率参考线.黑色曲线和*代表震源深度为0~5 km的事件，同理蓝色为5~10 km, 紫色为10~15 km，黄色为15~20 km, 红色为20~25 km，绿色为25~30 km的事件) Fig. 12 Normalized and averaged power density spectrum of Pb, Pg, Pn, and PmP phases for different epicentml distances in the plain region. Vertical black dashed line is 1.0 Hz reference, and red dashed line is 1;he observed minimum visible frequency. Color coding indicates different fbcal depths: black, 0~5 km; blue, 5~10 km; purple, 10~15 km; yellow, 15~20 km; red, 20~25 km; and green, 25~30 km

图 13

Fig. 13

图 13 山区不同震中距范围内的Pb, Pg, Pn, PmP震相平均归一化功率密度分布图(参考线和震源深度分类与图 12同) Fig. 13 Same as Fig. 12 but for mountainous regions

Pb相对于Pg，PmP，Pn的到时顺序与震源深度、震中距有很大关系, 在Pb走时曲线与Pg，Pn接近或者相交的震中距范围，Pb是可能为初至波的.考虑到地震定位精度，最优定位结果误差，在水平向只能做到小于3 km，在垂直向约束更差，不同定位结果深度相差10 km也是很常见.鉴于实际观测数据定位结果震源深度误差大，分不同深度范围去拟合震相在不同震中距的到时顺序分辨率不会高.本研究识别Pb震相的数据时，为了保证数据可靠，根据两个震相到时较为接近时不选，同时参考确定Pn出现之后Pb震相的波形特征的原则，这可能导致了Pb走时曲线与Pg，Pn接近或者相交的震中距范围的数据缺失或者拟合分辨率不高.为了分析在不同震中距范围内, Pb的相对到时的顺序是什么的问题，考虑到定位深度误差，以及本研究所涉及的大多数地震事件震源深度在10 km左右(见图 11—13)，我们不考虑震源深度，将本研究所校核与识别的所有的Pg，Pb，PmP，Pn震相标到相应归一化的波形记录上，并以相对到时为横坐标(发震时刻为零点)，震中距为纵坐标，将所有记录进行排序，并对相应的震相进行多项式拟合(图 14).发震时刻由于定位误差，可能有0.5 s左右的到时差，但是数据是多个数据的集合，因此不影响拟合结果.在Pn出现的震中距范围，已有研究者提出依次相应的顺序为Pn，Pb，Pg(如傅淑芳, 1991；郭永霞, 2010)，但是对于Pn未出现的震中距范围，Pb的到时顺序鲜有报道.从图 14可以看出，我们所识别的Pg，Pb，Pn，PmP多项式拟合结果具有较好的线性关系，而且数据收敛性较好.在Pn出现的范围，到时顺序为Pn，Pb，Pg, 与前人的研究结果一致.在Pg作为初至波的距离范围，相应震相的到时顺序为Pg, Pb, PmP.震中距约在120~140 km，顺序较为复杂，有Pg, Pb, PmP, 也有Pn，Pb, Pg.说明在此震中距范围内震相复杂.在校核震相时，我们考虑了发生在不同区域的地震Pn出现的相对距离，根据地震发生的区域，将数据分为平原区和山区两类，并分别进行多项式拟合分析.图 15为平原区数据的拟合结果，Pn出现距离约为120 km, 可以看出，拟合线性度很好.图 16为山区数据进行多项式拟合的结果，Pn出现距离约为140 km, 数据相对较少，拟合时曲线相对弯曲，但依然可以从其拟合结果看出与图 14和15相同的趋势.在Pn出现范围，本研究涉及震相依次到时顺序为Pn, Pb, Pg；在Pg作为初至波的范围，到时顺序为Pg，Pb, PmP.在Pb走时曲线与Pg，Pn接近或者相交的震中距范围, 由于观测误差的影响，实际观测数据分辨率并不高(图 14—16)，而这段震中距范围内的Pb相对位置需要结合震源深度误差、理论到时(见5.4节)、波形特征等来确定Pb是否为初至波.

图 14

Fig. 14

图 14 首都圈及邻区Pn, Pb, Pg, PmP在不同震中距相对到时顺序(红色线为Pg到时，黄绿色线为Pb到时，黄色线为PmP到时，紫色线为Pn到时) Fig. 14 Arriving sequence of Pn, Pb, Pg, and PmP phases at different epicentral distances in the Capital Circle region around Beijing and adjacent areas. Color lines show arrival times of Pg (red), Pb (yellow-green), PmP (yellow), and Pn (purple), respectively

图 15

Fig. 15

图 15 平原地区Pn, Pb, Pg, PmP在不同震中距相对到时顺序(拟合彩色线代表同图 14) Fig. 15 Same as Fig. 14 but for plain regions

图 16

Fig. 16

图 16 山区Pn, Pb, Pg, PmP在不同震中距相对到时顺序(拟合彩色线代表同图 14) Fig. 16 Same as Fig. 14 but for mountainous regions

利用观测到的Pb震相数据，我们可以根据台站到震中的射线路径在空间上的覆盖情况来判断康拉德界面在横向的分布情况，通过最小二乘法拟合所有Pb数据的走时，可以粗略估计研究区域的Pb传播速度.同时利用本研究所得的Pg, Pb, Pn, PmP震相数据，根据联合反演方法得到该区域的一维速度模型，并得到康拉德界面的深度.有了速度模型，我们则可以反过来计算震相相应的理论走时去验证分析的观测数据是否准确.

图 17是首都圈及邻区Pb的射线路径图，其中，黑色直线为Pb射线路径，红色实心三角形为识别出观测到Pb震相的台站，黄色空心三角形为没有识别出观测到Pb震相的台站.我们分析了186个台站的数据，其中157个台站我们识别出观测到了Pb震相, 29个台站没有识别出是否观测到Pb震相.从图中可以看出，首都圈及邻近地区大部分台站都识别出记录到了Pb，只有在平原沉积层较厚的地区或者海边的少数台站没有识别到记录Pb的情况.在沉积层较厚区域，或者海边的台站，由于沉积层场地响应以及海浪影响，常常导致观测记录噪声很大，增加了震相识别的困难.射线密度受到北京市地震局测震台网接入台站和河北省地震台网接入台站分布，以及地震震中分布的影响.地震主要分布在张家口—渤海断裂带，山西断陷带，以及邢台等地区.在华北平原中东部地区，地震相对较少，所以中东部地区射线密度相对稀疏.在台网边缘，接入台站相对少，再受到震中分布的影响，射线也相对稀疏.但是从总体上看，考虑到震中及台站分布的影响，首都圈及邻区Pb射线路径整体覆盖较好，可以得出，首都圈及邻区的康拉德界面横向分布是连续的.

图 17

Fig. 17

图 17 首都圈及邻区Pb射线路径覆盖情况图 Fig. 17 Pb ray path coverage in the Capital Circle region around Beijing and adjacent areas

本研究得到了1153个Pb的走时数据，通过利用最小二乘法拟合，我们可以初略地得出Pb的实际观测速度.图 18是在数据分布相对均匀的震源距范围内利用1143个Pb走时进行拟合, 得到Pb的拟合速度为7.0443 km·s^-1, 残差为0.82853.拟合残差受到定位的震中误差、震源深度误差、以及到时读取误差等的影响，精度可能相对较低，但是根据数据分布及其拟合几乎为一条直线的情况，我们可以得出，首都圈及邻区Pb的平均传播速度在7.0 km·s^-1左右.

图 18

Fig. 18

图 18 Pb速度线性拟合结果(*代表不同事件的到时) Fig. 18 Pb velocity linear fitting results. Stars show the arrivals of individual events

图 19是我们利用本研究重新校核过的Pg, Pb, Pn, PmP等震相走时结果，并对所有地震事件进行了重新定位，采用基于遗传算法的联合反演方法(万永革和李鸿吉，1996)，反演首都圈及邻区一维速度模型，从而确定康拉德界面的深度和相对速度.在联合反演时，考虑到首都圈地区地震台站分布较密，地震定位得到的发震时刻、深度的精度问题，我们设定经纬度允许变化范围为±0.05°, 深度允许变化范围为±10 km，发震时刻允许变化范围为±1.5 s.根据祝志平等(1997)物探的结果，我们将地壳分为5层的初始模型(见表 1).通过20000次迭代，得到了最佳模型(表 2，图 19d)，残差为0.004102.图 19a为迭代种群变化及平均残差分布，基本在迭代1万次以后拟合相对稳定了.图 19b为搜索发震时刻变化分布图，变化主要集中在±0.7 s.图 19c为重定位后地震深度分布图，优势分布在5~20 km.图 19d为最佳模型速度分布图.根据反演得到的模型结果，以及实际拟合的Pb速度结果，我们可以得出首都圈及邻区康拉德界面是第三层与第四层的分界面，平均深度在23 km左右，康拉德界面上P波速度反演结果为7.06 km·s^-1，与直接用Pb观测走时拟合到的速度7.04 km·s^-1非常接近.考虑到观测误差并结合人工地震剖面得到的速度结果(祝志平等，1997)，我们认为康拉德界面的平均速度应在6.8~7.0 km·s^-1之间.需要特别说明的是，在反演过程中我们不考虑康拉德界面的倾斜，所以得到的结果为区域的平均结果.

图 19

Fig. 19

图 19 基于遗传算法联合反演速度模型结果 (a)迭代种群变化及平均误差分布；(b)搜索发震时刻变化分布；(c)重定位后地震深度分布；(d)最佳模型速度分布. Fig. 19 Velocity models from joint inversion based on the genetic algorithm (a) Variation of iteration seed groups and average errors; (b) Variations of occurrence times searched; (c) Depths of relocated events; (d) Optical velocity model.

表 1 (Table 1)

表 1 用于联合反演的首都圈及邻区初始速度模型 Table 1 Initial velocity model used in joint inversion for the Capital Circle region around Beijing and adjacent areas 厚度变化范围(km) P波速度变化范围(km·s-1) 平均沉积层 0~4 2.0~5.4 第二层 5~15 5.5~6.2 第三层 5~15 6.2~6.4 第四层 5~15 6.4~7.3 上地幔 1000- 7.9~8.3

表 1 用于联合反演的首都圈及邻区初始速度模型 Table 1 Initial velocity model used in joint inversion for the Capital Circle region around Beijing and adjacent areas

表 2 (Table 2)

表 2 基于联合反演得到的首都圈及邻区最佳速度模型 Table 2 Optical velocity model derived from joint inversion for the Capital Circle region around Beijing and adjacent areas 厚度(km) P波速度(km·s-1) 平均沉积层 0.5 3.60 第二层 11.1 6.11 第三层 11.5 6.26 第四层 10.0 7.06 上地幔 — 7.94

表 2 基于联合反演得到的首都圈及邻区最佳速度模型 Table 2 Optical velocity model derived from joint inversion for the Capital Circle region around Beijing and adjacent areas

根据联合反演得到的P波最佳模型，结合经验公式V_S=V_P/1.732给出S波速度，以及Crust 2.0(Bassin et al., 2000)给出的相对密度建立模型.我们用Pyrocko(Heimann et al., 2017)方法模拟Pg, Pb, PmP, Pn的传播路径(图 2).分别用震源深度为5 km，10 km，15 km，20 km, 23.5 km计算走时曲线.结果表明：震源深度为5 km时，Pb, Pn, Pg走时曲线的交叉点几乎一致，即在震中距约为140 km处，震中距约＜140 km时，震相到时顺序为Pg, Pb, PmP，震中距约＞140 km时，震相到时顺序为Pn，Pb，Pg.震源深度为10 km时(见图 2和图 20)，Pb与Pg走时曲线交叉点约在震中距为120 km处，与Pn的交叉点约在135 km处，这说明震中距约＜120 km时，震相到时顺序为Pg, Pb, PmP；震中距约＞120 km且＜135 km时，震相到时顺序为Pb, Pg，PmP；震中距约＞135 km时，震相到时顺序为Pn, Pb, Pg.震源深度为15 km时，Pb与Pg走时曲线交叉点约在震中距为110 km处，与Pn走时曲线交叉点约在震中距为130 km处.说明震中距约＜110 km时，震相到时顺序为Pg，Pb, PmP；震中距约＞110 km且＜130 km时，震相到时顺序为Pb, Pg，PmP；震中距约＞130 km时，震相到时顺序为Pn，Pb, Pg.震源深度为20 km时，Pb与Pg走时曲线交叉点在震中距约为80 km处，与Pn走时曲线交叉点约在125 km处，即震中距约＜80 km时，震相到时顺序为Pg，Pb, PmP；震中距约＞80 km且＜125 km时，震相到时顺序为Pb, Pg, PmP；震中距约＞125 km时，震相到时顺序为Pn, Pb, Pg.震源深度在下地壳顶部时, 没有Pg震相，震中距约＜112 km时，Pb作为初至波，震中距约＞112 km时，Pn作为初至波.图 20中不同颜色的“*”表示实际观测到的不同深度和不同震中距的Pg, Pb, PmP, Pn等震相到时数据，与“*”对应颜色的曲线为震源深度为10 km时相应震相的理论走时曲线.从图 20理论走时曲线可以看出，在震中距约为120~135 km时，理论震相到时顺序为Pb, Pg, PmP，而我们的观测数据在这段距离内的到时顺序较为复杂，有震源深度的影响，也有定位精度的影响等，说明分辨率不高.在其他震中距范围与实际观测的结果基本一致，进一步说明了我们结果的可靠性.由于受到震相到时识别误差和震源位置定位误差等的影响，实测数据有相对离散性.实际区域速度结构在三维空间是不均匀的，在震中距比较近时，速度结构的各向异性对实际走时影响较大.我们的模型是一维的，所以在震中距比较近时理论走时曲线和实际观测数据稍微有些偏离，且我们的理论走时仅是震源深度为10 km的计算结果，而实测数据是包含了各种震源深度的结果，所以两者对比会有一定的偏差.理论计算和实测结果表明，震源深度大于5 km时，在Pb走时曲线与Pg，Pn接近或者相交的震中距范围，我们需要结合波形特征、震中距、震源深度等情况来判定Pb是否为初至震相.

图 20

Fig. 20

图 20 不同震源深度和震中距离的实测Pg，Pb，Pn，PmP到时与基于联合反演所得的4层速度模型计算得到的Pg，Pb，Pn，PmP四种震相在震源深度为10 km时的理论走时 Fig. 20 Observed arrival times of Pg, Pb, Pn and PmP phases compared with the calculated theoretical travel times with a 10 km deep source and the 4-layer crustal velocity model from the joint inversion of this study

需要特别说明的是，大概在120 km之后(Pn出现的范围)PmP理论走时依然存在，但是稍微比Pb慢，随着震中距的增加二者逐步趋于一致.但是实际观测中，在Pn出现距离内PmP不发育，很难识别，基本没有相应的观测报道.根据PmP的出射角度，以及实际观测的PmP震相的振幅和频率大多数相对比Pg大或高，这段距离内如果有PmP出现，同理它的特征即振幅和频率大多数也应该相对比Pg大或高.但实际情况几乎没有观测到，这也从另外一个侧面说明我们的Pb识别数据是可靠的.

关于Pb震相的识别与时频特征的相关文献较少.由于Pb震相识别比较困难，在平时区域地震台网编目工作中较少有台网进行系统分析.也许有些台站记录到了Pb，但是由于震源性质、传播路径、台基噪声等方面的影响，也会被直接淹没在噪声中，我们无法进行识别.尤其是Pb的振幅在比Pn或者Pg振幅更小的情况下，即使有了频率的稍微改变，也很难识别.只有频率改变相对明显，或者振幅差异相对比较明显的情况下，我们才能更好地保证Pb震相读取的可靠性.Pb震相和Pn震相一样，有相应的盲区和出现范围，即使在出现的范围内，也并不是每个台都能记录到.本研究中为了保证震相识别的准确性，我们选择只有P波内至少出现3种震相时的数据，才进行Pb震相的挑选和识别工作.同时结合不同的构造分区，对Pn出现的最小震中距区别对待.如山西和内蒙古地区山区采用Pn出现的最小震中距为140 km, 平原地区Pn出现的最小震中距为120 km, 即使基于此，在地震波传播路径既经过山区，又经过平原的记录，也许在震中距为120~140 km内的少部分震相分析数据可能还会存在误差.通过不同类型的震源，不同台站记录同一事件的Pb震相，以及同一台站记录到的不同事件的Pb震相在时频域内的特征可以看出，我们通常容易识别的Pb的振幅或频率高频部分相对比Pn和Pg大或高(如焦煜媛等，2017)，也有振幅变小或者频率变化不明显的情况，这可能与震源机制、台站方位、场地响应、仪器类型等方面有关.对于Pb震相并未表现出十分明显的强振幅(如图 4c)，而是有较为明显的频率特征变化情况，根据孙伟家等(2018)对莫霍面以及岩石圈内部界面(MLD)的研究结果，认为康拉德界面在这些地区的速度对比并不强，而可能是一个物质和速度渐变的界面.对于可以看到有明显特征的Pb震相(如图 4f)，可能意味着这条射线路径经过的康拉德界面速度对比较强.焦煜媛等(2017)认为强烈的康拉德界面的反射波产生于具有高阻抗差的康拉德界面，并根据观测数据提出甘东南地区的上、下地壳界面具有高阻抗差.通过功率谱密度和肉眼识别分析结果表明，P波的主要能量集中在低频部分.即使所得的功率谱密度在低频的主频部分包含了噪声，但是它也叠加了地震信号，比主频峰值大的主频部分代表震相的主频可靠性还是强的，小于峰值的部分，可能受噪声的影响会强一些，可靠性差.通过功率谱密度可以看出都属于P波的Pn, Pb, Pg, PmP四种震相有很强的共性，区别就是在传播路径上的不同，从而带来了频率或观测记录的周期上的小幅度的差别.

在Pn出现的震中距范围，本研究结果与之前有研究指出主要震相的到时顺序依次为Pn，Pb，Pg(如傅淑芳和刘宝诚，1991；郭永霞等，2010)相同，但是在Pn的盲区，Pb到时顺序鲜有报道，本研究利用多项式拟合实际观测数据相对走时曲线，同时基于理论计算，综合分析结果表明：鉴于震源深度对Pb到时顺序影响较大，以及深度定位结果精确度差的情况，我们得考虑区域地壳厚度横向不均匀、震中距、震源深度等情况并结合波形特征(如Pb震相的频率比Pg的频率高)，来判定Pb是否为初至震相.震源深度在康拉德界面之下时，在Pn未出现的震中距范围，Pb将直接为初至波，没有Pg震相.通过本研究联合反演所得模型计算理论走时结果和实际观测结果基本一致进一步证明了我们结果的可靠性.

关于华北地区康拉德界面的P波传播速度、康拉德界面深度、以及康拉德界面横向分布，前人有不少研究成果.滕吉文等(1979)用工程爆破和海城地震资料得到京津唐张地区的综合四层地壳模型并指出康拉德界面P波速度为6.83 km·s^-1，平均深度为19.73 km；万永革和李鸿吉(1996)利用联合反演的遗传算法反演了京津唐张地区的速度模型并指出康拉德界面P波速度为6.6216 km·s^-1，康拉德界面平均深度为18.36 km；赵仲和(1983)使用1979年记录的43个地震，共836个P波及S波震相到时，联合测定出了北京地区的速度模型MDBJ81，并指出康拉德界面平均P波速度为7.03 km·s^-1，平均深度为24.9 km(研究区涉及到了山区和平原区).IASP91模型(Kennett, 1993)在没有考虑沉积层厚度和速度的情况下指出康拉德界面的P波传播速度为6.5 km·s^-1，深度为20 km.BLOC86模型(来自中国地震台网中心赵仲和)同样在没有考虑沉积层厚度和速度的情况下，指出康拉德界面的P波传播速度为6.51 km·s^-1，深度为15 km.目前区域台网测震监测系统Jopens(广东省地震局开发)下华北地区使用的区域模型(郑秀芬等，1999)同样未考虑沉积层厚度和速度，康拉德界面的P波传播速度为6.51 km·s^-1, 深度为15 km，与BLOC86模型一致，但是康拉德界面之上的一层速度比BLOC86 (来自中国地震台网中心赵仲和)模型高0.31 km·s^-1.李松林等(2001)利用人工地震测深剖面联合反演首都圈地壳三维速度结构指出，首都圈地区康拉德界面最浅部位于东部的三河、宝坻一带, 深度为20~21 km, 北京附近为22 km, 至张家口北西加厚到23~24 km.在北京—涿鹿间深度较深, 为24 km.易县至蔚县北为一康氏界面隆起区, 深度为21 km左右.由蔚县向涞源方向康拉德界面迅速变深, 至研究区边缘约为26 km.其他地区康拉德界面大致在22~24 km之间变化, 且存在着局部起伏.郑天愉等(2015)指出了华北地区康拉德界面深度在16~26 km范围内变化，平原区康拉德界面相对较浅，山区较深.本研究使用369个M_L≥2.5地震事件的波形资料，重新分析震相，识别出1153条Pb震相，用最小二乘法拟合出康拉德界面的P波速度为7.0443 km·s^-1.在没有考虑康拉德界面倾斜的情况下，结合重新校正的Pg, Pn, PmP, Pb震相到时，通过联合反演遗传算法(万永革和李鸿吉, 1996)反演得到的康拉德界面的P波速度为7.0622 km·s^-1，平均深度为23 km.无论是直接拟合的康拉德界面P波传播速度，还是通过联合反演最终得到的P波速度，考虑到误差因素，与赵仲和(1983)、滕吉文(1979)的结果相对一致，比万永革和李鸿吉(1996)、李松林等(2001)、王丽婵(2016)等人的结果相对偏高，可能是因为我们使用了更多的Pb数据，研究区域的范围不同或者是由于人工地震射线传播路径与天然地震传播路径不一致等原因引起的.范玉兰等(1990)分析天然地震和人工爆破的基础上，得出华南地区康拉德界面的平均P波速度为6.88 km·s^-1.本文所得首都圈地区的结果为7.0 km·s^-1左右.考虑到地区接壤以及相关误差，华北和华南康拉德界面的平均P波速度可能具有一定的连续性或相似性，我们的结果或许有一定的合理性.康拉德界面深度则与前人结果(如赵仲和, 1983；滕吉文等，1979)等基本一致.

康拉德界面并不是一个连续的界面, 在大陆有些地区存在, 有些地区不存在.不存在的地区也可能是因为康拉德界面是一个较厚的过渡区而导致难以观测到.这种情况可能与Sun和Kennett(2017)认为华北克拉通西部(鄂尔多斯)下方的中岩石圈不连续面实质上是岩石圈到软流圈的过渡区域的起始位置的结果类似.虽然郑天愉等(2015)给出了华北地区康拉德界面的深度等值线分布, 具有重要的科学意义，但是他们的结果是通过人工地震剖面插值后得到的连续面.王丽婵等(2016)利用96个Pb数据的分布情况指出，康拉德界面在河北地区是不连续的.本研究通过1153个Pb震相的射线路径直接给出了康拉德界面在首都圈及邻区的分布情况，用观测数据直接证明了康拉德界面在首都圈地区的分布是连续的.

在结合前人的研究结果的基础上, 本研究通过分析大量的可靠数据首次详细地给出了首都圈地区Pb震相典型时频特征，为区域台网人员如何识别该震相提供了技术标准，对丰富首都圈地区台网观测报告震相产出，提高地震定位精度和确定康拉德界面等相关研究具有重要意义.

我们通过分析得出：在首都圈及邻区，从肉眼可识别时频域内的特征来看，通常容易识别的Pb的振幅或频率高频部分相对比Pn和Pg大或高，也有振幅变小或者频率变化不明显的情况，这可能与震源机制、台站方位、场地响应、仪器类型等方面有关.通过时频分析，功率谱密度和肉眼识别分析结果表明，P波的主要能量集中在相对低频部分，都属于P波的Pn，Pb，Pg，PmP四种震相有很强的共性，区别就是在传播路径上的不同，从而带来了频率或观测记录的周期上的小幅度的差别.在Pg作为初至波时，Pb震相的低频主频部分与Pg震相的低频主频部分差不多(受到包含Pg震相的影响)，但是高频主频部分频率更高(图 12a, f).Pn作为初至波时，Pn震相低频主频部分带宽比Pb宽，但是Pb高频主频部分频率相对更高(图 12b，c，d，e，k，l，m，n).通过多项式拟合实际观测数据相对走时曲线，与理论计算走时曲线分析表明：首都圈及邻区，震源深度约≤5 km时，Pb只能作为续至波出现，即震中距约＜140 km时，震相到时顺序为Pg, Pb，PmP，震中距约＞140 km时，震相到时顺序为Pn，Pb，Pg.随着震源深度从5 km增至20 km左右，Pb理论走时曲线与Pg，Pn接近或者相交的震中距范围由震中距约以140 km处为起点，逐渐向震中距减小的方向移动并扩大范围，即增加至最大震中距约在80~125 km范围内.在这个随深度变化而变化的震中距范围内，Pb理论震相作为初至震相, 到时顺序为Pb，Pg，PmP.鉴于Pb相对到时顺序与震源深度的关系，以及深度定位结果精确度差的情况，在震中距约在80~140 km范围内时，我们得考虑区域地壳厚度横向不均匀、震中距、震源深度等情况并结合波形特征(如大多数Pb震相的频率比Pg的频率高)，来判定Pb是否为初至震相.震源深度在下地壳顶部时, 没有Pg震相，震中距约＜112 km时，Pb作为初至波，震中距约＞112 km时，Pn作为初至波.Pb震相平均传播速度约在7.0 km·s^-1，康拉德界面平均深度约为23 km；Pb射线的分布情况直接证明了康拉德界面在首都圈的分布是连续的(图 17).基于本研究联合反演所得模型的基础上计算的理论走时结果和实际观测对比分析结果也说明了我们结果的可靠性(图 20).

在对分析出来的首都圈地区丰富的Pb震相到时数据的利用上，本研究通过分析和反演仅给出了康拉德界面下地壳的平均的速度及其平均深度，结果相对粗糙，下一步我们将致力于利用Pb观测数据，进一步研究并给出康拉德界面在该区域的三维分布.