一种基于时变滤波的旋转机械信号时频分析方法

1.本发明涉及旋转机械装备的状态监测及故障诊断技术领域，特别是指一种基于时变滤波的旋转机械信号时频分析方法。


背景技术：

2.旋转机械广泛应用在交通、能源、工业等领域，通过振动、电流、应力应变等动态运行信号对其进行状态监测和故障诊断，对设备运行维护至关重要。旋转机械中各个旋转部件(包括轴、齿轮、轴承等)的特征频率通常与输入转频成比例变化。旋转机械的运行工况往往存在较强的非稳定性，导致其信号中的各频率成分具有时变特征，且在频域中相互重叠，此时信号无法简单通过时域或频域分析提取特征。时频分析方法可以在时间、频率、幅值(或功率)三方面全面揭示非平稳信号的时变特征，对旋转机械状态监测与故障诊断具有重要意义。
3.传统时频分析方法(例如，短时fourier变换、连续小波变换等)存在时频分辨率较低、存在交叉项干扰(wigner-ville分布)等缺陷，在信号成分频率接近或噪声较强时难以准确显示微弱瞬态特征。希尔伯特-huang变换通过自适应模式分解方法将多分量信号分解为多个单分量信号，再通过分别计算单分量信号的瞬时频率和包络线来构造原信号的时频分布。此类时频分布具有很好的时频分辨率，且不受交叉项干扰，能够反映信号的瞬态特性，适合复杂多分量非平稳信号的时频分析。常用的自适应模式分解方法包括经验模式分解(emd)、局部均值分解(lmd)、本质时间尺度分解(itd)、变分模式分解(vmd)等。然而，在分解频率变化范围较大或存在频率中断的信号时，自适应模式分解方法存在模式混淆问题，具有相同物理意义的同一频率分量可能被分解到不同的单分量信号中，使得分解得到的信号无法进一步揭示原信号的时频分布。模式混淆问题导致基于自适应模式分解的时频分析方法鲁棒性较差，限制了其在工程问题中的应用。


技术实现要素：

4.本发明实施例提供了基于时变滤波的旋转机械信号时频分析方法，能够克服时频模糊、交叉项干扰、模式混淆问题，得到具有较高时频分辨率的时频分布。
5.本发明实施例提供的基于时变滤波的旋转机械信号时频分析方法，包括：
6.获取非平稳工况下目标旋转机械的时域信号作为原信号；
7.确定原信号中各单分量信号的时频曲线；
8.通过时变滤波在时域提取各时频曲线对应的单分量信号；
9.计算提取的各单分量信号的希尔伯特时频分布，将各希尔伯特时频分布叠加得到原信号的时频分布。
10.其中，所述时域信号包括：振动、位移、声音、应变、压力或电信号。
11.其中，所述确定原信号中各单分量信号的时频曲线包括：
12.获取目标旋转机械的输入转频；
13.计算旋转机械各特征频率相对输入转频的特征阶次；
14.将输入转频与特征阶次相乘，得到原信号中各单分量信号的时频曲线。
15.其中，所述确定原信号中各单分量信号的时频曲线包括：
16.对原信号进行短时傅里叶变换，得到原信号的传统时频分布图；
17.基于得到的传统时频分布图，估计各时频脊线峰值对应的时频曲线，得到原信号中各单分量信号的时频曲线。
18.其中，所述确定原信号中各单分量信号的时频曲线包括：
19.获取目标旋转机械的输入转频；
20.对原信号进行短时傅里叶变换，得到原信号的传统时频分布图；
21.对于与输入转频成正比的频率成分，计算特征阶次，将输入转频与特征阶次相乘估算其时频曲线；
22.对于其它与输入转频无关的频率成分，基于得到的传统时频分布图，估计各时频脊线峰值对应的时频曲线。
23.其中，所述通过时变滤波在时域提取各时频曲线对应的单分量信号包括：
24.a1，根据得到的时频曲线，确定时变滤波的上截止频率f
cu
(t)和下截止频率f
cl
(t)：
[0025][0026]
其中，t表示时间，sk(t)表示第k个特征阶次的时频曲线，b(t)表示时变滤波带宽；
[0027]
a2，根据得到的时变滤波的上截止频率f
cu
(t)和下截止频率f
cl
(t)，确定局部滤波器长度ln：
[0028]
其中，ρ
max
≠0
[0029][0030]
其中，n＝0,1,
…
，n-1，n为时域信号的长度，l表示滤波器长度参数，α表示变化系数，ρ(n)表示局部截止频率的导数平方和，ρ
max
表示ρ(n)的最大值，f
cu
(n)表示时变滤波的上截止频率，f
cl
(n)表示时变滤波的下截止频率，t(n)表示第n个时间点，表示向下取整；
[0031]
其中，l表示滤波器长度参数，α表示变化系数，n＝0,1,
…
，n-1，n为时域信号的长度；
[0032]
a3，根据得到的上截止频率f
cu
(t)、下截止频率f
cl
(t)和局部滤波器长度ln，构造滤波器矩阵w；
[0033]
a4，使用滤波器矩阵对原信号进行时变滤波得到单分量信号：
[0034]
y＝wx
[0035]
其中，x为原信号，y为单分量信号；
[0036]
a5，重复步骤a1-a4，得到各时频曲线对应的单分量信号。
[0037]
其中，滤波器矩阵w表示为：
[0038][0039]
其中，w
(n)
表示t(n)处的局部滤波器，且表示t(n)处的局部滤波器，且为局部滤波器w
(n)
的滤波器系数，局部滤波器长度等于2ln+1。
[0040]
其中，所述计算提取的各单分量信号的希尔伯特时频分布，将各希尔伯特时频分布叠加得到原信号的时频分布包括：
[0041]
对各单分量信号进行希尔伯特变换，将实信号变为复信号；
[0042]
基于复信号计算各单分量信号的幅值包络和瞬时频率；
[0043]
基于得到的各单分量信号的幅值包络和瞬时频率，构造各单分量信号的希尔伯特时频分布；
[0044]
将各单分量信号的希尔伯特时频分布叠加，得到原信号的时频分布。
[0045]
其中，希尔伯特变换的公式为：
[0046][0047]
其中，yk(t)为通过时变滤波提取的第k个特征阶次的单分量信号，h[yk(t)]表示yk(t)的希尔伯特变换复信号；
[0048]
单分量信号的希尔伯特时频分布tfrk(t,f)表示为：
[0049]
tfrk(t,f)＝ak(t)δ[f(t)-fk(t)]
[0050][0051][0052]
其中，yk(t)为通过时变滤波提取的第k个特征阶次的单分量信号，tfrk(t,f)表示单分量信号yk(t)的希尔伯特时频分布，ak(t)为单分量信号yk(t)的幅值包络，fk(t)为单分量信号yk(t)的瞬时频率，δ[
·
]为狄拉克函数，f(t)表示t时刻的频率自变量，t和f分别为时间和频率自变量；
[0053]
原信号的时频分布tfr(t,f)表示为：
[0054][0055]
本发明实施例提供的技术方案带来的有益效果至少包括：
[0056]
(1)本实施例采用的时变滤波可以直接在时域提取非平稳工况下时域信号中的时变单分量信号，从而保留时域信号的瞬时相位、频率、幅值等关键信息，同时避免了自适应分解方法的模式混淆问题，也避免了角域重采样方式引入的计算误差，提高了后续时频分布的准确性；
[0057]
(2)与基于积分变换的传统时频分析方法相比，本实施例中基于单分量信号和希尔伯特变换，能够得到具有较高时频分辨率的时频分布，以揭示信号瞬态特征，有利于被监
测设备的瞬态特征提取，不受海森堡不确定原理的限制，且不受交叉项干扰，具有较高的时间和频率分辨率，能够清晰地展示各密集分布成分的瞬态特征。
附图说明
[0058]
为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0059]
图1为本发明实施例提供的基于时变滤波的旋转机械信号时频分析方法的流程示意图；
[0060]
图2为本发明实施例2提供的基于时变滤波的旋转机械信号时频分析方法的流程示意图；
[0061]
图3为本发明实施例3提供的基于短时傅里叶变换的时频分布示意图；
[0062]
图4为本发明实施例3提供的基于经验模式分解的时频分布示意图；
[0063]
图5为本发明实施例3提供的基于时变滤波的时频分布示意图。
具体实施方式
[0064]
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。
[0065]
实施例一
[0066]
如图1所示，本发明实施例提供了一种基于时变滤波的旋转机械信号时频分析方法，包括：
[0067]
s101，获取非平稳工况下目标旋转机械的时域信号作为原信号；
[0068]
本实施例中，获取目标旋转机械的时域信号x(t)，时域信号x(t)长度为n；所述时域信号包括：振动、位移、声音、应变、压力或电信号。
[0069]
s102，确定原信号中各单分量信号的时频曲线；
[0070]
本实施例中，单分量信号是指只有一个瞬时频率成分的信号。
[0071]
本实施例中，提供了3种方式来确定原信号中各单分量信号的时频曲线；其中，方式一包括以下步骤：
[0072]
获取目标旋转机械的输入转频s(t)；
[0073]
计算旋转机械各特征频率(具体为：各齿轮的特征频率)相对输入转频的特征阶次，得到特征阶次序列pk，其中，k＝1,2,
…
k，k表示理论存在的特征频率数量；
[0074]
将输入转频s(t)与特征阶次pk相乘，得到原信号中各单分量信号的时频曲线sk(t)；其中，sk(t)表示k个特征阶次的时频曲线。
[0075]
方式二包括以下步骤：
[0076]
对原信号进行短时傅里叶变换，得到原信号的传统时频分布图；
[0077]
基于得到的传统时频分布图，估计各时频脊线峰值对应的时频曲线，得到原信号中各单分量信号的时频曲线。
[0078]
方式三包括以下步骤：
[0079]
获取目标旋转机械的输入转频；
[0080]
对原信号进行短时傅里叶变换，得到原信号的传统时频分布图；
[0081]
对于与输入转频成正比的频率成分，计算特征阶次，将输入转频与特征阶次相乘估算其时频曲线；
[0082]
对于其它与输入转频无关的频率成分，基于得到的传统时频分布图，估计各时频脊线峰值对应的时频曲线。
[0083]
s103，通过时变滤波在时域提取各时频曲线对应的单分量信号；具体可以包括以下步骤：
[0084]
a1，根据得到的时频曲线，确定时变滤波的上截止频率f
cu
(t)和下截止频率f
cl
(t)：
[0085][0086]
其中，t表示时间，sk(t)表示第k个特征阶次的时频曲线，b(t)表示时变滤波带宽；
[0087]
a2，根据得到的时变滤波的上截止频率f
cu
(t)和下截止频率f
cl
(t)，确定局部滤波器长度ln：
[0088]
其中，ρ
max
≠0
[0089][0090]
其中，n＝0,1,
…
，n-1，n为时域信号的长度，l表示滤波器长度参数，α表示变化系数，α取值范围(0,1)，ρ(n)表示局部截止频率的导数平方和，ρ
max
表示ρ(n)的最大值，f
cu
(n)表示时变滤波的上截止频率，f
cl
(n)表示时变滤波的下截止频率，t(n)表示第n个时间点，表示向下取整，此公式所求ln为局部滤波器长度，其取值范围被限制在[l－αl,l]内。
[0091]
本实施例中，n和t的本质都表示时间，当时间可表示为连续值时用t表示，无法表示为连续值时，则用n表示，n＝0,1,
…
，n-1。
[0092]
a3，根据得到的上截止频率f
cu
(t)、下截止频率f
cl
(t)和局部滤波器长度ln，构造滤波器矩阵w；
[0093]
本实施例中，根据得到的上截止频率f
cu
(t)、下截止频率f
cl
(t)和局部滤波器长度ln，使用窗函数法设计局部带通滤波器w
(n)
，并生成滤波器矩阵w：
[0094][0095]
其中，w
(n)
表示t(n)处的局部滤波器，且其中，为局部滤波器w
(n)
的滤波器系数，局部滤波器长度等于2ln+1。
[0096]
a4，使用滤波器矩阵对原信号进行时变滤波得到单分量信号：
[0097]
y＝wx
[0098]
其中，x为原信号，y为单分量信号；
[0099]
其中，为输入信号x(t)，x与x(t)同为输入信号的不同表达形式；为单次时变滤波输出信号，即所提取的第k个特征阶次的单分量信号yk(t)，(x、y分别为x(t)和yk(t)的n行1列的矩阵格式)；
[0100]
a5，重复步骤a1-a4，得到各时频曲线对应的单分量信号。
[0101]
s104，计算提取的各单分量信号的希尔伯特时频分布，将各希尔伯特时频分布叠加得到原信号的时频分布，具体可以包括以下步骤：
[0102]
b1，对各单分量信号进行希尔伯特变换，将实信号变为复信号；其中，希尔伯特变换的公式为：
[0103][0104]
其中，yk(t)为通过时变滤波提取的第k个特征阶次的单分量信号，h[yk(t)]表示yk(t)的希尔伯特变换复信号；
[0105]
b2，基于复信号计算各单分量信号的幅值包络ak(t)和瞬时频率fk(t)；其中，ak(t)和fk(t)分别表示为：
[0106][0107][0108]
其中，yk(t)为通过时变滤波提取的第k个特征阶次的单分量信号，ak(t)表示此单分量信号的幅值包络，fk(t)表示此单分量信号的瞬时频率，j表示虚数；
[0109]
b3，基于得到的各单分量信号的幅值包络和瞬时频率，构造各单分量信号的希尔伯特时频分布tfrk(t,f)：
[0110]
tfrk(t,f)＝ak(t)δ[f(t)-fk(t)]
[0111]
其中，δ[
·
]为狄拉克函数，tfrk(t,f)表示单分量信号yk(t)的希尔伯特时频分布，ak(t)表示此单分量信号yk(t)的幅值包络，fk(t)表示此单分量信号yk(t)的瞬时频率，f(t)表示t时刻的频率自变量，t和f分别为时间和频率自变量；
[0112]
b4，将各单分量信号的希尔伯特时频分布叠加，得到原信号的时频分布tfr(t,f)：
[0113][0114]
本实施例中，通过希尔伯特变化计算各单分量信号的瞬时频率和幅值，得到各单分量信号的希尔伯特时频分布，将各个单分量信号的时频分布矩阵叠加，得到原信号的时频分布。
[0115]
本发明实施例所述的基于时变滤波的旋转机械信号时频分析方法，对于非平稳工
况下旋转机械的时域信号，提取各单分量信号并构造其时频分布，克服时频模糊、交叉项干扰、模式混淆问题，得到具有较高时频分辨率的时频分布，以揭示信号瞬态特征。
[0116]
综上，本发明实施例提供的基于时变滤波的旋转机械信号时频分析方法至少具有以下有意效果：
[0117]
(1)本实施例采用的时变滤波可以直接在时域提取非平稳工况下时域信号中的时变单分量信号，从而保留时域信号的瞬时相位、频率、幅值等关键信息，同时避免了自适应分解方法的模式混淆问题，也避免了角域重采样方式引入的计算误差，提高了后续时频分布的准确性；
[0118]
(2)与基于积分变换的传统时频分析方法相比，本实施例中基于单分量信号和希尔伯特变换，能够得到具有较高时频分辨率的时频分布，以揭示信号瞬态特征，有利于被监测设备的瞬态特征提取，不受海森堡不确定原理的限制，且不受交叉项干扰，具有较高的时间和频率分辨率，能够清晰地展示各密集分布成分的瞬态特征。
[0119]
实施例2
[0120]
本实施例中，如图2所示，获取的时域信号x(t)为四个非线性调幅调频信号的叠加，其解析式为：
[0121][0122]
其中，x(k)表示第k个单分量信号，其调幅成分为0.5sin(2π*0.5t)，调频成分为2πk∫(-100t2+250t+50)dt，采样频率为2000hz，时长为1s，信号长度n为2000。此仿真信号可模拟时变多分量的各类信号类型，包括但不限于振动、声音、压力、电信号等物理量。
[0123]
对时域信号x(t)进行短时傅里叶变换，并基于时频脊线估计各单分量信号的时频曲线sk(t)；
[0124]
以得到的时频曲线sk(t)为中心频率，并设置时变滤波带宽b(t)，计算时变滤波的上截止频率f
cu
(t)和下截止频率f
cl
(t)；
[0125]
基于得到的上截止频率f
cu
(t)和下截止频率f
cl
(t)生成滤波器矩阵w，并执行时变滤波y＝wx，得到各单分量信号yk(t)；
[0126]
利用希尔伯特变换计算各单分量信号的瞬时频率fk(t)和幅值包络ak(t)，并生成各单分量信号的希尔伯特时频分布tfrk(t,f)；
[0127]
将各单分量信号的希尔伯特时频分布进行叠加，得到原信号的时频分布tfr(t,f)；
[0128]
图2中各子图的含义依次为：输入多分量信号波形，基于传统时频分布的脊线估计，各分量的时变截止频率，通过时变滤波提取的各个单分量信号，各个单分量信号的时频分布，重构的原始信号的时频分布。
[0129]
实施例3
[0130]
本实施例中，获取的时域信号x(t)为某水轮机转子的径向位移信号；
[0131]
计算水轮机转子的特征阶次pk，k＝1,2,
…
,5；
[0132]
将实施例3中水轮机转子的输入转频s(t)与特征阶次相乘，得到5个估算的时频曲线sk(t)，k＝1,2,
…
,5；
[0133]
计算时域信号x(t)的短时傅里叶变换，并基于时频脊线估计除上述5个时频曲线
外的其它单分量信号的时频曲线sk(t)，k＝6,
…
,k；
[0134]
以得到的时频曲线sk(t)为中心频率，并设置时变滤波带宽b(t)，计算时变滤波的上截止频率f
cu
(t)和下截止频率f
cl
(t)；
[0135]
基于得到的上截止频率f
cu
(t)和下截止频率f
cl
(t)生成滤波器矩阵w，并执行时变滤波y＝wx，得到各单分量信号yk(t)；
[0136]
利用希尔伯特变换计算各单分量信号的瞬时频率fk(t)和幅值包络ak(t)，并生成各单分量信号的希尔伯特时频分布tfrk(t,f)；
[0137]
将各单分量信号的希尔伯特时频分布进行叠加，得到原信号的时频分布tfr(t,f)。
[0138]
本实施例中，根据图3-图5可知，本发明实施例提供的基于时变滤波的旋转机械信号时频分析方法，能够得到具有较高时频分辨率的时频分布。
[0139]
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。