基于改進核平滑輔助粒子濾波的失效預測方法

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

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

基于改進核平滑輔助粒子濾波的失效預測方法

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

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

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

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

0 引 言

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

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

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

1 問題描述

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

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

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

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

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

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

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

2 改進的核平滑輔助粒子濾波

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

2.1 輔助粒子濾波

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

步驟5 歸一化權值wjk；

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

3 IKS-APF預測及其不確定性管理

基于IKS-APF的失效預測框架主要包括3個階段。假設當前時刻為λ，已知觀測序列z1∶λ＝｛z1，z2，…，zλ｝。

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

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

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

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

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

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

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

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

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

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

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

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

3.2.2 短期估計（λ－p*＋1≤k≤λ）

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

3.3 第3階段（k＞λ）：迭代預測與RUL計算

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

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

4 實驗與分析

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

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

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

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

圖3 實際裂紋增長曲線及其觀測值

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

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

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

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

圖5 k＝55ΔN時文獻［12］方法裂紋預測結果

圖6 k＝55ΔN時IKS-APF方法（加速）裂紋預測結果

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

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

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

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

式中，tTTF為真實的失效時間

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

RA指標用于表征預測的相對準確度，即

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

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

圖7 IKS-APF方法與文獻［12］方法預測性能指標比較

表1 IKS-APF方法與文獻［12］方法預測性能指標對比

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

5 結 論

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

［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．（孫博，康銳，謝勁松．故障預測與健康管理系統(tǒng)研究和應用現(xiàn)狀綜［J］．系統(tǒng)工程與電子技術，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．（張琪，胡昌華，喬玉坤，等．基于權值選優(yōu)粒子濾波器的故障預測算法［J］．系統(tǒng)工程與電子技術，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．（張磊，李行善，于勁松，等．一種基于高斯混合模型粒子濾波的故障預測算法［J］．航空學報，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-），男，博士研究生，主要研究方向為預測與健康管理技術、設備失效預測與剩余使用壽命估計。

E-mail：cxzbuaa＠163．com

于勁松（1968-），通訊作者，男，副教授，博士，主要研究方向為預測與健康管理技術、自動測試技術。

E-mail：yujs＠buaa．edu．cn

唐荻音（1986-），女，博士研究生，主要研究方向為預測與健康管理技術、設備退化建模與視情維修策略。

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

李行善（1938-），男，教授，博士研究生導師，主要研究方向為自動測試系統(tǒng)、虛擬儀器技術、預測與健康管理技術。

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)絡優(yōu)先出版日期：2014- 08- 20。

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

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