師培峰，邱 偉，吳瑞斌，王曉斌

（北京強度環(huán)境研究所 北京 100076）

引言

在各類飛行試驗中，遙測數(shù)據(jù)為制定飛行器環(huán)境試驗條件、檢驗系統(tǒng)工作性能和故障定位提供了重要依據(jù)，對遙測數(shù)據(jù)進行準確、有效處理［1］具有極其重要的意義。由于飛行環(huán)境、電磁干擾及傳感器零位漂移等原因，遙測數(shù)據(jù)往往存在虛假趨勢，并混雜有周期性干擾和噪聲干擾。為改善數(shù)據(jù)品質(zhì)，需要對數(shù)據(jù)進行預處理。此外，由于遙測數(shù)據(jù)具有非線性、非平穩(wěn)（時變）的特點，基于傅里葉變換的傳統(tǒng)頻域分析方法往往無法表述信號的時頻局部特征，而這種性質(zhì)恰恰是非平穩(wěn)信號［2］最根本和最關(guān)鍵的性質(zhì)?；贖HT［3］（Hilbert-Huang Transform）技術(shù)的遙測信號處理方法，可以有效完成信號預處理和時頻分析，為高效處理遙測數(shù)據(jù)提供了有力工具。

1 HHT基本理論

HHT利用經(jīng)驗模態(tài)分解［4］（EMD）將信號自適應分解為一系列的固有模態(tài)函數(shù)（IMF）信號分量，對每個IMF信號進行Hilbert變換得到信號的瞬時頻率和每個IMF信號分量的時頻序列。在HHT理論中，IMF信號分量滿足以下條件:

①整個數(shù)據(jù)序列中，極值點的數(shù)量與過零點的數(shù)量相等或至多差一個；

②信號上任意一點，由局部極大值點確定的包絡(luò)線和局部極小值點確定的包絡(luò)線的均值為零，即信號關(guān)于時間軸對稱。

IMF信號分量通過EMD計算得到。與經(jīng)典的Fourier變換理論中有關(guān)信號的基本組成分量不同，EMD方法認為非線性、非平穩(wěn)信號的基本組成分量不是正余弦信號，而是IMF信號分量。EMD是以IMF為基函數(shù)的時域分解方法。采用EMD方法分解任意一個信號s（t）的具體步驟有如下七步。

①提取信號s（t）的所有極大值和極小值極點；

②利用三次樣條函數(shù)分別基于所有極大值點與極小值點擬合信號s（t）的上包絡(luò)和下包絡(luò)，并用它們近似地表示信號s（t）真實的上、下包絡(luò)線；

③求上包絡(luò)和下包絡(luò)的均值m1（t），作為信號s（t）的真實均值包絡(luò)曲線；

④用信號s（t）減去均值包絡(luò)m1（t）得到新的信號h1（t），即

⑤將h1（t）視為新的s（t），對h1（t）重復步驟①～④。假設(shè)通過k次迭代（k≥2）得到滿足IMF信號分量兩個條件的h1k（t），則h1k（t）是第一個IMF信號分量，即

用字符C（t）表示IMF信號分量，Ci（t）表示第i個IMF信號分量，令

則C1（t）就是信號s（t）的第一個IMF信號分量，它包含了原始信號里最高頻部分的信號。

⑥用原始信號減去第一個IMF信號分量，得到第一個殘余信號R1（t），即

⑦對第一個殘余信號R1（t）重復步驟①～⑥，可以得到每個IMF信號分量Ci（t）及殘余信號Ri（t），即

若殘余信號Ri（t）滿足:Ri（t）或Ci（t）小于設(shè)定的誤差或者Ri（t）為單調(diào)函數(shù)，則EMD分解終止。

在實際的數(shù)據(jù)處理過程中，步驟⑤中迭代次數(shù)k的值通過迭代閾值（標記為SD）來確定，SD的計算如式（6）所示。實踐表明，SD的取值在0.2～0.3之間為宜，既可保證IMF信號分量的線性和穩(wěn)定性，又可使IMF信號分量具有相應的物理意義［5］。

式（6）中，t為信號采樣時刻點，T為總采樣時間。

原始信號s（t）通過EMD分解得到n個IMF信號分量（C1（t），…，Cn（t））和殘余信號Rn（t），根據(jù)EMD方法原理可知

所有的IMF信號分量包含了原始信號中頻率由高到低的部分。每次EMD分解產(chǎn)生的IMF信號的頻率是逐次降低的。

對單個IMF信號分量Ci（t）進行HT（Hilbert變換），可以得到Ci（t）的解析信號，即

其中P.V.表示柯西主值積分。構(gòu)造Ci（t）的解析函數(shù)為

其中幅值函數(shù)和相位函數(shù)分別為

在式（11）的基礎(chǔ)上定義瞬時頻率為

從式（12）可以看出，HT變換后的瞬時頻率是時間的函數(shù)，能夠?qū)θ我庑盘柸我鈺r刻頻率變化進行分析和度量，尤其是非平穩(wěn)信號。通過對每個IMF信號分量進行HT變換，就可以得到每個IMF信號分量的頻率隨時間的變化曲線。

2 基于HHT的信號處理

在對遙測信號進行數(shù)據(jù)分析之前，需要對信號進行預處理。遙測信號的預處理主要包括信號虛假趨勢項排除［6］和信號去噪［7，8］兩方面的工作。

2.1 趨勢項排除

目前趨勢項排除的主要方法有:①高通濾波［9］，濾波截止頻率應小于下限分析頻率的1/2；②最小二乘多項式擬合，提取出趨勢項信號的多項式；③滑動平均，根據(jù)趨勢項特征選擇階次；④小波分析［10］，利用小波分解去除趨勢項后進行波形重構(gòu)。上述四種方法都有特定的使用條件:①高通濾波的性能與阻帶衰減特性、通帶平坦度、窗類型及階數(shù)有關(guān)；②最小二乘多項式擬合的階數(shù)受到限制，而且高階多項式的擬合效果也不理想，需謹慎使用。③滑動平均的趨勢提取效果與點數(shù)相關(guān)，采用不同點數(shù)進行滑動平均差異較大；④小波分析往往需要選擇合適的小波類型和分解層數(shù)，才能取得較好的效果。這四種方法要選擇合適的參數(shù)類型往往比較困難，而且針對特定的數(shù)據(jù)才有較好的趨勢提取與排除效果。本文基于HHT理論，對遙測信號進行EMD分解，得到從高頻到低頻有序排列的多個IMF信號分量。這本質(zhì)上是一個從高頻到低頻的多分辨率自適應濾波［11］過程，分解過程中不存在信號失真問題，具有完備性和可重構(gòu)性。通過對高頻IMF信號分量線性疊加進行信號重構(gòu)，實現(xiàn)低頻趨勢項排除。圖1為一實測遙測振動信號時間歷程圖（g=9.8m/s2），通過EMD分解得到11個相互正交的IMF信號分量，分別記為imf-1、…imf-11，如圖2和圖3所示（t表示time，a表示acceleration）。

圖1 實測遙測振動信號時間歷程圖Fig.1 Time-domain curve of telemetry vibration signal

圖2 EMD分解后的前8組IMF信號Fig.2 The former 8 IMF signals by EMD decomposition

根據(jù)式（7），由EMD分解后的前8個IMF信號得到重構(gòu)后的信號，即排除了趨勢項后的數(shù)據(jù)。分別采用本文方法、高通濾波法和小波變換法進行趨勢項排除，結(jié)果對比如圖4所示。

從圖4中可以看出，從排除趨勢項的效果來看，對于初始值不在零線的數(shù)據(jù)，采用高通濾波法在初始位置會出現(xiàn)擾動，而且人為設(shè)置高通截止頻率具有較大的隨意性；小波變換方法和本文方法都能夠有效去除趨勢項，但小波變換法處理后數(shù)據(jù)的加速度均方根（RMS）值為4.23g，本文方法處理后數(shù)據(jù)的RMS值為4.92g，由于小波方法中小波基的選擇依據(jù)個人經(jīng)驗，因此在排除趨勢項的過程中，會丟失部分真實的高頻趨勢項數(shù)據(jù)信息，而本文方法通過自適應分解實現(xiàn)，在排除虛假趨勢項的同時保留了真實的高頻細節(jié)趨勢項數(shù)據(jù)，更真實地反映了遙測振動的物理過程。

圖3 EMD分解后具有趨勢項的后3組IMF信號Fig.3 The later 3 IMF signals containing trend by EMD decomposition

圖4 三種方法排除趨勢項效果對比Fig.4 Trend extraction and elimination contrast of three methods

2.2 去噪濾波

在飛行遙測試驗中，由于工頻、幀碼、采樣頻率和彈上記錄校準信號的干擾和影響，許多遙測數(shù)據(jù)存在某些頻率及其諧波的干擾周期分量，尤其是在低頻參數(shù)中，因此必須抑制干擾頻率及其諧波頻率的影響。目前抑制工頻及其諧波干擾的方法主要是設(shè)計帶阻濾波器和小波變換分析方法降噪。濾波器法是使信號幅值在需要消除的干擾頻率點上為零。由于干擾頻率往往是一個點頻率，加上濾波器的過渡帶衰減、通帶波動、相位延遲等特性的影響，濾波器法在實際操作中比較復雜，而且很難消除工頻及其諧波頻率的干擾。小波變換作為一種信號去噪方法，其基本思想是首先尋找一個滿足一定條件的母波函數(shù)，然后通過母波函數(shù)的平移和伸縮構(gòu)成小波基，再利用小波基去逼近所要研究的信號，進行時頻局部分析，從而達到去噪的目的。小波變換也具有多分辨特性，但小波變換本質(zhì)上是一種窗口可調(diào)的傅里葉變換，沒有從根本上擺脫傅里葉分析的局限。與上述兩種方法不同的是，經(jīng)EMD分解后的IMF信號分量都是平穩(wěn)的，由于EMD分解是從信號本身的特征時間尺度出發(fā)對信號進行分解，沒有固定的先驗基底，是自適應的［11］，因此得到的IMF信號分量具有明顯的物理意義，表現(xiàn)了信號的真實物理過程。本文基于HHT原理，通過以下四步驟實現(xiàn)濾波去噪:

①對原始信號s（t）進行EMD分解，得到M個IMF信號分量，記為IMFsi（t），i=1…M。

②對第i個IMF信號分量IMFsi（t）進行HT變換，得到IMFsi（t）的時間-頻率曲線，記為Fi（t）。

③從Fi（t）中標記出第i個IMF信號分量的干擾頻率及其諧波頻率，將Fi（t）中其他頻率對應的IMF信號分量序列標記為IMF（1）si（t），i=1…N，N＜M。

④去除干擾頻率后的信號s（t）通過IMF（1）si（t）實現(xiàn)信號重構(gòu)，即

通過疊加不同的IMF（1）si（t）信號分量，可以得到所需類型的濾波器（低通、高通、帶通、帶阻等），從而實現(xiàn)對遙測信號的自適應濾波和去噪分析。以實測信號為例，信號采樣頻率為1024Hz，信號成分主要是40Hz以內(nèi)的低頻部分，在采集過程中有較大的背景噪聲，并夾雜有50Hz及其倍頻諧波的干擾。時域波形及頻譜分析如圖5所示。

圖5 遙測振動信號時間歷程及頻譜Fig.5 Time-domain curve and frequency spectrum of telemetry vibration signal

對圖5（a）所示的實測遙測振動信號進行譜分析，得到頻譜圖，如圖5（b）所示。從圖5（b）中可以看出，振動信號中夾雜著50Hz及其倍頻諧波的干擾。分別采用數(shù)字濾波器法、小波變換法和本文方法進行濾波去噪處理，并分兩種情況進行討論。

①只保留40Hz以內(nèi)頻帶的信號。

三種方法濾波去噪效果如圖6所示。

圖6 保留40Hz以內(nèi)頻帶信號頻譜Fig.6 Frequency spectrum by keeping the signal within 40Hz

②只消除50Hz及其倍頻諧波的干擾。

三種方法濾波去噪效果如圖7所示，局部細節(jié)如圖8（a）、8（b）所示。

圖7 去除50Hz及諧波干擾后的信號頻譜Fig.7 Frequency spectrum after eliminating the interference of 50Hz and its harmonic frequency

從圖6對比結(jié)果可以看出，本文方法能夠有效消除頻率大于40Hz的信號。IIR濾波器由于過渡帶衰減斜率的影響，無法完全消除50Hz工頻信號的干擾。小波變換方法由于小波基不是自適應的，也無法完全消除50Hz及其諧波信號的干擾。

從圖7和圖8對比結(jié)果可以看出:采用帶阻濾波器方法，在50Hz及其諧波頻率處出現(xiàn)局部信號失真，表現(xiàn)為50Hz及其諧波頻率處頻譜出現(xiàn)谷值，而且主頻范圍（40Hz以內(nèi)）內(nèi)信號幅值明顯降低，說明有效信號的部分信息丟失；小波分析方法與本文方法的頻譜圖在主頻范圍內(nèi)保持較好的一致性，但是小波分析方法無法完全消除50Hz及其諧波頻率信號的干擾，而本文方法在消除50Hz及其諧波頻率干擾的同時，保證了主頻范圍內(nèi)信號不丟失。

圖8 局部細節(jié)展開圖Fig.8 Curve details of figure 7

3 結(jié)束語

本文研究了HHT在遙測信號處理中的應用。仿真實驗結(jié)果表明，基于HHT的數(shù)據(jù)處理技術(shù)可以自適應對信號進行分解，它結(jié)合了EMD分解和IMF信號分量的特性，有效地實現(xiàn)了信號趨勢項提取和排除以及自適應信號濾波降噪處理。本文方法應用前景廣闊，可為各種信號的數(shù)據(jù)處理提供一種新的有力工具。

［1］陳以恩，等.遙測數(shù)據(jù)處理［M］.北京:北京國防工業(yè)出版社，2002.

［2］科恩 L.時-頻分析:理論與應用［M］.白居憲，譯.西安:西安交通大學出版社，2000:22～40.

［3］Huang N E，Shen Z，Long SR，et al.The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non stationary Time Series and Analysis［J］.Pro.Roy.Soc.London.A，1998，454:903 ～995.

［4］Peng Z K，Tse PeterW，Chu F L.An Improved Hilbert-Huang Transform and Its Applicaition in Bibration Signal Analysis［J］.Journal of Sound and Vibration，2005，286:187 ～205

［5］Rilling G，F(xiàn)landrin P.On Empirical Mode Decomposition and Its Algorithm［C］//IEEE-ERUASIP Workshop on Nonlinear Signal Image Processing.NSIP-03，Gradp（I），June 2003.

［6］陳 燕，劉 哲，鄭 賓.基于LabVIEW的測試信號預處理方法研究［J］.國外電子測量技術(shù)，2008，27（10）:4～5，16.Chen Yan，Liu Zhe，Zheng Bin.Study on Test Signal Pre-processing Method Based on LabVIEW［J］.Foreign Electronic Measurement Technology，2008，27（10）:4～5，16.

［7］肖方煜，湯 偉，傅 娜.自尋優(yōu)閾值小波去噪方法［J］.信號處理，2012，28（4）:577～586.Xiao Fangyu，Tang Wei，F(xiàn)u Na.Wavelet Based De-noising Self-optimizing Method［J］.Signal Processing，2012，28（4）:577～586.

［8］趙瑞珍.小波理論及其在圖像、信號處理中的算法研究［D］.西安:西安電子科技大學，2001.Zhao Ruizhen.Study on Wavelet Theory and Its Algorithms for Image and Signal Processing［D］.Xi'an:Xidian University，2001.

［9］程佩青.數(shù)字信號處理［M］.北京:清華大學出版社，2001.

［10］劉明才.小波分析及其應用［M］.北京:清華大學出版社，2005.

［11］蓋 強，張海勇，徐曉剛.Hilbert-Huang變換的自適應頻率多分辨分析研究［J］.電子學報，2005，33（3）:563～566.Ge Qiang，Zhang Haiyong，Xu Xiaogang.Study of Adaptive Frequency Multiresolution Analysis of the Hilbert-Huang Transform［J］.Acta Electronica Sinica，2005，33（3）:563 ～566.