基于多属性的P波到达时间自动拾取方法*

辛 平， 李志华

(江南大学 物联网工程学院 计算机科学系，江苏 无锡 214122)

0 引 言

在地震震相检测中，首先要做的是精确拾取震相到时，即精确拾取P波到时。地震波的准确拾取在许多与地震相关的研究中起着关键作用。目前，自动拾取方法已占据主导地位[1]。文献[2]用简单的振幅和能量阈值来确定P波到时;由于该算法只用了振幅和能量是否大于阈值判定P波到时，所以其拾取的精确度不高；文献[3]采用长短时平均(short term average/long term average,STA/LTA)方法确定地震震相到时,可以快速拾取震相到时，但其拾取震相到时的精确度受噪声影响比较大，高噪音时出现误拾甚至拾取不到P波。此外经验阈值的选取直接影响STA/LTA拾取的准确度；文献[4]采用最大似然(maximum likelihood,ML)估计方法来确定地震P波到时，虽然能够比较精确确定P波到时，但ML方法认为噪声干扰是一成不变的，忽略了突发人为或者其他非地震因素所引起的噪声，这显然是不合理的；文献[5]通过计算P波偏振度来确定P波到时,该方法提高了识别可靠度，但计算量大；文献[6]将STA/LTA、偏振分析、AR分析以及瞬时频率分析这4种方法结合起来确定P波到时，虽然精确度上有了极大的提高，但算法的计算量很大；文献[7]采用自回归方法来拾取P波到时，能够比较精确确定P波到时；文献[8]提出了一种基于局部极大值分布的P波到时自动识别方法，算法将能量和频率特征结合起来，能够在高噪音时比较精确拾取P波。

本文提取地震信号的15维特征作为数据样本，并通过特征选择来对数据样本进行降维处理，最后利用支持向量机(support vector machine,SVM)进行分类预测，确定P波到时。

1 基于改进的局部线性嵌入和核分类方法的P波拾取

1.1 多维特征的构造

为了获取丰富的地震信号特征，本文通过对实测地震信号进行变换和处理，从时域和频域两方面综合选取15维特征利于数据深度挖掘。15维特征包括震相、振幅、偏振度、分形维数、STA/LTA，ML、振幅比、能量比、曲线长度比、峭度、偏斜度、均方根振幅、平均能量、振幅峰态，频域参数中选取频率。具体如下:

1)震相到时由STA/LTA确定。

2)振幅从信号中直接读取。

3)偏振度F采用文献[5]中的方法计算，F=m1/m2,m1和m2分别为由地震信号数据形成的3×3阶矩阵的最大特征值和次最大特征值。其矩阵构造为

(1)

式中x，y，z为地震信号3分量数据，cov(x,y)为x和y的协方差

(2)

式中μx和μy为x和y的平均值。

4)分形维数D的计算为D=1-S,S为不同步长和其对应的曲线长度拟合的直线斜率。

5)i时刻振幅比为

(3)

式中x为地震数据，m为地震数据时窗的起点，n为地震数据时窗的终点。

6)i时刻能量比为

(4)

7)i时刻曲线长度比为

(5)

式中 Δt为地震数据的采样间隔。

8)峭度为

(6)

式中s为时窗大小。

1.1.1 改进的STA/LTA方法

传统STA/LTA计算为

(7)

式中m为STA窗口大小,n为LTA窗口大小。

为提高震相拾取的准确度，对STA/LTA方法进行改进,即

(8)

1.1.2 改进的ML方法

由于震相初至前后地震信号的振幅有很大的变化，因此可以用ML估计来确定震相的初至时间。由于地震信号呈随机分布，所以本文假设地震信号变化前后的分布都服从正态分布，则ML计算为

(9)

式中σ1为震相初至前噪音方差，σ2为震相初至后方差。

(10)

改进的ML计算方法概括如下：

3)采用式(10)计算窗口内的似然估计值L。

1.2 基于改进的LLE特征抽取方法

LLE计算序列间距离采用欧氏距离，考虑每个区域上各个点分布是不均匀的，为了使每个区域上各点的分布整体呈均匀化，LLE中距离的计算改进为

(11)

式中M(i)，M(j)分别为与其点距离的均值。

改进的LLE算法具体描述如下:

1)计算每个样本点xi的k个近邻点xij(j=1,2,…,k)。

(12)

3)计算xi和xj的低维映射值yi和yj，使得ε(Y)最小

(13)

式中ε(Y)为损失函数，式(14)满足

(14)

M=(I-W)T(I-W)

(15)

1.3 面向时间序列的核分类方法

由于传统的SVM[10]距离度量采用欧氏距离，计算时间复杂度低(O(n2))，但其无量纲且只适用于等长的时间序列[11]。由于所研究的地震信号是时间序列，具有多维、离散的特性，且每一维都有其特定的量纲，为了使SVM更好地适用于地震时间序列数据，本文提出面向时间序列的SVM(time series SVM,TS-SVM)分类方法，在TS-SVM方法中，时间序列间相似性度量采用动态时间弯曲距离[12]。描述如下:

假设有时间序列X和Y，则其相似性度量测度为

(16)

(17)

式中 < >为空序列,Rest(X)={x2,x3,…，xn},Rest(Y)={y2,y3,…,yn}。

TS-SVM方法在核函数选取上采用高斯核函数，用DTW距离代替欧氏距离，核函数定义为

(18)

得到最优分类函数为

f(x)=sgn{(w*,x)+b*}

(19)

新提出的TS-SVM方法描述如下：

1)由计算时间序列X的动态弯曲距离DTW。

2)计算时间序列的内核函数，求解二次规划问题

得到最优解a*={a1,a2,…，ak}。

4)计算f(x)的值，实现对时间序列X的分类。

动态弯曲距离的时间复杂度为O(n2)，欧氏距离的时间复杂度也为O(n2)，TS-SVM并未改变传统SVM的时间复杂度。

1.4 基于TS-SVM的 P波拾取方法AMPAT

AMPAT流程如图1所示。

图1 基于TS-SVM的P波拾取方法流程

AMPAT方法概括如下：

1)数据预处理：计算15维特征构造特征空间集，考虑到特征空间的每一维特征数量级不同，需要进行归一化处理。对归一化处理的特征空间用LLE进行降维，选取最佳特征降维后的维数。选取合理数量的样本数据，50 %样本数据作为训练集，50 %作为测试集。

2)TS-SVM参数寻优:在训练样本集上进行有监督学习的训练，确定最佳参数。

3)构造P波拾取模型：对LLE降维后的特征数据，采用合理的TS-SVM寻优参数进行测试，拾取P波到时。

2 实验与结果分析

2.1 实验数据

本文数据由IRIS(Incorporated Research Institutions for Seismology)提供，数据格式为SAC。以美国东部海岸 2016年 3月份地震数据为样本，利用MATLAB进行仿真实验来验证本文提出的方法的准确性。样本的采样频率均为40Hz，记录了地震信号在某一时间段内的加速度变化情况。

2.2 结果分析

根据专家经验对样本信号进行人工拾取，其中P波到时如图2中直线。拾取得P波到时在[161.65,161.675]s之间，精确P波到时为161.652 8 s。

图2 人工拾取P波

对样本信号采用STA/LTA和ML方法处理,P波到时分别为161.618 7 s和161.631 5 s,如图3。

图3 2种方法处理结果

AMPAT方法得到P波到时为161.632 6 s。另取IRIS 2016年3月份美国中部海岸的10组地震数据，用这4种方法拾取P波到时如表1。可知，以人工拾取结果为参考，STA/LTA方法拾取误差比较大，最大为0.338 8 s，ML拾取结果比STA/LTA方法好，最大误差为0.061 4 s,AMPAT拾取结果比STA/LTA和ML结果要好，最大误差为0.055 4 s。进一步比较，以人工拾取方法结果为参照，其他3种方法计算平均误差为：STA/LTA，0.073 7 s；ML，0.0331 s;AMPAT，0.025 1 s。可知，AMPAT的平均误差要小于ML方法和STA/LTA方法的，ML拾取误差较好于STA/LTA方法。

表1 拾取10组地震数据的P波初至时间 s

综上，AMPAT在最大误差和平均误差都优于STA/LTA和ML方法，其对P波拾取的精度高于这2种方法。

3 结 论

针对地震数据单个特征受噪声等其他因素影响准确拾取P波到时难度之大这一问题，本文提取地震数据的15维特征并通过基于LLE降维和SVM拾取P波，实验对比证明，AMPAT方法获取的P波到时更加精确。由于AMPAT方法选取了15维特征，导致计算量较多，如何有效降低计算复杂度是下一阶段的研究目标。