0 引言

因控制需求，微差爆破被广泛应用于露天爆破、拆除爆破和地下爆破等工程.精确的微差时间，是达到优良爆破效果的前提和保证^[1]，但微差雷管本身存在着一定误差和合格率问题，相邻段别雷管甚至存在“跳段”现象，致使起爆顺序混乱，爆破效果和爆破安全达不到设计要求^[2-4].

众多学者对微差爆破延时识别进行研究，取得了诸多成果，先后出现了八线示波器和金属拉线方法^[5]、小波模极大值法^[6]和小波时-能密度等方法^[7].但前者操作复杂，精度有限；小波类方法借助实测爆破振动信号进行识别，取得了不错的效果，但存在分解度和小波基选择问题^[8-9]，不同的分解尺度和小波基，识别结果亦有差异.为此，采用HHT方法对微差爆破进行延时识别研究.

1 HHT方法简介 1.1 EMD分解

HHT（Hilbert-Huang Transform，简称HHT）是N.E.Huang等^[10]提出的信号处理方法，由EMD（Empirical Mode Decomposition，简称EMD）分解和Hilbert变换组成.针对爆破振动类非平稳数据，EMD分解无需先验基，能随信号的变化自适应分解成有限若干IMF分量，其原理如下^[11-12]：

步骤1：原始信号为S（t）（t=1，2，…，n，为时间采样点数），设定IMF分量的判断条件：①极点数和零点数差值不超过1；②上、下包络线以时间轴为准，局部对称，也即均值为0.

步骤2：确定S（t）所有极值点，通过3次样条插值获取S（t）的上包络线S_max（t）和下包络线S_min（t），并求取上、下包络线的均值曲线m（t）：

$m(t) = \frac{1}{2}\left[ {{S_{\max }}(t) + {S_{\min }}(t)} \right]$

将m（t）从S（t）中输出，则剩余值h₁（t）为：

${h_1}(t) = S(t) - m(t)$

步骤3：如果h₁（t）符合IMF分量判断条件，将其作为第一个IMF分量输出，否则代替S（t）重复筛选过程，直至k（k=1，2，…，n）次迭代后，剩余值h_k（t）成为一个IMF，即IMF₁（t）=h_k（t）.

将IMF₁（t）输出后，第一阶段的剩余信号R₁（t）为：

${R_1}(t) = S(t) - IM{F_1}(t)$

步骤4：对R₁（t）继续重复步骤2、步骤3的筛选工作，从高频到低频依次输出IMF₂（t）、IMF₃（t）、…、IMF_i（t），直至最后的残差：

${R_i}(t) = {R_{i - 1}}(t) - IM{F_i}(t)$

当R_i（t）再也无法分解出新的IMF分量时，终止筛选工作.那么信号S（t）为IMF₁（t）、IMF₂（t）、…、IMF_i（t）及余量R_i（t）之和：

$S(t) = \sum\limits_{i = 1}^n {IM{F_i}(t)} + {R_i}(t)$

1.2 分量包络和Hilbert谱

EMD分解后获取若干个IMF分量，对各分量作Hilbert变换：

$H\left[ {IMF(t)} \right] = \frac{1}{{\rm{ \mathsf{ π} }}}PV\int_{ - \infty }^{ + \infty } {\frac{{IMF(t')}}{{t - t'}}} {\rm{d}}t'$

其中，PV为柯西主值（Cauchy Principal Value），因此，构造信号z（t）：

$z(t) = IMF(t) + jH\left[ {IMF(t)} \right] = a(t){e^{j\varphi (T)}}$

其中，a（t）为z（t）的幅值函数，也即，分量IMF（t）的包络；φ（t）为相位函数：

$\begin{array}{l} a(t) = \sqrt {IM{F^2}(t) + {H^2}\left[ {IMF(t)} \right]} \\ \varphi (t) = \arctan \frac{{H\left[ {IMF(t)} \right]}}{{IMF(t)}} \end{array}$

于是，瞬时频率可表示为：

$f(t) = \frac{{\omega (t)}}{{2{\rm{ \mathsf{ π} }}}} = \frac{1}{{2{\rm{ \mathsf{ π} }}}} \times \frac{{{\rm{d}}\varphi (t)}}{{{\rm{d}}(t)}}$

如果将趋势分量R（t）忽略，Re表示取实部，则S（t）可表示为：

$S(t) = \mathit{Re}\sum\limits_{i = 1}^n {{a_i}(t)} {e^{j{\varphi _i}(t)}} = \mathit{Re}\sum\limits_{i = 1}^n {{a_i}(t)} {e^{j\int {{\omega _i}(t)} {\rm{d}}(t)}}$

表达了在时间和频率平面上瞬时振幅的分布规律，将其定义为Hilbert谱：

$H(\omega ,t) = \mathit{Re}\sum\limits_{i = 1}^n {{a_i}(t)} {e^{j\int {{\omega _i}(t)} {\rm{d}}(t)}}$

1.3 瞬时能量

依据巴什瓦（Parseval）定理^[13]，可将|IMF_i（t）|²视作是各IMF分量的能量密度，如此，H²（ω，t）同样具备能量密度的物理意义，将其称作Hilbert能量谱，则在Hilbert-Huang变换前后，信号的能量应该守恒，也即：

${\int_{ - \infty }^{ + \infty } {\left| {S(t)} \right|} ^2}{\rm{d}}t = \int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } H } (\omega ,t){\rm{d}}\omega {\rm{d}}t$

由此定义：

$IE(t) = \int_\omega {{H^2}(\omega ,t)} {\rm{d}}\omega $

表达了信号的频带能量随时间的分布情况.

2 MATLAB信号仿真

某矿采用水平深孔阶段矿房法采矿，对其某2次采场崩矿进行爆破振动监测，见图 1，信号S₁和S₂分别为第1次和第2次监测所得，详细爆破信息见表 1.

以S₁为例，对其进行EMD分解，获取了9个按照频率高低依次排列的IMF分量和一个余量R，见图 2（篇幅有限，仅列前5个IMF分量）.其中，IMF₁的频率最大，0.24~0.4 s存在衰减，但整个时间坐标轴分散着少量噪声，是相对高频部分；IMF₂的波形衰减明显，幅值最大，是与真实信号相关性最大的分量；IMF₃~IMF₅分量的幅值已大幅减小，携带真实信号的信息已逐渐减小.

3 微差延时识别

在微差爆破中，每段雷管的起爆就意味着炸药能量的释放，引发测点爆破振动信号速度或者能量突增^[14-16].由此可知，获取爆破主振分量的幅值包络线或瞬时能量曲线，寻找包络线或能量突峰点并计算对应的起爆时刻，就能获取各段雷管的实际延期时间，分别定义为微差爆破延时的幅值包络线识别法和瞬时能量识别法.

经EMD分解和识别，信号S₁和S₂均是IMF₂幅值大，波形衰减明显，为主振分量，携带了爆破的大部分信息，图 3和图 4分别是信号S₁和S₂的IMF₂分量幅值包络线.图 5和图 6则分别是信号S₁和S₂的瞬时能量.

图 3~图 6中突峰点计算可知，由幅值包络线法和瞬时能量法识别的信号S₁的两段雷管起爆时刻相同，分别为0.267 6 s和0.322 3 s，延时时间为0.054 7 s；由幅值包络线法识别的信号S₂的两段雷管起爆时刻分别为0.267 1 s和0.363 8 s，延时时间为0.096 7 s，由瞬时能量法识别的信号S₂的两段雷管的起爆时刻分别为0.263 2 s和0.363 8 s，延时时间为0.100 6 s.实际爆破设计使用的是国产第二系列毫秒延期雷管，对比微差雷管的理论延时和实测间隔时间（见表 2），可知幅值包络线法和瞬时能量法识别的延时时间都在理论值内，该批次雷管精度可靠，矿山可放心使用，也说明这2种方法用于识别微差爆破延时可行.

4 结论

借助MATLAB工具，基于某矿崩矿爆破实测数据，对微差爆破延时时间进行识别研究，得出了以下结论：

1）EMD分解具有自适应性，经Hilbert变换后，可以有效方便地获取幅值包络线和瞬时能量，为微差爆破延时时间识别提供基础.

2）幅值包络线法和瞬时能量法识别的微差爆破延时时间都在理论值内，证明该批次雷管安全可靠，也说明这2种方法识别微差爆破延时可行，为微差爆破延时识别提供了方法.

3）虽然幅值包络线法和瞬时能量法均可用于微差爆破延时识别，但识别结果有所差异（如信号S₂）.当信号的信噪比较低时，瞬时能量法识别易受噪声能量干扰，出现错误突峰点或突峰点被噪声淹没无法识别现象，因此信号应先进行去噪处理；EMD分解也会出现多个接近的分量，爆破信息被分散，没有绝对主振分量，此时重构多个相近分量形成爆破振动优势分量，然后再提取幅值包络线进行延时识别.