艾斯煜，杜博軍，王 軍,3

(1.蘇州科技大學(xué), 江蘇 蘇州 215009; 2.中國(guó)人民解放軍63850部隊(duì)， 吉林 白城 137001；3.中國(guó)科學(xué)院 長(zhǎng)春光學(xué)精密機(jī)械與物理研究所, 長(zhǎng)春 130033)

一般在靶場(chǎng)測(cè)試系統(tǒng)中，爆炸沖擊波的檢測(cè)對(duì)象主要是空氣超壓。而氣壓檢測(cè)系統(tǒng)易受爆炸所產(chǎn)生的光輻射和熱輻射影響。就本文中僅涉及的地面爆炸沖擊波檢測(cè)系統(tǒng)[1-2]而言，可以通過(guò)地下淺埋加速度傳感器檢測(cè)沖擊波的加速度值實(shí)現(xiàn)。

在爆炸沖擊波檢測(cè)領(lǐng)域，研究重點(diǎn)多為傳感器數(shù)據(jù)采樣，尤其是結(jié)合了前端數(shù)據(jù)處理的高精度采樣問(wèn)題。沖擊波中有效信號(hào)拾取是研究的關(guān)鍵，其方法主要有滑動(dòng)窗口法和相關(guān)耦合法?；瑒?dòng)窗口法是指在每個(gè)滑動(dòng)產(chǎn)生的窗口中計(jì)算窗內(nèi)時(shí)域波形的特征函數(shù)，根據(jù)函數(shù)的變化特征情況來(lái)標(biāo)示沖擊波的存在，并通過(guò)分析起突變情況以拾取沖擊波初至。而相關(guān)耦合法主要基于相鄰信道之間的波形相似性，對(duì)沖擊波波形信號(hào)上的部分或者整體波形截取并將其時(shí)間延遲，最終與待測(cè)沖擊波信道上的記錄對(duì)比以獲得沖擊波信號(hào)的存在性并確定波峰到達(dá)時(shí)間。

滑動(dòng)窗口法的最新應(yīng)用有學(xué)者時(shí)偉等[3]提出的雙約束變換時(shí)窗能量比法，該方法計(jì)算步驟簡(jiǎn)單，初至拾取效率高；左國(guó)平等[4]提出的滾動(dòng)能量比法，提高了計(jì)算精度和準(zhǔn)確度；鄭江龍等[5]提出的能量差法，將計(jì)算能量和比值轉(zhuǎn)換為能量差比值，提高了時(shí)窗檢測(cè)靈敏度。該類(lèi)方法雖然對(duì)多路信號(hào)獨(dú)立處理，無(wú)需考慮信號(hào)關(guān)聯(lián)性問(wèn)題，運(yùn)算簡(jiǎn)單、數(shù)據(jù)處理效率高，但對(duì)高頻低信噪比信號(hào)獲取效果較耦合法不明顯。

相關(guān)耦合法的最新研究成果有北京大學(xué)喻志超等[6]提出的利用波形相似性進(jìn)行初至?xí)r刻拾取的時(shí)窗選擇法,解決了沖擊波的微變化數(shù)據(jù)提取問(wèn)題；魏夢(mèng)祎等[7]學(xué)者提出了一種結(jié)合STA/LTA方法的基于波形互相關(guān)的沖擊波初至拾取方法，聯(lián)合多路數(shù)據(jù)為一個(gè)整體處理,實(shí)現(xiàn)初至獲取結(jié)果全局優(yōu)化；Akram等[8]國(guó)外學(xué)者基于迭代的互相關(guān)算法并在初至?xí)r刻獲取方面應(yīng)用，有效識(shí)別了地震沖擊波的微變化。但該類(lèi)方法普遍存在對(duì)高頻非線性、關(guān)聯(lián)性較弱的信號(hào)拾取效果不明顯等問(wèn)題。

本文提出了一種結(jié)合了變窗長(zhǎng)峰值濾波和耦合量互信息特性的改進(jìn)時(shí)窗法。以建模仿真、實(shí)驗(yàn)測(cè)量相結(jié)合的方式，完成了對(duì)爆炸沖擊波等實(shí)驗(yàn)數(shù)據(jù)的采集和效果偏差等后端數(shù)據(jù)處理。通過(guò)與傳統(tǒng)時(shí)窗法[9]對(duì)比，驗(yàn)證了該爆炸沖擊波檢測(cè)方法的可行性。

1 方法原理

1.1 時(shí)窗能量對(duì)比法拾取沖擊波信號(hào)

在三軸加速度傳感器組成的沖擊波檢測(cè)系統(tǒng)中[10-11]，沖擊波初至?xí)r刻是一個(gè)非常特殊的時(shí)刻，即此刻之前所記錄只有噪聲，而后所記錄則是噪聲與沖擊波的迭加信號(hào)。由于初至前后兩時(shí)窗積分能量有較大差異，故利用該特性來(lái)判斷初至?xí)r刻的發(fā)生點(diǎn)。時(shí)窗能量比E(i)反映了峰值時(shí)窗與前一時(shí)窗的變化程度,其表達(dá)式為：

(1)

其中，ck表示加速度傳感器測(cè)量值(m/s2)；M表示下一個(gè)時(shí)窗大小(s)；N表示當(dāng)前時(shí)間窗大小(s)；ε表示測(cè)試環(huán)境中的背景噪聲系數(shù)；i表示當(dāng)前時(shí)刻相鄰最大整數(shù)值。

傳統(tǒng)時(shí)窗能量對(duì)比法受固定窗長(zhǎng)制約，信號(hào)的高、低信噪比部分均采用同窗長(zhǎng)采樣，致使結(jié)果中高信噪比部分混疊較多噪聲信號(hào)，故而感知能力不佳，拾取精度不高。

1.2 互信息量概念

將隨機(jī)變量考慮成通信系統(tǒng)，兩信道即為2個(gè)相關(guān)變量?；バ畔⒘縄(x，y)表示兩個(gè)相關(guān)變量之間的統(tǒng)計(jì)學(xué)依存度，取值范圍為0～1。其中xi表示該系統(tǒng)的輸入量，yi表示系統(tǒng)的輸出量，通過(guò)互信息量這一概念建立兩路采集信號(hào)相關(guān)性。

互信息量表達(dá)式為：

(2)

式中：p(xi，yj)為2個(gè)變量的聯(lián)合分布；p(xi)、p(yj)為邊緣概率分布；o和s分別是X和Y的發(fā)生數(shù)量。

1.3 改進(jìn)后的時(shí)窗法拾取沖擊波信號(hào)

由于傳統(tǒng)時(shí)窗法中窗長(zhǎng)固定且對(duì)信號(hào)的高頻低信噪比部分感知效果差，進(jìn)而導(dǎo)致采樣精度低。本文結(jié)合變窗長(zhǎng)峰值濾波和互信息量概念對(duì)傳統(tǒng)時(shí)窗法進(jìn)行了改進(jìn)。

變窗長(zhǎng)峰值濾波是一種信號(hào)增強(qiáng)方法，操作步驟如下：

步驟1：加速度傳感器所直接采集輸入時(shí)域信號(hào)x(t)通過(guò)頻率調(diào)制得到解析信號(hào)z(t)：

(3)

在實(shí)際研究中原信號(hào)x(t)可表示為：

x(t)=s(t)+v1(t)

(4)

其中：s(t)是A通道加速度有效值；v1(t)是該通道的任意噪聲。

步驟2：對(duì)z(t)進(jìn)行Wigner-Ville轉(zhuǎn)換，獲得的時(shí)頻域信號(hào)Wz(t，f)為：

(5)

其中：h(τ)為時(shí)窗調(diào)節(jié)函數(shù)，表示時(shí)窗長(zhǎng)度且會(huì)隨所測(cè)數(shù)據(jù)做出調(diào)整；τ為時(shí)間變量(s);f為自變量頻率(Hz)。

步驟3：求取s(t)期望值，求取過(guò)程為:

(6)

其中：μ表示縮放比，由實(shí)際測(cè)試信號(hào)的強(qiáng)度設(shè)置。式(6)表示在對(duì)應(yīng)頻率f取值能讓W(xué)z(t，f)最大時(shí)，所對(duì)應(yīng)時(shí)域函數(shù)按一定比例縮放。

步驟4：選取時(shí)窗?？紤]到震動(dòng)模型與Ricker小波信號(hào)模型(如圖1所示)相似，時(shí)窗長(zhǎng)WL的選取參照Ricker小波信號(hào)模型常用式(7)，α為時(shí)窗選擇比，默認(rèn)0.707。針對(duì)沖擊波不同頻段，fs為波峰對(duì)應(yīng)頻率(Hz)；fd為有效信號(hào)帶寬(Hz)； 對(duì)應(yīng)時(shí)段為T(mén)s、Td，單位都為s。

(7)

圖1 Ricker小波信號(hào)模型

令2個(gè)地震信道記錄A和B采樣信號(hào)為x(t)和y(t)，根據(jù)沖擊波一般觀測(cè)方法和傳播特性，信號(hào)滿(mǎn)足以下模型：

x(t)=s(t)+v1(t)

y(t)=Q[s(t-D)]+v2(t)

(8)

B通道中Q[s(t-D)]是非線性時(shí)不變函數(shù)，它主要描述了傳播過(guò)程中沖擊波的波形變化程度。D是一個(gè)延遲參數(shù)，v2(t)為B通道觀測(cè)噪聲。通過(guò)檢測(cè)y(t)中不同延遲點(diǎn)D處的沖擊波形，并綜合通道A觀測(cè)信號(hào)x(t)進(jìn)而對(duì)沖擊波估計(jì)分析。

互信息量反映隨機(jī)變量間的統(tǒng)計(jì)關(guān)聯(lián)度，且不依賴(lài)于線性關(guān)系，可以有效避免非線性函數(shù)帶來(lái)的不利影響。在測(cè)量模型中x(t)為沖擊波信號(hào)，x(L)是波峰附近選擇的一段時(shí)窗，由于傳播通道間存在時(shí)延，可以將通道B看成A的輸出關(guān)聯(lián)量y(L+d)，即輸出量為輸入時(shí)窗長(zhǎng)加一段函數(shù)延遲d。在時(shí)窗大小不斷改變和移動(dòng)的情況下，求取對(duì)應(yīng)互信息量I(x，y)，I(x，y)最大值時(shí)關(guān)聯(lián)性達(dá)到最強(qiáng)，即此窗口時(shí)刻為該通道的初至?xí)r刻。

獲取A、B通道數(shù)據(jù)步驟如下：

1) 通過(guò)A、B通道加速度傳感器測(cè)量得到?jīng)_擊波數(shù)據(jù)x(t)，y(t+d)，計(jì)算得到其互信息量I(x，y)。

圖2給出了改進(jìn)算法的流程。

1.4 爆炸沖擊波超壓等效換算

圖2 改進(jìn)算法流程框圖

1.4.1加速度數(shù)據(jù)壓強(qiáng)轉(zhuǎn)換公式

通過(guò)實(shí)際測(cè)量和查閱相關(guān)資料[15-16]，壓電式加速度傳感器測(cè)量值與施加于其上壓強(qiáng)可以通過(guò)工程公式進(jìn)行轉(zhuǎn)換，公式如下：

(9)

反向推導(dǎo)得：

(10)

式中：R為檢測(cè)裝置與炸點(diǎn)中心間距(m)；d為傳感器探頭部分厚度(m)；amax為加速度測(cè)量最大值(m/s2)； Δpmax為接觸表面壓強(qiáng)最大變化量(MPa)；A1為常數(shù)，由剪切模量、材料密度等性能決定；α1、α2為加速度傳感器經(jīng)驗(yàn)參數(shù)(m/s2)，需要廠商進(jìn)行提供。

1.4.2薩道夫斯基公式

在1.4.1中，通過(guò)實(shí)際測(cè)量值換算得到?jīng)_擊波超壓，進(jìn)而通過(guò)薩道夫斯基算法，可以對(duì)其與TNT當(dāng)量進(jìn)行換算，薩道夫斯基公式如下：

(11)

2 系統(tǒng)搭建與測(cè)試結(jié)果

2.1 測(cè)試設(shè)備埋設(shè)布局

圖3所示為測(cè)試時(shí)炸點(diǎn)區(qū)及傳感器埋設(shè)布局示意圖，中心爆炸區(qū)為紅色半徑10 m圓形圍成的區(qū)域，三角形所示為布設(shè)的加速度傳感器，同心圓間隔2 m半徑,夾角約120°。圖4是測(cè)試系統(tǒng)框圖。

圖3 測(cè)試設(shè)備埋設(shè)布局示意圖

圖4 測(cè)試系統(tǒng)框圖

2.2 測(cè)試數(shù)據(jù)與分析

2.2.1仿真測(cè)試結(jié)果對(duì)比

通過(guò)3組Ricker小波信號(hào)進(jìn)行爆炸沖擊波信號(hào)仿真，包括主頻35 Hz、25 Hz和20 Hz三組反射信號(hào)和仿真頻率衰減[19]，其中增加了0.3 W高斯白噪聲。生成仿真波形如圖5所示。

圖5 源信號(hào)與不同方法提取的有效信號(hào)波形

圖5(a)是添加了白噪聲的源信號(hào)。圖5(b)和圖5(c)是濾波提取信號(hào)，三主頻對(duì)應(yīng)3個(gè)波峰。

圖5(b)為改進(jìn)時(shí)窗法提取的信號(hào)，采樣時(shí)窗長(zhǎng)為0.05s-0.1s。黑色線為信號(hào)波形，縱軸為傳感器的路數(shù)，橫軸為對(duì)應(yīng)時(shí)刻，亮色為所選時(shí)窗。圖5(c)是傳統(tǒng)時(shí)窗法所處理信號(hào)并選取的時(shí)窗，該方法取定長(zhǎng)時(shí)窗0.13 s。

對(duì)比圖5(b)和(c)可知，第二、三峰值附近波形相似，僅第一峰值差異較大。圖5(b)信號(hào)在第一小波峰即高頻部分保持效果較好。而圖5(c)信號(hào)中進(jìn)行了濾波去噪，對(duì)原信號(hào)高頻有效部分有明顯削弱。由此可知，相較于傳統(tǒng)時(shí)窗法而言，改進(jìn)時(shí)窗法不但能夠?qū)π盘?hào)低信噪比部分有效數(shù)據(jù)進(jìn)行獲取，也對(duì)高頻部分保留和增強(qiáng)。

2.2.2沖擊波實(shí)測(cè)效果對(duì)比

對(duì)11路沖擊波信號(hào)采用兩種方法采樣，頻率為1 000 Hz，采樣時(shí)長(zhǎng)為6 s。采樣曲線如圖6、圖7所示，圖6、圖7分別為傳統(tǒng)時(shí)窗法和改進(jìn)時(shí)窗法的解析數(shù)據(jù)與時(shí)窗選擇結(jié)果。 根據(jù)互信息量特性選取[2 2.5 3]、[3.5 4.5 5]、[5 5.5 6]三個(gè)初至?xí)r段。在該時(shí)段內(nèi)，圖5為1 s固定采樣時(shí)窗所取波形，圖6為改進(jìn)時(shí)窗法分別選取0.7 s、1.5 s、0.8 s窗長(zhǎng)解析所得波形。

圖6 傳統(tǒng)時(shí)窗法采樣曲線

圖7 改進(jìn)時(shí)窗法采樣曲線

圖6中沖擊波信號(hào)包括隨機(jī)噪聲，集中于框1、2、3中。噪聲在整個(gè)圖中分布使有效信號(hào)被抹去。從圖7的數(shù)據(jù)中可得出結(jié)論，在1.3中利用改進(jìn)時(shí)窗法去除高頻噪聲后，高頻低信噪比有效信號(hào)得以保留，信號(hào)得到增強(qiáng)。

2.2.3實(shí)際測(cè)試數(shù)據(jù)

如表1所示是采樣、換算所得數(shù)據(jù)。

表1 測(cè)量數(shù)據(jù)

在有效半徑范圍內(nèi)分別投放3 kg和6 kgTNT裝藥量測(cè)試彈，經(jīng)過(guò)前端采集和上位機(jī)換算得出上表數(shù)據(jù)。通過(guò)以上數(shù)據(jù)表，可以看出在10 m、12 m、14 m、16 m位置使用傳統(tǒng)時(shí)窗法的效果偏差均不小于7%，小部分甚至大于25%。

相較于傳統(tǒng)時(shí)窗法，改進(jìn)時(shí)窗法獲取沖擊波信號(hào)效果偏差降到10%以?xún)?nèi)，采樣精度提升明顯。

3 結(jié)論

1) 改進(jìn)后時(shí)窗法通過(guò)在采集前端增加了峰值濾波，降低了信號(hào)失真水平。

2) 在時(shí)窗選擇上，改進(jìn)后時(shí)窗法結(jié)合Ricker小波信號(hào)實(shí)現(xiàn)變窗長(zhǎng)選擇時(shí)窗，實(shí)現(xiàn)有效信號(hào)保持和增強(qiáng)。

3) 系統(tǒng)將采樣加速度有效值轉(zhuǎn)換為T(mén)NT當(dāng)量，滿(mǎn)足監(jiān)測(cè)數(shù)據(jù)直觀性要求。