爆破振动信号反问题研究初探

2020-11-14李立峰

金属矿山 2020年10期

李立峰

（武汉理工大学资源与环境工程学院，湖北武汉430070）

炸药用于岩石爆破历史悠久，可以追溯到10世纪左右，当时中国发明了黑火药。最初，岩石是被雷管和炸药瞬间引爆产生的能量破碎。到20世纪40年代的时候，毫秒延期雷管的使用产生了更大的影响［1-2］。在早期阶段，延期起爆的爆破设计主要侧重于岩石破碎，而不太考虑其副作用，包括空气冲击波和爆破振动。这种做法在当时是可以接受的，因为爆破工程主要集中在偏远地区。然而，随着采矿业的发展，爆破活动有时会接近有人居住的地区，这种转变使爆破活动的副作用越加明显。因此，爆破振动控制已成为爆破设计的必要部分。

为了控制爆破作业的不利影响，人们进行了许多相关的研究，并制定了爆炸安全规程。最初，爆破安全规程集中在一个变量上：质点峰值振速（PPV）［3-6］。为了预测现有爆破设计产生的爆破振动是否能得到控制，预测爆破振动是一个行之有效的方法。其中一种传统方法是使用比例距离法［7-8］，其常用表达式为

式中，R为爆心距；W为最大单段药量；K、β为与场地条件有关的参数。

还有其他研究者［9-11］采用3次方根的形式取代二次方根

另一种预测爆破振动的方法是利用人工智能。这种方法需要建立一个人工神经网络，然后根据爆炸参数通过学习过程预测爆破振动水平［12-13］。这种方法在现有文献中一般只预测PPV和主振频率，而不是整个振动波形。此外，人工神经网络需要一定的学习成本，如果学习过程中不考虑岩体性质或地质条件，预测结果的有效性可能是不确定的。然而，人工智能仍然是未来地面振动研究的一种有前景的方法，只是目前仍未充分发挥它的作用。

与仅预测爆破振动水平（如PPV）相比，爆破振动波形预测是一种更全面的方法，可以描述质点在整个爆破过程中的振动过程。该种方法提供的信息包括但不限于PPV和主振频率。如前所述，首次使用毫秒延期爆破是在20世纪40年代。在此期间，Thoenen和 Windes（1942）［14］开始注意到，可以使用某些适当的延期时间来减少爆破振动。由于当时雷管的技术有限，这个想法仍然只停留在理论阶段。直到20世纪80年代，Anderson（1983，1985）等［15，16］，Hinzen（1987）等［17］，Crenwelge（1988）［18］才开始利用信号与系统的理论，基于卷积计算发展出特征孔方法（signature hole method）。Anderson（2008）［19］对该方法进行了全面的综述。

特征孔法自提出以来，在爆破振动预测中发挥了重要作用。然而，传统的特征孔法由于其假设条件的限制而排除了许多随机因素。事实上，由于整个爆炸过程的内在随机性（如雷管、炸药、钻孔装药、岩石破坏、地质条件、波的传播路径等的随机性），实际的爆破效果并不能完全符合预期，各次爆破之间也不可能是一致的。Blair（1993）［20］和 Silva（2012）［21］改进了特征孔法，将随机因素和蒙特卡罗方法引入到爆破振动的预测中，得到了更实际的结果。在他们的方法中，每个炮孔产生的爆破振动波形是变化的，这种情况下一般的卷积计算不再适用，只能将其归入线性叠加。不论是排除随机因素的卷积模型，还是考虑随机因素的线性叠加模型，都是爆破振动传播现象的正问题。而本研究感兴趣的是，从监测到的爆破振动信号中，反向得到单孔爆破振动波形，这是一类典型的反问题。

1 爆破振动反问题

反问题是一个相对于正问题的概念。如果2个问题的表述都需要对另一个问题的充分或部分了解，那么这两个问题就互为反问题［22］。正问题和反问题之间没有绝对的区别。如果这2个问题中有一个已经在前面或更详细地研究过了，它就被称为正问题，而另一个则成为反问题。对存在于任何现象、过程或物理系统中反问题，可以通过时间、空间或因果顺序获得合理的理解。如果某些问题的解决遵循一定的时间、空间或因果顺序，它们就被称为正问题。相反，反问题是试图通过现象得到本质，或通过结果找到原因。具体到爆破振动现象，正问题和反问题在图1中从系统、输入和输出3个方面进行得到了阐释。

在图1中，中间的方框代表地震波传播的系统。对于爆破振动，系统是爆炸事件与爆破振动测点之间的岩土介质。系统的信息可以通过岩体的性质、地面介质的地质条件或反映地质条件的经验格林函数来揭示［24］。经验格林函数是在地震研究中广泛使用的方法，通常被称为小地震事件的地震测量［25-28］。在爆破工程中，单孔爆破也可以看作是经验格林函数［29］，这意味着它在一定程度上包含了地质条件，因此可以看作是整个爆破的特征孔，并用来进行卷积运算以预测爆破振动。

通过接收激励或输入，系统将产生反应或输出。输入是由炮孔发出的地震能量，输出则是在某一位置测量到的爆破振动。如果有一系列炮孔作为系统的输入（输入1到输入n），传感器测量所得的信号则是各炮孔相应输出的叠加。对于多个炮孔的爆破，通常是按照一定的时间顺序起爆，因此系统输入的信息也包括每个炮眼的起爆时间。

根据图1，爆破振动的监测和预测属于正问题。当爆破设计和一些系统信息已知时，可以预测距爆破中心一定距离处的爆破振动。但是，如果感兴趣的是在爆炸或波传播过程中系统或输入的信息，它就成为一个反问题。正问题和反问题通常是互补的。例如，利用预先测量的单孔波形作特征孔预测生产爆破的地面振动波形是一个正问题。相反，单孔波形的获取就是在解决一个相反的问题。与正问题相比，反问题通常更复杂，因为它们通常是不适定的问题。

Hadamard（1902）［30］定义的适定问题有如下性质:①问题有一个解；②解是唯一的；③解是数据的一个连续函数。不满足上述任何一个条件的问题，就成为不适定问题。

通常，一个正问题也是一个适定问题［31］。但反问题通常不满足以上3个条件中的1个或多个，这种情况增加了反问题求解的难度。因此，需要一些先验信息或对问题的附加限制。

2 单孔爆破振动信号的获取

2.1 爆破振动测试数据

本研究使用的数据来自于西弗吉尼亚州Guyan露天煤矿进行的试验，包含1个6孔爆破试验和1个83孔爆破试验，其爆破振动测试结果如图2所示。

对于6孔爆破测试，测点距爆炸中心210 m，其中包含了1个孔间延期为5 ms的6孔爆破和2孔爆破，以及3个单孔爆破测试；对于83孔爆破测试，测点距爆炸中心54.56 m，其中包含1个孔间延期为20 ms的83孔爆破和1个单孔爆破测试。

2.2 频域商反卷积法［23］

2.2.1 计算方法

如前所述，特征孔法预测爆破振动卷积运算可以表述为

式中，y[n]为爆破振动监测结果；d[n]为炮孔起爆时间序列；g[n]为特征孔信号；δ[n-ti]为带有时延ti单位脉冲响应；ai为幅度系数。

通过傅里叶变换将式（2）从时域变换到频域的式（3），就可以进行频域反卷积运算。

如果d[n]是已知的，也就是说如果炮孔起爆时间序列是已知的，那么特征孔g[n]的傅里叶变换G(f)就可以通过下述公式求出：

实际的点火时间ti，可以使用一个特殊的测量装置测量［32-33］。但如果在爆破中使用电子雷管，实际点火时间与标称延期时间相差很小，因此在这种情况下所设计延期序列可以近似地用作实际点火时间。为简单起见，可以假设ai为1。将G()f作傅里叶逆变换，就可以估计出单孔爆破振动波形，即特征孔波形g[n]。

2.2.2 零点不稳定性

在频域商方法的计算中，有2个重要问题:①分母为零的不稳定性；②在较高频率范围内的干扰成分。

以6孔爆破为例，其起爆延期序列由6个单位冲击响应组成，成为一个梳齿函数，如图3（a）所示。经过傅里叶变换，就生成如图3（b）所示的幅度谱。频域商的方法要求实测爆破振动信号的傅里叶变换和对应的梳齿函数的傅里叶变换具有相同的长度。因此，在进行傅里叶变换之前，需要补零至与实测爆破振动信号相同的长度。这样结果就是频谱中会出现零点。在图3（b）中，可以看到有1个零点出现在D(f)的512 Hz处。因为D(f)是式（1）中的分母，因此就有可能使解不稳定，产生无穷大的结果。

但是，从结果来看，零点的数目是有限的。因此，可以使用其他数值来代替这些零点。这里采用的方法是使用该零点前后2个数据点中较小的数据点来取代此零点。需要注意的是，零点的替换是复数的替换，而不仅仅是幅值的替换，见图4。

2.2.3 不同延期时间和孔数对频域商结果的影响

在解决了零点问题之后，就可以对y[n]和d[n]进行傅里叶变换，并利用公式（4）来进行频域商反卷积计算。从图5中可以看出，由于分母中在某些频率点存在极小值，因此会在商的结果中出现一些极大值，污染有用的成分。通常这一极大值可以使用一个低通滤波器去除，但是如果这一极大值在低频范围内，低通滤波器也会把有用成分过滤掉，使反卷积结果失去意义。因此，有必要考查不同的延期时间和孔数对商中的极大值所出现的频率范围的影响。而这又取决于分母D(f)中极小值点所出现的频率范围。通过对不同延期时间和孔数所形成的梳齿函数进行傅里叶变换，就可以看到这些极小值点的位置（如图6所示）。

从图6中可以看出，随着延期时间的增加以及孔数的增加，梳齿函数频谱中的极大值点有往低频区域移动的趋势。由此可以看出，频域商反卷积的方法对6孔爆破比较合适，但对文中83孔爆破的情况，就无能为力了。因此，文中仅列出6孔爆破情况下的结果，在不失一般性的情况下，文中只列出径向的结果（如图7所示）。

2.3 维纳滤波反卷积法

2.3.1 维纳滤波反卷积的计算

维纳滤波器本质上是一个最小二乘滤波器，其基本思想如下：

在图8中，期望输出是一个脉冲函数，而维纳滤波器的目的就是要使一个脉冲序列经过计算，尽量地和期望输出接近。将期望输出和实际输出以最小二乘法进行比较，形成一个函数I[n]用以计算最小误差。维纳滤波器的理论与分析已经有许多文献进行过详细的阐述［34-36］。本文只将维纳滤波器反卷积所涉及的必要公式列出。假设维纳滤波器以如下函数表示：

这一函数用来将输入的脉冲序列d[n]转变为实际输出：

而期望输出通过下式表示：

根据最小二乘理论，b[n]和c[n]之间的误差平方和I应该取最小值。将式（6）代入式（8），并使I对于f[n]的偏导数等于零，得到式（9）。

将式（9）做进一步的变换，得到

脉冲序列d[n]的自相关函数为

期望输出b[n]和实际输出d[n]的互相关函数为

结合式（10）～式（12），可得

式（13）若表示成矩阵形式，则为

上式就是著名的Wiener-Hopf公式，通过求解这一公式，就可计算出维纳滤波器f[n]的各个系数。对于维纳滤波器的效果如何，则可用滤波器性能参数P来表示［34-36］。

在式（15）中，P=0表示实际输出c[n]和期望输出b[n]没有相关性；而P=1意味着维纳滤波器可以使实际输出c[n]和期望输出b[n]完全一致。理论上讲，如果滤波器可以无限长，性能参数P可以完全等于1。但是一个数字滤波器一定是有限长的，所以只能适当地增加滤波器的长度或者调整期望输出的脉冲位置来提高P值。

2.3.2 维纳滤波反卷积在爆破振动信号中的应用［43］

爆破振动的卷积模型仍如式（2）所示，维纳滤波器就是要去除脉冲序列d[n]的影响，而估计出单孔爆破振动信号g[n]，如式（16）和式（17）所示。

以6孔爆破为例，其脉冲序列仍然如图3（a）所示，可以写成一个27×1向量：

首先令滤波器的长度等于d[n]的长度，即式（5）中K=L=26。因此，f[n]也是一个27×1的向量。根据式（6），滤波器的实际输出c[n]是一个53×1的向量。而式（7）中的期望输出b[n]则写为

根据式（11）和式（12）计算出d[n]的自相关函数rdd[t-s]以及d[n]和b[n]的互相关函数rbd[t]，并将其代入式（14），解之可得维纳滤波器f[n]。式（14）可以使用莱文森递归法求解［34，37］，其结果如图9所示。

将计算所得的维纳滤波器应用于实测的爆破振动信号，就能估计出单孔爆破振动信号，即特征孔信号gˆ[n]。估计的gˆ[n]对比实测的单孔信号g[n]，在振幅上偏小，因此可以用一个标量修正估计结果的幅值。图10仍然只列出了径向的计算结果，并且与3个实测的单孔爆破振动信号对比。

理论上讲，维纳滤波器避免了频域商反卷积因为分母的零点和极小值点而使解无意义的缺点，因而可以用于任何孔数和延期的爆破振动信号反卷积计算。然而，将前述同样的方法应用于文中的83孔数据后，却发现维纳滤波器仍难以给出满意的答案。

从图11中可以看出，估计所得的特征孔信号在波形包络上与实测的单孔爆破振动信号在时域和频域上均有很大的不同，与实际情况不相符。因此，仍然需要寻找适应性更加广的特征孔信号估计方法。

2.4 基于群延迟的单孔爆破振动信号合成

2.4.1 爆破振动信号与群延迟的关系

不同于式（2）的卷积模型，从更一般的情况考虑，实测的爆破振动信号其实就是每个单孔爆破振动信号的线性叠加，如式（19）所示。

式中，y表示实测的爆破振动信号；gi(i=0,1,…,D)表示D+1个炮孔分别产生的单孔爆破振动信号；ti(i=0,1,…,D)表示各个炮孔的起爆时间。

将式（19）作傅里叶变换，得到爆破振动信号的频域表达式：

式中，Y(ω)是实测爆破振动信号的傅里叶变换；Gi(ω)(i=0,1,…,D)是各个单孔爆破振动信号的傅里叶变换。

将式（21）中的欧拉公式代入式（20），改写成幅度和相位的形式：

式中，AY为整体爆破振动信号的幅度谱；为各个单孔爆破振动信号的幅度谱；θY为整体爆破振动信号的相位；ωti(i=0,1,…,D)为各个单孔爆破振动信号的相位，其中ωti表示相移。

式（22）告诉我们，爆破振动信号在频域由幅度谱和相位谱构成。若能够在幅度谱和相位谱上建立整体爆破振动信号和单孔爆破振动信号之间的关系，就能够根据实测的整体爆破振动信号估计出单孔爆破振动信号了。考虑到各个炮孔产生的爆破振动信号有所不同，在幅度谱和相位中引入随机变量，由此可以得到一系列各不相同的单孔爆破振动信号。

2.4.2 单孔爆破振动信号和多孔爆破振动信号幅度谱的关系

图12中给出了5个不同孔数的爆破情况下，多孔爆破振动和单孔爆破振动信号在幅度谱上的对比关系。可以看出，单孔爆破振动信号的幅度谱相比起多孔爆破振动信号的幅度谱不仅幅值比较小，曲线也更加平滑。因此，可以在经过平滑处理后的多孔爆破振动信号幅度谱曲线上乘以一个修正系数，用以估计单孔爆破振动信号的幅度谱，如式（23）所示。

考虑到每个炮孔产生爆破振动的随机性，要产生与炮孔数相等的一系列单孔爆破振动信号幅度谱，就需要再加入1个随机变量该随机变量服从分布

2.4.3 单孔爆破振动信号和多孔爆破振动信号相位的关系

这里不直接比较单孔与多孔爆破之间的相位，而是在相位对频率的导数上进行观察。相位对频率的导数称为群延迟［38］：

式中，τ(ω)为信号的群延迟；θ为信号的相位谱；ω为角频率。

群延迟的一个重要特点，就是其直方图（或者概率分布）形状与信号的时域波形的包络线是相似的，并且群延迟的均值点对应于信号包络的极值点［39-41］。如图13所示，若把爆破振动信号分为主体部分和尾部，群延迟的直方图形状和信号的主体部分包络线是相似的。根据这一特点，基于实测的爆破振动信号的群延迟的概率分布特征，来估计单孔爆破振动信号群延迟的概率分布参数，最终来估计单孔爆破振动信号的相位谱。

2.4.4 计算结果

当单孔爆破振动信号幅度谱和相位谱均估计出以后，就可以利用反傅里叶变换或者傅里叶级数的方法，合成一系列单孔爆破振动信号时域信号［42］。对于6孔爆破和83孔爆破2种情况，其计算结果如图14、图15所示。