基于混沌理論的往復(fù)式壓縮機(jī)故障診斷

張來(lái)斌，陳敬龍，段禮祥

(1.中國(guó)石油大學(xué)機(jī)械與儲(chǔ)運(yùn)工程學(xué)院北京 102249;2.中原石油勘探局勘察設(shè)計(jì)研究院，濮陽(yáng) 457001)

基于混沌理論的往復(fù)式壓縮機(jī)故障診斷

張來(lái)斌1，陳敬龍2，段禮祥1

(1.中國(guó)石油大學(xué)機(jī)械與儲(chǔ)運(yùn)工程學(xué)院北京 102249;2.中原石油勘探局勘察設(shè)計(jì)研究院，濮陽(yáng) 457001)

實(shí)測(cè)某往復(fù)壓縮機(jī)氣缸在正常、活塞體松動(dòng)及拉缸工況下的振動(dòng)信號(hào)，計(jì)算信號(hào)的關(guān)聯(lián)維數(shù)、Kolmogorov熵及最大Lyapunov指數(shù)，并證明信號(hào)具有非線性。用奇異值分解(SVD)降噪法對(duì)信號(hào)進(jìn)行降噪，通過(guò)奇異熵增量曲線選擇降噪階次，用互信息法求最佳延時(shí)，并用假近鄰域法求最佳嵌入維數(shù)，同時(shí)用G-P算法、小數(shù)據(jù)量法求出關(guān)聯(lián)維數(shù)、Kolmogorov熵及最大Lyapunov指數(shù)。計(jì)算結(jié)果表明，3種工況下的往復(fù)壓縮機(jī)氣缸振動(dòng)信號(hào)均為非線性混沌信號(hào)，故障越嚴(yán)重，信號(hào)表現(xiàn)出的混沌特性越明顯，可依據(jù)混沌指標(biāo)對(duì)往復(fù)壓縮機(jī)氣缸的狀態(tài)進(jìn)行監(jiān)測(cè)。

往復(fù)壓縮機(jī);混沌理論;非線性信號(hào);故障診斷

往復(fù)壓縮機(jī)是油田廣泛使用的設(shè)備，但其易損部件多。活塞和缸壁發(fā)生碰磨可導(dǎo)致氣缸損傷，嚴(yán)重時(shí)會(huì)引起活塞桿斷裂、連桿斷裂及曲軸變形，甚至導(dǎo)致整臺(tái)機(jī)組報(bào)廢，造成巨大的經(jīng)濟(jì)損失。對(duì)往復(fù)壓縮機(jī)的氣缸進(jìn)行狀態(tài)監(jiān)測(cè)非常重要，其中通過(guò)振動(dòng)信號(hào)對(duì)氣缸進(jìn)行故障診斷是較常用的方法。但是，往復(fù)壓縮機(jī)振源多、氣缸振動(dòng)信號(hào)具有強(qiáng)非線性，用常規(guī)的頻譜分析難以對(duì)其做出準(zhǔn)確的診斷?；煦缋碚搶儆诜蔷€性科學(xué)，適合于處理非線性信號(hào)［1-2］，可將混沌理論用于往復(fù)壓縮機(jī)的故障診斷，通過(guò)混沌指標(biāo)對(duì)往復(fù)壓縮機(jī)進(jìn)行狀態(tài)監(jiān)測(cè)。往復(fù)壓縮機(jī)氣缸振動(dòng)信號(hào)是否為非線性信號(hào)必須加以證明，否則計(jì)算信號(hào)的混沌特性將無(wú)意義。目前尚無(wú)學(xué)者研究往復(fù)壓縮機(jī)的氣缸振動(dòng)信號(hào)是否為非線性，因此筆者采用替代數(shù)據(jù)法對(duì)其進(jìn)行證明。

1 信號(hào)降噪

設(shè)信號(hào)X=(x1，x2，…，xN)，將其用Takens嵌入定理進(jìn)行重構(gòu)［3］:

式中，t=1，2，…，m;m=N－(n－1)τ;n為嵌入維數(shù);τ為時(shí)間延遲。

根據(jù)奇異值分解理論，A滿足以下關(guān)系式:

矩陣S的非對(duì)角元素全為0，其對(duì)角元素λ1≥λ2≥…≥λl≥0，λi(i=1，2，…，l)是矩陣A的奇異值。將較小的λi置為0，構(gòu)成一個(gè)新的對(duì)角矩陣S'，用S'代替S，代入式(2)中算出一個(gè)新的矩陣A'，通過(guò)構(gòu)造相空間的逆過(guò)程，從A'中可獲取降噪后的信號(hào)。通過(guò)奇異熵法可確定降噪階次，奇異熵的計(jì)算公式［4］為

式中，Ek為階次為k時(shí)的奇異熵;ΔEi為奇異熵在階次為i時(shí)的增量。

當(dāng)ΔEi曲線開(kāi)始下降并趨向于一個(gè)較小的穩(wěn)定值時(shí)，選擇此時(shí)的i作為降噪階次。

為驗(yàn)證奇異值分解降噪的效果，將其用于仿真信號(hào)的降噪。在Lorenz信號(hào)上疊加白噪聲信號(hào)，圖1為含噪信號(hào)的偽相圖。用奇異值分解降噪對(duì)信號(hào)進(jìn)行降噪處理，圖2為降噪后信號(hào)的偽相圖?？煽闯龌煦缥幼兊们逦？梢?jiàn)，奇異值分解降噪具有很好的降噪效果。

圖1 含噪信號(hào)的偽相圖Fig.1 Pseudo-phase portrait of original signal

圖2 SVD降噪后信號(hào)的偽相圖Fig.2 Pseudo-phase portrait for noise reduction by SVD

實(shí)測(cè)的某往復(fù)式壓縮機(jī)2號(hào)氣缸3種工況下的加速度振動(dòng)信號(hào)如圖3所示。采樣頻率為16 kHz，采樣長(zhǎng)度為10240個(gè)點(diǎn)。工況1為正常，工況2為活塞體松動(dòng)，工況3為活塞體與缸壁碰磨。用奇異值分解降噪對(duì)原始信號(hào)進(jìn)行降噪處理。

圖3 3種工況下氣缸的原始振動(dòng)信號(hào)Fig.3 Original cylinder vibration signals under three work conditions

2 最佳延時(shí)及嵌入維數(shù)

2.1 最佳延時(shí)的計(jì)算

自相關(guān)法是計(jì)算最佳延時(shí)的常用方法［5］，但自相關(guān)函數(shù)只能表征信號(hào)中的線性關(guān)系?；バ畔⒛鼙碚餍盘?hào)中的非線性，故筆者選用互信息法［6-7］計(jì)算最佳延時(shí)。

延時(shí)τ=kΔt，Δt為相鄰采樣點(diǎn)的時(shí)間間隔。令［s，q］=［x(n)，x(n+k)］，在(s，q)平面上用一個(gè)矩形劃出所求的所有點(diǎn)，s方向的長(zhǎng)度為Δs，q方向的長(zhǎng)度為Δq。等間隔劃分Δs和Δq，s方向格子數(shù)為M，格子長(zhǎng)度為εs;q方向格子數(shù)為M'，格子長(zhǎng)度為εq。(a，b)是所劃區(qū)域的起點(diǎn)，做如下判斷:若(i－1)εs≤s－a＜iεs，則s在第i個(gè)格子中，對(duì)Ni做一次記錄;若(j－1)εq≤q－b＜jεq，則q在第j個(gè)格子中，對(duì)Nj做一次記錄;若(i－1)εs≤s－a＜iεs且(j－1)εq≤q－b＜jεq，則(s，q)在標(biāo)號(hào)為(i，j)的格子中，對(duì)Nij做一次記錄。Ntotal為矩形區(qū)域中的總點(diǎn)數(shù)，可得ps(i)=Ni/Ntotal，pq(j)=Nj/Ntotal，psq(i，j)= Nij/Ntotal。

互信息的計(jì)算公式為

其中

互信息曲線第一次下降到極小值所對(duì)應(yīng)的延時(shí)為最佳延時(shí)。格子數(shù)取為80×80，經(jīng)計(jì)算，3種工況下的最佳延時(shí)分別為50、40和32。

2.2 最佳嵌入維數(shù)的計(jì)算

選用假近鄰域法［8］計(jì)算最佳嵌入維數(shù)，傳統(tǒng)的假近鄰域法有多處需要人為選取閾值，不利于計(jì)算，故用改進(jìn)的假近鄰域法［9］計(jì)算最佳嵌入維數(shù)。

對(duì)時(shí)間序列(x1，x2，…，xN)進(jìn)行相空間重構(gòu):

式中，m為嵌入維數(shù);k為時(shí)間延遲。

計(jì)算距離矩陣R，共有M(M－1)/2個(gè)元素，其中

求取R的最小值Rmin和最大值Rmax，則

式中，p0為閾值參數(shù);p0∈［0，1］，本文中選取p0=0.1。

令m=m1，統(tǒng)計(jì)鄰域點(diǎn)總數(shù)，再令m=m1+1，計(jì)算m=m1時(shí)的鄰域點(diǎn)是否為真鄰域點(diǎn)，統(tǒng)計(jì)非鄰域點(diǎn)總數(shù)，并計(jì)算非鄰域點(diǎn)總數(shù)與鄰域點(diǎn)總數(shù)的比值。不斷重復(fù)這個(gè)過(guò)程，直到比值小于0.05，或不再隨m的增加而增加，取此時(shí)的m為最佳嵌入維數(shù)。經(jīng)計(jì)算3種工況下氣缸振動(dòng)信號(hào)的最佳嵌入維數(shù)為13、16和20。

3 關(guān)聯(lián)維數(shù)及Kolmogorov熵

3.1 關(guān)聯(lián)維數(shù)的計(jì)算

Grassberger和Procacciay提出了G-P算法［10］，可用其計(jì)算關(guān)聯(lián)維數(shù)。對(duì)時(shí)間序列進(jìn)行相空間重構(gòu)后，計(jì)算向量的距離。計(jì)算有關(guān)聯(lián)的向量對(duì)數(shù)，它在一切可能的M(M－1)種配對(duì)中所占的比例稱為關(guān)聯(lián)積分表達(dá)式，表示為

式中，θ(·)為Heaviside單位函數(shù)，表示為

根據(jù)式(6)可計(jì)算關(guān)聯(lián)積分Cn(r)。當(dāng)r在某一范圍時(shí)，關(guān)聯(lián)積分Cn(r)與r存在以下關(guān)系:

圖4為3種工況下氣缸振動(dòng)信號(hào)的關(guān)聯(lián)積分曲線，圖5為相應(yīng)的局部斜率曲線。

圖4 3種工況下氣缸振動(dòng)信號(hào)的關(guān)聯(lián)積分曲線Fig.4 Correlation integral curves of cylinder vibration signals under three work conditions

圖5 局部斜率曲線Fig.5 Local slope curves

根據(jù)局部斜率曲線，用最小二乘法對(duì)關(guān)聯(lián)積分曲線的直線部分進(jìn)行計(jì)算，直線斜率即為關(guān)聯(lián)維數(shù)。3種工況下的關(guān)聯(lián)維數(shù)分別為1.1853、1.195 5和1.2951。隨著故障程度的加劇，關(guān)聯(lián)維數(shù)增大。這是因?yàn)楣收显絿?yán)重，系統(tǒng)的耗散能越大，引起能量變化的作用力也就越大。

3.2 Kolmogorov熵

Grassberger和Procacciay提出了計(jì)算Kolmogorov熵的算法［11］。Kolmogorov熵的定義如下:設(shè)一個(gè)動(dòng)力系統(tǒng)的自由度為F，將F-維相空間劃分為一系列的盒子，盒子的尺寸為εF，假設(shè)相空間中存在一個(gè)吸引子，每隔τ對(duì)系統(tǒng)的狀態(tài)進(jìn)行測(cè)量。令p(i1，i2，…，id)為在盒子i1中在盒子i2中，…，且在盒子id中的聯(lián)合概率，則Kolmogorov熵的表達(dá)式為

文獻(xiàn)［11］中已證明，在某些情況下K2在數(shù)值上接近于K。計(jì)算K是非常困難的，但計(jì)算信號(hào)的K2較容易，因此可以用K2估計(jì)K。K2的計(jì)算公式為

當(dāng)d→∞，m→∞，ε→0時(shí)，K2→K。

用G-P算法計(jì)算關(guān)聯(lián)積分時(shí)，可同時(shí)算出Kolmogorov熵。3種工況下氣缸振動(dòng)信號(hào)的Kolmogorove熵分別為7.68×10－4、8.32×10－4和8.53×10－4。隨著故障程度的加劇，Kolmogorove熵也增大，這是因?yàn)橐鹣到y(tǒng)能量耗散的沖擊力增大。

4 最大Lyapunov指數(shù)

最大Lyapunov指數(shù)是混沌特性的一個(gè)重要的指標(biāo)，如果信號(hào)有最大Lyapunov指數(shù)且為正值，則說(shuō)明該系統(tǒng)是混沌的。因?qū)嶋H獲取的信號(hào)長(zhǎng)度有限，故選用小數(shù)據(jù)量法［12］計(jì)算最大 Lyapunov指數(shù)。

對(duì)信號(hào)進(jìn)行相空間重構(gòu)，計(jì)算相空間每個(gè)點(diǎn)的最短初始距離，并限制短暫分離。

其中p為平均周期，可通過(guò)傅里葉變換估算出來(lái)。

式中，fi(i=1，2，…，L)為0到奈奎斯特頻率之間均勻分布的頻率點(diǎn);Ai為對(duì)應(yīng)于fi的幅值。

對(duì)相空間中的每個(gè)點(diǎn)，計(jì)算該鄰點(diǎn)對(duì)的i個(gè)離散時(shí)間步后的距離di(j)。

式中，i=1，2，…，min(M－j，M。對(duì)每個(gè)i，求出所有j的lgdi(j)的平均值y(i)。

式中，h為非零di(j)的數(shù)目。

用最小二乘法對(duì)曲線y(i)的直線部分做線性回歸，直線的斜率就是最大Lyapunov指數(shù)λ1。

圖6為求取氣缸振動(dòng)信號(hào)的最大Lyapunov指數(shù)曲線，工況1～工況3的最大Lyapunov指數(shù)分別為4.319 5、4.444 8和4.613 1。3種工況下的最大Lyapunov指數(shù)都是正值，且故障越嚴(yán)重，最大Lyapunov指數(shù)越大。可見(jiàn)，隨著故障程度的加劇，系統(tǒng)的混沌特性越明顯。

圖6 求取最大Lyapunov指數(shù)的曲線Fig.6 The largest Lyapunov exponent curves

5 氣缸振動(dòng)信號(hào)的非線性證明

Theiler等［13］提出了用替代數(shù)據(jù)法檢驗(yàn)時(shí)間序列中的非線性，雷敏［14］對(duì)該算法進(jìn)行了改進(jìn)。

對(duì)時(shí)間序列進(jìn)行傅里葉變換。

對(duì)相位進(jìn)行隨機(jī)化，新生成的相位φ(f)∈［－π，π］。當(dāng)信號(hào)長(zhǎng)度為奇數(shù)時(shí)，φ(f1=0)，φ(fi)=－φ(fk)，i=1～(N+1)/2，k=N～(N+1)/2+1。當(dāng)信號(hào)長(zhǎng)度為偶數(shù)時(shí)，φ(f1=0)，φ(fN/2+1)=0，φ(fi)=－φ(fk)，i=2～N/2，k=N～N/2+2。引入以下判據(jù)［15］:

式中，Dorig為原始數(shù)據(jù)的關(guān)聯(lián)維數(shù);ˉDsurr為所有替代數(shù)據(jù)關(guān)聯(lián)維數(shù)的平均值;σsurr為所有替代數(shù)據(jù)關(guān)聯(lián)維數(shù)的標(biāo)準(zhǔn)差。若取顯著水平α=0.05，則當(dāng)S≥1.96時(shí)，原始數(shù)據(jù)以95%的置信水平為非線性序列;當(dāng)S＜1.96時(shí)，原始數(shù)據(jù)以95%的置信水平為隨機(jī)序列。

各工況下都生成10組替代數(shù)據(jù)，計(jì)算其關(guān)聯(lián)維數(shù)，并計(jì)算各工況下關(guān)聯(lián)維數(shù)的均值和標(biāo)準(zhǔn)差，結(jié)果見(jiàn)表1。將表1數(shù)據(jù)代入式(12)，可算得各工況下的S。

表1 3種工況下關(guān)聯(lián)維數(shù)的統(tǒng)計(jì)量Table 1 Statistics of correlation dimensions under three work conditions

3種工況下的判據(jù)S都大于1.96，可見(jiàn)，往復(fù)壓縮機(jī)的氣缸振動(dòng)信號(hào)是非線性混沌信號(hào)。

6 結(jié)論

(1)計(jì)算的3種工況下某往復(fù)壓縮機(jī)氣缸振動(dòng)信號(hào)的混沌指標(biāo)包括關(guān)聯(lián)維數(shù)、Kolmogorov熵和最大Lyapunov指數(shù)。故障越嚴(yán)重，這3個(gè)指標(biāo)的值越大，且3種工況下的最大Lyapunov指數(shù)都為正值。往復(fù)壓縮機(jī)的氣缸振動(dòng)信號(hào)具有混沌特性，且故障越嚴(yán)重，混沌特性越明顯?；煦缰笜?biāo)對(duì)往復(fù)壓縮機(jī)氣缸的工況很敏感，可通過(guò)混沌指標(biāo)對(duì)往復(fù)壓縮機(jī)氣缸進(jìn)行狀態(tài)監(jiān)測(cè)。

(2)用替代數(shù)據(jù)法證明了3種工況下往復(fù)壓縮機(jī)氣缸振動(dòng)信號(hào)是非線性混沌信號(hào)，可用混沌理論對(duì)其進(jìn)行分析。

［1］吳光強(qiáng)，盛云.混沌理論在汽車非線性系統(tǒng)中的應(yīng)用進(jìn)展［J］.機(jī)械工程學(xué)報(bào)，2010，46(10):81-87.

WU Guang-qiang，SHENG Yun.Review on the application of chaos theory in automobile nonlinear system［J］.Journal of Mechanical Engineering，2010，46(10):81-87.

［2］韓敏.混沌時(shí)間序列預(yù)測(cè)理論與方法［M］.北京:中國(guó)水利水電出版社，2007:97-150.

［3］SHINK，HAMMONDJK，WHITE P R.Iterative SVD method for noise reduction of low dimension chaotic time series［J］.Mechanical Systems and Signal Processing，1999，13(1):115-124.

［4］YANG Wen-xian，PETERW.Development of an advanced noise reduction method for vibration analysis based on singular value decomposition［J］.NDT and E International，2003，36(6):419-432.

［5］陳鏗，韓伯棠.混沌時(shí)間序列分析中的相空間重構(gòu)技術(shù)綜述［J］.計(jì)算機(jī)科學(xué)，2005，32(4):67-70.

CHENKeng，HANBo-tang.A survey of state spaee reconstruction of chaotic time series analysis［J］.Computer Science，2005，32(4):67-70.

［6］LIEBERT W，SCHUSTERH G.Proper choice of the time delay for the analysis of chaotic time series［J］.Phys Lett A，1989，142:107-111.

［7］楊志安，王光瑞，陳式剛.用等間距分格子法計(jì)算互信息函數(shù)確定延遲時(shí)間［J］.計(jì)算物理，1995，12(4):442-448.

YANG Zhi-an，WANG Guang-rui，CHENShi-gang.Determination of delay time by calculating mutual information with equally distant space cells［J］.Chinese Journal of Computational Physics，1995，12(4):442-448.

［8］KENNELMathew B，BROWNRAbarbanel H D I.Determining embedding dimension for phase-space reconstruction using a geometrical construction［J］.Phy Rev A，1992，45:3403-3411.

［9］岳毅宏，韓文秀，王健.基于相點(diǎn)距離集的相空間重構(gòu)嵌入維數(shù)確定法［J］.機(jī)械工程學(xué)報(bào)，2005，41(10):35-38.

YUE Yi-hong，HANWen-xiu，WANG Jian.Embedding dimension in phase-space reconstruction based on distance set of phase points［J］.Chinese Journal of Mechanical，2005，41(10):35-38.

［10］GRASSBERGERP，PROCACCIA I.Characterization of strange attractors［J］.Phys Rev LettA，1983，50(5): 346-349.

［11］GRASSBERGERP，PORCACCIA I.Estimation of the Kolmogorov entropy from a chaotic signal［J］.Phys Rev A，1983，28(4):2591-2593.

［12］MICHAELT Rosenstein，JAMESJCollins，CARLO JDe Luca.A practical method for calculating largest Lyapunov exponents from small data sets［J］.Physical D，1993，65:117-134.

［13］THEILERJ，EUBANK S，LONGTINA，et al.Testing for nonlinearity in time series:the method of surrogate data［J］.Physica D，1992，58:77-94.

［14］MINLei，WANG Zhi-zhong，F(xiàn)ENG Zheng-jin.Detecting nonlinearity of action surface EMG signal［J］.Physics Letters A，2001，290:297-303.

［15］ROMBOUTSS，KEUNENR.Investigation of nonlinear structure in multichannel EEG［J］.Phys Lett A，1995，202:352-358.

Fault diagnosis of reciprocating compressor based on chaos theory

ZHANG Lai-bin1，CHENJing-long2，DUANLi-xiang1

(1.College of Mechanical and Transportation Engineering in China University of Petroleum，Beijing102249，China; 2.Survey and Design＆Research Institute，Zhongyuan Petroleum Exploration Bureau，Puyang457001，China)

The vibration signals of a cylinder were measured in three work conditions，including normal operation，loosening of piston and piston-liner wear，then chaotic indexes were calculated and the nonlinear characteristics of these signals were proved.Firstly，the original signals were decomposed using singular value decomposition(SVD)，and a reasonable order for noise reduction was selected according to the singular entropy of singular spectrum，and the best delay was calculated using mutual information method，and the best embedding dimension was calculated using fault near neighbour method.Then correlation dimensions and Kolmogorov entropies were calculated using G-P algorithm simultaneously，and the largest Lyapunov exponents were calculated using small-data method.The results show that the cylinder vibration signals are nonlinear chaotics signal，and the more serious the fault is，the more apparent the chaotic character is.The vibration condition monitoring for the reciprocating compressor can be implemented using chaotic indexes.

reciprocating compressor;chaos theory;nonlinear signal;fault diagnosis

TH 17

A

10.3969/j.issn.1673-5005.2012.01.019

1673-5005(2012)01-0112-05

2011－05－23

國(guó)家“863”計(jì)劃項(xiàng)目(2008AA06Z209);中國(guó)石油天然氣集團(tuán)公司創(chuàng)新基金項(xiàng)目(07E1005)

張來(lái)斌(1961－)，男(漢族)，安徽銅陵人，教授，博士，博士生導(dǎo)師，現(xiàn)從事石油機(jī)械工程方面的教學(xué)與科研工作。

(編輯 沈玉英)