基于輔助粒子濾波的疲勞裂紋擴展預測研究

楊偉博， 袁慎芳， 邱 雷， 陳 健(南京航空航天大學 機械結構力學及控制國家重點實驗室, 南京 210016)

疲勞裂紋是當前飛行器的主要損傷形式之一，實現疲勞裂紋的擴展預測，對飛行器延壽，提升飛行安全、降低維護成本具有重要意義。然而結構疲勞裂紋擴展是一個隨機過程，具有非線性非高斯的復雜動態特性，它包含各種不確定因素，如材料參數的變化、環境參數的變化、測試數據的分散性、建模的局限性等，如何克服此類不確定性的因素的影響，實現疲勞裂紋這類非線性非高斯動態系統的預測，一直是學者們研究的重點。

目前，可用于實現裂紋擴展的壽命預測方法的主要有失效物理(Physic of Failure, PoF)方法、數據驅動(Data-Driven)方法和融合(Hybrid)算法[1-3]。其中，融合算法意在融合前兩種算法的優點，使壽命預測算法具備更低的計算消耗和復雜度、需要更少的實驗與訓練數據，因此具備更好的發展前景。而作為融合算法的一種，基于粒子濾波(Particle Filter, PF)[4]的疲勞裂紋壽命預測方法結合了Paris規則的裂紋擴展規律，有限元建模和粒子濾波算法，與其他兩種方法相比，其建模相對較容易，并不需要大量的實驗與訓練數據，且與Lamb波主動監測方法相結合，可實現在線預測，因此更適合于工程應用。

近年來，Orchard等[5-6]利用粒子濾波算法，建立了一種裂紋擴展在線預測系統，對UH-60行星齒輪進行了壽命預測，其預測結果相比傳統的卡爾曼濾波等方法具有更高的預測精度，且能實時修正監測過程中的預測誤差，但觀測方程采用應變作為監測手段，觀測方程的精度還有進一步的提高空間；Cadini等[7-8]結合文獻中的仿真數據，利用粒子濾波算法實現了疲勞裂紋擴展的預測；袁慎芳等[9]以Paris規則作為狀態方程，以Lamb波主動健康監測的在線監測技術建立觀測方程，利用標準粒子濾波算法實現了中心裂紋鋁板上的裂紋擴展在線預測，有效地消除了裂紋擴展參數及觀測手段等多種不確定性的影響，獲得了較高的預測精度。

從上述粒子濾波的裂紋損傷預測研究可以看出，目前相關領域的研究都是基于最簡單的標準粒子濾波算法利用仿真數據或在簡單結構中開展的。然而由于采用了次優估計[10]，其預測結果并未考慮當前觀測更新的影響，因而依賴于建模的精確性。當模型誤差較大時，將有大量粒子處于低相似區中而在重采樣時被浪費，從而加劇了粒子多樣性匱乏現象，即原有粒子群中很多粒子由于預測偏差大而權值太小沒有“后代”，而少數幾個權值較高的粒子則有很多相同的“后代”，重采樣以后的粒子群由大量重復的粒子構成，使粒子群失去了多樣性，粒子多樣性匱乏現象的加劇，將導致預測精度下降。

解決上述問題的方法之一是增加粒子個數，然而這將影響計算效率。另一種方法是結合當前觀測更新選擇重要性函數。本文提出基于輔助粒子濾波(Auxiliary Particle Filter, APF)[11-12]的裂紋擴展壽命預測方法，通過Pairs規則建立狀態方程，以基于壓電主動Lamb波在線監測方法作為觀測方程，構建描述孔邊裂紋擴展的狀態空間模型。最后，基于孔邊裂紋的APF與PF壽命預測結果表明，APF算法充分結合裂紋在線觀測值，使重采樣后的粒子向似然函數的高值區移動，提高了先驗估計與似然函數之間重疊性，能有效緩解多樣性匱乏現象，克服PF算法過于依賴模型的缺點，具備更高的預測精度。

1 粒子濾波理論

結構疲勞裂紋擴展是一個動態變化的過程，它存在諸多的不確定性，其數學描述可用狀態空間模型來表征。狀態空間模型包含狀態方程和觀測方程

(1)

式中:xt表示t時刻裂紋狀態值；f(·)為系統的狀態方程，表征t時刻裂紋向t+1時刻擴展的規律；zt+1表示t+1時刻的觀測值，此處指通過健康監測手段所獲取的損傷特征參數；g(·)為系統的觀測方程；ωt為t時刻的過程噪聲，用于表征結構上的各種不確定影響因素；υt為t時刻的觀測噪聲，用于表征健康監測手段的不確定性監測誤差。

粒子濾波具備平滑、濾波和預測的功能，在裂紋擴展預測的實現，是根據0到t時刻的裂紋狀態值和1到t時刻的裂紋觀測值，結合狀態空間模型和t+1時刻的最新觀測值，去預測t+1時刻及其以后的裂紋狀態值，最后再根據裂紋參數計算結構是否可以使用，確定結構的壽命。上述過程實質是通過Bayes估計來實現裂紋擴展過程的無偏估計

(2)

式中狀態值x的上標“^”表征后驗估計，右上標“-”表征先驗估計，p(xt+1|z1:t+1)表征觀測值為z1,z2,…zt+1時，出現狀態值xt+1的概率密度函數，由于它是通過t+1時刻的觀測值zt+1修正的，因此常稱為后驗概率密度函數。p(x0:t+1|z1:t+1)為p(xt+1|z1:t+1)的邊緣概率密度。

假設裂紋擴展滿足一階Markov過程，則后驗概率密度函數p(xt+1|z1:t+1)可通過Bayes估計迭代計算獲得

(3)

然而式(3)中在迭代過程中包含多重積分，一般情況下并無解析解，因此采用Monte Carlo仿真方法，從后驗概率p(x0:t+1|z1:t+1)中抽取N個獨立同分布的樣本，當N～+∞時，樣本的均值可逼近于數學期望，即

(4)

但在實際問題中，往往很難從p(x0:t+1|z1:t+1)中進行抽樣，因此引入了一個容易采樣的重要性函數q(x0:t+1|z1:t+1)，該參考分布已知且與p(x0:t+1|z1:t+1)同分布，則最終可將裂紋擴展的無偏估計式(2)化簡為

(5)

式中：

(6)

以上過程即為序貫重要性采樣(Sequential Importance Sampling, SIS)，最終通過裂紋先驗估計與其權值的加權，可實現裂紋擴展過程的無偏估計。然而SIS算法存在嚴重的退化現象，因其最優解q(xt+1|xt,z1:t+1)=p(xt+1|xt,z1:t+1)在實際問題中常無法得到且對應積分不易求解。為此，Gordon提出使用次優解q(xt+1|xt,z1:t+1)=p(xt+1|xt)與重采樣相結合以解決退化問題，即序貫重要性采樣重采樣 (Sequential Importance Resampling, SIR)，又稱粒子濾波算法。

2 輔助粒子濾波

SIR算法解決了退化問題且抽樣簡單，易于實現，然而由于重要性函數使用了次優解p(xt+1|xt)，對比最優解p(xt+1|xt,z1:t+1)，次優解忽略了觀測值z1:t+1的影響，更依賴于狀態方程xt+1=f(xt,ωt)的建模精度，因此SIR(PF)算法存在過于依賴模型的特點。當重要性函數與真實分布非常接近時，SIR算法能達到很好的預測效果，但如果二者偏差較大，則預測效果將明顯下降，甚至導致預測失敗。

為克服上述缺陷，Pitty等提出了輔助粒子濾波算法，通過引入輔助變量l，將其與最新觀測值的聯合概率密度函數作為重要性函數，預先對原有粒子群進行抽樣處理，使抽樣后的粒子群向真實分布的高似然區移動，以此來克服算法過于依賴模型的缺點，提高預測精度。

假設t時刻p(xt|z1:t)已知，則t+1時刻的“經驗”后驗概率分布為

(7)

定義結合了輔助變量的聯合概率密度函數

(8)

(9)

(10)

(11)

對式(9)進行邊緣概率積分可得

(12)

輔助粒子濾波算法執行步驟如下：

3 疲勞裂紋擴展的狀態空間模型

3.1 狀態方程的建立

金屬結構裂紋的擴展存在三個階段：萌生階段、穩定擴展階段和快速擴展階段。其中萌生階段和快速擴展階段的持續時間短，因此工程上進行結構壽命預測時，常采用穩定擴展階段，此階段又稱為Paris規則區，裂紋擴展速率的Paris公式如式(13)所示

(13)

式中：a為表征裂紋長度；Nf為疲勞載荷周期數，C與m為材料參數，ΔK為應力強度因子幅，可通過有限元離線計算獲得，此處采用奇異單元法計算獲得。

假設兩次擴展之間的時間間隔ΔNf足夠小，則根據Paris規則可建立如式(14)所示的狀態方程，其中為使過程噪聲包含更大的預測空間，狀態噪聲采用加性的零均值高斯白噪聲[13]，而非相乘的對數正態分布噪聲。

at+1=at+Ct(ΔK)mΔNf+ω(t)

(14)

式中:at+1為t+1時刻的裂紋長度；Ct為t時刻的材料及環境相關的參數，針對特定的試件，C0可通過實驗或查找數據庫獲得，此處采用實驗手段，并結合曲先強驗證所總結材料參數分布假設[14]，令lgC0～N(-14.307 5, 0.174 0)，m取值為4.038 3，用于初始化狀態方程的不確定性；ω(t)～N(0,0.5)為高斯白噪聲，用于表征狀態方程預測過程中的不確定性。

式(14)存在兩個不確定分布Ct與ω(t)，它們分別分于初始化表征Paris規則的不確定性以及預測過程中的不確定性，二者的結合充分考慮了疲勞裂紋擴展的隨機因素，提升了狀態方程的魯棒性。

3.2 觀測方程的建立

觀測方程在裂紋擴展預測中的意義是將裂紋擴展尺寸與觀測手段的某種特征關聯起來，定量地描述裂紋擴展與觀測手段之間的規律。Lamb波主動健康監測技術具有對小損傷敏感的特性，因此目前已被廣泛應用于裂紋擴展的監測中[15-17]。該方法通過監測裂紋擴展時結構中傳播的Lamb波的變化，量化分析裂紋尺寸與Lamb波響應信號間的規律，如圖1所示為Lamb波響應信號隨裂紋尺寸的變化。

圖1 響應信號隨裂紋尺寸的變化

由圖1可知，金屬疲勞裂紋擴展對Lamb波信號傳播的影響主要表現為相位及幅值的變化，對應后續所開展的孔邊裂紋驗證實驗，主要為裂紋尖端的擴展所導致的Lamb波直達波相位延遲及幅值衰減現象。其中相位延遲現象在本實驗中占主導地位，因此根據無損檢測中飛行時間(Time of Fight)的概念，定義如式(15)所示的時間延遲損傷因子，用于表征結構損傷程度與相位延遲之間的關系。

(15)

式中:T為裂紋擴展過程中直達波的飛行時間;T0為健康狀態下直達波的飛行時間。實驗準備初期，先通過幾組實驗獲得DI與裂紋觀測值，再通過多項式擬合，建立損傷因子與裂紋擴展之間的觀測方程

DIt+1=g(at+1)+υ(t)

(16)

聯立狀態方程與觀測方程，描述疲勞裂紋擴展的狀態空間模型如式(17)所示

(17)

4 疲勞裂紋擴展預測的輔助粒子濾波實現

基于式(17)所建立的狀態空間模型，為解決PF過于依賴模型的缺點，緩解多樣性匱乏現象，將APF算法應用于裂紋擴展預測中。裂紋擴展的輔助粒子濾波實現流程如下，流程圖如圖2所示。

步驟1 初始化。由Lamb波在線監測獲得裂紋的萌生尺寸初始化裂紋狀態值a0及其對應循環周期數t0；由n組重復性疲勞實驗通過割線法初始化材料參數C0與m的分布；初始化過程噪聲ω(t)和觀測噪聲υ(t)分布；將a0和C0代入狀態方程中，初始化裂紋預測粒子群{a0}N，權值{w0}N=1/N；

(18)

(19)

注：t與t+1的時間間隔為ΔNf。

圖2 APF預測流程圖

5 實驗驗證

選用孔邊裂紋試驗件，結合狀態空間模型，對輔助粒子濾波改進算法的正確性進行驗證。試件如圖3所示，尺寸為240 mm×100 mm×2 mm，材料為LY12鋁合金，共計6件，編號為T1-T6，傳感器對稱分布于裂紋兩側，間距80 mm，孔邊預制3 mm裂紋。實驗加載設備采用MTS810，試件T1先進行靜載破壞實驗，確定所能承受的最大拉力為45 kN，在此基礎上，試件T2-T5采用正弦載荷等幅譜進行疲勞實驗，載荷峰值為15 kN，應力比為0.1，載荷頻率為10 Hz。疲勞實驗過程中，裂紋長度每增加2 mm，即保持載荷至5 kN，采集Lamb波響應信號，并記錄循環周期數，監測系統采用南京航空航天大學結構健康監測與預測研究中心自主研發的PHM掃查系統[18-19]，激勵頻率為290 kHz，采樣率為10 MHz。

DI=-6.234×10-5a3+4.854×10-3a2-

4.069×10-2a+9.231×10-2+υ(t)

(20)

作為驗證的T6試件，其裂紋擴展過程如表1所示，圖4所示為T2-T5的觀測方程擬合曲線與T6實驗損傷因子的對比，可見T6的實驗值均落于擬合曲線的95%的置信區間內，其最大誤差為0.144 3，驗證了觀測手段與觀測擬合方程的正確性。

表1 疲勞裂紋擴展過程

圖4 T6實驗損傷因子與觀測方程擬合曲線的對比

選取粒子數N為500，載荷循環迭代間隔ΔNf為40，分別利用PF與APF算法，在實時獲取T6試件的損傷因子更新時，實現T6試件的實時監測與疲勞壽命預測。APF與PF的預測結果如圖5所示，以預測值與真實值的均方根誤差(Root Mean Square Error, RMSE)作為預測精度對比指標，RMSE越小，則精度越高，從圖5中可以看出，APF的RMSE為0.438 5，低于PF的0.712 9，因此結合了最新Lamb波DI觀測值的APF具備比PF更優秀的預測效果。

圖5 PF與APF的預測結果

剃除重采樣后的重復粒子(材料參數Ct)，對比APF與PF算法的粒子多樣性匱乏程度，以重采樣后具有不同數值的材料參數Ct粒子個數作為多樣性指標，則APF與PF在迭代過程中粒子多樣性變化如圖6所示。

從圖 6可以看出，相比PF算法，APF算法的粒子多樣性匱乏現象得到了緩解，預測精度得到提升。此外，由于進行了二次重采樣，運算時間由0.276 1 s增至0.299 3 s(計算采用MATLAB R2013a的64位4G內存Win7系統的筆記本電腦)，雖然增加了部分計算量，但預測精度獲得了提升，且這樣的計算消耗，對在線監測的實現影響不大。

圖6 APF與PF多樣性的對比

將PF算法的中粒子數N由500增加至2 500，則APF與PF算法在不同粒子數下的預測結果如圖7所示。從圖7中可以看出，隨著粒子數的增加PF算法的預測精度增加，RMSE由0.712 9降至0.415 2，相比于粒子數為500的APF算法，其RMSE=0.438 5，PF算法的預測精度稍顯優勢。然而PF算法由于粒子數的增加，運算時間由0.276 1 s增至0.804 6 s。且若將此部分運算消耗提供給APF算法，APF將獲得更高的預測精度，因此APF算法的壽命預測性能要優于PF算法。

圖7 APF與PF不同粒子數的預測結果

6 結論

APF算法結合了Lamb波主動健康監測技術獲得的最新裂紋觀測值，克服了傳統PF算法過于依賴模型的特點，在孔邊裂紋結構的壽命預測中，有效地緩解了多樣性匱乏現象，雖然算法計算時間存在小幅度的增加，但并不影響在線預測的實現，預測精度相比PF得到了整體提升，其RMSE由0.712 9降低至0.438 5，更適用于結合了Lamb波在線監測的裂紋擴展實時壽命預測。

[1] ZHANG H, KANG R, PECHT M. A hybrid prognostics and health management approach for condition-based maintenance[C]∥2009 IEEE International Conference on Industrial Engineering and Engineering Management. Hongkong, China: IEEE, 2009.

[2] PECHT M, GU J. Physics-of-failure-based prognostics for electronic products[J]. Transactions of the Institute of Measurement and Control, 2009, 31(3/4): 309-322.

[3] SHEPPARD J W, KAUFMAN M A, WILMER T J. IEEE standards for prognostics and health management[J]. Aerospace and Electronic Systems Magazine, 2009, 24(9): 34-41.

[4] GORDON N J, SALMOND D J, SMITH A F M. Novel approach to nonlinear/non-Gaussian Bayesian state estimation[J].IEE Proceedings F (Radar and Signal Processing). IET Digital Library, 1993, 140(2): 107-113.

[5] ORCHARD M, VACHTSEVANOS G. A particle filtering approach for on-line failure prognosis in a planetary carrier plate[J]. International Journal of Fuzzy Logic and Intelligent Systems, 2007, 7(4): 221-227.

[6] ORCHARD M E, VACHTSEVANOS G J. A particle-filtering approach for on-line fault diagnosis and failure prognosis[J]. Transactions of the Institute of Measurement and Control, 2009, 31(3-4): 221-246.

[7] CADINI F, ZIO E, AVRAM D. Monte Carlo-based filtering for fatigue crack growth estimation[J]. Probabilistic Engineering Mechanics, 2009, 24(3): 367-373.

[8] ZIO E, PALONI G. Particle filtering prognostic estimation of the remaining useful life of nonlinear components[J]. Reliability Engineering & System Safety, 2011, 96(3): 403-409.

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

YUAN Shenfang, ZHANG Hua, QIU Lei, et al. A fatigue crack growth prediction method based on particle filter[J]. Acta Aeronautica et Astronautic Sinica, 2013, 34(12): 2740-2747.

[10] DOUCET A, GODSILL S, ANDRIEU C. On sequential Monte Carlo sampling methods for Bayesian filtering[J]. Statistics and Computing, 2000, 10(3): 197-208.

[11] PITT M K, SHEPHARD N. Filtering via simulation: auxiliary particle filters[J]. Journal of the American Statistical Association, 1999, 94(446): 590-599.

[12] PITT M K, SHEPHARD N. Auxiliary variable based particle filters[M]∥Sequential Monte Carlo Methods in Practice. New York: Springer, 2001: 273-293.

[13] SUN J, ZUO H, WANG W, et al. Prognostics uncertainty reduction by fusing on-line monitoring data based on a state-space-based degradation model[J]. Mechanical Systems and Signal Processing, 2014, 45(2): 396-407.

[14] 曲先強, 馬永亮, 崔洪斌, 等．Paris 公式中材料參數的統計特性分析[C]//2008 全國MTS斷裂測試研討會論文集, 2008: 26-31.

[15] 曹俊, 袁慎芳, 蔡建, 等．疲勞裂紋擴展的實時健康監測[J]．壓電與聲光, 2008, 30(6): 776-778.

CAO Jun, YUAN Shenfang, CAI Jian, et al. Real-time structural health monitoring for fatigue crack growth[J]. Piezoelectrics and Acoustooptics, 2008, 30(6): 776-778.

[16] 袁慎芳, 邱雷, 王強, 等．壓電-光纖綜合結構健康監測系統的研究及驗證[J]．航空學報, 2009, 30(2): 348-356.

YUAN Shenfang, QIU Lei, WANG Qiang, et al. Application research of a hybrid piezoelectric-optic fiber integrated structural health monitoring system[J]. Acta Aeronautica et Astronautica Sinica, 2009, 30(2): 348-356.

[17] 陸希,孟光,李富才．基于Lamb波的薄壁槽狀結構損傷檢測研究[J]．振動與沖擊, 2012, 31(12):63-67.

LU Xi, MENG Guang, LI Fucai, et al. Lamb wave-based damage detection for a channel-like thin-wall structure[J]. Journal of Vibration and Shock, 2012, 31(12):63-67.

[18] 邱雷, 袁慎芳．集成壓電健康監測掃查系統的研制及其應用[J]．壓電與聲光, 2008, 30(1): 39-41.

QIU Lei, YUAN Shenfang. Research and application of integrated health monitoring scanning system based on PZT sensor array[J]. Piezoelectrics & Acoustooptics, 2008,30(1): 39-41.

[19] QIU L, YUAN S F, WANG Q, et al. Design and experiment of PZT network-based structural health monitoring scanning system[J]. Chinese Journal of Aeronautics, 2009, 22(5): 505-512.

[20] DONG M, WANG N. Adaptive network-based fuzzy inference system with leave-one-out cross-validation approach for prediction of surface roughness[J]. Applied Mathematical Modelling, 2011, 35(3): 1024-1035.