一种新的基于地震波形指纹的特定地区核爆炸事件快速检测方法

弓妮 ,  商杰 ,  唐伟 ,  刘哲函 ,  王海军 ,  黄立洪 ,  韩守诚 ,  江宇

地球科学 ›› 2026, Vol. 51 ›› Issue (01) : 130 -145.

PDF (7919KB)
地球科学 ›› 2026, Vol. 51 ›› Issue (01) : 130 -145. DOI: 10.3799/dqkx.2025.234

一种新的基于地震波形指纹的特定地区核爆炸事件快速检测方法

作者信息 +

A Novel Method of Detecting Underground Nuclear Explosion at Specific Site Based on Seismic Waveform Fingerprint

Author information +
文章历史 +
PDF (8108K)

摘要

核爆炸监测是禁核试核查的关键技术.为监测全球可能发生的核试验,全面禁止核试验条约规定了一套严格的核查机制.其中国际监测系统(International Monitoring System, IMS)的波形数据实时传输至国际数据中心(International Data Centre, IDC)进行处理和分析,分别在大约1 h、4 h和6 h给出三个不同阶段的自动处理结果.对于特定地区的核爆监测,直接依赖IDC结果存在响应滞后和误检率高的问题.本文提出了一种基于地震波形指纹的快速检测方法Seisprint.该方法借鉴音频指纹识别思想,将历史核爆波形作为模板,利用滑动窗口与特征提取将连续波形压缩为多个二进制指纹,通过快速相似性匹配与聚类实现核爆事件自动检测并实时报警.采用朝鲜周边两个IMS地震台站和我国东北地区4个地震台站记录的朝鲜6次地下核试验以及历史天然地震事件数据对Seisprint进行测试.Seisprint生成的指纹可以有效区分核爆与非核爆信号,且具有较强的抗噪性;可在数分钟内完成多个地震台站24小时连续波形数据的处理,实现核爆事件的快速准确检出.结果表明,Seisprint可提高特定地区核爆事件监测的时效性和准确性.

Abstract

Nuclear explosion monitoring is a critical technology for nuclear test ban verification. In order to monitor potential nuclear tests worldwide, the Comprehensive Nuclear-Test-Ban Treaty (CTBT) establishes a rigorous verification regime. Within this framework, waveform data from the International Monitoring System (IMS) are transmitted in real time to the International Data Centre (IDC) for processing and analyzing. The IDC delivers automated processing results at three stages: approximately 1 h, 4 h and 6 h after data acquisition. For region-specific nuclear explosion monitoring, direct reliance on IDC results faces challenges of delayed response and high false detection rates. To address these limitations, we propose Seisprint, a rapid detection method based on seismic waveform fingerprints. Inspired by audio fingerprinting techniques, Seisprint utilizes historical nuclear explosion waveforms as templates. Continuous seismic waveforms are compressed into multiple binary fingerprints through sliding-window feature extraction. Automated nuclear event detection and real-time alerts are achieved via rapid similarity matching and clustering. The method was tested using data from two IMS seismic stations near North Korea and four seismic stations in Northeast China, covering six historical underground nuclear tests and natural seismic events in North Korea. The fingerprints generated by Seisprint effectively distinguish nuclear explosions from non-nuclear signals and demonstrate strong robustness against noise. The method processes an entire day of continuous data from multiple seismic stations within just a few minutes, enabling rapid detection of underground nuclear explosion events. Results confirm that Seisprint promote the timeliness and accuracy of underground nuclear explosion detection in specific area.

Graphical abstract

关键词

CTBT / 核爆炸监测 / 地震波形指纹 / 局部敏感哈希 / 自动处理 / 地震学.

Key words

CTBT / nuclear explosion monitoring / seismic waveform fingerprint / locality⁃sensitive hashing / automatic data processing / seismology

引用本文

引用格式 ▾
弓妮,商杰,唐伟,刘哲函,王海军,黄立洪,韩守诚,江宇. 一种新的基于地震波形指纹的特定地区核爆炸事件快速检测方法[J]. 地球科学, 2026, 51(01): 130-145 DOI:10.3799/dqkx.2025.234

登录浏览全文

4963

注册一个新账户 忘记密码

0 引言

1996年9月24日,全面禁止核试验条约(Comprehensive Nuclear⁃Test⁃Ban Treaty, CTBT)在纽约联合国总部颁布,CTBT禁止任何国家或个人进行任何形式的核爆炸试验(Asada, 2002Bauer and O’Reilly, 2016).为了监测全球可能发生的违约核试验,CTBT给出了一套严格的核查机制,其中包括由337个核查设施组成的国际监测系统(International Monitoring System, IMS),以及负责IMS数据处理和分析的国际数据中心(International Data Centre, IDC).IMS中的地震、水声和次声台站监测数据实时发送到IDC进行处理和分析,分别在事件发生后大约1 h、4 h和6 h给出三个不同阶段的自动处理公报(SEL1, SEL2, SEL3).随后,IDC分析员对自动处理结果SEL3进行校正、完善和质量审核,形成审核事件公报(Reviewed Event Bulletin, REB),上述自动处理结果和人工分析结果作为IDC的数据处理产品对授权用户发布(Bobrov et al., 2012).

相对于人工分析结果,IDC波形数据自动处理系统中事件的误检率大约是50%,事件漏检率接近30%(Arora, 2012).作为全面禁止核试验条约组织重要的技术部门,IDC一直通过优化系统参数、集成新技术新方法来不断改进数据处理系统的性能.Arora et al.(2009)在科学技术大会ISS09期间首次提出NET⁃VISA(Network Processing Vertically Integrated Seismic Analysis),这种方法利用历史事件的特征参数和关联台站的信号检测参数建立概率统计模型,通过贝叶斯理论筛选最佳台站组合,系统地给出了一种新的类似于机器学习的台网事件关联算法.经过十余年的不断完善和测试(Arora, 2012Arora et al., 2013),NET⁃VISA技术已于2018年1月正式集成应用在了IDC数据处理系统中(Le Bras et al., 2021).

Ali et al. (2025)的评估结果表明,NET⁃VISA的应用使得扫描事件总数增加了17.90%,平均每日新增 7 个事件;相比于传统依赖于信号特征参数和震相传播走时规律进行聚类的全球关联(Global Association, GA)(李健,2021)算法,NET⁃VISA显著降低了漏检率.但在实际应用中,NET⁃VISA需要足够多的历史事件和人工分析结果来训练事件和信号特征参数模型.针对特定地区,当可用历史事件数量较少时,NET⁃VISA算法性能下降.此外,IDC自动处理系统主要面向全球地震事件监测,为关联更多台站数据通常需等待约1 h才输出检测结果,其时效性难以满足特定地区实时监测任务的需求.鉴于特定地区的地震事件在相同台站记录中往往具有相似波形特征,一些学者提出了基于波形相似性的检测方法,例如波形互相关算法和指纹与相似性阈值检测法(Fingerprint and Similarity Threshold,FAST)(Yoon et al., 2015)等.波形互相关(Cross⁃Correlation)或模板匹配(Matched Filtering)是地震学中常用的波形检测和相似性分析方法,通过计算归一化相关系数并选择适当阈值,可以在连续波形数据中识别与模板波形相匹配的部分(Senobari et al., 2018).互相关方法已成功用于朝鲜核试验及大型余震序列的检测(Gibbons and Ringdal, 2012),提高了震相拾取和事件构建的准确性.当前IDC通过预定义高质量的主事件作为波形模板,寻找相似信号,并对检测出的信号进行关联形成事件列表XSEL,以补充SEL中的漏检事件(Bobrov et al., 2014,2016).然而,单台站互相关在噪声或复杂环境下误检率较高(Li and Zhan, 2018Muir et al., 2023),对长时间的连续数据或海量模板做滑动相关需要巨大的计算量和 I/O,在大规模应用上成为瓶颈(Gibbons and Ringdal, 2006).

FAST是斯坦福大学提出的一种基于区域台网的重复地震事件快速检测技术(Yoon et al., 2015Bergen and Beroza, 2019),将波形转换为紧凑指纹,通过哈希索引技术快速在数据库中进行相似性搜索,存储和检索都不需加载原始波形;相较于互相关算法,内存占用和I/O 成本显著降低,且计算复杂度为近线性复杂度.FAST凭借计算效率高、可扩展性强的优势,在连续数据处理时能够自动检测到重复出现的地震事件,适用于实时地震监测.Rong et al.(2018)将FAST检测管道应用到加州地区,显著提高了重复事件的检出效率,尤其能发现大量之前未记录的微小地震.然而,FAST算法只能检测到特定时间段内重复发生的事件,而无法检测到仅发生一次的新事件.核爆事件往往多年一次,在长时间窗内缺乏同源模板,使用FAST 算法存在漏检风险.

为此,本文提出了一种新的基于地震波形指纹的快速检测方法Seisprint(Seismic Waveform Fingerprint),将其应用于特定地区核爆炸事件监测中.该方法在FAST算法的基础上,借鉴通信领域中的帧同步理论(Sugita and Ito, 2015),将历史核爆事件的波形作为帧头,通过时域拼接嵌入待处理的连续地震记录数据流前端.在通信系统中,帧头用于标记数据帧的起始并携带识别信息;在Seisprint方法中,帧头可作为核爆事件标识,只有当待处理的连续地震记录中出现与该帧头波形相匹配的信号(即帧头时间戳重现)时,才会被初步判定为候选信号.为了提高事件检测的准确性,实际应用时采取多台站联合判据,只有当多个台站记录中同时检测到与帧头相匹配的波形时才形成候选事件.Seisprint将FAST与模板匹配相结合,是一种模板驱动的波形指纹相似性检测法,无需像 NET⁃VISA 那样依赖大规模历史事件库预训练,只需要少量模板即可实现快速检测.本文利用朝鲜核试验场周边的两个IMS地震台站(韩国KSRS台阵的子台KS31、俄罗斯USRK台阵的子台USA0B)和我国东北地区的4个国内地震台站(MDJ、CN2、KDN、YNB)记录到的朝鲜6次核试验以及发生在试验场周边的历史地震数据,对Seisprint方法进行测试评估,验证其在特定地区核爆事件检测中的有效性.

1 地震数据集

本文采用的朝鲜地区地震数据集来源于朝鲜核试验场区周边6个不同方位的地震台站(CN2、KDN、YNB、MDJ 、USA0B、KS31)的观测记录.数据集包含两部分:(1)离线事件库,用于评估指纹的鉴别性与抗噪性;(2)连续波形数据,用于检验Seisprint算法在实际流式数据中的端到端检测效果.

在构建离线事件库时,本文选取了朝鲜2006-2017年间进行的6次地下核试验事件(用NKT1⁃NKT6表示)、2018-2024年发生在朝鲜试验场周边(40°~42°N、128°~130°E)的134个天然地震事件(源自CTBTO REB目录,https://swp.ctbto.org/web/swp/reb,震级分布见图1,位置分布见图2),以及经人工校验的150段背景噪声记录,构建出包含核爆、天然地震和背景噪声的综合事件库.所有事件波形均为台站的垂直向通道记录(Z通道),截取自P波初至前10 s至后140 s的时窗(总长150 s);背景噪声数据则从无事件时段随机截取等长片段.

连续波形数据部分的选取了共24天连续波形,包含以下时间段:(1)NKT3、NKT5、NKT6三次核试验当天的24小时连续波形记录,用于验证对已知核爆事件的检测能力;(2)2016年9月1日至8日、9月10日至20日共19天的连续波形记录,用于评估一般背景下的事件检测能力;(3)NKT6后发生的3次余震事件(用PEV1⁃PEV3表示,时间分别为NKT6事件后约 8 min、2017年9月23日和 10月12日)当日24小时的连续波形记录,用于测试Seisprint方法对特定区域内余震事件的检测能力.需要说明的是,由于 NKT6 与PEV1 发生在同一天,二者共用同一段连续波形数据.6次核试验和3次余震事件参数见表1(来源于CTBTO REB目录).

2 Seisprint方法

Seisprint算法是将FAST和模板匹配相结合的一种模板驱动的波形指纹相似性检测方法,处理流程主要包括5个步骤,分别是:预处理、指纹提取、相似性搜索、后处理和事件筛选.各流程相互衔接,协同完成核爆事件的自动检测,其工作流程如图3所示.

2.1 预处理

Seisprint相较于FAST的一项关键改进在于历史核爆波形作为帧头驱动核爆炸事件检测.地下核爆炸在固定场地重复进行时,其震源机制具有高度一致性(各向同性膨胀源),即震源矩张量以各向同性分量为主,对应源区体积的瞬时增加,理论辐射图呈球对称,仅辐射等幅同极性的 P 波而不直接激发 S 波(Richards and Kim, 2005),假设传播路径的地壳结构相对稳定,爆炸产生的地震波在相同地震台站记录中具有相似性.每次核试验的当量不同(2006 年首次试验当量不足1 kt,而 2017 年试验推测超过100 kt)(赵连锋等,2017),震源埋深也有所变化,这直接影响地震波的绝对振幅和高频能量含量:较大当量和较深埋深的试验通常产生更强的低频能量,且高频部分受到不同程度的衰减,从而导致同一台站记录的核爆炸波形有所差异(图4),但是整体呈现出核爆炸的典型地震记录特征.其相同的波形特征为:(1)P波起始陡峭、Lg波相对较弱:爆炸源在起爆瞬间释放巨大能量,产生高频冲击波,导致P波初至尖锐且振幅陡增;(2)P/S波能量分配异常:P波高频成分相较于S波更丰富,且P波幅值大于S波幅值;(3)短周期Rayleigh(Rg)波发育:Rg为Rayleigh基模在短周期端的导波,常在区域‒近区距离内表现为显著能量包.在MDJ台站记录的6次核爆事件的波形中,Rg波主能量稳定落在3~5 s的周期带,显示出典型频散:群速度随周期变化,与区域尺度的基模Rayleigh波频散规律一致(Kafka, 1990).

Seisprint算法中,首先对输入的待检测连续波形数据进行预处理:(1)模板拼接:选取已知历史核爆事件作为模板,将其在多个台站垂直向通道记录的模板波形窗口WP=Parrival-Tbefore:Parrival+Tafter拼接至待检测的连续波形数据开头.其中,Parrival为P波到时,TbeforeTafter分别为P波到时前后截取的窗口长度.(2)时间同步:由于添加帧头,需要将待检测波形的起始时间统一向前调整Tbefore+Tafter(假设模板波形窗口长度为150 s,待检测数据为2016年9月9日的24小时连续波形,那么将拼接后的波形数据起始时间向前调整150 s,设置为2016年9月8日23:57:30).

2.2 指纹提取

指纹提取是Seisprint检测流程的核心环节,图5是KS31台记录的2021年6月13日目标区域内发生的一次天然地震事件与2009年5月25日朝鲜第二次核爆事件(NKT2)的波形(图5a、5b)及其指纹提取过程.

(1)首先通过短时傅里叶变换(STFT)将原始地震波形转换为时频谱图,将一维波形映射到“时间‒频率”平面,捕捉震相能量随时间演化的局部谱结构,计算公式如下所示:

         STFTxt,ω=-+xτWτ-te-jωτdτ,  

其中,t表示分析时刻,ω为角频率,xτ表示τ时刻的时域地震信号,Wτ-t表示地震信号在τ时刻的滑动窗函数,e-jωτ为傅里叶基函数.STFT 的核心思想是在时间轴上“挪动”一个窗,把以t为中心的一段信号拿出来做局部傅里叶变换,因此STFTxt,ω表示“在时刻t附近的频谱”.图5c、5d展示了天然地震与核爆的时频谱对比,直观反映了二者在不同频段的能量分布差异.

(2)使用滑动窗口将时频谱图分割生成重叠频谱片段,以增强对时间漂移与小尺度噪声的鲁棒性;然后对每个片段进行二维离散Haar小波变换提取特征.Haar小波变换作为一种对函数进行层次分解的数学工具,能够将图像描述为整体形状与逐级细化细节的叠加,在图像检索领域应用广泛(Mallat, 2009),详细计算过程见附件1.通过小波变换,把时频图分解为整体轮廓与逐级细节,如图5e、5f所示.

(3)使用中位数绝对偏差(MAD)对小波系数进行标准化,使得对异常脉冲与台站噪声更稳健:

         aj=medianiLji   bj=medianimedianiLji-aj    L'jiLji-ajbj

其中,median表示计算中位数,Lji为第i个指纹的第j个小波系数,L'ji为标准化后的小波系数;ajbj分别表示第j个小波系数在全部指纹集中分布的中心位置和离散度度量.

(4)为了增强对噪声的鲁棒性,对小波图像进行压缩,选取绝对幅值最大的K个标准化小波系数,保留系数的符号(正为+1,负为-1),其余系数置零.如图5g、5h所示,其中黑色表示-1、白色表示1、灰色表示0,可以看出,小波图像中能量不突出的部分均变成了灰色.选Top⁃K系数并仅保留符号的目的是为了强制保留“最具判别力”的结构;将高维连续系数压缩为稀疏符号图,显著降低存储与匹配计算量,同时减弱了幅值尺度差异对匹配的影响.

(5)对压缩后的小波图像进行二进制编码生成指纹向量(+1→01, -1→10, 0→00),以支持近线性复杂度的大规模相似度搜索.生成的指纹如图5i、5j所示,其中黑色表示0、白色表示1.

2.3 相似性搜索

稀疏的二进制指纹可以视为一个集合:取值为 1的系数索引构成集合 A,那么另外一个指纹对应集合B,使用 Jaccard 相似度度量两个集合的重叠程度,定义为:

         J(A,B)=ABAB.

基于Jaccard相似度,采用局部敏感哈希(Locality Sensitive Hashing, LSH)(Gionis et al., 1999)算法实现快速指纹匹配.LSH是一类为相似度搜索设计的哈希方法:让“越相似的对象在同一哈希桶中碰撞的概率越高、越不相似的碰撞概率越低”,从而用少量桶筛出候选近邻.

本研究中,指纹提取步骤生成的指纹具有以下特性:相似波形的指纹相似度显著高于非相似波形的指纹相似度.因此,采用基于LSH的相似性搜索算法,构成一种基于哈希空间的近似聚类机制,即相似的波形指纹倾向于多次碰撞,聚集于相同哈希桶中,形成稠密的指纹簇,而不相似的指纹则由于哈希特征离散化,在哈希空间中分布分散,难以跨越相似度阈值,最终被排除在候选列表之外.具体实现分为两步.

(1)Min⁃Hash降维:首先利用Min⁃Hash生成哈希签名,Min⁃Hash具有重要的LSH属性,即两个指纹fp1fp2的哈希值相等的概率等于其Jaccard相似度:

         Prhfp1=hfp2=Jfp1,fp2,    

其中,hfp1表示生成指纹fp1的哈希值.Min⁃Hash具体算法细节见Cohen et al. (2001),用p个独立 Min⁃Hash 函数H=h1x,h2x,,hix把高维稀疏指纹压缩成长度为p的哈希签名(p个整数组成的数组).

(2)LSH索引与搜索:构建哈希索引库进行近似相似性搜索.将长度为p的签名切成l个子向量(band),每个band长度为m=p/l,把band“完全相同”的指纹索引映射到相同哈希桶中.如果两个指纹落在同一哈希桶的哈希表个数vv为相似度阈值),就会判定其相似.匹配搜索时,用相同切分与哈希把查询项落入这l个哈希表的若干桶中,并收集桶内候选指纹.

图6展示了将波形指纹存储到哈希数据库的过程.哈希数据库包含多个哈希表,每个哈希表包含多个哈希桶.如果两个指纹的band相同,就会被哈希到相同哈希桶中,如图中哈希表2、哈希表(l-1)和哈希表l所示,否则被哈希到不同的哈希桶中,如哈希表1所示.图中两个指纹在3个哈希表中被哈希到相同的哈希桶,通过比较3与v的大小来判断这两个指纹是否相似.

相似性搜索步骤输出的是波形指纹对构成的相似度矩阵,如图7所示,由于连续地震波形会生成许多个指纹,因此相似指纹对数量也很庞大,需要通过后处理步骤进行指纹对的聚类.

2.4 后处理

后处理模块通过三级结构化流程将海量指纹对提升至事件层级:

(1)单台站事件提取:在相似度矩阵中识别沿对角线分布的高相似度簇——该模式表征具有连续时间关联性的相似指纹序列.基于事件持续时间的物理约束,设定时间邻近阈值合并重复检测窗口,生成初步台站级事件对列表.

(2)多台站伪事件关联:利用重复地震事件的固定发震时刻差特性(同一台站的波形走时差异可忽略),设置台站检测阈值(记录到同一事件的最小台站数).仅当多个台站独立检测到具有一致发震时刻差的事件对时,才关联构建候选事件.

(3)事件解析与重建:融合同一物理事件的多重检测结果,基于综合置信度指标解析最终事件清单,显著压缩数据规模,提升了检测可靠性,克服了单指纹对仅依赖相似度指标的局限,且有效抑制了单台误检.

2.5 事件筛选

为提升检测结果的可信度,在最终处理阶段,基于核爆模板的先验信息及台站空间分布,引入时空约束与波形互相关算法对候选事件进行筛选.

(1)时空约束:利用赤池信息准则(Akaike Information Criterion, AIC)方法优化检测信号初至(赵大鹏等, 2012),该方法通过最小化时间序列分段模型的信息损失实现高精度到时检测:

         AICk=kln varx1,k+N-k-                              1ln varxk+1,N,             (5)

其中k表示假设的初至位置(第k个采样点),x为地震波形数据,N为窗口长度, varx1,k计算从第一个采样点到第k个采样点的方差,表示初至前背景噪声段的能量水平;varxk+1,N计算从第k+1个采样点到第N个采样点的方差,表示包含地震信号的能量水平,P波初至时刻对应AICk的全局最小值点kmin.通过对比震前和震后段信号的统计特征,AIC 方法能有效识别震相初至的突变点,实现高精度自动到时拾取.依据区域速度模型计算监测区域内事件到周边各台站的理论到时差,设定时空约束条件:事件到各台站的实际观测到时差需在理论到时差允许的误差范围内,剔除不满足该条件的事件.

(2)波形互相关验证:以历史核爆事件在各台站的波形记录作为模板,计算候选事件波形与对应台站模板波形的归一化互相关系数(Cross⁃Correlation Coefficient, CC),计算公式如下:

        CC=i=1nTi-T¯Ci-C¯i=1nTi-T¯2i=1nCi-C2,

其中,Ti为模板波形的第i个采样点,Ci为候选事件的第i个采样点,T¯C¯为对应波形的平均幅值,n为波形采样点总数.保留互相关系数高于阈值的事件,筛除与模板互相关系数偏离较大的误判事件.

3 结果分析

Seisprint关键参数设置见表2,其中由于MDJ台站的采样率为20 Hz,为了保持所有台站的参数一致,其余5个台站统一重采样到20 Hz.接下来对影响算法性能的重点参数进行详细说明.

稀疏度参数K的选择:K表示稀疏度参数,即压缩后的非零小波系数的个数.K决定了指纹的稀疏程度与判别细节量,随着K增大,指纹更稠密,存储紧凑性下降,相似度搜索步骤的运行时间增加;随着K减小,指纹更稀疏,存储与检索更高效,但可能牺牲对于弱信号或细粒度差异的敏感度.Bergen and Beroza(2019)指出,对于FAST而言,在4 096维二进制指纹中,取K=300~500(对应稀疏度7%~12%)通常最优,可在效率与灵敏度之间取得良好折中.本研究中的Seisprint与FAST的目标存在区别,FAST更强调“检出微弱地震”;Seisprint除了要将核爆信号从噪声中稳健检出,还强调其与天然地震信号的鉴别性.为保留更多鉴别性细节,将稀疏度提高至约20%,即K=800(4 096维),该取值在指纹判别力与检索速度之间取得了较好的平衡(经数据集验证,见图8).

哈希表个数l与哈希签名长度p的选择:LSH 的精度和召回率主要通过这两个关键参数调节.由2.3小节可知,长度为m的band映射到l个独立哈希表;两个指纹在m个随机选取的非零位置上完全一致,才能被哈希到相同的哈希桶中.因此m越大,碰撞概率越低,候选项减少,召回率下降;m越小,碰撞概率越大,召回率更高,但精度下降,且运行时间增加.为了在召回不下降的前提下提高精度,可增加哈希表个数l,但代价是更高的内存与索引维护成本.Rong et al. (2018)指出哈希表个数l通常取100,m的经验建议为:短时长数据集(天-周)时m=4;较长的数据集(月-年)m=5;最长数据集(5~10 年)m=6.在本研究的数据集中,最长为1天的连续波形,因此取m=4 (p=400).

相似度阈值v的选择:为尽量保留潜在相似对,在相似性搜索阶段建议将初始阈值设得偏低,常用v=2作为起点,生成更充足的候选指纹集合;再于后处理中进一步修剪.

3.1 离线数据库测试

核爆事件自动检测的关键问题是需要具备鲁棒性和鉴别性的双重优化性能.鲁棒性表示为在较低信噪比下将信号检出的能力,而鉴别性表示将核爆信号与非核爆信号区分的能力.为了验证Seisprint指纹提取算法的鉴别性,在离线事件库中,提取核爆、天然地震和背景噪声波形的指纹,选取2016年1月6日核爆事件(NKT4)波形作为模板,计算其余指纹与模板指纹相似度分布,如图8所示.其中,由于核爆事件和天然地震事件不一定可以同时关联到6个台站(例如NKT1仅关联到了MDJ台站),因此图8中不同台站的核爆事件数和天然地震事件数会出现差异.

显然,核爆指纹与模板指纹的相似度得分最高,天然地震指纹处于中间区间,与核爆指纹之间存在稳定的可分界;而背景噪声指纹得分普遍最低,基本集中在低相似度区域,形成明显的三层区分带;其中,天然地震与噪声的指纹区分度不高,甚至有的重叠在一起(如CN2、KDN台站),主要源于部分地震事件(图1)震级小于2.5,信噪比低,地震信号被背景噪声淹没,导致其指纹与噪声指纹差异不显著.相反,震级较大且波形清晰的事件可保留更明显的特征,其指纹与噪声指纹差异明显.尽管各台站的匹配得分数值分布存在差异(可能源自于仪器响应、路径效应或背景噪声差异),但核爆、天然地震和背景噪声三类指纹在不同台站的分层趋势均一致.指纹空间中的聚类验证了以下特性:对于特定地区的事件,应用Seisprint指纹提取算法,相同震源类型(6次核爆炸)的波形生成高度相似的指纹特征,而不同震源类型(核爆炸与天然地震)则产生显著差异化的指纹,这说明Seisprint指纹提取算法对区分来自相同区域的不同震源类型的信号具有较好的鉴别性.

为评估Seisprint指纹提取算法对噪声干扰的鲁棒性,本研究采用向模板波形添加随机噪声的策略,并计算加噪指纹与模板指纹间的相似度.假设模板波形为xRM,在时域中用M个样本表示,噪声波形n(j)RM取自同台站的噪声数据集,F(x)0,1D表示指纹提取操作,通过缩放因子α控制信噪比(SNR):

           SNR=10lg(αx2n2),                                     

指纹精度Accuracy定义为加噪指纹与模板指纹的Jaccard相似度:

           Accuracyi=JFαix+ni,Fx.

基线相似度Baseline定义为噪声指纹与模板指纹间的Jaccard相似度:

           Baselinej=JFnj,Fx.

考虑到朝鲜地区历史核爆炸事件当量较大且选用近台数据,添加噪声对原始波形的影响可能不显著.为验证Seisprint算法对小当量核爆信号的检测鲁棒性,实验选用KDN台站记录的一次小震级天然地震事件波形为模板(事件参数与多台站波形见附件2),对其添加不同信噪比的噪声进行测试.在SNR=2 dB时(图9),添加的噪声几乎淹没了有用信号,即使采用1~2 Hz带通滤波,人工分析依旧存在困难.图10是应用Seisprint算法提取模板波形、加噪波形(SNR=2 dB)与噪声的指纹,获得指纹精度和基线相似度的分布.在低信噪比条件下,指纹精度与基线相似度的分布依然呈现出较为明显的分界,位置大约在Jaccard相似度值为0.2附近.结果显示,对于低信噪比的地震波形记录,Seisprint提取的指纹能够有效区分噪声和地震信号.

在指纹提取的基础上,结合相似性搜索算法,进一步采用Precision⁃Recall(PR)曲线对Seisprint算法的检测性能进行综合评估.对于给定的相似度检测阈值,将模板指纹作为查询对象,利用LSH算法在包含加噪指纹和噪声指纹的数据库中进行快速相似性搜索.定义真阳性(True Positive, TP)为被正确识别为与模板指纹相似的加噪指纹,假阳性(False Positive, FP)为被错误识别为与模板指纹相似的噪声指纹.基于TP和FP的统计可给出精度 (Precision)与召回率(Recall)这两个指标:

           Precision=TPTP+FP,                                      
           Recall=TPNtotal,                                                      

其中,Ntotal表示测试集中加噪指纹的总数量(即所有待检出的正样本数量).精度衡量系统返回结果中正确匹配的比例,召回率衡量系统找出所有真正匹配项的能力.

PR曲线描绘了当检测阈值在从最严格(相似度要求高)到最宽松(相似度要求低)的范围内连续变化时,系统所对应的Precision与Recall之间的关系.统计了不同SNR条件下的PR曲线,结果如图11所示.实验表明,当SNR5 dB时,在满足召回率大于80%的条件下,能够实现90%以上的检测精度.选择合适的检测阈值,系统能够较好地实现误检与漏检之间的权衡,这充分体现了Seisprint算法对噪声干扰的强鲁棒性.

3.2 连续波形数据测试

为验证Seisprint系统在实时监测场景下的端到端检测效率与可靠性,开展基于连续地震波形数据的测试.选取离线数据库中2016年1月6日朝鲜核爆事件的波形片段,构建“帧头”,将其拼接至对应台站的连续波形数据起始端,以驱动后续的核爆指纹匹配检测流程.如图3的检测流程,根据是否与帧头匹配,可以将事件分为波形相似的天然地震事件和候选事件;对于候选事件,算法启动事件筛选程序,满足筛选条件的则触发警报,表明监测区域内可能存在核爆炸事件,需由分析员进一步人工复核,否则不触发警报.本文利用Seisprint算法对24天的连续波形进行测试,共检测出22个事件,包括16次波形相似的天然地震事件和6次候选事件,经过事件筛选,候选事件中触发3次警报.

针对这3次报警事件和另外3次未触发警报的候选事件,绘制其在多台站的波形对比图(图12).采用AIC算法检测信号初至,报警事件在MDJ、CN2等5个台站的初至信号到时差不超过10 s,并且YNB台检测到的朝鲜核试验场的初至信号到时都早于其他台站(图2显示YNB台站距离监测区域最近,因此理论走时最短).相比之下,其他3个候选事件在这些地震台站记录上的初至信号到时差不具备上述特征,通过人工分析并与标准地震目录对比,确定候选事件均为监测区域外的天然地震事件,详细事件参数见附件3.

为进一步验证检测结果,对上述6个候选事件在多个台站分别截取检测时刻±60 s的波形段,与帧头波形进行互相关计算,并提取最大互相关系数CC表3).设定判别标准为:在6个参与台站中,需满足≥5个台站的CC  0.5,才会判定为核爆事件.结果显示,真实核爆事件全部满足判定标准,而候选事件的平均CC仅0.20~0.27,且所有台站数据的互相关系数CC<0.5.通过多台站的初至信号到时差与波形互相关分析,可确认3个候选事件不符合报警条件.

表4给出了Seisprint、NET⁃VISA(王娟等, 2021)和波形互相关方法对于NKT3、NKT5和NKT6当天连续波形数据的检测结果,并与CTBTO REB标准目录进行了对比.其中Seisprint和波形互相关均采用2016年1月6日朝鲜核爆炸事件(NKT4)的波形作为模板,波形互相关方法参数设置如下:带通滤波频段为 0.8~4 Hz,模板长度 20 s,单台站检测(MDJ 台站),采样率 20 Hz,相关阈值 0.5.检测结果显示,3种方法均成功检出3次核爆事件,NET⁃VISA将NKT5和NKT6均分裂成了2个独立事件,波形互相关的结果包含2次误检.

在时间误差方面,Seisprint相较于NET⁃VISA和波形互相关更大.这是因为Seisprint的检测流程中未包含震相拾取步骤.在后处理阶段的单台站事件提取步骤中,Seisprint设定将指纹合并的最大时间间隔为P波与S波的最大到时差(通常为15~40 s),并以相似度得分最高的指纹对应时间作为事件检测时间.在本研究中,指纹长度为10 s,因此理论上Seisprint的检测时间相较真实发震时刻的误差最大可达约50 s.

如果期望得到精确的发震时刻,在Seisprint流程之后可考虑添加震相拾取和定位步骤.例如,可截取检测时间前后各60 s的波形,利用匹配‒定位(M&L)方法或机器学习方法进行震相识别与定位,给出事件的精确时刻与位置.

对于NKT6后发生的3次余震事件PEV1⁃PEV3,Seisprint算法均未检出.已有研究(He et al., 2018)指出,这些事后事件在主频带与各相位(P、S、Rg/Lg)相对幅度关系上显著区别于核爆,提示其震源时间函数与机制不同;其中PEV1很可能为爆炸腔体坍塌或诱发构造地震,而PEV2、PEV3更可能属于构造型事件(Yao et al., 2018;刘俊清等, 2019).采用Seisprint方法对PEV1⁃PEV3的指纹提取过程与详细分析见附件4.因此,连续波形数据测试结果表明,Seisprint对核爆事件和天然地震事件的鉴别性较好,它能稳定检出核爆事件的同时不会将同一区域内的天然地震事件错误检出.

3.3 计算效率

为了评估Seisprint方法在实际应用中的计算效率,对其内存需求与运行时间进行了分析.其内存占用主要由以下三部分构成:(1)哈希表指针开销Ol·M,其中M为指纹数量,l为哈希表数量.每个哈希表中存储指向子指纹的指针,每个指针占用 4 字节.(2)哈希桶结构开销Ol·b,其中b为每个哈希表中哈希桶的数量.采用数组实现哈希表,每个哈希桶存储一个指向其内容链表的指针(4字节)和一个记录桶内子指纹数量的整数(4字节),共计8字节.(3)子指纹存储开销Op·M,其中p为哈希签名长度.

在本文的测试场景中,参数设定如下:每秒生成一个指纹,处理24小时波形数据并拼接 150 s帧头,总指纹数量M=86 550.哈希表数量l=100,子指纹长度(哈希签名长度)p=400,每个哈希表的哈希桶数量b=100 000.据此计算各部分内存开销:

         Ol·M:4 字节×100 ×86 550 =                            34 620 000 字节33 MB,              
         Ol·b:8 字节×100 ×100 000 =                         80 000 000 字节80 MB,                 
         Op·M:1 字节×400 × 86 550 =                              34 620 000 字节33 MB.

假设额外开销约为1 MB,那么内存需求总计约为147 MB.扩展至1年(365 d)的内存需求约为 53 GB,此内存规模在普通服务器或集群环境中可行.Seisprint的计算复杂度与FAST相似,特征提取与数据持续时间呈线性关系,相似性搜索为近线性复杂度O(M1.36).经统计,在Linux环境(Intel(R) Xeon(R) Gold 6248R CPU @ 3.00GHz,x86⁃64, 96 CPUs,1.5 MB缓存)下,处理6个台站的24小时连续波形数据的端到端耗时均在1 min 30 s以内.

当模板数增加时,Seisprint算法只需要将新的模板继续拼接至帧头,运行时间只与新增模板自身长度相关;相对于大规模连续波形,模板极短,因此带来的开销可忽略不计;而每增加一个模板,波形互相关算法的运行时间都需要增加一倍,当扩展到大规模连续数据时,Seisprint算法的计算效率要远远高于波形互相关.

4 结论

为监测特定地区的核试验活动,本研究提出了一种基于波形指纹的核爆炸事件快速检测方法Seisprint.该方法的核心流程包括:选取历史核爆炸波形作为模板帧头,拼接至连续波形数据起始处;将波形数据压缩为稀疏二进制指纹,借助哈希数据库实现快速相似性搜索;再结合时空约束与互相关阈值对候选事件进行严格筛选,仅当多个台站均检测到与模板波形相似且符合判别条件的事件时才触发警报.利用朝鲜核试验场周边6个台站记录的6次历史核爆及天然地震数据进行验证,结果表明,Seisprint可在2 min内完成多台站24小时连续数据的处理,实现了对该地区核爆炸事件的快速准确识别.与国际数据中心现有波形处理算法相比,本方法在特定地区的应用中表现出显著优势:运行效率高、内存占用低、具有良好的可扩展性,同时保持较高的检测准确性.通过替换模板帧头,该方法还可推广至其他地区的核爆炸或重复地震事件的检测中.

然而,Seisprint在指纹提取过程中对信号时频谱进行了多尺度边缘与纹理压缩,将丰富的能量分布信息转换为二值化轮廓,导致一定程度的信息损失.区域内不同类型地震事件因震源机制差异而呈现不同波形特征,当前指纹表征方式在区分细微特征方面尚存在局限.下一步改进将聚焦于更精细的波形指纹特征提取方法,并结合更高效的模式匹配算法,以提升对相似事件的区分能力与检测可靠性.

附件见https://doi.org/10.3799/dqkx.2025.234.

参考文献

[1]

Ali, S. M., Qorbani, E., Le Bras, R., et al., 2025. Interactive Analysis of the Results of NET⁃VISA, a Bayesian Inference System, in CTBTO’s International Data Centre Bulletin Production. Acta Geophysica, 73(1): 71-81. https://doi.org/10.1007/s11600⁃024⁃01398⁃0

[2]

Arora, N. S., 2012. Model⁃Based Bayesian Seismic Monitoring (Dissertation). EECS Department, University of California, Berkeley.

[3]

Arora, N. S., Russell, S. J., Sudderth, E. B., 2009. Vertically Integrated Seismological Analysis II: Inference. In: American Geophysical Union (AGU) Fall Meeting 2009. Bayesian Logic, Inc., San Francisco.

[4]

Arora, N. S., Russell, S. J., Sudderth, E. B., 2013. NET⁃VISA: Network Processing Vertically Integrated Seismic Analysis. Bulletin of the Seismological Society of America, 103(2A): 709-729. https://doi.org/10.1785/0120120107

[5]

Asada, M., 2002. CTBT: Legal Questions Arising from Its Non⁃Entry⁃into⁃Force. Journal of Conflict and Security Law, 7(1): 85-122. https://doi.org/10.1093/jcsl/7.1.85

[6]

Bauer, S., O’Reilly, C., 2016. The Comprehensive Nuclear⁃Test⁃Ban Treaty Organization (CTBTO): Current and Future Role in the Verification Regime of the Nuclear⁃Test⁃Ban Treaty. In: Black⁃Branch, J., Fleck, D., eds., Nuclear Non⁃Proliferation in International Law. T.M.C. Asser Press, The Hague. https://doi.org/10.1007/978⁃94⁃6265⁃075⁃6_6

[7]

Bergen, K. J., Beroza, G. C., 2019. Earthquake Fingerprints: Extracting Waveform Features for Similarity⁃Based Earthquake Detection. Pure and Applied Geophysics, 176(3): 1037-1059. https://doi.org/10.1007/s00024⁃018⁃1995⁃6

[8]

Bobrov, D., Bormann, P., Coyne, J., et al., 2012. CTBTO: Goals, Networks, Data Analysis and Data Availability. In: Bormann, P., ed., New Manual of Seismological Observatory Practice 2 (NMSOP⁃2). Deutsches GeoForschungs Zentrum GFZ, Potsdam. https://doi.org/10.2312/GFZ.NMSOP⁃2_ch17

[9]

Bobrov, D., Kitov, I., Zerbo, L., 2014. Perspectives of Cross⁃Correlation in Seismic Monitoring at the International Data Centre. Pure and Applied Geophysics, 171(3): 439-468. https://doi.org/10.1007/s00024⁃012⁃0626⁃x

[10]

Bobrov, D. I., Kitov, I. O., Rozhkov, M. V., et al., 2016. Towards Global Seismic Monitoring of Underground Nuclear Explosions Using Waveform Cross Correlation. Part II: Synthetic Master Events. Seismic Instruments, 52(3): 207-223. https://doi.org/10.3103/S0747923916030038

[11]

Cohen, E., Datar, M., Fujiwara, S., et al., 2001. Finding Interesting Associations without Support Pruning. IEEE Transactions on Knowledge and Data Engineering, 13(1): 64-78. https://doi.org/10.1109/69.908981

[12]

Gibbons, S. J., Ringdal, F., 2006. The Detection of Low Magnitude Seismic Events Using Array⁃Based Waveform Correlation. Geophysical Journal International, 165(1): 149-166. https://doi.org/10.1111/j.1365⁃246X.2006.02865.x

[13]

Gibbons, S. J., Ringdal, F., 2012. Seismic Monitoring of the North Korea Nuclear Test Site Using a Multichannel Correlation Detector. IEEE Transactions on Geoscience and Remote Sensing, 50(5): 1897-1909. https://doi.org/10.1109/TGRS.2011.2170429

[14]

Gionis, A., Indyk, P., Motwani, R., 1999. Similarity Search in High Dimensions via Hashing. In: Proceedings of the 25th International Conference on Very Large Data Bases (VLDB’99). Morgan Kaufmann Publishers Inc., San Francisco, 518-529.

[15]

He, X., Zhao, L. F., Xie, X. B., et al., 2018. High⁃ Precision Relocation and Event Discrimination for the 3 September 2017 Underground Nuclear Explosion and Subsequent Seismic Events at the North Korean Test Site. Seismological Research Letters, 89(6): 2042-2048. https://doi.org/10.1785/0220180164

[16]

Kafka, A. L., 1990. Rg as a Depth Discriminant for Earthquakes and Explosions: A Case Study in New England. Bulletin of the Seismological Society of America, 80(2): 373-394. https://doi.org/10.1785/bssa0800020373

[17]

Le Bras, R., Arora, N., Kushida, N., et al., 2021. NET⁃VISA from Cradle to Adulthood. A Machine⁃Learning Tool for Seismo⁃Acoustic Automatic Association. Pure and Applied Geophysics, 178(7): 2437-2458. https://doi.org/10.1007/s00024⁃020⁃02508⁃x

[18]

Li, J., 2021. Key Technologies for Intelligent Processing of Seismic Data from Global Sparse Networks (Dissertation). Beijing University of Posts and Telecommunications, Beijing (in Chinese with English abstract).

[19]

Li, Z. F., Zhan, Z. W., 2018. Pushing the Limit of Earthquake Detection with Distributed Acoustic Sensing and Template Matching: A Case Study at the Brady Geothermal Field. Geophysical Journal International, 215(3): 1583-1593. https://doi.org/10.1093/gji/ggy359

[20]

Liu, J. Q., Liu, C., Zhang, Y., et al., 2019. The Moment Tensor of ML3.4 Earthquake on September 23, 2017, Democratic People’s Republic of Korea (DPRK). Chinese Journal of Geophysics, 62(2): 534-543 (in Chinese with English abstract).

[21]

Mallat, S., 2009. Wavelet Bases. In: Mallat, S., ed., A Wavelet Tour of Signal Processing (Third Edition). Academic Press, Boston, 263-376. https://doi.org/10.1016/B978⁃0⁃12⁃466606⁃1.X5000⁃4

[22]

Muir, J., Fernando, B., Barrett, E., 2023. False Positives Are Common in Single⁃Station Template Matching. Seismica, 2(2). https://doi.org/10.26443/seismica.v2i2.385

[23]

Richards, P. G., Kim, W. Y., 2005. Equivalent Volume Sources for Explosions at Depth: Theory and Observations. Bulletin of the Seismological Society of America, 95(2), 401-407. https://doi.org/10.1785/0120040034

[24]

Rong, K. X., Yoon, C. E., Bergen, K. J., et al., 2018. Locality⁃Sensitive Hashing for Earthquake Detection: A Case Study of Scaling Data⁃Driven Science. Proceedings of the VLDB Endowment,11(11): 1674-1687.https://doi.org/10.14778/3236187.3236214

[25]

Senobari, N. S., Funning, G. J., Keogh, E., et al., 2018. Super⁃Efficient Cross⁃Correlation (Sec⁃C): A Fast Matched Filtering Code Suitable for Desktop Computers. Seismological Research Letters, 90 (1): 322-334. https://doi.org/10.1785/0220180122

[26]

Sugita, D., Ito, Y., 2015. Study of Web Service Estimation Using only Ethernet Frame Header or IP Packet Header. In: Proceedings of the 2015 IEEE 4th Global Conference on Consumer Electronics (GCCE).IEEE,Osaka. https://doi.org/10.1109/GCCE.2015.7398648

[27]

Wang, J., Liu, J. M., Li, J., et al., 2021. Progress and Testing of NET⁃VISA Software. Progress in Earthquake Sciences, 51(6): 274-278 (in Chinese with English abstract).

[28]

Yao, J. Y., Tian, D. D., Lu, Z., et al., 2018. Triggered Seismicity after North Korea’s 3 September 2017 Nuclear Test. Seismological Research Letters, 89(6): 2085-2093. https://doi.org/10.1785/0220180135

[29]

Yoon, C. E., O’Reilly, O., Bergen, K. J., et al., 2015. Earthquake Detection through Computationally Efficient Similarity Search. Science Advances, 1(11): e1501057. https://doi.org/10.1126/sciadv.1501057

[30]

Zhao, D. P., Liu, X. Q., Li, H., et al., 2012. Detection of Regional Seismic Events by Kurtosis Method and Automatic Identification of Direct P⁃Wave First Motion by Kurtosis⁃AIC Method. Journal of Seismological Research, 35(2): 220-225, 295 (in Chinese with English abstract).

[31]

Zhao, L. F., Xie, X. B., He, X., et al., 2017. Seismological Discrimination and Yield Estimation of the 3 September 2017 Democratic People’s Republic of Korea(DPRK) Underground Nuclear Test. Chinese Science Bulletin, 62(35): 4163-4168 (in Chinese).

AI Summary AI Mindmap
PDF (7919KB)

0

访问

0

被引

详细

导航
相关文章

AI思维导图

/