石寶蘭 葉振信 楊小龍 付維賢

北京宇航系統(tǒng)工程研究所，北京 100076

固體導(dǎo)彈通常采用耗盡關(guān)機方式，當(dāng)下面級發(fā)動機燃燒室壓力下降到某一預(yù)裝訂數(shù)值時發(fā)出上面級發(fā)動機點火指令。考慮采用發(fā)動機后效段加表測得的視加速度數(shù)據(jù)作為上面級發(fā)動機點火判據(jù)，因此作為判據(jù)的視加速度數(shù)據(jù)的準(zhǔn)確度將直接影響級間分離發(fā)生的時刻，對下面級殘余推力、分離的相對速度、上面級起控均會產(chǎn)生較大影響。

導(dǎo)彈飛行過程中，視加速度由慣組中的加表測得。但原始測量數(shù)據(jù)跳動比較劇烈，直接作為判據(jù)使用將帶來很大的時間偏差，給級間分離帶來不利影響。因此必須對視加速度原始數(shù)據(jù)進行濾波，以減小數(shù)據(jù)跳動，提高級間分離時刻的準(zhǔn)確性。本文將從濾波平滑度、實時性和計算量3方面對4種濾波方法：動窗平均法、動窗擬合法、頻域濾波法及數(shù)字濾波法進行對比分析，摸索各自規(guī)律，以找到一種適合工程應(yīng)用的簡單有效的濾波方法。

1 問題的描述

某型導(dǎo)彈發(fā)動機后效段某秒內(nèi)視加速度原始數(shù)據(jù)經(jīng)歸一化處理后如圖1所示，圖中虛線為所有數(shù)據(jù)5階多項式擬合結(jié)果。顯然視加速度跳動劇烈，跳動幅值最高達0.3m/s2，無法直接作為判據(jù)使用，必須濾波。

圖1 視加速度原始數(shù)據(jù)

進行濾波前，先來分析導(dǎo)致視加速度數(shù)據(jù)跳動的原因，主要有2個：

1)發(fā)動機后效段燒蝕不穩(wěn)定

發(fā)動機推力下降段藥柱已基本消耗完全，主要燒蝕內(nèi)壁上的隔熱層，其燒蝕過程不穩(wěn)定，產(chǎn)生的推力本身就是跳動的。易證明推力與視加速度之間存在如下對應(yīng)關(guān)系[1]

2)加表存在測量噪聲

導(dǎo)彈飛行過程中存在風(fēng)干擾、結(jié)構(gòu)干擾等諸多不確定性外部干擾，加表無法辨識這些隨機干擾，因而存在系統(tǒng)噪聲。同時，加表內(nèi)部也存在多種電噪聲。

2 四種濾波方法

2.1 動窗平均法

動窗平均法是最簡單的數(shù)據(jù)濾波方法。其基本思想是以當(dāng)前時刻為基準(zhǔn)點，往前逐項獲取1組歷史數(shù)據(jù)，計算這組數(shù)據(jù)的代數(shù)平均值作為當(dāng)前時刻的濾波結(jié)果。到下一時刻則去掉最早的采樣點，補充當(dāng)前采樣點，即數(shù)據(jù)窗口隨時間往后移動，窗口內(nèi)采樣點數(shù)量固定不變[2]。

動窗平均法不可避免地會帶來時延問題。比如以2N+1個周期為1組數(shù)據(jù)，其代數(shù)平均值應(yīng)為采樣點中間時刻的近似值，以它為當(dāng)前時刻的濾波結(jié)果，必然會帶來N個周期的時延。圖2展示了窗口周期數(shù)為2，10，20時的濾波結(jié)果。

圖2 動窗平均法濾波結(jié)果局部放大圖

從上圖容易看出窗口周期數(shù)越大，時延越大。此外，還可以看出其濾波平滑度與窗口周期數(shù)亦正相關(guān)，即窗口周期數(shù)越大越能體現(xiàn)平均的濾波效果。工程應(yīng)用時必需根據(jù)具體數(shù)據(jù)合理選擇窗口周期數(shù)，周期數(shù)太大則時延導(dǎo)致的偏差大，周期數(shù)太小則濾波平滑度不好。對于這組數(shù)據(jù)，仿真結(jié)果表明：10周期平均的濾波結(jié)果較好，但濾波結(jié)果的跳動仍然不夠小，存在跳過預(yù)裝訂判據(jù)值的風(fēng)險。

2.2 動窗擬合法

動窗擬合法中的動窗概念與動窗平均法一樣，所不同的是對動窗內(nèi)數(shù)據(jù)的處理方法為基于最小二乘估計的多項式擬合法。利用Matlab提供的曲線擬合工具箱可以方便地實現(xiàn)多項式擬合，根據(jù)擬合多項式計算當(dāng)前時刻的擬合值作為濾波值。這樣的濾波結(jié)果是實時的，不存在時延。

2.2.1 設(shè)計參數(shù)

動窗擬合法的設(shè)計參數(shù)有2個：擬合階次與窗口數(shù)據(jù)量，需要根據(jù)原始數(shù)據(jù)特性仿真確定。

同樣其濾波平滑度與窗口數(shù)據(jù)量有關(guān)，而窗口數(shù)據(jù)量又決定了計算量。圖3是擬合階次固定為3，窗口數(shù)據(jù)量為40，60，80，100時的濾波結(jié)果。

圖3 動窗擬合法不同窗口數(shù)據(jù)量的濾波結(jié)果

仿真結(jié)果表明在相同的擬合階次下，窗口數(shù)據(jù)量低于80時濾波平滑度隨數(shù)據(jù)量增加而提高，高于80時平滑度隨數(shù)據(jù)量增加而提高的程度不甚明顯。因此，對于這組數(shù)據(jù)，窗口數(shù)據(jù)量選擇80比較合理。

現(xiàn)在固定窗口數(shù)據(jù)量為80，再來分析不同擬合階次對濾波結(jié)果的影響，如圖4所示。

如果取所有數(shù)據(jù)進行一次擬合，顯然階次越高，精度越好。但從濾波結(jié)果上來看階次越高、數(shù)據(jù)平滑度反而越差。這可能是因為多項式階次越高對系數(shù)的靈敏度越高。不過1階擬合的濾波結(jié)果偏差很大。從原始數(shù)據(jù)圖可以看出，視加速度下降趨勢近似于拋物線。因此，對于這組數(shù)據(jù)，擬合階次選擇2階比較合理。

2.2.2 預(yù)測擬合法

動窗擬合法的濾波結(jié)果較好，計算量也可接受，如果能進一步壓縮窗口數(shù)據(jù)量，將更有利于工程應(yīng)用。為此基于動窗擬合法提出預(yù)測擬合法，其基本思想是在動窗擬合法的基礎(chǔ)上利用當(dāng)前時刻的擬合多項式預(yù)測下一時刻的數(shù)據(jù)，當(dāng)數(shù)據(jù)更新后，取預(yù)測值與測量值的平均值作為當(dāng)前時刻的數(shù)據(jù)進行擬合，仍取擬合值為當(dāng)前時刻的濾波值。預(yù)測擬合法比動窗擬合法多了一個預(yù)測步，能夠利用歷史信息部分修正測量值，因此相同窗口數(shù)據(jù)量和擬合階次下，其濾波結(jié)果要優(yōu)于動窗擬合法。仿真結(jié)果也說明了這一點。

圖5 預(yù)測擬合法與其它方法濾波結(jié)果對比

圖5中的第1幅圖是擬合階次均為2階、窗口數(shù)據(jù)量均為40的預(yù)測擬合法和動窗擬合法的濾波結(jié)果對比，可以看出預(yù)測擬合法的濾波平滑度更好。第2幅圖是2階、40窗口數(shù)據(jù)量的預(yù)測擬合法和10周期平均法的濾波結(jié)果對比，可知當(dāng)窗口數(shù)據(jù)量取40時預(yù)測擬合法的濾波結(jié)果接近于10周期平均法。這說明預(yù)測擬合法能夠起到壓縮窗口數(shù)據(jù)量的作用。

2.3 頻域濾波法

以上濾波方法均為時域方法。加表得到的視加速度數(shù)據(jù)是時間的顯式函數(shù)，在時域進行處理是自然的。但是外部干擾與內(nèi)部噪聲也是以時間的顯式函數(shù)形式疊加在視加速度信號上，要想在時域?qū)⑺鼈儚氐追珠_是十分困難的。而在導(dǎo)彈飛行力學(xué)環(huán)境下，視加速度信號的頻譜與干擾信號或噪聲的頻譜一般是不同的。也就是說，在頻域可以實現(xiàn)視加速度信號與干擾信號或噪聲的分離，達到濾波目的。這樣就產(chǎn)生了頻域濾波的思想[3]。作為溝通時域和頻域的有效工具，傅立葉變換(FT)和反變換(IFT)使得實現(xiàn)頻域濾波的思想成為可能。但是離散傅立葉變換的計算量非常巨大，限制了它的工程應(yīng)用。在數(shù)字信號處理領(lǐng)域中廣泛使用的是一種被稱為快速傅立葉變換(FFT)的離散傅立葉變換的快速算法，它大大減少了運算量。

頻域濾波的基本過程包括3步，如圖6所示：

1)利用FFT將被噪聲污染的測量信號轉(zhuǎn)換到頻域；

2)頻譜估計，即：在頻域中確定有用信號和各類噪聲占據(jù)的頻段；

3)利用IFFT將有用信號轉(zhuǎn)換回時域。

圖6 頻域濾波的基本過程

以窗口周期數(shù)為5為例，頻域濾波結(jié)果如圖7所示。仿真結(jié)果表明，頻域濾波結(jié)果的平滑度和實時性均非常好。

圖7 頻域濾波結(jié)果

2.4 數(shù)字濾波法

數(shù)字濾波法的核心是數(shù)字濾波器的設(shè)計。數(shù)字濾波器是數(shù)字信號處理的一個重要技術(shù)分支，其實質(zhì)是一種用來描述離散系統(tǒng)輸入與輸出關(guān)系差分方程的計算或卷積計算[4]。數(shù)字濾波器的設(shè)計就是根據(jù)要求選擇系統(tǒng)h(n)或H(z)，當(dāng)原始數(shù)據(jù)通過系統(tǒng)時，對其波形和頻譜進行加工，獲得人們需要的信號。

按單位沖激響應(yīng)的樣值個數(shù)是否有限，數(shù)字濾波器可分為有限沖激響應(yīng)(FIR)和無限沖激響應(yīng)(IIR)2類；按實現(xiàn)方法和形式不同，數(shù)字濾波器可分為遞歸型、非遞歸型和快速卷積型3類。本文的目的是尋找一種適合工程應(yīng)用的濾波方法，所以選擇的是FIR數(shù)字濾波器，采用非遞歸法實現(xiàn)，即輸出僅與輸入有關(guān)，與歷史輸出無關(guān)。這樣得到的數(shù)字濾波器非常簡單，計算量小。

利用Matlab提供的濾波器設(shè)計箱進行FIR數(shù)字濾波器的設(shè)計。其系統(tǒng)函數(shù)為

即濾波值是當(dāng)前數(shù)據(jù)和前N個測量數(shù)據(jù)的加權(quán)平均值，因此也存在時延。容易看出當(dāng)取h(n)=1/(N+1)時，即退化為動窗平均法。

沿用動窗的思想，本文采用窗函數(shù)設(shè)計FIR數(shù)字濾波器，需要設(shè)計的指標(biāo)有：窗函數(shù)類型、濾波器階數(shù)和通帶截止頻率ωc(Hz)。

與動窗平均法一樣，數(shù)字濾波器階數(shù)越高，時延越大。另一方面，濾波器階數(shù)太低會影響濾波效果。經(jīng)多次仿真確定最佳濾波器階數(shù)為10階。接下來設(shè)計ωc。ωc的選擇應(yīng)該保證能夠區(qū)分視加速度信號和噪聲信號，這里有頻域濾波的影子。從5階多項式擬合結(jié)果可以看出視加速度信號的趨勢項近似為拋物線，頻率甚低，而噪聲信號頻率比較高，兩者容易區(qū)分。

圖8 不同ωc的數(shù)字濾波結(jié)果

圖8展示的是ωc取50，30，10，5和1時的10階Hamming窗FIR數(shù)字濾波器濾波結(jié)果?？梢钥闯靓豤越小濾波平滑度越好，減小到10以內(nèi)時，濾波平滑度已經(jīng)較好，幾乎不再提高。因此這里ωc取1。

當(dāng)采用Hamming窗時，濾波器階數(shù)取10，ωc取1的濾波效果較好。下面再來看不同窗函數(shù)類型的濾波結(jié)果。除了Hamming窗，常用的窗函數(shù)還有2種：Hann窗和Blackman窗。其它設(shè)計參數(shù)相同的情況下，Hamming窗數(shù)字濾波器的濾波效果略好于其余兩種，仿真結(jié)果的局部放大圖(見圖9)說明了這一點。

圖9 不同窗函數(shù)的濾波結(jié)果局部放大圖

綜上所述，最終設(shè)計的濾波器為10階Hamming窗FIR數(shù)字濾波器，ωc取1，濾波結(jié)果與11周期動窗平均法的結(jié)果對比如圖10所示。可以看出，在利用同樣多歷史數(shù)據(jù)的情況下，數(shù)字濾波法的濾波平滑度明顯優(yōu)于動窗平均法。

圖10 數(shù)字濾波法與動窗平均法的濾波結(jié)果對比

3 四種濾波方法的對比

前面通過對每種方法單獨分析，確定了各自的最優(yōu)設(shè)計參數(shù)。下面對4種方法的濾波結(jié)果進行對比并總結(jié)各自的特點。11周期動窗平均法、2階40周期預(yù)測擬合法、頻域濾波法和10階Hamming窗數(shù)字濾波法的仿真結(jié)果如圖11所示。

圖11 四種方法的濾波結(jié)果對比

從圖11可以看出：

1)頻域濾波結(jié)果的平滑度最好，數(shù)字濾波次之，動窗平均和預(yù)測擬合相對略差；

2)動窗平均法和數(shù)字濾波法的濾波結(jié)果存在時延，其余方法均為實時濾波；

3)頻域濾波結(jié)果在某些時段偏離5階多項式擬合結(jié)果較大，且偏差量難以進行補償。

工程應(yīng)用時，受彈上計算機運算速度和內(nèi)存容量限制，濾波方法的計算量不能太大，并且要求實時性好。因此選擇濾波方法時要綜合考慮在濾波平滑度、實時性和計算量等方面的表現(xiàn)。4種濾波方法的特點總結(jié)見表1。

表1 4種濾波方法的特點匯總表

需要特別指出的是，頻域濾波的計算量特別巨大，現(xiàn)有彈上計算機的運算速度和內(nèi)存容量暫時無法滿足其要求，在導(dǎo)彈飛行控制中還未見應(yīng)用實例。從綜合表現(xiàn)來看，數(shù)字濾波法簡單而有效，比動窗平均法更加適合于工程應(yīng)用。

4 數(shù)字濾波法的工程應(yīng)用

數(shù)字濾波法的計算量非常小，濾波平滑度比較

好，能夠很好地滿足工程應(yīng)用要求。但也存在不可避免的時延問題，且視加速度數(shù)據(jù)變化越快，時延造成的數(shù)據(jù)偏差越大。工程應(yīng)用時必須對時延造成的數(shù)據(jù)偏差進行補償，以提高級間分離時刻的準(zhǔn)確性。

級間分離判據(jù)諸元必須預(yù)先裝訂好。本質(zhì)上是以發(fā)動機燃燒室壓強為判據(jù)，根據(jù)推力與視加速度之間的關(guān)系可以計算出發(fā)出點火指令時刻對應(yīng)的視加速度理論值。考慮濾波方法帶來的時延，還要加上數(shù)據(jù)偏差才能作為視加速度判據(jù)諸元的裝訂值。以10階Hamming窗數(shù)字濾波法為例，其固有的時延為5周期。理論彈道計算可以得到視加速度下降速率，其與時延的乘積作為數(shù)據(jù)偏差，與視加速度判據(jù)理論值相加即可得到視加速度判據(jù)諸元裝訂值。

5 結(jié)論

本文從提高級間分離時刻的準(zhǔn)確性這一工程背景出發(fā)對發(fā)動機后效段視加速度數(shù)據(jù)濾波方法進行了研究，通過對動窗平均法、動窗擬合法、頻域濾波法和數(shù)字濾波法4種濾波方法的對比分析，得到如下結(jié)論：

1)數(shù)字濾波法的計算量非常小，濾波平滑度比較好，比動窗平均法更加適合于工程應(yīng)用；

2)數(shù)字濾波法的工程應(yīng)用必須對時延帶來的數(shù)據(jù)偏差進行補償，這種補償體現(xiàn)在視加速度判據(jù)諸元裝訂值上；

3)動窗擬合法濾波效果相對略差，計算量稍大，但實時性好，工程上可以考慮使用；

4)頻域濾波法濾波平滑度非常好，但計算量太大，工程上尚難以應(yīng)用。

[1] 王錚，胡永強．固體火箭發(fā)動機[M]．北京：中國宇航出版社，1993．

[2] 鄒洪，向大威，景永剛．多普勒計程儀的數(shù)據(jù)平滑方法[J]．聲學(xué)技術(shù)，2008，27(4)：507-510．(ZOU Hong，XIANG Dawei，JING Yonggang．Data Smoothing Methods for DVL[J]．Technical Acoustics，2008，27(4)：507-510．)

[3] 劉保中．隨機噪聲的頻域濾波方法[J]．戰(zhàn)術(shù)導(dǎo)彈控制技術(shù)，2008(1)：15-17．(LIU Baozhong．Filtering Method Based on Frequency Domain of Random Noise[J]．Control Technology of Tactical Missile，2008(1)：15-17．)

[4] 叢玉良，王宏志．?dāng)?shù)字信號處理原理及其MATLAB實現(xiàn)[M]．北京：電子工業(yè)出版社，2005．