CN104568435A - 一种基于时频能量谱提取风机状态信号调制成分的方法 - Google Patents

Abstract

本发明公开了一种基于时频能量谱提取风机状态信号调制成分的方法，用于获得目标振动信号中指定频率的调制成分，包括以下步骤：A.绘制目标振动信号的希尔伯特-黄谱,并采用瞬时无量纲频率归一化方法将希尔伯特-黄谱转换为无量纲频率信号能量谱；B.从所述无量纲频率信号能量谱中提取频率在目标振动信号中指定频率附近的多组振动信号的无量纲能量，并采用频谱分析法获得所述多组无量纲能量对应的调制信号频率，然后将这些调制信号频率的均值作为所述调制成分的调制频率。发明的方法可以估计为了进行信号分析所需要的状态信号采集的最短时间长度，能够保证在不丢失状态信号特征的同时获得较高的检测精度及较少的数据长度。

Description

一种基于时频能量谱提取风机状态信号调制成分的方法

技术领域

本发明涉及风电技术领域，特别是涉及一种基于时频能量谱提取风机状态信号调制成分的方法。

背景技术

随着风电场规模和数量的增加，风机运行及维护问题越来越重要，风机可靠性及故障预警问题受到越来越多的重视。风机的机械构成比较复杂，加上风能的随机性特点，风速的变化及负荷的波动很容易导致风机机械部分的损伤或故障，直接增加了风机运行及维护成本，所以，实施风机机械部件的状态监测和故障诊断对提高风机的经济效益和可靠性有重要意义。

振动监测是风机机械安全评估的重要手段，已经取得了很多有效的分析方法和理论成果。传统的振动监测系统能够在时空间内对振动的幅度和频率进行有限的监测，随着数字信号处理及计算机技术的发展，数据挖掘及时-频分析方法取得了长足的进步，成为基于振动信号分析的齿轮箱故障诊断中应用最广泛的一种，而瞬时频率的提取成为整个变转速振动故障诊断的关键环节。

小波方法是具有多尺度分辨率的信号处理方法，先进技术[1](参见LinC.H.,Wang C.H..Adapt ive Wavelet Networks For Power Qual ityDetection and Discrimination in a Power System[J].IEEE Trans.PowerDelivery,2006,21:1106–1113)利用自适应小波分析方法，根据转速的变化在线调整小波参数，可以消除时变转速对振动信号的影响，但需要具有小波基函数的选择经验。先进技术[2](参见Wu T.Y.,Chen J.C.,Wang C.C..Characterization of Gear Faults in Variable Rotating SpeedUsing Hilbert-Huang Transform and Instantaneous DimensionlessFrequency Normalization[J].Mechanical Systems and SignalProcessing,2012,30:103–122)采用希尔伯特-黄变换方法分离故障信号，基于这种方法T.Y.Wu采用频率归一化的无量纲希尔伯特-黄谱方法，消除转速变化在希尔伯特-黄谱中对信号能量分布的影响，从而分离出真实的故障信号。与傅里叶变换相比，边际谱能比较准确地反映信号的实际频率成分，由信号时频能量衍生的边际谱可以处理非平稳信号。

因为实际检测中都是采集一个时间区间的状态数据进行计算，如何选择适当时长的数据，才能够保证用较少的数据包含尽可能多的频率特征并获得较高的检测精度是风机状态实时有效性的基础。传统方法是通过建模得到信号变化周期的数值估计，并将其作为采样区间长度。但因为风速变化会引起传动部件的转速的变化，导致状态信号的频率特征在时域上是时变的，并且包含大量的谐波或边带信号，给基于信号时频特征分析的风机状态监测带来很多困难，所以振动信号周期特性的估计必须考虑风机转速变化对故障信号提取的影响，而传统频域分析方法在处理此类问题的时候效果并不理想。

发明内容

本发明的目的是提供一种基于时频能量谱提取风机状态信号调制成分的方法，可以避免风机转速变化对信号提取的影响，能够估计为了进行信号分析所需要的振动信号采集的最短时间长度，保证在不丢失振动特征的同时获得较高的检测精度及较少的数据长度。

为实现上述目的，本发明采用如下技术方案：

一种基于时频能量谱提取风机状态信号调制成分的方法，用于获得目标振动信号中指定频率的调制成分，包括以下步骤：A.采集目标振动信号，绘制目标振动信号的希尔伯特-黄谱,并采用瞬时无量纲频率归一化方法将希尔伯特-黄谱转换为无量纲频率信号能量谱；B.从所述无量纲频率信号能量谱中提取频率在目标振动信号中指定频率附近的多组振动信号的无量纲能量，并采用频谱分析法获得所述多组无量纲能量对应的调制信号频率，然后将这些调制信号频率的均值作为所述目标振动信号中指定频率的调制成分的调制频率。

进一步地，所述步骤具体包括：

A.采集目标振动信号x(t)；

用EMD方法将x(t)分解为公式(1)的形式：

$x\left(t\right)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{y}_k\left(t\right)-r_m\left(t\right)---\left(1\right)$

其中y_k(t)是x(t)的第k个IMF分量，并且满足单分量信号条件；r_m(t)是经过m次分解后的残差；

利用公式(2)对所述IMF分量y_k(t)进行循环归一化，直到第i次迭代过程中，IMF分量的所有时刻的瞬时幅值都小于单位值，

$y_k^i\left(t\right)=y_k^{i-1}\left(t\right)/e_k^{i-1}\left(t\right)---\left(2\right)$

其中，i∈[1,n]为迭代次数；为的包络函数；

于是，所述IMF分量被分解为独立的载波函数F_k(t)和调制函数A_k(t)，并有：

$F_k\left(t\right)=y_k^n\left(t\right)---\left(3\right)$

$A_k\left(t\right)=e_k^1\left(t\right)\·e_k^2\left(t\right)...e_k^n\left(t\right)---\left(4\right)$

基于所述载波函数F_k(t)将IMF分量的瞬时频率定义为：

${\mathrm{\ω}}_k\left(t\right)=\frac{1}{2\mathrm{\π}}\frac{{\mathrm{dF}}_k\left(t\right)}{\mathrm{dt}}---\left(5\right)$

其中ω_k(t)为第k个IMF分量的瞬时频率；

基于所述A_k(t)和ω_k(t)，得到所述振动信号的希尔伯特-黄谱H(ω,t)：

$H(\mathrm{\ω},t)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{A}_k\left(t\right)\cos (\∫2{\mathrm{\π\ω}}_k\left(t\right)\mathrm{dt})---\left(6\right)$

选择风机上某测量点的转动频率ω_ref(t)为参考，用公式(7)计算所述无量纲频率信号能量谱H(ω,t)：

$H(\mathrm{\ω},t)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{A}_k\left(t\right)\cos \left(2\mathrm{\π}\frac{{\mathrm{\ω}}_k\left(t\right)}{{\mathrm{\ω}}_{\mathrm{ref}}\left(t\right)}\mathrm{dt}\right)---\left(7\right)$

B.根据公式(10)计算目标振动信号中指定频率的调制成分的调制频率f_F：

$f_F\left(t\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{f}_H({\mathrm{\ω}}_j,t)/N---\left(10\right)$

其中N为提取的无量纲能量的组数，f_H(ω_j,t)为在ω_j处的振动信号能量H(ω_j,t)的包络线的瞬时频率，H(ω_j,t)为无量纲频率信号能量。

进一步地，所述步骤A中选择风机传动链中高速轴或低速轴的转动频率ω_ref(t)为参考。

进一步地，所述目标振动信号为风机齿轮箱高速轴振动信号，所述步骤B中的指定频率为风机齿轮箱的啮合频率。

进一步地，所述风机齿轮箱的啮合频率f_m根据公式(9)计算获得：

f_m＝f_i·n_i＝f_i+1·n_i+1 (9)

其中f_i和f_i+1分别为以啮合方式连接的两个传动轴的转动频率，n_i和n_i+1分别为其相应的齿轮齿数。

由于采用上述技术方案，本发明至少具有以下优点：

通过本发明的基于时频能量谱提取风机状态信号调制成分的方法，可以在振动信号周期特性的估计过程中，消除风机转速变化对故障信号提取的影响；可以估计为了进行信号分析所需要的状态信号采集的最短时间长度，能够保证在不丢失状态信号特征的同时获得较高的检测精度及较少的数据长度。

附图说明

上述仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，以下结合附图与具体实施方式对本发明作进一步的详细说明。

图1A、1B分别是高速轴实测振动信号的振动加速度及转速示意图。

图2是部分实测振动信号的IMF成分示意图。

图3是部分F₅(t)分量示意图。

图4是部分A₅(t)分量示意图。

图5是实测振动信号的部分希尔伯特-黄能量谱示意图。

图6是实测振动信号的部分无量纲时频能量谱示意图。

图7是本发明的基于时频能量谱提取风机状态信号调制成分的方法流程图。

具体实施方式

本发明提供了一种基于时频能量谱提取风机状态信号调制成分的方法，用于获得目标振动信号中指定频率的调制成分，包括以下步骤：A.采集目标振动信号，绘制目标振动信号的希尔伯特-黄谱,并采用瞬时无量纲频率归一化方法将希尔伯特-黄谱转换为无量纲频率信号能量谱；B.从所述无量纲频率信号能量谱中提取频率在目标振动信号中指定频率附近的多组振动信号的无量纲能量，并采用频谱分析法获得所述多组无量纲能量对应的调制信号频率，然后将这些调制信号频率的均值作为所述目标振动信号中指定频率的调制成分的调制频率。

该方法是分析目标振动信号中指定频率的调制成分的方法。如果该方法用于诊断齿轮箱断齿故障，则“目标振动信号”为采集于齿轮箱的振动信号，“指定频率”可以是故障齿轮的固有啮合频率，“指定频率的调制成分”则指在该啮合频率上包含的振动幅值的变化。

如图7所示，本发明基于信号时频特性提出了一种利用希尔伯特-黄谱提取信号的周期性特征的计算方法，并根据周期性特征确定状态监测最短有效采样时间的方法。为分析非线性、非平稳信号，首先利用经验模态分解(EMD)将信号分解为若干个固有模态函数(intrinsic mode function，IMF)分量之和，然后对每个IMF分量进行希尔伯特变换(Hilberttransform，HT)，得到瞬时频率和瞬时幅值，从而得到信号的Hilbert谱。为了分析变转速条件下的振动信号特征，通过瞬时无量纲频率归一化(instantaneous DLF normalization)方法，将时频域的希尔伯特能量谱转换成无量纲频率信号能量谱的时间分布，以消除时变的转速对Hilbert谱的能量分布的影响。

以下步骤1～4介绍归一化能量谱的绘制过程：

步骤1：设原始信号为x(t)，用EMD方法将x(t)分解为下面的形式：

$x\left(t\right)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{y}_k\left(t\right)-r_m\left(t\right)---\left(1\right)$

其中y_k(t)是x(t)的第k个IMF分量，并且满足单分量信号条件，r_m(t)是经过m次分解后的残差。

步骤2：将上一步得到的IMF分量y_k(t)标准化，得到解耦的代表目标解析信号的瞬时幅值特征。对IMF分量y_k(t)和其包络函数e_k(t)，利用下面的公式对y_k(t)进行循环归一化，直到第i次迭代过程中，y_k(t)的所有时刻的瞬时幅值都小于单位值，并将此时的y_k(t)记为

$y_k^i\left(t\right)=y_k^{i-1}\left(t\right)/e_k^{i-1}\left(t\right)---\left(2\right)$

其中，i∈[1,n]代表迭代次数。于是，原IMF分量被分解为相互独立的载波函数F_k(t)和调制函数A_k(t)，并有：

$F_k\left(t\right)=y_k^n\left(t\right)---\left(3\right)$

$A_k\left(t\right)=e_k^1\left(t\right)\·e_k^2\left(t\right)...e_k^n\left(t\right)---\left(4\right)$

IMF的瞬时角频率ω(t)可基于F(t)分量定义为：

${\mathrm{\ω}}_k\left(t\right)=\frac{1}{2\mathrm{\π}}\frac{{\mathrm{dF}}_k\left(t\right)}{\mathrm{dt}}---\left(5\right)$

步骤3：基于上面的A_k(t)和角频率ω(t)的定义，则原始测量信号的能量谱为：

$H(\mathrm{\ω},t)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{A}_k\left(t\right)\cos (\∫2{\mathrm{\π\ω}}_k\left(t\right)\mathrm{dt})---\left(6\right)$

其中A_k(t)为第k个IMF信号的模分量，ω_k(t)为第k个IMF信号的瞬时频率。至此，原始测量信号的能量谱绘制完成。

步骤4：能量谱的归一化：振动信号的时频空间特性是对齿轮箱故障进行识别和归类的依据，所以对振动信号进行时频空间的分析对诊断齿轮箱故障类型十分重要。因为风能的波动性特点，迫使风机机械传动系统以时变的速度运转，转速的变化进一步影响振动信号能量的时频分布，增加了故障诊断的难度。选择传动链中能够反映风机状态的周期性信号点作为被测点，如高速轴、低速轴等，将被测点的转动频率ω_ref(t)作为参考，用公式(7)计算振动信号无量纲能量谱，则无量纲能量谱代表信号在时频域上相对于参考转动频率ω_ref(t)的能量分布，方便分析振动信号的固有频率和啮合频率附近的幅频特征。

$H(\mathrm{\ω},t)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{A}_k\left(t\right)\cos \left(2\mathrm{\π}\frac{{\mathrm{\ω}}_k\left(t\right)}{{\mathrm{\ω}}_{\mathrm{ref}}\left(t\right)}\mathrm{dt}\right)---\left(7\right)$

步骤5：齿轮箱啮合频率计算。因为生产精度不高、安装不当、润滑不充分等原因，齿轮箱振动信号包含丰富的频率特征，包括各个齿轮固有的频率特征、啮合频率特征以及故障频率特征，并且故障频率总是在固有频率及啮合频率附近出现，设以啮合方式连接的两个传动轴的转动频率分别为f_i和f_i+1，其相应的齿轮齿数分别为n_i和n_i+1，则两个传动轴的频率和齿数之间有关系：

$\frac{f_i}{f_{i+1}}=\frac{n_{i+1}}{n_i},---\left(8\right)$

啮合频率f_m为：

f_m＝f_i·n_i＝f_i+1·n_i+1， (9)

步骤6：振动信号调制频率计算。正常齿轮由于刚度的影响，故障振动信号经常是以某个频率为载波信号分量，以故障所在位置的固有振动为调制信号的调幅振动。其波形为周期性的衰减波形，且低频信号具有近似正弦波的啮合波形。若确定故障位置及故障类型，必须首先确定该位置调幅振动信号的调制频率。为减小计算误差，提取啮合频率附近的几组振动能量并分别计算其调制信号的频率，然后将这些调制信号频率的均值作为最终的调制信号的频率。假设提取了N组振动能量，令H(ω_j,t)代表振动信号在ω_j处的能量(j＝1,…,N)，C_H(ω_j,t)为H(ω_j,t)的包络线，f_H(ω_j,t)为该包络线的瞬时频率，则故障信号的调制频率f_F由式(10)确定：

$f_F\left(t\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{f}_H({\mathrm{\ω}}_j,t)/N---\left(10\right)$

以下选择两级行星轮系加一级平行轴式的齿轮箱为对象说明本发明的具体实施方法。

采集于高速轴的振动信号样本数据中共有100000个采样点，采样频率为6250Hz，时间长度为16秒。图1A中显示前50000个点，时间长度为8秒。与图1A中的振动数据同步采集的高速轴转速如图1B所示。其包含了齿轮和轴承的固有频率信号、啮合频率信号和大量的噪声信号，为了辨识异常振动，首先根据步骤1～4绘制能量谱。

步骤1：设原始信号为x(t)，根据公式(1)，利用EMD方法将x(t)分解。因为利用能量谱进行故障特征提取的时候着重于振动能量在时频空间的分布，即各个IMF函数的能量在时频空间的累积，所以不需要考虑EMD分解后各个IMF的混叠。设加入的噪声水平为0.1，ensemble次数为5，经EMD分解后得到的IMF如图2所示。

步骤2：从图2可以看出，步骤1共得到了12个IMF分量。令第k个IMF分量y_k(t)的包络函数为e_k(t)，利用公式(2)对y_k(t)进行循环归一化，直到第i次迭代过程中，信号的所有时刻的瞬时幅值都小于单位值，此时该IMF分量的F_k(t)和A_k(t)可利用公式(3)和公式(4)得出，以k＝5为例，该IMF分量的F₅(t)和A₅(t)分量分别为图3和图4所示。

步骤3：基于上面的A_k(t)和角频率ω(t)的定义，根据公式(6)，则原始测量信号的能量谱如图5所示。

步骤4：能量谱的归一化：因为风能的波动性特点，迫使风机机械传动系统以时变的速度运转，转速的变化进一步影响振动信号能量的时频分布，增加了故障诊断的难度。选择传动链中高速轴的转动频率为参考转速ω_ref(t)，用公式(7)计算振动信号无量纲能量谱(如图6所示)，消除转速变化对振动的影响，则无量纲能量谱代表信号在时频域上相对于参考转动频率ω_ref(t)的能量分布，方便分析振动信号的固有频率和啮合频率附近的幅频特征。

步骤5：啮合频率计算。实施方案针对的是采用两级行星轮系加一级平行轴式齿轮箱，齿轮箱参数如表1所示：

表1 齿轮箱参数表

因为振动信号采集于齿轮箱高速轴，所以振动频率主要由高速轴转速频率及倍频、平行级大小齿轮啮合频率构成。根据本实验的齿轮箱参数及机械振动原理可以得到下表中的齿轮箱各级传动机构的啮合频率：高速轴转速为r＝25转/秒，n₁、n₂分别为一级和一级行星齿轮的内齿圈齿数，n₃为平行级齿轮中小齿轮的齿数，并令G₁₂₃＝g₁·g₂·g₃，则根据齿轮箱参数及式(9)，可按照表2计算各级齿轮的啮合频率：

表2 齿轮箱各级机构啮合频率

步骤6：调制频率提取：如果用振动信号幅值极值的拟合曲线作为振动信号的包络函数，不仅不能保证振动信号在时间空间上的所有幅值都小于单位值，而且将会丢失原始信号中大量的谐波信息，而这些谐波信息正是风机齿轮箱故障信息的重要组成部分，因此根据公式(4)用多个包络函数完成振动信号的归一化，然后利用全部包络函数的乘积还原振动信号的幅值。

正常齿轮的信号反映在功率上，有啮合频率及其谐波分量，且以啮合频率成分为主，其高次谐波依次减小；同时，在低频处有齿轮轴旋转频率及其高次谐波。齿轮发生故障时，通常会使其正弦波式的啮合波形遭到破坏。从图6中可以看出，信号中含有明显的低频周期变化成分，频率上表现为以齿轮啮合频率及其谐波为载波频率，齿轮所在轴转频及其倍频为调制频率的啮合频率调制。其振动能量有较大幅度的波动，这是因为断齿故障的时域表现为幅值很大的冲击型振动，其振动频率等于发生断齿的轴的转频。

啮合振动信号的幅值调制反映的是振动信号的低频特征，通过频谱分析，能够容易地得到调制信号的成分构成。为得到调制信号的时频特性，根据表中齿轮箱各特征频率，按照表2中列出的振动信号成分，利用式(8)和式(9)介绍的特征频率计算方法，计算能量谱在这些特征频率附近的时间和频率特性。在图6中，在平行级齿轮的啮合频率，即高速轴转速的25倍频率附近，选取m组特定频率的能量信号，用频谱分析法分析其时间分布情况，分别得到m个调制信号的频率f_H(ω_j,t)，j＝1,…,m，然后根据式(10)计算调制信号的平均周期。取m＝3，表3是各调制频率分别在Hilbert能量谱中24、25、26倍频处的幅值，表示相应倍频处有该调制频率的周期性谐振出现，取表中的最低频率0.69，即为该振动信号的最低调制频率，其倒数即为数据采集的最短时间长度。

表3 基于HHT的调制频率分析结果

由于采用了以上技术方案，通过本发明的基于时频能量谱提取风机状态信号调制成分的方法，可以在振动信号周期特性的估计过程中，消除风机转速变化对故障信号提取的影响；可以估计为了进行信号分析所需要的状态信号采集的最短时间长度，能够保证在不丢失状态信号特征的同时获得较高的检测精度及较少的数据长度。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，本领域技术人员利用上述揭示的技术内容做出些许简单修改、等同变化或修饰，均落在本发明的保护范围内。

Claims (5)

1.一种基于时频能量谱提取风机状态信号调制成分的方法，用于获得目标振动信号中指定频率的调制成分，其特征在于，包括以下步骤：

A.采集目标振动信号，绘制目标振动信号的希尔伯特-黄谱,并采用瞬时无量纲频率归一化方法将希尔伯特-黄谱转换为无量纲频率信号能量谱；

B.从所述无量纲频率信号能量谱中提取频率在目标振动信号中指定频率附近的多组振动信号的无量纲能量，并采用频谱分析法获得所述多组无量纲能量对应的调制信号频率，然后将这些调制信号频率的均值作为所述目标振动信号中指定频率的调制成分的调制频率。

2.根据权利要求1所述的基于时频能量谱提取风机状态信号调制成分的方法，其特征在于，所述步骤具体包括：

A.采集目标振动信号x(t)；

用EMD方法将x(t)分解为公式(1)的形式：

$x\left(t\right)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{y}_k\left(t\right)-r_m\left(t\right)---\left(1\right)$

其中y_k(t)是x(t)的第k个IMF分量，并且满足单分量信号条件；r_m(t)是经过m次分解后的残差；

利用公式(2)对所述IMF分量y_k(t)进行循环归一化，直到第i次迭代过程中，IMF分量的所有时刻的瞬时幅值都小于单位值，

$y_k^i\left(t\right)=y_k^{i-1}\left(t\right)/e_k^{i-1}\left(t\right)---\left(2\right)$

其中，i∈[1,n]为迭代次数；为的包络函数；

于是，所述IMF分量被分解为独立的载波函数F_k(t)和调制函数A_k(t)，并有：

$F_k\left(t\right)=y_k^n\left(t\right)---\left(3\right)$

$A_k\left(t\right)=e_k^1\left(t\right)\·e_k^2\left(t\right)...e_k^n\left(t\right)---\left(4\right)$

基于所述载波函数F_k(t)将IMF分量的瞬时频率定义为：

${\mathrm{\ω}}_k\left(t\right)=\frac{1}{2\mathrm{\π}}\frac{{\mathrm{dF}}_k\left(t\right)}{\mathrm{dt}}---\left(5\right)$

其中ω_k(t)为第k个IMF分量的瞬时频率；

基于所述A_k(t)和ω_k(t)，得到所述振动信号的希尔伯特-黄谱H(ω,t)：

$H(\mathrm{\ω},t)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{A}_k\left(t\right)\cos (\∫2\mathrm{\π}{\mathrm{\ω}}_k\left(t\right)\mathrm{dt})---\left(6\right)$

选择风机上某测量点的转动频率ω_ref(t)为参考，用公式(7)计算所述无量纲频率信号能量谱H(ω,t)：

$H(\mathrm{\ω},t)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{A}_k\left(t\right)\cos (\∫2\mathrm{\π}\frac{{\mathrm{\ω}}_k\left(t\right)}{{\mathrm{\ω}}_{\mathrm{ref}}\left(t\right)}\mathrm{dt})---\left(7\right)$

B.根据公式(10)计算目标振动信号中指定频率的调制成分的调制频率f_F：

$f_F\left(t\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{f}_H({\mathrm{\ω}}_j,t)/N---\left(10\right)$

其中N为提取的无量纲能量的组数，f_H(ω_j,t)为在ω_j处的振动信号能量H(ω_j,t)的包络线的瞬时频率，H(ω_j,t)为无量纲频率信号能量。

3.根据权利要求2所述的基于时频能量谱提取风机状态信号调制成分的方法，其特征在于，所述步骤A中选择风机传动链中高速轴或低速轴的转动频率ω_ref(t)为参考。

4.根据权利要求1所述的基于时频能量谱提取风机状态信号调制成分的方法，其特征在于，所述目标振动信号为风机齿轮箱高速轴振动信号，所述步骤B中的指定频率为风机齿轮箱的啮合频率。

5.根据权利要求4所述的基于时频能量谱提取风机状态信号调制成分的方法，其特征在于，所述风机齿轮箱的啮合频率f_m根据公式(9)计算获得：

f_m＝f_i·n_i＝f_i+1·n_i+1 (9)

其中f_i和f_i+1分别为以啮合方式连接的两个传动轴的转动频率，n_i和n_i+1分别为其相应的齿轮齿数。