李文峰，許愛強，王　豐，張懷遠

(1.海軍航空工程學院 科研部, 山東 煙臺　264001； 2.92635部隊, 山東 青島　266000；3.91206 部隊, 山東 青島　266000； 4.91359部隊, 北京　102443)

?

基于改進確定性采樣濾波的疲勞裂紋擴展RUL預測

李文峰1,2，許愛強1，王豐3，張懷遠4

(1.海軍航空工程學院 科研部, 山東 煙臺264001； 2.92635部隊, 山東 青島266000；3.91206 部隊, 山東 青島266000； 4.91359部隊, 北京102443)

摘要：針對確定性采樣濾波在進行狀態(tài)估計預測時隨維數(shù)增加時出現(xiàn)計算量增加且精度不高的問題，提出一種確定性采樣濾波的算法并將其應用到疲勞裂紋擴展RUL預測當中去；首先，闡述了確定性采樣濾波的基本原理；其次，從多維數(shù)值積分的角度分析確定性采樣濾波所需計算的數(shù)學期望，根據(jù)完全對稱積分公式計算積分節(jié)點值、節(jié)點個數(shù)和權重；最后，將改進后的確定性采樣濾波器應用到構件疲勞裂紋損傷擴展中去，并與無跡卡爾曼濾波算法、容積卡爾曼濾波算法進行比較，提升了裂紋擴展RUL預測的精度，實例仿真分析驗證了該方法的可行性和有效性。

關鍵詞：確定采樣型濾波器；疲勞裂紋擴展；多維數(shù)值積分；剩余使用壽命預測

0引言

由于擴展卡爾曼濾波(EKF)方法存在需要求解雅可比矩陣、數(shù)值穩(wěn)定性較差和濾波精度不高等明顯缺點，已經不能滿足非線性系統(tǒng)的濾波要求。近年來，隨著對非線性系統(tǒng)統(tǒng)計濾波領域研究的日益深入，相關專家學者提出了無跡卡爾曼濾波(UKF)、高斯厄米特濾波(GHF)和容積卡爾曼濾波(CKF)等確定采樣性濾波方法。這些濾波器都是以確定的數(shù)學解析式完成采樣，統(tǒng)一稱為確定性采樣濾波器或確定采樣型濾波器[1-4](deterministic sampling filter， DSF)，不同的確定性采樣濾波方法區(qū)別在對濾波方法中均值和方差的計算，其中UKF采用的是UT變換方法[5]，GHF采用的是Gauss-Hermit積分方法計算采樣點[6]，CKF采用的是球面容積法則[7]。確定性采樣濾波方法的提出，奠定了其在信號處理、金融、導航和估計預測等領域應用的理論基礎。確定采樣性濾波方法在保證濾波精度的同時，需要大量的計算，這在實際應用時比較困難。

本文針對確定采樣型濾波器中存在的問題，對確定采樣型濾波器改進計算的基礎上，在至少不增加計算量的同時，提升確定采樣濾波器的估計精度。同時，結構疲勞損傷擴展受環(huán)境、材料等多種參數(shù)的不確定性影響，本身是一個隨機過程，利用濾波對疲勞裂紋擴展進行預測，非線性濾波方法進行疲勞裂紋擴展預測時可以消除這種不確定性。因此，針對構件疲勞裂紋擴展預測的特點，本文將改進后的確定采樣型濾波器應用于構件疲勞裂紋擴展剩余使用壽命(remaining useful life，RUL)預測中。

1確定性采樣濾波方法基本原理

確定采樣性濾波方法是基于線性最小方差估計框架，利用系統(tǒng)狀態(tài)的一階矩和二階矩(均值和方差)，依據(jù)不同的方法，確定狀態(tài)先驗概率的均值和方差構造的采樣點及相對應的權值，利用采樣點對狀態(tài)非線性函數(shù)的均值和方差進行近似估計。不同的確定采樣性濾波方法區(qū)別在對濾波方法中均值和方差的計算。

通常，確定采樣性濾波多假設狀態(tài)后驗分布為高斯型。

對于非線性系統(tǒng)

(1)

其中：xk和zk分別為狀態(tài)變量和量測向量，fk(xk)和hk(xk)為已知函數(shù)，wk和vk分別為隨機系統(tǒng)噪聲和隨機量測噪聲，且相互獨立。

(2)

對于式(2)，無法直接求解，只能通過一些數(shù)值算法近似計算，確定采樣性濾波方法的估計精度取決于公式(2)的計算精度，通過不同的數(shù)值近似方法計算公式(2)將衍生出不同的濾波方法。下面通過利用數(shù)值積分容積方法，推導出相應的確定采樣性濾波方法。

2基于完全對稱積分公式的多維數(shù)值積分

對確定的非線性函數(shù)y=f(x)，假設p(x)為x的概率密度函數(shù)，則式(2)需計算積分

(3)

若用函數(shù)g(x) 表示被積函數(shù)的其他部分，則計算其數(shù)學期望：

(4)

對公式(4)，其積分的計算精度與選取的積分點有關。計算出該積分，即可計算出隨機變量的均值和方差，從而完成對結果分布的高斯近似。

通常解決實際問題時，遇到的是多維隨機變量x1,…,xn，其函數(shù)g(x1,…,xn)的數(shù)學期望為：

(5)

利用數(shù)值積分方法獲得式(5)的近似解。數(shù)值積分的計算公式為：

(6)

其中:σj和ωj分別表示積分節(jié)點和權值，與濾波理論中的采樣點和權值相對應。

確定節(jié)點個數(shù)、節(jié)點和權值，才能用Q[g]近似I[g]。節(jié)點個數(shù)的多少對計算量的大小會產生比較大的影響。節(jié)點個數(shù)越少，計算量就越小。在節(jié)點個數(shù)相同的情況下，選取不同的節(jié)點和權值，積分精度也會有較大的不同。

為了計算公式(5)，相關專家學者作了大量的研究。利用高斯厄米特數(shù)值積分方法，選取高斯點和相應權值，提高系統(tǒng)狀態(tài)的均值和方差估計的精度，誕生了高斯厄米特濾波方法[6]。但是，在高維情況下，該方法的計算負荷將呈指數(shù)增長。McNamee和Stenger定義了完全對稱積分公式[8]，通過求解非線性方程組的方法來得到采樣點。該方法的計算負荷遠小于高斯厄米特數(shù)值積分方法的計算負荷。

完全對稱積分區(qū)域包括無窮區(qū)域Rn和有界區(qū)域如超球體和超球面。

對積分(5)，其積分區(qū)域為完全對稱，如果隨機變量之間為獨立同分布，那么權值函數(shù)也是全對稱的，以生成元[9]的方式表示式(6)，其完全對稱積分公式：

(7)

若代數(shù)精度為d=2s-1，則其完全對稱積分公式，所需要的最少節(jié)點個數(shù)為：

(8)

完全對稱積分公式[2]的特點為：

1)完全對稱積分公式能夠準確積分任何包含奇次冪的單項式，如果g(x)在某個方向上含有奇次冪項，那么I=0。

2)如果g(x)只包含偶次冪項，那么積分的結果與生成元中元素的排列無關。因此，建立完全對稱積分公式的精度為d=2k+1，只需要使積分公式對如下單項式能準確成立：

(9)

即

(10)

針對相應的精度，生成元的選取與單項式相關。生成元中非零元素的個數(shù)與單項式中階次非零的變量個數(shù)相等。

當d≥5時，會出現(xiàn)多個變量的單項式，此時選取的生成元組合會出現(xiàn)不同的情況。對d=5時的生成元為[0]，[u]1=[(u1,0,…,0)]，[u]2=[(u2,u2,0,…,0)]，根據(jù)式(8)積分節(jié)點個數(shù)為n2+n+1，根據(jù)式(10)有：

(11)

當隨機變量服從標準正態(tài)分布時，I0=I2=I2,2=1，I4=3。據(jù)式(11)，有

(12)

3實例仿真分析

采用完全對稱積分公式對多維數(shù)值積分的積分節(jié)點的采樣策略進行改進后，給出基于改進確定性采樣濾波(IDSF)的剩余使用壽命預測步驟：

步驟1：在建立狀態(tài)空間模型時，最佳途徑是能夠通過設備或部件的故障機理建立物理損傷數(shù)學模型，對于復雜的故障機理不能建立數(shù)學模型時，可根據(jù)已知的狀態(tài)退化數(shù)據(jù)，根據(jù)ARMA模型、Gamma過程、HMM模型等數(shù)學統(tǒng)計模型描述部件的退化過程。

經狀態(tài)方程中非線性函數(shù)fk(·)將ξi,k傳遞為γi,k+1|k，即γi,k+1|k=fk(ξi,k)。

(13)

(14)

(15)

(16)

(17)

步驟5：給定新的量測信息后，執(zhí)行步驟3、4，將計算結果代入貝葉斯估計框架下所需的5個數(shù)學期望進行濾波更新，得到下一時刻的狀態(tài)估計。

步驟6：在狀態(tài)估計的基礎上，執(zhí)行后期的計算處理，實現(xiàn)剩余使用壽命預測，并給出預測結果的概率密度分布。

給定基于改進確定性采樣濾波的預測方法、步驟后，在實例仿真分析中，主要從數(shù)據(jù)準備、狀態(tài)估計、評價指標以及RUL預測4個方面實現(xiàn)改進求積公式的確定性采樣濾波RUL預測，選用應用較廣泛的UKF及近年來新提出并逐漸成為研究熱點的CKF與本文所提的IDSF進行仿真比較。

3.1數(shù)據(jù)準備

(18)

(19)

其中:ω1(t)為高斯噪聲，均值為0.045，標準差為0.116；ω2(t)為零均值高斯噪聲，標準差為0.01，初始協(xié)方差矩陣P0|0=[0.10;00.1]。其他參數(shù)參照文獻[10-11]。

量測方程為：

(20)

其中:量測噪聲ν(t)為零均值高斯噪聲，標準差為0.074。

3.2狀態(tài)估計

結構疲勞損傷擴展受環(huán)境、材料等多種參數(shù)的不確定性影響，本身是一個隨機過程，利用濾波對疲勞裂紋擴展進行預測，后期得到的概率密度函數(shù)表達了這種不確定性。仿真時間100個載荷循環(huán)周期，在相同初始條件下采用CKF、UKF、IDSF分別進行了仿真，仿真結果如圖1所示。

圖1基于UKF、CKF以及IDSF的疲勞裂紋擴展狀態(tài)估計

從圖1狀態(tài)估計曲線可以看出，UKF和CKF在48個循環(huán)周期后對真實狀態(tài)的跟蹤性能變差，但UKF的性能在第80個循環(huán)周期后稍優(yōu)于CKF，這是因為當系統(tǒng)狀態(tài)維數(shù)不大于3時，UKF相比CKF，更適合解決非線性狀態(tài)估計問題；當系統(tǒng)狀態(tài)維數(shù)達到20以上時，UKF因參數(shù)選擇經常出現(xiàn)濾波發(fā)散的問題，相比UKF，CKF在采樣和濾波過程中，不需要進行參數(shù)的選擇計算且不受狀態(tài)維數(shù)高低的影響，權值始終保持為正值,算法應用范圍較廣。但是從仿真分析來看，UKF和CKF都不如本文提出的IDSF。

3.3評價指標

采用平均相對誤差作為估計預測評價指標，公式為：

(21)

表1對比了UKF、CKF以及本文提出的IDSF方在進行疲勞裂紋擴展RUL預測時的仿真時間和平均相對誤差，從表1中可以看出，從起始裂紋長度為6.93 mm開始，本文提出的IDSF方法的仿真時間大于CKF的仿真時間，而小于UKF的仿真時間，這說明UKF的計算復雜度要高于IDSF和CKF，在平均相對誤差方面，UKF的平均相對誤差稍小于CKF的平均相對誤差，而本文提出的IDSF方法平均相對誤差要小于UKF和CKF二者，這說明將3種方法用于疲勞裂紋擴展RUL預測中當屬本文提出的IDSF方法精度最高。

表1　3種方法的仿真時間和平均相對誤差

3.4RUL預測

上述分析可以看出，IDSF在進行疲勞裂紋擴展狀態(tài)估計時，估計精度優(yōu)于UKF和CKF，因此選擇本文方法進行疲勞裂紋擴展預測，在疲勞裂紋擴展狀態(tài)估計的基礎上，對數(shù)據(jù)做后期進一步處理，在此過程中，IDSF濾波消除了剩余使用壽命中的不確定性，從圖2中可以看出，在給定疲勞裂紋擴展長度失效閾值為40mm后，確定了RUL預測的預測區(qū)間，同時給出了預測結果的概率密度分布，從圖2中可以得到，以疲勞裂紋長度為6.93 mm為當前起始點，則疲勞裂紋擴展的剩余使用壽命為83個循環(huán)周期通過計算，疲勞裂紋擴展的預測區(qū)間為[70,100]，置信度為90%。

圖2　疲勞裂紋擴展的IDSF剩余使用壽命預測

4結論

本文針對確定性采樣濾波所需計算的數(shù)學期望，采用多維數(shù)值積分求積公式近似求取數(shù)學期望，采用完全對稱積分公式求取積分節(jié)點的位置、個數(shù)和權重，改進采樣節(jié)點的選取策略，提高剩余使用壽命預測的精度，通過疲勞裂紋擴展的算法仿真對比，本文提出的IDSF方法精度高于UKF和CKF，驗證了該方法的有效性和可行性。

參考文獻：

[1] 劉偉, 楊峰, 張洪才,等. 狀態(tài)估計中確定性采樣濾波器的比較分析[J]. 系統(tǒng)仿真學報， 2007, 19(18): 4265-4269.

[2] 叢源林，李文峰，王豐.基于多維數(shù)值積分的高階確定采樣型濾波方法[J]. 海軍航空工程學院學報， 2014, 29(6): 547-551.

[3] 王小旭, 潘泉, 黃鶴,等. 非線性系統(tǒng)確定采樣型濾波算法綜述[J]. 控制與決策. 2012, 27(6): 801-812.

[4] 劉濤, 解永春. 一種自適應確定性采樣濾波方法[J]. 信息與控制, 2010, 39(6): 673-680.

[5] 張勇剛，黃玉龍，武哲民,等.一種高階無跡卡爾曼濾方法[J]. 自動化學報，2014, 40(5): 838-847.

[6] Ito K, Xiong K. Gaussian filters for nonlinear filtering problems[J]. IEEE Transactions on Automatic Control, 2000, 45(5): 910-927.

[7] Arasaratnam I, Haykin S. Cubature Kalman Filters[J]. IEEE Transactions on Automatic Control, 2009, 54(6): 1254-1269.

[8] McNamee J, Stenger F. Construction of fully symmetric numerical integration formulas[J]. Numer. Math., 1967(10)：327-344.

[9] Stroud A H. Secrest D. Approximate integration formulas for certain spherically symmetric regions[J]. Math. Comput., 1963，17：105-135.

[10] Baraldi P, Compare M, Sauco S, et al. Ensemble neural network-based particle filtering for prognostics[J]. Mechanical System and Signal Processing, 2013, 41: 288-300.

[11] 袁慎芳, 張華, 邱雷,等. 基于粒子濾波算法的疲勞裂紋擴展預測方法[J]. 航空學報, 2013, 34(12): 2740-2747.

Research on RUL Prediction of Improved Deterministic Sampling Filter for Fatigue Crack Propagation

Li Wenfeng1,2, Xu Aiqiang1, Wang Feng3, Zhang Huaiyuan4

(1.Department of Scientific Research, Naval Aeronautical and Astronautical University, Yantai264001,China; 2.92635 Troop, Qingdao266000,China; 3.91206 Troop, Qingdao266000,China; 4.91359 Troop, Beijing102443,China)

Abstract:In view of the problem of increasing the accuracy of the computation of the deterministic sampling filter in the state estimation, a new algorithm is proposed and applied to the prediction of fatigue crack growth. First of all, the basic principle of determining the sampling filter is described. Secondly, from the point of view of multidimensional numerical integration, the mathematical expectation of deterministic sampling filter is analyzed, and the integral node value, the number of nodes and the weights are calculated according to the fully symmetrical integral formula. Finally, the improved model is applied to the fatigue crack propagation of the members, and the accuracy of the crack propagation prediction is improved by comparing with unscented Kalman filtering algorithm and cubature Kalman filtering algorithm. The feasibility and effectiveness of the proposed method are verified by simulation analysis.

Keywords:deterministic sampling filter; fatigue crack propagation; multidimensional numerical integration; remaining useful life prediction

文章編號：1671-4598(2016)02-0205-04

DOI：10.16526/j.cnki.11-4762/tp.2016.02.057

中圖分類號：TP273

文獻標識碼：A

作者簡介：李文峰(1983-)，男，山東榮成人，博士生，主要從事航空裝備故障預測與綜合保障研究。許愛強(1964-)，男，山東即墨人，教授，博士生導師，主要從事自動測試與裝備綜合保障研究。

基金項目：總裝武器裝備預研基金(9140A27020214JB14436)。

收稿日期：2015-08-16；修回日期：2015-09-17。