一种精确计算振动或压力脉动信号频率及幅值的方法
【专利摘要】本发明涉及一种精确计算振动或压力脉动信号频率及幅值的方法，其特征在于它包括以下步骤：1)采用二次曲线插值法获得相对准确的主频频率及其幅值，其步骤如下：(1)选择精确计算的对象，确定初步分析的主频频及其幅值；(2)采用二次插值法确定小山包峰值对应的主频频率fI及其幅值AI；2)采用拟合三次多项式精确计算主频频率f及其幅值A，其步骤如下：(1)确定计算频率间隔Δf；(2)计算δfi和ε；(3)分别按δfi＝0、δfi>0和δfi<0三种情况，精确计算振动或脉动信号的主频频率f及其幅值A，将ε代入相应的拟合三次多项式中，得到精确的主频频率f及其幅值A。本发明提高了计算振动或压力脉动信号频率及幅值的精度。
【专利说明】
一种精确计算振动或压力脉动信号频率及幅值的方法
技术领域
[0001] 本发明涉及一种计算频率及幅值的方法，特别是关于一种在水电站、水栗站运行 的水力机械及其他旋转机械的稳定性能测试中使用的精确计算振动或压力脉动信号频率 及幅值的方法。
【背景技术】
[0002] 目前，在水电站、水栗站运行的水力机械及其他旋转机械的性能测试中，为了解其 稳定性而测量其振动值，许多水力机械还要测试过机流体的压力脉动，看其振动或压力脉 动的幅值是否超标，检验其运行是否稳定;分析频谱特性，确定其主要频率成分，有针对性 的提出避振或减振措施，避免共振发生，降低压力脉动，减轻振动，确保运行稳定。在进行振 动或压力脉动信号的频谱分析时，主要方法是快速傅里叶变换(FFT)，把2N个离散数据看作 周期函数，用〇~(N-I)阶的调谐波进行合成。在转换分析出各阶傅里叶系数之后，再根据其 实部和虚部的大小计算出各阶调谐波的幅值及相位，完成对振动或压力脉动信号的频谱分 析。
[0003] 但是，在快速傅立叶变换中，受采样周期T影响，其分析出的频率不是连续的，是有 间隔的。例如，采样周期T = Is，其分析出的频率间隔Δ f = IHz ;也就是说，频率f是IHz的整 数倍，其不可能分析出〇. IHz、11.3Hz或100.5Hz等频率成分;在实际的振动和压力脉动测量 中，实际信号的频率不一定是A f的整数倍，FFT就会用该频率两侧的△ f整数倍频率来合成 该频率的信号波形。如某正弦波的频率f = 6.5Hz，采样周期T= Is，相应的频率间隔Δ f = IHz，则分析出的最高频率是6Hz，其余依次是7Hz、5Hz、8Hz等等，形成小山包信号，如图1、图 2所示。但是，同样是f = 6.5Hz，如采样周期T = IOs，Af = 〇.IHz，分析的最高幅值频率(或称 主频)是f max = 6.5Hz，和标准值非常一致，如图3、图4所示。
[0004] 通过对大量调谐波信号的分析对照，我们发现:如果该波形的频率不是分析频率 间隔Af的整数倍，则不能准确的分析出频率值及幅值。这就带来两个问题，一是分析频率 不准确，还凭空捏造出一些相近频率;其次是幅值分析不准确，分析获得的一组频率分量的 幅值都低于真实振动或压力脉动信号的幅值。现在对振动或压力脉动的评价多采用时域曲 线的"峰峰值"，也有采用有效值的，但都不太理想，有部分专家希望采用主频的幅值进行评 价，这就要求该主频的幅值必须分析准确，否则无法进行评价。为解决该问题，曾有专家提 出过解决办法，主要是延长采样周期T，从而减小FFT分析出的频率间隔Af，试图以此种方 式增加分析频率及其幅值的精度。这样分析出的频率更接近实际，但在某些特殊情况下可 能会出现更严重的小山包现象，幅值比实际值低得更多。如某正弦波频率为6.05Hz，其采样 周期T=ls，其波形及频谱分析结果如图5、图6所示，分析的主频f = 6Hz，主频幅值amax = 0.9919。但是，同样频率为6.05Hz的正弦波，如采样周期T = IOs，分析的主频同样是6Hz，主 频幅值却只有0.6392,如图7、图8所示。显然，延长采样周期能起到增加频率分析精度的作 用，但在提高幅值分析精度方面作用不确定，个别情况下还可能使分析精度降低，偏离真值 的幅度可能更大。
【发明内容】
[0005] 针对上述问题，本发明的目的是提供一种精确计算振动或压力脉动信号频率及幅 值的方法，其提高了振动或压力脉动信号分析分频频率及其幅值的精度。
[0006] 为实现上述目的，本发明采取以下技术方案:一种精确计算振动或压力脉动信号 频率及幅值的方法，其特征在于它的具体步骤如下：1)第一次修正，采用二次曲线插值法获 得相对准确的主频频率及其幅值：（1)在FFT频谱分析曲线中，选择峰值比较高和主频不突 出，形状类似于小山包的频率组合作为精确计算的对象，确定初步分析主频频率及其幅值； (2)第一次修正，采用二次曲线插值法确定小山包峰值对应的主频频率fLS其幅值AI;2)第 二次修正，采用拟合的三次多项式获得精确主频频率f及其幅值A。
[0007] 优选地，所述步骤2)中，第二次修正的具体步骤如下：（1)确定被分析振动或压力 脉动信号的分析频率间隔Δ f = fi+i-fi; (2)定义Sfi和ε，即= I 5fi/A f I ; (3)分 三种情况精确计算振动或脉动信号的主频频率f及其幅值A。
[0008] 优选地，所述步骤(3)中，三种情况精确计算振动或脉动信号的主频频率f及其幅 值A，其具体步骤如下：（3.1)如果Sf 1 = O,则二次曲线插值法获得的主频频率及其幅值准 确，不必进行第二次修正，即f = fi，A=Ai; (3.2)如果5fi>〇,则主频频率f在fi右侧，即f>fi， 将ε均代入下面公式，分别得到f右侧部分频率的拟合三次多项式y，f右侧部分频率对应幅 值的拟合三次多项式z;
[0009]
[0010] 其中，b〇 = 0.0535648 ;bi = 4.37226 ;b2 = -16.1524 ;b3 = 19.2103; CO = O .996647; Cl =-2.26355; c2 = 6.773067; C3 = -7.43096;将y和z分别代入下面公式，得到二次修正后振动 或压力脉动信号的主频频率f及其幅值A;
[0011] f = fi+y · Af,A=Ai/z；
[0012] (3.3)如果δ?^ο,则主频频率ftEfi左侧，即Kf1，将￡均代入下面公式，分别得到f 左侧部分频率的拟合三次多项式R，f左侧部分频率对应幅值的拟合三次多项式h;
[0013]
[0014] 其中，do = -。· 0767339 ;di = 3.74446; Cl2 = 12.75489 ;d3 = 14.5654; e〇 = 0.994020; ei = 2 · 13805; e2 = 5 · 858876; e3 = 5 · 79772;将g和h分别代入下面公式，得到二次修正后振动 或压力脉动信息的主频频率f及其幅值A;
[0015] f = fi_g · Δ f，A=Ai/h〇
[0016] 优选地，所述步骤1)的步骤（1)中，小山包需具备如下特征：（a)最高值两侧的频率 的幅值Ai-4PAi+1不低于最高值A i的l/10;(b)幅值最高频率fi和前第2个频率fV2对应的幅值 Ai-2需小于A^1，幅值最高频率f i和后第2个频率f i+2对应的幅值Ai+2需小于Ai+1。
[0017] 优选地，所述步骤1)的步骤(2)中，采用二次插值法确定小山包峰值对应的第一次 修正主频频率及其幅值，其具体步骤如下：（a)从该频率组合内选择幅值最高的频率为h， 该频率的幅值为Ai，该频率前后的频率分别为fi-1和fi+i，其对应的幅值分别为Ai-1和Ai+i，分 别将fVi、fi和fi+i命名为xi、X2和X3,将Ai-i、Ai和Ai+i命名为yi、y2和y3; (b)假设二次曲线为y = f(x)=ao+ai · x+a2 · X2过(11，71)、（12,72)和(13,73)3个点，则将3个点分别代入该方程，得 到如下3个等式：
[0018] "
'
[0019] (c)联合解上面三个方程式可得:4:fl =为-q·.' -α2?〇7'，
[0020]
[0021 ] (d)对y = a〇+ai · x+a2 · X2求一阶导数，并令其等于0得:y' =ai_2a2 · x = 0，
[0022] 得第一次修正后该"小山包"最高点对应的频率
[0023] (e)将该频率fi = x代入y = ao+ai · x+a2 · X2,得第一次修正后该"小山包"最高点幅 {直 Al - Ymax 〇
[0024] 优选地，所述步骤(3)中，y、z、g、h的4个拟合三次多项式公式是对步骤1)中步骤 (2)二次曲线插值结果的再修正和优化，这些拟合三次多项式公式和第一次修正的二次曲 线插值是一一对应的。
[0025]本发明由于采取以上技术方案，其具有以下优点：1、本发明采用精确计算振动或 压力脉动信号频率及幅值的方法，有针对性的提出二次修正方法，能大幅度提高振动或压 力脉动信号的频率及幅值分析精度，对需要精确分析频率和幅值来避免产生共振的仪器研 发有很大帮助。2、本发明采用精确计算振动或压力脉动信号频率及幅值的方法，发现影响 幅值计算精度的不是频率间隔A f，而是实际频率f与相距最近的计算频率A之差的绝对值 与频率间隔Δ f之比:b = I ^f11 + Δ f。3、本发明采用精确计算振动或压力脉动信号频率及 幅值的方法，提供的三次多项式进行修正，幅值分析误差可降低到3%以下，频率分辨率可 提高一个数量级。
【附图说明】
[0026]图1是现有技术中频率f = 6.5Hz、采样周期T = IS的调谐波时域曲线示意图；
[0027]图2是现有技术中频率f = 6.5Hz、采样周期T = Is的调谐波频谱分析结果示意图；
[0028I图3是现有技术中频率f = 6.5Hz、采样周期T = I Os的调谐波时域曲线示意图；
[0029]图4是现有技术中频率f = 6.5Hz、采样周期T = IOs的调谐波频谱分析结果示意图； [0030]图5是现有技术中频率f = 6.05Hz、采样周期T = IS的调谐波时域曲线示意图；
[0031]图6是现有技术中频率f = 6.05Hz、采样周期T = Is的调谐波频谱分析结果示意图； [0032 ]图7是现有技术中频率f = 6.05Hz、采样周期T = I Os的调谐波时域曲线示意图。 [0033]图8是现有技术中频率f = 6.05Hz、采样周期T = IOs的调谐波频谱分析结果示意 图。
【具体实施方式】
[0034]下面结合附图和实施例对本发明进行详细的描述。
[0035]本发明提供一种精确计算振动或压力脉动信号频率及幅值的方法，其具体步骤如 下：
[0036] 1)第一次修正:采用二次曲线插值法获得相对准确的主频频率及其幅值；
[0037] (1)在FFT频谱分析曲线中，选择峰值比较高和主频不突出，形状类似于小山包的 频率组合作为精确计算的对象；
[0038] (2)采用二次曲线插值法确定小山包峰值对应的主频频率h及其幅值六:，其具体步 骤如下：
[0039 ] (a)从该频率组合内选择幅值最高的频率为fi，该频率的幅值为Ai，该频率前后的 频率分别为fi-1和fi+ι，其对应的幅值分别为Ai-I和Ai+ι，分别将fi-i、fi和fi+ι命名为X1、X2和 叉3，将厶1-1^和八1+1命名为71、72和73;
[0040] (b)假设二次曲线为y = f(x)=a〇+ai · x+a2 · X2过（xi，yi)、（X2，y2WP(X3，y3)3f 点，则将3个点分别代入该方程。得到如下3个等式：
[0048] (d)对y = ao+ai · x+a2 · X2求一阶导数，并令其等于0:[0049] =ai~2a2 · x = 0 (7)
[0041 ] (1)
[0042] (2)
[0043] (3)
[0044]
[0045] (4)
[0046] (.5) 则 (§)
[0058] (8)
[0059] (9)
[0060] 其中，b〇 = 0.0535648 ;bi = 4.37226 ;b2 = -16.1524 ;b3 = 19.2103; Co = O .996647; ci =-2.26355; C2 = 6.773067; c3 = -7.43096。
[0061] 将y和z分别代入公式（10)和（11)，得到二次修正后振动或压力脉动信号的主频频 率f及其幅值A;
[0062] f = fi+y · Δ f (10)
[0063] A=Ai/z (11)
[0064] (3.3)如果6心〈〇，则主频频率€在乜左侧，8时〈心，将￡代入公式（12)，得到心左侧部 分频率的拟合三次多项式8;将 ￡代入公式（13)，得到fj左侧部分频率对应幅值的拟合三次 多项式h:
[0065] (12)
[0066] (13)
[0067] 其中，d〇 = -〇 .0767339 ;di = 3.74446; Cl2 = 12.75489 ;d3 = 14.5654; e〇 = 0.994020; ei = 2.13805 ； e2 = 5.858876 ；e3 = 5.79772〇
[0068] 将g和h分别代入（14)和（15)，得到二次修正后振动或压力脉动信息的主频频率f 及其幅值A。
[0069] f = fi-g · Af (14)
[0070] A=Ai/h (15)
[0071] 上述步骤1)中步骤(1)中，小山包需具备如下特征：
[0072] (a)最高值两侧的频率的幅值Ai-dPAi+1不低于最高值Ai的1/10;
[0073] (b)幅值最高频率(或称初步分析主频)fi和前第2个频率fV2对应的幅值Ai- 2需小 于Au，幅值最高频率fi和其后第2个频率fi+2对应的幅值Ai+2需小于A i+1。
[0074] 上述步骤2)的步骤(3)中，y、z、g、h的4个拟合三次多项式公式是对步骤1)二次曲 线插值结果的再修正和优化，这些拟合三次多项式公式和第一次修正的二次曲线插值是一 一对应的，不可更改。
[0075]上述各实施例仅用于说明本发明，各部件的结构、尺寸、设置位置及形状都是可以 有所变化的，在本发明技术方案的基础上，凡根据本发明原理对个别部件进行的改进和等 同变换，均不应排除在本发明的保护范围之外。
【主权项】
1. 一种精确计算振动或压力脉动信号频率及幅值的方法，其特征在于：它的具体步骤 如下： 1) 第一次修正，采用二次曲线插值法获得相对准确的主频频率及其幅值： (1) 在FFT频谱分析曲线中，选择峰值比较高和主频不突出，形状类似于小山包的频率 组合作为精确计算的对象，确定初步分析主频频率及其幅值； (2) 第一次修正，采用二次曲线插值法确定小山包峰值对应的主频频率fi及其幅值Ai; 2) 第二次修正，采用拟合的Ξ次多项式获得精确主频频率f及其幅值A。2. 如权利要求1所述的一种精确计算振动或压力脉动信号频率及幅值的方法，其特征 在于:所述步骤2)中，第二次修正的具体步骤如下： (1) 确定被分析振动或压力脉动信号的分析频率间隔Af = fw-fi; (2) 定义 Sfi 和 ε，即 Sfi = f 广 fi，e= |Sfi/Af I ; (3) 分Ξ种情况精确计算振动或脉动信号的主频频率f及其幅值A。3. 如权利要求2所述的一种精确计算振动或压力脉动信号频率及幅值的方法，其特征 在于:所述步骤(3)中，Ξ种情况精确计算振动或脉动信号的主频频率f及其幅值A，其具体 步骤如下： (3.1) 如果Sfi = 〇,则二次曲线插值法获得的主频频率及其幅值准确，不必进行第二次 修正，即 f = fi，A=Ai; (3.2) 如果δfl〉0，则主频频率f在fl右侧，即f乂l，将ε均代入下面公式，分别得到f右侧 部分频率的拟合Ξ次多项式y，f右侧部分频率对应幅值的拟合Ξ次多项式z;其中，6〇 = 0.0535648;61 = 4.37226;62 = -16.1524;63=19.2103;。〇 = 0.996647;。1 =- 2.26355; C2 = 6.773067; C3 = -7.43096; 将y和Z分别代入下面公式，得到二次修正后振动或压力脉动信号的主频频率f及其幅 值A; f = fi+y · A f， A=Ai/z； (3.3) 如果δfl<0，则主频频率f在fl左侧，即f<fl，将ε均代入下面公式，分别得到f左侧 部分频率的拟合Ξ次多项式g，f左侧部分频率对应幅值的拟合Ξ次多项式h;其中，d〇 = -〇. 0767339 ;di = 3.74446 ;d2 = 12.75489 ;d3 = 14.5654; e〇 = 0.994020 ;ei = 2.13805 ； Θ2 = 5.858876 ； es = 5.79772 ； 将g和h分别代入下面公式，得到二次修正后振动或压力脉动信息的主频频率f及其幅 值A; f = f广邑· A f， A=Ai/h。4. 如权利要求1所述的一种精确计算振动或压力脉动信号频率及幅值的方法，其特征 在于:所述步骤1)的步骤(1)中，小山包需具备如下特征： (a) 最高值两侧的频率的幅值Ai-i和Ai+i不低于最高值Ai的1/10; (b) 幅值最高频率fi和前第2个频率fi-拥应的幅值Ai-2需小于Ai-i，幅值最高频率fi和后 第2个频率fi+拥应的幅值Ai+2需小于Ai+i。5. 如权利要求1所述的一种精确计算振动或压力脉动信号频率及幅值的方法，其特征 在于:所述步骤1)的步骤(2)中，采用二次插值法确定小山包峰值对应的第一次修正主频频 率及其幅值，其具体步骤如下： (a) 从该频率组合内选择幅值最高的频率为fi，该频率的幅值为Ai，该频率前后的频率 分别为fi-l和fi+l，其对应的幅值分别为Ai-l和Ai+l，分别将和fi+l命名为X1、X2和X3,将 Ai-i、Ai 和 Ai+i 命名为 yi、y2 和 y3; (b) 假设二次曲线为y = f(x)=ao+ai · x+a2 · X2过(又1，71)、佔，72)和(又3,73)3个点，则将 3个点分别代入该方程，得到如下3个等式：(C)联合解上面Ξ个方程式可得：(d) 对y = ao+ai · x+a2 · X2求一阶导数，并令其等于0得： y' =a广2曰2 · x = 0， 得第一次修正后该"小山包"最高点对应的频率(e) 将该频率fi = x代入y = ao+ai · x+a2 . X2,得第一次修正后该"小山包'最高点幅值Ai -y max 〇6. 如权利要求1-5任一项所述的一种精确计算振动或压力脉动信号频率及幅值的方 法，其特征在于:所述步骤(3)中，y、z、g、h的4个拟合Ξ次多项式公式是对步骤1)中步骤(2) 二次曲线插值结果的再修正和优化，运些拟合Ξ次多项式公式和第一次修正的二次曲线插 值是--对应的。
【文档编号】G06F17/50GK105938508SQ201610318735
【公开日】2016年9月14日
【申请日】2016年5月13日
【发明人】徐洪泉, 吴剑, 王万鹏, 廖翠林, 李铁友, 陆力, 瞿军, 易艳林, 赵立策, 王武昌, 范小付
【申请人】中国水利水电科学研究院, 中国水利学会, 北京中水科水电科技开发有限公司