使用无监督机器学习技术寻找2015年科托帕希火山爆发的可能前兆

摘要

科托帕希火山自2015年4月以来活动增加，最终演变为2015年8月的轻度喷发。在这项工作中，我们利用位于喷口不足4公里的宽带地震台站的记录，包括2015年4月至12月的数据，来检测和研究低频地震事件。我们应用无监督学习方案对可能的前兆低频地震家庭进行分组和识别。为了找到这些族，我们采用了一个两阶段的过程，首先通过应用k-means算法对信号的谱密度向量进行频率含量分离，然后通过Correntropy和Dynamic Time Warping进一步分离事件的波形。因此，通过探索其时间分布和估计其事件的位置，我们发现了一个与火山活动状态相关的特殊家族。

1.介绍

火山监测是潜在和当前活火山的一项关键任务。特别是，监测火山活动状态之间的过渡是很有意义的，因为这对风险管理和缓解至关重要。监测可以通过几种方式进行，但在所有使用的地球物理信号中，地震活动是最能充分代表火山活动的信号[1］.与构造地震活动相比，火山地震活动由于震源不同的物理过程而有很大的差异，因此事件的分类和术语也各不相同[2］.一般认为，在火山环境中有5种主要的事件类型:火山构造事件(VT)、长周期事件(LP)、极长周期事件(VLP)、地震事件(TR)和混合事件(HB)。虽然VTs与火山大厦内部的裂缝有关，但LP事件作为可能的喷发前兆一直受到高度关注[3.］.脂多糖与火山内部流体运动引起的结构共振有关[4]，其大部分能量在低频段，通常为0.3至3赫兹[5或高达5hz [3.］.周期为2 - 100 s的VLPs与质量运动直接相关，可出现在火山活动状态的转换中[6］.

迅速评估火山活动的变化对于近实时地探测和识别这些地震-火山事件是必要的。然而，在人员有限或危机期间没有完全自动化监测系统的火山观测站，这种认识可能不可能或不可靠，因为大量的地震数据可能会压倒监测。这个问题可以通过使用以前记录的地震事件的目录来解决，然后训练一个系统来检测和分类这些事件。这种方法被称为监督机器学习(SML)，世界各地的许多天文台已经或正在过渡到这种操作方式。研究/实施这种方法的一些火山是加莱拉火山[7和内华达德鲁伊兹[8]在哥伦比亚、圣克里斯托瓦尔和泰利卡[9]在尼加拉瓜，科利马[10]在墨西哥，默拉皮[10]在印度尼西亚，埃特纳和斯特龙博利[9在意大利，Piton de la Fournaise [10]和欺骗岛[11在南极洲。

SML技术在火山监测系统中的应用已被证明对地震事件的可靠检测和分类是有用的，但它们需要高质量的火山地震事件目录才能稳健地实施。因此，如果火山的地震目录质量很差，或者根本没有地震目录，那么SML实现可能会很困难或不可能。质量较差的地震目录是常见的，原因有:地震台网较新、活动性较低、操作员分类较差、仪器不正确或不存在、路径或场地反褶积。

在缺乏高质量目录的情况下，可以采用其他机器学习技术来缓解火山危机、仪器和系统故障或其他物流困难。快速识别和区分地震事件的一种有效替代方法是进行无监督分类。许多研究使用无监督分类来解决一些情况，作为SML的补充工具，或者只是在数据中寻找新的模式。例如，已经有了一些工作，在一个干净的数据集下分析特定类型的事件，例如地震[12- - - - - -15]及长期事件[16]或在没有事先进行清理或标记的情况下，直接对不同类型的火山地震事件进行分类[17］.无监督学习的优点是它不需要事先带安全标签的数据(数据预处理是快得多,甚至在某些情况下删除人类偏见),也可以应用于搜索特定的模式(或事件)的家庭在数据由数据本身。缺点是，在寻找这些家庭时，它的可靠性远远低于监督学习，并且在分类完成后需要人类验证。在这项工作中，我们对2015年4月至12月发生在厄瓜多尔科托帕希火山的低频率事件进行了无监督学习过程，其中包括火山短暂的苏醒。

科托帕希火山(0.68056S, 78.4378W)位于厄瓜多尔安第斯山脉中部，是一座高约5897 m的活动层状火山。这座火山可能是美国最危险的火山，因为它靠近居民区，历史上具有破坏性的火山泥流，以及最近历史上经常爆发的火山。在过去的几十年里，火山呈现出活动增加的时期，包括地震速率增加、变形和持续的排气。2015年4月初，该火山表现出以LPs为主的一般地震活动的增加。6月4日，地震活动从瞬变事件(主要是LPs)来回转换为地震(然而这些地震不是重复LPs的合并)[18］.8月13日，短暂信号再次以加速的方式出现在14日爆发的几个小时前，这是70多年平静后的第一次爆发。这次喷发是活动的高峰，活动逐渐减弱，直到2015年11月最后一次明显的气体排放。到2016年仍有少量活动[19］.

自1989年以来，在科托帕希发现了低频率事件[20.]，并与可能的岩浆侵入有关[21］.由于在整个2015年科托帕希觉醒中不断发现低频事件，我们以两阶段的方式在这些事件中寻找家庭。首先，我们通过比较它们的功率谱(谱密度向量)来寻找相似点，然后，在得到频率含量相似的组后，我们用动态时间翘曲和Correntropy两种方法来比较波形。最后，应用这些程序，我们发现了一组与火山活动状态有关的低频事件，它们可能是前兆，而无需事先选择或分类事件，也无需对信号进行严重过滤。

2.数据处理与方法

2006年年中，在科多帕希火山两侧安装了五个地震监测站(见[22];看到图6（左）。这些台站配备Guralp CMG-40T宽带地震仪，其平坦响应介于0.02和60秒之间，并连接到Smart24D数字化仪，该数字化仪以50秒的频率对信号进行采样 Hz，并将数据发送至基多地理信息网络研究所[22，23］.我们只处理BREF站的数据，因为它是离峰顶最近(~ 3km)的站，信噪比最高。特别是，我们对频率较低的事件(但不限于LPs或VLPs，这意味着它们也可能包含高频，如混合事件)感兴趣，并记录在BREF站的垂直分量上，时间为2015年4月至12月。为了识别低频内容的事件，我们首先对每日地震道应用0.05 ~ 1hz的Butterworth四阶带通滤波器。信号滤波后，对滤波后的日地震迹线进行经典STA-LTA，通过探测发现的检测阈值为4.5的特征函数得到检测触发点，得到1655个低频含量的事件。通过上述程序，我们确保检测到的信号具有低频率内容;然而，正如后面将解释的，我们也对这些事件的高频率的应用无监督学习算法感兴趣。从这个过程中，在整个时间段内共检测到1655个事件，如图所示1．

对于检测到的信号，我们将频率分析扩展到0-10 Hz，因为高频可以揭示断裂过程[24最终，更高的频率可以帮助我们更好地区分事件。使用更宽频带的逻辑是，同时具有更高和更低频率的事件，可能与混合火山过程有关，甚至可能是浅的LPs [1，20.，并应作为流体运动的可能迹象加以考虑。聚类的第一步是根据0- 10hz频带内事件的频率相似性对事件进行分组。为了对这些事件进行分组，我们计算了它们的功率谱并对它们进行了采样，同时还计算了每个间隔的最大频率(0.2 Hz)。由于0.2 Hz的步长相对较小，因此保留了足够的频率信息，作为事件之间的区分标准。在图中，我们用不同的颜色标出了这些向量2．然后将两个向量连接起来，得到长度为100的特征向量。类似大小的特征向量已经在其他工作中构建，以无监督的方式识别事件的家族或子家族[14- - - - - -17]，并成功地确定了事件的类别[17］.信号没有被转换成速度序列，因为在执行无监督学习时，处理不会产生任何效果。另一种可能的方法是在间隔上选择平均频率功率而不是每0.2 Hz的样本。虽然这种方法给略好结果聚类后执行(稍后讨论),95%的巧合分组进行比较后发现,第一个过程是保持更好的比赛以来频率需要类似事件分组(平均可以作为一个过滤器的频率向量)。

我们通过k-means方案比较和分组特征向量[25，设科数(k)为3。在本研究中，这个数字的选择是基于两个标准，戴维斯-博尔丁指数和肘关节标准。在第一个准则中，k = 3的指数比k = 2或k = 4的指数要小，除了在前面讨论的替代方法中选择平均超频率间隔而不是0.2 Hz采样的某些情况。在第一种方法中，DB-index略高，但在规模为3的家庭中始终较慢，这支持选择第一种方法。肘关节标准也表明科数可以选择3或4，因此根据DB-index选择第一种。家庭数量越多，就会自动产生频率内容更相似的群体，但这可能掩盖了一些事件可能在总体频谱内容上不同，但可能有类似的低频信号的事实;因此，我们选择更广泛的家庭。从这1655个事件中，1075个事件以有趣的时间分布和信号组合在一起(图)1)．这个家族拥有频率范围很广的事件，但在较低的频率上有重要的内容(图)3.)，以下简称家族1。此外，家族1代表了65%的总事件，并允许稳健的统计量化。

在时域中对属1的事件进行进一步分类，因为波形特征也可以用于基于相似事件形状识别相似的物理过程。在第二个聚类阶段，我们使用两种技术(为了验证目的)比较信号:动态时间翘曲(dynamic Time Warping, DTW)和Correntropy。我们使用这两种方法是因为它们解释了数据中的(不同)相似性，这些相似性可能不是与信号之间简单的互相关相反的线性。因此，这些技术不是用于检测目的，而是用于定义之前由STA/LTA检测并通过K-means分组的信号之间的相似性。

DTW算法在语音识别中得到了广泛的应用，它通过对数据进行“翘曲”来寻找时间序列之间的最优匹配[26］.该算法背后的想法是找到最好的方法来改变一个时间序列的数据点，以产生其他的数据点，并限制序列的端点重合，如图所示4．这种扭曲的代价定义了时间序列之间不相似性的度量:序列越难匹配，两个序列就越不相似。这种非线性非相似性度量在其他情况下的监督分类方案中使用[27]以及地震学，表现优于相互关联[28］.

相关系数是互相关的一种推广，通常称为广义相关函数[30.］.而互相关只考虑数据的二阶矩(协方差)，采用高斯核，Correntropy考虑数据的所有偶数矩，并推导出可用于监督和非监督学习的不同度量[31］.诱导度规(1的定义Correntropy (Correlation Induced Metric)与高斯核 在哪里X和Y两个向量有大小吗N和σ称为内核带宽，它与CIM作为L0、L1或L2范数的空间范围有关。在以0.1为步长(该参数的大值会导致Correntropy的行为类似Correlation)进行范围从0.1到1.0的测试后，本工作将内核带宽参数设置为0.1，并在视觉上验证CIM的大值在事件对之间显示出大致不同的形状。

利用这两种方法，计算了一族1中每对信号之间的异同点。由于需要计算相同大小的向量的Correntropy，信号就其开始(自动挑选)进行对齐，然后将零填充为相同大小。DTW也适用于相同的信号进行性能比较。由于我们主要对数据的低频内容感兴趣，因此我们计算了原始数据以及0.3 - 1 Hz之间过滤的数据的不同之处。此外，对每个信号的最大振幅进行归一化，只考虑波形形状的相似性而不是事件振幅大小。最后，得到时域比较差分矩阵，并将其用于层次聚类。

采用Ward的方法对差异矩阵进行分层聚类，因为它们验证了正确定性[32］.在家族1中获得了几个集群，但有一个突出(图)5)在2015年期间,它出现在类似的方式从Correntropy和DTW差别矩阵,以下简称集群A事件集群出现94%的次程序和他们将进一步验证时间出现正值关键活动政权科多帕希火山火山在2015年的喷发(18］.A组由308个事件组成。

3.分析和结果

不同方法之间的高符合性家族分析是通过不同的方式完成的:通过观察事件的时间分布，也通过观察事件本身的形状。

用两种方法研究A类事件:时间分布和波形特征。从波形特征来看，这些事件的低频内容在2015年几乎相同(图)5)．这构成了在火山中随时间而保持的可能来源的证据，并可能表明火山活动的基本水平。另一个亮点是，在6月初到8月中旬的地震活动增加期间，A星团的事件几乎消失，但在喷发前几天以加速的方式重新出现。8月13日，就在火山喷发前几个小时，这些事件每小时出现5次以上。值得注意的是，尽管2015年5月下旬事件的日发生率高于8月，但在一小时尺度上，事件发生的速度在喷发之前要快得多。

根据A簇中事件的时间分布，我们决定通过粒子运动计算源位置来进一步探究这些事件的本质(图)6)．为此，我们只使用了在0.1 ~ 0.6 Hz频段内具有较高信噪比的25个事件。利用BREF和BTAM的交叉极化方向估计了位置。图中红色部分显示了火山爆发前发生的事件的位置6，而喷发后的事件则用蓝色表示。

在构建的特征向量上进行K-means聚类只需要几秒钟的时间，而DTW和Correntropy在描述的过程中各需要大约30分钟。

4.讨论

A群事件的时间分布及其源位置揭示了2015年8月14日科托帕希火山喷发前可能有岩浆运动到地表。从时间分布来看，这些事件持续出现在喷发前的四个海拔高度，并且事件的出现在喷发开始前几小时加速。在图中所示的逐步过程中1,我们看到的是一个“缩小”的家庭,虽然时间distribution-wise每日频率似乎足够的监控火山活动,从过去的家庭出现不稀疏集群相比,发现只有分类使用频率,以及信号本身更相似。该亚家族出现在与火山活动有关的关键时刻:就在[18就在火山爆发之前。虽然每天的速率计数不能很好地区分5 - 6月可能的假阳性活动，但每小时的日期显示了它在喷发前的高峰时间(一个小时内该家族有5个快速连续事件，而5 - 6月的最大速率为3个)。由于低频率事件与火山内部的流体运动有关，这可能表明岩浆和气体在喷发前几个月在火山中稳步上升，甚至在喷发开始时向地表供应的岩浆加速上升。在事件的空间分布上，观察到一个显著的行为:与爆发后的事件相比，爆发前的事件是高度聚集的。这可能表明火山爆发后由于源的浅化或由于火山大厦内部的物理条件发生了变化而导致源特征的变化。值得一提的是，这些地点与另一项研究中发现的VLPs是一致的，该研究显示早在2002年就有类似的事件发生[20.]。特别是，308次事件中，共有71次与2009年之前确定和研究的特定VLP具有高于0.7的互相关，该VLP与科托帕西火山可能的气体释放过程有关[21］.此外，其中38个事件与IG人员进行的LPs/VLPs目视识别是一致的，可以忽略为衰减过程的一部分。不幸的是，由于在研究过程中一些气象站关闭，没有使用更精确的方法来定位这些事件。这些反复发生的事件在火山爆发后几乎消失了。在更长的时间内研究这些事件是很重要的，因为它们可能被用作喷发的前兆，或者它们可以帮助我们更好地了解科托帕希火山的基本活动水平。

5.结论

无监督机器学习方案分两个阶段应用于2015年科托帕西火山的重新唤醒，以发现可能的前兆事件。低频信号首先根据其频谱内容进行鉴别，然后根据其波形相似性进一步分离。为了测量波形相似性，两个非线性应用程序通过动态时间扭曲和Correntic实现Roach，以形成具有层次聚类的族。这两种方法在发现集群A中的事件时表现出极大的一致性，然后根据时间事件分布和源位置进行探索。值得注意的是，Correntic近似地发现了集群A，而不需要过滤signals。此外，集群A中的事件显示出与火山活动关键阶段相关的明确差异。这项工作是将无监督学习应用于之前未标记的事件数据集的一个示例，可以为识别喷发前兆提供见解。虽然程序不完善，但可能会引入假阳性，通过快速检查属于该家族的事件，人们至少可以识别与可能的火山喷发相关的信号。需要进一步研究，因为在地震群和嘈杂的环境中，该程序可能需要进一步完善。然而，它仍然可以提供一些通过进一步探索参数，可以在以前数据有限的观测站获得验证时间很少的信息，并进行更好的调整。

数据可用性

这项工作中使用的数据的所有权属于厄瓜多尔Instituto Geofísico-EPN和JICA。该数据是部分公开的，要访问它，任何请求都应向IG-EPN提交。

利益冲突

作者声明他们没有利益冲突。

致谢

这项工作完全由Pontificia Universidad Católica del Ecuador提供资金，作为与IPGH合作的“Caracterización y Estadística de Señales Volcánicas”项目的一部分。它还获得了Instituto Geofísico-EPN和JICA合作机构的数据支持。作者要感谢参与资助和发展这项工作的机构。

工具书类

1. S. R.麦克纳特，“火山地震学和火山爆发的监测”国际地震和工程地震学手册， H. Kanamori, P. Jennings，和W. Lee, Eds。，Academic Press, California, Calif, USA, 2002.视图:谷歌学者
2. S. R.麦克纳特，《火山地震学》，地球与行星科学年度回顾，第33卷，第461-491页，2005。视图:出版商的网站|谷歌学者
3. B. Chouet， "长周期火山地震活动:其来源及其在火山喷发预测中的应用"自然第380卷第2期6572页，309 - 316,1996。视图:出版商的网站|谷歌学者
4. H. Kumagai和B. A. Chouet，“作为火山下流体成分探测器的长周期地震事件的复杂频率，”国际地球物理杂志第138卷第1期2，页F7-F12, 1999。视图:出版商的网站|谷歌学者
5. r·特朗布利地震分析的概率贡献:火山地震的类型，第3章，火山爆发的预测，2006。
6. B. Chouet，《火山地震学》，纯地球物理学与应用地球物理学号，第160卷。3-4，第739-788页，2003。视图:出版商的网站|谷歌学者
7. M. Bicego, C. Acosta-Munoz, M. Orozco-Alzate，“基于隐马尔可夫模型的生成嵌入地震火山信号分类”，地球科学与遥感学报第51卷第1期6、pp. 3400-3409, 2013。视图:出版商的网站|谷歌学者
8. M. Orozco-Alzate, M. E. García, R. P. Duin，和C. Castellanos，“内华达德鲁伊斯火山地震信号的不同分类”，地球科学研究杂志，第10卷，第5期。2，页57-66,2006。视图:谷歌学者
9. l·古铁雷斯Sistema de Detección y Clasificación de señales sísmico-volcánicas utilizdo Modelos Ocultos de Markov (HMMs): aplicación a volcano activos de尼加拉瓜和意大利[博士论文]，格拉纳达大学，格拉纳达，España, 2013。
10. a . Boue数据挖掘与火山喷发预测[博士论文]，法国格勒诺布尔大学，格勒诺布尔，2015。
11. C. Benitez, J. Ramirez, J. Segura等人，“南极洲欺骗岛基于hmm的火山持续监测”，地球科学与遥感学报，第45卷，第1期，第138-146页，2007年。视图:谷歌学者
12. H. Langer和S. Falsaperla，“斯特龙博利火山(意大利)的火山震颤的长期观察:概要”，纯地球物理学与应用地球物理学，第147卷，第1期，第57-82页，1996年。视图:谷歌学者
13. A. Messina和H. Langer，“用KKAnalysis-A无监督分类软件程序对埃特纳火山(意大利)地震数据的模式识别”，电脑与地球科学，第37卷，第2期7, pp. 953-961, 2011。视图:出版商的网站|谷歌学者
14. H. Langer, S. Falsaperla, M. Masotti, R. Campanini, S. Spampinato，和A. Messina，“用于意大利埃特纳火山地震数据的监督和非监督模式分类技术概要”，国际地球物理杂志，第178卷，第2期2, pp. 1132-1144, 2009。视图:出版商的网站|谷歌学者
15. R. Carniel, L. Barbui, and A. D. Jolly，“通过自组织地图(SOM)分析探测动力机制:以2006年3月新西兰克马德克弧线拉乌尔岛潜水喷发为例”，意大利意大利大地之舞第54卷第5期1, pp. 39-52, 2013。视图:谷歌学者
16. A. M. Esposito, F. Giudicepietro, L. D’auria等，“使用自组织地图对斯特龙博利火山非常长周期事件的无监督神经分析”，美国地震学会公报第98卷第1期5，第2449-2459页，2008。视图:出版商的网站|谷歌学者
17. J. C. Anzieta和J. C. Jiménez，“厄瓜多尔通古拉瓦火山的火山-地震事件的非监督分类”，刊于IEEE第二届厄瓜多尔技术分会会议(ETCM)论文集，第1-6页，加利福尼亚州，美国，2017。视图:谷歌学者
18. S. Hidalgo, J. Battaglia, S. Arellano等，“通过综合地球化学和地震观测揭示2015年科托帕希火山喷发的演化”，地球化学、地球物理学、呈规则第19卷第2期7, pp. 2087-2108, 2018。视图:出版商的网站|谷歌学者
19. J.B.Johnson，M.C.Ruiz，H.D.Ortiz等人，“由科托帕西火山口产生的次声托尼洛，”《地球物理研究快报第45卷第5期11, pp. 5436-5444, 2018。视图:出版商的网站|谷歌学者
20. M.Ruiz、B.Guillier、J.-L.Chatelain、H.Yepes、M.Hall和P.Ramon，“厄瓜多尔科托帕西火山观测到的地震活动的可能原因，”《地球物理研究快报，第25卷，第2期13，页2305-2308,1998。视图:出版商的网站|谷歌学者
21. I.Molina、H.Kumagai、A.GarcíA-Aristizábal、M.Nakano和P.Mothes，“厄瓜多尔科托帕西火山伴随长周期信号的长周期事件的源过程，”火山学与地热研究第176期1，页119-133,2008。视图:出版商的网站|谷歌学者
22. H. Kumagai, H. Yepes, M. Vaca等人，“加强厄瓜多尔火山监测能力”，美国地球物理学会学报第88期23，页245-246,2007。视图:谷歌学者
23. H. Kumagai, M. Nakano, T. Maeda等，“使用确定性和随机方法的活火山宽带地震监测”，地球物理研究第115卷第1期8, 2010。视图:谷歌学者
24. p•鲍曼Ed。新地震台实践手册(NMSOP-2)德国地球科学研究中心，2012。
25. T. Hastie, R. Tibshirani, J. Friedman，统计学习的元素，数据挖掘，推断和预测，施普林格，柏林，德国，第二版，2009。视图:MathSciNet
26. D. J. Berndt和J. Clifford，“利用动态时间扭曲来寻找时间序列的模式”，发表于KDD研讨会的会议记录，第10卷，第359 - 370,1994年。视图:谷歌学者
27. R. J. Kate，“使用动态时间扭曲距离作为改进时间序列分类的特征”，数据挖掘和知识发现，第30卷，第2期2, pp. 283-312, 2016。视图:出版商的网站|谷歌学者|MathSciNet
28. T. D. Mikesell, A. E. Malcolm, D. Yang, M. M. Haney，“估计地震相位延迟方法的比较:尾波干涉测量的数值例子”，国际地球物理杂志，第202卷，第1期，第347-360页，2015年。视图:出版商的网站|谷歌学者
29. A. Yurtman和B. Barshan，“使用可穿戴运动传感器单元动态时间扭曲的物理治疗练习的检测和评估”信息科学与系统学报， pp. 305-314，施普林格国际出版，2013。视图:谷歌学者
30. I. Santamaría, P. P. Pokharel，和J. C. Principe，“广义相关函数:定义、性质和对盲目均衡的应用”，IEEE信号处理汇刊第54卷第5期6 I, pp. 2187-2197, 2006。视图:出版商的网站|谷歌学者
31. R. Chalasani和J. C. Principe，“自组织映射与相关性诱导度量”，在国际神经网络联合会议论文集(IJCNN’10)，第1-6页，巴塞罗那，2010。视图:谷歌学者
32. 宫本S.，安倍R.，远藤Y.和j - i。“非欧几里得相似性度量的分层聚类的Ward方法”，刊于第七届软计算与模式识别国际会议论文集， pp. 60-63, IEEE，福冈，2015年11月。视图:谷歌学者

版权

版权所有©2019 Juan C. Anzieta等人。这是一篇发布在知识共享署名许可协议，允许在任何媒介上不受限制地使用、传播和复制，但必须正确引用原作。