裝備剩余壽命預測平行仿真中模型動態演化方法*

葛承壟，朱元昌，邸彥強，崔浩浩

(陸軍工程大學石家莊校區， 河北 石家莊 050003)

裝備平行仿真首次完整提出是在2016年亞仿/秋季仿真大會上[1-2]。裝備平行仿真旨在將實際裝備和仿真系統結合在一起，仿真系統利用傳感器采集的裝備信息演化仿真模型，實際裝備受益于仿真系統的仿真結果，從而提高裝備的運用效能。以這種模式運行的仿真系統稱之為平行仿真系統。

裝備平行仿真中的一個重要概念是實時數據驅動的模型演化，即在裝備實時數據驅動下，仿真模型形態或參數進行適應性調整的過程，包括模型適宜性選擇和模型參數演化兩個方面的內涵，這被認為是其區別于以往仿真技術的典型特征。以往仿真技術中，仿真模型側重于一次性構建，仿真運行后模型形態和模型參數不再發生變化，并不存在模型演化過程。在裝備平行仿真中，模型演化的目的是為了滿足動態變化的仿真需求，通過基于實時裝備數據的動態仿真，動態提高仿真逼真度和仿真預測準確度，為仿真預測和輔助決策提供數據支持。仿真模型及其演化屬于裝備平行仿真的模型理論范疇，也是裝備平行仿真研究中的基礎問題。裝備平行仿真還具有虛實共生、平行運行、數據驅動等技術特征。目前，裝備平行仿真主要包括兩個研究分支：一是面向戰場指揮決策領域的平行仿真，代表性的研究學者包括邱曉剛[3]、周芳[4]、竇林濤[5]等，此類平行仿真已初步建立演化建模框架，包括模型動態匹配、模型類別修正、模型參數校準和實體模型行為意圖更新等，相關演化方法的研究正在逐步展開；二是裝備維修保障領域中面向裝備剩余壽命(Remaining Useful Life，RUL)預測的平行仿真，已完成模型構建機理研究和模型演化方法設計。本文關注的是面向裝備RUL預測的平行仿真中模型演化方法研究。

在裝備維修保障領域中，由于外部沖擊、磨損、疲勞、腐蝕等原因，裝備的性能將不可避免地發生退化，最終引起裝備故障甚至造成嚴重事故，因此需要視情對裝備進行健康維護。近年來，預測與健康管理技術是實現裝備健康維護的主要方法[6]，它通過預測裝備的RUL，對裝備進行合理的維修與管理，保證設備運行的安全性、可靠性與經濟性，其關鍵內容包括RUL預測與健康管理兩個方面。其中，RUL預測是預測與健康管理技術的核心內容[7]。RUL預測主要是指解算RUL分布，即概率密度函數(Probability Density Function，PDF)，它表征了壽命預測的不確定性，是維修決策的重要依據。

當前，RUL預測方法可分為失效物理模型法、基于統計的方法和人工智能法[8]。對于復雜裝備來說，其失效機理很難獲得，因此后兩種方法得到更多關注。基于統計的方法和人工智能法分別通過統計模型、機器學習進行數據擬合，從而預測RUL。然而基于統計的方法和人工智能法仍存在以下典型不足：第一，在預測過程中通常假定裝備退化模式在整個預測周期中是固定不變的，利用單個模型進行預測，然而退化模式卻可能呈現多種退化模式混合的情況，其中連續退化與未知離散沖擊混合是一種典型情形；第二，模型參數往往通過極大似然估計等離線方法獲得，不能隨著退化數據的積累而在線更新；第三，難以得到RUL概率密度函數的解析表達式，制約預測方法的實時性，尤其是人工智能法，它還需要大量訓練樣本，這在實際中往往無法滿足。因此，針對現實中普遍存在的連續退化與未知離散沖擊混合的退化過程，迫切需要在RUL預測中同時考慮模型適宜性選擇和參數在線演化。依據其技術原理和典型特征，裝備平行仿真為解決此類RUL預測問題提供了新思路。本文研究涉及面向混合退化裝備RUL預測的平行仿真中仿真模型構建、模型動態演化以及基于平行仿真的RUL實時預測等問題。

1 平行仿真建模

1.1 建模分析

文獻[1,9]指出，構建裝備性能退化狀態空間模型(State Space Model，SSM)是面向裝備RUL預測平行仿真的建模方向。退化狀態估計是RUL預測的前提，裝備性能退化SSM兼顧了退化過程的動態性和時變性，易于進行退化狀態估計，有利于RUL預測。特別地，考慮到退化狀態是經平行仿真得到的，故可稱之為仿真退化狀態。

由于基于統計的方法更容易得到RUL概率密度函數的解析表達式，本文將隨機過程方法與SSM結合起來構建平行仿真模型。隨機過程適宜描述退化過程的隨機性和不確定性，Wiener過程和Gamma過程是最常用的隨機過程，但Gamma過程適合于描述具有嚴格單調特征的退化過程，適用條件過于苛刻，而Wiener過程具有適用范圍廣、首達時(First Hitting Time，FHT)分布明確等建模優勢。因此，在平行仿真建模中將其與SSM建模法結合起來，構建Wiener狀態空間模型(Wiener SSM，WSSM)是一種合理選擇[9]。特別地，為描述帶未知離散沖擊的混合退化過程，宜建立多態WSSM作為平行仿真模型，它包括連續退化模型和含未知離散沖擊的退化模型兩種形態。

1.2 多態WSSM構建

多態WSSM構建涉及裝備混合退化狀態方程和觀測方程。裝備混合退化狀態方程包括兩種形態，一種是連續退化狀態方程，另一種是帶未知離散沖擊的退化狀態方程。首先，利用Wiener過程構建連續退化狀態方程，有

x(t)=x(0)+ηt+σB(t)

(1)

式中：{x(t),t≥0}是由標準Brownian運動{B(t),t≥0}驅動的連續退化過程，且有B(t)～N(0,t)；x(0)為初始退化狀態；η、σ分別是標準Brownian運動的漂移系數和擴散系數。對式(1)進行Euler離散化，得到在離散時間點tk(k=1,2,…)上的不考慮離散沖擊的狀態方程

(2)

隨機離散沖擊的未知特性是指沖擊時刻未知，其對裝備性能造成的損傷可用一個沖擊變量表示。根據沖擊的離散特性，將沖擊造成的損傷融入式(2)中，得到帶未知離散沖擊的退化狀態方程，有

(3)

式中，D為沖擊造成的損傷。假設沖擊的到達服從Poisson過程[10]，即有

(4)

式中，ρ為沖擊到達率，M(tk)表示從初始時刻至k時刻沖擊出現的總數量，n的取值為0或1[11]。為簡便起見，這里假定Poisson過程與Brownian運動相互獨立。

裝備退化觀測數據y(t)與x(t)的隨機關系可由觀測方程描述，即

y(t)=x(t)+π(t)

(5)

式中，π(t)～N(0,φ2)，φ2為測量噪聲的方差，并假設π(t)與Brownian運動B(t)相互獨立。對式(5)進行Euler離散化，得到在離散時間點tk(k=1,2,…)上的觀測方程為

yk=xk+φζk

(6)

根據式(1)～(6)，可以得到在離散時間點tk(k=1,2,…)上的多態WSSM為

(7)

2 多態WSSM演化

2.1 模型演化分析

在SSM建模方式下，仿真退化狀態可通過實際退化觀測值驅動多態WSSM運行而獲得，仿真退化數據與實際退化觀測值構成完全退化數據。在SSM建模框架下，仿真退化狀態的實質是一種隱含退化狀態。為實現多態WSSM動態演化，提出如下演化機制：一是基于交互多模型(Interactive Multiple Model，IMM)濾波[12]的模型軟切換，即利用觀測數據修正仿真預測結果，并動態計算不同仿真模型的概率，實現仿真輸出和實際觀測數據的同化，得到仿真退化狀態估計；二是在實時退化觀測數據驅動下在線估計模型參數。模型軟切換和參數在線估計不是孤立執行的，而是相互迭代，通過二者的不斷迭代，動態校正模型預測輸出，從而動態提高仿真預測準確度。

在多態WSSM中，沖擊的未知特性使得無法通過觀測數據得知是否有沖擊到達，這就導致無法得知仿真模型切換的判斷條件，因此需要利用IMM濾波計算不同模型形態的概率。由于裝備性能退化過程受離散沖擊影響，具有退化軌跡突變、非平穩的特點，考慮在IMM濾波中采用強跟蹤濾波器[13]。基于IMM強跟蹤濾波(IMM Strong Tracking Filter，IMMSTF)的模型軟切換包括5個階段，即模型輸入交互、模型預測、濾波、模型概率計算和模型輸出交互。多態WSSM是一種含有隱含狀態的狀態空間模型，故考慮利用期望最大化(Expectation Maximum，EM)算法[14]實現模型參數演化，EM算法能有效估計此類SSM的模型參數。綜合以上分析，裝備平行仿真中模型動態演化示意如圖1所示。

圖1 裝備平行仿真中模型動態演化Fig.1 Model dynamic evolution of equipment parallel simulation

2.2 基于IMMSTF的模型軟切換

(8)

(9)

(10)

(11)

(12)

(13)

(14)

則仿真模型u的概率為

(15)

(16)

(17)

2.3 模型參數演化

在多態WSSM中，模型參數θ包括μ0、Σ0、η、σ、D和φ，其中μ0、Σ0表示初始退化狀態x0的均值和協方差。根據EM算法，在監測時刻k、EM算法第j步時，模型參數θ={μ0,Σ0,η,σ,D,φ}的估計可由式(18)獲得。

(18)

(19)

2.3.1 計算聯合對數似然函數的數學期望

(20)

根據多態WSSM可知

p(yi|xi,θ)=N(yi;xi,φ2)

p(x0|θ)=N(x0;μ0,Σ0)

(21)

(22)

IMMBS算法包括后向濾波和模型狀態融合兩個階段。后向濾波是指從最新的測量值開始后向濾波，其運算流程與IMMSTF類似，但是二者也存在明顯差異。IMMSTF先執行輸入交互再執行一步預測，而后向濾波則是先執行一步預測再執行輸入交互。后向濾波可以分為后向一步預測、后向輸入交互、后向濾波更新、后向模型概率計算和后向輸出融合5個步驟。特別地，后向濾波的后向一步預測方程為：

(23)

(24)

(25)

(26)

(27)

(28)

2.3.2 最大化聯合對數似然函數的數學期望

時刻k處、期望最大化算法第j步時，θ的估計可利用對聯合對數似然函數的數學期望取偏導數并令偏導數為0求得，即

(29)

求解可得參數θ的在線估計值，即

(30)

3 剩余壽命分布實時解算

根據FHT的概念，裝備RULT定義為x(t)首次通過失效閾值w的時間[17]，即

T(w)=inf{t:x(t)≥w|x(0)

(31)

鑒于多態維納狀態空間模型在某特定時刻可能存在式(2)、式(3)兩種仿真退化狀態方程，無法直接得到壽命分布的PDF，需對x(t)進行轉化。根據文獻[11]可知，維納過程可改寫為

x(t)=xk+η(t-tk)+σB(t-tk)

(32)

其中，B代表布朗運動。將式(2)、式(3)進行融合，即將n次Poisson沖擊造成的損傷融入式(32)中，有

x(t)=xk+nD+η(t-tk)+σB(t-tk)

(33)

根據式(31)中RUL的定義和維納過程的FHT分布性質可知，多態WSSM中以n次Poisson沖擊、xk、θ和Yk為條件的RUL概率密度函數服從逆高斯分布，即

(34)

其中，Tk為監測時刻k處的剩余壽命，其概率密度函數可以寫為

f(Tk|n,xk,θ,Yk)

(35)

引理1 設Ω～N(γ,ξ2)，A、B、C為常數，則有

(36)

引理的證明見文獻[18]。

定理1對于混合退化過程{x(t),t≥0}，時刻k處裝備剩余壽命概率密度函數為

(37)

fT(Tk|θ,Yk)

□

根據數學期望的定義，平行仿真系統利用數值積分計算RUL的數學期望值，即

(38)

根據k時刻RUL概率密度函數和數學期望的解析表達式，平行仿真系統能夠以在線、實時的方式計算RUL的概率密度函數和期望值，為裝備維修決策提供數據支撐。

4 實例研究

4.1 數據介紹

機械裝備中軸承的性能退化是典型的含有沖擊特性的混合退化過程，本文采用某軸承全壽命試驗數據進行模型演化方法驗證。該數據由法國FEMTO-ST研究所提供[19]，全壽命試驗在PRONOSTIA平臺上進行，近年來被廣泛應用于可靠性領域的方法驗證。試驗過程工況條件以及試驗數據相關信息詳見文獻[11]，本文仍以軸承1_3為例進行實例研究。試驗過程中采集的是軸承1_3的振動加速度信號，該信號的均方根(Root Mean Square，RMS)值是常用的退化特征量，計算公式為

(39)

式中，N是采樣點數，這里N=2560，xi為第i個采樣點對應的振動加速度信號。軸承1_3的RMS如圖2所示。根據圖2可知，軸承1_3的性能退化過程沖擊特性明顯，適宜用于驗證本文方法。根據文獻[11]可知，失效閾值設為4.714 5，即w=4.714 5。

圖2 軸承1_3的均方根值Fig.2 RMS of bearing 1_3

4.2 多態WSSM動態演化與RUL實時預測

多態WSSM參數初始設置為x0=0.2、η=0.02、σ=0.5、D=0.02、ρ=0.5、τ=1、φ=0.1，Markov模型轉移概率矩陣P=[0.50.5;0.60.4]，模型初始概率分別為0.6和0.4，即μ0=[0.60.4]T。得到的退化軌跡對比如圖3所示，仿真退化軌跡和實際退化軌跡差異較小，表明在實時退化數據驅動下，仿真退化軌跡能有效逼近實際退化軌跡。利用均方根誤差(Root Mean Square Error，RMSE)來量化退化軌跡對比結果，其計算公式為

(40)

其中，r=842為監測時間點數目。經計算，兩種退化軌跡的RMSE僅為3.497%，充分說明本文提出的模型動態演化方法能有效實現軸承1_3性能退化過程的建模與仿真。

圖3 退化軌跡對比Fig.3 Comparison of the degraded traces

模型概率如圖4所示。由圖4可知，在t1500至t1765監測時間段內，由于Poisson沖擊特性并不顯著，此時線性退化特性較為明顯，線性退化模型(模型1)的概率明顯高于帶未知離散Poisson沖擊退化模型(模型2)的概率，線性退化模型的概率保持在0.74左右，帶未知離散Poisson沖擊退化模型的概率保持在0.26左右，此時處于主導地位的模型是線性退化模型。但隨著時間的推移，Poisson沖擊特性越發顯著，突出表現在t1766、t1827、t1877、t2130、t2234等時刻，帶未知離散Poisson沖擊退化模型的概率總體呈現動態上升的趨勢，反之，線性退化模型的概率呈現動態下降的趨勢。在退化后期，帶未知離散Poisson沖擊退化模型的概率超過線性退化模型的概率，說明帶離散Poisson沖擊退化模型更適宜描述當前的退化過程。從以上分析可以看出，本文方法能有效實現仿真模型的“軟”切換，使之滿足壽命預測對模型適宜性的需求。值得注意的是，由于模型庫中僅存在兩種模型，在模型軟切換條件下，二者模型概率之和為1，因此兩個模型概率曲線關于概率μ=0.5對稱。

圖4 模型概率Fig.4 Model probability

在軸承退化數據的動態驅動下，模型參數θ的在線估計結果如圖5所示。由圖5可知，漂移系數η在0.004附近動態波動，波動區間為[0,0.012]，尤其是在t1766、t1827、t1877等沖擊特性明顯的時刻，η的數值波動較大，反映出軸承1_3的退化速率加快；擴散系數σ能較快收斂，在沖擊特性明顯的時刻，擴散系數發生較大波動，并達到新的收斂狀態，擴散系數的收斂有利于獲得穩定的剩余壽命概率密度函數，同時擴散系數的收斂值較小，有利于獲得更為狹窄的剩余壽命概率密度函數，提高剩余壽命預測的準確性；Poisson沖擊造成的損傷D在初始時刻附近波動較大，后在區間[0.03,0.07]內動態波動，將損傷D考慮到剩余壽命預測中，能有效減少甚至避免“欠維修”的發生；測量誤差的標準差φ收斂速度較快，收斂值約為0.01，反映出測量誤差的波動逐漸穩定。

(a) 參數η在線演化(a) Online evolution of parameter η

(b) 參數σ在線演化(b) Online evolution of parameter σ

(c) 參數D在線演化(c) Online evolution of parameter D

(d) 參數φ在線演化(d) Online evolution of parameter φ圖5 多態WSSM參數在線演化Fig.5 Online evolution of polymorphic WSSM parameters

在每一個監測點，隨著模型演化的執行，平行仿真系統利用式(37)和式(38)預測軸承1_3的剩余壽命，包括剩余壽命概率密度函數及其數學期望。以100個監測時間點長度為預測間隔，可得到軸承1_3在t1600，t1700，…，t2300處共計8個監測點的剩余壽命預測結果，如圖6所示。

圖6 軸承1_3在不同監測點處的剩余壽命預測結果Fig.6 RUL prognostic results of bearing 1_3 at different monitoring time

根據圖6可知，在t1600，t1700，…，t2300等監測時間點，RUL概率密度曲線均能有效覆蓋實際RUL，并且隨著軸承1_3退化數據的動態注入，RUL概率密度函數逐漸收窄，右偏特性減弱，正態特性越發明顯，說明通過模型軟切換和參數在線估計，模型匹配度逐漸提高，RUL預測的不確定性逐漸減小。此外，RUL期望值與真實RUL之間的誤差較小，RUL期望值接近PDF峰值對應的RUL，表明PDF的不確定性較小，有利于輔助維修方案的制訂。

4.3 比較研究

圖7 不同方法剩余壽命預測對比結果Fig.7 Comparative RUL prediction results of different methods

圖8 第1900個監測時間點處對比結果Fig.8 Comparative results at the 1900th monitoring time

為量化預測結果，給出兩個評價指標：平均相對精度(Mean Relative Accuracy, MRA)和總均方差(Total Mean Square Error, TMSE)，它們的定義分別為：

(41)

(42)

經計算得到的方法評價結果如表1所示。根據表1可知，較之不考慮模型軟切換的演化方法，本文方法的MRA明顯較大，說明經過模型動態演化后剩余壽命預測的相對預測精度更高；本文方法的TMSE明顯較小，說明剩余壽命預測的誤差更小。

表1 方法評價結果

5 結論

以連續退化與未知離散沖擊混合的混合退化裝備剩余壽命預測為背景，提出一種裝備平行仿真中模型演化方法。該方法構建了多態Wiener狀態空間模型，利用交互多模型強跟蹤濾波動態計算模型概率，實現模型軟切換，利用期望最大化算法在線更新模型參數。通過模型軟切換和參數在線估計的迭代，動態提高了仿真逼真度，得到了首達時意義下的剩余壽命分布解析表達式。通過模型動態演化，有效提高了剩余壽命預測準確度，減小了預測不確定性。該方法還具有可擴展性，可為解決更為復雜條件下平行仿真在壽命預測領域中的應用問題提供模型演化框架。