CN110196448A - 一种滑坡次声信号识别方法 - Google Patents

Abstract

本发明提供一种滑坡次声信号识别方法，包括以下步骤：采集多种不同类型的来自滑坡体的土壤，获取滑坡次声信号；通过野外现场采集环境干扰次声信号；分别分析所述滑坡次声信号和干扰次声信号在时、频域上的相关特征以及特征差异，提炼出一系列具有表征性的关键特征；根据所述关键特征，结合KNN分类算法，建立滑坡次声信号智能识别模型，对监测到的信号进行判别。本发明能够排除环境噪声干扰，实现滑坡次声信号的准确识别。

Description

一种滑坡次声信号识别方法

技术领域

本发明涉及地质灾害监测技术领域，具体的说是一种滑坡次声信号识别方法。

背景技术

滑坡体从开始发育到最后滑动过程中，会因坡体变形、断裂和摩擦而释放出具有一定频率的、可在空气中稳定传播且蕴含重要地质信息的次声波。滑坡体在灾变过程中产生的次声波信号在时、频域上有其独特特征，次声波的传播速度与空气密度、压力等均无关，仅与温度有关，且其声发射源是球面波，无方向性，向四面八方进行传播。此外，次声具有穿透力强、受空气和水的粘滞作用弱、在传播过程中不易发生畸变且能量衰减极小等特点，能够为滑坡次声波源远距离监测提供有利条件。因此，利用次声信号的特点，次声传感器可以接收滑坡体在灾变过程中产生的次声信号。因而，可以通过对滑坡体灾变过程中产生的次声波进行监测获得一定的提前量以实现警报。

利用次声监测滑坡灾害的现有技术主要是：通过部署一个或多个次声传感器来监测某个滑坡易发区，从而判断该区域内是否将要发生滑坡灾害。然而，这类技术存在的问题是没有排除环境干扰噪声的影响(自然界很多现象或事件都会产生次声波信号)，致使识别准确率不高。

现有技术《泥石流次声信号筛选方法、发生定位方法、路径监控方法》(刘敦龙，2014)公开了一种利用次声对泥石流次声信号进行筛选、发生地定位和运动路径监控的方法。该方法在泥石流易发区域内选取3个易于接收泥石流次声信号的监测地点，每个点安装1个次声传感器，构成一个三角形阵，阵元间距为1～2Km。根据获取的次声信号的特征，判别接收到的信号是否来自泥石流，若来自泥石流，则计算各阵元接收到同一次声信号的时间差，根据基于声达时间差的声源定位算法估算泥石流的发生位置。

该方法若应用于滑坡监测预警，存在的缺陷是：滑坡次声信号特征与泥石流次声信号特征完全不同，无法利用泥石流次声信号识别方法对滑坡次声信号进行判别。

发明内容

针对现有技术中存在的上述不足之处，本发明要解决的技术问题是提供一种滑坡次声信号识别方法。

本发明为实现上述目的所采用的技术方案是：一种滑坡次声信号识别方法，包括以下步骤：

采集多种不同类型的来自滑坡体的土壤，获取滑坡次声信号；

通过野外现场采集环境干扰次声信号；

分别分析所述滑坡次声信号和干扰次声信号在时、频域上的相关特征以及特征差异，提炼出一系列具有表征性的关键特征；

根据所述关键特征，结合KNN分类算法，建立滑坡次声信号智能识别模型，对监测到的信号进行判别。

对采集到多种不同类型的土壤开展滑坡灾变室内模拟实验，从而获取滑坡次声信号。

通过野外现场采集的环境干扰次声信号包括以下信号的任意组合：大风、雷电、爆破以及飞机和汽车的引擎产生的次声信号。

所述分析所述滑坡次声信号和干扰次声信号在时、频域上的相关特征及特征差异，提炼出具有表征性的关键特征，包括以下步骤：

通过切比雪夫II型滤波器对采集到的滑坡次声信号进行滤波降噪处理，以去除随机噪声的影响，同时获取信号的最大振幅，通过信号振幅的变化来计算短时过零率；

通过海明窗对去除随机噪声的滑坡次声信号进行分帧加窗；

通过快速傅里叶变换和短时傅里叶变换对分帧加窗处理后的滑坡次声信号进行频率分布情况分析，提取出峰值频率；

通过一阶矩二阶矩谱估计分析并提取分帧加窗处理后的滑坡次声信号的中心频率和谱线宽度。

所述关键特征包括短时过零率、峰值频率、中心频率、谱线宽度以及最大有效声压。

滑坡次声信号/干扰次声信号的特征向量是[短时过零率，峰值频率，中心频率，谱线宽度]。

所述根据所述相关特征建立滑坡次声信号智能识别模型，对监测到的信号进行判别，包括以下步骤：

滤波降噪处理：利用切比雪夫II型滤波器对待识别信号S₁进行低通滤波降噪，滤除高频部分，只保留次声部分；

分析有效声压：分析该段待识别信号S₁的有效声压P₁，若P₁＜5Pa，则继续判别，否则视为非滑坡次声信号；

分析特征值：利用快速傅里叶变换以及一阶矩和和二阶矩谱估计方法分析该段待识别信号S₁的短时过零率、峰值频率、中心频率和谱线宽度这四个特征值，分别为g1、f1、m1、w1；

建立特征向量：利用滑坡、风、爆破、雷电、引擎以及S1的短时过零率、峰值频率、中心频率和谱线宽度四个特征值建立各次声源的特征向量；

特征向量归一化：将待识别信号S₁与已知类别声源的短时过零率、峰值频率、中心频率和谱线宽度四个特征值分别进行2-范数归一化处理；

待识别信号分类：将待识别信号S₁作归一化处理后的特征向量分别与滑坡、风、爆破、雷电、引擎作归一化处理后的特征向量进行K-近邻分类，记录与待识别信号S₁距离最近的声源类型，即为S₁的声源类型，从而实现待识别信号S₁的类别判断。

所述各次声源的特征向量包括：

滑坡：[0.01，2.0，3.6，3.3]，风：[0.11，2.7，5.6，3.9]，爆破：[0.06，4.5，8.3，6.4]，雷电：[0.12，2.8，4.8，3.8]，引擎：[0.22，5.2，10.4，6.7]；S1：[g1，g1，m1，w1]。

所述K-近邻分类取K＝1，采用欧式距离。

还包括对收集的滑坡次声信号和环境干扰次声信号进行划分，按照4∶1的比例，随机分为训练集和测试集，对训练集执行分析所述滑坡次声信号和干扰次声信号在时、频域上的相关特征以及特征差异步骤，所述测试集用于验证：利用各类型声源剩余的若干段信号对所述滑坡次声信号智能识别模型进行验证。

本发明具有以下优点及有益效果：

1、能够排除环境噪声干扰，实现滑坡次声信号的准确识别。

2、能够进行实时现场监测与判别。

3、能够为滑坡监测预警提供可靠的信息依据。

附图说明

图1为本发明方法流程图；

图2为本发明的土质滑坡灾变室内模拟试验装置的结构图；

图3为本发明实际案例1中的一段滑坡的次声信号图；其中(a)为时序波形图，(b)是时频联合域分析图，(c)是功率谱密度图；

图4为本发明实际案例1中的同类别特征2-范数归一化图；

图5为本发明实际案例2中的一段飞机的次声信号图；其中(a)为时序波形图，(b)是时频联合域分析图，(c)是功率谱密度图；

图6为本发明实际案例3中的一段雷电的次声信号图；其中(a)为时序波形图，(b)是时频联合域分析图，(c)是功率谱密度图。

具体实施方式

下面结合附图及实施例对本发明做进一步的详细说明。

本发明的整体步骤如图1所示。滑坡按组成物质分为土质滑坡和岩质滑坡两类，本发明的实施例主要介绍土质滑坡，该方法同样适用于岩质滑坡。本发明的实施例如下：

(1)对采集的多种(至少有10种类型的土壤)不同类型的土壤(均来自滑坡体)进行多次滑坡灾变室内模拟试验(每种土壤根据压实程度和颗粒级配不同，至少进行20组试验，每组试验根据信号质量和持续时间，可提取出若干段有效的滑坡次声信号)，大量收集不同土质滑坡在灾变过程中产生的次声信号(至少收集1000组有效样本数据，每组3000个样本数值)。同时，对风、雨、雷电以及人类活动如飞机、引擎和矿山爆破等山区常见环境干扰次声进行大量观测采样，均收集大量的次声信号数据样本(每种干扰噪声至少收集1000组有效样本数据，每组3000个样本数值)；

滑坡灾变室内模拟试验是通过在野外的真实滑坡体上采集的各种类型的土壤，然后利用如图2所示的土质滑坡灾变模拟试验装置(该装置的技术内容记载在专利申请2017103173576中)来进行滑坡实验，主要目的是为了获取不同类型土壤、不同固结度的滑坡体在整个灾变过程中释放出的次声信号，为下面的特征分析与总结提供稳定的数据来源。因为在野外真实滑坡不经常发生，难以获得大量的次声数据，只能通过模拟实验来获取。

本步骤主要是获取目标次声信号(滑坡次声信号)和环境干扰次声信号。因为在滑坡监测过程中，自然界中很多事件或现象都会产生次声波，所以即使监测到了次声信号，也未必是滑坡次声信号。为了排除噪声干扰，大幅降低误报率，我们就先收集各种声源的次声信号，对其特征进行分析和总结，就是为了能通过特征差异，排除环境噪声的干扰。

(2)利用切比雪夫II型滤波器、快速傅里叶变换(FFT)、短时傅里叶变换(STFT)以及一阶矩和和二阶矩谱估计分析等方法分析与总结这些不同类型与不同固结度的土质滑坡产生的次声信号在时、频域上的相关特征，如短时过零率、峰值频率、中心频率、谱线宽度以及最大有效声压等(均可通过Matlab编程实现)。其中，短时过零率表征信号样本在时域上的波形有无明显突变现象；峰值频率是指功率谱密度最大值所对应的频率；中心频率与谱线宽度可通过一阶矩和和二阶矩谱估计分析方法获得，可表征信号样本的主要能量集中的频响范围E_W(中心频率-谱线宽度/2.0≤E_W≤中心频率+谱线宽度/2.0)；有效声压是指在一定时间间隔中，瞬时声压对时间取均方根值称为有效声压，可表征信号样本的能量强度，如下公式：

式中，“P_e”代表有效值，T代表取平均的时间间隔，它可以是一个周期或比周期大得多的时间间隔，P是指瞬时声压。

滑坡次声的特征提取与各个环境干扰次声的方式相同，下面以滑坡次声的特征提取为例：

通过切比雪夫II型滤波器对训练集中的滑坡次声信号进行滤波降噪处理，以去除随机噪声的影响，以获取较为纯净的滑坡次声信号；海明窗用来对信号进行分帧加窗，以防止能量泄露，用于快速傅里叶变换；快速傅里叶变换和短时傅里叶变换用来分析滑坡次声信号的频率分布情况，提取出峰值频率；一阶矩二阶矩谱估计用来分析与提取滑坡次声信号的中心频率和谱线宽度；短时过零率是指一小段时间内，信号幅值跨越所设置的过零门限(默认为0)的次数，可通过信号幅值(振幅)的变化来计算；小波分解，通过设置分解层数，可获取每层的小波系数。然而，小波系数在这里并没有使用，因为这些信号的小波分解系数相差不大。本实施例中可以不需要小波分解，小波分解的作用多提取一些特征，然后对比一下，看哪些特征的差异比较大，我们就选取这几个特征进行组合。

(3)根据各种类型的土壤，进行固结程度不同以及颗粒级配不同情况下的多组室内滑坡灾变模拟实验(20组)。共筛选出1000段有效次声信号(每段有3000个样本数值)用于分析和测试(800段用于分析获取特征，200段用于后续的测试)。随机选取其中的800段进行分析，总结了每种类型土壤的滑坡次声信号特征如下：

图1各种类型土壤的滑坡次声信号特征表

分析上表可以看出，虽然土壤类型、压实程度与颗粒级配不同，但其相关特征均在一定范围内变化。故，将它们的短时过零率、峰值频率、中心频率和谱线宽度的平均值作为土质滑坡的相关特征值，分别为0.01，2.0，3.6，3.3。有效声压与同时处于灾变状态的滑坡体方量有关，但一般不会超过3Pa。

(4)分别各随机选取山区常见环境干扰次声(风、雷电、爆破、汽车与飞机的引擎)800个样本(剩余的200样本用于测试)，利用(2)中方法手段处理，分析与总结它们的短时过零率、峰值频率、中心频率和谱线宽度以及最大有效声压。具体如下表：

表2山区常见环境干扰次声的信号特征表

根据上述技术方案分析出的滑坡次声与环境干扰次声的相关特征，根据它们的相关特征采用K-近邻算法的思想建立滑坡次声信号智能识别模型，对监测到的信号进行判别。以下为一个实施例。

相关特性介绍：次声监测终端设备部署在滑坡体附近，与监控中心的上位机实时通信，采样频率为1000Hz，即1s时间内采样1000次，得到1000个数据值。监控中心记录监测信号并进行处理分析。该模型每3s进行一次判别，即每次判别有3000个数据值，作为一段待识别信号，记为S1。判别方案如下：

1)滤波降噪处理：利用切比雪夫II型滤波器对S₁进行低通滤波降噪滤除高频部分，只保留次声部分，即＜20Hz的部分；

2)分析有效声压：分析该段待识别信号S₁的有效声压P₁，若P₁＜5Pa，则继续判别，否则视为非滑坡次声信号；

3)分析特征值：利用快速傅里叶变换(FFT)以及一阶矩和和二阶矩谱估计方法分析该段待识别信号S₁的短时过零率、峰值频率、中心频率和谱线宽度这四个特征值，分别为g₁、f₁、m₁、w₁；

4)建立特征向量：利用滑坡、风、爆破、雷电、引擎以及S₁的短时过零率、峰值频率、中心频率和谱线宽度四个特征值建立各次声源的特征向量，即：滑坡：[0.01，2.0，3.6，3.3]，风：[0.11，2.7，5.6，3.9]，爆破：[0.06，4.5，8.3，6.4]，雷电：[0.12，2.8，4.8，3.8]，引擎：[0.22，5.2，10.4，6.7]；S₁：[g₁，f₁，m₁，w₁]

5)特征向量归一化：将待识别信号段S₁与滑坡、风、爆破、雷电、引擎五种已知类别声源的短时过零率、峰值频率、中心频率和谱线宽度四个特征值分别进行2-范数归一化处理。即分别将短时过零率、峰值频率、中心频率和谱线宽度单独做归一化处理。S₁信号段，含有3000个数值(每3s采集的信号作为一个信号段进行识别)。

6)待识别信号分类：将S₁作归一化处理后的特征向量[g₂，f₂，m₂，w₂]分别与滑坡、风、爆破、雷电、引擎作归一化处理后的特征向量进行K-近邻分类(取K＝1，采用欧式距离)，记录与S₁距离最近的声源类型，即为S₁的声源类型，从而实现待识别信号S₁的类别判断。

实际案例1：

一段来自滑坡的次声信号如图3所示：

(1)分析其相关特征，得到短时过零率、峰值频率、中心频率、谱线宽度和有效声压分别是0.008、1.36、3.4、3.1和0.026；

(2)该段信号的有效声压＜5Pa，则继续判别；

(3)该段信号的特征向量：[0.008，1.36，3.4，3.1]；

(4)将该信号段与滑坡、风、爆破、雷电、引擎五种已知类别声源的特征向量分别进行2-范数归一化处理(处理方式如图4，按列进行归一化，然后按行的方式进行判别)，得到处理之后的特征向量分别是：滑坡：[0.03565431，0.242056260.22501318，0.28297267]风：[0.39219737，0.32677595，0.35002051 0.33442224]；爆破：[0.21392584，0.54462658，0.5187804，0.54879547]；闪电：[0.42785167，0.33887876，0.30001758，0.32584731]；引擎：[0.78439473，0.62934627 0.65003809，0.57452026]；待识别信号：[0.02852344，0.16459826，0.21251245 0.26582281]。

(5)欧式距离计算：该信号段特征向量与滑坡的欧式距离约为0.081，与风的欧氏距离约为0.43，与闪电的欧氏距离约为0.45，与爆破的欧氏距离约为0.59，与引擎的欧氏距离约为1；

(6)该信号段特征向量与滑坡的欧式距离最小，则判断该段信号来自滑坡，与实际情况吻合。

实际案例2：

一段来自飞机的次声信号如图5所示：

(1)分析其相关特征，得到短时过零率、峰值频率、中心频率、谱线宽度和有效声压分别是0.25、2.67、10.4、7.0和0.64；

(2)该段信号的有效声压＜5Pa，则继续判别；

(3)该段信号的特征向量：[0.25，2.67，10.4，7.0]；

(4)将该信号段与滑坡、风、爆破、雷电、引擎五种已知类别声源的特征向量分别进行2-范数归一化处理，得到处理之后的特征向量分别是：滑坡：[0.02662174，0.23320717，0.19172548，0.24917905]风：[0.29283918，0.31482968，0.29823963，0.29448433]；爆破：[0.15973046，0.52471614，0.44203373，0.48325634]；闪电：[0.31946092，0.32649004，0.25563397，0.28693345]；引擎：[0.58567835，0.60633865，0.5538736，0.50590898]；待识别信号：[0.66554359，0.31133158，0.5538736，0.52856162]。

(5)欧式距离计算：该信号段特征向量与滑坡的欧式距离约为0.79，与风的欧氏距离约为0.51，与闪电的欧氏距离约为0.52，与爆破的欧氏距离约为0.56，与引擎的欧氏距离约为0.31；

(6)该信号段特征向量与引擎的欧式距离最小，则判断该段信号来自飞机或汽车引擎，与实际情况吻合。

实际案例3：

一段来自雷电的次声信号如图6所示：

(1)分析其相关特征，得到短时过零率、峰值频率、中心频率、谱线宽度和有效声压分别是0.1、4、4.7、3.9和5.5；

(2)该段信号的有效声压＞5Pa，则可判断为非滑坡次声信号，停止判断，判断结果与实际情况吻合。

Claims (9)

1.一种滑坡次声信号识别方法，其特征在于，包括以下步骤：

采集多种不同类型的来自滑坡体的土壤，获取滑坡次声信号；

通过野外现场采集环境干扰次声信号；

分别分析所述滑坡次声信号和干扰次声信号在时、频域上的相关特征以及特征差异，提炼出一系列具有表征性的关键特征；

根据所述关键特征，结合KNN分类算法，建立滑坡次声信号智能识别模型，对监测到的信号进行判别。

2.根据权利要求1所述的一种滑坡次声信号识别方法，其特征在于，对采集到多种不同类型的土壤开展滑坡灾变室内模拟实验，从而获取滑坡次声信号。

3.根据权利要求1所述的一种滑坡次声信号识别方法，其特征在于，通过野外现场采集的环境干扰次声信号包括以下信号的任意组合：大风、雷电、爆破以及飞机和汽车的引擎产生的次声信号。

4.根据权利要求1所述的一种滑坡次声信号识别方法，其特征在于，所述分析所述滑坡次声信号和干扰次声信号在时、频域上的相关特征及特征差异，提炼出具有表征性的关键特征，包括以下步骤：

通过切比雪夫II型滤波器对采集到的滑坡次声信号进行滤波降噪处理，以去除随机噪声的影响，同时获取信号的最大振幅，通过信号振幅的变化来计算短时过零率；

通过海明窗对去除随机噪声的滑坡次声信号进行分帧加窗；

通过快速傅里叶变换和短时傅里叶变换对分帧加窗处理后的滑坡次声信号进行频率分布情况分析，提取出峰值频率；

通过一阶矩二阶矩谱估计分析并提取分帧加窗处理后的滑坡次声信号的中心频率和谱线宽度。

5.根据权利要求1所述的一种滑坡次声信号识别方法，其特征在于，所述关键特征包括短时过零率、峰值频率、中心频率、谱线宽度以及最大有效声压。

滑坡次声信号/干扰次声信号的特征向量是[短时过零率，峰值频率，中心频率，谱线宽度]。

6.根据权利要求1所述的一种滑坡次声信号识别方法，其特征在于，所述根据所述相关特征建立滑坡次声信号智能识别模型，对监测到的信号进行判别，包括以下步骤：

滤波降噪处理：利用切比雪夫II型滤波器对待识别信号S₁进行低通滤波降噪，滤除高频部分，只保留次声部分；

分析有效声压：分析该段待识别信号S₁的有效声压P₁，若P₁＜5Pa，则继续判别，否则视为非滑坡次声信号；

分析特征值：利用快速傅里叶变换以及一阶矩和和二阶矩谱估计方法分析该段待识别信号S₁的短时过零率、峰值频率、中心频率和谱线宽度这四个特征值，分别为g1、f1、m1、w1；

建立特征向量：利用滑坡、风、爆破、雷电、引擎以及S1的短时过零率、峰值频率、中心频率和谱线宽度四个特征值建立各次声源的特征向量；

特征向量归一化：将待识别信号S₁与已知类别声源的短时过零率、峰值频率、中心频率和谱线宽度四个特征值分别进行2-范数归一化处理；

待识别信号分类：将待识别信号S1作归一化处理后的特征向量分别与滑坡、风、爆破、雷电、引擎作归一化处理后的特征向量进行K-近邻分类，记录与待识别信号S₁距离最近的声源类型，即为S₁的声源类型，从而实现待识别信号S₁的类别判断。

7.根据权利要求6所述的一种滑坡次声信号识别方法，其特征在于，所述各次声源的特征向量包括：

滑坡：[0.01，2.0，3.6，3.3]，风：[0.11，2.7，5.6，3.9]，爆破：[0.06，4.5，8.3，6.4]，雷电：[0.12，2.8，4.8，3.8]，引擎：[0.22，5.2，10.4，6.7]；S1：[g1，g1，m1，w1]。

8.根据权利要求6所述的一种滑坡次声信号识别方法，其特征在于，所述K-近邻分类取K＝1，采用欧式距离。

9.根据权利要求6所述的一种滑坡次声信号识别方法，其特征在于，还包括对收集的滑坡次声信号和环境干扰次声信号进行划分，按照4∶1的比例，随机分为训练集和测试集，对训练集执行分析所述滑坡次声信号和干扰次声信号在时、频域上的相关特征以及特征差异步骤，所述测试集用于验证：利用各类型声源剩余的若干段信号对所述滑坡次声信号智能识别模型进行验证。