CN113565584A - 一种叶端定时信号时频滤波方法 - Google Patents
一种叶端定时信号时频滤波方法 Download PDFInfo
- Publication number
- CN113565584A CN113565584A CN202111017893.7A CN202111017893A CN113565584A CN 113565584 A CN113565584 A CN 113565584A CN 202111017893 A CN202111017893 A CN 202111017893A CN 113565584 A CN113565584 A CN 113565584A
- Authority
- CN
- China
- Prior art keywords
- frequency
- aliasing
- data
- index
- time
- 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.)
- Granted
Links
- 238000001914 filtration Methods 0.000 title claims abstract description 59
- 238000000034 method Methods 0.000 title claims abstract description 39
- 238000006243 chemical reaction Methods 0.000 claims abstract description 50
- 238000006073 displacement reaction Methods 0.000 claims abstract description 37
- 239000011159 matrix material Substances 0.000 claims abstract description 22
- 238000005070 sampling Methods 0.000 claims abstract description 17
- 239000013598 vector Substances 0.000 claims abstract description 16
- 238000010183 spectrum analysis Methods 0.000 claims abstract description 6
- 238000004458 analytical method Methods 0.000 claims abstract description 4
- 238000010586 diagram Methods 0.000 claims description 14
- 238000012546 transfer Methods 0.000 claims description 12
- 229910052704 radon Inorganic materials 0.000 claims description 6
- SYUHGPGVQRZVTB-UHFFFAOYSA-N radon atom Chemical compound [Rn] SYUHGPGVQRZVTB-UHFFFAOYSA-N 0.000 claims description 6
- 238000004364 calculation method Methods 0.000 claims description 4
- 238000001228 spectrum Methods 0.000 claims description 4
- 239000000306 component Substances 0.000 description 8
- 238000004088 simulation Methods 0.000 description 7
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000003745 diagnosis Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000012423 maintenance Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000005215 recombination Methods 0.000 description 2
- 230000006798 recombination Effects 0.000 description 2
- 239000008358 core component Substances 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000009659 non-destructive testing Methods 0.000 description 1
- 239000013307 optical fiber Substances 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 239000000523 sample Substances 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 238000013519 translation Methods 0.000 description 1
Images
Classifications
-
- F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
- F01—MACHINES OR ENGINES IN GENERAL; ENGINE PLANTS IN GENERAL; STEAM ENGINES
- F01D—NON-POSITIVE DISPLACEMENT MACHINES OR ENGINES, e.g. STEAM TURBINES
- F01D21/00—Shutting-down of machines or engines, e.g. in emergency; Regulating, controlling, or safety means not otherwise provided for
- F01D21/003—Arrangements for testing or measuring
-
- F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
- F01—MACHINES OR ENGINES IN GENERAL; ENGINE PLANTS IN GENERAL; STEAM ENGINES
- F01D—NON-POSITIVE DISPLACEMENT MACHINES OR ENGINES, e.g. STEAM TURBINES
- F01D25/00—Component parts, details, or accessories, not provided for in, or of interest apart from, other groups
-
- 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
- Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
- Y02E10/00—Energy generation through renewable energy sources
- Y02E10/70—Wind energy
- Y02E10/72—Wind turbines with rotation axis in wind direction
Landscapes
- Engineering & Computer Science (AREA)
- Mechanical Engineering (AREA)
- General Engineering & Computer Science (AREA)
- Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
本发明公开了一种叶端定时信号时频滤波方法，方法包括以下步骤：利用叶端定时传感器获取旋转叶片的时间脉冲，然后根据叶片半径R和转速n将实际达到时间和理论达到时间差Δt转换为位移数据；根据采样频率选择分析数据，使用频谱分析方法得到混叠频率幅值向量，将不同采样频率下的混叠频率幅值向量进行组合，形成二维幅值矩阵，由转频序列、混叠频率序列、二维幅值矩阵得到转频‑混叠频率图；识别RAF中的直线，得到直线的表达式，即斜率与截距参数；根据直线表达式计算不同采样频率下的混叠频率，根据该混叠频率产生带通滤波频带范围，使用该带通滤波器对信号进行滤波。本发明方法具有可解释性，滤波能力强的特点，可用于对叶端定时信号滤波。
Description
技术领域
本发明属于叶片无损检测领域，具体涉及一种叶端定时信号时频滤波方法。
背景技术
涡轮机械广泛应用于工业中，同时旋转叶片又是涡轮机械中核心部件。但由于长期承受气流冲击、高速旋转所引起的离心力，外物冲击，叶片成为涡轮机械中较为薄弱的部件，叶片的损伤将会导致设备失效甚至是安全事故，所以对叶片进行故障诊断十分有必要。传统的维修方式的定期维修，严重依赖于维修经验，且容易维修过剩和维修不及时。测量叶片的振动，从中分析叶片的健康状态称为了流行的方法，单传统的接触式测量方法，如应变片测量，耗时，寿命短且效率低，很难应用到实际的工业环境下，叶端定时(Blade TipTiming，BTT)是一种非接触式方法，它使用安装在静止机匣中的探头(电容式、光纤式、电涡流式等)，记录叶片的达到时间脉冲，将不考虑叶片振动情况下的理论达到时间和实际达到时间之差转换成叶端位移，并以此计算叶片末端的位移。叶端定时信号中存在着噪声，包括背景造成，转频及其倍频所产生的振动噪声，为了得到准确的故障诊断信息，有必要对叶端定时信号进行滤波，但叶端定时信号具有欠采样的特性，这也给滤波工作带来了困扰，传统的叶端定时滤波方式往往是通过多项式滤波、小波降噪等方法，效果十分有限，所以需要一种更先进的叶端定时滤波方法。
发明内容
鉴于以上内容，本发明的目的在于提出一种叶端定时信号时频滤波方法，对严重欠采样的叶端定时数据进行滤波。。
为实现上述发明目的，本发明采用以下技术方案：
一种叶端定时信号时频滤波方法包括以下步骤：
第一步骤，利用叶端定时传感器获取叶片的时间脉冲，基于所述时间脉冲生成转速n，然后根据叶片半径R和转速n将实际达到时间和理论达到时间差Δt转换为位移数据；
第二步骤，基于所述位移数据选择转频范围，基于频谱分析得到混叠频率幅值向量，将不同采样频率下的混叠频率幅值向量进行组合，形成二维幅值矩阵，基于转频频率序列、混叠频率序列、二维幅值矩阵得到转频-混叠频率图，其中，转频频率序列经由转频频率步长和转频范围确定，混叠频率序列经由混叠频率分辨率和最高转频确定；
第三步骤，识别转频-混叠频率图中的直线，得到直线的表达式fh＝afr+b，a，b为参数，
第四步骤，基于表达式计算不同采样频率下的混叠频率，根据所述混叠频率产生带通滤波频带范围，基于带通滤波频带范围对位移数据进行滑动滤波，所述位移数据包括位移时域信号。
所述的叶端定时信号时频滤波方法中，第一步骤中，
使用叶端定时传感器采集变转速情况下叶盘叶片的达到时间脉冲t，利用转速传感器获取转轴的转频fr，计算在无振动情况下，叶片达到第k个传感器的理论达到时间
其中，
表示第i个叶片在第p圈到达第k个传感器的理论时间，θi表示第i个叶片的角度，αk表示第k个传感器的安装角度，fr(p)表示第p圈时的转频，
表示第i个叶片在第p-1圈达到第k个传感器的实际时间。
根据叶片半径R和转频fr将实际达到时间与理论达到时间差转换为位移数据：
其中，
表示第i个叶片在
时刻的位移，其中
表示第i个叶片在第p圈到达第k个传感器的实际时间与理论时间之差，
得到了叶片末端的位移信号d，其中，d表示的是位移序列，及位移向量，包含了多个位移值。
所述的叶端定时信号时频滤波方法中，第二步骤中，位移数据的转频范围
基于旋转频率间隔Δfr生成转频序列Fr，确定分析数据截取的中心位置：
Fh＝[0，Δfh，2Δfh，3Δfh，…，(M-1)Δfh]1×M，其中N是整数，满足
由此来确定N的值，M是整数，满足
由此来确定M的值，
根据混叠频率分辨率Δfh和转速计算转频序列Fr中每一个转频对应的数据截取长度：
根据第i个理想的截取数据中心位置的转频频率Fr(i)计算实际最接近该转频的数据索引index：
其中，
表示求使得后式取到最小值的index值，
以index(i)为截取数据中心，以wlen(i)为数据长度，截取位移数据d[index(i)-round(wlen(i)/2)：index(i)+round(wlen(i)/2)]，选取位移数据中索引为index(i)-round(wlen(i)/2)到index(i)+round(wlen(i)/2)，其利用快速傅里叶变换进行频谱估计，
其中，d(n)为采样得到的信号，i是虚数符号，
Nfft是所选择的数据的长度；n是一个迭代数，从index(i)-round(wlen(i)/2)遍历到index(i)+round(wlen(i)/2)，取遍d中的所有元素，k是一个0到N-1整数，D(k)表示离散傅里叶变换后的第k个数据，将D取绝对值，即得到了混叠频率幅值向量AD，
生成零矩阵AN×M，其中N为转频序列长度，M为混叠频率序列长度，
将由转频序列Fr中第i个转频数据Fr(i)确定的数据索引范围index(i)-round(wlen(i)/2)到index(i)+round(wlen(i)/2)所截取的数据频谱幅值估计结果ADi填入矩阵AN×M中的第i列，以此类推，当选取转频序列中的所有转频元素，重复上述操作后，得到了完整的幅值矩阵AN×M，由转频序列Fr、混叠频率序列Fh、二维幅值矩阵AN×M得到转频-混叠频率图。
所述的叶端定时信号时频滤波方法中，第三步骤中，利用直线拟合或者Radon变换提取直线参数。
所述的叶端定时信号时频滤波方法中，第三步骤中，转频-混叠频率图的混叠频率方向进行遍历搜索，找到幅值最大的连续n个点，n为2-5，
其中I(k)为使得
取到最大值的索引位置，
所得到的点的坐标为：
(Fr(k)，Fh(I(k)))，(Fr(k)，Fh(I(k)))，…，(Fr(k)，Fh(I(k)))，
记录下坐标点的横坐标的采样频率和纵坐标的混叠频率，选定直线线段范围，对该线段范围内的坐标点进行最小二乘拟合，得到直线的斜率：
其中
d＝[Fh(index(1)，Fh(index(2)，…，Fh(index(np)]T，a是直线斜率，b是直线截距。
所述的叶端定时信号时频滤波方法中，第三步骤中，Radon变换中，将转频-混叠频率图上的值沿着各个直线方向投影，设定角度步长Δβ和距离步长ΔL以产生角度序列β和距离序列L：
β＝[-90°，-90°+Δθ，-90°+2Δθ，…，90°]，
其中H、W分别为RAF图像的宽，高像素尺寸，
在Radon变换后，转频-混叠频率图的峰值坐标(β，L)如下转换得到斜率与截距：
通过转频推算出混叠频率fh＝kfr+b。
所述的叶端定时信号时频滤波方法中，第四步骤中，设定滑动滤波时的窗长len，重叠率kp，形成带通滤波滑动窗口，第i个带通滤波滑动窗口的包含的索引为1+(i-1)·kp·len：(i-1)·kp·len+len，对数据d(1+(i-1)·kp·len)，…，d((i-1)·kp·len+len)进行带通滤波，将带通滤波后的数据首尾重组，得到滤波后的数据，其中带通滤波的滤波频带为[fh(i)-Δh1，fh(i)+Δh2]，使得频率在fh(i)-Δh1～fh(i)+Δh2范围内的频率通过，其中
本发明方法是一种叶端定时信号时频滤波方法，基于对转频-混叠频率图中的直线进行特征识别，由此根据转频产生响应的带通频带，然后使用带通滤波器对信号进行滑动滤波，本发明专利的滤波方法具有可解释性，滤波能力强的特点，可用于对叶端定时信号滤波。
附图说明
参照下述说明，结合附图，可以对本发明有最佳的理解。
图1为本专利提出的一种叶端定时信号时频滤波方法；
图2为单个叶端定时传感器原始仿真信号的时域图；
图3为单个叶端定时传感器加噪后仿真信号的时域图；
图4为仿真信号的转频-混叠频率(RAF)图；
图5为同一数据范围下原始信号、加噪信号、滤波信号的频域图；
图6为同一数据范围下原始信号、加噪信号、滤波信号的时域图。
具体实施方式
为更加清楚地表明本发明的目的、技术方案和优点，现结合附图1至图6及示例性实例，进一步详细地说明本发明。应当理解，此处所描述的示例性实例仅用以解释本发明，并不限定本发明的适用范围。
本发明提出一种叶端定时信号时频滤波方法，包括下述步骤：
(1)利用端定时传感器获取旋转叶片的时间脉冲，首先将其转换成转速，然后根据叶片半径R和转速n将实际达到时间和理论达到时间差Δt转换为位移数据。
在本示例性实例中，根据叶端定时的原理产生仿真信号，叶片的固有频率设定为875Hz，叶片末端振动位移包含875Hz的固有频率，另外还包含转频，转频的倍频，具体仿真参数设定如表1所示：
表1 叶端定时仿真信号参数说明
由此产生叶端位移，并通过单个叶端定时传感器采集转换成位移，其中转频及其倍频成分及添加的噪声都被视为噪声，将固有频率成分信号视为原始信号(不含噪)，将含转频及其倍频成分以及背景噪声的信号视为含噪信号，获取到的原始信号如图2所示，加噪信号如图3所示。
(2)根据采样频率选择分析数据，使用频谱分析方法得到混叠频率幅值向量，将不同采样频率下的混叠频率幅值向量进行组合，形成二维幅值矩阵，由转频序列、混叠频率序列、二维幅值矩阵得到转频-混叠频率(Rotating Aliasing Frequency，RAF)图。
在一个实施例中，转频频率序列和叠频率序列是等差序列。
在本示例性实例中，选择分析的数据的转频范围为91Hz到191Hz，采样频率步长Δfr为0.4Hz，混叠频率步长Δfh为0.5Hz，由此生成转频序列Fr和混叠频率步长序列Fh：
Fr＝[91，91.4，91.8，92.2，…，189]1×246
Fh＝[0，0.5，1，1.5，2，2.5，…，94.5]1×190
选择转频序列Fr中第1个转频Fr(1)＝91Hz，计算截取长度：
根据第1个理想的截取数据中心位置的转频频率Fr(1)＝91Hz，计算实际最接近该采样频率的数据索引index：
以92为截取数据中心，以182为数据长度，截取位移数据d(1)～d(182)。对该部分数据进行频谱分析，本实例中选用快速傅里叶变换：
其中d(n)为采样得到的信号，i是虚数符号，
n是一个迭代数，从1遍历到128，即取遍d中的所有元素，k是一个1到128整数，D(k)表示离散傅里叶变换后的第k个数据；D1＝[D(1)，D(2)，D(3)，…，D(182)]，对D1取绝对值，即得到了混叠频率幅值向量AD1。
选择转频序列Fr中第2个转频Fr(2)＝91.5Hz，重复上述操作，得到AD2，对转频序列中的所有元素进行遍历后，将最终得到的，混叠频率幅值向量，AD1～AD246横向拼接，即按列组合，形成一个二维矩阵，其中长度不足的地方用0填满，由此得到了二维幅值矩阵A190×246。
由转频序列Fr、混叠频率序列Fh、二维幅值矩阵A190×246即可得到转频-混叠频率(RAF)图，如图4所示。
(3)识别RAF图中的直线，得到直线的表达式fh＝afr+b，即a，b参数。
在本实例中，选用直线拟合的方法得到直线表达式，为了避免低频成分的干扰，将RAF图中混叠频率低于4Hz的部分舍去。对处理后的RAF图的混叠频率方向进行遍历搜索，找到幅值最大的连续3个点，选择A190×246的第1列，进行搜索：
注意这里是索引52是对于去掉混叠低频4Hz及其以下部分后的剩下部分的索引，即是新的混叠序列
山的索引为52的位置。
所得到的点的坐标为：
(91，30)，(91，30.5)，(91，31)
对A190×246中的所有列进行遍历后，共得到738个坐标，选定直线线段范围，对该范围内的坐标点进行最小二乘拟合，得到直线的斜率：
其中
a是直线斜率，b是直线截距。
对所有直线拟合，拟合结果如表2所示：
表2 直线拟合结果
| 参数 | 斜率a | 截距b |
| 直线1 | -8.86 | 875.88 |
| 直线2 | 8.86 | -874.58 |
| 直线3 | -8.37 | 872.01 |
| 直线4 | 7.98 | -869.02 |
| 直线5 | -6.99 | 872.35 |
| 直线6 | 7.01 | -874.85 |
| 直线7 | -6.00 | 874.80 |
| 直线8 | 5.94 | -866.05 |
| 直线9 | -4.99 | 872.75 |
| 直线10 | 4.99 | -872.75 |
(4)根据直线表达式计算不同采样频率下的混叠频率，根据该混叠频率产生带通滤波频带范围，使用该带通滤波器对信号进行滤波。
在本示例性实例中，设定滑动滤波所用的窗长len＝64，重叠率kp＝0.25，以第8条直线所在位置为例，展示滤波效果，第1个滑动窗口索引为：7000∶7063，计算混叠频率：
设定Δh1＝2，Δh2＝2，形成带通滤波频带为：16.21～20.21Hz，对信号进行带通滤波，将滤波后的时域数据首位重组，得到滤波后的时域数据，原始信号，加噪信号和滤波后信号相同索引位置下数据的频域图如图5所示，可以看到，噪声信号被大量过滤，仅保留了与混叠后的固有频率相关的频率分量。滤波后信号的时域图如图6所示，可以看到与原始信号的时域波形更为接近，对叶端定时信号的去噪将有利于后续的信号处理。
本发明专利方法是一种叶端定时信号时频滤波方法，基于对转频-混叠频率图中的直线进行特征识别，由此根据转频产生响应的带通频带，然后使用带通滤波器对信号进行滑动滤波，本发明专利的滤波方法具有可解释性，滤波能力强的特点，可用于对叶端定时信号滤波。
【应用实例】
根据叶端定时的原理产生仿真信号，叶片的固有频率设定为875Hz，叶片末端振动位移包含875Hz的固有频率，另外还包含转频，转频的倍频，具体仿真参数设定如表1所示：
表1 叶端定时仿真信号参数说明
由此产生叶端位移，并通过单个叶端定时传感器采集转换成位移，其中转频及其倍频成分及添加的噪声都被视为噪声，将固有频率成分信号视为原始信号(不含噪)，将含转频及其倍频成分以及背景噪声的信号视为含噪信号，获取到的原始信号如图2所示，加噪信号如图3所示。
选择分析的数据的转频范围为91Hz到191Hz，采样频率步长Δfr为0.4Hz，混叠频率步长Δfh为0.5Hz，由此生成转频序列Fr和混叠频率步长序列Fh：
Fr＝[91，91.4，91.8，92.2，…，189]1×246
Fh＝[0，0.5，1，1.5，2，2.5，…，94.5]1×190
选择转频序列Fr中第1个转频Fr(1)＝91Hz，计算截取长度：
根据第1个理想的截取数据中心位置的转频频率Fr(1)＝91Hz，计算实际最接近该采样频率的数据索引index：
以92为截取数据中心，以182为数据长度，截取位移数据d(1)～d(182)。对该部分数据进行频谱分析，本实例中选用快速傅里叶变换：
其中d(n)为采样得到的信号，i是虚数符号，
n是一个迭代数，从1遍历到128，即取遍d中的所有元素，k是一个1到128整数，D(k)表示离散傅里叶变换后的第k个数据；D1＝[D(1)，D(2)，D(3)，…，D(182)]，对D1取绝对值，即得到了混叠频率幅值向量AD1。
选择转频序列Fr中第2个转频Fr(2)＝91.5Hz，重复上述操作，得到AD2，对转频序列中的所有元素进行遍历后，将最终得到的，混叠频率幅值向量，AD1～AD246横向拼接，即按列组合，形成一个二维矩阵，其中长度不足的地方用0填满，由此得到了二维幅值矩阵A190×246。
由转频序列Fr、混叠频率序列Fh、二维幅值矩阵A190×246即可得到转频-混叠频率(RAF)图，如图4所示。
选用直线拟合的方法得到直线表达式，为了避免低频成分的干扰，将RAF图中混叠频率低于4Hz的部分舍去。对处理后的RAF图的混叠频率方向进行遍历搜索，找到幅值最大的连续3个点，选择A190×246的第1列，进行搜索：
注意这里是索引52是对于去掉混叠低频4Hz及其以下部分后的剩下部分的索引，即是新的混叠序列
中的索引为52的位置。
所得到的点的坐标为：
(91，30)，(91，30.5)，(91，31)
对A190×246中的所有列进行遍历后，共得到738个坐标，选定直线线段范围，对该范围内的坐标点进行最小二乘拟合，得到直线的斜率：
其中
a是直线斜率，b是直线截距。
对所有直线拟合，拟合结果如表2所示：
表2 直线拟合结果
| 参数 | 斜率a | 截距b |
| 直线1 | -8.86 | 875.88 |
| 直线2 | 8.86 | -874.58 |
| 直线3 | -8.37 | 872.01 |
| 直线4 | 7.98 | -869.02 |
| 直线5 | -6.99 | 872.35 |
| 直线6 | 7.01 | -874.85 |
| 直线7 | -6.00 | 874.80 |
| 直线8 | 5.94 | -866.05 |
| 直线9 | -4.99 | 872.75 |
| 直线10 | 4.99 | -872.75 |
设定滑动滤波所用的窗长len＝64，重叠率kp＝0.25，以第8条直线所在位置为例，展示滤波效果，第1个滑动窗口索引为：7000∶7063，计算混叠频率：
设定Δh1＝2，Δh2＝2，形成带通滤波频带为：16.21～20.21Hz，对信号进行带通滤波，将滤波后的时域数据首位重组，得到滤波后的时域数据，原始信号，加噪信号和滤波后信号相同索引位置下数据的频域图如图5所示，可以看到，噪声信号被大量过滤，仅保留了与混叠后的固有频率相关的频率分量。滤波后信号的时域图如图6所示，可以看到与原始信号的时域波形更为接近，对叶端定时信号的去噪将有利于后续的信号处理。
本发明专利方法是一种叶端定时信号时频滤波方法，基于对转频-混叠频率图中的直线进行特征识别，由此根据转频产生响应的带通频带，然后使用带通滤波器对信号进行滑动滤波，本发明专利的滤波方法具有可解释性，滤波能力强的特点，可用于对叶端定时信号滤波。
以上所述仅为本发明的较佳实例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
Claims (7)
1.一种叶端定时信号时频滤波方法，其特征在于，其包括以下步骤：
第一步骤，利用叶端定时传感器获取叶片的时间脉冲，基于所述时间脉冲生成转速n，然后根据叶片半径R和转速n将实际达到时间和理论达到时间差Δt转换为位移数据；
第二步骤，基于所述位移数据选择转频范围，基于频谱分析得到混叠频率幅值向量，将不同采样频率下的混叠频率幅值向量进行组合，形成二维幅值矩阵，基于转频频率序列、混叠频率序列、二维幅值矩阵得到转频-混叠频率图，其中，转频频率序列经由转频频率步长和转频范围确定，混叠频率序列经由混叠频率分辨率和最高转频确定；
第三步骤，识别转频-混叠频率图中的直线，得到直线的表达式fh＝afr+b，a，b为参数，
第四步骤，基于表达式计算不同采样频率下的混叠频率，根据所述混叠频率产生带通滤波频带范围，基于带通滤波频带范围对位移数据进行滑动滤波。
2.根据权利要求1所述的叶端定时信号时频滤波方法，其中，优选的，第一步骤中，
使用叶端定时传感器采集变转速情况下叶盘叶片的达到时间脉冲t，利用转速传感器获取转轴的转频fr，计算在无振动情况下，叶片达到第k个传感器的理论达到时间
其中，
表示第i个叶片在第p圈到达第k个传感器的理论时间，θi表示第i个叶片的角度，αk表示第k个传感器的安装角度，fr(p)表示第p圈时的转频，
表示第i个叶片在第p-1圈达到第k个传感器的实际时间，
根据叶片半径R和转频fr将实际达到时间与理论达到时间差转换为位移数据：
其中，
表示第i个叶片在
时刻的位移，其中
表示第i个叶片在第p圈到达第k个传感器的实际时间与理论时间之差，
得到了叶片末端的位移信号d，其为位移向量。
3.根据权利要求1所述的叶端定时信号时频滤波方法，其中，第二步骤中，位移数据的转频范围
基于旋转频率间隔Δfr生成转频序列Fr，确定分析数据截取的中心位置：
Fh＝[0，Δfh，2Δfh，3Δfh，…，(M-1)Δfh]1×M，其中N是整数，满足
由此来确定N的值，M是整数，满足
由此来确定M的值，
根据混叠频率分辨率Δfh和转速计算转频序列Fr中每一个转频对应的数据截取长度：
根据第i个理想的截取数据中心位置的转频频率Fr(i)计算实际最接近该转频的数据索引index：
其中，
表示求使得后式取到最小值的index值，
以index(i)为截取数据中心，以wlen(i)为数据长度，截取位移数据d[index(i)-round(wlen(i)/2)：index(i)+round(wlen(i)/2)]，选取位移数据中索引为index(i)-round(wlen(i)/2)到index(i)+round(wlen(i)/2)，其利用快速傅里叶变换进行频谱估计，
其中，d(n)为采样得到的信号，i是虚数符号，
Nfft是所选择的数据的长度；n是一个迭代数，从index(i)-round(wlen(i)/2)遍历到index(i)+round(wlen(i)/2)，取遍d中的所有元素，k是一个0到N-1整数，D(k)表示离散傅里叶变换后的第k个数据，将D取绝对值，即得到了混叠频率幅值向量AD，
生成零矩阵AN×M，其中N为转频序列长度，M为混叠频率序列长度，
将由转频序列Fr中第i个转频数据Fr(i)确定的数据索引范围index(i)-round(wlen(i)/2)到index(i)+round(wlen(i)/2)所截取的数据频谱幅值估计结果ADi填入矩阵AN×M中的第i列，以此类推，当选取转频序列中的所有转频元素，重复上述操作后，得到了完整的幅值矩阵AN×M，由转频序列Fr、混叠频率序列Fh、二维幅值矩阵AN×M得到转频-混叠频率图。
4.根据权利要求1所述的叶端定时信号时频滤波方法，其中，第三步骤中，利用直线拟合或者Radon变换提取直线参数。
5.根据权利要求4所述的叶端定时信号时频滤波方法，其中，第三步骤中，转频-混叠频率图的混叠频率方向进行遍历搜索，找到幅值最大的连续n个点，n为2-5：
所得到的点的坐标为：
(Fr(k)，Fh(I(k)))，(Fr(k)，Fh(I(k)))，…，(Fr(k)，Fh(I(k)))，
其中I(k)为使得
取到最大值的索引位置，
记录下坐标点的横坐标的采样频率和纵坐标的混叠频率，选定直线线段范围，对该线段范围内的坐标点进行最小二乘拟合，得到直线的斜率：
其中
d＝[Fh(index(1)，Fh(index(2)，…，Fh(index(np)]T，a是直线斜率，b是直线截距。
6.根据权利要求4所述的叶端定时信号时频滤波方法，其中，第三步骤中，Radon变换中，将转频-混叠频率图上的值沿着各个直线方向投影，设定角度步长Δβ和距离步长ΔL以产生角度序列β和距离序列L：
β＝[-90°，-90°+Δθ，-90°+2Δθ，…，90°]，
其中H、W分别为RAF图像的宽，高像素尺寸，
在Radon变换后，转频-混叠频率图的峰值坐标(β，L)如下转换得到斜率与截距：
通过转频推算出混叠频率fh＝afr+b。
7.根据权利要求1所述的叶端定时信号时频滤波方法，其中，第四步骤中，设定滑动滤波时的窗长len，重叠率kp，形成带通滤波滑动窗口，第i个带通滤波滑动窗口的包含的索引为1+(i-1)·kp·len：(i-1)·kp·len+len，对数据d(1+(i-1)·kp·len)，…，d((i-1)·kp·len+len)进行带通滤波，将带通滤波后的数据首尾重组，得到滤波后的数据，其中带通滤波的滤波频带为[fh(i)-Δh1，fh(i)+Δh2]，使得频率在fh(i)-Δh1～fh(i)+Δh2范围内的频率通过，其中
Δh1和Δh2分为带通滤波上、下边界阈值，fr(i)表示第i个带通滤波滑动窗口中心数据对应的转频，fh(i)表示由表达式
计算得到的混叠频率，其中a，b为转频-混叠频率图中直线的斜率和截距。
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110911788 | 2021-08-10 | ||
| CN2021109117881 | 2021-08-10 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113565584A true CN113565584A (zh) | 2021-10-29 |
| CN113565584B CN113565584B (zh) | 2022-08-09 |
Family
ID=78173454
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111017893.7A Active CN113565584B (zh) | 2021-08-10 | 2021-08-31 | 一种叶端定时信号时频滤波方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113565584B (zh) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114858270A (zh) * | 2022-05-06 | 2022-08-05 | 南京方一测控科技有限公司 | 一种用于核电汽轮机末级叶片低频振动检测系统及方法 |
| CN115114740A (zh) * | 2022-06-02 | 2022-09-27 | 西安交通大学 | 基于脉冲序列生成的非接触式测量的校核方法与系统 |
| CN115586260A (zh) * | 2022-09-29 | 2023-01-10 | 西安交通大学 | 基于连续压缩感知的叶端定时信号无网格频谱估计方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105043386A (zh) * | 2015-06-11 | 2015-11-11 | 北京航空航天大学 | 光纤陀螺类盲发变滤波滑窗长度的异步通信数据传输方法 |
| WO2015196735A1 (zh) * | 2014-06-23 | 2015-12-30 | 华南理工大学 | 基于啮合频率和频谱校正技术的风电齿轮箱阶次跟踪方法 |
| CN109682601A (zh) * | 2019-03-04 | 2019-04-26 | 北京天泽智云科技有限公司 | 一种变转速工况下滚动轴承的早期故障识别方法 |
| CN111879508A (zh) * | 2020-07-28 | 2020-11-03 | 无锡迈斯德智能测控技术有限公司 | 基于时频变换的旋转机械瞬时转速估计方法、装置和存储介质 |
| WO2021063294A1 (zh) * | 2019-09-30 | 2021-04-08 | 华能四川水电有限公司 | 一种转子的中心偏移量检测方法、装置、存储介质和设备 |
-
2021
- 2021-08-31 CN CN202111017893.7A patent/CN113565584B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2015196735A1 (zh) * | 2014-06-23 | 2015-12-30 | 华南理工大学 | 基于啮合频率和频谱校正技术的风电齿轮箱阶次跟踪方法 |
| CN105043386A (zh) * | 2015-06-11 | 2015-11-11 | 北京航空航天大学 | 光纤陀螺类盲发变滤波滑窗长度的异步通信数据传输方法 |
| CN109682601A (zh) * | 2019-03-04 | 2019-04-26 | 北京天泽智云科技有限公司 | 一种变转速工况下滚动轴承的早期故障识别方法 |
| WO2021063294A1 (zh) * | 2019-09-30 | 2021-04-08 | 华能四川水电有限公司 | 一种转子的中心偏移量检测方法、装置、存储介质和设备 |
| CN111879508A (zh) * | 2020-07-28 | 2020-11-03 | 无锡迈斯德智能测控技术有限公司 | 基于时频变换的旋转机械瞬时转速估计方法、装置和存储介质 |
Non-Patent Citations (2)
| Title |
|---|
| 张西宁等: "分数阶域滤波在启停车过程转频振动分量提取中的应用", 《西安交通大学学报》 * |
| 滕伟等: "基于时频滤波的汽轮机半速涡动故障成分提取", 《振动与冲击》 * |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114858270A (zh) * | 2022-05-06 | 2022-08-05 | 南京方一测控科技有限公司 | 一种用于核电汽轮机末级叶片低频振动检测系统及方法 |
| CN115114740A (zh) * | 2022-06-02 | 2022-09-27 | 西安交通大学 | 基于脉冲序列生成的非接触式测量的校核方法与系统 |
| CN115114740B (zh) * | 2022-06-02 | 2024-03-19 | 西安交通大学 | 基于脉冲序列生成的非接触式测量的校核方法与系统 |
| CN115586260A (zh) * | 2022-09-29 | 2023-01-10 | 西安交通大学 | 基于连续压缩感知的叶端定时信号无网格频谱估计方法 |
| CN115586260B (zh) * | 2022-09-29 | 2024-05-28 | 西安交通大学 | 基于连续压缩感知的叶端定时信号无网格频谱估计方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113565584B (zh) | 2022-08-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113565584B (zh) | 一种叶端定时信号时频滤波方法 | |
| CN110567574B (zh) | 一种旋转叶片叶端定时振动参数辨识方法与系统 | |
| Wang et al. | Rolling element bearing fault diagnosis via fault characteristic order (FCO) analysis | |
| CN111089726B (zh) | 一种基于最优维数奇异谱分解的滚动轴承故障诊断方法 | |
| CN110763462B (zh) | 一种基于同步压缩算子的时变振动信号故障诊断方法 | |
| Wang et al. | Bearing fault diagnosis under time-varying rotational speed via the fault characteristic order (FCO) index based demodulation and the stepwise resampling in the fault phase angle (FPA) domain | |
| CN113702044B (zh) | 一种轴承故障检测方法及系统 | |
| CN110046476B (zh) | 滚动轴承故障的三元二进分形小波稀疏诊断方法 | |
| CN104215456B (zh) | 一种基于平面聚类和频域压缩感知重构的机械故障诊断方法 | |
| Lin et al. | A review and strategy for the diagnosis of speed-varying machinery | |
| CN111879508B (zh) | 基于时频变换的旋转机械瞬时转速估计方法、装置和存储介质 | |
| CN111256993A (zh) | 一种风电机组主轴承故障类型诊断方法及系统 | |
| CN105628176A (zh) | 一种旋转机械扭转振动信号采集分析方法 | |
| CN116361733A (zh) | 一种故障诊断方法、装置、系统以及存储介质 | |
| CN113586177B (zh) | 基于单叶端定时传感器的叶片固有频率识别方法 | |
| CN108317052B (zh) | 齿轮的损伤因子的检测方法及装置、风力发电机组 | |
| CN111323233B (zh) | 一种用于低速旋转机械故障诊断的局部均值分解方法 | |
| Sabbatini et al. | Data acquisition and processing for tip timing and operational modal analysis of turbomachinery blades | |
| CN112597969A (zh) | 一种滚动轴承故障诊断方法、系统及介质 | |
| CN109916624B (zh) | 一种基于希尔伯特黄的滚珠丝杠副疲劳失效诊断方法 | |
| CN114486252B (zh) | 一种矢量模极大值包络的滚动轴承故障诊断方法 | |
| CN114061746B (zh) | 旋转机械故障诊断中的重复瞬变信号提取方法 | |
| CN113504310B (zh) | 基于单个叶端定时传感器的叶片固有频率识别方法 | |
| CN110779723B (zh) | 一种基于霍尔信号的变速工况电机轴承精确故障诊断方法 | |
| CN113565585B (zh) | 单叶端定时传感器变工况旋转叶片固有频率提取方法 |
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 |