劉 忠，劉 振，鄒淑云，陳 瑩，蔣 盈

(1．長(zhǎng)沙理工大學(xué) 能源與動(dòng)力工程學(xué)院，長(zhǎng)沙 410114；2．中國(guó)水利水電科學(xué)研究院，北京 100038)

空化會(huì)破壞水流正常的流動(dòng)狀態(tài)，降低水輪機(jī)的水力性能，嚴(yán)重時(shí)會(huì)導(dǎo)致過(guò)流部件損壞，造成機(jī)組嚴(yán)重振動(dòng)，產(chǎn)生噪聲，威脅水電站的安全經(jīng)濟(jì)運(yùn)行[1]。水輪機(jī)空化過(guò)程中大量微型射流或沖擊波作用于葉片和管壁，形成包含0.02～1 MHz的高頻聲發(fā)射(Acoustic Emission，AE)信號(hào)，并沿著水力、機(jī)械系統(tǒng)進(jìn)行傳播[2]。對(duì)伴隨空化產(chǎn)生的AE信號(hào)進(jìn)行有效檢測(cè)與特征提取，有望實(shí)現(xiàn)對(duì)空化狀態(tài)的監(jiān)測(cè)和辨別。國(guó)內(nèi)外學(xué)者對(duì)水輪機(jī)空化產(chǎn)生的非平穩(wěn)、非周期的AE信號(hào)進(jìn)行了諸多研究。Husin等[3]證實(shí)了在氣液兩相流動(dòng)條件下將聲發(fā)射技術(shù)應(yīng)用于氣泡的在線(xiàn)檢測(cè)是可行的。Rus等[4]在軸流式低水頭閉環(huán)試驗(yàn)臺(tái)上進(jìn)行試驗(yàn)，揭示了AE信號(hào)的葉片流道調(diào)制水平相對(duì)值與空化系數(shù)的關(guān)系。劉忠等[5-6]采用小波包方法和改進(jìn)經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)對(duì)混流式模型水輪機(jī)空化聲發(fā)射信號(hào)進(jìn)行了特性分析，揭示了各典型頻段小波包能量和各本征模態(tài)函數(shù)(Intrinsic Mode Function，IMF)特征分量隨空化系數(shù)的變化規(guī)律。但小波包在變換算法上并沒(méi)有取得明顯突破，仍然受限于小波基和分解層數(shù)的選擇，而EMD及其改進(jìn)算法(EEMD)仍然存在因極值點(diǎn)包絡(luò)估計(jì)誤差易出現(xiàn)模態(tài)混疊和端點(diǎn)效應(yīng)，參數(shù)難以選擇，重構(gòu)信號(hào)中存在殘余噪聲等問(wèn)題[7]。

變分模態(tài)分解(Variational Mode Decomposition，VMD)[8]方法通過(guò)迭代搜尋變分模型最優(yōu)解來(lái)確定每個(gè)分量的頻率中心及帶寬，具有更堅(jiān)實(shí)的理論基礎(chǔ)，因此得到了廣泛應(yīng)用。但VMD算法中分解層數(shù)和懲罰因子需要預(yù)先設(shè)定，極大地影響了分解效果。李華等[9]采用峭度最大值優(yōu)化VMD分解層數(shù)的方法對(duì)滾動(dòng)軸承早期故障進(jìn)行了特征提取，但僅考慮了分解層數(shù)優(yōu)化而將懲罰因子設(shè)為定值。Shi等[10]采用獨(dú)立優(yōu)化VMD參數(shù)的方法對(duì)風(fēng)力機(jī)進(jìn)行了狀態(tài)監(jiān)測(cè)，卻忽略了分解層數(shù)與懲罰因子的內(nèi)在聯(lián)系，易使迭代算法陷入局部最優(yōu)。為了實(shí)現(xiàn)全局尋優(yōu)，蝗蟲(chóng)優(yōu)化算法[11]、多目標(biāo)粒子群算法[12]等已被用于VMD參數(shù)組合優(yōu)化。灰狼優(yōu)化(Grey Wolf Optimization，GWO)[13]算法是一種模擬灰狼種群捕食行為機(jī)制和種群等級(jí)制度的智能元啟發(fā)式優(yōu)化算法，具有收斂性較強(qiáng)、設(shè)置參數(shù)較少和易實(shí)現(xiàn)等優(yōu)點(diǎn)，但這一方法容易陷入局部最優(yōu)解。布谷鳥(niǎo)搜索(Cuckoo Search，CS)[14]算法是模擬布谷鳥(niǎo)的寄生育雛行為與萊維飛行搜索機(jī)制相結(jié)合而形成的一種有效求解最優(yōu)化問(wèn)題的算法，具備很強(qiáng)的全局搜索能力。筆者結(jié)合灰狼優(yōu)化算法與布谷鳥(niǎo)搜索算法，對(duì)VMD的分解層數(shù)和懲罰因子參數(shù)進(jìn)行全局組合搜尋，形成了GWO和CS混合優(yōu)化VMD算法(簡(jiǎn)稱(chēng)優(yōu)化VMD算法)，并將其應(yīng)用于水輪機(jī)空化AE信號(hào)的特征提取中。

1 基本理論

1.1 VMD算法

VMD 算法實(shí)際上是一個(gè)變分問(wèn)題的求解過(guò)程。其將信號(hào)f分解為k個(gè)模態(tài)函數(shù)uk，且每個(gè)IMF具有固定的中心頻率ωk，使得每個(gè)模態(tài)函數(shù)的估計(jì)帶寬之和最小。其建立的約束變分模型為：

(1)

(2)

式中：δ(t)為沖擊函數(shù)。

為了將上述約束性變分問(wèn)題轉(zhuǎn)變?yōu)榉羌s束性變分問(wèn)題，引入懲罰因子α和拉格朗日乘子λ(t)，懲罰因子α可以在有高斯噪聲的情況下保證信號(hào)的重構(gòu)精度，而拉格朗日乘子λ(t)則保證了約束條件的嚴(yán)格性，故擴(kuò)展的拉格朗日表達(dá)式L為：

(3)

利用乘子交替方向算法不斷更新各uk,n+1、ωk,n+1和λk,n+1，求得式(3)的鞍點(diǎn)即為原問(wèn)題的最優(yōu)解。利用Parseval/Plancherel傅里葉等距變換，將式(3)變換到描述頻率ω特性的頻域求解，模態(tài)uk和中心頻率ωk可分別通過(guò)式(4)和式(5)進(jìn)行更新。

(4)

(5)

VMD算法的具體步驟如下：

(1) 進(jìn)行參數(shù)初始化，{uk,1}、{ωk,1}。

(2) 根據(jù)式(4)和式(5)在頻域內(nèi)更新uk和ωk。

1.2 灰狼優(yōu)化算法和布谷鳥(niǎo)搜索算法

1.2.1 灰狼優(yōu)化算法

在GWO算法中，為模擬灰狼種群的社會(huì)等級(jí)制度，定義狼群中的當(dāng)前最優(yōu)解為η，次最優(yōu)解為β，第三最優(yōu)解為φ，其他解均為v。

灰狼追捕包圍獵物的數(shù)學(xué)模型為：

(6)

式中：D為灰狼與獵物間的距離；l為當(dāng)前迭代次數(shù)；Xp為獵物的位置向量；X為灰狼的位置向量；A和C為協(xié)同系數(shù)向量。

其中,A和C的計(jì)算式為：

(7)

式中：r1、r2為[0,1]之間的隨機(jī)向量；a值隨迭代次數(shù)的增加從2到0之間線(xiàn)性減小。

在狩獵過(guò)程中，若η、β和φ灰狼知道獵物的位置，則灰狼群體就可以利用這三者的位置來(lái)判斷獵物所在的方位。其數(shù)學(xué)表達(dá)式為：

(8)

(9)

式中：Xη、Xβ和Xφ為η，β和φ狼的當(dāng)前位置；Dη、Dβ和Dφ為η，β和φ狼分別與獵物的距離。

1.2.2 布谷鳥(niǎo)搜索算法

CS算法的基本原理是把布谷鳥(niǎo)所寄生的鳥(niǎo)巢位置映射為算法的解，以寄生巢位置的優(yōu)劣作為算法的適應(yīng)度值，尋找鳥(niǎo)巢的位置是隨機(jī)或者類(lèi)似隨機(jī)的方式。

布谷鳥(niǎo)尋找鳥(niǎo)巢路徑和位置的更新表達(dá)式為：

xi(l+1)=xi(l)+γ⊕L(θ),i=1,2,…,n

(10)

式中：xi(l)為第i個(gè)鳥(niǎo)巢在第l次迭代后的位置；γ為步長(zhǎng)控制量；⊕表示點(diǎn)對(duì)點(diǎn)乘積；L(θ)為列維分布函數(shù)，服從L～u=l-θ(1<θ≤3 )。

1.2.3 灰狼和布谷鳥(niǎo)混合優(yōu)化算法

在GWO算法中，狼群始終根據(jù)η、β和φ狼的位置信息來(lái)更新自身位置，使其迅速向最優(yōu)解聚集，但這容易導(dǎo)致整個(gè)狼群過(guò)早地聚集于當(dāng)前最優(yōu)位置的某一鄰域內(nèi)，從而陷入局部最優(yōu)解。CS算法通過(guò)萊維飛行機(jī)制可隨機(jī)更新鳥(niǎo)巢位置，跳出當(dāng)前區(qū)域，全局搜索能力強(qiáng)。因此，筆者將GWO算法和CS算法進(jìn)行混合優(yōu)化，利用CS算法對(duì)適應(yīng)度值最好的3只灰狼和整個(gè)種群進(jìn)行位置擾動(dòng)，讓狼群在逼近獵物的過(guò)程中根據(jù)是否被發(fā)現(xiàn)的概率p隨機(jī)地更新位置，使算法具備更強(qiáng)的全局搜索能力和更快的收斂速度，適應(yīng)性更強(qiáng)。

GWO和CS混合優(yōu)化算法流程圖見(jiàn)圖1。

圖1 GWO和CS混合優(yōu)化算法流程圖Fig.1 Flow chart of GWO and CS hybrid optimization algorithm

2 GWO和CS混合優(yōu)化VMD算法

采用GWO和CS混合優(yōu)化算法對(duì)VMD算法中的分解層數(shù)K以及懲罰因子α同時(shí)進(jìn)行尋優(yōu)，具體步驟如下：

(1) 對(duì)GWO和CS混合優(yōu)化算法進(jìn)行參數(shù)設(shè)置，其中，灰狼種群規(guī)模為M，最大迭代次數(shù)為T(mén)max，優(yōu)化參數(shù)個(gè)數(shù)為d。

(2) 設(shè)定優(yōu)化過(guò)程中的適應(yīng)度函數(shù)。計(jì)算VMD分解的K個(gè)模態(tài)分量的包絡(luò)熵E1,E2,…,EK，求得K個(gè)模態(tài)分量的平均包絡(luò)熵S，同時(shí)定義各個(gè)模態(tài)包絡(luò)熵與平均包絡(luò)熵的差方和為包絡(luò)熵差異系數(shù)B，表達(dá)式為：

(11)

式中：Ei為第i個(gè)IMF分量的包絡(luò)熵；pj為a(j)的歸一化形式；a(j)為經(jīng)Hilbert解調(diào)后得到的包絡(luò)信號(hào)；N為數(shù)據(jù)長(zhǎng)度。

將各模態(tài)分量與原始信號(hào)之間的關(guān)聯(lián)程度定義為互相關(guān)系數(shù)H。

(12)

包絡(luò)熵差異系數(shù)B代表VMD分解模態(tài)的差異程度，B值越大表明分解的效果越好；互相關(guān)系數(shù)H表示分解模態(tài)與原始信號(hào)的相關(guān)性，H值越大表明分解模態(tài)與原始信號(hào)越相關(guān)，包含更多的信號(hào)信息。為了實(shí)現(xiàn)最終分解模態(tài)的差異最大化，同時(shí)保留更多的含有原始信號(hào)信息的模態(tài)分量，故定義包絡(luò)熵差異互相關(guān)系數(shù)Q為：

Q=-max(B·H)

(13)

為了獲得全局最優(yōu)值，將包絡(luò)熵差異互相關(guān)系數(shù)Q作為適應(yīng)度值，以其最小值作為尋優(yōu)目標(biāo)，利用GWO和CS混合優(yōu)化算法對(duì)VMD的參數(shù)組合進(jìn)行尋優(yōu)。

(3) 循環(huán)迭代直至最大迭代次數(shù)后結(jié)束計(jì)算，得出包絡(luò)熵差異互相關(guān)系數(shù)的最小值，其對(duì)應(yīng)的K和α即為VMD最佳分解層數(shù)Kopt和懲罰因子αopt。其算法流程見(jiàn)圖2。

圖2 GWO和CS混合優(yōu)化VMD算法流程圖Fig.2 Flow chart of GWO and CS hybrid optimized VMD

3 仿真信號(hào)分析

為了檢驗(yàn)優(yōu)化VMD算法的有效性，構(gòu)造高頻仿真信號(hào)x(t)，如式(14)所示。取采樣點(diǎn)數(shù)4 096，時(shí)域波形如圖3所示。

設(shè)置優(yōu)化VMD算法的參數(shù)M=10，Tmax=20，d=2，α的搜索范圍為[200,5 000]，K的搜索范圍為[2,10]。

圖3 仿真信號(hào)波形Fig.3 Simulation signal waveform

(14)

經(jīng)過(guò)優(yōu)化尋優(yōu)，得到適應(yīng)度函數(shù)的收斂曲線(xiàn)(見(jiàn)圖4)。由圖4可知，全局最佳的適應(yīng)度值即最小包絡(luò)熵差異互相關(guān)系數(shù)Q=-1.653 1。對(duì)應(yīng)的VMD最佳參數(shù)組合為[1 303,3]。

圖4 仿真信號(hào)迭代收斂圖Fig.4 Iterative convergence diagram of the simulation signal

將優(yōu)化VMD算法處理結(jié)果與常規(guī)VMD和EEMD算法進(jìn)行對(duì)比，其時(shí)域波形和頻譜如圖5所示。

從圖5可以看出，3種算法均分解出了仿真信號(hào)3個(gè)主要分量的頻率為60 kHz、30 kHz和8 kHz。但是，相對(duì)于原始信號(hào)分量的波形，常規(guī)VMD和EEMD算法處理得到的IMF2和IMF3分量在波形和幅值上均出現(xiàn)了失真；其對(duì)應(yīng)的頻譜圖中出現(xiàn)了模態(tài)混疊，如圖5(b)中IMF2的8 kHz成分、IMF3的30 kHz和60 kHz成分以及圖5(c)中IMF2的60 kHz成分。而優(yōu)化VMD算法得到的各IMF的時(shí)域波形和頻譜成分更接近原始信號(hào)分量。為進(jìn)一步比較3種算法的分解效果，采用常用的反映測(cè)量數(shù)據(jù)偏離真實(shí)值程度的均方差J。

(15)

(a) 優(yōu)化VMD算法

(b) 常規(guī)VMD算法

(c) EEMD算法

式中：Ji為分解得到的第i個(gè)分量的信號(hào)誤差；xi(j)為對(duì)應(yīng)分量的原始信號(hào)；Ii(j)為分解的模態(tài)分量。

Ji越小，表示分解得到的第i個(gè)分量與第i個(gè)原始信號(hào)的偏差越小，分解效果越好，反之，則越差。

上述3種算法處理得到的各分量的均方差如表1所示。由表1可知，優(yōu)化VMD算法的各均方差均小于常規(guī)VMD和EEMD算法，表明優(yōu)化VMD算法的分解效果最佳。

表1 不同分解算法的均方差

4 試驗(yàn)信號(hào)分析

所研究的混流式水輪機(jī)空化試驗(yàn)信號(hào)來(lái)自于一座處于國(guó)內(nèi)領(lǐng)先水平、綜合精度<±0.2%的閉式水輪機(jī)模型試驗(yàn)臺(tái)。該試驗(yàn)臺(tái)配備有水輪機(jī)性能測(cè)試分析系統(tǒng)，能采集并保存試驗(yàn)數(shù)據(jù)，實(shí)時(shí)監(jiān)測(cè)模型水輪機(jī)水頭和壓力的變化，測(cè)量并計(jì)算單位流量、單位轉(zhuǎn)數(shù)等衡量水輪機(jī)效率水平的性能參數(shù)，了解水輪機(jī)空化狀態(tài)對(duì)模型水輪機(jī)運(yùn)行效率的影響。通過(guò)調(diào)整試驗(yàn)臺(tái)蓄水箱水位，控制尾水管壓力，構(gòu)造偏離設(shè)計(jì)工況的試驗(yàn)工況點(diǎn)，按從無(wú)到有、從弱到強(qiáng)的順序改變水輪機(jī)的空化狀態(tài)。水輪機(jī)尾水管段為透明有機(jī)玻璃，可采用閃頻儀記錄水輪機(jī)流態(tài)的變化。為了貼近空化聲源，減少信號(hào)在過(guò)長(zhǎng)流道和不同部件的損耗，完整地采集水輪機(jī)不同空化狀態(tài)的聲發(fā)射信號(hào)，分別在模型水輪機(jī)導(dǎo)葉拐臂和轉(zhuǎn)輪下環(huán)處布置了2套SR-150N聲發(fā)射傳感器，采用PCI-9846H多功能數(shù)據(jù)采集卡進(jìn)行信號(hào)采集。所分析的聲發(fā)射信號(hào)數(shù)據(jù)截取于水輪機(jī)導(dǎo)葉拐臂上不同空化狀態(tài)下對(duì)應(yīng)的一段，每段試驗(yàn)數(shù)據(jù)點(diǎn)數(shù)為4 096(見(jiàn)圖6)。

分別將3種不同空化狀態(tài)下的水輪機(jī)聲發(fā)射信號(hào)導(dǎo)入以灰狼和布谷鳥(niǎo)組合優(yōu)化的VMD算法中，得到了不同的分解層數(shù)和懲罰因子的最佳參數(shù)組合。優(yōu)化VMD算法分解試驗(yàn)信號(hào)的初始參數(shù)設(shè)置同仿真信號(hào)一致。處理后的聲發(fā)射信號(hào)各階IMF的時(shí)域波形及其頻譜如圖7所示。

(a) 無(wú)空化信號(hào)

圖7 無(wú)空化、空化初生和臨界空化狀態(tài)下的IMF及其頻譜圖Fig.7 IMF and the spectrum under the condition without cavitation, with initial cavitation and with critical cavitation

以臨界空化狀態(tài)為例，優(yōu)化VMD算法的收斂曲線(xiàn)如圖8所示。全局最佳的包絡(luò)熵差異互相關(guān)系數(shù)Q=-1.585 3，其對(duì)應(yīng)的分解層數(shù)K=8，懲罰因子α=987。其他工況因篇幅限制，不再詳述。優(yōu)化VMD算法參數(shù)的尋優(yōu)結(jié)果如表2所示。

不同空化狀態(tài)的聲發(fā)射信號(hào)分解得到的IMF不同，對(duì)應(yīng)的各頻帶能量也不同，可作為反映空化程度的特征參數(shù)。設(shè)xi為分解得到的第i個(gè)模態(tài)分量IMF，其對(duì)應(yīng)的能量為Fi，則有：

(16)

圖8 臨界空化狀態(tài)信號(hào)迭代收斂圖Fig.8 Iterative convergence diagram of critical cavitation signals

表2 優(yōu)化VMD算法參數(shù)尋優(yōu)結(jié)果

計(jì)算各階IMF能量，保留能量值較大的前5層IMF，其隨空化系數(shù)變化的關(guān)系如圖9所示。從圖9可以看出，除IMF1外，其他IMF能量值均隨空化系數(shù)的減小逐漸增大，具有明顯的規(guī)律性。結(jié)合圖7中各IMF的頻譜圖可以看出，IMF1所處的頻率段較低，對(duì)應(yīng)的信號(hào)成分應(yīng)為水輪機(jī)運(yùn)行的背景噪聲[5]，因此與其他IMF的能量變化規(guī)律不同。在尾水真空壓力不斷降低且尚未低至飽和汽化壓力時(shí)，空化還未形成，水流流態(tài)相對(duì)簡(jiǎn)單穩(wěn)定，此時(shí)采集的信號(hào)主要是由水體流動(dòng)以及水流與水輪機(jī)的相互作用形成的。隨著壓力繼續(xù)降低至飽和汽化壓力時(shí)，空泡開(kāi)始出現(xiàn)，形成氣液兩相流動(dòng)，水流流態(tài)開(kāi)始變得復(fù)雜，此時(shí)為空化初生狀態(tài)，聲發(fā)射信號(hào)能量較無(wú)空化時(shí)增大。隨著水中空泡數(shù)量的急劇增多，潰滅沖擊程度加劇，氣液兩相流的流動(dòng)更加紊亂，此時(shí)采集的聲發(fā)射信號(hào)各頻段能量更大。因此，聲發(fā)射信號(hào)分解得到的主要IMF的能量隨空化系數(shù)變化的關(guān)系，可以敏銳地反映混流式水輪機(jī)從無(wú)空化、空化初生到臨界空化的狀態(tài)變化。

圖9 IMF能量與空化系數(shù)的關(guān)系Fig.9 IMF energy vs. cavitation coefficient

5 結(jié) 論

(1) 優(yōu)化VMD算法以包絡(luò)熵差異互相關(guān)系數(shù)為適應(yīng)度函數(shù)，自適應(yīng)地確定了影響VMD處理效果的最佳分解層數(shù)Kopt和懲罰因子αopt，在保證各模態(tài)之間差異化最大的同時(shí)，又選擇了與原始信號(hào)相關(guān)度最為緊密的模態(tài)分量，減少了人為因素造成的主觀(guān)誤差。

(2) 將優(yōu)化VMD算法用于水輪機(jī)空化聲發(fā)射信號(hào)的特征提取，得到位于空化主要特征頻率段各階IMF能量隨空化系數(shù)的減小不斷增大，直接反映了水輪機(jī)空化狀態(tài)加劇，水流更加紊亂的過(guò)程，驗(yàn)證了優(yōu)化VMD算法用于水輪機(jī)空化狀態(tài)識(shí)別的有效性。