技术领域
[0001] 本
发明涉及地震
信号处理领域，具体为一种基于波形特征的地震信号检测方法。
背景技术
[0002] 近年来，台阵的出现逐渐取代了使用单个地震台收集数据的方法，收集到的地震数据量变得更加庞大，传统的使用人工来进行地震信号分析的方法效率低下，检测结果也容易受到主观判断的影响，已经不再适用于数据量剧增的时代。
[0003] 因此，地震研究人员提出了多种地震检测
算法，大致分为两类，一类是基于单一特征的检测方法，另一类是基于信号整体特征或多个特征来完成地震信号的检测。如长短时窗比值法(STA/LTA)，F检测算法等，STA/LTA算法能够快速的识别地震信号的初至时刻，分形维理论用来判断
地震波的到时，相关法通过分离信号和噪声实现地震波的初至拾取，
自回归模型能够对地震信号的震相进行拾取，F检测算法是基于方差分析的F统计算法，能够检测出地震信号。这些算法能够实时有效地对地震信号进行分析处理，为事件判断和任务部署争取时间。但是，长短时窗比值法不适用于低
信噪比的信号，F检测算法更适用于相关性强的远震信号。所以在加强对阵列自动处理功能研究的同时，改进地震信号检测算法，快速准确地检测出地震信号并降低虚假检测率是研究的重点。
发明内容
[0004] 本发明为了解决地震信号的检测问题，提出一种基于波形特征的地震信号检测方法，基于贝叶斯的理论知识，对信号的STA/LTA波形提取特征，然后利用高斯过程对于历史事件的特征进行训练获得特征模型，以此来区分地震事件和噪声，提升检测事件的性能，改善实际地震事件信号的检测效率，降低误检测率。
[0005] 基于波形特征的地震信号检测方法，具体步骤如下：
[0006] 步骤一、选取台阵收集的历史事件中的地震信号和噪声信号作为数据集，并进行筛选；
[0007] 筛选具体是：丢弃与历史观测明显不符的样本；
[0008] 步骤二、提取每条地震信号和每条噪声信号中的幅值特征α，比值特征ρ，特定频段
能量均值特征γ并进行归一化，以及能量和特征λ并进行归一化；
[0009] 将每条地震信号和每条噪声信号分别经过STA/LTA处理，提取各信号的幅值特征和比值特征；幅值特征的提取方法是对每条信号的STA/LTA波形取最大值；
[0010] 比值特征的提取方法为：
[0011] 首先，针对每条信号的STA/LTA波形，通过预设
阈值筛选所有大于等于阈值的尖峰值，去除STA/LTA波形中的零分量。
[0012] 然后，统计所有峰值点对应的
位置以及个数，设置拾取峰值宽度的参数大小，剔除将一个峰值识别为多个峰值的情况；
[0013] 最后，将最大的峰值(即幅值特征α对应的峰值)剔除出去，对所有剩余的尖峰值求和取平均值β，计算比值特征
[0014] 特定频段能量均值特征的提取方法是：对每条地震信号和每条噪声信号分别进行
频谱分析(FFT变换)，提取每条信号的频谱的特定频段(1-2Hz)能量，然后求和取平均作为能量均值特征γ，并进行对数归一化，得到归一化后的特定频段能量均值特征log(γ)；
[0015] 能量和特征的提取方法是：对每条地震信号和每条噪声信号分别进行EMD(Empirical mode decomposition)分解，得到固有模态函数IMF1分量的能量和特征λ，并进行对数归一化，得到归一化后的能量和特征log(λ)；
[0016] 步骤三、将所有的地震信号划分为训练样本和测试样本，将所有的噪声信号定义为测试样本；
[0017] 步骤四、将训练样本中所有地震信号的幅值特征α组成幅值矩阵，将所有比值特征ρ组成比值矩阵，将所有归一化后的能量均值特征log(γ)组成能量均值矩阵，将所有归一化后的能量和特征log(λ)λ组成能量和矩阵；
[0018] 步骤五、将各特征矩阵带入高斯函数，并利用
梯度下降法优化，得到各特征分别对应的三个最优超参数；
[0019] 幅值特征α，比值特征ρ，归一化后的能量均值特征log(γ)，以及归一化后的能量和特征log(λ)中每个特征分别对应三个超参数：边际方差 尺度l以及噪声方差[0020] 高斯函数公式如下：
[0021]
[0022] L(θ)表示对数边际似然函数；X表示地震波信号组成的训练集；y表示四种特征矩阵中的一种矩阵；θ表示三个超参数组成的向量，即 Ky表示含高斯噪声的观测值变量；n表示训练集中地震波信号的数目；
[0023] 利用梯度下降法优化高斯函数，当L(θ)达到最大时，得到最优的三个超参数，具体公式如下：
[0024]
[0025] θj依次等于训练集中地震信号的超参数 l；
[0026] 步骤六、利用四个特征分别对应的各自的超参数，计算各特征高斯过程的后验均值和协方差；
[0027] 后验均值 计算公式如下：
[0028] X*表示一个待测信号形成的1*1的矩阵；
[0029] 协方差 计算公式如下：
[0030] 步骤七、利用每个特征的后验均值和协方差，得到高斯分布的关于该特征的模型，共得到四个特征模型；
[0031] 步骤八、利用包含地震信号和噪声信号的测试样本，将每条测试样本中的四个波形特征分别输入各自对应的特征模型中进行验证；
[0032] 步骤九、用贝叶斯的思想预测新的事件发生的概率，根据事件发生概率的大小来判断是否为地震事件。
[0033] 具体为：
[0034] 首先、实时采集一条真实的数据，提取该数据信号中的幅值特征α，比值特征ρ，归一化后的特定频段能量均值特征log(γ)，以及归一化后的能量和特征log(λ)；输入四种特征模型中得到四个特征各自的概率值，将四个独立特征各自的概率值进行相乘，得到该数据信号的整体概率值；
[0035] 然后、判断该数据信号的整体概率值与设定的阈值之间的关系，判断该数据信号为地震信号还是噪声信号。
[0036] 设定的阈值包括两部分：
[0037] 第一部分是将测试集中的所有噪声信号样本分别输入到幅值特征α的高斯分布模型中，得到各样本的幅值概率值；设置噪声的虚警率为5％，提取位于此时的幅值概率值为第一个
门限thre1；
[0038] 第二部分是将利用幅值特征训练得到的高斯分布的均值作为第二个门限thre2；
[0039] 当新的实时采集的数据信号的幅值特征α小于第二门限thre2，且该数据信号的四个特征各自的概率值相乘后的整体概率值低于第一门限thre1时，判定该数据信号为噪声。否则，该数据信号为地震信号。
[0040] 本发明的优点在于：
[0041] 1)、一种基于波形特征的地震信号检测方法，相比于F检测算法，不要求地震事件存在强相关性，适用性更强。
[0042] 2)、一种基于波形特征的地震信号检测方法，相比于传统的STA/LTA检测算法，能够有效地提升正确检测率，降低误判率。
附图说明
[0043] 图1是某个区域发生的历史事件散点图；
[0044] 图2是本发明基于波形特征的地震信号检测方法的
流程图；
[0045] 图3是
实施例中地震信号经过STA/LTA处理后的波形图；
[0046] 图4是实施例中噪声信号经过STA/LTA处理后的波形图；
[0047] 图5是实施例中距离与协方差的关系图。
具体实施方式
[0048] 下面将结合附图和实施例对本发明作进一步的详细说明。
[0049] 本发明考虑的场景是针对某个经常发生地震的地理区域。对于地震频发的地带，地震台阵收集到大量的地震数据，如图1所示，如何更加有效的利用这些数据研究地震波包含的丰富信息是一项重要而艰巨的任务。
[0050] 地震波是地球内部传播的一种弹性波，分为纵波、横波和面波三种。纵波传播速度较快，破坏性较小，又称P波。横波又称S波，紧随P波后到达，破坏性较大。面波是混合波，又称L波，有横纵波在地表相遇造成的，对地面的
建筑物破坏
力较大。同时，采集到的地震信号会受到噪声的影响，噪声的来源有多种，如自然源(海洋微震、
风等)引起的环境振动；人造振动(来自工业、交通)；信号散射；重力效应(如大气
牛顿引力、地表倾斜引起的
水平
加速度)；
传感器技术
缺陷或老化产生的信号；采集设备不可避免的自噪声；
数据处理的产物。因此，如何将有效信号与噪声进行区分以及抑制噪声一直是数据处理的重要工作。
[0051] 本发明利用地球上发生的历史地震事件，产生一个联合概率分布，包括过去、现在以及未来的所有信号。在实践中，模型的许多方面是根据历史数据进行调整的，而这些数据来自于过去发生的事件：
[0052] p(wordfurure|worldpast)∝p(wordfurure,worldpast)
[0053] 因此，需要选取大量的历史事件作为数据集，并对其进行筛选，丢弃与历史观测明显不符的样本，利用筛选后的数据构建模型。此模型用于描述事件与产生信号之间的关系，使得现实发生的事件与之对应。即定义了未观察事件的后验分布，当模型的效果很好时，此时后验分布将主要集中于实际发生的事件上。
[0054] 地震信号的原始波形中包含带大量有用的信息，如初至到时，幅值，尾波等特征，这些都能为研究
震源以及地球的内部结构提供帮助。基于这种思想，在处理和观察大量的地震事件和噪声事件的
基础上，本发明具体提出了基于波形特征的地震信号检测方法，其流程如图2所示，具体步骤如下：
[0055] 步骤一、选取台阵收集的历史事件中的地震信号和噪声信号作为数据集，并进行筛选和预处理；
[0056] 筛选具体是：丢弃与历史观测明显不符的样本；
[0057] 预处理是指：将这个区域的数据通过随机的赋值来进行0或1的编号，选取0.5作为门限，将地震
波数据分为训练集和测试集，相当于数据的均匀化。
[0058] 步骤二、提取每条地震信号和每条噪声信号中的幅值特征α，比值特征ρ，进行对数归一化后的特定频段能量均值特征log(γ)，以及进行对数归一化后的能量和特征log(λ)；
[0059] 将每条地震信号和每条噪声信号分别经过STA/LTA处理，提取各信号的幅值特征和比值特征；本实施例中地震信号经过STA/LTA处理后的波形如图3所示，噪声信号经过STA/LTA处理后的波形如图4所示，根据波形的变化曲线选取的形状参数特征。
[0060] 幅值特征的提取方法是对每条信号的STA/LTA波形取最大值；选取峰值α作为特征，观察可知大部分地震信号的峰值要大于噪声信号的峰值。
[0061] 比值特征的提取方法为：
[0062] 首先，针对每个STA/LTA波形，通过预设阈值筛选所有大于等于阈值的尖峰值，去除STA/LTA波形中的零分量。
[0063] 然后，统计所有峰值点对应的位置以及个数，设置拾取峰值宽度的参数大小，剔除将一个峰值识别为多个峰值的情况；
[0064] 最后，将最大的峰值(即幅值特征α对应的峰值)剔除出去，对所有剩余的尖峰值求和取平均值β，计算比值特征 观察可知大部分的地震信号峰值前后波形变化的程度小于噪声信号的变化程度。
[0065] 特定频段能量均值特征的提取方法是：对每条地震信号和每条噪声信号分别进行频谱分析(FFT变换)，提取每条信号的频谱的特定频段(1-2Hz)能量，然后求和取平均并归一化作为能量均值特征log(γ)；
[0066] 能量和特征的提取方法是：对每条地震信号和每条噪声信号分别进行EMD(Empirical mode decomposition)分解，得到固有模态函数IMF1分量的能量和，并进行归一化后得到特征log(λ)；
[0067] 步骤三、将所有的地震信号划分为训练样本和测试样本，将所有的噪声信号定义为测试样本；
[0068] 步骤四、将训练样本中所有地震信号的幅值特征α组成幅值矩阵，将所有比值特征ρ组成比值矩阵，将所有归一化后的能量均值特征log(γ)组成能量均值矩阵，将所有归一化后的能量和特征log(λ)组成能量和矩阵；
[0069] 步骤五、将各特征矩阵带入高斯函数，并利用梯度下降法优化，得到各特征分别对应的三个最优超参数；
[0070] 幅值特征α，比值特征ρ，归一化后的能量均值特征log(γ)，以及归一化后的能量和特征log(λ)中每个特征分别对应三个超参数：边际方差 尺度l以及噪声方差[0071] 对STA/LTA的波形进行提取特征，进行GP建模，使用历史事件的信息作为先验，使用高斯过程处理先验信息，协方差函数是高斯过程预测值中的关键因素，因为它对想要学习的函数的进行假设编码。GP协方差是一个Matern核，其超参数是由噪声方差 边际方差和特征尺度l组成，训练过程通过最大限度地利用历史数据的边际可能性，为每个
站点、
相位和形状参数选择超参数。
[0072] 高斯过程核函数的参数称为超参数(hyperparameters)，超参数决定了GP表达函数的先验信息。Matern核的超参数包括边际方差 尺度l以及(在噪声观察下)噪声方差使用向量θ来表示超参数。通过
训练数据集(X,y)，最大化对数边际似然L(θ)决定超参数，其中L(θ)通过高斯(对数)
密度得出：
[0073]
[0074] L(θ)表示对数边际似然函数；X表示地震波信号组成的训练集；y表示四种特征矩阵中的一种矩阵；θ表示三个超参数组成的向量，即 Ky表示含高斯噪声的观测值变量；n表示训练集中地震波信号的数目；
[0075] 通常情况下，L(θ)对于Ky是可微的，而Ky对于θ也是可微的，因此可以通梯度优化来寻找最优的超参数。由于这种问题通常的非凸的，所以可能需要多初始化。
[0076] 为了能够最大化边际似然来设置合适的超参数，需要进行梯度优化，此过程需要边际似然的偏导数，偏导数如下：
[0077]
[0078] θj依次等于训练集中地震信号的超参数 l；
[0079] 超参数的
选定是通过梯度下降寻找，具体采用的是Adam算法来进行梯度优化，它利用梯度的一阶矩阵和二阶矩阵估计，可以动态调整每个参数的学习率，使得参数更新比较平稳。
[0080] 超参数的训练过程如下：
[0081] (1)参数初始化；
[0082] (2)计算偏导数矩阵，θj依次等于 l；
[0083] (3)然后根据偏导数计算此时的边际似然值；
[0084] (4)根据学习率更新超参数；
[0085] (5)判断边际似然是否已经最大化，满足条件则输出超参数的结果，不满足则继续执行，直到最大化边际似然值。
[0086] 步骤六、利用四个特征分别对应的各自的超参数，计算各特征高斯过程的后验均值和协方差；
[0087] 在测试点矩阵X*上的函数值f*＝f(X*)是通过(含噪声)训练来预测的。在这种情况下，条件分布呈封闭高斯形式，
[0088]
[0089] 后验均值 计算公式如下：
[0090] X*表示一个待测信号形成的1*1的矩阵；
[0091] 协方差 计算公式如下：
[0092] 协方差函数对想要学习的函数进行假设编码，是高斯过程预测值中的关键因素。在
监督学习中，数据点之间的相似性概念是至关重要的，这是一个基本假设，即接近输入x的点很可能具有相似的目标值，因此接近测试点的训练点应该能够提供关于该点的预测信息。
[0093] 步骤七、利用每个特征的后验均值和协方差，得到高斯分布的关于该特征的模型，共得到四个特征模型；
[0094] 步骤八、利用包含地震信号和噪声信号的测试样本，将每条测试样本中的四个波形特征分别输入各自对应的特征模型中进行验证；
[0095] 步骤九、用贝叶斯的思想预测新的事件发生的概率，根据事件发生概率的大小来判断是否为地震事件。
[0096] 具体为：
[0097] 首先、实时采集一条真实的数据，提取该数据信号中的幅值特征α，比值特征ρ，归一化后的特定频段能量均值特征log(γ)，以及归一化后的能量和特征log(λ)；输入四种特征模型中得到四个特征各自的概率值，将四个独立特征各自的概率值进行相乘，得到该数据信号的整体概率值；
[0098] 单个形状参数利用事件位置和深度进行高斯过程建模。
[0099] 对于每个站j和每个相位k，假定参数遵循事件起源 的未知函数，利用高斯过程处理带噪声的观测值，对提取出的各特征参数分别进行建模，以幅值特征α为例：
[0100] 第i个新信号的STA/LTA波形的幅值αi为：αi＝fα(si)+ε
[0101] 对实际观测信号(含噪声ε)幅值的数学表示为：
[0102] 其他特征参数的建模与上式类似，其中ε表示每个观测信号中的高斯噪声，并被合并到GP协方差中，得到新信号为地震信号的概率为：
[0103]
[0104] 其中 和 是利用测试点X*和训练点E的GP后验求出的均值和协方差，N表示高斯分布，E代表历史事件，α*表示测试点X*的幅值。
[0105] 地震信号的概率判断公式如下：
[0106] p(si|E)＝p(αi|E)p(ρi|E)p(logγi|E)p(logλi|E)
[0107] 其中，p(si|E)代表地震信号发生的概率，E代表在某个区域上发生的历史地震事件，si代表第i个新的信号，αi代表第i个新信号的STA/LTA波形的幅值，ρi代表第i个新信号的STA/LTA波形的比值，γi代表第i个新信号的频域能量均值；p(αi|E)代表在已知历史事件幅值的前提下，新信号幅值是地震事件幅值的概率大小；同理，p(ρi|E)代表是地震事件比值的概率；p(logγi|E)代表是地震事件特定频段能量均值的概率，p(logλi|E)代表是地震事件IMF1分量的能量和的概率。
[0108] 然后、判断该数据信号的整体概率值与设定的阈值之间的关系，判断该数据信号为地震信号还是噪声信号。
[0109] 设定的阈值包括两部分：
[0110] 第一部分是将测试集中的所有噪声信号样本分别输入到幅值特征α的高斯分布模型中，得到各样本的幅值概率值；设置噪声的虚警率为5％，提取位于此时的幅值概率值为第一个门限thre1；
[0111] 第二部分是将利用幅值特征训练得到的高斯分布的均值作为第二个门限thre2；
[0112] 当新的实时采集的数据信号的幅值特征α小于第二门限thre2，且该数据信号的四个特征各自的概率值相乘后的整体概率值低于第一门限thre1时，判定该数据信号为噪声。否则，该数据信号为地震信号。
[0113] 利用贝叶斯的思想通过训练先验的信息来预测新的事件发生的概率，设置阈值，根据事件发生概率的大小来判断是否为地震事件。如图5所示，四条线中除了实线型线外，其余三条线分别表示地震发生位置位于(136°，34°)，(135.4°，37°)，(139.6°，36°)的三个不同的历史事件，横坐标对应的是预测新事件发生的位置，纵坐标表示不同位置用于预测新事件分布的协方差，纵坐标对应的值越大，代表发生新事件发生在历史事件附近的可能性越大，由图5可以看出新事件可能发生的位置依次是(136.11°，34.06°)，(135.3°，37°)，(139.66°，35.98°)。实线型为选取了一个未曾发生过事件的位置(137°，35°)，从图像可以看出，协方差在所有位置都为0，即不能有效地区分信号与噪声。
[0114] 在训练完成后，对地震事件和噪声分别进行检测，在相同误警率的条件下，将本发明基于波形特征的信号检测算法得出的测试结果与传统的STA/LTA的测试结果进行对比，根据测试结果表明，本发明相对于传统的STA/LTA算法，可以有效地提升地震事件检测的性能，降低误检测事件的数量。
