基于时窗熵的冲击波到时提取方法研究

刘晓佳，李 剑，刘代劲， 魏晓曼， 孔庆珊，金 艳

(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方法对于信号变化幅度小而发生漏拾的情况，避免了不必要的提取误差，实现了更高精度的到达时间提取，有利于后续的炸点定位以及冲击波场重建，有利于重建模型的正确构建。