一种改进的基于实倒频谱分析的旋转机械调制频率提取方法与流程

本发明涉及信号的频率成分提取领域，尤其涉及一种改进的基于实倒频谱分析的旋转机械调制频率提取方法。

背景技术：

信号的频率成分提取属于信号频域分析领域的内容，现有的应用于实际信号的频域分析方法主要有快速傅里叶变换、短时傅里叶变换以及小波变换等。快速傅里叶变换是利用wn因子的周期性和对称性，构造的一种快速、高效的算法，能够较为完美地描绘平稳信号；短时傅里叶变换是将信号加窗，将加窗后的信号进行快速傅里叶变换，能够描绘信号的二维时频图；小波变换也是对信号进行加窗分析，但其窗函数形状可以改变，能够保证信号的频域分析结果在低频处有较高的频率分辨率，在高频处有较高的时间分辨率，从而实现对信号的自适应分析。

这些方法能够将实际信号的时域信息转换成频域中的信息，但是针对于旋转机械的实际振动信号来说，以上三种方法有时候并不能识别出旋转机械振动中的一些基本调制频率成分，而这些基本调制频率包括了转子轴频、叶片叶频等极其关键的重要信息，因此非常有必要被识别出来。

《科技视界》(机械与电子)2014年第11期94-103页“倒频谱在振动测试中的应用”的文献公开了一种运用倒频谱分析技术识别频域调制信号的边频成分的方法。倒频谱具有谱线定位准确、幅值突出，可较好地识别频域调制信号的边频成分的特点。其缺点在于运用倒频谱定义对实际信号进行处理时，所获取的图像信息不够明显，曲线变化相对平缓，象征着调制频率的波峰也不够突出，难以有效地获得所需信息。

技术实现要素：

本发明提供了一种改进的基于实倒频谱分析的旋转机械调制频率提取方法，能够克服在旋转机械实际工作的大部分情况下，存在的调制频率在频谱图中不明显，不容易被识别的缺陷，能够在振动时域信号的实倒频谱分析图中简单、直接地识别旋转机械的调制频率，得到的图像信息更加突出鲜明，曲线变化更明显，波峰更突出。

一种改进的基于实倒频谱分析的旋转机械调制频率提取方法，包括：

步骤1，对旋转机械振动信号进行采集，获得离散的振动信号时序序列；

步骤2，对离散的振动信号时序序列进行实倒频谱处理，得到每个点所对应的实倒频谱分析值；

步骤3，以离散的振动信号时序序列每个点的时间倒数即对应的频率为横坐标，以每个点所对应的降噪后的倒频谱处理值为纵坐标，绘图得到振动时域信号的实倒频谱处理图；

步骤4，根据得到的实倒频谱处理图，识别出旋转机械的基本调制频率。

步骤1中，在对旋转机械振动信号进行采集时，振动信号采样频率与采样时间间隔的关系为：

式中：

fs为振动信号的采样频率；ts为振动信号的采样时间间隔；

振动信号采样总时长与信号频率分辨率的关系为：

式中：

f0为振动信号的频率分辨率；n为振动信号的采样总点数；t为振动信号的采样总时长。

步骤2中，在对离散的振动信号时序序列进行实倒频谱处理的基础上，对实倒频谱分析值随时间变化的曲线进行降噪处理，得到时序序列的每个点所对应的降噪后的实倒频谱处理值。

所述的降噪处理，在matlab中的语句为：

c(t)＝wden(c(t),′sqtwolog′,′s′,′mln′,3,′db2′)；

式中：

c(t)为振动信号时间序列上每个点所对应的降噪后的实倒频谱处理值；c(t)为振动信号时间序列上每个点所对应的快速傅里叶变换之后再取绝对值，再对其取自然对数，接着对其作快速傅里叶变换逆变换，然后再对其进行取实部操作得到的实倒频谱分析值；wden为小波降噪；sqtwolog代表采用固定阈值原则；s代表阈值函数选择方式为软阈值；mln代表根据每一层小波分解的噪声水平估计进行调整；3和db2代表利用db2的小波对信号进行3层分解。

步骤2中，所述的实倒频谱处理的具体步骤包括：

步骤2-1，对离散的振动信号时序序列进行快速傅里叶变换；

步骤2-2，在快速傅里叶变换的基础上，对所得结果取绝对值；

步骤2-3，在取绝对值的基础上，对所得结果取自然对数；

步骤2-4，在取自然对数的基础上，对所得结果进行快速傅里叶逆变换；

步骤2-5，在快速傅里叶逆变换的基础上，对其进行取实部的操作，得到振动信号时间序列上每个点所对应的实倒频谱分析值。

其中：步骤2-1中，所述的快速傅里叶变换，在matlab中的语句为：

fft-t＝fft(t)

式中：

fft-t为振动信号时间序列上每个点所对应的快速傅里叶变换结果；t为离散的原始振动信号时序序列；fft表示快速傅里叶正变换。

步骤2-2中，所述的取绝对值，在matlab中的语句为：

abs-fft-t＝abs(fft-t)

式中：

abs-fft-t为振动信号时间序列上每个点所对应的快速傅里叶变换之后再取绝对值的结果；abs表示取绝对值。

步骤2-3中，所述的取自然对数，在matlab中的语句为：

ln-abs-fft-t＝log(abs-fft-t)

式中：

ln-abs-fft-t为振动信号时间序列上每个点所对应的快速傅里叶变换之后再取绝对值，然后再对其取自然对数的结果；log表示取自然对数。

步骤2-4中，所述的快速傅里叶变换的逆变换，在matlab中的语句为：

ifft-ln-abs-fft-t＝ifft(ln-abs-fft-t)

式中：

ifft-ln-abs-fft-t为振动信号时间序列上每个点所对应的快速傅里叶变换之后再取绝对值，再对其取自然对数，然后再对其作快速傅里叶变换逆变换的结果；ifft表示快速傅里叶逆变换。

步骤2-5中，所述的取实部的操作在matlab中的语句为：

c(t)＝real(ifft-ln-abs-fft-t)

式中：

c(t)为振动信号时间序列上每个点所对应的快速傅里叶变换之后再取绝对值，再对其取自然对数，接着对其作快速傅里叶变换逆变换，然后再对其进行取实部操作得到的实倒频谱分析值；real表示对复数进行取实部的操作。

本发明提供的技术方案采用实倒频谱处理方法，只保留了信号的频谱幅值信息，摒弃了相位信息，实现了对旋转机械的转子轴频的有效提取，通过对实倒频谱处理值进行适当降噪处理，可以很大程度上消除实际振动信号因为噪声带来的干扰，解决了现有的快速傅里叶变换、短时傅里叶变换以及小波变换等方法不能显著识别出旋转机械振动中的一些基本调制频率成分的缺陷。

应用本发明能够在振动时域信号的实倒频谱分析图中简单、直接地识别旋转机械的调制频率，得到的图像信息更加突出鲜明，曲线变化更明显，波峰更突出。结合旋转机械的机器运行规律，能够进一步有效地识别出旋转机械的轴频、叶频，对于旋转机械的性能分析和状态监测具有重要意义。

附图说明

图1为原始的仿真信号的实倒频谱处理图；

图2为改进的仿真信号的实倒频谱处理图；

图3为加信噪比为20db的白噪声的仿真信号的实倒频谱处理图；

图4a为加15db噪声仿真信号的实倒频谱处理图；

图4b为加12db噪声仿真信号的实倒频谱处理图；

图4c为加9db噪声仿真信号的实倒频谱处理图；

图4d为加6db噪声仿真信号的实倒频谱处理图；

图4e为加3db噪声仿真信号的实倒频谱处理图；

图4f为加0db噪声仿真信号的实倒频谱处理图；

图5为小风机50hz正常工作的振动信号时域图；

图6为小风机10hz工作的实倒频谱处理图；

图7为小风机20hz工作的实倒频谱处理图；

图8为小风机30hz工作的实倒频谱处理图；

图9为小风机40hz工作的实倒频谱处理图；

图10为小风机50hz工作的实倒频谱处理图；

图11为小风机附重10克工作频率10hz的实倒频谱处理图；

图12为小风机附重10克工作频率20hz的实倒频谱处理图；

图13为小风机附重20克工作频率10hz的实倒频谱处理图；

图14为小风机附重20克工作频率20hz的实倒频谱处理图。

具体实施方式

为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步地详细描述。

本发明提出的改进的基于实倒频谱分析的旋转机械调制频率提取方法，分别对仿真信号和实际信号进行处理分析，具体步骤描述如下：

(1)根据实际需要选定振动信号的采样频率fs，进而确定采样时间间隔ts，选取恰当的采样总时长t，在此条件下，对旋转机械的振动信号进行测取。

振动信号采样频率与采样时间间隔的关系如公式(1)所示，

式中：

fs为振动信号的采样频率；

ts为振动信号的采样时间间隔；

振动信号采样总时长与信号频率分辨率的关系如公式(2)所示，

式中：

f0为振动信号的频率分辨率；

n为振动信号的采样总点数；

t为振动信号的采样总时长。

(2)在得到离散的振动信号时序序列的基础上，对其进行快速傅里叶变换，得到振动信号时间序列上每个点所对应的快速傅里叶变换结果。

matlab(用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境)中的语句为：

fft-t＝fft(t)

式中：

fft-t为振动信号时间序列上每个点所对应的快速傅里叶变换结果；

t为离散的原始振动信号时序序列；fft表示快速傅里叶正变换。

(3)在得到振动信号时间序列上每个点所对应的快速傅里叶变换结果的基础上，对其取绝对值，得到每个点所对应的值。

matlab中的语句为：

abs-fft-t＝abs(fft-t)

式中：

abs-fft-t为振动信号时间序列上每个点所对应的快速傅里叶变换之后再取绝对值的结果；abs表示取绝对值。

(4)在上一步计算得到振动信号时间序列上每个点所对应的值的基础上，对其取自然对数，得到每个点所对应的值。

matlab中的语句为：

ln-abs-fft-t＝log(abs-fft-t)

式中：

ln-abs-fft-t为振动信号时间序列上每个点所对应的快速傅里叶变换之后再取绝对值，然后再对其取自然对数的结果；log表示取自然对数。

(5)在上一步计算得到振动信号时间序列上每个点所对应的值的基础上，对其进行快速傅里叶变换的逆变换，得到每个点所对应的值。

matlab中的语句为：

ifft-ln-abs-fft-t＝ifft(ln-abs-fft-t)

式中：

ifft-ln-abs-fft-t为振动信号时间序列上每个点所对应的快速傅里叶变换之后再取绝对值，再对其取自然对数，然后再对其作快速傅里叶变换逆变换的结果；ifft表示快速傅里叶逆变换。

(6)在上一步计算得到振动信号时间序列上每个点所对应的值的基础上，对其进行取实部的操作，得到振动信号时间序列上每个点所对应的实倒频谱分析值。

matlab中的语句为：

c(t)＝real(ifft-ln-abs-fft-t)

式中：

c(t)为振动信号时间序列上每个点所对应的快速傅里叶变换之后再取绝对值，再对其取自然对数，接着对其作快速傅里叶变换逆变换，然后再对其进行取实部操作得到的实倒频谱分析值；real表示对复数进行取实部的操作。

(7)在上一步计算得到振动信号时间序列上每个点所对应的实倒频谱分析值的基础上，对其进行适当的降噪处理，得到降噪后每个点所对应的实倒频谱处理值。

matlab中的语句为：

c(t)＝wden(c(t),′sqtwolog′,′s′,′mln′,3,′db2′)；

式中：

c(t)为振动信号时间序列上每个点所对应的降噪后的实倒频谱处理值；wden是小波降噪，对c(t)进行小波降噪得到c(t)；sqtwolog代表采用固定阈值原则；s代表阈值函数选择方式为软阈值；mln代表根据每一层小波分解的噪声水平估计进行调整；3和db2代表利用db2的小波对信号进行3层分解。

(8)在计算得到振动信号时间序列上每个点所对应的降噪后的实倒频谱处理值的基础上，以振动信号时间序列上每个点的对应时间的倒数为横坐标，以上一步得到的振动信号时间序列上每个点所对应的降噪后的实倒频谱处理值为纵坐标，一一对应，进行绘图，得到振动信号的实倒频谱分析图(以频率为横坐标)，从而可以直观简单地识别、提取出振动信号中的调制频率。

为了凸显本发明的优越性，本实例将对原始的仿真信号、改进的仿真信号、加噪声的仿真信号以及实际信号、降噪后的实际信号进行基于实倒频谱技术的处理，将其中的调制频率成分识别出来。

(1)仿真信号为：

y(t)＝(cos(2π·10t)+cos(2π·20t)+cos(2π·30t))·sin(2π·100t)

这是一个调制基频为10hz的信号，同时还有10hz的二倍频、三倍频的调制频率，这些频率一同对100hz的频率进行调制，得到原始的仿真信号。

对原始的仿真信号进行实倒频谱处理后，得到图1，可以发现在10hz处出现了第一个较为明显的波峰，而10hz就是原始的仿真信号的调制基频。

(2)多项调制信号的检测

对(1)中仿真信号进行改进，使其调制信号中含有多种调制基频，并在改进的仿真信号上加上一个信噪比为20db的白噪声，再对其进行实倒频谱处理。改进后的仿真信号为：

y(t)＝(cos(2π·10t)+cos(2π·20t)+cos(2π·30t)+cos(2π·7t)+cos(2π·14t)+cos(2π·21t)+cos(2π·28t))·sin(2π·100t)

可以看出，对比于原来的仿真信号，调制信号在原来的10hz、20hz、30hz的基础上增加了7hz、14hz、21hz、28hz的频率成分，由10hz及其谐波和7hz及其谐波组成。对改进的仿真信号进行实倒频谱处理后，得到图2，可以发现，在7hz和10hz附近出现了两个最明显的波峰，这证明了本方法不仅能够识别单一的调制基频，还能同时识别出多个不同的调制基频。

(3)在原始的仿真信号基础上，加上一个信噪比为20db的白噪声，对加噪声的仿真信号进行实倒频谱处理后，得到图3，可以发现在10hz处存在明显的波峰，能够有效识别调制基频。

为了验证所加噪声的信噪比大小对于本发明识别调制频率的影响，在原始的仿真信号基础上，依次分别加上信噪比为15db、12db、9db、6db、3db、0db的白噪声，在对这些加噪声的仿真信号进行实倒频谱处理后，得到图4a～图4f，可以发现在10hz处依旧存在较为明显的波峰，证明本发明在存在一定的噪声干扰的情况下，依旧能够有效识别调制基频，但是随着所加噪声信噪比的降低，基本调制频率处的波峰逐渐变得不再是最明显的波峰，这也说明了识别基本调制频率变得越来越困难，但是即便在所加噪声的信噪比为0db的情况下，依旧能够对基本调制频率进行有效识别，说明了本发明能够有效提取出低信噪比信号中的调制频率。

(4)小风机正常工作的振动信号处理

对于旋转机械的实际数据分析处理，选取了小风机作为典型的旋转机械，确定采样频率fs为5120hz，采样总时长t为15s，对小风机以10hz、20hz、30hz、40hz、50hz的频率正常工作时的振动信号进行测取，其中小风机以50hz正常工作时的实际振动信号时域图如图5所示：

分别对小风机以10hz、20hz、30hz、40hz、50hz频率工作的振动信号进行实倒频谱处理后，得到图6至图10。

可以发现在图6和图7中，10hz和20hz左右处出现了最为明显的波峰，而此时小风机的工作频率就是10hz和20hz。换而言之，可以通过寻找实倒频谱图中出现的波峰，找到起主要调制作用的频率成分，而在小风机以较低的工作频率工作时，主要调制频率就是其转子频率。

观察图9和图10，可以发现当小风机以40hz或者50hz的频率工作时，对应的振动信号倒频谱图中，最明显波峰不再位于对应的转子频率处，而是出现在转子频率的二分之一倍频处，如在图9中，较为明显的第一个波峰出现在横坐标19.39hz处，恰好是工作频率40hz的一半左右，在图10中，较为明显的第一个波峰出现在横坐标24.98hz处，恰好是工作频率50hz的一半左右。这也就说明，当小风机以40hz、50hz等相对较高的工作频率运转时，其振动信号中对主要频率成分起到调制作用的基本调制频率不再是转子频率，而是转子频率的二分之一倍频。

将小风机分别以10hz、20hz、30hz、40hz、50hz频率工作的振动信号实倒频谱处理结果图进行横向对比，可以发现：当小风机以相对较低的工作频率运转时，如以10hz、20hz的工作频率工作时，可以通过实倒频谱识别其基本调制频率，此时其基本调制频率就是转子频率。随着工作频率的逐渐提高，基本调制频率逐渐地由转子频率变成转子频率的二分之一倍频，如以40hz、50hz的工作频率工作时，倒频谱识别出来的最明显的调制频率已经变成转子频率的二分之一倍频了。而以中间频率30hz工作时，如图8所示，在实倒频谱处理图上既能够找到代表转子频率的波峰，也能够在转子频率的二分之一倍频处找到波峰，但是此时两个波峰都不是特别明显，并且在还有一些其他高度和陡峭程度相当的波峰，这证明此时调制频率的主要成分比较分散，既包括了转子频率，还包括了转子频率的二分之一倍频以及其他整数谐波、分数倍频等多种频率成分，不是集中于某一主要频率成分。

(5)小风机转子不平衡的振动信号处理

转子不平衡故障的实质就是旋转部件质量的偏心，由于旋转部件的质心偏离旋转的轴线，当转子旋转时就会产生一个不平衡的离心惯性力，从而使得机械设备产生周期性的振动。为了实现转子不平衡的故障工况，在风机的某个叶片上增加质量已知的附重，使得转子发生质量偏心的情况，附重质量有10克和20克两种情况，风机的工作频率则选择了10hz和20hz两种情况，之所以没有选择更大的附重质量和更宽的工作频率范围，是因为在继续增加附重质量和调高工作频率的情况下，就会发生将附重甩出的危险状况，这是由于附重是由可以吸附在叶片上的已知重量的磁铁所营造，一旦磁铁所受到的离心力超过叶片对其的吸附力，就会将磁铁甩出。

叶片附重的转子不平衡实验中，共进行了四组实验，分别是附重10克工作频率10hz，附重10克工作频率20hz、附重20克工作频率10hz和附重20克20hz等四种情况，对其振动信号进行实倒频谱处理得到图11至图14。

对比图11和图12，图13和14可以发现，在附重相同的情况下，随着转子频率由10hz变为20hz，实倒频谱图中最明显的波峰也由10hz变为20hz，通过实倒频谱依旧能够直观、有效地读取出此时的基本调制频率——转子频率。然而也可以发现随着转子频率的增大，转子频率所代表的主波峰高度下降，其他次级波峰的数目和幅值都有所增加，这与小风机在正常情况下工作频率逐渐增大时的变化趋势也是一致的。

对比图11和图13，图12和图14可以发现，在频率相同时，随着附重质量由10克变成20克，实倒频谱中最明显的波峰横坐标位置依旧保持基本不变，但是主波峰的高度显著下降，与周围其他次级波峰的差距也变得不太明显，这就说明了调制频率的成分也变得不再主要集中于转子频率处，而是出现了更多的其他调制频率成分，这与工作频率增大带来的影响是比较接近的。

根据以上情况进行总结，在转子不平衡的工况下，随着旋转机械工作频率的增大、不平衡质量的增大，会导致调制频率由只存在一个主要的调制频率成分变为同时存在多种较显著的调制频率成分，而通过实倒频谱分析，可以有效地识别出旋转机械中这些主要的调制频率及其变化过程。

以上所述仅为本发明的优选实施方式，本发明的保护范围并不仅限于上述实施方式，凡是属于本发明原理的技术方案均属于本发明的保护范围。对于本领域的技术人员而言，在不脱离本发明的原理的前提下进行的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。