火炮射击冲击加速度信号提取方法研究*

张宇曦，贾彦斌，张凤云，闫 琦，白 羽

（1.北方自动控制技术研究所，太原 030006；2.解放军63961 部队，北京 100012）

0 引言

火炮射击过程中，身管受到的冲击载荷具有输入间断性、作用时间极短、包含大量的高幅值高频成分等特点，会使火炮零部件结构发生形变和损伤，严重影响火炮结构的安全、寿命和性能。火炮结构在发生形变和损伤后，在大多数情况下表现为特定频率上的冲击，这种冲击相对火炮发射的主要冲击特征很弱，容易被淹没，提取该信号可以用于结构损伤识别、故障诊断等领域。

加速度信号的响应最能反映火炮发射过程的非线性、非平稳特征［1］。而加速度信号的强度和频率都会随时间变化，所以需选取时频分析方法研究火炮发射过程的信号变化情况。

传统的时频分析方法需预置参数，缺乏自适应特性［2］。针对传统方法缺乏自适应性的问题，Gilles提出了基于经验模态分解（EMD）的自适应性和小波分析的理论框架的经验小波变换（EWT）［3］。EWT通过对信号频谱的自适应划分，构造合适的正交小波滤波器组，提取具有紧支撑傅里叶频谱的AM-FM 成分，然后对分解出的信号进行Hilbert 变换得到瞬时幅值和瞬时频率，以表征信号的时频分布情况。EWT 解决了经验模态分解方法模态分量正交性无理论支撑、易出现模态混叠、虚假模态、计算量大、计算时间长等问题，相比经验模态分解方法分解效率和准确性更高［4］。

然而EWT 也需要预置分解的模态数，模态数的选择直接影响傅里叶频谱的分割结果。另外傅里叶频谱对冲击特征相对较弱的成分表现力较差，也会影响频谱分割。频谱的分割结果最终影响EWT信号分析效果。本文针对频谱分割结果影响EWT分析效果的问题，提出了一种根据能量算子得到的能量谱分割频谱的方法，对EWT 进行改进，最后通过分析叠加已知信号后的火炮发射加速度数据，验证该方法提取冲击特征相对较弱信号的有效性。

1 经验小波变换原理

EWT 首先要对傅里叶频谱进行划分，根据预置的分解模态数N，取傅里叶频谱由大到小排序为N的极大值所对应的频率为中心频率，通过N+1 条边界线将范围为［0，π］的傅里叶支撑域分为N 个各包含1 个中心频率的连续区域。边界频率定义以ωn为中心、宽度为的过渡段，如图1 所示的阴影区。

图1 频谱分割图

傅里叶频谱分割完成后，根据Meyer 小波的构造方法构建经验尺度函数φn（t）和经验小波函数ψn（t），其傅里叶变换分别如下：

其中，

经验小波函数和尺度函数与信号的内积分别为N 个细节系数和1 个近似系数，N 个细节系数和经验尺度函数的卷积为N 个细节经验模态，1 个近似系数和经验尺度函数的卷积为1 个近似经验模态，这N+1 个经验模态的总和为重新构建的信号。最后对每个经验模态进行Hilbert 变换，得到信号频率和幅值随时间变化的情况。

2 能量算子原理

离散信号的瞬时幅值A（n）和瞬时频率ω（n）估计值如下：

信号在噪声干扰下，傅里叶变换得到的幅频特性结果突出表现了干扰频率，而无信号包含信息频率［9-11］。能量算子对窄带信号的瞬时能量可进行追踪和分析，能够很好地反应信号的瞬时能量变化，对瞬态冲击成分的表现有增强效果［12］，且不同于傅里叶变换的积分运算，是基于求导的运算，计算复杂性低。

3 基于经验小波变换和能量算子的火炮冲击加速度信号信息提取方法

传统经验小波变换中，傅里叶频谱分割容易忽略冲击响应特征相对微弱的成分，而能量算子对瞬态冲击成分的表现有增强效果。针对上述问题，提出了一种基于经验小波变换和能量算子的火炮冲击加速度信号信息提取方法。首先对待测信号进行能量算子解调得到能量谱，再对能量谱进行阈值修剪，能量谱中可识别的幅值极大值对应的频率为主频率分量，用于确定EWT 的模态数和中心频率。具体步骤如下：

1）获取火炮冲击信号x（t），将x（t）信号用Teager 能量算子进行解算，得到信号能量谱，令A（ω）代表能量谱幅值，ω 为频率。

2）查找能量谱局部的极大值和极小值。

3）计算阈值T［6］：

其中，Am和AM分别为能量谱中的最大值和最小值，r 与信噪比成反比，常数C 控制r 取值在［0，1］范围内。

4）使用阈值T 截取能量谱中大于T 的极大值和极小值，并定位其对应的频率。极大值的个数即为经验小波变换分解的模态数N，极大值对应的频率定义为该模态的中心频率。极小值对应的频率则作为模态的边界频率，对频谱进行分割。

5）根据频谱分割结果对x（t）作经验小波变换分解。

对火炮冲击加速度信号信息提取流程如图2所示。

图2 火炮冲击加速度信号信息提取流程图

4 实验验证

本实验模拟火炮发射造成结构故障的情况，原始信号为一次火炮发射时的实测冲击加速度数据，将其叠加信号以模拟产生的冲击响应特征较弱的信号，对叠加信号分析以验证基于能量算子和EWT方法的有效性。叠加的信号y（t）为20 ms 后发生的幅值为10 g、频率为1 500 Hz 的正弦函数和26 ms 后发生的幅值为20 g、频率为100 Hz 的正弦函数。火炮载体坐标系定义如图3 所示，实验中的实测火炮数据x（1t）为MEMS 加速度传感器测试x 轴方向上的一次冲击加速度信号，采样频率为9 kHz，记录时长为40 ms。实验信号的时域波形如图4 所示。

图3 火炮载体坐标系示意图

图4 实验信号时域波形

根据MEMS 加速度传感器采集的冲击加速度x1（t）叠加y（t）后的时域波形分析，火炮冲击加速度在大约8ms 内无明显波动，之后10.028 ms 剧烈变化至-480.3 g，然后在10.080 ms 迅速变化至472.8 g，17.238 ms 后加速度开始逐渐稳定。本次实验数据记录时长为40 ms，但在40 ms 时冲击加速度还未稳定为0 g。

对x2（t）作能量算子解算，得到的能量序列结果如图5 所示。

图5 能量算子解算结果

由解算结果分析，前8 ms 能量有小幅震荡，8 ms～12 ms 能量激增至峰值后迅速回落，12 ms 以后能量振荡且幅度较大，20 ms 后能量逐渐减少，能量主要集中在8 ms～12 ms。对比x2（t）波形，能量谱更明显地反映出了前8 ms 和后20 ms 信号中的瞬态冲击成分。

对x2（t）进行FFT 和阈值计算。假设实验选取的加速度传感器信噪比大于85 dB，令C=10，则r 取0.1。傅里叶分析结果和阈值计算结果如图6。从图中可以看出1 500 Hz 的信号小于阈值，被其他频率成分淹没。

图6 傅里叶频谱图

对能量序列作FFT 得到能量谱，然后按流程计算阈值，最后分割频谱，得到的结果如图7，其中黑色虚线代表阈值，红色竖线代表边界频率。图7 能量谱相比傅里叶频谱图突显了更多相对弹丸主要冲击响应特征较弱的频率分量。

图7 能量谱分割结果

图7 中，8 条边界将能量谱分割为9 个部分，即9 个模态，其中能量谱中大于阈值的7 个极大值为7 个主模态的中心频率。根据上述能量谱分割结果对x2（t）进行经验小波分解，分解结果如下页图8 所示。

图8 改进EWT 分解结果

图9 x1（t）传统EWT 时频谱

图10 x2（t）传统EWT 时频谱

图11 x1（t）改进EWT 时频谱

图12 x2（t）改进EWT 时频谱

将x1（t）和x2（t）的信号分别进行传统EWT 和改进EWT 分解，得到其时频谱如图9～图12 所示。对比传统EWT 分解后的时频图，能看出加y（t）后信号的时频图有新的频率成分出现，但是无法辨别其幅值、频率和出现时刻。通过对图11 和图12 的对比分析，改进EWT 时频谱能更清晰地体现y（t）的成分。为了清晰地分析新的频率成分，将图12 局部放大如图13 可以观察到，20 ms 开始出现频率约为1 500 Hz、幅值基本稳定为10 g 的信号，26 ms 左右开始出现频率约为100 Hz、幅值基本稳定为20 g 的信号，并且这两种频率成分出现后持续存在于频谱中。图9 和图11 对比，传统方法的时频图中幅值差异过大，瞬时幅值最大值对应的频率成分和时间点极少，且瞬时幅值大部分在250 g 以下，较小幅值的瞬态信号观测困难。综上所述，传统EWT 方法时频谱整体观测难度更大，基本无法观测实验中叠加y（t）的变化情况，说明改进EWT 比传统EWT 对冲击加速度信号特征的提取效果更好。

5 结论

图13 改进EWT 时频谱细节图

本文采用一种通过能量算子输出指导频谱分割的改进EWT 方法对一次叠加已知信号的实测火炮射击冲击加速度信号进行分析，解决了傅里叶频谱中相对火炮发射主要冲击特征很弱的故障频率成分被淹没的问题。改进EWT 方法与传统EWT 方法相比，能更精确、清晰地表现信号的时频特征，提取信号特征的能力更强。能量算子能加强冲击特征，但同时对噪声敏感，所以需要对信号去噪方法进一步研究。