劉曉佳，李 劍，劉代勁， 魏曉曼， 孔慶珊，金 艷

(1.中北大學(xué) 信息探測與處理山西省重點(diǎn)實(shí)驗(yàn)室，太原 030051；2.山東省軍區(qū)數(shù)據(jù)信息室，濟(jì)南 250099； 3.湖南云箭集團(tuán)有限公司， 長沙 410000)

0 引言

近年來，隨著國際形勢以及高新技術(shù)的發(fā)展，世界各國加強(qiáng)了在軍事防御、武器裝備等國防安全領(lǐng)域的研發(fā)力度，各類新型武器裝備層出不窮[1]。與此同時，爆炸波動理論、沖擊波測試技術(shù)、場重建技術(shù)等研究也發(fā)揮著越來越重要的作用[2]，各類武器彈藥主要是憑借沖擊波和破片來達(dá)到毀傷作用，在其研制、驗(yàn)收等過程中，如何準(zhǔn)確的掌握沖擊波超壓在爆炸區(qū)域的分布特性一直是研究的重點(diǎn)[3]。爆炸是在極短時間內(nèi)，釋放出大量能量，產(chǎn)生高溫，同時破壞性極強(qiáng)，難以通過直接觀察的方式研究。隨著戰(zhàn)斗部研制向著智能化，一體化發(fā)展，對于爆炸場的壓力毀傷能力等參量研究具有重要意義。沖擊波場重建的方法有走時層析法、插值法、數(shù)值模擬法等，其中，基于走時層析成像的方法是沖擊波場重建的重要方式[14]，沖擊波的到達(dá)時間是實(shí)現(xiàn)走時層析成像的關(guān)鍵參數(shù)，隨著戰(zhàn)斗部威力的不斷加大，對其毀傷威力的評估區(qū)域發(fā)生了變化，沖擊波場重建作為毀傷評估的重要方式，沖擊波場重建的區(qū)域隨著戰(zhàn)斗部威力的增大而擴(kuò)大，場重建區(qū)域由原來的小區(qū)域向大區(qū)域發(fā)展，在近場時，沖擊波信號是一個非平穩(wěn)的瞬時脈沖信號，現(xiàn)有的沖擊波拾取到達(dá)時間的方法已經(jīng)可以滿足近場沖擊波信號到達(dá)時間的拾取，但是對于中遠(yuǎn)場區(qū)域的沖擊波信號而言，沖擊波信號由非平穩(wěn)的瞬時信號變?yōu)樯仙鼐徛男盘?，并且隨著爆炸過程本身具有底噪大，強(qiáng)干擾、高溫等綜合影響，現(xiàn)有的提取沖擊波到達(dá)時間的方法不適用，針對這個問題，本文針對中遠(yuǎn)場區(qū)域的復(fù)雜條件下的微弱沖擊波信號的到達(dá)時間提取問題，提出了基于時窗熵的沖擊波到達(dá)時間提取方法。

沖擊波到達(dá)時間提取是沖擊波信號處理和成像的關(guān)鍵步驟，人工拾取到時信息，既耗時又不準(zhǔn)確，為了提高到達(dá)時間拾取的效率和速度，通常根據(jù)某些特定指標(biāo)使拾取過程實(shí)現(xiàn)自動化提取，隨著越來越多的研究，自動檢測到達(dá)時間的方法已經(jīng)成為陣列化數(shù)據(jù)處理和沖擊波數(shù)據(jù)處理的重要組成部分[22]。

多年來，研究學(xué)者提出了很多關(guān)于到達(dá)時間提取的方法，最常見的是使用半自動/自動提取算法，再通過視覺檢查交互式地校正結(jié)果，而且在這一個過程中的某些信號組成復(fù)雜的位置需要重復(fù)幾次，因此，當(dāng)數(shù)據(jù)量大且信噪比低時，初至波拾取的過程耗時較長。傳統(tǒng)的到達(dá)時間拾取方法有分形錐數(shù)法、能量比算法、AIC法、相關(guān)法、偏振分析法和基于熵等方法。其中，能量比算法因?yàn)橛?jì)算方法簡單高信噪比情況下精度較高而被廣泛應(yīng)用，標(biāo)準(zhǔn)短期平均值/標(biāo)準(zhǔn)長期平均值(STA/LTA,sort term average/long term average)又稱為長短時窗法是能量比算法的代表，最常用的特征函數(shù)是能量比，通過計(jì)算兩個窗口內(nèi)的能量比，根據(jù)能量比最大的時間位置判定為是沖擊波到達(dá)的時間，Coppens F[6]指出，當(dāng)信噪比高的時候，基于能量比的方法所得到的沖擊波到達(dá)時間是比較準(zhǔn)確的；但是長短時窗法的抗噪性差，沖擊波的到達(dá)時間提取的過程，大多為低信噪比信號的信號到達(dá)時間提取，諸多學(xué)者對于長短時窗法進(jìn)行了改進(jìn)，使其在低信噪比條件下也可以精確拾取沖擊波的到達(dá)時間，韓浩等[7]提出使用互相關(guān)約束來進(jìn)行長短時窗發(fā)的到達(dá)時間提取方法；林凡生[8]等先用長短時窗法確定信號到達(dá)時間的大致范圍，再利用沖擊波多個信號互相關(guān)精確提取沖擊波到達(dá)時間；何先龍等[9]將到達(dá)時間拾取轉(zhuǎn)化為能量變化率的提取，并用二次方自回歸對能量變化率曲線進(jìn)行精確提取到達(dá)時間；金澤龍[10]將極小時窗的協(xié)方差矩陣作為特征函數(shù)引入到STA/LTA中; 譚玉陽等[11]在單一的窗長識別因子的基礎(chǔ)上，引入了偏振度函數(shù)與邊緣檢測函數(shù)；喻志超等[12]利用互相關(guān)對各通道信號進(jìn)行時差校正，將所有信號的疊加道進(jìn)行全局互相關(guān)得到初至相對校正量，最后的初至?xí)r間由初始到達(dá)時間與相對校正量疊加得到；孟娟等[13]先進(jìn)性長短時窗到達(dá)時間前后2～3 s，進(jìn)行變分模態(tài)分解(VMD)后，提取每個分量的到時信息，進(jìn)行加權(quán)得到二次拾取的到時信息；趙揚(yáng)鋒等[15]用疊加的信噪指數(shù)代替信噪比，將信噪指數(shù)較高的信號作為輸入信號來提高信號到達(dá)時間提取精度。

針對上述問題，本文將信息論理論引入到了信號的到達(dá)時間拾取中，首先，為驗(yàn)證本文提出的拾取方法受噪聲影響較小，本文比較了不同噪聲條件下信號的拾取結(jié)果，其次，通過STA/LTA方法和本文方法對同一噪聲條件下的沖擊波進(jìn)行到達(dá)時間拾取，結(jié)果表明了基于時窗熵的沖擊波到達(dá)時間提取方法具有魯棒性，大大降低了噪聲對于信號初至波拾取的影響，本文提出的基于時窗熵的初至波拾取方法，既降低了噪聲對于信號到達(dá)時間拾取的影響，與STA/LTA方法相比又可以提高到達(dá)時間提取精度，解決了STA/LTA方法對于信號變化幅度小而發(fā)生漏拾的情況，避免了不必要的提取誤差，實(shí)現(xiàn)了更高精度的到達(dá)時間提取。

1 STA/LTA方法原理

長短時窗比值(STA/LTA)方法通過計(jì)算隨時間移動的長短時窗的平均能量比，選擇大于定義的閾值的峰值作為初至波的到達(dá)時間，STA與LTA的比率可以表示短時間內(nèi)信號的振幅、頻率和能量的異常變化。圖1描述了短時間平均值對長時間平均值方法的確定原理。

其表達(dá)式為式(1):

(1)

其中：i是當(dāng)前時間；STA(i)表示短時間窗信號的平均值；LTA(i)表示長時間窗信號的平均值；CF(xi)表示信號xi的特征函數(shù)；N代表短時間窗口的長度(即短時間窗口中信號數(shù)據(jù)點(diǎn)的總數(shù))；M是長時間窗口的長度(即長時間窗口中信號數(shù)據(jù)點(diǎn)的總數(shù)(M=5N))，窗長取決于信號的頻率，設(shè)置太短可能因?yàn)樵肼晫?dǎo)致錯誤拾取，設(shè)置太長可能會錯過初至點(diǎn)，窗長的選取一般為短時窗窗長為信號周期的2～3倍，長時窗的窗長一般為短時窗窗長的5倍。特征函數(shù)有信號的幅值、能量、頻率的變化，計(jì)算方式如下所示：

CF1=|z(i)|

CF2=|z(i)i|

CF3=|z(i)-z(i-1)|

CF4=|z(i)2-z(i-1)×z(i+1)|

(2)

上述式子中，i為采樣點(diǎn)數(shù)，前兩個特征函數(shù)表征時域特征，后兩個特征函數(shù)體現(xiàn)信號幅值與頻率的變化關(guān)系，一般在時域進(jìn)行初至波拾取，選取信號的能量特征作為特征函數(shù)。表達(dá)式如下式：

(3)

長短時窗可能會發(fā)生將噪聲誤認(rèn)為是突變點(diǎn)，為了降低噪聲對于初至波拾取的影響，更加精確得識別沖擊波初至點(diǎn)，劉曉紅等[4]在單獨(dú)的長短時窗的基礎(chǔ)上，加入了等窗長的識別因子，表達(dá)式如下：

(4)

結(jié)合兩種窗長的識別因子，有了一個新的初至波識別因子，其表達(dá)式如下:

(5)

P(i)確定初至點(diǎn)的范圍，Q(i)精確定位初至波到達(dá)時間。但是這種方法也知識能夠粗略估計(jì)初至波到達(dá)的時間，賀銘[5]為了放大信號突變情況，將高階統(tǒng)計(jì)量引入到初至波拾取中，得到了更加精確的時間信息。

不論哪個特征函數(shù)的長短時窗方法拾取得到的初至點(diǎn)都與真實(shí)的初至波到達(dá)時間存在一定誤差，所以在STA/LTA基礎(chǔ)上引入高階統(tǒng)計(jì)量，高階統(tǒng)計(jì)量在隨機(jī)變量中有高階矩和高階累積量，在隨機(jī)過程中有高階矩、高階累積量、高階譜，高階統(tǒng)計(jì)量因?yàn)榭梢垣@得一階、二階所不能獲得的信息或者不滿意的結(jié)果而被廣泛應(yīng)用于雷達(dá)、聲納、通信、地球物理、故障診斷等領(lǐng)域的信號處理問題中，典型的問題有自適應(yīng)估計(jì)濾波、陣列信號處理、時間序列分析、圖像處理等。高階統(tǒng)計(jì)量在微弱信號識別方面具有優(yōu)勢，高階統(tǒng)計(jì)量在信號處理方面的應(yīng)用，多為三階偏斜度和四階峰度，放大了信號的突變信息，其中，偏斜度是信號對于均值偏離的度量，體現(xiàn)了信息偏斜的程度，偏斜度的表達(dá)式如下式：

(6)

峰度表示信號中分布集中的程度，在峰度曲線中值越大的部分代表信號突變情況越明顯，峰度的表達(dá)式如下式：

(7)

將式(5)的識別因子與高階統(tǒng)計(jì)量的峰度結(jié)合起來，對于識別信號的突變更加敏銳，從而達(dá)到提高初至波拾取精度的目的，得到一個新的長短時窗識別因子，公式如下所示：

R=P(i)·Q(i)·K(i)

(8)

2 基于時窗熵的沖擊波到時提取理論

STA/LTA方法簡單快捷，但是STA/LTA在低信噪比情況下容易出現(xiàn)誤拾、漏拾的情況，拾取精度不高。沖擊波初至波到時拾取影響沖擊波超壓場重建精度，STA/LTA方法對于低信噪比的沖擊波信號到時提取方面拾取精度較低，會對后續(xù)應(yīng)用造成較大失誤，所以本文引入了信息論理論。在信息論理論中，信號與噪聲看作兩個獨(dú)立的變量，二者的互信息為零，因此在利用信號信息量進(jìn)行初至波拾取時，可以忽略噪聲對拾取結(jié)果的影響。

2.1 信息論理論

熵的概念最初被在熱力學(xué)中被提出[19]，是熱力學(xué)系統(tǒng)中表征物質(zhì)體系混亂程度的度量，熱熵越大，表明物理系統(tǒng)可能的圍觀狀態(tài)越多，從微觀上看，系統(tǒng)變化多端，沒有秩序。香農(nóng)為表示信息中的不確定性的度量，解決信息的量化度量問題，將熵的概念引入到信息論中，一般稱為香農(nóng)熵或者信息熵，信息熵[20]表示信息中的不確定性，信息中包含的有用信息量越多，信息熵越大。

在信息論中，隨機(jī)變量X的信息量可以表示為：

I(x)=-logp(x)

(9)

式中，I(x)是一個隨機(jī)變量，p(x)表示隨機(jī)變量的概率分布，采用信息熵來衡量變量集{x1,x2,x3,...,xn}的平均信息量，信息熵[21]的定義是為：

(10)

另一種信息度量的形式互信息I(X,Y)的定義為：

(11)

其中：fXY是隨機(jī)變量X和另外一個隨機(jī)變量Y的聯(lián)合概率密度函數(shù)，fX(x)、fY(y)是變量X和變量Y的邊際概率密度函數(shù)互信息可以表示為：

I(X:Y)=H(Y)-H(Y|X)=H(X)+H(Y)-H(XY)

(12)

其中：H(X,Y)為兩個隨機(jī)信號的交叉熵，H(X|Y)和H(Y|X)分別為信號X對于信號Y的條件熵和信號Y對于信號X的條件熵，H(X)、H(Y)分別為信號X和信號Y的信息熵，I(X;Y)為信號X和信號Y的互信息；由圖幾可以看出，互信息表示兩個隨機(jī)信號所攜帶的信息的共同部分，當(dāng)兩個信號相互獨(dú)立時，互信息為零。在沖擊波到時提取的問題中，沖擊波和背景噪聲可以看作是兩個相互獨(dú)立的變量，互信息為零同樣適用，所以無論采集到的沖擊波信號帶有多大的噪聲，都不會影響沖擊波到達(dá)時間拾取。

2.2 時窗熵理論

信息熵用來衡量變量集的不確定性，針對中遠(yuǎn)場沖擊波信號微弱，底噪大，上升沿緩慢等特點(diǎn)，本文參照STA/LTA原理結(jié)合互信息理論，提出了基于時窗熵的沖擊波到時提取方法，該方法解決了STA/LTA方法對于信號變化幅度小而發(fā)生漏拾的情況，避免了不必要的提取誤差，實(shí)現(xiàn)了更高精度的沖擊波到達(dá)時間的提取，時窗熵反映的是時窗內(nèi)信號的不確定性，兩個時窗內(nèi)信息越不一致，兩時窗的互信息越小。時窗內(nèi)的信息包含沖擊波信號和噪聲，沖擊波和背景噪聲可以看作是兩個相互獨(dú)立的變量，互信息為零，所以噪聲信號不會對到時提取的結(jié)果造成干擾。該方法的原理圖如圖2所示。

圖3 沖擊波信號及添加不同噪聲后的信號

圖2 基于時窗熵的沖擊波到時提取示意圖

其表達(dá)式為式(13):

h(i)=-∑i

(13)

圖4 沖擊波到達(dá)時間提取結(jié)果

式中，i為第i個時間點(diǎn)，N為時窗長度，pj為短時窗的信息量，pk為長時窗的信息量。當(dāng)短時窗滑動到初至點(diǎn)時，h(i)開始減小，當(dāng)長時窗到達(dá)初至點(diǎn)時，h(i)達(dá)到最大值，第i個點(diǎn)即為本文要提取的到達(dá)時間。

基于時窗熵的沖擊波到達(dá)時間提取方法將STA/LTA方法里面的不等長窗長換為了兩個等窗長的時窗，計(jì)算兩個時窗內(nèi)信號的互信息熵，時窗a位于式(13)中i點(diǎn)的左面，長度為N，時窗b位于式(13)中i點(diǎn)的右面，長度為N，隨著時窗的滑動，當(dāng)i點(diǎn)到達(dá)沖擊波到達(dá)時間的起點(diǎn)后，兩個時窗內(nèi)信號的互信息互不相關(guān)，互信息熵達(dá)到最小，這個最小的點(diǎn)就是本文要提取的到時信息。

基于時窗熵的沖擊波到達(dá)時間提取的步驟如下。

步驟一：輸入采集到的沖擊波信號;

步驟二：對沖擊波信號進(jìn)行頻譜分析，得到?jīng)_擊波信號的瞬時頻率;

步驟三：根據(jù)沖擊波信號的順勢頻率來確定時窗的長度;

步驟四：計(jì)算每一時刻的時窗熵;

步驟五：提取沖擊波的到達(dá)時間。

3 實(shí)驗(yàn)仿真與結(jié)果分析

首先，為驗(yàn)證本文提出的基于時窗熵的沖擊波到達(dá)時間提取的方法，生成當(dāng)量為20 kg，距離炸點(diǎn)3.7 m處的沖擊波信號，該沖擊波信號頻率為3.5 kHz，對沖擊波信號添加信噪比分別為30 dB、60 dB和75 dB的高斯白噪聲，按照上述步驟分別對沖擊波信號和加噪信號進(jìn)行基于時窗熵的沖擊波到達(dá)時間進(jìn)行提取，來驗(yàn)證本文提出的基于時窗熵的沖擊波到達(dá)時間的提取結(jié)果具有魯棒性，沖擊波信號和添加不同噪聲后的三組信號如圖3，先對圖3的四組信號進(jìn)行基于時窗熵的到達(dá)時間提取，提取結(jié)果如圖4。

圖5 沖擊波到達(dá)時間提取結(jié)果

圖4分別為沖擊波信號和添加30 dB、60 dB、75 dB噪聲信號的到達(dá)時間提取結(jié)果。

對添加不同噪聲的信號進(jìn)行基于時窗熵的信號到時提取可以得出其拾取誤差，添加不同噪聲的沖擊波到達(dá)時間提取的誤差如表1所示。

表1 不同噪聲條件下本文方法的提取結(jié)果及誤差

由表1可以看出，本文方法對于未加噪信號的到達(dá)時間提取結(jié)果為0.167 ms；添加30 dB噪聲的信號本文方法的拾取結(jié)果為0.166 9 ms，拾取誤差為0.65%；添加60 dB噪聲的信號本文方法的拾取結(jié)果為0.166 7 ms，拾取誤差為0.77%；添加75 dB噪聲的信號本文方法的提取結(jié)果為0.166 7 ms，拾取誤差為0.77%；本文通過比較沖擊波信號添加不同噪聲的信號通過基于時窗熵的到達(dá)時間的提取誤差，可以看出拾取結(jié)果的誤差精度基本保持在0.77%，可得出噪聲對基于時窗熵的到達(dá)時間提取結(jié)果干擾較小，基于時窗熵的沖擊波到達(dá)時間提取的方法具有魯棒性。

其次，為驗(yàn)證基于時窗熵的沖擊波到達(dá)時間提取方法在沖擊波到達(dá)時間拾取精度方面同樣具有優(yōu)勢，本文比較了同一噪聲條件下的信號分別采用STA/LTA方法和本文方法對與沖擊波到達(dá)時間提取精度，同樣對未添加噪聲的沖擊波信號和分別添加30 dB、60 dB和75 dB的高斯白噪聲的沖擊波信號使用STA/LTA方法進(jìn)行沖擊波到達(dá)時間提取的仿真實(shí)驗(yàn)。

圖5為使用STA/LTA方法對于沖擊波信號和添加30 dB、60 dB、75 dB噪聲信號的到達(dá)時間提取結(jié)果。

對添加不同噪聲的信號進(jìn)行STA/LTA方法的信號到時提取可以得出其拾取誤差，添加不同噪聲的沖擊波到達(dá)時間提取的誤差如表2所示。

表2 2種方法在同一噪聲下的拾取結(jié)果

由表2可以看出，在添加噪聲為30 dB條件下，STA/LTA方法對于沖擊波到達(dá)時間提取結(jié)果的誤差為2.8%，本文提出的基于時窗熵的沖擊波到達(dá)時間提取方法的拾取結(jié)果的誤差為0.65%，在添加噪聲為60 dB條件下，STA/LTA方法對于沖擊波到達(dá)時間提取結(jié)果的誤差為4.8%，本文提出的基于時窗熵的沖擊波到達(dá)時間拾取方法的拾取結(jié)果的誤差為0.77%，在添加噪聲為75 dB條件下，STA/LTA方法對于沖擊波到達(dá)時間拾取結(jié)果的誤差為6.6%，本文提出的基于時窗熵的沖擊波到達(dá)時間拾取方法的拾取結(jié)果的誤差為0.77%。

結(jié)合表1和表2，可以明顯看出本文提出的基于時窗熵的沖擊波到達(dá)時間提取方法在抗噪干擾方面以及拾取結(jié)果精度方面都具有明顯的改進(jìn)，在后續(xù)沖擊波場重建時的先驗(yàn)信息為沖擊波的到達(dá)時間信息，由于沖擊波傳播速度快，不同位置處的沖擊波的到達(dá)時間相差僅有幾毫秒，所以提高沖擊波到達(dá)時間提取的精確度至關(guān)重要，本文提出的基于時窗熵的沖擊波到達(dá)時間提取方法，該方法在提取到達(dá)時間的結(jié)果具有魯棒性，既降低了噪聲對于沖擊波信號到達(dá)時間提取的干擾，與STA/LTA方法相比又可以實(shí)現(xiàn)對于沖擊波到達(dá)時間的高精度提取，從而實(shí)現(xiàn)了更高精度的沖擊波到達(dá)時間的提取，使提取結(jié)果具有更高的提取精度，能夠?yàn)榇髤^(qū)域的沖擊波超壓場高精度重建提供有效的到達(dá)時間特征參數(shù)，在高價值彈藥毀傷效能參數(shù)中具有一定的理論意義和工程使用價值。

4 結(jié)束語

為滿足大威力炸藥的全區(qū)域重建，根據(jù)中遠(yuǎn)場沖擊波底噪大，信號微弱等特點(diǎn)，本文結(jié)合信息論與STA/LTA方法，提出基于時窗熵的沖擊波到達(dá)時間的提取方法，根據(jù)信息論理論，沖擊波信號和噪聲相互獨(dú)立，兩者的互信息為零，通過仿真實(shí)驗(yàn)結(jié)果表明，在不同噪聲條件下，本文提出的沖擊波到達(dá)時間拾取方法具有魯棒性，噪聲對于沖擊波到達(dá)時間的提取結(jié)果影響較小，提取誤差保持在0.77%，避免了對于沖擊波到達(dá)時間的提取過程中噪聲對于提取結(jié)果的干擾；其次，在同一噪聲條件下，對比了STA/LTA方法和本文提出的基于時窗熵的到達(dá)時間提取方法在進(jìn)行沖擊波到達(dá)時間提取的結(jié)果，結(jié)果表明本文提出的沖擊波到達(dá)時間提取方法對于提取結(jié)果的誤差精度方面高于STA/LTA方法，提取精度高于STA/LTA約6%，解決了STA/LTA方法對于信號變化幅度小而發(fā)生漏拾的情況，避免了不必要的提取誤差，實(shí)現(xiàn)了更高精度的到達(dá)時間提取，有利于后續(xù)的炸點(diǎn)定位以及沖擊波場重建，有利于重建模型的正確構(gòu)建。