一种低信噪比声发射信号的到时拾取方法与流程

本发明属于信号处理技术领域，尤其是涉及一种低信噪比声发射信号的到时拾取方法。

背景技术：

近年来，微震或声发射监测技术已经被广泛运用到室内岩石试验和各类岩石工程应用中，比如水压致裂开采油气、大型水利水电工程地下洞室围岩稳定性监测，深部岩体开挖及岩爆预警测等。微震或声发射监测技术包括波形识别、声波到时拾取、岩石破裂定位、破裂机制解释和灾害预测等，其中声波到时拾取是事件定位和灾害预测最为基础和重要的一步，其拾取精度直接影响到岩石破裂定位的精度并对破裂机制和灾害预测产生重要影响。

目前，声发射到时主要依靠人工拾取，但人工拾取主要依赖于拾取者的经验，易受个人因素影响且拾取效率低，更严重的是在噪音环境下声发射信号到时很难精确拾取。因此国内外学者提出了一系列含噪音信号的到时自动拾取方法。gaci(2014)提出了基于离散波形转换(dwt)的修改的能量比(mer)法和长短时窗均值比(sta/lta)，在信噪比是10到15时，基于dwt的mer方法精确地拾取了地震信号的到时，其精度高于基于dwt的sta/lta方法，但是当信噪比在3以下时，前一种方法拾取的精度却不如后一种方法。liandshang等(2017)基于dwt和sta/lta的改进峰度法拾取了微震信号的p相到时，这种方法在各种信噪比的条件下都可以精确拾取信号初至时刻，但sta/lta方法采用阀值触发点作为到时点，阀值较小时,声发射信号容易被噪音提前触发；阀值较大时,声发射信号可能不被触发或触发较晚,致使误差较大,由此很难找到一个适合不同波形的阀值。zhang等(2003)提出了一种基于dwt的赤池信息准则(aic)的方法来拾取低信噪比的p波到时，这种方法克服了dwt-sta/lta方法阀值选取的问题，但是dwt-aic方法在某些信噪比的条件下拾取精度不如li等(2017)提出的基于经验模态分解(emd)的aic的方法，然而emd在去除声发射信号噪音的同时会损失大量有用信息，致使拾取结果精度降低。

可见现有的低信噪比声发射或微震信号到时拾取方法存在很大的局限,需要研究一种拾取精度高、稳定性好的自动拾取方法。

技术实现要素：

本发明所要解决的技术问题是提供一种低信噪比声发射信号的到时拾取方法,解决aic法拾取低信噪比、尾部震荡信号不稳定的技术问题,该声发射信号到时拾取方法适用性强、准确性高。

针对aic方法对高信噪比声发射信号拾取精度高，但对低信噪比和尾部震荡信号拾取不稳定的技术问题，本发明提出一种新方法，该方法借助mer方法能很好的克服尾部震荡的影响及多尺度(ms)理论既能平滑噪音同时能够较好保持声发射信号特征的粗粒化思想，极大地提高了低信噪比声发射信号到时的准确性。

为实现上述技术目的，本发明采用如下技术方案实现：

一种低信噪比声发射信号的到时拾取方法，包括以下步骤：

步骤1：导入低信噪比声发射信号；

导入低信噪比声发射信号x(t)t＝1,2,…,n，其中采样频率为3-6mhz，采样起始点为声发射信号的起始时刻，采样点的个数为n；

步骤2：依据mer算法按照以下公式确定低信噪比声发射信号x(t)t＝1,2,…,n的大致到时点k1；

mer(i)＝(er(i)|x(i)|)³(2)

其中，x(j)为待拾取信号x(t)第j个采样点的振幅，ne为能量窗口的长度，er(i)为第i点处的能量比，mer(i)为第i点处修改的能量比；

计算全部采样点的修改的能量比，并将最大值所对应的采样点作为声发射信号的大致到时点k1；

基于mer方法拾取信号的大致到时点k1，且取信号初始时刻到大致到时点k1的时间段作为更新后的低信噪比声发射信号，从而能够排除声发射信号尾部震荡的影响；

步骤3：重新选择低信噪比声发射信号x(t)t＝1,2,…,n的长度；

选择大致到时点k1作为更新后的低信噪比声发射信号x(t)t＝1,2,…,k1的长度；

步骤4：采用多尺度粗粒化方法对更新后的低信噪比声发射信号x(t)t＝1,2,…,k1进行处理，分别得到尺度因子q不同取值时的粗粒化时间序列x^q，粗粒化公式为：

其中，q表示尺度因子，x^q(j)为在尺度因子为q时粗粒化时间序列x^q在第j点的粗粒化振幅；[k1/q]为最靠近k1/q的整数；

通过利用ms理论粗粒化信号对排除尾部震荡后的声发射信号分解成多个尺度因子的粗粒化时间序列，在不同的尺度因子下使得信号平滑，从而排除低信噪比的影响。

步骤5：利用aic方法按照以下公式拾取粗粒化时间序列x^q在尺度因子q不同取值时的的到时点mq：

aic(k)＝klog(var{x(1,k)})+([k1/q]-k-1)log(var{x(k+1,[k1/q])})

(4)；

其中,var{x(1,k)}为粗粒化时间序列的第1点至第k点的方差，var{x(k+1,[k1/q])}为粗粒化时间序列的第k+1点至第[k1/q]点的方差，aic(k)为第k点的声发射aic值；

对每个尺度因子q的粗粒化时间序列x^q，计算每个时间点的声发射aic值，并将最小值所对应的时间点作为当前尺度因子q取值时的粗粒化时间序列x^q的到时点mq；

通过利用aic技术拾取每个尺度因子所对应的粗粒化时间序列的到时，并把所有到时与当前尺度因子的乘积中的最小值，作为信号的最终精确到时。

步骤6：根据尺度因子q与该尺度因子q取值时的粗粒化时间序列x^q的到时点mq，按以下公式确定声发射信号的精确到时l：

lq＝q×mq(5)

l＝min[l1,l2,…lq](6)

其中，lq为在尺度因子为q的低信噪比声发射信号x(t)的到时点。

本方案综合使用mer、ms和aic技术，有效排除了低信噪比、尾部震荡的影响，极大地提高了声发射信号到时拾取精度。

进一步地，能量窗口的长度ne的取值大小为：2t≤ne≤3t,t为低信噪比声发射信号的周期。

进一步地，所述步骤4中，尺度因子q的取值为1至20的整数值。

进一步地，所述步骤1中，低信噪比声发射信号x(t)t＝1,2,…,n的采样频率为3mhz，总采样点数n的取值范围为：3000≤n≤20000。

有益效果

本发明提供了一种低信噪比声发射信号到时自动拾取方法，包括如下步骤：导入低信噪比声发射信号；依据mer算法按照以下公式确定低信噪比声发射信号的大致到时点k1；选择大致到时点k1作为更新后的低信噪比声发射信号的长度；采用多尺度粗粒化方法对更新后的低信噪比声发射信号进行处理，分别得到尺度因子q不同取值时的粗粒化时间序列x^q；利用aic方法拾取尺度因子q不同取值时的粗粒化时间序列x^q的到时点mq；选择尺度因子q与该尺度因子的到时点mq的乘积的最小值作为低信噪比声发射信息的精确到时l。本发明在拾取有尾部噪音的声发射信号时，基于mer方法能够排除声发射信号的尾部震荡的影响；拾取低信噪比声发射信号时，利用ms理论粗粒化信号序列，在不同的尺度下使得信号得到平滑，从而排除低信噪比的影响。由此有效排除了低信噪比、尾部震荡的影响，极大地提高了声发射信号到时拾取精度。

附图说明

图1是本发明所述方法流程图。

图2是本发明对低信噪比声发射信号进行到时拾取的步骤图，其中(a)表示步骤1导入的低信噪比声发射信号；(b)表示步骤2依据mer算法确定大致到时点k1＝3089；(c)表示步骤3重新确定长度后的得到的更新后的低信噪比声发射信号；(d)表示步骤4采用多尺度粗粒化方法且尺度因子分别取值1-20所得到的一组粗粒化时间序列x^q；(e)表示步骤5利用aic方法拾取粗粒化时间序列x^q在尺度因子分别取值1-20时的到时mq；(f)表示步骤6得到的尺度因子q与每个尺度下的到时mq的乘积最小所对应的粗粒化时间序列x¹⁸。

图3是手动、aic、本发明方法拾取声发射信号典型案例图，其中实线、箭头、数值是每种方法的拾取结果。(a)表示手动拾取方法案例图；(b)表示aic拾取方法案例图；(c)表示本发明方法案例图，当q＝19时，q和该尺度下拾取结果62的乘积最小，最终拾取到时为1178。

图4是aic方法和本发明方法分别和手动拾取方法拾取结果对比残差图。(a)表示aic方法和手动拾取方法的拾取结果对比残差图；(b)表示本发明方法和手动拾取方法的拾取结果对比残差图。

具体实施方式

下面结合附图1～4，对本发明提出的一种低信噪比声发射信号的到时拾取方法作进一步说明。

如图1所示，本发明提供的一种低信噪比声发射信号的到时拾取方法，包括以下步骤：

步骤1：导入低信噪比声发射信号；

导入低信噪比声发射信号x(t)t＝1,2,…,n，其中采样频率为3mhz，采样起始点为声发射信号的起始时刻，采样点的个数n＝8000，如图2(a)所示。

步骤2：依据mer算法按照以下公式确定低信噪比声发射信号x(t)t＝1,2,…,n的大致到时点k1；

mer(i)＝(er(i)|x(i)|)³(2)

其中，x(j)为待拾取信号x(t)第j个采样点的振幅，ne为能量窗口的长度，er(i)为第i点处的能量比，mer(i)为第i点处修改的能量比；

计算全部采样点的修改的能量比，并将最大值所对应的采样点作为声发射信号的大致到时点k1，计算得到的k1＝3089，如图2(b)所示。

步骤3：重新选择低信噪比声发射信号x(t)t＝1,2,…,n的长度；

选择大致到时点k1＝3089作为更新后的低信噪比声发射信号x(t)t＝1,2,…,k1的长度，即从导入的低信噪比声发射信号中截取信号起始点到大致到时点3089处的长度区间内的信号，作为更新后的低信噪比声发射信号x(t)t＝1,2,…,k1，如图2(c)所示。

步骤4：采用多尺度粗粒化方法对更新后的低信噪比声发射信号x(t)t＝1,2,…,k1进行处理，分别得到尺度因子q不同取值时的粗粒化时间序列x^q，其粗粒化公式为：

其中，q表示尺度因子，x^q(j)为在尺度因子为q时粗粒化时间序列x^q在第j点的粗粒化振幅；尺度因子q分别取值1至20，得到如图2(d)所示的20个粗粒化时间序列x¹、x²、……、x²⁰。

步骤5：利用aic方法按照以下公式拾取粗粒化时间序列x^q在尺度因子q不同取值时的的到时点mq：

aic(k)＝klog(var{x(1,k)})+([k1/q]-k-1)log(var{x(k+1,[k1/q])})(4)

其中,var{x(1,k)}为粗粒化时间序列的第1点至第k点的方差，var{x(k+1,[k1/q])}为粗粒化时间序列的第k+1点至第[k1/q]点的方差，aic(k)为第k点的声发射aic值；

对每个尺度因子q的粗粒化时间序列x^q，计算每个时间点的声发射aic值，并将最小值所对应的时间点作为当前尺度因子q取值时的粗粒化时间序列x^q的到时点mq。在尺度因子取值1-20的到时结果如图2(e)所示。

步骤6：根据尺度因子q与该尺度因子q取值时的粗粒化时间序列x^q的到时点mq，按以下公式确定声发射信号的精确到时l：

lq＝q×mq(5)

l＝min[l1,l2,…lq](6)

其中，lq为在尺度因子为q的低信噪比声发射信号x(t)的到时点。

经过计算发现，当尺度因子q＝18时，粗粒化时间序列x¹⁸的到时点m18＝142，尺度因子和在该尺度下的到时的乘积最小，其值为2556，因此选择该值作为最终的信号精确到时点。

实施实例1：

图3(a)～(c)分别是手动、aic、本发明方法拾取声发射信号典型案例图，由图可知：aic拾取低信噪比信号时，由于受尾部震荡的干扰，出现多个局部极小值(1227、2775所在点)，全局最大值对应的时间点远大于手动拾取的到时；另外在到时附近aic时间序列尽管有局部极小值，但其局部极小值对应的时间点(1227)和手动到时(1172)相差55个点，明显受噪音影响。由此对于声发射信号，尤其是尾部噪音比较严重的，aic方法拾取精度较差；而本发明方法对于低信噪比信号或者尾部噪音严重的声发射信号拾取精度非常高，成功解决了aic方法低信噪比信号到时拾取不稳定和精度较差的技术问题。

实施实例2：

图4是aic方法和本发明方法分别与手动拾取结果对比残差图。100组声发射信号了从声发射监测系统ds5-16c中记录的岩石破裂事件随机抽取得到。本发明方法计算参数如下：声发射信号采样频率3mhz，采样点个数在8000-30000之间，mer方法的滑动窗口长度为30个采样点。

由图4可知aic方法拾取误差在10个采样点以内、10-20、20-30、30-40、40-50以及大于50占总事件的比例分别为17％、5％、8％、1％、0％、69％；本发明方法对应上述拾取误差划分的占比分别为47％、18％、21％、5％、3％、6％。由此可知，本发明方法对于含有噪音的信号拾取精度较低的比例极低，极大地提高了含噪声发射信号的到时自动拾取的准确性。

以上所述仅为本发明的实施例而已,并不用以限制本发明,凡在本发明精神和原则之内,所作任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。