CN110319995A - 火工冲击响应数据时频谱图分析方法 - Google Patents

Abstract

本发明的实施例提出一种火工冲击响应数据时频分析方法，包括以下步骤：获取火工冲击响应离散加速度信号；计算火工冲击响应数据VMD‑Rihaczek时频分布函数，分析火工冲击激励非平稳时变特征；利用非负矩阵分解方法求解高维火工冲击响应时频图像矩阵的低维特征信息，量化火工冲击响应数据非平稳时变特征；基于冲击响应谱与冲击响应时频分布特征信息，利用D‑S证据融合理论建立火工冲击响应数据特征信息与火工冲击环境之间的映射关系。根据本发明的实施例提供的分析方法，在反映冲击损伤效应的同时进行火工冲击响应激励的精细辨识，并提出更加全面的火工冲击环境评价指标，从而为地面模拟试验的开展提供参考。

Description

火工冲击响应数据时频谱图分析方法

技术领域

本发明涉及数据处理技术领域，特别是涉及对火工冲击响应数据进行处理的方法。

背景技术

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

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

发明内容

为了克服上述现有技术的不足，本发明提出一种火工冲击响应数据时频谱图分析方法，该方法能够同时反映冲击响应的能量时频分布特征，更直观地表征冲击激发的结构模态响应分布规律，弥补了冲击响应谱分析方法的局限性。

根据本发明的实施例，其提供了一种火工冲击响应数据时频谱图分析方法，包括以下步骤：S1，获取火工冲击响应离散加速度信号；S2，计算火工冲击响应数据VMD-Rihaczek时频分布函数，分析火工冲击激励非平稳时变特征；S3，利用非负矩阵分解方法求解高维火工冲击响应时频图像矩阵的低维特征信息，量化火工冲击响应数据非平稳时变特征；S4，基于冲击响应谱与冲击响应时频分布特征信息，利用D-S证据融合理论建立火工冲击响应数据特征信息与火工冲击环境之间的映射关系。

通过以上步骤，将冲击响应谱关键信息和冲击响应时频分布特征信息有效结合，在反映冲击损伤效应的同时进行火工冲击响应激励的精细辨识，为地面模拟试验的开展提供参考，从而，实现航天产品真实火工冲击环境的精确模拟、测试和分析。

附图说明

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

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

具体实施方式

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

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

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

如图1所示，根据本发明的实施例提出一种火工冲击响应数据时频谱图分析方法，其包括以下步骤：S1，获取火工冲击响应离散加速度信号；S2，计算火工冲击响应数据VMD-Rihaczek时频分布函数，分析火工冲击激励非平稳时变特征；S3，利用非负矩阵分解方法求解高维火工冲击响应时频图像矩阵的低维特征信息，量化火工冲击响应数据非平稳时变特征；S4，基于冲击响应谱与冲击响应时频分布特征信息，利用D-S证据融合理论建立火工冲击响应数据特征信息与火工冲击环境之间的映射关系。

在步骤S1中，根据本发明优选的实施例，火工冲击响应离散加速度信号可以由振动加速度传感器测量得到。原始信号来自于单机冲击试验火箭端测点Z方向的数据，通过对加速度信号进行采样获取离散加速度信号。也可以采用其他方式获得信号，例如，声发射测量或者多普勒激光测量等。

在步骤S2中，根据本发明优选的实施例，在S2步骤中，基于变分模态分解方法和Rihaczek时频分布方法构造火工冲击响应数据VMD-Rihaczek时频分布。

在步骤S2中，进一步的，利用变分模态分解方法将火工冲击响应离散加速度信号分解为一系列从高频到低频分布排列的单模态冲击响应信号分量；计算各单模态火工冲击响应信号分量的Rihaczek分布函数并对结果进行线性叠加，表征到二维时频平面，进而得到火工冲击响应数据时频分布。

变分模态分解(variational mode decomposition，VMD)方法将火工冲击响应离散加速度信号分解为一系列从高频到低频分布排列的单模态冲击响应信号分量(Intrinsic Mode Function，IMF)。

初始信号x(t)经VMD分解得到K个单模态冲击响应信号分量x_IMF(t)，可表示为：

计算各单模态冲击响应信号分量的Rihaczek分布函数；

非线性时频分析方法是以维格纳变换(Wigner-Ville Distribution，WVD)为基础发展来的。信号x(t)的WVD定义为

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

对WVD进行时频平滑(即对时频分布求解进行加窗)，进而抑制交叉干扰项的时频分析方法形式如下：

其中φ(τ,θ)代表不同的核函数，选择不同的核函数，可以得到不同的交叉项抑制效果。当核函数为指数函数exp(iπωθ)时，又可以得到里哈契克时频分布(RihaczekDistribution，RD)

VMD分解所得各个IMF分量的Rihaczek分布可用公式描述为：

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

进一步地，将待分析的信号经VMD分解成一组单分量信号，对各个单分量信号单独进行Rihaczek分析和处理，消除在时域方向上的交叉干扰项，再将结果线性叠加。信号x(t)的时频分布定义为：

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

所得VMD-RD函数即为火工冲击响应数据时频分布函数。该方法充分利用了VMD方法对不同频率分量的解析能力和Rihaczek分布对时域交叉项的抑制能力，具有较高的时频聚集性，并有效消除交叉项成分的干扰。

如上所述，根据本发明优选的实施例，初始火工冲击响应离散加速度信号x(t)经变分模态分解方法分解得到K个单模态冲击响应信号分量x_IMF(t)，表示为

在一些实施例中，变分模态分解方法分解所得各个单模态冲击响应信号分量的Rihaczek分布具有如下关系：

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

在一些实施例中，将待分析的信号经变分模态分解方法分解成一组单分量信号，对各个单分量信号单独进行Rihaczek分析和处理，消除在时域方向上的交叉干扰项，再将结果线性叠加；

火工冲击响应信号x(t)的VMD-Rihaczek时频分布函数为：

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

在一些实施例中，利用分块并行火工冲击时频谱图非负矩阵分解方法获取时频图像的局部特征信息。

在一些实施例中，将火工冲击时频非负观测矩阵V_m×n(VMD-RD_x(t,ω))按列平均分解为块矩阵的形式V_m×n＝[V₁ V₂ … V_b]，其中，b为分块数，n₀＝n/b；

矩阵的分块并行火工冲击时频谱图非负矩阵分解形式表示为

利用分块并行火工冲击时频谱图非负矩阵分解算法(block non-negativematrix factorization,BNMF)对V进行分解，得到特征矩阵W，将火工冲击时频谱图向特征矩阵W进行投影，得到相应的特征系数。

根据本发明的一些实施例，可以利用非负矩阵分解方法(Non-negative MatrixFactorization，NMF)抽取时频图像主要特征；给定时频图像矩阵V_m×n，利用NMF算法可将矩阵V_m×n近似分解为两个矩阵W_m×r与H_r×n的乘积：

V_m×n＝W_m×r·H_r×n (6)

其中，称W_m×r为基矩阵，称H_r×n为系数矩阵，W、H均要求非负，每一个样本通过向基矩阵投影，可以得到对应的特征编码(或特征系数)，该特征编码可看做时频图像的特征参量，用于模式识别。进一步地，利用欧式距离作为V与W·H间的近似误差：

其对应优化问题为：

minE(V||WH),s.t.W,H≥0

相应迭代规则为：

NMF使用不同数量和不同标记的基向量的组合表示数据，所以可以抽取数据的多线性结构，具有一定的非线性数据分析能力。

进一步地，在本发明的一些实施例中，基于D-S证据理论综合考虑冲击响应谱与冲击响应时频分布特征信息；首先进行基本可信度函数的构造，在D-S证据理论中并没有给出基本可信度函数的一般形式，要根据具体问题构造具体形式；然后确定网络不确定度分配；

根据本发明的一些实施例，采用如下方法来求取网络的不确定度分配函数：确定距离和贴近度；相关性度量；基本概率分配。

通过以上步骤，将冲击响应谱关键信息和冲击响应时频分布特征信息有效结合，在反映冲击损伤效应的同时进行火工冲击响应激励的精细辨识，从而为地面模拟试验的开展提供参考。

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

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

S2，计算火工冲击响应数据高精度VMD-Rihaczek时频分布函数，分析火工冲击激励非平稳时变特征；

S3，利用非负矩阵分解方法求解高维火工冲击响应时频图像矩阵的低维特征信息，量化火工冲击响应数据非平稳时变特征；

S4，综合考虑冲击响应谱与冲击响应时频分布特征信息，利用D-S证据融合理论建立火工冲击响应数据特征信息与火工冲击环境之间的映射关系。在反映冲击损伤效应的同时进行火工冲击响应激励的精细辨识，并提出更加全面的火工冲击环境评价指标，从而为地面模拟试验的开展提供参考。

利用VMD-Rihaczek分析法进行火工冲击响应离散加速度信号的时频分析；利用变分模态分解(variational mode decomposition，VMD)方法将信号分解为一系列从高频到低频分布排列的单模态冲击响应信号分量；以VMD分解所得分量信号的功率谱熵值大小为评价指标，判别VMD分解的彻底性。

初始信号x(t)经VMD分解得到K个单模态冲击响应信号分量x_IMF(t)，如公式(1)所示。

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

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

将指数函数exp(-jω_it)加入各变分模态分量x_IMF(t)之中，调整分量对应的中心频率函数：

利用L²范数对解调信号进行高斯平滑，估计各模态带宽：

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

计算各单模态冲击响应信号分量(Intrinsic Mode Function，IMF)的Rihaczek分布函数。Rihaczek分布中的指数核函数能够有效消除时域上的交叉项干扰，时频分布中可以保持较高的时域分辨率，但是频域聚集性很差。VMD与Rihaczek分布的结合恰好起到了取长补短的效果，既能够消除时、频域上的交叉干扰项，又能够保持较高的时频聚集性。

对各个单模态冲击响应信号分量单独进行Rihaczek分析和处理，消除在时域方向上的交叉干扰项，再将结果线性叠加，得到的VMD-RD函数即为火工冲击响应数据时频分布函数。该方法充分利用了VMD方法对不同频率分量的解析能力和Rihaczek分布对时域交叉项的抑制能力，具有较高的时频聚集性，并有效消除交叉项成分的干扰。

在本发明的实施例中，可以采用如下方式，获得原始信号，进一步获得数据：随着火工品起爆t＝0.01s时，产生了极宽频带[0-10000]Hz下的宽频激励。紧随其后t＝0.012时，分离螺母的应变能释放产生了二次冲击，应变能释放的冲击能量集中在[3000-5500]Hz内。两次激励之后，结构的模态振动被激发，冲击响应迅速衰减。

在本发明的实施例中，步骤S3，利用非负矩阵分解方法(Non-negative MatrixFactorization，NMF)抽取时频图像主要特征；给定时频图像矩阵V_m×n，利用NMF算法可将矩阵V_m×n近似分解为两个矩阵W_m×r与H_r×n的乘积，设定n为数据样本的维数；m为集合中数据样本的个数；r为特征维数，一般情况下，r远小于n且满足r(m+n)＜mn。

利用欧式距离作为V与W·H间的近似误差，计算出非负矩阵分解编码的投影基矩阵。将原始图像矩阵向基矩阵投影，可获得特征编码。原始彩色时频图像转化为灰度图像后的维度为420×560，NMF编码值被限定于0～255的区间内，所选择的特征维度值r为64，样本个数为30。编码矩阵每一行为一个数据样本，由编码矩阵可见所提取特征向量类内差异小，可有效表征时频图像特征状态。

在本发明的一些实施例中，在步骤S4，基于D-S证据理论综合考虑冲击响应谱与冲击响应时频分布特征信息；

D-S证据理论假设有n个互相排斥的原始子命题存在，其类型为a₁,a₂,…,a_n，可记为集合{a₁,a₂,…,a_n}，这个命题集组成了整个假想空间，称之为辨识框架，用Θ表示。Bayes推理是对其中的n个元素运算，而D-S理论则是对2ⁿ个命题(包括原始子命题所有可能的并命题)进行运算。

首先进行基本可信度函数的构造，在D-S证据理论中并没有给出基本可信度函数的一般形式，要根据具体问题构造具体形式；

设第i个特征量由诊断参数进行局部诊断的可靠性为A_i，那么它对应的在本证据的基础上对j状态的可信度分配为：

m_i(Θ)＝1-A_i(i＝1,2,…,p) (14)

式中，m_i(j)代表第i个证据对j状态的可信度分配，m_i(Θ)表示根据第i个证据不能确定的可信度分配，即不能判断哪种状态发生的可能性(也就是各自状态都有可能发生)。然后，可进一步求取信度函数Bel和似真度函数pl，根据Bel和pl的值来判断系统所在状态的可能性。

然后确定网络不确定度分配：

确定距离和贴近度，设网络的输出有M个节点，即对应M种状态(正常与故障情况)，{Y_j}为网络M种状态对应的理想输出Y_j＝{y_j1,y_j2,…,y_jM},(j＝1,2,…,M)，{X_i}为网络对已知的M种状态验证样本经网络计算的实际输出X_i＝{x_i1,x_i2,…,x_iM},(i＝1,2,…,M)，二者存在一一对应关系。取{X_i}中一个实际输出X_k(k＝1,2,…,M)。则实际输出与标准理想输出的Manhattan距离为：

距离表示了实际输出与理想输出的贴近程度。距离越大，则贴近度越低；反之，则贴近程度越高；

相关性度量；

定义实际输出与理想输出的相关系数如下：

基本概率分配；根据相关性的定义求取X_k的基本可信度分配m_k(F_j)以及不确定描述m_k(Θ)。

通过以上方步骤，将冲击响应谱关键信息和冲击响应时频分布特征信息有效结合，在反映冲击损伤效应的同时进行火工冲击响应激励的精细辨识，从而为地面模拟试验的开展提供参考。

根据本发明实施例的火工冲击响应数据智能时频谱图分析方法，利用所提出的先进时频分析和图像处理技术，得到火工冲击响应数据时变规律的定量化精细描述，与冲击响应谱相结合，提出能够全面反映火工冲击环境的评价指标，提高地面试验模拟真实火工冲击环境的精准性，具有较强的工程背景，对航天器重量减轻、缓冲设计和冲击验收试验标准制定有较大的参考意义。

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

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

Claims (8)

1.一种火工冲击响应数据时频谱图分析方法，其特征在于，包括以下步骤：

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

S2，计算火工冲击响应数据VMD-Rihaczek时频分布函数，分析火工冲击激励非平稳时变特征；

S3，利用非负矩阵分解方法求解高维火工冲击响应时频图像矩阵的低维特征信息，量化火工冲击响应数据非平稳时变特征；

S4，基于冲击响应谱与冲击响应时频分布特征信息，利用D-S证据融合理论建立火工冲击响应数据特征信息与火工冲击环境之间的映射关系。

2.根据权利要求1所述的火工冲击响应数据时频谱图分析方法，其特征在于：

在S2步骤中，基于变分模态分解方法和Rihaczek时频分布方法构造火工冲击响应数据VMD-Rihaczek时频分布。

3.根据权利要求2所述的火工冲击响应数据时频谱图分析方法，其特征在于：

在S2步骤中，利用变分模态分解方法将火工冲击响应离散加速度信号分解为一系列从高频到低频分布排列的单模态冲击响应信号分量；

计算各单模态火工冲击响应信号分量的Rihaczek分布函数并对结果进行线性叠加，表征到二维时频平面，进而得到火工冲击响应数据时频分布。

4.根据权利要求3所述的火工冲击响应数据时频谱图分析方法，其特征在于：初始火工冲击响应离散加速度信号x(t)经变分模态分解方法分解得到K个单模态冲击响应信号分量x_IMF(t)，表示为

5.根据权利要求3所述的火工冲击响应数据时频谱图分析方法，其特征在于：

变分模态分解方法分解所得各个单模态冲击响应信号分量的Rihaczek分布具有如下关系：

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

6.根据权利要求3所述的火工冲击响应数据时频谱图分析方法，其特征在于：将待分析的信号经变分模态分解方法分解成一组单分量信号，对各个单分量信号单独进行Rihaczek分析和处理，消除在时域方向上的交叉干扰项，再将结果线性叠加；

火工冲击响应信号x(t)的VMD-Rihaczek时频分布函数为：

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

7.根据权利要求1所述的火工冲击响应数据时频谱图分析方法，其特征在于：利用分块并行火工冲击时频谱图非负矩阵分解方法获取时频图像的局部特征信息。

8.根据权利要求7所述的火工冲击响应数据时频谱图分析方法，其特征在于：将火工冲击时频非负观测矩阵V_m×n(VMD-RD_x(t,ω))按列平均分解为块矩阵的形式V_m×n＝[V₁ V₂ …V_b]，其中b为分块数，n₀＝n/b；

矩阵的分块并行火工冲击时频谱图非负矩阵分解形式表示为

利用分块并行火工冲击时频谱图非负矩阵分解算法对V进行分解，得到特征矩阵W，将火工冲击时频谱图向特征矩阵W进行投影，得到相应的特征系数。