CN111931288A - 一种火工冲击响应时频域分布相似性测度方法 - Google Patents

Abstract

本发明提出一种火工冲击响应时频域分布相似性测度方法，涉及火工冲击信号处理技术领域。该方法分别获取真实火工冲击响应离散加速度信号和每种待进行相似性测度模拟环境下火工冲击响应离散加速度信号，并对每组信号进行降噪处理；然后获取降噪后的信号对应的时频分布并计算该时频分布对应的Rényi熵。比较真实火工冲击响应Rényi熵序列与各火工冲击模拟环境下火工冲击响应Rényi熵的动态时间规整距离，其中距离最小值对应的模拟环境即为与真实火工冲击响应相似度最高的环境。该方法能够量化不同火工冲击响应时频谱的相似性，更全面地比较火工冲击环境模拟实验等效程度，弥补了传统冲击响应谱分析方法的局限性。

Description

一种火工冲击响应时频域分布相似性测度方法

技术领域

本发明涉及火工冲击信号处理技术领域，特别是涉及一种火工冲击响应时频域分布相似性测度方法。

背景技术

火工冲击是由航天器星箭分离、部组件展开等工作过程中的火工品起爆引起的作用于结构上的瞬态冲击响应，具有瞬态、高频、高量级的特点，是航天器在全生命周期内经历的最苛刻的力学环境之一。为了提高航天器的工作性能和可靠性，需要在地面精确模拟航天器火工冲击环境，考核航天产品来对火工品爆炸冲击的耐受性。目前，对于火工品爆炸引起的高频结构响应，国内外普遍采用冲击响应谱分析方法来鉴定冲击响应量级、规定结构对冲击环境的承受能力。但是，由于冲击响应谱仅仅考虑了冲击信号作用在单自由度系统上响应的峰值信息，忽略了冲击信号的时变特征，冲击响应谱与时域信号难以一一对应，也就造成地面试验与实际火工环境不完全等效的问题，在指导地面试验时具有一定的局限性。为了更精准地对应冲击信号中的有效信息，以便更好地开展地面试验，有必要提出一种适合火工冲击响应数据等效性分析方法，更全面的反映冲击响应随时间的变化细节。

发明内容

本发明的目的是为克服现有技术的不足之处，提出一种火工冲击响应时频域分布相似性测度方法。该方法能够量化不同火工冲击响应时频谱的相似性，更全面地比较火工冲击环境模拟实验等效程度，弥补了传统冲击响应谱分析方法的局限性。

本发明提出一种火工冲击响应时频域分布相似性测度方法，其特征在于，包括以下步骤：

1)获取真实火工冲击响应的离散加速度信号；对该加速度信号进行降噪处理，得到降噪后的真实火工冲击响应离散加速度信号x(t)；

2)利用连续小波变换CWT方法将步骤1)得到的降噪后的真实火工冲击响应的离散加速度信号转化为对应的真实火工冲击响应的时频分布；

其中，采用复小波函数ψ(t)作为连续小波变换的基本小波，表达式如下：

式中f_B为带宽，f_C为小波中心频率，t为时间变量，j为虚数；

通过对连续小波变换基本小波的扩展和变换得到连续小波变换基波函数ψ_a，b(t)，表达式如下：

式中，a为展开变量、b为平移变量；

降噪后的真实火工冲击响应离散加速度信号s(t)经过CWT变换后得到对应的时频分布C_x(a，t)：

其中，

表示基波函数ψ_a，b(t)的共轭函数，时间变量t、展开变量a、平移变量b为式中的连续变量；

3)计算C_x(a，t)中每个采样时刻t_i对应时频分布C_x(a，t_i)的Rényi熵H_αx(t_i)，表达式如下：

其中，α为Shannon熵的均值松弛系数；n为真实火工冲击响应离散加速度信号采样点数；

对式(4)计算得到的结果进行排列得到真实火工冲击响应Rényi熵序列：

x＝[H_αx(t₁)，H_αx(t₂)，…，H_αx(t_n)] (5)

4)从待进行相似性测度的火工冲击模拟环境中任意选取一种模拟环境；在该环境下进行火攻冲击环境模拟，对该模拟环境下火工冲击响应离散加速度信号进行采样；然后对采样得到的加速度信号行降噪处理，得到降噪后的模拟环境下火工冲击响应离散加速度信号s(t)；

5)利用连续小波变换CWT方法将步骤4)得到的降噪后的模拟环境下火工冲击响应离散加速度信号转化为对应的模拟环境下火工冲击响应的时频分布；

降噪后的模拟环境下火工冲击响应离散加速度信号s(t)经过CWT变换后得到对应的时频分布C_s(a，t)：

6)计算C_s(a，t)中每个采样时刻t_i对应时频分布C_s(a，t_i)的Rényi熵H_αs(t_i)，表达式如下：

其中，m为该模拟环境下火工冲击响应离散加速度信号采样点数；

对式(9)计算得到的结果进行排列得到该模拟环境下火工冲击响应Rényi熵序列：

y＝[H_αs(t₁)，H_αs(t₂)，…，H_αs(t_m)] (10)

7)计算真实火工冲击响应Rényi熵序列x与步骤4)选取的模拟环境下火工冲击响应Rényi熵序列y的动态时间规整距离d_DTW(x，y)：

d_DTW(x，y)＝γ(x_n，y_m) (11)

其中，γ(x_i，y_j)是x中的第i个元素到y中的第j个元素之间的累积距离，i＝1，2，…，n，j＝1，2，…，m，计算表达式如下：

其中，|x_i-y_j|表示x中的第i个元素到y中的第j个元素之间的欧氏距离；当x_i或y_j的下标为0时，γ(x_i，y_j)按0计算；

8)重复步骤4)到步骤7)，得到真实火工冲击响应Rényi熵序列与每种待进行相似性测度的火工冲击模拟环境下火工冲击响应Rényi熵序列的动态时间规整距离；对所有动态时间规整距离结果按从小到大进行排序，其中距离最小值对应的模拟环境即为与真实火工冲击响应相似度最高的环境。

本发明的特点及有益效果在于：

相对于冲击响应谱分析，时频分析方法能够表征信号能量和频率随时间的变化规律，可用于分析冲击响应数据的时变特性。时频分析是在时频联合域上将一维信号转化为二维时频图像进行分析的方法，是分析时变非平稳信号的有力工具。相比于时域分析与频域分析的方法，时频分析方法能够更加全面地反映火工冲击响应数据信号在时域和频域的分布特性。本发明从能量角度对火工冲击响应进行时频分析，建立基于Rényi熵函数和动态时间规整的火工冲击响应时频域分布相似性测度方法，以用于衡量不同环境火工冲击响应时频域的等效程度，弥补单一使用冲击响应谱进行环境评价的不足。该方法能够量化不同火工冲击响应时频谱的相似性，更全面地比较火工冲击环境模拟实验等效程度，为地面冲击试验条件的制定提供参考。

附图说明

图1是本发明实施例的整体流程图；

图2是本发明实施例的真实火工冲击响应时域信号图；

图3是本发明实施例的真实火工冲击响应数据CWT时频分布图；

图4是本发明实施例的火工冲击响应数据CWT-Rényi熵的DTW相似性测度示意图。

具体实施方式

本发明提出一种火工冲击响应时频域分布相似性测度方法，下面结合附图和具体实施例对本发明进一步详细说明如下

本发明提出一种火工冲击响应时频域分布相似性测度方法，整体流程如图1所示，包括以下步骤：

1)获取真实火工冲击响应的离散加速度信号；对该加速度信号进行初步的降噪处理，得到降噪后的真实火工冲击响应离散加速度信号x(t)，记信号采样点数为n。

真实火工冲击响应的离散加速度信号来源于火箭级间分离时爆炸螺栓爆炸产生的箭体振动，为箭上实测数据，具体测量要求根据火箭总体设计确定。

2)利用连续小波变换(CWT)方法将步骤1)得到的降噪后的真实火工冲击响应的离散加速度信号转化为对应的真实火工冲击响应的时频分布；

其中，本发明采用复小波函数ψ(t)作为连续小波变换的基本小波，表达式如下：

式中f_B为带宽，f_C为小波中心频率。t为时间变量，j为虚数。

通过对连续小波变换基本小波的扩展和变换得到连续小波变换基波函数ψ_a,b(t)，表达式如下：

式中，a为展开变量、b为平移变量。

降噪后的真实火工冲击响应离散加速度信号x(t)经过CWT变换后得到对应的时频分布C_x(a,t)：

其中，

表示基波函数ψ_a,b(t)的共轭函数，时间变量t、展开变量a、平移变量b为式中的连续变量。

利用CWT时频分析方法将一维火工冲击响应离散加速度信号转换为二维时频联合域能量分布函数，全面反映火工冲击响应信号的时频域非平稳时变特性。初始信号x(t)经CWT变换得到时频分布C_x(a,t)。

CWT可以部分定位时间，并且时间窗的大小是变化的，适用于频率波动大的非平稳信号，是一种自适应的时频分析方法，可以进行多分辨率分析。本发明中CWT的实现依赖于复合Morlet小波函数，所得到的火工冲击响应时频分布能够保持较高的时频聚集性。

3)降噪后的真实火工冲击响应离散加速度信号的时频分布C_x(a,t)可以看作是时域和频域的二维能量密度分布，降噪后的真实火工冲击响应离散加速度信号采样点数为n，计算C_x(a,t)中每个采样时刻t_i对应火工冲击响应时频分布C_x(a,t_i)的Rényi熵H_αx(t_i)，表达式如下：

通过Shannon熵的均值松弛系数α>0，得到了参数化的Rényi熵。三阶Rényi熵可以度量时频分布的信号复杂度，即α＝3；

对式(4)计算得到的结果进行排列得到真实火工冲击响应Rényi熵序列：

x＝[H_αx(t₁),H_αx(t₂),L,H_αx(t_n)] (5)

利用Rényi熵函数计算固定时间步长火工冲击响应信号时频分布能量密度信息，由少量基本时频成分组成的时频分布产生的Rényi熵值较小，而较复杂信号的时频分布产生的Rényi熵值较大。

图2是本发明实施例所参照的的真实火工冲击响应时域信号图，本发明实施例的真实火工冲击响应离散加速度信号由振动加速度传感器测量得到。冲击信号来自于单机冲击试验火箭端测点Z方向的数据，通过对加速度信号进行采样获取离散加速度信号，信号的采样率为10000Hz，采样点数为5201。对应的真实火工冲击响应时频分布如图3所示。

4)从待进行相似性测度的火工冲击模拟环境中任意选取一种模拟环境；在该环境下进行火攻冲击环境模拟，对该模拟环境下火工冲击响应离散加速度信号进行采样；然后对采样得到的加速度信号进行初步的降噪处理，得到降噪后的模拟环境下火工冲击响应离散加速度信号。其中，加速度信号的采样频率fs应高于冲击的理论频率，信号的采样时长应包含一整个冲激响应从激励到衰减的全部振动，记信号采样点数为m。

本实施例中采用火攻爆炸式、摆锤冲击、气枪冲击、电动振动或电液振动等任一种方式进行火攻冲击环境模拟，采用振动加速度传感器测量、声发射测量或者多普勒激光测量得到模拟环境下火工冲击响应离散加速度信号，利用小波阈值降噪的方法滤除获取的加速度信号中的高斯白噪声，得到降噪后的模拟环境下火工冲击响应离散加速度信号s(t)；

5)利用连续小波变换(CWT)方法将步骤4)得到的降噪后的模拟环境下火工冲击响应离散加速度信号转化为对应的模拟环境下火工冲击响应的时频分布；其中，本发明采用复小波函数ψ(t)作为连续小波变换的基本小波，表达式如下：

式中f_B为带宽，f_C为小波中心频率。t为时间变量，j为虚数。

通过对连续小波变换基本小波的扩展和变换得到连续小波变换基波函数ψ_a,b(t)，表达式如下：

式中，a为展开变量、b为平移变量。

降噪后的模拟环境下火工冲击响应离散加速度信号s(t)经过CWT变换后得到对应的时频分布C_s(a,t)：

其中

表示基波函数ψ_a,b(t)的共轭函数，时间变量t、展开变量a、平移变量b为式中的连续变量。

6)计算C_s(a,t)中每个采样时刻t_i对应火工冲击响应时频分布C_s(a,t_i)的Rényi熵H_αs(t_i)，表达式如下：

通过Shannon熵的均值松弛系数α＞0，得到了参数化的Rényi熵。三阶Rényi熵可以度量时频分布的信号复杂度，即α＝3

对式(9)计算得到的结果进行排列得到该模拟环境下火工冲击响应Rényi熵序列：

y＝[H_αs(t₁)，H_αs(t₂)，…，H_αs(t_m)] (10)

7)利用动态时间规整方法对非时齐条件下的时频分布相似性进行测量，通过递归计算真实火工冲击响应Rényi熵序列x与模拟环境下火工冲击响应Rényi熵序列y的动态时间规整距离d_DTW(x，y)：

d_DTW(x，y)＝γ(x_n，y_m) (11)

其中，γ(x_i，y_j)(i＝1，2，…，n，j＝1，2，…，m)是x中的第i个元素到y中的第j个元素之间的累积距离，计算表达式如下：

其中，|x_i-y_j|表示x中的第i个元素到y中的第j个元素之间的欧氏距离，当x_i或y_j的下标为0时，γ(x_i，y_j)按0计算；

8)重复步骤4)到步骤7)，得到真实火工冲击响应Rényi熵序列与每种待进行相似性测度的火工冲击模拟环境下火工冲击响应Rényi熵序列的动态时间规整距离；对所有动态时间规整距离结果按从小到大进行排序，其中距离最小值对应的模拟环境即为与真实火工冲击响应相似度最高的环境。

图4是本发明实施例的火工冲击响应数据CWT-Rényi熵与DTW的相似性测度示意图。结合图4，动态时间翘曲首先被提出作为一种比较时间序列的方法。DTW是两个有序数列之间的弹性距离测度。一个迭代过程检查可能的畸变序列，以尽量减少他们之间的总距离。设Rényi熵值向量时间函数为x和y，式中，通过局部优化的方法实现加权距离总和最小，即最优匹配对应于一个序列与另一个序列的最短翘曲路径。利用DTW进行距离测量的优势在于，不同的火工冲击响应难以做到严格时齐，DTW在测量相似度之前，利用动态路径规划的方法将其中一个(或者两个)序列在时间轴进行规整，以达到更好的对齐。相比于其他距离测度，更适用于Rényi熵值向量相似度的测量。

对于任意的模拟火工冲击环境，根据模拟环境火工冲击响应离散加速度信号与真实火工冲击响应离散加速度信号对应时频分布的DTW距离，距离值越小，则说明该模拟环境与真实环境下所得到的冲击响应的时频特征越接近。对于实际的工程应用可作为火工冲击环境模拟程度的一种度量方式。

本发明采用现代谱分析的理论方法，从能量角度对火工冲击响应进行分析，建立以为三阶Rényi熵函数和动态时间规整为基础的火工冲击响应时频域分布相似性测度分析方法。相比于冲击响应谱，CWT时频分析方法能够同时反映冲击响应的能量时频分布特征，能够更清晰地反映冲击响应的触发时间和次数，更直观地表征冲击激发的结构模态响应分布规律。在进行地面冲击试验条件的制定时，采用本发明实施例的火工冲击响应时频域分布相似性测度指标可在保证试验条件制定方便、高效的同时，保证模拟试验方法与真实火工冲击在频域和时域上均具有较好的一致性，能够提高地面试验模拟真实火工冲击环境的精准性，具有较强的工程背景，对航天器重量减轻、缓冲设计和冲击验收试验标准制定有较大的参考意义。

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

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

Claims (1)

1.一种火工冲击响应时频域分布相似性测度方法，其特征在于，包括以下步骤：

1)获取真实火工冲击响应的离散加速度信号；对该加速度信号进行降噪处理，得到降噪后的真实火工冲击响应离散加速度信号x(t)；

2)利用连续小波变换CWT方法将步骤1)得到的降噪后的真实火工冲击响应的离散加速度信号转化为对应的真实火工冲击响应的时频分布；

其中，采用复小波函数ψ(t)作为连续小波变换的基本小波，表达式如下：

式中f_B为带宽，f_C为小波中心频率，t为时间变量，j为虚数；

通过对连续小波变换基本小波的扩展和变换得到连续小波变换基波函数ψ_a,b(t)，表达式如下：

式中，a为展开变量、b为平移变量；

降噪后的真实火工冲击响应离散加速度信号s(t)经过CWT变换后得到对应的时频分布C_x(a,t)：

其中，

表示基波函数ψ_a,b(t)的共轭函数，时间变量t、展开变量a、平移变量b为式中的连续变量；

3)计算C_x(a,t)中每个采样时刻t_i对应时频分布C_x(a,t_i)的Rényi熵H_αx(t_i)，表达式如下：

其中，α为Shannon熵的均值松弛系数；n为真实火工冲击响应离散加速度信号采样点数；

对式(4)计算得到的结果进行排列得到真实火工冲击响应Rényi熵序列：

x＝[H_ax(t₁)，H_αx(t₂)，…，H_αx(t_n)] (5)

4)从待进行相似性测度的火工冲击模拟环境中任意选取一种模拟环境；在该环境下进行火攻冲击环境模拟，对该模拟环境下火工冲击响应离散加速度信号进行采样；然后对采样得到的加速度信号行降噪处理，得到降噪后的模拟环境下火工冲击响应离散加速度信号s(t)；

5)利用连续小波变换CWT方法将步骤4)得到的降噪后的模拟环境下火工冲击响应离散加速度信号转化为对应的模拟环境下火工冲击响应的时频分布；

降噪后的模拟环境下火工冲击响应离散加速度信号s(t)经过CWT变换后得到对应的时频分布C_s(a,t)：

6)计算C_s(a,t)中每个采样时刻t_i对应时频分布C_s(a,t_i)的Rényi熵H_αs(t_i)，表达式如下：

其中，m为该模拟环境下火工冲击响应离散加速度信号采样点数；

对式(9)计算得到的结果进行排列得到该模拟环境下火工冲击响应Rényi熵序列：

y＝[H_αs(t₁)，H_αs(t₂)，…,H_αs(t_m)] (10)

7)计算真实火工冲击响应Rényi熵序列x与步骤4)选取的模拟环境下火工冲击响应Rényi熵序列y的动态时间规整距离d_DTW(x,y)：

d_DTW(x,y)＝γ(x_n,y_m) (11)

其中，γ(x_i,y_j)是x中的第i个元素到y中的第j个元素之间的累积距离，i＝1,2,…,n,j＝1,2,…,m，计算表达式如下：

其中，|x_i-y_j|表示x中的第i个元素到y中的第j个元素之间的欧氏距离；当x_i或y_j的下标为0时，γ(x_i,y_j)按0计算；

8)重复步骤4)到步骤7)，得到真实火工冲击响应Rényi熵序列与每种待进行相似性测度的火工冲击模拟环境下火工冲击响应Rényi熵序列的动态时间规整距离；对所有动态时间规整距离结果按从小到大进行排序，其中距离最小值对应的模拟环境即为与真实火工冲击响应相似度最高的环境。

Images