發動機后效段加表數據濾波方法研究

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

北京宇航系統工程研究所，北京 100076

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

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

1 問題的描述

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

圖1 視加速度原始數據

進行濾波前，先來分析導致視加速度數據跳動的原因，主要有2個：

1)發動機后效段燒蝕不穩定

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

2)加表存在測量噪聲

導彈飛行過程中存在風干擾、結構干擾等諸多不確定性外部干擾，加表無法辨識這些隨機干擾，因而存在系統噪聲。同時，加表內部也存在多種電噪聲。

2 四種濾波方法

2.1 動窗平均法

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

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

圖2 動窗平均法濾波結果局部放大圖

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

2.2 動窗擬合法

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

2.2.1 設計參數

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

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

圖3 動窗擬合法不同窗口數據量的濾波結果

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

現在固定窗口數據量為80，再來分析不同擬合階次對濾波結果的影響，如圖4所示。

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

2.2.2 預測擬合法

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

圖5 預測擬合法與其它方法濾波結果對比

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

2.3 頻域濾波法

以上濾波方法均為時域方法。加表得到的視加速度數據是時間的顯式函數，在時域進行處理是自然的。但是外部干擾與內部噪聲也是以時間的顯式函數形式疊加在視加速度信號上，要想在時域將它們徹底分開是十分困難的。而在導彈飛行力學環境下，視加速度信號的頻譜與干擾信號或噪聲的頻譜一般是不同的。也就是說，在頻域可以實現視加速度信號與干擾信號或噪聲的分離，達到濾波目的。這樣就產生了頻域濾波的思想[3]。作為溝通時域和頻域的有效工具，傅立葉變換(FT)和反變換(IFT)使得實現頻域濾波的思想成為可能。但是離散傅立葉變換的計算量非常巨大，限制了它的工程應用。在數字信號處理領域中廣泛使用的是一種被稱為快速傅立葉變換(FFT)的離散傅立葉變換的快速算法，它大大減少了運算量。

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

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

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

3)利用IFFT將有用信號轉換回時域。

圖6 頻域濾波的基本過程

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

圖7 頻域濾波結果

2.4 數字濾波法

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

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

利用Matlab提供的濾波器設計箱進行FIR數字濾波器的設計。其系統函數為

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

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

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

圖8 不同ωc的數字濾波結果

圖8展示的是ωc取50，30，10，5和1時的10階Hamming窗FIR數字濾波器濾波結果。可以看出ωc越小濾波平滑度越好，減小到10以內時，濾波平滑度已經較好，幾乎不再提高。因此這里ωc取1。

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

圖9 不同窗函數的濾波結果局部放大圖

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

圖10 數字濾波法與動窗平均法的濾波結果對比

3 四種濾波方法的對比

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

圖11 四種方法的濾波結果對比

從圖11可以看出：

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

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

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

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

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

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

4 數字濾波法的工程應用

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

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

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

5 結論

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

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

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

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

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

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

[2] 鄒洪，向大威，景永剛．多普勒計程儀的數據平滑方法[J]．聲學技術，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]．戰術導彈控制技術，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] 叢玉良，王宏志．數字信號處理原理及其MATLAB實現[M]．北京：電子工業出版社，2005．