CN116010806B - 旋转机械时变多分量信号的时频分析方法 - Google Patents
旋转机械时变多分量信号的时频分析方法 Download PDFInfo
- Publication number
- CN116010806B CN116010806B CN202310313811.6A CN202310313811A CN116010806B CN 116010806 B CN116010806 B CN 116010806B CN 202310313811 A CN202310313811 A CN 202310313811A CN 116010806 B CN116010806 B CN 116010806B
- Authority
- CN
- China
- Prior art keywords
- time
- frequency
- instantaneous
- component
- signal
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000004458 analytical method Methods 0.000 title claims abstract description 81
- 239000011159 matrix material Substances 0.000 claims abstract description 65
- 238000000034 method Methods 0.000 claims abstract description 44
- 238000001228 spectrum Methods 0.000 claims abstract description 33
- 238000000605 extraction Methods 0.000 claims abstract description 20
- 238000009826 distribution Methods 0.000 claims abstract description 15
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 11
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 10
- 238000012545 processing Methods 0.000 claims abstract description 10
- 230000009467 reduction Effects 0.000 claims abstract description 8
- 230000001360 synchronised effect Effects 0.000 claims description 15
- 230000002441 reversible effect Effects 0.000 claims description 10
- 230000001052 transient effect Effects 0.000 claims description 6
- 238000006073 displacement reaction Methods 0.000 claims description 3
- 238000005457 optimization Methods 0.000 claims description 3
- 230000008707 rearrangement Effects 0.000 claims description 3
- 230000002776 aggregation Effects 0.000 abstract description 5
- 238000004220 aggregation Methods 0.000 abstract description 5
- 238000010586 diagram Methods 0.000 description 17
- 101001095088 Homo sapiens Melanoma antigen preferentially expressed in tumors Proteins 0.000 description 5
- 102100037020 Melanoma antigen preferentially expressed in tumors Human genes 0.000 description 5
- 230000006835 compression Effects 0.000 description 5
- 238000007906 compression Methods 0.000 description 5
- 230000009466 transformation Effects 0.000 description 5
- 230000008569 process Effects 0.000 description 4
- 238000004364 calculation method Methods 0.000 description 2
- 238000013016 damping Methods 0.000 description 2
- 238000003745 diagnosis Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 230000005654 stationary process Effects 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000012805 post-processing Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000010349 pulsation Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 239000000523 sample Substances 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
- 238000012706 support-vector machine Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
本发明属于信号处理领域，尤其涉及一种旋转机械时变多分量信号的时频分析方法，包括：1）获取待分析的旋转机械时变多分量信号；2）基于经验模态分解的多尺度小波阈值降噪算法处理原始信号；3）基于短时傅里叶变换对降噪后的信号进行初始化处理；4）通过EDSE的时频脊线提取方法重构初始时频分布结果，基于重构时频矩阵和最值搜索方法提取时频脊线并获得其瞬时频率估计值；5）基于Vold‑Kalman滤波器和瞬时频率估计值从降噪后的信号中分离谐波分量，计算谐波分量轴心轨迹的瞬时全谱参数；6）对旋转机械时变多分量信号的时频进行分析。本发明具有分辨率高、时频能量聚集性高以及在时频图谱上可观测到进动方向的优点。
Description
技术领域
本发明属于信号处理领域，涉及一种信号分析处理方法，尤其涉及一种旋转机械时变多分量信号的时频分析方法。
背景技术
由于旋转机械运行条件的多变性，其信号中普遍存在非平稳特性，而时频分析方法是表征非平稳信号时变特征的有效工具，经典时频分析方法包括短时傅里叶变换、小波变换和Wigner-Ville分布。高传昌等人基于经典时频分析方法分析了结构参数对压力脉动的影响，然而受海森堡测不准原理和交叉项的限制，经典时频分析方法的时频分辨率较低，时频谱的能量聚集性有限，难以实现非平稳信号的精确描述。近年来涌现出一系列后处理方法来提高这些经典方法的时频分辨率，如重分配方法、同步压缩变换等。
针对旋转机械信号特征复杂，提出应用同源多传感器信号的联合信息更加全面地提供相关旋转机械健康状况重要信息，吴峰崎等通过将全谱方法扩展并结合级联全谱和支持向量机实现了旋转机械多种故障的分析和诊断。温广瑞等将 Vold-Kalman 滤波阶次跟踪技术和分数阶傅里叶变换与全息谱结合实现转子启停车故障特征的提取。王丽雅等将全矢谱技术与短时傅里叶变换相结合提出了短时矢谱的概念，分析了压缩机喘振故障时振动矢量信号短时能量时频分布的规律。综上所述，时频分析方法结合以全频谱和全息谱为代表的信息融合技术可以实现转子系统非平稳振动信号分析，但大多数方法都是以短时窗稳定为假设，难以实现对机组非平稳过程的准确描述。
发明内容
为了解决背景技术中存在的上述技术问题，本发明提供了一种分辨率高、时频能量聚集性高以及在时频图谱上可观测到进动方向的旋转机械时变多分量信号的时频分析方法。
为了实现上述目的，本发明采用如下技术方案：
一种旋转机械时变多分量信号的时频分析方法，其特征在于：所述旋转机械时变多分量信号的时频分析方法包括以下步骤：
1）获取待分析的旋转机械时变多分量信号，所述待分析的旋转机械时变多分量信号是原始信号；
2）基于经验模态分解的多尺度小波阈值降噪算法处理原始信号，得到降噪后的信号；
3）基于短时傅里叶变换对降噪后的信号进行初始化处理，获得初始时频分布结果；
4）通过EDSE的时频脊线提取方法对步骤3）获取得到的初始时频分布结果重构时频矩阵，随后通过最值搜索方法提取时频脊线并获得其瞬时频率估计值；
5）基于Vold-Kalman滤波器和步骤4）所得到的瞬时频率估计值从步骤2）所得到的降噪后的信号中分离谐波分量，计算谐波分量轴心轨迹的瞬时全谱参数；
6）根据步骤5）所得到的谐波分量的瞬时全谱参数对旋转机械时变多分量信号的时频进行分析。
作为优选，本发明所采用的步骤1）中待分析的旋转机械时变多分量信号通过采用两只互相垂直安装的非接触式位移传感器对旋转机械轴系振动测量获得，所述原始信号包括x(t)和y(t)，其中，t表示全局时间变量。
作为优选，本发明所采用的步骤2）的具体实现方式是：
2.1）通过经验模态分解算法将原始信号x(t)分解为本征模态分量IMF k （k=1…K）和残差分量r(t)；
2.2）结合小波阈值降噪方法完成对本征模态分量以及残差分量的降噪处理，得到降噪后的信号，记作，
，
其中：表示降噪后的信号，/>表示降噪后的本征模态分量，/>表示降噪后的残差分量；
2.3）重复步骤2.1）以及步骤2.2），获取原始信号y(t)的降噪后的信号。
作为优选，本发明所采用的步骤4）的具体实现方式是：
4.1）基于步骤3）获得的初始时频分布结果，采用同步提取算子重构时频矩阵，获得初始时频矩阵重分配结果的视频矩阵；
4.2）根据基于欧氏距离迭代初始时频矩阵重分配结果的时频矩阵，获得优化重排后的重构时频矩阵；
4.3）通过最值搜索方法对步骤4.2）获取得到的重构时频矩阵进行提取，获得时频脊线和瞬时频率估计值。
作为优选，本发明所采用的步骤4.1）的具体实现方式是：
将步骤3）获得的初始时频分布结果经同步提取算子压缩得到待处理的时频分析结果/>，即得
，
其中：表示迪利克雷函数；ω表示频率；ω 0(t, ω)表示谐波分量瞬时频率值；
将待处理的时频分析结果离散化处理得到时频矩阵TFM，该矩阵列代表时间步长，行表示频率，时频矩阵TFM如下所示；
。
作为优选，本发明所采用的步骤4.2）的具体实现方式是：
根据时频矩阵TFM坐标值与根据欧式距离计算的惩罚值重构坐标值，具体是：
4.2.1）取时频矩阵TFM的；
4.2.2）更新：对最后一列应用惩罚产生最后一列的最小值/>，将最后一列中的最小值与当前索引值相加/>；
，
，
其中：是时频矩阵中（x 1，y 1）和（x 2，y 2）两个点之间欧氏距离；
表示当前索引最小值；/>表示更新前坐标值；/>表示更新后坐标值；/>为惩罚系数；
4.2.3）重复步骤4.2.1）以及步骤4.2.2）直至对时频矩阵TFM的全部更新，得到基于EDSE方法的重构时频矩阵ST 1，重构时频矩阵ST 1表达式为：
。
作为优选，本发明所采用的步骤4.3）的具体实现方式是：
基于最值搜索方法，对重构时频矩阵ST 1的延频域方向进行搜索，找出最小值，跟踪最小值坐标索引，最小值坐标索引构成了第一条时频脊线，所述时频脊线离散化的点表示t时刻瞬时频率估计值；
，
式中：表示第i条时频脊线；/>表示第i次更新后的时频矩阵；表示第i次重构的时频矩阵。
作为优选，本发明所采用的步骤5）的具体实现方式是：
基于Vold-Kalman滤波器和步骤4.3）获取得到的瞬时频率估计值，在降噪后信号中分离谐波分量，并计算谐波分量的瞬时全谱参数，所述谐波分量的瞬时全谱参数包括谐波分量的瞬时轨道的长轴、谐波分量的瞬时轨道的短轴、瞬时轨道的正向分量、瞬时轨道的反向分量以及进动方向系数；
设转子在互相垂直方向上的运动方程为：
，
式中：和/>分别是谐波分量；/>（i=1，2，···）是若干时变的谐波分量瞬时频率；
和/>分别是谐波/>在两方向上的幅值；
和/>分别是互相垂直方向上谐波/>的测量信号相位；
所述谐波分量的瞬时轨道的长轴表示为：
，
所述谐波分量的瞬时轨道的短轴表示为：
，
，
其中：、/>、/>和/>；
和/>分别是谐波/>在两方向上的幅值；
和/>分别是互相垂直方向上谐波/>的测量信号相位；
所述瞬时轨道的正向分量以及瞬时轨道的反向分量表示为：
，
所述进动方向系数SDI表示为：
，
依据全谱理论分析可知，轴心瞬时轨迹的进动方向由瞬时轴心轨迹的瞬时特征正向进动分量和反向进动分量决定，当>/>时，轴心瞬时轨迹进动方向为正向，1>SDI>0；反之，则轴心瞬时轨迹进动方向为反向, 此时0>SDI>-1；此外还有两种特殊状态，当/>为零时，谐波轴心瞬时轨迹为反向进动的正圆，此时SDI=-1；当/>为零时，谐波轴心瞬时轨迹为正向进动的正圆，此时SDI=1。
作为优选，本发明所采用的步骤6）的具体实现方式是：
6.1）获取步骤5）得到的瞬时轨道的正向分量以及瞬时轨道的反向分量；
6.2）将瞬时轨道的正向分量以及瞬时轨道的反向分量投影到时频谱三维坐标轴，得到谐波时频分析结果，所述谐波时频分析结果表达式是：
，
其中：瞬时频率使用的是谐波分量/>和谐波分量/>的瞬时频率的均值；
6.3）将步骤6.2）获取得到的各谐波分量时频谱图叠加，得到旋转机械时变多分量信号的分析结果，具体表示为：
，
其中：表示旋转机械时变多分量信号时频分析结果；/>表示谐波时频分析结果。
与现有技术相比，本发明的优点在于：
本发明在传统时频分析方法基础上，引入了欧几里得距离，重构时频矩阵，提高信号的时频聚集性，而且使时频矩阵具有更好的极值点连贯性避免极值点偏移，进而更准确地检出时频脊线并估计信号谐波分量瞬时频率值。
为提取谐波信号分量提供精准的转子速度信号，并基于谐波分量瞬时特征参数构建转子轴心高分辨率时频分析。该方法实现了转子时变复杂信号的谐波分量提取，同时计算转子谐波信号的瞬时特征参数，可用于描述非平稳过程中转子谐波分量瞬时振动状态。
同时高分辨率的时频图提供了转子谐波分量的进动方向信息，清晰地展现转子不对中二倍频信号的特征信息，为旋转机械振动信号分析提供更直观的特征，有利于非平稳过程中旋转机械转子的状态监测。经本发明处理的信号，实现谐波分量的高精度重构，重构后的谐波分量与理论谐波分量信号相比，还原度高，本发明是一种旋转机械振动信号高分辨率时频分析方法，在本发明对信号分析方法的基础上，可进行振动趋势预测、故障诊断等其他科学研究。本发明公开了一种用于旋转机械时变多分量信号时频分析方法，所述方案包括：获取旋转机械轴系振动信号；基于短时傅里叶变换对轴系振动信号进行处理，获取初始时频分布结果；通过引入欧几里得距离和同步提取算子重构时频矩阵，进而更准确地检出时频脊线并估计信号谐波分量瞬时频率值；采用Vold-Kalman滤波器在时域信号中分离谐波分量，并计算谐波分量的瞬时全谱参数，构建旋转机械时变多分量信号的高分辨率时频分析。本发明提供了一种瞬时频率估计方法，估计值精度高。而且给出一种新的旋转机械振动信号时频分析方法，不仅提高时频分析结果的时频聚集性，而且提供了轴心轨迹进动方向信息。
附图说明
图1为本发明的流程图。
图2为本发明使用的仿真信号时域波形及频域波形。
图3为降噪后信号分析图。
图4为瞬时频率估计值对比图。
图5为1倍谐波分量时域分析图。
图6为1.5倍谐波分量时域分析图。
图7为真实谐波分量与拟合分量对比图。
图8为进动方向系数图。
图9为信号时频分析对比图。
图10为转子不对中信号分析图。
图11为一倍频谐波分量时域分析图。
图12为二倍频谐波分量时域分析图。
图13为三倍频谐波分量时域分析图。
图14为进动方向系数图。
图15为信号时频分析对比图。
具体实施方式
下面结合附图及具体实施例对本发明做进一步的介绍。
参见图1，本发明提供了一种旋转机械时变多分量信号的时频分析方法，该方法包括以下步骤：
1）获取一组待分析的旋转机械时变多分量信号，此组信号通过采用两只互相垂直安装的非接触式位移传感器对旋转机械轴系振动测量获得，记作原始信号x(t)和y(t)，t表示全局时间变量；
2）基于经验模态分解的多尺度小波阈值降噪算法分别处理原始信号x(t)和y(t)，首先通过经验模态分解算法将信号分解为一组本征模态分量IMF k (k=1…K)和一个残差分量r(t)，然后结合小波阈值降噪方法完成对信号本征模态分量和残差分量的降噪处理，得到降噪后的信号，记作和/>。以信号x(t)为例，具体表示如下：
，
表示降噪后的信号，/>表示降噪后的本征模态分量，/>表示降噪后的残差分量；
3）基于短时傅里叶变换对降噪后的信号进行处理，获得初始时频分布结果；
4）基于获得的初始时频分布结果，基于欧氏距离的同步提取方法（Euclideandistance based Synchronous extraction method，EDSE）重构时频矩阵，获得初始时频矩阵重分配结果；
基于欧氏距离的同步提取方法，迭代初始时频矩阵，获得优化重排后的重构时频矩阵，基于重构时频矩阵和最值搜索方法，获得时频脊线和瞬时频率估计值，具体表示：
以初始时频分布第一列为起点，逐列计算更新矩阵列向量，然后根据当前列更新值，继续计算后续列，直到计算时频矩阵完成。然后从新矩阵每一列，找出最小值其坐标就是时频脊线坐标，进而在时频矩阵中提取最大能量脊。以三阶矩阵为例详细阐述基于EDSE的时频脊线提取方法，算法步骤如下：
4.1）首先将初始时频分析结果经同步提取算子压缩得到待处理的时频分析结果/>，即得
（1）
表示迪利克雷函数，ω表示频率，ω 0 (t, ω)表示谐波分量瞬时频率值，将待处理的时频分析结果/>离散化处理得到时频矩阵TFM，该矩阵列代表时间步长，行表示频率，时频矩阵TFM如下所示：
。
4.2）根据时频矩阵坐标值与根据欧式距离计算的惩罚值重构坐标值，以达到寻优和限制频率跳跃的目的，以更新为例，对最后一列应用惩罚产生最后一列的最小值，将最后一列中的最小值与当前索引值相加/>，
（2）
（3）
式中：时频矩阵中（x 1，y 1）和（x 2，y 2）两个点之间欧氏距离；/>表示当前索引最小值；/>表示更新前坐标值； />表示更新后坐标值；/>为惩罚系数。
4.3）继续更新第1列中其余元素的值，如下所示：
4.4）重复步骤4.2）和步骤4.3），获取后续列的更新值，得到基于EDSE方法的重构时频矩阵，矩阵表达式为：
基于最值搜索方法，对初始时频矩阵重分配结果延频域方向进行搜索，找出最小值，跟踪最小值坐标索引，它们构成了组成时频脊线，时频脊线离散化的点表示t时刻瞬时频率估计值。
（4）
其中：表示第i条时频脊线；/>表示第i次更新后的时频矩阵；表示第i次重构的时频矩阵。
基于Vold-Kalman滤波器和前述步骤得到的瞬时频率估计值，在降噪后信号中分离谐波分量，并计算谐波分量的瞬时全谱参数，构建旋转机械时变多分量信号的高分辨率时频分析。
设转子在互相垂直方向上的运动方程为：
（5）
式中： 和/>信号谐波分量；信号 />（i=1，2，···）是若干时变的谐波分量瞬时频率；/>和/>分别是谐波/>在两方向上的幅值；/>和/>分别是互相垂直方向上谐波/>的测量信号相位；t是时间变量。
根据全频谱理论分析，式（5）可以写成
（6）
其中，
（7）
其中、/>、/>和/>。
以转子简谐运动轴心轨迹分析为例，根据式（7）三角函数关系与简谐运动轨迹联立可得：
（8）
即得
（9）
假设转子的振动由给定时刻t处的瞬时参数决定且暂态稳定，则时刻t处的第i次谐波瞬时轨道定义为，表示为
（10）
的最大值和最小值分别位于椭圆轨道的长轴和短轴的末端。
将式（9）带入式（10）得表示为
（11）
在函数一阶导数为0且二阶导数为正数时取到长轴值，二阶导数为负数时取到短轴值计算函数对/>的导数
（12）
（13）
椭圆长轴和短轴的坐标点应是上式极值点处取到，即：
（14）
得到函数在极值点的值为
（15）
轨道长轴如下所示
（16）
椭圆的长轴垂直于短轴，短轴表示为
（17）
瞬时轴轨道正向进动分量和反向进动分量表示为
（18）
进动方向系数SDI表示为：
（19）
据全谱理论分析可知，轴心瞬时轨迹的进动方向由瞬时轴心轨迹的瞬时特征正向进动分量和反向进动分量决定，当>/>时，轴心瞬时轨迹进动方向为正向，1>SDI>0；反之，则轴心瞬时轨迹进动方向为反向, 此时0>SDI>-1；此外还有两种特殊状态，当/>为零时，谐波轴心瞬时轨迹为反向进动的正圆，此时SDI=-1；当/>为零时，谐波轴心瞬时轨迹为正向进动的正圆，此时SDI=1。
基于全频谱的轴心轨迹瞬时特征参数构建转子信号谐波分量的高分辨率时频分析
，
需要注意的是，此处瞬时频率使用的是谐波分量/>和谐波分量/>瞬时频率的均值，由于两只探头测量的是同一转轴的振动，各个阶次的瞬时频率理论上保持一致。
所述构建旋转机械时变多分量信号的高分辨率时频分析采用各谐波分量时频谱图叠加方式求得，具体表示为：
（20）
其中表示旋转机械时变多分量信号时频分析结果，/>表示谐波时频分析结果。
实施例1：
参见图1，一种旋转机械时变多分量信号的时频分析方法，包括以下步骤：
旋转机械时变多分量模拟信号包含两个非平稳分量，它们模拟旋转机械停机过程中转子的振动响应，用数学表达式为：
，
式中：x表示为转子横轴振动信号；y表示为转子纵轴振动信号；信号x谐波分量阻尼系数被设置为0.6；信号y谐波分量阻尼系数/>被设置为0.5；/>、，分别表示两个分量的瞬时频率；x和y的相位差由/>表示，分别设置、/>。x和y具有相同的瞬时频率分量。产生两个合成信号，采样频率为2048 Hz，采样时间为1s，在信号中人工加入高斯噪声，信号的信噪比为0.78dB。噪声信号如图2所示，图2中的（a）图和（b）图是信号x和y的时域波形图，由快速傅里叶变换获得的x和y的频域波形如图2中的（c）图和（d）图所示，可以看出噪声很强。
首先采用基于经验模态分解的多尺度小波阈值降噪算法对信号进行处理，分别对信号x、y进行经验模态分解得到一组本征模态分量和残差分量，然后对各个分量进行小波分解，采用软阈值方法对各个分量进行处理，得到降噪后的信号的时域波形如图3中的（a）图和（b）图所示，由频域波形图3中的（c）图和（d）图观察可以看出高频噪声被明显抑制，通过短时傅里叶变换得到的x、y的时频表示如图3中的（e）图和（f）图所示，可以观察到信号由两个谐波分量构成，然而每个谐波分量轨迹方向信息无法确定。从图3中的（g）轴心轨迹图和图3中的（h）全谱图中无法观测到更多信息。
在信号时频图谱中通过基于EDSE的时频脊线提取方法提取时频脊线并获取瞬时频率估计值如图4所示（实线表示实际值，细虚线表示EDSE提取结果，粗虚线表示同步压缩变换提取结果），精度较同步压缩变换有明显优势，与实际值较为接近。
将瞬时频率估计值作为转速信号应用于Vold-Kalman滤波器提取信号1倍和1.5倍谐波分量，图5中的（a）图表示信号x1倍谐波分量C1，图5中的（b）图表示信号y1倍谐波分量C1，图5中的（c）图表示信号倍谐波分量轴心轨迹图，图6中的（a）图表示信号x1.5倍谐波分量C2，图6中的（b）图表示信号y1.5倍谐波分量C2，图6中的（c）图表示1.5倍谐波分量轴心轨迹图，将提取到的谐波分量波形与理论谐波波形对比进行对比，清晰观察出提取到的谐波分量分量与理论谐波分量拟合度较高，表明本文所提出的方法对谐波信号的提取精度较高，噪声鲁棒性高，提取到的谐波分量与理论谐波分量对比图如图7所示。
进一步计算各阶谐波分量的瞬时轨道的长轴、短轴、瞬时轨道的正反向分量。
为验证信号时频表示能够展现旋转机械信号谐波分量的进动方向，定义振动信号的进动方向系数计算公式为，其中/>、/>表示x、y谐波信号相位，进动方向系数结果如图8所示，左侧为含噪信号方向系数图，右侧为不含噪信号方向系数图。为验证方法的性能，通过定义特征平均绝对百分比误差（Mean AbsolutePercentage Error， MAPE）和均方根误差（Root Mean Square Error，RMSE）指标对含噪各分量进动方向系数的各项指标进行评价。
，
，
式中：N为信号长度，下标r表示实际值，下标e为估计值。
含噪信号1倍谐波MAPE为5.5%、RMSE为0.0057，1.5倍谐波分量MAPE为11.3%、RMSE为0.0496，不含噪信号1倍谐波MAPE为0.37%、RMSE为2.38×e-4，1.5倍谐波分量MAPE为0.39%、RMSE为0.0017，通过特征指标表明该方法重构的信号谐波分量轴心轨迹进动方向系数与理论值较为一致，能精准表达旋转机械振动状态信息。
最后，将各个谐波分量轴轨道的正向分量和反向分量振幅投影到时频图三维坐标中，得到谐波分量轴轨道时频分析。原始信号的转子信号高分辨率时频分析由两个谐波分量的高分辨率时频分析叠加得到。为验证方法的准确性和有效性，本文将短时傅里叶变换如图9中的（a）图、同步压缩变换的时频分析如图9中的（b）图、理论上的时频分析如图9中的（c）图所示和本文所提出方法的结果如图9中的（d）图进行对比，本文所提方法的时频分析具有良好的时频分辨率，C1和C2分量的时频图可以清晰地分辨出来，与理论时频分析非常接近，并且从时频谱图上可以清晰地观察出各谐波分量的进动方向信息与SDI计算结果一致。
实施例2：本实例是为了验证本方法在旋转机械轴心现场应用
选取了某旋转机械发生转子不对中时转子振动信号进行分析说明，该信号x和y采集时转子转速为3180r/min（53Hz），信号x和y时域分析图如图10中的（a）图和（b）图所示。实施方法与实施例1相同。从图10中的（b）图和（c）图信号频谱分析，信号主要以一倍频和二倍频分量为主，在x和y的频谱上二倍频分量幅值大于一倍频。从由短时傅里叶变换得到的时频图（d）和（f）中，可以判断信号x和y分别由三个分量。从全频谱图10（h）直接看出信号一倍频分量是反进动，三倍频分量是正进动除此之外无法直观获取更多有价值信息，从图10中的（g）轨迹图无法判断转子任何状态信息。
接下来，使用本文所提出的算法对上述信号进行分析，图11中的（a）图表示信号x一倍谐波分量C1，图11中的（b）图表示信号y一倍谐波分量C1，图11中的（c）图表示信号一倍谐波分量轴心轨迹图，图12中的（a）图表示信号x二倍谐波分量C2，图12中的（b）图表示信号y二倍谐波分量C2，图12中的（c）图表示二倍谐波分量轴心轨迹图，图13中的（a）图表示信号x三倍谐波分量C2，图13中的（b）图表示信号y三倍谐波分量C2，图13中的（c）图表示三倍谐波分量轴心轨迹图。从图中可以观察到，谐波分量可以很好的分离出来，一倍频和三倍频谐波分量是椭圆形，且信号二倍频谐波分量的轴心轨迹符合转子不对中信号二倍频谐波分量特征。
计算谐波分量瞬时特征参数，进一步通过计算进动方向系数（如图14所示），获取谐波分量的进动方向系数，C3始终保持正进动，C1始终保持负进动，符合全频谱所观察到的信息，C2进动方向不确定符合转子不对中二倍频谐波分量特征。接下来，构建高分辨率时频分析谱图（如图15中的（d）图所示），并给出传统时频分析方法短时傅里叶变换如图15中的（a）图、同步压缩变换的时频分析如图15中的（b）图、小波变换的时频分析如图15中的（c）图所示作相比，此图清晰地展示了信号高分辨率时频表示谱图，并且提供了谐波分量C1、C2和C3的进动方向及信号时频域的变化信息。
Claims (1)
1.一种旋转机械时变多分量信号的时频分析方法，其特征在于：所述旋转机械时变多分量信号的时频分析方法包括以下步骤：
1)获取待分析的旋转机械时变多分量信号，所述待分析的旋转机械时变多分量信号是原始信号，其中待分析的旋转机械时变多分量信号通过采用两只互相垂直安装的非接触式位移传感器对旋转机械轴系振动测量获得，所述原始信号包括x(t)和y(t)，其中，t表示全局时间变量；
2)基于经验模态分解的多尺度小波阈值降噪算法处理原始信号，得到降噪后的信号；具体实现方式是：
2.1)通过经验模态分解算法将原始信号x(t)分解为本征模态分量IMFk和残差分量r(t)，k＝1…K；
2.2)结合小波阈值降噪方法完成对本征模态分量以及残差分量的降噪处理，得到降噪后的信号，记作x′(t)，
其中：x′(t)表示降噪后的信号，IMFk′表示降噪后的本征模态分量，r(t)′表示降噪后的残差分量；
2.3)重复步骤2.1)以及步骤2.2)，获取原始信号y(t)的降噪后的信号y′(t)；
3)基于短时傅里叶变换对降噪后的信号进行初始化处理，获得初始时频分布结果Ge(t,ω)；
4)通过基于欧氏距离的同步提取方法EDSE的时频脊线提取方法对步骤3)获取得到的初始时频分布结果Ge(t,ω)重构时频矩阵，随后通过最值搜索方法提取时频脊线并获得其瞬时频率估计值；具体实现方式是：
4.1)基于步骤3)获得的初始时频分布结果Ge(t,ω)，采用同步提取算子重构时频矩阵，获得初始时频矩阵重分配结果的时频矩阵，具体是：
将步骤3)获得的初始时频分布结果Ge(t,ω)经同步提取算子压缩得到待处理的时频分析结果TFse(t,ω)，即得
TFse(t,ω)＝Ge(t,ω)·δ[ω-ω0(t,ω)]，
其中：δ[·]表示迪利克雷函数；ω表示频率；ω0(t,ω)表示谐波分量瞬时频率值；
将待处理的时频分析结果TFse(t,ω)离散化处理得到时频矩阵TFM，该矩阵列代表时间步长，行表示频率，时频矩阵TFM如下所示；
4.2)根据基于欧氏距离的迭代初始时频矩阵重分配结果的时频矩阵，获得优化重排后的重构时频矩阵，具体是：
根据时频矩阵TFM坐标值与根据欧式距离计算的惩罚值重构坐标值，具体是：
4.2.1)取时频矩阵TFM的a11；
4.2.2)更新a11：对最后一列应用惩罚产生最后一列的最小值tip11，将最后一列中的最小值与当前索引值相加得到a′11；
其中：
是时频矩阵中(x1，y1)和(x2，y2)两个点之间欧氏距离；tip11表示当前索引最小值；a11表示更新前坐标值；a′11表示更新后坐标值；p为惩罚系数；
4.2.3)重复步骤4.2.1)以及步骤4.2.2)直至对时频矩阵TFM的全部更新，得到基于EDSE方法的重构时频矩阵ST1，重构时频矩阵ST1表达式为：
4.3)通过最值搜索方法对步骤4.2)获取得到的重构时频矩阵进行提取，获得时频脊线和瞬时频率估计值，具体是：
基于最值搜索方法，对重构时频矩阵ST1的沿频域方向进行搜索，找出最小值，跟踪最小值坐标索引，最小值坐标索引构成了时频脊线，所述时频脊线离散化的点表示t时刻瞬时频率估计值；
式中：fi(t)表示第i条时频脊线；TFMi(t,f)表示第i次更新后的时频矩阵；STi(t,f)表示第i次重构的时频矩阵；
5)基于Vold-Kalman滤波器和步骤4)所得到的瞬时频率估计值从步骤2)所得到的降噪后的信号中分离谐波分量，计算谐波分量轴心轨迹的瞬时全谱参数，具体实现方式是：
基于Vold-Kalman滤波器和步骤4.3)获取得到的瞬时频率估计值，在降噪后信号中分离谐波分量，并计算谐波分量的瞬时全谱参数，所述谐波分量的瞬时全谱参数包括谐波分量的瞬时轨道的长轴、谐波分量的瞬时轨道的短轴、瞬时轨道的正向分量、瞬时轨道的反向分量以及进动方向系数；
设转子在互相垂直方向上的运动方程为：
式中：xi和yi分别是谐波分量；ωi(t)是若干时变的谐波分量瞬时频率，i＝1，2，···；ai(t)和bi(t)分别是谐波ωi(t)在两方向上的幅值；ɸxi(t)和ɸyi(t)分别是互相垂直方向上谐波ωi(t)的测量信号相位；
所述谐波分量的瞬时轨道的长轴表示为：
所述谐波分量的瞬时轨道的短轴表示为：
其中：ait＝ai(t)、bit＝bi(t)、ɸxit＝ɸxi(t)和ɸyit＝ɸyi(t)；
所述瞬时轨道的正向分量以及瞬时轨道的反向分量表示为：
所述进动方向系数SDI表示为：
SDI＝sin(ɸxit-ɸyit)，
依据全谱理论分析可知，轴心瞬时轨迹的进动方向由轴心瞬时轨迹的瞬时特征正向进动分量和反向进动分量决定，当Rit+>Rit-时，轴心瞬时轨迹进动方向为正向，1>SDI>0；反之，则轴心瞬时轨迹进动方向为反向,此时0>SDI>-1；此外还有两种特殊状态，当Rit+为零时，谐波轴心瞬时轨迹为反向进动的正圆，此时SDI＝-1；当Rit-为零时，谐波轴心瞬时轨迹为正向进动的正圆，此时SDI＝1；
6)根据步骤5)所得到的谐波分量的瞬时全谱参数对旋转机械时变多分量信号的时频进行分析，具体实现方式是：
6.1)获取步骤5)得到的瞬时轨道的正向分量以及瞬时轨道的反向分量；
6.2)将瞬时轨道的正向分量以及瞬时轨道的反向分量投影到时频谱三维坐标轴，得到谐波时频分析结果，所述谐波时频分析结果表达式是：
TFRi(t，ω)＝Rit+δ[ω-ωi(t)]+Rit-δ[-ω-ωi(t)]
其中：瞬时频率ωi(t)使用的是谐波分量xi和谐波分量yi的瞬时频率的均值；
6.3)将步骤6.2)获取得到的各谐波分量时频分析结果叠加，得到旋转机械时变多分量信号的分析结果，具体表示为：
其中：TFR(t,ω)表示旋转机械时变多分量信号时频分析结果；TFRi(t,ω)表示谐波时频分析结果。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310313811.6A CN116010806B (zh) | 2023-03-28 | 2023-03-28 | 旋转机械时变多分量信号的时频分析方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310313811.6A CN116010806B (zh) | 2023-03-28 | 2023-03-28 | 旋转机械时变多分量信号的时频分析方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116010806A CN116010806A (zh) | 2023-04-25 |
| CN116010806B true CN116010806B (zh) | 2023-09-01 |
Family
ID=86025294
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310313811.6A Active CN116010806B (zh) | 2023-03-28 | 2023-03-28 | 旋转机械时变多分量信号的时频分析方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116010806B (zh) |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110926594A (zh) * | 2019-11-22 | 2020-03-27 | 北京科技大学 | 一种旋转机械信号时变频率特征提取方法 |
| CN111879508A (zh) * | 2020-07-28 | 2020-11-03 | 无锡迈斯德智能测控技术有限公司 | 基于时频变换的旋转机械瞬时转速估计方法、装置和存储介质 |
| WO2021068939A1 (zh) * | 2019-10-12 | 2021-04-15 | 中科新松有限公司 | 一种基于多分量信号分解的机械臂关节振动识别方法 |
| CN113125179A (zh) * | 2021-03-11 | 2021-07-16 | 同济大学 | 一种旋转机械转速波动无键相阶次跟踪方法 |
| CN113702043A (zh) * | 2021-08-09 | 2021-11-26 | 大连理工大学 | 一种基于povmd和fdtw的时变转速下行星齿轮箱故障诊断方法 |
| CN113988125A (zh) * | 2021-10-25 | 2022-01-28 | 西安交通大学 | 一种基于改进同步压缩变换的扭振信号瞬时频率提取方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20200051419A1 (en) * | 2017-10-11 | 2020-02-13 | Analog Devices Global Unlimited Company | Cloud-based machine health monitoring |
| EP3702796A1 (en) * | 2019-02-28 | 2020-09-02 | Julius-Maximilians-Universität Würzburg | System for the one-sided generation of magnetic fields for the multidimensional encoding of magnetic particles and method of operation thereof |
-
2023
- 2023-03-28 CN CN202310313811.6A patent/CN116010806B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2021068939A1 (zh) * | 2019-10-12 | 2021-04-15 | 中科新松有限公司 | 一种基于多分量信号分解的机械臂关节振动识别方法 |
| CN110926594A (zh) * | 2019-11-22 | 2020-03-27 | 北京科技大学 | 一种旋转机械信号时变频率特征提取方法 |
| CN111879508A (zh) * | 2020-07-28 | 2020-11-03 | 无锡迈斯德智能测控技术有限公司 | 基于时频变换的旋转机械瞬时转速估计方法、装置和存储介质 |
| CN113125179A (zh) * | 2021-03-11 | 2021-07-16 | 同济大学 | 一种旋转机械转速波动无键相阶次跟踪方法 |
| CN113702043A (zh) * | 2021-08-09 | 2021-11-26 | 大连理工大学 | 一种基于povmd和fdtw的时变转速下行星齿轮箱故障诊断方法 |
| CN113988125A (zh) * | 2021-10-25 | 2022-01-28 | 西安交通大学 | 一种基于改进同步压缩变换的扭振信号瞬时频率提取方法 |
Non-Patent Citations (1)
| Title |
|---|
| 付波等.基于EDM-VKF的转子时变复杂信号时频表示方法.《中国农村水利水电》.2023,第1-14页. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116010806A (zh) | 2023-04-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Yu et al. | Multisynchrosqueezing transform | |
| Li et al. | An improved local mean decomposition method based on improved composite interpolation envelope and its application in bearing fault feature extraction | |
| Yu | A multisynchrosqueezing-based high-resolution time-frequency analysis tool for the analysis of non-stationary signals | |
| CN104034412B (zh) | 一种基于分数阶全息原理的旋转机械故障特征提取方法 | |
| CN106501602B (zh) | 一种基于滑窗频谱分离的基波参数测量方法 | |
| CN114992033B (zh) | 基于nlm-ceemdan的水电机组信号去噪方法 | |
| CN113310684B (zh) | 基于尺度空间和改进稀疏表示的齿轮箱故障特征提取方法 | |
| CN109117896B (zh) | 一种基于ksvd字典学习的滚动轴承故障特征提取方法 | |
| Deng et al. | An improved spline-local mean decomposition and its application to vibration analysis of rotating machinery with rub-impact fault | |
| He et al. | Local maximum synchrosqueezes from entropy matching chirplet transform | |
| Chen et al. | Adaptive scale decomposition and weighted multikernel correntropy for wheelset axle box bearing diagnosis under impact interference | |
| CN111582128B (zh) | 一种基于狼群参数化联合字典的机械故障稀疏表示方法 | |
| Liu et al. | Two-step adaptive chirp mode decomposition for time-varying bearing fault diagnosis | |
| Xu et al. | Match-extracting chirplet transform with application to bearing fault diagnosis | |
| Dong et al. | Time–frequency-multisqueezing transform | |
| Lv et al. | Generalized synchroextracting-based stepwise demodulation transform and its application to fault diagnosis of rotating machinery | |
| CN112229627A (zh) | 基于短时稀疏傅里叶变换的旋转机械瞬时转速估计方法 | |
| CN116010806B (zh) | 旋转机械时变多分量信号的时频分析方法 | |
| Ding et al. | Slope synchronous chirplet transform and its application to tacho-less order tracking of rotating machineries | |
| Wang et al. | A novel time-frequency analysis method for fault diagnosis based on generalized S-transform and synchroextracting transform | |
| CN113281047A (zh) | 一种基于变尺度Lempel-Ziv的轴承内外圈故障定量趋势诊断方法 | |
| Liu et al. | An approximate maximum likelihood estimator for instantaneous frequency estimation of multicomponent nonstationary signals | |
| Zhong et al. | Measurement while drilling mud pulse signal denoising and extraction approach based on particle-swarm-optimized time-varying filtering empirical mode decomposition | |
| CN114098690A (zh) | 一种轻量级脉搏波信号噪声去除及心率检测方法 | |
| Zhou et al. | A study on fault diagnosis of rotating machinery combined wavelet transform with VMD |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |