基於改進粒子濾波算法的設備剩餘壽命預測方法與流程
2023-05-21 09:13:31 3
本發明涉及機電設備的故障預測領域，尤其涉及基於改進粒子濾波算法的設備剩餘壽命預測方法。
背景技術：
隨著人類文明進步和科技加速發展，諸如飛機、船舶、發電裝備、軌道交通設備等重大機電裝備的安全服役對於我國國民經濟發展和國防安全具有重大意義。然而由於實際工況複雜多變、環境惡劣，設備在長期的運行過程中會逐漸老化，剩餘壽命會逐漸下降，易產生惡性事故，造成重大的財產損失或人員傷亡；而若盲目進行維修更換則會維修過剩，造成巨大的浪費。因此，準確的壽命預測，對於保障設備安全、高效運行，提高經濟效益具有重要意義。
近年來，出現了很多可用於機電設備剩餘使用壽命預測的方法，如自回歸(AR,Auto Regressive)時序模型方法、人工神經網絡(ANN,Artificial Neural Networks)、支持向量機(SVM,Support Vector Machine)、相關向量機(RVM,Relevance Vector Machine)等基於數據驅動的預測方法，以及直接對產品具體零部件解析建模或仿真進行預測的基於物理模型方法等。
其中，粒子濾波(PF,Particle Filter)是一種基於模型與數據混合的設備壽命預測方法，相較於傳統的數據驅動預測方法，該方法能提供預測結果的不確定性表達，這對於決策者來說十分重要。另外，其基於的模型不同於一般的物理模型，而是一種經驗退化模型。該模型從數據角度出發，尋求採集數據隨時間變化的規律或尋求系統內部狀態的遞推關係，這種數學模型相比更容易獲取，適用範圍也更廣。同時，採用粒子濾波方法結合經驗退化模型實現剩餘壽命預測，對於非高斯、非線性退化過程具有良好的適應能力。目前，基於粒子濾波的壽命預測方法在鋰電池、滾動軸承、齒輪箱等得到了成功應用。
儘管如此，該類基於粒子濾波算法的壽命預測方法，目前存在兩個主要的缺陷，即粒子退化與樣本枯竭，導致設備剩餘壽命預測精度不高。
技術實現要素：
本發明要解決的技術問題是：提供一種基於改進粒子濾波算法的設備剩餘壽命預測方法，提高設備壽命預測的精度。
為解決上述問題，本發明採用的技術方案是：基於改進粒子濾波算法的設備剩餘壽命預測方法，所採用的粒子濾波算法包括重要性採樣階段和重採樣階段，其特徵在於，在重要性採樣階段引入無跡卡爾曼濾波(UKF,Unscented Kalman Filter)方法更新粒子，生成建議分布，從而抑制粒子退化；在重採樣階段增加馬爾科夫鏈蒙特卡洛(MCMC,Markov Chain Monte Carlo)步驟，抑制樣本枯竭。
進一步的，生成建議分布的步驟包括：
根據前一時刻粒子的狀態值及粒子的初始方差計算前一時刻粒子的sigma點集；
將前一時刻粒子的sigma點集中的點帶入粒子狀態方程，得到當前時刻的粒子sigma點集；
對當前時刻的粒子sigma點集加權求和，從而得到當前時刻粒子的期望；
做當前時刻粒子的期望與當前時刻粒子的sigma點集的方差，並對所得方差進行加權求和，得到當前時刻粒子的方差的時間更新值；
將當前時刻粒子的sigma點集中的點以及前一時刻粒子的sigma點集中的點帶入粒子狀態方程，從而得到當前時刻粒子的觀測值；
對當前時刻的粒子的觀測值加權求和，得到當前時刻粒子的觀測值的更新值；
基於當前時刻的粒子sigma點集、當前時刻粒子的期望、當前時刻的粒子的觀測值、當前時刻粒子的觀測值的更新值，計算得到當前時刻的粒子的方差的量測更新值；
計算當前時刻粒子的方差的修正係數；
基於當前時刻粒子的方差的時間更新值、當前時刻的粒子的方差的量測更新值、當前時刻粒子的方差的修正係數，計算得到當前時刻粒子的方差的修正值；
用當前時刻粒子的期望、當前時刻粒子的方差的修正值作為高斯分布的期望與方差，生成建議分。
進一步的，馬爾科夫鏈蒙特卡洛步驟包括：
抽取閾值u～U[0,1]，從建議分布中抽樣得到粒子的候選狀態利用粒子的候選狀態計算其接受概率A，若u≤A，則接受即：否則丟棄進行重採樣並保留重採樣的粒子即：
具體的，本發明的完整步驟包括：
(1)建立預測模型；
(2)對退化性能評估，並選取敏感特徵作為退化特徵；
(3)對選取的原始退化特徵進行預處理；
(4)在預測模型的基礎上，從退化特徵的歷史數據中獲取先驗信息；
(5)根據退化模型分布及模型中各參數的初始信息，對改進粒子濾波算法進行初始化；
(6)利用改進粒子濾波算法依次對各個時刻進行迭代估計和外推預測兩個過程，當退化狀態預測值小於或等於失效閾值時，預測結束；
(7)對所有粒子樣本的狀態進行統計，得到最終的預測結果。
進一步的，步驟(1)的預測模型中，預測對象的系統狀態更新方程與量測方程分別為：
xk＝fk(xk-1,vk-1)
zk＝hk(xk,nk)
其中，xk表示當前時刻的系統狀態值，zk表示當前時刻含有加性噪聲的量測值，vk-1為前一時刻的過程噪聲，nk為當前時刻的量測噪聲，k∈N，N粒子個數，且過程噪聲與量測噪聲均為平穩噪聲序列。
進一步的，步驟(3)對原始特徵進行平滑和單調化的預處理。
進一步的，先驗信息包括退化特徵隨變量的分布圖及預測模型中各參數的初始值。
本發明的有益效果是：本發明通過在重要性採樣階段引入UKF方法，對標準粒子濾波算法的重要密度函數進行改進，抑制粒子退化，在重採樣階段增加MCMC步驟，抑制樣本枯竭，相較於標準粒子濾波預測方法，預測精度更高，尤其能提供更加準確的早期預測。
附圖說明
圖1是本發明所改進粒子濾波算法流程圖；
圖2是先驗信息獲取流程圖；
圖3是本發明所提出的監測目標對象壽命預測流程圖；
圖4-7分別是電池A1-A4容量退化真實數據與曲線擬合；
圖8是已知A4前18組測量數據時改進粒子濾波預測結果；
圖9是已知A4前30組測量數據時改進粒子濾波預測結果
圖10是已知A4前18組測量數據時標準粒子濾波預測與改進粒子濾波預測比較結果；
圖11是已知A4前30組測量數據時標準粒子濾波預測與改進粒子濾波預測比較結果。
圖中編號：Un為失效閾值，Tn迭代更新周期，PDF為實施例的概率密度函數，PDF＇為標準粒子濾波的概率密度函數。
具體實施方式
下面從兩個部分——改進算法部分和設備剩餘壽命預測部分，對本發明進行詳細闡述。(一)改進粒子濾波算法
該改進粒子濾波算法通過在重要性採樣階段引入UKF方法，對標準粒子濾波算法的重要密度函數進行改進，抑制粒子退化，在重採樣階段增加MCMC步驟，抑制樣本枯竭。算法的流程圖如圖1所示，其具體步驟為：
(1)粒子初始化
對於時刻k＝0，i＝1,2,…,N，j＝1,2,…,2l，從先驗概率p(x0)中隨機抽取N個粒子得到並令粒子權值
計算均值：
擴維處理：
計算方差：
擴維得到：
(2)重要性採樣
(一)、對於k＝1,2,…，的時刻利用UKF更新粒子，生成建議分布。具體步驟如下：
①計算粒子sigma點集。即，根據前一時刻粒子的狀態值及粒子的初始方差計算前一時刻粒子的sigma點集公式如下：
式中，λ＝α2(n+κ)-n為一比例係數；α的大小決定了選取的樣本圍繞粒子均值的分布情況，通常取一小的正數；κ為另一比例係數，可以取為0；參數β反映關於粒子x的先驗分布，對於高斯分布，β＝2；表示矩陣平方根的第i行。
②時間更新：
a.將前一時刻粒子的sigma點集中的點帶入粒子狀態方程，得到當前時刻的粒子sigma點集公式如下：
b.對當前時刻的粒子sigma點集加權求和，從而得到當前時刻粒子的期望公式如下：
c.做當前時刻粒子的期望與當前時刻粒子的sigma點集的方差，並對所得方差進行加權求和，得到當前時刻粒子的方差的時間更新值公式如下：
d.將當前時刻粒子的sigma點集中的點以及前一時刻粒子的sigma點集中的點帶入粒子狀態方程，從而得到當前時刻粒子的觀測值公式如下：
e.對當前時刻的粒子的觀測值加權求和，得到當前時刻粒子的觀測值的更新值公式如下：
③量測更新(融入最新觀測值)：
i.基於當前時刻的粒子sigma點集當前時刻粒子的期望當前時刻的粒子的觀測值當前時刻粒子的觀測值的更新值計算得到當前時刻的粒子的方差的量測更新值公式如下：
ii.計算當前時刻粒子的方差的修正係數Kk，公式如下：
iii.基於當前時刻粒子的方差的時間更新值當前時刻的粒子的方差的量測更新值當前時刻粒子的方差的修正係數Kk，計算得到當前時刻粒子的方差的修正值公式如下：
④用當前時刻粒子的期望當前時刻粒子的方差的修正值作為高斯分布的期望與方差，最終生成建議分布
(二)採樣：從建議分布獲得採樣
重要性權值更新為：
歸一化權值更新為：
(3)重採樣：
MCMC步驟，抽取閾值u～U[0,1]，從建議分布中抽樣得到候選狀態計算接受概率若u≤A，則接受即：否則丟棄按照公式(16)進行重採樣並保留重採樣的粒子即：
(4)輸出結果：
經過上面算法的迭代，則時刻k的系統狀態估計與方差估計分別為：
若k≤T0(T0為已知測量值的個數)，則令k＝k+1並返回步驟(2)進行重要性權值更新，繼續執行後續程序。否則，結束程序，輸出結果。
(二)設備壽命預測
利用上述改進粒子濾波算法進行設備剩餘壽命預測，是一種基於模型和數據的壽命預測方法，其預測過程如下：
(1)建立預測模型
在應用改進粒子濾波算法之前，首先要建立預測對象的退化模型。一般可以根據已知的系統物理知識來建立模型，也可根據經驗來建立模型。建立模型最重要的一點是模型必須是健壯的，能夠準確快速響應系統的變化。在通用預測模型中，預測對象的系統狀態更新方程與量測方程分別為：
xk＝fk(xk-1,vk-1) (19)
zk＝hk(xk,nk) (20)
其中，xk表示當前時刻的系統狀態值，zk表示當前時刻含有加性噪聲的量測值，vk-1為前一時刻的過程噪聲，nk為當前時刻的量測噪聲，k∈N，N粒子個數，且過程噪聲與量測噪聲均為平穩噪聲序列。
(2)退化性能評估與敏感特徵選取
退化特徵的「質量」會直接影響預測結果的準確程度。所謂「質量」指的是用於預測的退化特徵反映目標系統退化趨勢的能力。通常應根據設備運行狀態信息，進行性能衰退評估，選取能較好反映設備性能退化的敏感特徵用於剩餘壽命預測。
(3)原始敏感特徵預處理
在實際預測中，通常只關心退化特徵反映系統退化的趨勢，而不關心退化特徵所包含的其他細節信息。由於噪聲等外部調節的影響，對於原始退化特徵，其退化範圍往往難以確定。另外，設備的退化過程是單調不可逆的。所以，對初步提取的退化特徵進行平滑和單調化的預處理，能夠更好的適用於粒子濾波預測。
(4)歷史數據的先驗信息獲取
在上述預測模型的基礎上，從退化特徵的歷史數據中可以獲取兩方面的先驗信息。一是退化特徵隨變量(時間或周期等)的變化情況，用於完善建立的通用預測模型。二是預測模型中各參數的初始值，與改進粒子濾波預測樣本的初始化有關。先驗信息的獲取流程如圖2所示:
(5)監測對象的預測
在獲取所需的先驗信息後，就能用改進的粒子濾波算法對已有監測數據的目標對象進行預測。首先，根據退化模型分布及模型中各參數的初始信息，對算法進行初始化。然後，依次進行迭代估計和外推預測兩個過程，當退化狀態預測值小於或等於失效閾值時，預測結束。最後，對所有粒子樣本的狀態進行統計，得到最終的預測結果。其預測流程如圖3所示。
實施例
以下結合具體實例——鋰離子電池系統狀態估計與壽命預測，來闡述本發明的具體實施方式。由於鋰離子電池工作原理複雜，其內部狀態參數難以測量，一般通過採集電池充放電過程中的電流、電壓等數據，建立容量與循環次數間的數學模型來預測其剩餘壽命。其具體過程如下：
(1)建立通用預測模型
建立目標預測模型，預測對象的系統狀態更新方程與量測方程分別為式(19)，(20)所示。
(2)敏感特徵選取
同大多數二次電池一樣，鋰離子電池在循環使用過程中，其容量會隨著循環次數的增加而越來越小，一旦下降到預設的臨界值，電池將視為失效，不能再作為可靠的能量供應。因此，可選取鋰離子電池的容量值作為壽命預測的敏感特徵。
(3)容量數據預處理
在不影響鋰離子電池容量退化趨勢的前提下，對其容量退化數據進行平滑單調化預處理，處理後的容量退化數據如圖4-7中圓點所示。
(4)歷史數據的先驗信息獲取
首先，對於鋰離子電池而言，其容量的經驗退化模型為：
Qk＝a·exp(b·k)+c·exp(d·k) (21)
式中：a,b,c,d為模型參數；k是模型變量，代表充放電循環次數。
從圖4-7曲線擬合結果可以看出，退化模型的曲線能很好地擬合容量真實退化數據。因此，可以利用其容量退化的經驗模型作為預測的量測方程，即：
Qk＝akexp(bk·k)+ckexp(dk·k)+v,v～N(0,σv) (22)
其中，v為高斯白噪聲。
鋰離子電池系統的狀態轉移模型和方程分別為：
xk＝[ak bk ck dk]T (23)
式中：wa,wb,wc,wd為高斯白噪聲，其標準差可根據模型初值確定。
這裡將電池A1，A2，A3的數據作為已知歷史數據用於獲取信息的先驗信息。A4作為檢測電池用於壽命預測。首先藉助Matlab曲線擬合工具箱對A1，A2，A3數據進行擬合，得到模型參數a,b,c,d(見表1)；為了減少特定電池參數對電池A4的影響，這裡選取a,b,c,d的平均值作為電池A4退化模型的初始值。
表1 模型參數值
根據A4電池模型參數a,b,c,d初值，以及噪聲相對初值
可以確定模型噪聲為：wa～N(0,1e-10)，wb～N(0,1e-6)，wc～N(0,1e-5)，wd～N(0,1e-7)，v～N(0,1e-3)。
(5)監測對象的預測
取採樣粒子數為N＝1000，失效閾值為電池額定容量的80％，按照改進粒子濾波算法對已知部分容量數據的監測目標(A4)進行預測。圖8為已知電池A4前18個測量數據的壽命預測結果，圖9為已知A4前30個容量數據時的壽命預測結果。從圖8可以看出預測的壽命為48個循環，預測誤差為1個周期。PDF區域表示剩餘壽命預測的概率密度函數(PDF,Probability distribution function)，其最大寬度越窄表面預測精度越高，預測結果更具參考價值。圖9的最終預測結果為49個循環，誤差為0，說明隨著監測對象的已知測量數據越多，剩餘壽命分布越集中，預測結果也越準確。詳細結果見表3。
表3 預測結果詳細值
為了說明實施例相較於基於傳統標準粒子濾波壽命預測方法的優勢，這裡在同一張圖上給出了二者在相同條件下的預測結果，如圖10和圖11所示，表3給出了詳細預測結果。從對比圖中可以看出，在相同條件下，實施例的預測精度明顯要高於標準粒子濾波的預測精度。並且，實施例在監測對象已知數據有限的情況下仍能提供較準確的預測結果，為故障的預防提供可靠的依據。
需要指出的是，上面所述只是說明本發明的一些原理，由於對相同技術領域的普通技術人員來說是很容易在此基礎上進行若干修改和改動的。因此，本說明書並非是要將本發明局限在所示和所述的具體結構和適用範圍內，故凡是所有可能被利用的相應修改以及等同物，均屬於本發明所申請的專利範圍。