基于時窗熵的沖擊波到時提取方法研究

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

(1.中北大學 信息探測與處理山西省重點實驗室，太原 030051；2.山東省軍區數據信息室，濟南 250099； 3.湖南云箭集團有限公司， 長沙 410000)

0 引言

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

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

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

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

1 STA/LTA方法原理

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

其表達式為式(1):

(1)

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

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為采樣點數，前兩個特征函數表征時域特征，后兩個特征函數體現信號幅值與頻率的變化關系，一般在時域進行初至波拾取，選取信號的能量特征作為特征函數。表達式如下式：

(3)

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

(4)

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

(5)

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

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

(6)

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

(7)

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

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

(8)

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

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

2.1 信息論理論

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

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

I(x)=-logp(x)

(9)

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

(10)

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

(11)

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

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

(12)

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

2.2 時窗熵理論

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

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

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

其表達式為式(13):

h(i)=-∑i

(13)

圖4 沖擊波到達時間提取結果

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

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

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

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

步驟二：對沖擊波信號進行頻譜分析，得到沖擊波信號的瞬時頻率;

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

步驟四：計算每一時刻的時窗熵;

步驟五：提取沖擊波的到達時間。

3 實驗仿真與結果分析

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

圖5 沖擊波到達時間提取結果

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

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

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

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

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

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

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

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

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

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

4 結束語

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