余發(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)分解為少量原子的線性組合，通過(guò)原子可提取當(dāng)前信號(hào)的波形成分. 與傅里葉變換和小波變換相比，稀疏分解不再保證分解基的正交性，具有更強(qiáng)的靈活性，得到了廣泛的應(yīng)用.

自1993年匹配追蹤(matching pursuit，MP)[1]算法提出以來(lái)，稀疏分解發(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)長(zhǎng)度較大時(shí)，內(nèi)積運(yùn)算開(kāi)銷會(huì)使分解速度快速下降，影響了該方法在實(shí)際中的應(yīng)用. 利用人工智能算法提高最佳原子搜索速度是解決該問(wèn)題的有效途徑[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]，采用量子比特編碼，通過(guò)量子相位旋轉(zhuǎn)和門操作進(jìn)行染色體更新，豐富了種群個(gè)體的多樣性和密度，結(jié)合了量子比特運(yùn)算的并行性和遺傳算法優(yōu)勢(shì)，具備很強(qiáng)的搜索能力和速度. 近年來(lái)，被快速應(yīng)用到多個(gè)領(lǐng)域，如：神經(jīng)網(wǎng)絡(luò)優(yōu)化[10]，支持向量機(jī)[11]等. 然而，QEA采用二進(jìn)制的編碼方式需要一定的編解碼開(kāi)銷；另外，采用固定的旋轉(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)換過(guò)程中的編解碼耗時(shí)，并引入編碼因子，增強(qiáng)種群個(gè)體多樣性和密度；利用簡(jiǎn)化形式的梯度進(jìn)化操作和逐代縮減的變異操作，在保障有效進(jìn)化和變異的前提下節(jié)約時(shí)間開(kāi)銷，具有更強(qiáng)的運(yùn)算速度和自適應(yīng)性. 提出的基于IQEA的稀疏特征提取方法，利用IQEA強(qiáng)大的搜索能力和運(yùn)算速度，解決以L0-范數(shù)為模型的稀疏分解過(guò)程中內(nèi)積運(yùn)算開(kāi)銷大、搜索速度慢的問(wè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)

簡(jiǎn)寫(xiě)為[cosθsinθ]T.

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

(4)

式中：m為種群規(guī)模；n為個(gè)體染色體長(zhǎng)度；θ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為染色體長(zhǎng)度.

證明設(shè)某優(yōu)化問(wèn)題的最優(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è)體更新操作.

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

(10)

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

(11)

式中θ0為當(dāng)前最優(yōu)解對(duì)應(yīng)的種群個(gè)體中量子比特相位角. 旋轉(zhuǎn)角Δθ大小采用簡(jiǎn)化形式的梯度運(yùn)算加以確定，以減少算法的開(kāi)銷，如式(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的稀疏特征提取

通過(guò)稀疏分解將信號(hào)表征為一組原子的線性組合，由于信號(hào)中不同成分具有不同的結(jié)構(gòu)特征，因此反映在原子上也具有不同的結(jié)構(gòu)，借此可通過(guò)原子的重新組合進(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è)過(guò)完備字典，則據(jù)稀疏分解理論可將f在D表示成

(16)

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

由于x和y具有不同的結(jié)構(gòu)特征，所以稀疏分解過(guò)程中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中提取出來(lái)[16].

2.2 提取算法

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

(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)的長(zhǎng)度，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)長(zhǎng)度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)一步說(shuō)明在相同的分解結(jié)束條件下，所提方法在每次稀疏分解過(guò)程中篩選出的匹配原子比其他兩種方法更精確.

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í)別和診斷軸承常見(jiàn)的內(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所示)，僅從此圖無(wú)法判斷該軸承的故障類型.

采用所提方法對(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ì)算出該軸承常見(jiàn)故障的理論特征頻率如表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ì)算開(kāi)銷；以待提取信號(hào)與量子編碼的原子字典內(nèi)積為適應(yīng)度函數(shù)，篩選出每次稀疏分解的最佳原子. 通過(guò)對(duì)含噪仿真信號(hào)的稀疏特征提取表明，所提方法在相同噪聲環(huán)境下單次稀疏分解篩選出最佳原子較其他方法更精確，提取的特征成分與原始信號(hào)更接近；對(duì)故障軸承振動(dòng)信號(hào)的稀疏特征提取表明，所提方法可有效提取故障瞬態(tài)信息，為分析實(shí)際工程信號(hào)提供了一種有效的特征提取方法.