胡啟國(guó)， 何奇， 曹歷杰

(1. 重慶交通大學(xué) 機(jī)電與車輛工程學(xué)院， 重慶 400074；2. 川慶鉆探工程公司 安全環(huán)保質(zhì)量監(jiān)督檢測(cè)研究院， 四川 廣漢 618300)

滾動(dòng)軸承作為機(jī)械設(shè)備的典型零部件，是決定機(jī)械標(biāo)準(zhǔn)件故障預(yù)測(cè)和健康管理的關(guān)鍵部件之一，剩余壽命對(duì)于衡量軸承故障信息有著重要價(jià)值[1].隨著機(jī)器學(xué)習(xí)的發(fā)展，基于數(shù)據(jù)驅(qū)動(dòng)的軸承剩余壽命預(yù)測(cè)逐漸成為國(guó)內(nèi)外研究的主流.但軸承原始信號(hào)中不可避免地會(huì)有大量噪聲，影響了軸承壽命預(yù)測(cè)的精確度，因此，需要對(duì)原始信號(hào)進(jìn)行降噪處理.為了進(jìn)行退化特征的提取和壽命預(yù)測(cè)，Yan等[2]提出一種使用支持向量機(jī)和混合退化跟蹤模型，預(yù)測(cè)軸承的剩余使用壽命；Liu等[3]提出一種磷蝦群算法，優(yōu)化基于混合核函數(shù)-支持向量回歸(HKF-SVR)對(duì)軸承性能退化趨勢(shì)的預(yù)測(cè)；葛陽等[4]基于t-分布領(lǐng)域嵌入算法(t-SNE)和長(zhǎng)短期記憶(LSTM)神經(jīng)網(wǎng)絡(luò)，對(duì)滾動(dòng)軸承的剩余壽命進(jìn)行研究；Jiang等[5]提出一種基于時(shí)間序列的多通道卷積神經(jīng)網(wǎng)絡(luò)的軸承剩余使用壽命模型.上述方法為軸承退化特征的提取和壽命預(yù)測(cè)提供了有效借鑒，但仍有考慮不周的地方，如對(duì)原始信號(hào)噪聲沒有進(jìn)行有效處理，壽命預(yù)測(cè)只能預(yù)測(cè)短期的退化趨勢(shì)且預(yù)測(cè)結(jié)果較為保守.Lu等[6]提出最小二乘支持向量機(jī)(LSSVM)模型，該模型具有在預(yù)測(cè)壽命過程中預(yù)測(cè)精度高和在非線性樣本數(shù)據(jù)中泛化能力好的優(yōu)點(diǎn)，其預(yù)測(cè)精度主要受正規(guī)化參數(shù)C和核參數(shù)δ取值不同的影響.

針對(duì)以上問題，本文以提高軸承壽命預(yù)測(cè)精度為目標(biāo)，提出一種新的壽命預(yù)測(cè)方法.首先，用集合經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD)方法對(duì)原始振動(dòng)信號(hào)進(jìn)行去噪，再根據(jù)相關(guān)性系數(shù)和峭度值對(duì)固有模態(tài)函數(shù)(IMF)分量信號(hào)進(jìn)行篩選并重構(gòu)數(shù)據(jù)；然后，提取重構(gòu)數(shù)據(jù)中的時(shí)域、頻域、小波包能量譜等22個(gè)特征指標(biāo)作為特征指標(biāo)集，并對(duì)其進(jìn)行核主成分分析，選取累計(jì)貢獻(xiàn)率達(dá)到要求的核主成分作為性能退化特征，建立改進(jìn)的哈里斯鷹優(yōu)化-最小二乘支持向量機(jī)(IHHO-LSSVM)模型.

1 集合經(jīng)驗(yàn)?zāi)B(tài)分解原理

Wu等[7]針對(duì)經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)在使用過程中存在噪聲或脈沖產(chǎn)生模態(tài)混疊的現(xiàn)象，提出一種EEMD算法，有效地抑制了EMD產(chǎn)生的模態(tài)混疊現(xiàn)象.EEMD具體分解過程如下.

1) 在信號(hào)中加入一定幅值的白噪聲序列，得到新的信號(hào)x(t).

2) 新信號(hào)x(t)經(jīng)過EEMD分解，獲得h個(gè)IMF分量cj(t)和1個(gè)余量r(t)，即

(1)

3) 將分解得到的各個(gè)IMF分量和剩余余量的均值作為最后結(jié)果，最大限度地消除上、下包絡(luò)線擬合誤差.每個(gè)IMF分量具有穩(wěn)態(tài)性能和相互獨(dú)立的特點(diǎn).分解后的IMF表示為

(2)

2 軸承壽命預(yù)測(cè)模型的建立

2.1 最小二乘支持向量機(jī)預(yù)測(cè)模型

最小二乘支持向量機(jī)由改進(jìn)支持向量機(jī)(SVM)演化而來，其原理是用等式約束替換SVM中的不等式約束，將不等式約束的二次規(guī)劃問題轉(zhuǎn)化成線性矩陣求解問題，與標(biāo)準(zhǔn)的SVM相比，LSSVM的求解速度更快，精度更高.LSSVM的函數(shù)估計(jì)問題為

(3)

s.t.yi-ωΦ(xj)-b=0.

(4)

式(3)，(4)中：ω為權(quán)向量；e為單位矩陣；ei為松弛變量，表示第i個(gè)數(shù)據(jù)的預(yù)測(cè)輸出和實(shí)際輸出的誤差值；γ為懲罰參數(shù)；Φ=[φ(x1)+φ(x2)，…，φ(xn)]T；b為偏差向量.

式(3)和式(4)對(duì)應(yīng)的拉格朗日函數(shù)為

(5)

式(5)中：α=[α1，α2，…，αn]T；αi為L(zhǎng)agrange乘子.依據(jù)卡魯什-庫恩-塔克(KKT)條件進(jìn)行求解，對(duì)式(5)中的ω，b，e，αi求偏導(dǎo)數(shù)，并令其偏導(dǎo)數(shù)為零.消掉變量ω，e，最終得到的方程組為

(6)

式(6)中:I=[1，1，…，1]T；y=[y1，y2，…，yn]T.

根據(jù)Mercer條件，選擇徑向基核函數(shù)(RBF)作為核函數(shù)，其公式為

(7)

式(7)中:σ為核函數(shù)寬密度.

最終得到的LSSVM預(yù)測(cè)回歸模型為

(8)

2.2 改進(jìn)的哈里斯鷹優(yōu)化(IHHO)算法

哈里斯鷹優(yōu)化(HHO)算法[8]是Heidari受栗翅鷹捕食獵物過程中的群體合作和追逐方式啟發(fā)，提出的一種群智能優(yōu)化算法，可根據(jù)場(chǎng)景的動(dòng)態(tài)性和獵物的逃跑模式揭示多種追逐模式，由搜索階段、轉(zhuǎn)化階段、開發(fā)階段組成.但原算法易陷入局部最優(yōu)值，因此，對(duì)搜索階段和能量調(diào)控機(jī)制進(jìn)行改進(jìn).

改進(jìn)的搜索階段中，哈里斯鷹棲息在[lb，ub]的空間某位置觀察獵物，迭代時(shí)根據(jù)所在位置更新概率p.當(dāng)p<0.5時(shí)，其他個(gè)體捕食獵物被認(rèn)為是模仿哈里斯鷹曲折運(yùn)動(dòng)的，采用Levy函數(shù)(LF)飛行模擬；當(dāng)p≥0.5時(shí)，哈里斯鷹個(gè)體移動(dòng)在內(nèi)部是隨機(jī)的，因此，通過哈里斯鷹與其他3只個(gè)體初始值的平均位置來更新位置[9]，即

(9)

式(9)中：Xj(t)，Xk(t)，Xl(t)為哈里斯鷹隨機(jī)初始選擇的位置；X0(t)為獵物位置；X(t+1)，X(t)分別為第t+1，t次迭代時(shí)個(gè)體的位置；Xm(t)為第t次迭代時(shí)N只個(gè)體的平均位置；D為問題維度.

HHO算法根據(jù)獵物的逃逸能量E在搜索和不同的開發(fā)行為之間轉(zhuǎn)換，但E的變換是單周期性遞減的，無法描述多個(gè)輪次共同追捕獵物的情況[10].因此，提出一種調(diào)控能量周期性的機(jī)制，即多個(gè)輪次數(shù)以獵物逃逸能量E的遞變周期數(shù)表示，實(shí)現(xiàn)“全局+局部”的尋優(yōu)搜索能力.以余弦函數(shù)刻畫E的周期遞變性，其定義式為

(10)

式(10)中：E0為獵物的初始能量，E0∈[-1，1]；T為最大迭代數(shù)；kp為獵物能量的周期遞減數(shù)，kp=0，1，2，….當(dāng)|E|≥1時(shí)，HHO進(jìn)入全局搜索階段；當(dāng)|E|<1時(shí)，HHO算法進(jìn)入局部開發(fā)過程.

在HHO算法中，通過因子β∈[0，1]描述獵物是否成功逃脫，當(dāng)β<0.5時(shí)，表示獵物逃脫成功，反之則失效.在開發(fā)階段的尋優(yōu)過程中包含了4種情形[11].

1) 軟圍攻.當(dāng)β≥0.5且0.5≤|E|≤1.0時(shí)，獵物的能量E充足并逃脫，但最終被捕獲，位置更新為

X(t+1)=ΔX0(t)-E|JpX0(t)-X0(t)|.

(11)

式(11)中：ΔX0(t)為第t次迭代時(shí)鷹與獵物的位置偏差；Jp為獵物逃跑時(shí)的隨機(jī)跳躍能力，Jp=2(1-r5)，r5是(0，1)上的一個(gè)隨機(jī)數(shù).

2) 硬圍攻.當(dāng)β≥0.5且|E|<0.5時(shí)，獵物因能量E較低而被鷹直接捕獲，位置更新為

X(t+1)=X0(t)-|EΔX0(t)|.

(12)

3) 漸進(jìn)式快速俯沖圍捕.當(dāng)β<0.5且0.5<|E|<1.0時(shí)，鷹以最優(yōu)方向俯沖軟圍捕獵物，位置為

(13)

式(13)中：S是D維隨機(jī)向量.

4) 漸進(jìn)式快速俯沖硬捕.當(dāng)β<0.5且|E|<0.5時(shí)，獵物能量較低，鷹俯沖硬捕獵物，位置更新為

(14)

2.3 IHHO-LSSVM模型流程圖

選取軸承退化特征作為輸入向量，以預(yù)測(cè)壽命均方誤差(MSE)最小化為IHHO尋優(yōu)目標(biāo).以核參數(shù)δ和正規(guī)化參數(shù)C為優(yōu)化參數(shù)的IHHO-LSSVM模型流程圖，如圖1所示.

圖1 IHHO-LSSVM模型流程圖Fig.1 Flow chart of IHHO-LSSVM model

圖2 軸承實(shí)驗(yàn)臺(tái)Fig.2 Bearing test bench

3 基于軸承全壽命預(yù)測(cè)的驗(yàn)證

3.1 實(shí)驗(yàn)裝置

采用Gousseau等[12]得到的滾動(dòng)軸承全壽命周期加速軸承性能退化數(shù)據(jù).軸承實(shí)驗(yàn)臺(tái)，如圖2所示.實(shí)驗(yàn)過程中，每個(gè)軸承承受26 670 N的徑向載荷，每個(gè)軸承的橫向和縱向各安裝了一個(gè)壓電式加速度傳感器，采樣間隔為10 min，共采集984個(gè)樣本，采樣頻率為20 kHz，共采集20 152 320個(gè)數(shù)據(jù).軸承運(yùn)行了一周，直到軸承退化失效.原始信號(hào)圖，如圖3所示.圖3中：a為加速度；N為數(shù)據(jù)樣本數(shù)量；f為頻率.

(a) 原始時(shí)域圖 (b) 原始頻譜圖圖3 原始信號(hào)圖Fig.3 Original signal diagram

3.2 基于EEMD的軸承信號(hào)降噪

采集的原始信號(hào)中一般含有噪聲，因此，需要對(duì)原始數(shù)據(jù)進(jìn)行降噪處理.分解得到21個(gè)IMF分量，其相關(guān)系數(shù)和峭度值，如表1所示.各個(gè)IMF分量，如圖4所示.

表1 各個(gè)IMF的相關(guān)系數(shù)和峭度值Tab.1 Correlation coefficient and kurtosis value of each IMF

由圖4可知：IMF1～I(xiàn)MF12的振動(dòng)沖擊信號(hào)較多，而IMF13～I(xiàn)MF20缺少故障的退化信息，因此，對(duì)前12個(gè)IMF分量進(jìn)行篩選.分解結(jié)果中，IMF分量容易產(chǎn)生虛假的分量，故引入相關(guān)系數(shù)和峭度值以辨別和剔除虛假分量.根據(jù)相關(guān)系數(shù)和峭度值的計(jì)算結(jié)果可知：相關(guān)系數(shù)越接近1，表明該IMF與原始信息關(guān)聯(lián)越緊密，對(duì)應(yīng)的距離越??；峭度值偏離越大，說明包含的軸承故障信息越多[13].因此，根據(jù)表1選取IMF1，IMF3～I(xiàn)MF9作為重構(gòu)信號(hào)分量.

(a) IMF1 (b) IMF2 (c) IMF3 (d) IMF4

3.3 驗(yàn)證結(jié)果與分析

由于時(shí)域、頻域各個(gè)特征的敏感度皆不相同，篩選對(duì)故障敏感度較高的特征，剔除退化趨勢(shì)不明顯的特征.在時(shí)域退化特征中，剔除脈沖指標(biāo)、裕度指標(biāo)、歪度等幾乎無任何變化的退化特征，篩選出8個(gè)時(shí)域特征，如圖5所示.在13個(gè)單一頻域特征中，選取能夠很好地呈現(xiàn)軸承全壽命退化過程的6個(gè)頻域特征指標(biāo)，如圖6所示.圖6中：P1為頻域幅值均值；P2為頻域標(biāo)準(zhǔn)差；P4為頻域主頻帶位置；P5，P7，P13為頻譜分散或集中程度.

(a) 均方根 (b) 標(biāo)準(zhǔn)差 (c) 最大值

由圖5可知：8個(gè)時(shí)域特征均能反映軸承性能的退化趨勢(shì)，但其敏感程度各有不同；均方根、標(biāo)準(zhǔn)差和絕對(duì)平均值在第500個(gè)數(shù)據(jù)點(diǎn)后呈現(xiàn)良好的上升趨勢(shì)；而峰值、方差、最大值則在第600個(gè)數(shù)據(jù)點(diǎn)后逐漸上升.因此，單一的時(shí)域特征可能存在評(píng)估能力不足的問題.

(a) P1 (b) P2 (c) P4

由圖6可知：P1在第500個(gè)數(shù)據(jù)點(diǎn)之后出現(xiàn)了小幅度的變化，隨著故障程度的加深，其特征指標(biāo)值也不斷變大；P2，P4，P5，P7，P13也都反映了軸承退化的趨勢(shì).

小波包分解克服了小波分解的不足，提高了信號(hào)高頻帶的頻率分辨率和低頻帶的時(shí)間分辨率[14]，其原理是將原始信號(hào)分解到不同頻帶，生成頻帶能量譜.選用db3小波基函數(shù)對(duì)重構(gòu)信號(hào)進(jìn)行3層的小波包分解，得到8個(gè)能夠反映軸承全壽命的小波包能量譜(S3，1～S3，8)，如圖7所示.S3，i中的3表示分解層數(shù)，i表示8個(gè)頻段位置，i=1，2，…，8.

(a) S3，1 (b) S3，2 (c) S3，3 (d) S3，4

由圖7可知：8個(gè)小波包能量譜能有效增強(qiáng)故障和退化特征，完整保留信號(hào)頻率信息.在非線性信號(hào)中，相比時(shí)域、頻域的單一特征，小波包能量譜具有分辨率高、精度細(xì)化高的優(yōu)點(diǎn).小波包能量譜與時(shí)域、頻域共同構(gòu)成多維特征集，更能反映軸承全壽命過程的運(yùn)行狀況.共選取22個(gè)退化特征，利用核主成分分析(KPCA)將特征退化指標(biāo)加權(quán)融合.將KPCA第一主成分標(biāo)準(zhǔn)化處理，結(jié)果如圖8所示.

圖8 KPCA第一主成分趨勢(shì)圖Fig.8 First principal component trend chart of KPCA

由圖8可知：KPCA第一主成分在第500個(gè)數(shù)據(jù)點(diǎn)(早期故障)后開始呈現(xiàn)優(yōu)良的上升趨勢(shì).單一的退化特征則不具備這一優(yōu)勢(shì)，利用KPCA加權(quán)融合提取主成分呈現(xiàn)出的軸承退化性能趨勢(shì)比單一的特征更全面，該方法使維數(shù)得到約簡(jiǎn)，計(jì)算復(fù)雜程度減小.將KPCA各主成分的貢獻(xiàn)率從大到小依次排列，各主成分的貢獻(xiàn)率，如表2所示.

表2 各主成分的貢獻(xiàn)率Tab.2 Contribution rate of each principal component

選取累計(jì)貢獻(xiàn)率達(dá)到85%的前幾個(gè)主成分[15].由表2可知：前3個(gè)主成分的貢獻(xiàn)率相加得到的累計(jì)貢獻(xiàn)率為86.64%，已超過85%.故選取軸承退化期的前3個(gè)主成分和對(duì)應(yīng)的剩余壽命序列組成訓(xùn)練樣本，建立IHHO-LSSVM壽命預(yù)測(cè)模型.將第500～703個(gè)數(shù)據(jù)點(diǎn)的退化特征和每一個(gè)樣本對(duì)應(yīng)的剩余壽命時(shí)間序列作為訓(xùn)練樣本，等間隔選取第704～984個(gè)數(shù)據(jù)點(diǎn)間的7個(gè)樣本作為預(yù)測(cè)樣本，代入剩余壽命預(yù)測(cè)模型[16]中，將LSSVM的核參數(shù)δ和正規(guī)化參數(shù)C通過改進(jìn)的哈里斯鷹算法進(jìn)行尋優(yōu).對(duì)比HHO算法與IHHO算法的適應(yīng)度(適應(yīng)度值選擇均方誤差)曲線，如圖9所示.圖9中：n為迭代次數(shù).由圖9可知：IHHO算法的收斂精度比HHO算法提高了許多.

圖9 HHO算法與IHHO算法的適應(yīng)度曲線對(duì)比Fig.9 Comparison of fitness curves between HHO algorithm and IHHO algorithm

文中方法是利用EEMD去除噪聲，分解、重構(gòu)信號(hào)，再用KPCA融合的前3個(gè)主成分作為軸承性能退化指標(biāo)，最后采用IHHO-LSSVM作為壽命預(yù)測(cè)模型進(jìn)行模擬.為了評(píng)估軸承壽命的預(yù)測(cè)效果，將IHHO-LSSVM模型與采用EEMD-KPCA處理的LSSVM模型和IHHO-LSSVM模型的壽命預(yù)測(cè)結(jié)果進(jìn)行對(duì)比，如圖10所示.圖10中：tr為剩余壽命.由圖10可知：3種模型預(yù)測(cè)的軸承剩余壽命與真實(shí)壽命變化趨勢(shì)一致.

(a) EEMD-KPCA處理的LSSVM模型

3種模型的預(yù)測(cè)誤差對(duì)比，如表3所示.對(duì)比圖10和表3可知：EEMD-KPCA處理的IHHO-LSSVM模型的均方根誤差、平均絕對(duì)誤差和均方誤差均小于IHHO-LSSVM模型和EEMD-KPCA處理的LSSVM模型，說明EEMD分解重構(gòu)具有消除噪聲的作用，KPCA消除了時(shí)域、頻域等各特征指標(biāo)間相關(guān)性導(dǎo)致的冗余數(shù)據(jù)；IHHO-LSSVM模型的預(yù)測(cè)準(zhǔn)確率相比其他支持向量機(jī)模型有所提高.

表3 3種模型的預(yù)測(cè)誤差對(duì)比Tab.3 Comparison of prediction errors of three models

4 結(jié)論

1) 提出一種新的EEMD-KPCA方法，對(duì)軸承進(jìn)行特征處理與提取，該方法既能去噪、又能采用核主成分分析方法剔除無關(guān)特征，提高了數(shù)據(jù)特征的敏感性，為建立壽命預(yù)測(cè)模型提供了良好的基礎(chǔ).

2) 提出一種改變搜索階段和能量調(diào)控機(jī)制的改進(jìn)的哈里斯鷹優(yōu)化算法.與原算法相比，該方法更容易找到LSSVM預(yù)測(cè)模型的核參數(shù)δ和正規(guī)化參數(shù)C的最優(yōu)解.

3) 建立一種采用EEMD-KPCA處理的IHHO-LSSVM滾動(dòng)軸承壽命預(yù)測(cè)模型，與IHHO-LSSVM預(yù)測(cè)模型和EEMD-KPCA處理的LSSVM模型相比，所提模型的預(yù)測(cè)精度得到了提高，為滾動(dòng)軸承壽命預(yù)測(cè)提供了一種新的方法.