刘飞跃, 杨天鸿, 邓文学, 张鹏海
东北大学 资源与土木工程学院, 辽宁 沈阳 110819
收稿日期：2017-08-23
基金项目：国家自然科学基金资助项目(51604062)；国家重点基础研究发展计划项目(2016YFC0801602)。
作者简介：刘飞跃(1996-)，男，安徽阜阳人，东北大学博士研究生;
杨天鸿(1968-)，男，辽宁抚顺人，东北大学教授，博士生导师。

摘要：岩石破裂过程中的声发射(AE)监测技术已广泛应用于岩石工程领域, 声发射信号的到时提取是深入开展岩石力学声发射研究的基础和关键.为了获取更为准确的声发射信号P波到时, 在前人研究的基础上, 提出了一种改进的Pphase-Picker声发射P波到时自动提取算法并对其原理进行详细描述.选取1 000个典型声发射波形, 对比分析了几种到时自动提取算法的计算结果, 表明改进的算法不仅可以获取更为精确的P波到时, 还可以获取到时提取误差的估计值.
关键词：P波到时提取Pphase-Picker算法置信度误差估计声发射
An Improved Automatic Algorithm of Acoustic Emission P-Wave Arrival Time Picking
LIU Fei-yue, YANG Tian-hong, DENG Wen-xue, ZHANG Peng-hai
School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China
Corresponding author: YANG Tian-hong, E-mail: yangtianhong@mail.neu.edu.cn
Abstract: Acoustic emission(AE)monitoring technique has been widely applied in the area of geotechnical engineering. Automated arrival time picking of P-wave is the fundament and key in the further rock mechanical AE study. In order to get the more accurate arrival time of AE single, an improved Pphase-Picker P-wave arrival time picking algorithm was proposed and its principle was described in detail. Several automatic arrival time picking algorithm were compared in 1 000 typical AE wave, the result indicated that the improved algorithm can not only get a more accurate P-wave arrival time rapidly, but also can get an error estimation of P-wave arrival time.
Key words: P-wave arrival time pickingPphase-Picker algorithmcertain degreeerror estimationacoustic emission
声发射技术是研究岩石力学性质的最有效工具之一, 该方法中关键的一步就是P波到时的准确提取.声发射监测过程中记录得到的海量波形数据文件远远超出人工所能处理的范围, 国内外学者提出多种到时自动提取方法.其中有以下两种方法应用最为广泛:1)长短时均值比方法(STA/LTA法).Stevenson^[1]首次提出了用短时窗信号平均值(STA)和长时窗信号平均值(LTA)之比反映信号振幅或能量的变化, 由此确定到时; 在此基础上Allen^[2]提出了用特征函数的方式来灵敏反映并增强初动到达时其振幅和频率变化特征, 以此更加精确确定到时; Munro^[3]提出了多窗口算法, 克服了STA/LTA算法受STA长度的影响而不能准确给出到时以及难以在高振幅噪音中检测出初动的缺陷.2)时间序列自回归模型法(AR-AIC法).Ozaki等^[4]假定波形数据在初动到达前和到达后各自为一定常态, 用自回归模型模拟地震波, 将振幅和频率的信息反映在自回归模型的系数上并根据模型的变化来检测初动; Maeda^[5]确定了AR-AIC算法的置信度来衡量到时提取的准确性; Carpinteri等^[6]利用设置时间窗口和小波滤波来提高到时提取的计算速度和准确度.此外, 基于波形偏斜度和峰度的PAI-S/K方法^[7]、从时频域分解出发的ITD算法^[8]、神经网络法^[9]、分形分维方法^[10]等方法在声发射P波到时自动提取上也有应用.
上述研究方法各有其特点, 从不同角度丰富和改进了到时提取方法, 并在实际应用中获得了较为良好的效果.但是上述算法多是从统计学角度出发, 根据初动到达前后波形的差别构建特征函数加以判断, 缺乏明确的物理学意义, 从声发射力学机理上着手研究的算法较少.对此，Kalkan^[11]提出了Pphase-Picker算法, 该算法将质点的振动简化成一个单自由度阻尼振荡器, 利用阻尼能在初动前后的差异, 实现P波到时的自动提取, 因此Pphase-Picker算法具有对各种波形很强的适应性.但在实际应用中, Pphase-Picker算法存在部分参数选取困难、精度不稳定等问题.为解决上述问题, 本文从多角度出发, 对Pphase-Picker算法进行了改进.选取1 000个典型声发射波形, 将多种到时提取算法分别应用到声发射波形到时提取上并进行比较, 结果表明改进的Pphase-Picker算法不仅具有更好的计算精度, 而且可以获取到时提取误差的估计值, 验证了改进算法的优越性.
1 Pphase-Picker到时提取方法原理岩石在外部载荷作用下能量快速释放引起质点振动并产生弹性波向外传播的现象称为声发射.如图 1所示, Kalkan用单自由度阻尼振荡器模拟质点的振动, 图 1a为移动支座单自由度振荡器, 图 1b为固定支座单自由度振荡器.对于声发射波形, 质点振动的方式与图 1b中固定支座单自由度振荡器一致, 其阻尼能计算公式可由移动支座单自由度振荡器变形推导而来, Pphase-Picker算法的原理简要叙述如下.
图 1(Fig. 1)

图 1 理想状态下黏性阻尼单自由度振荡器数学模型Fig.1 Mathematical models of idealized single- degree-of-freedom oscillators with viscous damping (a)—绝对能量模型；(b)—相对能量模型.

对于两种不同的单自由度振荡器, 考虑到力平衡, 均有

(1)

式中:m为质量；c为阻尼项；k为刚度项；u_t(u_t=u+u_g)为绝对位移, u_g为地面位移；u为相对位移.对于图 1b中固定支座单自由度振荡器, 式(1)可改写为

(2)

将式(2)等号两侧分别对u进行积分, 可得

(3)

式中t为时间.可将式(3)写成简便形式, 表示不同的能量项, 有

(4)

式中:E′_K为相对动能；E_ζ为阻尼能；E_S为应变能; 它们的和为相对输入能量E′_I.
在固定支座阻尼单自由度振荡器中, 阻尼能对时间的导数(阻尼能率)可以表示为

(5)

式中:ζ为阻尼率；ω_D为圆频率, 定义为2π/T_D, T_D是阻尼振动的自然周期, 与无阻尼振动的自然周期T_n的关系为

(6)

如图 2所示, 当初动尚未到达时, 阻尼能维持在一个较低水平, 初动到达后, 阻尼能急剧上升, 基于此原理对到时进行提取.在求解到时的过程中, 应用直方图法.直方图法求解到时原理如图 2所示.
图 2(Fig. 2)

图 2 直方图法求解到时示意图Fig.2 Illustration of arrival time solved by histogram method

2 改进的Pphase-Picker算法在实际应用中, 考虑多因素对声发射波形造成的影响, 为获取更佳的计算精度, 从多角度对Pphase-Picker算法进行改进.
2.1 声发射波形降噪在实际监测过程中, 采集到的声发射波形必然会受到各种噪声的干扰, 对到时提取造成影响.通常情况下, 噪声与声发射信号在频率和振幅上有明显区别, 因此通过数字滤波可以有效消除噪音干扰.声发射信号频带较宽, 且频率成分丰富, 为更好保存有用信号, 采用通带平缓的Butterworth无限脉冲响应数字滤波器(IIR)带通滤波器, 其幅度特性可以表示为

(7)

式中:H_a(jΩ)为滤波器单位脉冲响应；Ω_c为幅值衰减到-3dB时的频率；N为滤波器的阶数.
已有较多学者对声发射信号频谱进行研究, 结果表明声发射信号主要集中在1~100 kHz, 本文中设置高通转角频率为10 kHz, 截止频率为1 kHz, 低通转角频率200 kHz, 截止频率500 kHz, 截止频率处幅值至少衰减到-30 dB, 确定滤波器阶数为4阶.
图 3给出了一组原始含噪波形和Butterworth带通滤波器去噪后的声发射波形, 可以看出经过滤波后信噪比有了显著的提高.利用该滤波方法对声发射波形进行预处理, 可以有效提高波形的信噪比, 达到很好的信噪分离效果.
图 3(Fig. 3)

图 3 Butterworth滤波Fig.3 Butterworth filtering

2.2 时间窗口的设定对于一个声发射波形序列, 波形两端存在着较长的背景噪音段, 噪音段的存在使得计算速度减慢, 且波前噪音对到时提取的精确度有影响.通过设置时间窗口选取初动前典型噪音与初动后典型信号波形段来解决上述问题, 为此需首先确定到时的估计值.为减少人为参数设定, 提高到时提取的自动化程度, 参照文献[6], 当某时刻后10个数据点的绝对平均值首次大于先前所有数据点绝对平均值的4倍时, 认定该时刻为波形到时的估计值k₀.

(8)

时间窗口长度过长, 参与到时寻找的数据点过多, 导致计算速度减慢; 时间窗口长度过短, 初动前后波形特性不能充分显现, 影响到时提取精度.由于声发射波形的持续时间为10²~10³ μs级别, 考虑到声发射监测系统采样频率8 Msps(million samples per second, 简称Msps, 表示每秒采样百万次)、采样长度5 k(表示每个声发射波形采集数据点个数为5 120)与预触发160(表示触发时刻来到之前, 所采集的时长为160 μs)等参数的设定, 在本文中将时间窗口长度设定为1 001, 步长为1个采样点时间间隔, 即参与到时寻找的数据点范围为[k₀-500, k₀ + 500], 图 3中红色线框内即为设定的时间窗口.
2.3 置信度的确定应用直方图法求取P波到时, 为衡量到时提取的准确性, 引入置信度的概念.由于定义P波到时为低阻尼能水平至高阻尼能水平的起跳点, 对于同一波形, 到时前后阻尼能大小差异越明显, 即E_high/E_low比值越大, 到时提取结果越准确.考虑到不同波形振幅上存在差异, 将高子直方图平均阻尼能与低子直方图平均阻尼能的比值R_high/R_low作为基准, 定义置信度DD:

(9)

其中：δk为置信度计算窗口长度, 本文选取30个采样点长度；E_low和E_high分别表示在[E(k－δk), E(k)]和[E(k), E(k+δk)]上的阻尼能平均值.从下文室内试验中选取1 000个典型声发射波形, 定义到时提取误差为自动算法计算到时与人工提取到时之间的差值, 以人工提取到时误差为零, 计算得到置信度与到时提取误差的关系如图 4所示, 图中数据点的颜色表示数据的集中程度.可以看出置信度与到时提取误差具有良好的相关性, 置信度越高, 到时提取的精确性越好, 利用置信度可以有效衡量到时提取的准确性.
图 4(Fig. 4)

图 4 置信度与到时提取误差Fig.4 Confidence and arrival time picking error (a)—单一直方图元；(b)—优化直方图元.

置信度可为直方图法求解到时提供参数选取的依据, 同时也可以为到时提取误差估算提供依据.
2.4 直方图元数的确定在直方图法确定到时的过程中, 直方图元数对到时提取的精确性有着很大的影响, 对于不同波形, 其平均阻尼能差异较大, 给定一个固定直方图元数是不可取的.经验表明, 直方图元数应该有一个上限与下限, 直方图元数应该足够大以达到足够的分辨率; 直方图元数也不应该太大以至于单元宽度Δy小于幅度量化位数.为此考虑将直方图元数进行迭代优化以期获取误差的到时提取结果.随着直方图元数的增大, 到时提取误差为先下降后上升的阶跃函数, 置信度为先上升后下降的单谷函数, 且当置信度最高时对应的到时提取误差最小, 对应的直方图元数为最佳直方图元数.据此原理, 将求解最佳直方图元数的实际问题简化为求单谷函数极值问题的数学问题.在本文中采用粒子群方法求解函数极值, 该方法具有搜索速度快、需调整参数少、结构简单易于实现等优点, 其具体原理本文中不赘述.
如图 4所示, 采用优化直方图元数较单一直方图元数置信度有了明显上升, 对应的到时提取误差也有明显降低.采用粒子群方法优化得到的直方图元数可以有效降低到时提取误差.
2.5 到时提取误差的获取使用相关性数来考查1 000个典型声发射波形到时提取获取的置信度DD与到时提取误差Δt之间的相关性.相关性系数的计算公式如式(10)所示, 其中Cov(DD, Δt)为置信度、到时提取误差之间的协方差; D(DD)，D(Δt)分别为置信度、到时提取误差的方差.计算得到置信度与到时提取误差的相关性系数达0.851, 具有良好的相关性.

(10)

等间距选取到时提取误差-置信度的25个均值点进行曲线拟合, 如图 4b所示, 置信度与到时提取误差之间的定量关系近似为负指数函数, 拟合优度达0.967.如式(11)所示, 可以通过到时提取算法得到的置信度DD计算到时提取误差Δt.

(11)

值得注意的是, 该方法获取的到时提取误差Δt与真实到时提取误差存在差异, 只能作为真实到时提取误差的估计值.但该方法获取的到时提取误差在无到时提取误差实测情况下, 可一定程度上定量反映到时提取的准确性, 也可为源定位误差分析等进一步的声发射研究提供依据.
3 到时提取算法的验证3.1 声发射波形的获取岩石试件为边长100 mm的立方体花岗岩, 岩石试件均匀致密, 无明显缺陷, 如图 5a所示.声发射系统的6个传感器布置方式如图 5b所示, 用胶皮带固定至岩石表面并以凡士林为耦合剂增加传感器与岩石表面的耦合效果, 自检声学特性幅值均能达到80(最大100)以上, 耦合效果较好.试验采用加载速率为0.05 mm/min的单轴位移加载方式, 人工选取1 000个典型声发射波形.
图 5(Fig. 5)

图 5 声发射试验Fig.5 Acoustic emission test (a)—试件照片；(b)—传感器布置.

3.2 多种到时提取算法的对比改进的Pphase-Picker算法由质点的加速度计算阻尼能, 而试验中采用的Nano30传感器为速度传感器, 其电压幅值V_s正相关于接收点速度v, 即V_s∝v, 因此首先需要将原始波形信号电压幅值转换为速度, 在声发射监测频带范围内传感器灵敏度S=62 dB, 单位:V/(m·s^-1).可以通过式(12)得到测点的速度波形, 再对速度波形进行一次数值求导, 得到加速度波形, 然后进行到时提取.

(12)

为比较多种到时提取算法的精度, 运用Matlab软件, 将长短时均值比方法^[2]、时间序列自回归模型法^[6], Pphase-Picker算法^[11], 改进Pphase-Picker算法分别对声发射试验获取的1 000个典型波形进行处理, 考察各算法的计算速度与计算精度, 将各算法计算得到P波到时与人工提取结果进行对比得到误差.如表 1所示, 改进的Pphase-Picker算法较其他算法精度高, 误差分别为STA/LTA算法的20.1%, AR-AIC算法的53.5%, Pphase-Picker算法的55.2%.可以看出本文提出的改进Pphase-Picker算法获取的误差估计与实际较为接近, 是一种优良的声发射P波到时提取方法.
表 1(Table 1)

表 1 多种到时提取算法计算结果Table 1 Results obtained by different arrival time picking algorithms μs 误差 STA/LTA算法 AR-AIC算法 Pphase-Picker算法 改进Pphase-Picker算法 误差估计 最大误差 8.35 3.80 5.48 2.96 3.07 平均误差 4.22 1.59 1.54 0.85 0.92

表 1 多种到时提取算法计算结果 Table 1 Results obtained by different arrival time picking algorithms

4 结论1) 本文结合Butterworth带通滤波、时间窗口设定、粒子群方法优化直方图元数以及置信度与到时提取误差之间的定量关系, 对原有Pphase-Picker声发射P波到时自动提取算法进行了改进, 并且通过不同算法对比验证了改进Pphase-Picker算法的优越性.
2) 与传统的STA/LTA算法、AR-AIC算法以及原有Pphase-Picker算法相比, 改进的Pphase-Picker算法显著提高了到时提取精度.
3) 改进的Pphase-Picker算法可以获得对到时提取误差的估计值, 为无误差实测条件下误差的量化提出了初步解决方法.
参考文献

[1]	Stevenson R. Micro earthquakes at Flathead Lake, Montana:A study using automatic earthquake processing[J].Bulletin of the Seismological Society of America, 1976, 66(1): 61–80.
[2]	Allen R. Automatic phase pickers:their present use and future prospects[J].Bulletin of the Seismological Society of America, 1982, 72(6): 225–242.
[3]	Munro K.Automatic event detection and picking of P-wave arrivals[R].Calgary: CREWES Research Report, 2004, 16: 1-10.
[4]	Ozaki T, Tong H. On fitting of non-stationary autoregressive models in time series analysis[J].Japanese Journal of Industrial Health, 1975, 34(7): 642–647.
[5]	Maeda N. A method for reading and checking phase times in auto processing system of seismic wave data[J].Journal of the Seismological Society of Japan, 1985, 38(2): 365–379.
[6]	Carpinteri A, Xu J, Lacidogna G, et al. Reliable onset time determination and source location of acoustic emissions in concrete structures[J].Cement & Concrete Composites, 2012, 34: 529–537.
[7]	Saragiotis C D. PAI-S/K:A robust automatic seismic P phase arrival identification scheme[J].IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(6): 1395–1404.DOI:10.1109/TGRS.2002.800438
[8]	Zhang R H, Zhang L H. Method for identifying micro-seismic P-arrival by time-frequency analysis using intrinsic time-scale decomposition[J].Acta Geophysica, 2015, 63(2): 468–485.DOI:10.1515/acgeo-2015-0007
[9]	Gentili S, Bragato P. A neural-tree-based system or automatic location of earthquakes in northeastern Italy[J].Journal of Seismology, 2006, 10(1): 73–89.DOI:10.1007/s10950-005-9001-z
[10]	Boschetti F, Dentirh M, List R. A fractal-based algorithm for detecting first arrivals on seismic traces[J].Geophysics, 1996, 61(4): 1095–1102.DOI:10.1190/1.1444030
[11]	Kalkan E. An automatic P-phase arrival-time picker[J].Bulletin of the Seismological Society of America, 2016, 106(3): 971–986.DOI:10.1785/0120150111