融合ESPRIT和FFT的动叶片振动参数辨识方法研究

Research on Vibration Parameter Identification Method of Rotating Blade Based on ESPRIT and FFT Algorithm

摘要

叶端定时（Blade Tip Timing， BTT）是当前研究叶片振动问题最有效的测试技术之一， 被广泛应用于动叶片的状态监测与故障诊断。然而， BTT技术固有的非均匀、 欠采样特性引发的叶片振动参数辨识困难已成为当前BTT技术亟需解决的难题。首先， 研究旋转不变子空间（Estimating Signal Parameter via Rotational Invariance Techniques， ESPRIT）和快速傅里叶变换（Fast Fourier Transform， FFT）辨识叶片振动参数的特征； 然后， 对叶片振动的频谱分析结果进行扩展； 最后， 提出一种融合ESPRIT和FFT的动叶片振动参数辨识方法。ESPRIT-FFT方法利用扩展频谱寻找叶片振动频率的估计值， 通过ESPRIT算法确定频率所属的频带， 进而辨识叶片振动参数。仿真研究证明了ESPRIT-FFT方法的可行性和抗噪性。离心压缩机叶轮振动试验显示ESPRIT-FFT方法辨识叶片振动频率的绝对误差为1.4 Hz， 相对误差为0.22%。

Abstract

Currently，blade tip timing（BTT）stands as one of the most effective monitoring techniques for studying blade vibration issues，and it is extensively utilized in condition monitoring and fault diagnosis of rotating blades. However，the difficulties in identifying blade vibration parameters caused by the inherent non-uniformity and under sampling characteristics of the BTT technique have become an urgent challenge that needs to be addressed. This paper initially investigates the attributes of blade vibration parameters identified by the estimating signal parameter via rotational invariance techniques（ESPRIT）and fast Fourier transform（FFT） algorithms. Subsequently， it extended the spectral analysis results of blade vibration，and ultimately proposed a method for identifying vibration parameters of the rotating blade that integrates ESPRIT and FFT algorithms. The ESPRIT-FFT method employed the extended spectrum to estimate the blade vibration frequency， determined the frequency band via the ESPRIT algorithm， and subsequently identified the blade vibration parameters. Simulation studies demonstrate the feasibility and noise resistance of the ESPRIT-FFT method. The vibration test of the centrifugal compressor impeller reveals that the absolute error of the ESPRIT-FFT method in identifying blade vibration frequency is 1.4 Hz，with a relative error of 0.22%.

Graphical abstract

关键词

叶端定时 / 动叶片 / 旋转不变子空间 / 快速傅里叶变换 / 频谱扩展 / 参数辨识

Key words

blade tip timing / rotating blade / estimating signal parameter via rotational invariance techniques / fast Fourier transform / extended spectrum / parameter identification

0 引　言

动叶片是航空发动机、 大型压缩机、 连续风洞等重大装备完成能量转换的核心部件， 其安全运行是重大装备机组健康管理和智能运维的重要组成^［1］。当前， 叶端定时（Blade Tip Timing， BTT）是研究叶片振动问题最有效的测试技术之一， 被广泛应用于动叶片状态监测与故障诊断^［2］。

BTT技术起源于间断相位法和脉冲调制法， 其测试原理为通过在机组上安装的定时传感器连续感知叶片的实际到达时间， 然后， 结合转速同步信号计算叶片的理论到达时间， 最后， 结合两者的时间差获取叶片到达定时传感器的振动位移^［3］。对于单支定时传感器， 转子旋转一圈， 传感器对单个叶片的振动位移实现一次测量， 单支定时传感器的采样频率近似等于转频^［4］。然而， 实际叶片的振动频率远超转子旋转频率。因此， 基于BTT技术获取的叶片振动信号具有严重的欠采样特性^［5］， 且现有的频谱分析方法无法有效辨识叶片的振动参数。

叶片振动参数辨识是通过特定算法对欠采样振动信号进行特征提取， 包括幅值、 相位、 频率等。当前， 国内外学者对叶片振动参数辨识算法进行了广泛探索， 主要分为两类： 时域类算法和空间类算法。时域类算法基于叶片单频或多频响应模型， 利用非线性最小二乘法拟合叶片振动波形， 进而辨识叶片振动参数， 主要包括单参数法^［6］、 双参数法^［7-8］、 自回归法^［9］、 正弦拟合法^［10-11］等， 一般用于解决变转速工况下叶片单频振动响应的参数辨识问题。空间类算法是对叶片振动信号的自相关矩阵进行特征分解， 然后， 利用叶片振动参数与特征值、 特征向量之间的映射关系实现叶片振动参数辨识。空间类算法主要包括多信号分类（Multiple Signal Classification， MUSIC）算法和旋转不变子空间（Estimating Signal Parameter via Rotational Invariance Techniques， ESPRIT）算法。MUSIC算法根据叶片振动信号向量与特征向量之间的正交性对叶片振动频率进行搜索估计^{［12‐14］}， 其通过频率搜索辨识叶片振动参数， 对定时传感器布局有严格的要求。ESPRIT算法是将两支定时传感器构成有延迟和无延迟通道， 根据信号所携带的叶片振动相位差信息， 利用空间类算法解模糊， 得到振动频率的估计值， 该算法通过相关矩阵特征分解能够直接辨识叶片振动频率^［15-16］。ESPRIT算法通过两组信号的相位差辨识叶片振动频率， 当信号信噪比较大时， 相位差辨识精度低， 进而会导致叶片振动参数辨识误差。此外， 基于伪谱估计^［17］、 压缩感知^［18］、 叶片振动主动混叠^［19］等方法也可用于对叶片振动参数进行辨识。

本文在ESPRIT算法原理的基础上， 首先， 研究该算法辨识叶片振动参数的特征， 揭示ESPRIT算法在辨识叶片振动频带方面的优越性； 然后， 对叶片振动欠采样信号进行快速傅里叶变换（Fast Fourier Transform， FFT）， 研究叶片振动频谱的混叠特性， 提出叶片振动频谱扩展方法； 最后， 结合ESPRIT算法和频谱扩展， 提出一种融合ESPRIT和FFT的动叶片振动参数辨识方法（ESPRIT-FFT）， 有效地解决了叶片振动参数辨识精度低的问题， 为基于BTT技术的动叶片状态监测与故障诊断提供理论支撑。

1 ESPRIT算法

1.1 振动信号相关矩阵特征值分解

假设叶片振动信号y\left(t\right)由复正弦信号和白噪声组成， 其振动信号的表达式为

y\left(t\right)={\displaystyle \sum _{k=i}^K}{A}_i{\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_{i}t\mathrm{j}+{\phi }_i}+v\left(t\right),（1）

式中： A_i、 f_i和{\phi }_i分别为叶片振动的幅值、 频率和相位； i为振动分量； \mathrm{j}为复数单位； K为振动分量数量； v\left(t\right)为叶片振动信号噪声。

通过定时传感器连续采集M个数据点， 第q次采样的信号向量可以定义为

\mathit{y}\left(q\right)={\left[\begin{array}{llll}y\left(t_q\right)& y\left(t_{q+\mathrm{1}}\right)& \cdots & y\left(t_{q+M-\mathrm{1}}\right)\end{array}\right]}^{\mathrm{T}}\mathrm{。}（2）

忽略噪声信号v\left(t\right)， 式（2）中复正弦信号可以展开为

\mathit{y}\left(q\right)=\left[\begin{array}{c}{A}_{\mathrm{1}}{\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_{\mathrm{1}}{t}_q\mathrm{j}+{\phi }_{\mathrm{1}}}+A_{\mathrm{2}}{\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_{\mathrm{2}}{t}_q\mathrm{j}+{\phi }_{\mathrm{2}}}+\cdots +A_K{\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_K{t}_q\mathrm{j}+{\phi }_K}\\ A_{\mathrm{1}}{\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_{\mathrm{1}}{t}_{q+\mathrm{1}}\mathrm{j}+{\phi }_{\mathrm{1}}}+A_{\mathrm{2}}{\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_{\mathrm{2}}{t}_{q+\mathrm{1}}\mathrm{j}+{\phi }_{\mathrm{2}}}+\cdots +A_K{e}^{\mathrm{2}\mathrm{\pi }{f}_K{t}_{q+\mathrm{1}}\mathrm{j}+{\phi }_K}\\ ⋮\\ A_{\mathrm{1}}{\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_{\mathrm{1}}{t}_{q+M-\mathrm{1}}\mathrm{j}+{\phi }_{\mathrm{1}}}+A_{\mathrm{2}}{\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_{\mathrm{2}}{t}_{q+M-\mathrm{1}}\mathrm{j}+{\phi }_{\mathrm{2}}}+\cdots +A_K{\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_K{t}_{q+M-\mathrm{1}}\mathrm{j}+{\phi }_K}\end{array}\right]=

\left[\begin{array}{cccc}\mathrm{1}& \mathrm{1}& \cdots & \mathrm{1}\\ {\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_{\mathrm{1}}{\tau }_{\mathrm{1}}\mathrm{j}}& {\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_{\mathrm{2}}{\tau }_{\mathrm{1}}\mathrm{j}}& \cdots & {\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_K{\tau }_{\mathrm{1}}\mathrm{j}}\\ ⋮& ⋮& ⋮& ⋮\\ {\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_{\mathrm{1}}{\tau }_{M-\mathrm{1}}\mathrm{j}}& {\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_{\mathrm{2}}{\tau }_{M-\mathrm{1}}\mathrm{j}}& \cdots & {\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_K{\tau }_{M-\mathrm{1}}\mathrm{j}}\end{array}\right]\left[\begin{array}{c}{A}_{\mathrm{1}}{\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_{\mathrm{1}}{t}_q\mathrm{j}+{\phi }_{\mathrm{1}}}\\ A_{\mathrm{2}}{\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_{\mathrm{2}}{t}_q\mathrm{j}+{\phi }_{\mathrm{2}}}\\ ⋮\\ A_K{\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_K{t}_q\mathrm{j}+{\phi }_K}\end{array}\right],（3）

式中： {\tau }_i为相邻采样点的时间差值， {\tau }_i=t_{q+i}-t_q。

进一步， 式（1）可以写成矩阵形式为

\mathit{y}\left(q\right)=\mathit{A}\mathit{s}\left(q\right)+\mathit{v}\left(q\right),（4）

式中：

\mathit{A}=\left[\begin{array}{cccc}\mathit{a}\left(f_{\mathrm{1}}\right)& \mathit{a}\left(f_{\mathrm{2}}\right)& \cdots & \mathit{a}\left(f_K\right)\end{array}\right]=

\left[\begin{array}{cccc}\mathrm{1}& \mathrm{1}& \cdots & \mathrm{1}\\ {\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_{\mathrm{1}}{\tau }_{\mathrm{1}}\mathrm{j}}& {\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_{\mathrm{2}}{\tau }_{\mathrm{1}}\mathrm{j}}& \cdots & {\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_K{\tau }_{\mathrm{1}}\mathrm{j}}\\ ⋮& ⋮& ⋮& ⋮\\ {\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_{\mathrm{1}}{\tau }_{M-\mathrm{1}}\mathrm{j}}& {\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_{\mathrm{2}}{\tau }_{M-\mathrm{1}}\mathrm{j}}& \cdots & {\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_K{\tau }_{M-\mathrm{1}}\mathrm{j}}\end{array}\right],（5）

其中， 向量\mathit{a}\left(f\right)、 \mathit{s}\left(q\right)和\mathit{v}\left(q\right)的表达式为

\mathit{a}\left(f\right)={\left[\begin{array}{cccc}\mathrm{1}& {\mathrm{e}}^{\mathrm{2}\mathrm{\pi }f{\tau }_{\mathrm{1}}\mathrm{j}}& \cdots & {\mathrm{e}}^{\mathrm{2}\mathrm{\pi }f{\tau }_{M-\mathrm{1}}\mathrm{j}}\end{array}\right]}^{\mathrm{T}},（6）

\mathit{s}\left(q\right)={\left[\begin{array}{cccc}{A}_{\mathrm{1}}{\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_{\mathrm{1}}{t}_q\mathrm{j}+{\phi }_{\mathrm{1}}}& A_{\mathrm{2}}{\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_{\mathrm{2}}{t}_q\mathrm{j}+{\phi }_{\mathrm{2}}}& \cdots & A_K{\mathrm{e}}^{\mathrm{2}\mathrm{\pi }{f}_K{t}_q\mathrm{j}+{\phi }_K}\end{array}\right]}^{\mathrm{T}},（7）

\mathit{v}\left(q\right)={\left[\begin{array}{cccc}v\left(t_q\right)& v\left(t_{q+\mathrm{1}}\right)& \cdots & v\left(t_{q+M-\mathrm{1}}\right)\end{array}\right]}^{\mathrm{T}}\mathrm{。}（8）

叶片振动信号\mathit{y}\left(q\right)的自相关矩阵计算公式为

\mathit{R}=E\left\{\mathit{y}\left(q\right){\mathit{y}}^{\mathrm{H}}\left(q\right)\right\}\mathrm{。}（9）

这里， 假设噪声信号是均值为0， 方差为{\sigma }_v^{\mathrm{2}}的白噪声。因此， 式（9）的自相关矩阵可以表示为

\mathit{R}=\mathit{A}\mathit{P}{\mathit{A}}^{\mathrm{H}}+{\sigma }_v^{\mathrm{2}}\mathit{E}，（10）

式中： \mathit{E}为单位矩阵， \mathit{P}为正定的对角矩阵，

\mathit{P}=diag\left\{\begin{array}{cccc}{\left|A_{\mathrm{1}}\right|}^{\mathrm{2}},& {\left|A_{\mathrm{2}}\right|}^{\mathrm{2}},& \cdots ,& {\left|A_K\right|}^{\mathrm{2}}\end{array}\right\}\mathrm{。}（11）

因此， 可得矩阵\mathit{A}\mathit{P}{\mathit{A}}^{\mathrm{H}}共有K个非零特征值。对\mathit{A}\mathit{P}{\mathit{A}}^{\mathrm{H}}进行特征值分解， 设{\overline{\lambda }}_{\mathrm{1}}、 {\overline{\lambda }}_{\mathrm{2}}、 …、 {\overline{\lambda }}_M为特征值， 则{\mathit{u}}_{\mathrm{1}}、 {\mathit{u}}_{\mathrm{2}}、 …、 {\mathit{u}}_M表示对应的特征向量。相关矩阵 R 共有M个特征值， 对相关矩阵 R 进行特征值分解， 可得特征值的表达式为

\left\{\begin{array}{ll}{\lambda }_i={\overline{\lambda }}_i+{\sigma }_v^{\mathrm{2}},& i=\mathrm{1},\cdots ,K，\\ {\lambda }_i={\sigma }_v^{\mathrm{2}},& i=K+\mathrm{1},\cdots ,M\mathrm{。}\end{array}\right.（12）

自相关矩阵 R 中有K个特征值与信号有关， 其余M-K个特征值仅与噪声有关。根据以上信息， 有以下重要定义：

1） 信号子空间{\mathit{E}}_{\mathrm{s}}： {\mathit{E}}_{\mathrm{s}}是由{\lambda }_{\mathrm{1}}、 {\lambda }_{\mathrm{2}}、 …， {\lambda }_K对应的特征向量{\mathit{u}}_{\mathrm{1}}、 {\mathit{u}}_{\mathrm{2}}、 …、 {\mathit{u}}_K生成的子空间， 记为

{\mathit{E}}_{\mathrm{s}}=span\left\{\begin{array}{cccc}{\mathit{u}}_{\mathrm{1}}& {\mathit{u}}_{\mathrm{2}}& \cdots & {\mathit{u}}_K\end{array}\right\}\mathrm{。}（13）

2） 噪声子空间{\mathit{E}}_{\mathrm{N}}： {\mathit{E}}_{\mathrm{N}}是由{\lambda }_{K+\mathrm{1}}、 …、 {\lambda }_M对应的特征向量{\mathit{u}}_{K+\mathrm{1}}、 …、 {\mathit{u}}_M生成的子空间， 记为

{\mathit{E}}_{\mathrm{N}}=span\left\{\begin{array}{cccc}{\mathit{u}}_{K+\mathrm{1}}& {\mathit{u}}_{K+\mathrm{2}}& \cdots & {\mathit{u}}_M\end{array}\right\}\mathrm{。}（14）

1.2 ESPRIT算法原理

在实际试验过程中， 两支定时传感器可以构成有延迟的双通道。ESPRIT算法根据两组叶片振动信号所携带的相位差信息， 通过频率解模糊， 对叶片振动频率进行估计。根据式（2）， 两支定时传感器第q次采样的叶片振动信号可以定义为

{\mathit{y}}_{\mathrm{1}}\left(q\right)={\left[\begin{array}{cccc}{y}_{\mathrm{1}}\left(t_q\right)& y_{\mathrm{1}}\left(t_{q+\mathrm{1}}\right)& \cdots & y_{\mathrm{1}}\left(t_{q+M-\mathrm{1}}\right)\end{array}\right]}^{\mathrm{T}}，（15）

{\mathit{y}}_{\mathrm{2}}\left(q\right)={\left[\begin{array}{cccc}{y}_{\mathrm{2}}\left(t_q\right)& y_{\mathrm{2}}\left(t_{q+\mathrm{1}}\right)& \cdots & y_{\mathrm{2}}\left(t_{q+M-\mathrm{1}}\right)\end{array}\right]}^{\mathrm{T}}\mathrm{。}（16）

这里， 假设两支定时传感器之间的安装夹角为\mathrm{\Delta }\beta， 因此， 两支定时传感器通道之间的延迟时间为

{\tau }^{†}=\frac{\mathrm{\Delta }\beta }{\mathrm{2}\mathrm{\pi }{f}_s}，（17）

式中： f_s为转轴旋转频率， 单位为Hz。

进一步， 式（16）可以改写为

{\mathit{y}}_{\mathrm{2}}\left(q\right)=\begin{array}{cc}\left[y_{\mathrm{1}}\left(t_q+{\tau }^{†}\right)\right.& y_{\mathrm{1}}\left(t_{q+\mathrm{1}}+{\tau }^{†}\right)\end{array}

{\left.\cdots y_{\mathrm{1}}\left(t_{q+M-\mathrm{1}}+{\tau }^{†}\right)\right]}^{\mathrm{T}}\mathrm{。}（18）

根据式（4）， 两支定时传感器采集的振动信号可以表示为

{\mathit{y}}_{\mathrm{1}}\left(q\right)=\mathit{A}\mathit{s}\left(q\right)+{\mathit{v}}_{\mathrm{1}}\left(q\right)，（19）

{\mathit{y}}_{\mathrm{2}}\left(q\right)=\mathit{A}{\mathit{\Lambda }}^{†}\mathit{s}\left(q\right)+{\mathit{v}}_{\mathrm{2}}\left(q\right)，（20）

式中：

{\mathit{\Lambda }}^{†}=\left[\begin{array}{cccc}{\mathrm{e}}^{-\mathrm{2}\mathrm{\pi }{f}_{\mathrm{1}}{\tau }^{†}\mathrm{j}}& & & \\ & {\mathrm{e}}^{-\mathrm{2}\mathrm{\pi }{f}_{\mathrm{2}}{\tau }^{†}\mathrm{j}}& & \\ & & \ddots & \\ & & & {\mathrm{e}}^{-\mathrm{2}\mathrm{\pi }{f}_K{\tau }^{†}\mathrm{j}}\end{array}\right]\mathrm{。}（21）

将向量{\mathit{y}}_{\mathrm{1}}\left(q\right)和{\mathit{y}}_{\mathrm{2}}\left(q\right)合并成一个向量\mathit{Y}\left(q\right)， 可以表示为

\mathit{Y}\left(q\right)=\left[\begin{array}{c}{\mathit{y}}_{\mathrm{1}}\left(q\right)\\ {\mathit{y}}_{\mathrm{2}}\left(q\right)\end{array}\right]=\left[\begin{array}{c}\mathit{A}\\ \mathit{A}{\mathit{\Lambda }}^{†}\end{array}\right]\mathit{s}\left(q\right)+\left[\begin{array}{c}{\mathit{v}}_{\mathrm{1}}\left(q\right)\\ {\mathit{v}}_{\mathrm{2}}\left(q\right)\end{array}\right]\mathrm{。}（22）

向量\mathit{Y}\left(q\right)的自相关矩阵计算公式为

\mathit{R}=E\left\{\mathit{Y}\left(q\right){\mathit{Y}}^{\mathrm{H}}\left(q\right)\right\}=\overline{\mathit{A}}\mathit{P}{\overline{\mathit{A}}}^{\mathrm{H}}+{\sigma }_v^{\mathrm{2}}\mathit{E}\mathrm{。}（23）

然后， 对自相关矩阵 R 进行特征值分解， {\mathit{E}}_{\mathrm{s}}为前K个大特征值对应的特征向量组成的信号子空间， {\mathit{E}}_{\mathrm{s}}=span\left\{{\mathit{u}}_{\mathrm{1}},{\mathit{u}}_{\mathrm{2}},\cdots ,{\mathit{u}}_K\right\}， {\mathit{E}}_{\mathrm{N}}为2M-K个较小特征值对应的特征向量组成的噪声子空间， {\mathit{E}}_{\mathrm{N}}=span\left\{{\mathit{u}}_{K+\mathrm{1}},{\mathit{u}}_{K+\mathrm{2}},\cdots ,{\mathit{u}}_{\mathrm{2}M}\right\}。这时， 存在唯一的、 非奇异的M\times M维满秩矩阵 T 使下式成立

{\mathit{E}}_{\mathrm{s}}=\overline{\mathit{A}}\mathit{T}\mathrm{。}（24）

根据矩阵的时移不变特性， 可以将{\mathit{E}}_{\mathrm{s}}分解为两个矩阵， 即{\mathit{E}}_{y_{\mathrm{1}}}\in {\mathbf{C}}^{M\times K}和{\mathit{E}}_{y_{\mathrm{2}}}\in {\mathbf{C}}^{M\times K}， 可以得到如下关系式

{\mathit{E}}_{\mathrm{s}}=\left[\begin{array}{c}{\mathit{E}}_{y_{\mathrm{1}}}\\ {\mathit{E}}_{y_{\mathrm{2}}}\end{array}\right]=\left[\begin{array}{c}\mathit{A}\mathit{T}\\ \mathit{A}{\mathit{\Lambda }}^{†}\mathit{T}\end{array}\right]\mathrm{。}（25）

根据式（25）， 可以得到

{\mathit{T}}^{-\mathrm{1}}{\mathit{\Lambda }}^{†}\mathit{T}={\left({\mathit{E}}_{y_{\mathrm{1}}}^{\mathrm{H}}{\mathit{E}}_{y_{\mathrm{1}}}\right)}^{-\mathrm{1}}{\mathit{E}}_{y_{\mathrm{1}}}^{\mathrm{H}}{\mathit{E}}_{y_{\mathrm{2}}}，（26）

式中： {\left({\mathit{E}}_{y_{\mathrm{1}}}^{\mathrm{H}}{\mathit{E}}_{y_{\mathrm{1}}}\right)}^{-\mathrm{1}}{\mathit{E}}_{y_{\mathrm{1}}}^{\mathrm{H}}{\mathit{E}}_{y_{\mathrm{2}}}的特征值为{\mathit{\Lambda }}^{†}的对角元素。对{\left({\mathit{E}}_{y_{\mathrm{1}}}^{\mathrm{H}}{\mathit{E}}_{y_{\mathrm{1}}}\right)}^{-\mathrm{1}}{\mathit{E}}_{y_{\mathrm{1}}}^{\mathrm{H}}{\mathit{E}}_{y_{\mathrm{2}}}进行特征值分解， 通过计算特征值{\lambda }_{\mathrm{1}}、 {\lambda }_{\mathrm{2}}、 …、 {\lambda }_K携带的信号延迟相位差得到振动频率估计值为

P_{\mathrm{e}\mathrm{s}\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{t}}\left(k\right)=\left|\frac{\mathrm{\angle }{\lambda }_k}{\mathrm{2}\mathrm{\pi }{\tau }^{†}}\right|，（27）

式中： \mathrm{\angle }\left(\cdot \right)表示取相位。

ESPRIT算法根据两支传感器通道所携带的相位差信息， 通过两次特征值分解， 对叶片振动频率进行辨识。然而， 在实际工程应用中， ESPRIT算法对叶片振动信号的噪声比较敏感， 通过ESPRIT算法对BTT获取的欠采样信号进行频率辨识存在较大误差。

1.3 ESPRIT算法仿真分析

基于MATLAB软件模拟叶片振动信号， 并通过仿真研究ESPRIT算法的误差特征。模拟叶片单频异步振动， 振动幅值为1 mm， 振动频率为378 Hz。模拟叶片振动全采样， 根据香农采样定理， 全采样的采样频率设置为1 000 Hz。模拟基于BTT技术的叶片振动测试， 仿真2支BTT传感器， 分别命名为Tip1和Tip2， 传感器安装角度分别设置为0°和30°。模拟叶片转子的旋转频率为100 Hz， 即BTT技术单探针的采样频率为100 Hz。叶片振动采样时长为100 s。模拟叶片振动测试过程中的非确定性因素， 在叶片振动信号中添加白噪声， 信噪比SNR=10 dB。图 1 展示了0.2 s的叶片振动信号， 图中折线表示叶片振动全采样信号， 散点图表示叶片振动欠采样信号。

对叶片振动全采样信号和BTT欠采样信号进行FFT分析， 不同信号的频谱分析结果如图 2 所示。可以发现： 1） 全采样信号的频谱分析结果为378 Hz和622 Hz， 其中622 Hz是由于FFT双边谱特性导致的； 2） BTT欠采样信号的频谱分析结果为22 Hz和78 Hz， 这是由于在欠采样状态下， 频率混叠造成的。BTT欠采样信号的频谱分析结果与仿真假设完全不同。

利用ESPRIT算法对叶片振动信号进行参数辨识。由上述内容可知， 仿真信号的数据长度为1×10⁴， 设置每次连续采集的数据点长度M=100， 因此， 根据仿真信号可以对叶片振动频率重复估计100次。为了将ESPRIT算法的结果进行图形可视化， 对每次频率估计结果进行赋值， 其计算公式为

\begin{array}{cc}Q\left(i\right)=\frac{F_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\mathrm{100}},& i=\mathrm{1,2},\cdots ,\mathrm{100}\end{array}，（28）

式中： F_{\mathrm{m}\mathrm{a}\mathrm{x}}为扩展频谱幅值的最大值； Q\left(i\right)为无量纲参数； i为估计的次数。

接着， 将每一次频率估计结果作为横坐标， 对应的Q\left(i\right)作为纵坐标绘制ESPRIT算法结果的散点图， 如图 3 所示。

图 3 显示， ESPRIT算法获取的叶片振动频率分布在实际频率378 Hz附近， 平均值为373.83 Hz， 频率辨识绝对误差和相对误差分别为4.17 Hz和1.10%。

2 ESPRIT⁃FFT方法

2.1 基于ESPRIT增强FFT的方法

叶片振动模拟频率为378 Hz， BTT采样频率为100 Hz， 因此， 可以得到叶片振动的混叠频率为

Mod\left(\mathrm{378,100}\right)=\mathrm{78}，（29）

式中： Mod\left(\cdot \right)为求余函数。

根据图 2（b）， 对BTT欠采样振动信号进行FFT分析， 虽然无法准确辨识叶片振动频率， 但可以准确辨识叶片振动的混叠频率。另外， 图3显示ESPRIT算法的估计结果分布在真实频率附近， 辨识精度较低， 但ESPRIT算法可以准确估计叶片振动频率所在的频带。因此， 可以通过将FFT结果进行倍频扩展， 然后， 结合ESPRIT算法的结果， 实现叶片振动频率的准确辨识。本文将该方法称为融合ESPRIT和FFT的动叶片振动参数辨识方法， 简称为ESPRIT-FFT方法。

2.2 ESPRIT⁃FFT方法仿真分析

首先， 将图2（b）的2组叶片振动频谱进行平均处理。然后， 将处理后的频谱按照8倍频进行扩展。接着， 将ESPRIT算法的估计结果绘制到扩展后的频谱。ESPRIT-FFT方法的辨识结果如图 4 所示， 可以发现ESPRIT-FFT方法辨识结果与仿真假设完全一致， 能够准确辨识叶片振动频率。

接着， 进一步考察ESPRIT-FFT方法的抗噪性。将叶片振动信号中的白噪声设置为SNR=-5 dB。利用ESPRIT-FFT方法辨识叶片振动频率， 结果如图 5 所示。图 5 显示ESPRIT算法的辨识结果方差较大， 辨识频率的平均值为401.19 Hz， 绝对误差达到23.19 Hz， 与仿真值偏差较大。

为了更准确地辨识叶片振动频率， 本文以1/2倍采样频率为基准， 对扩展后的频谱进行频段划分， 然后， 统计不同频段内ESPRIT算法辨识结果所占的比例， 结果如图 6 所示。图 6 中黄色矩形条表示ESPRIT算法辨识结果所占的比例， 根据条形图可以准确定位叶片振动频率所在的频带， 进而可以确定叶片振动频率为378 Hz， 与仿真假设完全一致。

ESPRIT-FFT方法结合了ESPRIT算法和FFT算法频率辨识的优势， 通过扩展FFT算法的频谱预测叶片振动可能出现的频率， 然后， 利用ESPRIT算法确定叶片振动频率的频带， 最后， 确定叶片振动频率。基于ESPRIT-FFT方法的叶片振动频率辨识方法如图 7 所示。仿真结果显示ESPRIT-FFT方法的辨识结果与仿真假设完全一致， 可以证明ESPRIT-FFT方法的可行性和可靠性。

3 离心压缩机试验及结果分析

3.1 离心压缩机试验台

在沈鼓集团Φ800型号离心压缩机模型级试验台上开展基于BTT技术的叶轮叶片振动试验。离心压缩机试验台主要由进口导叶、 离心叶轮、 扩压器、 回流通道等组成， 离心压缩机及叶轮实物图见图 8 所示。离心压缩机叶轮额定转速可达9 000 r/min。

本文利用BTT技术和应变片法开展离心压缩机叶轮叶片振动测试与分析。在离心压缩机进口的型环上安装2支BTT传感器， 分别命名为Tip1和Tip2， 2支传感器的夹角为10°。同时， 在叶轮叶片表面粘贴应变片， 对叶片振动应力进行测试， 用于验证本文研究方法在工业应用中的准确性。离心压缩机叶轮振动测试传感器布局如图 9 所示。

离心压缩机叶轮叶片振动测试系统包括应变测试系统和BTT测试系统。应变信号传输采用德国MANNER公司的遥测系统， 遥测系统的转子部件固定在压缩机的主轴端， 静子部件固定在压缩机进口管路中。应变信号经NI-9231采集卡将模拟信号传输到计算机端的DASP采集系统， 进行信号的采集和分析， 应变信号的采样频率设置为10 240 Hz。BTT测试系统采用美国HOOD公司的激光发生器， 当叶片经过定时传感器时， 定时传感器生成方波信号并发送到PCI-6602采集卡， 进而记录试验过程中叶片的BTT信号。离心压缩机叶轮叶片振动测试系统如图 10 所示。

3.2 离心压缩机试验结果分析

离心压缩机在恒速工况下运行， 压缩机运行转速为5 628 r/min， 待转速稳定后， 开始叶轮叶片振动测试。基于实测的BTT信号计算叶片振动位移， 其中7号叶片振动测量结果如图 11 所示。对叶片的振动信号进行FFT分析， 7号叶片振动的频谱计算结果如图 12 所示。

根据ESPRIT-FFT方法的原理， 首先， 将图 12 的频谱进行10倍扩展； 然后， 基于叶片振动测量结果， 建立叶片振动自相关矩阵， 利用ESPRIT算法对叶片振动频率进行辨识， 并统计各频带内辨识结果的占比； 最后， 绘制ESPRIT-FFT方法的频率辨识结果， 如图 13 所示， 可以发现， 离心压缩机叶轮7号叶片振动频率为635.3 Hz。

为了进一步验证ESPRIT-FFT方法的准确性， 本文对7号叶片的应变信号进行分析， 取应变信号数据长度为10 s， 应变信号的时域图如图 14（a） 所示， 并对应变信号进行频谱分析， 结果如图 14（b） 所示， 可以发现叶片振动频率为633.9 Hz。与应变片测试结果相比， ESPRIT-FFT方法叶片振动频率的绝对误差为1.4 Hz， 相对误差为0.22%。因此， 可以证明ESPRIT-FFT方法在实际工程应用中的可靠性。

4 结　论

本文通过理论研究和仿真分析揭示了ESPRIT和FFT算法辨识叶片振动参数的优点和缺陷。ESPRIT算法的频率辨识精度低， 但能够准确辨识叶片振动频率的频带， 通过对FFT的频谱进行扩展能够有效获取叶片振动频率的估计值。然后， 本文结合两种算法的优点， 提出一种融合ESPRIT和FFT的动叶片振动参数辨识方法。仿真分析证明了ESPRIT-FFT方法的可行性和抗噪性。在离心压缩机试验台开展叶轮叶片振动试验， 与应变片法的结果相比， ESPRIT-FFT方法频率辨识的绝对误差为1.4 Hz， 相对误差为0.22%， 进一步证明了ESPRIT-FFT方法在实际工程中应用的价值。然而， ESPRIT-FFT方法只适用于解决频带内叶片单频振动的问题， 对于动叶片多频振动参数辨识， 还需在未来的工作中继续研究。

参考文献

[1] 贺长波. 离心式压缩机叶片状态在线监测与故障预警方法研究［D］. 大连： 大连理工大学， 2019.

[2] 欧阳涛. 基于叶尖定时的旋转叶片振动检测及参数辨识技术［D］. 天津： 天津大学， 2011.

[3] FAN Z， LI H， DONG J， et al. An improved multiple per revolution-based blade tip timing method and its applications on large-scale compressor blades［J］. Mechanical Systems and Signal Processing， 2022， 167： 108538.

[4] DONG J， LI H， FAN Z， et al. Time-frequency sparse reconstruction of non-uniform sampling for non-stationary signal［J］. IEEE Transactions on Vehicular Technology， 2021， 70（11）： 11145-11153.

[5] CHEN Z， SHENG H， XIA Y， et al. A comprehensive review on blade tip timing-based health monitoring： status and future［J］. Mechanical Systems and Signal Processing， 2021， 149： 107330.

[6] 郭浩天， 段发阶， 汪猛. 基于叶尖定时单参数法的叶片共振倍频数辨识［J］. 天津大学学报（自然科学与工程技术版）， 2016， 49（9）： 951-955.

[7] GUO Haotian， DUAN Fajie， WANG Meng. Engine order identification based on blade tip-timing single parameter method［J］. Journal of Tianjin University （Science and Technology）， 2016， 49（9）： 951-955. （in Chinese）

[8] HEATH S. A new technique for identifying synchronous resonances using tip-timing［J］. Journal of Engineering for Gas Turbines and Power， 2000， 122（2）： 219-225.

[9] LI H， FAN Z， DONG J， et al. An improved blade vibration difference-based two-parameter plot method for synchronous vibration parameter identification of rotating blades［J］. Measurement， 2023， 207： 112397.

[10] CARRINGTON I B， WRIGHT J R， COOPER J E， et al. A comparison of blade tip timing data analysis methods［J］. Proceedings of the Institution of Mechanical Engineers， Part G： Journal of Aerospace Engineering， 2001， 215（5）： 301-312.

[11] KAŹMIERCZAK K， PRZYSOWA R. Standard sine fitting algorithms applied to blade tip timing data［J］. Journal of KONBiN， 2014， 30（1）： 21-30.

[12] FAN Z， LI H， HUANG J， et al. Blade vibration difference-based circumferential Fourier fitting algorithm for synchronous vibration parameter identification of rotation blades［J］. Sensors， 2024， 24（24）： 8083.

[13] WANG Z， YANG Z， WU S， et al. An improved multiple signal classification for nonuniform sampling in blade tip timing［J］. IEEE Transactions on Instrumentation and Measurement， 2020， 69（10）： 7941-7952.

[14] WANG Z， YANG Z， LI H， et al. Subspace dimension reduction for faster multiple signal classification in blade tip timing［J］. IEEE Transactions on Instrumentation and Measurement， 2021， 70： 3509510.

[15] LIU Z， DUAN F， NIU G， et al. Reconstruction of blade tip-timing signals based on the MUSIC algorithm［J］. Mechanical Systems and Signal Processing， 2022， 163： 108137.

[16] 李孟麟， 段发阶， 欧阳涛， 等. 基于叶尖定时的旋转机械叶片振动频率辨识ESPRIT方法［J］. 振动与冲击， 2010， 29（12）： 18-21.

[17] LI Menglin， DUAN Fajie， OUYANG Tao， et al. ESPRIT method for blade vibration frequency identification in rotating machinery based on blade tip-timing measurement［J］. Journal of Vibration and Shock， 2010， 29（12）： 18-21. （in Chinese）

[18] 贺长波， 李宏坤， 赵新维， 等. 基于总体最小二乘准则旋转不变子空间法的叶尖定时欠采样信号分析［J］. 机械工程学报， 2019， 55（19）： 103-111.

[19] HE Changbo， LI Hongkun， ZHAO Xinwei， et al. Analysis method for under-sampled blade tip-timing signal based on the rotational invariance technique with total least squares principle［J］. Journal of Mechanical Engineering， 2019， 55（19）： 103-111. （in Chinese）

[20] VERCOUTTER A， LARDIES J， BERTHILLIER M， et al. Improvement of compressor blade vibrations spectral analysis from tip timing data： aliasing reduction［C］//ASME Turbo Expo 2013： Turbine Technical Conference and Exposition， 2013： 1-8.

[21] DONG J， LI H， FAN Z， et al. Characteristics analysis of blade tip timing signals in synchronous resonance and frequency recovery based on subspace pursuit algorithm［J］. Mechanical Systems and Signal Processing， 2023， 183： 109632.

[22] CAO J， YANG Z， TIAN S， et al. Biprobes blade tip timing method for frequency identification based on active aliasing time-delay estimation and dealiasing［J］. IEEE Transactions on Industrial Electronics， 2023， 70（2）： 1939-1948.

基金资助

国家自然科学基金资助项目(52405141)

山西省基础研究计划资助项目(202403021212143)

山西省高等学校科技创新资助项目(2024L177)

中国博士后科学基金资助项目(2024M752992)