基于MPE-CEEMD算法的轴流压缩机旋转失速状态识别*

姚金鑫 杨树华 王培林 太兴宇 李宏坤

（1.大连理工大学机械工程学院;2.沈阳鼓风机集团股份有限公司；3.闽西职业技术学院信息与制造学院）

0 引言

早期的研究认为压气机进入旋转失速是由本质上二维、线性、低频、小幅值的模态波扰动诱发，其以20%～50%转子旋转速度沿周向传播。其后，Day在研究中发现压气机在失速前会出现尖脉冲型扰动信号，其本质上是一种大幅值、三维、高频率、非线性的扰动，其一般以60%～80%的转子旋转速度沿周向传播[1]。西北工业大学的强冠杰等人以单级低速轴流风扇为研究对象，采用动态探针测量捕捉到了“尖峰型”失速先兆波及其特征，并且测量得到的扰动波的初始转速均约为70%转速[2]。因此，可以通过分析压缩机内流道声压信号的特征频率，来进行旋转失速的试验研究。

当前国内外对压缩机旋转失速状态识别的主要方法有时域分析、频域分析、极坐标图、时频分析法和基于混沌理论的分析方法等。本文采用了一种基于多尺度排列熵（MPE）值检测CEEMD集成平均分量随机性的改进的CEEMD算法，通过分析降噪后的声压信号傅里叶变换后各个频率对应的幅值来判断压缩机是否处于旋转失速状态。

1 MPE-CEEMD算法

补充的经验模态分解（Complementary Ensemble Empirical Mode Decomposition，CEEMD）是通过在原始信号中多次加入大小相等，方向相反的高斯白噪声，最后以集成平均的方式来消除白噪声对原始信号的影响[3]。事实上，在信号分解过程中，引起模态混叠的成分以及高斯白噪声信号会首先被分解出来，在分解出异常信号之后，信号逐渐平稳，极值点分布较为均匀。因此，郑近德等人提出了基于排列熵检测IMF分量随机的MEEMD算法，但并没有考虑各个尺度下的排列熵值。本文提出了一种基于多尺度排列熵（MPE）值检测CEEMD集成平均分量随机性的改进的CEEMD算法，与传统的排列熵相比，各个尺度的排列熵值均被考虑在内，具有更好的抗干扰能力和适应性，能更有效的分析时间序列信息。

1.1 MPE-CEEMD算法分解步骤

MPE-CEEMD算法剔除了EEMD，CEEMD算法中合成信号自适应分解之后的异常分量，避免了不必要的集成平均，在保证分解的完备性的同时，减小了EEMD和CEEMD的计算量，使得到的IMF分量更具有意义。对于非平稳信号S(t)，MPE-CEEMD算法的分解步骤如下：

图1 算法流程图Fig.1 Algorithm flow chart

1）在原始信号S(t)中，加入大小相等方向相反的2组白噪声，加入的白噪声均方根值应接近于原始信号噪声，或不超过0.3倍。即：

2）用EMD对Ne加入噪声后的信号依次做经验模态分解，从第一对开始，依次得到第一阶固有模态函数，

如果集成分量多尺度排列熵的平均值MMPE值大于某个阈值，那么该集成分量的随机性较强，被认为是异常分量。

4）如果I1(t)是异常分量，则继续执行步骤（1），直到EMD分解后的集成平均分量Ip(t)不是异常信号。

5）将已分解的前 p-1个分量从原始信号中分离出来，即：

并对剩余信号r(t)进行EMD分解，并将所有得到IMF分量按照高频到低频排列。

1.2 多尺度排列熵（MPE）

排列熵是一种检测一维离散时间序列动力学特征的指标，能够准确有效的提取出信号的微弱特征，在海洋、大气、机械故障诊断和脑电信号分析等领域得到了广泛的应用。单独的排列熵只能获取信号在某个尺度上的动力学信息[4]，为了更充分地研究样本的特征信息，本文将多尺度排列熵引入轴流压缩机旋转失速特征识别中，并将压缩机旋转失速信号CEEMD分解后的集成平均分量在尺度[1,5]范围内的多尺度排列熵的平均值用来作为筛选IMF分量的标准。多尺度排列熵基本思想是首先将时间序列进行多尺度粗粒化，然后计算其排列熵[5]。假设长度为L的时间序列xL={ }x1,x2,…xL的粗粒化过程如下：

对于任意一个重构后的时间序列Yk，都可以表示为符号序列s(l)={ j1,j2,…,jm} ，其中，{ j1,j2,…,jm} 表示相空间重构后时间序列中的各个元素原始位置索引，其中，l=1,2,3…K，且K≤m!。所以，相空间重构后的时间序列为：

{ j1,j2,…,jm}共有m!种排列可能，s(l)只是其中的一种，若P1,P2,…,PK表示每种序列出现的概率，此时根据香农熵的定义，离散序列X的k种符号序列具有的香农熵Hp为：

当Pj=1m!时，Hp达到最大值ln(m !)，对尺度s下的排列熵进行归一化可得到归一化的多尺度排列熵：Hp=Hp(m )ln(m !)。Hp值反映了时间序列的随机性，放大了序列的微小变化[4]。对于尺度l(l =1,2,…,s) ，计算时间序列X={x (i),i=1,2,3,…,n}在所有尺度下的排列熵MPE(l)={Hp(1),Hp(2),…,Hp(s)}，定义多尺度排列熵的最小值：

将MMPE值作为筛选集成平均分量的指标，所有尺度下的排列熵值都被考虑了进来，将MMPE值大于阈值的集成平均分量继续分解，直至满足终止条件为止。

1.3 多尺度排列熵参数的选取

在多尺度排列熵参数的选择中，有三个参数需要确定，尺度因子s，嵌入维数m，时间延迟τ和多尺度排列熵阈值（MMPE）。嵌入维数m和实践延迟τ决定相空间重构后时间序列的长度，若选取的过大，则会出现片段细节信息丢失，无法反映时间序列的细微变化；若选取的过小，那么相空间重构的片段过短，无法反映整体的片段信息。综上所述，选取的嵌入维数m=5。时间延迟τ对时间序列的计算影响较小，本文选择τ=1。从图3中可以看出，在尺度因子超过4之后，声压仿真信号的多尺度排列熵（MPE）值随着尺度的增加而越来越小，因此尺度因子区间选择的过大对于剔除集成平均分量中的随机分量是不利的，因此尺度因子的选择区间为[1,5]。

风洞压缩机在运行的过程中，会产生大量的非线性，非平稳的噪声信号，主要由四部分组成，空气动力噪声信号；机壳、管路、电动机轴承等辐射的机械式噪声信号；电动机的电磁噪声信号；压缩机振动通过基础辐射的固体声信号[7]。故而构造以下几个代表性信号：

x1(t)为均值100，标准差为20的，长度为1 024的白噪声。

仿真信号x2(t)表示信号的主要成分-空气动力噪声信号，主要包含轴频，叶通频，旋转失速频率成分，将电动机的电磁噪声构造为调幅信号x3(t)；将压缩机振动通过基础辐射的固体声信号构造为调幅调频信号x4(t)，高频噪声信号表示为x5(t)。故而构造的压缩机实际运行状态下的仿真信号为：s(t)=x1(t)+x2(t)+x3(t)+x4(t)+x5(t)，其中空气动力噪声信号，也就是压缩机出口机壳处的压力脉动信号构造为s1(t)=x2(t)，噪声仿真信号构造为s2(t)=x1(t)+x3(t)+x4(t)+x5(t)。仿真信号的时域图如图2所示。

图2 仿真信号时域图Fig.2 Time domain diagram of simulation signal

分别计算压缩机实际运行仿真信号s(t)、声压信号s1(t)和噪声仿真信号s2(t)在尺度区间[1,5]范围内变化的多尺度排列熵值，并展示在图3中。

图3 仿真信号多尺度排列熵值计算Fig.3 The multi-scale entropy calculation diagram of the simulation signal

从图中可以看出，在加入噪声仿真信号后，压缩机实际运行状态的仿真信号的多尺度排列熵值和噪声仿真信号的多尺度排列熵值基本一致，且大于声压仿真信号的多尺度排列熵值的最大值0.6。通过将尺度[1,5]范围内多尺度排列熵均值大于0.6的集成平均分量剔除出去，实现对信号的随机性检测，这个方法可以应用于轴流压缩机旋转失速特征提取中，且多尺度排列熵的阈值可以确定为0.6。

2 实验介绍

实验对象为沈鼓集团某型风洞压缩机，该型压缩机有四级动叶，四级静叶，每级26个叶片，叶片材料碳纤维，主轴电机额定功率80 000kW/h。实验地点位于营口透平试验基地，实验设备布置在厂房的低台位区。

测试系统采用课题组自主开发的风洞压缩机压力脉动在线监测系统，其中有6个声压通道，能实现数据的在线采集与离线分析。4级静叶后和出口管道截面内呈120°各布置3个声压传感器。声压传感器将采集到的信号通过电荷放大器转化为低阻抗电压信号后，再通过光纤交换机传输至上位机，上位机对信号进行实时的处理。采样频率10.24kHz，采样时长10s，取其中5s作为研究对象，数据点数为51 200个。数据采集系统的参数配置如表1所示。

表1 数据采集系统配置Tab.1 Data acquisition system configuration

NI CRIO基于LabView软件开发架构包含三部分，上位机LabVIEW、机箱RT实时控制器LabVIEW Real-Time、机箱背板FPGA。上位机实时控制LabVIEW RT程序运行于实时操作系统平台，与FPGA端程序进行数据交互[8]，LabVIEW FPGA模块对采集到的数据进行简单操作后，通过FIFO将采集到的数据传输到LabVIEW RT模块，RT端通过TCP协议将数据传输到上位机[9]，然后将采集到的数据进行简单的数据处理后，实现数据的实时采集功能。

为了方便对现场采集到的数据进行更为精确地分析，在软件中还加入了离线分析模块，该模块的主要功能是对数据进行时域分析、频域分析、包络谱分析和倍频程分析，另外，还可以将各通道每个工况下采集到的电压信号转化为工程单位，经A计权后用能量平均法转化为连续等效声级，用于评估压缩机蜂窝降噪装置的降噪性能。

在实验过程中，首先调整压缩机动叶至某一个角度，逐渐关小阀门，测量由低到高各个压比的试验点，直至发生旋转失速现象；此时应立即打开阀门，降低负荷，阀门开度以试验回路不发生旋转失速现象，恢复正常为宜。然后，调节动叶至下一个角度，之后增大阀门开度，直至能找到试验回路能循环的最佳压比点，此时可以开始进行下一个叶片角度的测试。最后，依次重复上面的操作，直至完成各马赫数（即不同转速）下各动叶角度的性能曲线测试。

3 实验数据分析

压缩机处于旋转失速状态时，旋转失速团的特征频率为低频，其频率低于转速频率。当压缩机组发生旋转失速时，该低频信号产生，并且随着流量的减小，其幅值逐渐增大。基于声压信号频谱分析，未发生旋转失速时其频谱主要频率成分为转速频率和叶片通过频率及其相关倍频，而发生旋转失速时会出现频率较低的旋转失速团特征频率[10]。本文通过分析压缩机内流道声压信号的特征频率，来进行旋转失速的试验研究。

3.1 失速频率的确定

图4 560r/min旋转失速频率识别Fig.4 560r/min rotating stall frequency identification

图4展示了压缩机性能试验过程中声压信号频谱的演化过程，由于变频器转速的误差，导致实际采集得到的频谱图中的频率成分会有微小的波动。从流量点1到流量点3流量依次减小，分别为754m3/s,719m3/s,629m3/s，压比依次升高，分别为1.057 0,1.107 1,1.205 4。从图4可以明显的看出转子的转频9.2Hz和9.4Hz，转频二次谐波频率18.6Hz和转频三次谐波频率28Hz。在压缩机流量点进入754m3/s时，该压缩机已经进入旋转失速状态。

对声压信号的频谱分析结果表明，在声压信号的频谱中，包含了强烈的5.6Hz和5.8Hz的频率成分，它们可以认为是在同一频段内，且该频率成分对应的幅值随着流量的变小而增加，如表2所示，所以旋转失速频率成分既包含5.6Hz，又包含5.8Hz的频率成分。C.Freeman在Viper MK522型高速单轴轴流式涡喷发动机上所做的试验和I.J.Day在剑桥4级轴流式压气机上所做的试验表明：对于同一类型的压气机，失速团的旋转速度与失速时压气机转速有一基本恒定比例关系[11]。性能试验逼喘时压缩机进入失速时的转速在560r/min，转子转动频率为9.2Hz和9.4Hz，即失速团以压缩机转速的60%恒定转速旋转。

表2 560r/min失速工况频率统计Tab.2 Statistics of frequency 560r/min stall condition

结合图5和表3可以看出在未失速声压信号的频谱图中，其频谱图中的主要频率成分为转频和叶通频及其相关倍频，另外，在小于2Hz的频率范围内存在1.6Hz和1.8Hz的特征频带，这些频率成分同样存在于压缩机机械运转时未发生旋转失速的声压信号频谱图中（如图5所示），因此这些频谱图中的低频成分可以认为是噪声干扰，应予以剔除。

图5 500r/min,550r/min,600r/min未失速频谱图Fig.5 Spectra of 500r/min,550r/min,and 600r/min without rotating stall

表3 未发生旋转失速频率统计Tab.3 Data statistics of frequency without rotating stall occurred

3.2 数据分析

当压缩机处于旋转失速状态时，在频谱图上存在小于轴频，且随着流量的减小，而幅值增大的特征频率成分。选取压缩机气流处于1.3Ma，叶片角度20°，压比1.107 1下的数据进行傅里叶变换，此时的流量为719m3/s。从图6的时域部分和表4中可以看出，由于噪声污染，从时域波形中无法观察出压缩机是否处于旋转失速的状态。在图6的频域部分和表4中可以看出，小于轴频范围内较明显的特征频率有1.4Hz，3.6Hz和5.8Hz，其中，1.4Hz的特征频率为低频噪声干扰，3.6Hz特征频率为未知频率，5.8Hz属于旋转失速特征频率。

表4 实测信号FFT分析的频率成分Tab.4 Frequency components of measured signal FFT analysis

图6 实测信号FFT分析Fig.6 FFT analysis of measured signal

对风洞压缩机内流道声压信号做MPE-CEEMD分解，首先将多次加入白噪声的原始声压信号逐级分解为多个不同特征时间尺度相对平稳的固有模态函数（IMF），集成平均后计算所有集成平均分量的尺度为[1,5]区间内多尺度排列熵的平均值，检索随机性程度较大的集成平均分量并将其中的间歇信号，脉冲信号和干扰信号等异常成分一次性剔除。最后，对剩余信号再次做EMD分解，得到7个固有模态函数IMF，如图7所示。

图7 信号分解图Fig.7 Decomposition diagram of measured signal

图8 实测信号MPE-CEEMD分析Fig.8 MPE-CEEMD analysis of measured signals

表5 实测信号MPE-CEEMD分析频率成分Tab.5 Frequency components of measured signal MPECEEMD analysis

重构后信号的时域图和频域图如图8所示。从时域部分可以看出，信号存在明显的脉动成分，能明显的识别出压缩机旋转失速故障。从图8的频域部分和表5可以看出，在频谱图上小于轴频的频率成分中，1.4Hz的低频干扰成分已经消失，仅存在3.6Hz的未知频率和5.8Hz的旋转失速频率成分。

4 结论

1）通过对传统的单尺度排列熵进行改造，建立以多尺度排列熵均值（MMPE）为筛选CEEMD分解后的集成平均分量的体系，考虑了尺度[1,5]范围内的每一个集成分量的多尺度排列熵值。通过将CEEMD与多尺度排列熵结合，在保证小剩余噪声的情况下，不仅可以降低重构次数，节省计算时间，还可以最大限度的剔除信号中的间歇信号，脉冲信号和干扰信号等异常成分。本文提出的信号处理方法抑制低频噪声效果明显，具有一定的工程应用价值。

2）通过使用压缩机出口机壳处的声压信号来实现旋转失速的特征频率，不需要额外布置失速探头，既简化了测试系统的结构，又为企业降低了成本。

3）相比于MEEMD，MPE-CEEMD算法将各个尺度的排列熵值均被考虑在内，这虽然增加了计算时间，但总体上增加的时间并不多，这对于压缩机旋转失速状态的在线监测具有一定的意义。