基于改进量子进化算法的稀疏特征提取方法

余发军， 刘义才

(1.中原工学院 电子信息学院，河南，郑州 450007；2.武汉商学院 机电工程与汽车服务学院，湖北，武汉 430056)

对信号进行分解可以观察其内部结构信息，为提取其特征信息提供了一种有效途径. 例如，傅里叶变换理论将信号分解为一系列具有不同频率的正余弦基的线性组合，可提取信号中不同频谱信息；小波变换理论将信号分解为一系列具有不同伸展和尺度的小波基的线性组合，可提取信号中不同成分的伸展和尺度信息；稀疏分解理论将信号分解为少量原子的线性组合，通过原子可提取当前信号的波形成分. 与傅里叶变换和小波变换相比，稀疏分解不再保证分解基的正交性，具有更强的灵活性，得到了广泛的应用.

自1993年匹配追踪(matching pursuit，MP)[1]算法提出以来，稀疏分解发展起以L0-范数为模型的贪婪算法[2-3]和以L1-范数为模型的基追踪算法(basis pursuit，BP)[4]及其变体等. 前者利用信号与原子的内积筛选出最佳原子，具有稀疏度高和重构性强的优点. 然而，当原子数量和信号长度较大时，内积运算开销会使分解速度快速下降，影响了该方法在实际中的应用. 利用人工智能算法提高最佳原子搜索速度是解决该问题的有效途径[5]，如：遗传算法(genetic algorithm，GA)[6]、粒子群优化(particle swarm optimization，PSO)[7]、人工蜂群优化(artificial bee colony optimization，ABCO)[8]等，这些人工智能算法可以使稀疏分解的速度得到有效提高.

量子进化算法(quantum evolutionary algorithm，QEA)[9]，采用量子比特编码，通过量子相位旋转和门操作进行染色体更新，丰富了种群个体的多样性和密度，结合了量子比特运算的并行性和遗传算法优势，具备很强的搜索能力和速度. 近年来，被快速应用到多个领域，如：神经网络优化[10]，支持向量机[11]等. 然而，QEA采用二进制的编码方式需要一定的编解码开销；另外，采用固定的旋转相位角或适应度梯度更新量子比特相位角作为进化策略和NOT门或其他门操作作为变异策略，在搜索优化参数复杂情况下，缺乏一定的自适应性.

本文提出的改进型量子进化算法(IQEA)，采用量子比特概率幅的编码方式[12-13]，避免二进制数与实际参数转换过程中的编解码耗时，并引入编码因子，增强种群个体多样性和密度；利用简化形式的梯度进化操作和逐代缩减的变异操作，在保障有效进化和变异的前提下节约时间开销，具有更强的运算速度和自适应性. 提出的基于IQEA的稀疏特征提取方法，利用IQEA强大的搜索能力和运算速度，解决以L0-范数为模型的稀疏分解过程中内积运算开销大、搜索速度慢的问题. 将所提方法应用于仿真信号和故障轴承的振动信号的特征提取实验中，结果验证了所提方法的有效性和优越性.

1 改进量子进化算法(IQEA)

1.1 量子编码

在量子计算中，量子比特表达为

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

(1)

(2)

利用某个角度的余弦和正弦作为量子比特的概率幅，则该角度称为量子比特相位角，记为θ，量子比特就表达为

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

(3)

简写为[cosθsinθ]T.

设种群个体为

(4)

式中：m为种群规模；n为个体染色体长度；θij(i=1,2,…,m;j=1,2,…n)为第i个个体的第j个量子比特相位角，可取[0，2π)范围内任意值.

为了保证量子进化算法的搜索全局收敛性和高效性，要求种群个体的多样性和高密度性. 利用增强型的量子比特编码方式[14]，将编码因子λ(≥1)引入种群染色体中，即

(5)

定理1在量子比特相位角取值范围相同的情况下，编码因子λ引入种群染色体后，种群中所含最优解的密度为引入前的λn倍，其中n为染色体长度.

证明设某优化问题的最优解映射到单位空间In=[-1, 1]n后为X=[x1x2…xn]，则由文献[10]可知，对于∀xi∈X，(i=1,2,…n)，在[0，2π)范围内均存在4个量子比特相位角与之对应，即

(6)

(7)

式中θci和θsi分别为最优解xi对应的余弦相位角和正弦相位角.

编码因子λ引入种群染色体中后，当θ仍然在[0，2π)范围内取值时，λθ则在[0，2πλ)范围取值. 因为λ≥1，可将[0，2πλ)分解为[0，2π)∪[2π,4π)∪…∪[2π(λ-1),2πλ)，则对∀xi∈X，(i=1,2,…n)，对应的最优量子比特相位角为

(8)

(9)

由此可知，编码因子λ后增加了最优解的个数，提高了寻优概率. 实际应用中，可设置λ=2、θ取值范围为[0，π)，以保证θ取值在小范围时仍具有很强的寻优性能.

1.2 交叉进化-变异操作

考虑到生物体代与代之间即存在“小变化”的进化，也存在“大变化”的变异几率. 基于此在生成下一代种群时，提出交叉进化-变异的种群个体更新操作.

首先通过旋转量子比特相位角进行进化操作，即

(10)

旋转角Δθ的方向由式(11)确定[12]：

(11)

式中θ0为当前最优解对应的种群个体中量子比特相位角. 旋转角Δθ大小采用简化形式的梯度运算加以确定，以减少算法的开销，如式(12).

(12)

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

由式(12)确定的旋转相位角大小与当前量子比特与最优量子比特的相位角差值的关系如图1所示. 可以看出：由于指数函数ex的斜率随x的增大而增大，当当前量子比特与最优量子比特的相位角接近时，由式(12)计算的Δθ值较小；反之，计算的Δθ值较大.

变异操作时，考虑种群个体的优良性，当其具备良好的环境适应能力和生存能力时，变异几率小；相反，则变异几率大. 种群个体的优良性用η表示.

(13)

式中Fs和F分别为设定的适应度阈值和该种群个体的当前适应度.

以种群个体优良性作为发生变异的几率基础， 提出的变异操作表达为[15]

(14)

式中G和Ngen分别为设定的最大代数和当前代数.

由式(14)可以看出：当种群个体达到“优良”前，即η=0，得Δθ=π/2，相当于NOT门操作，发生“大变化”的变异；而当种群个体达到“优良”后，其变异的幅度随代的增加而逐渐缩减，直到减小到0.

更新下一代种群个体整体采用这种交叉进化-变异操作，使得生成的子代具有更好的适应度.

2 基于IQEA的稀疏特征提取

通过稀疏分解将信号表征为一组原子的线性组合，由于信号中不同成分具有不同的结构特征，因此反映在原子上也具有不同的结构，借此可通过原子的重新组合进行信号不同成分的特征提取.

2.1 提取模型

设一维信号f∈RN中包含两种不同结构的信号成分x∈RN和y∈RN，三者满足叠加关系：

f=x+y.

(15)

若D=[d1d2…dM]∈RN×M(M>>N,‖di‖2=1)为一个过完备字典，则据稀疏分解理论可将f在D表示成

(16)

式中：〈,〉为内积运算；ε为信号的残差.

由于x和y具有不同的结构特征，所以稀疏分解过程中x与di表现为强相关性，而与dj表现为弱相关性；同理y与dj为强相关性，而与di为弱相关性. 若要提取出信号中x成分，则可建立与其结构相似的原子字典，当将信号f用匹配追踪类算法在该字典上稀疏分解时，x成分的投影系数大，故首先得到分解，而y成分的投影系数小，故后得到分解. 选取恰当的稀疏分解终止条件，就可有效将x成分从信号f中提取出来[16].

2.2 提取算法

首先利用所提IQEA对原子字典进行量子编码. 下面以参数化Gabor字典为例说明编码过程，原子定义为

(17)

式中：参数γ=(s,u,v,w)=(2q,p,2πα/N,πβ/6)，0≤q≤lbN，0≤p≤N-1，0≤α≤N-1，0≤β≤12，N为信号的长度，g(t)=exp(-πt2)为高斯窗函数.

利用量子比特概率幅表达原子的参数(q,p,α,β)值，将其线性映射到单位空间[-1, 1]上. 若量子比特为[cosλθij， sinλθij]T，则原子参数(q,p,α,β)编码为

(18)

编码完成后，采用正交匹配追踪算法(OMP)将信号f在编码后的原子字典上稀疏分解. 将IQEA的适应度函数fit(•)设为

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

(19)

式中：k为稀疏分解的次数；Rk-1f为第k-1次稀疏分解后信号的残差，当k=1时，R0f=f.

在每次稀疏分解前，采用随机数初始化染色体个体的每个量子比特相位角，利用适应度函数筛选出本代最优个体基因位；再进行种群个体的交叉进化-变异操作，最后利用适应度函数选出最优个体基因位，直到最优适应度函数值满足式(20)结束本次稀疏分解.

(20)

式中：fit*()表示最优适应度值.

对于整个稀疏分解结束的条件，可根据所提取特征成分的信息进行设置. 若提取含噪信号的特征成分，则根据稀疏分解的信号残差是否小于噪声功率作为结束条件. 若信号噪声功率难以估计，则可根据近两次稀疏分解信号残差变化率是否小于设定的阈值作为结束条件.

总结基于IQEA的稀疏特征提取算法步骤如下.

步骤1：建立恰当的原子字典和一定规模的量子种群；

步骤2：初始化种群个体的每个量子比特相位角，并对原子字典进行量子编码；

步骤3：利用OMP算法将信号残差在量子编码后的原子字典上进行稀疏分解，筛选最佳原子和对应种群个体及量子比特相位；

步骤4：利用交叉进化-变异操作进行量子种群更新；

步骤5：重复步骤3～4，直到满足式(20)结束本次稀疏分解；

步骤6：重复步骤2～5，直到满足信号残差变化率小于设定阈值，结束整体稀疏分解；

步骤7：利用每次稀疏分解时筛选出的最佳原子进行稀疏重构(如式(21))，即为提取的特征成分.

(21)

3 实验与分析

3.1 仿真信号的稀疏特征提取

已知仿真信号长度N=128，其中包含3个特征成分.

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). 现加入不同噪声强度的高斯白噪声，分别利用所提方法、基于遗传算法的正交匹配追踪方法(GA-OMP)和基于双链量子遗传算法的正交匹配追踪方法(DCQGA-OMP)，对含噪信号进行稀疏特征提取，3种进化方法的参数设置如表1所示.

表1 3种进化方法参数设置Tab.1 Parameter settings of the three evolution methods

将3种方法稀疏分解的结束条件均设定为

(23)

式中：Rkf为第k稀疏分解残差信号；Pn为噪声功率. 提取结果的优劣性采用均方根误差(RMSE)衡量为

(24)

式中f(t)和f′(t)分别为原信号和提取的特征成分.

3种方法提取结果的RMSE随信噪比SNR变化关系如图2所示. 可以看出,本文方法在相同强度噪声环境下提取还原出特征分量的失真度比GA-OMP方法和DCQGA-OMP方法更小. 这主要由于所提IQEA在匹配具有一定结构特征的信号成分时，寻优筛选出的Gabor原子参数(q,p,α,β)取值是连续的，使得匹配进一步提高；另外，由于随机噪声不具有一定的结构特征，与Gabor原子做内积运算得到的适应度小，所以在设定一定的稀疏分解结束条件下，噪声得不到稀疏分解，使所提方法在一定噪声强度下具有较强的提取原信号能力.

3种方法提取特征成分所需的稀疏分解次数与SNR关系如图3所示. 可以看出，所提方法在相同强度噪声环境下稀疏提取特征成分所需的分解次数比GA-OMP方法和DCQGA-OMP方法明显少，这进一步说明在相同的分解结束条件下，所提方法在每次稀疏分解过程中筛选出的匹配原子比其他两种方法更精确.

3.2 故障轴承振动信号的稀疏特征提取

滚动轴承是旋转机械的关键部件之一，由于工况的复杂性其极易发生故障. 为避免经济损失和发生事故，需在轴承发生故障的早期阶段及时识别并诊断出故障类型，以便及时更换. 利用安装在滚动轴承周围的加速度传感器采集其振动信号，并提取出振动信号中的早期故障特征，再根据故障特征的频谱分布，可有效识别和诊断轴承常见的内环、外环、滚动体和保持架等故障类型. 本实验利用所提方法对故障轴承的振动信号进行稀疏特征提取，再进行包络谱分析确定故障类型，进一步验证所提方法对实际工程信号的有效性.

实验平台采用型号为QPZZ-Ⅱ旋转机械振动故障试验设备，该设备由调速电机、传动轴、轴承支架、模拟负载等部件组成. 用激光机将型号为N205EM轴承外环跑道内壁加工一个凹点瑕疵点，然后将其安装在轴承支架上，其外环固定、内环随传动轴转动，设置电机转速为1 800 r/min. 振动信号由安装在该轴承外环径向周围的压电加速度传感器采集得到，采样频率为12 kHz. 任取采集的1 s时间段内的振动信号波形(如图4所示)，仅从此图无法判断该轴承的故障类型.

采用所提方法对该段轴承振动信号进行故障特征提取，量子进化参数设置如表1，稀疏分解的结束条件设置为

(25)

提取结果如图5所示，对提取结果进行Hilbert包络谱分析，其频谱分布如图6所示. 可以看出，图5重构信号中冲击成分周期性很明显，图6中频谱峰值处的频率约为145.6 Hz. 而型号为N205EM的轴承其参数：外径为52 mm、内径为25 mm、滚动体数为12、滚动体直径为7.5 mm、接触角为0°. 当转速为1 800 r/min时，根据轴承故障机理可分别计算出该轴承常见故障的理论特征频率如表2所示. 而图6中峰值处的频率恰与外环故障理论特征频率相吻合，由此可诊断出该轴承外环存在缺陷点，这一诊断结论与事实相符.

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

故障类型理论特征频率值/Hz外环故障144.95内环故障215.06滚动体故障74.08保持架故障12.08

4 结 论

提出的基于IQEA的稀疏特征提取方法，将编码因子参数引入到量子比特相位角中，增强寻优能力；采用交叉进化-变异操作更新种群个体，在保障有效更新前提下减小计算开销；以待提取信号与量子编码的原子字典内积为适应度函数，筛选出每次稀疏分解的最佳原子. 通过对含噪仿真信号的稀疏特征提取表明，所提方法在相同噪声环境下单次稀疏分解筛选出最佳原子较其他方法更精确，提取的特征成分与原始信号更接近；对故障轴承振动信号的稀疏特征提取表明，所提方法可有效提取故障瞬态信息，为分析实际工程信号提供了一种有效的特征提取方法.