李美威1 謝小鵬1 馮 偉 賀石中

(1.華南理工大學(xué)機(jī)械與汽車(chē)工程學(xué)院 廣東廣州 510640； 2.廣州機(jī)械科學(xué)研究院有限公司設(shè)備狀態(tài)檢測(cè)研究所 廣東廣州 510700)

在線(xiàn)油液監(jiān)測(cè)是指在設(shè)備不停機(jī)的情況下，通過(guò)各種在線(xiàn)傳感器對(duì)油液理化指標(biāo)和鐵磁磨粒等進(jìn)行實(shí)時(shí)監(jiān)測(cè)，根據(jù)監(jiān)測(cè)參數(shù)來(lái)判定設(shè)備的運(yùn)行情況，并進(jìn)行相應(yīng)的故障診斷[1]。早期的異常檢測(cè)能夠避免更多的嚴(yán)重后果,保證系統(tǒng)性能和效率并減少維修費(fèi)用。在整個(gè)在線(xiàn)監(jiān)測(cè)過(guò)程中，當(dāng)系統(tǒng)嚴(yán)重偏離其正常情形時(shí)就會(huì)導(dǎo)致異常發(fā)生[2]，異常檢測(cè)就是識(shí)別出監(jiān)測(cè)數(shù)據(jù)偏離正常情況的程度。對(duì)在線(xiàn)油液監(jiān)測(cè)數(shù)據(jù)進(jìn)行趨勢(shì)分析和預(yù)測(cè)，以提前發(fā)現(xiàn)異常并報(bào)警，可有效避免故障發(fā)生。目前，常用的趨勢(shì)分析方法有回歸擬合、統(tǒng)計(jì)分析、支持向量機(jī)、時(shí)間序列分析、灰色理論和人工神經(jīng)網(wǎng)絡(luò)等。趨勢(shì)預(yù)測(cè)的方法有很多，但是適用于油液監(jiān)測(cè)系統(tǒng)的較少[3]。

針對(duì)油液檢測(cè)特征的趨勢(shì)分析，國(guó)內(nèi)外學(xué)者都展開(kāi)了大量研究。林麗等人[4]基于油液在線(xiàn)監(jiān)測(cè)的磨粒信息對(duì)齒輪箱磨損狀態(tài)進(jìn)行了診斷和預(yù)測(cè)。張紅和龔玉[5]以光譜分析為例，運(yùn)用灰色理論對(duì)磨損趨勢(shì)進(jìn)行了預(yù)測(cè)。高經(jīng)緯等[6]對(duì)內(nèi)燃機(jī)潤(rùn)滑油中的各元素運(yùn)用時(shí)序模型進(jìn)行了趨勢(shì)分析，以判斷內(nèi)燃機(jī)的磨損狀態(tài)。TOBON-MEJIA等[7]使用小波分解技術(shù)和基于高斯混合隱馬爾科夫來(lái)估計(jì)軸承的壽命。SOUALHI等[8]運(yùn)用基于向量回歸機(jī)的方法對(duì)軸承進(jìn)行了監(jiān)測(cè)分析。上述研究主要是通過(guò)離線(xiàn)油液分析或是對(duì)一個(gè)周期內(nèi)某個(gè)特征進(jìn)行實(shí)驗(yàn)擬合分析再運(yùn)用于在線(xiàn)，這往往需要大量的時(shí)間去收集足夠全面的數(shù)據(jù)，而且用某個(gè)周期內(nèi)的特征變化代表所有情況是不準(zhǔn)確的。

在線(xiàn)監(jiān)測(cè)數(shù)據(jù)往往是不平衡的，能監(jiān)測(cè)到的故障數(shù)據(jù)比正常數(shù)據(jù)通常要少得多，這使得許多傳統(tǒng)的診斷方法并不適用于這種問(wèn)題[9]。而各種設(shè)備故障的發(fā)生往往會(huì)有一個(gè)潛在過(guò)程。針對(duì)在線(xiàn)數(shù)據(jù)的這些特點(diǎn)，本文作者對(duì)設(shè)備運(yùn)行平穩(wěn)期的監(jiān)測(cè)數(shù)據(jù)構(gòu)建ARMA模型，然后對(duì)整個(gè)時(shí)期的模型殘差進(jìn)行分析，將監(jiān)測(cè)數(shù)據(jù)特征量在故障發(fā)生之前的一段時(shí)間內(nèi)的殘差劃分為平穩(wěn)期和故障潛伏期，設(shè)定故障潛伏期殘差界限值，一旦殘差越界即可提前報(bào)警。

1 ARMA模型

1.1 模型定義

自回歸滑動(dòng)平均(Auto Regression Moving Average，ARMA)模型是研究平穩(wěn)隨機(jī)過(guò)程的典型方法。ARMA模型認(rèn)為一個(gè)時(shí)間序列的相互依存關(guān)系表現(xiàn)在原始數(shù)據(jù)的延續(xù)性上，在某時(shí)刻的值受到歷史值和噪聲的影響[10-11]。對(duì)于一個(gè)離散的時(shí)間序列{x1,x2,......,xn,......}，ARMA數(shù)學(xué)模型表示為

xt=μ+φ1xt-1+φ2xt-2+......+φpxt-p+εt-

θ1εt-1-......-θqεt-q

(1)

式中：xt是當(dāng)前值;μ是常數(shù)項(xiàng);p和q是模型階數(shù)；φi和θi是模型參數(shù)；{εt}是白噪聲序列。

當(dāng)q=0時(shí)，ARMA(p,q)模型就退化成自回歸AR(p)模型，即

xt=μ+φ1xt-1+φ2xt-2+......+φpxt-p+εt

(2)

當(dāng)p=0時(shí)，ARMA(p,q)模型就退化成移動(dòng)平均MA(q)模型，即

xt=μ+εt-θ1εt-1-......-θqεt-q

(3)

1.2 模型定階和參數(shù)估計(jì)

在使用ARMA(p,q)模型對(duì)時(shí)間序列進(jìn)行擬合時(shí)，首先需要對(duì)時(shí)間序列進(jìn)行平穩(wěn)非白噪聲檢驗(yàn)?？墒褂脝挝桓鶛z驗(yàn)和Ljung-Box檢驗(yàn)進(jìn)行序列的平穩(wěn)非白噪聲檢驗(yàn)。ADF檢驗(yàn)通過(guò)檢驗(yàn)序列中是否存在單位根來(lái)判斷序列的平穩(wěn)性[12]，如果序列存在單位根則為非平穩(wěn)序列，否則就是平穩(wěn)序列。Ljung-Box是對(duì)時(shí)間序列是否存在滯后相關(guān)的一種檢驗(yàn)方法[13]，如果序列存在滯后相關(guān)則為非白噪聲序列，否則為白噪聲序列。

檢驗(yàn)完序列為平穩(wěn)非白噪聲后，就可使用ARMA(p,q)模型進(jìn)行分析。首先需要確定模型中的p和q值。ARMA(p,q)模型的階數(shù)主要是根據(jù)序列的自相關(guān)系數(shù)和偏自相關(guān)系數(shù)的拖尾性和截尾性來(lái)判斷。其對(duì)應(yīng)的模型如表1所示。

表1 ARMA模型定階基本原則

確定模型階數(shù)后，估計(jì)模型中的參數(shù)。使用最小二乘估計(jì)法對(duì)模型參數(shù)進(jìn)行估計(jì)[14]，最小二乘法能充實(shí)使用序列值，精度較高。此時(shí)，記

a=(μ,φ1,φ2,......,φp,θ1,θ2,......,θq)T

(4)

Ft(a)=μ+φ1xt-1+φ2xt-2+......+φpxt-p-

θ1εt-1-......-θqεt-q

(5)

計(jì)算殘差：

εt=xt-Ft(a)

(6)

計(jì)算殘差平方和：

φpxt-p+θ1εt-1+......+θqεt-q)2

(7)

使式(7)達(dá)到最小的參數(shù)值a即為參數(shù)的估計(jì)值。由于白噪聲序列的值無(wú)法確定，所以需要運(yùn)用迭代的方法進(jìn)行計(jì)算。

1.3 模型檢驗(yàn)與優(yōu)化

在模型構(gòu)建完成以后，就可以得到擬合殘差序列。如果殘差序列為白噪聲，則認(rèn)為擬合模型充分提取了原始序列中的所有信息，這樣模型才顯著有效。所以需要使用Ljung-Box方法對(duì)殘差序列進(jìn)行白噪聲檢驗(yàn)以判斷模型的有效性。

在進(jìn)行模型定階時(shí)只是主觀判斷相關(guān)系數(shù)的拖尾性和截尾性來(lái)確定p和q的值，但該模型不一定是最優(yōu)模型。選用最小信息量準(zhǔn)則(Akaike Information Criterion，AIC)來(lái)進(jìn)行模型優(yōu)化，AIC準(zhǔn)則通過(guò)模型的似然函數(shù)值和參數(shù)個(gè)數(shù)來(lái)衡量模型的擬合效果[15]，即

(8)

通過(guò)比較多個(gè)模型的AIC值，即可選取AIC最小的模型作為最優(yōu)模型。

1.4 故障潛伏期界限值設(shè)定

K均值是一種聚類(lèi)方法。對(duì)于一維K均值，其算法如下：

(1)隨機(jī)選取序列中的2個(gè)數(shù)作為初始中心點(diǎn)k1、k2；

(2)分別計(jì)算所有點(diǎn)到k1和k2的距離d1和d2，若d1≤d2，則記為1類(lèi)，否則為2類(lèi)；

(3)更新中心點(diǎn)，將所有1類(lèi)點(diǎn)的均值賦給k1，所有2類(lèi)點(diǎn)的均值賦給k2；

(4)重復(fù)步驟(2)、(3)，直至k1和k2不再變化。

由于數(shù)據(jù)是一維時(shí)間序列，即可使用兩中心點(diǎn)的均值作為界限值，即

thresholdε=(k1+k2)/2

(9)

2 實(shí)驗(yàn)驗(yàn)證與分析

2.1 實(shí)驗(yàn)數(shù)據(jù)

以某水電站水輪機(jī)組為研究對(duì)象，用實(shí)時(shí)獲取的推力油槽油液含水量進(jìn)行模型驗(yàn)證和分析。該機(jī)組安裝了廣研檢測(cè)的在線(xiàn)油液監(jiān)測(cè)儀，可實(shí)時(shí)采集機(jī)組油液的各項(xiàng)特征，采樣周期為1 h。該設(shè)備在運(yùn)行過(guò)程中遭受了水污染，油液含水量大幅上升。水分會(huì)使油液乳化，降低油液黏度和油膜厚度，還會(huì)促使油品氧化，加速腐蝕，惡化油質(zhì)，所以應(yīng)保證油液中的含水量盡可能地低。文中選取一次含水量超標(biāo)的數(shù)據(jù)作為樣本數(shù)據(jù)進(jìn)行分析，總共454條數(shù)據(jù)，根據(jù)油液標(biāo)準(zhǔn)設(shè)定的含水量界限值為1×10-4(質(zhì)量分?jǐn)?shù)，下同),含水量實(shí)時(shí)趨勢(shì)如圖1所示。

圖1 原始數(shù)據(jù)時(shí)序Fig 1 Original time series

2.2 平穩(wěn)序列建模

使用ADF檢驗(yàn)序列平穩(wěn)性，其檢驗(yàn)結(jié)果如表2所示。

表2 序列檢驗(yàn)結(jié)果

表3 序列檢驗(yàn)結(jié)果

2.2.2 模型定階與參數(shù)估計(jì)

圖2 序列自相關(guān)函數(shù)Fig 2 Time autocorrelation function (ACF)

圖3 序列偏相關(guān)函數(shù)Fig 3 Time partial autocorrelation function (ACF)

從圖2和圖3中可看出，其自相關(guān)表現(xiàn)出明顯的拖尾性，偏相關(guān)在9階之后取值均在置信區(qū)間內(nèi)，可認(rèn)為偏相關(guān)9階截尾，根據(jù)表1可判定為AR(9)模型。再分別對(duì)AR(7)、AR(8)、AR(9)、AR(10)、AR(11)進(jìn)行分析，比較它們的AIC，選取最優(yōu)模型。5種模型的AIC值如表4所示。

表4 各模型的AIC值對(duì)比

從表4可知，AR(9)模型的AIC值最小，則選取AR(9)模型進(jìn)行序列擬合，并運(yùn)用最小二乘法對(duì)AR(9)模型進(jìn)行參數(shù)估計(jì)，其參數(shù)估計(jì)結(jié)果如表5所示

表5 參數(shù)估計(jì)結(jié)果

則模型的表達(dá)式為

xt=26.288 9+0.639 8(xt-1-26.288 9)+

0.123 6(xt-2-26.288 9)-0.004 1(xt-3-26.288 9)-

0.020 6(xt-4-26.288 9)+0.003 8(xt-5-26.288 9)+

0.027 3(xt-6-26.288 9)+0.378 8(xt-7-26.288 9)-

0.243 5(xt-8-26.288 9)-0.160 6(xt-9-26.288 9)+

εt=6.716 8+0.639 8xt-1+0.123 6xt-2-0.004 1xt-3-

0.020 6xt-4+0.003 8xt-5+0.027 3xt-6+0.378 8xt-7-

0.243 5xt-8-0.160 6xt-9+εt

圖4 AR(9)模型擬合結(jié)果和原始數(shù)據(jù)比較Fig 4 Comparison of fitting results of AR(9) model and original data

2.2.3 模型檢驗(yàn)

驗(yàn)證擬合模型的殘差序列是否為白噪聲序列。同理，使用Ljung-Box檢驗(yàn)方法對(duì)殘差序列的3、6、9期和12期滯后進(jìn)行自相關(guān)分析，其結(jié)果如表6所示。

表6 殘差Ljung-Box檢驗(yàn)結(jié)果

從表6中可知?dú)埐钚蛄械?、6、9和12期滯后的檢驗(yàn)概率遠(yuǎn)大于0.05，并且均接近于1，則可認(rèn)為殘差序列為白噪聲序列，模型有效。

3 殘差序列分析

圖5 原始時(shí)間序列及其擬合結(jié)果Fig 5 Original time series and its fitting results

3.1 對(duì)殘差序列進(jìn)行統(tǒng)計(jì)分析

圖6 不同時(shí)期殘差統(tǒng)計(jì)直方圖Fig 6 Residual statistical histogram of different period

從圖6中可以看出，平穩(wěn)期殘差是一個(gè)白噪聲序列，大致符合正態(tài)分布。而故障潛伏期的殘差分布發(fā)生了較大變化，對(duì)故障潛伏期的殘差序列進(jìn)一步分析劃分界限值。

3.2 對(duì)殘差進(jìn)行K均值聚類(lèi)

圖7 原始序列報(bào)警示意圖Fig 7 Alarming schematic diagram of original time sequence

從圖7中可以看出，油液含水量的界限值為1×10-4， 2019年2月14日4∶22時(shí)的實(shí)測(cè)值為9.685×10-5，尚未達(dá)到報(bào)警界限值，而且此時(shí)下一步的預(yù)測(cè)值為7.83×10-5，也未達(dá)到報(bào)警界限值，即不會(huì)提前報(bào)警。只有在2019年2月14日5∶22時(shí)的實(shí)測(cè)值為1.157 2×10-4，超過(guò)報(bào)警界限值，系統(tǒng)及時(shí)報(bào)警。即在這種情況下，系統(tǒng)在測(cè)量值超過(guò)報(bào)警界限值時(shí)會(huì)及時(shí)報(bào)警，但無(wú)法提前報(bào)警。

圖8 殘差報(bào)警示意圖Fig 8 Alarming schematic diagram of residual

從圖8中可以看出，根據(jù)3sigma原則確定殘差界限值2.643×10-5，在2019年2月14日4∶22時(shí)的殘差值為2.987×10-5，超過(guò)殘差界限值，系統(tǒng)發(fā)出警報(bào)，即系統(tǒng)提前1 h發(fā)出警報(bào)。而根據(jù)K均值得到的界限值1.271×10-5，在2019年2月14日2∶22時(shí)的殘差值為1.877×10-5，超過(guò)殘差界限值，系統(tǒng)發(fā)出警報(bào)。相對(duì)于3sigma原則，使用K均值可在相同的采樣周期下大大增加提前報(bào)警的時(shí)間，而且報(bào)警時(shí)的含水量相對(duì)更小，可有效地防止設(shè)備故障發(fā)生和進(jìn)一步惡化。

4 結(jié)論

(1)利用ARMA模型對(duì)設(shè)備平穩(wěn)運(yùn)行時(shí)期的在線(xiàn)油液監(jiān)測(cè)數(shù)據(jù)進(jìn)行建模，利用該模型可及時(shí)掌握監(jiān)測(cè)特征量的變化趨勢(shì)，運(yùn)用殘差來(lái)表示設(shè)備偏離平穩(wěn)狀態(tài)的程度。

(2)針對(duì)ARMA模型對(duì)非平穩(wěn)期的預(yù)測(cè)精度較低且只能進(jìn)行一步預(yù)測(cè)的特點(diǎn)，利用模型擬合殘差判定設(shè)備處于正常期還是故障潛伏期，并通過(guò)K均值聚類(lèi)將設(shè)備發(fā)生故障前的殘差進(jìn)行二分類(lèi)，將兩類(lèi)中心點(diǎn)的均值作為正常期與故障潛伏期間的界限值。

(3)通過(guò)實(shí)際的監(jiān)測(cè)數(shù)據(jù)進(jìn)行模型驗(yàn)證和分析，結(jié)果表明，通過(guò)此方法可提前系統(tǒng)警報(bào)時(shí)間，有效控制油液監(jiān)測(cè)特征量的值并且防止設(shè)備故障發(fā)生。