基于模糊聚类和赤池信息准则的微震信号到时拾取方法

1.本发明属于微震监测技术领域，具体涉及一种基于模糊聚类和赤池信息准则的微震信号到时拾取方法，用于识别微震信号并拾取初至时刻。


背景技术：

2.随着我国煤炭开采深度的增加，面临的煤矿动力灾害更加频繁，因此冲击地压、顶板灾害等灾害的有效预警是目前亟需解决的工程问题。微震监测技术利用检波器采集大规模岩体失稳、岩层滑移等产生的微震信号，并从微震信号中提取出有效信息，推测出煤岩体的力学状态，能够有效地预防煤矿动力灾害。
3.微震p波初至时刻的拾取是微震定位的基础，其拾取精度的高低直接影响微震监测的预警效果。国内外研究人员提出了许多p波到时拾取方法，其中长短时窗均值比法(sta/lta)与赤池信息量准则法(aic)是目前应用最广泛的方法。sta/lta方法基于波形振幅的突然变化来判断是否为微震信号，算法简单，计算速度快，但拾取精度不高，并且当信噪比较低时难以拾取。aic方法通过求取局部波形段的aic函数最小值所对应的时间作为p波初至时刻，该方法合理应用的前提是已知信号初至时刻的大致位置；而传统aic方法是利用sta/lta方法初步确定微震信号的初至时刻，对信噪比水平较低的微震信号不敏感，拾取效果不稳定，甚至会漏掉微震p波，并依赖于经验选取阈值。
4.为满足微震监测系统的精确性、实时性和稳定性，亟需一种计算效率高、拾取精度高、抗噪声能力强的自动拾取方法。


技术实现要素：

5.针对现有aic拾取方法中的不足，本发明所要解决的技术问题就是提供一种基于模糊聚类和赤池信息准则的微震信号到时拾取方法，它能增强抗噪声能力，能初选更准确的p波初至时间窗，对连续的微震数据进行可靠分析，能提高aic拾取法对低信噪比信号的拾取效果。
6.本发明所要解决的技术问题是通过这样的技术方案实现的，它包括以下步骤：
7.步骤1、输入微震序列数据，对原始微震数据进行滤波；
8.步骤2、计算用于聚类的特征函数的时间序列，该特征函数的时间序列包括功率e(i)、峰度k(i)和sta/lta值r(i)；
9.步骤3、通过模糊c均值算法得到聚类结果：
10.模糊c均值算法设置有聚类的个数c和一个用来度量全局最优的目标函数jm，在计算目标函数的过程中，给予并更新每一个数据独有的隶属度u，并通过隶属度u的更新得到每个类别的聚类中心c；通过反复迭代训练，降低目标函数的误差值，当目标函数接近收敛时，得出聚类中心c，并根据数据的隶属度u判定数据所归属的聚类；
11.步骤4、确定包含微震p波到时的到时区间t0；
12.步骤5、精确找出微震信号初至时刻：
13.计算到时区间t0内的aic函数值，并作出aic曲线，再拾取到时区间内aic曲线的最小值点作为精确的微震信号初至时刻。
14.本发明的技术效果是：
15.本发明在步骤2、3使用模糊c均值算法(fcm)改进了现有的aic拾取法，将传统sta/lta方法替换为模糊聚类方法，模糊聚类具有更强的抗噪声能力且计算并不复杂，作为微震信号初拾的方法更为稳健，特别是对低信噪比的信号，解决了现有aic拾取法拾取精度高但计算效率低、低信噪比信号拾取不稳定的问题。
附图说明
16.本发明的附图说明如下：
17.图1为本发明的流程图；
18.图2是微震信号原始波形图；
19.图3是特征函数图；
20.(a)是峰度函数；(b)是功率函数；(c)是sta/lta值；
21.图4是隶属度函数图；
22.虚线表示设置的阈值为0.6，左边的括号表示微震p波初动点的位置，右边的括号表示微震p波结束点的位置，括号内的数字分别对应fcm法拾取点和该点的隶属度；
23.图5是aic曲线图；
24.虚线表示到时区间的范围，括号中的数字分别对应aic法拾取点和该点的aic函数值。
具体实施方式
25.下面结合附图和实施例对本发明作进一步说明：
26.聚类属于无监督机器学习方法，不需要大量的训练数据，而是根据对象本身的属性把给定样本划分为不同的类，与神经网络等受监督机器学习方法相比更加灵活简便，在p波到时拾取工作中表现出计算简便、抗噪声能力较强的优势。目前，p波到时拾取方法有allen拾取法、baer-kradolfer拾取法、高阶统计量拾取法等只计算单一特征函数，而模糊聚类算法基于多个特征函数来划分信号模型和噪声模型，对波形变化识别更加敏感，具有更强的抗噪声能力。
27.如图1所示，本发明具体包括以下步骤：
28.步骤1、对原始微震数据进行简单滤波
29.如图2所示一种微震信号原始波形，利用带通滤波器对微震数据进行预处理，在不破坏有用信号的前提下利用简单温和的滤波器去除少量噪声，得到微震序列；
30.步骤2、计算用于聚类的特征函数的时间序列
31.在一个移动时窗内，利用下列公式计算微震序列的功率e(i)、峰度k(i)和sta/lta值r(i)，得到特征函数的时间序列。
32.功率e(i)计算公式如下：
[0033][0034]
短时窗sta反映了微震信号的平均值，长时窗lta反映了背景噪声的平均值，现有sta/lta法以sta与lta的比值的突变作为判定条件，当比值超过设定阈值时，代表有微震事件发生。本步骤使用微震系统中常用的递归sta/lta算法计算sta/lta值：
[0035][0036][0037][0038]
上式中，n是滑动时间窗内的样本数，x(i)表示微震序列的波形数据，i＝1,2,
…
,n，sta和lta分别为信号在i时刻的短时平均值和长时平均值，w表示滑动时间窗长度的一半，nsta和nlta分别表示短时间窗和长时间窗。
[0039]
峰度k(i)是描述分布峰峭程度和尾部粗细的一个高阶统计量，当非高斯信号即微震信号到达时峰度会快速增加。
[0040][0041]
式中，μ和σ分别是有限序列的算术平均值和标准差，m4是有限序列的四阶中心矩。正态分布的随机变量的峰度为3，峰度k(i)与3的差值的正负表示分布的加宽或变窄。
[0042]
如图3所示，上述微震信号经本步骤处理后获得的特征函数。
[0043]
步骤3、计算隶属度函数
[0044]
在计算出峰度k(i)、功率e(i)和sta/lta比值r(i)三个特征函数的时间序列后，通过模糊c均值聚类分析可以得到聚类结果，即隶属度函数。
[0045]
模糊c均值算法(fcm)是目前应用最广泛的模糊聚类方法，可以允许一个数据属于两个或以上的类。模糊c均值算法是由普通c均值算法改进而来，普通c均值算法对于数据的划分是硬性的，而fcm则是一种柔性的模糊划分，聚类结果更加灵活，能更准确地划分噪声和微震信号。
[0046]
fcm需要指定聚类的个数c，还需要一个用来度量全局最优的聚类标准，即目标函数jm，fcm在计算目标函数的过程中，给予并更新每一个数据独有的隶属度u，并通过隶属度u的更新得到每个类别的聚类中心c。通过反复迭代训练，降低目标函数的误差值，当目标函数接近收敛时，得出聚类中心c，并根据数据的隶属度u判定数据所归属的类别。
[0047]
对于含有n个样本点的样本数据集x＝{x1,x2,
…
,xn}，设聚类中心的集合为c＝{c1,c2,
…
,cj}，模糊聚类是基于目标函数jm的最小化：
[0048][0049]
式中，n是微震序列的样本点数，c是聚类的个数，也是中心点的个数，m》1是加权指
数，u
ij
是xi在聚类j上的隶属度，xi表示第i个d维的样本数据(d是xi的维数)，cj是聚类的中心，‖*‖表示xi和cj之间的欧氏距离，其中xi和cj具有相同的维数。
[0050]
计算隶属度函数按以下四个步骤：
[0051]
(1)、初始化隶属度矩阵u＝[u
i,j
]，记为u0；
[0052]
(2)、在k步骤中：用uk计算中心向量ck＝[cj]，使用以下公式：
[0053][0054]
(3)、迭代计算uk，u
k+1
，所用公式为：
[0055][0056]
(4)判断条件：如果||u
k+1-uk||《ε或者达到最大迭代次数100次，则迭代停止，然后使用得到的隶属度函数检测微震事件；否则返回步骤(2)。
[0057]
通常固定m＝2，这里聚类的数量c＝2，特征向量xi的维数为3。最终得到两类隶属度结果：一类代表属于噪声的概率u1，另一类代表属于微震信号的概率u2。
[0058]
步骤4、确定包含微震p波到时的到时区间t0：
[0059]
设置初始阈值a，当微震信号的隶属度u2≥a时，代表一个微震事件的出现，确定u2≥a曲线段，得到微震p波初动点a0和结束点a1以及微震p波的估计时长l0，若l0所包含的样本数小于250，则令df表示微震信号的采样率。
[0060]
如图4所示，u2≥0.6，产生了微震，微震p波初动点a0为34.275s，结束点a1为36.915s，估计时长l0为2.64s，l0包含528个样本点。
[0061]
以初动点a0为基点，在a0前取一段长度为l0的时间窗，在a0后取一段长度为l0的时间窗，得到到时区间t0为[a
0-l0，a0+l0]。
[0062]
将到时区间t0作为var-aic拾取法的计算时窗，var-aic方法是一种快速的无需计算ar自回归系数、直接根据微震记录计算得到aic函数值的方法。
[0063]
当微震信号的隶属度u2＜a时，无微震信号输出。
[0064]
步骤5、精确找出微震信号初至时刻
[0065]
利用赤池信息量准则法，根据下列计算公式得到到时区间t0内的aic函数值，并作出aic曲线，再拾取到时区间内aic曲线的最小值点作为精确的微震信号初至时刻。
[0066]
微震序列x(i)的aic(s)值计算公式如下：
[0067]
aic(s)＝slog(var(x[1,s]))+(n-s-1)log(var(x[s+1,n])
[0068]
式中，n为采样点的总个数，分割点s的取值范围从(a
0-l0)*df到(a0+l0)*df，与到时区间[a
0-l0,a0+l0]对应，var(x[1,s])和var(x[s+1,n])分别表示微震序列在两个时间段内的方差。
[0069]
利用上述计算公式计算aic(s)值，以s为横坐标，aic(s)为纵坐标，并作出aic(s)曲线如图5所示，初至时刻位于33.255s。

技术特征：
1.基于模糊聚类和赤池信息准则的微震信号到时拾取方法，其特征是，包括以下步骤：步骤1、输入微震序列数据，对原始微震数据进行滤波；步骤2、计算用于聚类的特征函数的时间序列，该特征函数的时间序列包括功率e(i)、峰度k(i)和sta/lta值r(i)；步骤3、通过模糊c均值算法得到聚类结果：模糊c均值算法设置有聚类的个数c和一个用来度量全局最优的目标函数j
m
，在计算目标函数的过程中，给予并更新每一个数据独有的隶属度u，并通过隶属度u的更新得到每个类别的聚类中心c；通过反复迭代训练，降低目标函数的误差值，当目标函数接近收敛时，得出聚类中心c，并根据数据的隶属度u判定数据所归属的聚类；步骤4、确定包含微震p波到时的到时区间t0；步骤5、精确找出微震信号初至时刻：计算到时区间t0内的aic函数值，并作出aic曲线，再拾取到时区间内aic曲线的最小值点作为精确的微震信号初至时刻。2.根据权利要求1所述的微震信号到时拾取方法，其特征是：在步骤2中，计算特征函数的时间序列的公式为：的时间序列的公式为：的时间序列的公式为：的时间序列的公式为：的时间序列的公式为：式中，n是滑动时间窗内的样本数，x(i)表示微震序列的波形数据，i＝1,2,
…
,n，sta和lta分别为信号在i时刻的短时平均值和长时平均值，w表示滑动时间窗长度的一半，nsta和nlta分别表示短时间窗和长时间窗内的样本数；μ和σ分别是有限序列的算术平均值和标准差，m4是有限序列的四阶中心矩。3.根据权利要求2所述的微震信号到时拾取方法，其特征是：在步骤3中，计算隶属度函数为：对于含有n个样本点的样本数据集x＝{x1,x2,
…
,x
n
}，设聚类中心的集合为c＝{c1,c2,
…
,c
j
}，模糊聚类是基于目标函数j
m
的最小化：式中，n是微震序列的样本点数，c是聚类的个数，也是中心点的个数，m>1是加权指数，u
ij
是x
i
在聚类j上的隶属度，x
i
表示第i个d维的样本数据，c
j
是聚类的中心，‖*‖表示x
i
和c
j
之间的欧氏距离。4.根据权利要求3所述的微震信号到时拾取方法，其特征是：在步骤3中，计算隶属度函数按以下步骤：(1)、初始化隶属度矩阵u＝[u
i,j
]，记为u0；(2)、在k步骤中：用u
k
计算中心向量c
k
＝[c
j
]：(3)、迭代计算u
k
：(4)判断条件：如果||u
k+1-u
k
||<ε或者达到最大迭代次数，则迭代停止，否则返回步骤(2)。5.根据权利要求4所述的微震信号到时拾取方法，其特征是：在步骤3中，取m＝2，聚类的数量c＝2，特征向量x
i
的维数为3，得到两类隶属度结果：一类代表属于噪声的概率u1，另一类代表属于微震信号的概率u2。6.根据权利要求5所述的微震信号到时拾取方法，其特征是：在步骤4中，设置初始阈值a，确定u2≥a曲线段，得到微震p波初动点a0与结束点a1以及微震p波的估计时长l0，令df表示微震信号的采样率；以初动点a0为基点，在a0前取一段长度为l0的时间窗，在a0后取一段长度为l0的时间窗，得到到时区间t0为[a
0-l0，a0+l0]。

技术总结
本发明公开了基于模糊聚类和赤池信息准则的微震信号到时拾取方法，包括以下步骤：步骤1、输入微震序列数据，对原始微震数据进行简单滤波；步骤2、计算用于聚类的特征函数的时间序列的功率E(i)、峰度K(i)和STA/LTA值R(i)；步骤3、通过模糊C均值算法得到聚类结果；步骤4、确定包含微震P波到时的到时区间T0；步骤5、计算到时区间T0内的AIC曲线，再拾取局部窗口内AIC曲线的最小值点作为精确的微震信号初至时刻。本发明的技术效果是：用模糊聚类方法替换传统STA/LTA方法，具有更强的抗噪声能力，作为微震信号初拾的方法更为稳健，解决了现有AIC拾取法计算效率低、低信噪比信号拾取不稳定的问题。问题。问题。


技术研发人员：罗永江 王思成 张振宇 郭坤勇 胡千庭
受保护的技术使用者：重庆大学
技术研发日：2023.06.26
技术公布日：2023/9/20