本發(fā)明涉及一種基于高斯權(quán)值-混合粒子濾波的疲勞裂紋擴展預測方法，屬于故障預測與健康管理技術(shù)領域。

背景技術(shù)：

隨著現(xiàn)代工程系統(tǒng)的復雜程度和綜合化程度越來越高，傳統(tǒng)的事后維修和定期維護策略已經(jīng)無法滿足其保障維護的要求。近年來，故障預測與健康管理(Prognostics and Health Management，PHM)技術(shù)可以通過傳感器在線獲取系統(tǒng)的實際健康狀態(tài)，并預測其退化情況，從而制定最佳的運行和維護策略，以實現(xiàn)視情維護。因此PHM技術(shù)受到越來越多的關注。工程結(jié)構(gòu)作為系統(tǒng)的核心組成部分，其服役環(huán)境復雜并承受各種交變載荷，容易產(chǎn)生疲勞裂紋。疲勞裂紋的存在和擴展將嚴重削弱結(jié)構(gòu)承載能力，甚至導致災難性的事故的發(fā)生。因此對結(jié)構(gòu)的疲勞裂紋進行在線準確地預測具有重要的理論意義和工程應用價值。

然而疲勞裂紋擴展過程是一個包含各種不確定性的過程，比如材料性質(zhì)的不確定性，服役環(huán)境的不確定性，以及載荷的不確定性等。同時基于結(jié)構(gòu)健康監(jiān)測方法的在線裂紋監(jiān)測也受到各種不確定性因素的影響。這些不確定性因素嚴重影響了疲勞裂紋擴展預測的準確性。近年來，基于貝葉斯理論的概率方法通過結(jié)合裂紋擴展模型與實際的結(jié)構(gòu)健康監(jiān)測數(shù)據(jù)以消除這些不確定性的影響，受到越來越多的關注。在結(jié)構(gòu)健康監(jiān)測方法中，基于導波的結(jié)構(gòu)健康監(jiān)測方法具有小損傷敏感，監(jiān)測范圍廣等優(yōu)點，被認為是非常有前景的方法之一。

在基于貝葉斯理論的方法中，粒子濾波方法由于不需要滿足線性和高斯過程假設，被認為非常適用于解決非線性非高斯的疲勞裂紋擴展預測問題。但將粒子濾波方法應用于疲勞裂紋擴展預測時，常常難以準確地定義先驗裂紋擴展狀態(tài)方程，也就是描述當前結(jié)構(gòu)疲勞裂紋擴展的模型，這會加劇粒子濾波算法的粒子退化和多樣性匱乏問題，使得裂紋擴展預測結(jié)果存在較大的誤差。

技術(shù)實現(xiàn)要素：

本發(fā)明針對先驗定義的疲勞裂紋擴展狀態(tài)方程通常與結(jié)構(gòu)實際的裂紋擴展存在較大偏差的問題，提出了一種基于高斯權(quán)值-混合粒子濾波的疲勞裂紋擴展預測方法，結(jié)合基于主動導波的裂紋監(jiān)測方法實現(xiàn)結(jié)構(gòu)疲勞裂紋擴展的準確預測。

本發(fā)明為解決其技術(shù)問題采用如下技術(shù)方案：

一種基于高斯權(quán)值-混合粒子濾波的疲勞裂紋擴展預測方法，包括如下步驟：

(1)在結(jié)構(gòu)處于服役條件下，采用基于主動導波的結(jié)構(gòu)健康監(jiān)測方法對結(jié)構(gòu)的疲勞裂紋進行在線監(jiān)測，通過獲取的導波信號計算損傷因子，并以損傷因子作為裂紋長度的觀測值；同時根據(jù)結(jié)構(gòu)的實際服役情況定義結(jié)構(gòu)的疲勞裂紋擴展狀態(tài)空間模型；

(2)在結(jié)構(gòu)處于服役條件下，每當獲得一個新的裂紋長度觀測值，按照定義的混合系數(shù)λ，從均勻分布U(0,1)中隨機采樣得到隨機數(shù)，如果該隨機數(shù)大于λ，從觀測概率密度中隨機采樣得到一個表征裂紋長度的粒子，否則從先驗轉(zhuǎn)移概率密度中隨機采樣得到粒子，重復該過程的得到確定數(shù)目粒子組成的粒子集；

(3)對于從先驗轉(zhuǎn)移概率密度中采樣得到的粒子，通過觀測似然概率密度加權(quán)；而對于從觀測概率密度中采樣得到的粒子，根據(jù)先驗估計定義一個高斯權(quán)值分布，通過粒子與先驗估計之間的高斯權(quán)值作為粒子的權(quán)值，然后對所有粒子的權(quán)值進行歸一化，得到粒子的歸一化權(quán)值；

(4)通過粒子集和相應的歸一化權(quán)值表征裂紋長度的后驗概率密度，并計算裂紋長度的后驗估計；在獲得的后驗概率密度基礎上，將粒子向未來時刻投影得到當前時刻的疲勞裂紋擴展預測結(jié)果；

(5)根據(jù)每個粒子的歸一化權(quán)值，使用系統(tǒng)重采樣算法進行重采樣；

(6)重復上述步驟(2)，(3)，(4)，(5)，(6)實現(xiàn)在線的疲勞裂紋擴展預測。

步驟(2)所述混合系數(shù)λ取值為0.5；此外，觀測概率密度定義為p(x_k|y_k)，式中y_k為k時刻獲得的裂紋長度觀測值，x_k為k時刻的裂紋長度；觀測概率密度的概率密度形式與表征觀測不確定性的概率密度一致，以觀測值y_k對應的裂紋長度為均值，k時刻先驗粒子集的標準差作為觀測概率密度的標準差。

步驟(3)中定義了一個高斯權(quán)值分布，對從觀測概率密度中采樣得到的粒子加權(quán)，權(quán)值計算方式如下，

式中，是從觀測概率密度p(x_k|y_k)中采樣得到的第i個粒子，為第i個粒子在k-1時刻的權(quán)值；為第i個粒子在k時刻的權(quán)值；σ_v為表征觀測不確定性的概率密度標準差；為定義的高斯權(quán)值分布，其均值為裂紋長度的先驗估計標準差σ_p為經(jīng)驗設定值。

本發(fā)明的有益效果如下：

本發(fā)明提出了一種基于高斯權(quán)值-混合粒子濾波的疲勞裂紋擴展預測方法，以觀測概率密度和先驗轉(zhuǎn)移概率密度的線性組合作為粒子濾波的重要性密度函數(shù)。從重要性密度函數(shù)中采樣得到的粒子集同時考慮了觀測信息和先驗信息，降低了對準確定義疲勞裂紋擴展狀態(tài)方程的依賴，同時結(jié)合了基于主動導波的裂紋監(jiān)測方法，能有效地用于結(jié)構(gòu)的疲勞裂紋擴展預測。

附圖說明

圖1為本發(fā)明提出的疲勞裂紋擴展預測方法流程圖。

圖2為實施例中的結(jié)構(gòu)尺寸圖。

圖3為實施例中的結(jié)構(gòu)疲勞裂紋擴展圖。

圖4為實施例中結(jié)構(gòu)健康監(jiān)測方法獲得的導波損傷因子圖。

圖5為實施例中裂紋長度后驗估計對比圖。

圖6為實施例中預測的失效循環(huán)載荷數(shù)的相對誤差對比圖。

具體實施方式

下面結(jié)合附圖對本發(fā)明創(chuàng)造的技術(shù)方案進行詳細說明。

圖1所示，為一種基于高斯權(quán)值-混合粒子濾波的疲勞裂紋擴展預測方法的流程圖，方法步驟如下，

(1)在結(jié)構(gòu)處于服役條件下，采用基于主動導波的結(jié)構(gòu)健康監(jiān)測方法對結(jié)構(gòu)的疲勞裂紋進行在線監(jiān)測。在結(jié)構(gòu)的關鍵部位粘貼壓電傳感器，通過壓電傳感器激勵和采集結(jié)構(gòu)中的導波信號。利用結(jié)構(gòu)健康監(jiān)測系統(tǒng)實現(xiàn)導波信號的在線采集和調(diào)理，對獲取的導波信號，計算其損傷因子作為裂紋長度的觀測值；

(2)根據(jù)結(jié)構(gòu)的實際情況定義疲勞裂紋擴展狀態(tài)空間模型，包括狀態(tài)方程和觀測方程兩個部分，其中狀態(tài)方程由結(jié)構(gòu)的形狀及其服役條件確定，觀測方程由基于導波的結(jié)構(gòu)健康監(jiān)測方法確定，如下所示。

狀態(tài)方程：x_k＝x_k-1+Δx_k-1(x_k-1)·exp(ω)

觀測方程：y_k＝g(x_k)+v

式中，k是離散的時間；x_k是k時刻的裂紋長度；x_k-1是k-1時刻的裂紋長度；Δx_k-1(x_k-1)是裂紋擴展增量，為x_k-1的函數(shù)，由疲勞裂紋擴展模型定義，比如Paris模型，NASGRO模型；ω是一個隨機變量，服從高斯分布表征疲勞裂紋擴展的不確定性，其均值選擇為以保證exp(ω)的期望為1，指數(shù)項exp(ω)保證裂紋擴展增量非負；y_k為k時刻獲得的導波損傷因子；g(·)是觀測映射，表征裂紋長度和損傷因子之間的關系；v是一個隨機變量，表征觀測的不確定性。

(3)定義狀態(tài)空間模型后，通過結(jié)構(gòu)健康監(jiān)測方法在線獲取結(jié)構(gòu)當前時刻k的裂紋長度觀測值y_k，每當獲得一個新的裂紋長度觀測值，從如下所示的混合建議分布中采樣，

式中，k是離散的時間；x_k是k時刻的裂紋長度；x_k-1是k-1時刻的裂紋長度，y_k是k時刻的裂紋長度觀測值；i為粒子序號，i＝1,...,N_s，N_s為粒子總數(shù)；為重要性密度函數(shù)；為先驗轉(zhuǎn)移概率密度；λ為混合系數(shù)取為λ＝0.5；p(x_k|y_k)為觀測概率密度，其形式與表征觀測不確定性的概率密度一致。

對于觀測概率密度，其均值為觀測值y_k對應的裂紋長度g^-1(y_k)；由上一時刻粒子集通過狀態(tài)方程一步轉(zhuǎn)移得到先驗粒子集以的標準差作為觀測概率密度的標準差。

從混合建議分布中采樣得到粒子集的方式如下：

(a)令i＝1

(b)從均勻分布U(0,1)中隨機采樣得到u⁽ⁱ⁾，

如果u⁽ⁱ⁾＞λ，從p(x_k|y_k)中隨機采樣得到

如果u⁽ⁱ⁾≤λ，從中隨機采樣得到

(c)i＝i+1，并且重復步驟(a)(b)，直到i＝N_s

(4)采樣得到粒子集后，分別計算每個粒子的權(quán)值。從先驗概率分布中采樣得到的樣本，其權(quán)值為，

式中，為第i個粒子在k-1時刻的權(quán)值；為第i個粒子在k時刻的權(quán)值；為觀測似然概率密度。

如果該粒子是從觀測概率密度p(x_k|y_k)中采樣得到的，權(quán)值計算如下，

式中，為第i個粒子在k-1時刻的權(quán)值；為第i個粒子在k時刻的權(quán)值；為定義的高斯權(quán)值分布，其均值為裂紋長度的先驗估計標準差為σ_p，該條件概率密度表征粒子與先驗估計之間的高斯權(quán)值。先驗估計由先驗粒子集加權(quán)求和得到，比例項(σ_p/σ_v)保證兩種不同加權(quán)方式得到的粒子權(quán)值在同一數(shù)量級，式中σ_v為表征觀測不確定性的概率密度標準差。

計算得到所有粒子的權(quán)值后，對所有粒子的權(quán)值進行歸一化處理，如下式所示

當前時刻裂紋長度的后驗概率密度由粒子集以及對應的歸一化權(quán)值表征如下，

式中，p(x_k|y_1:k)為當前時刻裂紋長度的后驗概率密度；δ是狄利克雷函數(shù)，其表達式如下，

當前時刻裂紋長度的后驗估計由粒子加權(quán)求和得到，

(5)在獲得的粒子集和歸一化權(quán)值的基礎上，通過裂紋擴展狀態(tài)空間模型預測每個粒子在未來時刻的疲勞裂紋擴展，得到粒子集里面所有粒子的裂紋擴展預測結(jié)果，其表達式如下，

式中為第i個粒子在未來時刻k+d的裂紋長度，為第i個粒子在未來時刻k+d-1的裂紋長度(d＝1,2,3,......)。最后得到預測的裂紋長度概率密度以及裂紋長度預測值如下，

式中，為k時刻的粒子權(quán)值；p(x_k+d|y_1:k)為預測的k+d時刻裂紋長度概率密度，粒子為預測的k+d時刻裂紋長度

(6)根據(jù)歸一化權(quán)值對所有的粒子進行重采樣，采用系統(tǒng)重采樣算法，步驟如下：

(a)計算歸一化權(quán)值的累積數(shù)列C，共N_s個元素。

(b)從均勻分布U(0,1)中隨機采樣得到u，

(c)令i＝1

計算u⁽ⁱ⁾＝u+(i-1)/N_s

順序查找數(shù)列C中第一個大于u⁽ⁱ⁾的元素下標j，

令重采樣得到的粒子并且權(quán)值

(d)i＝i+1，并重復步驟(c)，直到i＝N_s

(7)系統(tǒng)重采樣算法完成后，將重采樣得到的粒子和相應的權(quán)值作為當前值，繼續(xù)通過結(jié)構(gòu)健康監(jiān)測方法獲得新的損傷因子，k＝k+1，算法循環(huán)迭代。

本實施例中以仿真的AL2024-T6鋁板的中心裂紋擴展為例來具體說明本發(fā)明方法的具體實施過程。

如圖2所示的中心裂紋鋁板，其厚度為3mm，尺寸為b＝200mm,h＝400mm。所承受疲勞載荷為常幅正弦載荷，載荷峰值S_max＝60Mpa，應力比為R＝0。通過下式所述的NASGRO模型仿真結(jié)構(gòu)真實的疲勞裂紋擴展，

式中，k為離散的時刻，x_k為k時刻的裂紋長度，x_k-1為k-1時刻的裂紋長度。ΔK為應力強度因子幅；K_max為應力強度因子最大值；ΔN為載荷循環(huán)的步長，C，m，p，q為經(jīng)驗常數(shù)；K_C為斷裂韌度；ΔK_th為應力強度因子幅閾值；f為裂紋張開系數(shù)；應力強度因子式中形狀參數(shù)b為矩形板寬度，S為矩形板承受的拉伸應力。

取ΔN＝50個載荷循環(huán)；f＝0.3958；C＝exp(-26.76)；m＝3.2；p＝0.25；q＝1；初始裂紋長度x₀＝3mm，得到疲勞裂紋擴展軌跡如圖3所示。

假設導波損傷因子與實際的裂紋長度存在如下關系，

y＝-5×10^-5(x-x₀)³+2×10^-3(x-x₀)²+1.5×10^-2(x-x₀)+v

式中，x為裂紋長度；y為導波損傷因子；x₀為初始裂紋長度；觀測過程受到高斯白噪聲v的影響，其服從高斯分布N(0,0.03²)；并且每隔15000個載荷循環(huán)進行裂紋長度觀測，得到的表征裂紋長度的損傷因子如圖4所示。

首先根據(jù)具體情況定義該結(jié)構(gòu)的狀態(tài)空間模型，包括狀態(tài)方程和觀測方程。其中狀態(tài)方程定義如下，

式中，x_k為k時刻的裂紋長度；x_k-1為k-1時刻的裂紋長度；參數(shù)C選取為C＝exp(-26.6)以體現(xiàn)先驗定義的裂紋擴展狀態(tài)方程不準確的情形，其余參數(shù)保持不變；表征裂紋擴展不確定性的隨機變量ω標準差為σ_ω＝1.3。

觀測方程定義根據(jù)觀測映射定義，如下式所示。

y_k＝-5×10^-5(x_k-x₀)³+2×10^-3(x_k-x₀)²+1.5×10^-2(x_k-x₀)+v

式中，y_k為導波損傷因子；x₀為初始裂紋長度；觀測不確定性v服從高斯分布N(0,0.03²)，其標準差σ_v＝0.03。

接下來對高斯權(quán)值-混合粒子濾波算法進行初始化，選取粒子數(shù)N_s＝500，高斯權(quán)值分布的標準差σ_p＝0.4。在初始時刻k＝0所有的粒子的裂紋長度都初始化為x₀。

在結(jié)構(gòu)處于服役過程中，通過基于主動導波的結(jié)構(gòu)健康監(jiān)測方法對結(jié)構(gòu)的裂紋長度進行在線監(jiān)測，序貫獲得導波損傷因子。在當前時刻，一旦獲得新的損傷因子，粒子濾波方法首先計算裂紋長度的后驗概率密度，以及裂紋長度的后驗估計，然后在該后驗概率密度的基礎上對未來的疲勞裂紋擴展進行預測。圖5所示為結(jié)合了所有觀測值后的裂紋長度后驗估計結(jié)果，其中先驗模型表征的是裂紋擴展狀態(tài)方程的裂紋擴展期望值。可以看到，先驗模型的誤差很大。而本發(fā)明的方法在裂紋長度觀測點的后驗估計均方根誤差僅為0.28mm。另一方面，假設以裂紋長度18.69mm為臨界裂紋長度，超過該裂紋長度的循環(huán)載荷數(shù)被認為是失效循環(huán)載荷數(shù)。本發(fā)明提出的方法預測失效循環(huán)載荷數(shù)的相對誤差如圖6所示，可以看到，本文提出的方法可以較為準確的預測結(jié)構(gòu)疲勞裂紋擴展。