余發(fā)軍， 劉義才

(1.中原工學(xué)院 電子信息學(xué)院，河南，鄭州 450007；2.武漢商學(xué)院 機(jī)電工程與汽車服務(wù)學(xué)院，湖北，武漢 430056)

對(duì)信號(hào)進(jìn)行分解可以觀察其內(nèi)部結(jié)構(gòu)信息，為提取其特征信息提供了一種有效途徑. 例如，傅里葉變換理論將信號(hào)分解為一系列具有不同頻率的正余弦基的線性組合，可提取信號(hào)中不同頻譜信息；小波變換理論將信號(hào)分解為一系列具有不同伸展和尺度的小波基的線性組合，可提取信號(hào)中不同成分的伸展和尺度信息；稀疏分解理論將信號(hào)分解為少量原子的線性組合，通過原子可提取當(dāng)前信號(hào)的波形成分. 與傅里葉變換和小波變換相比，稀疏分解不再保證分解基的正交性，具有更強(qiáng)的靈活性，得到了廣泛的應(yīng)用.

自1993年匹配追蹤(matching pursuit，MP)[1]算法提出以來，稀疏分解發(fā)展起以L0-范數(shù)為模型的貪婪算法[2-3]和以L1-范數(shù)為模型的基追蹤算法(basis pursuit，BP)[4]及其變體等. 前者利用信號(hào)與原子的內(nèi)積篩選出最佳原子，具有稀疏度高和重構(gòu)性強(qiáng)的優(yōu)點(diǎn). 然而，當(dāng)原子數(shù)量和信號(hào)長度較大時(shí)，內(nèi)積運(yùn)算開銷會(huì)使分解速度快速下降，影響了該方法在實(shí)際中的應(yīng)用. 利用人工智能算法提高最佳原子搜索速度是解決該問題的有效途徑[5]，如：遺傳算法(genetic algorithm，GA)[6]、粒子群優(yōu)化(particle swarm optimization，PSO)[7]、人工蜂群優(yōu)化(artificial bee colony optimization，ABCO)[8]等，這些人工智能算法可以使稀疏分解的速度得到有效提高.

量子進(jìn)化算法(quantum evolutionary algorithm，QEA)[9]，采用量子比特編碼，通過量子相位旋轉(zhuǎn)和門操作進(jìn)行染色體更新，豐富了種群個(gè)體的多樣性和密度，結(jié)合了量子比特運(yùn)算的并行性和遺傳算法優(yōu)勢，具備很強(qiáng)的搜索能力和速度. 近年來，被快速應(yīng)用到多個(gè)領(lǐng)域，如：神經(jīng)網(wǎng)絡(luò)優(yōu)化[10]，支持向量機(jī)[11]等. 然而，QEA采用二進(jìn)制的編碼方式需要一定的編解碼開銷；另外，采用固定的旋轉(zhuǎn)相位角或適應(yīng)度梯度更新量子比特相位角作為進(jìn)化策略和NOT門或其他門操作作為變異策略，在搜索優(yōu)化參數(shù)復(fù)雜情況下，缺乏一定的自適應(yīng)性.

本文提出的改進(jìn)型量子進(jìn)化算法(IQEA)，采用量子比特概率幅的編碼方式[12-13]，避免二進(jìn)制數(shù)與實(shí)際參數(shù)轉(zhuǎn)換過程中的編解碼耗時(shí)，并引入編碼因子，增強(qiáng)種群個(gè)體多樣性和密度；利用簡化形式的梯度進(jìn)化操作和逐代縮減的變異操作，在保障有效進(jìn)化和變異的前提下節(jié)約時(shí)間開銷，具有更強(qiáng)的運(yùn)算速度和自適應(yīng)性. 提出的基于IQEA的稀疏特征提取方法，利用IQEA強(qiáng)大的搜索能力和運(yùn)算速度，解決以L0-范數(shù)為模型的稀疏分解過程中內(nèi)積運(yùn)算開銷大、搜索速度慢的問題. 將所提方法應(yīng)用于仿真信號(hào)和故障軸承的振動(dòng)信號(hào)的特征提取實(shí)驗(yàn)中，結(jié)果驗(yàn)證了所提方法的有效性和優(yōu)越性.

1 改進(jìn)量子進(jìn)化算法(IQEA)

1.1 量子編碼

在量子計(jì)算中，量子比特表達(dá)為

|φ〉=α|0〉+β|1〉.

(1)

(2)

利用某個(gè)角度的余弦和正弦作為量子比特的概率幅，則該角度稱為量子比特相位角，記為θ，量子比特就表達(dá)為

|φ〉=cosθ|0〉+sinθ|1〉.

(3)

簡寫為[cosθsinθ]T.

設(shè)種群個(gè)體為

(4)

式中：m為種群規(guī)模；n為個(gè)體染色體長度；θij(i=1,2,…,m;j=1,2,…n)為第i個(gè)個(gè)體的第j個(gè)量子比特相位角，可取[0，2π)范圍內(nèi)任意值.

為了保證量子進(jìn)化算法的搜索全局收斂性和高效性，要求種群個(gè)體的多樣性和高密度性. 利用增強(qiáng)型的量子比特編碼方式[14]，將編碼因子λ(≥1)引入種群染色體中，即

(5)

定理1在量子比特相位角取值范圍相同的情況下，編碼因子λ引入種群染色體后，種群中所含最優(yōu)解的密度為引入前的λn倍，其中n為染色體長度.

證明設(shè)某優(yōu)化問題的最優(yōu)解映射到單位空間In=[-1, 1]n后為X=[x1x2…xn]，則由文獻(xiàn)[10]可知，對(duì)于?xi∈X，(i=1,2,…n)，在[0，2π)范圍內(nèi)均存在4個(gè)量子比特相位角與之對(duì)應(yīng)，即

(6)

(7)

式中θci和θsi分別為最優(yōu)解xi對(duì)應(yīng)的余弦相位角和正弦相位角.

編碼因子λ引入種群染色體中后，當(dāng)θ仍然在[0，2π)范圍內(nèi)取值時(shí)，λθ則在[0，2πλ)范圍取值. 因?yàn)棣恕?，可將[0，2πλ)分解為[0，2π)∪[2π,4π)∪…∪[2π(λ-1),2πλ)，則對(duì)?xi∈X，(i=1,2,…n)，對(duì)應(yīng)的最優(yōu)量子比特相位角為

(8)

(9)

由此可知，編碼因子λ后增加了最優(yōu)解的個(gè)數(shù)，提高了尋優(yōu)概率. 實(shí)際應(yīng)用中，可設(shè)置λ=2、θ取值范圍為[0，π)，以保證θ取值在小范圍時(shí)仍具有很強(qiáng)的尋優(yōu)性能.

1.2 交叉進(jìn)化-變異操作

考慮到生物體代與代之間即存在“小變化”的進(jìn)化，也存在“大變化”的變異幾率. 基于此在生成下一代種群時(shí)，提出交叉進(jìn)化-變異的種群個(gè)體更新操作.

首先通過旋轉(zhuǎn)量子比特相位角進(jìn)行進(jìn)化操作，即

(10)

旋轉(zhuǎn)角Δθ的方向由式(11)確定[12]：

(11)

式中θ0為當(dāng)前最優(yōu)解對(duì)應(yīng)的種群個(gè)體中量子比特相位角. 旋轉(zhuǎn)角Δθ大小采用簡化形式的梯度運(yùn)算加以確定，以減少算法的開銷，如式(12).

(12)

式中θmax和θmin分別取0.1π和0.005π.

由式(12)確定的旋轉(zhuǎn)相位角大小與當(dāng)前量子比特與最優(yōu)量子比特的相位角差值的關(guān)系如圖1所示. 可以看出：由于指數(shù)函數(shù)ex的斜率隨x的增大而增大，當(dāng)當(dāng)前量子比特與最優(yōu)量子比特的相位角接近時(shí)，由式(12)計(jì)算的Δθ值較??；反之，計(jì)算的Δθ值較大.

變異操作時(shí)，考慮種群個(gè)體的優(yōu)良性，當(dāng)其具備良好的環(huán)境適應(yīng)能力和生存能力時(shí)，變異幾率?。幌喾?，則變異幾率大. 種群個(gè)體的優(yōu)良性用η表示.

(13)

式中Fs和F分別為設(shè)定的適應(yīng)度閾值和該種群個(gè)體的當(dāng)前適應(yīng)度.

以種群個(gè)體優(yōu)良性作為發(fā)生變異的幾率基礎(chǔ)， 提出的變異操作表達(dá)為[15]

(14)

式中G和Ngen分別為設(shè)定的最大代數(shù)和當(dāng)前代數(shù).

由式(14)可以看出：當(dāng)種群個(gè)體達(dá)到“優(yōu)良”前，即η=0，得Δθ=π/2，相當(dāng)于NOT門操作，發(fā)生“大變化”的變異；而當(dāng)種群個(gè)體達(dá)到“優(yōu)良”后，其變異的幅度隨代的增加而逐漸縮減，直到減小到0.

更新下一代種群個(gè)體整體采用這種交叉進(jìn)化-變異操作，使得生成的子代具有更好的適應(yīng)度.

2 基于IQEA的稀疏特征提取

通過稀疏分解將信號(hào)表征為一組原子的線性組合，由于信號(hào)中不同成分具有不同的結(jié)構(gòu)特征，因此反映在原子上也具有不同的結(jié)構(gòu)，借此可通過原子的重新組合進(jìn)行信號(hào)不同成分的特征提取.

2.1 提取模型

設(shè)一維信號(hào)f∈RN中包含兩種不同結(jié)構(gòu)的信號(hào)成分x∈RN和y∈RN，三者滿足疊加關(guān)系：

f=x+y.

(15)

若D=[d1d2…dM]∈RN×M(M>>N,‖di‖2=1)為一個(gè)過完備字典，則據(jù)稀疏分解理論可將f在D表示成

(16)

式中：〈,〉為內(nèi)積運(yùn)算；ε為信號(hào)的殘差.

由于x和y具有不同的結(jié)構(gòu)特征，所以稀疏分解過程中x與di表現(xiàn)為強(qiáng)相關(guān)性，而與dj表現(xiàn)為弱相關(guān)性；同理y與dj為強(qiáng)相關(guān)性，而與di為弱相關(guān)性. 若要提取出信號(hào)中x成分，則可建立與其結(jié)構(gòu)相似的原子字典，當(dāng)將信號(hào)f用匹配追蹤類算法在該字典上稀疏分解時(shí)，x成分的投影系數(shù)大，故首先得到分解，而y成分的投影系數(shù)小，故后得到分解. 選取恰當(dāng)?shù)南∈璺纸饨K止條件，就可有效將x成分從信號(hào)f中提取出來[16].

2.2 提取算法

首先利用所提IQEA對(duì)原子字典進(jìn)行量子編碼. 下面以參數(shù)化Gabor字典為例說明編碼過程，原子定義為

(17)

式中：參數(shù)γ=(s,u,v,w)=(2q,p,2πα/N,πβ/6)，0≤q≤lbN，0≤p≤N-1，0≤α≤N-1，0≤β≤12，N為信號(hào)的長度，g(t)=exp(-πt2)為高斯窗函數(shù).

利用量子比特概率幅表達(dá)原子的參數(shù)(q,p,α,β)值，將其線性映射到單位空間[-1, 1]上. 若量子比特為[cosλθij， sinλθij]T，則原子參數(shù)(q,p,α,β)編碼為

(18)

編碼完成后，采用正交匹配追蹤算法(OMP)將信號(hào)f在編碼后的原子字典上稀疏分解. 將IQEA的適應(yīng)度函數(shù)fit(?)設(shè)為

fit(k)=〈Rk-1f,D〉.

(19)

式中：k為稀疏分解的次數(shù)；Rk-1f為第k-1次稀疏分解后信號(hào)的殘差，當(dāng)k=1時(shí)，R0f=f.

在每次稀疏分解前，采用隨機(jī)數(shù)初始化染色體個(gè)體的每個(gè)量子比特相位角，利用適應(yīng)度函數(shù)篩選出本代最優(yōu)個(gè)體基因位；再進(jìn)行種群個(gè)體的交叉進(jìn)化-變異操作，最后利用適應(yīng)度函數(shù)選出最優(yōu)個(gè)體基因位，直到最優(yōu)適應(yīng)度函數(shù)值滿足式(20)結(jié)束本次稀疏分解.

(20)

式中：fit*()表示最優(yōu)適應(yīng)度值.

對(duì)于整個(gè)稀疏分解結(jié)束的條件，可根據(jù)所提取特征成分的信息進(jìn)行設(shè)置. 若提取含噪信號(hào)的特征成分，則根據(jù)稀疏分解的信號(hào)殘差是否小于噪聲功率作為結(jié)束條件. 若信號(hào)噪聲功率難以估計(jì)，則可根據(jù)近兩次稀疏分解信號(hào)殘差變化率是否小于設(shè)定的閾值作為結(jié)束條件.

總結(jié)基于IQEA的稀疏特征提取算法步驟如下.

步驟1：建立恰當(dāng)?shù)脑幼值浜鸵欢ㄒ?guī)模的量子種群；

步驟2：初始化種群個(gè)體的每個(gè)量子比特相位角，并對(duì)原子字典進(jìn)行量子編碼；

步驟3：利用OMP算法將信號(hào)殘差在量子編碼后的原子字典上進(jìn)行稀疏分解，篩選最佳原子和對(duì)應(yīng)種群個(gè)體及量子比特相位；

步驟4：利用交叉進(jìn)化-變異操作進(jìn)行量子種群更新；

步驟5：重復(fù)步驟3～4，直到滿足式(20)結(jié)束本次稀疏分解；

步驟6：重復(fù)步驟2～5，直到滿足信號(hào)殘差變化率小于設(shè)定閾值，結(jié)束整體稀疏分解；

步驟7：利用每次稀疏分解時(shí)篩選出的最佳原子進(jìn)行稀疏重構(gòu)(如式(21))，即為提取的特征成分.

(21)

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

3.1 仿真信號(hào)的稀疏特征提取

已知仿真信號(hào)長度N=128，其中包含3個(gè)特征成分.

f(t)==gγ1(t)+gγ2(t)+gγ3(t).

(22)

式中：γ1=(4,33,4π,5π/6)；γ2=(8,100,2π,3π/7)；γ3=(6,64,3π,4π/9). 現(xiàn)加入不同噪聲強(qiáng)度的高斯白噪聲，分別利用所提方法、基于遺傳算法的正交匹配追蹤方法(GA-OMP)和基于雙鏈量子遺傳算法的正交匹配追蹤方法(DCQGA-OMP)，對(duì)含噪信號(hào)進(jìn)行稀疏特征提取，3種進(jìn)化方法的參數(shù)設(shè)置如表1所示.

表1 3種進(jìn)化方法參數(shù)設(shè)置Tab.1 Parameter settings of the three evolution methods

將3種方法稀疏分解的結(jié)束條件均設(shè)定為

(23)

式中：Rkf為第k稀疏分解殘差信號(hào)；Pn為噪聲功率. 提取結(jié)果的優(yōu)劣性采用均方根誤差(RMSE)衡量為

(24)

式中f(t)和f′(t)分別為原信號(hào)和提取的特征成分.

3種方法提取結(jié)果的RMSE隨信噪比SNR變化關(guān)系如圖2所示. 可以看出,本文方法在相同強(qiáng)度噪聲環(huán)境下提取還原出特征分量的失真度比GA-OMP方法和DCQGA-OMP方法更小. 這主要由于所提IQEA在匹配具有一定結(jié)構(gòu)特征的信號(hào)成分時(shí)，尋優(yōu)篩選出的Gabor原子參數(shù)(q,p,α,β)取值是連續(xù)的，使得匹配進(jìn)一步提高；另外，由于隨機(jī)噪聲不具有一定的結(jié)構(gòu)特征，與Gabor原子做內(nèi)積運(yùn)算得到的適應(yīng)度小，所以在設(shè)定一定的稀疏分解結(jié)束條件下，噪聲得不到稀疏分解，使所提方法在一定噪聲強(qiáng)度下具有較強(qiáng)的提取原信號(hào)能力.

3種方法提取特征成分所需的稀疏分解次數(shù)與SNR關(guān)系如圖3所示. 可以看出，所提方法在相同強(qiáng)度噪聲環(huán)境下稀疏提取特征成分所需的分解次數(shù)比GA-OMP方法和DCQGA-OMP方法明顯少，這進(jìn)一步說明在相同的分解結(jié)束條件下，所提方法在每次稀疏分解過程中篩選出的匹配原子比其他兩種方法更精確.

3.2 故障軸承振動(dòng)信號(hào)的稀疏特征提取

滾動(dòng)軸承是旋轉(zhuǎn)機(jī)械的關(guān)鍵部件之一，由于工況的復(fù)雜性其極易發(fā)生故障. 為避免經(jīng)濟(jì)損失和發(fā)生事故，需在軸承發(fā)生故障的早期階段及時(shí)識(shí)別并診斷出故障類型，以便及時(shí)更換. 利用安裝在滾動(dòng)軸承周圍的加速度傳感器采集其振動(dòng)信號(hào)，并提取出振動(dòng)信號(hào)中的早期故障特征，再根據(jù)故障特征的頻譜分布，可有效識(shí)別和診斷軸承常見的內(nèi)環(huán)、外環(huán)、滾動(dòng)體和保持架等故障類型. 本實(shí)驗(yàn)利用所提方法對(duì)故障軸承的振動(dòng)信號(hào)進(jìn)行稀疏特征提取，再進(jìn)行包絡(luò)譜分析確定故障類型，進(jìn)一步驗(yàn)證所提方法對(duì)實(shí)際工程信號(hào)的有效性.

實(shí)驗(yàn)平臺(tái)采用型號(hào)為QPZZ-Ⅱ旋轉(zhuǎn)機(jī)械振動(dòng)故障試驗(yàn)設(shè)備，該設(shè)備由調(diào)速電機(jī)、傳動(dòng)軸、軸承支架、模擬負(fù)載等部件組成. 用激光機(jī)將型號(hào)為N205EM軸承外環(huán)跑道內(nèi)壁加工一個(gè)凹點(diǎn)瑕疵點(diǎn)，然后將其安裝在軸承支架上，其外環(huán)固定、內(nèi)環(huán)隨傳動(dòng)軸轉(zhuǎn)動(dòng)，設(shè)置電機(jī)轉(zhuǎn)速為1 800 r/min. 振動(dòng)信號(hào)由安裝在該軸承外環(huán)徑向周圍的壓電加速度傳感器采集得到，采樣頻率為12 kHz. 任取采集的1 s時(shí)間段內(nèi)的振動(dòng)信號(hào)波形(如圖4所示)，僅從此圖無法判斷該軸承的故障類型.

采用所提方法對(duì)該段軸承振動(dòng)信號(hào)進(jìn)行故障特征提取，量子進(jìn)化參數(shù)設(shè)置如表1，稀疏分解的結(jié)束條件設(shè)置為

(25)

提取結(jié)果如圖5所示，對(duì)提取結(jié)果進(jìn)行Hilbert包絡(luò)譜分析，其頻譜分布如圖6所示. 可以看出，圖5重構(gòu)信號(hào)中沖擊成分周期性很明顯，圖6中頻譜峰值處的頻率約為145.6 Hz. 而型號(hào)為N205EM的軸承其參數(shù)：外徑為52 mm、內(nèi)徑為25 mm、滾動(dòng)體數(shù)為12、滾動(dòng)體直徑為7.5 mm、接觸角為0°. 當(dāng)轉(zhuǎn)速為1 800 r/min時(shí)，根據(jù)軸承故障機(jī)理可分別計(jì)算出該軸承常見故障的理論特征頻率如表2所示. 而圖6中峰值處的頻率恰與外環(huán)故障理論特征頻率相吻合，由此可診斷出該軸承外環(huán)存在缺陷點(diǎn)，這一診斷結(jié)論與事實(shí)相符.

Tab.2 Theoretical characteristic frequencies of N205EM bearing at rotation speed of 1 800 r/min

故障類型理論特征頻率值/Hz外環(huán)故障144.95內(nèi)環(huán)故障215.06滾動(dòng)體故障74.08保持架故障12.08

4 結(jié) 論

提出的基于IQEA的稀疏特征提取方法，將編碼因子參數(shù)引入到量子比特相位角中，增強(qiáng)尋優(yōu)能力；采用交叉進(jìn)化-變異操作更新種群個(gè)體，在保障有效更新前提下減小計(jì)算開銷；以待提取信號(hào)與量子編碼的原子字典內(nèi)積為適應(yīng)度函數(shù)，篩選出每次稀疏分解的最佳原子. 通過對(duì)含噪仿真信號(hào)的稀疏特征提取表明，所提方法在相同噪聲環(huán)境下單次稀疏分解篩選出最佳原子較其他方法更精確，提取的特征成分與原始信號(hào)更接近；對(duì)故障軸承振動(dòng)信號(hào)的稀疏特征提取表明，所提方法可有效提取故障瞬態(tài)信息，為分析實(shí)際工程信號(hào)提供了一種有效的特征提取方法.