付晓强， 俞 缙， 戴良玉， 秦双双

(1.三明学院 建筑工程学院，福建 三明 365004； 2.华侨大学 福建省隧道与城市地下空间工程技术研究中心，福建 厦门 361021； 3.三明科飞产气新材料股份有限公司，福建 三明 365500)

钻爆法是现阶段土石方开挖中应用最为广泛的施工方法之一[1-2]。爆破振动监测作为爆破损伤效应最行之有效的评价手段，其时频分析对于爆破方案优化和参数调整具有至关重要的作用。受测试环境、仪器性能和地质条件等因素的影响，爆破信号中普遍存在不相干分量，致使其时频谱中均不同程度存在交叉项干扰，从而对信号分析精度产生不利影响。而传统的以傅里叶变换为基础的信号分析存在表示结果不稀疏，不具备自适应性且时频分辨率差的缺陷，未能有效解决交叉项对信号辨识度的影响[3]。

近年来，在改变传统爆破信号表示方面取得了很大的进展，一定程度上丰富和推动了爆破信号时频分析理论的发展。如关晓磊等[4]采用极值点对称延拓，消除了希尔伯特-黄变换方法中的端点效应缺陷，改善了算法中瞬时频率出现负值问题，提高了爆破信号时频分析精度。郭涛等[5]将频率切片小波变换应用于爆破信号分析过程中，通过细化信号频率区间，提取到更为真实的信号时频特征，一定程度上抑制了交叉项干扰。杨仁树等[6]对隧道爆破信号进行集合经验模态分解，利用各模态信号分形特征确定优势分量并重构，采用平滑伪魏格纳分布(smooth pseudo Wigner-Ville distribution,SPWVD)对重构信号进行了时频分析，削弱了低频交叉项干扰的影响。同时应注意到，由于受爆破复杂环境和仪器参数设置等的影响，信号中含有的趋势项和噪声干扰成分会进一步污染其时频平面上的能量分布。趋势项和噪声通常为时频谱产生交叉项干扰的主要来源，趋势项的剔除对信号时频细节的表征具有明确的物理意义。如龙源等[7]分别采用最小二乘法、小波及经验模态分解(empirical mode decomposition,EMD)3种方法对爆破信号进行趋势项消除，并对3种方法的分析效果进行了对比，初步探讨了EMD算法在信号趋势项处理方面的优良特性，但EMD算法对信号端点处理不当却会引起模态混叠现象。上述诸多算法对爆破信号分析均取了良好的效果，但参数设置及处理过程较为复杂，面对海量的监测数据，算法缺乏自适应性严重制约了脑力时间的合理分配，甚至无法穷尽时间去精细刻画信号的时频特征。

针对上述难题，相关学者将匹配追踪(matching pursuit，MP)算法引入到爆破信号分析领域[8-9]。该算法实现了信号灵活表达，其基函数采用被称之为原子字典的过完备冗余函数(过完备原子库或原子字典)集合取代，可将信号表示为原子字典中少数向量的线性组合。这样，信号在过完备字典上的分解结果将是稀疏的，信号信息或能量集中在少数原子上，从而便于对信息的进一步处理。本文采用匹配追踪-平滑伪魏格纳分布(matching pursuit-smooth pseudo Wigner Ville distribution, MP-SPWVD)组合方法，对实测爆破振动信号时频谱中交叉项干扰消除难题进行了分析与实践，验证了该算法用于多因素导致的爆破信号时频交叉项消除的效果。

1 基本算法

1.1 MP算法

MP算法是在信号处理和函数逼近领域独立发展起来的，是信号稀疏表示最重要的算法之一[10]。它将信号按照字典原子逐步分解，通过多次迭代并选取原子库中匹配度高的原子近似逼近，最终获得最稀疏且具有明确物理意义的信号表示，其分解过程如图1所示：① 选择与信号Xn相似度高的原子Ψn并计算其投影值ɑn和差值信号Xn+1；② 对残差信号Xn+1重新选择原子进行匹配，获取其在最相近原子Ψn+1上的投影，得到差值信号Xn+2；③ 重复上述过程，当残差信号能量小于设定的阈值则计算终止。

图1 信号投影过程示意图

设D为有限维Hilbert空间H中的一个超完备字典，假定f为待分解信号(f∈H)，数据长度为N，则D满足[11]

(1)

式中:gγ为子波算子;Γ为伽马函数。

(2)

(3)

(4)

(5)

1.2 MP-SPWVD算法实现过程

Cohen类非线性时频分布中的魏格纳分布(Wigner-Ville distribution，WVD)算法会产生交叉项，伪魏格纳分布(pseudo Wigner-Ville distribution，PWVD)算法可一定程度上削弱交叉项，而相对WVD和PWVD分布，SPWVD时频解析精度更高，具体算法详见文献[14]。对于任意给定爆破振动信号X(t)，利用Hilbert变换可预先提取其瞬时频率和相位信息，从而很大程度上减少信号原子库的匹配数量，提高运算效率。具体步骤如下[15-17]。

步骤1选择Gabor原子并对其扩展，形成匹配过完备原子数据库Di(i=1,2,…,I)。

步骤2令初始差值信号r0等于待分解的初始信号X(t)，并将其代入过完备原子库中，筛选并预先确定子波分解原子库范围。

步骤3经过m次迭代后，从确定的原子库内找出与差值信号rm局部特征最为匹配的原子，也就是内积最大的原子dmi，将匹配到的原子从差值信号中减去，获取新的差值信号rm+1。

步骤4重复步骤2和步骤3，直至差值信号满足规定的条件。最终，X(t)被分解为

(6)

步骤5剔除信号中干扰成分R，分别对通过反复迭代获取的各分量信号cmidmi进行SPWVD运算并叠加，从而获得在时频平面的最优化特征分布。信号MP-SPWVD算法流程，如图2所示。

图2 MP-SPWVD信号处理流程

2 仿真信号分析

为了说明算法的优良特性，采用以子波瞬时频率变化表示信号信息的频率调制方式建立含4个高斯分量复杂调制信号，如图3(a)所示。各分量归一化频率中心分别为0.15,0.25,0.35和0.45，峰值均为2 cm/s，调制信号由各分量信号实部采用线性叠加形成，尽可能模拟多段别雷管起爆多频率多振型特点,如图3(b)所示。其中x1,x4波形起始于32 ms，x2,x3初振时刻为96 ms，为了提高运算效率，设定分量波形振动时长均为128 ms。调制信号WVD时频分布，如图3(c)所示。

(a) 调制信号

可以看出：信号时频面上除了4个真实高斯分量的能量分布外，亦包含任意两个分量时频中心连线处出现的交叉项，信号4个有效分量共产生6个交叉项，分别为交叉项1(由分量x1,x4产生)、交叉项2(由分量x3,x4产生)、交叉项3(由分量x2,x3产生)、交叉项4(由分量x1,x2产生)、交叉项5(由分量x1,x3产生)和交叉项6(由分量x2,x4产生)，且对角交叉项5、交叉项6互相叠加产生了重叠。信号PWVD时频分布将交叉项个数减少至两个，分别为x1,x4(交叉项1)和x2,x3(交叉项2)产生,如图4(a)所示。PWVD一定程度上减少了交叉项的个数，但仍未能从根本上抑制交叉项的产生。对调制信号进行MP迭代分解后，可准确获得4个高斯分量x1～x4的波形曲线，如图4(b)所示。对各分量分别求取SPWVD后叠加，便得到信号真实的时-频-能分布，如图4(c)所示。

(a) PWVD时频谱

调制信号WVD和PWVD在时频域本不应有能量存在的局部出现交叉项。而MP-SPWVD算法将信号分解为相互独立且包含局部结构特征的有限个原子的线性组合，由于原子间相互独立，则其SPWVD分布亦是相互独立且不存在交叉项，通过相互独立原子SPWVD分布叠加得到信号时频表达。因此，MP-SPWVD可有效避免交叉项的影响并改善多分量信号时频分辨率，进而获得更为真实的时频谱分布。

3 实例分析

3.1 隧道爆破信号

图5(a)为典型的隧道多段别爆破信号，该信号波峰值为3.36 cm/s,波谷值为2.77 cm/s,主频为180.22 Hz。隧道掘进断面41.2 m2，爆破采用MS1～MS15段别雷管跳段使用，其中，掏槽孔采用MS1、MS3段，单孔装药量0.8 kg；辅助孔采用MS5、MS7、MS9、MS11及MS13段，单孔装药量0.6 kg；周边孔采用MS15段，单孔装药量0.4 kg；底孔为MS15，单孔装药量0.6 kg。采用为采用“短进尺+弱爆破”方案，掏槽眼深度1.5 m，其余炮眼1.2 m，单循环总装药量为56.2 kg。

图5(b)功率谱表明：信号在低频段存在明显的趋势项，在起始频率0～3 Hz附近产生0.02 dB/Hz的奇异值，而在大于500 Hz的频带出现显著的低幅高斯白噪声，信号频谱中存在明显的伪信息，其在信号时频谱上产生无法解释的交叉项干扰，严重影响了对信号特征的解读和判别。图5(c)中MP重构信号功率谱有效克服了低频趋势项对信号真实频谱的影响，获得了具有明确物理意义的信号频谱特征。

(a) 隧道爆破信号

3.2 信号MP分解与重构

在信号频谱分析中，Gaussian函数具有良好的时频聚集性[18]。因此，文中建立的Gabor原子库中的每个原子都是由Gaussian函数经伸缩、平移和调制得到的，如图6(a)所示。

通过引入本征时间尺度分解得到信号一系列不同频率段的PR(proper rotation)分量[19-20]。对各PR分量进行Hilbert变换，获取信号的瞬时优势频率和相位，随后将其代入MP算法中，可大大减少程序循环步数。采用1.1节的MP算法对原始信号进行处理，建立与分析信号相似度高，长度为10 ms的过完备稀疏Gabor原子字典库。实践证明，迭代次数满足10n(n≤3且为正整数)指数关系时能取得较为理想的分析效果，n较小时(n=1)存在分解不彻底缺陷，而n较大时(n=3)易产生过分解，同时运算时长激增。根据爆破信号特点及子波分量频带限制，此处设置迭代次数为100(n=2)，从而将原信号分解为7个子波分量，如图6(b)所示。由于MP算法是以信号时频能量高低对其进行分解的，因此将迭代获得的前4个能量占比高的分量进行线性叠加，得到重构后的真实信号如图6(c)所示。

从图6中可知，重构信号很好地继承了原始信号的波动特征，同时对信号波形局部进行了平滑拟合。为了便于分析，分别标识原始信号、MP处理信号为X1,X2，采用文献[21]中的交叉小波变换方法，求取两者的相关性凝聚谱,如图7所示。

(a) Gobor原子库

图7 信号相关性谱

相关性凝聚谱X2和X1在时域和主振频域上均为显著正相关。X2与X1在低频部分(<32 Hz)和高频部分(>512 Hz)为负相关，有效滤除了信号中的低频趋势项和高频噪声且相关性在时间轴上连续统一，显著性检验基本贯穿小波影响椎范围内的整个时频域。定义互相关系数为

λs1,s2=max[Rs1(τ)]、max[Rs2(τ)]，

λs1,s2∈(0,1)

(7)

式中:Rs1(τ),Rs2(τ)为分别为X1,X2的自相关值;λs1,s2值越大，则两信号间互相关性越强。计算X2与X1之间的互相关系数为0.968，交叉小波变换清晰描述X2与X1在时域和频域不同层次的波动性和共振相位相关度。MP算法在信号处理过程中对低频趋势项进行了有效校正，同时对高频噪声分量起到了良好的抑制效果，体现了算法信号处理优势。

3.3 时频分析

为了分析爆破信号交叉项干扰的抑制效果，分别采用WVD，PWVD和MP-SPWVD算法求取信号的时频谱，如图8所示。隧道爆破信号WVD时频谱交叉项干扰严重，时频谱交叉项出现在原本不存在能量的位置，时频分布的可辨性差。同时，段间延期时间越短，交叉项越明显。PWVD谱削弱并剔除了信号时频面上的高频交叉项，但对于低频交叉项分布却有一定程度放大，表现在时频面上的分布密度增强且分布形态趋于一致。MP-SPWVD时频谱能量分布清晰可辨，实现了多段别微差爆破起爆能量的有效分割。爆破产生的能量主要聚集在掏槽段起爆0～35 ms时程内，能量峰值位于归一化频率0.05～0.07处(采样频率为3 000 Hz，约150～210 Hz内)，低频交叉项位于0～35 ms及归一化频率0～0.02处(约0～60 Hz)，而高频交叉项主要位于波形时程尾部边缘，且集中在归一化频率0.15～0.25处(约450～750 Hz)。MP-SPWVD算法时频分布能够随着被分析信号局部特征自适应地调节原子基函数的特征参数，如时宽和带宽，避免了传统方法中窗效应的影响，还可很好地匹配信号的局部特征，得到的时频谱对于低频漂零和高频边缘交叉项均有优良的抑制效果。

(a) WVD时频谱

4 可靠性验证

为了验证1.2节组合方法的有效性，选取煤矿冻结立井某次爆破井壁振动信号，如图9(a)所示。信号采用文献[22]中井壁传感器预埋法采集，测点位于爆破掌子面距离为6 m的井壁处。由于测点位于爆破近区，该信号含有明显的偏离基线的趋势项，导致其在时间轴上0.1 s后出现了“甩尾”(基线漂零)现象。爆破采用雷管为MS1～MS5 5个段别，其中，周边眼采用的MS5段雷管的起爆及误差区间为(110±15) ms，与趋势项在时间轴上的出现区间重叠。因此，趋势项引起的时频谱交叉项干扰对其频谱特征的解读产生误判，特别是0.1 s后的时频平面上的能量分布尤其值得关注。

对图9(a)中井壁振动信号分别采用小波方法和MP算法进行信号恢复，其结果见图9(b)和图9(c)。图9对比表明：小波方法信号重构过程中，阈值选取不当易使信号产生了微幅的高频振荡噪声，同时信号幅值大幅降低，损失了信号部分有用信息。MP算法校正并剔除了信号中的趋势项分量，信号特征辨识度高。

(a) 立井爆破信号

求取原始信号WVD,PWVD和MP-SPWVD时频谱分别如图10(a)～图(c)所示。

(a) WVD时频谱

从图10可知：WVD时频谱窗函数不可调，致使其时频分辨率是固定的，存在严重的交叉项干扰。PWVD时频谱虽然改善了交叉项干扰的影响程度，但对于时频谱在频率边界处的交叉项干扰有一定程度的放大。MP-SPWVD时频谱是通过相互独立的原子重构子信号SPWVD相互叠加构成的，可有效抑制交叉项的产生，其更加符合微差爆破信号的局部特征，能够更清晰揭示立井信号的时频分布特性。

同时应注意到，MP-SPWVD算法也有一定不足之处，即当待提取振动信号数据较长时，原子数量将急剧增加，这种情形下会严重降低运算速度。选取文献[23]中高斯调制信号、三正弦和典型爆破信号在时间轴上分别进行人为延拓至1 s,3 s,5 s，来比较不同信号长度处理所耗机时，具体如表1所示。

由表1可知，随着信号长度和复杂性的增大，分析所耗机时也会增加，尤其对于采样点密集或持续时长较大(大于5 s)的爆破信号计算量过大，严重时甚至导致内存溢出。因此求解过程需耗费很长机时，对求解器配置和性能有特殊要求。

表1 不同信号CPU运行机时

5 结 论

(1) 爆破信号时频谱提供了信号在时域和频域的联合分布，清晰描述了能量在时频域上的变化关系。受测试环境及仪器自身原因的影响，爆破信号时频谱中普遍存在交叉项干扰，交叉项的存在严重影响了对爆破信号信息特征的判别，致使时频域能量分布在局部物理意义不明确，影响后续爆破参数调整优化和效果评价。

(2) MP算法实现了爆破信号的稀疏化分解和降噪重构，MP-SPWVD分布可有效抑制了时频谱交叉项，降低不确定因素对信号谱特征的干扰，精细刻画爆破信号能量分布的非平稳时变特征。交叉项抑制效果与雷管段别及延期时间密切相关。

(3) MP-SPWVD算法虽是一种高精度的信号恢复和时频表征方法，但受信号时长的限制较为严重，尤其对于采样密集或持续时长大于5 s以上的信号数据，求解速度受到一定的制约，所耗机时激增，对求解器配置和性能有特殊要求。后续研究过程中，应重点优化长时信号计算速度，提高算法的高效性。