陳雄姿，于勁松，2，唐荻音，李行善

（1．北京航空航天大學(xué)自動化科學(xué)與電氣工程學(xué)院，北京100191；2．先進(jìn)航空發(fā)動機(jī)協(xié)同創(chuàng)新中心，北京100191）

基于改進(jìn)核平滑輔助粒子濾波的失效預(yù)測方法

陳雄姿1，于勁松1，2，唐荻音1，李行善1

（1．北京航空航天大學(xué)自動化科學(xué)與電氣工程學(xué)院，北京100191；2．先進(jìn)航空發(fā)動機(jī)協(xié)同創(chuàng)新中心，北京100191）

針對系統(tǒng)模型存在多個(gè)未知參數(shù)的情況，提出了一種基于改進(jìn)核平滑輔助粒子濾波（improved kernel smoothing auxiliary particle filtering，IKS-APF）的失效預(yù)測方法。首先，在已有核平滑輔助粒子濾波基礎(chǔ)上引入增益因子和加速因子，使其具有參數(shù)方差雙向調(diào)節(jié)能力和更快的參數(shù)估計(jì)收斂速度。然后，使用ISK-APF進(jìn)行狀態(tài)和參數(shù)的聯(lián)合估計(jì)，為確保參數(shù)估計(jì)的準(zhǔn)確性同時(shí)減少參數(shù)的不確定性，設(shè)計(jì)了方差監(jiān)視和短期預(yù)測誤差匹配相結(jié)合的自適應(yīng)粒子方差控制方案。最后，使用最新估計(jì)到的狀態(tài)和參數(shù)粒子進(jìn)行迭代預(yù)測，并通過統(tǒng)計(jì)狀態(tài)粒子首達(dá)失效狀態(tài)空間的時(shí)間計(jì)算出剩余使用壽命（remaining useful life，RUL）。仿真結(jié)果證明了本文方法的有效性和優(yōu)越性。

失效預(yù)測；核平滑；輔助粒子濾波；參數(shù)估計(jì)；方差控制；不確定性管理

0 引 言

航空航天和武器裝備領(lǐng)域應(yīng)用預(yù)測與健康管理（prognostics and health management，PH M）技術(shù)具有提高對象系統(tǒng)的安全性、可靠性和降低維護(hù)成本的顯著意義，該技術(shù)自21世紀(jì)初被提出以來，越來越受到行業(yè)內(nèi)的重視［1］。故障／失效預(yù)測技術(shù)作為PHM的核心技術(shù)，一直是研究的熱點(diǎn)和難點(diǎn)。現(xiàn)有的失效預(yù)測方法主要可以分為統(tǒng)計(jì)方法、人工智能方法和基于模型的方法［2］，前兩者又可以統(tǒng)一歸為數(shù)據(jù)驅(qū)動的方法?；谀Ｐ偷姆椒ň哂蓄A(yù)測準(zhǔn)確度高、計(jì)算復(fù)雜度小的優(yōu)點(diǎn)，在能夠獲得系統(tǒng)失效模型的前提下，往往是預(yù)測方法的首選。

基于狀態(tài)空間演化模型的粒子濾波（particle filtering，PF）方法，是當(dāng)前一種較為流行的基于模型的預(yù)測方法。它一方面適用于非線性／非高斯系統(tǒng)，另一方面具有概率性的輸出，有利于預(yù)測不確定性的表示和管理。文獻(xiàn)［3］較早提出了一個(gè)完整的基于PF的在線故障預(yù)測框架；文獻(xiàn)［4］將PF成功應(yīng)用到鋰電池的健康預(yù)測中；文獻(xiàn)［5］為解決樣本貧化問題，提出了一種權(quán)值選優(yōu)的故障預(yù)測算法；文獻(xiàn)［6］提出了一種免重采樣的高斯混合模型PF故障預(yù)測算法；文獻(xiàn)［7］設(shè)計(jì)了3種不同權(quán)值的剩余使用壽命（remaining useful life，RUL）PF有偏估計(jì)器；文獻(xiàn)［8］針對模型參數(shù)持續(xù)時(shí)變的情況，提出了一種PF與最小二乘支持向量回歸相融合的失效預(yù)測框架；文獻(xiàn)［9］使用無跡PF來預(yù)測鋰電池的剩余壽命。然而，上述文獻(xiàn)中使用到的系統(tǒng)模型一般要求其所有參數(shù)都是已知的［3，5，7］，或者即使包含動態(tài)可變參數(shù)，也一般假設(shè)參數(shù)的初始值和隨機(jī)游走模型參數(shù)的方差是已知的［4，6，8-9］。系統(tǒng)物理模型是基于模型的預(yù)測方法的根本。在實(shí)際應(yīng)用中，由于設(shè)備工作負(fù)載、工作環(huán)境、工作時(shí)長、個(gè)體屬性等差異，同一類設(shè)備不同個(gè)體的模型也不盡相同，部分模型參數(shù)很多時(shí)候只能獲悉其大概的取值范圍。當(dāng)模型參數(shù)的不確定性變大時(shí)，預(yù)測結(jié)果的不確定性也會隨之變大。因此，研究模型包含多個(gè)未知參數(shù)情況下如何準(zhǔn)確估計(jì)模型參數(shù)并盡可能減少其不確定性，從而提升設(shè)備RUL預(yù)測的質(zhì)量具有重要意義。

為解決上述問題，本文提出了一種基于改進(jìn)核平滑輔助粒子濾波（improved kernel smoothing auxiliary particle filtering，IKS-APF）的失效預(yù)測方法。使用IKS-APF算法聯(lián)合估計(jì)系統(tǒng)的狀態(tài)和參數(shù)，并設(shè)計(jì)了基于方差監(jiān)視和短期預(yù)測誤差控制相結(jié)合的參數(shù)方差雙階段自適應(yīng)調(diào)節(jié)策略，在準(zhǔn)確估計(jì)參數(shù)的同時(shí)減少其不確定性；最后，基于最新的狀態(tài)和參數(shù)估計(jì)進(jìn)行迭代實(shí)現(xiàn)設(shè)備RUL的預(yù)測。通過裂紋增長仿真實(shí)驗(yàn)對所提出的方法進(jìn)行了驗(yàn)證。

1 問題描述

假設(shè)系統(tǒng)的模型在tk＝kΔt時(shí)間點(diǎn)（簡稱k時(shí)刻）可以表示為

式中，xk∈Rnx是狀態(tài)向量，又可稱為健康指征，一般與系統(tǒng)的退化直接相關(guān)；θk∈Rnθ是參數(shù)向量，假設(shè)參數(shù)的初始值和變化方式未知，只知道參數(shù)的大概取值范圍，同時(shí)假設(shè)參數(shù)是靜態(tài)的或緩慢變化的，這種情況在實(shí)際應(yīng)用中也最為常見；uk∈Rnu是輸入向量，一般表示系統(tǒng)的外部工作條件；nk∈Rnn是過程噪聲向量；fk是狀態(tài)方程；zk∈Rnz是觀測向量；hk是觀測方程；vk∈Rnv是觀測噪聲向量。hk和fk可以是非線性函數(shù)，nk和vk可以是非高斯噪聲，且都是已知的。

k時(shí)刻的失效預(yù)測就是基于傳感器觀測序列z1∶k＝｛z1，z2，…，zk｝，估計(jì)當(dāng)前的狀態(tài)xk和參數(shù)θk值，并迭代預(yù)測未來時(shí)刻的狀態(tài)xk＋n（n＝1，2，…）直至其達(dá)到失效閾值，該時(shí)刻就是系統(tǒng)的失效時(shí)間。根據(jù)貝葉斯準(zhǔn)則，同時(shí)估計(jì)狀態(tài)和參數(shù)就是要求解它們的聯(lián)合概率分布

為了同時(shí)估計(jì)狀態(tài)和參數(shù)，常采用聯(lián)合估計(jì)方法［4，6，8-9］，將系統(tǒng)狀態(tài)和參數(shù)合并為一個(gè)增廣向量x－k＝［xkTθkT］T，模型式（1）轉(zhuǎn)變?yōu)?/p>

這里假設(shè)未知參數(shù)θk服從一個(gè)隨機(jī)游走模型，模型參數(shù)ξk～N（0，Σk），即ξk服從均值為0、方差為Σk的高斯隨機(jī)分布，Σk是ξk對應(yīng)的協(xié)方差矩陣，簡稱為超參數(shù)。PF方法得到的參數(shù)粒子越接近其真實(shí)值，分配得到的權(quán)值越大，并最終收斂于真實(shí)值。超參數(shù)Σk的大小決定了參數(shù)θk收斂的速度和收斂后的估計(jì)性能。若Σk設(shè)置較大，收斂的速度較快，但是收斂后，θk的方差較大，即不確定性較大；若Σk設(shè)置過小，收斂速度變慢，但是一旦收斂，參數(shù)的方差比較小，跟蹤估計(jì)精度較高。

為了減少預(yù)測的不確定性，文獻(xiàn)［10- 11］基于短期預(yù)測誤差的外環(huán)校正來自適應(yīng)調(diào)節(jié)超參數(shù)Σk，預(yù)測誤差滿足閾值條件，縮小超參數(shù)，反之則放大超參數(shù)。不過該方法有兩個(gè)明顯的不足：首先，它沒有明確地給出啟用外環(huán)校正的條件，如果從起始時(shí)就啟用，在參數(shù)初始值未知，超參數(shù)初始設(shè)置不合適的時(shí)候，短期預(yù)測必然不能滿足誤差要求，這時(shí)就將陷入不斷擴(kuò)大超參數(shù)的惡性循環(huán)；其次，該方法只應(yīng)用在了單一未知參數(shù)的情況，當(dāng)存在多個(gè)未知參數(shù)時(shí)，該方法根據(jù)短期預(yù)測效果對所有超參數(shù)進(jìn)行完全相同的操作，顯然是不合適的。文獻(xiàn)［12］為了處理多失效模式并存的情況，提出基于參數(shù)θk粒子群方差的多階段自適應(yīng)控制方法，各階段均包括3個(gè)重要參數(shù)：理想的方差值、進(jìn)入下一階段的方差閾值和用于控制自適應(yīng)調(diào)節(jié)程度的比例增益。該方法可以處理多未知參數(shù)的情況，然而各階段參數(shù)的設(shè)置并無統(tǒng)一方法，特別是比例增益因子，很難確定；另外，該方法只能控制方差變小，一旦估計(jì)出現(xiàn)偏差，并不能擴(kuò)大超參數(shù)，重新達(dá)到收斂，因此其魯棒性并無保證。

為了解決上述問題，本文提出一種基于IKS-APF且具有超參數(shù)自適應(yīng)調(diào)節(jié)能力的新的失效預(yù)測框架。

2 改進(jìn)的核平滑輔助粒子濾波

文獻(xiàn)［13］采用APF與參數(shù)核平滑近似相結(jié)合來估計(jì)狀態(tài)和參數(shù)的聯(lián)合后驗(yàn)分布，該方法是一種聯(lián)合估計(jì)系統(tǒng)狀態(tài)和靜態(tài)未知參數(shù)的有效方法，簡稱為核平滑輔助粒子濾波方法。雖然預(yù)測問題中系統(tǒng)的未知參數(shù)可能是小幅波動或者緩慢變化的，但在任意特定需要進(jìn)行參數(shù)估計(jì)的短暫時(shí)間片內(nèi)，完全可以將它們當(dāng)作靜態(tài)參數(shù)來對待。然而，若直接將該方法應(yīng)用到失效預(yù)測問題，參數(shù)粒子的方差會持續(xù)縮小至0，將不能動態(tài)跟蹤參數(shù)的變化和表示參數(shù)估計(jì)的不確定性。因此，本節(jié)對原有參數(shù)核平滑方法進(jìn)行了改進(jìn)，引入了增益因子和加速因子，并進(jìn)一步得到IKSAPF算法。

2.1 輔助粒子濾波

最常用的PF算法是采樣重要性重采樣（sampling importance resampling，SIR）算法，它使用先驗(yàn)分布p（xk|）作為重要性密度函數(shù)（其中｛，i＝1，2，…，Ns｝是k－1時(shí)刻對應(yīng)權(quán)值為｛，i＝1，2，…，Ns｝的粒子群，Ns為粒子數(shù)目），抽樣過程簡單易求，但該密度函數(shù)獨(dú)立于當(dāng)前時(shí)刻的觀測值zk，僅僅根據(jù)粒子運(yùn)動和先前狀態(tài)盲目抽樣，使其對異常值過于靈敏，丟失大量低權(quán)重粒子。文獻(xiàn)［14］提出采用一種更加可靠的重要性密度函數(shù)

式中，μik＝E［xk|xik－1］，i是粒子索引，同時(shí)也是導(dǎo)出重要性密度的輔助變量。μik的引入，使得k時(shí)刻的粒子xik是基于k－1時(shí)刻的先驗(yàn)粒子和k時(shí)刻的觀測值共同產(chǎn)生的，更加接近于真實(shí)的狀態(tài)值。相應(yīng)的權(quán)值為

歸一化此權(quán)值后重采樣得到具有等權(quán)值的粒子集｛xjk，j＝1，2，…，Ns｝。另外根據(jù)貝葉斯準(zhǔn)則有

由式（4）和式（6）可以得到

這就是APF，它的重要性密度函數(shù)考慮了最新的觀測信息，更接近真實(shí)狀態(tài)。理論上，在過程噪聲較小的條件下，可以獲得更好的估計(jì)性能。

2.2 參數(shù)核平滑及其改進(jìn)

該分布的均值為θˉk－1，協(xié)方差矩陣為Vk－1＋Σk。由此可見，隨機(jī)游走模型致使協(xié)方差不斷增大。于是，文獻(xiàn)［13］采用核平滑的方法來減小協(xié)方差，引入平滑因子h（0＜h＜1）

此時(shí)，超參數(shù)Σk＝h2Vk－1，核心mik－1用收縮規(guī)則迫使參數(shù)粒子向其均值靠攏

可以直接證明式（9）中混合分布的均值為θˉk－1，協(xié)方差矩陣為Vk－1，和原來的大小保持一致。不過，由于核縮的作用，參數(shù)粒子的方差將持續(xù)減小直至為0。在設(shè)備失效預(yù)測的過程中，由于參數(shù)的動態(tài)特性、估計(jì)的不準(zhǔn)確性（PF本身就是一種近似估計(jì)）以及各種噪聲的影響，不可能也無需獲得精確的參數(shù)估計(jì)值，而應(yīng)該保有一定的參數(shù)粒子方差來跟蹤其動態(tài)和表征它的不確定性。因此，本文對原有的核平滑方法進(jìn)行改進(jìn)，引入加速因子q和增益因子γ，式（9）和式（10）分別變?yōu)槭剑?1）和式（12）：

為了處理第1節(jié)中提到的參數(shù)初始值未知、僅知大概取值范圍條件下的參數(shù)估計(jì)問題，對式（10）中mik－1的定義做了改進(jìn)，引入了加速因子q（0＜q＜1），加塊收斂速度。同時(shí)，超參數(shù)Σk*中添加了增益因子γ（γ≥1），它的引入使得核平滑不再僅能實(shí)現(xiàn)方差的縮小，通過設(shè)置它也能實(shí)現(xiàn)方差的放大，這種雙向調(diào)節(jié)能力對于參數(shù)方差自適應(yīng)調(diào)節(jié)是非常重要的。

為了避免持續(xù)核縮導(dǎo)致參數(shù)方差減小為0，當(dāng)參數(shù)粒子方差達(dá)到理想的閾值時(shí)，應(yīng)停止核縮，回到使用傳統(tǒng)隨機(jī)游走模型的跟蹤估計(jì)階段。此時(shí)式（12）中參數(shù)的值變?yōu)?/p>

2.3 改進(jìn)核平滑輔助粒子濾波聯(lián)合估計(jì)算法

將改進(jìn)的核平滑方法融入到第2.1節(jié)中的APF算法，可得到IKS-APF狀態(tài)和參數(shù)聯(lián)合估計(jì)算法，算法步驟如下：

步驟1 初始化：取k＝0，按p（x0）和p（θ0）分別抽取Ns個(gè)樣本粒子xi0和θi0，i＝1，2，…，Ns；

步驟2 對于i＝1，2，…，Ns的每一個(gè)粒子，根據(jù)式（4）計(jì)算，根據(jù)式（12）或者式（13）計(jì)算m和Σ，再根據(jù)式（5）計(jì)算；

步驟4 j＝1，2，…，Ns，采樣N（）得到，采樣pk（x|，） 得到， 計(jì)算權(quán)值

步驟5 歸一化權(quán)值wjk；

步驟6 設(shè)置k＝k＋1，當(dāng)下一次量測值到來時(shí)，繼續(xù)重復(fù)執(zhí)行步驟2～步驟6。

3 IKS-APF預(yù)測及其不確定性管理

基于IKS-APF的失效預(yù)測框架主要包括3個(gè)階段。假設(shè)當(dāng)前時(shí)刻為λ，已知觀測序列z1∶λ＝｛z1，z2，…，zλ｝。

3.1 第1階段（k=0）：粒子初始化

PF算法首先需要生成等權(quán)值的初始狀態(tài)空間粒子，假設(shè)狀態(tài)初始值x0已知。對于模型參數(shù)θ，若先驗(yàn)分布已知，直接取樣即可；若只知其大致取值區(qū)間，可以在區(qū)間上進(jìn)行均勻隨機(jī)采樣，設(shè)粒子數(shù)目為Ns。

3.2 第2階段（1≤k≤λ）：狀態(tài)與參數(shù)聯(lián)合估計(jì)

3.2.1 參數(shù)方差自適應(yīng)調(diào)節(jié)（1≤k≤λ－p*）

根據(jù)初始粒子和觀測序列，在各時(shí)刻運(yùn)用IKS-APF對狀態(tài)和參數(shù)進(jìn)行聯(lián)合估計(jì)。由于核縮的作用，參數(shù)方差逐步縮小，趨于收斂。值得注意的是：上述方法如無約束條件，參數(shù)方差將減小至0。一方面，需要減小參數(shù)方差來提高預(yù)測的精度；另一方面，由于各種噪聲的存在，PF作為一種離散近似估計(jì)，即使模型參數(shù)是靜態(tài)的，參數(shù)估計(jì)過程也存在不確定性，更何況參數(shù)經(jīng)常會隨著負(fù)載和環(huán)境的變化而緩慢變化或小幅波動，需要保有一定方差維持動態(tài)更新。因此，這里參數(shù)的精度通過監(jiān)視各參數(shù)粒子的方差來控制，參數(shù)估計(jì)的準(zhǔn)確性通過短期預(yù)測誤差來判斷（p*為短期預(yù)測的步數(shù)）。參數(shù)方差自適應(yīng)調(diào)節(jié)策略如圖1所示。

圖1 參數(shù)方差雙階段自適應(yīng)調(diào)節(jié)策略圖

參數(shù)的方差采用文獻(xiàn)［12］中的相對中值絕對偏差（relative median absolute deviation，RMAD）來度量，一方面它是一個(gè)魯棒統(tǒng)計(jì)量，另一方面它是個(gè)相對量，相同類參數(shù)可以等同設(shè)置。其定義如下：

假設(shè)未知參數(shù)θ的維數(shù)是d，為了能夠?qū)Ω鱾€(gè)參數(shù)進(jìn)行獨(dú)立調(diào)節(jié)，每個(gè)參數(shù)可以設(shè)有獨(dú)立的平滑因子h、增益因子γ和加速因子q，分別是｛，…，｝、｛，，…，γ｝和｛q10，q20，…，qd0｝。雙階段自適應(yīng)調(diào)節(jié)策略如下：

第1階段 設(shè)定各未知參數(shù)RMAD高限閾值：｛RT1up，RT2up，…，RTdup｝，任意k時(shí)刻，使用IKS-APF方法估計(jì)狀態(tài)和參數(shù)后，計(jì)算各參數(shù)的RMAD值Rjk（j＝1，2，…，d），并調(diào)整下一時(shí)刻的平滑因子和加速因子

式中，i＝1，2，…，Ns；j＝1，2，…，d。若在k*時(shí)刻，所有參數(shù)的RMAD都小于對應(yīng)的高限閾值，進(jìn)入下一階段。

式中，i＝1，2，…，Ns；j＝1，2，…，d。

上述參數(shù)自適應(yīng)調(diào)節(jié)策略的思想是：在第1階段，采用一個(gè)比較粗獷的方差閾值（即高限閾值）作為參考，運(yùn)用IKS-APF方法估計(jì)參數(shù)并使其快速收斂；一旦某個(gè)參數(shù)方差滿足閾值條件，保持其相應(yīng)的超參數(shù)不變，當(dāng)所有的參數(shù)都滿足閾值條件時(shí)，可進(jìn)入第2階段。第1階段調(diào)節(jié)是為第2階段的精細(xì)調(diào)節(jié)做準(zhǔn)備，一般在所有參數(shù)滿足高限閾值要求的條件下，短期預(yù)測誤差滿足條件有保證。在第2階段，在參數(shù)估計(jì)準(zhǔn)確的基礎(chǔ)上，盡可能地縮小參數(shù)的方差，減少其不確定性。如果短期預(yù)測誤差滿足要求，某個(gè)參數(shù)方差不滿足低限閾值要求，繼續(xù)對其進(jìn)行核縮；若兩個(gè)要求都滿足，停止核縮，保持現(xiàn)有超參數(shù)轉(zhuǎn)為普通的隨機(jī)游走；如果預(yù)測誤差不滿足要求，則將所有參數(shù)的加速因子設(shè)回1，增益因子和平滑因子調(diào)為設(shè)定的值，緩慢地增大各個(gè)參數(shù)的超參數(shù)，直至再次收斂。最后一種情況出現(xiàn)的原因有兩個(gè)：核縮過度或者系統(tǒng)中的某些參數(shù)出現(xiàn)幅度較大的變化，兩種情況都使得參數(shù)粒子集覆蓋不到參數(shù)的真實(shí)值。由于不能確定哪些參數(shù)的變化導(dǎo)致預(yù)測誤差的增加，因此需要對所有參數(shù)同時(shí)進(jìn)行超參數(shù)放大，重新達(dá)到新的平衡。由此可見，該調(diào)節(jié)方法具有魯棒性。各個(gè)參數(shù)的RMAD高低限閾值根據(jù)參數(shù)的動態(tài)特性來定，變化或者波動的幅度越大，對應(yīng)的閾值也越高。另外，兩個(gè)參數(shù)調(diào)節(jié)階段是順序執(zhí)行的，一旦進(jìn)入階段2，不用再返回到階段1。

3.2.2 短期估計(jì)（λ－p*＋1≤k≤λ）

從λ－p*＋1時(shí)刻起，不論短期失效閾值和RMAD是否滿足條件，保持所有參數(shù)的超參數(shù)不變，γ和q的值均為1，持續(xù)進(jìn)行狀態(tài)和參數(shù)的更新。

3.3 第3階段（k＞λ）：迭代預(yù)測與RUL計(jì)算

基于當(dāng)前λ時(shí)刻狀態(tài)和參數(shù)的聯(lián)合估計(jì)＾p（xλ，θλ|z1∶λ）迭代預(yù)測p次，可得到λ＋p時(shí)刻的狀態(tài)和參數(shù)的粒子分布

圖2 二維狀態(tài)下系統(tǒng)狀態(tài)粒子軌跡演化示意圖

4 實(shí)驗(yàn)與分析

本文使用飛機(jī)壁板的裂紋增長仿真實(shí)驗(yàn)來驗(yàn)證所提方法的有效性。裂紋增長和疲勞載荷關(guān)系密切，一般使用Paris模型［15］來描述裂紋增長與循環(huán)載荷周期之間的關(guān)系，該模型經(jīng)常被用于裂紋的增長預(yù)測［7，1617］。模型如下：

式中，a是裂紋尺寸；N是周期數(shù)；C和m是與材料和環(huán)境相關(guān)的損傷模型參數(shù)；ΔK是應(yīng)力強(qiáng)度因子幅；Δσ是應(yīng)力幅。狀態(tài)方程為

裂紋尺寸可直接測得，因此觀測方程為

假設(shè)Δσ＝78 MPa，每ΔN＝30個(gè)周期，根據(jù)觀測值zk估計(jì)一次模型參數(shù)Ck和mk。由于制造工藝和材料屬性等的不同，飛機(jī)壁板個(gè)體之間損傷模型參數(shù)會存在差異，但是可以通過實(shí)驗(yàn)獲得這些參數(shù)值范圍的上下界［18］。針對7075-T651型號的鋁合金結(jié)構(gòu)材料，本文假設(shè)已知的C和m的初始分布分別服從（5×10－11，5×10－10）和（3，6）的均勻分布。設(shè)定真實(shí)參數(shù)值為ln Ctrue＝ln（1.5×10－10）＝－22.620 4，mtrue＝3.8，在仿真過程中分別添加方差為0.04和0.01的高斯白噪聲。能觀測到的初始裂紋長度是0.01 m，過程噪聲na是標(biāo)準(zhǔn)方差為3×10－4的高斯噪聲。首先依靠式（22）中的損傷增長方程得到真實(shí)的裂紋增長數(shù)據(jù)，再添加標(biāo)準(zhǔn)方差為0．001 m的高斯噪聲va，得到觀測數(shù)據(jù)zk。設(shè)定失效閾值為0.046 3 m，實(shí)際裂紋增長曲線及其觀測值如圖3所示，真實(shí)失效時(shí)間為第88ΔN。

圖3 實(shí)際裂紋增長曲線及其觀測值

首先根據(jù)觀測數(shù)據(jù)序列，分別使用文獻(xiàn)［10- 11］方法、文獻(xiàn)［12］方法和本文提出的IKS-APF方法對系統(tǒng)的參數(shù)Ct和mt進(jìn)行估計(jì)。其中，文獻(xiàn)［10- 11］方法的短期預(yù)測步數(shù)設(shè)為5，誤差閾值為0.1，方差增大系數(shù)為1.02，縮小系數(shù)為0.95；文獻(xiàn)［12］方法采用二階段控制，v1*＝［3%，0．6%］，v2*＝［1%，0．2%］，T1＝［10%，0］，T2＝［3%，0］，經(jīng)過反復(fù)調(diào)整，比例增益因子P1＝P2＝［2×10－2，5×10－2］時(shí)效果較理想；IKS-APF方法的短期預(yù)測步數(shù)同樣設(shè)為5，誤差閾值為0．1，兩個(gè)參數(shù)的RMAD高低限閾值分別與文獻(xiàn)［12］方法的T1和T2一致，另外兩個(gè)參數(shù)的初始平滑因子均為h0＝0.2，初始增益因子均為γ0＝2，兩個(gè)參數(shù)的加速因子也設(shè)為相同，但是分q0＝1（不加速）和q0＝0.95（加速）兩種情況來實(shí)驗(yàn)。前兩種方法濾波都采用SIR算法。另外4種方法粒子數(shù)目均設(shè)為2 000，初始時(shí)刻ξC與ξm的超參數(shù)ΣC和Σm均為h20與該時(shí)刻粒子方差的乘積，它們的參數(shù)估計(jì)結(jié)果如圖4所示（參數(shù)C用對數(shù)形式表示）。

圖4 4種方法參數(shù)估計(jì)效果對比

由圖4（a）可以看出，文獻(xiàn)［10- 11］方法估計(jì)得到的參數(shù)粒子方差不但沒有減小，反而逐步增大。因?yàn)槠鹗紩r(shí)刻參數(shù)粒子并未收斂，短期預(yù)測誤差很難滿足誤差要求，調(diào)節(jié)程序進(jìn)一步擴(kuò)大方差，進(jìn)入惡性循環(huán)，這種情況下并不適用。圖4（b）表明在比例增益因子設(shè)置合理時(shí)文獻(xiàn)［12］方法能夠有效減小參數(shù)mt的方差，但是收斂速度并不快，另外對參數(shù)Ct無明顯效果。究其原因，根據(jù)第2.2節(jié)分析，該方法雖然可以減小超參數(shù)的值，但隨機(jī)游走模型得到的方差總體上是變大的，該方法減小粒子方差主要還是依靠PF的濾波能力（給靠近真實(shí)值的粒子更大權(quán)值），因此方差減小的速度是有限的，特別對于觀測值影響不是很大的參數(shù)（例如Ct）很難起作用。由圖4（c）可知，IKS-APF方法在不加速的條件下參數(shù)估計(jì)的效果和文獻(xiàn)［12］方法比較接近。圖4（d）是該方法在加速條件下的參數(shù)估計(jì)結(jié)果，相比前面的3種方法，它能夠更快地將參數(shù)方差收斂到一個(gè)較低的水平，當(dāng)達(dá)到對應(yīng)的RMAD低限時(shí)，動態(tài)地維持穩(wěn)定。

在k＝55ΔN時(shí)，分別使用文獻(xiàn)［12］方法和IKS-APF方法（加速）來預(yù)測未來裂紋增長軌跡與RUL的概率分布，結(jié)果分別如圖5和圖6所示。

圖5 k＝55ΔN時(shí)文獻(xiàn)［12］方法裂紋預(yù)測結(jié)果

圖6 k＝55ΔN時(shí)IKS-APF方法（加速）裂紋預(yù)測結(jié)果

由圖5可以看出，隨著預(yù)測的推進(jìn)，文獻(xiàn)［12］方法預(yù)測得到的裂紋粒子增長軌跡越來越發(fā)散，預(yù)測不確定性迅速增大。雖然各時(shí)刻粒子的期望值與裂紋真實(shí)值比較接近（由于是指數(shù)增長關(guān)系，一個(gè)大于閾值的粒子可以平衡多個(gè)小于閾值的粒子），但是此時(shí)的RUL預(yù)測期望為ERUL＝43．6ΔN，與真實(shí)壽命33ΔN有一定距離，且RUL的概率分布范圍較寬，預(yù)測的置信度并不高。

由圖6（a）可知，IKS-APF（加速）的預(yù)測期望軌跡與真實(shí)軌跡同樣吻合較好，與文獻(xiàn)［12］方法結(jié)果的一個(gè)顯著不同是：由于估計(jì)得到的模型參數(shù)具有更小的方差，在預(yù)測推進(jìn)過程中它的裂紋粒子增長軌跡發(fā)散的速度要慢很多。此時(shí)，RUL預(yù)測期望為39.2ΔN，更接近真實(shí)壽命33ΔN。得到的RUL概率分布如圖6（b）所示，分布較陡峭，預(yù)測置信度較高。

為了更好地進(jìn)行比較，分別在k＝45，50，60，65，70時(shí)重復(fù)上述實(shí)驗(yàn)，并采用文獻(xiàn)［19］中的PH指標(biāo)、α－λ指標(biāo)以及RA指標(biāo)來進(jìn)行預(yù)測性能評估，各指標(biāo)定義如下。

PH指標(biāo)值表示預(yù)測結(jié)果首次滿足指定的預(yù)測準(zhǔn)確度要求的時(shí)間點(diǎn)與設(shè)備失效時(shí)間點(diǎn)之間的時(shí)間差值。α－λ指標(biāo)用來判斷λ時(shí)刻是否預(yù)測準(zhǔn)確度在此時(shí)真實(shí)RUL的α× 100%的范圍內(nèi)。這里將兩個(gè)指標(biāo)相結(jié)合，有

式中，tTTF為真實(shí)的失效時(shí)間

表示預(yù)測結(jié)果首次以大于β的概率位于真實(shí)RUL±α準(zhǔn)確度區(qū)域的時(shí)間索引。其中，l為所有執(zhí)行預(yù)測的時(shí)間索引。π［r（λ）］為預(yù)測得到的RUL概率分布在±α區(qū)域內(nèi)的概率質(zhì)量和，可用于度量預(yù)測的精度，α和β可根據(jù)情況設(shè)定，這里取α＝0.1，β＝0.4。

RA指標(biāo)用于表征預(yù)測的相對準(zhǔn)確度，即

式中，TRUL（λ）為λ時(shí)刻系統(tǒng)的真實(shí)RUL；ERUL（λ）為λ時(shí)刻的預(yù)測RUL期望。

當(dāng)預(yù)測所得PH值、π值和RA值越大，預(yù)測方法越優(yōu)越，詳細(xì)比較如圖7和表1所示。

圖7 IKS-APF方法與文獻(xiàn)［12］方法預(yù)測性能指標(biāo)比較

表1 IKS-APF方法與文獻(xiàn)［12］方法預(yù)測性能指標(biāo)對比

由圖7（a）可知，文獻(xiàn)［12］方法隨著預(yù)測起始點(diǎn)的增大，RUL預(yù)測逐步靠近真實(shí)值，但是直到第70ΔN，文獻(xiàn)［12］的預(yù)測結(jié)果都沒有進(jìn)入±α的準(zhǔn)確區(qū)域。圖7（b）中，IKSAPF方法（加速）隨著預(yù)測起始點(diǎn)的增大，預(yù)測準(zhǔn)確度提高，并在第60ΔN時(shí)以概率π＝0．489（＞β）首次進(jìn)入±α的準(zhǔn)確區(qū)域，所以PH＝28ΔN。圖7的一個(gè)顯著特征是IKSAPF方法（加速）獲得的RUL分布與文獻(xiàn)［12］方法得到的分布相比更加集中在一段較小的區(qū)域內(nèi)，分布曲線更加陡峭，這反映出本文方法的預(yù)測結(jié)果置信度水平更高，表1中兩種方法π值的比較也證實(shí)了這一點(diǎn)。由表1可知，IKSAPF方法（加速）的π值和RA值要全面優(yōu)于文獻(xiàn)［12］方法，該方法估計(jì)得到的模型參數(shù)具有更小的不確定性，保證了其預(yù)測結(jié)果具有更好的準(zhǔn)確度和置信度。

5 結(jié) 論

本文提出了一種基于IKS-APF的失效預(yù)測方法，該方法可用于解決系統(tǒng)模型中存在多個(gè)未知參數(shù)情況下的失效預(yù)測問題，并具有以下優(yōu)點(diǎn)：①在原有核平滑輔助粒子濾波基礎(chǔ)上引入了增益因子和加速因子，參數(shù)估計(jì)的魯棒性更好，收斂速度更快；②具有參數(shù)自適應(yīng)方差調(diào)節(jié)能力，通過短期預(yù)測誤差來度量參數(shù)估計(jì)的準(zhǔn)確度，同時(shí)監(jiān)控參數(shù)粒子方差來調(diào)整參數(shù)估計(jì)的精度，融合了文獻(xiàn)［10- 11］方法和文獻(xiàn)［12］方法兩者的優(yōu)點(diǎn)，有利于減少預(yù)測的不確定性；③該方法只需知道參數(shù)的大致范圍，并能夠在一定程度上對多個(gè)參數(shù)進(jìn)行獨(dú)立的調(diào)節(jié)。裂紋增長仿真實(shí)驗(yàn)結(jié)果表明了本文方法比已有方法具有更好的預(yù)測準(zhǔn)確度和置信度。

［1］Sun B，Kang R，Xie J S．Research and application of the prognostic and health management system［J］．Systems Engineering and Electronics，2007，29（10）：1762- 1767．（孫博，康銳，謝勁松．故障預(yù)測與健康管理系統(tǒng)研究和應(yīng)用現(xiàn)狀綜［J］．系統(tǒng)工程與電子技術(shù)，2007，29（10）：1762- 1767．）

［2］Jardine A K S，Lin D，Banjevic D．A review on machinery diagnostics and prognostics implementing condition-based maintenance［J］．Mechanical Systems and Signal Processing，2006，20（7）：1483- 1510．

［3］Orchard M，Wu B，Vachtsevanos G．A particle filtering framework for failure prognosis［C］∥Proc.of the World Tribology Congress III，2005.

［4］Saha B，Goebel K，Poll S，et al．Prognostics methods for battery health monitoring using a Bayesian framework［J］．IEEE Trans.on Instrumentation and Measurement，2009，58（2）：291- 296．

［5］Zhang Q，Hu C H，Qiao Y K，et al．Fault prediction algorithm based on weight selected particle filter［J］．Systems Engineering and Electronics，2009，31（1）：221- 224．（張琪，胡昌華，喬玉坤，等．基于權(quán)值選優(yōu)粒子濾波器的故障預(yù)測算法［J］．系統(tǒng)工程與電子技術(shù)，2009，31（1）：221- 224．）

［6］Zhang L，Li X S，Yu J S，et al．A fault prognostic algorithm based on Gaussian mixture model particle filter［J］．Acta Aeronautica et Astronautica Sinica，2009，30（2）：319- 324．（張磊，李行善，于勁松，等．一種基于高斯混合模型粒子濾波的故障預(yù)測算法［J］．航空學(xué)報(bào)，2009，30（2）：319- 324．）

［7］Zio E，Peloni G．Particle filtering prognostic estimation of the remaining useful life of nonlinear components［J］．Reliability Engineering and System Safety，2011，96（3）：403- 409．

［8］Chen X Z，Yu J S，Tang D Y，et al．A novel PF-LSSVR-based framework for failure prognosis of nonlinear systems with timevarying parameters［J］．Chinese Journal of Aeronautics，2012， 25（5）：715- 724．

［9］Miao Q，Xie L，Cui H J，et al．Remaining useful life prediction of lithium-ion battery with unscented particle filter technique［J］．Microelectronics Reliability，2013，53（6）：805- 810．

［10］Orchard M，Kacprzynski G，Goebel K，et al．Advances in uncertainty representation and management for particle filtering applied to prognostics［C］∥Proc.of the International Conference on Prognostics and Health Management，2008：1- 6.

［11］Orchard M，Tobar F，Vachtsevanos G．Outer feedback correction loops in particle filtering-based prognostic algorithms：statistical performance comparison［J］．Studies in Informatics and Control，2009，18（4）：295- 304.

［12］Orchard M J，Goebel K．Model-based prognostics with concurrent damage progression processes［J］．IEEE Trans.on Systems，Man，and Cybernetics-Part A：Systems and Humans，2013，43（3）：535- 546.

［13］Liu J，West M．Combined parameter and state estimation in simulation-based filtering［M］∥Doucet A，de Freitas J F G，Gordon N J．Sequential Monte Carlo methods.New York：Springer，2001：197- 223．

［14］Pitt M K，Shephard N．Filtering via simulation：auxiliary particle filters［J］．Journal of the American Statistical Association，1999，94（446）：590- 599.

［15］Paris P C，Erdogan F．A critical analysis of crack propagation laws［J］．Journal of Basic Engineering，1963，85（4）：528- 533.

［16］Zhang B，Sconyers C，Orchard M，et al．Fault progression modeling：an application to bearing diagnosis and prognosis［C］∥Proc.of the American Control Conference，2010：6993- 6998．

［17］An D，Choi J H，Kim N H．Prognostics 101：a tutorial for particle filter-based prognostics algorithm using matlab［J］．Reliability Engineering and System Safety，2013，115：161- 169．

［18］Newman J C，Phillips E P，Swain M H．Fatigue-life prediction methodology using small-crack theory［J］．International Journal of Fatigue，2009，21（2）：109- 119．

［19］Saxena A，Celaya J，Saha B，et al．On applying the prognostic performance metrics［C］∥Proc.of the Annual Conference of Prognostics and Health Management Society，2009：1- 16.

Failure prognostics using improved kernel smoothing auxiliary particle filtering

CHEN Xiong-zi1，YU Jin-song1，2，TANG Di-yin1，LI Xing-shan1
（1.School of Automation Science and Electrical Engineering，Beihang University，Beijing 100191，China；2.Co-Innovation Center for Advanced Aero-Engine，Beijing 100191，China）

A failure prognostics approach based on improved kernel smoothing auxiliary particle filtering（IKS-APF）is proposed for systems with multiple unknown parameters．Firstly，a gain factor and an acceleration factor are employed in the kernel smoothing APF to bi-directionally control the parameter variance and accelerate the parameter convergence．Secondly，the IKS-APF method is used to jointly estimate the states and parameters．In order to ensure the accuracy of parameter estimation and reduce its uncertainty，an adaptive control scheme for the particle variance of parameters is presented，combining the variance monitoring and the short-term prediction errors．Finally，iterative prediction is implemented based on the latest estimated state and parameter particles，and then the remaining useful life（RUL）is calculated by the time of each propagated state particle first entering the failure zone．Simulation results demonstrate the effectiveness and superiority of the proposed approach．

failure prognostics；kernel smoothing；auxiliary particle filtering；parameter estimation；variance control；uncertainty management

TP 206

A

10．3969／j．issn．1001-506X．2015．01．17

陳雄姿（1985-），男，博士研究生，主要研究方向?yàn)轭A(yù)測與健康管理技術(shù)、設(shè)備失效預(yù)測與剩余使用壽命估計(jì)。

E-mail：cxzbuaa＠163．com

于勁松（1968-），通訊作者，男，副教授，博士，主要研究方向?yàn)轭A(yù)測與健康管理技術(shù)、自動測試技術(shù)。

E-mail：yujs＠buaa．edu．cn

唐荻音（1986-），女，博士研究生，主要研究方向?yàn)轭A(yù)測與健康管理技術(shù)、設(shè)備退化建模與視情維修策略。

E-mail：amytdy＠asee．buaa．edu．cn

李行善（1938-），男，教授，博士研究生導(dǎo)師，主要研究方向?yàn)樽詣訙y試系統(tǒng)、虛擬儀器技術(shù)、預(yù)測與健康管理技術(shù)。

E-mail：lixingshan＠263．net

1001-506X（2015）01-0101-08

網(wǎng)址：www．sys-ele．com

2014- 04- 25；

2014- 06- 19；網(wǎng)絡(luò)優(yōu)先出版日期：2014- 08- 20。

網(wǎng)絡(luò)優(yōu)先出版地址：http：／／w ww．cnki．net／kcms／detail／11．2422．TN．20140820．1733．006．html

航空科學(xué)基金（20100751010）；北京市自然科學(xué)基金（4113073）資助課題