爆破振動信號反問題研究初探

李立峰

（武漢理工大學資源與環境工程學院，湖北武漢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所示。

3 討論與結論

（1）單孔爆破振動信號（即特征孔信號）的合成屬于爆破工程中的一類反問題，該問題的研究具有重要的理論意義和實際應用價值。

（2）電子雷管與普通雷管相比，在延期時間上具有較高的精度，這為本研究中假設延遲時序為已知提供了足夠的支持。

（3）在反卷積法和基于群延遲的信號合成方法中，目前所研究的反卷積法對于較大規模的爆破情況，無能為力；相比較下，基于群延遲的方法，不僅可以隨機地合成一系列單孔爆破振動信號，并且對于較大規模的爆破振動，仍能合成具有合理包絡線的單孔爆破振動信號。

（4）現有的爆破振動反問題研究方法，仍有很大的改進空間，目前僅僅處于初探階段。下一階段的研究，將會集中在從多孔爆破振動信號中分離出對應于各個炮孔的單孔爆破振動信號。