CN110991055A - 一种旋转类机械产品剩余寿命预测系统 - Google Patents

一种旋转类机械产品剩余寿命预测系统 Download PDF

Info

Publication number
CN110991055A
CN110991055A CN201911246253.6A CN201911246253A CN110991055A CN 110991055 A CN110991055 A CN 110991055A CN 201911246253 A CN201911246253 A CN 201911246253A CN 110991055 A CN110991055 A CN 110991055A
Authority
CN
China
Prior art keywords
module
data
unit
signal
analog 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.)
Granted
Application number
CN201911246253.6A
Other languages
English (en)
Other versions
CN110991055B (zh
Inventor
刘浩阔
边智
孟理华
杨俊斌
刘漪纹
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Aero Polytechnology Establishment
Original Assignee
China Aero Polytechnology Establishment
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China Aero Polytechnology Establishment filed Critical China Aero Polytechnology Establishment
Priority to CN201911246253.6A priority Critical patent/CN110991055B/zh
Publication of CN110991055A publication Critical patent/CN110991055A/zh
Application granted granted Critical
Publication of CN110991055B publication Critical patent/CN110991055B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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
    • Y02P—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN THE PRODUCTION OR PROCESSING OF GOODS
    • Y02P90/00—Enabling technologies with a potential contribution to greenhouse gas [GHG] emissions mitigation
    • Y02P90/30—Computing systems specially adapted for manufacturing

Landscapes

  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)

Abstract

本发明提供一种旋转类机械产品剩余寿命预测系统，其包括数据载入模块、模拟信号生成模块、数据预处理模块、频域分析模块、时频域分析模块、特征提取模块以及寿命预测模块，数据载入模块用于载入TDMS文件和mat文件；模拟信号生成模块用于生成正弦波、三角波、方波、锯齿波和白噪声；数据预处理模块包括去趋势单元、加窗单元、滤波单元以及奇异值分解单元，去趋势单元可以对模拟信号进行去趋势处理，消除信号中的直流分量。本发明的系统覆盖了信号预处理、频域分析、时频域分析、特征提取以及寿命预测的整个数据处理流程，方法具有通用性，可以在旋转类机械产品信号分析中使用。

Description

一种旋转类机械产品剩余寿命预测系统
技术领域
本发明涉及机械装备可靠性与寿命预测领域，具体地涉及一种旋转类机械产品剩余寿命预测系统。
背景技术
随着机械产品的可靠性要求越来越高，针对旋转类机械产品的状态监测与健康管理也变得越来越重要。在传统的机械产品信号分析、故障诊断与预测手段中，往往采用底层编程语言逐项进行分析，不仅工作量大，准确性也无法保证。针对旋转类机械产品来说，对采集的数据的分析处理、故障特征的特征提取、退化模型建立与预测等具体应用存在诸多通用方法。
现有技术中，尚没有类似的旋转类机械产品剩余寿命预测系统，在进行寿命预测时容易出现预测不准确的问题。
发明内容
针对现有方案的不足，本发明提供一种旋转类机械产品剩余寿命预测系统，实现了旋转类机械产品的数据载入、模拟信号生成、数据预处理、频域分析、时频域分析、特征提取、寿命预测等功能。
具体地，本发明提供一种旋转类机械产品剩余寿命预测系统，其包括数据载入模块、模拟信号生成模块、数据预处理模块、频域分析模块、时频域分析模块、特征提取模块以及寿命预测模块，数据载入模块、模拟信号生成模块、数据预处理模块、频域分析模块、时频域分析模块、特征提取模块以及寿命预测模块相互之间通讯连接；
所述数据载入模块用于载入TDMS文件和mat文件格式的产品数据信息；
所述模拟信号生成模块用于根据产品数据信息生成正弦波、三角波、方波、锯齿波或白噪声的模拟信号；
数据预处理模块对模拟信号进行处理，所述数据预处理模块包括去趋势单元、加窗单元、滤波单元以及奇异值分解单元，所述去趋势单元对模拟信号进行去趋势处理，消除信号中的直流分量；
所述加窗单元用于给模拟信号增加窗口，所述窗口的类型包括矩形窗、离散汉宁窗、离散哈明窗和布莱克曼窗，所述加窗单元能够显示原始波形、加窗后波形、加窗后波形时域图和加窗后波形频域图，模拟信号由加窗单元导出之后形成加窗后波形，
其中长度为N的离散化矩形窗时域表达式为：
Figure BDA0002307761460000021
长度为N的离散汉宁窗时域表达式为：
Figure BDA0002307761460000022
长度为N的离散哈明窗时域表达式为：
Figure BDA0002307761460000023
长度为N的布莱克曼窗时域表达式为：
Figure BDA0002307761460000024
所述滤波单元用于对模拟信号进行滤波，所述滤波单元能够显示模拟信号的原始时域信号、原始信号幅频谱、原始信号相频谱、滤波后时域信号、滤波后信号幅频谱和滤波后信号相频谱，所述滤波模块中能够设置滤波器类型、纹波、衰减、阶次、下限截止频率和上限截止频率，所述奇异值分解单元用于对模拟信号进行奇异值分解和重构；
所述奇异值分解单元能够对模拟信号进行奇异值分解和重构，对于长度小于或等于10000的数据进行奇异值分解，分解成功后绘制奇异值差分结果，输入重构阶次后，对分解后的信号进行重构，从而达到主成分提取和滤波的目的，奇异值分解具体为：
假设M为m*n阶矩阵，其中的矩阵元素均属于实数域或复数域，则存在一个分解使得
M＝U∑V*
其中，U是m*m阶酉矩阵，∑是m*n阶非负实数对角矩阵；而V*即V的共轭转置，是n*n阶酉矩阵；此分解称作M的奇异值分解；
所述频域分析模块能够自动绘制模拟信号的频谱图、周期图法功率谱和间接法功率谱，在细化谱设置模块中输入频率上限、频率下限和快速傅里叶变换点数，绘制模拟信号的细化谱，从而进行频域分析；
所述时频域分析模块包括经验模态分解单元和希尔伯特-黄变换单元，所述时频域分析模块能够对模拟信号进行重构，得到模拟信号各阶次的瞬时频率图，进行时频域分析，在经验模态分解单元中，使用经验模态分解方法将自动对模拟信号进行分解，绘制采集的数据及其内涵模态分量的时域图，以及采集的数据及其内涵模态分量的频谱图，设置重构信号的起始内涵模态分量和终止内涵模态分量，并对起始内涵模态分量和终止内涵模态分量中间的所有分量进行重构；
经验模态分解的具体流程如下：
步骤一：假设原始模拟信号为s(t)，找出原始模拟信号s(t)中的所有局部极大值以及局部极小值，接着利用三次样条，分别将局部极大值串连成上包络线，将局部极小值串连成下包络线；
步骤二：求出上包络线和下包络线的平均值，得到均值包络线m1(t)；
步骤三：将原始模拟信号s(t)与均值包络线m1(t)相减，得到第一个分量h1(t)；
h1(t)＝s(t)-m1(t)
步骤四：检查h1(t)是否符合IMF的条件，如果不符合，则回到步骤1并且将h1(t)当作原始模拟信号，进行第二次的筛选，亦即
h2(t)＝h1(t)-m2(t)
重复筛选k次，
hk(t)＝hk-1(t)-mk(t)
直到hk(t)符合IMF的条件，得到第一个IMF分量c1(t)，亦即
c1(t)＝hk(t)
步骤五：原始讯号s(t)减去c1(t)得到剩余量r1(t)，表示如下式：
r1(t)＝s(t)-c1(t)
步骤六：将r1(t)当做新的原始模拟信号，重复执行步骤一直步骤五，得到新的剩余量r2(t)，如此重复n次
r2(t)＝r1(t)-c2(t)
r3(t)＝r2(t)-c3(t)
…
rn(t)＝rn-1(t)-cn(t)
当第n个剩余量rn(t)已成为单调函数，将无法再分解IMF时，整个EMD的分解过程完成，原始讯号s(t)可以表示成n个IMF分量与一个平均趋势分量rn(t)的组合，亦即
Figure BDA0002307761460000041
所述希尔伯特-黄变换单元使用希尔伯特-黄变换方法，得到模拟信号各阶次的瞬时频率图，
希尔伯特-黄变换是在EMD分解的基础上，将IMF进行希尔伯特变换得到的，函数u(t)的希尔伯特变换公式是：
Figure BDA0002307761460000051
其中，x(t)为连续的时间信号，y(t)为经过希尔伯特变换后的信号；
所述特征提取模块包括线性判别分析单元和盲源分离单元，所述特征提取模块用于对模拟信号进行特征提取，获取模拟信号的故障征向量，并得到故障诊断结果，给出置信度，同时获得模拟信号的源信号；
线性判别分析单元利用线性判别分析提取模拟信号数据中的主成分，通过对监督数据的训练，形成故障特征向量，实现对未知试验数据的故障模式识别，线性判别分析单元通过载入不同故障模式的数据，计算特征评价指标值，选择指标值最大的特征值，在VMD分解后训练数据，导入待分类数据，得到故障诊断结果，并给出置信度；盲源分离单元利用盲源分离在没有任何先验知识的前提下，将源信号从混叠信号中提取和分离出来；
所述寿命预测模块使用自回归滑动平均模型进行退化模型的建立和预测，寿命预测模块通过载入时间序列数据、模型参数辨识、参数估计和模型检验步骤实现建模过程，并在模型建立完成后对退化特征的趋势进行预测，得到给定失效阈值下的机械产品寿命评估值。
优选地，线性判别分析的具体方法为：
假设C个类中每一个类都有均值μi和相同的协方差∑，那么，类间的变化通过类均值的协方差来定义：
Figure BDA0002307761460000052
这里μ是各类均值的均值，在
Figure BDA0002307761460000053
的方向区分类由下式给出：
Figure BDA0002307761460000054
如果
Figure BDA0002307761460000061
是∑-1∑b的特征向量，等同于用对应的特征值进行分类；
如果∑-1∑b是可对角化矩阵，特征之间的变化性就会被保留在C-1个最大特征值对应的特征向量构成子空间内；
优选地，盲源分离的具体方法为：
将盲源分离算法做出如下假定：具有m个独立的信号源s1(t),…,sm(t)和n个独立的观察量x1(t),…,xm(t)，观察量和信号源具有如下的关系：
x(t)＝As(t)
其中，x(t)＝[x1(t)…,xm(t)]T，s(t)＝[s1(t)…,sm(t)]T，A是一个n*m的系数矩阵，原问题变成了已知x(t)和s(t)的独立性，求对s(t)的估计问题，假设有如下公式：
y(t)＝Wx(t)
其中y(t)是对s(t)的估计，W是一个m*n系数矩阵；
优选地，所述寿命预测模块建模的具体方法如下：
自回归滑动平均模型中包含了p个自回归项和q个移动平均项，自回归滑动平均模型表示为：
Figure BDA0002307761460000062
其中εt是白噪声序列，p和q都是非负整数，c为常数，φi和θi为待定系数；
使用贝叶斯信息量准则对模型参数进行辨识，贝叶斯信息量定义为：
BIC＝ln(n)k-2ln(L)
其中，k为模型参数个数，n为样本数量，L为似然函数；
使用极大似然估计法对自回归滑动平均模型中的参数进行估计，假设x1,x2,…,xn为独立同分布的采样，θ为模型参数,f为我们所使用的模型，遵循上述的独立同分布假设，参数为θ的模型f产生上述采样表示为：
f(x1,x2,…,xn|θ)＝f(x1|θ)×f(x2|θ)…,f(xn|θ)
似然定义为：
Figure BDA0002307761460000071
其中lnL(θ|x1,x2,…,xn)为对数似然，极大似然函数表示为：
Figure BDA0002307761460000072
模型检验对模型估计的残差值进行检验，使用自相关和偏自相关检验残差的自相关和偏自相关特性，
其中自相关公式为：
Figure BDA0002307761460000073
偏自相关公式为：
Figure BDA0002307761460000074
优选地，数据载入模块在TDMS文件载入单元中读取tdms数据后，能够显示数据中的通道组以及通道信息，在点击不同通道后能够显示该通道的属性值、图和数据表，在mat文件载入单元中读取mat数据后，能够显示数据中包含的所有数据变量名，在点击数据变量名后能够显示该数据波形图。
优选地，所述模拟信号生成模块能够根据需要设置正弦波、三角波、方波、锯齿波和白噪声的幅值、频率、数据长度、偏置和采样率，在波形叠加模式中，并能够将基础波形进行相加、相乘运算，进而模拟生成各种复杂形式的仿真信号。
优选地，所述滤波器类型包括低通滤波器、高通滤波器、带通滤波器以及带阻滤波器。
优选地，所述滤波单元的滤波器设计方法为Butterworth，Butterworth滤波器幅频响应关系如下:
Figure BDA0002307761460000081
式中，Ω、Ωc、N分别表示频率、转折频率、系统阶数。
本发明的效果如下：
1.本发明的系统覆盖了信号预处理、频域分析、时频域分析、特征提取以及寿命预测的整个数据处理流程，方法具有通用性，可以在旋转类机械产品信号分析中使用。
2.本发明相较于其他信号分析手段，该软件系统更具有针对性，系统集成了信号分析与预测最为常用的分析方法，并包含专用工具箱，降低了学习成本，提高了分析效率。
3.本发明支持不同格式的数据文件，可以实现与MATLAB、Labview等软件的联合分析，灵活性更强。
附图说明
图1是本发明的结构示意框图；
图2是本发明的各个模块的结构示意框图；
图3是本发明的实施例中的截段后的采集的数据；
图4是本发明的实施例中的时域特征变化；
图5是本发明的实施例中的残差检验结果；
图6是本发明的实施例中的寿命预测结构。
具体实施方式
以下，参照附图对本发明的实施方式进行说明。
具体地，本发明提供一种旋转类机械产品剩余寿命预测系统，实现了旋转类机械产品的数据载入、数据预处理、频域分析、时频域分析、特征提取、寿命预测等功能。
本发明提供一种旋转类机械产品剩余寿命预测系统，如图1及图2所示，其包括数据载入模块1、模拟信号生成模块2、数据预处理模块3、频域分析模块4、时频域分析模块5、特征提取模块6以及寿命预测模块7。
所述数据载入模块用于载入TDMS文件和mat文件。数据载入模块在TDMS文件载入单元中读取tdms数据后，显示数据中的通道组以及通道信息，在点击不同通道后能够显示该通道的属性值、图和数据表，在mat文件载入单元中读取mat数据后，能够显示数据中包含的所有数据变量名，在点击数据变量名后能够显示该数据波形图。
所述模拟信号生成模块用于生成正弦波、三角波、方波、锯齿波和白噪声。所述模拟信号生成模块能够根据需要设置正弦波、三角波、方波、锯齿波和白噪声的幅值、频率、数据长度、偏置和采样率，在波形叠加模式中，能够将基础波形进行相加、相乘等运算，进而模拟生成各种复杂形式的仿真信号。
数据预处理模块包括去趋势单元、加窗单元、滤波单元以及奇异值分解单元，所述去趋势单元可以对模拟信号进行去趋势处理，消除信号中的直流分量。
所述加窗单元用于给模拟信号增加窗口，所述窗口的类型包括矩形窗、汉宁窗、哈明窗和布莱克曼窗，所述加窗单元能够显示原始波形、加窗后波形、加窗后波形时域图和加窗后波形频域图，模拟信号由加窗单元导出之后将被替换为加窗后波形，
其中长度为N的离散化矩形窗为：
Figure BDA0002307761460000101
长度为N的离散汉宁窗时域表达式为：
Figure BDA0002307761460000102
长度为N的离散哈明窗时域表达式为：
Figure BDA0002307761460000103
长度为N的布莱克曼窗时域表达式为：
Figure BDA0002307761460000104
所述滤波单元能够实现针对模拟信号的滤波，所述滤波单元能够显示原始时域信号、原始信号幅频谱、原始信号相频谱、滤波后时域信号、滤波后信号幅频谱和滤波后信号相频谱，所述滤波模块能够设置滤波器类型、滤波器设计方法、纹波、衰减、阶次、下限截止频率和上限截止频率，所述奇异值分解单元用于对模拟信号进行奇异值分解和重构。
优选地，所述滤波器类型包括低通滤波器、高通滤波器、带通滤波器以及带阻滤波器。
优选地，所述滤波单元的滤波器设计方法为Butterworth，Butterworth滤波器幅频响应关系如下:
Figure BDA0002307761460000105
式中，Ω、Ωc、N分别表示频率、转折频率、系统阶数。
所述奇异值分解单元能够对模拟信号进行奇异值分解和重构，对于长度小于或等于10000的数据进行奇异值分解，分解成功后绘制奇异值差分结果，输入重构阶次后，对分解后的信号进行重构，从而达到主成分提取和滤波的目的，奇异值分解具体为：
假设M为m*n阶矩阵，其中的矩阵元素均属于实数域或复数域，则存在一个分解使得
M＝U∑V*
其中，U是m*m阶酉矩阵，∑是m*n阶非负实数对角矩阵；而V*，即V的共轭转置，是n*n阶酉矩阵；此分解称作M的奇异值分解。
奇异值阈值的选取决定了奇异值分解降噪的效果，阈值数值过大会丢失重要的冲击信号，阈值数值过小则会使得降噪效果不明显。因此选择合适的阈值，在噪声降低的情况下不丢失冲击特征，是准确提取出冲击特征的关键之一，本专利利用奇异值差分谱的奇异值阈值确定方法来寻找合适的奇异值阈值。
所述频域分析模块，能够自动绘制模拟信号的频谱图、周期图法功率谱和间接法功率谱，在细化谱设置模块中输入频率上限、频率下限和NFFT(快速傅里叶变换点数)，绘制模拟信号的细化谱。
所述时频域分析模块包括经验模态分解单元和希尔伯特-黄变换单元，在经验模态分解单元中，使用EMD(经验模态)分解方法将自动对模拟信号进行分解，绘制采集的数据及其内涵模态分量的时域图，以及采集的数据及其内涵模态分量的频谱图，设置重构信号的起始IMF(内涵模态分量)和终止IMF，并对起始IMF和终止IMF中间的所有分量进行重构。
经验模态分解流程如下：
步骤一：假设原始信号为s(t)，找出s(t)中的所有局部极大值以及局部极小值，接着利用三次样条(cubic spline)，分别将局部极大值串连成上包络线与局部极小值串连成下包络线；
步骤二：求出上下包络线之平均，得到均值包络线m1(t)；
步骤三：原始信号s(t)与均值包络线相减，得到第一个分量h1(t)；
h1(t)＝s(t)-m1(t)
步骤四：检查h1(t)是否符合IMF的条件，如果不符合，则回到步骤1并且将h1(t)当作原始讯号，进行第二次的筛选，亦即
h2(t)＝h1(t)-m2(t)
重复筛选k次，
hk(t)＝hk-1(t)-mk(t)
直到hk(t)符合IMF的条件，即得到第一个IMF分量c1(t)，亦即
c1(t)＝hk(t)
步骤五：原始讯号s(t)减去c1(t)可得到剩余量r1(t)，表示如下式：
r1(t)＝s(t)-c1(t)
步骤六：将r1(t)当做新的原始信号，重复执行步骤一直步骤五，得到新的剩余量r2(t)，如此重复n次
r2(t)＝r1(t)-c2(t)
r3(t)＝r2(t)-c3(t)
…
rn(t)＝rn-1(t)-cn(t)
当第n个剩余量rn(t)已成为单调函数(monotonic function)，将无法再分解IMF时，整个EMD的分解过程完成，原始讯号s(t)可以表示成n个IMF分量与一个平均趋势(meantrend)分量rn(t)的组合，亦即
Figure BDA0002307761460000121
希尔伯特-黄变换单元中，使用希尔伯特-黄变换方法，得到各阶次的瞬时频率图，
希尔伯特-黄变换是在EMD分解的基础上，将IMF进行希尔伯特变换得到的，函数u(t)的希尔伯特变换是：
Figure BDA0002307761460000122
所述特征提取模块包括线性判别分析单元和盲源分离单元。
线性判别分析单元利用线性判别分析提取数据中的主成分，通过对监督数据的训练，形成故障特征向量，实现对未知试验数据的故障模式识别，线性判别分析单元通过载入不同故障模式的数据，计算特征评价指标值，选择指标值最大的特征值，在VMD分解后训练数据，导入待分类数据，得到故障诊断结果，并给出置信度，其具体方法为：
假设C个类中每一个类都有均值μi和相同的协方差∑。那么，类间的变化通过类均值的协方差来定义：
Figure BDA0002307761460000131
这里μ是各类均值的均值，在
Figure BDA0002307761460000132
的方向区分类由下式给出：
Figure BDA0002307761460000133
如果
Figure BDA0002307761460000134
是∑-1∑b的特征向量，等同于用对应的特征值进行分类。
如果∑-1∑b是可对角化矩阵，特征之间的变化性就会被保留在C-1个最大特征值对应的特征向量构成子空间内。
盲源分离单元利用盲源分离在没有任何先验知识的前提下，将源信号从混叠信号中提取和分离出来，盲源分离的具体方法为：
将盲源分离算法做出如下假定：具有m个独立的信号源s1(t),…,sm(t)和n个独立的观察量x1(t),…,xm(t)，观察量和信号源具有如下的关系：
x(t)＝As(t)
其中，x(t)＝[x1(t)…,xm(t)]T，s(t)＝[s1(t)…,sm(t)]T，A是一个n*m的系数矩阵，原问题变成了已知x(t)和s(t)的独立性，求对s(t)的估计问题，假设有如下公式：
y(t)＝Wx(t)
其中y(t)是对s(t)的估计，W是一个m*n系数矩阵；
寿命预测模块中，使用自回归滑动平均模型进行退化模型的建立和预测，寿命预测模块通过载入时间序列数据、模型参数辨识、参数估计和模型检验步骤实现建模过程，并在模型建立完成后进行时间序列预测，具体建模方法如下：
自回归滑动平均模型中包含了p个自回归项和q个移动平均项，自回归滑动平均模型表示为：
Figure BDA0002307761460000141
其中εt是白噪声序列，p和q都是非负整数，c为常数，φi和θi为待定系数；
使用贝叶斯信息量准则对模型参数进行辨识，贝叶斯信息量定义为：
BIC＝ln(n)k-2ln(L)
其中，k为模型参数个数，n为样本数量，L为似然函数；
使用极大似然估计法对自回归滑动平均模型中的参数进行估计，假设x1,x2,…,xn为独立同分布的采样，θ为模型参数,f为我们所使用的模型，遵循上述的独立同分布假设，参数为θ的模型f产生上述采样表示为：
f(x1,x2,…,xn|θ)＝f(x1|θ)×f(x2|θ)…,f(xn|θ)
似然定义为：
Figure BDA0002307761460000142
其中lnL(θ|x1,x2,…,xn)为对数似然，极大似然函数表示为：
Figure BDA0002307761460000151
模型检验对模型估计的残差值进行检验，使用自相关和偏自相关检验残差的自相关和偏自相关特性，
其中自相关公式为：
Figure BDA0002307761460000152
偏自相关公式为：
Figure BDA0002307761460000153
具体实施例
以一种齿轮箱为例，在数据载入模块中，可以载入TDMS文件和mat文件。在TDMS文件载入子模块中读取tdms数据后，可以显示数据中的通道组、通道信息，在点击不同通道后可以显示该通道的属性值、图和数据表。在mat文件载入子模块中读取mat数据后，可以显示数据中包含的所有数据变量名，在点击数据变量名后可以显示该数据波形图。
利用数据预处理模块对数据信号进行预处理，之后数据进入频域分析模块。
在频域分析模块模块中，可以自动绘制模拟信号的频谱图、周期图法功率谱、间接法功率谱。在细化谱设置模块中输入频率上限、频率下限和NFFT，可以绘制模拟信号的细化谱，之后数据进入时频域分析模块。
在时频域分析模块中进行经验模态分解和希尔伯特-黄变换，之后数据进入特征提取模块。
在特征提取模块中，利用线性判别分析和盲源分离进行特征提取。
寿命预测模块中，使用自回归滑动平均模型(ARMA)进行退化模型的建立和预测。该模块通过载入时间序列数据、模型参数辨识、参数估计、模型检验等步骤实现建模过程。并在模型建立完成后可以时间序列预测。
具体地寿命预测包括以下几个步骤：
步骤一：利用数据载入模块载入产品的数据；
步骤二：数据预处理步骤：数据载入模块总共采集约2T数据，在具体预测时，将采集的数据进行截段，将每30分钟的采集的数据分别截取前10秒钟的采集的数据作为该采集的数据段的分析对象，按照这个方法，所有采集的数据总共被分为2324个数据段。
图3所示为被截取后的采集的数据在整个实施例的时域变化过程，图3中，横坐标为数据的点数。
步骤三：进行特征提取步骤：采集的数据中含有大量的噪声信号与工频干扰信号，必须使用奇异值分解与重构的方法实现原始模拟信号的滤波与主成分提取。
奇异值分解降噪的关键在于选定阈值，奇异值阈值过大会丢失重要的冲击信号，奇异值阈值过小则会降低降噪效果。选择合适的阈值，在滤除噪声的前提下不丢失冲击信号，是准确提取出冲击特征信号的一个关键步骤，本实施例中，利用奇异值差分谱的奇异值阈值确定方法来寻找合适的奇异值阈值。
以10秒为一个单元对原始信号进行奇异值分解和重构，再将重构后的信号单元拼接，实现原始信号的主成分的提取。
步骤四：进行频域分析和时频域分析：在本实施例中，选取均方根、方差、标准差、峭度、偏度、波形因子、峰值因子、脉冲因子、裕度因子和方根幅值作为特征提取的时域指标。以10秒为一个数据段对经过奇异值分解与重构后的信号进行时域特征值计算，得到由不同时域指标表征的产品性能退化曲线，如图4所示。
图4中所表示的各通道不同时域特征的退化过程中，输入轴电机侧轴承y方向的退化趋势较为明显，这是受到了信号的传递路径、信号特性等多种因素共同影响的结果。最终选择方根幅值作为表征性能退化的特征参数。
步骤五：退化数据建模与预测，其还包括以下几个子步骤：
步骤51：使用一阶或者多阶差分方法使序列变得平稳，并通过单位根方法进行检验。
步骤52：确定p和q的值，构建ARIMA模型：
进行一阶差分之后，序列变为平稳序列，此时的模型为ARIMA(p，1，q)。本实施例中，使用贝叶斯信息准则(BIC)的方法，通过加入模型复杂度的惩罚项来避免由p和q的值过大导致的过拟合，其表达式为：
BIC＝klnn-2lnL
式中，k为模型参数个数，L为似然函数，n为样本数量。从一组可供选择的模型中选择最佳模型时，选择BIC最小的模型。通过计算得到p＝3，q＝2。
步骤53：模型残差检验：对模型进行残差检验，图5为残差检验的结果图，主要是查看残差是否接近正态分布，接近正态分布即为理想的残差；ACF和PACF检验残差的自相关和偏自相关，上述图片结构能够表明证明，残差是接近正态分布且相互独立的，因此步骤51的建模是符合要求的。
步骤54：寿命预测步骤：在寿命预测模型确定之后，对故障特征值的趋势进行预测，预测结果如图6所示。根据建模结果绘制故障特征值预测曲线，特征值达到失效阈值时，该产品的寿命值为3973小时。
本发明的效果如下：
1.本发明的系统覆盖了信号预处理、频域分析、时频域分析、特征提取以及寿命预测的整个数据处理流程，方法具有通用性，可以在旋转类机械产品信号分析中使用。
2.本发明相较于其他信号分析手段，该软件系统更具有针对性，系统集成了信号分析与预测最为常用的分析方法，并包含专用工具箱，降低了学习成本，提高了分析效率。
3.本发明支持不同格式的数据文件，可以实现与MATLAB、Labview等软件的联合分析，灵活性更强。
以上所述的实施例仅是对本发明的优选实施方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明权利要求书确定的保护范围内。

Claims (8)

1.一种旋转类机械产品剩余寿命预测系统，其特征在于：其包括数据载入模块、模拟信号生成模块、数据预处理模块、频域分析模块、时频域分析模块、特征提取模块以及寿命预测模块，数据载入模块、模拟信号生成模块、数据预处理模块、频域分析模块、时频域分析模块、特征提取模块以及寿命预测模块相互之间通讯连接；
所述数据载入模块用于载入TDMS文件和mat文件格式的产品数据信息；
所述模拟信号生成模块用于根据产品数据信息生成正弦波、三角波、方波、锯齿波或白噪声的模拟信号；
数据预处理模块对模拟信号进行处理，所述数据预处理模块包括去趋势单元、加窗单元、滤波单元以及奇异值分解单元，所述去趋势单元对模拟信号进行去趋势处理，消除信号中的直流分量；
所述加窗单元用于给模拟信号增加窗口，所述窗口的类型包括矩形窗、离散汉宁窗、离散哈明窗和布莱克曼窗，所述加窗单元能够显示原始波形、加窗后波形、加窗后波形时域图和加窗后波形频域图，模拟信号由加窗单元导出之后形成加窗后波形，
其中长度为N的离散化矩形窗时域表达式为：
Figure FDA0002307761450000011
长度为N的离散汉宁窗时域表达式为：
Figure FDA0002307761450000012
长度为N的离散哈明窗时域表达式为：
Figure FDA0002307761450000013
长度为N的布莱克曼窗时域表达式为：
Figure FDA0002307761450000021
所述滤波单元用于对模拟信号进行滤波，所述滤波单元能够显示模拟信号的原始时域信号、原始信号幅频谱、原始信号相频谱、滤波后时域信号、滤波后信号幅频谱和滤波后信号相频谱，所述滤波模块中能够设置滤波器类型、纹波、衰减、阶次、下限截止频率和上限截止频率，所述奇异值分解单元用于对模拟信号进行奇异值分解和重构；
所述奇异值分解单元能够对模拟信号进行奇异值分解和重构，对于长度小于或等于10000的数据进行奇异值分解，分解成功后绘制奇异值差分结果，输入重构阶次后，对分解后的信号进行重构，从而达到主成分提取和滤波的目的，奇异值分解具体为：
假设M为m*n阶矩阵，其中的矩阵元素均属于实数域或复数域，则存在一个分解使得
M＝U∑V*
其中，U是m*m阶酉矩阵，∑是m*n阶非负实数对角矩阵；而V*即V的共轭转置，是n*n阶酉矩阵；此分解称作M的奇异值分解；
所述频域分析模块能够自动绘制模拟信号的频谱图、周期图法功率谱和间接法功率谱，在细化谱设置模块中输入频率上限、频率下限和快速傅里叶变换点数，绘制模拟信号的细化谱，从而进行频域分析；
所述时频域分析模块包括经验模态分解单元和希尔伯特-黄变换单元，所述时频域分析模块能够对模拟信号进行重构，得到模拟信号各阶次的瞬时频率图，进行时频域分析，在经验模态分解单元中，使用经验模态分解方法将自动对模拟信号进行分解，绘制采集的数据及其内涵模态分量的时域图，以及采集的数据及其内涵模态分量的频谱图，设置重构信号的起始内涵模态分量和终止内涵模态分量，并对起始内涵模态分量和终止内涵模态分量中间的所有分量进行重构；
经验模态分解的具体流程如下：
步骤一：假设原始模拟信号为s(t)，找出原始模拟信号s(t)中的所有局部极大值以及局部极小值，接着利用三次样条，分别将局部极大值串连成上包络线，将局部极小值串连成下包络线；
步骤二：求出上包络线和下包络线的平均值，得到均值包络线m1(t)；
步骤三：将原始模拟信号s(t)与均值包络线ml(t)相减，得到第一个分量h1(t)；
h1(t)＝s(t)-m1(t)
步骤四：检查h1(t)是否符合IMF的条件，如果不符合，则回到步骤1并且将h1(t)当作原始模拟信号，进行第二次的筛选，亦即
h2(t)＝h1(t)-m2(t)
重复筛选k次，
hk(t)＝hk-1(t)-mk(t)
直到hk(t)符合IMF的条件，得到第一个IMF分量c1(t)，亦即
c1(t)＝hk(t)
步骤五：原始讯号s(t)减去c1(t)得到剩余量r1(t)，表示如下式：
r1(t)＝s(t)-c1(t)
步骤六：将r1(t)当做新的原始模拟信号，重复执行步骤一直步骤五，得到新的剩余量r2(t)，如此重复n次
r2(t)＝r1(t)-c2(t)
r3(t)＝r2(t)-c3(t)
…
rn(t)＝rn-1(t)-cn(t)
当第n个剩余量rn(t)已成为单调函数，将无法再分解IMF时，整个EMD的分解过程完成，原始讯号s(t)可以表示成n个IMF分量与一个平均趋势分量rn(t)的组合，亦即
Figure FDA0002307761450000041
所述希尔伯特-黄变换单元使用希尔伯特-黄变换方法，得到模拟信号各阶次的瞬时频率图，
希尔伯特-黄变换是在EMD分解的基础上，将IMF进行希尔伯特变换得到的，函数u(t)的希尔伯特变换公式是：
Figure FDA0002307761450000042
其中，x(t)为连续的时间信号，y(t)为经过希尔伯特变换后的信号；
所述特征提取模块包括线性判别分析单元和盲源分离单元，所述特征提取模块用于对模拟信号进行特征提取，获取模拟信号的故障征向量，并得到故障诊断结果，给出置信度，同时获得模拟信号的源信号；
线性判别分析单元利用线性判别分析提取模拟信号数据中的主成分，通过对监督数据的训练，形成故障特征向量，实现对未知试验数据的故障模式识别，线性判别分析单元通过载入不同故障模式的数据，计算特征评价指标值，选择指标值最大的特征值，在VMD分解后训练数据，导入待分类数据，得到故障诊断结果，并给出置信度；盲源分离单元利用盲源分离在没有任何先验知识的前提下，将源信号从混叠信号中提取和分离出来；
所述寿命预测模块使用自回归滑动平均模型进行退化模型的建立和预测，寿命预测模块通过载入时间序列数据、模型参数辨识、参数估计和模型检验步骤实现建模过程，并在模型建立完成后对退化特征的趋势进行预测，得到给定失效阈值下的机械产品寿命评估值。
2.根据权利要求1所述的旋转类机械产品剩余寿命预测系统，其特征在于：线性判别分析的具体方法为：
假设C个类中每一个类都有均值μi和相同的协方差∑，那么，类间的变化通过类均值的协方差来定义：
Figure FDA0002307761450000051
这里μ是各类均值的均值，在
Figure FDA0002307761450000052
的方向区分类由下式给出：
Figure FDA0002307761450000053
如果
Figure FDA0002307761450000054
是∑-1∑b的特征向量，等同于用对应的特征值进行分类；
如果∑-1∑b是可对角化矩阵，特征之间的变化性就会被保留在C-1个最大特征值对应的特征向量构成子空间内。
3.根据权利要求2所述的旋转类机械产品剩余寿命预测系统，其特征在于：盲源分离的具体方法为：
将盲源分离算法做出如下假定：具有m个独立的信号源s1(t),…,sm(t)和n个独立的观察量x1(t),…,xm(t)，观察量和信号源具有如下的关系：
x(t)＝As(t)
其中，x(t)＝[x1(t)…,xm(t)]T，s(t)＝[s1(t)…,sm(t)]T，A是一个n*m的系数矩阵，原问题变成了已知x(t)和s(t)的独立性，求对s(t)的估计问题，假设有如下公式：
y(t)＝Wx(t)
其中y(t)是对s(t)的估计，W是一个m*n系数矩阵。
4.根据权利要求1所述的旋转类机械产品剩余寿命预测系统，其特征在于：所述寿命预测模块建模的具体方法如下：
步骤a、自回归滑动平均模型中包含了p个自回归项和q个移动平均项，自回归滑动平均模型表示为：
Figure FDA0002307761450000055
其中εt是白噪声序列，p和q都是非负整数，c为常数，φi和θi为待定系数；
使用贝叶斯信息量准则对模型参数进行辨识，贝叶斯信息量定义为：
BIC＝ln(n)k-2ln (L)
其中，k为模型参数个数，n为样本数量，L为似然函数；
步骤b、使用极大似然估计法对自回归滑动平均模型中的参数进行估计，假设x1,x2,…,xn为独立同分布的采样，θ为模型参数,f为我们所使用的模型，遵循上述的独立同分布假设，参数为θ的模型f产生上述采样表示为：
f(x1,x2,…,xn|θ)＝f(x1|θ)×f(x2|θ)…,f(xn|θ)
似然定义为：
Figure FDA0002307761450000061
其中lnL(θ|x1,x2,…,xn)为对数似然，极大似然函数表示为：
Figure FDA0002307761450000062
步骤c、模型检验对模型估计的残差值进行检验，使用自相关和偏自相关检验残差的自相关和偏自相关特性，
其中自相关公式为：
Figure FDA0002307761450000063
偏自相关公式为：
Figure FDA0002307761450000064
步骤d、根据步骤a确定的自回归滑动平均模型，对故障特征值的趋势展开预测，绘制故障特征值预测曲线，分析当达到失效阈值时，对应的寿命值。
5.根据权利要求1所述的旋转类机械产品剩余寿命预测系统，其特征在于：数据载入模块在TDMS文件载入单元中读取tdms数据后，能够显示数据中的通道组以及通道信息，在点击不同通道后能够显示该通道的属性值、图和数据表，在mat文件载入单元中读取mat数据后，能够显示数据中包含的所有数据变量名，在点击数据变量名后能够显示该数据波形图。
6.根据权利要求1所述的旋转类机械产品剩余寿命预测系统，其特征在于：所述模拟信号生成模块能够根据需要设置正弦波、三角波、方波、锯齿波和白噪声的幅值、频率、数据长度、偏置和采样率，在波形叠加模式中，并能够将基础波形进行相加、相乘运算，进而模拟生成各种复杂形式的仿真信号。
7.根据权利要求1所述的旋转类机械产品剩余寿命预测系统，其特征在于：所述滤波器类型包括低通滤波器、高通滤波器、带通滤波器以及带阻滤波器。
8.根据权利要求7所述的旋转类机械产品剩余寿命预测系统，其特征在于：所述滤波单元的滤波器设计方法为Butterworth，Butterworth滤波器幅频响应关系如下:
Figure FDA0002307761450000071
式中，Ω、Ωc、N分别表示频率、转折频率、系统阶数。
CN201911246253.6A 2019-12-08 2019-12-08 一种旋转类机械产品剩余寿命预测系统 Active CN110991055B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911246253.6A CN110991055B (zh) 2019-12-08 2019-12-08 一种旋转类机械产品剩余寿命预测系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911246253.6A CN110991055B (zh) 2019-12-08 2019-12-08 一种旋转类机械产品剩余寿命预测系统

Publications (2)

Publication Number Publication Date
CN110991055A true CN110991055A (zh) 2020-04-10
CN110991055B CN110991055B (zh) 2022-09-16

Family

ID=70091208

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911246253.6A Active CN110991055B (zh) 2019-12-08 2019-12-08 一种旋转类机械产品剩余寿命预测系统

Country Status (1)

Country Link
CN (1) CN110991055B (zh)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111476220A (zh) * 2020-06-03 2020-07-31 中国南方电网有限责任公司超高压输电公司大理局 一种换流阀空冷器故障定位方法
CN111854920A (zh) * 2020-07-24 2020-10-30 贵州电网有限责任公司 一种基于dvs振动监测信号的预处理方法及系统
CN112053009A (zh) * 2020-09-30 2020-12-08 华人运通(江苏)技术有限公司 一种故障预测方法、装置、系统及存储介质
CN112798279A (zh) * 2020-12-30 2021-05-14 杭州朗阳科技有限公司 一种新的诊断电机轴承故障的检测方法
CN113375933A (zh) * 2021-05-31 2021-09-10 中国矿业大学 一种刮板输送机故障诊断系统和诊断方法
CN113836756A (zh) * 2021-11-29 2021-12-24 山东华尚电气有限公司 一种立体卷铁心变压器的退火工艺智能监测方法及系统
CN117332214A (zh) * 2023-11-30 2024-01-02 中国航发上海商用航空发动机制造有限责任公司 基于小波变换与频域相干函数融合式喘振报警方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106909756A (zh) * 2017-03-29 2017-06-30 电子科技大学 一种滚动轴承剩余寿命预测方法
CN108875841A (zh) * 2018-06-29 2018-11-23 国家电网有限公司 一种抽蓄机组振动趋势预测方法
CN109211548A (zh) * 2018-08-31 2019-01-15 沃德传动（天津）股份有限公司 一种机械故障诊断方法
WO2019043600A1 (en) * 2017-08-29 2019-03-07 Emerson Electric (Us) Holding Corporation (Chile) Limitada ESTIMATOR OF REMAINING USEFUL LIFE ESTIMATOR
US20190272354A1 (en) * 2019-03-06 2019-09-05 Beihang University Method for degradation modeling and lifetime prediction considering recoverable shock damages

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106909756A (zh) * 2017-03-29 2017-06-30 电子科技大学 一种滚动轴承剩余寿命预测方法
WO2019043600A1 (en) * 2017-08-29 2019-03-07 Emerson Electric (Us) Holding Corporation (Chile) Limitada ESTIMATOR OF REMAINING USEFUL LIFE ESTIMATOR
CN108875841A (zh) * 2018-06-29 2018-11-23 国家电网有限公司 一种抽蓄机组振动趋势预测方法
CN109211548A (zh) * 2018-08-31 2019-01-15 沃德传动（天津）股份有限公司 一种机械故障诊断方法
US20190272354A1 (en) * 2019-03-06 2019-09-05 Beihang University Method for degradation modeling and lifetime prediction considering recoverable shock damages

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
肖顺根等: "基于EEMD和PCA滚动轴承性能退化指标的提取方法", 《江南大学学报(自然科学版)》 *
边智 等: "综合加速寿命试验与仿真的密封圈性能退化规律研究", 《中国测试》 *

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111476220A (zh) * 2020-06-03 2020-07-31 中国南方电网有限责任公司超高压输电公司大理局 一种换流阀空冷器故障定位方法
CN111854920A (zh) * 2020-07-24 2020-10-30 贵州电网有限责任公司 一种基于dvs振动监测信号的预处理方法及系统
CN112053009A (zh) * 2020-09-30 2020-12-08 华人运通(江苏)技术有限公司 一种故障预测方法、装置、系统及存储介质
CN112053009B (zh) * 2020-09-30 2023-08-01 华人运通(江苏)技术有限公司 一种故障预测方法、装置、系统及存储介质
CN112798279A (zh) * 2020-12-30 2021-05-14 杭州朗阳科技有限公司 一种新的诊断电机轴承故障的检测方法
CN113375933A (zh) * 2021-05-31 2021-09-10 中国矿业大学 一种刮板输送机故障诊断系统和诊断方法
CN113836756A (zh) * 2021-11-29 2021-12-24 山东华尚电气有限公司 一种立体卷铁心变压器的退火工艺智能监测方法及系统
CN117332214A (zh) * 2023-11-30 2024-01-02 中国航发上海商用航空发动机制造有限责任公司 基于小波变换与频域相干函数融合式喘振报警方法
CN117332214B (zh) * 2023-11-30 2024-04-12 中国航发上海商用航空发动机制造有限责任公司 基于小波变换与频域相干函数融合式喘振报警方法

Also Published As

Publication number Publication date
CN110991055B (zh) 2022-09-16

Similar Documents

Publication Publication Date Title
CN110991055B (zh) 一种旋转类机械产品剩余寿命预测系统
US11022633B2 (en) Enhanced system and method for conducting PCA analysis on data signals
Yu et al. Time-reassigned multisynchrosqueezing transform for bearing fault diagnosis of rotating machinery
Manhertz et al. STFT spectrogram based hybrid evaluation method for rotating machine transient vibration analysis
Kantz et al. Dimension estimates and physiological data
Yen et al. Wavelet packet feature extraction for vibration monitoring
Yang et al. Component extraction for non-stationary multi-component signal using parameterized de-chirping and band-pass filter
US7401008B2 (en) Method, computer program, and system for intrinsic timescale decomposition, filtering, and automated analysis of signals of arbitrary origin or timescale
US5251151A (en) Method and apparatus for diagnosing the state of a machine
Srinivasan et al. A modified empirical mode decomposition (EMD) process for oscillation characterization in control loops
Hao et al. Underdetermined source separation of bearing faults based on optimized intrinsic characteristic-scale decomposition and local non-negative matrix factorization
KR20200075148A (ko) 문제소음 발음원 식별을 위한 소음데이터의 인공지능 장치 및 전처리 방법
CN109635864B (zh) 一种基于数据的容错控制方法及装置
CN111307438A (zh) 一种基于信息熵的旋转机械振动故障诊断方法及其系统
Chen et al. Nonstationary signal denoising using an envelope-tracking filter
CN115954017A (zh) 一种基于hht的发动机小样本声音异常故障识别方法及系统
CN111898644A (zh) 一种无故障样本下航天液体发动机健康状态智能识别方法
Antoni et al. Standalone extraction of tonal components from aeroacoustic signals
Arundale Studying change over time: Criteria for sampling from continuous variables
Dang et al. A fault diagnosis method for one-dimensional vibration signal based on multiresolution tlsDMD and approximate entropy
Cicone et al. Jot: a variational signal decomposition into jump, oscillation and trend
CA2390384A1 (en) Improved spectrum analysis method and apparatus
CN110222390B (zh) 基于小波神经网络的齿轮裂纹识别方法
Reidy An introduction to random processes for the spectral analysis of speech data
Deng et al. Soft fault feature extraction in nonlinear analog circuit fault diagnosis

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