本發(fā)明涉及一種基于數(shù)學(xué)形態(tài)學(xué)和IFOA-SVR的滾動軸承可靠度預(yù)測方法，涉及滾動軸承可靠度預(yù)測技術(shù)領(lǐng)域。

背景技術(shù)：

滾動軸承是旋轉(zhuǎn)機械中的關(guān)鍵部件，一旦出現(xiàn)故障將會造成大量的經(jīng)濟損失甚至危害人的生命安全^[1-2]。因此，準(zhǔn)確預(yù)測滾動軸承在下一階段的工作狀態(tài)是合理制定機械設(shè)備維修計劃的前提基礎(chǔ)^[3-4]。

目前，滾動軸承振動信號的特征提取方法研究受到學(xué)者的廣泛關(guān)注。文獻[5]提出了基于形態(tài)分量分析和包絡(luò)譜的滾動軸承故障診斷方法，該方法可有效提取滾動軸承振動信號中的故障特征。文獻[6]提出局部均值分解與形態(tài)學(xué)分形維數(shù)相結(jié)合的滾動軸承故障診斷方法，可有效地對滾動軸承進行故障診斷。本發(fā)明將Hilbert包絡(luò)解調(diào)與形態(tài)學(xué)分形維數(shù)結(jié)合，對滾動軸承全壽命數(shù)據(jù)提取出數(shù)學(xué)形態(tài)學(xué)分形維數(shù)，并將其作為性能退化狀態(tài)特征。

在建立支持向量回歸(Support vector regression,SVR)模型方面，文獻[7]利用粒子群算法(Particle swarm optimization,PSO)優(yōu)化最小二乘SVR對滾動軸承的性能退化趨勢進行了20點預(yù)測，平均絕對百分誤差(Mean absolute percentage error,MAPE)為15.82％。文獻[8]利用交叉驗證法對C，g尋優(yōu)建立SVR模型，對滾動軸承故障趨勢進行了7點預(yù)測，MAPE為0.2％。文獻[9]通過經(jīng)驗選取C，g參數(shù)建立SVR模型，并對電廠汽輪發(fā)電機組轉(zhuǎn)子的振動幅值進行了30點預(yù)測，平均絕對誤差(Mean absolute error,MAE)為0.6424，MAPE為1.48％。文獻[10]利用PSO對SVR模型中的C，g以及ε尋優(yōu)，預(yù)測航空飛行器的剩余壽命。文獻[7]到文獻[9]預(yù)測步長較短，文獻[8]和[9]只對SVR模型中的C和g進行選取，并沒有提出對不敏感誤差ε的取值問題。同時，憑經(jīng)驗選取C和g有很大的盲目性，通過交叉對比測試確定C和g的值則需要進行大量的實驗，而PSO相對復(fù)雜，并且容易陷入局部最優(yōu)，導(dǎo)致對模型參數(shù)的賦值不夠準(zhǔn)確。文獻[11]提出一種果蠅優(yōu)化算法(Fruit fly optimization algorithm,FOA)。文獻[12]提出遞減步長果蠅優(yōu)化算法(Diminishing step fruit fly optimization algorithm,DS-FOA)實現(xiàn)了全局搜索能力和局部尋優(yōu)能力的平衡。但DS-FOA在二維空間進行搜索，不能真實反映果蠅的覓食行為。SVR模型內(nèi)部參數(shù)的選取會對SVR的預(yù)測精度以及預(yù)測步長產(chǎn)生較大影響，本發(fā)明采用在三維空間搜索的改進果蠅優(yōu)化算法(Improved fruit fly optimization algorithm,IFOA)對SVR模型中的3個參數(shù)C、g以及ε同時進行參數(shù)尋優(yōu)。并且在進行預(yù)測時，將MAE、均方根誤差(Root mean square error,RMSE)、歸一化均方誤差(Normalized mean square error,NMSE)以及MAPE的和作為適應(yīng)度函數(shù)，相比于任一誤差作為適應(yīng)度函數(shù)具有更強的尋優(yōu)能力。

在可靠度模型方面，可靠度理論已被應(yīng)用于機械工程以及航空航天等各個領(lǐng)域^[13]。比例故障率模型的優(yōu)勢在于將軸承的性能退化特征與可靠度理論相結(jié)合，這樣，軸承的可靠度評估就能夠在設(shè)備狀態(tài)監(jiān)測數(shù)據(jù)的基礎(chǔ)上得到更新。文獻[14]提出將峭度、均方根作為特征指標(biāo)，利用威布爾比例故障率模型(Weibull proportional hazard model,WPHM)作為可靠度評估模型，并利用fminsearch優(yōu)化函數(shù)求解極大似然方程組確定WPHM的待定參數(shù)，有效地對鐵路機車輪對軸承進行可靠度評估。文獻[15]在估計WPHM的待定參數(shù)時，利用牛頓迭代法對極大似然方程組進行求解。但上述兩種求解方法都需要根據(jù)經(jīng)驗設(shè)定搜索初值，并且運算時間較長。本發(fā)明將極大似然方程組中每個方程的絕對值的和作為IFOA的適應(yīng)度函數(shù)，向適應(yīng)度函數(shù)最小的方向優(yōu)化，可以快速求出方程組的解，確定WPHM的待定參數(shù)。

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

本發(fā)明為了保證滾動軸承運行狀態(tài)可靠度的預(yù)測精度同時增加預(yù)測步長，提出一種數(shù)學(xué)形態(tài)學(xué)分形維數(shù)結(jié)合改進果蠅優(yōu)化算法-支持向量回歸(Improved fruit fly optimization algorithm-support vector regression,IFOA-SVR)的滾動軸承可靠度預(yù)測方法，即提出一種基于數(shù)學(xué)形態(tài)學(xué)和IFOA-SVR的滾動軸承可靠度預(yù)測方法。

本發(fā)明為解決上述技術(shù)問題采取的技術(shù)方案是：

一種基于數(shù)學(xué)形態(tài)學(xué)和IFOA-SVR的滾動軸承可靠度預(yù)測方法，所述方法的實現(xiàn)過程為：

步驟一、求出包絡(luò)信號的數(shù)學(xué)形態(tài)學(xué)分形維數(shù)D_M；

步驟二、從獲取的數(shù)學(xué)形態(tài)學(xué)分形維數(shù)D_M選取訓(xùn)練樣本對，基于所述訓(xùn)練樣本對利用IFOA對SVR模型中的參數(shù)C、g以及ε同時進行尋優(yōu)，建立IFOA-SVR預(yù)測模型：

構(gòu)建IFOA的過程：

(1)初始化算法參數(shù)：設(shè)置果蠅種群規(guī)模Sizepop，最大覓食代數(shù)Maxgen，并隨機初始化果蠅群體位置坐標(biāo)(X₀,Y₀,Z₀)；

(2)果蠅個體利用嗅覺隨機搜索的方向和距離可以通過式(12)獲得

式中，i＝1,2,...,Sizepop，L₀為初始步長值，gen為當(dāng)前覓食代數(shù)；

(3)由于無法確定食物源的具體位置，所以需要通過式(13)估計第i個果蠅個體的當(dāng)前位置與坐標(biāo)原點間的距離Dist_i，之后計算出味道濃度判定值S_i

S_i＝1/Dist_i (14)

(4)將S_i代入味道濃度判定函數(shù)，計算出果蠅個體當(dāng)前位置的味道濃度

Smell_i＝function(S_i) (15)

(5)當(dāng)前果蠅群體中具有最高味道濃度的個體，可由式(16)獲得

[bestSmell,bestIndex]＝max(Smell_i) (16)

式中，bestSmell表示果蠅群體中具有最高味道濃度的個體的味道濃度值，bestIndex表示果蠅群體中具有最高味道濃度的個體的位置；

(6)保留果蠅群體中最佳味道濃度值和與其對應(yīng)的個體坐標(biāo)，同時果蠅群體利用自身的視覺對食物源進行定位，然后飛往食物源所在的位置；

(7)進入迭代尋優(yōu)過程，重復(fù)步驟(2)-(5)，并判斷當(dāng)前味最高道濃度是否好于前一迭代味道濃度，且gen<Maxgen；若成立，則執(zhí)行步驟(6)；

建立IFOA-SVR模型，其過程為：

(1)初始化IFOA參數(shù)，包括Sizepop、Maxgen以及果蠅個體初始位置；對SVR中的3個參數(shù)：懲罰系數(shù)C、核函數(shù)參數(shù)g以及不敏感誤差ε進行尋優(yōu)，所以初始坐標(biāo)為(X₀¹,Y₀¹,Z₀¹)，(X₀²,Y₀²,Z₀²)和(X₀³,Y₀³,Z₀³)；

(2)附與每個果蠅個體隨機方向和飛行距離，并用遞減搜索步長代替固定步長，得到(X_i¹,Y_i¹,Z_i¹)，(X_i²,Y_i²,Z_i²)，(X_i³,Y_i³,Z_i³)，計算當(dāng)前果蠅個體與原點間距離的倒數(shù)，得到味道濃度判定值S_i¹，S_i²和S_i³；

(3)確定SVR中參數(shù)C、g和ε的范圍，即C∈[2^-14,2¹⁴]，g∈[2^-14,2¹⁴]，ε∈[2^-14,2¹⁴]；

(4)將所述訓(xùn)練樣本對輸入到SVR模型中，對模型進行訓(xùn)練，將平均絕對誤差MAE、均方根誤差RMSE、歸一化均方誤差NMSE以及平均絕對百分誤差MAPE的和作為適應(yīng)度函數(shù)，即

(5)找到適應(yīng)度函數(shù)Fitness的最小值對應(yīng)的果蠅個體，開始迭代尋優(yōu)，并判斷最小Fitness是否低于前一代最小Fitness；如低于，則保留最小Fitness值及其對應(yīng)的坐標(biāo)，并將其賦給初始坐標(biāo)；如高于，則返回步驟(2)；

(6)找到C、g和ε的最佳值，建立IFOA-SVR預(yù)測模型；

步驟三、從所述性能退化狀態(tài)特征中截取一定長度Z作為WPHM的響應(yīng)協(xié)變量，結(jié)合極大似然估計獲得似然函數(shù)方程組，把方程組中每個方程絕對值的和作為IFOA的適應(yīng)度函數(shù)，求出方程組的解，確定WPHM的待定參數(shù)，進而得到可靠度模型R(t,Z)，過程如下：

基于WPHM建立滾動軸承性能退化狀態(tài)特征與可靠度之間的數(shù)學(xué)關(guān)系，WPHM的表達式為

式中β為形狀參數(shù)，η為尺度參數(shù)，μ為協(xié)變量回歸參數(shù)，t表示時間，Z為截取一定長度的性能退化狀態(tài)特征，即從數(shù)學(xué)形態(tài)學(xué)分形維數(shù)D_M中截取一定長度數(shù)據(jù)；

h(t,Z)和可靠度函數(shù)R(t,Z)之間的關(guān)系為

從而，R(t,Z)可表示為

利用極大似然估計法得到關(guān)于待定參數(shù)β、η、μ的方程組并利用IFOA對該方程組進行求解，確定WPHM的表達式h(t,Z)的待定參數(shù)β、η、μ的最終取值；將參數(shù)β、η、μ代入式(21)即可建立可靠度模型；

步驟四、將性能退化狀態(tài)特征的最后一部分?jǐn)?shù)據(jù)作為IFOA-SVR預(yù)測模型的輸入，采用長期迭代預(yù)測法獲取性能退化狀態(tài)特征預(yù)測結(jié)果；

步驟五、將IFOA-SVR預(yù)測模型得出的預(yù)測結(jié)果作為響應(yīng)協(xié)變量嵌入到可靠度模型R(t,Z)中，即可計算出所述性能退化狀態(tài)特征預(yù)測結(jié)果所對應(yīng)的可靠度，實現(xiàn)滾動軸承可靠度的預(yù)測。

在步驟一中，通過Hilbert變換計算出狀態(tài)監(jiān)測數(shù)據(jù)的包絡(luò)信號；

步驟一一、獲取滾動軸承的狀態(tài)監(jiān)測數(shù)據(jù)：滾動軸承的振動信號；

步驟一二、通過Hilbert變換計算出狀態(tài)監(jiān)測數(shù)據(jù)的包絡(luò)信號；

步驟一三、利用式(4)求出每次采集振動信號的包絡(luò)信號的數(shù)學(xué)形態(tài)學(xué)分形維數(shù)D_M，將D_M作為滾動軸承性能退化狀態(tài)特征；

式中：D_M為信號的Minkowski-Bouligand維數(shù)即數(shù)學(xué)形態(tài)學(xué)分形維數(shù)，c為常數(shù)，對和log(1/λ)進行最小二乘線性擬合即可得到對信號數(shù)學(xué)形態(tài)學(xué)分形維數(shù)D_M；λ＝1,2,…,λ_max，λ_max為最大分析尺度；

A_Se(λ)為定義在分析尺度λ下f(n)關(guān)于Se(n)膨脹和腐蝕運算的覆蓋面積，f(n)表示滾動軸承振動信號的包絡(luò)信號；Se(n)為結(jié)構(gòu)元素，為一維離散向量；Se(n)的長度及λ_max的值通過控制變量法求得。

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

本發(fā)明提出的一種基于數(shù)學(xué)形態(tài)學(xué)分形維數(shù)結(jié)合改進果蠅算法-支持向量回歸(Improved fruit fly optimization algorithm-support vector regression,IFOA-SVR)的滾動軸承可靠度預(yù)測方法，充分發(fā)揮IFOA-SVR預(yù)測模型的優(yōu)勢，對提取的滾動軸承狀態(tài)特征數(shù)據(jù)進行預(yù)測。將預(yù)測結(jié)果嵌入到已建立的可靠度模型中進行滾動軸承可靠度預(yù)測，在保證滾動軸承運行狀態(tài)可靠度預(yù)測精度的同時增加了預(yù)測步長。

該方法首先提取振動信號的包絡(luò)信號，計算該包絡(luò)信號的數(shù)學(xué)形態(tài)學(xué)分形維數(shù)，將其作為滾動軸承性能退化狀態(tài)特征；其次，利用IFOA對SVR中的參數(shù)C，g以及ε同時進行尋優(yōu)，建立IFOA-SVR預(yù)測模型。同時，利用極大似然估計結(jié)合IFOA建立威布爾比例故障率模型(Weibull proportional hazard model,WPHM)，進而得到可靠度模型；最后，將退化狀態(tài)特征作為IFOA-SVR預(yù)測模型的輸入，采用長期迭代預(yù)測法獲取特征預(yù)測結(jié)果，并將該結(jié)果嵌入到可靠度模型中，從而預(yù)測出軸承運行狀態(tài)的可靠度。實驗結(jié)果表明，利用所提方法對滾動軸承可靠度進行預(yù)測，能在保證預(yù)測精度的前提下增加預(yù)測步長。

具體而言，本發(fā)明將Hilbert包絡(luò)解調(diào)與形態(tài)學(xué)分形維數(shù)結(jié)合，對滾動軸承全壽命數(shù)據(jù)提取出數(shù)學(xué)形態(tài)學(xué)分形維數(shù)，并將其作為性能退化狀態(tài)特征。本發(fā)明采用在三維空間搜索的改進果蠅優(yōu)化算法(Improved fruit fly optimization algorithm,IFOA)對SVR模型中的3個參數(shù)C、g以及ε同時進行參數(shù)尋優(yōu)。并且在進行預(yù)測時，將MAE、均方根誤差(Root mean square error,RMSE)、歸一化均方誤差(Normalized mean square error,NMSE)以及MAPE的和作為適應(yīng)度函數(shù)，相比于任一誤差作為適應(yīng)度函數(shù)具有更強的尋優(yōu)能力。本發(fā)明將極大似然方程組中每個方程的絕對值的和作為IFOA的適應(yīng)度函數(shù)，向適應(yīng)度函數(shù)最小的方向優(yōu)化，可以快速求出方程組的解，確定WPHM的待定參數(shù)。

附圖說明

圖1是本發(fā)明中的果蠅三維空間覓食示意圖，圖2是本發(fā)明的可靠度預(yù)測方法流程圖；

圖3為本發(fā)明方法在應(yīng)用中的第860次采集數(shù)據(jù)的時域波形圖，圖4是本發(fā)明方法在應(yīng)用中對應(yīng)的第860次采集數(shù)據(jù)的包絡(luò)信號波形圖，圖5為本發(fā)明在應(yīng)用中第1點至第860點的可靠度曲線圖，圖6為本發(fā)明在應(yīng)用中GA優(yōu)化SVR模型得出的預(yù)測結(jié)果圖，圖7為本發(fā)明在應(yīng)用中PSO優(yōu)化SVR模型得出的預(yù)測結(jié)果圖，圖8為本發(fā)明在應(yīng)用中DS-FOA優(yōu)化SVR模型得出的預(yù)測結(jié)果圖，圖9為本發(fā)明在應(yīng)用中IFOA優(yōu)化SVR模型得出的預(yù)測結(jié)果圖，圖10為本發(fā)明在應(yīng)用中50點預(yù)測數(shù)據(jù)的可靠度趨勢曲線圖；

圖11為由包絡(luò)求分形維數(shù)與直接求分形維數(shù)的對比圖，本發(fā)明采用由由包絡(luò)求分形維數(shù)。

具體實施方式

具體實施方式一：如圖1至11所示，本實施方式針對所述的基于數(shù)學(xué)形態(tài)學(xué)和IFOA-SVR的滾動軸承可靠度預(yù)測方法的實現(xiàn)過程具體描述如下：

1基于數(shù)學(xué)形態(tài)學(xué)的分形維數(shù)

數(shù)學(xué)形態(tài)學(xué)包括兩種基本算子，分別是膨脹運算和腐蝕運算。設(shè)原始信號f(n)和結(jié)構(gòu)元素Se(n)分別為定義在集合F＝{0,1,…,N-1}和集合G＝{0,1,…,M-1}上的兩個一維離散函數(shù)，且N≥M。在每一個分析尺度λ下，令Se(n)對f(n)進行一次膨脹和腐蝕運算，即：

式中：表示膨脹運算，Θ表示腐蝕運算，λ＝1,2,…,λ_max，λ_max為最大分析尺度。

定義在尺度λ下f(n)關(guān)于Se(n)膨脹和腐蝕運算的覆蓋面積A_Se(λ)為

Maragos證明，A_Se(λ)滿足式(4)

式中D_M為信號的Minkowski-Bouligand維數(shù)即數(shù)學(xué)形態(tài)學(xué)分形維數(shù)，c為常數(shù)，對和log(1/λ)進行最小二乘線性擬合即可得到對信號數(shù)學(xué)形態(tài)學(xué)分形維數(shù)的估計。

2支持向量回歸機

SVR算法一般用于解決數(shù)據(jù)回歸預(yù)測問題，通過建立一個最優(yōu)超平面，使?fàn)顟B(tài)空間內(nèi)的各個數(shù)據(jù)點距離該超平面最近，并將該超平面作為回歸模型進行預(yù)測。

訓(xùn)練樣本對(x₁,y₁)，(x₂,y₂)，…，(x_r,y_r)，x_i∈Rⁿ，y_i∈R，x_i為輸入樣本，y_i為輸出樣本，i＝1,2,...,r。SVR的回歸函數(shù)表達式為

f(x)＝<ω·x_i>+b (5)

式中<ω·x_i>為ω和x_i的內(nèi)積。系數(shù)ω和b通過求解式(6)的最小值優(yōu)化問題獲得

由于擬合曲線必然會存在誤差，但要把誤差控制在一定的可容許范圍內(nèi)，ε為不敏感誤差，式(6)中L_ε為不敏感損失函數(shù)，其表達式為

對于回歸錯誤的數(shù)據(jù)點，需要引入松弛變量ξ和ξ^*。將不敏感損失函數(shù)Lε代入式(6)可得

滿足約束條件

求解上式時，一般采用對偶理論將其轉(zhuǎn)化為二次規(guī)劃問題。對于非線性數(shù)據(jù)，引入非線性映射函數(shù)Φ，建立拉格朗日方程，通過化簡，可得式(8)的對偶式

式中α，α^*為拉格朗日乘子。

令K(x_i,x_j)＝Φ(x_i)^TΦ(x_j)為特征空間的內(nèi)積，稱K(x_i,x_j)為核函數(shù)，選擇徑向基函數(shù)作為核函數(shù)。根據(jù)KKT定理，可求得變量α_i，α_i^*，b，最終可得支持向量機的回歸函數(shù)為

3基于IFOA的SVR尋優(yōu)

3.1 IFOA算法

FOA是一種基于果蠅覓食行為推演出尋求全局優(yōu)化的算法。利用在三維空間搜索的IFOA對SVR模型進行優(yōu)化，可以有效地獲取模型中的最優(yōu)參數(shù)。

如圖1所示為果蠅三維空間覓食示意圖。

IFOA步驟：

(1)初始化算法參數(shù)。設(shè)置果蠅種群規(guī)模Sizepop，最大覓食代數(shù)Maxgen，并隨機初始化果蠅群體位置坐標(biāo)(X₀,Y₀,Z₀)。

(2)果蠅個體利用嗅覺隨機搜索的方向和距離可以通過式(12)獲得

式中，i＝1,2,...,Sizepop，L₀為初始步長值，gen為當(dāng)前覓食代數(shù)。

(3)由于無法確定食物源的具體位置，所以需要通過式(13)估計第i個果蠅個體的當(dāng)前位置與坐標(biāo)原點間的距離Dist_i，之后計算出味道濃度判定值S_i

S_i＝1/Dist_i (14)

(4)將S_i代入味道濃度判定函數(shù)，計算出果蠅個體當(dāng)前位置的味道濃度

Smell_i＝function(S_i) (15)

(5)當(dāng)前果蠅群體中具有最高味道濃度的個體，可以由式(16)獲得

[bestSmell,bestIndex]＝max(Smell_i) (16)

(6)保留果蠅群體中最佳味道濃度值和與其對應(yīng)的個體坐標(biāo)，同時果蠅群體利用自身的視覺對食物源進行定位，然后飛往食物源所在的位置。

(7)進入迭代尋優(yōu)過程，重復(fù)步驟(2)-(5)，并判斷當(dāng)前味最高道濃度是否好于前一迭代味道濃度，且gen<Maxgen；若成立，則執(zhí)行步驟(6)。

3.2 IFOA-SVR模型

(1)初始化IFOA參數(shù)。包括Sizepop、Maxgen以及果蠅個體初始位置。由于本發(fā)明需要對SVR中的3個參數(shù)C、g以及ε進行尋優(yōu)，所以初始坐標(biāo)為(X₀¹,Y₀¹,Z₀¹)，(X₀²,Y₀²,Z₀²)和(X₀³,Y₀³,Z₀³)。

(2)附與每個果蠅個體隨機方向和飛行距離，并用遞減搜索步長代替固定步長，得到(X_i¹,Y_i¹,Z_i¹)，(X_i²,Y_i²,Z_i²)，(X_i³,Y_i³,Z_i³)。計算當(dāng)前果蠅個體與原點間距離的倒數(shù)，得到味道濃度判定值S_i¹，S_i²和S_i³。

(3)確定SVR中參數(shù)C、g和ε的范圍，即C∈[2^-14,2¹⁴]，g∈[2^-14,2¹⁴]，ε∈[2^-14,2¹⁴]。

(4)將訓(xùn)練樣本對輸入到SVR模型中，對模型進行訓(xùn)練，將MAE、RMSE、NMSE以及MAPE的和作為適應(yīng)度函數(shù)，即

(5)找到Fitness的最小值對應(yīng)的果蠅個體，開始迭代尋優(yōu)，并判斷最小Fitness是否低于前一代最小Fitness。如低于，則保留最小Fitness值及其對應(yīng)的坐標(biāo)，并將其賦給初始坐標(biāo)。如高于，則返回步驟(2)。

(6)找到C、g和ε的最佳值，建立IFOA-SVR預(yù)測模型。

4滾動軸承運行狀態(tài)可靠度

WPHM建立了設(shè)備運行狀態(tài)特征指標(biāo)與可靠度之間的數(shù)學(xué)關(guān)系，表達式為

式中β為形狀參數(shù)，η為尺度參數(shù)，μ為協(xié)變量回歸參數(shù)，Z為振動特征指標(biāo)數(shù)據(jù)。

h(t,Z)和可靠度函數(shù)R(t,Z)之間的關(guān)系為

從而，R(t,Z)可表示為

利用極大似然估計法對待定參數(shù)β、η、μ進行估計，確定待定參數(shù)最終取值。

5滾動軸承可靠度預(yù)測方法及流程

在WPHM參數(shù)確定的情況下，建立可靠度模型，將特征預(yù)測結(jié)果嵌入已建立的可靠度模型中，即可預(yù)測未來一段時間內(nèi)軸承運行的可靠度，具體的預(yù)測流程見圖2。

(1)通過Hilbert變換計算出狀態(tài)監(jiān)測數(shù)據(jù)的包絡(luò)信號，利用式(4)，求出包絡(luò)信號的數(shù)學(xué)形態(tài)學(xué)分形維數(shù)D_M。

(2)選取適當(dāng)?shù)挠?xùn)練樣本對，并利用IFOA對SVR模型中的參數(shù)C、g以及ε同時進行尋優(yōu)，建立IFOA-SVR預(yù)測模型。

(3)選取一定長度的振動特征指標(biāo)數(shù)據(jù)作為WPHM的響應(yīng)協(xié)變量，結(jié)合極大似然估計獲得似然函數(shù)方程組，把方程組中每個方程絕對值的和作為IFOA的適應(yīng)度函數(shù)，求出方程組的解，確定WPHM的待定參數(shù)，進而得到可靠度模型。

(4)將振動特征指標(biāo)數(shù)據(jù)作為IFOA-SVR預(yù)測模型的輸入，采用長期迭代預(yù)測法^[16]獲取振動特征預(yù)測結(jié)果。

(5)將IFOA-SVR預(yù)測模型得出的預(yù)測結(jié)果作為響應(yīng)協(xié)變量嵌入到可靠度模型中，即可計算出這些特征預(yù)測結(jié)果所對應(yīng)的可靠度，實現(xiàn)滾動軸承可靠度的預(yù)測。

針對本明提出的方法進行具體應(yīng)用與分析

滾動軸承全壽命數(shù)據(jù)來自于Cincinnati大學(xué)^[17]，每隔10min進行一次數(shù)據(jù)采集，共采集984次，每次采樣長度為20480點，采樣頻率為20KHz。圖3為第860次采集數(shù)據(jù)的時域波形，圖4為其對應(yīng)的包絡(luò)信號波形。

在每次采集的數(shù)據(jù)中，取中間的4096點進行Hilbert變換，求出包絡(luò)信號，進而計算出數(shù)學(xué)形態(tài)學(xué)分形維數(shù)，并作為一個特征數(shù)據(jù)點。通過控制變量法實驗對比證明，當(dāng)λ_max＝11，Se(n)＝(0,0,0,0,0,0,0,0,0,0,0,0)時，利用式(4)求出包絡(luò)信號的數(shù)學(xué)形態(tài)學(xué)分形維數(shù)能夠作為刻畫滾動軸承性能退化過程的特征。

將滾動軸承性能退化狀態(tài)特征作為響應(yīng)協(xié)變量，取600點到860點之間的特征指標(biāo)值，采用極大似然估計，利用IFOA求得WPHM模型的待定參數(shù)，如表1所示。

表1 WPHM待定參數(shù)的估計結(jié)果

將待定參數(shù)代入式(19)，即可確定當(dāng)數(shù)學(xué)形態(tài)學(xué)分形維數(shù)作為響應(yīng)協(xié)變量時的WPHM，進而建立可靠度模型，按照式(22)可以計算出滾動軸承在任意時刻下的可靠度

式中m＝1,2,…,984，984表示數(shù)據(jù)采集總次數(shù)；

每個數(shù)據(jù)點等于10分鐘，數(shù)據(jù)點與時間存在對應(yīng)關(guān)系。

將第1點至第860點之間對應(yīng)的特征指標(biāo)值嵌入可靠度模型即可計算出在這段時間內(nèi)軸承運行的可靠度，如圖5所示。

由圖5可知，滾動軸承在525點左右出現(xiàn)早期故障，在525點以前軸承對應(yīng)的可靠度為1。之后可靠度開始逐漸減小，到860點時滾動軸承可靠度已經(jīng)降到0.25附近，此時軸承處于故障程度比較嚴(yán)重的狀態(tài)。

利用IFOA對SVR中的參數(shù)C，g以及ε進行尋優(yōu)，L₀＝30，Sizepop＝20，Maxgen＝100。訓(xùn)練數(shù)據(jù)取第680點到第860點，輸入向量維數(shù)為20，建立IFOA-SVR預(yù)測模型。表2為利用不同的優(yōu)化算法對SVR進行尋優(yōu)得到的參數(shù)值。

表2利用不同算法尋優(yōu)得到的3個參數(shù)值

預(yù)測數(shù)據(jù)從第861點開始，共進行了50次迭代運算，通過長期迭代預(yù)測法得到最終預(yù)測結(jié)果。如果只優(yōu)化SVR模型中的C和g，預(yù)測誤差會很大并且預(yù)測步長較短。

利用GA對SVR模型中的C、g以及ε同時進行尋優(yōu)，預(yù)測結(jié)果如圖6所示。

利用PSO對SVR模型中的C、g以及ε同時進行尋優(yōu)，預(yù)測結(jié)果如圖7所示。

利用DS-FOA對SVR模型中的C、g以及ε同時進行尋優(yōu)，預(yù)測結(jié)果如圖8所示。

利用IFOA對SVR模型中的C、g以及ε同時進行尋優(yōu)，并將MAE、RMSE、NMSE以及MAPE的和作為適應(yīng)度函數(shù)，預(yù)測結(jié)果如圖9所示。

實際退化曲線在890點以后呈上升趨勢，由圖6可以看出，預(yù)測曲線在890點以后，明顯偏離實際退化曲線，預(yù)測效果不好；由圖7可以看出，預(yù)測曲線在890點以后上升趨勢不明顯；由圖8可以看出，預(yù)測曲線與實際退化曲線較為接近，但是在893至901點，與實際退化曲線相差較大；對比圖6至圖9可以看出，IFOA優(yōu)化SVR所得出的預(yù)測曲線與實際性能退化曲線更為接近，預(yù)測效果更好。

為了定量地評估預(yù)測結(jié)果，分別計算所預(yù)測50個點的MAE、RMSE、NMSE和MAPE，計算結(jié)果如表3所示。

表3 50點預(yù)測誤差分析

由表3可知，在本發(fā)明對SVR中3個參數(shù)均優(yōu)化的情況下，進行50點預(yù)測時，與GA、PSO以及DS-FOA優(yōu)化的SVR預(yù)測模型相比，IFOA-SVR模型的4類預(yù)測誤差均最小，其4類誤差和也最小。證明IFOA具有更好的空間尋優(yōu)能力，能夠準(zhǔn)確地找到SVR模型中參數(shù)的最佳值，體現(xiàn)了IFOA-SVR模型在滾動軸承性能退化趨勢預(yù)測方面的優(yōu)勢。

與文獻[7]到[9]對比分析可知，本發(fā)明所得實驗結(jié)果MAPE低于文獻[7]，預(yù)測步長增加了30點；與文獻[8]相比，MAPE在同一數(shù)量級內(nèi)稍高，但預(yù)測步長增加了43點；相比于文獻[9]，在降低了MAE和MAPE的同時預(yù)測步長增加了20點。

將預(yù)測所得50點數(shù)據(jù)嵌入可靠度模型中，可以計算出這50點的可靠度值，可靠度曲線如圖10所示。

從圖10中可以看出，第861點的可靠度值在0.25附近，并且這50點的可靠度呈下降趨勢，到910點時可靠度已經(jīng)下降到0.08左右。這表明軸承運行狀態(tài)很差，已經(jīng)出現(xiàn)非常嚴(yán)重的故障，應(yīng)該做好及時更換軸承的準(zhǔn)備。

通過對本發(fā)明方法的應(yīng)用得出如下結(jié)論：

(1)提出Hilbert包絡(luò)解調(diào)與數(shù)學(xué)形態(tài)學(xué)分形維數(shù)相結(jié)合的特征提取方法。通過控制變量法得出最大分析尺度λ_max以及結(jié)構(gòu)元素Se這兩個參數(shù)的最佳設(shè)置，計算出包絡(luò)信號的出數(shù)學(xué)形態(tài)學(xué)分形維數(shù)并作為刻畫滾動軸承性能退化過程的特征。

(2)將MAE、RMSE、NMSE以及MAPE的和作為IFOA的適應(yīng)度函數(shù)，對SVR模型中的3個參數(shù)進行尋優(yōu)，建立了IFOA-SVR模型。利用該模型對滾動軸承的性能退化趨勢進行了50點預(yù)測，與GA，PSO和DS-FOA優(yōu)化的SVR模型相比預(yù)測誤差最小。同時，與相關(guān)文獻進行了橫向?qū)Ρ?，在保證預(yù)測精度的前提下，增加了預(yù)測步長。

(3)利用IFOA求解WPHM的待定參數(shù)，建立了可靠度模型。將IFOA-SVR模型所得的特征預(yù)測結(jié)果嵌入到可靠度模型中，得到可靠度預(yù)測曲線，能夠預(yù)測滾動軸承在未來50點的運行狀態(tài)可靠度。

本發(fā)明中涉及的參考文獻：

[1]何正嘉,曹宏瑞,訾艷陽,等.機械設(shè)備運行可靠性評估的發(fā)展與思考[J].機械工程學(xué)報,2014,50(2)：171-186.

HE Zhengjia,CAO Hongrui,ZI Yanyang,et al.Developments and Thoughts on Operational Reliability Assessment of Mechanical Equipment[J].Journal of Mechanical Engineering,

2014,50(2)：171-186.

[2]WANG Zhonglai,HUANG Hongzhong,LI Yanfeng.An approach to reliability assessment under degradation and shock process[J].IEEE Transactions on Reliability,2011,60(4)：852-863.

[3]何群,李磊,江國乾,等.基于PCA和多變量極限學(xué)習(xí)機的軸承剩余壽命預(yù)測[J].中國機械工程,2014,25(7)：984-988+989.

HE Qun,LI Lei,JIANG Guoqian,et al.Residual Life Predictions for Bearings Based on PCA and MELM[J].China Mechanical Engineering,2014,25(7)：984-988+989.

[4]申中杰,陳雪峰,何正嘉,等.基于相對特征和多變量支持向量機的滾動軸承剩余壽命預(yù)測[J].機械工程學(xué)報,2013,49(2)：183-189.

SHEN Zhongjie,CHEN Xuefeng,HE Zhengjia,et al.Remaining Life Predictions of Rolling Bearing Based on Relative Features and Multivariable Support Vector Machine[J].Journal of Mechanical Engineering,2013,49(2)：183-189.

[5]陳向民,于德介,李蓉.基于形態(tài)分量分析和包絡(luò)譜的軸承故障診斷[J].中國機械工程,2014,25(8)：1047-1052+1053.

CHEN Xiangmin,YU Dejie,LI Rong.Fault Diagnosis of Rolling Bearings Based on Morphological Componenet Analysis and Envelope Spectrum[J].China Mechanical Engineering,2014,25(8)：1047-1052+1053.

[6]張亢,程軍圣,楊宇.基于局部均值分解與形態(tài)學(xué)分形維數(shù)的滾動軸承故障診斷方法[J].振動與沖擊,2013,32(9)：90-94.

ZHANG Kang,CHENG Junsheng,YANG Yu.Roller bearing fault diagnosis based on local mean decomposition and morphological fractal dimension[J].Journal of Vibration and Shock,2013,32(9)：90-94.

[7]陳昌,湯寶平,呂中亮.基于威布爾分布及最小二乘支持向量機的滾動軸承退化趨勢預(yù)測[J].振動與沖擊,2014,33(20)：52-56.

CHEN Chang,TANG Baoping,LU Zhongliang.Degradation trend prediction of rolling bearings based on Weibull distribution and least squares support vector machine[J].Journal of Vibration and Shock,2014,33(20)：52-56.

[8]宋梅村,蔡琦.基于支持向量回歸的設(shè)備故障趨勢預(yù)測[J].原子能科學(xué)技術(shù),2011,45(8)：972-976.

SONG Meicun,CAI Qi.Fault Trend Prediction of Device Based on Support Vector Regression[J].Atomic Energy Science and Technology,2011,45(8)：972-976.

[9]朱霄珣.基于支持向量機的旋轉(zhuǎn)機械故障診斷與預(yù)測方法研究[D].北京：華北電力大學(xué),2013.

ZHU Xiaoxun.Research on Rotating Machine Fault Diagnosis and State Prediction Method Based on Support Vector Machine[D].Beijing：North China Electric Power University,2013.

[10]P.J.García Nieto,E.García-Gonzalo,F.Sánchez Lasheras et al.Hybrid PSO-SVM-based method for forecasting of the remaining useful life for aircraft engines and evaluation of its reliability[J].Reliability Engineering and System Safety,2015,138：219-231.

[11]PAN W T.A new fruit fly optimization algorithm:Taking the financial distress model as an example[J].Knowledge-Based Systems,2012,26(2)：69-74.

[12]寧劍平,王冰,李洪儒,等.遞減步長果蠅優(yōu)化算法及應(yīng)用[J].深圳大學(xué)學(xué)報理工版,2014,31(4)：367-373.

NING Jianping,WANG Bing,LI Hongru,et al.Research on and application of diminishing step fruit fly optimization algorithm[J].Journal of Shenzhen University Science and Engineering,2014,31(4)：367-373.

[13]張義民,孫志禮.機械產(chǎn)品的可靠性大綱[J].機械工程學(xué)報,2014,50(14)：14-20.

ZHANG Yimin,SUN Zhili.The Reliability Syllabus of Mechanical Products[J].Journal of Mechanical Engineering,2014,50(14)：14-20.

[14]丁鋒,何正嘉,訾艷陽,等.基于設(shè)備狀態(tài)振動特征的比例故障率模型可靠性評估[J].機械工程學(xué)報,2009,45(12)：89-94.

DING Feng,HE Zhengjia,ZI Yanyang,et al.Reliability Assessment Based on Equipment Condition Vibration Feature Using Proportional Hazards Model[J].Journal of Mechanical Engineering,2009,45(12)：89-94.

[15]陳昌.基于狀態(tài)振動特征的空間滾動軸承可靠性評估方法研究[D].重慶：重慶大學(xué),2014.

CHEN Chang.Reliability Assessment Method for Space Rolling Bearing Based on Condition Vibration Feature[D].Chongqing：Chongqing University,2014.

[16]董紹江.基于優(yōu)化支持向量機的空間滾動軸承壽命預(yù)測方法研究[D].重慶：重慶大學(xué),2012.

DONG Shaojiang.Research on Space Bearing Life Prediction Method Based on Optimized Support Vector Machine[D].Chongqing：Chongqing University,2012.

[17]Rexnord technical services,bearing data set,IMS,University of Cincinnati,NASA Ames Prognostics Data Repository,NASA Ames,Moffett Field,CA,2007[2014-04].http://ti.arc.nasa.gov/project/prognostics-data-repository.