基于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值選取提供參考。