一种工程尺度下微震信号及p波初至自动识别算法的制造方法与工艺

本发明涉及微震监测技术领域。具体涉及一种工程尺度下微震信号及p波初至自动识别算法，该方法可广泛用于矿业工程、水利水电工程、石油工程、岩土工程以及地下工程。

背景技术：
微震源定位是微震监测与灾害预警的重要组成部分，而微震信号的识别及P波初至拾取是微震源定位的关键。现有的微震信号识别方法大部分来源于地震领域，方法很多，主要有：根据在时间域能量和能量变化构建特征函数时间域的STA/LTA算法；根据地震信号到达前后地震波形数据统计性质的差别，如AIC算法；此外还有神经网络法、高阶统计法、以及基于小波理论的小波变换法等。STA/LTA及其改进算法，由于算法简单、计算效率高、适于实时处理，在地震领域被广泛应用。这些方法用于地震信号的处理具有较好的适用性，但工程尺度下产生的岩石破裂信号与天然地震信号不同：1)地震信号频率一般为几十赫兹以下，而微震信号的频率一般为几十到几百赫兹，有的高达几千赫兹；2)地震信号持续时间较长，一般大于1秒，微震信号持续时间较短，一般不超过0.1秒；3)微震与天然地震相比，一般震级小，较小的微震事件不易被拾取；4)工程环境下由于机械施工、电气干扰、车辆行驶等因素造成背景噪音复杂，使得监测到的大多数微震信号信噪比较低。因此，将地震领域的信号识别及P波到时拾取方法用于工程尺度微震信号是不合适的，主要表现为：1)难以识别低信噪比的微弱信号；2)P波初至自动拾取误差较大。为此，针对上述不足，发明一种可以有效识别工程尺度下微震信号及其p波初至的自动识别的方法，提高工程尺度下微震信号及其P波初至的拾取准确率，进而提高岩爆、矿震、塌方等地质灾害预警的及时性与准确性，是非常必要的。

技术实现要素：
本发明的目的在于克服现有技术存在的问题，提供了一种工程尺度下微震信号及p波初至自动识别算法，提高工程尺度下微震信号及其P波初至的拾取准确率，进而提高岩爆、矿震、塌方等地质灾害预警的及时性与准确性。一种工程尺度下微震信号及p波初至自动识别算法，包括以下步骤：步骤1、读取振动传感器实时监测到的设定时窗内波形的采样点数据，建立平面直角坐标系，其中X轴对应波形的采样点编号，Y轴对应波形的幅值，在平面直角坐标系中对时窗内波形的采样点进行对称校准；步骤2、计算采样点个数N，设置零交叉点个数上限M1，零交叉点个数下限M2；设置微震信号时间判断阀值T；步骤3、设定第A个采样点的加权因子K(A)和特征函数CF(A)；设定第A个采样点短时平均值STA(A)和长时平均值LTA(A)；设定第A个采样点的动态阈值r(A)；其中A为采样点个数的编号；步骤4、以平面直角坐标系中编号最小的采样点为当前采样点，当前采样点的编号为i；步骤5、设置算法参数变量M、S、L、t，M、S、L、t的初值均设置为0，其中M为零交叉点个数；S为计数器；L的计算公式为L＝3+M/3；t为拾取出的微震信号的持续时间；步骤6、比较当前采样点的短时平均值STA(i)和动态阈值r(i)的大小；步骤7、若当前采样点的短时平均值STA(i)小于当前采样点的动态阀值r(i)，则当前采样点不是P波初至点，此时将当前采样点的下一个采样点作为当前采样点，返回到步骤(6)；若当前采样点的短时平均值STA(i)大于等于当前点的动态阀值r(i)，则将当前采样点的编号i赋值给k，计算当前采样点的幅值P1，将当前采样点的幅值大小P1赋值给缓存最大幅值P0，然后进行步骤(8)；步骤8、计算编号为k的采样点的下一个采样点k+1的幅值P2；步骤9、如果采样点的编号为k的采样点的下一个采样点k+1不是零交叉点，将缓存最大幅值P0和采样点k+1的幅值P2中较大值赋值给缓存最大幅值P0，将k+1赋值给k，当k的值小于波形采样点个数N时，返回步骤8；当k的值大于等于波形采样点个数N时，则将当前采样点的下一个采样点作为当前采样点，返回步骤5；如果采样点的编号为k的下一个采样点k+1是零交叉点，则零交叉点个数M值加1，将k+1赋值给k，然后进行步骤10；步骤10、若零交叉点个数M值大于零交叉点个数上限M1值，则当前采样点不为P波初至点，则将当前采样点的下一个采样点作为当前采样点，返回步骤5；若零交叉点个数M值小于等于零交叉点个数上限M1值，则计算当前采样点的判定参数δ＝LTA(i)×M，然后进行步骤11；步骤11、比较编号为k的采样点的短时平均值STA(k)与当前采样点的判定参数δ的大小：若编号为k的采样点的短时平均值STA(k)小于等于当前采样点的判定参数δ，计数器S值归零，将P2赋值给P0，此时将k加1，返回步骤8；若编号为k的采样点的短时平均值STA(k)大于当前采样点的判定参数δ，计数器S值加1，计算L＝3+M/3，进入步骤12；步骤12、比较计数器S值与L的大小：若计数器S值小于等于L，则将P2赋值给P0，将k加1，返回步骤8；若计数器S值大于L，计算从当前采样点i到编号为k的采样点的这段波形的持续时间t，t＝(k-i)/f，其中f为采样频率；步骤13、当t小于微震信号时间判断阀值T或零交叉点个数M小于零交叉点个数下限M2时，则当前采样点i与编号为k的采样点之间的采样点不是微震信号，将编号为k的采样点的下一个采样点作为当前采样点，清零M、P0、P1、P2、S、L的取值，进行步骤5；当t大于等于时间阀值T且零交叉点个数M大于等于零交叉点个数下限M2时，则当前采样点i与编号为k的采样点之间的采样点为微震信号，将当前采样点的编号i赋值给T1，将k赋值给T2，其中编号为T1的采样点为微震信号的P波初至点，编号为T2的采样点为信号结束点，进入步骤14；步骤14、将编号为k的采样点的下一个采样点作为当前采样点，清零M、P0、P1、P2、S、L的取值，进行步骤5，直至将设定时窗内波形的采样点数据分析完。如上所述的对称校准包括以下步骤：选取时窗中设定个数的连续的采样点作为校正样本，若校正样本中幅值为正的采样点个数与幅值为负的采样点个数的比值不在设定比值范围内，则设定时窗内波形的采样点数据需对称校准，计算校正样本幅值的平均值，将设定时窗内波形的采样点幅值与校正样本幅值的平均值相加，完成对称校准；若校正样本中幅值为正的采样点个数与幅值为负的采样点个数的比值不在设定比值范围内，则设定时窗内波形的采样点数据不需对称校准。如上所述的第A个采样点的加权因子K(A)基于以下公式：所述的第A个采样点的特征函数CF(A)基于以下公式：其中，A为对称校准后波形中采样点个数的编号，y(A)为对称校准后波形的第A个采样点的幅值，为y(A)一阶差分；所述的第A个采样点的短时平均值STA(A)基于以下公式：STA(A)＝STA(A-1)+C3×[CF(A)-STA(A-1)]所述的第A个采样点的长时平均值LTA(A)基于以下公式：LTA(A)＝LTA(A-1)+C4×[CF(A)-LTA(A-1)]所述的第A个采样点的动态阈值r(A)基于以下公式：r(A)＝C5×LTA(A)其中，C3为短时平均系数，C4为长时平均系数，C5为触发阀值。当P波到来时，将引起微震信号振幅或者频率等波形特征变化，若信号的信噪比越低，则P波初至越不明显。因而提高拾取算法对该类变化的敏感性，将能有利于提高工程尺度下微震信号的P波初至拾取能力。而算法中的加权因子K、特征函数CF、动态长短时平均比值ε均能反应波形的特征变化，直接影响拾取准确率，动态长短时平均比值ε为长时平均值LTA(A)与短时平均值STA(A)的比值：分析本发明的加权因子K、特征函数CF、动态长短时平均比值ε对信号振幅或者频率波形等特征变化的敏感性，分别如图2、3、4。对比结果分析本发明三种变量在频率或振幅变化处数值改变明显，达到峰值的迅速，因此在理论上本发明对信噪比低的信号及其P波初至的拾取效果明显。从锦屏引水隧洞实时监测的数据中随机选取555个微震信号为样本，定义两个拾取结果判断准则：1)信号成功拾取：能将样本信号自动识别出来而不产生漏波与误判。这样将有利于人工对自动拾取出来的信号部分进行进一步处理。信号成功拾取个数越多，信号拾取率越高。2)P波初至拾取准确：将自动算法拾取P波初至结果与人工拾取结果相比，若拾取误差不大于3个采样点，则认为该算法自动拾取结果准确。本发明信号拾取率为94.23％，P波拾取准确率为73.51％，拾取率完全满足微震领域目前定位误差。分析2015年6月15日至2015年7月15日锦屏极深地下实验室微震监测数据，一共524个微震信号，109个事件，将人工拾取P波、本发明拾取结果用于定位，结果分别如图6、7所示，其中7#实验室是岩石破裂高风险区域，可以看出本发明定位结果微震事件的聚集度较高，微震事件空间分布规律与人工拾取结果接近，与工程实际情况相符。将人工拾取过程中信噪比非常低、P波不明显、定位结果明显异常的13个事件，如图6中的三角形所示，用本发明拾取定位分析：1)本发明拾取的P波结果均能成功定位，如图7中的三角形所示；2)本发明该13个事件拾取定位效果相对于人工拾取改进明显，绝大部分在定位误差允许的范围内。有益效果：本发明针对工程尺度下微震信号特点，改进地震领域的识别算法，本发明能够准确的自动拾取工程尺度下的微震信号特别是性噪比较低的弱信号，同时提高了自动拾取P波初至的准确率，进而提高岩爆、矿震、塌方等地质灾害预警的及时性与准确性。附图说明图1为本发明流程图；图2为本发明加权因子K对振幅或频率变化敏感性分析图；如图2(a)当振幅发生变化时，本发明的加权因子振幅约增大2倍；如图2(b)当频率发生变化时，本发明K值变化幅度较为明显。图3为本发明特征函数CF对振幅或频率变化敏感性分析图；图3(a)模拟背景噪音中出现幅值较小的微震信号，即幅值由1突变到2时，本发明特征函数在突变点的幅值由1到16，突变前后波形对称轴由Y＝1变为Y＝16；本发明在突变前后幅值变化和对称轴偏移明显，特征函数曲线突变迅速，图3(b)模拟信号出现频率变化时，本发明特征函数在突变点的幅值由1到15.6，幅值变化明显，幅值大小增长迅速。图4为本发明ε比值对振幅或频率变化敏感性分析图；图4(a)所示：本发明ε变化曲线在波形振幅变化后3个采样内即达到峰值，曲线变化明显，；如图4(b)当频率变化时，人工能准确拾取出本发明频率改变位置，曲线变化明显。图5为本发明实例1拾取过程S-L值变化图；图6人工拾取微震信号P波定位效果左视图；图7本发明自动拾取微震信号P波定位效果左视图。具体实施方式为了使本发明的技术手段、创作特征、工作流程、使用方法达成目的与功效易于明白了解，下面结合具体实施例，进一步阐述本发明。本发明的保护范围不受以下实例的限制。锦屏地下实验室垂直岩石覆盖达2400米，是目前世界岩石覆盖最深的实验室，工程区的最大主应力可达63MPa，属典型的高应力区，其中7#、8#实验室工程施工时已遭遇到强～极强岩爆，给现场人员和设备安全造成了严重威胁和损失。现场微震监测过程中发现，由于传感器距离岩爆风险较高的7-8#实验室施工段距离较远，大概为400～500m，微震信号在传播过程中衰减严重，且现场监测背景噪音复杂，监测到的微震信号大多数信噪比低，能量幅值小，P波初至不明显，致使微震信号不易识别，P波不易拾取。利用该工程微震样本信号，将本发明拾取和人工拾取结果对比分析，利用粒子群算法优化出适于该工程的本发明最佳参数组：C3＝0.25，C4＝0.08，C5＝2，上限M1＝300，下限M2＝15，T＝0.01。本实例以分析该工程某时窗内实时监测采样点数据为例加以说明。实例1：一种工程尺度下微震信号及p波初至自动识别算法，流程如图1所示，其步骤如下：(1)、读取振动传感器实时监测到的设定时窗内波形的采样点数据，建立平面直角坐标系X0Y，其中X轴对应波形的采样点编号，从平面直角坐标系的原点到X轴正轴方向编号从1开始，按照编号从小到大排列采样点，采样点编号为正整数，Y轴对应波形的幅值，在平面直角坐标系中对时窗内波形的采样点进行对称校准，进入步骤(2)。波形对称校准：选取时窗中设定个数(在本例中为1000个)连续的采样点作为校正样本，若校正样本中幅值为正的采样点个数与幅值为负的采样点个数的比值不在设定比值范围内(在本例中为小于0.95或大于1.05)，则选取的时窗内波形的采样点数据需对称校准，计算校正样本幅值的平均值，将步骤(1)中设定时窗内波形的采样点幅值与校正样本幅值的平均值相加，完成对称校准，这样有利于准确统计波形中的零交叉点个数M，减小算法拾取误差。若校正样本中幅值为正的采样点个数与幅值为负的采样点个数的比值在设定比值范围内(在本例中为大于等于0.95小于等于1.05)，则步骤(1)中设定时窗内波形的采样点数据不需对称校准。(2)计算对称校准后平面直角坐标系中波形的采样点个数N。设置零交叉点个数上限M1＝300，零交叉点个数下限M2＝15；设置微震信号时间判断阀值T＝0.01，单位为秒。微震信号持续时间短，设置零交叉点个数上限M1和零交叉点个数下限M2，能有效避免持续时间较短的电气尖峰信号干扰和对大部分持续时间较长的非微震信号误判。(3)、利用对称校准后平面直角坐标系中的波形的采样点数据计算本发明的第A个采样点的加权因子K(A)：和第A个采样点的特征函数CF(A)：其中A为对称校准后波形中采样点个数的编号，y(A)为对称校准后波形的第A个采样点的幅值，为y(A)一阶差分。计算短时平均值STA(A)和长时平均值LTA(A)：STA(A)＝STA(A-1)+C3×[CF(A)-STA(A-1)]LTA(A)＝LTA(A-1)+C4×[CF(A)-LTA(A-1)]其中STA(A)为第A个采样点的短时平均值；C3为短时平均系数，实例1中C3取值为0.25；LTA(A)为第A个采样点长时平均值；C4为长时平均系数，实例1中C4取值为0.08。计算r(A)＝C5×LTA(A)，r(A)为第A个采样点的动态阀值，C5为触发阀值，实例1中C5值为取2。(4)、以平面直角坐标系中编号最小的采样点为当前采样点，当前采样点的编号为i。(5)、设置算法参数变量M、S、L、t，初值均设置为0，其中M为零交叉点个数；S为计数器，起计数作用；L的计算公式为L＝3+M/3；t为拾取出的微震信号的持续时间。(6)、比较当前采样点的短时平均STA(i)值和动态阀值r(i)的大小。(7)、若当前采样点的短时平均值STA(i)小于当前采样点的动态阀值r(i)，则当前采样点不是P波初至点，此时将当前采样点的下一个采样点作为当前采样点，返回到步骤(6)；若当前采样点的短时平均值STA(i)大于等于当前点的动态阀值r(i)，则当前采样点可能为P波初至点，将当前采样点的编号i赋值给k，即为k＝i，计算当前采样点的幅值大小P1，为找出当前采样点后第一个峰值，将当前采样点的幅值大小P1赋值给缓存最大幅值P0暂存为最大峰值，然后进行步骤(8)。(8)、计算对称校准后坐标系波形中采样点的编号为k的下一个采样点k+1的幅值P2。(9)、如果对称校准后波形中采样点的编号为k的下一个采样点k+1不是零交叉点，比较缓存最大幅值P0、采样点k+1的幅值P2的大小，将较大值赋值给缓存最大幅值P0，将k+1赋值给k，当k的值小于波形采样点个数N时，返回步骤(8)；当k的值大于等于波形采样点个数N时，则将当前采样点的下一个采样点作为当前采样点，即i＝i+1，返回步骤(5)。如果对称校准后坐标系波形中采样点的编号为k的下一个采样点k+1是零交叉点，则零交叉点个数M值加1，即M＝M+1，将k+1赋值给k，然后进行步骤(10)。(10)、判断若零交叉点个数M值大于零交叉点个数上限M1值，则当前采样点不为P波初至点，则将当前采样点的下一个采样点作为当前采样点，即i＝i+1，返回步骤(5)。若零交叉点个数M值小于等于零交叉点个数上限M1值，则计算当前采样点的判定参数δ＝LTA(i)×M，然后进行步骤(11)。(11)、比较编号为k的采样点的短时平均值STA(k)与当前采样点的判定参数δ的大小：若编号为k的采样点的短时平均值STA(k)小于等于当前采样点的判定参数δ，计数器S值归零，将P2赋值给P0，此时将k加1，即k＝k+1，返回步骤(8)；若编号为k的采样点短时平均值STA(k)大于当前采样点的判定参数δ，计数器S值加1，即S＝S+1，计算L＝3+M/3，进入步骤(12)。(12)、比较计数器S值与L的大小：若计数器S值小于等于L，则将P2赋值给P0，此时将k加1，即k＝k+1，返回步骤(8)；若计数器S值大于L，计算对称校准后坐标系波形中从当前采样点到编号为k的采样点的这段波形的持续时间t，t＝(k-i)/f，其中f为采样频率。(13)、当t小于微震信号时间判断阀值T，T为0.01，单位为秒，或零交叉点个数M小于零交叉点个数下限M2时，M2为15，则判断当前采样点i与编号为k的采样点之间的采样点不是微震信号，需对此时的采样点k后的信号继续拾取，此时将对称校准后波形坐标系中编号为k的采样点的下一个采样点作为当前采样点，清零M、P0、P1、P2、S、L的取值，进行步骤(5)；当t大于等于微震信号时间阀值T且零交叉点个数M大于等于零交叉点个数下限M2时，则判断当前采样点i与编号为k的采样点之间的采样点为微震信号，将i赋值给T1，k赋值给T2，其中编号为T1的采样点为微震信号的P波初至点，编号为T2的采样点为信号结束点，进入步骤(14)。(14)、继续对下一个信号自动分析，对称校准后波形坐标系中编号为k的采样点的下一个采样点作为当前采样点，清零M、P0、P1、P2、S、L的取值，进行步骤(5)，直至将设定时窗内波形的采样点数据分析完。以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。