VMD與多重分形奇異譜的往復(fù)壓縮機(jī)故障特征預(yù)測(cè)方法

劉 巖 ，王金東 ，李 穎

（1.東北石油大學(xué) 機(jī)械科學(xué)與工程學(xué)院，黑龍江 大慶 163318；2.齊齊哈爾大學(xué) 機(jī)電工程學(xué)院，黑龍江 齊齊哈爾 161006）

1 引言

對(duì)于天然氣與石化行業(yè)的大型往復(fù)式壓縮機(jī)故障診斷而言，軸承類故障因其監(jiān)測(cè)的困難性及后果的突發(fā)性與危害性，決定了對(duì)其開展有效的故障診斷的意義重大[1，6]。故障狀態(tài)及壽命預(yù)測(cè)作為故障診斷技術(shù)不可或缺的環(huán)節(jié)，在往復(fù)機(jī)軸承類故障診斷中其重要性尤為突出。

大型往復(fù)機(jī)作為典型的非線性系統(tǒng)，其預(yù)測(cè)不準(zhǔn)確的主要原因在于：系統(tǒng)動(dòng)力模型與實(shí)際的差異、初始條件的影響2方面[1]。從預(yù)測(cè)模型角度分析，確定性的參數(shù)類回歸模型由于模型影響參數(shù)過(guò)多、影響因素考慮不足以及模型修正過(guò)分依賴主觀性等不足，造成可靠度和精度不高而限制其應(yīng)用；非參數(shù)回歸方法是一種近年來(lái)在預(yù)測(cè)學(xué)中較熱門建模方法，特別適合于不確定性、非線性動(dòng)態(tài)系統(tǒng)。另外，對(duì)于初始條件的影響，帕科德（Packard）在1980年提出用時(shí)間序列相空間重構(gòu)理論，從靜態(tài)點(diǎn)與時(shí)間關(guān)系結(jié)合的角度再現(xiàn)了系統(tǒng)動(dòng)態(tài)演化軌跡；為復(fù)雜非線性系統(tǒng)演化趨勢(shì)分析提供有效途徑[1-2]。學(xué)者們基于KNN及與相空間重構(gòu)方法結(jié)合的改進(jìn)算法方面開展了大量研究，如文獻(xiàn)[3]在股票價(jià)格預(yù)測(cè)中提出一種基于徑向基網(wǎng)絡(luò)的KNN預(yù)測(cè)模型；文獻(xiàn)[4]提出EMD分解與KNN結(jié)合的股票價(jià)格改進(jìn)預(yù)測(cè)模型；文獻(xiàn)[5]將改進(jìn)KNN算法和相空間重構(gòu)結(jié)合在風(fēng)電功率預(yù)測(cè)中取得較好的效果。

然而，往復(fù)機(jī)振動(dòng)信號(hào)表現(xiàn)出的類周期震蕩和復(fù)雜的多重分形特征—即長(zhǎng)期不可預(yù)測(cè)性，使得對(duì)該類振動(dòng)信號(hào)預(yù)測(cè)而言，局限在振動(dòng)時(shí)序離散點(diǎn)值的預(yù)測(cè)，而忽略因初始條件敏感性決定的預(yù)測(cè)可信時(shí)間的方法顯然不可取，同時(shí)預(yù)測(cè)結(jié)果受熵值復(fù)雜度和類周期振動(dòng)的制約，有限個(gè)振動(dòng)幅值因無(wú)法展現(xiàn)系統(tǒng)整周期包含的細(xì)節(jié)特征而失去意義。

基于以上分析，筆者將最大預(yù)測(cè)可信時(shí)間引入預(yù)測(cè)算法，對(duì)時(shí)序信號(hào)等間隔周期采樣，利用VMD算法提取主分量并提取各采樣時(shí)段的多重分形奇異譜特征參數(shù)構(gòu)造預(yù)測(cè)序列，計(jì)算各模態(tài)分量的最大預(yù)測(cè)可信時(shí)間以確定預(yù)測(cè)可信時(shí)段，基于相空間重構(gòu)的改進(jìn)KNN創(chuàng)建初始建模域，以C_C算法變參數(shù)尋優(yōu)還原系統(tǒng)演化趨勢(shì)；算法將度量預(yù)測(cè)時(shí)間的熵值和基于模態(tài)分解的多重分形譜特征有效結(jié)合，將系統(tǒng)動(dòng)力學(xué)模型和初始條件影響因素作為整體統(tǒng)一在算法中。

2 系統(tǒng)預(yù)測(cè)可信時(shí)間

在香農(nóng)（Shannon）信息論中，對(duì)于M維時(shí)間序列的相空間點(diǎn)對(duì)，若將其劃分為k個(gè)尺寸為rM的單元，令單元序列si的生成概率為pi；其信息值Ik它正比于以精度r確定系統(tǒng)單元區(qū)域序列si（或稱軌道）所需信息；即Ik+1-Ik度量了系統(tǒng)從時(shí)間（k-1）τ到kτ的信息損失[1，6]。

系統(tǒng)的最大預(yù)測(cè)可信時(shí)間可用K（柯爾莫哥洛夫，Kolmogorov）熵表示，設(shè)t時(shí)刻信息量為I（t），經(jīng)過(guò)Δt后的信息量為I（t+Δt），則有：I（t+Δt）=I（t）-KΔt；

令I(lǐng)（t）=1，當(dāng)I（t+Δt）=0時(shí)，即可得到系統(tǒng)的最大預(yù)測(cè)可信時(shí)間T=1/K。由于Renyi熵具有性質(zhì)Kq≤Kq+1，可用Renyi熵的下界K2作為K熵的近似估計(jì)[1]。式（1）可近似用式（2）表示：

式中：Cm（r）—m維相空間距離小于r的點(diǎn)對(duì)占的比例，即C_C算

法的相似容限，故有：

式（4）中當(dāng)嵌入維數(shù)m飽和時(shí)，可得T≈1/K2。具體算法流程，如圖1所示。

圖1 預(yù)測(cè)可信時(shí)間計(jì)算流程Fig.1 The Flowchart of Maximum Dependable Prediction Time

3 變分模態(tài)分解

變分模態(tài)分解（VMD）是由Dragomiretskiy和Zosso于2014年提出的一種新的自適應(yīng)信號(hào)分解方法。通過(guò)求解變分框架中的約束變分模型，自適應(yīng)地將復(fù)雜的多分量振動(dòng)信號(hào)分解成一組有限帶寬的模態(tài)函數(shù)，實(shí)現(xiàn)各信號(hào)分量的頻率分離[7]。

對(duì)任意信號(hào)x（t），分解步驟可分為兩部分：

3.1 建立變分模型：

（1）定義有限帶寬的內(nèi)稟模態(tài)函數(shù)（Band-Limited Intrinsic Mode Function，BLIMF）；

（2）構(gòu)建解析信號(hào)獲取單邊譜，然后混頻得：

（3）估計(jì)帶寬：通過(guò)高斯平滑（即L2范數(shù)梯度的平方根）對(duì)信號(hào)解調(diào)獲得各模態(tài)函數(shù)帶寬；

（4）建立優(yōu)化模型：引入約束條件，構(gòu)造優(yōu)化變分模型如下：

式中：K—BLIMF 分量｛uk｝=｛u1，u2，L，uK｝的數(shù)目；

｛ωk｝=｛ω1，ω2，L，ωK｝—uk（t）中心頻率。

3.2 變分模型求解

（1）為獲得上述變分模型的最優(yōu)解，VMD通過(guò)引入二次罰因子 α 和拉格朗日乘數(shù) λ（t）構(gòu)建增廣拉格朗日函數(shù) L（｛uk｝，｛ωk｝，λ），把約束問(wèn)題轉(zhuǎn)化為無(wú)約束問(wèn)題。

其中，二次懲罰因子在高斯噪聲的存在下保證了信號(hào)重構(gòu)的準(zhǔn)確性，而拉格朗日乘子可以保證模型約束的剛度；

（2）用乘數(shù)交替方向法更新｛uk｝，｛ωk｝和 λ，并尋求增廣拉格朗日函數(shù)的鞍點(diǎn)，變分模型的最優(yōu)值顯示在方程（8）里；由此，輸入信號(hào)x（t）被分解成K個(gè)BLIMF分量。

4 基于PSR的KNN預(yù)測(cè)模型

傳統(tǒng)的KNN預(yù)測(cè)算法原理是對(duì)原n個(gè)數(shù)據(jù)構(gòu)造m個(gè)序列，找到其最鄰近的k個(gè)值，然后依據(jù)k個(gè)點(diǎn)的均值作為第（n+1）個(gè)數(shù)據(jù)的預(yù)測(cè)值[8]。算法認(rèn)為，系統(tǒng)的內(nèi)在復(fù)雜的映射關(guān)系均包含在歷史數(shù)據(jù)中，它既沒(méi)有針對(duì)歷史數(shù)據(jù)建立訓(xùn)練模型，也沒(méi)有將歷史數(shù)據(jù)作處理，特別適合于非線性問(wèn)題的快速分類[9]。

基于C_C算法對(duì)序列相空間重構(gòu)（Phasespacereconstruction，PSR），改進(jìn)了KNN算法對(duì)模型的決策性，由于空間重構(gòu)與原動(dòng)力系統(tǒng)方程微分同胚，能使相點(diǎn)分布表現(xiàn)出最優(yōu)空間支撐；將構(gòu)造的空間矩陣行向量作為KNN的建模域，每次計(jì)算得到的k個(gè)鄰近點(diǎn)均值作為狀態(tài)空間的下一個(gè)演化軌跡點(diǎn)；以預(yù)測(cè)值滾動(dòng)更新建模域，并重新計(jì)算嵌入維度m和延遲時(shí)間τ實(shí)現(xiàn)變參數(shù)預(yù)測(cè)，依此循環(huán)直至序列最大可信時(shí)間尺度。

具體預(yù)測(cè)步驟如下：（1）等間隔抽取測(cè)點(diǎn)的等長(zhǎng)短時(shí)周期振動(dòng)信號(hào)，構(gòu)造長(zhǎng)度為L(zhǎng)的振動(dòng)時(shí)序X0；（2）計(jì)算X0分解尺度并VMD分解，計(jì)算各BLIMF的最大預(yù)測(cè)可信時(shí)間尺度向量Ti=［T1，T2…Tn］；（3）提取各主模態(tài)不同時(shí)段的多重分形奇異譜特征值，形成特征序列；（4）創(chuàng)建不同模態(tài)特征序列初始建模域，通過(guò)C_C算法計(jì)算嵌入維度m和時(shí)延；（5）以該重構(gòu)相空間點(diǎn)列為初始待分類序列，KNN算法計(jì)算該相型分布的k個(gè)鄰近；（6）以K個(gè)鄰近的加權(quán)均值為該奇異譜特征預(yù)測(cè)值，建模域滾動(dòng)前移再預(yù)測(cè)；（8）重復(fù)（2）～（7）步；直至各模態(tài)的min｛Ti｝為截止時(shí)間尺度，獲得特征值預(yù)測(cè)序列。因此對(duì)于多重分形奇異譜特征向量而言，基于VMD和PSR（相空間重構(gòu)）的改進(jìn)KNN預(yù)測(cè)算法（稱為VMD＆PSR_KNN法）計(jì)算框圖，如圖2所示。

圖2 VMD＆PSR_KNN預(yù)測(cè)流程Fig.2 The Flowchart of VMD＆PSR_KNN Method

5 實(shí)例分析

結(jié)合以上分析與算法流程，課組對(duì)大慶天然氣公司壓氣站2D12型往復(fù)式壓縮機(jī)為研究對(duì)象；其中，壓縮機(jī)主要參數(shù)為：軸功率500kW，排氣量70m3/min，活塞行程240mm，電機(jī)轉(zhuǎn)速496r/min。采集設(shè)備大修前30天的大頭軸瓦敏感測(cè)點(diǎn)（缸體外曲軸箱處）振動(dòng)信號(hào)；數(shù)據(jù)采集系統(tǒng)為湖北優(yōu)泰公司的uT3416型16通道數(shù)據(jù)采集儀；壓縮機(jī)結(jié)構(gòu)與測(cè)點(diǎn)位置，如圖3所示。其中，黑色三角形表示二級(jí)軸瓦敏感測(cè)點(diǎn)位置。

5.1 算法參數(shù)確定

由于振動(dòng)信號(hào)內(nèi)部的分形本質(zhì)限定了預(yù)測(cè)只能是有限時(shí)間長(zhǎng)度的可信值；而基于故障識(shí)別的要求，預(yù)測(cè)得到的信號(hào)要展現(xiàn)至少一個(gè)周期內(nèi)的振動(dòng)特征，以滿足對(duì)其進(jìn)一步分析和處理的需求；同時(shí)，過(guò)分強(qiáng)調(diào)增加采樣間隔，似乎能增到預(yù)測(cè)可信時(shí)間，但必然會(huì)與復(fù)雜系統(tǒng)動(dòng)力學(xué)演化的長(zhǎng)期不可預(yù)測(cè)性特性相悖。

為保證采樣的有效長(zhǎng)度，對(duì)長(zhǎng)期采樣序列每間隔1h取得2周期約2410個(gè)采樣點(diǎn)（采樣頻率10kHz），獲得30h的間隔采樣數(shù)據(jù)，為便于程序化計(jì)算，采用MRMR（最大相關(guān)最小冗余）法計(jì)算各個(gè)采樣時(shí)段的VMD預(yù)分解個(gè)數(shù)均為4個(gè)，計(jì)算各模態(tài)序列的最大預(yù)測(cè)可信時(shí)間間隔向量為：T=［43Vt 33Vt 35Vt 37Vt］，可知在計(jì)算點(diǎn)后33h內(nèi)的預(yù)測(cè)值是可信的；KNN算法中取鄰近值K=3。

圖3 2D12往復(fù)式壓縮機(jī)結(jié)構(gòu)如圖Fig.3 The Structure of 2D12 Reciprocating Compressor

多重分形奇異譜特征量中，譜中心值α0（fmax＝f（α0））代表信號(hào)長(zhǎng)程相關(guān)性，其值越大相關(guān)性越強(qiáng)，其穩(wěn)定性雖好而可分性并不好[5-6]；通過(guò)構(gòu)造包含對(duì)稱性和奇異性的譜差比△αdr=（αmax-α0）/（α0-αmin）；以物理意義明確的α0和△αdr形成穩(wěn)定、可分的特征向量作為預(yù)測(cè)指標(biāo)，擬將其用于相空間重構(gòu)改進(jìn)的KNN預(yù)測(cè)模型加以預(yù)測(cè)。多重分形奇異譜形態(tài)參數(shù)示意，如圖4所示。

圖4 多重分形奇異譜形態(tài)參數(shù)示意Fig.4 Schematic Diagram of Shape Parameter for MSS

5.2 特征值預(yù)測(cè)結(jié)果

基于上述算法參數(shù)計(jì)算與分析，將該狀態(tài)時(shí)間序列經(jīng)VMD分解的4個(gè)BLIMF的譜參數(shù)α0和Δαdr計(jì)算結(jié)果與實(shí)測(cè)結(jié)果，如圖5、圖6所示。

圖5 各主模態(tài)奇異譜中心值Fig.5 Central Valueof MSS for Each BLIMF

圖6 各主模態(tài)奇異譜譜差比Fig.6 Difference Ratioof MSS for Each BLIMF

從各個(gè)模態(tài)分量的特征參數(shù)預(yù)測(cè)趨勢(shì)可見，譜寬和中心值兩個(gè)參數(shù)其幅值和波動(dòng)趨勢(shì)均與實(shí)測(cè)數(shù)據(jù)偏差較小，特別是前20個(gè)預(yù)測(cè)可信時(shí)間間隔以內(nèi)，預(yù)測(cè)值的吻合度較高。

5.3 方法評(píng)估

為證明算法有效性，課組采用線性擬合回歸的方式[4]；構(gòu)造方程為Ypredict=R×Xactual+Err，Ypredict和Xactual分別表示預(yù)測(cè)值和實(shí)測(cè)值，Err為擬和誤差，R為曲線斜率，其值越接近1則曲線預(yù)測(cè)效果越好，用相同數(shù)據(jù)分別對(duì)EMD_KNN、VMD_KNN和VMD&PSR_KNN三種預(yù)測(cè)方法做對(duì)比分析（其中在EMD分解中取前4個(gè)主模態(tài)，其信號(hào)能量比率已達(dá)到80%，已包含了高頻的特征信息[10]），兩特征參數(shù)回歸比較結(jié)果（圖略）。

回歸分析可見，EMD_KNN算法由于模態(tài)混疊等因素造成特征值預(yù)測(cè)結(jié)果不佳，而VMD法因模態(tài)分解方面由于EMD而使預(yù)測(cè)精度明顯提高，特別是結(jié)合空間重構(gòu)變參數(shù)建模后，算法考慮了時(shí)序熵值的變化的影響，預(yù)測(cè)回歸結(jié)果最優(yōu)。

表1 各方法a0預(yù)測(cè)值誤差比較Tab.1 Comparison of Predicted Errors for a0by Different Methods

表2 各方法Δadr預(yù)測(cè)值誤差比較Tab.2 Comparison of Predicted Errors for Δadrby Different Methods

如表1、表2所示對(duì)三種方法得到的特征參數(shù)預(yù)測(cè)誤差用標(biāo)準(zhǔn)差（STD）和均方差（RMSE）匯總并比較，兩種基于VMD分解的預(yù)測(cè)算法得到的預(yù)測(cè)特征值與實(shí)測(cè)序列計(jì)算的特征值誤差較小；可以證實(shí)，文中多重分型譜參數(shù)KNN預(yù)測(cè)算法的有效性；更為重要的是，從基于多重分形譜形參數(shù)衰減指標(biāo)程度判定，該故障已達(dá)到重度磨損閾值，這與停機(jī)大修檢查的實(shí)際結(jié)果是一致。

6 結(jié)論

（1）從熵值變化反映復(fù)雜系統(tǒng)初始敏感性角度將最大預(yù)測(cè)可信時(shí)間算法引入預(yù)測(cè)程序，優(yōu)化預(yù)測(cè)時(shí)間尺度；（2）從相空間重構(gòu)與原系統(tǒng)狀態(tài)方程微分同胚原理出發(fā)，利用VMD算法改進(jìn)KNN非參數(shù)預(yù)測(cè)模型缺陷，結(jié)合相空間重構(gòu)變尺度構(gòu)造建模域，再現(xiàn)系統(tǒng)動(dòng)力學(xué)演化趨勢(shì)；（3）基于多重分形理論，選擇可分性和魯棒性較好的奇異譜特征參數(shù)，結(jié)合VMD分解方法提取表征精細(xì)分形的特征向量，提高非參數(shù)預(yù)測(cè)算法理論的可解釋性。通過(guò)特征提取與實(shí)測(cè)分析證實(shí)，相比于EMD和傳統(tǒng)KNN算法，基于VMD多重分形奇異譜的變參數(shù)相空間重構(gòu)預(yù)測(cè)算法為設(shè)備健康狀態(tài)評(píng)估提供新的途徑和理論支撐。

［1］石博強(qiáng)，申焱華．機(jī)械故障診斷的分形方法：理論與實(shí)踐［M］．北京：冶金工業(yè)出版社，2001：288-300．（Shi Bo-qiang，Shen Yan-hua．Fractal Methods for Mechanical Fault Diagnosis：Theory and Practice［M］.BeiJing：Metallurgical Industry Press，2001：288-300.）

［2］劉巖，王金東，郭建華．基于參量?jī)?yōu)化的多變量相空間重構(gòu)方法［J］．機(jī)械傳動(dòng)，2013（6）：5-7．（Liu Yan，Wang Jin-dong，Guo Jian-hua．Method of Phase Space Reconstruction of Multivariate Data based on Parameters Optimization［J］.Journal of Mechanical Transmission，2013（6）：5-7．）

［3］周廣旭．一種新的時(shí)間序列分析算法及其在股票預(yù)測(cè)中的應(yīng)用［J］．計(jì)算機(jī)應(yīng)用，2005（9）：2179-2181．（Zhou Guang-xu．RBF-based tmi e-series forecasting［J］.Journal of Computer Applications，2005（9）：2179-2181．）

［4］Lin Ai-jing，Shang Peng-jian，F(xiàn)eng Guo-chen.Application of EMD combined with KNN approach in financial time series forecasting［J］.Fluctuation Noise Letters，2012（11）：1250018.

［5］楊茂，賈云彭，穆鋼.基于改進(jìn)KNN算法的風(fēng)電功率實(shí)時(shí)預(yù)測(cè)研究［J］.電測(cè)與儀表，2014（24）：38-43.（Yang Mao，Jia Yun-peng，Mu Gang.Wind Power Real-time Prediction Research Based on the Improved KNN Algorithm［J］.Electrical Measurement&Instrumentation，2014（24）：38-43.）

［6］J.Wang，W.Wei，H.Li.Multifractal features extraction of bearing faults for reciprocating compressor［J］.Journal of Vibration and Shock，2008，27（S）：313-315.

［7］Dragomiretskiy K and Zosso D．Variational mode decomposition IEEE Trans［C］.Signal Process，2014，62（3）：531-544.

［8］Bagui S C，Bagui S，Pal K.Breast Cancer Detection using Nearest Neighbor Classification Rules［J］.Pattern Recognition，2003（36）：25-34.

［9］Zeng Y，Yang Y，Zhou L.Pseudo Nearest Neighbor Rule for Pattern Recognition ［J］.Expert Systems with Applications，2009 （36）：3587-3595.

［10］張維強(qiáng)．基于MSKPCA和SVM的轉(zhuǎn)子故障診斷模型及應(yīng)用［J］.機(jī)械設(shè)計(jì)與制造，2015（10）：4-8.（Zhang Wei-qiang.Rotor fault diagnosis model and application based on multi-scale kernel principal component analysis and support vector machine［J］.Machinery Design&Manufacture，2015（10）：4-8.）