基于分形维数的高速铁路地震预警系统地震P波震相识别方法

杨长卫，张凯文，吴东升，张志方，张 良，瞿立明

(1.西南交通大学 土木工程学院，四川 成都 610031；2.中国铁路武汉局集团有限公司，湖北 武汉 430071；3.西南交通大学 地球科学与环境工程学院，四川 成都 611756)

截至2022年底,中国高铁运营里程突破4.2万km,占据世界高铁运营里程2/3以上[1]。高速铁路的运行很大程度地缓解了我国交通运营压力,为我国经济发展做出了巨大贡献。地震作为一种随机性强、破环性大的自然灾害,对高速铁路的运营造成了极大的安全隐患[2]。运行中的高速铁路列车在地震作用下可能有脱轨、碰撞等危害性极强的次生灾害发生,造成严重的人员伤亡及经济损失,许多研究学者研究了地震作用下的铁路桥梁、轨道的运行安全[3-6]。我国在2017年建成了高铁地震预警系统[7],该系统主要利用P波和S波、S波和电磁波的速度差来实现的[8-9]。地震P波自动识别作为高速铁路地震预警系统的首要工作,其识别地震P波的精度与速度直接影响了后续地震预警系统工作震中定位[10]及震级估算[11-12]的准确性。研究适应我国高速铁路地震预警的快速、精准地震P波震相识别算法十分必要,能够有效保障我国高速铁路的安全运营。

数十年来,大量研究学者研究了提高地震P波自动识别的速度与精度的方法,主要形成了基于地震波长时窗均值与短时窗均值之比(short time average/ long time average,STA/LTA)[13-15]与基于赤池信息法则(Akaike information criterion,AIC)[16-17]的地震波识别方法。此外,也有国内外研究学者提出了不同的地震P波震相初至识别方法,包括相邻道互相关法[18]、分形维数方法[19]、小波变换方法[20-21]、数字图像分割法[22]、CM及MCM算法[23]、人工智能算法[24-25]等。上述方法各有特点,但因存在局限性没有作为地震P波的震相识别主流算法。已有研究结果显示,STA/LTA方法在地震波振幅较大时的工作效果良好,AIC方法的识别的精准度较高。STA/LTA方法由于工作原理方法,在地震P波振幅较小时的识别精度较差;而AIC方法通过建立自回归模型,识别时间窗口长度内的函数最小值作为地震波的精准到达时刻,所需要的时间较长。据统计,仅在四川境内,2009年1月1日至2022年2月12日共发生3级以上地震 1 084次,其中5级以下的地震1 040次,占比95.94%。在小震频发的背景下,为确保高速铁路地震预警系统的良好工作及高速铁路的安全运营,研究快速精准的小震P波自动识别算法十分必要。分形理论应用在地震P波震相识别最早由Boshetti在1996提出的,后续不断有学者对分形理论应用在地震波识别的方法进行研究[26-28],提出STAFD/LTAFD、灰度边缘识别等方法,并得出微震震相识别分形理论具有更高精确度的结论。该理论由于分形维数计算较慢的缺点,许多研究学者将计算一次分形维数的时间间隔设置为5个采样点及以上以提高识别速度,识别精度仍有进一步提高的空间。

本文在已有研究学者的基础上,对分形维数应用在地震P波震相识别进行了进一步的研究,优化了分形维数的计算方法,提高了分形维数的计算速度,提升了识别精度。此外,本文在时效性和准确性上与STA/LTA方法及AIC方法进行对比。结果表明,在准确性方面,本文提出方法的平均识别误差达到了0.006 3 s,优于STA/LTA方法且与AIC方法基本持平(<0.01 s)。在时效性方面,本文提出方法的平均识别时长为0.16 s,优于AIC方法。这些结果满足我国高速铁路地震预警系统的要求,以期为我国高速铁路地震预警系统提供参考并作为补充算法。

1 基于分形斜率的地震P波识别算法

1.1 分形理论

分形理论最早由美国数学家Mandelbrot在1975年提出。自此,分形理论被广泛应用解决不同的实际问题。分形理论改变了传统的维数观念,如点是一维的,面是二维的,体是三维的。在分形理论中,曲线的维数被认为是1～2维的,而不规则平面的维数被认为是2～3维的,即引入线、面的“粗糙度”的概念。分形理论在地震波震相识别的核心理念是,通过识别地震波到达前后曲线粗糙度变化来识别地震波的精准到达时间。曲线分形维数的计算方法通常包括尺码法和网格法,已有研究成果表明尺码法更适合应用在地震P波自动震相识别[29],本文同样采用尺码法进行分形维数计算。尺码法计算流程见图1(a),其核心计算原理为对一段曲线,移动不同的尺码ri覆盖整条曲线,共移动了Ni次,则可得到曲线的近似长度Li,对不同的尺码ri和近似长度Li取对数进行一次拟合,见图1(c),得到斜率F,则分形维数D计算式为

D=1-F

(1)

式中:D为分形维数;F为拟合一次斜率。

在计算分形维数时,应该保证在移动尺码时,尺码的两端都在曲线上,见图1(b)。同时,为得到曲线的近似长度Li,需要多次改变尺码ri的长度,通常需要改变十余次尺码长度才可以完成一次分形维数计算,这一行为严重降低了分形维数的计算速度。为提高分形维数的计算速度,改变尺码ri长度的次数被降低至5次,减少了分形维数计算工作量,同时并没有影响分形维数曲线的性质。结果表明,在降低尺码变化次数后,一次分形维数的计算速度由0.27 s提高至0.05 s。

1.2 计算原理

地震预警系统对自动震相识别算法的最大要求在于快速及时准确。很多研究学者为实现分形理论在地震预警系统实现的可能,选择间隔数个采样点进行一次分形维数计算,然而这一做法损失了识别的精度。基于改进的分形维数算法,本文提出了连续采样点计算分形维数的方法。连续地震波分形维数具体计算方法是选择一定长度的时间窗口覆盖地震波,计算一次分形维数,并将分形维数标记在时间窗口的右侧时刻。随着时间窗口的逐步移动,可以在地震波的每个时刻标记上对应分形维数,即曲线的粗糙度,进而形成连续的分形维数曲线。该方法在不损失计算速度的情况下,得到连续分形维数曲线,充分反映了地震波的细节特征,从而得到更加精准的地震P波到达时间。分形维数连续与间断计算曲线见图2。图2中,显示了间隔5个采样点计算一次分形维数得到的分形维数曲线与连续分形维数曲线的区别。

图2 分形维数连续与间断计算曲线

计算分形维数的时间窗口沿着地震波曲线移动,地震波可以被划分为3个阶段,见图3(a),同样地,在分形维数曲线中3个阶段也有明显的特征,见图3(b)。第一阶段为地震波到达之前,时间窗口内的地震波信号仅有震前白噪声组成,因此分形维数曲线十分平坦,记该段为CFD1;第二阶段为地震波初至,时间窗口内的地震波信号由部分震前白噪声及地震波与白噪声组成,曲线的粗糙度快速上升,因此分形维数曲线也快速上升,记该段为CFD2;第三阶段为地震波完全到达,时间窗口内的地震波信号仅有地震波与白噪声组成,分形维数曲线反映了地震波变化的细节特征,开始无规则波动,记该段为CFD3。其中第一阶段与第二阶段的临界点记为点1,第二阶段与第三阶段的临界点记为点2。

图3 分形斜率识别原理

本文进一步引入了分形斜率Ki去分析分形维数曲线的变化趋势,其中斜率的计算式为

(2)

式中:Di为i时刻的分形维数;t为时间窗口长度。

分形斜率即一个时间窗口的左右两点的分形维数形成的斜率,反应了分形维数曲线随着时间窗口逐步移动的变化特征,进而形成分形斜率曲线,进一步反映了地震波的波动特点,分形斜率曲线可以划分为四个阶段,见图3(c)。第一阶段,分形斜率十分平坦,此阶段时间窗口的左点与右点均在CFD1上,两点的分形维数基本持平,从而分形斜率也基本持平;第二阶段,分形斜率快速上升,此阶段时间窗口的左点在CFD1上、右点在CFD2上,随着时间窗口的移动,左点保持不变,右点快速上升,分形斜率也快速上升并上升;第三阶段,分形斜率快速下降,此阶段时间窗口的左点在CFD2上、右点在CFD3上,随着时间窗口的移动,左点快速上升,而右点进入波动段,从而分形斜率快速下降;第四阶段,分形斜率不断波动,此阶段时间窗口的左点和右点均在CFD3上,反映了分形维数与地震波信号的变化特征。同样地,临界点1、临界点2在分形斜率曲线中也是第一阶段与第二阶段、第二阶段与第三阶段的临界点2。在分形维数斜率第二阶段的末尾时刻,时间窗口的左点刚好在临界点1上,右点刚好在临界点2上。若时间窗口继续推进,由于时间窗口左点落入到CFD2上,分形维数开始上升,则分形斜率开始下降,将该临界点2时刻定义为极值时刻,而分形斜率也有极值的表现。本文通过识别分形斜率极值判断地震波的精准到达时刻。若分形斜率极值时刻为i,则地震波的精准到达时刻为i-t。

本文通过设置极值时刻判断标准和引入标准差去提高识别精度。首先,分形斜率极值时刻应该大于相邻4个采样的分形斜率。其次,计算极值点时刻前15个采样点至215个采样点的分形斜率标准差,即CFD1阶段的标准差,评估该阶段的分形斜率曲线的波动程度。最后,要求极值时刻的分形斜率应该大于分形斜率第一阶段最大值数倍及以上。标准差的计算式为

(3)

根据震前白噪声的分形斜率波动情况,可以将其分为3个波动程度水平,分别为:标准差s小于0.002、大于0.02和0.002～0.02之间。为进一步减少识别误差,当标准差s小于0.002时,分形斜率的极值时刻应该大于分形斜率第一阶段最大值的6倍及以上;当标准差s在0.002～0.02之间时,分形斜率的极值时刻应该大于分形斜率第一阶段最大值的4倍及以上;当标准差s大于0.02时,分形斜率的极值时刻应该大于分形斜率第一阶段最大值的2倍及以上。本文提出的方法的具体识别流程见图4。

图4 识别流程

1.3 影响因素分析

时间窗长度的选择直接影响了分形维数地震P波震相识别工作的速度与精度。时间窗长度过长导致分形维数计算速度降低,时间窗长度过短则导致地震波加速度曲线波动过于敏感容易误触发。为保证高速铁路地震预警系统工作的时效性与准确度,有必要研究确定合适的时间窗长度。

本文选取了100组日本K-Net &Kik-net地震台网的地震数据进行分析,其地震事件的震级均小于M 5,震中距小于100 km,采样频率为100 Hz,地震波方向均为UD方向。在进行对比时,认为人工检测地震P波到达时间是精准的。

本文提出的地震P波自动识别算法的精准到达时刻为极值时刻减去时间窗长度,为了满足高速铁路地震预警系统的时效性,时间窗口长度的选择不亦过大,保证得到地震波精准到达时刻的速度满足Q/CR 634—2018《高速铁路地震预警监测系统技术条件》[30]的技术要求。总共4个时间窗口长度被选择进行分形维数、分形斜率计算并进行识别,包括5、10、15、20个采样点,即0.05、0.10、0.15、0.20 s。按照本文提出的识别方法,可以得到4种不同时间窗口自动识别得到的地震波震相到达时刻,其识别结果与人工检测的误差统计见图5。结果显示,不同时间窗口长度下的识别误差基本保持在±0.1 s以内,但也存在一定差异。不同时间窗口长度的识别平均误差均在±0.05 s以内,以不同时间窗口长度识别标准差作为标准能够更好的评价识别效果。显然,15个采样点的识别结果标准差更加优异。

图5 不同时间窗口长度识别误差直方图

2 案例分析

STA/LTA与AIC方法作为地震预警系统的主流地震P波震相识别方法,将本文提出的方法与这两种方法进行比较。其中,STA/LTA的计算参数按照马强提出的进行设置,短时间窗口的长度为50个采样点,长时间窗口的长度为3 000个采样点,阈值为10[31]。AIC的时间窗范围为STA/LTA方法识别得到的到达时刻前300个采样点及到达时刻后30个采样点,计算共计330个采样点的AIC函数,识别最小值作为地震波的精准到达时刻,其中,k时刻AIC函数fALC(k)的计算式为

fAIC(k)=k×lg{var[x(1,k)]}+(N-k-1)×lg{var[x(k+1,N)]}

(4)

式中:N为选定窗口中所有数据的数量;var[x(1,k)]、var[x(k+1,N)]均为两个窗口数据采样点的方差;x(i)为i时刻的地震波信号加速度值,i=1,2,…,N。

本文共选取368组地震波进行了鲁棒性测试,并与STA/LTA及AIC方法进行对比,其选择标准与之前相同。同样,人工检测的地震波到达时刻被认为是准确的。本文选择的对比标准是平均误差及标准差。三种不同方法的对比结果见表1,统计直方图如图6所示。根据结果显示,AIC算法的平均识别误差及标准差最低达到了0.003 4 s及0.025 6 s。STA/LTA算法的平均识别误差为0.090 3 s,标准差为0.063 0 s。而本文提出的分形斜率算法的平均识别误差为0.006 3 s、标准差为0.043 8 s,基本与AIC方法持平,优于STA/LTA算法。

表1 不同识别算法识别结果 s

图6 不同识别方法的识别误差统计

本文提出的地震P波识别算法,是通过识别地震波分形斜率极值时刻来判断地震波的精准到达时刻,即在地震波到达一个时间窗口后加上分形维数运算时间可以得到地震波的精准到达时刻。本文选取的时间窗口长度为0.15 s,而优化分形维数计算后的一次计算时间为0.05 s,平均误差为0.006 3 s,即平均在地震波到达0.16 s内可以得到地震波的精准到达时间。STA/LTA方法是识别均值触发阈值的时刻作为地震波的精准到达时刻,即平均在地震波到达后0.090 3 s可以得到地震波的精准到达时间。而AIC算法需要在利用STA/LTA方法识别震相时刻后,后推0.3 s后进行AIC函数计算,即地震波到达时刻后0.39 s后获得地震波的到达时刻,同样滞后于本文提出的分形斜率算法。本文提出的算法,在仅比STA/LTA算法延迟0.07 s的情况下,提高了0.1 s的识别精度,在与AIC算法相比提前了0.23 s获得了与AIC算法基本相同的精准度,为后续高速铁路地震预警系统的震中定位及震级估算工作的进行提供了参考,并且在识别速度上达到《高速铁路地震预警监测系统技术条件》[30]的技术标准。

3 结论

1)提出了一种基于分形维数的小震地震P波震相识别方法。与现有分形理论地震P波自动震相识别算法比较,本文改进了分形维数的计算方法,在不改变分形维数曲线的特征的情况下,提高了分形维数的计算速度。

2)在改进的分形维数计算方法的基础上,改变了以往间隔计算分形维数的算法,采用连续滑移时间窗口进行分形维数计算,得到连续分形维数曲线充分反映地震波的细节特征。此外,分形斜率被引入,并划分为4个阶段,分析临界点的分形斜率特征,通过引入标准差和判断标准,识别分形斜率极值点识别判断地震P波的精准到达时刻。

3)此外,研究了在震级小于5的地震事件中,分形维数时间窗口长度的选择的最优值,为0.15 s。本文提出的算法与STA/LTA、AIC算法相比,识别精度的平均误差达到了0.006 3 s,标准差达到了0.043 8 s,识别误差优于STA/LTA算法。在识别速度上,本文提出的算法平均可以在地震P波到达0.16 s后得到精准到达时刻,优于AIC算法,且可以满足高速铁路地震预警系统时效性与准确度的要求。