基于希尔伯特-黄与小波变换的压气机失速信号分析方法

向宏辉，侯敏杰，杨荣菲，葛宁，刘志刚，吴晨

（1.南京航空航天大学能源与动力学院，南京210016；2.中国燃气涡轮研究院航空发动机高空模拟技术重点实验室，四川江油621703）

基于希尔伯特-黄与小波变换的压气机失速信号分析方法

向宏辉1，2，侯敏杰2，杨荣菲1，葛宁1，刘志刚2，吴晨2

（1.南京航空航天大学能源与动力学院，南京210016；2.中国燃气涡轮研究院航空发动机高空模拟技术重点实验室，四川江油621703）

针对多级跨声速压气机失速过程壁面静压测量数据，分别采用希尔伯特-黄变换（HHT）和小波分析方法进行时频分析。相比于小波分析方法，HHT时频图中失速过程的频率特征变化更为明显，且利用HHT图中失速发生后的特征频率可判断出失速团的旋转速度。结合方差分析，进一步探索了获取失速先兆信号的可行性。结果表明，两种方法均可用于压气机失速先兆预测，其中HHT方法能够提供更高层次的分解信号，得到的失速先兆更为提前。

压气机；失速信号；失速先兆；希尔伯特-黄变换；小波变换；方差分析；时频分析

1　引言

压气机的失速与喘振直接影响航空发动机的稳定工作范围，主动控制技术是实现发动机扩稳的重要技术手段。如果在失速完全发生后进行主动控制，由于流体惯性作用，往往需要很大的输入能量及较长的控制时间才能成功退出失速。反之，如果能在失速发生时刻前及时启动主动控制系统，则能在较短时间以较小能量（如叶尖微射流技术）实现失速抑制，确保航空发动机安全可靠工作。因此，如何尽早地获取与识别失速先兆信号，是压气机失速主动控制研究的一个重要内容。

压气机失速信号分析通常包括信号选取和先兆捕捉两方面内容。其中，选取的信号可以是原始采集信号，也可以是滤波处理信号。原始信号由于混入大量背景噪声，相对较弱的失速信号容易淹没其中，影响失速先兆捕捉精度。因此，失速信号分析过程中，首先需对原始采集信号进行滤波降噪处理。小波分析是应用较广的一种降噪方法，其降噪结果依赖于小波基函数与分解层次的选取，压气机失速信号通常选用Daubechies正交小波系，采用4至5层信号分解［1-2］。除此外，黄锷等［3-4］提出的经验模态分解（EMD）也是一种有效的滤波降噪方法。将此方法与希尔伯特谱分析相结合，得到了著名的时频分析方法，即希尔伯特-黄变换（HHT）。该方法相较于小波分析，最大优势在于分解结果不依赖于基函数，能自适应处理带噪声的非线性非稳态信号，频率分辨率较高。基于此，HHT在多个工程领域得到广泛采用，如机械故障诊断［5］、地震信号分析［6］等。在压气机失速信号分析方面，张勇等［7］利用HHT时频分析功能研究了小波降噪后离心压缩机试验数据，得到了压缩机喘振频率；郭贵喜等［8］利用HHT分析了轴流压气机失速信号，得到了失速前、后频率的变化特征；李长征等［9］利用HHT分析了小型压气机出口总压信号，通过判断时域图中频率突变点发生时刻，得到了压气机失速频率及失速起始时间；而符娆等［10］利用HHT分析了三级轴流压气机失速过程的小波降噪信号，通过预先设定阈值来监测失速先兆频带与失速起始时刻。上述研究主要集中在时频分析结果方面，而在失速先兆捕捉方面，将HHT分析方法与失速判断方法相结合的研究还未见公开报道。

在失速先兆捕捉方面，可分为两类方法：利用降噪信号的阶跃点进行失速判别与利用失速特征参数的阶跃点进行失速判别。其中，失速特征参数包括统计参数（如方差［1，11］、自相关系数［12-13］、互相关系数［13-14］）和混沌特征参数（如关联积分［15-16］、最大Lyapunov指数［17］、模糊熵［18］等）。混沌特征参数计算需首先进行相空间重构，会增加失速计算过程的时间与复杂程度；而统计参数仅针对信号进行直接处理，计算时间短。马翼等［1］对小波降噪信号进行方差分析提前获得了压气机失速先兆；高曼等［11］对比分析了原始失速信号的方差分析结果与自相关系数分析结果，认为亚声速压气机中自相关分析法得到的失速起始时间早于方差分析结果；李继超等［13］指出互相关系数法在失速先兆捕捉方面优于自相关系数法。

为探索高负荷多级轴流压气机动态压力信号分析与失速先兆捕捉的适用方法，本文以某六级跨声速高压压气机失速试验数据为研究对象，首先对比HHT与小波时频分析结果，然后结合失速先兆捕捉方法，分析EMD降噪与小波降噪用于提前获取压气机失速先兆的可行性。本文研究结果有望为高负荷跨声速压气机失速先兆信号及时捕获与主动控制策略实施提供技术基础。

2　试验样本数据

动态压力试验样本数据来源于某带可调进口导叶的六级轴流高压压气机［19］。该压气机设计转速为10 512 r/min，第1～2级转/静子叶片数分别为47、68、57、82。试验时，压气机进口保持均匀来流条件，分别在第1级转/静子之间、第2级转/静子之间的机匣壁面沿周向布置2支动态压力传感器。从压气机进口往后看，沿逆时针方向，第1级壁面2点脉动静压周向安装角度分别为50°、140°（分别用S1-P1、S1-P2表示），第2级壁面2点脉动静压周向安装角度分别为42°、86°（分别用S2-P1、S2-P2表示），静压测点采样频率为5 kHz。

试验用传感器选用美国PCB公司生产的M101A05型压电式动态压力传感器（固有频率为250 kHz），并采用半无限长管测试技术解决管腔的驻波效应。动态数据采集系统为Genesis系统，采样率6.4 kS/s，分析带宽2.5 kHz。试验中，通过调节试验器排气节气门开度改变压气机出口反压，实时记录压气机95%设计转速下，从稳定工作状态连续变化到失速工作状态的机匣壁面脉动静压信号。

3　信号处理方法

3.1HHT方法

HHT方法主要包括EMD和希尔伯特谱分析两部分。其中EMD类似于小波分解，通过筛选算法将原始信号分解为一系列不同尺度的本征模态函数（IMF），原始信号等于所有IMF分量与残差之和。希尔伯特谱分析类似于小波时频分析中的谱分析，得到的是信号瞬时频率，而不是小波谱得到的一小段时间窗口的计算频率。相比小波分析，希尔伯特谱分析的频率分辨率更高。

HHT分析的核心是EMD，利用时间序列的上下包络线平均值对信号进行筛选。筛选的终止条件为：信号极值点的数量与零点数相等或最多相差1个，信号由极大值点和由极小值点构成的上下包络线平均值为零。EMD具体分解过程如下：

（1）对原始信号X（t）的极大值点、极小值点采用三次样条函数插值，构造原始信号的上下包络线XU（t）、XL（t），以及包络线平均值M1=（XU（t）+XL（t））/2。

（2）考察H1=X（t）-M1是否满足筛选条件，如果满足则转到下一步，否则将H1作为原始信号继续操作步骤（1），求得M11及H11=H1-M11；依此筛选，直至第k步H1k=H1k-1-M1k满足筛选条件，则第一个IMF分量C1=H1k。

（3）将第一个残留值R1=X（t）-C1作为原始信号，按照上述方法继续筛选，得到第二个IMF分量C2（t）；依此类推，得到一系列IMF分量，直至Rn为单调信号或只存在一个极点位置。

由上述分解过程可知，EMD得到的IMF信号为不同分解层次的降噪信号，类似于小波分解重构信号，理论上可用于压气机失速先兆分析。另外，对不同分解层次的IMF信号进行希尔伯特谱分析，还可得到频率-幅值随时间的变化。

3.2方差分析法

时间序列Cj（t）的方差计算公式见式（1）。其中，N为计算窗口内的数据长度，tN为失速监测时刻。将Cj（t）进行分段滑动，滑动宽度为ΔN，则得到方差随时间的变化关系。

针对试验样本数据中的S1-P1原始信号，选取N为100，ΔN为5，计算得到方差时间序列如图1所示（图中纵坐标V代表方差）。可见，方差在稳定工作状态时基本为常值，接近失速时出现明显突变。定义突变点发生时刻为失速起始时刻。本文定义当前方差值为前100个方差均值的2倍时压气机开始进入失速，图中虚线为根据此判据得到的失速起始时刻，可见该方法能有效捕捉失速起始时刻。

图1　原始信号的方差Fig.1 Variance of raw signal

图2　S1-P1测量点信号的时频分析结果Fig.2 Time frequency diagram of Hilbert-Huang transform and wavelet at S1-P1

4　基于HHT与小波方法的失速信号分析

4.1时频分析

采用HHT方法与小波方法对失速信号进行时频分析，可得到信号频率与幅值随时间的变化关系，图2给出了S1-P1测量点的时频分析结果。对比发现：HHT方法得到的频率分辨率高于小波方法，在失速发展过程中，HHT时频图存在幅值占优的一个或几个频率，而小波时频图中幅值占优的频率为一定范围，不能明确与失速过程相关的主要特征频率，不利于失速发展过程中干扰频率判断。HHT时频图中504.40 s左右出现一个约10 Hz的低频脉动，且随着时间推移，脉动频率逐渐增加并稳定在78 Hz左右，此频率增加过程对应着从小尺度低频扰动逐渐发展为大尺度中频失速团的演化过程。失速发生后，对比转子转动一周产生的166.5 Hz转动频率，HHT时频图中78 Hz频率占转动频率的46.8%，可认为是失速团通过频率，由此得到失速团旋转速度为46.8%转子转速。由图3中S1-P1、S1-P2测量点失速信号峰值点确定的失速团旋转速度为47.6%转子转速（图中纵坐标ΔP代表压力脉动幅值），与HHT图中失速点频率判断出的失速团旋转速度相近。因此，采用单点HHT分析即可获得失速团旋转速度，相较于传统测量方案中利用周向两点测量确定失速团旋转速度而言，可以减少测量点数，降低试验成本。

图3　不同周向测量位置原始信号的时域图Fig.3 Time domain diagram of raw signal at S1-P1 and S1-P2

4.2失速先兆监测

首先利用HHT与小波方法对采集信号进行多重分解，获得失速先兆计算的待分析信号。图4给出了S1-P1测量点两种方法的分解结果对比，其中C0表示原始信号，C1～C8分别代表分解得到的IMF分量或小波分解重构信号分量。可见，两种方法得到的分解信号在高频率时基本重合，随着分解层次的增加，信号中的突变点前移，如果将信号剧烈变化的时刻看成失速起始时刻，则信号分解后计算得到的失速先兆时间相较于原始信号提前。但当分解层次过高时，如第8层，分解信号逐渐变为定常时均信号，信号中不再包含低频周期性失速信号，由此计算得到的失速先兆可能与真实结果相矛盾。

进一步将上述分解信号与失速先兆捕捉方法相结合，计算不同分解信号失速先兆起始时间。图5给出了各分解信号失速先兆的方差判断结果，从中可以看出，对于两种信号分解方法，随着信号分解层次的增加，方差判断得到的失速先兆时间提前。对于低分解层次，高频信号预测得到的失速起始时间相同；对于中等分解层次得到的信号（如C4～C7），基于IMF分量预测得到的失速时间略微早于基于小波分解预测得到的失速时间；对于高层次分解信号（如C8～C11），两种方法预测得到的失速起始时间呈现出相反的变化趋势。由图4可知，C8或更高层次分解信号的脉动逐渐接近于时均信号，这会对具有低频信号特征的失速先兆预测带来误判。

图4　本征模态函数分量与小波分解重构信号对比Fig.4 Comparison of intrinsic mode function and wavelet decomposition signal

图5　S1-P1测量点处本征模态函数分量与小波分解重构信号的失速先兆Fig.5 Stall precursor obtained from intrinsic mode function and wavelet decomposition signal at S1-P1

按照相同的分析步骤，对其他测量点（S1-P2、S2-P1和S2-P2）进行失速先兆判断，得到不同位置测量点的失速先兆随信号分解层次的变化规律，如图6所示。可见，随着分解层次的增加，不同测量位置的失速信号采用HHT方法均能获得一个最低的失速起始时间，而采用小波分解分析方法获得的最低起始时间所对应的分解层次，会随着测量点位置的变化而变化。因此，采用小波分解方法进行失速判断时，建议取第4～6层分解信号进行失速分析，可避免测量点位置对失速先兆判断的影响；而采用HHT方法时，建议取第6～8层分解信号进行失速分析，可获得最早的失速起始时间。相比之下，HHT方法获得的失速起始时间，要早于小波分析获得的失速起始时间，该结论在其他转速下也成立。

图6　不同测量点位置处本征模态函数分量与小波分解重构信号的失速先兆Fig.6 Stall precursor obtained from intrinsic mode function and wavelet decomposition signal at different measurement location

5　结论

（1）HHT方法相对于小波方法频率分辨率更高，从HHT时频图中可清楚地显示出压气机从小尺度低频扰动信号向大尺度中频失速团发展过程的频率变化。

（2）HHT时频图中失速发生后产生的频率为失速团旋转频率，利用该方法判定失速团旋转速度仅需一个测量点即可，减少了试验测量布点数量。

（3）采用HHT分解与小波分解重构获得的不同分解层次信号，结合方差分析法，均能提前预测压气机的失速先兆。相比于小波方法，HHT分解信号得到的失速先兆时间更为提前。

［1］马翼，陆凤霞.基于Daubechies小波的压气机失速分析及预测［J］.机械工程与自动化，2011，40（4）：123—125.

［2］姜涛，贾智伟，李应红.基于小波分析的压气机失速故障检测［J］.推进技术，2003，24（1）：55—57.

［3］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 London A，1998，454：903—995.

［4］Wu H，Huang N E.On the filtering properties of the empirical mode decomposition［J］.Advances in Adaptive Data Analysis，2010，2（4）：397—414.

［5］Wu F J，Qu L S.Diagnosis of subharmonic faults of large rotating machinery based on EMD［J］.Mechanical Systems and Signal Processing，2009，23：467—475.

［6］杨培杰，印兴耀，张广智.希尔伯特-黄变换地震信号时频分析与属性提取［J］.地球物理学进展，2007，22（5）：1585—1590.

［7］张勇，张春梅.用HHT变换处理离心压缩机喘振试验数据［J］.流体机械，2012，40（1）：10—12.

［8］郭贵喜，王振华，王洪祥，等.基于HHT的航空发动机气动失稳信号时频分析方法［J］.航空发动机，2013，39（3）：62—67.

［9］Li C Z，Xiong B.Application of HHT for aerodynamic instability signal analysis in compressor［C］//.3rdInternational Conference on Signal Processing Systems.Singapore：2011.

［10］符娆，雷勇.基于HHT的压气机气动失稳检测技术研究［J］.测控技术，2010，29（S0）：71—75.

［11］高曼，楚武利，吴艳辉，等.亚声速轴流压气机失速监测方法试验［J］.推进技术，2012，33（3）：398—404.

［12］Nakajima T，Tahara N，Kurosaki M.Stall prediction methodforaxialflowcompressor：US，6715984B2［P］. 2003-01-09.

［13］李继超，童志庭，聂超群，等.基于互相关分析的前失速先兆检测分析［J］.航空学报，2013，24（1）：28—36.

［14］Salunkhe P B，Joseph J，Pradeep A M.Active feedback control of stall in an axial flow fan under dynamic inflow distortion［J］.Experimental Thermal and Fluid Science，2011，35（6）：1135—1142.

［15］Tahara N，Nakajima T，Kurosaki M，et al.Active stall control with practicable stall prediction system using autocorrelation［R］.AIAA 2001-3623，2001.

［16］吴云，李应红，张朴.某涡喷发动机压气机气动失稳过程的非线性分析［J］.航空动力学报，2005，20（5）：890—894.

［17］杨荣菲，侯安平，周盛，等.最大Lyapunov指数用于压气机局部性能的分析［J］.航空学报，2009，30（4）：622—624.

［18］刘志刚，杨荣菲，向宏辉.基于模糊熵的压气机失速信号分析研究［J］.测控技术，2014，33（2）：53—56.

［19］向宏辉，侯敏杰，葛宁，等.多级环境下径向总压畸变影响效应的试验评价［J］.航空动力学报，2012，27（6）：1270—1277.

Analysis method investigation of axial compressor stall signal based on Hilbert-Huang transform and wavelet transform

XIANG Hong-hui1，2，HOU Min-jie2，YANG Rong-fei1，GE Ning1，LIU Zhi-gang2，WU Chen2
（1.College of Energy and Power Engineering，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，China；2.Key Laboratory on Aero-Engine Altitude Simulation Technology，China Gas Turbine Establishment，Jiangyou 621703，China）

A comparative study of Hilbert-Huang transform（HHT）and wavelet was presented through analyzing the wall static pressure acquired from a transonic multistage compressor which operated from stable to stall point continuously.Firstly，the time-frequency chart was compared.It is shown that the changes of frequency get more recognizable by HHT than by the method of wavelet in the development process of compressor stall.And the rotation speed of stall cells can be calculated based on the frequency that taking from HHT time-frequency chart at the time of compressor undergoing stall.Secondly，feasibility study on the stall precursor capture method was carried out which includes two steps of signal decomposition and variances calculation.The results suggest that the decomposition signal by HHT and wavelet can capture stall inception earlier than raw signal.As HHT can provide high levels of decomposition signal for stall prediction，it obtained earlier stall precursor than wavelet decomposition.

compressor；stall signal；stall precursor；Hilbert-Huang transform；wavelet transform；variance analysis；time-frequency analysis

V231.3

A

1672-2620（2015）05-0001-05

2015-04-28；

2015-10-25

江苏省普通高校研究生科研创新计划项目（KYLX15_0260）

向宏辉（1979-），男，湖南沅陵人，高级工程师，博士研究生，主要从事叶轮机性能评定与试验技术研究工作。