CN104374939A - 基于振动信号同步压缩变换的旋转机械瞬时转速估测方法 - Google Patents

Abstract

本发明公开了一种基于振动信号同步压缩变换的旋转机械瞬时转速估测方法，包括以下步骤：1)获取旋转机械运行过程中的振动信号；2)对测得的振动信号进行频移处理；3)对频移后的振动信号进行同步压缩连续小波变换，得到频移后振动信号的时频分布；4)利用Viterbi算法从上述获得的时频分布中提取经频移处理后振动信号的一阶瞬时频率成分；5)利用提取的一阶瞬时频率恢复计算得到旋转机械的瞬时转速。本发明采用频移算法和同步压缩连续小波变换处理信号，实现振动信号瞬时频率的精确估计，Viterbi算法实现瞬时频率的精确提取，该方法可以实现对无法进行瞬时转速直接测量的旋转机械,通过振动信号准确地提取其瞬时转速。

Description

基于振动信号同步压缩变换的旋转机械瞬时转速估测方法

技术领域

本发明涉及一种旋转机械瞬时转速估测方法。

背景技术

旋转机械的瞬时转速中常常包含有大量自身运行状态信息，而这些信息对于设备的实时控制、状态检测和故障诊断等都具有十分重要的意义。如何能够准确的提取旋转机械的瞬时转速已经成为近年来的研究热点。

旋转机械的瞬时转速的直接测量方法是利用旋转编码器，对编码器产生的脉冲信号采用特定的信号处理方法进行处理，从而得到旋转机械的瞬时转速。但是由于一些诸如运行工况不允许、设计装配不合理以及成本的原因，使得对旋转设备瞬时转速的直接测量在很多时候无法实现。因此，许多基于旋转机械振动信号的瞬时转速提取方法被提出，诸如相位解调方法、基于信号模型的方法以及时频分布的方法等。这些方法当中，基于时频分布的方法由于其强大的时域信号分析能力，以及优异的抗噪声干扰特性而受到广泛的研究和应用。

由于旋转机械瞬时转速在工业应用中具有重要意义，国内许多学者对于瞬时转速的估测进行了较多研究。史振盛(史振盛，王桂华，等.基于瞬时频率估计的谐振分析及应用[J].内燃机与动力装置，2014,31(4)：11-15)提出了采用结合短时傅里叶变换和峰值搜索的基于时频分布的瞬时频率估计方法，并将其应用于谐振分析。曹书峰(曹书峰，等.基于时频融合的转速估计及轴承故障特征提取研究[J].振动与冲击，2013,32(8)：184-178)提出将Wigner-Vile分布与小波尺度谱进行融合，对融合后的时频分布用峰值搜索法估计转速，并将其运用于轴承故障特征提取中。胡爱军(胡爱军，朱瑜.基于改进峰值搜索法的旋转机械瞬时频率估计[J].振动与冲击，2013,32(7)：113-117)提出基于改进峰值搜索法的旋转机械瞬时频率估计方法，提出将瞬时频率中相邻两点一阶导数差值作为搜索峰值是否合理的判别条件，可避免传统峰值搜索法在干扰信号作用下提取虚假峰值，并将此应用与阶次跟踪。

申请号为201410104252.9的中国专利公开了一种基于拉格朗日插值多项式的瞬时转速估计方法，其特征在于,利用2阶拉格朗日插值多项式对电机位置轨迹进行函数拟合,然后把得到的位置函数对时间t求导得到其相应的速度函数,对指定的任一时刻,在速度函数中代入位置增量和时间间隔,即求得该时刻对应的瞬时速度。申请号为201310247876.1的中国专利公开了一种基于齿轮啮合振动的转速估计方法及装置，其特征在于计算齿轮啮合振动信号的Wigner-Ville分布和小波尺度谱；对所述Wigner-Ville分布和所述小波尺度谱进行时频融合,得到时频融合分布；根据所述时频融合分布,确定出齿轮啮合频率曲线；根据所述齿轮啮合频率曲线以及主动轮的齿数,确定出主动轮的转速曲线。

基于时频分布的瞬时转速估计主要包含两个过程：1)信号时频分布的获得；2)信号时频分布中瞬时频率的提取。而目前的研究多只考虑到这两个过程中的一个过程，而两者中任何一个的误差都将降低瞬时转速提取的精度，因此目前对于旋转机械瞬时转速估计的研究仍存在一定不足，提取的瞬时转速精度有待提高。

发明内容

本发明的目的在于利用同步压缩连续小波变换时频分布低频处的高频率分辨率，将高频成分平移到时频分布的低频位置以获得高的频率分辨率来处理信号，从而获得精确的高频成分瞬时频率估计值，然后利用比传统峰值检索方法更精确有效的Viterbi算法从时频分布中提取瞬时频率，进而获得精确的旋转机械的瞬时转速。

为达到以上目的，本发明是采取如下技术方案予以实现：

一种基于振动信号同步压缩变换的旋转机械瞬时转速估测方法，其特征在于，包括下述步骤：

(1)对从旋转机械运行过程中测得的振动信号s(t)构造其解析信号：

$\tilde{s}\left(t\right)=s\left(t\right)+\mathrm{jH}\left(s\left(t\right)\right)$

其中H(s(t))是s(t)的希尔伯特黄变换；

(2)根据分析信号所需要满足的频率分辨率Δf来选择频移频率ω₀，包含以下子步骤：

1)已知原始振动信号s(t)相邻两点之间的时间间隔为Δt，总数据长度为n,信号持续时间为T，构造变量L_f和H_f分别如下：

$\left\{\begin{array}{c}{L}_f={\log }_2\left(\frac{1}{T}\right)\\ H_f={\log }_2\left(\frac{1}{2\mathrm{\Δt}}\right)\end{array}\right.$

2)由L_f和H_f，得到同步压缩变换时频分布中成指数增长的离散频率序列：

$F_S\left(k\right)=2^{[L_f+\frac{H_f-L_f}{\mathrm{na}-1}(k-1)]},k=\mathrm{1,2},...,\mathrm{na}$

其中，na＝(log₂N-1)×nv，N是比n大的2的下一个指数次幂，nv为同步压缩变换中的一个参数；

3)对离散频率序列求差分，计算得到同步压缩变换时频分布中成指数增长的频率间隔序列，即相邻两个离散频率之间的差值：

$\mathrm{\Δf}\left(p\right)=2^{[L_f+\frac{H_f-L_f}{\mathrm{na}-1}p(p-1)]}(1-2^{-\frac{H_f-L_f}{\mathrm{na}-1}}),p=\mathrm{1,2},...,\mathrm{na}$

4)从频率间隔序列中找出所需频率分辨率Δf在频率间隔序列中所处的位置p₀，然后将p₀带入到子步骤2)所得的离散频率序列中，即可得到所需频率分辨率在时频分布中所对应的中心频率值ω₀；该ω₀即为频移频率。在获得频移频率ω₀之后，构造解析信号的频移信号：

${\tilde{s}}^{*}\left(t\right)=\tilde{s}\left(t\right)e^{{-\mathrm{j\ω}}_{0}t};$

(3)对步骤(2)获得的频移信号进行同步压缩变换，以得到频移信号的时频分布，所述同步压缩变换是基于连续小波变换来实现的，包含以下子步骤：

1)对原始振动信号s(t)进行连续小波变换：

$W_s(a,b)=\∫s\left(t\right)a^{-1/2}\stackrel{\‾}{\mathrm{\ψ}\left(\frac{t-b}{a}\right)\mathrm{dt}}$

其中，ψ是选取的母小波，a是小波变换的尺度因子，b是小波变换的时移因子；

2)对信号连续小波变换结果，利用帕萨瓦尔定理计算得到信号基于连续小波变换的瞬时频率：

${\mathrm{\ω}}_s(a,b)=-j{\left(W_s(a,b)\right)}^{-1}\frac{\∂}{\∂b}{W}_s(a,b)$

3)根据瞬时频率对信号连续小波变换时频分布进行映射变换，其离散表达式如下：

$T_s({\mathrm{\ω}}_l,b)={\left(\mathrm{\Δ\ω}\right)}^{-1}\underset{a_k:|\mathrm{\ω}(a_k,b)-{\mathrm{\ω}}_l|\≤\mathrm{\Δ\ω}/2}{\mathrm{\Σ}}{W}_s(a_k,b)a_k^{-3/2}{\left(\mathrm{\Δa}\right)}_k$

其中，Δω＝ω_l-ω_l-1，(Δa)_k＝a_k-a_k-1；经映射变换之后的连续小波变换时频分布即为同步压缩变换时频分布；

(4)对步骤(3)获得的时频分布，利用Viterbi算法提取时频分布中频移信号的第一阶瞬时频率f^*(t)；

(5)利用步骤(4)中提取得到的瞬时频率f^*(t)，根据f(t)＝f^*(t)+ω₀/2π计算原始振动信号s(t)的瞬时频率，从而得到旋转机械的瞬时转速。

上述方法中，步骤(4)所述从时频分布中提取瞬时频率的方法具体如下：

对于得到的同步压缩时频分布T_s(ω_l,b)，定义时间间隔n∈[n₁,n₂]内的所有路径都属于集合K，则信号s(t)从时刻n₁到n₂的瞬时频率估计是通过寻找最优路径，使得罚函数之和取得最小值来实现的，表达式如下：

$\widehat{\mathrm{\ω}}\left(n\right)=\arg \underset{k\left(n\right)\∈K}{\min }[\underset{n=n_1}{\overset{n_2-1}{\mathrm{\Σ}}}g(k\left(n\right),k(n+1))+\underset{n=n_1}{\overset{n_2}{\mathrm{\Σ}}}f\left(T_s(k\left(n\right),n)\right)]$

其中，g(x,y)和f(x)都是定义的不同罚函数；g(x,y)是一个线性形式的罚函数，表达式如下：

$g(x,y)=\left\{\begin{array}{cc}0,& |x-y|\≤\mathrm{\Δ}\\ c\left(\right|x-y|-\mathrm{\Δ}),& |x-y|>\mathrm{\Δ}\end{array}\right.$

其中，Δ为相邻两点之间的惩罚阈值，c为罚函数的权重；对于给定的时刻n，对获得的同步压缩时频分布中不同频率点对应的函数值进行降序排列，罚函数f(x)被定义为：

f(T_s(ω_l,n))＝q-1,q＝1,2,…,m

其中，q是T_s(ω_l,n)在降序序列中的所排的序号。

与现有技术相比，本发明具有下述优点：

1、由于本发明同时考虑并改善基于时频分布的瞬时转速估计方法中的两个过程。1)采用频移技术和同步压缩连续小波变换方法提高时频分布的时频凝聚性，从而提高信号瞬时频率的估计精度。2)采用Viterbi算法提取时频分布中的瞬时频率，改善传统峰值搜索方法在干扰噪声作用下提取虚假峰值的问题。从而可以实现更加精确的旋转机械瞬时转速的估计。

2、该方法对于噪声具有更好的抗干扰能力。

3、通过振动信号提取旋转机械的瞬时转速，不需要安装旋转编码器，简单易行，便于在工程实践中使用。

附图说明

图1为本发明方法流程示意图。

图2为转子振动位移信号时域波形。

图3为转子振动位移信号频谱布局放大图。

图4为本发明所得频移同步压缩时频分布。

图5为本发明提取的转子瞬时转速同转速计测得瞬时转速对比图。

具体实施方式

以下结合附图及一个具体实施例对本发明基于振动信号同步压缩变换的旋转机械瞬时转速估测的方法进行详细说明，同时验证本发明在工程应用中的有效性。该实施例以旋转机械振动位移信号为例进行说明，但本发明并不局限于使用位移信号，其他旋转机械振动信号如振动加速度信号等均可使用。

参考图1，本实施例在一个转子实验台上进行，利用一个位移传感器测量转子和定子之间的水平方向相对位移，同时利用转子试验台自带的转速计测量转子的实际瞬时转速，用测得的实际瞬时转速与利用本发明提取的瞬时转速进行对比以验证本发明在工程应用中的有效性。实验过程中，转子运行在4800rpm附近，通过手动调节使转子存在一定的瞬时转速波动，波动的幅值小于120rpm。

基于振动信号同步压缩变换的旋转机械瞬时转速估测方法，包括以下步骤：

(1)利用位移传感器测量电机转轴的振动位移信号s(t)；构造其解析信号：其中H(s(t))是s(t)的希尔伯特黄变换；

(2)根据分析信号所需要满足的频率分辨率Δf来选择频移频率ω₀，包含以下子步骤：

1)已知原始振动信号s(t)相邻两点之间的时间间隔为Δt，总数据长度为n,信号持续时间为T，构造变量L_f和H_f分别如下：

$\left\{\begin{array}{c}{L}_f={\log }_2\left(\frac{1}{T}\right)\\ H_f={\log }_2\left(\frac{1}{2\mathrm{\Δt}}\right)\end{array}\right.$

2)由L_f和H_f，得到同步压缩变换时频分布中成指数增长的离散频率序列：

$F_S\left(k\right)=2^{[L_f+\frac{H_f-L_f}{\mathrm{na}-1}(k-1)]},k=\mathrm{1,2},...,\mathrm{na}$

其中，na＝(log₂N-1)×nv，N是比n大的2的下一个指数次幂，nv为同步压缩变换中的一个参数；

3)对离散频率序列求差分，计算得到同步压缩变换时频分布中成指数增长的频率间隔序列，即相邻两个离散频率之间的差值：

$\mathrm{\Δf}\left(p\right)=2^{[L_f+\frac{H_f-L_f}{\mathrm{na}-1}p(p-1)]}(1-2^{-\frac{H_f-L_f}{\mathrm{na}-1}}),p=\mathrm{1,2},...,\mathrm{na}$

4)从频率间隔序列中找出所需频率分辨率Δf在频率间隔序列中所处的位置p₀，然后将p₀带入到子步骤2)所得的离散频率序列中，即可得到所需频率分辨率在时频分布中所对应的中心频率值ω₀；该ω₀即为频移频率。在获得频移频率ω₀之后，构造解析信号的频移信号：

${\tilde{s}}^{*}\left(t\right)=\tilde{s}\left(t\right)e^{{-\mathrm{j\ω}}_{0}t};$

本实施例中由于频率波动小于2Hz，选择0.2Hz的频率分辨率去分析信号，得到频移频率f₀＝55Hz，构造解析信号的频移信号

(3)对步骤(2)获得的频移信号进行同步压缩变换，以得到频移信号的时频分布，同步压缩变换是基于连续小波变换来实现的，包含以下子步骤：

1)对原始振动信号s(t)进行连续小波变换：

$W_s(a,b)=\∫s\left(t\right)a^{-1/2}\stackrel{\‾}{\mathrm{\ψ}\left(\frac{t-b}{a}\right)\mathrm{dt}}$

其中，ψ是选取的母小波，a是小波变换的尺度因子，b是小波变换的时移因子；

2)对信号连续小波变换结果，利用帕萨瓦尔定理计算得到信号基于连续小波变换的瞬时频率：

${\mathrm{\ω}}_s(a,b)=-j{\left(W_s(a,b)\right)}^{-1}\frac{\∂}{\∂b}{W}_s(a,b)$

3)根据瞬时频率对信号连续小波变换时频分布进行映射变换，其离散表达式如下：

$T_s({\mathrm{\ω}}_l,b)={\left(\mathrm{\Δ\ω}\right)}^{-1}\underset{a_k:|\mathrm{\ω}(a_k,b)-{\mathrm{\ω}}_l|\≤\mathrm{\Δ\ω}/2}{\mathrm{\Σ}}{W}_s(a_k,b)a_k^{-3/2}{\left(\mathrm{\Δa}\right)}_k$

其中，Δω＝ω_l-ω_l-1，(Δa)_k＝a_k-a_k-1；经映射变换之后的连续小波变换时频分布即为同步压缩变换时频分布；同步压缩时频分布较连续小波时频分布具有更好的时频凝聚性。本实施例得到的时频分布如图4所示。

(4)对步骤(3)获得的时频分布，利用Viterbi算法提取时频分布中频移信号的第一阶瞬时频率f^*(t)；具体如下：

对于得到的同步压缩时频分布T_s(ω_l,b)，定义时间间隔n∈[n₁,n₂]内的所有路径都属于集合K，则信号s(t)从时刻n₁到n₂的瞬时频率估计是通过寻找最优路径，使得罚函数之和取得最小值来实现的，表达式如下：

$\widehat{\mathrm{\ω}}\left(n\right)=\arg \underset{k\left(n\right)\∈K}{\min }[\underset{n=n_1}{\overset{n_2-1}{\mathrm{\Σ}}}g(k\left(n\right),k(n+1))+\underset{n=n_1}{\overset{n_2}{\mathrm{\Σ}}}f\left(T_s(k\left(n\right),n)\right)]$

其中，g(x,y)和f(x)都是定义的不同罚函数；g(x,y)是一个线性形式的罚函数，表达式如下：

$g(x,y)=\left\{\begin{array}{cc}0,& |x-y|\≤\mathrm{\Δ}\\ c\left(\right|x-y|-\mathrm{\Δ}),& |x-y|>\mathrm{\Δ}\end{array}\right.$

其中，Δ为相邻两点之间的惩罚阈值，c为罚函数的权重；对于给定的时刻n，对获得的同步压缩时频分布中不同频率点对应的函数值进行降序排列，罚函数f(x)被定义为：

f(T_s(ω_l,n))＝q-1,q＝1,2,…,m

其中，q是T_s(ω_l,n)在降序序列中的所排的序号。

(5)利用步骤(4)中提取得到的瞬时频率f^*(t)，根据f(t)＝f^*(t)+ω₀/2π计算原始振动信号s(t)的瞬时频率，从而得到旋转机械的瞬时转速。

本实施例测得的转子振动位移信号的时域波形如图2所示，位移信号频谱的局部放大图如图3所示，从图3局部放大的频谱中我们可以看到，位移信号存在一定的非平稳特性，基频80Hz附近频谱聚集性较差，在频谱中占据一定的宽度。频移同步压缩的时频分布如图4所示，从图4中，我们发现频移同步压缩具有很高的时频凝聚性，能够清晰的显示出信号的瞬时频率波动。本实施例提取的瞬时转速与转速计测量的瞬时转速对比如图5所示。图5表明本发明提取的瞬时转速同转速计测量得到的瞬时转速基本一致，两者之间的相对误差小于0.5％，说明本发明可以有效地从振动信号中提取精确的瞬时转速，验证了本发明所述方法的有效性。

Claims (2)

1.一种基于振动信号同步压缩变换的旋转机械瞬时转速估测方法，其特征在于，包括下述步骤：

(1)对从旋转机械运行过程中测得的振动信号s(t)构造其解析信号：

$\tilde{s}\left(t\right)=s\left(t\right)+\mathrm{jH}\left(s\left(t\right)\right)$

其中H(s(t))是s(t)的希尔伯特黄变换；

(2)根据分析信号所需要满足的频率分辨率Δf来选择频移频率ω₀，包含以下子步骤：

1)已知原始振动信号s(t)相邻两点之间的时间间隔为Δt，总数据长度为n,信号持续时间为T，构造变量L_f和H_f分别如下：

$\left\{\begin{array}{c}{L}_f={\log }_2\left(\frac{1}{T}\right)\\ H_f={\log }_2\left(\frac{1}{2\mathrm{\Δt}}\right)\end{array}\right.$

2)由L_f和H_f，得到同步压缩变换时频分布中成指数增长的离散频率序列：

$F_s\left(k\right)=2^{[L_f+\frac{H_f-L_f}{\mathrm{na}-1}(k-1)]},k=\mathrm{1,2},...,\mathrm{na}$

其中，na＝(log₂N-1)×nv，N是比n大的2的下一个指数次幂，nv为同步压缩变换中的一个参数；

3)对离散频率序列求差分，计算得到同步压缩变换时频分布中成指数增长的频率间隔序列，即相邻两个离散频率之间的差值：

$\mathrm{\Δf}\left(p\right)=2^{[L_f+\frac{H_f-L_f}{\mathrm{na}-1}(p-1)]}(1-2^{-\frac{H_f-L_f}{\mathrm{na}-1}}),p=\mathrm{1,2},...,\mathrm{na}$

4)从频率间隔序列中找出所需频率分辨率Δf在频率间隔序列中所处的位置p₀，然后将p₀带入到子步骤2)所得的离散频率序列中，即可得到所需频率分辨率在时频分布中所对应的中心频率值ω₀；该ω₀即为频移频率；在获得频移频率ω₀之后，构造解析信号的频移信号：

${\tilde{s}}^{*}\left(t\right)=\tilde{s}\left(t\right)e^{-j{\mathrm{\ω}}_{0}t};$

(3)对步骤(2)获得的频移信号进行同步压缩变换，以得到频移信号的时频分布，所述同步压缩变换是基于连续小波变换来实现的，包含以下子步骤：

1)对原始振动信号s(t)进行连续小波变换：

$W_s(a,b)=\∫s\left(t\right)a^{-1/2}\stackrel{\‾}{\mathrm{\ψ}\left(\frac{t-b}{a}\right)}\mathrm{dt}$

其中，ψ是选取的母小波，a是小波变换的尺度因子，b是小波变换的时移因子；

2)对信号连续小波变换结果，利用帕萨瓦尔定理计算得到信号基于连续小波变换的瞬时频率：

${\mathrm{\ω}}_s(a,b)=-j{\left(W_s(a,b)\right)}^{-1}\frac{\∂}{\∂b}{W}_s(a,b)$

3)根据瞬时频率对信号连续小波变换时频分布进行映射变换，其离散表达式如下：

$T_s({\mathrm{\ω}}_l,b)={\left(\mathrm{\Δ\ω}\right)}^{-1}\underset{a_k:|\mathrm{\ω}(a_k,b)-{\mathrm{\ω}}_l|\≤\mathrm{\Δ\ω}/2}{\mathrm{\Σ}}{W}_s(a_k,b)a_k^{-3/2}{\left(\mathrm{\Δa}\right)}_k$

其中，Δω＝ω_l-ω_l-1，(Δa)_k＝a_k-a_k-1；经映射变换之后的连续小波变换时频分布即为同步压缩变换时频分布；

(4)对步骤(3)获得的时频分布，利用Viterbi算法提取时频分布中频移信号的第一阶瞬时频率f^*(t)；

(5)利用步骤(4)中提取得到的瞬时频率f^*(t)，根据f(t)＝f^*(t)+ω₀/2π计算原始振动信号s(t)的瞬时频率，从而得到旋转机械的瞬时转速。

2.如权利要求1所述的基于振动信号同步压缩变换的旋转机械瞬时转速估测方法，其特征在于，步骤(4)所述从时频分布中提取瞬时频率的方法具体如下：

对于得到的同步压缩时频分布T_s(ω_l,b)，定义时间间隔n∈[n₁,n₂]内的所有路径都属于集合K，则信号s(t)从时刻n₁到n₂的瞬时频率估计是通过寻找最优路径，使得罚函数之和取得最小值来实现的，表达式如下：

$\widehat{\mathrm{\ω}}\left(n\right)=\arg \underset{k\left(n\right)\∈K}{\min }[\underset{n=n_1}{\overset{n_2-1}{\mathrm{\Σ}}}g(k\left(n\right),k(n+1))+\underset{n=n_1}{\overset{n_2}{\mathrm{\Σ}}}f\left(T_s(k\left(n\right),n)\right)]$

其中，g(x,y)和f(x)都是定义的不同罚函数；g(x,y)是一个线性形式的罚函数，表达式如下：

$g(x,y)=\left\{\begin{array}{cc}0,& |x-y|\≤\mathrm{\Δ}\\ c\left(\right|x-y|-\mathrm{\Δ}),& |x-y|>\mathrm{\Δ}\end{array}\right.$

其中，Δ为相邻两点之间的惩罚阈值，c为罚函数的权重；对于给定的时刻n，对获得的同步压缩时频分布中不同频率点对应的函数值进行降序排列，罚函数f(x)被定义为：

f(T_s(ω_l,n))＝q-1,q＝1,2,…,m

其中，q是T_s(ω_l,n)在降序序列中的所排的序号。