耿 宏，郗厚山

(中國民航大學 電子信息與自動化學院，天津 300300)

0 引 言

大型民航客機的飛行數(shù)據(jù)在飛行事故分析、視情維修、異常檢測、飛行品質(zhì)監(jiān)控等方面起到重要作用[1,2]。由于飛行環(huán)境、飛機自身運動狀態(tài)的復雜性及機載設備的電磁干擾等因素的影響，實際記錄的飛行數(shù)據(jù)存在干擾噪聲。因此，對飛行數(shù)據(jù)的有效去噪是其后續(xù)分析和各領域應用的重要前提。

傳統(tǒng)的飛行數(shù)據(jù)濾波方法如低通濾波、卡爾曼濾波[3]、滑動多項式回歸模型等，由于飛行數(shù)據(jù)具有非線性、非平穩(wěn)的特征，處理后的數(shù)據(jù)在瞬變信號處容易產(chǎn)生局部震蕩，去噪效果不佳?；谛〔ㄗ儞Q的飛行數(shù)據(jù)濾波目前已被較廣泛應用[4,5]。文獻[4]采取平穩(wěn)小波去噪的方式濾除燃油流量、發(fā)動機N1等飛行參數(shù)的野值和噪聲。基于小波變換的飛行數(shù)據(jù)濾波可以取得較好的去噪效果。但是需要依據(jù)經(jīng)驗或者進行大量的實驗以選取恰當?shù)男〔ɑ?、分解層?shù)，無自適應性。

自適應完備集合經(jīng)驗模態(tài)分解(complete ensemble empirical mode decomposition with adaptive noise，CEEMDAN)是一種基于數(shù)據(jù)自身時間尺度特征，將信號分解為固有模態(tài)函數(shù)(intrinsic mode function，IMF)的方法[6-8]。散布熵(dispersion entropy，DE)[9]可以用來分析時間序列的復雜性，它為判斷各IMF分量的噪聲含量提供了依據(jù)。故提出一種基于CEEMDAN和DE的飛行數(shù)據(jù)濾波方法，其具有自適應的優(yōu)點，無需數(shù)據(jù)的眾多先驗知識。經(jīng)該方法去噪后的飛行數(shù)據(jù)平滑性和相似性較好，滿足實際應用的需求。

1 算法原理

1.1 CEEMDAN分解

CEEMDAN算法的基礎是經(jīng)驗模態(tài)分解(empirical mode decomposition，EMD)，但它彌補了EMD算法具有模態(tài)混疊的不足，且完備性較好[10]。CEEMDAN算法的具體分解過程為：

(1)求解時間序列x(t)的第一階IMF和余項。向x(t)附加一系列自適應的輔助白噪聲，構(gòu)成新的時間序列xi(t)

xi(t)=x(t)+ω0ξi(t)，i∈(1,2,…,I)

(1)

式中：xi(t)為第i次添加白噪聲后的時間序列，ω0為噪聲系數(shù)，ξi(t)為白噪聲。I為添加輔助白噪聲的次數(shù)。

(2)

(3)

(2)向余項r1(t)中繼續(xù)添加若干輔助白噪聲，構(gòu)造一組新的時間序列r1i(t)

r1i(t)=r1(t)+ω1E1(ξi(t)),i∈(1,2,…,I)

(4)

采用EMD同(1)對r1i(t)進行處理，求得x(t)階次為2的IMF以及余項，如式(5)、式(6)所示。記Ej(·)為EMD分解得到的第j個IMF分量

(5)

(6)

(3)重復執(zhí)行式(7)、式(8)所示遞推關系式，直到最后的余項序列不能繼續(xù)進行EMD分解時停止，求得其余的IMF

(7)

(8)

最終時間序列信號x(t)被分解為

(9)

K為x(t)的最大分解階次，R(t)為x(t)經(jīng)CEEMDAN最后分解得到的余項。

1.2 飛行數(shù)據(jù)噪聲分析

民航飛機的飛行數(shù)據(jù)主要通過快速存取記錄器(quick access recorder，QAR)、飛行數(shù)據(jù)記錄器(flight data recorder，F(xiàn)DR)等機載設備進行記錄[1]。記錄器按時間順序采樣，記錄了起飛、爬升、巡航等各個飛行階段飛機實際的性能數(shù)據(jù)和狀態(tài)數(shù)據(jù)。但是，不同飛行階段的數(shù)據(jù)變化特征具有較大差異，例如在起飛階段空速、氣壓高度、發(fā)動機N1等參數(shù)變化快且幅度大，而在巡航階段其變化相對穩(wěn)定。這要求飛行數(shù)據(jù)濾波算法具有多尺度分析的能力。

此外，飛行參數(shù)包括大氣環(huán)境參數(shù)、發(fā)動機參數(shù)、氣動參數(shù)等不同類型的參數(shù)，這些參數(shù)中包含的噪聲類型和表現(xiàn)形式也不盡相同。例如：周期性噪聲常規(guī)情況下是由飛機各種電氣干擾及某些部件的周期性機械運轉(zhuǎn)導致，體現(xiàn)為數(shù)據(jù)中的窄譜峰；脈沖噪聲具有持續(xù)時間短和幅度大的特點，一般由電磁干擾等隨機性干擾引起，體現(xiàn)為數(shù)據(jù)中的野值；導致飛行數(shù)據(jù)中存在寬帶噪聲的誘因較多，比如大氣紊流的干擾、量化噪聲等，體現(xiàn)為短時間內(nèi)頻繁抖動的數(shù)據(jù)。因此，噪聲對飛行數(shù)據(jù)的干擾在數(shù)據(jù)中主要表現(xiàn)為以下兩方面：一是存在野值，即數(shù)據(jù)中的異常點；二是存在短時間內(nèi)頻繁抖動的數(shù)據(jù)，其變化趨勢呈現(xiàn)鋸齒狀。上述兩種現(xiàn)象為飛行數(shù)據(jù)的實際應用帶來較大阻礙。

總體上，飛行數(shù)據(jù)中包含的噪聲大多數(shù)為高頻干擾信號，而民航飛機作為一個六自由度的剛體動力學系統(tǒng)，其運動頻率相對較低，通常不會超過10 Hz。因此，對飛行數(shù)據(jù)的去噪相當于是一個低通濾波的過程，但是使用傳統(tǒng)的低通濾波器去噪，其截止頻率難以確定，需要針對不同的飛行參數(shù)進行大量實驗，且去噪效果不是很理想。利用小波閾值去噪的方式對飛行數(shù)據(jù)濾波同樣存在上述問題，在缺乏飛行數(shù)據(jù)先驗知識的情況下，小波基和分解層數(shù)的選取需要通過大量的實驗以及結(jié)果的比對才能確定。有時為了獲得更好的濾波效果，針對不同類型的飛行參數(shù)需要重新選取小波基和分解層數(shù)，實際應用的計算量較大。

從算法的原理可知，CEEMDAN算法可以基于飛行數(shù)據(jù)自身的時間尺度實現(xiàn)飛行數(shù)據(jù)的自適應分解，得到一組從高頻到低頻分布的IMF分量。這些IMF分別包含了飛行數(shù)據(jù)中不同頻率成分的信息，相當于對飛行數(shù)據(jù)按照頻率的高低進行了劃分。此外，CEEMDAN算法具有自適應性，可適用于不同類型的飛行參數(shù)，且不受飛行階段的影響。

1.3 基于CEEMDAN的飛行數(shù)據(jù)去噪原理

利用CEEMDAN算法能夠簡便且高效地獲得包含不同頻率信息的飛行數(shù)據(jù)IMF。依據(jù)對飛行數(shù)據(jù)的噪聲分析，數(shù)據(jù)中的干擾信號集中存在于頻率較高的IMF中，而有用信號大多數(shù)存在于頻率較低的IMF中。因此，CEEMDAN分解階次較低的IMF所含有的有用信息較少且容易被噪聲完全掩蓋。借鑒低通濾波的思想，將前N-1階頻率成分較高的IMF直接去除，使用其余階次的低頻IMF和余項疊加，獲得濾波后的飛行數(shù)據(jù)。在保留大部分有用信息的同時，實現(xiàn)去除飛行數(shù)據(jù)中野值和噪聲的目的。記Ik(t)為飛行數(shù)據(jù)經(jīng)CEEMDAN分解得到的第k階IMF。則此方式濾波后的飛行數(shù)據(jù)s(t)，其表達式如式(11)所示

(10)

(11)

但是，直接剔除前N-1階高頻IMF的方式存在兩個問題：①界限N需要主觀設定，存在不確定性；②直接剔除的方式可以有效去除飛行數(shù)據(jù)中的多數(shù)高頻干擾噪聲，但是仍有部分寬帶噪聲存在于包含有用信號分量中，影響最終的去噪效果。

故本文在使用CEEMDAN算法分解飛行數(shù)據(jù)后，借助散布熵會隨著時間序列復雜度的提高而增大的特性判斷飛行數(shù)據(jù)IMF的噪聲含量，據(jù)此確定界限N，優(yōu)化基于CEEMDAN的飛行數(shù)據(jù)濾波，消除其中的不確定性。其次，借助閾值法對界限N處的IMF濾波，解決常規(guī)CEEMDAN濾波后的飛行數(shù)據(jù)仍存有部分干擾噪聲的問題。

2 基于CEEMDAN和散布熵的飛行數(shù)據(jù)濾波

2.1 基于散布熵的IMF分量噪聲判斷

CEEMDAN去噪通常需要指標判斷IMF的噪聲含量，以避免盲目剔除IMF分量造成信號失真的問題，目前常用的指標有排列熵[8]、樣本熵等。本文選取散布熵作為評判指標，相對于樣本熵和排列熵，它具有計算速度快且兼顧時間序列幅值間關系的優(yōu)點。

散布熵的值會隨著時間序列復雜度的提高而增大，據(jù)此特性可以分析經(jīng)CEEMDAN分解獲得的飛行數(shù)據(jù)IMF分量。若飛行數(shù)據(jù)的IMF分量的散布熵值較大，則說明該IMF分量越不規(guī)則，復雜性越強，包含較多的噪聲信號；反之，則說明該IMF是原時間序列信號中的有用信息。因此，當相鄰IMF的散布熵值發(fā)生突變時，說明兩個IMF分量的噪聲含量急劇減少或增加，故本文以相鄰IMF分量散布熵比值最大處，作為界限判斷含噪IMF，克服了CEEMDAN去噪需要主觀確定界限N的困難。

散布熵的求解過程如下所示：

步驟1 給定長度為M的時間序列x={xi,i=1,2,…,M}，首先利用標準正態(tài)累積分布函數(shù)(normal cumulative distribution function，NCDF)將其映射為y={yi,i=1,2,…,M}

(12)

式中：σ和μ分別為時間序列x的標準差和期望，且yi∈(0,1)，然后將y中每個元素yi利用線性變換即式(13)映射為c個類別(1～c)，組成新的分類時間序列zc

(13)

步驟2 對序列zc進行相空間重構(gòu)，設嵌入維數(shù)取值為m，延遲時間取值為τ，則嵌入向量可表示為

(14)

步驟3 對每種散布模式，計算其概率p(πv0,v1,…,vm-1)，則按照Shannon信息熵的定義，散布熵的表達式定義為

(15)

De即為散布熵，常作為評判時間序列x不規(guī)則程度的指標。

2.2 閾值法去除IMF分量噪聲

對于界限N處的飛行數(shù)據(jù)IMF，其散布熵相對較小，故IMF中有用信號占主導。因而可以使用閾值法去除臨界IMF包含的少量干擾噪聲。閾值法的去噪思想是噪聲在飛行數(shù)據(jù)IMF分量中的系數(shù)幅值小于有效信號的系數(shù)幅值。通過選取適當?shù)拈撝?，判定大于該閾值的IMF系數(shù)幅值為有用信號，保持其值不變；對于低于該閾值的系數(shù)幅值判定為噪聲信號，將其減小至零，從而得到較為純凈的飛行數(shù)據(jù)IMF分量。

閾值法去除IMF分量噪聲的核心在于閾值和閾值函數(shù)的選取，為避免將有用信號誤當噪聲信號去除的問題，此處采取一個相對保守的策略，即基于Stein的無偏似然估計閾值法[11]計算待處理IMF分量的閾值。

對待處理的IMF分量Ii(t)，先對其每個采樣點的IMF系數(shù)幅值求絕對值的平方，其值按照由小到大的順序構(gòu)成序列 {w1,w2,…,wL}，L為Ii(t)序列的最大長度。若取 {w1,w2,…,wL} 第l個成員的平方根作為閾值，則該值作為閾值所造成的風險估計Risk(l)，如式(16)所示

(16)

假若令風險函數(shù)Risk(l)取得最小值的 {w1,w2,…,wL} 序列元素索引為lmin，那么最終確定的第i個待處理IMF信號的閾值λi，其表達式為

(17)

閾值函數(shù)的選取一般采取兩種形式，分別為硬閾值函數(shù)和軟閾值函數(shù)[12]。若選用硬閾值函數(shù)處理臨界IMF，會導致IMF出現(xiàn)不連續(xù)的間斷點；軟閾值函數(shù)雖然克服了硬閾值函數(shù)導致IMF分量產(chǎn)生間斷點的問題，但其處理后的IMF分量與原IMF之間存在較大誤差，會造成濾波后的飛行數(shù)據(jù)失真。目前，對閾值函數(shù)的改進較為成熟，有多種表達形式可供使用，且克服了上述兩種形式存在的問題，實際應用效果良好。采用文獻[12]得出的改進閾值函數(shù)，對IMF分量進行處理。改進的閾值函數(shù)的表達式為

(18)

式中：T=λe(1-α)(λ-|ωj,k|)，u=1-e-α(|ωj,k|-λ)2，0≤α≤1。

2.3 算法實現(xiàn)

首先從譯碼后的飛行數(shù)據(jù)中選取待處理數(shù)據(jù)，通常將其保存成CSV格式。程序?qū)崿F(xiàn)CEEMDAN算法，并設置輔助噪聲的標準差，最大篩選迭代次數(shù)，進而利用CEEMDAN算法對待處理飛行數(shù)據(jù)分解，得到其眾多IMF及余項。

其次，求解飛行數(shù)據(jù)各IMF的DE值即式(15)以及相鄰IMF分量的DE比值(取前項與后項的比值)，取DE比值最大處的IMF階次作為界限N。界限N之前的飛行數(shù)據(jù)IMF分量，其中噪聲信號占比較高，故直接將其去除；對包含少量噪聲的IMF分量，采取式(16)和式(17)計算其閾值，進而使用改進閾值函數(shù)即式(18)對該IMF分量處理，去除IMF中系數(shù)幅值較小的數(shù)據(jù)，保留飛行數(shù)據(jù)IMF分量中的有用信號。

最后，通過對有用信號分量、余項以及閾值法處理后的IMF累加，以獲得濾除干擾信號的飛行數(shù)據(jù)。算法流程如圖1所示。

圖1 CEEMDAN-DE飛行數(shù)據(jù)濾波算法流程

3 實驗與分析

3.1 仿真數(shù)據(jù)去噪

為檢驗所提飛行數(shù)據(jù)濾波方法的實際應用效果，采用Matlab提供的標準bumps信號作為原始信號，其波形如圖2(a)所示。向標準bumps信號添加白噪聲作為仿真飛行數(shù)據(jù)，其波形如圖2(b)所示，總采樣點數(shù)為2048個，信噪比設置為9 dB。

圖2 原始信號與含噪信號

從圖2(b)可以明顯看出該bumps信號的輪廓完全被噪聲掩蓋，包含眾多野值及劇烈抖動的數(shù)據(jù)。此外，采樣后該信號與飛行數(shù)據(jù)同樣為非線性離散時間序列，故可以作為驗證數(shù)據(jù)。利用CEEMDAN算法分解含噪bumps信號，可獲得其10個IMF及最后的余項，并計算全部IMF的散布熵。圖3直觀地展示了隨CEEMDAN分解次數(shù)的增加，各階次IMF散布熵的變動情況。圖4展示了信噪比為9 dB的bumps經(jīng)CEEMDAN算法分解得到的全部IMF的波形。

圖3 各階IMF分量的散布熵

圖4 仿真數(shù)據(jù)經(jīng)CEEMDAN分解后的IMF波形

通過計算，得到IMF1～IMF10的散布熵值為：[3.395,3.546,3.396,3.012,2.058,1.938,1.969,1.888,1.870,1.711]，以及低階IMF分量散布熵與相鄰高階IMF分量散布熵的比值為：[0.957,1.044,1.127,1.464,1.062,0.984,1.043,1.010,1.093]。IMF4與IMF5的散布熵比值最大，比值為1.464，說明在IMF4與IMF5分量之間噪聲含量驟減而有用信號增多，如圖3所示。因此，界限N取值為5，可以認為IMF1～IMF4分量中的有效信號完全被噪聲掩蓋，重構(gòu)信號時將其直接去除；IMF5～IMF10分量的散布熵值較小，故在這些分量中有用信號的系數(shù)幅值占主導。但IMF5分量處在有效信號占主導的臨界處，其包含有效信號的同時，不可避免地包含少量干擾信號，有待進一步處理。

通過觀察圖4(a)各階IMF分量的波形可以進一步驗證散布熵判斷的準確性。IMF1～IMF4分量的波形紊亂，存在劇烈震蕩，為噪聲信號；IMF5分量中有用信號的波形較為清晰，但系數(shù)幅值在等于零處存在頻繁地小幅度波動，且明顯小于有用信號的系數(shù)幅值。因此，可以利用閾值法處理IMF的方式，將IMF5中的部分干擾噪聲濾除。

為對比本文算法與EMD去噪、小波閾值去噪(小波基選為sym8小波，分解層數(shù)為5)以及傳統(tǒng)CEEMDAN去噪算法的降噪效果，采用信噪比(signal-to-noise ratio，SNR)和均方根誤差(root mean square error，RMSE)作為衡量指標，其結(jié)果見表1。

表1 各算法的SNR及RMSE

根據(jù)表1中各算法去噪后數(shù)據(jù)的信噪比和均方根誤差，可以表明文本去噪算法在提高數(shù)據(jù)信噪比和保持與原始數(shù)據(jù)相似性方面，效果最好。

上述各種算法的實際應用效果，如圖5所示，圖中小框內(nèi)的波形為信號第3個尖峰處的放大圖?？梢钥闯觯珽MD去噪算法對仿真數(shù)據(jù)的濾波效果不夠理想，依舊含有較多頻繁震蕩的數(shù)據(jù)；經(jīng)小波閾值去噪算法濾波后的數(shù)據(jù)雖然較為純凈，但在數(shù)據(jù)的第3個尖峰處過于平滑，丟失原bumps信號的部分細節(jié)信息，導致信號失真。傳統(tǒng)CEEMDAN去噪算法雖然效果優(yōu)于EMD去噪，但是未考慮臨界IMF分量中包含的少部分噪聲，導致去噪不徹底且波形平滑性較差。文中算法濾除噪聲的結(jié)果，如圖5(d)所示，其與原bumps信號波形的相似性最好，有效保持了在尖峰處的細節(jié)信息且波形相對平滑，無頻繁震蕩的數(shù)據(jù)。

圖5 去噪效果對比

3.2 實際飛行數(shù)據(jù)去噪

用于實驗的飛行數(shù)據(jù)來源于某航空公司譯碼后的QAR數(shù)據(jù)，本文選取其中某航班一次完整飛行過程中記錄的發(fā)動機N1數(shù)據(jù)進行去噪，以驗證本文算法對實際飛行數(shù)據(jù)的濾波效果。所選取的飛行數(shù)據(jù)長約8000 s，真實記錄了民航飛機在起飛、爬升、巡航等各個飛行階段的發(fā)動機N1參數(shù)。發(fā)動機N1參數(shù)由渦扇發(fā)動機風扇的抬輪速度提供，能夠表征發(fā)動機的推力大小，是用于民航發(fā)動機控制和監(jiān)控的關鍵參數(shù)之一。實際記錄的發(fā)動機N1數(shù)據(jù)如圖6所示。

圖6 實際飛行數(shù)據(jù)

如圖6所示，實際發(fā)動機N1數(shù)據(jù)在飛機巡航階段的變化較為穩(wěn)定，但仍受到氣流沖擊等隨機干擾的影響，導致存在小幅度頻繁震蕩；而在下降和進近階段，為了實現(xiàn)準確的飛行速度控制，發(fā)動機的推力會較為頻繁地改變。因此，發(fā)動機N1參數(shù)在這個階段的波動較為劇烈，更容易受到各種干擾噪聲的影響。故QAR實際記錄的N1數(shù)據(jù)中會出現(xiàn)較多的野值點以及局部頻繁震蕩的數(shù)據(jù)；另外，數(shù)據(jù)中存在較多離散的窄譜峰。利用EMD算法、小波閾值以及本文算法分別對QAR記錄的發(fā)動機N1數(shù)據(jù)進行濾波，其結(jié)果如圖7所示，圖7中小框內(nèi)的波形為信號最后一個尖峰處的放大圖。

圖7 實際飛行數(shù)據(jù)去噪效果對比

在使用EMD對N1數(shù)據(jù)分解時，共獲得15個IMF分量，通過多次實驗驗證，確定去除前5個被噪聲掩蓋的IMF可獲得最佳濾波效果。從圖7(a)所示波形得知，基于EMD的數(shù)據(jù)濾波雖有效剔除了N1數(shù)據(jù)中的野值，但因EMD存在模態(tài)混疊的問題，在原野值點處產(chǎn)生了震蕩數(shù)據(jù)，導致重構(gòu)后的N1數(shù)據(jù)波形不夠平滑。

圖7(b)為使用小波閾值法對發(fā)動機N1數(shù)據(jù)的濾波效果。在濾波過程中，為獲得較好的效果分別使用MATLAB提供的db系列以及sym系列小波基多次進行實驗，同時需要考慮分解層數(shù)的影響，濾波步驟相對復雜。最終，選取sym8小波作為小波基，分解層數(shù)設置為4，其閾值函數(shù)如式(18)所示，成功濾除了飛行數(shù)據(jù)中的野值。但是，對N1數(shù)據(jù)中眾多局部小幅度頻繁波動的數(shù)據(jù)，其濾除效果相對較差。

使用CEEMDAN算法對該發(fā)動機N1數(shù)據(jù)自適應分解，共獲得13個IMF分量，借助散布熵可以快速辨別出前4個IMF為被噪聲掩蓋的IMF，將其去除。其次，使用閾值法對階次為5的IMF進一步濾波，去除IMF系數(shù)幅值較小的數(shù)據(jù)，然后重構(gòu)發(fā)動機N1數(shù)據(jù)，其波形如圖7(c)所示。濾波結(jié)果能夠較為真實地保持發(fā)動機N1數(shù)據(jù)在不同飛行階段的變化特征，而且數(shù)據(jù)波形平滑、清晰，無毛刺、野值以及局部頻繁波動的數(shù)據(jù)。

實驗結(jié)果表明，本文濾波方法對實際飛行數(shù)據(jù)同樣具有良好的濾波效果，且不受飛行階段的影響。相比其它兩種方式，其在濾波步驟上得到簡化，且減少了局部震蕩數(shù)據(jù)的存在。

4 結(jié)束語

為簡便、有效地去除飛行數(shù)據(jù)中的高頻干擾，提高后續(xù)分析結(jié)果的準確度，本文綜合散布熵的特性以及CEEMDAN能自適應分解非線性時間序列信號的優(yōu)點，提出一種飛行數(shù)據(jù)濾波方法。在缺乏飛行數(shù)據(jù)先驗知識的狀況下，使用CEEMDAN算法，基于數(shù)據(jù)自身的時間尺度實現(xiàn)了對飛行數(shù)據(jù)的自適應分解；借助散布熵，簡單高效地識別出了包含噪聲的飛行數(shù)據(jù)IMF，無需大量實驗。從該方法實際應用于發(fā)動機N1數(shù)據(jù)的濾波效果可知，與QAR實際記錄的N1相比，濾波后的N1數(shù)據(jù)波形平滑，無野值，有效地降低了高頻干擾對飛行數(shù)據(jù)的影響。