振动信号包络线的稀疏重构最优化算法研究与应用

于岩君， 叶庆卫， 陆志华， 周 宇

(宁波大学 信息科学与工程学院，浙江 宁波 315211)

HHT(Hilbert-Huang Transform)技术是由Huang等[1-2]提出的处理非线性、非平稳信号的新方法，这是继传统傅里叶时频分析和小波分析方法后的具有重大突破性的应用数学方法。HHT需要经过经验模态分解(Empirical Mode Decomposition，EMD)和Hilbert变换两大步骤实现，其中EMD分解是其核心步骤。如今，EMD算法已在机械故障诊断、信号处理与分析、经济金融工程等领域[3-5]得到了广泛的应用，而包络线的提取好坏直接影响EMD分解的结果，因此越来越多的学者对此问题做了大量的研究。钟佑明[6]在经典三次样条插值理论下提出Akima插值算法和分段幂函数法，通过三者之间的对比实验，从拟合曲线的“光滑性”和“柔性”方面做了改进，但是仍无法消除端点飞翼和“过冲”、“欠冲”问题。广义检波滤波是除三次样条插值法外的提取包络的另一经典方法，此方法容易产生解调混频效应，张帆等[7]提出采用合适的采样频率或提高采样频率可使广义检波解调分析中的混频效应忽略不计，但是这种方法对于处理复杂信号效果不佳。张郁山等[8]提出了基于自回归模型的延拓结构,采用“边筛分边延拓”方法将延拓得到的极值点和信号本身的极值点联合起来构成包络，最后筛分出IMF分量，但是自回归模型需为预测的数据点个数设置一个较大的数值,影响筛选过程的精度,通用性比较差。

文献[9]提出了基于稀疏复原的包络线提取算法。稀疏复原本质是求解一个范数最小模型问题，它可以利用较少的观测值精确地恢复出原始信号，某些信号的小波变换、离散余弦变换等都具有一定的稀疏性[10]，所以稀疏重构具有普遍适用性。虽然文献[9]的稀疏复原方法在一定程度上抑制了端点效应，但是抗噪性能有待提高，而且在寻优改变DCT基频带宽度的变化因子m时，一方面需要人为挑选，降低了包络线提取精度；另一方面寻找到的m只具有局部最优的性质，不具有全局性，所以寻优结果存在一定的片面性。

针对以上问题，本文结合粒子群算法(Particle Swarm Optimization，PSO)的独特优势[11]进一步对基于稀疏复原提取信号包络线的算法做了改进。在文献[9]的基础上引入混合变异粒子群优化算法对变化因子m进行全局寻优，从而自适应地提取出最优的包络线；并且将包络线稀疏重构优化模型中的凹二次规划问题转换为凸二次规划问题，利用内点法求得该模型中的稀疏解，其中，采用混合变异粒子群有效地避免了变化因子m陷入局部最优。实验仿真与实际提取的斜拉索振动信号测试表明，基于PSO的信号包络线稀疏重构最优化算法同样可以克服端点飞翼问题，同时PSO的引入，可以自适应的找出最适合包络线平稳特性的DCT基，克服了人为主观确定的缺点，提高了包络线的提取精度与抗噪性能。

1 稀疏复原拟合包络线基本原理

包络线的提取贯穿于整个EMD分解过程，后续也会影响到瞬时频率的计算结果。但是至今在包络线提取算法上还没有系统的理论支撑，目前最常用的经验算法有三次样条插值、分段幂函数和广义检波滤波等。文献[9]已经对三次样条插值和包络检波原理做了详细介绍，在此不再赘述。

文献[9]基于压缩感知理论，提出了利用稀疏重构的思想提取包络线，下面以上包络线为例，详细说明此方法的基本原理。

用向量x(t)=[x1(t),x2(t),…,xN(t)]T表示所要求解的上包络线信号，其中N表示采样点的个数，用yd(t)=[xd1,xd2,…,xdk]表示从原始信号x(t)中提取的所有极大值点所构成的信号，其中K表示极大值点的个数。由于x(t)不是严格的K-稀疏，所以要将其投影在构造的基矩阵上进行稀疏表示，即：

x(t)=Ψs

(1)

s为稀疏基Ψ下的K稀疏信号，Ψ为N×N的余弦变换矩阵。稀疏基的构建是整个稀疏复原过程中的关键技术之一，根据包络线的低频特性和离散余弦变换的“能量集中”特性，经过离散余弦变换后信号的能量将集中在低频部分，所以选择了结构简单的离散余弦变换基(DCT基)。为了提高包络线的重构精度与效率，文献[9]定义变化因子m，主要用于改变稀疏基的频带宽度，减小重构信号的频率之差。其表达式如下

k=1,2,…,Nn=0,1,…,N-1

(2)

由于包络线信号x(t)都经过其极值点，所以将极大值点所组成的信号yd(t)作为信号x(t)压缩采样后得到的M维观测信号，其过程可描述为

yd(t)=Φ·x(t)

(3)

矩阵ΦM×N称为观测矩阵(又称测量矩阵)，其中M在数值上与极大值点的个数K是相等的。由于所求稀疏解s稀疏度越高恢复出的x(t)效果越好，所以由式(1)、(3)可得包络信号x(t)重构模型为

.t.yd(t)=Φ·Ψ(m)·s

(4)

基于信号包络线的低频特性，应该使恢复出的包络线足够光滑，从而根据频域特征提取更精确的信息，基于此设计了一阶微分矩阵D(N-1)，定义信号的平稳度函数E(t)来体现信号包络的变化程度。其中：

(xN-1-xN)2

(5)

因此提取的最佳信号包络线是满足E(t)最小的稀疏重构获得的曲线，综上所述，基于稀疏复原的信号包络线提取的数学模型为

(6)

2 基于粒子群的包络线稀疏复原算法

2.1 稀疏复原拟合包络线的改进模型

针对以上问题，本文对此做了两点改进。首先将式(6)中的非凸模型转化为凸二次规划问题，从而利用内点法求得稀疏解；其次为了得到全局最优m值，本文将DCT基中的m放在分子上，范围取[0,1]之间，利用混合变异粒子群算法对m自主寻优，避免了m陷入局部最优的问题。同时在寻优过程中，通过改变粒子群算法中的适应度函数，直接提取出误差最小的包络线信号。

针对式(6)中的l0范数最小化问题，匹配追踪类算法是目前比较成熟的求解算法，包括MP、OMP、ROMP等。Donoho在文献[12]提出满足RIP条件下可将l0模型转化为最小化l1模型，针对l1模型一般用凸松弛算法求解，包括BP、内点法等。因此式(6)中优化模型可改为

(7)

式(7)中的λ称为正则化参数，此模型所求解一般是稀疏的(其证明详见文献[14])，由此符合稀疏复原的要求。由于上述l1-正则化模型是不可微的，所以受文献[13]的启示，通过引入线性不等式约束将稀疏重构模型转化为一个凸二次规划问题，这样便可以用内点法求解此凸二次问题，重构模型修正为

(8)

若采用内点法求解上述模型中的稀疏解s，首先要构造惩罚函数，由此将不等式约束的优化问题转换为无约束优化问题，其中惩罚函数构造形式如下

(9)

然后利用惩罚函数构造增广目标函数，即无约束优化问题的目标函数

(10)

对于求解式(10)无约束优化问题本文选择了具有超线性收敛速度且结构简单的共轭梯度法，在每次迭代过程中根据凸二次函数f(s,u)的Hesse阵G产生的共轭方向以及由精确线搜索方法所确定的搜索步长来求解正定二次目标函数极小点，此极小点即为模型中的稀疏解s。

2.2 粒子群寻优变化因子m

本文利用粒子群算法对上述模型中的Ψ(m)进行寻优，将式(8)中的目标函数作为粒子群算法中的适应度函数，同时PSO的引入可以灵活改变稀疏复原方法，本文选择了内点法求得稀疏解s， 最后将寻找到的最优稀疏基Ψ(m)和与之对应的稀疏解s代入x(t)=Ψ(m)·s中即可求得误差最小的包络信号。

用混合变异粒子群算法寻优变化因子m，能使搜索算法很快地收敛到全局最优值。由于只对m寻优，则搜索空间大小为1×1，且m的维数为一维信号。设种群粒子数为N,第i(i=1,2,…,N)个粒子的位置为xi，飞行速度为vi，第i个粒子目前为止搜索到的最优位置记为pbesti(简记pbi)，整个种群中目前搜索到的最优位置记为gbesti(简记gbi)，每个粒子根据以下公式更新自己的飞行速度和位置

vi(t+1)=w·vi(t)+c1·r1·(pbi(t)-xi(t))+

c2·r2·(gbi(t)-xi(t))

(11)

xi(t+1)=xi(t)+vi(t+1)

(12)

其中t为当前进化代数，w为惯性因子(一般取0～1)，c1,c2为学习因子(本文取c1=1,c2=2)，r1,r2为区间[0,1]之间的均匀随机数。具体寻优步骤如下：

步骤1输入所需要的参数，除了以上所提到的w,c1,c2,r1,r2,N还要输入最大迭代次数maxD；

步骤2种群初始化，不同的粒子代表不同的变化因子m，且以随机的方式求出每个粒子的初始位置与初始速度；

步骤3构造适应度函数，每一个m值对应一个DCT基，调用内点法(其求解过程在2.3节中详细说明)子程序求出与DCT基对应的稀疏解s，由式(1)即可求出与之对应的包络线x(t)，进而由式(5)得到不同的m所对应的平稳度值，用以下公式作为判断每个粒子(即m)好坏的适应度函数

(13)

步骤4更新粒子速度，根据上面式子(13)中的适应度值找出局部最优值pbesti与全局最优值gbesti，然后根据式(11)、(12)更新每个粒子m的速度与位置；

步骤5按照文献[14]进行粒子变异，得到变异后的粒子群；

步骤6检验终止条件，当达到最大迭代次数或者全局最优位置满足最小界限时，则终止迭代；否则转到步骤3。

2.3 基于粒子群的稀疏复原算法总体流程

利用混合变异粒子群算法在每次迭代过程中都需要由当前m值所构造的DCT基和观测值yd(t)来重构包络信号x(t)，所以稀疏复原方法的选择对整体的算法流程至关重要。本文选择了内点法，下面详细介绍怎样利用内点法解决式(10)中的模型。

步骤1参数设置，输入初始点s0,u0，惩罚因子τ0>0，衰减系数ρ∈(0,1)，终止误差0≤ε≤1，迭代次数为k=1；

步骤2根据式(9)构造惩罚函数P(s,u)，以sk-1为初始点利用共轭梯度法求解式(10)中的无约束子问题f(s,u)，求得极小值点sk；

步骤3判断终止条件，若τk·p(s,u)≤ε，停算，输出s*≈sk作为稀疏解；

步骤4令τk+1=ρ·τk，k=k+1，返回步骤1。

综上所述，基于粒子群的包络线稀疏重构最优化算法总体流程如图1。

图1 粒子群算法的总体流程图

粒子群算法在每次迭代过程中通过调用上述内点法子程序对求得的m寻优，直到迭代终止时找到最优的变化因子m*,进而得到最适合包络线变化特性的DCT(m*)基，最后利用内点法复原出原始包络信号。内点法的引入提高了稀疏复原的速度，可用于复杂信号包络的提取，同时利用粒子群算法解决了局部寻优变化因子m的缺陷，较大程度上提高了包络线提取精度。

3 实验仿真

为了验证稀疏复原方法提取包络线的有效性，本文使用的仿真信号为

x(t)=1.8e-3tcos(2πf1t+1.3)+

0.6e-3tcos(2πf2t+0.7)+

0.1e-2tcos(2πf3t+0.7)

(14)

仿真信号中的模态频率分别取f1=0.8 Hz、f2=3 Hz、f3=10 Hz，时间t的取值范围为[0,3]s，信号采样点数为500个。为了突出本文改进算法的优越性，结合粒子群寻优算法对式(14)中的信号采用内点法复原包络信号，与三次样条插值算法和分段幂函数算法对比，其结果如图2，同时与文献[9]的结果(见图3)相对比。

图2 基于粒子群的内点法复原方法提取信号上包络

由实验结果可知，基于粒子群的内点法复原信号上包络线明显地抑制了端点效应问题，通过与分段幂函数法结果对比，明显的解决了端点处过冲与欠冲问题，而且自适应地找到改变DCT基频带宽度的变化因子m，从而自适应地寻找到一条符合包络线缓慢变化特性的最佳曲线。与文献[9]的结果相对比，本文提取包络更加平滑，其中dd=0.02，而文献[9]中dd=0.04，而且文献[9]中通过m=m+1人为筛选取得m=18时得到最优包络，而本文利用混合变异粒子群算法自适应地在[0,1]之间寻找全局最优m值(m=4.707 426×10-2)取其倒数大概是m=23，不仅减小了人为确定所带来的误差，解决了局部寻优的片面性，而且提高了提取精度。为了更好地观察粒子群寻优过程，图4给出了每次迭代过程中适应度函数所对应的值。

图3 基于稀疏复原方法(OMP)提取信号上包络线

图4 内点法复原包络时粒子群算法寻优m过程

从图4中可以看出，在0～10次迭代过程中，寻优速度较快，在10～25次时，开始逐渐靠近最优值，25～40次时，寻优值基本上保持不变。所以，利用粒子群对变化因子m寻优不仅可以精确地找到m值以提高包络线的提取精度，而且可以快速地找到变化因子m以提高包络线的提取速度。同时，从统计出的粒子群寻优程序运行的平均时间(见表1)来分析本文算法的运算复杂度。本数据是通过大量的实验仿真运行统计得到的，其运行环境如下：MATLAB 2014b，Intel Core i3处理器，3.3 GHz的CPU处理器，4 GB内存。

从表格中可以看到，当粒子数目一定时，随着迭代次数的增加，程序运行的平均时间逐渐减少，从而说明达到一定的迭代次数时可以找到最优的粒子，当超过这个迭代次数时，粒子继续保持最优值。通过大量实验，我们得到这个迭代次数大概为30，所以无需过多迭代次数即可得到全局最优值；还可以看到，随着粒子数目的增多，程序运行的平均时间逐渐增加，但是通过仿真数据观察到，当粒子数目达到一定数量时，提取的包络误差保持不变，从而说明并不是粒子数量越多越好，通过大量实验得到最佳粒子数目为20，所以本文中的仿真取粒子数N=20。

表1 粒子群寻优算法的时间统计

为了检验改进算法对信号处理的抗噪性能，仍然采用式(14)仿真信号x(t)，然后给x(t)添加不同信噪比的高斯白噪声n(t),得到加噪信号x1(t)，然后对x1(t)用本文提出的基于粒子群的稀疏优化算法提取包络线。图5中给出的是信噪比取15 dB时的提取结果,从图中可以看出，本文方法可以自适应地提取出最佳包络线，消除了端点效应问题，而且除端点外提取的结果几乎和三次样条插值、分段幂函数算法相吻合。与文献[9]中取同样信噪比的结果(见图6)相比，本文提出的方法不仅具有一定的抗噪性，而且提高了包络线的精度，明显减小了提取误差(dd=0.13)，同时对比图6中的m=2,本文对m进行了全局寻优。

图5 信噪比15 dB时基于粒子群的内点法复原上包络

Fig.5 The interior-point methods based on PSO for extraction the upper envelope of the signal added noise when SNR is 15 dB

当噪声强度增强时，选取信噪比为10 dB的加噪信号再次进行仿真，由图7、图8可以看出，本文方法仍具有一定的抗噪性能，不仅可以解决端点效应问题，而且误差仍可以控制在一定的范围之内，其中dd=0.22，而文献[9]中的dd=1.41。随着信噪比的降低，本文提取信号包络线的精度没有明显降低，相比文献[9]有了更强的抗噪性能。

图6 信噪比15 dB时基于稀疏复原(OMP)提取信号上包络

Fig.6 The OMP based on sparse recovery for extraction the upper envelope of the signal added noise when SNR is 15 dB

图7 信噪比10 dB时基于粒子群的内点法复原上包络

图8 信噪比10 dB时基于稀疏复原(OMP)复原上包络

4 实际斜拉索振动信号的应用测试

在上述实验仿真中，为验证本文算法提取的包络线的抗噪性能，人为添加了高斯白噪声，但是自然界中噪声模型很复杂，为了验证包络线提取在EMD算法中的应用，对实际采样信号进行了仿真与EMD分解。本文的实测信号采集于宁波市某一斜拉索大桥，此桥总长67 m，由102根直径0.15 m拉索构成支撑系统；采用WS-ZHT2振动设备和双传感器，并且将双传感器装置在拉索和梁端的铰支上[15-16]，如图9所示。

图9 传感器安装图

为了验证本文提出的改进后的基于稀疏复原提取包络线方法的有效性，本实验仿真中采样频率fs=100 Hz采样512个点。图10是采用本文改进算法对采集信号的包络线提取结果，并且与三次样条插值结果对比，由图看出，基于粒子群的稀疏优化方法仍可以精确地提取采集信号的包络线，并且数值保持在允许误差之内，再次验证了本文提出方法的可行性和实用性。

图10 采用基于粒子群的稀疏优化和三次样条插值、分段幂 函数法对实际采集信号提取上包络的比较

Fig.10 The comparison chart on extracted the envelope of actual acquisition signal based on sparse optimization and based on cubic spline interpolation method and subsection power function method

EMD分解将非平稳、非线性信号分解为有限个本征模函数IMF，这是Hilbert Huang变换(HHT)的关键步骤。图12是采用基于粒子群的稀疏优化方法提取采集信号包络线的EMD分解结果，其中IMFi表示本征模函数的阶数，随着阶数的增加，频率依次降低。与文献[9]中的结果(见图11)对比可以看出基于粒子群的稀疏复原方法经EMD分解得到的IMF分量在端点处更加精确，图形表现的更加平稳，更接近于正弦函数，这一点在低频分量上表现的更为明显。

图11 基于稀疏复原算法对实测信号EMD分解结果

图12 基于粒子群的稀疏复原对实测信号EMD分解结果

5 结 论

EMD分解的关键步骤是信号上下包络线的提取，本文将包络提取模型中的凹问题转换为凸二次规划问题，利用内点法求得模型中的稀疏解，并且将粒子群优化算法与稀疏复原算法相结合，利用混合变异粒子群算法，将改变DCT基频带宽度的变化因子m作为粒子属性，通过多次迭代获得全局最优m，最终自适应地提取出最优的包络线信号。经验证该方法可以有效解决端点效应问题，提高包络线的提取精度与抗噪性能，粒子群自适应寻优方法的引入，减少了人为选择所产生的提取误差，从实际桥梁振动信号的模态分解过程中，证明了本文方法的有效性和实用性。

[1] HUANG N E, SHEN Z, LONG S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Proceedings of the Royal Society A Mathematical Physical & Engineering Sciences, 1998, 454(1971):903-995.

[2] DELECHELLE E, LEMOINE J, NIANG O. Empirical mode decomposition: An analytical approach for sifting process[J]. IEEE Signal Processing Letters, 2005, 12(11):764-767.

[3] CHENG J, YU D, YU Y. The application of energy operator demodulation approach based on EMD in machinery fault diagnosis[J]. Mechanical Systems & Signal Processing, 2007, 21(2):668-677.

[4] FLEUREAU J, NUNES J C, KACHENOURA A, et al. Turning tangent empirical mode decomposition: a framework for mono-and multivariate signals[J]. IEEE Transactions on Signal Processing, 2011, 59(3):1309-1316.

[5] 杨云飞. 基于EMD分解技术的不同市场原油价格相关性分析及预测研究[D].武汉：华中科技大学,2011.

[6] 钟佑明. 希尔伯特-黄变换局瞬信号分析理论的研究[D]. 重庆：重庆大学,2002..

[7] 张帆, 丁康. 广义检波解调分析的三种算法及其局限性研究[J]. 振动工程学报, 2002, 15(2):243-248.

ZHANG Fan，DING Kang.Research on the three algorithms and limitations of generalized detection-filtering demodulation analysis[J].Journal of Vibration Enginering，2002, 15(2):243-248.

[8] 张郁山,梁建文,胡聿贤. 应用自回归模型处理EMD方法中的边界问题[J]. 自然科学进展, 2003,13(10):48-53.

ZHANG Yushan，LIANG Jianwen，HU Yuxian. The processing of end effects in EMD method by autoregressive model.[J]. Progresses in Nature Science，2003,13(10):48-53.

[9] 徐静妹，叶庆卫，王晓东，等. 斜拉索振动信号的包络线稀疏复原算法[J]. 振动与冲击，2015，34(23)：187-193.

XU Jingmei，YE Qingwei，WANG Xiaodong，et al. Envelope sparse recovery algorithm for stay cable vibration signals[J]. Journal of Vibration and Shock，2015，34(23)：187-193.

[10] CANDES E J, TAO T. Near-optimal signal recovery from random projections: universal encoding strategies?[J]. IEEE Transactions on Information Theory, 2007, 52(12):5406-5425.

[11] BOERINGER D W, WERNER D H. Particle swarm optimization versus genetic algorithms for phased array synthesis[J]. IEEE Transactions on Antennas & Propagation, 2004, 52(3):771-779.

[12] MALLAT S G, ZHANG Z. Matching pursuits with time-frequency dictionaries[J]. IEEE Transactions on Signal Processing, 1993, 41(12):3397-3415.

[13] KIM S J, KOH K, LUSTIG M,et al. An interior-point method for large-scalel1—regularized least squares[J]. IEEE Journal of Selected Topics in Signal Processing, 2007, 8(3):1519-1555.

[14] 陈君波, 叶庆卫, 周宇,等. 一种新的混合变异粒子群算法[J]. 计算机工程与应用, 2007, 43(7):59-61.

CHEN Junbo，YE Qingwei，ZHOU Yu，et al. New multi-mutation particle swarm optimization algorithm[J]. Computer Engineering and Applications, 2007, 43(7): 59-61.

[15] YE Q, FENG Z, YUAN D. Nonlinear vibration signal tracking of large offshore bridge stayed cable based on particle filter[J]. Polish Maritime Research, 2016, 22(4):70-75.

[16] YE Q, FENG Z, ZHAN Q. Vibration signal fundamental frequency estimation of bridge stayed cable based on era algorithm.[J].Journal of the Balkan Tribological Association 2016, 22(2):1335-1343.