CN114707558A - 冰崩次声特征提取及分类识别方法和介质 - Google Patents

冰崩次声特征提取及分类识别方法和介质 Download PDF

Info

Publication number
CN114707558A
CN114707558A CN202210423964.1A CN202210423964A CN114707558A CN 114707558 A CN114707558 A CN 114707558A CN 202210423964 A CN202210423964 A CN 202210423964A CN 114707558 A CN114707558 A CN 114707558A
Authority
CN
China
Prior art keywords
signal
infrasound
wavelet
entropy
scale
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202210423964.1A
Other languages
English (en)
Inventor
陈乔
韦方强
张议芳
刘敦龙
陈吉龙
刘阳
郭小载
贾东铭
王莉莎
马健荣
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Chongqing Institute of Green and Intelligent Technology of CAS
Original Assignee
Chongqing Institute of Green and Intelligent Technology of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Chongqing Institute of Green and Intelligent Technology of CAS filed Critical Chongqing Institute of Green and Intelligent Technology of CAS
Priority to CN202210423964.1A priority Critical patent/CN114707558A/zh
Publication of CN114707558A publication Critical patent/CN114707558A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/02—Preprocessing
    • G06F2218/04—Denoising
    • G06F2218/06—Denoising by applying a scale-space analysis, e.g. using wavelet analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G06F18/2411—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on the proximity to a decision surface, e.g. support vector machines
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N20/00—Machine learning
    • G06N20/10—Machine learning using kernel methods, e.g. support vector machines [SVM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/08—Learning methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08—Feature extraction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/12—Classification; Matching
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Artificial Intelligence (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Software Systems (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Mathematical Physics (AREA)
  • Computing Systems (AREA)
  • Computational Linguistics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Health & Medical Sciences (AREA)
  • Biomedical Technology (AREA)
  • Biophysics (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Molecular Biology (AREA)
  • Medical Informatics (AREA)
  • Signal Processing (AREA)
  • Complex Calculations (AREA)

Abstract

本发明公开了一种冰崩次声特征提取及分类识别方法和介质，冰崩次声特征提取方法包括：对原始信号进行清洗，得到次声信号；对次声信号的时域、频域进行分析，得到声压范围、峭度、偏度以及功率谱；基于所述次声信号，通过离散小波变换确定次声信号在各尺度上的能量分布，并进行特征值的提取，得到小波能量熵、尺度熵、奇异熵三类特征；基于所述次声信号，通过希尔伯特黄变换得到希尔伯特边际谱。本发明建立了原始次声数据库、视频数据库以及次声样本特征数据库，并可以基于次声样本特征数据库训练分类模型以识别冰崩事件的发生。

Description

冰崩次声特征提取及分类识别方法和介质
技术领域
本发明涉及次声处理技术领域，特别是涉及一种冰崩次声特征提取及分类识别方法和介质。
背景技术
冰崩是重力作用下冰川前端发生的冰体瞬间崩落事件，多发生在冰川末端，有时也会发生在冰川的中上部。近年来，受人为活动和气候变化等因素的影响，全球气候变暖趋势加剧，冰冻圈因此受到严重影响，尤其是青藏高原作为全球气候变化的敏感区，出现大面积的冰川消融退缩现象，导致冰川的失稳现象逐渐增加，同时受地形、地震、降水等内外作用力的影响，使得冰崩灾害的发生概率逐渐增加。
由于冰崩大多发生在偏远地区，并没有造成严重的经济损失和人员伤亡，加之受监测技术的限制，早期并没有受到广泛的关注和研究。而近年来频繁的冰崩灾害以及其形成的灾害链严重威胁当地居民的生命财产安全，同时也给高原地区的生态环境带来极大的隐患，加之随着遥感技术的发展，为冰崩灾害提供有力的监测数据，因而使得对冰崩灾害的研究逐渐增加。目前对于冰崩灾害的研究主要是基于前期现场勘察收集到的资料、借助遥感影像、地震仪以及进行数值模拟。由于冰崩灾害常发生在偏远山区，地势险恶，现场勘察困难且危险，因而实际收集到的资料及其有限；遥感技术的发展促进了对冰崩灾害研究的进展，遥感影像可以对冰川进行全方位的监测，但遥感影像依赖卫星的时空分辨率，同时无法进行实时监测；而地震仪目前主要用于记录冰崩的频次，并不对冰崩的整个动态过程进行记录；数值模拟可以对冰崩的主要影响因素进行过程模拟，而对其机理的模拟目前还存在一定的难度。
发明内容
本发明提供一种冰崩次声特征提取及分类识别方法和介质，拟开展冰崩次声监测，以构建次声特征数据库，为基于次声识别冰崩事件提供技术铺垫。
为了解决上述技术问题，本发明采用以下技术方案：
根据本发明的第一方面，提供一种冰崩次声特征提取方法，所述方法包括：对原始信号进行清洗，得到次声信号；对次声信号的时域、频域进行分析，得到声压范围、峭度、偏度以及功率谱；基于所述次声信号，通过离散小波变换确定次声信号在各尺度上的能量分布，并进行特征值的提取，得到小波能量熵、尺度熵、奇异熵三类特征；基于所述次声信号，通过希尔伯特黄变换得到希尔伯特边际谱。
根据本发明的第二方面，提供一种冰崩次声分类识别方法，其特征在于，基于本发明任一实施例所述的冰崩次声特征提取方法，所述方法包括：将提取到的声压范围、峭度、偏度、功率谱、小波能量熵、尺度熵、奇异熵和/或希尔伯特边际谱组成特征向量，并根据所述特征向量判断冰崩事件的发生。
根据本发明的第三方面，提供一种计算机可读存储介质，其上存储有计算机可读指令，当所述计算机可读指令被计算机的处理器执行时，使计算机执行本发明任一实施例所述的方法
与现有技术相比，本发明的有益效果在于：
本发明采用次声监测的方法对冰崩产生的次声进行直接监测，通过对采集的次声信号进行事件样本提取、滤波降噪，采用时-频分析、小波及小波包分析以及希尔伯特黄变换来提取各事件的次声特征，并对不同事件的特征进行分析，随后采用支持向量机和人工神经网络的方法对不同事件进行识别与分类，并将分类识别模型导出用于预测未知事件的类别，以实现对冰崩事件进行实时准确地预测，保护当地居民的人身财产安全。
附图说明
图1是本发明的一种冰崩次声特征提取方法的流程图。
图2a是根据本发明实施例的原始信号图。
图2b是根据本发明实施例的去直流偏移信号图。
图3a是根据本发明实施例的凯泽窗的幅值响应图。
图3b是根据本发明实施例的凯泽窗的相位响应图。
图4是根据本发明实施例的凯泽窗时域、频域响应图。
图5是根据本发明实施例的sym4小波基示意图。
图6是根据本发明实施例的sym2、sym4、sym8小波基时、频域对比图。
图7是根据本发明实施例的小波硬阈值去噪效果图。
图8是根据本发明实施例的小波软阈值去噪效果图。
图9是根据本发明实施例的时频分析图。
图10是根据本发明实施例的小波分解5层示意图。
图11是根据本发明实施例的小波包分解3层示意图。
图12是根据本发明实施例的小波分解每一层的波形图。
图13是根据本发明实施例的每层的能量占比图。
图14是根据本发明实施例的小波包每个节点的能量占比图。
图15是根据本发明实施例的奇异值分解占比图。
图16是根据本发明实施例的希尔伯特谱图。
图17是根据本发明实施例的希尔伯特边际谱图。
图18是根据本发明实施例的0～7Hz希尔伯特边际谱图。
图19是根据本发明实施例的冰崩事件21维特征向量示意图。
图20是根据本发明实施例的背景事件21维特征向量示意图。
图21是根据本发明实施例的2021-8-1-08:24:00冰崩事件的信号图。
图22是根据本发明实施例的2021-6-10-12:56:00冰崩事件的信号图。
图23是根据本发明实施例的2021-8-1-10:23:00冰崩事件的信号图。
图24是根据本发明实施例的一种冰崩次声分类识别方法的技术路线图。
图25是根据本发明实施例的次声信号的分类识别系统图。
具体实施方式
现在结合说明书附图对本发明做进一步的说明。
本发明实施例提供一种冰崩次声特征提取方法。如图1所示，所述方法包括：
步骤S100，对原始信号进行清洗，得到次声信号。采集的原始信号通常会存在大量的无用信息及高频噪音，因此在进行数据分析前需要对原始信号进行清洗，即通过去直流、滤波、降噪等处理，以提高信号的信噪比，去除无用的信号成分。
由于受采集器设备内部各电子元器件之间相互关系的影响，系统会存在直流漂移，或直流分量。直流分量即信号的平均值，或频率为零的信号分量，其特点是信号的幅值不随时间变化而发生改变，因此这部分信号不携带信息，但其携带能量，会对后续的数据处理带来能量误差，因此在进行信号处理之前应把这部分信号去除。如图2a和图2b所示，原始信号存在约0.2帕斯卡的直流漂移，因此通过yDC＝y(x)-mean[y(x)]即可去除直流漂移的信号。
实际的信号采集中，为获取更为准确的次声信号，传感器设置的采样频率为100Hz，即采集的声波信号存在大于20Hz的可闻声段，由于分析时仅考虑次声波段，因而在处理信号之前，需设置低通滤波器将大于20Hz的信号滤除掉。
合适的滤波器，能尽可能减少信号失真。滤波器按照信号处理的形式可以分为模拟滤波器和数字滤波器，在本发明实施例中选用的是数字滤波器，即通过数字滤波器软件对信号进行滤波处理；按照通频带可以分为高通滤波器、低通滤波器、带通滤波器、全通滤波器和带阻滤波器，由于是对次声波进行分析，因而选用的是低通滤波器，即低于20Hz频带的声波可以通过，而高于20Hz的频段被阻拦；按照滤波网络结构(运算结构)分类可分为无限长单位冲激响应滤波器(Infinite impulse response，IIR)和有限长单位冲激响应滤波器(Finite impulse response，FIR)，本发明中选用的是FIR型滤波器。FIR滤波器较IIR 滤波器具有严格的线性相位，能避免在信号处理过程中产生信号的相位失真，同时因为其单位冲激响应为有限长，没有递归结构，其系统更稳定。
滤波器采用设计相对简单，应用广泛的时域设计法，即窗函数设计法，选择适应性较强的凯泽窗(Kaiser)，其窗函数的形式如式(1)所示：
Figure BDA0003607736440000051
其中，N为凯泽窗的长度，n为滤波器阶数，I0为第一类零阶变型贝塞尔函数，RN(n)为矩形窗函数，β是调节窗形状的参数，可以同时调整主瓣宽度和旁瓣衰减，β越大，ω(n)窗越窄，窗谱的旁瓣衰减越大，但主瓣宽度有相应的增加，因此并不是β越大越好，通常β取值在4～9之间，如表1，主瓣宽度主要影响信号能量分布和频率分辨能力，即主瓣越宽，在同等情况下其频率分辨率越低，而旁瓣衰减率反映能量泄漏程度(频谱拖尾效应)，衰减越慢，说明能量泄漏、频谱拖尾越严重。
表1.凯泽窗参数对滤波器性能的影响
Figure BDA0003607736440000052
通过MATLAB2020b内置Filter Designer进行低通滤波器设计，得到相应的滤波器参数。示例性的，当选择FIR低通滤波器时，设置滤波器阶数(即谐波过滤的次数)为200，阶数越高，性能越好，但计算时间会延长。选择凯泽窗，对应β值取9，原始信号采样频率为100Hz，设置截止频率为20Hz。
由图3a和图3b可以看出，凯泽窗的幅值响应和相位响应都在一定程度上逼近理想线性相位滤波器，即幅度响应为矩形，相位响应呈线性。如图4所示，其在时域响应上，其泄露因子为零，在频域响应上，其主瓣宽度为0.017，相对旁瓣衰减为66.2dB，即信号能量主要集中于主瓣，越接近真实的频谱。因此在利用本发明实施例设计的滤波器能够有效滤除无用的信号频段。
为进一步提高信号的信噪比，提高所需信号的权重，对过滤的信号进一步进行降噪处理，以去除部分的高频部分。本发明实施例选用的是小波软阈值降噪法，调用MATLAB2020b内置小波降噪工具箱Wavelet Signal Denoiser进行小波降噪处理。降噪的关键是小波基、阈值、阈值函数的选择，本发明实施例中选取sym4小波基进行5层小波分解，采用无偏风险估计(SURE)软阈值降噪。
symlets小波具有很好的正交性、正则性、紧支撑性和近似对称性，其支撑长度为2N-1，消失矩为N，本研究中选用的sym4小波基，如图5，其支撑长度为7，消失矩为4，通常支撑长度在5～9之间较为合适，支撑长度过长，消失矩太高，会导致小波函数过于震荡，边界能量过于集中，支撑长度过短，消失矩太低，则又不利于能量集中，如图6所示。
常用的阈值有无偏风险估计阈值(SURE)、极大极小阈值(Minimaxi)、通用阈值(Universal threshold)、贝叶斯阈值(Bayes)等。由于在降噪之前已进行滤波处理，滤除大部分的高频噪音，因此本发明实施例中采用的是无偏风险估计阈值(SURE)，其去噪相对保守，适用于噪声在高频段分布较少的情况下，减少对有用信号的去除。
阈值函数分为硬阈值去噪法和软阈值去噪。如图7和图8所示，硬阈值去噪即给定一个阈值，当小波系数的绝对值小于阈值，则令其为零，大于阈值时则保持不变，软阈值去噪则当其大于阈值时，减去阈值。由此可知，硬阈值去噪会使去噪后的小波系数出现阶跃，会使重构的信号产生一定的震荡，但其更为准确的保留原始信号，而软阈值去噪相对能减少阶跃的出现，重构信号更为光滑，但对原始信号造成一定的压缩，存在一定的信号失真，因此本发明实施例中选用硬阈值去噪法。
在对原始信号进行处理后，即进行特征提取。特征提取部分是本发明的核心内容。特征提取就是从大量复杂的数据中提取出有效信息的过程，特征的提取要满足两个特点，一是同一类别的信息具有高度的相似性，二是不同类别的信息之间具有很大的差异性。特征提取在信号处理领域主要分为时域分析、频域分析、时频联合域分析以及其他一些非线性分析。通过对次声信号的特征提取，可以构建相应的次声特征数据库，为冰崩的次声研究提供数据基础。
因此，在步骤S200，对次声信号的时域、频域进行分析，得到声压范围、峭度、偏度以及功率谱。
信号的时域分析主要用于提取信号的声压、波形、事件持续的时间等特征信息，本研究中主要提取信号的声压范围、峭度和偏度，声压范围反应信号在时域上的振幅大小；峭度反应波形的平缓程度，峭度越高，波形越陡；偏度反应波形的偏斜程度，偏度为正值，表示右偏态，偏度为负值，表示左偏态。
声压范围表示：
Pk＝max(x)-min(x) 式(2)
峭度表示：
Figure BDA0003607736440000071
偏度表示：
Figure BDA0003607736440000072
式(2)-(4)中，x表示输入信号，μ表示输入信号的均值，σ表示输入信号的标准差，E表示输入信号的期望值。
信号的频域分析主要是通过傅里叶变换(FT)得到信号的频谱，再通过傅里叶变换的平方与信号长度的比值得到信号的功率谱等特征信息。频谱反应的是信号在频域上的分布情况，功率谱即功率谱密度函数，表征单位频率内信号功率分布。本发明实施例中主要采用计算更为便捷的快速傅里叶变换(FFT) 的方法求取信号的功率谱
Figure BDA0003607736440000081
并提取信号的功率谱范围。
Figure BDA0003607736440000082
式中，X(f)表示信号的快速傅里叶变换，N表示信号的长度。
根据上述提取到的相关信号特征，可以得到如图9所示的时频分析图。
在步骤S300，基于所述次声信号，通过离散小波变换确定次声信号在各尺度上的能量分布，并进行特征值的提取，得到小波能量熵、尺度熵、奇异熵三类特征。
小波分析属于时频分析的一种，其可以对时间窗和频率窗均进行改变，且都能很好的表征信号的局部特征。在具体分析时，时间窗和频率窗都可以根据需要作出调整。在对信号低频部分进行分析时，小波变换可以牺牲时间分辨率来换取较高的频率分辨率：在对信号高频部分进行分析时，小波变换又可以牺牲频率分辨率以提高时间分辨率。本发明实施例中，主要采用离散小波(包) 变换来进行信号的特征提取。小波包变换原理同小波变换，小波变换仅对每一层的近似部分进行再分解，小波包变换同时对每一层的近似部分和高频部分进行分解，以判断是否在高频部分也存在信号的分量。图10为小波分解示意图。图11为小波包分解示意图。
本发明实施例中，通过如下步骤来确定次声信号在各尺度上的能量分布：
第一步：假设函数ψ(t)满足ψ(t)∈L2(R)，即平方可积，同时其傅里叶变换
Figure BDA0003607736440000091
满足如式(6)所示条件，则ψ(t)为一个母小波。
Figure BDA0003607736440000092
第二步：对母小波进行伸缩平移得到一组小波基函数ψa,b(t)，如下式(7)，其中a,b∈R,a>0，a为尺度因子，b为平移因子，均为连续变化的值。
Figure BDA0003607736440000093
第三步：对尺度参数a和平移参数b进行离散化，即a＝2j，bj,k＝2jk， j,k∈Z，求得离散小波基函数ψj,k(t)。
Figure BDA0003607736440000094
第四步：对原始数据χ(t)进行预处理得到预处理后的次声信号S(t)。该步骤是对传感器采集到的相关数据进行处理，可以参照步骤S100所述的方法进行处理。若已执行步骤S100，则可以直接将步骤S100处理后的次声信号作为次声信号S(t)。
第五步：对次声信号S(t)进行离散小波变换。
Figure BDA0003607736440000095
第六步：计算不同尺度上的能量及能量占比。
仅作为示例，采用小波变换对原始信号进行5层分解，小波包变换进行3 层分解。如图12为小波分解每一层的波形，图13为每层的能量占比，1、2、3、 4、5、6对应的频段为0～0.79Hz、0.79～1.57Hz、1.57～3.13Hz、3.13～6.25Hz、 6.25Hz～12.5Hz、12.5～25Hz。图14为小波包每个节点的能量占比，图15为奇异值分解占比，1、2、3、4、5、6、7、8对应的频段为0～3.13Hz、3.13～6.25Hz、6.25～9.38Hz、9.38～12.5Hz、12.5～15.63Hz、15.63～18.75Hz、18.75～21.88Hz、 21.88～25Hz。
在通过上述步骤获得小波分解各尺度的能量分布情况后，将进行特征值的提取，本发明实施例中主要提取小波能量熵、尺度熵、奇异熵三类特征。
在已经确定小波分解后每一层的能量Ei及能量占比Pi，则能量熵为
Figure BDA0003607736440000101
其中，
Figure BDA0003607736440000102
同理，小波尺度熵即对小波包分解的所有节点小波系数Cn的香农熵
Figure BDA0003607736440000103
其中，
Figure BDA0003607736440000104
在求解奇异熵前需要对信号进行奇异谱分解，它根据信号的时间序列进行分解、重构，提取出信号中不同的成分。
X(t)＝Um×m∑m×nVn×n T 式(12)
X(t)表示原始信号构成的矩阵。本研究中对小波包分解的每一个节点求奇异值Sj，再求信号的奇异熵：
Figure BDA0003607736440000105
其中，
Figure BDA0003607736440000106
从小波分解和小波包分解的能量占比图以及奇异值占比图中可以看出，信号在低频段分布集中，且分布差异较大，因此在选取特征时，不仅提取小波的总能量熵、总尺度熵、总奇异熵，还对小波包分解的前两个节点和前四个节点分别取1/4能量熵、1/4尺度熵、1/4奇异熵和1/2能量熵、1/2尺度熵、1/2奇异熵作为信号的特征。
最后在步骤S400，基于所述次声信号，通过希尔伯特黄变换得到希尔伯特边际谱。信号经过EMD后被分解为若干单一频率的本征模函数(Intrinsic Mode Function,IMF),然后对每个IMF分量作Hilbert变换得其相应的瞬时频率，这些瞬时频率反映了原始信号的时频能量分布，即Hilbert谱，谱能精确反应信号的分布规律，包括信号的能量在时间和频率上的分布规律。
其具体包括如下步骤：
第一步：对原始数据χ(t)进行预处理得到预处理后的次声信号S(t)；
第二步：寻找次声信号S(t)的所有极大值点、极小值点和端点；
第三步：通过三次样条函数对所有标记的极值点进行插值计算，然后求得次声信号S(t)的上下包络序列Eup、Edown，并计算上下包络的均值M1(t)；
M1(t)＝(Eup+Edown)/2 式(14)
第四步：计算次声信号S(t)的IMF信号I(t)；
I1(t)＝S(t)-M1(t) 式(15)
第五步：判断I(t)是否满足IMF的两个条件，即信号波形的总的跨零点个数与全部极值点的个数相减的差值小于等于1，以及信号的上下包络的局部均值为零。若满足以上两个条件，则IMF1＝I1(t)，得到第一个IMF分量；若不满足条件，则将I1(t)当做原始信号重复上述第二到第四步的步骤，直到满足IMF 的两个条件，得到IMF1；
第六步：将第一个IMF分量去掉，得到剩余分量r1；
r1＝S(t)-IMF1 式(16)
第七步：判断r1是否符合EMD的三个条件，即一是信号的极值点个数大于等于2，即至少有一个极大值一个极小值点；二是极值点的时间尺度可以并且唯一确定原始信号的局部时域特性；三是如果一个信号没有极值点，则必须要有通过数次微分可以求得极值点的奇异点或者拐点。若满足上述三个条件，则 EMD分解结束，若不满足，则以r1为新的原始信号重复上述第二道第六步的步骤，得到其他IMFn分量；
第八步：当rn满足终止条件时，EMD分解结束，原始信号S(t)可以表示为多个IMF分量的和的形式；
Figure BDA0003607736440000121
第九步：对每一个IMF分量进行希尔伯特黄变换，得到希尔伯特谱H(w,t)，其中Re表示取IMF分量的实部；
Figure BDA0003607736440000122
第十步：对希尔伯特谱H(w,t)进行积分，得到边际谱h(w)；
Figure BDA0003607736440000123
基于上述步骤可以得到希尔伯特谱和希尔伯特边际谱，如图16和图17所示，本发明主要以希尔伯特边际谱作为不同样本信号的特征。希尔伯特边际谱可以反应出信号在较低频的局部特征，因此以信号在0～7Hz范围的边际谱作为信号的特征值。即如图18所示。
基于上述时域、频域分析，小波(包)变换，希尔伯特黄变换，目前共提取特征141维，其中时域频域特征共8维，小波(包)变换提取特征21维，希尔伯特黄变换提取特征112维。图19和图20为冰崩事件和背景事件的基于小波(包)变换提取的特征构成的特征向量。从图19和图20可以看出，冰崩事件和背景事件的特征向量存在明显不同，基于小波(包)变换提取的特征构成的特征向量可以作为判断是否为冰崩事件的特征向量。
在一些实施例中，本发明实施例通过野外监测点的设计与建设研究来获取原始信号。野外监测要考虑到设备对环境的适应性，背景噪声的干扰，待监测事件的发生概率及对下游的危害，数据传输的实时性，以及设备运输安装的便捷性等。综合上述要求，对设备进行环境适应性测试，确保设备能适应-30℃的环境温度。根据设备安装需要考虑太阳能补给、避雷、防水等问题，设计的安装立杆要保证具有上述的各模块或功能同时还要考虑立杆的高度和力的平衡性，最后对设备的传输方式进行改进，确保采集的次声数据能在线传输。
本发明实施例的数据采集来源于青藏高原地区。在数据采集前期的前期工作包括电源环境适应性测试、设备安装立杆设计、次声采集器在线传输测试以及设备安装选址。
电源环境适应性测试：青藏高原地区环境恶劣，冬季温度可达-30℃，对设备的环境适应性及其考验，因此在开展实际的设备安装之前，先对设备进行测试，将电源放入-30℃的冰箱中冷冻一周后再测试电源是否能正常工作，测试结果表明电源设备具有很好的抗冻性，能够适应青藏高原恶劣的环境。
设备安装立杆设计：开展设备安装立杆的设计，需要综合考量环境对立杆的影响，因此不能设计过高、也不能太低，电源线路不能暴露在空气中，不同的部件箍在立杆上要保证整个圆周受力稳定等因素，确定设计方案之后前往加工厂加工并组装测试，确保立杆的稳定性。
次声采集器在线传输测试：由于数据需要实时传输，因此对选用的次声采集器进行在线传输测试，确保其能在网络稳定时在线传输，网络异常时能本地保存。
设备安装选址：设备安装受多因素的影响，(1)由于选用的是次声监测的方式，因而要考虑到尽量少的背景噪音的干扰，嘉龙错上端冰川位于喜马拉雅山脉中段北坡，属于大陆性冰川，受印度洋季风气候影响较小，局部气候相对稳定，减少对次声信号的干扰；(2)2002年，嘉龙错发生过两次冰崩事件，导致冰崩体入湖引发冰湖溃决，形成泥石流、洪水，对下游的聂拉木县城甚至中尼友谊桥均造成一定的影响；(3)由于设置的次声采样频率为100Hz，数量较大，北斗通讯无法满足大量数据传输的要求，而嘉龙错附近网络覆盖，有利于数据的在线传输；(4)由于设备安装涉及到立杆、太阳能板、电源等较重的设备，嘉龙错可以驱车直达，减少人力畜力的运输成本。综合考量，最终选址位于西藏自治区日喀则市聂拉木县嘉龙错侧碛垄一平缓处，嘉龙错北纬28°12′ 24.55″～28°12′54.13″，东经85°51′00.63″～85°51′11.48″，海拔约4334 米。
在通过上述前期准备工作后，即采集到相应的声波信号数据集。为了保证不漏采次声信号，通常采样频率会设置在100Hz以下，并且数据实时采集，这将会造成数据量比较大，传统的人工识别进行样本提取会耗费大量的时间，因此本研究中采用Short timeaverage(STA)与Long time average(LTA)的比值，即STA/LTA时窗比值法来获取事件到达的初始时间，根据获取的初始时间结合视频监控判断是否为事件，并进行事件样本提取。
STA/LTA时窗比值法最早用于地震事件的识别，通过设置一个长时窗和一个短时窗，窗随时间序列移动，在移动的同时对长时窗和短时窗内的信号幅值取平均值，当两者的比值超过设定的阈值时，判定为事件发生。
Figure BDA0003607736440000141
式中Wlta、Wsta分别为长时窗和短时窗的长度，A(i)为信号的幅值。长时窗用于衡量环境噪音的强度，短时窗则体现为对事件信号的敏感度，因此长时窗不能设置过长，短时窗不能设置过短，阈值的大小也反映对事件的敏感度，阈值越低，对事件的识别越敏感，会存在事件的误识别，因此在设置阈值时，应结合周围环境。本研究长时窗和短时窗设置的长度为30s、0.5s，阈值设定为 2。如图21-23为识别到的三次冰崩事件。
本发明实施例还提供一种冰崩次声分类识别方法，请参考图24所示，该方法基于如上所述的方法构建的次声特征数据库，可以用于判断冰崩事件的发生。所述方法包括：基于时、频域分析、小波变换、希尔伯特黄变换提取的特征构成特征向量，构建分类模型，通过不断训练改进模型参数，达到最佳的分类效果；将训练得到的分类模型用于未知信号的识别，对信号进行类别判定，输出分类结果。
本发明实施例所采用的分类识别即有监督的模式识别，即在已知分类结果的情况下对信号进行分类训练，通过不断的学习测试过程，实现对未知信号的分类识别。一个完整的分类识别系统是由训练和识别两部分构成的，如图25，通过对样本数据进行训练构建识别模型，并不断对模型进行改进，使分类效果达到最佳。常见的分类识别算法有人工神经网络、支持向量机、贝叶斯分类器、 K-近邻法、决策树等。
本发明实施例采用有监督的模式识别方式对次声信号进行分类，选用的是支持向量机(Support Vector Machine，SVM)和人工神经网络(Artificial Neural Network，ANN)分类识别算法，将次声信号分为训练样本和测试样本，对训练样本进行预处理，基于时、频域分析、小波变换、希尔伯特黄变换提取的特征构成特征向量，构建分类模型，通过不断训练改进模型参数，达到最佳的分类效果；将训练得到的分类模型用于未知信号的识别，对信号进行类别判定，输出分类结果。
支持向量机的分类方式是将向量映射到高维空间，在高维空间中找到一个超平面，使得该平面两侧距平面最近的样本之间的距离最大。该算法运用到核函数的思想，因而可以解决非线性问题的分类，同时其适合小样本分类，即适合于对灾害类型进行分类。
人工神经网络通过数学、物理以及信息处理方法对人脑的神经网络进行模拟，并建立简化的模型，由输入层、隐藏层、输出层组成，其中隐藏层可以有若干层，每一层由大量简单处理元件连接构成高度并行的非线性系统，可实现信息的并行处理，同时其高度的非线性系统使得整个网络具有良好的容错性，即网络中部分神经元损坏不会对整个网络造成影响，能根据自身系统的性能来适应环境的变化。
特征提取结论基于79组冰崩样本及71组背景噪声样本提取的时域、频域、小波包变换的部分特征进行对比。分类识别结论基于70组冰崩样本及70组背景噪声样本提取的小波包特征，包括小波能量熵、小波尺度熵、小波奇异熵9 组，以及第三层前四个节点功率谱特征共12组，即共计21组特征构成的特征向量进行分类识别。最终，得出以下结论：
1)冰崩事件声压范围均值为0.0042帕斯卡，背景事件的声压范围为0.00074 帕斯卡，冰崩事件的声压范围约为背景事件的5.66倍，即冰崩声压较背景声压广；而冰崩事件的功率谱范围均值58.12dB，背景事件的功率谱范围均值为69.54 dB，冰崩功率谱范围较背景功率谱范围窄。
2)冰崩事件的峭度、偏度均值为6.52和﹣0.14，背景事件的峭度、偏度均值为8.59和﹣0.53，即背景事件具有更陡、左侧拖尾更长的波形特征。
3)基于小波能量分布可知，冰崩次声频段主要分布在0～7Hz，背景次声频段主要分布在0～3Hz，即两种类型的次声频段存在重叠，因此也表明在滤波和降噪处理不能直接去除背景噪声。
4)基于小波能量熵E、尺度熵SE、奇异熵SVD，如表2所示，冰崩事件的均熵值基本高于背景事件，表明冰崩事件所包含的信息的复杂度更大，也再次证明两种类型次声信号存在重叠频段。
表2.基于小波熵的冰崩特征与背景特征对比
Figure BDA0003607736440000161
5)基于支持向量机的分类识别算法的正确识别率为98.6％，其中冰崩事件的正确识别率为98.6％，背景事件的正确识别率为98.6％。基于人工神经网络算法的总正确识别率为99.3％，其中训练集和验证集的正确识别率均为100％，测试集的正确识别率为95.2％，其中冰崩事件的正确识别率为90.9％，背景事件的正确识别率为100％。由特征提取结论可知，目前的两类样本之间差异较大，特征较为明显。
以上描述旨在是说明性的而不是限制性的。例如，上述示例(或其一个或更多方案)可以彼此组合使用。例如本领域普通技术人员在阅读上述描述时可以使用其它实施例。另外，在上述具体实施方式中，各种特征可以被分组在一起以简单化本公开。这不应解释为一种不要求保护的公开的特征对于任一权利要求是必要的意图。相反，本发明的主题可以少于特定的公开的实施例的全部特征。从而，以下权利要求书作为示例或实施例在此并入具体实施方式中，其中每个权利要求独立地作为单独的实施例，并且考虑这些实施例可以以各种组合或排列彼此组合。本发明的范围应参照所附权利要求以及这些权利要求赋权的等同形式的全部范围来确定。

Claims (10)

1.一种冰崩次声特征提取方法，其特征在于，所述方法包括：
对原始信号进行清洗，得到次声信号；
对次声信号的时域、频域进行分析，得到声压范围、峭度、偏度以及功率谱；
基于所述次声信号，通过离散小波变换确定次声信号在各尺度上的能量分布，并进行特征值的提取，得到小波能量熵、尺度熵、奇异熵三类特征；
基于所述次声信号，通过希尔伯特黄变换得到希尔伯特边际谱。
2.根据权利要求1所述的方法，其特征在于，所述对原始信号进行清洗，得到次声信号，包括：
去除原始信号中的直流漂移的信号得到去直流偏移信号；
通过FIR滤波器对所述去直流偏移信号进行滤波；
所述FIR滤波器通过式(1)所示的窗函数进行设计：
Figure FDA0003607736430000011
其中，N为凯泽窗的长度，n为滤波器阶数，I0为第一类零阶变型贝塞尔函数，RN(n)为矩形窗函数，β是调节窗形状的参数，用于同时调整主瓣宽度和旁瓣衰减，β越大，ω(n)窗越窄，窗谱的旁瓣衰减越大，但主瓣宽度有相应的增加，β取值在4～9之间；
通过小波软阈值降噪法对过滤后的信号进一步进行降噪处理，以去除部分的高频部分；所述小波软阈值降噪法的阈值选择为无偏风险估计阈值。
3.根据权利要求1所述的方法，其特征在于，所述对次声信号的时域、频域进行分析，得到声压范围、峭度、偏度以及功率谱，包括：
通过如下式(2)-(4)计算所述声压范围、峭度、偏度以及功率谱：
Pk＝max(x)-min(x) 式(2)
Figure FDA0003607736430000021
Figure FDA0003607736430000022
式(2)-(4)中，Pk为声压范围，Ku为峭度，Sk为偏度，x表示输入信号，μ表示输入信号的均值，σ表示输入信号的标准差，E表示输入信号的期望值。
4.根据权利要求1所述的方法，其特征在于，所述基于所述次声信号，通过离散小波变换确定次声信号在各尺度上的能量分布，并进行特征值的提取，得到小波能量熵、尺度熵、奇异熵三类特征，包括：
假设函数ψ(t)满足ψ(t)∈L2(R)，即平方可积，同时其傅里叶变换
Figure FDA0003607736430000023
满足如式(6)所示条件，则ψ(t)为一个母小波；
Figure FDA0003607736430000024
对母小波进行伸缩平移得到一组如式(7)所示的小波基函数ψa,b(t)；
Figure FDA0003607736430000025
其中，a,b∈R,a>0，a为尺度因子，b为平移因子，均为连续变化的值
对尺度参数a和平移参数b进行离散化，即a＝2j，bj,k＝2jk，j,k∈Z，求得如式(8)所示的离散小波基函数ψj,k(t)：
Figure FDA0003607736430000026
根据式(9)对次声信号S(t)进行离散小波变换；
Figure FDA0003607736430000031
计算不同尺度上的能量及能量占比；
根据各尺度的能量分布情况，进行特征值的提取，得到小波能量熵、尺度熵、奇异熵三类特征。
5.根据权利要求4所述的方法，其特征在于，所述根据各尺度的能量分布情况，进行特征值的提取，得到小波能量熵、尺度熵、奇异熵三类特征，包括：
根据小波分解后每一层的能量Ei及能量占比Pi，通过式(10)计算能量熵：
Figure FDA0003607736430000032
其中，
Figure FDA0003607736430000033
通过式(11)计算小波尺度熵SE(p)：
Figure FDA0003607736430000034
其中，
Figure FDA0003607736430000035
Cn为所有节点小波系数；
根据信号的时间序列进行分解、重构，提取出信号中不同的成分：
X(t)＝Um×m∑m×nVn×n T 式(12)
其中，X(t)表示原始信号构成的矩阵；
对小波包分解的每一个节点求奇异值Sj，并根据式(13)计算信号的奇异熵：
Figure FDA0003607736430000036
其中，
Figure FDA0003607736430000037
6.根据权利要求1所述的方法，其特征在于，基于所述次声信号，通过希尔伯特黄变换得到希尔伯特边际谱之后，所述方法还包括：以次声信号在0～7Hz范围的边际谱作为次声信号的特征值。
7.一种冰崩次声分类识别方法，其特征在于，基于如权利要求1至6任一项所述的冰崩次声特征提取方法，所述方法包括：
将提取到的声压范围、峭度、偏度、功率谱、小波能量熵、尺度熵、奇异熵和/或希尔伯特边际谱组成特征向量，并根据所述特征向量判断冰崩事件的发生。
8.根据权利要求7所述的方法，其特征在于，所述将提取到的声压范围、峭度、偏度、功率谱、小波能量熵、尺度熵、奇异熵和/或希尔伯特边际谱组成特征向量，并根据所述特征向量判断冰崩事件的发生，包括：
根据所述特征向量，构建分类模型，通过实测数据集不断训练改进模型参数；
将训练得到的分类模型用于未知信号的识别，对信号进行类别判定，输出分类结果。
9.根据权利要求8所述的方法，其特征在于，通过如下方式获取实测数据集：
基于在冰川实时采集的信号，通过设置一个长时窗和一个短时窗，所述长时窗和一个短时窗随时间序列移动，在移动的同时对长时窗和短时窗内的信号幅值取平均值，当两者的比值超过设定的阈值时，判定为冰崩事件，当两者的比值未超过设定的阈值时，判定为背景事件。
10.一种计算机可读存储介质，其特征在于，其上存储有计算机可读指令，当所述计算机可读指令被计算机的处理器执行时，使计算机执行权利要求1至6或7至9中的任一项所述的方法。
CN202210423964.1A 2022-04-21 2022-04-21 冰崩次声特征提取及分类识别方法和介质 Pending CN114707558A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210423964.1A CN114707558A (zh) 2022-04-21 2022-04-21 冰崩次声特征提取及分类识别方法和介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210423964.1A CN114707558A (zh) 2022-04-21 2022-04-21 冰崩次声特征提取及分类识别方法和介质

Publications (1)

Publication Number Publication Date
CN114707558A true CN114707558A (zh) 2022-07-05

Family

ID=82174637

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210423964.1A Pending CN114707558A (zh) 2022-04-21 2022-04-21 冰崩次声特征提取及分类识别方法和介质

Country Status (1)

Country Link
CN (1) CN114707558A (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115236741A (zh) * 2022-09-26 2022-10-25 成都理工大学 一种基于地震动信号的高速远程冰岩崩灾害链预警方法
CN115422984A (zh) * 2022-11-04 2022-12-02 北京理工大学 一种基于时间尺度信号分解及熵特征的信号分类方法
CN115824481A (zh) * 2022-10-01 2023-03-21 同济大学 一种基于递归演化的实时索杆力识别方法

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115236741A (zh) * 2022-09-26 2022-10-25 成都理工大学 一种基于地震动信号的高速远程冰岩崩灾害链预警方法
CN115236741B (zh) * 2022-09-26 2022-12-09 成都理工大学 一种基于地震动信号的高速远程冰岩崩灾害链预警方法
CN115824481A (zh) * 2022-10-01 2023-03-21 同济大学 一种基于递归演化的实时索杆力识别方法
CN115422984A (zh) * 2022-11-04 2022-12-02 北京理工大学 一种基于时间尺度信号分解及熵特征的信号分类方法
CN115422984B (zh) * 2022-11-04 2023-01-24 北京理工大学 一种基于时间尺度信号分解及熵特征的信号分类方法

Similar Documents

Publication Publication Date Title
CN114707558A (zh) 冰崩次声特征提取及分类识别方法和介质
Özger et al. Long lead time drought forecasting using a wavelet and fuzzy logic combination model: a case study in Texas
CN106886044B (zh) 一种基于剪切波与Akaike信息准则的微地震初至拾取方法
CN111832506B (zh) 一种基于长时序植被指数的重建植被遥感判别方法
CN110398647B (zh) 变压器状态监测方法
CN108648764B (zh) 基于雨水敲击声识别的雨量测量系统及其测量方法
Kunii et al. Estimating the impact of real observations in regional numerical weather prediction using an ensemble Kalman filter
CN104714925A (zh) 一种基于分数阶傅里叶变换和支持向量机的齿轮传动噪声分析方法
Liu et al. An improved empirical mode decomposition method for vibration signal
Du et al. Study on Optical Fiber Gas-Holdup Meter Signal Denoising Using Improved Threshold Wavelet Transform
CN108061653B (zh) 基于谐波-冲击多普勒调制复合字典的列车轮对轴承轨边声信号分离方法
Hu et al. Performance analysis of a wavelet packet transform applied to concrete ultrasonic detection signals
CN103323853B (zh) 一种基于小波包和双谱的鱼类识别方法及系统
CN107290047A (zh) 一种拟合次声台站风噪声‑风速的方法
Dorian et al. Bi-class classification of humpback whale sound units against complex background noise with Deep Convolution Neural Network
CN114743562B (zh) 一种飞机声纹识别方法、系统、电子设备及存储介质
CN104463325A (zh) 一种极地探冰雷达原始数据噪声抑制方法
CN107328578A (zh) 一种用于列车轴承轨边声学故障检测的声源分离方法
Phusakulkajorn et al. Wavelet-transform based artificial neural network for daily rainfall prediction in Southern Thailand
CN113782054B (zh) 基于智能语音技术的闪电哨声波自动识别方法及系统
Yang et al. Simulation System of Lake Eutrophication Evolution based on RS & GIS Technology—a Case Study in Wuhan East Lake
Roshni et al. Modelling land surface temperature using gamma test coupled wavelet neural network
CN113409819B (zh) 一种基于听觉谱特征提取的直升机声信号识别方法
Anarde Transformation and Morphological Impact of Low-Frequency Waves During Hurricane Attack
CN114626417A (zh) 一种基于声信号特征分析的雨量检测方法和系统

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination