多源时频脊线提取方法与流程

本发明涉及机械设备信号处理领域，具体涉及一种基于快速路径最优搜索和动态基角度的多源时频脊线提取方法。

背景技术：

轴承作为旋转机械的关键零部件之一，其状态直接影响到旋转机械的运行状态，一旦轴承出现故障，将很可能带来巨大的经济损失，甚至导致人员伤亡。信号处理与分析是实现轴承故障诊断的有效途径。在机械故障特征提取方面，较早的研究是基于时域分析和频域分析的信号处理方法，大多适用于分析线性、平稳的动态信号。这些方法简单高效，在实际机械故障诊断中得到广泛应用。

传统技术存在以下技术问题：

但是，在旋转机械设备运行工况变化或关键部件出现局部故障时，采集的振动信号往往呈现非线性、非平稳特性，传统的时域、频域分析方法往往存在很大的局限性。

技术实现要素：

本发明要解决的技术问题是提供一种基于快速路径最优搜索和动态基角度的多源时频脊线提取方法，在用传统的短时傅里叶对信号进行分析的基础上，引入动态基角度策略，优化时频脊线的聚集性，从而可以提取得到更加准确的结果。此外，该发明利用快速路径最优搜索算法考虑整个始于中所有幅值极大值，建立前后时刻点所有幅值极大值之间所有可能的优化路径，保证了瞬时时频脊线的连续性和准确性。然后，在提取的顺势极限的基础上计算得出频率与目标相匹配的基函数，从而提高时频聚集性。

为了解决上述技术问题，本发明提供了一种基于快速路径最优搜索和动态基角度的多源时频脊线提取方法，包括：

步骤1：短时傅里叶变换及频带的选择。对振动信号中采用短时傅里叶分析，将轴承信号分为低频段和共振频带。

步骤2：计算出每个时刻点处所有幅值的局部极大值及其对应的频率。根据快速路径最优算法确定每个极值点与前后极值点之间可能的优化路径。根据优化方程的解确定最后一个时刻点处应被提取的频率值对应的极值，根据该极值点选取与其对应的那条快速最优路径，这条快速最优路径即为最终提取得到的瞬时频率脊线；

步骤3：应用一种动态基角度策略，结合快速路径最优搜索结果对原始时频图做后处理，即利用快速路径最优算法提取得到的瞬时频率脊线求取脊线各点的切角，得出频率与瞬时频率匹配的基函数；最后，由于低频段和共振频带中任意两条修正后的时频脊线之间的比值表征故障类型，据此进行轴承故障诊断。

在其中一个实施例中，步骤2具体包括：

步骤2.1：计算时频图上每个时刻点处幅值的局部极大值，得到该时刻点处所有幅值极大值对应的频率值

fm(τn)＝f,(d[x(τn,f)]/df＝0||d²[x(τn,f)]/df²＜0),m＝1,2,...,np(9)

其中，τn(n＝1,2,...,n)表示时频图上的时刻点，np(τn)表示时刻点τn处局部极大值的个数，fm(τn)是时刻点τn处第m个极值对应的频率值，其幅值记为spm(τn)。

步骤2.2：为准确提取时频脊线，提出一种路径优化算法确定每个时刻点处应被提取的幅值极大值，该极值对应的、fm(τn)就是最终应被提取的频率值，该路径优化算法表示如下：

其中，mc(τn)决定了时刻点τn处被提取的幅值极大值，f[]表示优化方程，根据该方程的解即可提取出最优的时频脊线。优化方程具体表示如下：

m[]＝perc0.5[],iqr[]＝perc0.75[]-perc0.25[](14)

其中，fd(τn-1)表示时刻点τn-1处的候选脊频率值，δfd是fd的导数，percp[f(t)]表示f(t)的第p分位数。

步骤2.3：优化方程的解可以表示为

其中，n＝1,2,...,n，m＝1,2,...,np(τn)，q(m,τn)用以确定时刻τn(n＝2,...,n)处每个幅值极大值应与前一时刻点处哪一个极值点构成优化路径，规划所有可能路径。u(m,τn)表示得到优化结果的中间向量，即选择时刻τn处u(m,τn)取最大值的幅值极大值点对应的频率值作为最终应被提取的瞬时频率，相应的由q(m,τn)确定的该极值点下的路径作为快速路径最优，至此提取出瞬时频率脊线。

在其中一个实施例中，所述步骤3包括：

步骤3.1：应用一种动态基角度时频图增强策略，考虑线性变换基的角度而非窗长，利用由快速路径最优算法从原始时频图中提取得到的瞬时频率脊线计算出切角。得出频率与目标信号频率相匹配的基函数，而后对时频图做后处理，以提高时频图聚集性。

步骤3.2：低频段和共振频带中任意两条修正后的时频脊线之间的比值表征着故障类型，据此可进行轴承故障诊断。从低频段中提取并修正得到两条时频脊线fr1和fr2，从共振频带中提取并修正得到两条时频脊线fr3和fr4，其两两之间的比值有θ1＝fr3/fr1，θ2＝fr3/fr2，θ3＝fr4/fr1，θ4＝fr4/fr2，它们应为故障特征系数或其倍数。诊断不同故障类型的轴承，计算其故障特征系数的均值时，应对这些比值θ1，θ2，θ3和θ4做同步化处理，再求取其均值，有θ＝(η1θ1+η2θ2+η3θ3+η4θ4)/4；

其中，η1,η2,η3,η4表示相应比值的同步化系数。至此，轴承故障类型可根据故障特征系数均值诊断得出。

在其中一个实施例中，步骤1中，低频段为[0,400]hz。

一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现任一项所述方法的步骤。

一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现任一项所述方法的步骤。

一种处理器，所述处理器用于运行程序，其中，所述程序运行时执行任一项所述的方法。

本发明的有益效果：

1、该发明首先应用快速路径最优搜索策略，得到连续和准确的瞬时频率脊线，然后在得到的脊线的基础上，计算对应时刻点的切角，得到频率与目标信号频率相匹配的基函数，适用于处理变转速非平稳状态下的滚动轴承振动信号。

2、该发明所提出的快速路径最优搜索和动态基角度结合的多源时频脊线提取方法与传统的方法相比，都能得到更加精确的结果。相比传统的区域峰值搜索，快速路径搜索可以减少信号分析的计算量并提升分析结果的准确性；相比传统的短时傅里叶变换，由于短时傅里叶变换采用固定的窗长分段对信号进行时频分析，结合与实际频率变化相对应的基函数，可以进一步增强时频聚集性。

附图说明

图1为采集实验信号的轴承试验台。

图2为本发明所提出的基于快速路径最优搜索和动态基角度的多源时频脊线提取方法的流程图。

图3为快速路径最优搜索算法的示意图，(a)计算局部极大值；(b)规划路径；(c)选取最优路径。

图4为基于本发明的故障轴承模拟信号快速路径最优时频脊线提取结果，其中，图(a)为低频段提取结果，图(b)为共振频带提取结果。

图5为本发明提出的基于快速路径最优搜索得到脊线后计算得出的故障轴承模拟信号脊线切角。

图6为本发明模拟信号低频段时频图比较结果，其中，(a)窗长为1200时的时频图；(b)窗长为1200时的脊线提取结果；(c)窗长为2161时的时频图；(d)窗长为2161时的脊线提取结果。

图7为本发明模拟信号共振频带时频图比较结果，其中，(a)窗长为5200时的时频图；(b)窗长为5200时的脊线提取结果；(c)窗长为9295时的时频图；(d)窗长为9295时的脊线提取结果。

图8为基于本发明的轴承故障实验信号快速路径最优时频脊线提取结果，其中，图(a)为低频段提取结果，图(b)为共振频带提取结果。

图9为本发明提出的基于快速路径最优搜索得到脊线后计算得出的故障轴承实际信号脊线切角。

图10为本发明模拟信号低频段时频图比较结果，其中，(a)窗长为4590时的时频图；(b)窗长为4590时的脊线提取结果；(c)窗长为5818时的时频图；(d)窗长为5818时的脊线提取结果。

图11本发明模拟信号共振频段时频图比较结果，其中，(a)窗长为2500时的时频图；(b)窗长为2500时的脊线提取结果；(c)窗长为3377时的时频图；(d)窗长为3377时的脊线提取结果。

具体实施方式

下面结合附图和具体实施例对本发明作进一步说明，以使本领域的技术人员可以更好地理解本发明并能予以实施，但所举实施例不作为对本发明的限定。

仿真信号模型的构造如下所示

x(t)＝x1(t)+x2(t)+η(t)(17)

其中，x1(t)表示故障特征相关成分，定义为

其中，[1+αcos(2πx2(t))]反应出由转速信号x2(t)解调的共振现象，α＝0.11为调制幅值，m表示信号长度，am表示第m个脉冲响应的幅值，β＝800表示与阻尼有关的系数，ωr表示被激励的共振频率，u(t)是一个单位阶跃函数，tm表示第m个脉冲发生的时间，定义为

其中，μm表示滑移率，在0.01～0.02之间变化，x2(t)表示旋转轴转速相关的成分，η(t)表示高斯白噪声(信噪比为-12db)。x2(t)被定义为

其中，ai＝1,2,3表示幅值，相应地取值为1.6，1.2和1，模拟信号的瞬时转轴频率被定义为故障特征频率设为3.7。采样频率为20khz，采样时间为5.1s。

实验数据为轴承外圈故障。轴承试验台设备如图1所示。

驱动电机(sj200-022nfu)通过联轴器链接主轴，主轴直径为1英寸。两个轴承(er16k)安装在主轴上的轴承座内，左侧滚动轴承存在外圈故障，一个质量为5.03kg的质量盘安装在两轴承之间以提供负载。驱动轴的右端链接到较小的滑轮上，齿轮箱轴链接较大滑轮，两个滑轮上的皮带链接驱动轴以驱动齿轮箱轴，两滑轮的直径比约为1:2.6。转速计和加速度传感器安装在试验台上，分别用以测量轴转速和振动信号。为确保变速箱的振动具有振幅优势，加速度计安装在变速箱附近。然后，通过数据采集卡收集振动信号并使用labview进行采样。轴承具体参数如表1所示：

表1外圈故障轴承参数

如图2所示，本发明包括以下步骤：

步骤1：短时傅里叶变换及频带的选择。对振动信号中采用短时傅里叶分析，将轴承信号分为低频段和共振频带。其中，低频段约为[0,400]hz，足以覆盖常见的轴转速及其倍频信息。

步骤2：计算出每个时刻点处所有幅值的局部极大值及其对应的频率。根据快速路径最优算法确定每个极值点与前后极值点之间可能的优化路径。根据优化方程的解确定最后一个时刻点处应被提取的频率值对应的极值，根据该极值点选取与其对应的那条快速最优路径，这条快速最优路径即为最终提取得到的瞬时频率脊线。

步骤3：为解决利用短时傅里叶变换得到的时频图聚集性受线性变换基函数频率影响较大的问题，应用一种动态基角度策略，结合快速路径最优搜索结果对原始时频图做后处理，即利用快速路径最优算法提取得到的瞬时频率脊线求取脊线各点的切角(即与水平线所夹的锐角，顺时针方向为正值，逆时针方程为负值)，得出频率与瞬时频率匹配的基函数。最后，由于低频段和共振频带中任意两条修正后的时频脊线之间的比值表征故障类型，据此进行轴承故障诊断。

进一步的，所述步骤具体包括以下步骤：

步骤1：对振动信号中采用短时傅里叶分析，将轴承信号分为低频段和共振频带。

所述步骤2包括：

步骤2.1：计算时频图上每个时刻点处幅值的局部极大值，得到该时刻点处所有幅值极大值对应的频率值

fm(τn)＝f,(d[x(τn,f)]/df＝0||d²[x(τn,f)]/df²＜0),m＝1,2,...,np(21)

其中，τn(n＝1,2,...,n)表示时频图上的时刻点，np(τn)表示时刻点τn处局部极大值的个数，fm(τn)是时刻点τn处第m个极值对应的频率值，其幅值记为spm(τn)。

步骤2.2：为准确提取时频脊线，提出一种路径优化算法确定每个时刻点处应被提取的幅值极大值，该极值对应的、fm(τn)就是最终应被提取的频率值，该路径优化算法表示如下

其中，mc(τn)决定了时刻点τn处被提取的幅值极大值，f[]表示优化方程，根据该方程的解即可提取出最优的时频脊线。优化方程具体表示如下

m[]＝perc0.5[],iqr[]＝perc0.75[]-perc0.25[](26)

其中，fd(τn-1)表示时刻点τn-1处的候选脊频率值，δfd是fd的导数，percp[f(t)]表示f(t)的第p分位数。

步骤2.3：优化方程的解可以表示为

其中，n＝1,2,...,n，m＝1,2,...,np(τn)，q(m,τn)用以确定时刻τn(n＝2,...,n)处每个幅值极大值应与前一时刻点处哪一个极值点构成优化路径，规划所有可能路径。u(m,τn)表示得到优化结果的中间向量，即选择时刻τn处u(m,τn)取最大值的幅值极大值点对应的频率值作为最终应被提取的瞬时频率，相应的由q(m,τn)确定的该极值点下的路径作为快速路径最优，至此提取出瞬时频率脊线。

总的来说，步骤2通过快速路径最优搜索算法考虑了整个时域中所有幅值极大值，建立了前后时刻点所有幅值极大值之间所有可能的优化路径，保证了瞬时时频脊线的连续性和准确性。

所述步骤3包括：

步骤3.1：应用一种动态基角度时频图增强策略，考虑线性变换基的角度而非窗长，利用由快速路径最优算法从原始时频图中提取得到的瞬时频率脊线计算出切角。得出频率与目标信号频率相匹配的基函数，而后对时频图做后处理，以提高时频图聚集性。

步骤3.2：低频段和共振频带中任意两条修正后的时频脊线之间的比值表征着故障类型，据此可进行轴承故障诊断。从低频段中提取并修正得到两条时频脊线fr1和fr2，从共振频带中提取并修正得到两条时频脊线fr3和fr4，其两两之间的比值有θ1＝fr3/fr1，θ2＝fr3/fr2，θ3＝fr4/fr1，θ4＝fr4/fr2，它们应为故障特征系数或其倍数。诊断不同故障类型的轴承，计算其故障特征系数的均值时，应对这些比值θ1，θ2，θ3和θ4做同步化处理，再求取其均值，有θ＝(η1θ1+η2θ2+η3θ3+η4θ4)/4。

其中，η1,η2,η3,η4表示相应比值的同步化系数。至此，轴承故障类型可根据故障特征系数均值诊断得出。

该策略由于基函数频率与目标脊线频率匹配，在一定程度上也降低了时频图可读性对窗长的敏感性，有效解决了传统短时傅里叶变换下基函数频率与目标信号频率不匹配和最优窗长难以选择的难题。

快速路径最优搜索算法的示意图如图3所示。

经过分析，轴承的故障仿真信号结果如图4所示，图4(a)表示低频段提取结果，图4(b)表示共振频带提取结果。然后引入平均相对误差，两种脊线提取算法的瞬时频率误差如表2所示。

表2模拟信号局部峰值和快速路径最优搜索结果比较

与局部峰值搜索算法提取结果比较，可明显看到瞬时频率脊线整体的波动性和脊线后半部分的突变有所改善。分析表格数据可知，无论对于低频段或共振频带，快速路径最优搜索算法相较局部峰值搜索算法准确度更高，最高可使瞬时频率脊线准确度提高74％。

为验证动态基角度策略的有效性，计算利用快速路径最优搜索算法提取得到的瞬时频率脊线在每个时刻点的切角，如图5所示，并据此选取得到模拟信号低频段和共振频带时频图。为验证该策略下窗长对时频图聚集性的影响，本章将变换窗长并利用快速路径最优搜索算法从不同窗长下的时频图中提取瞬时频率脊线，通过分析脊线的平均相对误差来验证。

对于低频段，图6(a)和6(c)分别表示窗长为1200和2161时利用动态基角度策略得到的时频图，图6(c)和6(d)分别表示利用快速路径最优搜索算法从中提取得到的瞬时频率脊线，与图4(a)所示的从短时傅里叶变换得到的时频图中利用快速路径最优搜索算法提取得到的瞬时频率脊线相比，可以看到脊线后半部分已无突变。利用快速路径最优搜索算法提取得到的瞬时频率脊线的平均相对误差都有所减少，最高可减少52％的误差，验证了该脊线提取算法的有效性。

然后，利用基于故障特征系数均值的故障诊断策略，利用快速路径最优搜索算法从动态基角度策略下的时频图中提取得到的瞬时频率脊线之间的两两比值如表3所示。求其均值为3.684，与设定值3.7之间的相对误差为0.4％，诊断结果进一步验证了快速路径最优搜索算法和动态基角度策略的有效性。

表3模拟信号低频段和共振频带中两两瞬时频率脊线之比

利用快速路径最优搜索算法从低频段和共振频带中提取多源瞬时频率脊线，提取结果分别如图8(a)和8(b)所示。与局部峰值搜索算法提取结果相比曲线整体波动性减少，尤其低频段中第一条脊线前端的突变情况得到较大改善。引入平均相对误差，表4表示两种脊线搜索算法下瞬时频率脊线的平均相对误差。分析表中数据可知，利用快速路径最优搜索算法提取得到的瞬时频率脊线的平均相对误差都有所减少，最高可减少52％的误差，验证了该脊线提取算法的有效性。

表4轴承外圈故障实验信号局部峰值和快速路径最优搜索结果比较

图9表示利用快速路径最优搜索算法提取得到的瞬时频率脊线在每个时刻点处的切角，据此得到外圈故障轴承信号低频段和共振频带时频图。窗长为4590和5815时利用动态基角度策略得到低频段时频图分别如图10(a)和10(c)所示，利用快速路径最优搜索算法从中提取得到的瞬时频率脊线分别如图10(c)和10(d)所示，与局部峰值搜索算法提取结果相比，脊线前端突变得到改善。脊线提取准确性提高了35％以上，验证了动态基角度策略相较于传统短时傅里叶变换能有效提高时频图聚集性和瞬时频率脊线提取准确性。

求取动态基角度策略下由快速路径最优搜索算法提取得到的瞬时频率脊线两两比值，结果如表5所示。利用上述提出的基于故障特征系数均值的轴承故障诊断策略，得到故障特征系数为3.533，与真实值相对误差为0.9％。进一步分析表格可知，每个系数更接近真实值，表明脊线提取更加准确，有效验证了快速路径最优搜索算法和动态基角度策略的可行性。

表5轴承外圈故障信号低频段和共振频带中两两瞬时频率脊线之比

针对局部峰值搜索算法和“一步法”对初始时刻点的依赖和由经验选取局部峰值搜索算法搜索范围的局限性问题，本发明应用快速路径最优搜索算法，根据快速路径最优算法规划前后时刻点所有幅值极大值点之间可能的快速路径最优，再从中选取最优解作为最终提取脊线，保证了瞬时频率脊线提取的准确性和连续性。针对变转速下短时傅里叶变换基函数的频率难以和时变的瞬时频率相匹配和窗长对时频图聚集性影响较大的问题，本发明应用动态基角度策略，利用快速路径最优算法提取得到的瞬时频率脊线，求取脊线各点的切角，得出频率与目标信号频率相匹配的基函数，再进行时频变换，进而得到整个时域上的时频图。在此基础上构建模拟信号进行分析验证，结果表明快速路径最优搜索算法和动态基角度策略能有效提高瞬时频率脊线提取的准确性。利用基于故障特征系数均值的诊断策略，基于提取得到的多源瞬时频率脊线对轴承进行故障诊断，诊断结果进一步验证了快速路径最优搜索算法和动态基角度策略的有效性。考虑线性变换基的角度，利用快速路径最优搜索结果的切角指导获取优化后的模拟信号低频段和共振频带时频图，再次利用快速路径最优搜索算法从中提取出的瞬时频率脊线误差可降低79％。快速路径最优搜索算法和动态基角度策略在提高瞬时频率脊线提取和轴承故障诊断准确度方面具有一定的有效性。

以上所述实施例仅是为充分说明本发明而所举的较佳的实施例，本发明的保护范围不限于此。本技术领域的技术人员在本发明基础上所作的等同替代或变换，均在本发明的保护范围之内。本发明的保护范围以权利要求书为准。