基于EDM-VKF 的转子时变复杂信号时频表示方法

付 波，李乐全，李超顺，王诗雯

（1. 湖北工业大学 太阳能高效利用及储能运行控制湖北省重点实验室，湖北 武汉 430068；2. 华中科技大学土木与水利工程学院，湖北 武汉 430074）

0 引 言

由于旋转机械运行条件的多变性，其信号中普遍存在非平稳特性［1，2］，而时频分析方法是表征非平稳信号时变特征的有效工具［3］，经典时频分析方法［4］包括短时傅里叶变换、小波变换和Wigner-Ville 分布。高传昌［5］等人基于经典时频分析方法分析了结构参数对压力脉动的影响，然而受海森堡测不准原理和交叉项的限制，经典时频分析方法的时频分辨率较低，难以实现非平稳信号的精确描述。近年来涌现出一系列后处理方法［6］来提高这些经典方法的时频分辨率，如重分配方法、同步压缩变换等。

上述方法是对单一信号进行分析，然而同源多传感器信号的联合信息可以更加全面地提供相关旋转机械健康状况重要信息［7-9］，吴峰崎［10］等通过将全谱方法扩展并结合级联全谱和支持向量机实现了旋转机械多种故障的分析和诊断。温广瑞等［11，12］将Vold-Kalman 滤波阶次跟踪技术和分数阶傅里叶变换与全息谱结合实现转子启停车故障特征的提取。王丽雅等［13］将全矢谱技术与短时傅里叶变换相结合提出了短时矢谱的概念，分析了压缩机喘振故障时振动矢量信号短时能量时频分布的规律。综上所述，时频分析方法结合以全频谱和全息谱为代表的信息融合技术可以实现转子系统非平稳振动信号分析，但大多数方法都是以短时窗稳定为假设，难以实现对机组非平稳过程的准确描述。

本文提出了一种转子时变复杂信号高分辨率时频表示方法，通过提取基于欧几里得距离的时频矩阵（Euclidean distance based the time-frequency matrix，EDM）的时频脊线获取精准的转子旋转速率，其次提取转子谐波信号的瞬时轨道特征，构建转子时变复杂信号高分辨率时频表示。

1 理论基础

1.1 Vold-Kalman滤波器

Vold-Kalman 滤波器［14-17］（Vold-Kalman Filter VKF）的基本思想是设置局部约束，声明未知相位分配的阶数是平滑的，并且这些阶数的和应该接近总振动信号。旋转机械的振动信号，包括第k次谐波分量，轴转速的第k次分量可写为：

式中：ak(t)表示幅度包络；fsf(t)是轴转速。

以ΔT为时间增量的上式的离散形式可以表示为：

在真实信号中，由于机械系统的惯性，与载波信号相比，信号变化缓慢。VKF流程函数被描述为：

其中εk(t)是高阶项，s表示微分调和。对于离散信号，上式可以修改为：

式中：εk(n)是非齐次项（在某种意义上是小的）；∇表示差分算子。

通过设置差阶s=2，上式变为：

考虑到所有样本，（5）的矩阵形式可以重写为：

其中M为N×N的矩阵。

对于由几个阶次的谐波分量组成的信号y（n），噪声可以表示为：

其中δ（n）是干扰噪声和误差。结果，上式的矩阵形式可以表示为：

其中C表示载波矩阵。

当通过时频分析给出或估计目标阶信号xk(t)的瞬时速度时，这意味着可以获得载波信号Θk(n)。为了求解上式以获得幅值包络，利用最小二乘法将非齐次项和噪声或误差项最小化：

式中：上标∗表示共轭转置；εT表示向量ε的转置；r是加权因子，用于通过平衡δ和ε来控制VKF带宽。

基于上式，给出目标分量的幅值包络矩阵A：

其中E表示单位矩阵。

在已知幅值包络矩阵A和载波矩阵C的情况下，目标谐波分量可以表示为：

1.2 基于全频谱的轴心轨迹瞬时特征参数

2 基于EDM-VKF 的转子时变复杂信号时频表示方法

2.1 基于EDM的时频脊线提取方法

转子旋转速度作为VKF 重要的输入参量，其准确性直接影响VKF 的提取能力，所以在转子多分量非平稳信号中精确提取瞬时频率非常重要。由此，本节提出基于EDM的时频脊线提取方法，通过引入欧几里得距离重构时频矩阵，提高时频聚集性，解决了传统方法无法高精度提取多分量复杂信号时频脊线的弊端。

此方法以时频矩阵第一列为起点，逐列计算更新矩阵列向量，然后根据当前列更新值，继续计算后续列，直到计算时频矩阵完成。然后从新矩阵每一列，找出最小值其坐标就是时频脊线坐标，进而在时频矩阵中提取最大能量脊。以三阶正弦波矩阵为例详细阐述基于EDM 的时频脊线提取方法，算法步骤如下：

（1）构建以3×3矩阵验证算法有效性（为简化计算将时频矩阵值都设置为负值），该矩阵列代表时间步长，行表示频率，设时频谱图矩阵TFM=-ST如下所示。

（2）以更新ST11为例时，保持第一个时间点的值不变。根据第一列中的值与根据欧式距离计算的惩罚值重构坐标值，以达到寻优和限制频率跳跃的目的。

式中：D代表两个点之间欧氏距离；tip表示当前索引最小值；value表示更新前坐标值；value′表示更新后坐标值；p为惩罚系数；本次以惩罚系数2为例计算。

（3）对最后一列应用惩罚产生最后一列的最小值tip11是8，将最后一列中的最小值与当前索引值相加value′11为10，进而更新第1列中其余元素的值，如下所示：

（4）使用步骤上述中相同的过程，用惩罚因子重新计算前一列的值，获取后序列的新值。最终矩阵为：

（5）从矩阵的最后一列开始，找出最小值，跟踪最小值坐标索引，它们构成了组成时频脊线的路径fi(t)。

在本例中，重构时频矩阵的列最小值坐标为（1，2），（2，2），（3，2），这与所构成矩阵的第2 行的正弦波能量路径相匹配，验证了方案的可行性。

2.2 基于EDM-VKF 的转子时变复杂信号时频表示方法

基于EDM-VKF 的转子时变复杂信号时频表示的构建由以下步骤组成，流程图如图1所示。

图1 基于EDM-VKF的转子时变复杂信号时频表示方法流程图Fig.1 Flowchart of EDM-VKF based time-frequency representation of rotor time-varying complex signals

（1）通过基于经验模态分解的多尺度小波阈值降噪算法处理原始信号，得到降噪后的转子振动信号。

（2）对于基频分量及其谐波分量，通过基于EDM 的时频脊线提取方法从时频矩阵中提取时频脊线并获得其瞬时频率估计值，并使用两个基频分量均值作为转速信号的估计。

（3）在已知谐波分量阶次和转速信号前提下，通过VKF 将基频和谐波阶次分量从原始信号中分离。通过希尔伯特变换得到所有单分量的幅值包络、瞬时相位和瞬时频率。进一步计算各阶谐波分量的瞬时轨道的长轴、短轴、瞬时轨道的正反向分量以及形状指向系数。

（4）将谐波分量轴轨道的正向分量和反向分量幅值投影到时频图三维坐标中，得到谐波分量高分辨率时频表示。所有谐波时频表示的叠加是原始信号的转子时频表示。

3 试验分析

3.1 模拟信号仿真分析

在仿真分析中，对仿真信号的参数进行了精确获知和精确控制，实现本方法揭示的信号特征与预定义的理论信号特征进行比较，检验该方法的性能。在本节中，将通过数值模拟信号分析来说明和测试所提出的方法。第一组模拟信号包含两个非平稳分量，它们模拟停机过程中转子的振动响应，用数学表达式为：

式中：x表示为转子横轴振动信号；y表示为转子纵轴振动信号；信号x基频分量阻尼系数χx被设置为0.6；信号y基频阻尼系数χy被设置为0.5；f1(t) = 160 - 220t+ 110t2、f2(t) = 1.5f1(t)，分别表示两个分量的瞬时频率；x和y的相位差由φi表示，分别设置φ1= π/3、φ2= 4π/3。

x和y具有相同的瞬时频率分量。产生两个合成信号，采样频率为2 048 Hz，采样时间为1 s，在信号中人工加入高斯噪声，信号的信噪比为2.78 dB。噪声信号如图1 所示。由快速傅里叶变换获得的x和y的频域波形如图2（c）和（d）所示，可以看出噪声很强，人们无法从频谱中找到有用的信息。

图2 原始信号时域波形及频域波形Fig.2 Original signal time domain waveform and frequency domain waveform

将原始信号通过经验模态分解算法将振动信号分解为本征模态分量，以x轴信号为例分解结果如图3所示。

图3 信号x的EMD分解图Fig.3 EMD diagram of signal x

通过选取“sym6”的小波基函数和三层小波分解层数，对本征模态分量进行小波分解，采用软阈值方法对信号分量进行处理，得到降噪后的信号，如图4降噪后信号分析图所示。由短时傅里叶变换获得的x1和y1的频谱如图4（c）和（d）所示，可以直观看出经过经验模态分解结合小波阈值降噪，有效滤除了部分高频噪声，信号信噪比提高了一倍同时保留了信号特征，根据皮尔逊系数计算，得到重构信号与无噪声信号相关性系数大于0.91，从时频谱图可以看出，信号由两个分量组成，然而，两个分量从0.4 s 开始出现时频线混叠迹象，难以直观区分，从原始转子轨迹图和全谱也无法获得更多信息。

图4 降噪后信号分析图Fig.4 Signal analysis after noise reduction

在信号时频图谱中通过基于EDM 的时频脊线提取方法提取时频脊线并获取瞬时频率估计值，如图5所示（蓝色线段为拟合值，红色实线为信号真实值），将其作为转速信号应用于VKF提取信号谐波分量，提取结果如图6、图7 所示。通过希尔伯特变换求取幅值包络、瞬时相位和瞬时频率，并构建转子谐波分量的瞬时轴心轨迹。

图7 1.5倍谐波分量时域分析图Fig.7 Time domain analysis diagram of 1.5X harmonic component

真实谐波分量与重构谐波分量对比如图8所示，（左侧为信号x的重构谐波信号，右侧为信号y的重构谐波信号）。为验证方法的性能，通过定义特征平均绝对百分比误差（Mean Absolute Percentage Error，MAPE）和均方根误差（Root Mean Square Error，RMSE）指标进行评价。

图8 真实谐波分量与拟合分量对比图Fig.8 Comparison of the real harmonic component and the fitted component

式中：N为信号长度；下标r表示实际值；下标e为估计值。计算含噪各分量的各项指标，如表1所示。

通过图8 所示，清晰观察出拟合分量与真实分量拟合度较高，表明本文所提出的方法对重构谐波信号估计精度较高，噪声鲁棒性高。

转子振动信号的进动方向系数计算公式为SDI= sin(φxit-φyit)，其中φxit、φyit表示x、y谐波信号相位，估计结果如图9所示，左侧为含噪信号方向系数图，右侧为不含噪信号方向系数图。计算量化指标如表1所示。

图9 进动方向系数图Fig.9 Shape and directivity index

进动方向系数实现了瞬时轨道方向性信息的定量表达，量化指标表示所提出的方法对瞬时轨迹特征估计值非常接近真实值。

最后，将各个谐波分量轴轨道的正向分量和反向分量振幅投影到时频图三维坐标中，得到谐波分量轴轨道时频表示。原始信号的转子信号高分辨率时频表示由两个谐波分量的高分辨率时频表示叠加得到。为验证方法的准确性和有效性，本文将短时傅里叶变换如图10（a）、同步压缩变换的时频表示如图10（b）、理论上的时频表示如图10（c）所示和本文所提出方法的结果如图10（d）进行对比，它具有良好的时频分辨率，C1 和C2分量的时频图可以清晰地分辨出来，与理论时频表示非常接近，并且从时频谱图上可以清晰地观察出各谐波分量的进动方向信息。

图10 信号时频表示对比图Fig.10 Signal time and frequency representation comparison chart

由本方法提出的时频表示结果可以观察出，C1分量反进动分量大于正进动分量，C2 分量正进动分量大于反进动分量，从信号的方向指向系数计算结果（如图9所示）可以证实高分辨率时频图所展示进动分量的正确性。

3.2 水轮机转子不对中信号分析

为进一步验证本文所提出方法的有效性，选取了某旋转机械发生转子不对中时转子振动信号进行分析说明，该信号采集时转子转速为3 180 r/min（53 Hz），信号分析图如图11所示。

图11 转子不对中信号分析图Fig.11 Rotor misalignment signal analysis diagram

从信号频谱分析，信号主要以一倍频和二倍频分量为主，在x的频谱上二倍频分量幅值大于一倍频。从全频谱直接看出信号一倍频分量是反进动，三倍频分量是正进动除此之外无法直观获取更多有价值信息。

接下来，使用本文所提出的算法对上述信号进行分析，所提出的一倍频、二倍频和三倍频谐波分量如图12、图13 和图14所示。从图中可以观察到，谐波分量可以很好的分离出来，三倍频谐波分量是椭圆形，并且信号二倍频谐波分量的轴心轨迹符合转子不对中特征。

图12 一倍频谐波分量时域分析图Fig.12 Time domain analysis diagram of 1X harmonic component

图13 二倍频谐波分量时域分析图Fig.13 Time domain analysis diagram of 2X harmonic component

图14 三倍频谐波分量时域分析图Fig.14 Time domain analysis diagram of 3X harmonic component

计算谐波分量瞬时特征参数，进一步通过计算进动方向系数（如图15 所示），获取谐波分量的进动方向系数，C3 始终保持正进动，C1始终保持负进动。

图15 进动方向系数Fig.15 Shape and directivity index

接下来，构建高分辨率时频谱图［如图16（d）所示］，并给出传统时频表示方法作相比，此图清晰地展示了信号时频谱图，并且提供了基频及谐波分量的进动方向及信号时频域的变化。

图16 信号时频表示对比图Fig.16 Signal time and frequency representation comparison chart

4 结 论

本文提出了引入欧几里得距离的重构时频矩阵的时频脊线提取方法，提高了时频聚集性和时频脊线连贯性，为VKF 提取谐波信号分量提供精准的转子速度信号，并基于谐波分量瞬时特征参数构建转子轴心高分辨率时频表示。该方法实现了转子时变复杂信号的谐波分量提取，同时计算转子谐波信号的瞬时特征参数，可用于描述非平稳过程中转子谐波分量瞬时振动状态。通过仿真实验与转子不对中数据分析，基于EDMVKF 的转子时变复杂信号时频表示方法对信号谐波分量重构精度高，同时高分辨率的时频图提供了转子谐波分量的进动方向信息，清晰地展现转子不对中二倍频信号的特征信息，为旋转机械振动信号分析提供更直观的特征，有利于非平稳过程中旋转机械转子的状态监测。