CN105938508A - 一种精确计算振动或压力脉动信号频率及幅值的方法 - Google Patents

Abstract

本发明涉及一种精确计算振动或压力脉动信号频率及幅值的方法，其特征在于它包括以下步骤：1)采用二次曲线插值法获得相对准确的主频频率及其幅值，其步骤如下：(1)选择精确计算的对象，确定初步分析的主频频及其幅值；(2)采用二次插值法确定小山包峰值对应的主频频率f_I及其幅值A_I；2)采用拟合三次多项式精确计算主频频率f及其幅值A，其步骤如下：(1)确定计算频率间隔Δf；(2)计算δf_i和ε；(3)分别按δf_i＝0、δf_i>0和δf_i<0三种情况，精确计算振动或脉动信号的主频频率f及其幅值A，将ε代入相应的拟合三次多项式中，得到精确的主频频率f及其幅值A。本发明提高了计算振动或压力脉动信号频率及幅值的精度。

Description

一种精确计算振动或压力脉动信号频率及幅值的方法

技术领域

本发明涉及一种计算频率及幅值的方法，特别是关于一种在水电站、水泵站运行的水力机械及其他旋转机械的稳定性能测试中使用的精确计算振动或压力脉动信号频率及幅值的方法。

背景技术

目前，在水电站、水泵站运行的水力机械及其他旋转机械的性能测试中，为了解其稳定性而测量其振动值，许多水力机械还要测试过机流体的压力脉动，看其振动或压力脉动的幅值是否超标，检验其运行是否稳定；分析频谱特性，确定其主要频率成分，有针对性的提出避振或减振措施，避免共振发生，降低压力脉动，减轻振动，确保运行稳定。在进行振动或压力脉动信号的频谱分析时，主要方法是快速傅里叶变换(FFT)，把2N个离散数据看作周期函数，用0～(N-1)阶的调谐波进行合成。在转换分析出各阶傅里叶系数之后，再根据其实部和虚部的大小计算出各阶调谐波的幅值及相位，完成对振动或压力脉动信号的频谱分析。

但是，在快速傅立叶变换中，受采样周期T影响，其分析出的频率不是连续的，是有间隔的。例如，采样周期T＝1s，其分析出的频率间隔Δf＝1Hz；也就是说，频率f是1Hz的整数倍，其不可能分析出0.1Hz、11.3Hz或100.5Hz等频率成分；在实际的振动和压力脉动测量中，实际信号的频率不一定是Δf的整数倍，FFT就会用该频率两侧的Δf整数倍频率来合成该频率的信号波形。如某正弦波的频率f＝6.5Hz，采样周期T＝1s，相应的频率间隔Δf＝1Hz，则分析出的最高频率是6Hz，其余依次是7Hz、5Hz、8Hz等等，形成小山包信号，如图1、图2所示。但是，同样是f＝6.5Hz，如采样周期T＝10s，Δf＝0.1Hz，分析的最高幅值频率(或称主频)是f_max＝6.5Hz，和标准值非常一致，如图3、图4所示。

通过对大量调谐波信号的分析对照，我们发现：如果该波形的频率不是分析频率间隔Δf的整数倍，则不能准确的分析出频率值及幅值。这就带来两个问题，一是分析频率不准确，还凭空捏造出一些相近频率；其次是幅值分析不准确，分析获得的一组频率分量的幅值都低于真实振动或压力脉动信号的幅值。现在对振动或压力脉动的评价多采用时域曲线的“峰峰值”，也有采用有效值的，但都不太理想，有部分专家希望采用主频的幅值进行评价，这就要求该主频的幅值必须分析准确，否则无法进行评价。为解决该问题，曾有专家提出过解决办法，主要是延长采样周期T，从而减小FFT分析出的频率间隔Δf，试图以此种方式增加分析频率及其幅值的精度。这样分析出的频率更接近实际，但在某些特殊情况下可能会出现更严重的小山包现象，幅值比实际值低得更多。如某正弦波频率为6.05Hz，其采样周期T＝1s，其波形及频谱分析结果如图5、图6所示，分析的主频f＝6Hz，主频幅值a_max＝0.9919。但是，同样频率为6.05Hz的正弦波，如采样周期T＝10s，分析的主频同样是6Hz，主频幅值却只有0.6392，如图7、图8所示。显然，延长采样周期能起到增加频率分析精度的作用，但在提高幅值分析精度方面作用不确定，个别情况下还可能使分析精度降低，偏离真值的幅度可能更大。

发明内容

针对上述问题，本发明的目的是提供一种精确计算振动或压力脉动信号频率及幅值的方法，其提高了振动或压力脉动信号分析分频频率及其幅值的精度。

为实现上述目的，本发明采取以下技术方案：一种精确计算振动或压力脉动信号频率及幅值的方法，其特征在于它的具体步骤如下：1)第一次修正，采用二次曲线插值法获得相对准确的主频频率及其幅值：(1)在FFT频谱分析曲线中，选择峰值比较高和主频不突出，形状类似于小山包的频率组合作为精确计算的对象，确定初步分析主频频率及其幅值；(2)第一次修正，采用二次曲线插值法确定小山包峰值对应的主频频率f_I及其幅值A_I；2)第二次修正，采用拟合的三次多项式获得精确主频频率f及其幅值A。

优选地，所述步骤2)中，第二次修正的具体步骤如下：(1)确定被分析振动或压力脉动信号的分析频率间隔Δf＝f_i+1-f_i；(2)定义δf_i和ε，即δf_i＝f_I-f_i，ε＝|δf_i/Δf|；(3)分三种情况精确计算振动或脉动信号的主频频率f及其幅值A。

优选地，所述步骤(3)中，三种情况精确计算振动或脉动信号的主频频率f及其幅值A，其具体步骤如下：(3.1)如果δf_i＝0，则二次曲线插值法获得的主频频率及其幅值准确，不必进行第二次修正，即f＝f_I，A＝A_I；(3.2)如果δf_i>0，则主频频率f在f_i右侧，即f>f_i，将ε均代入下面公式，分别得到f右侧部分频率的拟合三次多项式y，f右侧部分频率对应幅值的拟合三次多项式z；

$y=\underset{k=0}{\overset{3}{\&Sigma;}}(b_k\&CenterDot;{\mathrm{\&epsiv;}}^k),Z=\underset{k=0}{\overset{3}{\&Sigma;}}(c_k\&CenterDot;{\mathrm{\&epsiv;}}^k),$

其中，b₀＝0.0535648；b₁＝4.37226；b₂＝-16.1524；b₃＝19.2103；c₀＝0.996647；c₁＝-2.26355；c₂＝6.773067；c₃＝-7.43096；将y和z分别代入下面公式，得到二次修正后振动或压力脉动信号的主频频率f及其幅值A；

f＝f_I+y·Δf，A＝A_I/z；

(3.3)如果δf_i<0，则主频频率f在f_i左侧，即f<f_i，将ε均代入下面公式，分别得到f左侧部分频率的拟合三次多项式g，f左侧部分频率对应幅值的拟合三次多项式h；

$g=\underset{k=0}{\overset{3}{\&Sigma;}}(d_k\&CenterDot;{\mathrm{\&epsiv;}}^k),h=\underset{k=0}{\overset{3}{\&Sigma;}}(e_k\&CenterDot;{\mathrm{\&epsiv;}}^k),$

其中，d₀＝-0.0767339；d₁＝3.74446；d₂＝12.75489；d₃＝14.5654；e₀＝0.994020；e₁＝2.13805；e₂＝5.858876；e₃＝5.79772；将g和h分别代入下面公式，得到二次修正后振动或压力脉动信息的主频频率f及其幅值A；

f＝f_I-g·Δf，A＝A_I/h。

优选地，所述步骤1)的步骤(1)中，小山包需具备如下特征：(a)最高值两侧的频率的幅值A_i-1和A_i+1不低于最高值A_i的1/10；(b)幅值最高频率f_i和前第2个频率f_i-2对应的幅值A_i-2需小于A_i-1，幅值最高频率f_i和后第2个频率f_i+2对应的幅值A_i+2需小于A_i+1。

优选地，所述步骤1)的步骤(2)中，采用二次插值法确定小山包峰值对应的第一次修正主频频率及其幅值，其具体步骤如下：(a)从该频率组合内选择幅值最高的频率为f_i，该频率的幅值为A_i，该频率前后的频率分别为f_i-1和f_i+1，其对应的幅值分别为A_i-1和A_i+1，分别将f_i-1、f_i和f_i+1命名为x₁、x₂和x₃，将A_i-1、A_i和A_i+1命名为y₁、y₂和y₃；(b)假设二次曲线为y＝f(x)＝a₀+a₁·x+a₂·x₂过(x₁,y₁)、(x₂,y₂)和(x₃,y₃)3个点，则将3个点分别代入该方程，得到如下3个等式：

$y_1=a_0+a_1\&CenterDot;x_1+a_2\&CenterDot;x_1^2,y_2=a_0+a_1\&CenterDot;x_2+a_2\&CenterDot;x_2^2,y_3=a_0+a_1\&CenterDot;x_3+a_2\&CenterDot;x_3^2;$

(c)联合解上面三个方程式可得：

$a_1=\frac{y_1-y_2}{x_1-x_2}-a_2\&CenterDot;(x_1-x_2),a_2=\frac{(\frac{y_1-y_2}{x_1-x_2}-\frac{y_1-y_3}{x_1-x_3})}{x_2-x_3};$

(d)对y＝a₀+a₁·x+a₂·x²求一阶导数，并令其等于0得：y′＝a₁-2a₂·x＝0，

得第一次修正后该“小山包”最高点对应的频率

(e)将该频率f_I＝x代入y＝a₀+a₁·x+a₂·x²，得第一次修正后该“小山包”最高点幅值A_I＝y_max。

优选地，所述步骤(3)中，y、z、g、h的4个拟合三次多项式公式是对步骤1)中步骤(2)二次曲线插值结果的再修正和优化，这些拟合三次多项式公式和第一次修正的二次曲线插值是一一对应的。

本发明由于采取以上技术方案，其具有以下优点：1、本发明采用精确计算振动或压力脉动信号频率及幅值的方法，有针对性的提出二次修正方法，能大幅度提高振动或压力脉动信号的频率及幅值分析精度，对需要精确分析频率和幅值来避免产生共振的仪器研发有很大帮助。2、本发明采用精确计算振动或压力脉动信号频率及幅值的方法，发现影响幅值计算精度的不是频率间隔Δf，而是实际频率f与相距最近的计算频率f_i之差的绝对值与频率间隔Δf之比：b＝|f-f_i|÷Δf。3、本发明采用精确计算振动或压力脉动信号频率及幅值的方法，提供的三次多项式进行修正，幅值分析误差可降低到3％以下，频率分辨率可提高一个数量级。

附图说明

图1是现有技术中频率f＝6.5Hz、采样周期T＝1s的调谐波时域曲线示意图；

图2是现有技术中频率f＝6.5Hz、采样周期T＝1s的调谐波频谱分析结果示意图；

图3是现有技术中频率f＝6.5Hz、采样周期T＝10s的调谐波时域曲线示意图；

图4是现有技术中频率f＝6.5Hz、采样周期T＝10s的调谐波频谱分析结果示意图；

图5是现有技术中频率f＝6.05Hz、采样周期T＝1s的调谐波时域曲线示意图；

图6是现有技术中频率f＝6.05Hz、采样周期T＝1s的调谐波频谱分析结果示意图；

图7是现有技术中频率f＝6.05Hz、采样周期T＝10s的调谐波时域曲线示意图。

图8是现有技术中频率f＝6.05Hz、采样周期T＝10s的调谐波频谱分析结果示意图。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

本发明提供一种精确计算振动或压力脉动信号频率及幅值的方法，其具体步骤如下：

1)第一次修正：采用二次曲线插值法获得相对准确的主频频率及其幅值；

(1)在FFT频谱分析曲线中，选择峰值比较高和主频不突出，形状类似于小山包的频率组合作为精确计算的对象；

(2)采用二次曲线插值法确定小山包峰值对应的主频频率f_I及其幅值A_I，其具体步骤如下：

(a)从该频率组合内选择幅值最高的频率为f_i，该频率的幅值为A_i，该频率前后的频率分别为f_i-1和f_i+1，其对应的幅值分别为A_i-1和A_i+1，分别将f_i-1、f_i和f_i+1命名为x₁、x₂和x₃，将A_i-1、A_i和A_i+1命名为y₁、y₂和y₃；

(b)假设二次曲线为y＝f(x)＝a₀+a₁·x+a₂·x₂过(x₁,y₁)、(x₂,y₂)和(x₃,y₃)3个点，则将3个点分别代入该方程。得到如下3个等式：

$y_1=a_0+a_1\&CenterDot;x_1+a_2\&CenterDot;x_1^2---\left(1\right)$

$y_2=a_0+a_1\&CenterDot;x_2+a_2\&CenterDot;x_2^2---\left(2\right)$

$y_3=a_0+a_1\&CenterDot;x_3+a_2\&CenterDot;x_3^2---\left(3\right)$

(c)联合解(1)、(2)、(3)三个方程式可得：

$a_0=y_1-a_1\&CenterDot;x_1-a_2\&CenterDot;x_1^2---\left(4\right)$

$a_1=\frac{y_1-y_2}{x_1-x_2}-a_2\&CenterDot;(x_1-x_2)---\left(5\right)$

$a_2=\frac{(\frac{y_1-y_2}{x_1-x_2}-\frac{y_1-y_3}{x_1-x_3})}{x_2-x_3}---\left(6\right)$

(d)对y＝a₀+a₁·x+a₂·x²求一阶导数，并令其等于0：

y′＝a₁-2a₂·x＝0 (7)

可得该小山包最高点对应的第一次修正主频频率

(e)将该频率f_I＝x代入y＝a₀+a₁·x+a₂·x²，得该小山包最高点主频频率对应的幅值A_I＝y_max；

2)第二次修正：采用拟合的三次多项式获得精确主频频率及其幅值

(1)确定被分析振动或压力脉动信号的分析频率间隔Δf＝f_i+1-f_i；

(2)定义δf_i和ε，即δf_i＝f_I-f_i，ε＝|δf_i/Δf|；

(3)分三种情况精确计算振动或脉动信号的主频频率f及其幅值A：

(3.1)如果δf_i＝0，则二次曲线插值法获得的主频频率及其幅值准确，不必进行第二次修正，即f＝f_I，A＝A_I；

(3.2)如果δf_i>0，则主频频率f在f_i右侧，即f>f_i，将ε代入公式(8)，得到f右侧部分频率的拟合三次多项式y；将ε代入公式(9)，得到f_j右侧部分频率对应幅值的拟合三次多项式z：

$y=\underset{k=0}{\overset{3}{\&Sigma;}}(b_k\&CenterDot;{\mathrm{\&epsiv;}}^k)---\left(8\right)$

$z=\underset{k=0}{\overset{3}{\&Sigma;}}(c_k\&CenterDot;{\mathrm{\&epsiv;}}^k)---\left(9\right)$

其中，b₀＝0.0535648；b₁＝4.37226；b₂＝-16.1524；b₃＝19.2103；c₀＝0.996647；c₁＝-2.26355；c₂＝6.773067；c₃＝-7.43096。

将y和z分别代入公式(10)和(11)，得到二次修正后振动或压力脉动信号的主频频率f及其幅值A；

f＝f_I+y·Δf (10)

A＝A_I/z (11)

(3.3)如果δf_i<0，则主频频率f在f_i左侧，即f<f_i，将ε代入公式(12)，得到f_j左侧部分频率的拟合三次多项式g；将ε代入公式(13)，得到f_j左侧部分频率对应幅值的拟合三次多项式h：

$g=\underset{k=0}{\overset{3}{\&Sigma;}}(d_k\&CenterDot;{\mathrm{\&epsiv;}}^k)---\left(12\right)$

$h=\underset{k=0}{\overset{3}{\&Sigma;}}(e_k\&CenterDot;{\mathrm{\&epsiv;}}^k)---\left(13\right)$

其中，d₀＝-0.0767339；d₁＝3.74446；d₂＝12.75489；d₃＝14.5654；e₀＝0.994020；e₁＝2.13805；e₂＝5.858876；e₃＝5.79772。

将g和h分别代入(14)和(15)，得到二次修正后振动或压力脉动信息的主频频率f及其幅值A。

f＝f_I-g·Δf (14)

A＝A_I/h (15)

上述步骤1)中步骤(1)中，小山包需具备如下特征：

(a)最高值两侧的频率的幅值A_i-1和A_i+1不低于最高值A_i的1/10；

(b)幅值最高频率(或称初步分析主频)f_i和前第2个频率f_i-2对应的幅值A_i-2需小于A_i-1，幅值最高频率f_i和其后第2个频率f_i+2对应的幅值A_i+2需小于A_i+1。

上述步骤2)的步骤(3)中，y、z、g、h的4个拟合三次多项式公式是对步骤1)二次曲线插值结果的再修正和优化，这些拟合三次多项式公式和第一次修正的二次曲线插值是一一对应的，不可更改。

上述各实施例仅用于说明本发明，各部件的结构、尺寸、设置位置及形状都是可以有所变化的，在本发明技术方案的基础上，凡根据本发明原理对个别部件进行的改进和等同变换，均不应排除在本发明的保护范围之外。

Claims (6)

1.一种精确计算振动或压力脉动信号频率及幅值的方法，其特征在于：它的具体步骤如下：

1)第一次修正，采用二次曲线插值法获得相对准确的主频频率及其幅值：

(1)在FFT频谱分析曲线中，选择峰值比较高和主频不突出，形状类似于小山包的频率组合作为精确计算的对象，确定初步分析主频频率及其幅值；

(2)第一次修正，采用二次曲线插值法确定小山包峰值对应的主频频率f_I及其幅值A_I；

2)第二次修正，采用拟合的三次多项式获得精确主频频率f及其幅值A。

2.如权利要求1所述的一种精确计算振动或压力脉动信号频率及幅值的方法，其特征在于：所述步骤2)中，第二次修正的具体步骤如下：

(1)确定被分析振动或压力脉动信号的分析频率间隔Δf＝f_i+1-f_i；

(2)定义δf_i和ε，即δf_i＝f_I-f_i，ε＝|δf_i/Δf|；

(3)分三种情况精确计算振动或脉动信号的主频频率f及其幅值A。

3.如权利要求2所述的一种精确计算振动或压力脉动信号频率及幅值的方法，其特征在于：所述步骤(3)中，三种情况精确计算振动或脉动信号的主频频率f及其幅值A，其具体步骤如下：

(3.1)如果δf_i＝0，则二次曲线插值法获得的主频频率及其幅值准确，不必进行第二次修正，即f＝f_I，A＝A_I；

(3.2)如果δf_i>0，则主频频率f在f_i右侧，即f>f_i，将ε均代入下面公式，分别得到f右侧部分频率的拟合三次多项式y，f右侧部分频率对应幅值的拟合三次多项式z；

$y=\underset{k=0}{\overset{3}{\&Sigma;}}(b_k\&CenterDot;{\mathrm{\&epsiv;}}^k),$

$Z=\underset{k=0}{\overset{3}{\&Sigma;}}(c_k\&CenterDot;{\mathrm{\&epsiv;}}^k),$

其中，b₀＝0.0535648；b₁＝4.37226；b₂＝-16.1524；b₃＝19.2103；c₀＝0.996647；c₁＝-2.26355；c₂＝6.773067；c₃＝-7.43096；

将y和z分别代入下面公式，得到二次修正后振动或压力脉动信号的主频频率f及其幅值A；

f＝f_I+y·Δf，

A＝A_I/z；

(3.3)如果δf_i<0，则主频频率f在f_i左侧，即f<f_i，将ε均代入下面公式，分别得到f左侧部分频率的拟合三次多项式g，f左侧部分频率对应幅值的拟合三次多项式h；

$g=\underset{k=0}{\overset{3}{\&Sigma;}}(d_k\&CenterDot;{\mathrm{\&epsiv;}}^k),$

$h=\underset{k=0}{\overset{3}{\&Sigma;}}(e_k\&CenterDot;{\mathrm{\&epsiv;}}^k),$

其中，d₀＝-0.0767339；d₁＝3.74446；d₂＝12.75489；d₃＝14.5654；e₀＝0.994020；e₁＝2.13805；e₂＝5.858876；e₃＝5.79772；

将g和h分别代入下面公式，得到二次修正后振动或压力脉动信息的主频频率f及其幅值A；

f＝f_I-g·Δf，

A＝A_I/h。

4.如权利要求1所述的一种精确计算振动或压力脉动信号频率及幅值的方法，其特征在于：所述步骤1)的步骤(1)中，小山包需具备如下特征：

(a)最高值两侧的频率的幅值A_i-1和A_i+1不低于最高值A_i的1/10；

(b)幅值最高频率f_i和前第2个频率f_i-2对应的幅值A_i-2需小于A_i-1，幅值最高频率f_i和后第2个频率f_i+2对应的幅值A_i+2需小于A_i+1。

5.如权利要求1所述的一种精确计算振动或压力脉动信号频率及幅值的方法，其特征在于：所述步骤1)的步骤(2)中，采用二次插值法确定小山包峰值对应的第一次修正主频频率及其幅值，其具体步骤如下：

(a)从该频率组合内选择幅值最高的频率为f_i，该频率的幅值为A_i，该频率前后的频率分别为f_i-1和f_i+1，其对应的幅值分别为A_i-1和A_i+1，分别将f_i-1、f_i和f_i+1命名为x₁、x₂和x₃，将A_i-1、A_i和A_i+1命名为y₁、y₂和y₃；

(b)假设二次曲线为y＝f(x)＝a₀+a₁·x+a₂·x₂过(x₁,y₁)、(x₂,y₂)和(x₃,y₃)3个点，则将3个点分别代入该方程，得到如下3个等式：

$y_1=a_0+a_1\&CenterDot;x_1+a_2\&CenterDot;x_1^2,$

$y_2=a_0+a_1\&CenterDot;x_2+a_2\&CenterDot;x_2^2,$

$y_3=a_0+a_1\&CenterDot;x_3+a_2\&CenterDot;x_3^2;$

(c)联合解上面三个方程式可得：

$a_0=y_1-a_1\&CenterDot;x_1-a_2\&CenterDot;x_1^2,$

$a_1=\frac{y_1-y_2}{x_1-x_2}-a_2\&CenterDot;(x_1-x_2),$

$a_2=\frac{(\frac{y_1-y_2}{x_1-x_2}-\frac{y_1-y_3}{x_1-x_3})}{x_2-x_3};$

(d)对y＝a₀+a₁·x+a₂·x²求一阶导数，并令其等于0得：

y′＝a₁-2a₂·x＝0，

得第一次修正后该“小山包”最高点对应的频率

(e)将该频率f_I＝x代入y＝a₀+a₁·x+a₂·x²，得第一次修正后该“小山包”最高点幅值A_I＝y_max。

6.如权利要求1-5任一项所述的一种精确计算振动或压力脉动信号频率及幅值的方法，其特征在于：所述步骤(3)中，y、z、g、h的4个拟合三次多项式公式是对步骤1)中步骤(2)二次曲线插值结果的再修正和优化，这些拟合三次多项式公式和第一次修正的二次曲线插值是一一对应的。