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

背景技術(shù)：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

附圖說明

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

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

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

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

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

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

具體實施方式

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

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

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

(2)根據(jù)結(jié)構(gòu)的實際情況定義疲勞裂紋擴(kuò)展?fàn)顟B(tài)空間模型，包括狀態(tài)方程和觀測方程兩個部分，其中狀態(tài)方程由結(jié)構(gòu)的形狀及其服役條件確定，觀測方程由基于導(dǎo)波的結(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)是裂紋擴(kuò)展增量，為x_k-1的函數(shù)，由疲勞裂紋擴(kuò)展模型定義，比如Paris模型，NASGRO模型；ω是一個隨機(jī)變量，服從高斯分布表征疲勞裂紋擴(kuò)展的不確定性，其均值選擇為以保證exp(ω)的期望為1，指數(shù)項exp(ω)保證裂紋擴(kuò)展增量非負(fù)；y_k為k時刻獲得的導(dǎo)波損傷因子；g(·)是觀測映射，表征裂紋長度和損傷因子之間的關(guān)系；v是一個隨機(jī)變量，表征觀測的不確定性。

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

式中，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對應(yīng)的裂紋長度g^-1(y_k)；由上一時刻粒子集通過狀態(tài)方程一步轉(zhuǎn)移得到先驗粒子集以的標(biāo)準(zhǔn)差作為觀測概率密度的標(biāo)準(zhǔn)差。

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

(a)令i＝1

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

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

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

(c)i＝i+1，并且重復(fù)步驟(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)值分布，其均值為裂紋長度的先驗估計標(biāo)準(zhǔn)差為σ_p，該條件概率密度表征粒子與先驗估計之間的高斯權(quán)值。先驗估計由先驗粒子集加權(quán)求和得到，比例項(σ_p/σ_v)保證兩種不同加權(quán)方式得到的粒子權(quán)值在同一數(shù)量級，式中σ_v為表征觀測不確定性的概率密度標(biāo)準(zhǔn)差。

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

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

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

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

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

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

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

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

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

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

(c)令i＝1

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

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

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

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

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

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

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

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

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

假設(shè)導(dǎo)波損傷因子與實際的裂紋長度存在如下關(guān)系，

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

式中，x為裂紋長度；y為導(dǎo)波損傷因子；x₀為初始裂紋長度；觀測過程受到高斯白噪聲v的影響，其服從高斯分布N(0,0.03²)；并且每隔15000個載荷循環(huán)進(jì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)先驗定義的裂紋擴(kuò)展?fàn)顟B(tài)方程不準(zhǔn)確的情形，其余參數(shù)保持不變；表征裂紋擴(kuò)展不確定性的隨機(jī)變量ω標(biāo)準(zhǔn)差為σ_ω＝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為導(dǎo)波損傷因子；x₀為初始裂紋長度；觀測不確定性v服從高斯分布N(0,0.03²)，其標(biāo)準(zhǔn)差σ_v＝0.03。

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

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