基于時(shí)頻分析的運(yùn)載火箭故障特征提取技術(shù)

崔桂磊 唐永華

(宇航動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，西安 710043)

基于時(shí)頻分析的運(yùn)載火箭故障特征提取技術(shù)

崔桂磊 唐永華

(宇航動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，西安 710043)

為盡可能從運(yùn)載火箭飛行試驗(yàn)數(shù)據(jù)中獲取更多有價(jià)值的信息，提出了新的火箭振動(dòng)參數(shù)故障特征提取模型，用于火箭振動(dòng)問(wèn)題分析.構(gòu)建了三維傅里葉譜，用于高效地分析全程觀測(cè)信號(hào)頻譜的時(shí)變趨勢(shì)，研究了Wigner-Ville分布的標(biāo)準(zhǔn)算法與快速算法，用于精確刻畫(huà)振動(dòng)參數(shù)感興趣數(shù)據(jù)段的頻譜隨時(shí)間的分布情況.試驗(yàn)結(jié)果表明:模型及方法處理火箭遙測(cè)數(shù)據(jù)直觀、高效、精確，可以應(yīng)用于實(shí)測(cè)數(shù)據(jù)分析.研究時(shí)頻分析方法對(duì)試驗(yàn)數(shù)據(jù)進(jìn)行更深層次的分析處理，對(duì)設(shè)計(jì)高可靠、高性能的運(yùn)載火箭，具有重要意義.

運(yùn)載火箭;遙測(cè)數(shù)據(jù);故障定位;時(shí)頻分析;特征提取

飛行試驗(yàn)遙測(cè)數(shù)據(jù)反映運(yùn)載火箭的工作環(huán)境和各系統(tǒng)、各部件的工作狀態(tài)，對(duì)于分析飛行試驗(yàn)過(guò)程是否正常以及分析飛行中出現(xiàn)的各種異常現(xiàn)象有著重要作用［1］.目前，火箭故障診斷的研究主要包括四個(gè)方面的內(nèi)容:故障機(jī)理研究、狀態(tài)信號(hào)采集、信號(hào)分析處理和故障特征提取、故障狀態(tài)識(shí)別和診斷決策，其中的核心問(wèn)題是故障信號(hào)的特征提取，其處理結(jié)果是故障診斷的重要依據(jù).無(wú)論從地面試車狀態(tài)還是飛行的實(shí)際來(lái)看，火箭遙測(cè)數(shù)據(jù)中往往含有大量的時(shí)變、短時(shí)突發(fā)性質(zhì)的成分，其頻率在故障發(fā)生時(shí)刻是跳變的，這使得目前基于傳統(tǒng)線性和平穩(wěn)理論發(fā)展起來(lái)的遙測(cè)數(shù)據(jù)處理技術(shù)如傅里葉變換、現(xiàn)代譜估計(jì)等［2］，難以有效地揭示故障信號(hào)的非平穩(wěn)特征，因?yàn)樗鼈兌际菍?duì)測(cè)量數(shù)據(jù)的整體變換，即對(duì)信號(hào)的表征或是完全在時(shí)域，或是完全在頻域，只能對(duì)信號(hào)進(jìn)行時(shí)域或者頻域分析，作為頻域表示的功率譜并不能說(shuō)明其中某種頻率分量出現(xiàn)在什么時(shí)候及其變化情況，故無(wú)法滿足遙測(cè)故障信號(hào)的分析處理需要.在這種情況下，必須研究建立新的故障特征提取模型進(jìn)行飛行試驗(yàn)數(shù)據(jù)的處理，以提高故障信號(hào)特征提取的品質(zhì)，為進(jìn)一步提高對(duì)試驗(yàn)數(shù)據(jù)的認(rèn)識(shí)水平，進(jìn)行型號(hào)任務(wù)火箭狀態(tài)分析與故障定位提供保障.

1 兩種時(shí)頻分析方法

時(shí)頻分析被認(rèn)為是近年來(lái)對(duì)以傅里葉變換為基礎(chǔ)的線性和平穩(wěn)數(shù)據(jù)分析的一個(gè)重大突破，產(chǎn)生的原因主要是由于傅里葉變換在分析非平穩(wěn)信號(hào)方面存在局限性，不能有效分析信號(hào)的瞬時(shí)特性［3－4］.時(shí)頻分析提供了信號(hào)從時(shí)域到時(shí)間-頻率域的變換，通過(guò)二維和三維時(shí)頻分布圖，能夠在時(shí)頻平面上直觀表示出故障信號(hào)中各個(gè)分量的時(shí)間關(guān)聯(lián)譜特性，在每一時(shí)刻指示出信號(hào)在瞬時(shí)頻率附近的能量聚集情況，從而同時(shí)掌握信號(hào)的時(shí)域、頻域信息，清楚地了解信號(hào)的頻譜、能量譜是如何隨時(shí)間變化的［5－6］.作為一種高效、高分辨率的數(shù)據(jù)處理技術(shù)，時(shí)頻分析理應(yīng)在飛行試驗(yàn)數(shù)據(jù)的使用和挖掘方面得到更好的運(yùn)用.時(shí)頻分析技術(shù)的基本思想是設(shè)計(jì)時(shí)間和頻率聯(lián)合函數(shù)，將一維時(shí)域信號(hào)映射到二維時(shí)頻平面，在二維時(shí)間-頻域平面全面分析觀測(cè)信號(hào)的時(shí)頻聯(lián)合特征，實(shí)現(xiàn)對(duì)信號(hào)頻譜或能量譜的時(shí)間定位［7－9］.按照時(shí)頻聯(lián)合函數(shù)的不同可以分為線性時(shí)頻分析和非線性時(shí)頻分析兩類，其時(shí)頻分辨率和聚集性各不相同［10－11］.因此工程實(shí)際中要研究最優(yōu)方法來(lái)處理遙測(cè)信號(hào)，以獲得高效、高分辨率的火箭飛行環(huán)境分析效果.

1.1 三維傅里葉譜

結(jié)合運(yùn)載火箭遙測(cè)振動(dòng)參數(shù)的特性，提出了全程數(shù)據(jù)三維傅里葉譜處理方案:

1)將振動(dòng)參數(shù)在飛行試驗(yàn)的全程時(shí)域觀測(cè)數(shù)據(jù)以1 024點(diǎn)為一個(gè)數(shù)據(jù)幀，劃分為連續(xù)的若干幀信號(hào)，每個(gè)數(shù)據(jù)幀在時(shí)域占用的時(shí)間片長(zhǎng)度為Δt=1024/fsr，fsr為該參數(shù)的采樣率.

2)對(duì)每個(gè)數(shù)據(jù)幀進(jìn)行快速傅里葉變換FFT(Fast Fourier Transform)變換得到其頻譜結(jié)構(gòu)，頻譜的分辨率Δf=fsr/1024.

3)構(gòu)建三維傅里葉譜的原理如圖1所示，圖中X為時(shí)間序列，以起測(cè)時(shí)間T0為時(shí)間原點(diǎn)，以時(shí)間片長(zhǎng)度Δt為時(shí)間間隔逐幀增加;Y為頻率序列，以譜分析的下限頻率(一般為0到5Hz)為原點(diǎn)，以頻譜分辨率Δf為頻率間隔;Z為頻譜幅值.

1.2 Wigner-Ville分布

在所有非線性時(shí)頻表示中，Wigner-Ville分布(WVD)［12］是非常重要的一類二次時(shí)頻分布，WVD的時(shí)頻帶寬積達(dá)到了Heisenberg不確定性原理給出的下界，使得WVD比線性時(shí)頻分析具有更高的時(shí)頻分辨率，時(shí)頻聚集性比較好，并在許多方面如地震信號(hào)處理、故障特征提取、高光譜探測(cè)等［13－15］取得了成功的應(yīng)用.

圖1 三維傅里葉譜構(gòu)造原理

1.2.1 WVD 的原理

設(shè)觀測(cè)信號(hào)為s(t)，它的WVD表示為

式中，t為時(shí)間變量;f為頻率變量;x(t)是信號(hào)s(t)的解析信號(hào)，即 x(t)=s(t)+j H［s(t)］，其中H［］表示希爾伯特變換;*表示復(fù)共軛.

式(1)中，信號(hào)x(t)出現(xiàn)了兩次，即WVD以雙線性或非線性的形式將信號(hào)從時(shí)域轉(zhuǎn)換到時(shí)頻域，描述了信號(hào)能量譜隨時(shí)間、頻率的變化規(guī)律，且不含任何窗函數(shù)，這樣就避免了線性時(shí)頻分析中時(shí)間和頻率分辨率的相互牽制.

1.2.2 WVD標(biāo)準(zhǔn)算法與快速算法

1)WVD標(biāo)準(zhǔn)算法.

設(shè)觀測(cè)信號(hào)為 s(n)，n=1，2，…，N，其離散解析信號(hào)為 x(n)，n=1，2，…，N，則其離散 WVD 可以表示為

式中，WM=e－j2π/M，M 為窗口長(zhǎng)度.

分析式(2)可知，離散WVD可以通過(guò)FFT變換來(lái)實(shí)現(xiàn).

式(2)可以視為WVD標(biāo)準(zhǔn)算法.標(biāo)準(zhǔn)算法將長(zhǎng)度為N的觀測(cè)信號(hào)作為一個(gè)整體進(jìn)行WVD變換.作為信號(hào)的二次變換，WVD具有分辨率高、時(shí)頻聚集性好的優(yōu)點(diǎn)，但隨著N的增大，時(shí)效性變差.

2)WVD快速算法.

當(dāng)遙測(cè)信號(hào)長(zhǎng)度N較大時(shí)，應(yīng)用WVD標(biāo)準(zhǔn)算法進(jìn)行時(shí)頻分析，其處理效率往往不能滿足工程快速分析的要求，因此實(shí)際中提出了快速算法方案.具體方案如下:①將待處理數(shù)據(jù)x(n)按1024點(diǎn)為1幀劃分為n=N/1024幀，xi(n)為第i幀數(shù)據(jù);②對(duì)劃分好的數(shù)據(jù)逐幀進(jìn)行WVD變換Φ［xi(n)］=Φxi(n，k);③在時(shí)頻域?qū)⒚繋幚斫Y(jié)果進(jìn)行拼接.

快速算法可以有效提高長(zhǎng)數(shù)據(jù)序列的處理效率，與之對(duì)應(yīng)的是時(shí)頻分辨率有所降低，且頻譜在幀間分割處會(huì)出現(xiàn)斷痕，但一般不會(huì)影響對(duì)振動(dòng)參數(shù)的正確分析.

2 火箭振動(dòng)參數(shù)故障特征提取模型

根據(jù)測(cè)量目的不同，火箭上多個(gè)部位安置了不同類型的力學(xué)環(huán)境參數(shù)傳感器，一般可分為振動(dòng)參數(shù)、沖擊參數(shù)、噪聲參數(shù)和脈動(dòng)壓力參數(shù)等，其中振動(dòng)參數(shù)又包括POGO振動(dòng)、低頻振動(dòng)和高頻振動(dòng)［1］.由于運(yùn)載火箭發(fā)射成本相當(dāng)高昂，不可能過(guò)于頻繁地進(jìn)行試驗(yàn)，因此遙測(cè)數(shù)據(jù)具有很高的價(jià)值.環(huán)境參數(shù)處理的任務(wù)就是對(duì)測(cè)量數(shù)據(jù)進(jìn)行分析，描述它們?cè)跁r(shí)域、頻域的基本特性，為火箭的設(shè)計(jì)評(píng)測(cè)及故障診斷提供可靠的依據(jù).

對(duì)于長(zhǎng)度為N的火箭遙測(cè)數(shù)據(jù)，經(jīng)過(guò)大量數(shù)據(jù)試驗(yàn)，建立了如圖2所示的振動(dòng)參數(shù)故障特征提取模型，描述如下:

1)全程數(shù)據(jù)預(yù)處理.讀取待分析參數(shù)的全程數(shù)據(jù)文件，對(duì)全程數(shù)據(jù)進(jìn)行中心化、濾波等預(yù)處理工作，并繪制其時(shí)域波形圖.

圖2 火箭振動(dòng)參數(shù)故障特征提取模型

2)全程數(shù)據(jù)三維傅里葉譜分析.全程數(shù)據(jù)以每1024點(diǎn)為1個(gè)數(shù)據(jù)幀，逐幀進(jìn)行FFT變換，并在頻率-時(shí)間-幅值三維空間繪制全程三維傅里葉譜，觀測(cè)全程數(shù)據(jù)頻譜變化的大概趨勢(shì).

3)數(shù)據(jù)選段.根據(jù)遙測(cè)數(shù)據(jù)處理要求及全程三維傅里葉譜，選取感興趣的數(shù)據(jù)段落，長(zhǎng)度為1024×m，m為幀數(shù).對(duì)選段數(shù)據(jù)進(jìn)行傅里葉變換，觀察分析頻譜結(jié)構(gòu)，根據(jù)頻譜信息進(jìn)行時(shí)域數(shù)據(jù)濾波.

4)選段數(shù)據(jù)的WVD分析處理.根據(jù)幀數(shù)m是否小于等于4，選取WVD標(biāo)準(zhǔn)算法或快速算法.觀察WVD二維時(shí)頻譜、三維能量譜，分析故障信號(hào)特征頻率隨時(shí)間的分布情況.

5)進(jìn)一步分析處理.若需要獲取更細(xì)微的信息，則回到步驟3)，進(jìn)一步進(jìn)行數(shù)據(jù)選段并進(jìn)行WVD時(shí)頻分析.

3 火箭振動(dòng)參數(shù)故障特征提取實(shí)例

以某次試驗(yàn)任務(wù)某型號(hào)運(yùn)載火箭低頻振動(dòng)參數(shù)ZD01觀測(cè)數(shù)據(jù)為例，分析如下:

1)全程數(shù)據(jù)預(yù)處理.遙測(cè)參數(shù)ZD01測(cè)量時(shí)間段為 －29.10 ～200.03 s，全程數(shù)據(jù)經(jīng)中心化及5～100Hz濾波后，其時(shí)域波形如圖3所示.

圖3 振動(dòng)參數(shù)全程數(shù)據(jù)時(shí)域圖

圖4 振動(dòng)參數(shù)全程三維傅里葉譜

2)全程數(shù)據(jù)三維傅里葉譜分析.通過(guò)圖4所示的全程三維傅里葉譜，可觀察到120 s及150 s附近頻率成分極為豐富.同時(shí)可以大致判斷參數(shù)ZD01在0～100 s有較為明顯的40Hz左右的振動(dòng)頻率存在，這就是型號(hào)任務(wù)專家高度重視的低頻振動(dòng)問(wèn)題.但由于該方法為逐幀分析繪圖，每幀1024點(diǎn)，因此時(shí)間分辨率與頻率分辨率均較低，僅適用于遙測(cè)全程數(shù)據(jù)譜結(jié)構(gòu)的粗略分析.

3)數(shù)據(jù)選段.選取參數(shù) ZD01在 －1.23～101.15 s范圍內(nèi)的N=32768點(diǎn)(32幀)遙測(cè)數(shù)據(jù)進(jìn)行時(shí)頻分析處理，其時(shí)域波形如圖5所示.傳統(tǒng)傅里葉變換的頻譜如圖6所示，可以看出存在明顯的40Hz頻率信號(hào)，但是，圖6的傅里葉頻譜是對(duì)－1.23～101.15 s數(shù)據(jù)從時(shí)域到頻域的整體變換，它只能說(shuō)明在 －1.23 ～101.15 s中40Hz信號(hào)存在，卻不能進(jìn)一步給出關(guān)于40 Hz頻率出現(xiàn)的時(shí)間以及40Hz頻譜、能量譜隨時(shí)間的分布情況.

圖5 選段數(shù)據(jù)時(shí)域圖

圖6 選段數(shù)據(jù)頻譜圖

4)選段數(shù)據(jù)的WVD分析處理.對(duì)選段數(shù)據(jù)做WVD變換，將信號(hào)分析從時(shí)域轉(zhuǎn)換到時(shí)間-頻率域，在時(shí)頻域全面觀測(cè)信號(hào)的聯(lián)合時(shí)頻特性.由于選取的數(shù)據(jù)較長(zhǎng)，因此應(yīng)用WVD快速算法進(jìn)行時(shí)頻分析處理.圖7為選段數(shù)據(jù)WVD的二維時(shí)頻譜圖，可以清晰地看出，參數(shù)ZD01在10～60 s存在40Hz信號(hào);圖8為選段數(shù)據(jù)WVD的三維能量譜圖，進(jìn)一步描述了40 Hz頻率信號(hào)的能量隨時(shí)間的分布情況.

5)進(jìn)一步分析處理.選取參數(shù)ZD01在10～60 s范圍內(nèi)的遙測(cè)數(shù)據(jù)可以獲得40 Hz頻率更精確的時(shí)變特性;選取120 s與150 s附近的遙測(cè)數(shù)據(jù)，可以分別得到火箭一級(jí)關(guān)機(jī)與一二級(jí)分離時(shí)刻的振動(dòng)特性，本文不再贅述.

6)振動(dòng)參數(shù)ZD01的時(shí)變特性分析結(jié)論.由全程三維傅里葉譜、WVD的二維時(shí)頻譜及三維能量譜全面觀測(cè)信號(hào)的聯(lián)合時(shí)頻特性，可以準(zhǔn)確判斷出40Hz振動(dòng)頻率發(fā)生與消失的時(shí)間及其變化趨勢(shì).

圖7 選段數(shù)據(jù)WVD的二維時(shí)頻譜

圖8 選段數(shù)據(jù)WVD的三維能量譜

4 結(jié)論

1)著眼于長(zhǎng)征系列運(yùn)載火箭振動(dòng)問(wèn)題分析的需要，建立了新的故障特征提取模型，從聯(lián)合時(shí)頻分布角度直觀描述了振動(dòng)參數(shù)特征頻率的頻譜與能量譜隨時(shí)間的變換趨勢(shì)，提高了對(duì)試驗(yàn)數(shù)據(jù)的認(rèn)識(shí)水平.

2)針對(duì)遙測(cè)振動(dòng)參數(shù)時(shí)變特性，構(gòu)建了三維傅里葉譜，研究了Wigner-Ville分布標(biāo)準(zhǔn)算法及快速算法，實(shí)現(xiàn)了箭遙振動(dòng)參數(shù)從單純的時(shí)域或頻域分析到時(shí)頻域分析的轉(zhuǎn)變，提高了數(shù)據(jù)處理能力.

3)實(shí)測(cè)數(shù)據(jù)試驗(yàn)表明，新處理模型直觀、高效，可提供更精確、更豐富的診斷信息，其方法結(jié)果參考價(jià)值更高，能夠?yàn)樾吞?hào)部門進(jìn)行故障狀態(tài)識(shí)別和診斷決策提供更為科學(xué)的依據(jù).

References)

［1］陳以恩，李連周，唐永華，等.遙測(cè)數(shù)據(jù)處理［M］.北京:國(guó)防工業(yè)出版社，2002:3－20 Chen Yien，Li Lianzhou，Tang Yonghua，et al.Telemetry data processing［M］.Beijing:National Defence Industry Press，2002:3－20(in Chinese)

［2］GJB2238A—2004遙測(cè)數(shù)據(jù)處理［S］GJB2238A—2004 Telemetry data processing［S］(in Chinese)

［3］ Chaparro L F，Akan A.Special issue on time-frequency signal analysis and its applications［J］.Journal of the Franklin Institute，2000，337(4):287 －289

［4］ Samantaray L，Dash M，Panda R.A review on time-frequency，time-scale and scale-frequency domain signal analysis［J］.IETE Journal of Research，2005，51(4):287 －293

［5］ Stankovic S.Time-frequency analysis and its application in digital watermarking［J］.EURASIP Journal on Advances in Signal Processing，2010(2):1 －20

［6］ Sejdic E，Djurovic I，Jiang J.Time-frequency feature representation using energy concentration:an overview of recent advances［J］.Digital Signal Processing，2009，19(1):153 － 183

［7］ Ristic B，Boashash B.Kernel design for time-frequency signal analysis using the Radon transform［J］.IEEE Transactionson Signal Processing，1993，41(5):1996 －2008

［8］ Barkat B，Boashash B.A high-resolution quadratic time-frequency distribution for multicomponent signals analysis［J］.IEEE Transactions on Signal Processing，2001，49(10):2232 －2239

［9］ Shafi，I，Ahmad J，Shah，S I，et al.Techniques to obtain good resolution and concentrated time-frequency distributions:a review［J］.EURASIP Journal on Advances in Signal Processing，2009(2):1－43

［10］ Hlawatsch F，Boudreaux-Bartels G F.Linear and quadratic timefrequency signal representations［J］.IEEE Signal Processing，1992，9(2):21 －67

［11］ Jeffrey C，Patrick F.Virtues and vices of quadratic time-frequency distributions［J］.IEEE Transactions on Signal Processing，2000，48(9):2641 －2650

［12］ Cohen L.Time-frequency distributions-a review［J］.Proceedings of the IEEE，1989，77(7):941 －981

［13］ Boashash B，Black P J.An efficient real-time implementation of the Wigner-Ville distribution［J］.IEEE Transactions on Acoustics Speech and Signal Processing，1987，35(11):1611 －1618

［14］ Xu Y N，Zhao Y，Jin C F，et al.Salient target detection based on pseudo-Wigner-Ville distribution and Renyi entropy［J］.Optics Letters，2010，35(4):475 －477

［15］ Climente-Alarcon V，Antonino-Daviu JA，Riera M，et al.Wigner-Ville distribution for the detection of high-order harmonics due to rotor asymmetries［C］//IEEE International Symposium on Diagnostics for Electric Machines，Power Electronics and Drives.［s.l.］:IEEE，2009:119 －124

(編 輯:李 晶)

Rocket fault feature extraction based on time-frequency analysis

Cui Guilei Tang Yonghua

(State Key Laboratory of Astronautic Dynamics，Xian 710043，China)

To obtain more valuable information from rocket flight test data asmuch as possible，a fault feature extraction model was proposed to the need of rocket vibration analysis.3D Fourier spectrum was constucted to analyze the time-varying trend of entire observation signal efficiently.Standard and fast algorithm of Wigner-Ville distribution was studied to describe the frequency spectrum of interesting data segments accurately.Themodel and methodswere illustrated on rocket telemetry data.The results prove the model's direct-vision，high efficiency，and accuracy.It's of great significance for designing a rocket with high reliability and performance to analyze the flight test data more deeply，based on time-frequency analysis methods.

rocket;telemetry data;fault location;time-frequency analysis;feature extraction

V 557+.3;TJ 760.6+28;TN 911.72

A

1001-5965(2012)06-0741-04

2011-03-10;網(wǎng)絡(luò)出版時(shí)間:2012-06-15 15:42

www.cnki.net/kcms/detail/11.2625.V.20120615.1542.016.htm l

總裝預(yù)研項(xiàng)目(2009SY2603009)

崔桂磊(1976 －)，男，山東青島人，工程師，CGL11@mails.tsinghua.edu.cn.