李文峰，許愛強(qiáng)，王　豐，張懷遠(yuǎn)

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

?

基于改進(jìn)確定性采樣濾波的疲勞裂紋擴(kuò)展RUL預(yù)測(cè)

李文峰1,2，許愛強(qiáng)1，王豐3，張懷遠(yuǎn)4

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

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

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

0引言

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

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

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

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

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

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

(1)

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

(2)

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

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

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

(3)

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

(4)

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

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

(5)

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

(6)

其中:σj和ωj分別表示積分節(jié)點(diǎn)和權(quán)值，與濾波理論中的采樣點(diǎn)和權(quán)值相對(duì)應(yīng)。

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

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

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

對(duì)積分(5)，其積分區(qū)域?yàn)橥耆珜?duì)稱，如果隨機(jī)變量之間為獨(dú)立同分布，那么權(quán)值函數(shù)也是全對(duì)稱的，以生成元[9]的方式表示式(6)，其完全對(duì)稱積分公式：

(7)

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

(8)

完全對(duì)稱積分公式[2]的特點(diǎn)為：

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

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

(9)

即

(10)

針對(duì)相應(yīng)的精度，生成元的選取與單項(xiàng)式相關(guān)。生成元中非零元素的個(gè)數(shù)與單項(xiàng)式中階次非零的變量個(gè)數(shù)相等。

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

(11)

當(dāng)隨機(jī)變量服從標(biāo)準(zhǔn)正態(tài)分布時(shí)，I0=I2=I2,2=1，I4=3。據(jù)式(11)，有

(12)

3實(shí)例仿真分析

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

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

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

(13)

(14)

(15)

(16)

(17)

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

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

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

3.1數(shù)據(jù)準(zhǔn)備

(18)

(19)

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

量測(cè)方程為：

(20)

其中:量測(cè)噪聲ν(t)為零均值高斯噪聲，標(biāo)準(zhǔn)差為0.074。

3.2狀態(tài)估計(jì)

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

圖1基于UKF、CKF以及IDSF的疲勞裂紋擴(kuò)展?fàn)顟B(tài)估計(jì)

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

3.3評(píng)價(jià)指標(biāo)

采用平均相對(duì)誤差作為估計(jì)預(yù)測(cè)評(píng)價(jià)指標(biāo)，公式為：

(21)

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

表1　3種方法的仿真時(shí)間和平均相對(duì)誤差

3.4RUL預(yù)測(cè)

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

圖2　疲勞裂紋擴(kuò)展的IDSF剩余使用壽命預(yù)測(cè)

4結(jié)論

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

參考文獻(xiàn)：

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

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

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

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

[5] 張勇剛，黃玉龍，武哲民,等.一種高階無跡卡爾曼濾方法[J]. 自動(dòng)化學(xué)報(bào)，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] 袁慎芳, 張華, 邱雷,等. 基于粒子濾波算法的疲勞裂紋擴(kuò)展預(yù)測(cè)方法[J]. 航空學(xué)報(bào), 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

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

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

中圖分類號(hào)：TP273

文獻(xiàn)標(biāo)識(shí)碼：A

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

基金項(xiàng)目：總裝武器裝備預(yù)研基金(9140A27020214JB14436)。

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