CN114111997A - 一种基于叶端定时欠采样信号特性的叶片同步共振频率恢复方法 - Google Patents

Abstract

一种基于叶端定时欠采样信号特性的叶片同步共振频率恢复方法，步骤如下：1)随机布局安装叶端定时传感器，获取叶端定时信号的振动位移，对每个通道进行离散傅里叶变换形成低维矩阵；2)测量叶端定时传感器之间的相对角度，设置满足奈奎斯特采样定理需要安装的传感器个数，计算每个叶端定时传感器对应的位置参数；3)计算低维矩阵的补偿矩阵以及参数化的压缩矩阵，构成多个列向量的欠定方程组；4)迭代求解过程中选择多个原子计算最小二乘解的稀疏估计值，获得重构的高维矩阵；5)将高维矩阵的行向量按照顺序首尾相连为恢复的频谱值，获得叶片振动的恢复频率。通过与迭代过程选个单原子的方法相比，本发明恢复叶片振动频率的准确性更高。

Description

一种基于叶端定时欠采样信号特性的叶片同步共振频率恢复 方法

技术领域

本发明属于旋转机械叶片振动信号的在线监测领域，涉及一种基于叶端定时欠采样信号特性的叶片同步共振频率恢复方法。

背景技术

叶片广泛应用于压气机、燃气轮机等涡轮机械中，是设备运行中实现机组内能量转换的关键部件。目前叶片的失效模式很多，主要原因是高周疲劳振动引起的疲劳失效。因此，需要对叶片振动进行实时在线监测，以了解叶片的健康状态。准确识别叶片的固有频率对各种设备的安全运行具有重要意义。叶端定时法是非接触式方法测量叶片振动的研究热点，已经称为一种很有前景的叶片振动监测方法。

叶端定时法的采样频率取决于叶轮转速和传感器的个数，但由于空间约束，传感器的安装个数有限，叶片的共振频率远高于采样频率。因此，叶端定时信号在频率识别上存在严重的欠采样现象，不能真实有效地反应叶片的振动信息，需要从这些欠采样的信号中恢复出原始的叶片振动频率。已有的叶端定时信号的重构方法是基于迭代过程选择最优原子求解最小二乘问题的恢复算法。但由于叶片同步共振频率的特殊性，迭代选择单原子的方法会造成叶片振动频率的失效识别。因此需要一种基于叶端定时欠采样信号的高效识别叶片振动频率的方法，对旋转机械的叶片振动监测有着重要的意义。

发明内容

本发明针对叶端定时欠采样信号的混叠成分的特殊性质，提出了一种叶片同步共振频率的恢复方法。该方法在重构过程中寻找压缩矩阵的最佳原子时，同时选择多个原子，避免叶片同步共振频率贡献相同的两个频谱值因迭代次序不同而使其中一个被削弱的问题，有效识别叶片的振动频率。本发明对选择单原子的重构方法做了改进，考虑重构过程中单原子选择的计算误差，能够有效提高叶片振动频率恢复的准确性。

为了达到上述目的，本发明采用的技术方案为：

一种基于叶端定时欠采样信号特性的叶片同步共振频率恢复方法，包括以下步骤：

第一步：随机安装P个叶端定时传感器(按照叶轮旋转方向，随机选择一个叶端定时传感器标记顺序从1到P)在叶轮周围的圆周型壳体上，用来获取叶片的实际到达时间；将同步转速参考传感器固定在轴的附近，用来获取叶片的理想到达时间；计算第p个叶端定时传感器获得的第n圈的叶端定时振动位移信号：

s_p[n]＝2πf_shaftR(t_p，n-t′_p，n) (1)

其中，p∈{1，2，...P}，f_shaft为叶轮的旋转频率，R为叶尖到叶轮中心的距离，t_p，n为通过叶端定时传感器获取的叶片实际到达时间，t′_p，n为通过同步转速参考传感器获取的叶片理想到达时间。

通过公式(1)计算每个叶端定时传感器每圈的振动位移后，采用公式(2)对每个叶端定时传感器的振动位移序列进行离散傅里叶变换：

其中，k_P∈{0，1，...，N-1}，N表示叶轮旋转的圈数，也就是每个叶端定时传感器的采样点数为N。

将每个叶端定时传感器计算得到的离散傅里叶变换结果作为行向量形成矩阵M_BTT：

第二步：根据实际安装的叶端定时传感器之间的相对角度，计算每个叶端定时传感器的位置在频谱恢复过程中的参数值；记录叶端定时传感器之间的相对角度集合α：

α＝{α_(1，1)，…，α_(1，p)，…α_(1，P)} (4)

其中，α_(1，p)表示第p个叶端定时传感器和第一个叶端定时传感器之间的相对角度。

在理想状态下，为了满足奈奎斯特采样定理，假设需要在叶轮周围圆周型壳体上均匀安装的叶端定时传感器个数为Q，定义每个叶端定时传感器的位置序号集合为C＝{0，1，...，Q-1}。每相邻的两个叶端定时传感器之间的夹角为：

α_Q＝360°/Q (5)

根据相对角度集合α中叶端定时传感器之间的相对角度，采用公式(6)计算实际安装的每个叶端定时传感器的相对应的位置参数：

c_p＝round(α_(1，p)/α_Q) (6)

其中，函数round为四舍五入取整函数。

因α_(1，1)＝0°，最终得到实际安装的叶端定时传感器的位置参数集合：

C_P＝{0，...c_p，...c_P} (7)

第三步：计算压缩矩阵Φ和矩阵M_BTT的补偿矩阵B；定义构成压缩矩阵Φ和补偿矩阵B的元素分别为：

其中，p∈{1，2，...P}，q∈{1，2，...Q}，k_P∈{0，1，...，N-1}，P表示实际安装的叶端定时传感器个数，Q表示为满足奈奎斯特采样定理假设安装的叶端定时传感器个数，N表示叶轮旋转的圈数，c_p表示实际安装的叶端定时传感器的位置参数，

表示压缩矩阵Φ中第p行第q列元素，

表示补偿矩阵B中第p行第(k_P+1)列元素。

由公式(8)、(9)可知，压缩矩阵Φ和补偿矩阵B的维度分别为P×Q和P×N。

第四步：利用最小二乘法重构叶片振动频率。假设安装Q个叶端定时传感器，叶轮旋转N圈后，获取的整个振动位移序列为s[m]，而对于整个振动位移序列的离散傅里叶变换结果为：

其中，M＝Q·N，M表示假设安装Q个叶端定时传感器获取的整个振动位移序列的点数，k_Q∈{0，1，...，M-1}。

由于是假设安装Q个叶端定时传感器来满足奈奎斯特采样定理，因此S[k_Q]是未知的奈奎斯特采样序列的离散傅里叶变换结果。在已知的叶端定时采样序列的离散傅里叶变换结果S_p[k_P]和未知的奈奎斯特采样序列的离散傅里叶变换结果S[k_Q]之间存在如下关系：

其中，(q-1)N+k_P＝k_Q。

将M＝Q·N点未知的奈奎斯特采样序列的离散傅里叶变换结果S[k_Q]沿频率方向分割Q段，每段N点，每段作为行向量形成矩阵M_NYQ：

因此，公式(11)可以转化为矩阵的表达形式：

M＝ΦM_NYQ (13)

其中，

运算符号

为Hadamard积(矩阵对应元素相乘)，矩阵M_NYQ的维度为Q×N。

在工程应用中，由于压缩机组的空间约束，叶端定时传感器的安装个数受限，所以实际安装的叶端定时传感器个数P小于需要满足奈奎斯特采样定理的假设安装的叶端定时传感器个数Q。因此，矩阵M和矩阵M_NYQ中相对应的列向量是一种典型的欠定问题。可根据矩阵M中的低维列向量和矩阵M_NYQ中高维列向量的稀疏度K，采用最小二乘法经K次迭代获得M_NYQ中相对应的高维列向量的稀疏估计值，形成重构的矩阵M′_NYQ。

对于整个迭代过程：输入端为矩阵M中已知的低维列向量、压缩矩阵Φ和高维列向量的稀疏度K；输出端为矩阵M_NYQ中未知的高维列向量的稀疏估计值。

在每次迭代计算过程中，利用上一次迭代获得的残差向量(第一次迭代计算时的残差向量即为矩阵M中已知的低维列向量)与压缩矩阵Φ中的原子(矩阵的列向量)(需去除上一次迭代选择的原子)的内积值的模的最大值来寻找最佳原子，并根据最佳原子的位置确定高维向量中稀疏估计值的位置，可表示为公式(14)：

其中，r_t-1为第(t-1)次迭代获得的残差向量，r₀为M中已知的低维列向量，

为压缩矩阵Φ的第q个原子。

叶片的同步共振频率为叶轮旋转频率的整数倍δ，在叶端定时信号上会发生混叠效应，表现在频率为0 Hz上，即出现在S_p[k_P＝0]的位置上，此时可以得到β_p，0＝1。在进行第一次迭代时，r₀＝[QS₁[0]，...，QS_p[0]，...，QS_P[0]]^T，可以得到：

对于两个原子

和

有如下结果：

对比公式(16)和(17)，可以得到两个内积值的模相等的最佳原子：

λ_(1，δ+1)＝λ_(1，Q-δ+1) (18)

传统的计算方法是每次迭代选择单个原子放入最优原子集合，经K次迭代可重构获得高维向量中K个稀疏估计值。基于公式(18)的性质，本发明每次迭代同时选择K个原子放入最优原子集合，再经迭代过程获取稀疏估计值。叶片同步共振频率的离散傅里叶变换结果形成的两个频谱值位于同一个列向量，所以本方法的最小稀疏度为K＝2，需要同时选择2个原子放入最优集合。

需要注意的是：本发明方法针对的是叶片的同步共振频率，主要表现在矩阵M的第一个列向量，即r₀＝[QS₁[0]，...，QS_p[0]，...，QS_P[0]]^T，而引起公式(18)的主要原因主要是由于k_P＝0，因此，主要描述矩阵M的第一个列向量在第一次迭代计算过程中的相关问题。

对于矩阵M的第一个列向量[QS₁[0]，...，QS_p[0]，...，QS_P[0]]^T的迭代计算过程：输入端为r₀＝[QS₁[0]，...，QS_p[0]，...，QS_P[0]]^T、压缩矩阵Φ和高维列向量的稀疏度K＝2。

第一次迭代时，有λ_(1，δ+1)＝λ_(1，Q-δ+1)，第(δ+1)和第(Q-δ+1)原子被同时选入最优原子集合：

求解r₀＝Φ₁θ₁的最小二乘解：

然后计算迭代后的残差向量r₁：

第二次迭代计算的输入向量为第一次迭代后得到的残差向量r₁，计算过程与第一次迭代计算过程相同(第二次迭代中选择原子的过程需去除第一次迭代选择的第(δ+1)和第(Q-δ+1)原子)。经过2次迭代可获得作为输出端的稀疏估计值

对应最优原子集合中原子的位置获得重构的高维列向量的稀疏估计值。

矩阵M的其他列向量的计算过程与第一个列向量相同，获得重构的矩阵M′_NYQ。

以上所述的第一步至第四步是针对叶片同步共振下的叶端定时欠采样信号的特殊性质，实现了叶片同步共振频率的恢复。

本发明的有益效果：

(1)本发明属于欠采样信号的频率恢复方法。叶片同步共振的叶端定时信号具有特殊的性质，采用低维向量与压缩矩阵原子的内积值的模作为度量寻找最佳原子时，会出现两个值相等的情况。在其他的欠定方程的恢复方法中，每次迭代只选择单个最优原子。

(2)本发明同时选择多个原子，避免了选中一个原子后因更新残差使另一个原子的贡献被忽略的情况，能够显著提高叶端定时欠采样信号频率恢复的准确性。

附图说明

图1为应变信号的频谱分析；

图2为叶端定时法计算叶端定时信号的振动位移；

图3为不同传感器的单通道的傅里叶变换结果；

图4为第一次迭代时不同原子的内积值的模；

图5为选择单个原子后的第二次迭代的内积值的模；

图6为迭代过程选择单个原子的恢复频谱图；

图7为本发明方法的恢复频谱图；

图8为叶端定时欠采样信号的频率恢复方法的流程图。

具体实施方式

以下为结合技术方案和附图详细叙述本发明的实施例。实施例的数据来自于叶端定时测振试验台。该试验台主要由驱动电机、叶轮和带有12个激励磁铁的保护罩组成，叶轮直接安装在电机的输出轴上，通过变频器编程可改变旋转试验台的运行转速。叶端定时测振系统主要由光纤传感器、光电前放、信号调理模块和计数器组成。

首先确定叶片在激励磁铁作用下发生同步共振时的转速。采用应变片法获取叶片的振动频率，用来对比试验结果。应变信号的采样频率为8000Hz，根据应变数据分析，在转速为65.18Hz时叶片发生同步共振，其振动频率约为2346.48Hz，如图1为采用应变法获取叶片同步共振的频率。

按照稀疏度的重构要求，应安装4个叶端定时传感器，同步转速参考传感器安装在轴的附近。随机布局4个叶端定时传感器，标记传感器顺序为S1-S4，测量叶端定时传感器之间的相对角度记录如表1所示。

表1传感器相对角度

调整旋转试验台转速为65.18Hz，采集叶端定时信号，每个叶端定时传感器采集11923圈实际到达时间，根据同步转速参考传感器采集的理想到达时间，计算每个叶端定时传感器的11923点的振动位移，如图2为根据4个叶端定时传感器信号计算得到的叶片的振动位移序列。如图3为计算的4个叶端定时传感器的振动位移序列的离散傅里叶变换结果，然后形成矩阵M_BTT。

设置满足奈奎斯特采样定理的参数Q为77，计算假设安装的每两个相邻的叶端定时传感器之间的角度应为360°/77≈4.68°。根据表1中的相对角度计算实际安装的叶端定时传感器的位置参数集合为C₄＝{0，2，16，62}。根据参数Q和集合C₄计算压缩矩阵Φ和补偿矩阵B，然后计算矩阵M。

根据公式(13)，对矩阵M进行逐个列向量的重构，获得重构的矩阵M′_NYQ。由于叶片发生同步共振，在叶端定时信号的离散傅里叶变换结果上，主要表现在k_P＝0，即矩阵M的第一列向量。在计算该列向量和压缩矩阵Φ的原子的内积值的模时，如图4为第一次迭代时不同原子的内积值的模，得到结果λ_(1，max)＝λ_(1，37)＝λ_(1，42)＝518.24。

为了对比本发明方法的优势，首先采用传统方法中每次迭代过程中只选择一个原子的方法。根据稀疏度K＝2，需要两次迭代过程，如图5为选择单个原子后的第二次迭代的内积值的模。在第一次迭代选择第37列原子后，更新的残差用于第二次迭代时，出现λ_(2，max)＝λ_(2，16)＝293.66＞λ_(2，42)。如图6为迭代过程选择单个原子的恢复频谱图，恢复的两个频率977.7Hz和2346.48Hz不满足离散傅里叶变换后的对称关系，两者之和不等于预设参数下的奈奎斯特采样频率5018.86Hz，属于错误的恢复结果。

而本发明会同时选择多个原子放入最优原子集合，所以同时选择第37列和第42列原子放入最优原子集合，利用最小二乘法求解高维向量中的稀疏估计值。如图7为本发明方法同时选择多个原子的恢复频谱图，恢复的两个频率2346.48Hz和2672.38Hz满足离散傅里叶变换后的对称关系，两者之和满足预设参数下的奈奎斯特采样频率5018.86Hz。恢复的频率2346.48Hz和应变信号中的频率2346.94Hz满足试验误差要求。这样，利用叶端定时欠采样的信号恢复出叶片的同步共振频率。本发明方法针对叶片同步共振下的叶端定时欠采样信号的特殊性质，准确得到叶片同步共振的频率。

以上所述实施例仅表达本发明的实施方式，但并不能因此而理解为对本发明专利的范围的限制，应当指出，对于本领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些均属于本发明的保护范围。

Claims (3)

1.一种基于叶端定时欠采样信号特性的叶片同步共振频率恢复方法，其特征在于，包括以下步骤：

第一步：在叶轮周围的圆周型壳体上随机安装P个叶端定时传感器，用来获取叶片的实际到达时间；将同步转速参考传感器固定在轴附近，用来获取叶片的理想到达时间；计算第p个叶端定时传感器获得的第n圈的叶端定时振动位移信号：

s_p[n]＝2πf_shaftR(t_p，n-t′_p，n) (1)

其中，p∈{1，2，...P}，f_shaft为叶轮的旋转频率，R为叶尖到叶轮中心的距离，t_p，n为通过叶端定时传感器获取的叶片实际到达时间，t′_p，n为通过同步转速参考传感器获取的叶片理想到达时间；

通过公式(1)计算每个叶端定时传感器每圈的振动位移后，对每个叶端定时传感器的振动位移序列进行离散傅里叶变换；

将每个叶端定时传感器计算得到的离散傅里叶变换结果作为行向量形成矩阵M_BTT：

第二步：根据实际安装的叶端定时传感器之间的相对角度，计算每个叶端定时传感器的位置在频谱恢复过程中的参数值；记录叶端定时传感器之间的相对角度集合α：

α＝{α_(1，1)，…，α_(1，p)，…α_(1，P)} (4)

其中，α_(1，p)表示第p个叶端定时传感器和第一个叶端定时传感器之间的相对角度；

在理想状态下，为了满足奈奎斯特采样定理，假设需要在叶轮周围圆周型壳体上均匀安装的叶端定时传感器个数为Q，定义每个叶端定时传感器的位置序号集合为C＝{0，1，...，Q-1}；每相邻的两个叶端定时传感器之间的夹角为：

α_Q＝360°/Q (5)

根据相对角度集合α中叶端定时传感器之间的相对角度，采用公式(6)计算实际安装的每个叶端定时传感器的相对应的位置参数：

c_p＝round(α_(1，p)/α_Q) (6)

其中，函数round为四舍五入取整函数；

因α_(1，1)＝0°，最终得到实际安装的叶端定时传感器的位置参数集合：

C_P＝{0，...c_p，...c_P} (7)

第三步：计算压缩矩阵Φ和矩阵M_BTT的补偿矩阵B；定义构成压缩矩阵Φ和补偿矩阵B的元素分别为：

其中，p∈{1，2，...P}，q∈{1，2，...Q}，k_P∈{0，1，...，N-1}，P表示实际安装的叶端定时传感器个数，Q表示为满足奈奎斯特采样定理假设安装的叶端定时传感器个数，N表示叶轮旋转的圈数，c_p表示实际安装的叶端定时传感器的位置参数，

表示压缩矩阵Φ中第p行第q列元素，

表示补偿矩阵B中第p行第(k_P+1)列元素；

由公式(8)、(9)可知，压缩矩阵Φ和补偿矩阵B的维度分别为P×Q和P×N；

第四步：利用最小二乘法重构叶片振动频率；假设安装Q个叶端定时传感器，叶轮旋转N圈后，获取的整个振动位移序列为s[m]，对于整个振动位移序列的离散傅里叶变换结果为：

其中，M＝Q·N，M表示假设安装Q个叶端定时传感器获取的整个振动位移序列的点数，k_Q∈{0，1，...，M-1}；

在已知的叶端定时采样序列的离散傅里叶变换结果S_p[k_P]和未知的奈奎斯特采样序列的离散傅里叶变换结果S[k_Q]之间存在如下关系：

其中，(q-1)N+k_p＝k_Q；

将M＝Q·N点未知的奈奎斯特采样序列的离散傅里叶变换结果S[k_Q]沿频率方向分割Q段，每段N点，每段作为行向量形成矩阵M_NYQ：

因此，公式(11)可以转化为矩阵的表达形式：

M＝ΦM_NYQ (13)

其中，

运算符号

为Hadamard积(矩阵对应元素相乘)，矩阵M_NYQ的维度为Q×N；

根据矩阵M中的低维列向量和矩阵M_NYQ中高维列向量的稀疏度K，采用最小二乘法经K次迭代获得M_NYQ中相对应的高维列向量的稀疏估计值，形成重构的矩阵M′_NYQ；

对于整个迭代过程：输入端为矩阵M中已知的低维列向量、压缩矩阵Φ和高维列向量的稀疏度K；输出端为矩阵M_NYQ中未知的高维列向量的稀疏估计值；

在每次迭代计算过程中，利用上一次迭代获得的残差向量与压缩矩阵Φ中的原子的内积值的模的最大值来寻找最佳原子，其中，第一次迭代计算时的残差向量为矩阵M中已知的低维列向量，压缩矩阵Φ中的原子为矩阵的列向量，需去除上一次迭代选择的原子；并根据最佳原子的位置确定高维向量中稀疏估计值的位置，表示为公式(14)：

其中，r_t-1为第(t-1)次迭代获得的残差向量，r₀为M中已知的低维列向量，

为压缩矩阵Φ的第q个原子；

叶片的同步共振频率为叶轮旋转频率的整数倍δ，在叶端定时信号上发生混叠效应，出现在S_p[k_P＝0]的位置上，得到β_p，0＝1；在进行第一次迭代时，r₀＝[QS₁[0]，...，QS_p[0]，...，QS_P[0]]^T，可以得到：

对于两个原子

和

有如下结果：

对比公式(16)和(17)，得到两个内积值的模相等的最佳原子：

λ_(1，δ+1)＝λ_(1，Q-δ+1) (18)

基于公式，每次迭代同时选择2个原子放入最优原子集合，再经迭代过程获取稀疏估计值；矩阵M的其他列向量的计算过程与第一个列向量相同，获得重构的矩阵M′_NYQ，最终实现叶片同步共振频率的恢复。

2.根据权利要求1所述的一种基于叶端定时欠采样信号特性的叶片同步共振频率恢复方法，其特征在于，第四步迭代过程中，叶片同步共振频率主要表现在矩阵M的第一个列向量，即r₀＝[QS₁[0]，...，QS_p[0]，...，QS_P[0]]^T，主要描述矩阵M的第一个列向量在第一次迭代计算过程中的相关问题，具体如下：

对于矩阵M的第一个列向量[QS₁[0]，...，QS_p[0]，...，QS_P[0]]^T的迭代计算过程：输入端为r₀＝[QS₁[0]，...，QS_p[0]，...，QS_P[0]]^T、压缩矩阵Φ和高维列向量的稀疏度K＝2；

第一次迭代时，有λ_(1，δ+1)＝λ_(1，Q-δ+1)，第(δ+1)和第(Q-δ+1)原子被同时选入最优原子集合：

求解r₀＝Φ₁θ₁的最小二乘解：

然后计算迭代后的残差向量r₁：

第二次迭代计算的输入向量为第一次迭代后得到的残差向量r1，计算过程与第一次迭代计算过程相同，其中第二次迭代中选择原子的过程需去除第一次迭代选择的第(δ+1)和第(Q-δ+1)原子；经过2次迭代可获得作为输出端的稀疏估计值

对应最优原子集合中原子的位置获得重构的高维列向量的稀疏估计值。

3.根据权利要求1或2所述的一种基于叶端定时欠采样信号特性的叶片同步共振频率恢复方法，其特征在于，第一步中，所述的采用公式(2)对每个叶端定时传感器的振动位移序列进行离散傅里叶变换：

其中，k_P∈{0，1，...，N-1}，N表示叶轮旋转的圈数，也就是每个叶端定时传感器的采样点数为N。

Images