基于双时间尺度粒子滤波的机电系统剩余寿命预测方法与流程

本发明涉及预测设备的剩余寿命的领域，尤其是基于双时间尺度粒子滤波的机电系统剩余寿命预测方法。

背景技术：

随着科学技术和现代工业的快速发展，以及随着工业系统和大量机电设备复杂程度和综合性能的提高，准确快速的故障诊断和精确的剩余寿命预测变得越来越重要。有效的故障诊断和预测方法能够快速准确的定位系统或设备的故障源，并精确的预测出设备的剩余寿命，这些关键信息为系统或设备的维护提供必要条件，从而避免灾难性事故的发生以及不必要的经济损失。

目前常用的非线性系统的故障参数的估计方法包括粒子滤波算法、卡尔曼滤波算法以及最小二乘估计法，但是由于非线性系统的故障参数具有缓慢的时变特性，而其系统状态具有快速的时变特性，当系统出现随时间而缓慢变化的故障参数时，利用上述方法进行系统状态和故障参数的联合估计时会浪费大量的计算资源在故障参数基本不变的时间段，导致降低了故障诊断的效率，从而导致降低了剩余寿命预测的实时性。

技术实现要素：

为了克服上述现有技术中的缺陷，本发明提供基于双时间尺度粒子滤波的机电系统剩余寿命预测方法，解决了系统状态和故障参数在进行联合估计时会浪费大量的计算资源在故障参数基本不变的时间段的问题，提高了故障诊断的效率，保障了剩余寿命预测的实时性。

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

基于双时间尺度粒子滤波的机电系统剩余寿命预测方法，包括以下步骤：

s1，采用欧拉离散化方法对连续的系统状态空间方程进行离散化处理，得到双时间尺度下离散的系统状态空间方程；所述系统状态空间方程为机电系统的状态空间方程；

s2，利用双时间尺度粒子滤波算法对系统状态和故障参数进行联合估计，得到故障参数的估计值，并通过故障参数的估计值识别系统故障源；所述故障参数为故障集合中的参数；

s3，预测系统的退化模型，利用所述故障参数的估计值和故障失效阈值fend计算系统的剩余寿命。

步骤s1中，将故障集合中的故障参数视为系统的特殊状态对连续的系统状态空间方程进行增广，并在双时间尺度下对增广后的连续的系统状态空间方程进行离散化处理，得到双时间尺度下离散的系统状态空间方程；

所述连续的系统状态空间方程为：

其中，t表示实际时间，是系统状态的一阶导数，x(t)是系统状态，u(t)是系统输入，w是系统状态的过程噪声；y(t)是系统输出，v是系统输出的测量噪声；f(·)和g(·)表示非线性函数；

所述双时间尺度包括宏观时间尺度和微观时间尺度；微观时间尺度下的微观时间l的循环周期为从1开始到l结束，l＞1；当l＝l时，宏观时间尺度下的宏观时间k加1，即宏观时间由k变换为k+1，所述微观时间l重新开始循环，微观时间l由l变为1；

所述双时间尺度下离散的系统状态空间方程为：

其中，下标k表示宏观时间为k，下标l表示微观时间为l，下标l+1表示微观时间为l+1；下标k+1表示宏观时间为k+1；

xk，l表示宏观时间为k、微观时间为l的时刻的系统状态；uk，l是宏观时间为k、微观时间为l的时刻的系统输入；θk表示宏观时间为k时刻的故障参数，视为系统的特殊状态；ωk，l表示宏观时间为k、微观时间为l的时刻的系统状态的过程噪声；tk，l表示宏观时间为k、微观时间为l的时刻；xk，l+1表示宏观时间为k、微观时间为l+1时刻的系统状态；f(·)表示非线性函数；

ρk表示宏观时间为k时刻的故障参数的过程噪声；θk+1表示宏观时间为k+1时刻的故障参数；

yk，l表示宏观时间为k、微观时间为l的时刻的系统输出；vk，l表示宏观时间为k、微观时间为l的时刻的系统输出的测量噪声；h(·)表示非线性函数。

步骤s2中，包括以下具体步骤：

s21，基于双时间尺度下离散的系统状态空间方程，进行初始化设置：

l＝l0；k＝k0；x＝x0，0+ω0；θ＝θ0+ρ0；

其中，将微观时间l初始化为l0；将宏观时间k初始化为k0；将系统状态x初始化为初始的系统状态x0，0加上初始的系统状态的过程噪声ω0；将故障参数θ初始化为初始的故障参数θ0加上初始的故障参数的过程噪声ρ0；

上标i表示第i个粒子；表示系统状态的第i个粒子的初始权重；表示故障参数的第i个粒子的初始权重；每个系统状态和每个故障参数的总粒子数均为np；

s22，微观时间l在一个循环周期内，即微观时间l从1开始到l结束，进行系统状态的估计，得到系统状态的估计值包括以下具体步骤：

s221，对系统状态进行采样，得到系统状态随机变量，为其中，为宏观时间为k、微观时间为l的时刻采样得到的系统状态的第i个粒子；为重要性概率密度函数；

s222，分别对系统状态中的np个粒子的权重进行计算，其中，计算系统状态的第i个粒子的权重的方式为：

其中，为系统状态转移概率密度函数；为由测量函数决定的概率密度函数；为重要性概率密度函数；为宏观时间为k、微观时间为l的时刻的系统状态的第i个粒子的权重；若k＝0，l＝1，则为系统状态的第i个粒子的初始权重否则，为系统状态的第i个粒子在宏观时间为k、微观时间为l-1的时刻，即前一微观时刻进行归一化处理后得到的权重；

s223，分别对系统状态中的np个粒子的权重进行归一化处理，其中，对系统状态i个粒子的权重进行归一化处理的方式为：

s224，对系统状态中的np个粒子进行重采样判断：

其中，neff是决定系统状态是否需要重采样的重要指标，neff的值的大小反映粒子退化的严重程度，若neff＜nt时，nt是预先设定的一个阈值，则系统状态的所有粒子的权重均被重新设置为1/np，并进行重采样，重新执行步骤s221～s224；否则，执行步骤s225；

s225，计算宏观时间为k、微观时间为l的时刻的系统状态的估计值计算方式为：

s23，每当微观时间l＝l时，宏观时间k加1，即宏观时间k＝k+1，且在每当微观时间l＝l时，进行故障参数的估计，得到故障参数的估计值包括以下具体步骤：

s231，对故障参数进行采样，得到故障参数随机变量，为其中，为宏观时间为k、微观时间为l的时刻采样得到的故障参数的第i个粒子；q(θk|θk-1，yk)为重要性概率密度函数；yk为宏观时间为k时刻的系统输出；

s232，分别对故障参数中的np个粒子的权重进行计算，其中，计算故障参数的第i个粒子的权重的方式为：

其中，是状态转移概率密度函数；是由测量函数决定的概率密度函数；是重要性概率密度函数；为宏观时间为k时刻的故障参数的第i个粒子的权重；若k＝1，则为故障参数的第i个粒子的初始权重否则，为故障参数的第i个粒子在宏观时间为k时刻，即前一宏观时刻进行归一化处理后得到的权重；

s233，分别对故障参数中的np个粒子的权重进行归一化处理，其中对故障参数的第i个粒子的权重进行归一化处理的方式为：

s234，对故障参数的np个粒子进行重采样判断：

其中，neff是决定故障参数是否需要重采样的重要指标；neff的值的大小反映权重粒子退化的严重程度，当neff＜nt时，nt是预先设定的一个阈值，所有故障参数的粒子的权重均被重新设置为1/np，并进行重采样，重新执行步骤s221～s224；否则，执行步骤s225；

s235，计算宏观时间为k时刻的故障参数的估计值计算方式为：

s24，以此类推，循环上述步骤s22和步骤s23，直至终止条件满足；在仿真过程中，所述终止条件为宏观时间达到了用户设定的运行时间，在实际预测过程中，所述终止条件为电机达到停机状态。

s25，根据所述故障参数的估计值作出故障参数估计值的变化曲线，并根据所述故障参数估计值的变化曲线是否远离该故障参数的标称值来识别系统故障源，若远离，则该故障参数即为系统故障源。

步骤s25中，若故障参数得估计值的变化曲线远离该故障参数的标称值的10％以上时，则该故障参数即为系统故障源。

步骤s3中，包括以下具体步骤：

s31，根据系统故障源所对应的故障参数的估计值的变化曲线拟合指数退化模型，所述指数退化模型fexp：

其中，是系统故障源所对应的故障参数的标称值，ε是非0常数，代表系统的退化参数，t0是未知的故障发生时间；

s32，分别对系统故障源所对应的故障参数的np个粒子的剩余寿命进行计算，

其中，计算该故障参数的第i个粒子的剩余寿命rulⁱ的方式为：

根据拟合的指数退化模型，先计算出第i个粒子达到故障失效阈值fend的寿命终止时间eolⁱ，所述寿命终止时间eolⁱ的计算方式为：

再根据第i个粒子的寿命终止时间eolⁱ计算出第i个粒子的剩余寿命rulⁱ，所述剩余寿命rulⁱ的计算方式为：

rulⁱ＝eolⁱ-tε

其中，tε表示故障源被诊断出来的时间，即指数退化模型被预测出来的时间；

s33，根据系统故障源所对应的故障参数的np个粒子的剩余寿命计算系统的剩余寿命rulmean，所述系统的剩余寿命rulmean的计算方式为：

本发明的优点在于：

(1)本发明避免了进行联合估计时浪费大量的计算资源在故障参数基本不变的时间段，提高了系统状态和故障参数联合估计的效率，能够更好的将计算资源利用在剩余寿命的实时预测上。

(2)本发明通过双时间尺度粒子滤波算法有效的估计出故障参数的估计值，并及时准确的隔离出故障源。

(3)本发明通过拟合指数退化模型，精确的预测系统的剩余寿命，为系统或设备的维护提供有用信息。

附图说明

图1为本发明的方法流程图。

图2为本发明的系统中注入阶梯型跳变故障时，故障参数jm估计值的变化曲线。

图3为本发明的系统中注入阶梯型跳变故障时，故障参数fm估计值的变化曲线。

图4为本发明的系统中注入阶梯型跳变故障时，故障参数fec估计值的变化曲线。

图5为本发明的指数退化模型的预测结果。

图6为本发明的系统剩余寿命概率分布的预测结果。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本实施例中，输入信号后的5s时，先给系统注入故障参数fec的阶梯型跳变故障，再通过本发明的预测方法去判断故障源并预测剩余寿命，已验证本发明能够准确判断故障源并预测剩余寿命。

由图1所示，基于双时间尺度粒子滤波的机电系统剩余寿命预测方法，包括以下步骤：

s1，采用欧拉离散化方法对机电系统的连续的系统状态空间方程进行离散化处理，得到双时间尺度下离散的系统状态空间方程；

s2，利用双时间尺度粒子滤波算法对系统状态和故障参数进行联合估计，得到故障参数的估计值，并通过故障参数的估计值识别系统故障源；

s3，预测系统参数的退化模型，利用所述故障参数的估计值和故障失效阈值fend计算系统的剩余寿命。

步骤s1中，将故障集合中的故障参数视为系统的特殊状态对连续的系统状态空间方程进行增广，并在双时间尺度下对增广后的连续的系统状态空间方程进行离散化处理，得到双时间尺度下离散的系统状态空间方程。

所述连续的系统状态空间方程为：

其中，t表示实际的连续时间，x(t)是系统状态，是系统状态的一阶导数，u(t)是系统输入，y(t)是系统输出，w是系统状态的过程噪声，v是系统输出的测量噪声，f(·)和g(·)表示非线性函数；

所述双时间尺度包括宏观时间和微观时间，且微观时间l的循环周期为从1开始到l结束，l＞1，当l＝l时，所述宏观时间k加1，即宏观时间由k变换为k+1，所述微观时间l重新开始循环，微观时间l由l变为1；

所述双时间尺度下离散的系统状态空间方程为：

其中，下标k表示宏观时间为k，下标l表示微观时间为l，下标l+1表示微观时间为l+1；下标k+1表示宏观时间为k+1；

xk，l表示宏观时间为k、微观时间为l的时刻的系统状态；uk，l是宏观时间为k、微观时间为l的时刻的系统输入；θk表示宏观时间为k时刻的故障参数，视为系统的特殊状态；ωk，l表示宏观时间为k、微观时间为l的时刻的系统状态的过程噪声；tk，l表示宏观时间为k、微观时间为l的时刻；xk，l+1表示宏观时间为k、微观时间为l+1时刻的系统状态；f(·)表示非线性函数；

ρk表示宏观时间为k时刻的故障参数的过程噪声；θk+1表示宏观时间为k+1时刻的故障参数；

yk，l表示宏观时间为k、微观时间为l的时刻的系统输出；vk，l表示宏观时间为k、微观时间为l的时刻的系统输出的测量噪声；h(·)表示非线性函数。

本实施例中，输入信号后的5s时，先给系统注入故障参数fec的阶梯型跳变故障，再通过本发明的预测方法去判断故障源并预测剩余寿命，已验证本发明能够准确判断故障源并预测剩余寿命。

本实施例中，将隔离出的可能故障参数的集合{jm，fm&fec}视为系统的特殊状态对已知的的系统状态进行增广，并在双时间尺度下对增广后的连续的系统状态空间方程进行离散化处理，得到双时间尺度下离散的系统状态空间方程。

本实施例中，所述双时间尺度下离散的系统状态空间方程为：

其中，代表状态向量，代表输出向量，和分别代表过程噪声和测量噪声，q1、q2、q3表示参数的噪声协方差，α1、α2、α3是0到1之间均匀分布的随机数，tl是采样时间，vin是系统输入电压，k1为输入电压与电流的转换比，k2是电流与机械转矩的转换比，k为电机传动轴的刚度，jm是电机机械部分转动惯量，fm是电机机械部分粘摩擦系数，fec是电机机械部分库伦摩擦转矩，js是负载部分转动惯量，fs是负载部分粘摩擦系数，fsc是电机负载部分库伦摩擦转矩，θe和θs分别表示编码器测得的电机和负载转角，和表示电机和负载的角速度。

步骤s2中，包括以下具体步骤：

s21，基于双时间尺度下离散的系统状态空间方程，进行初始化设置：

l＝l0；k＝k0；x＝x0，0+ω0；θ＝θ0+ρ0；

其中，将微观时间l初始化为l0；将宏观时间k初始化为k0；将系统状态x初始化为初始的系统状态x0，0加上初始的系统状态的过程噪声ω0；将故障参数θ初始化为初始的故障参数θ0加上初始的故障参数的过程噪声ρ0；

上标i表示第i个粒子；表示系统状态的第i个粒子的初始权重；表示故障参数的第i个粒子的初始权重；每个系统状态和每个故障参数的总粒子数均为np；

s22，微观时间l在一个循环周期内，即微观时间l从1开始到l结束，进行系统状态的估计，得到系统状态的估计值包括以下具体步骤：

s221，对系统状态进行采样，得到系统状态随机变量，为其中，为宏观时间为k、微观时间为l的时刻采样得到的系统状态的第i个粒子；为重要性概率密度函数，概率密度函数的计算方式参见现有技术；

s222，分别对系统状态中的np个粒子的权重进行计算，其中，计算系统状态的第i个粒子的权重的方式为：

其中，为系统状态转移概率密度函数；为由测量函数决定的概率密度函数；为重要性概率密度函数；为宏观时间为k、微观时间为l的时刻的系统状态的第i个粒子的权重；若k＝0，l＝1，则为系统状态的第i个粒子的初始权重否则，为系统状态的第i个粒子在宏观时间为k、微观时间为l-1的时刻，即前一微观时刻进行归一化处理后得到的权重；

s223，分别对系统状态中的np个粒子的权重进行归一化处理，其中，对系统状态i个粒子的权重进行归一化处理的方式为：

s224，对系统状态中的np个粒子进行重采样判断：

其中，neff是决定系统状态是否需要重采样的重要指标，neff的值的大小反映粒子退化的严重程度，若neff＜nt时，nt是预先设定的一个阈值，则系统状态的所有粒子的权重均被重新设置为1/np，并进行重采样，重新执行步骤s221～s224；否则，执行步骤s225；

s225，计算宏观时间为k、微观时间为l的时刻的系统状态的估计值计算方式为：

s23，每当微观时间l＝l时，宏观时间k加1，即宏观时间k＝k+1，且在每当微观时间l＝l时，进行故障参数的估计，得到故障参数的估计值包括以下具体步骤：

s231，对故障参数进行采样，得到故障参数随机变量，为其中，为宏观时间为k、微观时间为l的时刻采样得到的故障参数的第i个粒子；q(θk|θk-1，yk)为重要性概率密度函数；

s232，分别对故障参数中的np个粒子的权重进行计算，其中，计算故障参数的第i个粒子的权重的方式为：

其中，是状态转移概率密度函数；是由测量函数决定的概率密度函数；是重要性概率密度函数；为宏观时间为k时刻的故障参数的第i个粒子的权重；若k＝1，则为故障参数的第i个粒子的初始权重否则，为故障参数的第i个粒子在宏观时间为k时刻，即前一宏观时刻进行归一化处理后得到的权重；

s233，分别对故障参数中的np个粒子的权重进行归一化处理，其中对故障参数的第i个粒子的权重进行归一化处理的方式为：

s234，对故障参数的np个粒子进行重采样判断：

其中，neff是决定故障参数是否需要重采样的重要指标；neff的值的大小反映权重粒子退化的严重程度，当neff＜nt时，nt是预先设定的一个阈值，所有故障参数的粒子的权重均被重新设置为1/np，并进行重采样，重新执行步骤s221～s224；否则，执行步骤s225；

s235，计算宏观时间为k时刻的故障参数的估计值计算方式为：

s24，以此类推，循环上述步骤s22和步骤s23，直至终止条件满足在仿真过程中，所述终止条件为宏观时间达到了用户设定的运行时间，在实际预测过程中，所述终止条件为电机达到停机状态；本实施例中，设定的运行时间为30s。

s25，根据所述故障参数的估计值作出故障参数估计值的变化曲线，并根据所述故障参数估计值的变化曲线是否远离该故障参数的标称值来识别系统故障源，若故障参数得估计值的变化曲线远离该故障参数的标称值的10％以上时，则该故障参数即为系统故障源。

由图2-图4所示，分别为故障参数jm、fm、fec估计值的变化曲线，纵坐标为故障参数的估计值，横坐标为步长，即为采样点数，且每个采样点之间的采样间隔为0.01s。

其中，由图2和图3所示，故障参数jm和故障参数fm的标称值分为0.0001和0.001，且故障参数jm和故障参数fm的估计值均未远离其对应的标称值；由图4所示，故障参数fec的标称值为0.08，且故障参数fec的估计值远离该故障参数的标称值且远离的范围超过标称值的10％，并呈现出阶梯型跳变故障，即故障参数fec为系统故障源。

根据步骤s2得出的故障参数的估计值和系统故障源，与给系统注入故障参数fec的阶梯型跳变故障一致，因此，本发明精准定位出了系统的故障源。

步骤s3中，包括以下具体步骤：

s31，根据系统故障源所对应的故障参数的估计值的变化曲线拟合指数退化模型，所述指数退化模型fexp：

其中，是系统故障源所对应的故障参数的标称值，ε是非0常数，代表系统的退化参数，t0是未知的故障发生时间；

s32，分别对系统故障源所对应的故障参数的np个粒子的剩余寿命进行计算，

其中，计算该故障参数的第i个粒子的剩余寿命rulⁱ的方式为：

根据拟合的指数退化模型，先计算出第i个粒子达到故障失效阈值fend的寿命终止时间eolⁱ，所述寿命终止时间eolⁱ的计算方式为：

再根据第i个粒子的寿命终止时间eolⁱ计算出第i个粒子的剩余寿命rulⁱ，所述剩余寿命rulⁱ的计算方式为：

rulⁱ＝eolⁱ-tε

其中，tε表示故障源被诊断出来的时间，即指数退化模型被预测出来的时间；

s33，根据系统故障源所对应的故障参数的np个粒子的剩余寿命计算系统的剩余寿命rulmean，所述系统的剩余寿命rulmean的计算方式为：

由图5所示，故障参数的标称值设置的故障失效阈值fend＝0.18，且根据图4的故障参数fec估计值的变化曲线在第二次跳变，得出退化参数ε＝0.06，且将此时的时间作为指数退化模型被预测出来的时间tε＝10s；并根据图4的故障参数fec估计值的变化曲线，得出检测到系统故障的时间td＝5.2035s，并将td近似的看作系统故障发生时间t0＝5s，与给系统注入故障参数fec的阶梯型跳变故障的时间一致，因此，本发明及时的精准定位出了系统的故障源。

由图6所示，系统剩余寿命真实值在置信区间内，且估计均值和真实值基本吻合，因此，本发明能精确的计算出系统剩余寿命。

以上仅为本发明创造的较佳实施例而已，并不用以限制本发明创造，凡在本发明创造的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明创造的保护范围之内。