CN107783938A - 一种旋转设备瞬时转速估计方法 - Google Patents

一种旋转设备瞬时转速估计方法 Download PDF

Info

Publication number
CN107783938A
CN107783938A CN201710781137.9A CN201710781137A CN107783938A CN 107783938 A CN107783938 A CN 107783938A CN 201710781137 A CN201710781137 A CN 201710781137A CN 107783938 A CN107783938 A CN 107783938A
Authority
CN
China
Prior art keywords
mrow
msub
frequency
signal
envelope
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201710781137.9A
Other languages
English (en)
Other versions
CN107783938B (zh
Inventor
彭志科
陈是扦
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shanghai Jiaotong University
Original Assignee
Shanghai Jiaotong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shanghai Jiaotong University filed Critical Shanghai Jiaotong University
Priority to CN201710781137.9A priority Critical patent/CN107783938B/zh
Publication of CN107783938A publication Critical patent/CN107783938A/zh
Application granted granted Critical
Publication of CN107783938B publication Critical patent/CN107783938B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02P—CONTROL OR REGULATION OF ELECTRIC MOTORS, ELECTRIC GENERATORS OR DYNAMO-ELECTRIC CONVERTERS; CONTROLLING TRANSFORMERS, REACTORS OR CHOKE COILS
    • H02P23/00—Arrangements or methods for the control of AC motors characterised by a control method other than vector control
    • H02P23/0077—Characterised by the use of a particular software algorithm
    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02P—CONTROL OR REGULATION OF ELECTRIC MOTORS, ELECTRIC GENERATORS OR DYNAMO-ELECTRIC CONVERTERS; CONTROLLING TRANSFORMERS, REACTORS OR CHOKE COILS
    • H02P23/00—Arrangements or methods for the control of AC motors characterised by a control method other than vector control
    • H02P23/14—Estimation or adaptation of motor parameters, e.g. rotor time constant, flux, speed, current or voltage

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Power Engineering (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

本发明提供了一种旋转设备瞬时转速估计方法，该方法通过精确估计转子振动信号瞬时频率，实现设备转速估计。本发明首先利用参数化时频分析方法粗略估计振动信号瞬时频率，然后用一种新的包络跟踪滤波器对初始估计结果进行精确修正，最终得到的信号瞬时频率直接反映设备的瞬时转速，该包络跟踪滤波器利用傅立叶级数模型刻画信号包络，采用最小二乘法对包络参数优化求解。本发明采用参数化时频变换与包络跟踪滤波器结合的瞬时转速估计方法，能够克服传统技术估计精度受时频分辨率限制的缺点，具有估计精度高，实现简单的优点。

Description

一种旋转设备瞬时转速估计方法
技术领域
本发明涉及信号处理领域，尤其是一种估计旋转设备瞬时转速的方法，具体是一种基于参数化时频变换与包络跟踪滤波的瞬时转速估计方法。
背景技术
在大型旋转设备状态监测、故障诊断等应用中，通常需要精确测量设备的瞬时转速。键相测量是目前常用的一种设备瞬时转速测量方法。键相测量方法通过探测键相脉冲信号产生的时刻，确定转子在一个旋转周期内的位置。这种方法得到的是转子的平均转速，当转速变化时，测量误差较大。同时，键相测量需要额外的硬件设备，测量成本较高。转子振动信号的频率信息直接与转子瞬时转速相关。因此，通过分析振动信号，可以精确估计旋转设备瞬时转速。
当旋转设备转速变化时，转子振动信号的频率成分呈现出时变特征。瞬时频率常用来刻画信号的这类瞬时变化特征。时频分析技术是估计信号瞬时频率的有效工具，根据分析参数是否与信号有关，时频方法可以分为非参数化时频方法和参数化时频方法。常见的短时傅里叶变换、小波变换、魏格纳威尔分布等方法属于非参数化时频方法。这些方法的缺点是：时频集中性差，交叉项干扰等。为了提高非参数化时频方法的分辨率，有学者提出了时频重排和同步压缩等方法。这些方法通过后处理手段，将时频分布上每一点的值搬移到信号能量重心，以压缩时频能量带，提高时频分辨率，但方法计算量大、抗噪性能差。另一类非参数化时频方法利用先进的数学优化方法提升信号瞬时频率估计精度，例如现有技术中提出的数据驱动时频分析方法。这类方法对优化初值敏感，当初值选择不当时，算法难以收敛。参数化时频方法采用与信号模型匹配的核函数，能够有效提高时频分辨率，减小瞬时频率估计偏差。线性调频小波变换是处理线性调频信号的一种参数化时频方法。已有现有技术扩展了线性调频小波变换，提出了多种处理非线性调频信号的时频分析方法。参数化时频分析方法虽然提升了时频分布的能量聚集性，但是其时频分辨率仍然有限，难以得到高精度的瞬时频率估计值。
至今为止，还没有采用参数化时频方法和包络跟踪滤波相结合进行旋转设备瞬时转速估计。
发明内容
为了克服传统技术估计精度受时频分辨率限制的缺点，本发明提供一种基于参数化时频变换与包络跟踪滤波结合的方法实现对旋转设备瞬时转速的高精度估计：首先采用以傅立叶级数为核函数的参数化时频变换粗略估计振动信号瞬时频率；以瞬时频率初始估计值为输入，利用一种新的包络跟踪滤波器估计信号复数包络；利用信号包络的相位信息精确修正瞬时频率估计结果，得到的振动信号瞬时频率可以直接反映设备的瞬时转速。本发明提供的包络跟踪滤波算法利用傅立叶级数模型描述信号复数包络，用正则化最小二乘法估计包络模型参数。该包络跟踪算法可以有效抑制噪声，精确提取各种复杂的信号包络信息。
本发明是根据以下技术方案实现的：
一种旋转设备瞬时转速估计方法，其特征在于，包括如下步骤：
步骤S1，利用参数化时频变换估计设备振动信号的瞬时频率；
步骤S2，利用包络跟踪滤波器估计信号复数包络；
步骤S3，得到的信号包络被用于精确修正瞬时频率初始估计结果，信号的瞬时频率直接反映设备瞬时转速。
上述技术方案中，步骤S1通过迭代拟合参数化时频分布脊线的方法估计信号瞬时频率，其过程包括：
S101：参数初始化，令核参数α(i)＝{0，…，0}，设定傅立叶级数阶次M，收敛阈值ε，迭代次数i＝1；
S102：以核参数α(i)计算参数化时频分布TF(t，f；α(i))；
S103：从时频分布TF(t，f；α(i))中提取时频脊线
S104：用M阶傅立叶级数似合更新核参数α(i+1)；
S105：计算迭代终止条件
S106：若ξ(i)＞ε，令i＝i+1并跳转至步骤S102；否则执行步骤S107；
S107：返回瞬时频率估计结果
上述技术方案中，步骤S1中的参数化时频变换具体包括：
定义：
其中
z(t)是由希尔伯特变换得到的振动信号的解析形式，gσ(t)为高斯窗函数，α＝为变换核参数，F0＝Fs/2N为傅立叶级数的频率分辨率(Fs为信号采样频率)，参数化时频变换利用算子降低目标信号的频率调制程度，以得到集中的时频表示。
上述技术方案中，步骤S2中的包络跟踪滤波器算法具体包括：
假设当旋转设备转速发生变化时，转子振动信号模型用调幅-调频模型描述：
其中t＝t0，…tN-1为采样时刻，A(t)、f(t)、分别表示振动信号基频分量的瞬时幅值、瞬时频率和初始相位，n(t)代表噪声和其他不相关的谐波成分。瞬时频率f(t)直接对应旋转设备的瞬时转速。
假设参数化时频变换得到的瞬时频率估计值为振动信号模型重新整理为：
其中为信号复数包络，a(t)的复数相位与瞬时频率估计误差有关，因此可以利用a(t)的相位信息精确修正瞬时频率估计结果。为了精确提取a(t)，用K阶傅立叶级数模型刻画a(t)：
将上式带入信号模型，得到如下回归方程：
z＝Ga+n
其中z＝[z(t0)…z(tN-1)]T，a＝[a0…aKb1…bK]T，n＝[n(t0)…n(tN-1)]T，G矩阵的维数为N×(2K+1)，其元素为：
其中
信号包络系数向量a通过正则化最小二乘法估计：
其中α为正则化参数，I代表单位矩阵，上标H代表矩阵共轭转置，由估计到的系数向量可以得到信号包络估计进而提取信号基频分量包络跟踪滤波器可视为时频滤波器，其中心频率为参数化时频变换得到的瞬时频率其带宽为2KF0。
上述技术方案中，步骤S3中瞬时频率修正过程具体包括：
S301：将步骤S1中估计到的瞬时频率作为步骤S2中包络跟踪滤波器输入，提取信号复数包络具体需建立信号包络傅立叶级数模型，并对模型参数优化求解；
S302：提取包络信号的复数相位对相位函数求导得到包络瞬时频率用傅立叶级数模型对拟合得到瞬时频率修正量
S303：修正瞬时频率估计
与现有技术相比，本发明具有如下的有益效果：
1.由于采用傅立叶级数瞬时频率模型，与传统的多项式模型相比，本发明能够处理更复杂的设备瞬时转速变化情况。
2.本发明采用包络跟踪滤波器修正瞬时频率估计结果，能够克服传统技术估计精度受时频分辨率限制的缺点，能够得到高精度的瞬时转速估计。
3.算法实现简单，性能稳定，能够适用于多种应用领域。
附图说明
通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：
图1为本发明的瞬时转速估计方法流程示意图；
图2为本发明的仿真信号的时频表示示意图；
图3为本发明的仿真信号瞬时频率估计结果示意图；
图4为本发明的不同信噪比下仿真信号瞬时频率估计误差示意图；
图5为本发明的转子试验台示意图；
图6为本发明的转子振动信号的时频表示示意图；
图7为本发明的转子瞬时转速估计结果示意图。
具体实施方式
下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变化和改进。这些都属于本发明的保护范围。
图1为本发明的瞬时转速估计方法整体流程示意图，如图1所示，其特征在于，包括如下步骤：
步骤S1，利用参数化时频变换估计设备振动信号的瞬时频率；
步骤S2，利用包络跟踪滤波器估计信号复数包络；
步骤S3，得到的信号包络被用于精确修正瞬时频率初始估计结果，信号的瞬时频率直接反映设备瞬时转速。
步骤S1在估计过程中迭代拟合参数化时频分布脊线，其过程包括：
S101：参数初始化，令核参数α(i)＝{0，…，0}，设定傅立叶级数阶次M，收敛阈值ε，迭代次数i＝1；
S102：以核参数α(i)计算参数化时频分布TF(t，f；α(i))；
S103：从时频分布TF(t，f；α(i))中提取时频脊线
S104：用M阶傅立叶级数拟合更新核参数α(i+1)；
S105：计算迭代终止条件
S106：若ξ(i)＞ε，令i＝i+1并跳转至步骤S102；否则执行步骤S107；
S107：返回瞬时频率估计结果
上述技术方案中，步骤S1中的参数化时频变换具体包括：
定义：
其中
z(t）是由希尔伯特变换得到的振动信号的解析形式，gσ（t)为高斯窗函数， 为变换核参数，F0＝Fs／2N为傅立叶级数的频率分辨率(Fs为信号采样频率)，参数化时频变换利用算子降低目标信号的频率调制程度，以得到集中的时频表示。
步骤S2中的包络跟踪滤波算法估计信号包络的过程具体包括：
假设当旋转设备转速发生变化时，转子振动信号模型可用调幅-调频模型描述：
其中t＝t0，…tN-1为采样时刻，A（t）、f（t)、分别表示振动信号基频分量的瞬时幅值、瞬时频率和初始相位，n(t)代表噪声和其他不相关的谐波成分。瞬时频率f(t)直接对应旋转设备的瞬时转速。
假设参数化时频变换得到的瞬时频率估计值为振动信号模型重新整理为：
其中为信号复数包络，a(t)的复数相位与瞬时频率估计误差有关，因此可以利用a(t)的相位信息精确修正瞬时频率估计结果。为了精确提取a(t)，用K阶傅立叶级数模型刻画a(t)：
将上式带入信号模型，得到如下回归方程：
z＝Ga+n
其中z＝[z(t0)…z(tN-1)]T，a＝[a0…aKb1…bK]T，n＝[n(t0)…n(tN-1)]T，G矩阵的维数为N×(2K+1)，其元素为：
其中
信号包络系数向量a通过正则化最小二乘法估计：
其中α为正则化参数，I代表单位矩阵，上标H代表矩阵共轭转置，由估计到的系数向量可以得到信号包络估计进而提取信号基频分量包络跟踪滤波器可视为时频滤波器，其中心频率为参数化时频变换得到的瞬时频率其带宽为2KF0。
步骤S3中瞬时频率修正过程具体包括：
S301：将步骤S1中估计到的瞬时频率作为步骤S2中包络跟踪滤波器输入，提取信号复数包络
S302：提取包络信号的复数相位对相位函数求导得到包络瞬时频率用傅立叶级数模型对拟合得到瞬时频率修正量
S303：修正瞬时频率估计
实施例1
图2为一仿真信号的时频表示，信噪比为0dB，采样频率为100Hz。设定参数化时频变换傅立叶核的阶次M为6，迭代收敛阈值ε为1e-3。参数化时频变换总共迭代3次达到收敛条件，每次迭代中瞬时频率估计相对误差分别为-43.9dB、-44.9dB、-45.0dB。利用包络跟踪滤波器对瞬时频率估计结果进行修正，设定滤波器带宽为1Hz(即傅立叶阶次K＝15)，正则化参数α为0.5，修正后的瞬时频率如图3所示，其相对误差为-64.9dB。图4为参数化时频变换和本发明提供的方法在不同信噪比条件下瞬时频率估计的相对误差。结果表明本发明在强噪声情况下也能得到精确的瞬时频率估计。
实施例2
本发明提供的方法将用于估计实际设备的瞬时转速。图5为一转子试验台示意图，其中，电机经过减速器调速后驱动负载转子运动，加速度计采集轴承处的振动信号用于分析处理。图6为试验台启停机阶段实测的转子振动信号的时频表示，采样频率为100Hz。设定参数化时频变换傅立叶核的阶次M为30，迭代收敛阈值ε为1e-3。参数化时频变换总共迭代3次。利用包络跟踪滤波器对瞬时频率估计结果进行修正，设定滤波器带宽为1Hz(即傅立叶阶次K＝15)，正则化参数α为0.5。最终估计到的转子瞬时转速如图7所示。
以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变化或修改，这并不影响本发明的实质内容。在不冲突的情况下，本申请的实施例和实施例中的特征可以任意相互组合。

Claims (5)

1.一种旋转设备瞬时转速估计方法，其特征在于，包括如下步骤：
步骤S1，利用参数化时频变换估计设备振动信号的瞬时频率；
步骤S2，利用包络跟踪滤波器估计信号复数包络；
步骤S3，得到的信号包络被用于精确修正瞬时频率初始估计结果，信号的瞬时频率直接反映设备瞬时转速。
2.根据权利要求1所述的一种旋转设备瞬时转速估计方法，其特征在于，步骤S1通过迭代拟合参数化时频分布脊线的方法估计信号瞬时频率，其过程包括：
S101：参数初始化，令核参数α(i)＝{0，…，0}，设定傅立叶级数阶次M，收敛阈值ε，迭代次数i＝1；
S102：以核参数α(i)计算参数化时频分布TF(t，f；α(i))；
S103：从时频分布TF(t，f；α(i))中提取时频脊线
S104：用M阶傅立叶级数拟合更新核参数α(i+1)；
S105：计算迭代终止条件
S106：若ξ(i)＞ε，令i＝i+1并跳转至步骤S102；否则执行步骤S107；
S107：返回瞬时频率估计结果
3.根据权利要求1所述的一种旋转设备瞬时转速估计方法，其特征在于，步骤S1中的参数化时频变换具体包括：
定义：
<mrow> <mi>T</mi> <mi>F</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>,</mo> <mi>f</mi> <mo>;</mo> <mi>&amp;alpha;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>&amp;Integral;</mo> <mover> <mi>z</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <msub> <mi>g</mi> <mi>&amp;sigma;</mi> </msub> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>-</mo> <mi>t</mi> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>j</mi> <mn>2</mn> <mi>&amp;pi;</mi> <mi>f</mi> <mi>&amp;tau;</mi> </mrow> </msup> <mi>d</mi> <mi>&amp;tau;</mi> </mrow>
其中
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mover> <mi>z</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>z</mi> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <msubsup> <mi>&amp;Psi;</mi> <mi>&amp;alpha;</mi> <mi>D</mi> </msubsup> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <msubsup> <mi>&amp;Psi;</mi> <mrow> <mi>t</mi> <mo>,</mo> <mi>&amp;alpha;</mi> </mrow> <mi>C</mi> </msubsup> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msubsup> <mi>&amp;Psi;</mi> <mi>&amp;alpha;</mi> <mi>D</mi> </msubsup> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>j</mi> <mn>2</mn> <msubsup> <mi>&amp;pi;&amp;Sigma;</mi> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </msubsup> <mrow> <mo>(</mo> <mfrac> <msub> <mi>&amp;alpha;</mi> <mi>k</mi> </msub> <mrow> <mn>2</mn> <msub> <mi>&amp;pi;kF</mi> <mn>0</mn> </msub> </mrow> </mfrac> <mi>s</mi> <mi>i</mi> <mi>n</mi> <mo>(</mo> <mrow> <mn>2</mn> <msub> <mi>&amp;pi;kF</mi> <mn>0</mn> </msub> <mi>&amp;tau;</mi> </mrow> <mo>)</mo> <mo>-</mo> <mfrac> <msub> <mover> <mi>&amp;alpha;</mi> <mo>&amp;OverBar;</mo> </mover> <mi>k</mi> </msub> <mrow> <mn>2</mn> <msub> <mi>&amp;pi;kF</mi> <mn>0</mn> </msub> </mrow> </mfrac> <mi>c</mi> <mi>o</mi> <mi>s</mi> <mo>(</mo> <mrow> <mn>2</mn> <msub> <mi>&amp;pi;kF</mi> <mn>0</mn> </msub> <mi>&amp;tau;</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> </mrow> </msup> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msubsup> <mi>&amp;Psi;</mi> <mrow> <mi>t</mi> <mo>,</mo> <mi>&amp;alpha;</mi> </mrow> <mi>C</mi> </msubsup> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>e</mi> <mrow> <mi>j</mi> <mn>2</mn> <msubsup> <mi>&amp;pi;&amp;tau;&amp;Sigma;</mi> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>&amp;alpha;</mi> <mi>k</mi> </msub> <mi>cos</mi> <mo>(</mo> <mrow> <mn>2</mn> <msub> <mi>&amp;pi;kF</mi> <mn>0</mn> </msub> <mi>t</mi> </mrow> <mo>)</mo> <mo>+</mo> <msub> <mover> <mi>&amp;alpha;</mi> <mo>&amp;OverBar;</mo> </mover> <mi>k</mi> </msub> <mi>sin</mi> <mo>(</mo> <mrow> <mn>2</mn> <msub> <mi>&amp;pi;kF</mi> <mn>0</mn> </msub> <mi>&amp;tau;</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> </mrow> </msup> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> </mrow>
z(t)是由希尔伯特变换得到的振动信号的解析形式，gσ(t)为高斯窗函数， 为变换核参数，F0＝Fs/2N为傅立叶级数的频率分辨率(Fs为信号采样频率)，参数化时频变换利用算子降低目标信号的频率调制程度，以得到集中的时频表示。
4.根据权利要求1所述的一种旋转设备瞬时转速估计方法，其特征在于，步骤S2中的包络跟踪滤波器具体包括：
当旋转设备转速发生变化时，转子振动信号模型用调幅-调频模型描述：
其中t＝t0，…tN-1为采样时刻，A(t)、f(t)、分别表示振动信号基频分量的瞬时幅值、瞬时频率和初始相位，n(t)代表噪声和其他不相关的谐波成分。瞬时频率f(t)直接对应旋转设备的瞬时转速。
假设参数化时频变换得到的瞬时频率估计值为振动信号模型重新整理为：
<mrow> <mi>z</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>a</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mi>j</mi> <mn>2</mn> <mi>&amp;pi;</mi> <msubsup> <mo>&amp;Integral;</mo> <mn>0</mn> <mi>t</mi> </msubsup> <mover> <mi>f</mi> <mo>~</mo> </mover> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <mi>d</mi> <mi>r</mi> </mrow> </msup> <mo>+</mo> <mi>n</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>,</mo> <mi>t</mi> <mo>=</mo> <msub> <mi>t</mi> <mn>0</mn> </msub> <mo>,</mo> <mo>...</mo> <msub> <mi>t</mi> <mrow> <mi>N</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> </mrow>
其中为信号复数包络，a(t)的复数相位与瞬时频率估计误差有关，因此可以利用a(t)的相位信息精确修正瞬时频率估计结果。为了精确提取a(t)，用K阶傅立叶级数模型刻画a(t)：
<mrow> <mi>a</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>a</mi> <mn>0</mn> </msub> <mo>+</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>K</mi> </munderover> <msub> <mi>a</mi> <mi>k</mi> </msub> <mi>c</mi> <mi>o</mi> <mi>s</mi> <mrow> <mo>(</mo> <mn>2</mn> <msub> <mi>&amp;pi;kF</mi> <mn>0</mn> </msub> <mi>t</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>b</mi> <mi>k</mi> </msub> <mi>s</mi> <mi>i</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>2</mn> <msub> <mi>&amp;pi;kF</mi> <mn>0</mn> </msub> <mi>t</mi> <mo>)</mo> </mrow> </mrow>
将上式带入信号模型，得到如下回归方程：
z＝Ga+n
其中z＝[z(t0)…z(tN-1)]T，a＝[a0…aKb1…bK]T，n＝[n(t0)…n(tN-1)]T，G矩阵的维数为N×(2K+1)，其元素为：
<mrow> <msub> <mrow> <mo>(</mo> <mi>G</mi> <mo>)</mo> </mrow> <mrow> <mi>c</mi> <mi>d</mi> </mrow> </msub> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <msup> <mi>e</mi> <mrow> <mi>j</mi> <mi>&amp;psi;</mi> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mrow> <mi>c</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </msup> <mo>&amp;CenterDot;</mo> <mi>c</mi> <mi>o</mi> <mi>s</mi> <mo>&amp;lsqb;</mo> <mn>2</mn> <mo>(</mo> <mi>d</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> <mi>&amp;pi;</mi> <msub> <mi>F</mi> <mn>0</mn> </msub> <msub> <mi>t</mi> <mrow> <mi>c</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>&amp;rsqb;</mo> <mo>,</mo> <mi>d</mi> <mo>&amp;Element;</mo> <mo>&amp;lsqb;</mo> <mn>1</mn> <mo>,</mo> <mi>K</mi> <mo>+</mo> <mn>1</mn> <mo>&amp;rsqb;</mo> </mtd> </mtr> <mtr> <mtd> <msup> <mi>e</mi> <mrow> <mi>j</mi> <mi>&amp;psi;</mi> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mrow> <mi>c</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </msup> <mo>&amp;CenterDot;</mo> <mi>sin</mi> <mo>&amp;lsqb;</mo> <mn>2</mn> <mo>(</mo> <mi>d</mi> <mo>-</mo> <mi>K</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> <mi>&amp;pi;</mi> <msub> <mi>F</mi> <mn>0</mn> </msub> <msub> <mi>t</mi> <mrow> <mi>c</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>&amp;rsqb;</mo> <mo>,</mo> <mi>d</mi> <mo>&amp;Element;</mo> <mo>&amp;lsqb;</mo> <mi>K</mi> <mo>+</mo> <mn>2</mn> <mo>,</mo> <mn>2</mn> <mi>K</mi> <mo>+</mo> <mn>1</mn> <mo>&amp;rsqb;</mo> </mtd> </mtr> </mtable> </mfenced> </mrow>
其中
信号包络系数向量a通过正则化最小二乘法估计：
<mrow> <mover> <mi>a</mi> <mo>~</mo> </mover> <mo>=</mo> <msup> <mrow> <mo>(</mo> <msup> <mi>G</mi> <mi>H</mi> </msup> <mi>G</mi> <mo>+</mo> <mi>&amp;alpha;</mi> <mi>I</mi> <mo>)</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <msup> <mi>G</mi> <mi>H</mi> </msup> <mi>z</mi> </mrow>
其中α为正则化参数，I代表单位矩阵，上标H代表矩阵共轭转置，由估计到的系数向量可以得到信号包络估计进而提取信号基频分量包络跟踪滤波器可视为时频滤波器，其中心频率为参数化时频变换得到的瞬时频率其带宽为2KF0。
5.根据权利要求1所述的一种旋转设备瞬时转速估计方法，其特征在于，步骤S3中瞬时频率修正过程具体包括：
S301：将步骤S1中估计到的瞬时频率作为步骤S2中包络跟踪滤波器输入，提取信号复数包络具体需建立信号包络傅立叶级数模型，并对模型参数优化求解；
S302：提取包络信号的复数相位对相位函数求导得到包络瞬时频率用权利要求3中的傅立叶级数模型对拟合得到瞬时频率修正量
S303：修正瞬时频率估计
CN201710781137.9A 2017-09-01 2017-09-01 一种旋转设备瞬时转速估计方法 Active CN107783938B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710781137.9A CN107783938B (zh) 2017-09-01 2017-09-01 一种旋转设备瞬时转速估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710781137.9A CN107783938B (zh) 2017-09-01 2017-09-01 一种旋转设备瞬时转速估计方法

Publications (2)

Publication Number Publication Date
CN107783938A true CN107783938A (zh) 2018-03-09
CN107783938B CN107783938B (zh) 2021-04-02

Family

ID=61438038

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710781137.9A Active CN107783938B (zh) 2017-09-01 2017-09-01 一种旋转设备瞬时转速估计方法

Country Status (1)

Country Link
CN (1) CN107783938B (zh)

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108871429A (zh) * 2018-05-09 2018-11-23 宜兴市华鼎机械有限公司 一种两相、三相离心机全息监测方法及***
CN109541590A (zh) * 2018-12-19 2019-03-29 北京科技大学 一种高炉料面点云成像的方法
CN110763462A (zh) * 2019-04-26 2020-02-07 武汉科技大学 一种基于同步压缩算子的时变振动信号故障诊断方法
CN111879508A (zh) * 2020-07-28 2020-11-03 无锡迈斯德智能测控技术有限公司 基于时频变换的旋转机械瞬时转速估计方法、装置和存储介质
CN112746875A (zh) * 2019-10-31 2021-05-04 中国航发商用航空发动机有限责任公司 航空发动机转子轴系复杂振动的主动控制***、方法
CN113252929A (zh) * 2021-07-05 2021-08-13 格创东智(深圳)科技有限公司 转速确定方法、装置、电子设备及计算机可读存储介质
CN114880627A (zh) * 2022-07-08 2022-08-09 西南交通大学 一种基于迭代解调时变滤波的自适应瞬时频率估计方法
CN116312623A (zh) * 2023-03-20 2023-06-23 安徽大学 一种鲸类信号重叠分量的方向脊线预测追踪方法及***

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103020479A (zh) * 2012-12-28 2013-04-03 上海交通大学 基于非线性调频小波变换的信号瞬时频率估计方法
CN105388012A (zh) * 2015-10-22 2016-03-09 兰州理工大学 基于非线性调频小波变换的阶次跟踪方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103020479A (zh) * 2012-12-28 2013-04-03 上海交通大学 基于非线性调频小波变换的信号瞬时频率估计方法
CN105388012A (zh) * 2015-10-22 2016-03-09 兰州理工大学 基于非线性调频小波变换的阶次跟踪方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
J. WANG AND Q. HE: "《Wavelet Packet Envelope Manifold for Fault Diagnosis of Rolling Element Bearings》", 《IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT》 *
Y. YANG ; Z. K. PENG ; G. MENG ; W. M. ZHANG: "《Spline-Kernelled Chirplet Transform for the Analysis of Signals With Time-Varying Frequency and Its Application》", 《IEEE TRANSACTIONS ON INDUSTRIAL ELECTRONICS》 *
杨扬: "《参数化时频分析理论、方法及其在工程信号分析中的应用》", 《中国博士学位论文全文数据库信息科技辑》 *

Cited By (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108871429A (zh) * 2018-05-09 2018-11-23 宜兴市华鼎机械有限公司 一种两相、三相离心机全息监测方法及***
CN109541590A (zh) * 2018-12-19 2019-03-29 北京科技大学 一种高炉料面点云成像的方法
CN109541590B (zh) * 2018-12-19 2020-07-10 北京科技大学 一种高炉料面点云成像的方法
CN110763462A (zh) * 2019-04-26 2020-02-07 武汉科技大学 一种基于同步压缩算子的时变振动信号故障诊断方法
CN110763462B (zh) * 2019-04-26 2023-09-26 武汉科技大学 一种基于同步压缩算子的时变振动信号故障诊断方法
CN112746875B (zh) * 2019-10-31 2022-08-19 中国航发商用航空发动机有限责任公司 航空发动机转子轴系复杂振动的主动控制***、方法
CN112746875A (zh) * 2019-10-31 2021-05-04 中国航发商用航空发动机有限责任公司 航空发动机转子轴系复杂振动的主动控制***、方法
CN111879508A (zh) * 2020-07-28 2020-11-03 无锡迈斯德智能测控技术有限公司 基于时频变换的旋转机械瞬时转速估计方法、装置和存储介质
CN111879508B (zh) * 2020-07-28 2022-06-10 无锡迈斯德智能测控技术有限公司 基于时频变换的旋转机械瞬时转速估计方法、装置和存储介质
CN113252929A (zh) * 2021-07-05 2021-08-13 格创东智(深圳)科技有限公司 转速确定方法、装置、电子设备及计算机可读存储介质
CN113252929B (zh) * 2021-07-05 2022-02-01 格创东智(深圳)科技有限公司 转速确定方法、装置、电子设备及计算机可读存储介质
CN114880627A (zh) * 2022-07-08 2022-08-09 西南交通大学 一种基于迭代解调时变滤波的自适应瞬时频率估计方法
CN114880627B (zh) * 2022-07-08 2022-09-23 西南交通大学 一种基于迭代解调时变滤波的自适应瞬时频率估计方法
CN116312623A (zh) * 2023-03-20 2023-06-23 安徽大学 一种鲸类信号重叠分量的方向脊线预测追踪方法及***
CN116312623B (zh) * 2023-03-20 2023-10-13 安徽大学 一种鲸类信号重叠分量的方向脊线预测追踪方法及***

Also Published As

Publication number Publication date
CN107783938B (zh) 2021-04-02

Similar Documents

Publication Publication Date Title
CN107783938A (zh) 一种旋转设备瞬时转速估计方法
Pan et al. Symplectic geometry mode decomposition and its application to rotating machinery compound fault diagnosis
CN102435844B (zh) 一种频率无关的正弦信号相量计算方法
CN102680948B (zh) 一种线性调频信号调频率和起始频率估计方法
CN106483374A (zh) 一种基于Nuttall双窗全相位FFT的谐波间谐波检测方法
CN107590317A (zh) 一种计及模型参数不确定性的发电机动态估计方法
CN104089774B (zh) 一种基于并行多字典正交匹配的齿轮故障诊断方法
US9645595B2 (en) Method for measuring frequency of phasor of power system
CN103399203B (zh) 一种基于复合迭代算法的谐波参数高精度估计方法
CN105738696A (zh) 全相位时移相位差频率估计方法及装置
CN108871742B (zh) 一种改进的无键相故障特征阶次提取方法
CN101509945B (zh) 正负序电量实时检测的方法
CN108037361A (zh) 一种基于滑动窗dft的高精度谐波参数估计方法
CN103245832A (zh) 基于快速s变换的谐波时频特性参数估计方法及分析仪
CN103618492A (zh) 一种基于时频变换的同步发电机参数辨识方法
CN203287435U (zh) 一种基于stm32f107vct6的微电网谐波与间谐波检测装置
CN107478896A (zh) 一种基于级联广义积分器的频率自适应谐波电流检测方法
CN107807278A (zh) 基于h∞扩展卡尔曼滤波的低频振荡信号参数辨识方法
CN108051189A (zh) 一种旋转机械故障特征提取方法及装置
CN103018044A (zh) 一种改进冲击字典匹配追踪算法的轴承复合故障诊断方法
CN107229597A (zh) 同步挤压广义s变换信号时频分解与重构方法
CN102221639A (zh) 正负序电流实时检测的方法
CN110346772A (zh) 一种高频雷达大幅度电离层相径扰动抑制方法
CN103018555A (zh) 一种高精度的电力参数软件同步采样方法
CN112362343A (zh) 基于调频字典的齿轮箱变转速下分布型故障特征提取方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant