安宇晨，鄭 陽，陳啟卷，席 慧，閆懂林，游仕豪

（武漢大學(xué)動(dòng)力與機(jī)械學(xué)院，武漢430072）

0 引 言

水電機(jī)組是電站的重要設(shè)備，其運(yùn)行的穩(wěn)定性對(duì)電站安全運(yùn)行影響重大，振動(dòng)則是影響水電機(jī)組穩(wěn)定性的重要因素。由于水電機(jī)組受到機(jī)械、水力和電氣因素的共同作用，其振動(dòng)信號(hào)具有明顯的非平穩(wěn)、非線性特點(diǎn)［1］。對(duì)于水電機(jī)組的振動(dòng)信號(hào)，通常從時(shí)域、頻域或時(shí)頻域的角度，通過幅值域分析［2］、傅里葉變換、小波變換、Hilbert-Huang 變換等方法進(jìn)行處理［3］。但是這些方法各有不足之處：傅里葉變換不能分析非平穩(wěn)信號(hào)的時(shí)域特征，小波變換雖然能夠分析信號(hào)的時(shí)頻特征，但是基函數(shù)的選擇對(duì)分析結(jié)果的影響較大，且缺乏統(tǒng)一理論依據(jù)，Hil?bert-Huang 變換雖然能夠?qū)崿F(xiàn)非平穩(wěn)信號(hào)的自適應(yīng)分析，但存在端點(diǎn)效應(yīng)和模態(tài)混疊現(xiàn)象，可能導(dǎo)致分解結(jié)果失真，甚至產(chǎn)生虛假的模態(tài)分量等情形，嚴(yán)重影響信號(hào)處理的效果，多位學(xué)者針對(duì)上述問題進(jìn)行了深入分析［2-5］。

研究表明，水電機(jī)組振動(dòng)信號(hào)呈現(xiàn)出分形特征，是典型的分形信號(hào)［6］，因此可使用多重分形去趨勢(shì)波動(dòng)分析（multifractal detrend fluctuation analysis，MFDFA）處理。該方法以分形幾何為基礎(chǔ)，研究整體與局部間自相似性與不規(guī)則性，具有不依賴系統(tǒng)模型即可實(shí)現(xiàn)系統(tǒng)狀態(tài)信息提取的特點(diǎn)［6，7］，從而避免了常見的時(shí)頻域信號(hào)分析方法所存在的不足，基于此即可實(shí)現(xiàn)信號(hào)的特征提取和分類識(shí)別。陳龍［8］分析了不同類型軸承的時(shí)間序列信號(hào)并進(jìn)行研究，證明了軸承振動(dòng)信號(hào)具有多重分形特性；林近山［9］將多重分形分析方法引入至齒輪箱的故障診斷中，證明了基于該方法提取的特征對(duì)齒輪箱狀態(tài)變化十分敏感，并根據(jù)多重分形譜特征的均值與均方差選擇特征，最終實(shí)現(xiàn)了對(duì)齒輪箱正常、輕度磨損、中度磨損和斷齒4 種狀態(tài)的識(shí)別；裴峻峰［10］基于聲信號(hào)的多重分形譜特征和改進(jìn)的SVM 實(shí)現(xiàn)了往復(fù)泵的故障診斷，提高了往復(fù)泵故障診斷準(zhǔn)確率；唐利平［11］提取刀具磨損實(shí)驗(yàn)聲發(fā)射信號(hào)的多重分形譜參數(shù)，對(duì)比不同狀態(tài)下多重分形譜特征的分布情況選擇特征，實(shí)現(xiàn)了刀具磨損狀態(tài)的識(shí)別。

目前雖已有學(xué)者從多重分形的角度實(shí)現(xiàn)了設(shè)備的故障診斷，但由于多重分形分析方法提取出的特征種類較為豐富，通常根據(jù)不同故障下各種特征數(shù)值的分布情況進(jìn)行特征的選擇［8-11］，或者直接將譜圖曲線輸入至分類器中實(shí)現(xiàn)分類［7］，因此亟須在故障特征選擇方面開展進(jìn)一步研究。

1 基于多重分形分析與BGSA-PNN 的狀態(tài)識(shí)別算法

本文提出的狀態(tài)識(shí)別算法主要由特征提取、特征選擇與分類器參數(shù)的聯(lián)合尋優(yōu)、狀態(tài)識(shí)別模型構(gòu)建三部分組成，涉及的理論和方法如下。

1.1 多重分形去趨勢(shì)波動(dòng)分析

多重分形去趨勢(shì)波動(dòng)分析方法（multifractal detrend fluctua?tion analysis，MFDFA）由Kantelhardt于2002年提出［12］，可以精確地計(jì)算出振動(dòng)信號(hào)的多重譜。多重分形譜圖的特征對(duì)振動(dòng)信號(hào)的微弱變化較為敏感，因此其可作為振動(dòng)信號(hào)的特征，用于設(shè)備的狀態(tài)監(jiān)測(cè)與故障診斷中［9，13-15］。多重分形去趨勢(shì)波動(dòng)分析的計(jì)算過程如下。

對(duì)于時(shí)間序列xk(k= 1，2，…，n)，計(jì)算其累積離差Y(i)：

將序列Y劃分成Ns個(gè)長(zhǎng)度為s的子序列，由于序列Y的長(zhǎng)度可能不是子序列長(zhǎng)度的整數(shù)倍，從反方向重復(fù)這一過程，共得到2Ns個(gè)子序列。對(duì)于每個(gè)子序列，使用最小二乘法擬合并計(jì)算其方差：

其中yv(i)為各個(gè)子序列的擬合函數(shù)。求取所有子序列方差的平均值以計(jì)算q階波動(dòng)函數(shù)：

式中：q可取任意整數(shù)。當(dāng)q取各個(gè)整數(shù)時(shí)，分析Fq(s)與s的雙對(duì)數(shù)圖，確定Fq(s)與s間的冪律關(guān)系：

式中：h(q)為廣義Hurst 指數(shù)，對(duì)于具有多重分形特征的序列，h(q)將隨q的變化而變化。對(duì)于多重分形體，標(biāo)度指數(shù)τ(q)可通過如下關(guān)系計(jì)算：

根據(jù)Legendre 變換可以計(jì)算奇異性指數(shù)α和多重分形譜f(α)：

對(duì)于多重分形體，其多重分形f(α) -α譜圖為上凸單峰函數(shù)，多重分形譜圖特征的含義如下［8，11，13，14］：多重分形譜寬度Δα=αmax-αmin反映信號(hào)多重分形特征的強(qiáng)弱，即信號(hào)波動(dòng)的劇烈程度，Δα越大，信號(hào)的波動(dòng)越劇烈，信號(hào)的能量越強(qiáng)，多重分形特征越強(qiáng)，反之則信號(hào)波動(dòng)越平均，多重分形特征越弱；Δf=f(αmin)-f(αmax)反映了振動(dòng)信號(hào)中大、小峰值所占有的比例，Δf> 0 則信號(hào)的奇異性較強(qiáng)，Δf< 0 則信號(hào)的奇異性較弱，其中f(αmin)反映信號(hào)的大波動(dòng)特征，f(αmax)反映了信號(hào)的小波動(dòng)特征；α0=α(f=fmax)描述了振動(dòng)信號(hào)的不規(guī)則性，該數(shù)值越大，則信號(hào)隨機(jī)性越強(qiáng)，信號(hào)越趨向于隨機(jī)振動(dòng)。

將多重分形譜圖的特征組成信號(hào)的特征向量輸入至分類器中，即可實(shí)現(xiàn)信號(hào)的模式識(shí)別，由于提取到的特征較多，可基于二進(jìn)制引力搜索算法對(duì)特征進(jìn)行降維。

1.2 二進(jìn)制引力搜索算法

二進(jìn)制引力搜索算法［16］（binary gravitational search algo?rithm，BGSA）由Reshedi 提出，適合解決離散和二進(jìn)制變量的優(yōu)化問題。在粒子屬性的計(jì)算方式上與常規(guī)的引力搜索算法類似，主要區(qū)別在于粒子每個(gè)維度下坐標(biāo)的取值范圍為0 或1［17，18］。該算法主要流程如下：

定義搜索空間中第i個(gè)粒子的質(zhì)量Mi(t)為：

式中：i= 1，2，…，N，N為粒子的個(gè)數(shù)為粒子在第l維下的位置坐標(biāo)，為第i個(gè)粒子在時(shí)間為t時(shí)相對(duì)于最佳和最差粒子的適應(yīng)度的占比，用于計(jì)算粒子的質(zhì)量；fiti(t)為時(shí)間為t時(shí)第i個(gè)粒子的適應(yīng)度，worst(t)和best(t)分別代表該時(shí)間下，所有粒子中最差和最佳適應(yīng)度的值。在時(shí)間為t時(shí)，第j個(gè)粒子對(duì)第i個(gè)粒子在第l維空間內(nèi)的引力為：

式中：Rij(t)為第j個(gè)粒子和第i個(gè)粒子間的距離；ε為保證分母不為0 的極小常量為引力常數(shù)，G0為初始引力值；α為衰減系數(shù)；T為最大時(shí)間即最大迭代次數(shù)。因此第i個(gè)粒子受到的合力粒子在合力的作用下將產(chǎn)生加速度迭代過程中粒子的位置按照以下公式更新：

在引力搜索算法中，粒子的位置坐標(biāo)即為待尋優(yōu)的參數(shù)信息，因此需要將待尋優(yōu)的參數(shù)轉(zhuǎn)換為粒子的位置坐標(biāo)。根據(jù)本文研究?jī)?nèi)容，將粒子的尋優(yōu)空間維數(shù)確定為17 維，每個(gè)維下的坐標(biāo)只有0 或1，其中第1～10 維對(duì)應(yīng)MFDFA 提取出的10 個(gè)特征，將坐標(biāo)為1 的特征組合成為降維后的特征集，將第11～17 維坐標(biāo)組合成的7位二進(jìn)制數(shù)轉(zhuǎn)換為十進(jìn)制后作為分類器參數(shù)。

1.3 概率神經(jīng)網(wǎng)絡(luò)

概率神經(jīng)網(wǎng)絡(luò)（Probabilistic Neural Network，PNN）由D.F.Specht 于1989年提出［19］，與BP 神經(jīng)網(wǎng)絡(luò)和徑向基神經(jīng)網(wǎng)絡(luò)（Radial Basis Function Neural Network）相比，具有網(wǎng)絡(luò)結(jié)構(gòu)簡(jiǎn)單、收斂速度快、分類準(zhǔn)確等優(yōu)點(diǎn)［20］，更加適用于分類問題，由輸入層、隱層、求和層和輸出層4層組成，每層的作用如下：

輸入層：接收輸入數(shù)據(jù)并傳至隱層，該層節(jié)點(diǎn)個(gè)數(shù)即為樣本的特征個(gè)數(shù)。

隱層：該層神經(jīng)元個(gè)數(shù)即為輸入樣本個(gè)數(shù)，該層通過連接權(quán)值與輸入層連接并接收數(shù)據(jù)，計(jì)算輸入數(shù)據(jù)與隱層中心的距離，該層輸入輸出之間的關(guān)系如下：

式中：i= 1，2，…，M，Φij為第i類神經(jīng)元的第j個(gè)神經(jīng)元對(duì)應(yīng)的隱層輸出；σ為平滑參數(shù)；d為樣本的特征維數(shù)；i為輸入樣本的類別；xij為第i類的第j個(gè)中心。

求和層：該層神經(jīng)元與隱層中對(duì)應(yīng)的神經(jīng)元連接，接收隱層數(shù)據(jù)，對(duì)隱層同一類神經(jīng)元的輸出結(jié)果求加權(quán)平均，并將計(jì)算結(jié)果輸出至輸出層。

輸出層：根據(jù)求和層計(jì)算結(jié)果，判斷該輸入所對(duì)應(yīng)的輸出。

1.4 基于多重分形分析與BGSA-PNN 的特征提取及狀態(tài)識(shí)別

基于多重分形分析與BGSA-PNN 的特征提取方法流程如圖1所示，其步驟為：

圖1 基于多重分形分析與BGSO-PNN的特征提取與識(shí)別流程Fig.1 Flowchart of state identification of vibration signal based on multifractal and BGSA-PNN

（1）特征提?。菏褂枚嘀胤中畏椒ǚ治鏊姍C(jī)振動(dòng)信號(hào)數(shù)據(jù)，繪制其多重分形譜圖并提取多重分形特征。

（2）初始化BGSA 算法粒子：運(yùn)行BGSA 算法，根據(jù)設(shè)定參數(shù)初始化粒子。

（3）解釋粒子信息：對(duì)于BGSA 算法中的某個(gè)粒子，根據(jù)粒子的信息進(jìn)行特征降維，并選擇PNN的參數(shù)。

（4）訓(xùn)練分類器：將降維后的數(shù)據(jù)按照80%和20%的比例劃分為訓(xùn)練集和測(cè)試集，使用訓(xùn)練集訓(xùn)練PNN網(wǎng)絡(luò)。

（5）輸出PNN 分類準(zhǔn)確度：將測(cè)試集數(shù)據(jù)輸入至經(jīng)過訓(xùn)練的PNN網(wǎng)絡(luò)中，輸出得到分類器分類精度，即BGSA算法中該粒子的適應(yīng)度函數(shù)值。

（6）重復(fù)步驟3至步驟5，得到每個(gè)粒子的適應(yīng)度函數(shù)值。

（7）判斷是否達(dá)到BGSA 算法停止的條件，若未達(dá)到停止條件，則更新粒子位置，并返回至步驟3，若已達(dá)到停止條件，則停止BGSO算法并結(jié)束。

2 軸承數(shù)據(jù)驗(yàn)證

軸承數(shù)據(jù)來源于美國凱斯西儲(chǔ)大學(xué)軸承數(shù)據(jù)中心網(wǎng)站，由正常、軸承內(nèi)圈故障、外圈故障、滾動(dòng)體故障4 種狀態(tài)下的數(shù)據(jù)組成。本文選擇4 種軸承狀態(tài)下的數(shù)據(jù)各50 條，基于多重分形分析與BGSA-PNN 方法提取并選擇數(shù)據(jù)的特征，并使用PNN進(jìn)行狀態(tài)識(shí)別，以驗(yàn)證本文提出的方法。

2.1 軸承信號(hào)特征提取

提取4 種狀態(tài)下的軸承信號(hào)，計(jì)算并繪制軸承信號(hào)的多重分形譜圖如圖2 所示。對(duì)比發(fā)現(xiàn)軸承信號(hào)在4 種狀態(tài)下的多重分形譜圖差異明顯，可使用多重分形譜的圖形特征進(jìn)行軸承狀態(tài)的識(shí)別。本文提取多重分形譜圖中的αmin、αmax、f(αmin)、f(αmax)、Δf=f(αmin)-f(αmax)、Δα=αmax-αmin、α0=α(fmax)、αL=α0-αmin、αR=αmax-α0以及曲線的二階矩μ作為軸承信號(hào)的MFDFA 特征，并基于此進(jìn)行特征降維工作。部分軸承信號(hào)的MFDFA特征數(shù)據(jù)見表1。

表1 部分軸承信號(hào)特征數(shù)據(jù)Tab.1 Feature data of the partial bearing signal

圖2 4種狀態(tài)下的軸承多重分形譜圖Fig.2 Multifractal spectra of bearings in four states

考慮到提取的特征量較多，構(gòu)成的特征向量維度較高，其中可能存在一些冗余特征，將增加算力的消耗，延長(zhǎng)運(yùn)算時(shí)間，甚至可能會(huì)影響分類器的識(shí)別效果，因此有必要進(jìn)行降維處理。

2.2 基于BGSA-PNN 聯(lián)合尋優(yōu)的特征降維與狀態(tài)識(shí)別

本文選擇PNN 作為分類器，使用二進(jìn)制引力搜索算法進(jìn)行特征與分類器參數(shù)的聯(lián)合優(yōu)化。以特征向量和PNN 網(wǎng)絡(luò)參數(shù)作為BGSA 算法的優(yōu)化對(duì)象，PNN 網(wǎng)絡(luò)模型作為適應(yīng)度函數(shù)，PNN 網(wǎng)絡(luò)輸出的準(zhǔn)確度作為其適應(yīng)度函數(shù)值。尋優(yōu)算法的參數(shù)為：群體規(guī)模N= 20，最大迭代次數(shù)max_it= 50。分類準(zhǔn)確度隨迭代次數(shù)的變化曲線如圖3所示。

圖3 分類準(zhǔn)確度隨迭代次數(shù)的變化曲線Fig.3 Accuracy curve of the iterations

為避免隨機(jī)性，重復(fù)進(jìn)行50 次實(shí)驗(yàn)并記錄，取實(shí)驗(yàn)結(jié)果的平均值作為PNN 網(wǎng)絡(luò)參數(shù)的最終確定值，選擇出現(xiàn)概率大于50%的特征作為最終降維結(jié)果。PNN 參數(shù)值最終確定為0.15，各特征出現(xiàn)頻次如圖4 所示，其中特征序號(hào)1～10 分別代表特征：αmin、αmax、Δα、Δf、α0、αL、αR、f(αmin)、f(αmax)、μ，根據(jù)圖4 選擇αmin、α0、αL、f(αmin)、μ作為特征降維結(jié)果，構(gòu)成特征向量。

圖4 特征出現(xiàn)頻次統(tǒng)計(jì)Fig.4 Feature appearance statistics

將降維后的特征向量和未經(jīng)過降維的MFDFA 特征輸入至同參數(shù)PNN 進(jìn)行識(shí)別。各重復(fù)50 次實(shí)驗(yàn)，統(tǒng)計(jì)其準(zhǔn)確度與耗時(shí)，結(jié)果如表2 所示。從表2 可以看出，對(duì)于軸承故障數(shù)據(jù)，降維前后的MFDFA 特征識(shí)別的平均準(zhǔn)確度均接近99%，證明了MFDFA 特征可應(yīng)用于信號(hào)的模式識(shí)別中，且經(jīng)過降維后，計(jì)算時(shí)間得到縮短。

表2 基于軸承數(shù)據(jù)的不同方法的識(shí)別準(zhǔn)確度與計(jì)算時(shí)間統(tǒng)計(jì)Tab.2 Statistics of different Recognition accuracy and calculation time of different methods of bearing data

3 目標(biāo)電站實(shí)測(cè)數(shù)據(jù)分析

實(shí)測(cè)數(shù)據(jù)采集自國內(nèi)某電站機(jī)組，該電站水輪機(jī)型號(hào)為HLA835i-LJ-305，發(fā)電機(jī)型號(hào)為SF55-30/6400，機(jī)組設(shè)計(jì)水頭69 m，最高水頭85 m，額定有功功率為55 MW，額定轉(zhuǎn)速200 r/min。分別選擇該電站機(jī)組在穩(wěn)定運(yùn)行區(qū)、小負(fù)荷工況區(qū)和渦帶運(yùn)行區(qū)3種工況下的數(shù)據(jù)各50條進(jìn)行分析。

機(jī)組在3種工況下的水導(dǎo)軸承X向擺度的波形如圖5所示，其多重分形譜圖如圖6 所示，3 種運(yùn)行工況下部分?jǐn)?shù)據(jù)的特征值如表3所示。

表3 部分水導(dǎo)軸承X向擺度的特征數(shù)據(jù)Tab.3 Feature data of the partial water guide bearing X-direction swing waveform

圖5 水導(dǎo)軸承X向擺度波形Fig.5 X-direction swing waveform of water guide bearing

圖6 3個(gè)工況下水導(dǎo)軸承X向擺度的多重分形譜圖Fig.6 Multifractal spectra of X-direction swing of water guide bearing under three working conditions

使用BGSA-PNN 進(jìn)行特征降維和PNN參數(shù)的聯(lián)合尋優(yōu)，其中二進(jìn)制引力搜索算法的參數(shù)設(shè)置為：群體規(guī)模N= 50，最大迭代次數(shù)max_it= 50，優(yōu)化算法中的適應(yīng)度函數(shù)值即分類準(zhǔn)確度隨迭代次數(shù)變化的曲線如圖7 所示。重復(fù)實(shí)驗(yàn)50 次，根據(jù)實(shí)驗(yàn)結(jié)果確定PNN 參數(shù)值為0.013，繪制各個(gè)特征的出現(xiàn)頻次如圖8 所示，由于αmin、Δα、αL、f(αmin)、μ五維特征出現(xiàn)的頻次明顯高于其他特征，故將其選擇為分類器特征降維后的特征集。

圖7 電站實(shí)測(cè)數(shù)據(jù)的分類準(zhǔn)確度隨迭代次數(shù)的變化曲線Fig.7 Classification accuracy curve of the measured data of the power station with iterations

圖8 電站實(shí)測(cè)數(shù)據(jù)特征出現(xiàn)頻次統(tǒng)計(jì)Fig.8 Feature appearance statistics of the measured data of the power station

經(jīng)過降維，水電機(jī)組振動(dòng)信號(hào)的MFDFA 特征維數(shù)從10 降至5。使用經(jīng)過降維的MFDFA 特征進(jìn)行一次識(shí)別的結(jié)果如圖9所示，其中類別1 表示小負(fù)荷運(yùn)行工況、類別2 表示渦帶工況、類別3表示穩(wěn)定運(yùn)行工況，由圖9可以看出，基于BGSA-PNN 降維的MFDFA 特征的識(shí)別結(jié)果與實(shí)際情況十分接近，幾乎能夠完全識(shí)別渦帶工況和穩(wěn)定運(yùn)行工況，對(duì)于小負(fù)荷運(yùn)行工況則能夠完全識(shí)別。

圖9 電站實(shí)測(cè)數(shù)據(jù)特征降維后的分類結(jié)果Fig.9 Classification results after dimension reduction of measured data features of power plants

作為對(duì)比，使用降維前后的MFDFA 特征和機(jī)組實(shí)測(cè)信號(hào)的EMD 模糊熵作為特征向量，輸入至同參數(shù)的PNN 進(jìn)行識(shí)別分類，重復(fù)進(jìn)行50 次實(shí)驗(yàn)并將其統(tǒng)計(jì)結(jié)果記錄于表4。對(duì)比可見，水電機(jī)組振動(dòng)信號(hào)的MFDFA 特征在識(shí)別準(zhǔn)確度與時(shí)間上優(yōu)于EMD 模糊熵方法，且經(jīng)過降維后，神經(jīng)網(wǎng)絡(luò)分類器對(duì)于機(jī)組狀態(tài)的識(shí)別效率得到了提升，對(duì)于水電機(jī)組的實(shí)測(cè)振動(dòng)信號(hào)識(shí)別效果較好。

表4 基于實(shí)測(cè)數(shù)據(jù)的不同方法的識(shí)別準(zhǔn)確度與計(jì)算時(shí)間統(tǒng)計(jì)Tab.4 Statistics of Recognition accuracy and calculation time of different methods of the measured data

4 結(jié) 論

本文提出了一種基于多重分形分析與BSGA-PNN 的水電機(jī)組信號(hào)特征提取方法。其特點(diǎn)在于通過多重分形分析繪制水電機(jī)組振動(dòng)信號(hào)的多重分形譜圖并提取其特征，并采用BG?SA-PNN 進(jìn)行聯(lián)合尋優(yōu)，對(duì)PNN 參數(shù)進(jìn)行尋優(yōu)的同時(shí)開展多維特征向量降維操作。根據(jù)4種故障狀態(tài)下的軸承數(shù)據(jù)和三種運(yùn)行工況下的水電機(jī)組實(shí)測(cè)數(shù)據(jù)進(jìn)行對(duì)比驗(yàn)證，證明了MFDFA特征能夠有效表征振動(dòng)信號(hào)中包含的信息，且經(jīng)過BGSA-PNN聯(lián)合尋優(yōu)后，識(shí)別效率得到了提高。經(jīng)與基于EMD-PNN 的傳統(tǒng)狀態(tài)識(shí)別方法進(jìn)行對(duì)比，發(fā)現(xiàn)本文方法在識(shí)別準(zhǔn)確度和計(jì)算時(shí)間上均優(yōu)于傳統(tǒng)方法。本文方法有望推廣應(yīng)用于更多類型的水電機(jī)組振動(dòng)信號(hào)狀態(tài)識(shí)別中，提高電站安全、穩(wěn)定運(yùn)行的水平。 □