基于最大譜的非平穩隨機振動數據分析處理方法研究

徐俊，傅耘，張建軍，孫建勇，王星皓

（中國航空綜合技術研究所，北京 100028）

基于最大譜的非平穩隨機振動數據分析處理方法研究

徐俊，傅耘，張建軍，孫建勇，王星皓

（中國航空綜合技術研究所，北京 100028）

目的 解決非平穩隨機振動環境實測數據無法采用功率譜密度分析方法的問題。方法 基于最小均方誤差準則，歸納總結FFT和數字濾波兩種最大譜分析方法。通過對比兩種方法在不同分辨率帶寬時計算的優缺點，并分析其誤差。對優選出來的最大譜分析方法進行舉例驗證。結果 理論分析和仿真結果均表明基于比例分辨率帶寬的FFT最大譜分析方法的誤差最小。將FFT最大譜分析方法應用于某型炸彈9個架次自由飛振動數據的處理和環境條件確定中，結果合理有效。結論 可用比例分辨率帶寬的FFT法分析非平穩隨機振動環境實測數據，為基于最大譜制定型號振動環境試驗條件提供指導。

非平穩隨機信號；最大譜； FFT；數字濾波

近年來，隨著裝備作戰需求的提高，各個航空航天型號所面臨的環境愈加復雜多樣，其環境適應性要求也越來越高，而GJB 150A/899A是基于美國裝備的特點制訂的標準，不完全適用于國內裝備產品，因而基于實測數據確定環境條件也顯得愈加重要。對于運載器、火箭、炸彈等裝備在飛行過程中，經歷熱、振和噪聲等綜合環境尤其復雜，并且采集到的振動樣本少、飛行時間短、狀態變化大，主要表現為非平穩隨機振動數據。因此，有必要選擇合理的數據分析方法對非平穩隨機振動數據進行處理，為振動環境條件確定提供基礎[1—3]。

目前，功率譜密度分析方法是處理平穩隨機振動數據的重要工具，但不適用于非平穩隨機振動環境實測數據的處理。對于非平穩隨機信號，大多通過時頻分析識別不同時刻的不同頻率成分。裝備環境適應性的驗證，主要是利用特定量值的等效平穩振動代替難以復現的非平穩隨機振動進行環境驗證試驗。因此需要對現有的時頻分析算法進行改進，根據相關設計規范和試驗標準的應用指南，推薦采用基于短時傅里葉變換（STFT）的最大譜作為非平穩隨機振動數據的主要分析處理方法。國內外關于最大譜的工程應用研究較少，難以指導設計和試驗人員合理分析非平穩隨機信號，并得出正確的環境條件。

文中首先結合GJB/Z 222—2005《動力學環境數據采集和分析指南》附錄B[4]中的非平穩隨機振動相關資料，歸納出兩種最大譜分析方法。其次詳細介紹了這兩種方法的分析流程，給出了最大譜分析方法的工程應用建議。最后針對某型號炸彈實測的非平穩隨機振動數據，選擇合理的分析方法計算最大譜，并確定其環境條件。

1　最大譜分析方法

1.1 分析流程

非平穩振動的最大譜通常可解釋為設計和試驗所需要的等效平穩載荷自譜，為非平穩事件中每個頻率分辨率帶寬下最大的譜值，因此采用最大譜的均方根值應該大于原始信號的均方根值。主要分析流程如圖1所示，首先剔除非平穩隨機信號中的瞬態事件；然后確定在最大譜分析過程中的參數，即平均時間T和分析帶寬B；最后根據選用的最大譜計算方法計算自譜，并且通過平均處理提高統計精度得到最大譜結果。

圖1 非平穩隨機信號最大譜分析流程Fig.1 The analysis process of maximum spectral fornonstationary random data

對于隨機信號x(t)的功率譜密度G（f）是x(t)通過中心頻率為f、分析帶寬為B的窄帶濾波器，當B趨于0和平均時間趨于無限時的均方值[5]，數學上可表示為：

式中：T為平均時間；B為分析帶寬；x2(f，t，B)為經窄帶濾波器后信號的瞬時平方值。

在Bendat和Piersol[6]的論文中已經表明，公式(1)在T→∞并且B→0，BT→∞的極限條件下得到了隨機振動的真實功率譜密度。在現實中這些極限條件無法實現，即對于一個實際的PSD估計，BT必須是有限的，并且帶寬B為非0，因此實際的PSD估計存在隨機（統計抽樣）誤差、時間分辨偏置誤差和頻率分辨偏置（頻域平滑）誤差[7—9]。為了將功率譜密度估計的總誤差最小化，需要確定最大譜分析的最佳參數（平均時間 T和分析帶寬B）。

1.2 參數確定

在計算最大譜的過程中，選取合理的參數能夠減小功率譜密度估計的總誤差，主要包括平均時間和頻率分辨率帶寬的確定。

1.2.1 平均時間的確定

在整個非平穩事件上，在所要求的頻率分析范圍內每一個分析帶寬上，通過滑動平均計算時變自譜[10]。對所測的非平穩隨機振動信號，推薦用下列方法選擇譜分析用的固定平均時間[11]：在用解析法或試湊法計算振動信號總均方根值時，對分段線性平均，近似取最佳平均時間T0；當計算整個頻帶的譜時，平均時間取2T0。

1.2.2 頻率分辨率帶寬的確定

1) 比例分辨率帶寬。為了使典型運載器發射時振動數據譜分析的均方誤差接近最小，相應的比例帶寬是1/12倍頻程帶寬（B0=0.058f），而當小于43 Hz時，修正為常帶寬（B0=2.5 Hz）[12]。

2) 固定分辨率帶寬。這種方法實質為常分辨率帶寬，在頻率范圍從20～2000 Hz的譜分析，要達到可接受的均方誤差，可用表1中三種常帶寬來實現[13]。

中心頻率/Hz 　　分析帶寬/Hz 20≤f＜100 　3 100≤f＜400 　10 400≤f≤2000 　30

表1 固定分辨率帶寬

Table 1 Fixed resolution bandwidth

式中：i=0，1，2，…，(n-1)；n為總段數；M為段長（每段之間的數據值數目）。

以上兩種分辨率帶寬能夠在一定程度上減小由于分析帶寬無法趨近于 0而帶來的頻率分辨偏置誤差。對于統計隨機誤差，則采樣分段線性平均的方法，對離散時間序列測量數據x(nΔt)，n=1，2，…，滑動平均的等效值可在N個數據點上(T=N·Δt)用分段的無加權（線性的）平均來計算[14—15]：

1.3 最大譜計算方法

基于最小均方誤差準則，通過功率譜密度的兩種表示形式，歸納出兩種最大譜分析方法，分別是FFT法和數字濾波法。

1.3.1 FFT最大譜分析法

非平穩隨機振動信號的時變自譜可以用 FFT法來估計：

式中：E[|X(f,T)|2]為傅里葉變換后平方值的期望。式(5)的計算框圖如圖2所示。

1.3.1.1 集合平均FFT法

圖2 時變自譜密度函數估計程序(FFT法)Fig.2 Procedure for estimating time-varying autospectral density function (FFT)

計算振動數據的自譜最常用的方法是FFT法。其主要思路是將時域信號分段，然后分別求各個子樣本的自譜，最后進行平均處理提高精度，如圖 3所示。計算步驟為：把感興趣的時間段分隔成 nd個不相連的數據段，每段長為T；計算每段的原始自譜；對原始譜進行集合平均，得到時間為T0=ndT的平均自譜。當采用滑動平均時，步長為 0.2ndTs或小于0.2ndN。

圖3 集合平均FFT法Fig.3 The drawing of conventional FFT computation

1.3.1.2 單次FFT法

單次FFT法進行譜分析如圖4所示，其主要思路是將整體樣本進行單次FFT，然后在頻域上滑動平均求自譜。計算步驟為：將原始信號進行單次FFT變換；在平均時間為T上，采用比例分辨率帶寬計算原始自譜，得到一系列頻率間隔為 Δf=1/T的離散譜分量（該計算應用矩形窗來完成，不應加為抑制旁瓣泄漏的加窗處理）；按要求把相鄰的譜分量一起平均，以得到頻率分辨率帶寬下譜分量的集合平均，它近似于 1/12倍頻程；最后從兩個或更多個的搭接段來計算頻率平均自譜，每段長為N·Δt=T

1.3.2 數字濾波最大譜分析法

非平穩隨機振動信號的時變自譜可以用濾波法來估計：

式中：x2(t,B,f)表示經過中心頻率為f、帶寬為B的窄帶濾波器后的信號瞬時平方值。

圖4 單次FFT法Fig.4 The drawing of single FFT computation

式(4)的計算框圖如圖5所示。

圖5 時變自譜密度函數估計程序(數字濾波法)Fig.5 Procedure for estimating time-varying autospectral density function (digital filter)

1.3.2.1 時域濾波法

使用數字濾波法進行最大譜分析，當采用時域濾波時，原始信號通過窄帶濾波器后，求其瞬時平方值，在時間T上平均，如圖6所示。計算步驟為：x(t) 通過中心頻率為f、帶寬為B的時域帶通濾波器，得到x(f, B, t)；對時域濾波后（采用1/12倍頻程濾波器）的信號進行平方，在時間T上平均，得到自譜φx(f, B, t)；將自譜除以帶寬B得到Gxx(f,t)。

1.3.2.2 頻域濾波法

圖6 時域濾波法Fig.6 The drawing of time-domain filter

使用頻域數字濾波法時，與時域濾波不同的是，首先將原始信號進行快速傅里葉變換，然后采取頻域濾波，得到各個中心頻率上的時域信號，如圖7所示。計算步驟為：將原始信號進行單次FFT變換，在頻域上進行濾波，采用固定分辨率帶寬；然后利用傅立葉逆變換得到時域信號，在平均時長T0上平均；除以帶寬。

圖7 頻域濾波法Fig.7 The drawing of frequency-domain filter

與 FFT最大譜分析法相比，數字濾波方法由于濾波器的頻率響應函數不可能是理想的矩形，在分析的過程中會額外引入靜態誤差，而 FFT法則不存在濾波誤差。因此在平均時間和分析帶寬等參數相同的條件下，選用 FFT最大譜分析法的誤差會相對較小。另外，集合平均 FFT法對總樣本的時間進行分段，計算每個頻率分辨率帶寬下最大的譜值。單次 FFT法則在計算過程中丟失了很多時域上的信息，只是對分析帶寬內的多個譜值進行平均，因此單次 FFT法計算的最大譜值偏小，而集合平均FFT法的最大譜值更為準確。

1.4 最大譜分析誤差

在非平穩隨機振動數據譜分析中，分辨率帶寬是可選的。在分析中，存在時間分辨率偏置誤差和隨機誤差，這些誤差可用下列歸一化形式來估算。

時間分辨率偏置誤差：

隨機誤差：

式中：b[?xx(f,t)]為 ?xx(f,t)估計的偏置誤差，b[?xx(f,t)]=E[?xx(f,t)]-Gxx(f,t)；σ[?xx(f,t)]為 ?xx(f,t)估計的標準差； B為譜估計的頻率分辨率帶寬；T為譜估計的時間分辨率（平均時間）。

為了評估最大譜分析過程中的時間分辨偏置誤差，選擇式(7)中的 d2[Gxx(f,t)]/dt2的值是非常必要的，而二次求導的值與實際樣本分析的自譜值密切相關。在分析樣本相同、平均時間相同的情況下，各方法計算的時間分辨偏置誤差在理論上基本相同。

分析式(6)可知，最大譜分析過程中的隨機誤差與平均時間和分析帶寬相關，通常各分析方法選取的平均時間是相同的，因此無論是 FFT還是數字濾波最大譜分析方法，隨機誤差主要與分析帶寬有關。對比上述幾種最大譜分析方法的隨機誤差，其關系見表2。以中心頻率為橫坐標，隨機誤差為縱坐標，繪制比例分辨率帶寬和固定分辨率帶寬下的隨機誤差如圖所示。

表2 FFT和數字濾波最大譜分析方法的隨機誤差對比Table 2 The comparison of random error for FFT and digital filter maximum spectral

圖8 不同分析帶寬下的隨機誤差Fig.8 The random error of different bandwidths

通過對比分析，當各分析方法的平均時間相同時，在大部分頻率范圍內，比例分辨率帶寬最大譜分析的隨機誤差較小。

綜上所述，從理論分析的角度判斷可知，選用比例分辨率帶寬的 FFT最大譜分析方法時，分析過程中產生的隨機誤差較小，且其誤差隨中心頻率的增大而不斷減小。

2　仿真實例

根據某型號炸彈在自由飛過程中采集的振動數據，結合其物理飛行過程，包括點火、舵機啟控等瞬態事件，設計如下非平穩隨機仿真數據，其振動曲線如圖所示。由于其中包括瞬態事件，因此截取30～49 s之間的數據進行最大譜分析。

圖9 某型炸彈實測非平穩隨機信號Fig.9 Non-stationary random data of a certain missile

經過計算，采用上述的最大譜分析方法進行處理，得到各最大譜分析方法計算結果見圖和表3。

圖10 各方法計算得到的最大譜圖Fig.10 The maximum spectrals

表3 各最大譜分析方法計算結果Table 3 The results of maximum spectrals computation

在非平穩事件期間，每個頻率上譜密度值達到最大值時，它的偏置誤差一般將是最大的，因此在最大譜的每個頻率上計算得到的誤差代表著最大誤差。各種方法的時間偏置誤差和隨機誤差如圖和圖12所示。

圖11 時間偏置誤差Fig.11 Time resolution bias error

圖12 隨機誤差Fig.12 Random error

對比分析結果可知：

1) 集合平均FFT法、單次FFT法和數字濾波法（比例分辨率帶寬）進行最大譜分析得到的均方根值較為接近，且隨機誤差均隨著中心頻率的增大而減小。在100～2000 Hz范圍內時間偏置誤差小于5%；在中心頻率較小時，時間偏置誤差和隨機誤差均較大。

2) 在所有的方法中，單次FFT法的計算結果與原始信號的均方根值最接近，因為單次 FFT法進行了頻率平均，并不是每個頻率帶寬下的最大的譜值，因此均方根值與原始信號最接近。

3) 數字濾波法中的頻域濾波計算的結果偏差最大，這是由于原始信號進行濾波時引入了靜態誤差。

4) 無論是FFT法還是數字濾波法，選用固定分辨率帶寬計算出的最大譜均方根偏大，因為在平均時間相同時，固定分辨率帶寬引入了更大的隨機誤差。

5) 通過對比時間偏置誤差和隨機誤差的大小，隨機誤差普遍較大，是仿真實例最大譜計算總誤差的主要組成部分。

綜上所述，通過仿真分析，各種最大譜分析方法的誤差與理論分析一致，推薦在工程應用中選用FFT比例分辨率帶寬的最大譜分析方法。

3　應用舉例

針對某型炸彈自由飛過程，共實測9個架次的數據，選用基于 FFT比例分辨率帶寬的最大譜分析方法進行計算，并確定其環境試驗條件。選取該型炸彈所有實測架次的后艙段y向實測數據，分別計算其最大譜并取包絡，得到最終的試驗規范譜如圖所示。

圖13 后艙段y向規范譜Fig.13 Normalized spectral of rear section in y

考慮到外場實測與實驗室模擬試驗的差別，以及實測型號的樣本數量有限，分別在實測規范譜譜型上乘以載荷分散系數1.2，安全系數1.2，得到后艙段y向試驗條件見表4。

表4 后艙段y向試驗條件Table 4 The value of normalized spectral of rear section in y

4　結論

通過對比 FFT和數字濾波兩種最大譜分析方法的優缺點并分析其誤差，可得以下結論。

1) 與數字濾波最大譜分析方法相比，選用 FFT（比例分辨率帶寬）最大譜分析方法的誤差最小。

2) 對比兩種方法在不同分辨率帶寬時計算的誤差，理論分析和仿真結果一致，驗證了最大譜分析方法的有效性。

3) 將 FFT最大譜分析方法應用于某型炸彈 9個架次自由飛振動數據的處理和環境條件確定中，獲得的自由飛振動試驗條件合理有效。

[1] 郭迅, 郭強嶺. 空空導彈振動試驗條件分析[J]. 裝備環境工程, 2012, 9(3): 99—103. GUO Xun, GUO Qiang-ling. Analysis of Vibration Test Condition of Air-to-Air Missile[J]. Equipment Environmental Engineering, 2012, 9(3): 99—103.

[2] 郭強嶺, 李立民. 空空導彈掛飛振動試驗條件探討[J].航空兵器, 2003(6): 21—23. GUO Qiang-ling, LI Li-min. Discussion of Vibration Test Condition for Air-to-Air Missile[J]. Aero Weaponry, 2003(6): 21—23.

[3] 郭強嶺. 空空導彈振動試驗條件確定方法[J]. 航空標準化與質量, 2006(3): 34—37. GUO Qiang-ling. Determination of Vibration Test Condition for Air-to-Air Missile[J]. Aeronautic Standardization & Quality, 2006(3): 34—37.

[4] GJB/Z 222—2005, 動力學環境數據采集和分析指南[S]. GJB/Z 222—2005, Guidelines for Dynamic Environmental Data Acquisition and Analysis[S].

[5] 施榮明. 飛機結構振動設計與試驗[M]. 北京: 航空工業出版社, 2014: 110—111. SHI Rong-ming. Aircraft Structure Design and Vibration Test[M]. Beijing: Aviation Industry Press, 2014: 110—111.

[6] BENDAT J S, PIERSOL A G. Random Data: Analysis and Measurement Procedures[M]. 2nd Ed. New York: Wiley, 1986.

[7] BENDAT J S, PIERSOL A G. Engineering Applications of Correlation and Spectral Analysis[M]. 2nd Ed. New York: Wiley, 1993.

[8] FORLIFER W R. The Effects of Filter Bandwidth in the Spectrum Analysis of Random Vibration[R]. Washington D C: Department of Defense,1964.

[9] SCHMIDT R. Resolution Bias Errors in Spectral Density, Frequency Response and Coherence Function Measurements, III: Application to Second-order Systems[J]. Journal of Sound and Vibration, 1985, 101: 377—404.

[10] 蘇成, 徐瑞. 非平穩激勵下結構隨機振動時域分析法[J].工程力學, 2010, 27(12): 77—83. SU Cheng, XU Rui. Random Vibration Analysis of Structures Subjected to Non-stationary Excitations by Time Domain Method[J]. Engineering Mechanics, 2010, 27(12): 77—83.

[11] 張賢達, 保錚. 非平穩信號分析與處理[M]. 北京: 國防工業出版社, 1998. ZHANG Xian-da, BAO Zheng. Analysis Procedures for Non-stationary Signal[M]. Beijing: National Defence Industry Press, 1998.

[12] 張海勇, 李勘. 非平穩隨機信號的參數模型分析方法[J].系統工程與電子技術, 2003, 25(3): 386—390. ZHANG Hai-yong, LI Kan. Analytical Methods for the Non-Stationary Random Signals Based on Parametric Models[J]. Systems Engineering and Electronics, 2003, 25(3): 386—390.

[13] ALLAN G P. Optimum Data Analysis Procedures for Titan IV and Space Shuttle Payload[M]. California: Jet Propulsion Laboratory, California Institute of Technology, 1991.

[14] HARRY H, ALLAN G P, JAMES H W, et al. Handbook of Dynamic Data Acquisition and Analysis[M]. USA: Institute of Environmental Sciences, 1994.

[15] 張海勇. 基于局域波法的非平穩隨機信號分析中的若干問題的研究[D]. 大連: 大連理工大學, 2001. ZHANG Hai-yong. Research on the Problems of NonStationary Random Signal Analysis and Processing Based on Local Wave Method[D]. Dalian: Dalian University of Technology, 2001.

The Research of Processing Method Based on Maximum Spectral for Non-stationary Random Vibration Data

XU Jun, FU Yun, ZHANG Jian-jun, SUN Jian-yong, WANG Xing-hao
（China Aero-Polytechnology Establishment, Beijing 100028, China）

Objective To solve the problem that the environmental measurement of non-stationary random vibration data can’t deal with power spectral density. Methods Two methods of maximum spectral analysis were summed up, namely FFT and digital filter which were based on minimum mean-square error criterion. The advantages, disadvantages and error of two methods with different bandwidth were discussed, and then examples were given to testify the better method. Results The theoretic analysis and simulation results showed that: the error was minimum with conventional FFT computations (proportional resolution bandwidth). The conventional FFT computations were applied to 9 sorties of non-stationary vibration data processing and environmental condition determination for a certain missile, and the result was reasonable and effective. Conclusion The conventional FFT is the best method to deal with non-stationary random vibration data, which provides guidance for formulating environmental condition based on maximum spectral.

non-stationary random vibration data; maximum spectral; FFT; digital filter

2016-08-16； Revised：2016-08-16

10.7643/ issn.1672-9242.2016.05.009

TJ86

A

1672-9242(2016)05-0054-07

2016-08-16；

2016-08-16

徐俊（1990—），男，湖北武漢人，碩士研究生，主要研究方向為裝備環境工程。

Biography：XU Jun（1990—）， Male， from Wuhan， Hubei， Maser graduate student， Research focus: equipment environmental engineering.