基于极限梯度提升树的实时P波初至自动拾取方法

A Real-Time Automatic Picking Method of P-Wave First-Arrival Based on Extreme Gradient Boosting Trees

摘要

针对传统P波初至自动拾取方法抗干扰能力弱、拾取精度低的问题，提出1种基于极限梯度提升树（XGBoost）的实时P波初至自动拾取方法。首先，选择有助于区分地震信号与背景噪声的4种特征参数作为模型的输入，以降低模型的复杂度；其次，构建P波初至自动拾取XGBoost模型，并对模型进行训练和测试；最后，通过与目前地震预警中常用的P波初至实时识别方法进行对比，验证模型的有效性。结果表明：所提方法在±0.5 s误差范围内的拾取样本占比为93.3%，优于能量周期双参数（EDP-Picker）方法和短时/长时平均比（STA/LTA）方法，二者拾取样本占比分别为91.9%和83.6%；当误差超出±0.5 s时，XGBoost方法的提前和滞后触发样本占比分别为4.27%和5.26%，而EDP-Picker的相应比例分别为5.03%和6.50%，STA/LTA的相应比例分别为5.39%和1.71%。相较于2种传统方法，XGBoost方法的综合性能显著提升，且具有更高的识别精度和更强的抗干扰能力，能够更稳定地适应复杂场景下的拾取需求。

Abstract

To address the issues of weak anti-interference capability and low picking accuracy in traditional automatic picking methods of P-wave first-arrival, a real-time automatic picking method of P-wave first-arrival based on Extreme Gradient Boosting Trees (XGBoost) is proposed. First, 4 characteristic parameters that were conducive to distinguishing seismic signals from background noise were selected as the model's input to reduce complexity of the model. Subsequently, an XGBoost model for automatic picking of P-wave first-arrival was constructed, trained, and tested. Finally, the model's effectiveness was validated by comparing it with commonly used real-time identification methods of P-wave first-arrival in current earthquake early warning systems. The results indicated that the proposed method achieved a picking sample proportion of 93.3% within an error range of ±0.5 s, outperforming both the Energy Periodic Dual Parameter Picking (EDP-Picker) method and the Short-Term Average/Long-Term Average (STA/LTA) ratio method, the picking sample proportion of which was 91.9% and 83.6%, respectively. When the error exceeded ±0.5 s, the proportion of early-triggered and lag-triggered samples of the XGBoost method was 4.27% and 5.26%, respectively, while the corresponding proportion of EDP-Picker was 5.03% and 6.50%, respectively, and that of STA/LTA was 5.39% and 1.71%, respectively. Compared with 2 traditional methods, the overall performance of XGBoost method was significantly enhanced, along with superior identification accuracy and stronger anti-interference capability, meeting picking demands with greater stability in complex scenarios.

Graphical abstract

关键词

现地预警 / 地震紧急处置 / P波初至自动识别 / 极限梯度提升树 / 机器学习

Key words

On-site early warning / Earthquake emergency disposal / Automatic identification of P-wave first-arrival / Extreme Gradient Boosting Trees (XGBoost) / Machine learning

P波初至的精确拾取是实时地震预警系统的关键技术环节，其准确性直接影响减灾效果。在地震预警系统中，P波初至拾取精度不仅决定了实时地震定位的可靠性^［1］，还制约着震级估计的准确性^［2］，微秒级的拾取误差都可能导致现地地震动预测出现显著偏差，进而影响列车紧急制动决策^［3-5］。目前主流的自动拾取方法可分为3类：基于信号能量突变特征^［6-7］、基于信号周期突变特征^［8］以及基于人工智能^［9-11］的P波初至自动识别方法。

基于能量突变特征的P波初至自动识别方法包括长短时平均（Short-Term Average/Long-Term Average，STA/LTA）方法^［12-13］、高阶矩法^［14-15］等。STA/LTA方法因计算效率高被广泛应用于实时地震预警系统。然而，该方法存在2个显著局限性：一是对高频干扰极为敏感，容易产生误触发；二是在大地震远场记录中，由于P波能量较弱且与背景噪声差异较小，STA/LTA比值变化不明显，往往导致P波初至拾取滞后。地震信号和背景噪声间的差异不仅体现在能量方面，还体现在信号周期方面。文献［8］借鉴实时P波卓越周期的计算方法，将卓越周期用于P波初至自动拾取。然而，由于该方法忽略了信号能量突变的特征，和基于能量的识别方法相比，二者的P波初至自动识别精度相当。为了克服基于单一特征的P波初至自动拾取精度不高的缺陷，文献［16］提出了1种基于能量与周期双参数融合的P波初至自动拾取方法（EDP-Picker），该方法首先利用周期和能量突变特征参数进行粗略拾取，然后结合AIC函数进行P波初至精确拾取，和基于能量或周期单参数的拾取方法相比精度有了明显提升。然而，尽管EDP-Picker方法具有较高的拾取精度，但参数阈值的选取、判别条件的建立仍然通过经验确定，在一些特殊工况（如小震远场记录、低频背景干扰等）下存在误触发较多的现象。

从背景噪声信息中拾取P波初至实质上属于分类问题。近年来，机器学习技术在解决分类问题时体现出了明显的优势，大量机器学习方法也被应用于P波初至自动识别研究中。如支持向量机、卷积神经网络、双向LSTM等方法被应用于P波初至自动拾取。然而，部分基于机器学习的拾取方法需要较长的数据时长，不适用于实时预警系统；部分方法由于使用的网络深度过大，会消耗实时预警系统的大量计算资源。因此，建立1种兼具时效性和精度的P波初至自动拾取模型仍然是目前需要解决的关键问题。文献［17］提出了1种基于极限梯度提升树（Extreme Gradient Boosting，XGBoost）的机器学习算法，在分类问题中具有更好的性能。此外，所提算法具有效率高、结构简单等特点，适合应用于实时系统。

本研究提出1种基于XGBoost的实时P波初至自动拾取方法。为降低模型复杂度，选取4种能够有效区分地震信号与背景噪声的特征参数作为输入，并采用日本K-NET台网的9 978条竖向强震动记录进行训练，同时利用KiK-net台网的6 726条记录进行对比验证，可为重大工程地震预警及紧急处置提供参考。

1 模型构建

1.1　XGBoost算法原理

XGBoost是1种以梯度提升框架为基础的算法，通过集成多个弱分类器，构建性能更优的强分类器，从而提升模型的整体表现。对于给定有监督学习数据集{\left\{\left({\mathit{x}}_i,y_i\right)\right\}}_{i=\mathrm{1}}^N，则XGBoost分类模型可表示为多个分类决策树的加权和，其中：{\mathit{x}}_i为输入特征；y_i为分类标签；N为样本数量。该模型的预测值{\widehat{y}}_i为

{\widehat{y}}_i={\displaystyle \sum _{k=\mathrm{1}}^K}{f}_k({\mathit{x}}_i)\forall f_k\in \mathrm{Q}（1）

式中：f_k为第k棵分类树的分类函数，其与树的结构和叶子权重值相关；K为分类树的数量；Q为所有分类树的集合。

为防止模型过于复杂导致的过拟合问题，XGBoost模型的目标函数L（\varphi）由损失函数和正则项组成

L(\varphi )={\displaystyle \sum _{i=\mathrm{1}}^N}l(y_i,{\widehat{y}}_i)+{\displaystyle \sum _{k=\mathrm{1}}^K}\Omega (f_k)（2）

式中：l（y_i,{\widehat{y}}_i）为损失函数，表示真实值y_i与预测值{\widehat{y}}_i的误差；\mathrm{\Omega }（f_k）为正则项，用于控制模型的复杂度以防止过拟合。

1.2　特征与标签

为使模型满足实时性要求，采用固定时窗长度（L=2 s）计算特征参数x，将其时窗末端与加速度记录的采样时刻对齐，标签则根据时窗末端样本的属性确定。若该样本为背景噪声则标签赋值为0，若为地震信号，则赋值为1，特征与标签提取示意图如图1所示。以加速度记录的采样间隔为时间步长滑动时窗，可获取任意1条加速度记录后L秒内各采样点对应的特征及标签值。此时，每个时刻的特征-标签对即构成1个训练样本。在模型应用阶段，当实时观测数据的长度达到时窗长度L时，即可对每个样本进行预测。为减少模型训练过程中的数据量，选取每条记录开始后2 s至P波初至后3 s时间窗内的特征和标签用于训练模型。

1.2.1　特征

Lu等^［16］的研究表明，基于能量和周期参数的P波初至自动识别方法具有较好的识别效果。该方法的能量特征采用1 s滑动窗口计算信号方差的对数能量，通过当前能量与历史最小能量的相对增量直接量化P波能量相对于噪声的绝对变化，避免传统STA/LTA等方法对长窗口的依赖。使用的周期特征可以捕捉P波到达时的频率突变特征。因此，选取3种特征参数，分别为：卓越周期T^{\mathrm{p}\mathrm{d}}；卓越周期梯度∆T^{\mathrm{p}\mathrm{d}}；加速度能量梯度∆E。为抑制高频干扰导致的误触发问题，通过实时积分将加速度记录转换为速度记录，并从中提取加速度能量梯度参数作为补充特征。各特征参数的计算方法如下。

1）卓越周期

Hildyard提出了1种实时算法来计算卓越周期^［8］，并将该算法应用于P波初至捡拾的研究中。Hildyard卓越周期T_i^{\mathrm{p}\mathrm{d}}的计算式为

T_i^{\mathrm{p}\mathrm{d}}=\mathrm{2}\mathrm{\pi }\sqrt[]{\frac{X_i}{D_i+D_{\mathrm{s}}}}（3）

其中，

X_i=\alpha X_{i-\mathrm{1}}+v_i^{\mathrm{2}}

D_i=\alpha D_{i-\mathrm{1}}+{\left(\partial v/\partial t\right)}_i^{\mathrm{2}}

式中：v_i为第i个样点的速度幅值；\mathrm{\alpha }为平滑系数，\alpha ={\mathrm{e}}^{\mathrm{l}\mathrm{n}\left(\mathrm{0.1}\right)/\left({\tau }_{\mathrm{w}}/∆t\right)}；D_{\mathrm{s}}为1个常数，D_{\mathrm{s}}=\frac{\mathrm{4}{\mathrm{\pi }}^{\mathrm{2}}〈v_i^{\mathrm{2}}〉}{{\tau }_{\mathrm{m}\mathrm{x}}^{\mathrm{2}}}，当D_i趋于无穷小时用于保证T^{\mathrm{p}\mathrm{d}}的取值在正常范围内；∆t为加速度记录的采样间隔；{\tau }_{\mathrm{w}}为调整T^{\mathrm{p}\mathrm{d}}的取值的平滑系数，{\tau }_{\mathrm{w}}值越大，T^{\mathrm{p}\mathrm{d}}曲线越平滑，此处取值为4.5；{\tau }_{\mathrm{m}\mathrm{x}}为控制D_{\mathrm{s}}变化的平滑系数，此处取值为0.019。

2）卓越周期梯度

卓越周期梯度∆T^{\mathrm{p}\mathrm{d}}是衡量卓越周期变化程度的参数，当P波卓越周期与背景噪声卓越周期存在较大差异时，该参数能够很好地反映这一变化过程。其计算式为

∆T_i^{\mathrm{p}\mathrm{d}}=\left|\frac{T_i^{\mathrm{p}\mathrm{d}}-T_{i-n_{\mathrm{0.05}}}^{\mathrm{p}\mathrm{d}}}{\mathrm{0.05}}\right|（4）

式中：n_{\mathrm{0.05}}为0.05 s内加速度记录的样点数。

3）加速度能量梯度

为凸显P波到达时刻地震动幅值或能量的变化特征，研究者提出STA/LTA及高阶矩等特征参数。然而，这些参数仅在P波到达时刻出现突变，在P波到达后又重新回到正常水平，从而导致漏检拾的可能。针对这一问题，引入加速度能量梯度参数∆E。该参数值在P波到达前保持在较低水平，而P波到达后保持在较高水平，从而避免漏检拾的发生。样本i的加速度能量E_i和加速度能量梯度∆E_i分别为

E_i=n_{\mathrm{w}}\mathrm{l}\mathrm{g}\left[\mathrm{v}\mathrm{a}\mathrm{r}\left(a_{i-n_w+\mathrm{1}},a_{i-n_w+\mathrm{2}},\cdots ,a_i\right)\right]（5）

∆E_i=\left|\frac{E_i-E_{\mathrm{m}\mathrm{i}\mathrm{n}}}{E_{\mathrm{m}\mathrm{i}\mathrm{n}}}\right|（6）

式中：n_{\mathrm{w}}为时窗长度，取值为2 s；a_i为样本i的加速度幅值，\mathrm{g}\mathrm{a}\mathrm{l}；\mathrm{v}\mathrm{a}\mathrm{r}（·）为方差算子；E_{\mathrm{m}\mathrm{i}\mathrm{n}}为样本i之前通过实时观测记录计算得到的加速度能量最小值。

4）三特征之和

由于上述3个参数在P波到达后均呈现增大趋势，将3个参数相加得到复合参数F_{\mathrm{s}\mathrm{u}\mathrm{m}}，该处理使得P波到达时刻的特征变化幅度更为显著。单条加速度记录以及上述参数的特征变化如图2所示，图中红色直线表示P波到达时刻。

1.2.2　标签

采用滑动时间窗法提取特征参数，时间窗以系统最小采样间隔（0.01 s）为步长进行滑动。每1组特征参数对应1个标签，将记录开始后2 s至P波初至前的时间窗标记为0（表示P波未到达），将P波初至时刻及之后3 s内的时间窗标记为1（表示P波已到达）。

1.3　特征参数优选

为剔除特征参数中可能存在的冗余，构建15个P波初至自动识别模型，每个模型均以上述4个特征参数组合作为输入，具体组合见表1。AUC值是机器学习中用于评估二分类模型性能的重要指标，其定义为 ROC曲线（Receiver Operating Characteristic Curve）下方的面积，该值越大表明模型的分类性能越好。以AUC值作为标准来选择最佳的参数组合，其随不同参数组合变化情况如图3所示。当4个参数全部被选用时，模型得到的AUC值最高，表明该组合为最优参数组合。

1.4　P波初至精确识别

由于XGBoost模型需要特征参数达到特定阈值才能判定P波初至，其识别结果通常会比真实P波初至延迟若干个样点。为获得精确的P波初至，在XGBoost模型触发后，在触发时刻至前3 s的时窗内采用AIC函数进行精确拾取，AIC最小值点即为精确的P波到时。AIC函数的计算式为

V_k^{\mathrm{A}\mathrm{I}\mathrm{C}}=k\mathrm{l}\mathrm{g}\left[\mathrm{v}\mathrm{a}\mathrm{r}\left({\mathit{Z}}_{\mathrm{1}}\right)\right]+\left(m-k-\mathrm{1}\right)\mathrm{l}\mathrm{g}\left[\mathrm{v}\mathrm{a}\mathrm{r}\left({\mathit{Z}}_{\mathrm{2}}\right)\right]（7）

其中，

{\mathit{Z}}_{\mathrm{1}}=\left[a_{\mathrm{1}}{T}_{\mathrm{1}}^{\mathrm{p}\mathrm{d}},a_{\mathrm{2}}{T}_{\mathrm{2}}^{\mathrm{p}\mathrm{d}},\cdots ,a_k{T}_k^{\mathrm{p}\mathrm{d}}\right]

{\mathit{Z}}_{\mathrm{2}}=\left[a_{k+\mathrm{1}}{T}_{k+\mathrm{1}}^{\mathrm{p}\mathrm{d}},a_{k+\mathrm{2}}{T}_{k+\mathrm{2}}^{\mathrm{p}\mathrm{d}},\cdots ,a_{k+m}{T}_{k+m}^{\mathrm{p}\mathrm{d}}\right]

式中：V_k^{\mathrm{A}\mathrm{I}\mathrm{C}}为第k个样点对应的AIC值；m为时间窗口内的总样本点数；a_k为第k个样点在垂直方向的地震动加速度幅值；T_k^{\mathrm{p}\mathrm{d}}为第k个样点所对应的卓越周期。

1.5　超参数优化

XGBoost模型的性能依赖于超参数的设置。因此，针对7个关键超参数进行系统优化，包括：决策树数量n_{\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}}、树的深度d_{\mathrm{m}\mathrm{a}\mathrm{x}}、叶节点最小权重w_{\mathrm{m}\mathrm{c}}、损失函数下降值\mathrm{\gamma }、特征随机采样比例R_{\mathrm{b}\mathrm{t}}、训练样本采样率S_{\mathrm{s}\mathrm{u}\mathrm{b}}以及学习率L_{\mathrm{r}}。参数优化过程采用Python开源库scikit-learn中的网格搜索方法实现。

在优化超参数的过程中，效率的提升需要综合评估超参数的属性和各自的取值域来实现。综合分析超参数的特性及其可能的取值范围，并学习相关文献中超参数的取值经验^［18］。首先，划定1个较大的优化范围，并以相对较大的步长进行初步的粗略调整；然后，以初步调整得到的粗调值为基础，逐步缩小优化范围，并采用较小的步长进行精细微调。具体的调参取值见表2。通过这种方法，可达到最佳的模型性能。

2 数据处理

2.1　数据来源

机器学习是1种数据驱动的智能算法，通过自动分析海量数据识别其中潜在的模式与规律，并基于学习到的知识实现对新数据的预测。因此，选取的数据集是否具有代表意义是决定机器学习模型是否普遍适用的重要因素。为了使模型能够适用于不同大小地震及不同距离处的台站，避免产生漏捡拾及误捡拾风险，选取日本K-NET强震动观测台网1996年—2017年间197次地震的9 978条竖向强震动观测记录作为模型的训练集和验证集，震级范围为2.0~9.0级，震中距范围为0~655 km。此外，选取日本KiK-net强震动台网1996年—2019年间299次地震的6 726条竖向强震动观测记录作为测试集，震级范围为1.7~6.1级，震中距范围为0~502 km。3种数据集共计496次地震记录、16 704条加速度记录。所选记录的台站分布与震源分布如图4所示。3种数据集的震级和震中距分布如图5所示。

2.2　P波初至人工提取

采用多阶段人工识别方法确定P波初至。首先，根据竖向加速度记录时程曲线的幅值及周期变化特征；然后，将竖向加速度记录按震中距和绝对记录时间绘制震中距-时间剖面图，并利用地壳速度结构（Crust 1.0）计算不同距离处的P波理论走时进行对比；最后，检查P波初至手动识别的正确性，对无法准确识别P波初至的地震记录予以剔除。

2.3　测试集有效性

信噪比是反映背景噪声大小和地震记录中数据质量的重要标准，也是衡量P波拾取算法能力的重要判别标准。每条记录的信噪比S_{\mathrm{P}\mathrm{N}}为

S_{\mathrm{P}\mathrm{N}}=\mathrm{l}\mathrm{g}\left(S_{\mathrm{P}}/S_{\mathrm{N}}\right)（8）

其中，

S_{\mathrm{P}}=\mathrm{v}\mathrm{a}\mathrm{r}\left(a_{\mathrm{P}}\right)

S_{\mathrm{N}}=\mathrm{v}\mathrm{a}\mathrm{r}\left(a_{\mathrm{n}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{e}}\right)

式中：\mathrm{v}\mathrm{a}\mathrm{r}（·）为方差计算函数；a_{\mathrm{P}}为P波到达后3 s内竖向加速度幅值的集合；a_{\mathrm{n}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{e}}为由竖向背景噪声加速度幅值构成的集合，其时间窗取自记录开始后2 s至P波到达前；S_{\mathrm{P}}为P波信号方差；S_{\mathrm{N}}为背景噪声方差。

3种数据集的信噪比分布情况如图6所示。由图6可知，3种数据集的信噪比分布接近，表明测试集数据可用于模型性能的客观评估。

3 试验结果与分析

3.1　模型评估

在模型评估过程中采用测试样本是否被正确或错误识别的百分比为评价指标。使用真正率T_{\mathrm{P}\mathrm{R}}、真负率T_{\mathrm{N}\mathrm{R}}和假正率F_{\mathrm{P}\mathrm{R}}衡量模型的预测效果。以地震信号为正样本，背景噪声为负样本，则T_{\mathrm{P}\mathrm{R}}为模型将所有地震信号样本正确识别的比例；T_{\mathrm{N}\mathrm{R}}为模型将背景噪声样本正确识别的比例；F_{\mathrm{P}\mathrm{R}}为模型将背景噪声样本错误识别为地震信号的比例。其定义式为

T_{\mathrm{P}\mathrm{R}}=\frac{T_{\mathrm{P}}}{T_{\mathrm{P}}+F_{\mathrm{N}}}\times \mathrm{100}\mathrm{\%}（9）

T_{\mathrm{N}\mathrm{R}}=\frac{T_{\mathrm{N}}}{T_{\mathrm{N}}+F_{\mathrm{P}}}\times \mathrm{100}\mathrm{\%}（10）

F_{\mathrm{P}\mathrm{R}}=\frac{F_{\mathrm{P}}}{F_{\mathrm{P}}+T_{\mathrm{N}}}\times \mathrm{100}\mathrm{\%}（11）

式中：T_{\mathrm{P}}为地震信号样本被模型正确识别的个数；T_{\mathrm{N}}为背景噪声样本被模型正确识别的个数；F_{\mathrm{P}}为背景噪声样本被模型错误识别为地震信号的个数；F_{\mathrm{N}}为地震信号样本被模型错误识别为背景噪声的个数。

模型评价指标见表3。由表3可知：所提模型的训练集真正率为92.34%，验证集真正率为93.13%，测试集真正率为91.27%，表明模型在3个数据集上的真正率接近，没有出现过拟合。该模型在3个数据集上的AUC-ROC曲线如图7所示。由图7可知，模型在3个数据集上的AUC数值均较高且接近，说明模型具有很好的分类性能。

为评估XGBoost模型在实时系统中的适用性，计算了XGBoost模型的运行耗时和内存占用情况。选用测试设备配置为Intel Xeon W-2223 CPU @ 3.60 GHz处理器和内存为16.0 GB的台式工作站。在不考虑模型加载、特征参数提取等处理过程的条件下，XGBoost模型预测1个样本所用的平均时间为1.496 ms，CPU占用率为65.60%。通常，地震观测设备的采样间隔为5~10 ms，可见所提模型具备在实时系统中使用的条件。

3.2　对比试验

为进一步展示P波识别精度，将所提方法与EDP-Picker及STA/LTA方法进行对比。其中，STA/LTA方法的特征函数C(k)和长短时窗特征函数比R_{\mathrm{S}\mathrm{L}\mathrm{T}\mathrm{A}}（j）为

C(k)=a_k^{\mathrm{2}}+{\left[a_k-a_{k-\mathrm{1}}\right]}^{\mathrm{2}}（12）

R_{\mathrm{S}\mathrm{L}\mathrm{T}\mathrm{A}}(j)=\alpha \frac{k_{\mathrm{2}}{\displaystyle \sum _{i=j-k_{\mathrm{1}}+\mathrm{1}}^j}C(i)}{k_{\mathrm{1}}{\displaystyle \sum _{i=j-k_{\mathrm{2}}+\mathrm{1}}^j}C(i)}（13）

式中：\alpha为平滑系数，通常取值为0.99；k_{\mathrm{1}}和k_{\mathrm{2}}为当前时刻j前某时刻点且。

基于日本K-NET和KiK-net台网数据特性（P波到达前噪声段长度≤15 s），通过测试集优化确定短时窗长度为0.4 s，长时窗长度为10.0 s，触发阈值为10。

为了便于在同等条件下对比3种方法的性能，选取测试集中人工识别P波到时不小于10 s的记录进行对比，不同模型误差分布对比如图8所示。由图8可知，当P波初至拾取误差在±0.5 s内时，XGBoost模型的样本占比为93.3%，EDP-Picker为91.9%，STA/LTA为83.6%，表明XGBoost模型具有更高的拾取精度。

测试集P波初至拾取误差随信噪比的分布情况如图9所示。为定量评估模型的拾取性能，以±0.5 s为阈值划分误差区域：蓝色背景标记提前触发区（误差<-0.5 s），橙色背景标记滞后触发区（误差>0.5 s）。由图9可知：XGBoost模型的误触发样本中，提前触发占比为4.27%，滞后触发占比为5.26%，总误触发率为9.53%；EDP-Picker方法的提前触发比例达5.03%，滞后触发比例为6.50%，总误触发率为11.53%，表明XGBoost方法在抗干扰和防止滞后触发方面均具有更好的性能；对于S_{\mathrm{P}\mathrm{N}}>4的记录，模型误触发的比例很低，而对于3的记录，模型的识别成功率显著降低。

为进一步分析模型拾取误差较大的样本与震级、震中距的关系，统计了测试集识别误差大于1 s的所有样本随震级、震中距分布情况，结果如图10所示。由图10可知：小震的误触发率较高而大震误触发率较低，6级以上地震的误触发主要出现在震中距200 km以外。

3.3　特征重要性

特征重要性反映了特征参数在P波初至拾取中的贡献，使用开源库scikit-learn计算所得各特征参数的重要性排序如图11所示。特征参数重要程度依次为：卓越周期、三特征之和、加速度能量梯度和卓越周期梯度，表明卓越周期是P波初至自动识别中最关键的特征参数。

为进一步分析不同特征参数值与拾取误差之间的关系，统计了4个特征参数的取值和误差分布，结果如图12所示。由图12可知，误触发主要集中在特征参数值较小的条件下，加速度能量梯度在0.8以上、卓越周期在0.5 s以上、卓越周期梯度在0.5以上以及三特征之和在1以上时误触发较少。

4 结论

（1）模型在3个数据集上的真正率、真负率、假正率等评价指标接近且均超过90%，表明模型具有高分类精度和强泛化能力的特点，没有出现过拟合；AUC值接近0.98，表明模型具备高分类精度与鲁棒性；采用XGBoost进行初判并结合AIC函数修正的双阶段拾取框架，可解决单一模型滞后触发问题。

（2）与现有的EDP-Picker方法及STA/LTA方法相比，XGBoost方法在±0.5 s误差范围内的拾取样本占比为93.3%，优于EDP-Picker（91.9%）和STA/LTA（83.6%）；所提模型提前触发样本占比为4.27%，优于EDP-Picker（5.03%）和STA/LTA（5.39%），表明该方法具有更高的识别精度和更强的抗干扰能力。

（3）单样本预测平均耗时1.496 ms（通常强震动观测设备的采样间隔在5~10 ms之间），CPU占用率为65.6%，表明所提方法满足实时地震预警系统的时效性需求。

参考文献

[1] 卢建旗，李山有，谢志南，等.实时演化地震定位方法及其不确定性分析［J］.自然灾害学报，2021，30（3）：52-62.

[2] LU Jianqi， LI Shanyou， XIE Zhinan， et al. Real-Time Evolutionary Earthquake Location Method and Its Uncertainty Analysis ［J］. Journal of Natural Disasters， 2021， 30 （3）： 52-62. in Chinese

[3] 王自法，廖吉安，王延伟，等.基于深层卷积神经网络的震级快速估算方法［J］.地球物理学报，2023，66（1）：272-288.

[4] WANG Zifa， LIAO Ji’an， WANG Yanwei， et al. A Fast Magnitude Estimation Method Based on Deep Convolutional Neural Networks ［J］. Chinese Journal of Geophysics， 2023， 66 （1）： 272-288. in Chinese

[5] 杨长卫，张凯文，吴东升，等.基于地震P波初至2s功率信息的高速铁路地震预警系统快速震级估算方法［J］.铁道学报，2024，46（7）：167-173.

[6] YANG Changwei， ZHANG Kaiwen， WU Dongsheng， et al. Fast Magnitude Estimation Method for High-Speed Railway Earthquake Early Warning System Based on P-Wave 2 s Power Information ［J］. Journal of the China Railway Society， 2024， 46 （7）： 167-173. in Chinese

[7] 王子珺，赵伯明.高速铁路地震预警震级快速准确预测方法［J］.中国铁道科学，2017，38（2）：127-132.

[8] WANG Zijun， ZHAO Boming. Rapid and Accurate Prediction Method for Magnitude in Earthquake Early Warning for High Speed Railway ［J］. China Railway Science， 2017， 38 （2）： 127-132. in Chinese

[9] 宋晋东，朱景宝，刘艳琼，等.基于支持向量机预测模型的高速铁路现地地震预警方法［J］.中国铁道科学，2022，43（5）：177-187.

[10] SONG Jindong， ZHU Jingbao， LIU Yanqiong， et al. On-Site Earthquake Early Warning Method for High-Speed Railway Based on Support Vector Machine Prediction Models ［J］. China Railway Science， 2022， 43 （5）： 177-187. in Chinese

[11] ROSS Z E， MEIER M， HAUKSSON E. P Wave Arrival Picking and First-Motion Polarity Determination with Deep Learning ［J］. Journal of Geophysical Research： Solid Earth， 2018， 123 （6）： 5120-5129.

[12] ANANT K S， DOWLA F U. Wavelet Transform Methods for Phase Identification in Three-Component Seismograms ［J］. Bulletin of the Seismological Society of America， 1997， 87 （6）： 1598-1612.

[13] HILDYARD M W， NIPPRESS S E， RIETBROCK A. Event Detection and Phase Picking Using a Time-Domain Estimate of Predominate Period Tpd ［J］. Bulletin of the Seismological Society of America， 2008， 98 （6）： 3025-3032.

[14] ROSS Z E， YUE Y， MEIER M A， et al. PhaseLink： a Deep Learning Approach to Seismic Phase Association ［J］. Journal of Geophysical Research： Solid Earth， 2019， 124 （1）： 856-869.

[15] ZHU W， BEROZA G C. PhaseNet： a Deep-Neural-Network-Based Seismic Arrival-Time Picking Method ［J］. Geophysical Journal International， 2019， 216 （1）： 261-273.

[16] 李健，王晓明，张英海，等.基于深度卷积神经网络的地震震相拾取方法研究［J］.地球物理学报，2020，63（4）：1591-1606.

[17] LI Jian， WANG Xiaoming， ZHANG Yinghai， et al. Research on Seismic Phase Picking Method Based on Deep Convolutional Neural Network ［J］. Chinese Journal of Geophysics， 2020， 63 （4）： 1591-1606. in Chinese

[18] ALLEN R. Automatic Phase Pickers： Their Present Use and Future Prospects ［J］. Bulletin of the Seismological Society of America， 1982， 72 （6B）： S225-S242.

[19] 马强，金星，李山有，等.用于地震预警的P波震相到时自动拾取［J］.地球物理学报，2013，56（7）：2313-2321.

[20] MA Qiang， JIN Xing， LI Shanyou， et al. Automatic P-Arrival Detection for Earthquake Early Warning ［J］. Chinese Journal of Geophysics， 2013， 56 （7）： 2313-2321. in Chinese

[21] SARAGIOTIS C D， HADJILEONTIADIS L J， PANAS S M. PAI-S/K： a Robust Automatic Seismic P Phase Arrival Identification Scheme ［J］. IEEE Transactions on Geoscience and Remote Sensing， 2002， 40 （6）： 1395-1404.

[22] GALIANA-MERINO J J， ROSA-HERRANZ J L， PAROLAI S. Seismic P Phase Picking Using a Kurtosis-Based Criterion in the Stationary Wavelet Domain ［J］. IEEE Transactions on Geoscience and Remote Sensing， 2008， 46 （11）： 3815-3826.

[23] LU J Q， LI S Y， XIE Z N， et al. Energy-and Predominant-Period-Dependent P-Wave Onset Picker （EDP-Picker） ［J］. Seismological Research Letters， 2020， 91 （4）： 2355-2367.

[24] CHEN T Q， GUESTRIN C. XGBoost： a Scalable Tree Boosting System ［C］// Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. San Francisco California， USA： ACM， 2016： 785-794.

[25] 王浩钰，付国永，李黎，等.一种基于Xgboost+GWR算法的中国地区降水量空间插值模型［J］.水电能源科学，2021，39（12）：7-11.

[26] WANG Haoyu， FU Guoyong， LI Li， et al. A Spatial Precipitation Distribution Prediction Model for Chinese Regions Based on Xgboost + GWR Algorithm ［J］. Water Resources and Power， 2021， 39 （12）： 7-11. in Chinese

基金资助

中国地震局工程力学研究所基本科研业务费专项资金资助项目(2024C27)

应急管理部重点科技计划项目(2024EMST040409)

国家重点研发计划项目(2024YFC3012802)