液体火箭发动机地面抽检试验是考核发动机工作状态和性能的重要手段。涡轮泵的技术指标是发动机健康寿命评估系统的重要参考指标之一，其工作状态对液体火箭发动机的正常运行起着举足轻重的作用^[1]。

火箭发动机涡轮泵是液体火箭发动机中最重要、工作条件最恶劣的组件之一。主要功能是将来自推进剂贮箱的低温液体推进剂进行增压，保证在一定压力和流量下将推进剂输送至推力室。

作为液体火箭发动机的主要组成部分，涡轮泵的工作状态对控制发动机推力和推进剂各组元的比例起着决定性作用。行业内多年故障统计数据显示：涡轮泵工作异常是发动机故障的主要来源^[2]。

涡轮泵在高速旋转中的剧烈振动是引发涡轮泵故障的重要原因之一。振动信号不仅包含涡轮泵工作的工频信息，还混杂着高速旋转所激起的振动混叠和噪声^[3]。通过分析试车中涡轮泵的振动数据，不仅能够研究其工作状态，为研究人员提供优化设计的依据，还能对故障进行定位，提高涡轮泵的使用寿命和可靠性，降低发动机试验成本，保障发射安全^[4]。

目前，对于非平稳型振动数据的分析方法主有快速Fourier变换(fast Fourier transform，FFT)、短时Fourier变换(short-time Fourier transform，STFT)、小波分析方法^[5]和Hilbert-Huang变换(Hilbert-Huang transform, HHT)^[6]方法等。FFT和STFT建立了频域和时域的相互转换关系，通过调整窗函数类型分析数据的频域特征。2种方法的本质都是Fourier变换，缺乏局部分析信号的能力^[5]。小波分析虽然实现了信号的多尺度、局部细化分析，但本质上还是Fourier变换和加窗分析，且分析过程中所选择的小波基函数对分析的效果影响较大。上述3种方法均并不具备自适应分析非平稳振动数据的能力^[7-8]。Huang^[9]提出了HHT理论。该理论的创新之处在于定义了瞬时频率的概念，引进了EMD和IMF的概念。通过对IMF分量进行Hilbert变换得到时频域的Hilbert能量分布频谱图，建立了谱分析与时域的对应关系^[10]，克服了FFT、STFT和小波分析的缺点，在工程应用领域得到了广泛的应用。薛志成等^[11]以HHT为基础，通过建立边际谱、边际谱频响函数和边际移动谱来识别建筑结构一阶模态频率时变特性。但由于离散型边际谱的自身特征和EMD分解过程的误差导致计算过程出现虚假峰值。虽然采取了加窗滤波的办法，但只能有效辨识结构的一阶频率，对本文的参考意义有限。文[12-13]利用HHT在时频分析领域的优点，分析信号的边际谱和能量谱，提取信号中的时频关键信息。文[14-16]采用HHT方法对设备运转原始信号进行分解，通过对比分析频谱信息与故障信息精准定位设备故障，较好地体现了HHT方法在故障定位方面的优势。

本文采用HHT方法对发动机试验涡轮泵振动数据进行分解，采用非闭合镜像延拓的方法解决了EMD过程中的端点飞翼问题。通过对Hilbert时频谱的分析实现涡轮泵故障定位，为型号试验人员和研究设计人员分析、定位涡轮泵故障提供了有力的依据。

1 HHT的基本理论、不足和解决方法 1.1 基本理论

基于HHT理论的数据分析方法主要以EMD和IMF的Hilbert变换为基础，采用瞬时频率表征信号变化。该方法主要分为2个基本步骤^[17]：

1) 通过EMD对信号进行预处理，将其分解为一些基本的IMF分量。

2) 对IMF基本分量进行Hilbert变换，获取时频平面Hilbert时频谱图。

在Hilbert-Huang变换理论中，定义瞬时频率为

$ \omega = \frac{{{\rm{d}}\theta \left( t \right)}}{{{\rm{d}}t}}. $

(1)

对于任意一个时间序列x(t)，总有Hilbert变换y(t)表示如下：

$ y\left( t \right) = \frac{1}{{\rm{ \mathsf{ π} }}}P\int_{ - \infty }^\infty {\frac{{x\left( {x' } \right)}}{{t - t' }}} {\rm{d}}t' . $

(2)

式中P为Cauchy原理值。当x(t)和y(t)为复共轭对时得到解析信号z(t)表示如下：

$ z\left( t \right) = x\left( t \right) + {\rm{i}}y\left( t \right) = a\left( t \right){{\rm{e}}^{{\rm{i}}\theta (t)}}. $

(3)

其中：

$ \alpha \left( t \right) = {\left[ {{x^2}\left( t \right) + {y^2}\left( t \right)} \right]^{1/2}}, \theta \left( t \right) = {\rm{arctan}}\left( {\frac{{y\left( t \right)}}{{x\left( t \right)}}} \right). $

分解得到的固有模态函数需要满足以下2个必要条件^[17]：

1) 在整个数据集中，极值点的个数和过零点的个数相等或者最多相差1个。

2) 在任意点，由局部极小值构成的下包络和局部极大值构成的上包络的均值为0。

条件1规定了对数据的限制，使得瞬时频率的假设能够成立。条件2通过局部极大、极小值信号组成的包络将传统数据分析领域中的全局限定调整为局部限定，有效避免了瞬时频率定义中因为包含有非对称的波形而导致的波动。

EMD分解信号生成IMF的步骤主要分为以下4步：

1) 识别时序函数x(t)的所有极值点。采用3次样条插值法拟合极大、极小值点的上、下包络线u(t)和v(t)。u(t)和v(t)满足关系式：v(t)≤x(t)≤u(t)。

2) 建立上下包络曲线的平均曲线m(t)。

$ m\left( t \right) = \frac{1}{2}\left( {v\left( t \right) + u\left( t \right)} \right). $

(4)

3) 用x(t)减去m(t)后剩余的部分就是IMF，用h₁(t)来表示。

$ {h_1}(t) = x(t) - m(t). $

(5)

4) 为了得到满足条件的IMF，进一步用h₁(t)代替x(t)，与h₁(t)相应的上、下包络曲线为u₁(t)和v₁(t)。重复上述过程，即

$ {m_{k - 1}}(t) = \frac{1}{2}({v_{k - 1}}(t) + {u_{k - 1}}(t)), $

(6)

$ {h_k}(t) = {h_{k - 1}}(t) - {m_{k - 1}}(t). $

(7)

直到h_k(t)满足IMF的2个必要条件时，得到第1个IMF，即C₁(t)和信号的剩余部分r₁(t)表示如下：

$ {C_1}(t) = {h_k}(t), $

(8)

$ {r_1}(t) = x(t) - {C_1}(t). $

(9)

对信号的剩余部分r₁(t)继续进行EMD分解，直到r_n(t)的值小于预先设定的值或者为某个单调信号时结束。从上述分析可以看出，IMF的瞬时频率和瞬时幅度均为时间的函数，原始信号x(t)可以看作是所有IMF分量以及余量之和表示如下：

$ x(t) = {C_1}(t) + {C_2}(t) + \cdots + {C_n}(t) + {R_n}(t). $

(10)

EMD的实质是对数据进行筛选。在筛选过程中不仅可以消除信号上的畸形波，还能对信号中高低不平的振幅进行平滑处理，使得IMF满足所设定的必要条件^[18]。

根据IMF的必要条件，在实际的应用中，很多信号本身很难满足第二个必要条件。Huang等^[19]提出限制标准差的大小来确定分离过程停止的标准，即

$ S=\sum\limits_{k=1}^{n}\left[\left|h_{1(k-1)}(t)-h_{1 k}(t)\right|^{2} / h_{1(k-1)}^{2}(t)\right]. $

(11)

在工程应用领域，一般将S的值限定在0.2~0.3^[20]。根据以往数据处理的经验，本文将S值设定为0.25。

IMF分量C₁，C₂，…，C_n的Hilbert变换表示如下：

$ x(t)=\sum\limits_{j=1}^{n} a_{j}(t) \exp \left(j \int \omega_{j}(t) \mathrm{d} t\right). $

(12)

由以上分析结论可以看出，HHT摆脱了FFT和小波分析的局限性，引入基于信号局部时间特性的IMF函数，获得了具有实际物理意义的瞬时频率表达式。IMF函数的Hilbert变换建立了联合时间频率和幅度分布之间的关系，能够在Hilbert谱图中精确描述信号中振动能量发生的时间，因此非常适合分析非平稳数据。

在工程应用中通常采用计算机编程的方式实现HHT数据分析。Hilbert-Huang变换的计算机实现流程如图 1所示。

1.2 HHT的不足

模态混叠和端点飞翼是HHT方法在应用中存在的2个主要问题^[17]。

1) 模态混叠。

模态混叠是指信号的时间尺度存在跳跃性的变化时，EMD方法很难从原始信号中分离出不同的固有模态函数，导致同一个IMF分量包括了多个模态且尺度差异较大的信号^[11]。在LRE试验过程中，试验过程均为时间连续型试验，若出现异常造成紧急关机，随即停止试验。模态混叠对数据分析的影响较小，本文不予讨论。

2) 端点飞翼。

对原始数据进行EMD分解得到数据的上下包络，再根据数据的极大、小值采用样条曲线进行逼近。由于所分析的信号长度通常是有限的，信号两个端点并不意味着就是极值点。因此采用3次样条插值法进行插值时，信号的上下包络会在信号两端附近产生严重的扭曲，严重影响信号中低频分量的EMD效果，并将误差传递到后续的分解过程，造成IMF函数失真，即端点飞翼问题。在对多分量复杂信号进行多次EMD分解时，端点飞翼的问题更为明显，有时甚至可以严重淹没信号的端部特征^[17]。

1.3 端点飞翼的解决方法

有多位研究人员对端点飞翼问题进行了研究并总结出了相应的解决办法。如去除端点数据、包络镜像闭合延拓、极值延拓、拟合多项式、固定端点延拓等^{[17, 21]}。但去除端点数据、极值延拓、拟合多项式和固定端点延拓的方法在处理数据过程中存在丢失数据、以偏概全的问题。液体火箭发动机试验数据往往包含发动机多个工况下的重要信息，非特殊情况一般不采取舍弃原始数据的方法。包络镜像闭合延拓法虽然增加了数据的处理量，但是保证了计算的精度，避免了数据丢失的问题。

镜像延拓，即参照镜子映射物体的方式，根据数据两端的数值分布特征，在数据两端选定某映射点，以该点为基准对某一段原始数据进行延伸。放“镜子”，即映射的位置通常选取放置在信号具有对称性极值所在的位置。延拓的方法分为闭合镜像延拓和非闭合镜像延拓^[17]。

例如某数据列长度范围为[0, 4T]，将数据列分为[0, T]、[T, 2T]、[2T, 3T]、[3T, 4T]4段。以数据列两端端点为映射位置，各延伸2T的长度。则映射后数据范围为[-2T, 6T]，得到2倍于原信号的数据序列，将其首尾相连构成一个环形闭合曲线。该方法称之为闭合镜像延拓。

若以数据列两端端点为映射位置，分别延伸长度小于1/2原始数据列长度，则称之为非闭合镜像延拓。例如在上述数据列中，在原始数据两端端点分别延伸T长度，得到映射数据范围为[-T, 5T]。

2种延拓的方法各有优势。闭合镜像延拓方法将原始数据处理成无端点闭合曲线，端点效应不影响EMD迭代分解，消除了人为影响。但缺点是处理的数据量是原始数据的2倍，大大增加了计算时间。非闭合镜像延拓方法不构造闭合曲线。优点是在延拓数据时减少了计算时间。缺点是该方法对映射点位置的选取要求较高，若映射点位置选取不合适将造成人为误差。

以液体火箭发动机试验振动数据为例，对比分析采用非闭合镜像延拓方法与不做延拓处理的EMD分解效果。

选取某切向振动参数的1 000个数据点，即长度为[1，1 000]。在数据列的两端分别以端点为映射点放置“镜子”，映射的数据长度为300个数据点，完成映射后的数据长度为1 600个数据点，即[-300，1 300]。对该映射数据列进行EMD分解。为便于观察，选取分解后IMF3数据列中[1, 100]数据进行画图分析。数据分析采用MATLAB编程实现。2种方法对比结果如图 2所示。

通过图 2的分析可以看出，非闭合镜像延拓方法与不做任何处理的数据EMD分解得到的IMF3分量波形在变化趋势上基本吻合。但是在数据列起点处，采用非闭合镜像延拓方法对数据列端点进行处理后，第1个波峰的幅值明显小于未采用镜像延拓方法的幅值。在第30个样本点之后，二者在幅值和变化趋势上基本重合。分析结果表明：镜像延拓方法能够有效解决端点飞翼的问题。

2 液体火箭发动机试验氧涡轮泵振动数据与故障分析

在某次试车中发现发动机喷管呈不规则收缩、放大状。试车结束后对发动机零部件进行检查发现燃气导管有细微裂缝。对破损导管进行探伤、裂纹及应力检测，排除了导管本身制造工艺、质量的问题。推断导致其产生裂纹的主要原因是外界应力、振动过大，超出了其承受极限。

对燃气导管上游的主要振动源--液氧涡轮泵及相关振动测点的数据进行分析。数据的采样率为10 000 Hz/s。各测点分析样本数据量均为1 000个样本点。发动机试车时间为150 s，分别对应工况1-工况3。3个工况的点火时间依次为20、60、70 s。发动机氧涡轮泵轴向、径向、切向振动测点名称分别为a1、a2、a3。推力室轴向、径向、切向振动测点名称分别为a4、a5、a6。

首先，采用FFT分析法在全频域范围对3个不同工况下氧涡轮泵振动数据进行频谱分析，发现氧涡轮泵在试验工况3的振动数据异常(工况3液氧泵转速频率为348.8 r/s)。该工况氧涡轮泵轴、径、切3个方向和推力室径向振动测点的FFT频率-幅值峰值数据分别如表 1-3所示。

分析表 1-3的数据可知，氧涡轮泵轴、径和切向振动测点在其工作频率(348.8 Hz)的倍频，即683.6 Hz(2倍频)、1 074.2 Hz(3倍频)、4 052.7 Hz(12倍频)、4 443.4 Hz(13倍频)、5 712.9 Hz(17倍频)、8 105.5 Hz(24倍频)处均存在较大的振动能量峰值。通过与以往同型号发动机相同工况下的试验数据进行对比分析，发现本次试验中工况3状态下氧涡轮泵径向和切向在4 000~8 000 Hz的频率区间内出现多个高幅值峰值。分析下游推力室径向a5振动测点的频谱数据发现该测点也存在相似的情况。由于FFT分析的局限性，无法定位振幅峰值点在时域中的位置。

根据行业内的研究统计数据，涡轮泵的故障多与转子质量不平衡有关。根据涡轮泵转子系统质量自由度信息，当转子系统质量失衡时，一般产生的是径向力^[22-23]。因此，以氧涡轮泵径向a2振动测点的数据为例，采用HHT方法对其进行分析。采用非闭合镜像延拓克服端点飞翼的问题。a2振动测点时域数据如图 3所示。

EMD分解所得到的8个IMF分量结果分别如图 4和5所示，第1-6层时间-频率特性曲线分别如图 6和7所示，EMD分解重构后的时间-幅值特性曲线如图 8所示。

从图 6和7可以看出，多层时间-频率曲线非常清晰地定位出时域内频率的变化情况，尤其是第2-6层的分辨率较高。在该时域内的0.15、0.24、0.34、0.43、0.49、0.67和0.68 s存在多个频率值突变。对比图 8中的时间幅值图数据，在频率变化点也基本对应振动幅值的峰值。第5层频率率曲线的均值基本对应了涡轮泵转子的一阶临界转速(689.320 Hz)，第3层频率率曲线的均值基本对应了涡轮泵转子的二阶临界转速(1 070.609 Hz)。这表明在该时间段内转子的临界转速已被激发。分析该测点后续时间段内的数据也存在相同的情况，且临界转速的数值基本不变。分析该振动测点启动段的时域数据，如图 9所示。从图 9中可知在启动段多个时刻存在多个高峰值冲击振动。

通过上述分析可以基本断定在该次试车中发生了涡轮泵转子次同步进动。在启动段，液氧涡轮泵的转速提高至临界转速时产生了自激振动。在试验工况3主级工作段时自激振动再次产生。

分析发动机试验工作机理，液氧从氧泵后分为2路。一路经氧主气蚀管及旁通的推进剂利用阀汇成一路，经氧主阀、燃烧室头腔进入燃烧室；另一路经氧调节器、气蚀管和副控阀进入燃发器。燃发器产生燃气分2路分流后分别驱动氢、氧涡轮后进入排气管。燃气分流的作用是调整氢氧涡轮泵的功率，进行转速匹配，从而保证发动机工作时推进剂流量和混合比^[24]。氧涡轮泵运转时的次同步进动所产生的高幅值自激振动向下游其他组件传递，从而破坏了燃气导管。导管燃气外泄并改变了喷管延伸段的压力场，使得喷管在压力场的作用下呈不规则呼吸状收缩、变大。

最后，由图 4-8还可以得到以下结论：

1) 由图 4和5可知，对原始振动数据进行EMD分解所得到的IMF函数图可以清晰地观察到各IMF分量所包含的频率成分、相应幅值大小和周期变化情况。信号各阶IMF两端较为平稳，不存在明显的端点效应。这表明非闭合镜像延拓克服端点飞翼的效果较好。

2) a2振动数据的能量主要包含在IMF1~IMF5分量中，信号中的突变成分主要分布在IMF1分量中。与其他分量的频率相比，IMF1分量的频率最高，频率分量主要涵盖了涡轮泵转子的工作频率分量。涡轮泵工作振动频率主要分布在IMF2分量中。IMF3分量包含了涡轮泵转子低频振动引起的频率分量。IMF2和IMF3代表了涡轮泵工作中的振动分量。IMF4和IMF5分量包含部分涡轮泵工作中的低频噪声分量及其附近的一些低频分量。IMF6、IMF7和一阶余量R主要为系统中一些其他的不相关低频分量，变化过程非常缓慢且无规律，不影响数据的分析结果。

3) EMD分解可以自适应地分离发动机氧涡轮泵振动信号的基频分量，并提取与各基频分量相对应的模态特征。对于振动数据所代表的的加速度特征而言，分解的结果具有明确的物理指向意义。

4) 图 6和7清晰地显示了原始振动振动数据中的频率成分以及频率随时间的动态变化规律。第1层时间-频率特性曲线以高频分量为主。随着时间的推移，主要频率在2 500~5 000 Hz间变化非常剧烈。第2层时间-频率特性曲线在1 000~1 500 Hz间变化。第1层和第2层时间频率特征曲线明显呈现振动分量的频率变化特征，尤其是第2层特性曲线清晰显示了该时间段内频率突变的时间点。第3~6层时间-频率特性曲线则非常清晰地显示了低频分量的变化相对较缓慢的特性，也能非常清晰定位对应的变化时间点。

5) 对比图 3与图 8可知，将IMF分量进行重构去除了原始数据中的部分噪声，提高了数据的分辨率。

3 结论

涡轮泵的工作状态直接影响液体火箭发动机的性能、寿命和可靠性。基于HHT的数据分析方法弥补了FFT的不足，能同时在时域和频域内对非平稳态数据进行分析。通过对某次试车涡轮泵振动的分析表明：HHT方法能够实现涡轮泵故障分析和定位。采用非闭合镜像延拓的方法能有效克服HHT理论中端点飞翼的缺点。HHT方法具有良好的自适应特性，在时间和频率2个维度的分析优势非常明显。该分析方法为型号设计人员和试验测量人员分析非平稳性数据提供有力的支撑，在液体火箭发动机、组合件以及及其他相关航空、航天型号试验数据分析中也具有重要的推广和实用价值。