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

李立峰

（武漢理工大學(xué)資源與環(huán)境工程學(xué)院，湖北武漢430070）

炸藥用于巖石爆破歷史悠久，可以追溯到10世紀(jì)左右，當(dāng)時(shí)中國發(fā)明了黑火藥。最初，巖石是被雷管和炸藥瞬間引爆產(chǎn)生的能量破碎。到20世紀(jì)40年代的時(shí)候，毫秒延期雷管的使用產(chǎn)生了更大的影響［1-2］。在早期階段，延期起爆的爆破設(shè)計(jì)主要側(cè)重于巖石破碎，而不太考慮其副作用，包括空氣沖擊波和爆破振動。這種做法在當(dāng)時(shí)是可以接受的，因?yàn)楸乒こ讨饕性谄h(yuǎn)地區(qū)。然而，隨著采礦業(yè)的發(fā)展，爆破活動有時(shí)會接近有人居住的地區(qū)，這種轉(zhuǎn)變使爆破活動的副作用越加明顯。因此，爆破振動控制已成為爆破設(shè)計(jì)的必要部分。

為了控制爆破作業(yè)的不利影響，人們進(jìn)行了許多相關(guān)的研究，并制定了爆炸安全規(guī)程。最初，爆破安全規(guī)程集中在一個(gè)變量上：質(zhì)點(diǎn)峰值振速（PPV）［3-6］。為了預(yù)測現(xiàn)有爆破設(shè)計(jì)產(chǎn)生的爆破振動是否能得到控制，預(yù)測爆破振動是一個(gè)行之有效的方法。其中一種傳統(tǒng)方法是使用比例距離法［7-8］，其常用表達(dá)式為

式中，R為爆心距；W為最大單段藥量；K、β為與場地條件有關(guān)的參數(shù)。

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

另一種預(yù)測爆破振動的方法是利用人工智能。這種方法需要建立一個(gè)人工神經(jīng)網(wǎng)絡(luò)，然后根據(jù)爆炸參數(shù)通過學(xué)習(xí)過程預(yù)測爆破振動水平［12-13］。這種方法在現(xiàn)有文獻(xiàn)中一般只預(yù)測PPV和主振頻率，而不是整個(gè)振動波形。此外，人工神經(jīng)網(wǎng)絡(luò)需要一定的學(xué)習(xí)成本，如果學(xué)習(xí)過程中不考慮巖體性質(zhì)或地質(zhì)條件，預(yù)測結(jié)果的有效性可能是不確定的。然而，人工智能仍然是未來地面振動研究的一種有前景的方法，只是目前仍未充分發(fā)揮它的作用。

與僅預(yù)測爆破振動水平（如PPV）相比，爆破振動波形預(yù)測是一種更全面的方法，可以描述質(zhì)點(diǎn)在整個(gè)爆破過程中的振動過程。該種方法提供的信息包括但不限于PPV和主振頻率。如前所述，首次使用毫秒延期爆破是在20世紀(jì)40年代。在此期間，Thoenen和 Windes（1942）［14］開始注意到，可以使用某些適當(dāng)?shù)难悠跁r(shí)間來減少爆破振動。由于當(dāng)時(shí)雷管的技術(shù)有限，這個(gè)想法仍然只停留在理論階段。直到20世紀(jì)80年代 ，Anderson（1983，1985）等［15，16］，Hinzen（1987）等［17］，Crenwelge（1988）［18］才開始利用信號與系統(tǒng)的理論，基于卷積計(jì)算發(fā)展出特征孔方法（signature hole method）。Anderson（2008）［19］對該方法進(jìn)行了全面的綜述。

特征孔法自提出以來，在爆破振動預(yù)測中發(fā)揮了重要作用。然而，傳統(tǒng)的特征孔法由于其假設(shè)條件的限制而排除了許多隨機(jī)因素。事實(shí)上，由于整個(gè)爆炸過程的內(nèi)在隨機(jī)性（如雷管、炸藥、鉆孔裝藥、巖石破壞、地質(zhì)條件、波的傳播路徑等的隨機(jī)性），實(shí)際的爆破效果并不能完全符合預(yù)期，各次爆破之間也不可能是一致的。Blair（1993）［20］和 Silva（2012）［21］改進(jìn)了特征孔法，將隨機(jī)因素和蒙特卡羅方法引入到爆破振動的預(yù)測中，得到了更實(shí)際的結(jié)果。在他們的方法中，每個(gè)炮孔產(chǎn)生的爆破振動波形是變化的，這種情況下一般的卷積計(jì)算不再適用，只能將其歸入線性疊加。不論是排除隨機(jī)因素的卷積模型，還是考慮隨機(jī)因素的線性疊加模型，都是爆破振動傳播現(xiàn)象的正問題。而本研究感興趣的是，從監(jiān)測到的爆破振動信號中，反向得到單孔爆破振動波形，這是一類典型的反問題。

1 爆破振動反問題

反問題是一個(gè)相對于正問題的概念。如果2個(gè)問題的表述都需要對另一個(gè)問題的充分或部分了解，那么這兩個(gè)問題就互為反問題［22］。正問題和反問題之間沒有絕對的區(qū)別。如果這2個(gè)問題中有一個(gè)已經(jīng)在前面或更詳細(xì)地研究過了，它就被稱為正問題，而另一個(gè)則成為反問題。對存在于任何現(xiàn)象、過程或物理系統(tǒng)中反問題，可以通過時(shí)間、空間或因果順序獲得合理的理解。如果某些問題的解決遵循一定的時(shí)間、空間或因果順序，它們就被稱為正問題。相反，反問題是試圖通過現(xiàn)象得到本質(zhì)，或通過結(jié)果找到原因。具體到爆破振動現(xiàn)象，正問題和反問題在圖1中從系統(tǒng)、輸入和輸出3個(gè)方面進(jìn)行得到了闡釋。

在圖1中，中間的方框代表地震波傳播的系統(tǒng)。對于爆破振動，系統(tǒng)是爆炸事件與爆破振動測點(diǎn)之間的巖土介質(zhì)。系統(tǒng)的信息可以通過巖體的性質(zhì)、地面介質(zhì)的地質(zhì)條件或反映地質(zhì)條件的經(jīng)驗(yàn)格林函數(shù)來揭示［24］。經(jīng)驗(yàn)格林函數(shù)是在地震研究中廣泛使用的方法，通常被稱為小地震事件的地震測量［25-28］。在爆破工程中，單孔爆破也可以看作是經(jīng)驗(yàn)格林函數(shù)［29］，這意味著它在一定程度上包含了地質(zhì)條件，因此可以看作是整個(gè)爆破的特征孔，并用來進(jìn)行卷積運(yùn)算以預(yù)測爆破振動。

通過接收激勵(lì)或輸入，系統(tǒng)將產(chǎn)生反應(yīng)或輸出。輸入是由炮孔發(fā)出的地震能量，輸出則是在某一位置測量到的爆破振動。如果有一系列炮孔作為系統(tǒng)的輸入（輸入1到輸入n），傳感器測量所得的信號則是各炮孔相應(yīng)輸出的疊加。對于多個(gè)炮孔的爆破，通常是按照一定的時(shí)間順序起爆，因此系統(tǒng)輸入的信息也包括每個(gè)炮眼的起爆時(shí)間。

根據(jù)圖1，爆破振動的監(jiān)測和預(yù)測屬于正問題。當(dāng)爆破設(shè)計(jì)和一些系統(tǒng)信息已知時(shí)，可以預(yù)測距爆破中心一定距離處的爆破振動。但是，如果感興趣的是在爆炸或波傳播過程中系統(tǒng)或輸入的信息，它就成為一個(gè)反問題。正問題和反問題通常是互補(bǔ)的。例如，利用預(yù)先測量的單孔波形作特征孔預(yù)測生產(chǎn)爆破的地面振動波形是一個(gè)正問題。相反，單孔波形的獲取就是在解決一個(gè)相反的問題。與正問題相比，反問題通常更復(fù)雜，因?yàn)樗鼈兺ǔＪ遣贿m定的問題。

Hadamard（1902）［30］定義的適定問題有如下性質(zhì):①問題有一個(gè)解；②解是唯一的；③解是數(shù)據(jù)的一個(gè)連續(xù)函數(shù)。不滿足上述任何一個(gè)條件的問題，就成為不適定問題。

通常，一個(gè)正問題也是一個(gè)適定問題［31］。但反問題通常不滿足以上3個(gè)條件中的1個(gè)或多個(gè)，這種情況增加了反問題求解的難度。因此，需要一些先驗(yàn)信息或?qū)栴}的附加限制。

2 單孔爆破振動信號的獲取

2.1 爆破振動測試數(shù)據(jù)

本研究使用的數(shù)據(jù)來自于西弗吉尼亞州Guyan露天煤礦進(jìn)行的試驗(yàn)，包含1個(gè)6孔爆破試驗(yàn)和1個(gè)83孔爆破試驗(yàn)，其爆破振動測試結(jié)果如圖2所示。

對于6孔爆破測試，測點(diǎn)距爆炸中心210 m，其中包含了1個(gè)孔間延期為5 ms的6孔爆破和2孔爆破，以及3個(gè)單孔爆破測試；對于83孔爆破測試，測點(diǎn)距爆炸中心54.56 m，其中包含1個(gè)孔間延期為20 ms的83孔爆破和1個(gè)單孔爆破測試。

2.2 頻域商反卷積法［23］

2.2.1 計(jì)算方法

如前所述，特征孔法預(yù)測爆破振動卷積運(yùn)算可以表述為

式中，y[n]為爆破振動監(jiān)測結(jié)果；d[n]為炮孔起爆時(shí)間序列；g[n]為特征孔信號；δ[n-ti]為帶有時(shí)延ti單位脈沖響應(yīng)；ai為幅度系數(shù)。

通過傅里葉變換將式（2）從時(shí)域變換到頻域的式（3），就可以進(jìn)行頻域反卷積運(yùn)算。

如果d[n]是已知的，也就是說如果炮孔起爆時(shí)間序列是已知的，那么特征孔g[n]的傅里葉變換G(f)就可以通過下述公式求出：

實(shí)際的點(diǎn)火時(shí)間ti，可以使用一個(gè)特殊的測量裝置測量［32-33］。但如果在爆破中使用電子雷管，實(shí)際點(diǎn)火時(shí)間與標(biāo)稱延期時(shí)間相差很小，因此在這種情況下所設(shè)計(jì)延期序列可以近似地用作實(shí)際點(diǎn)火時(shí)間。為簡單起見，可以假設(shè)ai為1。將G()f作傅里葉逆變換，就可以估計(jì)出單孔爆破振動波形，即特征孔波形g[n]。

2.2.2 零點(diǎn)不穩(wěn)定性

在頻域商方法的計(jì)算中，有2個(gè)重要問題:①分母為零的不穩(wěn)定性；②在較高頻率范圍內(nèi)的干擾成分。

以6孔爆破為例，其起爆延期序列由6個(gè)單位沖擊響應(yīng)組成，成為一個(gè)梳齒函數(shù)，如圖3（a）所示。經(jīng)過傅里葉變換，就生成如圖3（b）所示的幅度譜。頻域商的方法要求實(shí)測爆破振動信號的傅里葉變換和對應(yīng)的梳齒函數(shù)的傅里葉變換具有相同的長度。因此，在進(jìn)行傅里葉變換之前，需要補(bǔ)零至與實(shí)測爆破振動信號相同的長度。這樣結(jié)果就是頻譜中會出現(xiàn)零點(diǎn)。在圖3（b）中，可以看到有1個(gè)零點(diǎn)出現(xiàn)在D(f)的512 Hz處。因?yàn)镈(f)是式（1）中的分母，因此就有可能使解不穩(wěn)定，產(chǎn)生無窮大的結(jié)果。

但是，從結(jié)果來看，零點(diǎn)的數(shù)目是有限的。因此，可以使用其他數(shù)值來代替這些零點(diǎn)。這里采用的方法是使用該零點(diǎn)前后2個(gè)數(shù)據(jù)點(diǎn)中較小的數(shù)據(jù)點(diǎn)來取代此零點(diǎn)。需要注意的是，零點(diǎn)的替換是復(fù)數(shù)的替換，而不僅僅是幅值的替換，見圖4。

2.2.3 不同延期時(shí)間和孔數(shù)對頻域商結(jié)果的影響

在解決了零點(diǎn)問題之后，就可以對y[n]和d[n]進(jìn)行傅里葉變換，并利用公式（4）來進(jìn)行頻域商反卷積計(jì)算。從圖5中可以看出，由于分母中在某些頻率點(diǎn)存在極小值，因此會在商的結(jié)果中出現(xiàn)一些極大值，污染有用的成分。通常這一極大值可以使用一個(gè)低通濾波器去除，但是如果這一極大值在低頻范圍內(nèi)，低通濾波器也會把有用成分過濾掉，使反卷積結(jié)果失去意義。因此，有必要考查不同的延期時(shí)間和孔數(shù)對商中的極大值所出現(xiàn)的頻率范圍的影響。而這又取決于分母D(f)中極小值點(diǎn)所出現(xiàn)的頻率范圍。通過對不同延期時(shí)間和孔數(shù)所形成的梳齒函數(shù)進(jìn)行傅里葉變換，就可以看到這些極小值點(diǎn)的位置（如圖6所示）。

從圖6中可以看出，隨著延期時(shí)間的增加以及孔數(shù)的增加，梳齒函數(shù)頻譜中的極大值點(diǎn)有往低頻區(qū)域移動的趨勢。由此可以看出，頻域商反卷積的方法對6孔爆破比較合適，但對文中83孔爆破的情況，就無能為力了。因此，文中僅列出6孔爆破情況下的結(jié)果，在不失一般性的情況下，文中只列出徑向的結(jié)果（如圖7所示）。

2.3 維納濾波反卷積法

2.3.1 維納濾波反卷積的計(jì)算

維納濾波器本質(zhì)上是一個(gè)最小二乘濾波器，其基本思想如下：

在圖8中，期望輸出是一個(gè)脈沖函數(shù)，而維納濾波器的目的就是要使一個(gè)脈沖序列經(jīng)過計(jì)算，盡量地和期望輸出接近。將期望輸出和實(shí)際輸出以最小二乘法進(jìn)行比較，形成一個(gè)函數(shù)I[n]用以計(jì)算最小誤差。維納濾波器的理論與分析已經(jīng)有許多文獻(xiàn)進(jìn)行過詳細(xì)的闡述［34-36］。本文只將維納濾波器反卷積所涉及的必要公式列出。假設(shè)維納濾波器以如下函數(shù)表示：

這一函數(shù)用來將輸入的脈沖序列d[n]轉(zhuǎn)變?yōu)閷?shí)際輸出：

而期望輸出通過下式表示：

根據(jù)最小二乘理論，b[n]和c[n]之間的誤差平方和I應(yīng)該取最小值。將式（6）代入式（8），并使I對于f[n]的偏導(dǎo)數(shù)等于零，得到式（9）。

將式（9）做進(jìn)一步的變換，得到

脈沖序列d[n]的自相關(guān)函數(shù)為

期望輸出b[n]和實(shí)際輸出d[n]的互相關(guān)函數(shù)為

結(jié)合式（10）～式（12），可得

式（13）若表示成矩陣形式，則為

上式就是著名的Wiener-Hopf公式，通過求解這一公式，就可計(jì)算出維納濾波器f[n]的各個(gè)系數(shù)。對于維納濾波器的效果如何，則可用濾波器性能參數(shù)P來表示［34-36］。

在式（15）中，P=0表示實(shí)際輸出c[n]和期望輸出b[n]沒有相關(guān)性；而P=1意味著維納濾波器可以使實(shí)際輸出c[n]和期望輸出b[n]完全一致。理論上講，如果濾波器可以無限長，性能參數(shù)P可以完全等于1。但是一個(gè)數(shù)字濾波器一定是有限長的，所以只能適當(dāng)?shù)卦黾訛V波器的長度或者調(diào)整期望輸出的脈沖位置來提高P值。

2.3.2 維納濾波反卷積在爆破振動信號中的應(yīng)用［43］

爆破振動的卷積模型仍如式（2）所示，維納濾波器就是要去除脈沖序列d[n]的影響，而估計(jì)出單孔爆破振動信號g[n]，如式（16）和式（17）所示。

以6孔爆破為例，其脈沖序列仍然如圖3（a）所示，可以寫成一個(gè)27×1向量：

首先令濾波器的長度等于d[n]的長度，即式（5）中K=L=26。因此，f[n]也是一個(gè)27×1的向量。根據(jù)式（6），濾波器的實(shí)際輸出c[n]是一個(gè)53×1的向量。而式（7）中的期望輸出b[n]則寫為

根據(jù)式（11）和式（12）計(jì)算出d[n]的自相關(guān)函數(shù)rdd[t-s]以及d[n]和b[n]的互相關(guān)函數(shù)rbd[t]，并將其代入式（14），解之可得維納濾波器f[n]。式（14）可以使用萊文森遞歸法求解［34，37］，其結(jié)果如圖9所示。

將計(jì)算所得的維納濾波器應(yīng)用于實(shí)測的爆破振動信號，就能估計(jì)出單孔爆破振動信號，即特征孔信號g?[n]。估計(jì)的g?[n]對比實(shí)測的單孔信號g[n]，在振幅上偏小，因此可以用一個(gè)標(biāo)量修正估計(jì)結(jié)果的幅值。圖10仍然只列出了徑向的計(jì)算結(jié)果，并且與3個(gè)實(shí)測的單孔爆破振動信號對比。

理論上講，維納濾波器避免了頻域商反卷積因?yàn)榉帜傅牧泓c(diǎn)和極小值點(diǎn)而使解無意義的缺點(diǎn)，因而可以用于任何孔數(shù)和延期的爆破振動信號反卷積計(jì)算。然而，將前述同樣的方法應(yīng)用于文中的83孔數(shù)據(jù)后，卻發(fā)現(xiàn)維納濾波器仍難以給出滿意的答案。

從圖11中可以看出，估計(jì)所得的特征孔信號在波形包絡(luò)上與實(shí)測的單孔爆破振動信號在時(shí)域和頻域上均有很大的不同，與實(shí)際情況不相符。因此，仍然需要尋找適應(yīng)性更加廣的特征孔信號估計(jì)方法。

2.4 基于群延遲的單孔爆破振動信號合成

2.4.1 爆破振動信號與群延遲的關(guān)系

不同于式（2）的卷積模型，從更一般的情況考慮，實(shí)測的爆破振動信號其實(shí)就是每個(gè)單孔爆破振動信號的線性疊加，如式（19）所示。

式中，y表示實(shí)測的爆破振動信號；gi(i=0,1,…,D)表示D+1個(gè)炮孔分別產(chǎn)生的單孔爆破振動信號；ti(i=0,1,…,D)表示各個(gè)炮孔的起爆時(shí)間。

將式（19）作傅里葉變換，得到爆破振動信號的頻域表達(dá)式：

式中，Y(ω)是實(shí)測爆破振動信號的傅里葉變換；Gi(ω)(i=0,1,…,D)是各個(gè)單孔爆破振動信號的傅里葉變換。

將式（21）中的歐拉公式代入式（20），改寫成幅度和相位的形式：

式中，AY為整體爆破振動信號的幅度譜；為各個(gè)單孔爆破振動信號的幅度譜；θY為整體爆破振動信號的相位；ωti(i=0,1,…,D)為各個(gè)單孔爆破振動信號的相位，其中ωti表示相移。

式（22）告訴我們，爆破振動信號在頻域由幅度譜和相位譜構(gòu)成。若能夠在幅度譜和相位譜上建立整體爆破振動信號和單孔爆破振動信號之間的關(guān)系，就能夠根據(jù)實(shí)測的整體爆破振動信號估計(jì)出單孔爆破振動信號了?？紤]到各個(gè)炮孔產(chǎn)生的爆破振動信號有所不同，在幅度譜和相位中引入隨機(jī)變量，由此可以得到一系列各不相同的單孔爆破振動信號。

2.4.2 單孔爆破振動信號和多孔爆破振動信號幅度譜的關(guān)系

圖12中給出了5個(gè)不同孔數(shù)的爆破情況下，多孔爆破振動和單孔爆破振動信號在幅度譜上的對比關(guān)系?？梢钥闯?，單孔爆破振動信號的幅度譜相比起多孔爆破振動信號的幅度譜不僅幅值比較小，曲線也更加平滑。因此，可以在經(jīng)過平滑處理后的多孔爆破振動信號幅度譜曲線上乘以一個(gè)修正系數(shù)，用以估計(jì)單孔爆破振動信號的幅度譜，如式（23）所示。

考慮到每個(gè)炮孔產(chǎn)生爆破振動的隨機(jī)性，要產(chǎn)生與炮孔數(shù)相等的一系列單孔爆破振動信號幅度譜，就需要再加入1個(gè)隨機(jī)變量該隨機(jī)變量服從分布

2.4.3 單孔爆破振動信號和多孔爆破振動信號相位的關(guān)系

這里不直接比較單孔與多孔爆破之間的相位，而是在相位對頻率的導(dǎo)數(shù)上進(jìn)行觀察。相位對頻率的導(dǎo)數(shù)稱為群延遲［38］：

式中，τ(ω)為信號的群延遲；θ為信號的相位譜；ω為角頻率。

群延遲的一個(gè)重要特點(diǎn)，就是其直方圖（或者概率分布）形狀與信號的時(shí)域波形的包絡(luò)線是相似的，并且群延遲的均值點(diǎn)對應(yīng)于信號包絡(luò)的極值點(diǎn)［39-41］。如圖13所示，若把爆破振動信號分為主體部分和尾部，群延遲的直方圖形狀和信號的主體部分包絡(luò)線是相似的。根據(jù)這一特點(diǎn)，基于實(shí)測的爆破振動信號的群延遲的概率分布特征，來估計(jì)單孔爆破振動信號群延遲的概率分布參數(shù)，最終來估計(jì)單孔爆破振動信號的相位譜。

2.4.4 計(jì)算結(jié)果

當(dāng)單孔爆破振動信號幅度譜和相位譜均估計(jì)出以后，就可以利用反傅里葉變換或者傅里葉級數(shù)的方法，合成一系列單孔爆破振動信號時(shí)域信號［42］。對于6孔爆破和83孔爆破2種情況，其計(jì)算結(jié)果如圖14、圖15所示。

3 討論與結(jié)論

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

（2）電子雷管與普通雷管相比，在延期時(shí)間上具有較高的精度，這為本研究中假設(shè)延遲時(shí)序?yàn)橐阎峁┝俗銐虻闹С帧?/p>

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

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