在岩石钻孔爆破过程中，爆源近区岩石的冲击振动信号具有持时短、频带宽的特点^[1]。目前广泛使用的压电式加速度传感器在采集此类信号时，大多会出现“零漂”现象^[2-6]。造成压电传感器出现零漂误差的原因主要有^[2-3]：(1)受到高频冲击波作用时，传感器由于自身谐振反应导致内部压电晶体磁畴变化，这种变化反映到输出信号上，出现零漂现象；(2)高冲击作用下，传感器压电元件受到的作用力超过其在装配时受到的预压力，进而造成压电晶体与惯性质量块之间出现微量运动，造成零漂现象；(3)传感器的质量-弹簧系统受到高冲击作用时的动刚度不足会导致其幅值线性度变差，可能造成零漂；(4)受限于传感器压电材料和制造工艺，信号放大器电路易过早进入饱和堵塞状态，RC电路放电时间延长到冲击结束后，表现为零漂；(5)传感器基座与被测物体表面连接刚度不足造成零漂误差。由于爆破过程中岩石结构会产生不可逆的损伤，因而每次爆破信号的采集是单次性的。所以，如果爆破振动信号采集过程中出现零漂现象，采用合理方法对误差信号进行修正非常有必要^[4-6]。

近年来，不少学者对于零漂信号中产生的趋势项去除方法进行研究：龙源等^[6]基于经验模态分解(empirical mode decomposition, EMD)、小波法、最小二乘法分别研究了爆破震动测试信号中趋势项的去除算法，并对三种方法去除趋势项效果及重构信号的时频特征进行对比分析，指出多项式趋势项中的多项式阶数选取对于信号的修正效果有较大影响；王燕等^[7]采用EMD法和最小二乘法相结合，以侵彻结束后的零漂趋势作为侵彻过程中的零漂趋势，对侵彻不同介质的实测过载数据进行去除零漂的处理，通过对修正后的过载数据进行一次积分和二次积分，得到的侵彻速度和位移与实测结果吻合较好；王济^[8]研究发现滑动平均法应用较为简便，但其本质仅是为了曲线基线回归零点，在处理信号时有一定的盲目性；谢全民等^[9]认为小波法需要针对振动信号特点事先选取合适的小波基函数和分解层数；吴祖堂等^[10]通过理论分析和实验验证得出冲击结束后的压电式传感器信号呈指数形式衰减的结论，但是信号的零漂往往在冲击信号达到峰值前就已经形成，冲击结束后信号的指数衰减规律和冲击过程中的趋势项类型是否一致，目前尚无定论；Huang等人^[11]揭示EMD方法及其改进算法—集合经验模态分解（ensemble empirical mode decomposition, EEMD）方法的基函数都来自于其信号本身，对于信号的先验性信息要求较低，在处理爆破冲击这类非平稳信号时优势明显。实际上，EEMD方法在EMD方法基础上，通过添加在频谱上均匀分布的白噪声有效地解决了EMD方法分解得到的固有模态分量(intrinsic mode function, IMF)中存在的模态混叠和端点震荡问题^[12]，故认为EEMD方法处理零漂信号的效果优于EMD方法。但EEMD方法中需要给出添加白噪声的系数，目前没有确定的量化指标去评价零漂信号的修正效果，从而无法确定白噪声系数的合理取值范围。

基于此，本文先给出不同白噪声系数下EEMD结合高低频处理方法修正岩石爆破中加速度零漂信号的过程；再辅以冲击响应谱分析，提出表征修正前后频域平均偏差幅度的修正指数$\bar D$；在修正指数分析的基础上确定EEMD修正爆破零漂信号中白噪声系数的最优取值范围，再用加速度信号积分后的速度信号对比讨论最优噪声系数范围的合理性，最后对比EEMD法与滑动平均法，讨论二者的优点与不足。

1.   算法原理

EEMD是针对EMD算法中所存在的不足进行改进的算法，EMD算法实现过程如图1所示^[11]。EEMD主要利用高斯白噪声在频谱上均匀分布的特性来弥补原始信号中一些缺失的尺度，使得不同时间尺度的信号自动分解到其相适应的参考尺度上去，从而有效地抑制EMD分解信号时存在的模态混叠问题。由于添加的白噪声总体均值为零，所以将多次添加白噪声后所求得的IMF分量的集合平均值作为最终结果，以抵消白噪声对于真实IMF分量的影响。

图  1  EMD算法流程图

Figure  1.  Flow chart of EMD algorithm

下载: 全尺寸图片 幻灯片

EEMD算法的实现步骤如下^[12]。

（1） 将随机高斯白噪声信号n_j (t)添加到原始信号x(t)：

式中：x_j (t)是第j次加噪信号，k、m分别是添加白噪声的系数和总次数，σ是原始信号的标准差。

（2） 用EMD方法将x_j (t)分解成一系列IMF分量c_i,j 和余项δ_j :

式中：c_i,j和δ_j分别为第j次加入白噪声后信号分解得到的第i个IMF分量和余项，g_j为第j次加入白噪声后信号的分解尺度。

（3） 如果j＜m，则重复步骤（1）和（2）。

（4） 获取I=min(g₁, g₂,···, g_m)，计算I相应的分解的IMF的集合平均值作为最终结果:

2.   信号处理

2.1   信号介绍

本文所处理的零漂信号来源于花岗岩爆破试验，其数据采集过程和时程曲线分别如图2和图3所示。爆源为1.5 g太安炸药，爆速约为6 000 m·s⁻¹。图3中加速度信号A测点的爆心距约为60 cm。信号放大器的采样频率为5×10⁵ Hz，加速度传感器的频响范围是40～4×10⁴ Hz，设置低通3×10⁴ Hz。由图3可以看出峰值加速度超过了4×10⁴ m·s⁻²，同时加速度信号在冲击结束后没有迅速归零，而是持续一段时间后才缓慢归零，这是高速冲击环境下的典型零漂信号。

图  2  试验与数据采集

Figure  2.  Test and data collection

下载: 全尺寸图片 幻灯片

图  3  爆破加速度信号A

Figure  3.  Blast acceleration signal A

下载: 全尺寸图片 幻灯片

2.2   白噪声系数范围

Wu等^[12]在分析EEMD理论时，建议k近似取0.20；Jiang等^[13]用EEMD方法对振动信号分析时指出，如果信号主频段由高频分量组成，则k的取值范围建议为0～0.20。

岩石的爆破试验结果表明^[14-15]，岩石质点的振动主频率与炸药装药量成反比关系，与炸药爆速成正比关系；装药深度对振动主频率也有很大影响，深孔爆破产生的地震波主频较低，浅孔爆破产生的地震波主频较高。本文爆破试验的装药量小、深度浅、炸药爆速大且传播距离短，故认为其主频段主要由高频分量组成。因此，白噪声系数k选在0～0.20范围内讨论。

2.3   修正过程

2.3.1   EEMD分解

为分析白噪声系数对于EEMD方法修正零漂信号的影响，在0～0.20范围内选取0.05、0.10、0.15和0.20四个k值对应的白噪声添加到原始信号中进行EEMD分解，添加次数均为100次。图4给出k=0.10时分解得到的IMF分量和余项δ。

图  4  EEMD分解信号(k=0.10)

Figure  4.  Signals decomposed by EEMD method (k=0.10)

下载: 全尺寸图片 幻灯片

2.3.2   IMF分量频谱分析

将分解得到的IMF分量分别进行FFT变换，可得各IMF分量的频谱信息。不同白噪声系数下分解得到的各IMF分量的主频如表1所示。限于篇幅，仅给出部分IMF分量的频谱图，如图5所示。

表  1  不同白噪声系数下IMF分量主频

Table  1.  Dominant frequencies of IMF components under different white noise coefficients

IMF分量	主频/Hz					IMF分量	主频/Hz
	k=0.05	k=0.10	k=0.15	k=0.20			k=0.05	k=0.10	k=0.15	k=0.20
IMF1	18 860	18 740	18 760	18 760		IMF8	440	540	600	640
IMF2	18 880	18 800	18 840	1 800		IMF9	240	280	300	320
IMF3	7 500	8 420	8 480	8 500		IMF10	100	200	220	220
IMF4	5 100	7 520	7 520	7 520		IMF11	40	40	40	40
IMF5	2 280	3 180	5 160	5 160		IMF12	20	20	20	20
IMF6	1 240	1 360	1 560	2 340		IMF13	20	20	20	20
IMF7	660	740	800	840		δ	20	20	20	20

下载: 导出CSV 

| 显示表格

图  5  部分IMF分量频谱图

Figure  5.  Spectrograms of partial IMF components

下载: 全尺寸图片 幻灯片

从表1可以看出，白噪声系数对前两阶IMF分量的主频影响很小，前两阶IMF分量主频均在18 800 Hz左右。结合图5(a)可知，白噪声系数主要影响其主频幅值，系数越大，主频幅值越小；从IMF2至IMF11分量，随着分解阶数的增加，IMF分量主频迅速衰减，且不同白系数噪声下分解的同阶IMF分量主频与白噪声系数呈现正相关关系；IMF11之后的分量的主频降至低频段，且基本保持不变。同时，从图5可看出，相同白噪声系数下分解得到的IMF分量随着分解阶数增加，其频段降低，频带渐窄。

2.3.3   IMF分量处理

（1）低频处理。由表1知，IMF11至IMF13分量和余项δ的主频集中在40 Hz以下，低于试验中加速度传感器的频响下限 (40 Hz)，故认为其可信度很低。因此，将IMF11至IMF13分量和余项δ消除。

（2）高频处理。由于原始信号就含有大量高频噪声，故对不同白噪声系数下的IMF1至IMF10分量分别进行小波阈值去噪处理，阈值由heursure函数自适应获取，去噪方式采用软阈值处理^[16-17]。

分别对处理过的分量IMF1至IMF10进行重构，即可得到不同白噪声系数下EEMD修正后的爆破冲击信号。信号A修正后的四个加速度信号都很好地消除了零漂趋势，但由于彼此在加速度时程图上的波形差异很小导致难以区分，故仅给出白噪声系数为0.10时的修正信号，如图6所示。

图  6  EEMD修正前后的信号

Figure  6.  Signals before and after EEMD correction

下载: 全尺寸图片 幻灯片

3.   修正信号分析

3.1   频域分析

图7给出了修正前后信号的冲击响应谱。从冲击响应谱上来看，四个修正的信号与原始信号A的频域“差异”主要集中在4 000 Hz以下。在4 000 Hz以上的频段，修正前后的频域“差异”不明显，很好的保留了原始信号的峰值和变化规律。需要指出，这里的“差异”是指处理后的信号与原始信号在冲击响应谱幅值上的相对偏差度，而非绝对偏差值。由于40 Hz以下频段已经超出加速度传感器的频率响应范围，这部分原始信号认为是不可信的，采用的修正方法应将这部分的信号尽可能削弱，同时应该尽量保留中高频段可信信号的频谱特征。因此，为兼顾防止信号中高频段过度修正和尽量削弱低频失信段的目的，需要定义一个具体指标来评价信号的修正效果。

图  7  不同白噪声系数下EEMD修正信号的冲击响应谱

Figure  7.  Shock response spectra of signals modified by EEMD under different white noise coefficients

下载: 全尺寸图片 幻灯片

3.1.1   修正指数定义

设所选频段的频率序列为{f₁, f₂,···, f_l}，对应的冲击响应谱序列为{P₁, P₂, ···, P_l}，其中f为频率，P为每个频率对应的冲击响应谱幅值，l为频段的频率点总数。则反映信号修正前后冲击响应谱平均偏差幅度的修正指数$\bar D$定义如下：

式中：${P_i}$和$P_i'$分别是原始信号和修正信号在所选频段上的第i个频率点的冲击响应谱幅值。

3.1.2   指数应用

根据加速度传感器频响范围和图7中的加速度冲击响应谱，从整个频段中选出几个特殊频率，分别是传感器频响下限40 Hz、修正前后响应谱上存在明显差异频段的上限频率4×10³ Hz和响应谱峰值对应的频率1.9×10⁴ Hz。先将整个频段分为2个分段：失信频段0～40 Hz、置信频段40～3×10⁴ Hz，再将置信频段分为三个分段：40～4×10³ Hz、4×10³～1.9×10⁴ Hz和1.9×10⁴～3×10⁴ Hz，将这五个频段中的修正指数$\bar D$分别记为$\bar D_1$、$\bar D_2$、$\bar D_3$、$\bar D_4$和$\bar D_5$。

为了评价不同白噪声系数下信号在不同频段上的修正效果，在(0, 0.20]中以0.001为间隔取200个k值进行信号处理。原始信号A在这些白噪声系数下处理后，可得200个修正信号。然后，将这些修正信号与原始信号的冲击响应谱代入式(5)计算，即可得到不同白系数噪声下相应频段的修正指数$\bar D$，如图8所示。在给出的五个频段上，修正指数与白噪声系数都呈现近似幂指数的关系，可假定修正指数与白噪声系数的函数关系为：

图  8  不同白系数噪声下不同频段上的冲击响应谱修正指数

Figure  8.  Correction indices of shock response spectra on different frequency bands under different white noise coefficients

下载: 全尺寸图片 幻灯片

式中：a₀、a₁、和b为常数。

不同频段的拟合曲线及拟合函数如图8所示。拟合曲线的相关系数R²均大于0.95，拟合效果良好，可以认为修正指数与白噪声系数符合式(5)的幂指数关系。根据不同频段上修正指数及其拟合函数图，可得不同频段上修正指数$\bar D$的取值范围，如表2所示。其中，修正指数下限${\bar D_{\min }} = ({a_0} + b){\text{%}}$，修正指数上限${\bar D_{\max }} = b{\text{%}}$。

表  2  不同频段修正指数范围

Table  2.  Ranges of correction index on different frequency bands

频段/Hz	${\bar D_{\min }}/{\text{%}}$	${\bar D_{\max }}/{\text{%}}$		频段/Hz	${\bar D_{\min }}/{\text{%}}$	${\bar D_{\max }}/{\text{%}}$
(0, 40]	8.86	86.72		(4 000, 19 000]	0.79	1.88
(40, 30 000]	1.85	4.34		(19 000, 30 000]	0.37	0.85
(40, 4 000]	9.99	23.34

下载: 导出CSV 

| 显示表格

由于低频趋势项不仅会导致信号基线漂移，还使得加速度信号积分后的速度信号严重失真^[18-19]，故EEMD方法处理零漂信号时，首先是保证较好地去除其低频趋势项，其次是控制修正后信号与原始信号在高频段上偏差尽可能小。故针对信号中失信频段0～40 Hz的修正指数上限${\bar D _{\max }}$，将其值的90%±2%分别代入拟合函数中去，可得两个白噪声系数：k₁= 0.056，k₂ = 0.066。

结合图8与表2，对不同的噪声系数k取值范围作如下讨论：

(1) 0＜k＜k₁。在此区间内，失信频段上的修正指数随白噪声系数变化的幅度很大，拟合函数最大值小于相应的修正指数上限值的88%。故认为在此区间失信频段修正效果不够稳定，且对信号失信频段的削弱效果不足，不满足处理方法对于较好去除其低频趋势项的要求。

(2) k₁≤k≤k₂。在此区间内，失信频段上的修正指数随白噪声系数变化的幅度相比于前一区间减小很多。修正指数平均值达到80.63%，比失信频段的修正指数下限8.86%提高了810%，已经达到修正指数上限值的92.98%。虽然在置信频段及其三个分段上，修正指数会随着噪声系数取值区间的上升而增大，可能导致在这些频段上出现信号修正过度的问题，但在这四个频段上前后两个区间的修正指数平均值分别只相差0.69%、3.73%、0.30%和0.13%，修正指数上升的幅度不大。故认为这个区间内修正后信号较好地去除了低频信号，且在置信频段尤其是高频段的频域失真很小，是可以接受的。

(3) k₂＜k ≤0.20。这一区间与[k₁, k₂]区间相比，噪声系数的增大会使置信频段尤其是高频段的频域失真增大，但与此同时信号失信频段的削弱效果不能得到有效提升，白噪声系数增大的意义很小。

综上分析，可以确定信号A对应的最优的白噪声系数k取值区间为[0.256，0.066]。

3.1.3   修正指数验证

为了验证修正指数$\bar D$在EEMD修正零漂信号中确定最优白噪声系数范围方面的适用性，现选取另一花岗岩试样上测得的零漂信号B，采用2.3节所示的方法进行处理，得到不同噪声系数下不同频段上的修正指数分布，如图9所示。零漂信号B和白噪声系数k=0.10时的修正信号如图10所示。

图  9  不同白系数噪声下不同频段上的冲击响应谱修正指数

Figure  9.  Correction indices of shock response spectrum on different frequency bands under different white noise coefficients

下载: 全尺寸图片 幻灯片

图  10  爆破加速度信号B及其修正信号

Figure  10.  Blast acceleration signal B and its corrected signal

下载: 全尺寸图片 幻灯片

通过图9看出，随着噪声系数的提高，EEMD方法修正后信号与原始信号B在冲击响应谱上的差异呈上升趋势，修正指数与噪声系数依然呈现如式(5)的幂指数函数关系。但是当噪声系数k上升到0.10之后，EEMD方法对失信频段0～40 Hz上信号的频域削弱几乎不再增加，但是对于置信频段上的三个分段上信号的频域削弱依然在增加，这也就意味着在(0.10, 0.20)范围内，增大噪声系数不能进一步削弱失信频段信号，反而会导致修正后信号在中高频失真加大。如采用3.1.2节中的方法，同样可以确定一个大致的最优噪声系数取值范围，即为[0.098，0.114]。

3.2   时域分析

从图3、图6和图10中可以看出，仅从加速度零漂信号本身来说，EEMD处理零漂信号A和B后得到的信号均很好地去除了原始信号中的零漂趋势项，其波形在爆破冲击结束后都很快回归基线，保留了原始爆破信号的衰减特征。

事实上，零漂偏差在加速度信号上表现得不如加速度时域积分后得到的速度信号显著，这是因为不仅加速度信号基线漂移会导致积分后速度信号基线漂移严重，而且原始信号上下振幅的不对称也会加剧积分曲线漂移。因此，对于加速度信号A及其在不同白噪声系数下处理得到的修正信号进行积分，得到的速度信号如图11所示。

图  11  EEMD修正前后速度曲线

Figure  11.  Velocity curves before and after correcting by EEMD

下载: 全尺寸图片 幻灯片

图11表明EEMD方法对于速度信号的零漂失真有明显改善：随着白噪声系数的提高，修正后的速度曲线逐渐靠近零点。当k值大于0.60时，EEMD对于速度信号的零漂修正效果几乎没有变化，这也与图7(a)中失信频段的修正指数分布图相契合。同时，这也表明通过修正指数分析确定信号A的最优白噪声系数范围[0.056，0.066]是合理的：将白噪声系数范围上调已经不能明显改善速度曲线零漂趋势；而将白噪声系数范围下调，速度曲线明显随之下滑，零漂趋势加重。但是，图11也表明了EEMD修正后的速度曲线无法完全消除零漂趋势，优化的噪声系数也很难让速度曲线基线完全回归零点。

3.3   进一步讨论

以上的分析表明，EEMD对于积分后的速度信号的零漂趋势去除得不够完善，存在一定的局限性。仅是考虑让积分后速度信号基线完全恢复到零点上，滑动平均法是一种简单有效的方法。影响滑动平均法去除信号零漂趋势的因素有两个^[8]：平滑点数和滑动阶次。对加速度信号A进行平滑点数为20和滑动阶次为10的滑动平均法处理，处理后的信号积分得到速度曲线如图12所示。图11和图12对比显示，滑动平均法得到的速度曲线是优于EEMD的，更符合爆破速度曲线的振动特征。

图  12  滑动平均法修正后速度曲线

Figure  12.  Velocity curve modified by sliding average method

下载: 全尺寸图片 幻灯片

通过图13给出的原始信号以及两种方法处理后信号的冲击响应谱，发现滑动平均法相对于EEMD方法在各个频段上频域损伤均更大。而冲击响应谱的低频部分，往往反映了岩石自身大的位移或者应变，其主要是由于零漂误差还是岩石自身的动力响应造成，还难以界定。所以，尽管滑动平均法很好地去除加速度信号积分后的速度信号的零漂趋势，但这也会使得修正信号的频谱失真加重以及岩石自身的动力响应削弱。

图  13  两种方法处理后信号与原始信号的冲击响应谱

Figure  13.  Shock response spectra of original signals and processed signals by the two methods

下载: 全尺寸图片 幻灯片

4.   结　论

本文结合花岗岩的室外小型爆破试验，对EEMD修正爆破零漂信号中最优白噪声系数范围进行了探析，得出主要结论如下：（1）利用EEMD方法并结合高低频处理能够较好地去除爆破加速度信号的零漂趋势，是一种自适应且高效的处理方法；（2）基于冲击响应谱定义的修正指数$\bar D$，能够较好地反映信号修正前后在冲击响应谱不同频段上的平均偏差幅度；随着白噪声系数增大，不同频段上修正指数$\bar D$均不同程度上升，二者呈现相关性较高的幂指数函数关系；（3）通过修正指数$\bar D$并结合对现场测试传感器的频响范围分析，确定出加速度零漂信号对应的白噪声系数的最优取值范围，优化的噪声系数使得EEMD法在最大程度削弱信号零漂趋势的同时尽量保留信号的频谱信息；（4）优化EEMD法中的噪声系数只能改善加速度积分后速度信号的零漂趋势而无法将其彻底去除；相对应的滑动平均法在速度信号的零漂修正效果上表现更好，但同时也带来了修正前后频谱差异加大以及动力响应的过度削弱的问题。