基于多台站高频GPS的地震预警优化算法研究

张文浩，尹玲,2，胡文博

( 1. 上海工程技术大学 电子电气工程学院, 上海 201620；2. 中国地震局地质研究所, 北京 100029；3. 上海大学 通信与信息工程学院, 上海 200444 )

0 引 言

随着我国经济建设的发展，大型高层建筑、人口密集城市不断增多，我国作为地震活动最强烈的国家之一，地震造成的生命财产损失不计其数，如何减少地震造成的人员和财产损失，是目前迫切需要解决的问题. 由于地震的复杂性，当前的技术水平难以实现对地震的精准预测，而地震预警则是目前国际上公认的能够减轻地震灾害的有效手段之一[1]. 传统地震仪方法虽然比较成熟，能在地震发生后短时间内确定震级和震源，但对强震震级进行估算时会出现超出量程、积分结果带有误差等现象，导致震级饱和[2].2011年3月11日，日本东北部发生Mw9.0级地震，日本气象厅预警系统发布警报，但该预警系统的地震级别为Mw8.1，预测准确度有很大差距，同时确保正确预警并防止来自系统的误报也是不小的挑战.2016年8月1日，日本地震预警系统发出Mw9.1级地震警报，预估日本大面积地区将达到烈度7度，该误报造成了小田急电铁和首都圈大部分地铁紧急停运，一些重大工程紧急关停. 一旦发生误报，会引起民众恐慌，降低政府的信誉，造成不必要的损失.

地震预警系统(EEWs)基本可分为区域预警和单站预警[3]. 区域预警是基于密集的台站检测网络，收集并分析来自多个台站的数据估算相关源参数，用于远距离地点提供早期预警. 单站预警则是由震源区域附近的单个传感器组成，在某些系统中，为估计台站区域的震动强度，还需确定震级和震源[4]. 单站预警通常比区域预警更快，但容易发生误报，精确度较低[5-6]. 最合适的方式应该是既保留区域预警的准确度，又尽可能实现单站快速预警.

现有的EEWs多是基于实时强震网络框架建立的，系统中的强震仪以触发模式进行工作，当达到一定强度后记录数据并通过积分获得位移结果. 但强震仪有基线偏差、量程有限等问题，而现代预警技术将高频GPS技术引入地表位移监控中. 它既可观测静态位移也可记录由地震造成的动态位移，具有不限制量程，不积累误差，能直接获取位移数据等优点. 文献[7]阐述了高频GPS技术在地震学中的应用方法、实际意义及目前存在的问题. 文献[8]利用高频GPS数据使用最速下降法和OKADA模型进行反演，探讨了高频GPS的峰值地面加速度(PGA)、峰值地面速度(PGV)、峰值地面位移(PGD)计算地震烈度的可能性.

随着地震检测仪器在全球范围的大规模部署，数据规模增大，如何完成从大量的连续波形数据中实现到震拾取是一个棘手的问题. 近年来人们发明了各种震相识别算法，比较典型的有：基于短时窗平均/长时窗平均算法(STA/LTA)的方法[9]. 信号到达时，STA/LTA值会发生突变，当比值大于某个阈值后，则判定事件发生，该方法应用于信噪比(SNR)高的数据可达到较高准确率，SNR低的数据容易发生遗漏. 赤池信息量准则(AIC)方法[10]根据AIC寻求序列波形的极小点作为震相的到时点，该方法对震相和频率的变化较敏感，常被用于震相到时的精确拾取，与STA/LTA类似，AIC的结果也强烈依赖于SNR以及检测区间. 神经网络算法[11]可以利用神经元捕捉输入特征和输出值之间的非线性关系. 文献[12]利用神经网络对强震三分量数据进行到时提取，并加入噪声提高模型泛化能力，展现出神经网络在自动化震相拾取方面的优越性. 文献[13]针对地震事件和噪音分类这一问题训练卷积神经网络，对序列片段进行分类.

本文探索用深度学习的方法对高频GPS数据进行分析，将区域信息通过长短期记忆网络(LSTM)[14]融合到单站预警中，尽可能降低误报率同时保证快速反应. 分别使用融合区域特征的网络模型、未融合区域特征的网络模型、STA/LTA方法对新西兰地区多个台站无震高频GPS序列进行误报检测，比较每种方法的误报次数，另外为检测算法的有效性对新西兰凯库拉地震进行检出. 结果表明：融合区域特征后的模型不仅能够对地震事件进行检出，而且对比未融合区域特征的网络模型和STA/LTA能有效地减少误报次数，在应用于高频GPS时有很好地表现，可对当下的地震预警工作做出一定补充.

1 高频GPS数据获取与预处理

在选择台站时，尽可能选择离震中距离均匀、记录数据相对完整的台站，然后进行数据解算. 文章使用2016年11月13日凯库拉Mw7.8级地震[15]作为检出事件，由于临近该地震区域无震数据较少，本文选择周边区域记录大量无震数据的台站，利用6 h的无震数据去微调网络，达到模型迁移的效果.

1.1 实验数据

本文所使用的1 Hz高频GPS数据来源于新西兰地震信息网(Geonet)，该网络包括Strong motion sensor、检测区域形变的GPS sensor等台站，剔除掉部分观测数据质量较差和数据缺失严重的台站，选择2016年2月中数据较为完整的3个台站数据用于训练网络，处理时间范围为00:00:00—23:45:00 (UTC).选择临近地震区域3个台站数据用于迁移学习，处理时间为2016-11-13 T 00:00:00—06:00:00 (UTC)，为了验证本文所提出的方法在实际地震检出应用中的效果和可行性，本文对受凯库拉Mw7.8级地震影响的多个台站进行到震检出，该地震发生于2016-11-13 T 11:02:59 (UTC)的新西兰南岛中部地区，震中位置位于(173.064 7°E，42.724 5°S)，震源深度为15.0 km. 选取2016-11-13 T 06:00:00—09:00:00 (UTC)的多个台站无震时间序列作为测试数据. 表1列出了用于迁移学习的3个台站信息.

表 1 高频GPS台站信息表

1.2 数据处理方法

文章使用GAMIT track模块对新西兰南部地区1 Hz GPS数据进行解算，该模块采用差分定位模式，实际获取到的是目标站对于参考站的相对位置，所以参考站的选取会影响实际结果[15]. 对于震时数据的解算，如果基线距离过小，地震对参考站的影响会叠加在流动站位移序列上. 为尽可能避免同震干扰，选择距震中较远且观测质量较好的观测站作为参考站，并检查数据质量，准备精密星历文件，进行解算可得到所需台站的时间序列数据.

2 基于LSTM神经网络的地震预测模型

LSTM神经网络是一种特殊的循环神经网络(RNN)主要是解决长序列训练过程中的梯度消失和梯度爆炸问题[16]，在长序列中有更好地表现.

2.1 LSTM神经网络原理

LSTM与RNN相比不再是单一 t anh 层，而是多个门结构，分别是遗忘门、输入门、输出门，LSTM单元结构如图1所示.

图 1 LSTM单元结构图

2.1.1 遗忘门

遗忘门会使用当前时刻的输入xt和上一时刻的输出ht-1通 过σ函数计算上一时刻细胞状态Ct-1中信息保留的比重，0表示不保留，1表示保留. 其中Wf表示遗忘门的权重，bf为偏置项.

2.1.2 输入门

输入门可以利用当前时刻的输入xt和上一时刻的输出ht-1通 过σ函数决定更新哪些信息. 再利用xt和ht-1通 过 t anh 层得到临时细胞状态Wi表示输入门的权重，bi为偏置项，Wc表示细胞单元的权重，bc为偏置项.

将旧的细胞状态Ct-1更 新为新的细胞状态Ct.

2.1.3 输出门

输出门可以利用当前时刻的输入xt和上一时刻的输出ht-1通过σ函数判断输出细胞的哪些特征，并通过与 t anh 层的细胞状态Ct相乘得到单元输出ht.Wo表示输出门的权重，bo为偏置项.

2.2 网络结构

一个LSTM网络包含多个隐藏层，每个隐藏层由大量LSTM神经元组成，展开后等效于在不同层之间传递参数的前反馈神经网络. 每一层的输出作用于下一层的输入，所以它可以处理时间序列数据并提取特征，LSTM预测模型整体框架如图2所示. 为适应地震检出任务使用单层LSTM，使得网络能够学习高层次的时域特征，输入层为多个台站三分量数据，输出为单台站下一秒预测值.

图 2 LSTM时间序列预测框架

2.3 模型计算流程

在网络训练前需要对数据进行处理，主要包括划分数据集、将高频GPS数据转换成可以输入LSTM神经网络的格式，即样本、时间步、特征的三维数据.对输入的多台站数据进行归一化处理，将高频GPS时间序列转换化成监督学习数据，每个时间步前t时刻的数据作为训练样本x，当前t时刻的数据作为输出值y. 数据处理完毕后，由输入层输入处理好的多台站地震信号数据，隐藏层LSTM单元实现对地震信号进行时域分析，通过Dropout层避免过拟合，全连接层对每个时间步的特征进行全连接，再将每个时间步的输出结果经过激活函数后，由输出层及反归一化操作输出预测值，期间数据不断迭代，以降低损失函数值为目标直至收敛. 最后根据预测值动态生成阈值，通过真实值与阈值相比较，判断是否发生地震事件，算法流程如图3所示.

图 3 LSTM模型计算流程图

2.4 模型训练

训练预测模型的目的是让模型学习无震时的数据特征，并且用来预测未来时刻无震的数据值，通过预测值与观测值相比较来检出异常. 使用2016年2月中筛选的无震时间序列片段，最终选择使用YALD、VWMR、VBRM 3个台站数据构造数据集，选取85%的数据作为训练集，15%的数据作为测试集，隐藏层使用单层LSTM网络，并在每层之间加Dropout层，在训练时以抑制部分神经元的输出，避免在训练过程中过拟合. 最后一层是全连接层，使用线性激活函数，用于输出模型下一秒的预测值. 网络的优化算法选取为自适应运动估计算法(Adam)算法，该算法能使收敛过程更加平稳，可以为不同的参数自适应学习率. 使用均方根误差(RMSE)作为衡量模型性能的指标，度量真实值与预测值之间的相对偏差，其计算如式(8)所示. 其中yi(t) 是 在t时刻的实际观测值，(t) 是在t时刻模型输出的预测值，n是数据样本的总数.

训练所用硬件为AMD Ryzen 3900x 3.8 GHz CPU、16 GB内存、NVIDIA GeForce GTX2080ti 11 GB显存. 在选取batch size参数的过程中，batch size过小会导致随机性更大，难以达到收敛，但是能使模型学到一定噪声，提高模型的泛化能力. 过大的 batch size会增加迭代过程中的计算量，内存消耗大. batch size选择的影响如表2所示.

表 2 不同batch size对结果的影响

batch size大小对模型学习有很大影响，当batch size足够大时才可以对整个数据的梯度进行稳定的计算，但过大的batch size会在迭代的过程中增加计算消耗. 过小的batch size能学习到一定噪声，可以避免局部最优，但可能会由于样本不足，网络权重突然波动而导致收敛缓慢甚至不收敛，从而对计算造成负面影响. 由表2可以看出，在batch size小的情况下单次迭代时间长，batch size大的情况虽能加快计算，但是考虑到精度，选择batch size为64.

网络进行训练得到的Loss图像如图4所示，实线为训练集Loss值变化，虚线为验证集Loss值变化.训练Loss曲线从0.009开始减少，初始时收敛较快，经过约15次迭代后收敛速度降低，最终收敛于0.000 6.

图5为模型应用于测试数据集得出的YALD站预测图，计算得出RMSE为0.003 429，可以看到LSTM神经网络可以对台站时间序列进行有效预测，可以充分挖掘区域各个台站之间的相关性.

得到满足预测精度的模型后，用地震周边区域台站无震数据微调模型. 选择临近地震区域MRBL、WITH、WEST 3个台站6 h数据构建数据集用于迁移学习. 应用于测试集的预测结果如图6所示，实线为真实观测数据，虚线为来自预测模型的预测结果.由图可知预测模型可以准确预测到观测数据的变化趋势，最终测试集的RMSE为0.002 989.

图 4 LSTM网络模型Loss值变化曲线

图 5 YALD台站预测图

图 6 MRBL台站预测图

图 7 MRBL台站地震检出

图 8 WITH台站地震检出

图 9 WEST台站地震检出

3 模型应用

利用本文设计的地震检测模型可以对包含地震事件的GPS时间序列进行检出，并对多条无震长序列进行检测统计误报次数.

3.1 方法有效性和准确性

为验证本文模型可对异常事件进行检出，对包含凯库拉地震事件的时间序列进行检测，模型学到的特征为无震时的区域特征，而发生地震后序列会产生突变，突变若连续超出置信区间，则此次突变可作为异常实现检出. 在对多个受地震影响的台站进行检出，MRBL台站序列检出结果如图7所示，在第203时间步时观测值发生明显变化且超出置信区间，判定为异常实现地震事件，WITH台站序列检出结果如图8所示，在第141时间步时观测值发生异常，后续序列波动明显且多次超出置信区间，实现到震检出. WEST台站序列检出结果如图9所示，在第202时间步观测值异常变化超出置信区间，此时模型生成的阈值也跟随观测值发生一定突变，但模型结合区域之前的地表位移特征，使阈值变化不如实际值大，当实际值超出阈值时实现到震检出.

对无震时序进行误报检测时，正常情况下观测值不超出置信区间，若超出置信区间则判定为一次误报，文章对多个台站时间序列进行误报检测. 使用本文方法对MRBL台站10 800 s无震序列进行检测，检测所用时间为31.74 ms，最终未见误报；对WITH台站10 800 s无震序列进行检测，未见误报；对WEST台站10 800 s无震序列进行检测，出现1次误报. 另使用本文模型尝试对周边其他台站3 h无震序列进行检测，结果HANM、KAIK、LKTA台站均未发生误报，检测图如图10所示.

基于双轨迹构图可知本题有两解，解法1漏解.构造△ADB的外接圆圆心E，点E可在AB下方如图4；也可在AB上方，如图5；因此有必要思考圆心E还可在AB边上，如图6.

图 10 引入区域信息模型无震序列误报检测

3.2 对比及分析

用STA/LTA方法和未引入区域特征的单站方法对多段无震时间序列进行检测. STA/LTA方法对振幅变化敏感，但易受噪声干扰，STA记录的主要是振幅的瞬时变化，LTA的作用是衡量周围噪声的强度.当地震事件发生时，STA/LTA值明显发生变化，当STA/LTA值超过阈值，则判断此时刻为地震到时. 如果阈值设定过大则会遗漏地震事件，若阈值设定过小，则会因为噪声或其他原因发生误报，阈值不易设定.

由于高频GPS SNR较高，STA/LTA方法用于高频GPS时间序列数据效果并不理想，对多个台站时长10 800 s的无震时序进行检测，检测一次所用时间约为0.73 ms. 图11为STA/LTA方法检测图，虚线为触发检测器标记，由图11可知STA/LTA方法在多个台站上均误报较多，具体误报次数如表3所示.

未融入区域特征的单站模型数据复杂度较低，过深的LSTM网络在训练时不易收敛. 在搭建网络时减少每层网络的神经元数量，防止过拟合. 网络架构及相关参数如表4所示. 选取最优模型对长10 800 s的多个台站无震时间序列进行检测，检测一次所用时间约为23.59 ms，无震时间序列的检测结果如图12所示. 各个台站的检测结果均出现少量误报，较STA/LTA方法有所减少，具体检测结果如表3所示.

图 11 STA/LTA检出结果

图 12 未引入区域信息模型无震序列误报检测

由表3可以看出，两种神经网络方法在多个台站的无震序列上均减少了误报数量，而本文方法在引入区域特征后误报数量再次减少，虽然相比于STA/LTA方法，神经网络模型计算所需时间变长，考虑到神经网络的参数量，计算所需时间必然会增加.

表 3 不同方法检测所用时间及误报数

表 4 网络架构及参数表

4 结束语

近些年来越来越多的连续运行的GPS接收器将记录速率提升到1 Hz，高频GPS在地震领域的应用不断发展，GPS数据不论是应用于长期位移观测还是短期地壳瞬时形变都有不错的表现，在震级估算方面更为精确. 本文基于高频GPS数据提出了一种包含区域信息的神经网络模型，在应用于较高噪声的GPS数据时能得到比传统方法更好的效果.