一種改進(jìn)濾波器與LSTM融合的飛輪故障預(yù)測(cè)算法

摘" 要：""""" 針對(duì)基于濾波器的預(yù)測(cè)方法由于依賴(lài)模型精度而導(dǎo)致長(zhǎng)周期預(yù)測(cè)精度低的問(wèn)題， 提出了一種基于濾波器和長(zhǎng)短期記憶（LSTM）網(wǎng)絡(luò)的融合故障預(yù)測(cè)算法， 實(shí)現(xiàn)航天器飛輪執(zhí)行器緩變故障的預(yù)測(cè)。 首先， 分別設(shè)計(jì)了小批量標(biāo)準(zhǔn)化LSTM網(wǎng)絡(luò)和趨勢(shì)識(shí)別模塊， 二者串聯(lián)組成神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器提高對(duì)時(shí)序預(yù)測(cè)的準(zhǔn)確性。 然后， 利用遞歸最小二乘（RLS）參數(shù)估計(jì)原理改進(jìn)卡爾曼濾波器更新過(guò)程，以增強(qiáng)對(duì)時(shí)序預(yù)測(cè)誤差的魯棒性。 在此基礎(chǔ)上， 將神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器輸出的預(yù)測(cè)值與改進(jìn)后的濾波器相融合， 獲得未來(lái)時(shí)刻的預(yù)測(cè)殘差項(xiàng)實(shí)現(xiàn)迭代更新和預(yù)測(cè)， 克服了濾波器算法對(duì)模型的依賴(lài)， 提高了預(yù)測(cè)精度。 最后， 設(shè)計(jì)仿真實(shí)驗(yàn)比較了三種神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器的時(shí)序預(yù)測(cè)性能， 并考慮飛輪軸承性能退化故障， 利用所提融合預(yù)測(cè)算法判斷出飛輪在856 s時(shí)性能退化達(dá)到閾值， 預(yù)測(cè)時(shí)間誤差為36 s， 從而驗(yàn)證了預(yù)測(cè)算法對(duì)緩變故障的有效性。

關(guān)鍵詞："""" 故障預(yù)測(cè)； LSTM； 卡爾曼濾波器； RLS； 飛輪

中圖分類(lèi)號(hào)："""""" TJ760;V44

文獻(xiàn)標(biāo)識(shí)碼：""" A

文章編號(hào)："""" 1673-5048（2024）03-0129-08

DOI： 10.12132/ISSN.1673-5048.2023.0211

引用格式： 龍弟之， 陳辛， 魏炳翌， 等 ." 一種改進(jìn)濾波器與LSTM融合的飛輪故障預(yù)測(cè)算法［ J］. 航空兵器， 2024， 31（ 3）： 129-136.

Long Dizhi， Chen Xin， Wei Bingyi， et al. Flywheel Fault Prediction Algorithm Based on Improved Filter and LSTM Fusion［ J］. Aero Weaponry， 2024， 31（ 3）： 129-136.（ in Chinese）

0" 引" 言

故障預(yù)測(cè)是在系統(tǒng)早期異常而尚未產(chǎn)生明顯變化時(shí)， 根據(jù)歷史數(shù)據(jù)來(lái)估計(jì)系統(tǒng)發(fā)展趨勢(shì)從而判斷下一時(shí)刻是否會(huì)發(fā)生故障。 有效的預(yù)測(cè)方法能夠輔助監(jiān)測(cè)人員及時(shí)預(yù)報(bào)在軌航天器可能的潛在故障， 將傳統(tǒng)“檢測(cè)+應(yīng)對(duì)”的故障排除思路提升為“預(yù)測(cè)+干預(yù)”， 以最小的代價(jià)阻止故障影響擴(kuò)大， 顯著提高了航天系統(tǒng)的安全裕度［1-3］。

目前， 研究較為廣泛的是以卡爾曼濾波器（Kalman Filter， KF）和粒子濾波器（Particle Filter， PF）為代表的基于模型的預(yù)測(cè)方法。 在系統(tǒng)數(shù)學(xué)模型可用情況下， 利用KF ［4］、 擴(kuò)展KF［5］、 無(wú)跡KF［6］和PF［7］對(duì)含有高斯/非高斯噪聲的測(cè)量信號(hào)進(jìn)行處理， 從而實(shí)現(xiàn)系統(tǒng)狀態(tài)或故障參數(shù)的預(yù)測(cè)。 然而， 在預(yù)測(cè)過(guò)程中由于濾波器無(wú)法獲取未來(lái)時(shí)刻的系統(tǒng)測(cè)量輸出， 其修正項(xiàng)為零， 導(dǎo)致無(wú)法更新系統(tǒng)狀態(tài)量和協(xié)方差矩陣， 在模型不精確的情況下隨著預(yù)測(cè)步長(zhǎng)的增加， 預(yù)測(cè)誤差將會(huì)顯著增大。 文獻(xiàn)［8］兼顧預(yù)測(cè)性能和運(yùn)算時(shí)效， 提出了一類(lèi)融合預(yù)測(cè)的思想， 利用數(shù)據(jù)驅(qū)動(dòng)的方法預(yù)測(cè)系統(tǒng)的變化趨勢(shì)， 然后與基于模型的方法相融合開(kāi)展預(yù)測(cè)。 但融合預(yù)測(cè)算法在使用數(shù)據(jù)驅(qū)動(dòng)方法時(shí)， 其預(yù)測(cè)誤差也會(huì)傳遞到濾波器更新過(guò)程中。 當(dāng)進(jìn)行中長(zhǎng)期預(yù)測(cè)任務(wù)時(shí)， 預(yù)測(cè)精度不夠理想， 因此， 提高基于數(shù)據(jù)驅(qū)動(dòng)算法的預(yù)測(cè)性能是融合算法的核心需求之一。

隨著系統(tǒng)復(fù)雜性的提升以及樣本數(shù)據(jù)的多樣化， 相比傳統(tǒng)淺層神經(jīng)網(wǎng)絡(luò)， 基于深度神經(jīng)網(wǎng)絡(luò)的方法表現(xiàn)出更明顯的優(yōu)勢(shì)。 LSTM網(wǎng)絡(luò)是一種由多個(gè)非線(xiàn)性映射層疊加而成的深度循環(huán)神經(jīng)網(wǎng)絡(luò)模型， 克服了隨網(wǎng)絡(luò)層數(shù)遞增產(chǎn)生的梯度消失和過(guò)擬合問(wèn)題， 能夠?qū)r(shí)序數(shù)據(jù)之間的依賴(lài)性進(jìn)行建模［9］。 相比于其他預(yù)測(cè)模型， LSTM網(wǎng)絡(luò)對(duì)于時(shí)序數(shù)據(jù)的預(yù)測(cè)性能更佳［10-11］。

飛輪作為航天器常用執(zhí)行器之一， 長(zhǎng)期在軌工作， 零部件將發(fā)生損耗導(dǎo)致性能衰退。 由于變化幅度小且難以獲得精確模型， 采用傳統(tǒng)濾波器算法預(yù)測(cè)精度低， 無(wú)法有效預(yù)警故障。 考慮在軌航天器擁有海量遙測(cè)數(shù)據(jù)， 可將基于濾波器的預(yù)測(cè)算法與LSTM網(wǎng)絡(luò)相融合， 既彌補(bǔ)基于數(shù)據(jù)方法由于對(duì)樣本代表性和數(shù)量的高要求而魯

收稿日期： 2023-11-02

基金項(xiàng)目： 技術(shù)基礎(chǔ)研究項(xiàng)目（514010203-103）

*作者簡(jiǎn)介： 龍弟之（1995-）， 男， 重慶人， 博士。

棒性差的缺點(diǎn)， 又解決了濾波器算法對(duì)模型依賴(lài)的問(wèn)題。 因此， 如何發(fā)揮LSTM網(wǎng)絡(luò)與濾波器各自在融合預(yù)測(cè)方法中的優(yōu)勢(shì)并應(yīng)用于飛輪故障預(yù)測(cè)， 有待進(jìn)一步研究。

1" 姿控系統(tǒng)輸出時(shí)序預(yù)測(cè)

1.1" 小批量標(biāo)準(zhǔn)化LSTM網(wǎng)絡(luò)

為解決由于LSTM網(wǎng)絡(luò)結(jié)構(gòu)復(fù)雜度增加而出現(xiàn)的協(xié)變量漂移現(xiàn)象， 采用小批量標(biāo)準(zhǔn)化對(duì)輸入進(jìn)行處理。 每組小批量訓(xùn)練集對(duì)應(yīng)產(chǎn)生一組均值和方差的估計(jì)， 使得用于標(biāo)準(zhǔn)化的統(tǒng)計(jì)信息能夠參與到梯度反向傳播中［12］。 假設(shè)輸入節(jié)點(diǎn)的一組小批量訓(xùn)練集為

航空兵器" 2024年第31卷第3期

龍弟之， 等： 一種改進(jìn)濾波器與LSTM融合的飛輪故障預(yù)測(cè)算法

Seta={a1， a2， …， am}， ai∈RL（1）

式中： m為一組訓(xùn)練集包含的輸入特征序列數(shù)量； L為每個(gè)輸入特征序列的維度。 對(duì)所有維度進(jìn)行并行標(biāo)準(zhǔn)化處理， 以第j維輸入aji（i=1， 2， …， m； j=1， 2， …， L）為例， 過(guò)程如下。

首先計(jì)算訓(xùn)練集平均值和方差：

μj=1m∑mi=1aji，" σ2j=1m∑mi=1（aji－μj）2（2）

得到標(biāo)準(zhǔn)化后的輸入為

a^ji=（aji－μj）/σ2j+εj（3）

式中： εj是為了數(shù)值穩(wěn)定性而添加到小批量方差中的標(biāo)準(zhǔn)化參數(shù)。

為確保每層輸入的標(biāo)準(zhǔn)化不改變?cè)搶拥挠成浔碚鲀?nèi)容， 引入?yún)?shù)ρj和ηj對(duì)輸入進(jìn)行縮放和平移， 即

a-ji=ρja^ji+ηj（4）

式中的ρj和ηj在學(xué)習(xí)過(guò)程中將會(huì)不斷調(diào)整。

至此， 原始輸入aji轉(zhuǎn)變?yōu)閍-ji， 該標(biāo)準(zhǔn)化過(guò)程記作BN（ρ， η， a）： a1， …， m→a-1， …， m。 將小批量標(biāo)準(zhǔn)化用于LSTM網(wǎng)絡(luò)中， 分別對(duì)LSTM細(xì)胞元的狀態(tài)和輸出進(jìn)行處理。 為了避免過(guò)擬合和不必要的冗余， 簡(jiǎn)化LSTM網(wǎng)絡(luò)， 設(shè)置偏置和η為零， 得到帶窺視功能的小批量標(biāo)準(zhǔn)化LSTM網(wǎng)絡(luò)前向傳播表達(dá)式為

［z-Ti， t" z-iTi， tz-fTi， tz-oTi， t］T =Wxxi， t +BN（ρh， Wh·hi， t-1）+

""""Wpci

［zTi， tziTi， tzfTi， tzoTi，" t］T=

［Th（z-i， t ）T σ（z-ii， t）T σ（z-fi， t）T σ（z-oi， t）T］T

ci， t =zi， t ⊙zii， t+ci，t-1 ⊙zfi， t

hi， t =Th（BN（ρc， ci， t））⊙zoi， t（5）

式中： 下標(biāo)i表示第i組訓(xùn)練樣本； xi， t， ci， t和hi， t分別為當(dāng)前時(shí)刻細(xì)胞元的輸入、 狀態(tài)和輸出； ci， t－1和hi， t－1為前一時(shí)刻細(xì)胞元的狀態(tài)和輸出。 數(shù)據(jù)集經(jīng)過(guò)歸一化處理， 所以無(wú)需對(duì)xi， t標(biāo)準(zhǔn)化。 zi， t， zii， t， zfi， t和zoi， t分別為細(xì)胞元門(mén)控輸入、 輸入門(mén)、 遺忘門(mén)和輸出門(mén)。 σ和Th為sigmoid函數(shù)和tanh函數(shù)。 ρh=［ρh， zρh， iρh， fρh， o］和ρc分別為對(duì)應(yīng)門(mén)控結(jié)構(gòu)和狀態(tài)的縮放參數(shù)， 各項(xiàng)權(quán)值分別為

Wx=［WzWiWfWo］T∈R4N×M;

Wh=［RzRiRfRo］T∈R4N×N;

Wpci=［0pi⊙ci， t－1pf⊙ci， t－1po⊙ci， t］T∈R4N。

其中， M和N分別為輸入向量的維數(shù)和隱藏層細(xì)胞元數(shù)； 符號(hào)⊙表示元素對(duì)應(yīng)相乘。 定義下標(biāo)（·）為集合{z" i" f" o}中任意元素， 則W（·）∈RN×M， R（·）∈RN×N和p（·）∈RN分別為各門(mén)結(jié)構(gòu)的輸入權(quán)值矩陣、 循環(huán)權(quán)值矩陣和窺視孔權(quán)值向量。

LSTM網(wǎng)絡(luò)將細(xì)胞元輸出hi， t用于分類(lèi)或回歸， 根據(jù)不同任務(wù)選取不同激活函數(shù)。 用于回歸任務(wù)時(shí)， LSTM網(wǎng)絡(luò)可進(jìn)行時(shí)序預(yù)測(cè)， 輸出為y^i=g（Vo·hi， t）， 其中g(shù)（·）為特定激活函數(shù)。 給定預(yù)測(cè)期望值yi， 定義誤差損失函數(shù)為

L（y， y^）=12m∑mi=1（yi－y^i）2（6）

由于LSTM網(wǎng)絡(luò)在設(shè)計(jì)上具有時(shí)間關(guān)聯(lián)性， 其誤差反向傳播除了沿網(wǎng)絡(luò)結(jié)構(gòu)的反向傳播外， 還包括沿時(shí)間的反向傳播。 首先， 根據(jù)損失函數(shù)， 計(jì)算t時(shí)刻細(xì)胞元輸出、 門(mén)控輸入和三個(gè)門(mén)控節(jié)點(diǎn)對(duì)應(yīng)的誤差項(xiàng)分別為

Δi， t =Lhi， t" =-VTog′（Vo·hi， t ）/m

Δoi， t=Lz-oi， t=Δi，" t ⊙Th（c^i， t ）⊙σ′（z-oi， t）

Δfi， t=Lz-fi， t=Δi， t ⊙zoi， t⊙Th′（c^i， t ）⊙ci， t-1 ⊙σ′（z-fi， t）

Δii， t=Lz-ii， t=Δi， t ⊙zoi， t⊙Th′（c^i， t ）⊙zi， t ⊙σ′（z-ii， t）

Δzi， t=Lz-i， t =Δi， t ⊙zoi， t⊙Th′（c^i， t ）⊙zii， t⊙Th′（z-i， t ）（7）

式中： c^i， t=BN（ρc， ci， t）為標(biāo)準(zhǔn)化后的細(xì)胞元狀態(tài)。

然后沿時(shí)間反向傳播， 計(jì)算t－1時(shí)刻的誤差項(xiàng)：

Δi， t－1=Lhi， t－1=Lhi， t hi， thi， t－1=Δi， thi， thi， t－1（8）

由LSTM細(xì)胞結(jié)構(gòu)可知， hi， t－1分別與門(mén)控輸入和三個(gè)門(mén)控節(jié)點(diǎn)相連， 利用全導(dǎo)數(shù)公式可得

Δi， t－1=Δi， thi， tzoi， t zoi， tz-oi， t z-oi，" thi， t－1+Δi， thi， tci， t ci， tzfi， t zfi， tz-fi， t·

z-fi， thi， t－1+

Δi， thi， tci， t ci， tzii， t zii， tz-ii， t z-ii， thi， t－1+

Δi， thi， tci， t ci， tzi， t zi， tz-i， t z-i， thi， t－1=

Δoi， tz-oi， thi， t－1+

Δfi， tz-fi， thi， t－1+Δii， tz-ii， thi， t－1+Δzi， tz-i， thi， t－1（9）

根據(jù)前向傳播式（5）和小批量標(biāo)準(zhǔn)化式（2）～（4）， 可得z-oi， t關(guān)于hi， t－1的偏導(dǎo)數(shù)為

z-oi， thi， t－1=z-oi， th^i， t－1·1σ2h， o+εh， o+z-oi， tσ2h， o·

2（hi， t－1－μh， o）m+z-oi， tμh， o·1m（10）

式中： μh， o， σ2h， o和εh， o分別對(duì)應(yīng)標(biāo)準(zhǔn)化過(guò)程中的平均值、 方差和標(biāo)準(zhǔn)化參數(shù)； h^i， t－1為標(biāo)準(zhǔn)化后的細(xì)胞元輸出。 式（10）中各項(xiàng)偏導(dǎo)數(shù)分別為

z-oi， th^i， t－1=z-oi， tBN（ρh， o， Ro·hi， t－1）ρh， o

z-oi， tσ2h， o=z-oi， th^i， t－1（hi， t－1－μh， o）·－12（σ2h， o+εh， o）－32

z-oi， tμh， o=z-oi， th^i， t－1 －1σ2h， o+εh， o （11）

同理， 可得z-fi， t， z-ii， t和z-i， t分別關(guān)于hi， t－1的偏導(dǎo)數(shù)以及Th（c^i， t）關(guān)于c^i， t的偏導(dǎo)數(shù)。 至此， 得到小批量標(biāo)準(zhǔn)化LSTM網(wǎng)絡(luò)預(yù)測(cè)誤差沿時(shí)間反向傳播t－1時(shí)刻的表達(dá)式， 那么， 向前傳播到任意Ta時(shí)刻的誤差項(xiàng)為

Δi， Ta=∏t－1l=TaΔoi， l+1z-oi， l+1hi， l+Δfi， l+1z-fi， l+1hi， l+Δii， l+1z-ii， l+1hi， l+Δzi， l+1z-i， l+1hi， l（12）

利用誤差項(xiàng)求出權(quán)值矩陣的梯度， 其中循環(huán)連接權(quán)值矩陣和窺視孔權(quán)值向量需要將各時(shí)刻梯度相加得到最終梯度。 那么， 各權(quán)值的反向傳播梯度為

ΔW（·）=〈Δ（·）i， t， xi， t〉

ΔR（·）=∑tl=1〈Δ（·）i， l， ρhh^i， l－1〉

Δpi=∑tl=1ci， l－1⊙Δii， l

Δpf=∑tl=1ci， l－1⊙Δfi， l

Δpo=∑tl=1ci， l－1⊙Δoi， lci， l （13）

式中： 〈·， ·〉表示兩個(gè)向量的外積。

縮放參數(shù)在訓(xùn)練中的更新梯度為

Lρc=－∑tl=1Δi， l⊙zoi， l·Th′（BN（ρc， ci， l））c^i， l

Lρh， （·）=－∑tl=1Δ（·）i， l·Th′（BN（ρh， （·）， R（·）hi， l－1））·

R（·）h^i， l－1 （14）

1.2" 神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器

為實(shí)現(xiàn)時(shí)序預(yù)測(cè)功能， 設(shè)計(jì)以小批量標(biāo)準(zhǔn)化LSTM網(wǎng)絡(luò)為核心的預(yù)測(cè)神經(jīng)網(wǎng)絡(luò)， 結(jié)構(gòu)如圖1所示， 包括輸入層、 LSTM層、 線(xiàn)性全連接層和預(yù)測(cè)輸出層。

對(duì)于緩變故障情況， 輸入神經(jīng)網(wǎng)絡(luò)一組時(shí)間序列的整體變化趨勢(shì)對(duì)判斷后續(xù)變化有一定輔助作用。 基于該

特性， 在時(shí)間序列進(jìn)入預(yù)測(cè)神經(jīng)網(wǎng)絡(luò)輸入層前， 通過(guò)趨勢(shì)識(shí)別模塊， 將輸入的時(shí)間序列分為上升趨勢(shì)、 近似平緩趨勢(shì)或下降趨勢(shì)。 每一類(lèi)趨勢(shì)的時(shí)間序列用于訓(xùn)練對(duì)應(yīng)的預(yù)測(cè)神經(jīng)網(wǎng)絡(luò)， 那么神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器整體結(jié)構(gòu)如圖2所示。

為減小計(jì)算量， 基于最小二次線(xiàn)性回歸法設(shè)計(jì)趨勢(shì)識(shí)別模塊。 假設(shè)輸入神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器的一組時(shí)序數(shù)據(jù)為x1， x2， …， xL， 轉(zhuǎn)化為坐標(biāo)軸上的點(diǎn)（t1， x1）， （t2， x2）， …， （tM， xL）。 用線(xiàn)性方程x=kst+bs近似擬合， 二次擬合誤差為

e（ks， bs）=∑Li=1（ksi+bs－xi）2（15）

根據(jù)微積分求極值思想， 應(yīng)滿(mǎn)足

e（ks， bs）ks=0，nbsp; e（ks， bs）bs=0（16）

求解線(xiàn)性方程的斜率ks用于輸入時(shí)序數(shù)據(jù)的趨勢(shì)識(shí)別。 當(dāng)－π/2≤arctanks≤τd， 趨勢(shì)識(shí)別為下降； 當(dāng)τa≤arctanks≤π/2， 趨勢(shì)識(shí)別為上升； 當(dāng)τdlt;arctankslt;τa， 趨勢(shì)識(shí)別為平緩。 其中， τdlt;0和τagt;0為預(yù)先設(shè)置的常數(shù)。

對(duì)于航天器姿控系統(tǒng)輸出的時(shí)序預(yù)測(cè)， 采用6組神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器并聯(lián)的方式， 分別對(duì)三軸姿態(tài)角度和角速度進(jìn)行預(yù)測(cè)。 以俯仰角度為例， 假設(shè)當(dāng)前時(shí)刻為k， 給定一組輸入時(shí)序k， k－1， …， k－L+1， 其中L為時(shí)序的長(zhǎng)度。 為了預(yù)知故障信息， 開(kāi)展多步預(yù)測(cè)， 設(shè)置預(yù)測(cè)步長(zhǎng)為p， 采用迭代預(yù)測(cè)法對(duì)進(jìn)行預(yù)測(cè)， 具體步驟如下。

首先通過(guò)神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器多步輸入、 單步輸出生成第一步預(yù)測(cè)值：

^k+1=Pr（k， k－1， …， k+1－L）（17）

將前一步的輸出值作為已知量反饋至輸入中進(jìn)行下一步預(yù)測(cè)：

^k+2=Pr（^k+1， k， …， k+2－L）（18）

依次迭代計(jì)算出第p步預(yù)測(cè)值為

^k+p=Pr（^k+p－1， ^k+p－2， …， k+p－L）（19）

式中： ^k+i表示在當(dāng)前時(shí)刻k通過(guò)神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器產(chǎn)生的k+i時(shí)刻預(yù)測(cè)值， i=1， 2， …，p。 需要注意， 當(dāng)進(jìn)行中長(zhǎng)期預(yù)測(cè)時(shí)， 預(yù)測(cè)步長(zhǎng)p大于時(shí)序長(zhǎng)度L， 輸入神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器的時(shí)序數(shù)據(jù)均由預(yù)測(cè)值組成， 對(duì)預(yù)測(cè)模型的性能要求更高。

2" 融合預(yù)測(cè)算法

如圖3所示， 融合預(yù)測(cè)算法由神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器和KF構(gòu)成， 其預(yù)測(cè)精度很大程度上取決于神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器的性能和濾波器算法對(duì)時(shí)序預(yù)測(cè)誤差的魯棒性。 針對(duì)神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器， 通過(guò)小批量標(biāo)準(zhǔn)化處理、 加入趨勢(shì)識(shí)別模塊和超參數(shù)選取， 能夠有效提高其對(duì)時(shí)序預(yù)測(cè)的精度。 為解決時(shí)序預(yù)測(cè)誤差的魯棒性問(wèn)題， 在傳統(tǒng)KF算法中引入RLS參數(shù)估計(jì)取代增廣系統(tǒng)方式實(shí)現(xiàn)故障估計(jì)， 增強(qiáng)濾波器的魯棒性。

考慮執(zhí)行器故障， 航天器姿控系統(tǒng)線(xiàn)性離散狀態(tài)空間模型為［13］

x（k+1）=Adx（k）+Bdu（k）+Ff（k）+w（k）

y（k）=Cx（k）+v（k） （20）

式中： f（k）為執(zhí)行器故障； F=Bd； 假設(shè)過(guò)程噪聲w（k）和測(cè)量噪聲v（k）是相互獨(dú)立的均值為零的白噪聲， 協(xié)方差矩陣為Q（k）和R（k）， 且與初始狀態(tài)無(wú)關(guān)； Ad， Bd， C， F為已知適維矩陣。

針對(duì)一類(lèi)線(xiàn)性回歸問(wèn)題zr（k）=ΩT（k）hr+br， RLS參數(shù)估計(jì)算法為［14］

Λ（k）=［λI+ΩT（k）S（k－1）Ω（k）］－1

Γ（k）=S（k－1）Ω（k）Λ（k）

S（k）=1λS（k－1）－1λΓ（k）ΩT（k）S（k－1）

h^r（k）=h^r（k－1）+Γ（k）（zr（k）－ΩT（k）h^r（k－1）） （21）

式中： Λ（k）， Γ（k）和S（k）為遞歸過(guò)程中加入的輔助變量； λ為遺忘因子。

將hr（k）和br分別看作故障f（k）和新息y～（k）， 其他輔助變量和遺忘因子類(lèi)似地轉(zhuǎn)換到KF算法中。 假設(shè)當(dāng)前時(shí)刻為k， 得到改進(jìn)后的KF與神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器的融合預(yù)測(cè)算法如下， 其中預(yù)測(cè)步為

x^（k+1|k）=Adx^（k|k） +Bdu^（k）+Ff^（k）

P（k+1|k）=AdP（k|k）ATd+Q（k）（22）

式中： u^（k）=－diag（Kp， Kd）ylstm（k）； Kp和Kd為姿控系統(tǒng)輸出反饋PD控制參數(shù)。

根據(jù)神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器的輸出ylstm（k+1）， 更新步為

Σ（k+1）=CP（k+1|k）CT+R（k）

K（k+1）=P（k+1|k）CTΣ－1（k+1）

P（k+1|k+1）=［I－K（k+1）C］P（k+1|k）

γ（k+1）=［I－K（k+1）C］Adγ（k）+［I－K（k+1）C］F

Ω（k+1）=（CAdγ（k）+CF）T

Λ（k+1）=［λI+ΩT（k+1）S（k）Ω（k+1）］－1

Γ（k+1）=S（k）Ω（k+1）Λ（k+1）

S（k+1）=1λS（k）－1λΓ（k+1）ΩT（k+1）S（k）

y～（k+1）=ylstm（k+1）－Cx^（k+1|k）

f^（k+1）=f^（k）+Γ（k+1）y～（k+1）

x^（k+1|k+1）=x^（k+1|k）+K（k+1）y～（k+1）+

γ（k+1）［f^（k+1）－f^（k）］

（23）

式中： Ω（k）， Λ（k）， Γ（k）和S（k）與RLS算法中的變量一致； γ（k）為額外添加的輔助變量， 其余變量與傳統(tǒng)KF更新方程類(lèi)似。 根據(jù)KF的性質(zhì)， P（k|k）， Σ（k）， K（k）有界， 且Σ（k）gt;R（k）gt;0。

3" 仿真實(shí)驗(yàn)

為了驗(yàn)證所提算法的有效性， 建立以飛輪為執(zhí)行器的航天器姿控系統(tǒng)模型。 參考文獻(xiàn)［13］， 設(shè)置航天器主要參數(shù)如下：

主慣量矩陣Is=diag{18.73， 20.77， 23.63} kg·m2；

初始角度： ［0.084 7" －0.163 5" 0.124 8］T rad；

初始角速度： ［－0.041 6" 0.048 4" －0.055 6］T rad/s；

航天器軌道角速度ω0=0.001 2 rad/s；

陀螺常值漂移： 10－5I3 rad/s；

陀螺測(cè)量噪聲均方差： 3×10－5；

干擾力矩： 1.5×10－53cos（ω0t）+1

1.5sin（ω0t）+3cos（ω0t）3sin（ω0t）+1 N·m。

得到姿控系統(tǒng)狀態(tài)空間矩陣如下：

A=03×3 I3×3

I-1sΨ -I-1sΦ ， B=0I-1s，

Ψ=4ω20（Iy－Iz）00

03ω20（Ix－Iz）0

00ω20（Iy－Iz），

Φ=00－ω0（Ix－Iy+Iz）

000

ω0（Ix－Iy+Iz）00。

Ad和Bd分別為矩陣A和B在nk時(shí)刻的采樣值， Ad=expA（nk）T， Bd=Dd=∫0TexpA（nk）tdtB（nk）， Cd=Gd=C， 采樣周期T=0.1 s。

飛輪在工作時(shí)輸出的力矩為電磁驅(qū)動(dòng)力矩和飛輪軸承摩擦力矩合力矩的反作用力矩， 數(shù)學(xué)表示為

τw=－h(huán)·w=－（τm－τf）0 ""hw≤hmaxhwgt;hmax（24）

τm=Umτmaxsgn（Um） ""Um≤τmaxUmgt;τmax（25）

τf=τcsgn（Ωw）+τvΩw（26）

hw=IwΩw（27）

式中： τw為飛輪輸出力矩； hw為飛輪角動(dòng)量； τm為電磁驅(qū)動(dòng)力矩； τf為軸承摩擦力矩； hmax為飛輪最大角動(dòng)量； Um為控制指令； τc=0.002 N·m為軸承庫(kù)倫摩擦力矩； τv=3.60×10-4為粘性摩擦系數(shù)； Iw=0.007 8 kg·m2為飛輪的轉(zhuǎn)動(dòng)慣量； Ωw為飛輪角速度； sgn（·）為符號(hào)函數(shù)； 飛輪最大電磁驅(qū)動(dòng)力矩τmax=0.4 N·m。

3.1" 神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器性能分析

3.1.1" 數(shù) 據(jù) 集

訓(xùn)練數(shù)據(jù)集來(lái)自航天器姿控系統(tǒng)模型的輸出。 姿控系統(tǒng)大約30 s后進(jìn)入穩(wěn)定狀態(tài)， 為保證數(shù)據(jù)集的代表性， 選擇在平穩(wěn)運(yùn)行20 s后， 即50 s的時(shí)刻作為采樣起點(diǎn)。 假設(shè)預(yù)測(cè)的起始時(shí)刻為T(mén)p， 則利用50 s至Tp時(shí)間段內(nèi)的系統(tǒng)輸出構(gòu)建數(shù)據(jù)集。 航天器實(shí)際在軌時(shí)， 可通過(guò)地面遙測(cè)得到的某一時(shí)間段的姿態(tài)信息構(gòu)建數(shù)據(jù)集。 實(shí)驗(yàn)中， 設(shè)置仿真時(shí)間Tf=500 s， 采樣步長(zhǎng)T=0.1 s， 將采集到的時(shí)序數(shù)據(jù)經(jīng)過(guò)歸一化處理后表示為yTs/T， yTs/T+1， …， yTp/T， 其中Ts=50 s為采樣開(kāi)始的時(shí)刻。 在此基礎(chǔ)上分別構(gòu)建輸入數(shù)據(jù)集和對(duì)應(yīng)的期望輸出數(shù)據(jù)集：

P1=［yTs/T， yTs/T+1， …， yTs/T+L－1］→T1=yTs/T+L

P2=［yTs/T+1， yTs/T+2， …， yTs/T+L］→T2=yTs/T+L+1

Pn=［yTs/T+n－1， yTs/T+n， …， yTs/T+n+L－2］→Tn=yTs/T+n+L－1

式中： Pn為輸入樣本； Tn為期望輸出樣本； L為時(shí)序長(zhǎng)度； n為所采集數(shù)據(jù)產(chǎn)生的樣本數(shù)。 輸入數(shù)據(jù)集和期望輸出數(shù)據(jù)集需要依次輸入到神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器， 設(shè)置趨勢(shì)識(shí)別模塊的閾值τd=－0.5°， τa=0.5°， 然后通過(guò)趨勢(shì)識(shí)別模塊分別訓(xùn)練預(yù)測(cè)神經(jīng)網(wǎng)絡(luò)。

3.1.2" 超參數(shù)選取

預(yù)測(cè)神經(jīng)網(wǎng)絡(luò)的收斂性和預(yù)測(cè)精度受學(xué)習(xí)率、 LSTM層細(xì)胞元數(shù)量和小批量樣本數(shù)量的影響。

學(xué)習(xí)率用于控制神經(jīng)網(wǎng)絡(luò)權(quán)值更新的步長(zhǎng)， 考慮可變學(xué)習(xí)率， 設(shè)置初始學(xué)習(xí)率為0.005， 每進(jìn)行50次迭代乘衰減因子0.2來(lái)降低學(xué)習(xí)率。

通過(guò)對(duì)比40， 60， 80， 100四種不同細(xì)胞元數(shù)訓(xùn)練收斂速度， 最終選取60作為上升趨勢(shì)和下降趨勢(shì)對(duì)應(yīng)的兩個(gè)預(yù)測(cè)網(wǎng)絡(luò)LSTM層的細(xì)胞元數(shù)， 選取100作為平緩趨勢(shì)對(duì)應(yīng)的預(yù)測(cè)網(wǎng)絡(luò)LSTM層細(xì)胞元數(shù)。

小批量樣本數(shù)量通常設(shè)置為2的n次冪。 考慮網(wǎng)絡(luò)結(jié)構(gòu)的復(fù)雜度和樣本數(shù)， 設(shè)置上升趨勢(shì)和下降趨勢(shì)對(duì)應(yīng)兩個(gè)LSTM網(wǎng)絡(luò)的小批量樣本數(shù)量為16， 平緩趨勢(shì)對(duì)應(yīng)LSTM網(wǎng)絡(luò)的小批量樣本數(shù)量為64。

3.1.3" 時(shí)序預(yù)測(cè)性能比較

三種方法在預(yù)測(cè)段的預(yù)測(cè)結(jié)果及局部放大圖、 誤差曲線(xiàn)分別如圖4～7所示。 根據(jù)預(yù)測(cè)誤差計(jì)算三種方法輸出結(jié)果的RMSE分別為3.001 1×10-6， 9.164 0×10-7， 2.163 9×10-8。 結(jié)合三種方法的預(yù)測(cè)結(jié)果， 能夠直觀(guān)地看出兩種基于LSTM網(wǎng)絡(luò)的預(yù)測(cè)器相比基于傳統(tǒng)RNN的預(yù)測(cè)器在時(shí)序預(yù)測(cè)方面具備更出色的性能。 但單一LSTM網(wǎng)絡(luò)預(yù)測(cè)器短期內(nèi)預(yù)測(cè)精度較高， 在420 s后， 預(yù)測(cè)步長(zhǎng)

大于預(yù)測(cè)器輸入時(shí)序長(zhǎng)度， 預(yù)測(cè)誤差開(kāi)始逐漸增大。 而

所提方法能夠持續(xù)平穩(wěn)地預(yù)測(cè)輸出時(shí)序，" 在480" s后性

能開(kāi)始衰退， 但誤差仍比單一LSTM網(wǎng)絡(luò)預(yù)測(cè)器小一個(gè)量級(jí)， 處于可接受的范圍內(nèi)。 綜合以上實(shí)驗(yàn)結(jié)果得出結(jié)論， 利用所提出的神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器可以獲得更為持續(xù)和準(zhǔn)確的預(yù)測(cè)值， 有利于對(duì)基于KF的融合預(yù)測(cè)提供輔助更新作用。

3.2" 故障預(yù)測(cè)實(shí)驗(yàn)

故障預(yù)測(cè)主要針對(duì)緩變故障， 根據(jù)當(dāng)前或歷史數(shù)據(jù)特征對(duì)系統(tǒng)故障未來(lái)的發(fā)展趨勢(shì)進(jìn)行預(yù)測(cè)、 分析和判斷， 以便提前消除故障隱患［15］。 文獻(xiàn)［16］表明潤(rùn)滑劑不足是導(dǎo)致飛輪失效的關(guān)鍵因素， 故障表現(xiàn)為軸承粘性摩擦系數(shù)的異常。 因此， 為了驗(yàn)證融合預(yù)測(cè)算法對(duì)緩變故障的有效性， 在仿真中假設(shè)俯仰軸飛輪在100 s時(shí)發(fā)生潤(rùn)滑劑泄露故障， 設(shè)置軸承粘性摩擦系數(shù)τv的增量τ′v， 其表達(dá)式為

τ′v（t）=0""""""""" t≤100

0.1（exp10－5（t－100）－1）tgt;100 （28）

實(shí)驗(yàn)參數(shù)設(shè)置為： P（0）=0.1In， Q=0.1In， R=0.1· In， γ（0）=0， S（0）=0.05In， λ=0.95。 初值f^（0）=［0" 0" 0］T， x^（0）=［0" 0" 0" 0" 0" 0］T。 在仿真400 s預(yù)測(cè)算法介入前， 利用改進(jìn)KF算法對(duì)航天器姿控系統(tǒng)狀態(tài)量和故障值（即飛輪軸承摩擦力矩）進(jìn)行更新和估計(jì)； 400 s后， 使用融合預(yù)測(cè)算法對(duì)二者未來(lái)變化進(jìn)行預(yù)測(cè)。 仿真結(jié)果如圖8～9所示， 圖中均對(duì)預(yù)測(cè)誤差進(jìn)行了局部放大處理。

根據(jù)仿真結(jié)果得出結(jié)論： 在估計(jì)段， 使用改進(jìn)KF能夠快速跟蹤姿控系統(tǒng)的六個(gè)狀態(tài)量和飛輪軸承摩擦力矩； 在預(yù)測(cè)段， 使用融合預(yù)測(cè)算法可以根據(jù)歷史數(shù)據(jù)準(zhǔn)確預(yù)測(cè)系統(tǒng)未來(lái)的變化趨勢(shì)。 通過(guò)對(duì)未來(lái)100 s飛輪軸承摩擦力矩的分析可以看出， 俯仰軸飛輪摩擦力矩逐漸增大， 但由于軸承粘性摩擦系數(shù)的變化非常緩慢， 且摩擦力矩還與飛輪轉(zhuǎn)速相關(guān)， 因此難以通過(guò)其早期異常導(dǎo)致的摩擦力矩微小變化直接反映出故障。

為了更好地判斷飛輪是否將發(fā)生故障以及發(fā)生故障的時(shí)間， 根據(jù)圖3預(yù)警框架進(jìn)一步對(duì)飛輪軸承摩擦系數(shù)進(jìn)行預(yù)測(cè)分析。 首先使用所設(shè)計(jì)的預(yù)測(cè)網(wǎng)絡(luò)對(duì)飛輪轉(zhuǎn)速Ωw進(jìn)行預(yù)測(cè)， 然后利用融合預(yù)測(cè)算法得到的軸承摩擦力矩τf， 根據(jù)式（26）即可計(jì)算得到軸承粘性摩擦系數(shù)τv， 最后對(duì)τv進(jìn)行閾值分析實(shí)現(xiàn)故障情況的判斷。 設(shè)定τv超過(guò)標(biāo)準(zhǔn)值3.60×10－4的20%則為故障情況， 由此設(shè)置閾值Dτv=4.32×10－4， 持續(xù)預(yù)測(cè)τv直到數(shù)值超過(guò)閾值后停止。 實(shí)驗(yàn)結(jié)果如圖10～11所示。

根據(jù)預(yù)測(cè)得到的軸承粘性摩擦系數(shù)變化曲線(xiàn)可以判斷飛輪軸承性能發(fā)生了退化， 預(yù)計(jì)在856 s達(dá)到閾值， 此時(shí)將對(duì)航天器姿態(tài)穩(wěn)定產(chǎn)生明顯影響。 而根據(jù)式（28）設(shè)定的粘性摩擦系數(shù)故障增量， 飛輪將在820 s表現(xiàn)出故障， 預(yù)測(cè)時(shí)間與實(shí)際相差36 s， 結(jié)果有效。 綜合以上實(shí)驗(yàn)結(jié)果得出結(jié)論， 在發(fā)生緩變故障后， 神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器能夠通過(guò)系統(tǒng)歷史信息來(lái)識(shí)別故障特征， 并利用融合預(yù)測(cè)算法有效預(yù)測(cè)系統(tǒng)短期和長(zhǎng)期的變化趨勢(shì)， 通過(guò)獲得更多的故障信息實(shí)現(xiàn)對(duì)飛輪軸承性能退化故障的預(yù)警。

4" 結(jié)" 論

本文提出一種基于濾波器和LSTM的融合故障預(yù)測(cè)算法用于航天器飛輪執(zhí)行器故障的預(yù)測(cè)。 主要貢獻(xiàn)在于：

（1） 設(shè)計(jì)了小批量標(biāo)準(zhǔn)化LSTM網(wǎng)絡(luò)解決協(xié)變量漂移問(wèn)題， 并與所設(shè)計(jì)的趨勢(shì)識(shí)別模塊串聯(lián)， 有效提升對(duì)姿控系統(tǒng)輸出時(shí)序預(yù)測(cè)的精度。

（2） 在傳統(tǒng)KF更新步中引入RLS參數(shù)估計(jì)取代擴(kuò)展?fàn)顟B(tài)向量的方式實(shí)現(xiàn)故障估計(jì)， 提高了濾波器輸出對(duì)時(shí)序預(yù)測(cè)誤差的魯棒性。

（3）" 將所設(shè)計(jì)的神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)器與改進(jìn)后的濾波器相結(jié)合， 動(dòng)態(tài)調(diào)整濾波器的更新過(guò)程以修正預(yù)測(cè)結(jié)果， 并通過(guò)飛輪軸承性能退化故障仿真結(jié)果驗(yàn)證了所提算法的有效性。

參考文獻(xiàn)：

［1］ ElDali M，" Kumar K D. Fault Diagnosis and Prognosis of Aerospace Systems Using Growing Recurrent Neural Networks and LSTM［C］∥ IEEE Aerospace Conference， 2021： 1-20.

［2］ 林燕. 航空航天領(lǐng)域故障預(yù)測(cè)與健康管理技術(shù)研究綜述［J］. 中國(guó)航天， 2019（5）： 21-26.

Lin Yan. Summary of Research on Fault Prediction and Health Management Technology in Aerospace Field［J］. Aerospace China， 2019（5）： 21-26. （in Chinese）

［3］ Kordestani M， Saif M， Orchard M E， et al. Failure Prognosis and Applications—A Survey of Recent Literature［J］. IEEE Transactions on Reliability， 2021， 70（2）： 728-748.

［4］ Son J， Zhou S Y， Sankavaram C， et al. Remaining Useful Life Prediction Based on Noisy Condition Monitoring Signals Using Constrained Kalman Filter［J］. Reliability Engineering amp; System Safety， 2016， 152： 38-50.

［5］ Kordestani M， Samadi M F， Saif M， et al. A New Fault Prognosis of MFS System Using Integrated Extended Kalman Filter and Baye-sian Method［J］. IEEE Transactions on Industrial Informatics， 2018： 1.

［6］ Jun C， Lu T. Gyro Fault Prediction Algorithm Based on UKF［C］∥ International Conference on Image Analysis and Signal Processing， 2012： 1-4.

［7］ Rahimi A， Kumar K D， Alighanbari H. Failure Prognosis for Satellite Reaction Wheels Using Kalman Filter and Particle Filter［J］. Journal of Guidance Control Dynamics， 2020， 43（3）： 585-588.

［8］ Huang C G， Yin X H， Huang H Z， et al. An Enhanced Deep Learning-Based Fusion Prognostic Method for RUL Prediction［J］. IEEE Transactions on Reliability， 2020， 69（3）： 1097-1109.

［9］ 王鑫， 吳際， 劉超， 等. 基于LSTM循環(huán)神經(jīng)網(wǎng)絡(luò)的故障時(shí)間序列預(yù)測(cè)［J］. 北京航空航天大學(xué)學(xué)報(bào)， 2018， 44（4）： 772-784.

Wang Xin， Wu Ji， Liu Chao， et al. Exploring LSTM Based Recurrent Neural Network for Failure Time Series Prediction［J］. Journal of Beijing University of Aeronautics and Astronautics， 2018， 44（4）： 772-784. （in Chinese）

［10］ Meng X， Song H. Fault Prediction of Satellite Attitude Control System Based on Neural Network［J］. Journal of System Simulation， 2019， 31（11）： 10.

［11］ Sirajul Islam M， Rahimi A. Use of a Data-Driven Approach for Time Series Prediction in Fault Prognosis of Satellite Reaction Wheel［C］∥IEEE International Conference on Systems， Man， and Cybernetics （SMC）， 2020： 3624-3628.

［12］ Wang J R， Li S M， An Z H， et al. Batch-Normalized Deep Neural Networks for Achieving Fast Intelligent Fault Diagnosis of Machines［J］. Neurocomputing， 2019， 329： 53-65.

［13］ 賈慶賢. 基于學(xué)習(xí)觀(guān)測(cè)器的衛(wèi)星姿態(tài)控制系統(tǒng)故障重構(gòu)研究［D］. 哈爾濱： 哈爾濱工業(yè)大學(xué)， 2015.

Jia Qingxian. Research on Learning Observer-Based Fault Reconstruction for Satellite Attitude Control Systems［D］.Harbin： Harbin Institute of Technology， 2015. （in Chinese）

［14］ Zhang Q H. Adaptive Kalman Filter for Actuator Fault Diagnosis［J］. Automatica， 2018， 93： 333-342.

［15］ 鄭秀娟. 基于狀態(tài)估計(jì)的多方法融合的故障預(yù)測(cè)算法研究［D］. 武漢： 華中科技大學(xué)， 2016.

Zheng Xiujuan. Research on Fault Prediction Algorithm Based on State Estimation and Multi-Methods Integration［D］.Wuhan： Huazhong University of Science and Technology， 2016. （in Chinese）

［16］ Li Y F， Huang M， Bai S， et al. Reliability Modeling and Analysis of Multi-Degradation of Momentum Wheel Based on Copula Function［J］. Applied Sciences， 2021， 11（23）： 11563.

Flywheel Fault Prediction Algorithm Based on

Improved Filter and LSTM Fusion

Long Dizhi*， Chen Xin， Wei Bingyi， Shi Chao

（Beijing Aerospace Automatic Control Institute， Beijing 100854， China）

Abstract： To solve the problem that the poor accuracy of long-term prediction relying on model accuracy for the filter-based prediction method， a fusion fault prediction algorithm based on" filter and long short-term memory （LSTM） network is proposed to achieve the prediction of spacecraft flywheel slowly growing faults. Firstly， a mini-batch normalization LSTM network and a trend recognition module are designed， which are connected in series to form a neural network predictor to improve the time series prediction accuracy. Then， the Kalman filter update process is improved by the recursive least square （RLS） parameter estimation principle to enhance the robustness for time series prediction error. On this basis， the predicted values output by the neural network predictor are fused with the improved filter. Prediction residual can be obtained for iterative updating and prediction， overcoming the dependence of the filter algorithm on the model and improving prediction accuracy. Finally， the time series predictive performance of three neural network predictors are analyzed by simulation experiments. And considering the degradation fault of flywheel bearing performance， the proposed fusion prediction algorithm is used to determine that the flywheel degrades to a threshold at 856 s， with a prediction time error of 36 s. The simulation results verify the effectiveness of the proposed algorithm for slowly growing fault.

Key words：" fault prediction; LSTM; Kalman filter; RLS; flywheel