考虑传感器感度的声发射事件震级获取方法、系统及可读存储介质与流程

本发明属于声发射技术领域，尤其涉及考虑传感器感度的声发射事件震级获取方法、系统及可读存储介质。

背景技术：

声发射技术被广泛应用于岩石压缩、无损检测等领域，而声发射震级是研究破裂过程的一个重要依据，其可为声发射事件大小做定量评价，并为声发射b值和分形参数计算的重要参数。声发射震级计算主要参照了gutenberg-richter地震震级经验公式lgn＝a-bml，其中ml为地震震级，n为相应地震震级的频次，a、b为回归参数，b即为b值。richter震级ml与最大振幅amax成正比，即ml∝lgamax。对于声发射事件，gutenberg-richter震级公式可写为lgn＝a-badb/20，其中adb＝20lgamax(colombo等，2003)。由此，声发射震级可记为m＝lgamax。可知，对于一个声发射事件，每个触发的传感器均可得到一个震级。实验发现直接采用每个声发射信号最大振幅计算得到的震级差异性较大。另一方面，传感器在接收信号时，由于凡士林、贴合程度等影响，其感度是不是一样的，由此造成的声发射信号衰减程度不同。

声发射定位通常基于p波到时差，目前，人工p波初至拾取受信噪比、人为情绪、经验等影响，且其处理数据较慢。而长短时窗平均值比法(sta/lta法)、高阶统计量法(pai-s/k法)、池赤准则法(aic法)、互相关法、分形维数法和人工神经网络法等拾取稳定性和精度较人工拾取低。上述方法只利用了单个传感器，而全波形互相关法借助多个声发射信号全波形互相关最大值作为相对到时差，在地震波领域得到了广泛应用。然而，声发射信号的尾波差异性较大，由全波形计算的互相关到时差可靠性很低。

常用的震源定位方法大多从geiger法发展而来，这类方法的思路是利用射线走时模型，建立理论到时与观测到时的残差函数，并在空间中搜索最佳坐标以使残差函数最小。例如，waldhauser和ellsworth(2000)提出了双重残差定位法，其利用相近事件同时反演，消除了路径差异造成的误差，aki和lee(1976)和aki等(1997)分别利用近震和远震数据联合反演震源参数和台站下岩体的速度参数，由此开辟了地震成像领域。在矿山微震领域，董陇军等(2011)的无需测速法将速度看作未知数，免去了测定速度的过程，李健等将单纯形法与无需测速法相结合，提高了求解精度，李楠等(2014)基于l1范数统计和单纯形法解决求解系统发散、定位精度低和定位受微震台网影响大等问题。最近，王泽伟等(2016)提出了一种有别于传统定位思路的方法—虚拟场定位法，其对拾取误差不灵敏的信号源定位方法具有很大的实际应用价值。然而，王泽伟等的方法采用了牛顿迭代算法，其于声发射这种小尺度定位稳定不太好。

综上知，现有全波形互相关p波相对到时拾取不能用于尾波变化较大的信号，已有虚拟场定位采用的牛顿迭代法存在计算不稳定的问题，现有的声发射震级计算只考虑了单个传感器的最大振幅，且未考虑传感器感度的影响。

技术实现要素：

本发明所要解决的技术问题是提供一种考虑传感器感度的声发射事件震级获取方法、系统及可读存储介质，该声发射震级获取方法具有抗噪性强、易于得到全局最优、震级稳定等特点。

一种考虑传感器感度的声发射事件震级获取方法，包括以下步骤：

步骤1：将采集得到的声发射事件的信号集，依据信号时间差，将声发射信号划分至不同的声发射事件的信号集中，对每一个声发射事件的所有声发射信号按照触发先后排序，对排序后的声发射信号进行加窗补长至振幅为0的点，获得各声发射事件的加窗声发射信号集；

步骤2：基于互相关法计算各声发射事件的加窗声发射信号集中，每个加窗声发射信号相对第一个加窗声发射信号的p波走时差(t1,t2,...,tn)，n为声发射信号的数目；

步骤3：基于p波走时差(t1,t2,...,tn)，采用虚拟场法计算发射事件所在空间的各网格点虚拟场势值，取虚拟场势值最大的网格点作为声发射事件位置，进而利用声发射事件位置与传感器之间的距离计算传感器的感度系数；

所述虚拟场法采用已公开的专利一种用于均匀速度场的信号源定位方法进行确定；

步骤4：利用传感器的感度系数，对声发射事件的所有原始声发射信号的最大振幅进行修正；

当一个声发射事件发生后，其通过波的形式传到传感器，传感器将其转换为电压信号，即为声发射信号；

步骤5：依据修正后的声发射信号的最大振幅，确定每个声发射信号的震级，以所有声发射信号的震级中位数作为声发射事件的震级；

其中，mi表示第i个声发射信号的震级，为第i个声发射信号的修正后的最大振幅，单位为微伏，k为震级调整参数，取值范围为[-1,0]。

进一步地，对声发射事件的信号进行加窗补长的具体过程如下：

首先，拾取一个p波初至前的时间t1，再拾取首个p波峰值后的时间t2，接着自动延伸至下一个振幅为0的时间t3，获得加窗信号在区间[t3-max(t3-t1),t3]的时间序列，完成加窗补长；

其中，t3-max(t3-t1)用于统一加窗信号的长度，max(t3-t1)为所有加窗信号中t3-t1的最大值。

在进行互相关拾取前进行加窗处理，加窗后的声发射信号不会受尾部信号的影响，从而提高p波相对到时差的拾取精度。

加窗只需拾取p波初至前时间和首个峰值后且零点前的时间，再按照相同方法对信号自动补长，而无需像人工拾取p波初至需要对信号进行放大和很仔细的拾取，通过加窗补长处理后，利用互相关法能快速得到p波走时差；声发射信号受电流、背景噪音等影响，人工拾取p波初至较为困难，而该加窗处理后只需大致拾取p波初至前的时间，p波首个峰值后且零点前时间。因此，经过此处理后，可以使得所提的数据具有较强的抗噪性；

进一步地，采用变网格搜索法计算发射事件所在空间各网格点的虚拟场势值，其首先采用l×l×lmm³的网格，搜索整个岩样空间范围，计算得到声发射事件的初步位置，再以初步位置为中心生成一个l×l×lmm³的立方体，利用l×l×lmm³的细网格搜索所述l×l×lmm³立方体，得到声发射事件的准确位置；

其中，l的取值范围为[3,8]mm，l的取值范围为[0.3,0.8]mm。

进一步地，所述虚拟场势值采用指数衰减函数进行计算：

其中，σ为常数，取值范围为10≤σ≤200，d与声发射时间误差相关的参数：声发射时间误差越大，d值越大。

大拾取误差的贡献较小，采用指数衰减函数的虚拟场势值能降低大拾取误差的影响，d的计算方法采用已有专利一种用于均匀速度场的信号源定位方法。

进一步地，sci表示第i个传感器的感度系数。

基于传感器感度修正了声发射信号最大振幅，修正后的最大振幅明显较未修正的最大振幅更为一致；

进一步地，采用断铅实验计算每个传感器的感度系数。

进一步地，在10cm×10cm×10cm的方形岩石试样前后左右四个面各布置两个传感器，以左面和后面两个面交线段的底点作为(0,0,0)，水平向东记为正x轴，水平向南记为正y轴，垂直向上记为z轴，1～8号传感器坐标为(5,10,8.5)，(0,5,8.5)，(5,0,8.5)，(10,5,8.5)，(5,10,1.5)，(0,5,1.5)，(5,0,1.5)，(10,5,1.5)cm；在顶面划分1cm×1cm的网格进行断铅实验，并计算每次断铅实验时各个传感器的感度系数，取各个传感器多次断铅的感度系数均值作为每个传感器的感度系数(sc1,sc2,...,scn)。

进一步地，每个传感器的感度系数计算方法为sc＝ap·d/f(n,d)；

其中，ap为p波初至后的首个峰值振幅，d为声发射事件与传感器的距离，d为声发射事件与传感器连线方向的单位向量，n为传感器接收信号方向的单位向量，r为p、s波的速度比值，s＝n·d。

一种考虑传感器感度的声发射事件震级获取系统，包括：

信号加窗补长单元：将采集得到的声发射事件的信号集，依据信号时间差，将声发射信号划分至不同的声发射事件的信号集中，对每一个声发射事件的所有声发射信号按照触发先后排序，对排序后的声发射信号进行加窗补长至振幅为0的点，获得各声发射事件的加窗声发射信号集；

信号p波走时差计算单元：基于互相关法计算各声发射事件的加窗声发射信号集中，每个加窗声发射信号相对第一个加窗声发射信号的p波走时差(t1,t2,...,tn)，n为声发射信号的数目；

声发射事件定位单元：基于p波走时差(t1,t2,...,tn)，采用变网格搜索法计算声发射事件所在空间的各网格点虚拟场势值，取虚拟场势值最大的网格点作为声发射事件的定位结果；

所述虚拟场势值采用已公开的专利一种用于均匀速度场的信号源定位方法进行确定；

声发射信号最大振幅修正单元：利用传感器的感度系数，对声发射事件的所有原始声发射信号的最大振幅进行修正；

当一个声发射事件发生后，其通过波的形式传到传感器，传感器将其转换为电压信号，即为声发射信号；

声发射事件震级确定单元：依据修正后的声发射信号的最大振幅，确定每个声发射信号的震级，以所有声发射信号的震级的中位数作为声发射事件的震级；

其中，mi表示第i个声发射信号的震级，为第i个声发射信号的修正后的最大振幅，单位为微伏，k为震级调整参数，取值范围为[-1,0]。

一种计算机可读存储介质，包括计算机程序指令，所述计算机程序指令被处理终端执行时使所述处理终端执行一种考虑传感器感度的声发射震级获取方法。

《一种用于均匀速度场的信号源定位方法》为申请号：201510973875.4申请日：2015-12-22公开日：2016-05-04的中国专利。

有益效果

本发明提供了一种考虑传感器感度的声发射震级获取方法、系统及可读存储介质，主要解决现有互相关p波相对到时拾取不能用于尾波变化较大的信号，已有虚拟场定位采用的牛顿迭代法计算不稳定，现有的声发射震级计算只考虑了单个传感器的最大振幅，且未考虑传感器感度影响的问题。本方法包括如下步骤：首先利用加窗互相关法计算每个声发射信号相对第一个到达的声发射信号p波走时差；再采用变网格搜索的方式计算虚拟场势值，进而确定声发射事件发生位置；对声发射系统开展断铅实验，计算每个传感器的感度系数；计算每个声发射信号的最大振幅，借助传感器的感度系数得到修正后的最大振幅；计算每个传感器修正后最大振幅对应的震级，再取所有震级的中位数作为该声发射事件震级。该加窗方法只需大致拾取p波初至前时间和首个峰值后时间，再依靠程序自动补长，而无需像人工拾取p波初至需要对信号进行放大和很仔细的拾取，因此该方法能快速得到p波走时差；声发射信号受电流、背景噪音等影响，人工拾取p波初至较为困难，而该方法只需大致拾取p波初至前的时间，p波首个峰值后时间，因此该专利具有较强的抗噪性；虚拟场势值函数采用了指数衰减的形式，大拾取误差的贡献较小，因此该方法能降低大拾取误差的影响。此外，该方法基于传感器感度修正了声发射最大振幅，修正后的最大振幅明显较未修正的最大振幅更为一致(图5)，取其中位数作为声发射震级较传统单一传感器得到的震级更为可靠。此方法具有抗噪性强、易于得到全局最优、震级稳定等特点。

附图说明

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

图2是声发射监测时传感器布置图；

图3是加窗互相关拾取实例图，其中，(a)为信号原始波形，(b)为加窗补长后的波形，(c)为加窗互相关后的时间序列；

图4是变网格虚拟场定位实例图；

图5是声发射震级确定实例图。

具体实施方式

下面将结合附图1～5，对本发明提出的一种考虑传感器感度的声发射事件震级获取方法作进一步说明。

本发明算法思想的描述如下：采用加窗互相关拾取p波相对到时差，降低声发射信号尾波的影响；采用变网格法实现虚拟场定位，可克服牛顿迭代法计算不稳定的问题；进而提出一种考虑传感器感度的声发射震级计算方法，得到稳定的声发射震级。

如图1所示，一种考虑传感器感度的声发射事件震级获取方法，包括以下步骤：

步骤1：加窗互相关法计算声发射信号p波走时差；

将触发时间差小于40us的声发射信号移动至同一个文件夹，构成一个声发射事件的信号集(s1,s2,...,sn)，n为声发射信号的数目。然后按照声发射信号触发时间先后排序，排序后的信号集记为(x1,x2,...,xn)。进而对(x1,x2,...,xn)展开加窗补长等操作，得到的加窗信号集记为(u1,u2,...,un)。最后采用互相关法计算每个加窗声发射信号相对第一个加窗声发射信号的p波走时差(t1,t2,...,tn)。

步骤2：变网格搜索虚拟场法声发射事件定位；

以步骤1中的p波走时差(t1,t2,...,tn)作为定位数据，采用变网格搜索计算《一种用于均匀速度场的信号源定位方法》中的虚拟场势值，取势值最大的网格点作为声发射事件的定位结果。

步骤3：断铅实验计算每个传感器的感度系数；

在10cm×10cm×10cm的方形岩石试样前后左右四个面各布置两个传感器，传感器坐标为(5,10,8.5)，(0,5,8.5)，(5,0,8.5)，(10,5,8.5)，(5,10,1.5)，(0,5,1.5)，(5,0,1.5)，(10,5,1.5)cm。在顶面划分1cm×1cm的网格进行断铅实验，并计算每次断铅实验时各个传感器的感度系数，取各个传感器多次断铅的感度系数均值作为每个传感器的感度系数(sc1,sc2,...,scn)。

步骤4：计算修正后的声发射信号最大振幅；

计算每个声发射信号的最大振幅(a1,a2,...,an)，将每个传感器最大振幅除以其感度系数，得到修正后的最大振幅

步骤5：依据修正后最大振幅计算声发射震级；

采用公式计算第i个声发射信号的震级，取所有震级的中位数作为该事件的最终震级m。其中，的单位为微伏，k为震级调整参数。

实施例

图2是声发射监测时传感器布置图，其中黑圆点为可见面上的传感器，灰圆点为不可见面上的传感器。其特征在于，在10cm×10cm×10cm的方形岩石试样前后左右四个面各布置两个传感器，传感器坐标为(5,10,8.5)，(0,5,8.5)，(5,0,8.5)，(10,5,8.5)，(5,10,1.5)，(0,5,1.5)，(5,0,1.5)，(10,5,1.5)cm。

图3是加窗互相关拾取实例图。图3(a)为声发射信号全波形，垂直竖线为手动拾取p波初至，声发射信号1的p波初至较为明显，而声发射信号2的p波初至不明显，手动拾取p波初至可靠性较低。另一方面，p波初至初始阶段声发射信号波形较为相似，而后半段由于衰减等影响，波形差异性很大，常规的全波形互相关确定相对到时误差较大。

由此，提出了一种加窗互相关技术拾取p波相对到时差。其特征在于，首先大致拾取一个p波初至前的时间t1(图3a中浅灰色框的左边界)，再大致拾取首个p波峰值后的时间t2(图3a中浅灰色框的右边界)，进而自动延伸至下一个振幅为0的时间t3(图3a中灰色框的右边界)，最终加窗信号为区间[t3-max(t3-t1),t3]的时间序列(图3b中带框部分)。其中，t3-max(t3-t1)用于统一加窗信号的长度，max(t3-t1)为所有加窗信号中t3-t1的最大值。加窗后信号如图3(b)所示。

进而对加窗信号与第一个加窗信号展开互相关计算，得到加窗信号互相关序列图3(c)。加窗互相关最大值对应声发射信号相对到时差。加窗互相关序列存在明显的最大值，表明加窗互相关拾取相对到时差具有很好的效果。

图4是变网格虚拟场定位实例图，图中圆圈表示搜索网格点，背景色为虚拟场势值云图。将触发时间差小于40us的声发射信号，构成一个声发射事件的信号集(s1,s2,...,sn)，n为声发射信号的数目。然后按照声发射信号触发时间先后排序，排序后的信号集记为(x1,x2,...,xn)。进而对(x1,x2,...,xn)展开加窗补长等操作，得到的加窗信号集记为(u1,u2,...,un)。最后采用互相关法计算每个加窗声发射信号相对第一个加窗声发射信号的p波走时差(t1,t2,...,tn)。

以p波走时差(t1,t2,...,tn)作为定位数据，采用变网格搜索计算虚拟场势值。变网格的特征在于，首先采用3mm×3mm×3mm的网格计算声发射事件的大致位置，再以上述位置为中心生成一个3mm×3mm×3mm的立方体，利用0.5mm×0.5mm×0.5mm的细网格搜索上述立方体，见图4中的圆圈点。背景色为网格点的虚拟场势值云图，可知声发射定位收敛性很好(左下角整体云图)，对其细网格附近区域放大，可知声发射定位为(45,36)mm。

图5是声发射震级计算实例图。首先基于断铅实验计算每个传感器的感度系数，计算方法为sc＝ap·d/f(n,d)。其中，ap为p波初至后的首个峰值振幅，d为声发射事件与传感器的距离，d为声发射事件与传感器连线方向的单位向量，n为传感器接收信号方向的单位向量，r为p、s波的速度比值，s＝n·d。

进而计算每个声发射信号的最大振幅(a1,a2,...,an)，将每个传感器最大振幅除以其感度系数，得到修正后的最大振幅采用公式计算第i个声发射信号的震级。其中，的单位为微伏，k为震级调整参数。

最大振幅震级和修正后最大振幅震级如图5所示，其特征参数如表1所示。可知基于最大振幅的震级差异性较大，最大值和最小值分别为6.30和4.61，标准差为0.62，震级差值达到了1.69。而修改后最大振幅的震级差异性较小，最大值和最小值分别为5.93和5.54，标准差为0.11，震级差值仅为0.39。修正后最大振幅震级的中位数作为该事件的最终震级m，即该声发射事件的震级为5.72。

表1声发射震级特征参数

一种考虑传感器感度的声发射事件震级获取系统，包括：

信号加窗补长单元：将采集得到的声发射事件的信号集，依据信号时间差，将声发射信号划分至不同的声发射事件的信号集中，对每一个声发射事件的所有声发射信号按照触发先后排序，对排序后的声发射信号进行加窗补长至振幅为0的点，获得各声发射事件的加窗声发射信号集；

信号p波走时差计算单元：基于互相关法计算各声发射事件的加窗声发射信号集中，每个加窗声发射信号相对第一个加窗声发射信号的p波走时差(t1,t2,...,tn)，n为声发射信号的数目；

声发射事件定位单元：基于p波走时差(t1,t2,...,tn)，采用变网格搜索法计算声发射事件所在空间的各网格点虚拟场势值，取虚拟场势值最大的网格点作为声发射事件的定位结果；

所述虚拟场势值采用已公开的专利一种用于均匀速度场的信号源定位方法进行确定；

声发射信号最大振幅修正单元：利用传感器的感度系数，对声发射事件的所有原始声发射信号的最大振幅进行修正；

当一个声发射事件发生后，其通过波的形式传到传感器，传感器将其转换为电压信号，即为声发射信号；

声发射事件震级确定单元：依据修正后的声发射信号的最大振幅，确定每个声发射信号的震级，以所有声发射信号的震级的中位数作为声发射事件的震级；

其中，mi表示第i个加窗声发射信号的震级，为第i个加窗声发射信号的修正后的最大振幅，单位为微伏，k为震级调整参数，取值范围为[-1,0]0。

应当理解，本发明各个实施例中的功能单元模块可以集中在一个处理单元中，也可以是各个单元模块单独物理存在，也可以是两个或两个以上的单元模块集成在一个单元模块中，可以采用硬件或软件的形式来实现。

本发明实施例还提供一种计算机可读存储介质，包括计算机程序指令，所述计算机程序指令被处理终端执行时使所述处理终端执行上述的一种考虑传感器感度的声发射事件震级获取方法。其有益效果参见方法部分的有益效果，在此不再赘述。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人员应当理解：依然可以对本发明的具体实施方式进行修改或者等同替换，而未脱离本发明精神和范围的任何修改或者等同替换，其均应涵盖在本发明的权利要求保护范围之内。