火工冲击响应数据时频分析方法与流程

本发明涉及火工冲击信号处理技术领域，特别是涉及火工冲击响应数据时频分析方法。

背景技术：

航天器火工冲击环境是由星箭分离、部组件展开等工作过程中的火工品起爆引起的作用于结构上的瞬态冲击响应，具有瞬态、高频、高量级的特点，是航天器在全生命周期内经历的最苛刻的力学环境之一。为了提高航天器的工作性能和可靠性，需要精确模拟火工冲击环境，考核航天产品对火工冲击的耐受性。火工冲击环境模拟方法分为数值仿真和地面试验两大类。目前，对于火工品爆炸引起的高频结构响应，国内外尚不具备成熟、有效的分析技术，无法满足工程需求。

传统的冲击响应谱分析方法广泛应用于航天产品火工冲击环境地面模拟试验中。在鉴定冲击响应量级、规定结构对冲击环境的承受能力或确定设备级模拟冲击试验的输入谱时，它可以代替冲击的时域响应，常被认为是冲击破坏能力的衡量标准。但是，由于冲击响应谱仅仅考虑了冲击信号作用在单自由度系统上响应的峰值信息，忽略了冲击信号的时变特征，冲击响应谱与时域信号难以一一对应，也就造成地面试验与实际火工环境不完全等效的问题，在指导地面试验时具有一定的局限性。为了更精准地提取冲击信号中的有效信息，为地面试验的开展提供更好的理论依据，有必要提出一种适合火工冲击响应的数据分析方法，能刻画冲击响应随时间的变化细节。相对于冲击响应谱分析，时频分析方法能够表征信号能量和频率随时间变化的规律，可用于分析冲击响应数据的时变特性。

技术实现要素：

为了克服上述现有技术的不足，本发明实施例提供了一种基于变分模态分解和rihaczek分布的火工冲击响应数据时频分析方法，该方法能够精细刻画冲击响应的能量时频分布特征，更直观地表征冲击激发的结构模态响应分布规律，弥补了冲击响应谱分析方法的局限性。

本发明所要解决的技术问题是单纯依靠冲击响应谱分析的火工试验不完全等效，辅助实现航天产品真实火工冲击环境的精确模拟、测试和分析。

根据本发明实施例，提出一种火工冲击响应数据时频分析方法，包括以下步骤：s1，获取火工冲击响应离散加速度信号；s2，利用变分模态分解方法将信号分解为一系列从高频到低频分布排列的单模态冲击响应信号分量；s3，利用rihaczek函数计算各单模态冲击响应信号分量的时频分布函数，抑制单模态冲击响应信号分量在时域与频域方向上的交叉项干扰；s4，线性叠加各单模态冲击响应信号分量的时频分布函数并表征到二维时频联合域，得到火工冲击响应数据时频分布。

与现有技术相比，本发明具有以下有益效果中的至少一个：

(1)本发明的火工冲击时频分析方法，能够同时反映冲击响应的能量时频分布特征，能够更清晰地反映冲击响应的触发时间和次数，因而更直观地表征冲击激发的结构模态响应分布规律；

(2)有助于弥补冲击响应谱的不足，实现模拟试验方法与真实火工冲击在频域和时域上均具有较好的一致性；

(3)有助于提高地面试验模拟真实火工冲击环境的精准性，对航天器重量减轻、缓冲设计和冲击验收试验标准制定等提供参考意义。

附图说明

通过参照附图详细描述本发明的实施例，本发明将变得更加清楚，多个实施例被图示在附图中以用于说明性目的，并且决不应该被理解为限制实施例的范围。另外，不同的公开的实施例的各种特征可以组合以形成额外实施例，额外实施例是本公开的一部分，其中：

图1是根据本发明实施例的火工冲击响应数据时频分析方法结构框图；

图2是本发明实施例的原子仿真信号vmd分解结果图，其中图2(a)对应时域imf分量信号，图2(b)对应频域imf分量信号；

图3是本发明实施例的原子仿真信号vmd-rd时频分布图。

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。在说明书中，相同或相似的附图标号指示相同或相似的部件。下述参照附图对本发明实施方式的说明旨在对本发明的总体发明构思进行解释，而不应当理解为对本发明的一种限制。

在本发明的任何实施例被具体地描述之前，应该理解本文中公开的概念不将它们的应用限制到在以下描述中阐述或在以下附图中图示的构造的细部和构件的布置。在这些实施例中图示的概念能够以各种方式实践或执行。本文中使用的具体短语和术语为了便于描述，并且不应该被认为是限制性的。

下面将结合实施例中的附图，对实施例中的技术方案进行具体、清楚、完整地描述。

参照图1，根据本发明一个实施例的火工冲击响应数据时频分析方法，包括以下步骤：

s1，获取火工冲击响应离散加速度信号；

s2，利用变分模态分解方法将信号分解为一系列从高频到低频分布排列的单模态冲击响应信号分量；

s3，利用rihaczek函数计算各单模态冲击响应信号分量的时频分布函数，抑制单模态冲击响应信号分量在时域与频域方向上的交叉项干扰；

s4，线性叠加各单模态冲击响应信号分量的时频分布函数并表征到二维时频联合域，得到火工冲击响应数据时频分布。

其中，在s1步骤中，可以采用例如振动加速度传感器测量、声发射测量或者多普勒激光测量获取火工冲击响应离散加速度信号。

接着进行s2步骤，利用变分模态分解(variationalmodedecomposition，vmd)方法将步骤s1获得的信号分解为一系列从高频到低频分布排列的单模态冲击响应信号分量(intrinsicmodefunction，imf)；

例如，初始火工冲击响应信号设为x(t)，经模态分解得到k个单模态冲击响应信号分量ximf(t)，用下式表示为：

进一步地，在s3步骤中，计算各单模态冲击响应信号分量的rihaczek分布函数，可基于以下方法实现：

非线性时频分析方法是以维格纳变换(wigner-villedistribution，wvd)为基础发展来的，可将上述信号x(t)的wvd定义为

其中，*表示共轭，t表示时间，ω表示角速度，τ表示时间变化量。式中不包含任何的窗函数(为了减少频谱能量泄漏，可采用不同的截取函数对信号进行截断，截断函数称为窗函数)，然而，其他的非线性时频分析方法都可以看作是wvd的加窗形式(加窗可以实现时域局域化，还可以修正谱泄露问题)。为了正确描述火工冲击响应信号的局部能量分布，要求时频分布具有理想的时频局部聚集性(即在时域和频域分布上具有良好的集中性)。wvd时频分布的聚集性较高，但其瓶颈在于严重的交叉干扰项。wvd时频分布中每两个信号分量之间就会产生一个交叉项，交叉项也是传统非线性时频分析方法中不可避免的。

为抑制时频分析方法中的交叉项干扰，可对wvd进行时频平滑(即对时频分布求解进行加窗)，表示为下式：

其中φ(τ,θ)代表不同的核函数，选择不同的核函数，可以得到不同的交叉项抑制效果。当核函数为指数函数exp(iπωθ)时，又可以得到里哈契克时频分布(rihaczekdistribution，rd)，vmd分解所得各个imf分量信号的rihaczek分布可用下式描述为：

其中，*表示共轭，t表示时间，ω表示角速度，τ和θ分别表示时间和角度变化量。

进一步地，在s4步骤中，将各个单分量信号进行rihaczek分析和处理后消除了在时域方向上的交叉项干扰之后，将一系列分量结果线性叠加，可将初始火工冲击响应信号x(t)的时频分布定义为：

其中，*表示共轭，t表示时间，ω表示角速度，τ和θ分别表示时间和角度变化量。

所得vmd-rd函数即为火工冲击响应数据时频分布函数。根据上述实施例的时频分析方法，其充分利用了vmd方法对不同频率分量的解析能力和rihaczek分布对时域交叉项的抑制能力，能够有效提高时频聚集性，并消除交叉项成分的干扰。

参照图1，根据本发明一个较优实施例的火工冲击响应数据时频分析方法，在s1步骤中，采用振动加速度传感器测量获取火工冲击响应离散加速度信号，其中，冲击信号来自于单机冲击试验火箭端测点z方向的数据，通过对加速度信号进行采样获取离散加速度信号，采样率为10000hz，采样点数为5201。

进一步地，在s2步骤中，在利用变分模态分解(vmd)方法对信号进行分解的过程中，定义分解所得分量信号的功率谱熵，并将功率谱熵值大小作为评价指标，来判别vmd分解的彻底性。其中，

vmd算法的实现主要分为变分问题的构造和求解两部分；使用vmd算法对被分析信号x(t)进行分解，实质上就是求取k个单模态冲击响应信号分量ximf(t)，并保证各模态估计带宽之和最小；求解该变分问题过程中所使用的约束条件是：所有模态之和等于被分析信号。

求取各变分模态分量ximf(t)的hilbert变换函数，进而求得各ximf(t)分量对应的单边频谱，如公式(6)：

将指数函数exp(-jωit)加入各变分模态分量ximf(t)之中，调整分量对应的中心频率函数，如公式(7)：

利用l²范数对解调信号进行高斯平滑，估计各模态带宽，如公式(8)：

为了求解上述约束性变分问题的最优解，引入二次惩罚因子和拉格朗日乘法算子将约束性变分问题变为非约束性变分问题。采用admm算法迭代搜索求取上述扩展的拉格朗日函数的鞍点来解决上述非约束性变分问题。

vmd能够有效地抑制模态混叠现象，需要对分量个数k及惩罚函数α进行预先设定；为了实现这两个参数的自动寻优，引入shannon熵作为vmd分解效果的评价指标，因而提出冲击响应信号分量功率谱熵的概念，即对于每一个imf分量信号ximf(t)，其功率谱为

由于信号在时、频域变换的过程中能量是守恒的，因此，对于长度为n的imf分量信号，simf＝{s1,s2,…,sn}可以看作对信号的一种划分方式，分量信号对应的功率谱熵为

其中，qi为第i个功率谱值在整个谱中所占的百分比。

进一步地，利用粒子群优化(particleswarmoptimization，pso)方法对vmd分解效果影响较大的分量个数k及惩罚函数α两个参数进行自动寻优，包括以下步骤：

1)初始化所有粒子，设定种群数目n、迭代终止条件(最大迭代次数)和学习因子c1、c2；

2)对每个粒子的qi进行评价，以所求分量信号函数功率谱熵最小值为优化目标，直接利用目标函数值作为个体的适应度，计算其适应度函数值p(qi)，并与该粒子历史最优位置适应度函数值p^*(qi)进行比较；如果当前位置的适应度函数值比其历史最优位置适应度函数值更小，则令p^*(qi)＝p(qi)，相反，则保持p^*(qi)不变；

3)更新每个粒子的速度v(t)与位置s(t)，按如下公式：

vi(t+1)＝vi(t)+c1·rand·(pi(t)-si(t))+c2·rand·(pi^*(t)-si(t))(11)

si(t+1)＝vi(t+1)+si(t)(12)

4)判断是否达到该算法运行时的最大迭代数，满足迭代终止条件时结束该算法；否则，转到步骤2)。

为了验证参数优化策略的有效性，建立一个多分量仿真信号x(t)，信号长度为256，由5个gaussian原子复合而成，5个原子分别对应归一化频率的3个频率成分即0.1、0.25和0.4，在时域方向上原子分量中心分别位于采样点数60、130和200位置上，运用上述方法对仿真信号进行vmd分解，最优分量个数k取3，惩罚函数α取2000时，信号的imf分量功率谱熵值存在最小值。仿真信号包含3个频率分量，分别为10hz、25hz、40hz，目标函数在k＝3处取得最小值，k＜3时vmd分解不完全，imf分量中信号频率成分不够单一；k＞3时，随着k值增加，目标函数值随之增加，vmd出现了过分解现象。参数优选结果与仿真信号的实际情况吻合，证明了参数优化方法的有效性。vmd分解所得imf分量的时域、频域方向的波形如图2所示，由于仿真原子信号幅值已进行归一化处理，各分图(a)和(b)中纵坐标数值表示信号在时域或频域的相对幅值。

将该方法应用于实际火工冲击响应数据对信号x(t)进行vmd分解时，例如最优分量个数k取6，惩罚函数α取1800时，信号的imf分量功率谱熵值存在最小值。

进一步地，在s3步骤中，对上述s2步骤中的仿真信号分解得到的gaussian原子，计算其rihaczek分布函数，其得到的时频分布如图3所示，频率为10hz和40hz的imf分量中各包含2个原子分量，若将3个imf分量的wvd分布进行叠加，得到的时频分布中交叉项仍然无法消除。要进一步消除时域上的交叉项干扰，需要借助合适的核函数。rihaczek分布中的指数核函数能够有效消除wvd分布时域上的交叉项干扰，且时频分布中可以保持较高的时域分辨率，但是频域聚集性较差。而wvd分布使不同频率的原子分量更加清晰，时频局部聚集性较高。因此，vmd与rihaczek分布的结合恰好起到了取长补短的效果，既能够消除时、频域上的交叉项干扰，又能够保持较高的时频聚集性。

因而，在s4步骤中，将各个单分量信号进行rihaczek分析和处理后消除了在时域方向上的交叉项干扰之后，再将一些列分量结果线性叠加，便可得初始火工冲击响应信号的时频分布，其对应的时频分布函数如上式(5)所示。

根据本发明实施例的时频分析方法，从能量角度对火工冲击响应进行分析，建立基于vmd和rihaczek分布的火工冲击时频分析方法，该时频分析方法能够同时反映冲击响应的能量时频分布特征，更清晰地反映冲击响应的触发时间和次数，以及更直观地表征冲击激发的结构模态响应分布规律。在进行地面冲击试验时，能够使模拟试验方法与真实火工冲击在频域和时域上均具有较好的一致性。

利用上述时频分析技术，得到火工冲击响应数据时变规律的定量化描述，能够提高地面试验模拟真实火工冲击环境的精准性，对航天器重量减轻、缓冲设计和冲击验收试验标准制定等具有较重要的参考意义。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。

尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。