基于STFT-AOK 的星箭力学环境等效频谱分析

邹元杰，朱卫红，陈小旺，冯志鹏，刘绍奎，庞世伟

（1.北京空间飞行器总体设计部，北京 100094； 2.北京科技大学 机械工程学院，北京 100083）

0 引言

卫星和运载火箭在实际发射过程（如点火、关机、级间分离、星箭分离）以及航天器在着陆过程中要经历复杂的瞬态力学环境，研究人员从实测数据中获取的通常是航天器力学环境参数的时域信息。然而，航天器研制过程中的地面验证试验通常是在频域内开展的，如正弦振动、随机振动和混响噪声等试验都采用频域试验方法。为了制定合理的频域试验条件，必须得到合理的力学环境频谱幅值。因此，需要将用于试验控制的物理量的时域信息转换为相应的频域量（即等效频谱），再对频谱进行一定的包络处理。

在将低频（一般在100 Hz 以内）瞬态响应由时域转换至频域时，传统的傅里叶变换无法得到合理的幅值特性。因此，工程上往往利用冲击响应谱（shock response spectrum, SRS）变换方法得到正弦振动等效频谱[1-3]，即先计算加速度的冲击响应谱，再除以稳态放大系数Q。冲击响应谱是目前使用最为广泛的瞬态冲击信号的频域表征方法[4-5]。该方法在理论上对于任意平稳或非平稳信号均适用，而不受傅里叶变换的信号平稳性假设约束[6]，计算也比较简便。但有研究表明，采用工程方法获取的等效频谱随着Q取值的增大，幅值明显减小[3]。可见，传统的等效频谱分析工程方法受Q取值的影响较大，尚有不合理之处。

基于上述原因，本文尝试基于信号的时频分析构造合理的等效频谱表达。传统时频分析方法难以有效分析同时包含周期成分和瞬态成分的复杂非平稳信号[7]：以短时傅里叶变换（short-time Fourier transform, STFT）、连续小波变换为典型代表的线性时频分析方法受到Heisenberg 不确定性制约，无法同时具备高时间分辨率和高频率分辨率。双线性时频分析方法（如Wigner-Ville 分布）具备理想的时频分辨率，但是当待分析信号具有多个频率成分时，将不可避免地产生交叉项及自项干扰，妨碍频率结构的准确识别，且时频分布的幅值不与真实幅值成比例变化。如何同时满足高时频分辨率要求、避免干扰项且追求幅值与真实值一致，是时频分析研究领域的难题[8-10]。

针对该难题，本文提出一种基于短时傅里叶变换-自适应优化核函数（STFT-AOK）融合技术的等效频谱分析方法，拟充分利用线性时频分布STFT 幅值分析准确以及自适应优化核函数（adaptive optimal kernel, AOK）方法时频分辨率高的优点，构造较为理想的时频分布，从而获取瞬态信号合理的频谱特性。

1 基于冲击响应谱的等效频谱分析工程方法

1.1 等效频谱分析的工程方法

冲击响应谱是根据一定形式的外激励（通常为加速度基础激励）给出最大响应与固有频率的关系曲线[3]。对于单自由度系统，加速度基础激励对应的冲击响应谱为

式中：y(t)和x(t)分别为弹簧振子和基础的绝对位移；tm为最大幅值对应的时刻；ωn为系统的固有圆频率。

为了确定星箭力学环境等效频谱特性，工程上往往借助冲击响应谱的概念换算出基础加速度激励的等效频谱。比如，在星箭耦合分析中，首先利用瞬态响应分析方法得到星箭结构上任意节点的加速度时域响应；然后以该响应作为若干不同固有频率的虚拟单自由度弹簧振子的基础激励，求解这些振子的最大响应（即冲击响应谱）；最后将该响应谱除以稳态放大系数Q（Q=1/(2ζ)）作为该节点的加速度频谱，用于制定卫星或部组件的振动试验条件。即，工程上定义ẍ(t)的等效频谱为

1.2 工程方法的特点和局限性

冲击响应谱所假想的一系列线性单自由度弹簧阻尼系统可等效为一系列卷积滤波器[4-5]。然而，该方法所构造的不同固有频率的单自由度系统的幅频特性和相频特性[5]存在差异，导致冲击响应谱的不同频段理论上不具备统一量度，从而带来一系列问题。

比如，给定一个由10 个正弦波叠加而成的时域加速度信号，正弦波幅值和相位随机生成，即

其 中：[A1,A2,···,A10]=[66.1, 93.2, 34.5, 42.2, 54.3,68.7, 29.1, 34.0, 27.1, 77.7] m/s2；[φ1,φ2,···,φ10]=[5.2, 3.2, 2.3, 1.4, 3.4, 1.8, 0.4, 0.5, 0.4, 2.6] rad；fi=10×iHz。按照式(1)～式(3)，该信号不同Q取值下所对应的等效频谱如图1 所示。

图1 某时域加速度信号不同Q 取值下的等效频谱Fig.1 Equivalent frequency spectra of a given time-domain acceleration signal calculated with different Q values

从图1 可看出，采用工程方法获取等效频谱，受Q取值影响非常大：随着Q取值的增大，频谱幅值明显减小，甚至当Q趋于无穷大时幅值会无限趋于0（详见文献[3]）；与真实值相比，工程方法获得的等效频谱幅值或偏大或偏小。因此，利用工程方法获得等效频谱，进而制定力学环境试验条件，同时存在“过试验”和“欠试验”2 种风险；并且在星箭耦合分析等效频谱变换时，Q值的选取（目前工程上通常取20）与星箭耦合动力学分析中阻尼比的选择并无直接关系，主要依据工程经验。

2 基于STFT-AOK 的等效频谱分析方法

为了更加准确地分析瞬态载荷信号的频率结构、判断不同频率特征的强弱，本文从信号特征提取的角度开展研究。鉴于传统的频域分析方法不适合于分析瞬态非平稳信号，本文重点研究基于时频分析的载荷信号分析及频谱构造方法。

2.1 基于短时傅里叶变换的频谱分析

短时傅里叶变换（STFT）是最为常用的时频分析方法之一。其核心是将信号表示为一系列加权后的三角函数基函数的总和[7]。对于任意有限能量信号s(t)，其短时傅里叶变换可表示为

式中：w(·)为窗函数；τ为时间变量。通过设计合适的窗函数形式，可将时变的信号特征在时间-频率-幅值三维空间中表现出来。从时频分布中，可以准确识别变化的瞬时频率曲线以及对应的时变幅值。将时频矩阵每一行的最大幅值绘制在频域图中，可得到原始信号的频谱，

基于STFT 的时频分析和频谱构造具有计算量小、幅值较为准确的优点；但该方法的窗函数固定且受到Heisenberg 不确定性制约，不能很好地匹配复杂多分量信号的不同时变特征。当基函数尺寸与目标特征不匹配时，将造成时频模糊现象，并影响构造频谱的准确性。

2.2 基于自适应优化核函数方法的时频分析

自适应优化核函数（AOK）方法利用最优核函数抑制模糊域中的干扰项，可实现较理想的时频分辨率[8]。核函数为径向高斯函数

式中σ(θ)为高斯函数沿径向角度θ=arctan (τ/υ)的方差。在极坐标系中，可通过求解优化问题

获得最优核函数，其中：A(r,θ)为极坐标系中的模糊函数；r为高斯函数的模；c为高斯函数的体积。所求得的最优核函数类似低通滤波器，当信号的自项成分聚集在模糊平面坐标轴和坐标原点周围，且干扰项分布在远离坐标轴区域时，干扰项能够被滤除，从而获得分辨率较高且干扰项被抑制的时频分布AOKs(t,f)。

2.3 基于STFT-AOK 的等效频谱分析

由2.1 节和2.2 节可知，STFT 以基函数变换为核心，通过设计合适的基函数，可以在构造的时频分布中准确揭示特征频率的幅值大小，有利于准确构造信号频谱；但固定的时频分辨率使其产生时频模糊，不适用于分析包含瞬态冲击特征的信号。AOK方法能够根据不同信号特征变换核函数，从而匹配不同频率特征以及瞬态冲击，可较为准确地还原信号特征；但由于瞬时自相关函数的计算和核函数滤波，实际信号的幅值未被还原，无法直接用于信号构造频谱。为了同时实现准确的幅值识别和高分辨率特征提取，本文提出一种融合技术，将STFT 与AOK 方法所得时频分布相结合，同时保留STFT 的准确幅值和AOK 方法的准确时频结构。STFT-AOK方法的流程如图2 所示，具体包括：

图2 STFT-AOK 方法流程示意Fig.2 Flowchart of STFT-AOK method

1）计算待分析信号的AOK 时频分布；

2）通过设定阈值，从时频分布中定位幅值较大的主导时频区域Q1；

3）采用适合的窗函数（窗宽通常取1/4 信号长度，对于冲击信号取1/10 信号长度），计算待分析信号的STFT 时频分布；

4）将AOK 时频分布中主导时频区域Q1 对应的幅值修正为STFT 时频分布中对应值；

5）从AOK 时频分布中定位冲击时频区域Q2；

6）采用窗宽较短的窗函数，计算待分析信号的STFT 时频分布；

7）将AOK 时频分布中冲击时频区域Q2 对应的幅值修正为STFT 时频分布中对应值；

8）得到STFT 与AOK 融合后的时频分布；

9）取所得融合时频分布每行的最大值，构成对应频率处的幅值，构造信号的频谱。

3 已知仿真信号的等效频谱分析

开展已知仿真信号的数值分析，以说明STFTAOK 方法相比于传统时频分析方法和SRS 方法的瞬态非平稳特征表达优势。以一个包含时变频率成分和瞬态冲击的非平稳加速度信号s(t)（单位为g）为例，其表达式为

该仿真信号的时域波形和Fourier 频谱分别如图3(a)和图3(b)所示。从图3(b)中仅能找到100 Hz的谱峰，对应瞬态冲击振荡频率；此外，该谱峰对应的幅值为1.28g，远小于实际仿真信号在该频率处设置的幅值50g。

图3 仿真信号的时域波形和Fourier 频谱Fig.3 Time-domain waveform and Fourier frequency spectrum of the simulated signal

取Q=10，仿真信号的SRS 如图4 所示。对比图4 和图3(b)可知，SRS 对于线性变化的时变频率成分幅值提取较为准确；但对于瞬态冲击成分而言，SRS 的幅值与真实幅值相比差异较大。

图4 仿真信号的SRS 及真实幅值Fig.4 SRS and true amplitudes of the simulated signal

以下分别采用STFT、连续小波变换（continuous wavelet transform, CWT）、Wigner-Ville 分布（Wigner-Ville distribution, WVD)、AOK 以及STFT-AOK 方法，构造仿真信号（信号长度为1000 ms）的时频分布及频谱。

仿真信号的STFT 分析结果如图5 所示。采用高斯窗且设窗宽（wl）为60 ms 所得到的时频分布如图5(a)所示，从图中可准确识别时变频率成分和瞬态冲击成分。通过提取时频矩阵中每一行对应的最大幅值构造的频谱如图5(b)所示。对比构造的频谱和真实值可以看出：在线性变化频率范围内，STFT 能够较为准确地提取频率成分对应的幅值；但对于瞬态冲击，所得到的结果仍小于真实幅值。将STFT 中的窗宽设为10 ms 所得到的时频分布及频谱如图5(c)和图5(d)所示。对比可看出：当取较窄的窗宽时，冲击成分的幅值能近似地被识别，但所得线性变化频率成分的幅值不准确；不同窗函数取值下，STFT 可以分别提取变化较慢的频率成分特征和瞬态冲击特征。但STFT 应用时的窗宽是固定的，无法同时实现2 类特征的幅值提取；此外，在Heisenberg 不确定原理的限定下，STFT 的时频分辨率低，使得构造的频谱呈现宽带分布。

图5 仿真信号的STFT 分析Fig.5 STFT analysis of the simulated signal

CWT 也是一种常用的时频分析方法，相比于STFT 而言，具有变化的时频分辨率。仿真信号的Morlet 小波尺度谱如图6(a)所示，可以看出，虽然冲击成分被较好地体现出来，但由于频率分辨率降低，线性变化的时变频率成分无法得到准确揭示。此外，基于CWT 构造的频谱的幅值不与实际幅值成固定比例关系，CWT 时频分布的幅值除以比例系数ratio=1×103后，所得频谱无法准确反映实际信号频率成分的强弱，且由于线性变化频率成分在时频分布中能量泄漏严重，该特征在频谱中几乎得不到体现（见图6(b)）。

图6 仿真信号的CWT 分析Fig.6 CWT analysis of the simulated signal

WVD 是一种典型的二次型时频分布，该方法的特点是无须适用基函数内积变换，而是利用信号自身的瞬时自相关函数进行频域转换。仿真信号的WVD 时频分布如图7(a)所示，可以看出，WVD具有高时频分辨率，但存在除瞬态冲击和时变频率成分外的虚假频率成分，即交叉项干扰。基于该时频分布构建的频谱如图7(b)所示，可以看出：因为干扰项的存在，产生了额外的谱峰；且因为瞬时自相关函数的二次型计算，频谱幅值与真实幅值不成固定比例关系。

图7 仿真信号的WVD 分析Fig.7 WVD analysis of the simulated signal

基于AOK 的时频分布如图8(a)所示，可以看出，AOK 方法构造的时频分布相比STFT 而言，具有显著的高时频分辨率优势；但所构造的频谱（如图8(b)所示）的幅值不与真实幅值成固定比例关系，AOK 时频分布的幅值除以比例系数ratio=2×105后，所得频谱无法准确反映实际信号频率成分的强弱。

图8 仿真信号的AOK 分析Fig.8 AOK analysis of the simulated signal

采用STFT-AOK 方法构造的时频分布如图9(a)所示，可见，STFT-AOK 方法不但能够利用高分辨率AOK 时频分布作为时频区域划分依据，而且综合了不同窗函数下STFT 时频分布的准确幅值提取能力。STFT-AOK 方法所构造的频谱如图9(b)所示，可见采用该方法能够同时准确提取时变频率成分和瞬态冲击成分的频率及幅值信息。

图9 仿真信号的STFT-AOK 分析Fig.9 STFT-AOK analysis of the simulated signal

4 星箭力学环境数据的等效频谱分析

针对运载火箭助推器分离后状态和一、二级分离后状态，由星箭耦合分析得到典型星箭界面加速度时间历程ẍ1(t)（见图10(a)）和ẍ2(t)（见图10(b)）；分别采用本文提出的STFT-AOK 方法和传统的基于冲击响应谱的工程方法分析这2 种状态的加速度数据，得到的等效频谱分别见图11(a)和图11(b)，其中，基于冲击响应谱的工程方法中Q值分别取为10、20 和30。

图11 STFT-AOK 和SRS 方法得到的ẍ1(t)和ẍ2(t)的等效频谱对比Fig.11 Comparison of equivalent spectrum for ẍ1(t) and ẍ2(t)by STFT-AOK and SRS methods respectively

比较图11 中2 种方法的结果可以看出，用传统的工程方法得到的等效频谱幅值与STFT-AOK方法得到的有明显区别：

1）当Q=10 时，用传统工程方法得到ẍ1(t)的频谱峰值为1.59 m/s2，比STFT-AOK 方法得到的1.44 m/s2高10%；而用传统工程方法得到ẍ2(t)的频谱峰值为5.16 m/s2，比STFT-AOK 方法得到的6.86 m/s2低25%。

2）当Q=20 时，用传统工程方法得到ẍ1(t)的频谱峰值为1.28 m/s2，ẍ2(t)的频谱峰值为3.44 m/s2，分别比STFT-AOK 方法得到的低10%和50%。如果只关注2 个算例频谱曲线中最大的2 个峰值，那么STFT-AOK 方法得到的等效频谱幅值大概位于传统工程方法下Q=10 和Q=20 的结果之间。

目前我国星箭力学环境领域采用工程方法计算加速度等效频谱时通常取Q=20，而从上面的2 个算例分析结果看，传统工程方法得到的频谱峰值相对于STFT-AOK 方法的偏差可达50%以上，这可能对力学环境条件制定产生一定影响。

5 结束语

本文首先针对基于冲击响应谱的等效频谱分析工程方法的局限性，提出了基于STFT-AOK 融合技术的星箭力学环境等效频谱分析方法；而后，利用含瞬态冲击的复杂非平稳仿真信号开展等效频谱分析，对多种方法的分析结果进行对比研究。结果表明：利用STFT-AOK 方法，在基于AOK 时频分布准确定位频率成分和瞬态特征的时频区域基础上，结合STFT 时频分布的准确幅值，构造分辨率高且幅值准确的时频分布，可实现力学环境等效频谱的准确获取，适用于复杂的瞬态星箭力学环境。最后针对2 组典型星箭界面加速度信号，采用STFTAOK 方法和传统工程方法获得信号的等效频谱并进行对比，结果显示：用传统工程方法得到的等效频谱幅值与STFT-AOK 方法得到的有明显区别。

STFT-AOK 方法将参数视为瞬态信号，不受冲击响应谱对应的基础加速度激励响应分析理论限制，可用于各类星箭力学环境参数（如力、力矩、应力、应变、位移、速度和角速度等）分析。后续需要进一步简化STFT-AOK 方法的流程，形成便于操作的标准规范；对更多星箭耦合分析数据和飞行遥测数据开展分析，充分对比STFT-AOK 方法和传统工程方法获得的等效频谱，为传统工程方法的Q值选取提供参考。