CN104506162A - 基于ls-svr建模的高阶粒子滤波器的故障预示方法 - Google Patents

Abstract

本发明提供一种基于LS-SVR建模的高阶粒子滤波机械故障预示方法，目的是通过提取从正常状态演化到故障状态的信号特征，用数据驱动的方式建立状态方程，并能够实现在线实时的预测。本发明对于提供的机械振动信号，提取相应的信号的特征，例如：时域的6种统计特性并应用LLE方法对其进行低维特征提取，然后利用LS-SVR构建m阶HMM模型，并用低维特征对其进行训练，最够将其融入到粒子滤波算法框架中构成m阶粒子滤波，最后对机械故障进行预示。本发明和现有预示方法相比较具有预示精度高，适用范围广的特点。

Description

基于LS-SVR建模的高阶粒子滤波器的故障预示方法

技术领域

本发明涉及故障诊断方法领域，具体地说是用最小二乘支持向量机(LS-SVM)通过学习的方式建立高阶粒子滤波器中的状态方程，用于设备故障预示的方法。

背景技术

基于状态的维修(CBM)取代了传统的按照计划表的故障维护，成为各种工程***实施维修的优选方案，以保证它们的可靠性、安全性和可用性，CBM利用运行时间数据，以确定/预测机器状态，从而确定当前/未来的故障状态。通过预示机器状态来避免机器失效和突发性故障的发生。CBM实现技术主要包括传感和监测、信息处理、故障诊断和故障预示算法，能够准确及时检测早期故障和预测其余发生故障的组件的使用寿命。预示能力是其中的重要组成部分，可以准确地预测组件或子***未来状况和剩余使用寿命。

机械故障预示方法可以被分为两大类：基于模型的方法和基于数据驱动的方法。基于模型的方法是用数学模型来预测故障的演化趋势。给定一个***的模型，基于模型的方法能提供准确的预测估计。但是通常在实际情况下，难以设计准确的模型，尤其是当故障传播的过程是复杂的或并没有完全理解。另一种基于数据驱动的方法是采用收集的状态数据来建立故障传播模型。有很多数据驱动的方法，例如：马尔科夫模型、递归向神经网络、自适应模糊神经推理***、支持向量回归和最小二乘支持向量回归(LS-SVR)等，都是在机械故障预示领域中流行的预测器。

在数据驱动的方法里，模糊神经网络和支持向量回归已经在机械状态退化预测中取得了成功的应用。但是在机械状态随时间动态实时变化的情况下，如果在预测过程没有考虑到状态的动态变化，那么模糊神经网络和支持向量回归方法将不能做出准确的预测。

近几十年来，在统计学和各种工程领域中，许多专家学者都致力于动态***实时估计问题的研究和分析，发源于17世纪英国牧师T.R.Bayes的贝叶斯理论为动态***的状态估计问题提供了严格的理论框架，它利用所有的已知信息来构造***状态变量的后验概率密度函数，即用***模型预测状态的后验概率，再利用最新的观测值进行修正。状态的各种统计值如均值、方差等都可从后验概率密度函数中计算获得。20世纪90年代初，随着计算机计算能力和存储量的迅速提升，一种基于递推贝叶斯估计和蒙特卡洛方法结合的实时在线仿真算法——粒子滤波器，逐渐受到人们的重视。粒子滤波器是一种基于仿真的方法，它利用状态空间中一组带权值的随机样本粒子逼近目标状态变量的概率密度函数，每个样本代表目标的一个可能状态，综合所有的粒子状态可以得到目标状态的最小方差估计。该算法不受模型线性、高斯假设的约束，适用于任意非线性非高斯动态***。

在目前国内外公开的文献中，Zio E,Peloni G.Particle filtering prognostic estimationof the remaining useful life of nonlinear components.Reliability Engineering and SystemSafety 2011；96(3):403-409.中提出了基于粒子滤波的机械部件状态估计方法。该方法应用基于模型的状态方程和观测方程，结合蒙特卡洛仿真技术估计退化部件的每个状态的先验概率密度函数，避免了卡尔曼滤波技术对***简单的线性和高斯噪声的假设。提供了一个鲁棒性更强的状态预测框架。该方法被应用于断裂故障的剩余使用寿命的估计中，得出了满意的结果。Zhang L,Li X S,Yu J S,et al.A fault prognostic algorithmbased on Gaussian mixture model particle filter.Acta Aeronautica et Astronautica Sinica2009；30(2):319-324中提出一种预测算法，该算法的状态估计阶段，采用联合估计和粒子滤波同时估计对象***故障演化模型状态和未知参数的后验分布。在算法的状态预测阶段，采用了两种不同的计算方法：一种方法是对状态变量当前时刻的后验分布进行迭代采样，从而获得未来时刻的状态变量的先验分布；另一种方法是采用数据驱动的方法预测未来一段时间内对象***的量测信息，从而将未来时刻状态变量的先验分布的预测问题转化为一个求解后验分布的估计问题。采用高斯混合模型近似随机变量分布密度，从而将两种方法的计算结果在一个统一的预测框架之下进行有效交互，进一步提高了预测的准确性和可靠性。在算法的决策阶段，在获取的故障演化模型状态变量分布基础上，结合一定的故障判据近似计算出对象***剩余寿命分布。故障预测仿真实验结果证明了所提算法的有效性。

但上述两种预测方法有几点不足：

(1)基于模型的方法需要对研究对像进行数学建模，当故障传播的过程是复杂的或并没有完全理解的情况下，难以得出正确的故障状态模型。

(2)由于一阶马尔科夫形式利用了少量的历史信息，即下一时刻的状态值只和当前时刻状态值有关系，只能用于描述故障的生长过程。但状态的演化不仅与前一刻的状态有关，而且和之前的多个状态有关。所以应用一阶马尔科夫模型并不能真正的描述状态预示模型。

发明内容

为了克服现有技术的不足，本发明提供一种基于LS-SVR建模的高阶粒子滤波机械故障预示新方法，以有效提高故障预示精度。

本发明的技术方案为：

所述一种基于LS-SVR建模的高阶粒子滤波器的故障预示方法，其特征在于：包括以下步骤：

步骤1：对一组从正常状态演化到故障状态的信号通过滑动时间窗截取信号，截取的每一段信号作为一个样本，计算每一个样本的6种时域统计特征，时域统计特征包括峰值、平均幅值、均方根幅值、偏斜度指标和峭度指标；然后建立样本特征空间，并用局部线性嵌入方法对构建的样本特征空间进行低维特征提取，得到低维特征集；

步骤2：将步骤1中得到的低维特征集作为样本训练和测试预测器，其中低维特征集中的奇数点作为训练集，偶数点作为测试集，训练集是用来得到预测模型，测试集是用来测试预测模型的预测的效果；

步骤3：采用训练集训练最小二乘支持向量机，得到最小二乘支持向量机预测模型，将最小二乘支持向量机预测模型加上模型噪声组合成高阶隐马尔科夫模型：

$x_k={\hat{x}}_k+v_{k+1}$

${\hat{x}}_{k}f(x_{k-1},x_{k-2},......,x_{k-m})$

其中f()为最小二乘支持向量机预测模型，m为阶数，x_k为在k时刻的状态值，为在k时刻的状态估计值，x_k-1,x_k-2,......,x_k-m为k时刻之前1时刻到前m时刻的状态值，v_k-1为k-1时刻的模型噪声；

步骤4：根据步骤3得到的高阶隐马尔科夫模型产生粒子集粒子集中的粒子个数为N_s，第i粒子表示0到k时刻的状态集，粒子集中各粒子对应权值为权值由得到，其中由k时刻目标状态的后验概率密度得到，z_k表示k时刻的量测值；并将权值归一化为 $\underset{i}{\mathrm{\Σ}}{\mathrm{\ω}}_k^i=1;$

步骤5：根据步骤4得到的权值判定粒子退化程度，若粒子退化到达设定条件，则进行粒子重采样，并计算重采样后粒子权值；

步骤6：由粒子权值和粒子值得到状态估计值：

${\hat{x}}_k=\underset{i=1}{\overset{\mathrm{Ns}}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^i{x}_k^i$

再将状态估计值代入预测模型得到状态值；

步骤7：根据步骤4至步骤6进行下一时刻运算，并循环进行直至机械故障预示完成。

进一步的优选方案，所述一种基于LS-SVR建模的高阶粒子滤波器的故障预示方法，其特征在于：步骤6中，由粒子权值和粒子值通过以下多步状态预测得到状态估计值：

${\hat{x}}_k=\underset{i=1}{\overset{\mathrm{Ns}}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^i{x}_k^i$

$x_k^i=f(x_{k-1}^i,x_{k-2}^i,...,x_{k-m}^i).$

有益效果

本发明实现了一种机械故障预示的方法，该方法可以将高维特征进行降维，在低维特征空间中并应用m阶隐马尔科夫模型(HMM模型)训练状态演化模型，并结合粒子滤波器构成高阶粒子滤波器对非高斯、非线性状态进行预示，和现有预示方法相比较具有预示精度高，适用范围广的特点。

本发明之所以具有上述的有益效果其原因在于：用局部线性嵌入方法(LLE)降维，并不仅仅是纯粹的降维，更是一种有效地特征提取方式。它并非是通过从高维特征向量中选取部分分量而组成新的特征向量，而是高维特征向量在低维空间的一个映射。而通过这样一个映射，可以得到新的低维的特征向量，它既能保持原有数据的拓扑结构和局部关联性，而且新的向量的维数可以足够得低，大大减少了计算量。利用最小二乘支持向量机(LS-SVR)建立状态方程，相对于神经网络和SVR具有更强的泛化能力和更少的计算复杂度，更适宜应用在长时间的预测。在故障预示中，故障的状态演化不仅和前一时刻有关系，而且和前几时刻有关系，应用m阶HMM模型用来描述状态的演化模型能更好的反映***演化的真实情况。LS-SVR融入到粒子滤波算法框架中，设m＝3，构成3阶粒子滤波。和现有的方法相比较，其对非高斯、非线性状态进行预示具有更高的预示精度，适用范围更广。

附图说明

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

具体实施方式

下面结合附图和实施例对本发明进一步说明。

附图1展示了本发明所实现LS-SVR建模的高阶粒子滤波器的机械故障预示的总流程，该总流程图包含了实现最终预示所需的各个主要步骤。本发明的目的是通过提取从正常状态演化到故障状态的信号特征，用数据驱动的方式建立状态方程，并能够实现在线实时的预测。预测器通过程序来实现，对于提供的机械振动信号，提取相应的信号的特征，例如：时域的6种统计特性并应用LLE方法对其进行低维特征提取，然后利用LS-SVR构建m阶HMM模型，并用低维特征对其进行训练，最够将其融入到粒子滤波算法框架中构成m阶粒子滤波，最后对机械故障进行预示。

下面是具体的实现步骤

步骤1：对一组从正常状态演化到故障状态的信号通过滑动时间窗截取信号，截取的每一段信号作为一个样本，计算每一个样本的6种时域统计特征，时域统计特征包括峰值、平均幅值、均方根幅值、偏斜度指标和峭度指标；然后建立样本特征空间，并用局部线性嵌入方法对构建的样本特征空间进行低维特征提取，得到低维特征集。

局部线性嵌入方法(LLE)是本领域的公知方法，下面给出具体的LLE特征提取的过程：

原始特征空间是由6种统计特征向量组成的，分别有峰值、平均幅值、均方根幅值、偏斜度指标和峭度指标。用LLE算法提取6维特征空间中嵌入的d维特征。

LLE算法可以归结为三步：(1)寻找每个样本的k个近邻点；(2)由每个样本点的近邻点计算出该样本点的局部重建权值矩阵；(3)由该样本点的局部重建权值矩阵和其近邻点计算出该样本的输出值。

算法的第一步是计算出每个样本点的k个近邻点。把相对于所求样本点距离最近的k个样本点规定为所求样本点的k个近邻点。k是一个预先给定值。

LLE算法的第二步是计算出样本点的局部重建权值矩阵(W)。这里定义一个误差函数，如下所示：

$\min \mathrm{\ϵ}\left(W\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{|x_i-\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{w}_j^i{x}_{\mathrm{ij}}|}^2---\left(1\right)$

其中N表示样本个数，x_ij(j＝1,2,...,k)为x_i的k个近邻点，是x_i与x_ij之间的权值。且要满足条件：这里求取W矩阵，需要构造一个局部协方差矩阵Qⁱ。

$Q_{\mathrm{jm}}^i=(x_i-x_{\mathrm{ij}})(x_i-x_{\mathrm{im}})---\left(2\right)$

将上式与相结合，并采用拉格朗日乘子法，即可求出局部最优化重建权值矩阵：

$w_j^i=\frac{\underset{m=1}{\overset{k}{\mathrm{\Σ}}}{\left(Q^i\right)}_{\mathrm{jm}}^{-1}}{\underset{p=1}{\overset{k}{\mathrm{\Σ}}}\underset{q=1}{\overset{k}{\mathrm{\Σ}}}{\left(Q^i\right)}_{\mathrm{pq}}^{-1}}---\left(3\right)$

在实际运算中，Qⁱ可能是一个奇异矩阵，此时必须正则化Qⁱ，如下所示：

Qⁱ＝Qⁱ+rI  (4)其中r是正则化参数，I是一个k×k的单位矩阵。

LLE算法的最后一步是将所有的样本点映射到低维空间中。映射条件满足如下所示：

$\min \mathrm{\ϵ}\left(Y\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{|y_i-\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{w}_j^i{y}_{\mathrm{ij}}|}^2---\left(5\right)$

其中，ε(Y)为损失函数值，y_i是x_i的输出向量，y_ij(j＝1,2,...,k)是y_i的k个近邻点，且要满足两个条件，即：

$\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{y}_i=0,\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{y}_i{y}_i^T=I---\left(6\right)$

其中I是m×m的单位矩阵。这里的可以存储在N×N的稀疏矩阵W中，当x_j是x_i的近邻点时，否则，W_i,j＝0。则损失函数可重写为：

$\min \mathrm{\ϵ}\left(Y\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{M}_{i,j}{y}_i^T{y}_j---\left(7\right)$

其中M是一个N×N的对称矩阵,其表达式为:

M＝(I-W)^T(I-W)  (8)要使损失函数值达到最小,则取Y为M的最小m个非零特征值所对应的特征向量。在处理过程中，将M的特征值从小到大排列，第一个特征值几乎接近于零，那么舍去第一个特征值。通常取第2～m+1间的特征值所对应的特征向量作为输出结果。

在本发明中，嵌入为维数为d＝1，k值选取根据经验得到为k＝13。

步骤2：将步骤1中得到的低维特征集作为样本训练和测试预测器，其中低维特征集中的奇数点作为训练集，偶数点作为测试集，训练集是用来得到预测模型，测试集是用来测试预测模型的预测的效果。

步骤3：采用训练集训练最小二乘支持向量机，得到最小二乘支持向量机预测模型，将最小二乘支持向量机预测模型加上模型噪声组合成高阶隐马尔科夫模型：

$x_k={\hat{x}}_k+v_{k+1}$

${\hat{x}}_{k}f(x_{k-1},x_{k-2},......,x_{k-m})$

其中f()为最小二乘支持向量机预测模型，m为阶数，x_k为在k时刻的状态值，为在k时刻的状态估计值，x_k-1,x_k-2,......,x_k-m为k时刻之前1时刻到前m时刻的状态值，v_k-1为k-1时刻的模型噪声。

下面给出具体的LS-SVR的训练学习建立状态方程过程：

经典的支持向量机最终是解一个凸二次规划问题，其计算方法相对比较复杂，有较高的计算量。而这里所用到的LS-SVR是在经典SVR的基础上加入最小二乘方法，将二次规划问题转换为一个线性方程组的问题，简化其求解的过程，同时能够更新矩阵，使其在保证预测精度的情况下适用于动态预测。

小二乘支持向量机(LS-SVR)是比利时数学家Suykens等人提出对经典SVM算法的一种改进。他将标准的SVR算法中的不等式约束条件改为等式约束条件，以误差的平方项做为训练的损失函数，把一个凸二次规划问题转化为解一个二次线性方程组。对于线性的我们只要取

LSSVM回归问题可以表达为：

$\min J(\mathrm{\ω},\mathrm{\ξ})=\frac{1}{2}{\left|\right|\mathrm{\ω}\left|\right|}^2+\frac{1}{2}\mathrm{\γ}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{\mathrm{\ξ}}_i^2---\left(10\right)$

其中ω为权值向量,为松弛变量,γ是用于平衡拟合误差和模型复杂度的正则化参数,b为偏置。

对于优化问题，引入拉格朗日函数，最后得到以下方程

公式中a_i(i＝1,......,l)是拉格朗日乘子。由最优化条件，分别就ω,b,ξ,a，求偏导数，并设为0，可以得到：

$\frac{\∂L}{\∂b}=0\⇒\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{\mathrm{\α}}_i=0---\left(12\right)$

$\frac{\∂L}{\∂e_k}=0\⇒{\mathrm{\α}}_k={\mathrm{Ce}}_k,k=1,...,l$

为了求最优α和b，由KKT条件可以得到

$\left[\begin{array}{cc}0& I^T\\ I& \mathrm{\Ω}+{\mathrm{\γ}}^{-1}I\end{array}\right]\left[\begin{array}{c}b\\ \mathrm{\α}\end{array}\right]=\left[\begin{array}{c}0\\ y\end{array}\right]---\left(13\right)$

其中，y＝[y₁,......,y_l]；α＝[α₁,......,α_l]^l；I＝[1,......,1]_l；Ω为核矩阵，第i行j列的元素为最后将线性方程组(5)解出来的α和b，可以表出最小二乘回归函数：

$f\left(x\right)=\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{\mathrm{\α}}_{i}K(x,x_i)+b---\left(14\right)$

其中K(x,y)为核函数，

算法训练过程如下：

(1)输入N＝n个训练样本点，n是输入样本个数，对N个样本训练最小二乘支持向量机；

(2)对|α_i|(i＝1,.....,N)进行排序，剪枝M个|α_i|最小的训练样本(一般来说，M＝N×5％)；

(3)设N＝N-M，将剩余的N个训练样本组成新的训练样本；

(4)用最小二乘支持向量机训练削减过后的样本集；

(5)然后，回到第二步，反复迭代，直到最后的运算结果变差为止。

建立LS-SVR预测模型

1、确定核函数

多项式核函数：K(x,y)＝((x·y)+c)^d，其中c≥0，d是任意正整数。

高斯径向基核函数： $K(x_i,y)=\exp \left\{\frac{-{\left|\right|y-x_i\left|\right|}^2}{2{\mathrm{\σ}}^2}\right\}i=\mathrm{1,2},...,l.$

傅里叶核函数：其中，x,y∈R,0＜q＜1。

任何函数只要满足Mercer条件，都可以用来作为支持向量机的核函数，采用不同的函数为核函数，可以构造实现输入空间中不同类型的非线性决策面的学习机器。最常用的核函数是RBF高斯径向基核函数。

2、建立模型的方差，求解参数α和b，确定支持向量，并剪枝。

3、拟合曲线：根据求出的支持向量，拟合出趋势曲线，并通过误差分析得出最后的结论，也就是通过最小二乘法求解支持向量，然后“剪枝”一部分多余的支持向量，用剩下的支持向量去拟合成曲线在满足一定误差的情况下去做预测。

高阶隐马尔科夫(HMM)模型：

用某种适当的模型来描述一个机械振动***，对分析、研究该***非常重要。在实际工程应用中，经常采用动态空间模型来描述其中的问题。动态空间模型是一个很重要的统计分析工具，如Kalman滤波采用的高斯-马尔科夫线性模型就是一个很好的例子，它用***方程来描述状态随时间演变的过程，而测量方程来描述与状态有关的噪声变量。同样的，只要将高斯-马尔科夫线性模型写成一般的数学映射，就可以用两个方程来描述更一般的动态***。

机械故障预示问题是时变问题，分析一般动态状态空间模型分为状态转移模型p(x_k|x_k-1)和量测模型p(z_k|x_k),其中，代表***在时间k的状态变量(隐含变量或参数)，为***在时间k的测量值。对非线性、非高斯过程，其模型可表示为：

$\begin{array}{c}{x}_k=f(x_{k-1},v_{k-1})\\ z_k=h(x_k,n_k)\end{array}---\left(15\right)$

式中，和分别为过程噪声和量测噪声，并且是相互独立、协方差分别为Q_k和R_k的零均值加性噪声序列。

故障预示中，故障的状态演化不仅和前一时刻有关系，而且和前几时刻有关系，应用m阶HMM模型用来描述状态的演化模型能更好的反映***演化的真实情况。则需要m阶马尔科夫模型。m阶马尔科夫模型如下：

$\begin{array}{c}{x}_k={\hat{x}}_k+v_{k-1}\\ {\hat{x}}_k=f(x_{k-1},x_{k-2},......,x_{k-m})\end{array}---\left(16\right)$

x_k为在k时刻的状态值，为在k时刻的状态估计值，x_k-1,x_k-2,......,x_k-m为k时刻之前1时刻到前m时刻的状态值。

步骤4：根据步骤3得到的高阶隐马尔科夫模型产生粒子集粒子集中的粒子个数为N_s，第i粒子表示0到k时刻的状态集，粒子集中各粒子对应权值为权值由得到，其中由k时刻目标状态的后验概率密度得到，z_k表示k时刻的量测值；并将权值归一化为 $\underset{i}{\mathrm{\Σ}}{\mathrm{\ω}}_k^i=1.$

步骤5：根据步骤4得到的权值判定粒子退化程度，若粒子退化到达设定条件，则进行粒子重采样，并计算重采样后粒子权值。

步骤6：由粒子权值和粒子值得到状态估计值：

${\hat{x}}_k=\underset{i=1}{\overset{\mathrm{Ns}}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^i{x}_k^i$

再将状态估计值代入预测模型得到状态值。

步骤7：根据步骤4至步骤6进行下一时刻运算，并循环进行直至机械故障预示完成。

下面给出本实施例中基于LS-SVR建模的高阶粒子滤波预示算法的具体推导过程：

粒子滤波的基本思想是：首先依据***状态向量的经验条件分布在状态空间产生一组随机样本集合，称这些样本为粒子，然后根据量测并不断调整粒子的权重和位置，通过调整后的粒子信息修正最初的经验条件分布。其实质是利用粒子及其权重组成的离散随机测度近似相关的概率分布，并且根据算法递推更新离散随机测度。当样本容量很大时，这种蒙特卡罗描述就近似于状态变量真实的后验概率密度函数。这种技术适用于任何能用状态空间模型表述的非高斯背景的非线性随机***，精度可以逼近最优估计，是一种很有效的非线性滤波技术。

已知动态***的状态先验条件概率p(x₀),利用描述k时刻目标状态x_k的后验概率分布是对应权值为的粒子集，δ(·)是采样函数，第i粒子表示0到k时刻的状态集。k时刻目标状态的后验概率分布可离散地加权为：

$p\left(x_{0:k}\right|z_{1:k})\≈\underset{i=1}{\overset{N_s}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^i\mathrm{\δ}(x_{0:k}-x_{0:k}^i)---\left(17\right)$

其中，权值通过重要采样法确定：

${\mathrm{\ω}}_k^i\∝\frac{p\left(x_{0:k}^i\right|z_{1:k})}{q\left(x_{0:k}^i\right|z_{1:k})}---\left(18\right)$

对应的q(x_0:k|z_1:k)为重要密度函数，重要密度函数分解为：

$\begin{array}{c}q\left(x_{0:k}\right|z_{1:k})=q\left(x_k\right|x_{0:k-1},z_{1:k})q\left(x_{0:k-1}\right|z_{1:k-1})\\ =q\left(x_k\right|x_{k-m:k-1},z_{1:k})q\left(x_{k-m:k-1}\right|z_{1:k-1})\end{array}---\left(19\right)$

则通过由(19)得到粒子和由q(x_0:k-1|z_1:k-1)得到的粒子集可以得到新的粒子集

由于后验概率密度函数可以表述为：

$\begin{array}{c}p\left(x_{0:k}\right|z_{1:k})=\frac{p\left(z_k\right|x_{0:k},z_{1:k-1})p\left(x_{0:k}\right|z_{1:k-1})}{p\left(z_k\right|z_{1:k-1})}\\ =\frac{p\left(z_k\right|x_k)p\left(x_k\right|x_{k-m:k-1})}{p\left(z_k\right|z_{1:k-1})}p\left(x_{k-m:k-1}\right|z_{1:k-1})\\ \∝p\left(z_k\right|x_k)p\left(x_k\right|x_{k-m:k-1})p\left(x_{k-m:k-1}\right|z_{1:k-1})\end{array}---\left(20\right)$

则可得到重要性权值更新公式为：

${\mathrm{\ω}}_k^i\∝\frac{p\left(z_k\right|x_k^i)p\left(x_k^i\right|x_{k-m:k-1}^i)p\left(x_{k-m:k-1}^i\right|z_{1:k-1})}{q\left(x_k^i\right|x_{k-m:k-1}^i,z_{1:k})q\left(x_{k-m:k-1}^i\right|z_{1:k-1})}={\mathrm{\ω}}_{k-1}^i\frac{p\left(z_k\right|x_k^i)p\left(x_k^i\right|x_{k-m:k-1}^i)}{q\left(x_k^i\right|x_{k-m:k-1}^i,z_{1:k})}---\left(21\right)$

在高阶粒子滤波算法中选择最易于实现的先验概率密度作为重要密度函数，即

$q\left(x_k^i\right|x_{k-m:k-1}^i,z_{1:k})=p\left(x_k^i\right|x_{k-m:k-1}^i)---\left(22\right)$

将式(22)代入(21)，可将重要性权值简化为：

${\mathrm{\ω}}_k^i\∝{\mathrm{\ω}}_{k-1}^{i}p\left(z_k\right|x_k^i)---\left(23\right)$

将权值归一化，即

${\mathrm{\ω}}_k^i={\mathrm{\ω}}_k^i/\underset{i=1}{\overset{N_s}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^i---\left(24\right)$

粒子滤波最大的问题是粒子匮乏现象，即随着时间的增加，重要性权值有可能集中到少数粒子上，仅靠这些例子已不能有效表达出后验概率密度函数，为了避免这种退化现象GORDON提出了重采样方法。在重采样中，首先计算然后判断N_eff＜N_th，N_th为设定阈值，如果N_eff＜N_th，则粒子将被重采样，权值大的粒子被复制，权值小的被舍弃。粒子被选择复制的次数与它的权值大小成正比。重采样后的粒子拥有相等的权值1/N_s。经过这个步骤后得到的粒子代表了k时刻目标状态的后验分布。

而后验概率密度p(x_k|z_1:k)可表示为:

$p\left(x_k\right|z_{1:k})\≈\underset{i=1}{\overset{N_s}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^i\mathrm{\δ}(x_k-x_k^i)---\left(25\right)$

当N_s→∞时，有大数定理即可保证上式可逼近真实后验概率p(x_k|z_1:k)。

本实施例把LS-SVR预测模型整合到高阶粒子滤波框架中，取粒子滤波的阶数为m＝3时，描述故障生长趋势的3阶HMM可以表示为：

$\begin{array}{c}{x}_k={\hat{x}}_k+v_{k-1}\\ {\hat{x}}_k=f(x_{k-1},x_{k-2},x_{k-3})\end{array}---\left(26\right)$

其中，f(x_k-1,x_k-2,x_k-3)是由LS-SVR建立的预测模型，用当前时刻状态值，和前两个时刻的状态值去预测下一时刻状态值。如果预测p步状态，则对公式(26)进行递归运算，则可以得到p步的状态预测值。在高阶粒子滤波中每一状态的粒子数为100，对应每一个粒子的初始权值为0.01。

根据上述理论推导，本实施例可以归纳如下：

1、对一组从正常状态演化到故障状态的信号通过滑动时间窗截取信号然后建立样本的特征空间，并用局部线性嵌入方法(LLE)对构建的高维的样本特征空间进行低维特征提取，得到低维特征集。

2、用已有的状态数据训练LS-SVR建立机械故障演化模型。

3、用训练好的加上过程噪声的故障演化模型式(26)和3阶粒子滤波产生粒子集。根据公式(23)和当前的观测值计算每个粒子的权值，然后判定粒子退化程度是否需要重采样，最后根据公式(25)和(27)可以分别计算出概率密度函数和状态估计值。根据粒子值和当前的权值，一步状态预测可以用下式计算出：

${\hat{x}}_k=\underset{i=1}{\overset{\mathrm{Ns}}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^i{x}_k^i---\left(27\right)$

多步状态预测可以通过以下公式连续递推计算得到：

${\hat{x}}_k=\underset{i=1}{\overset{\mathrm{Ns}}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^i{x}_k^i$

$x_k^i=f(x_{k-1}^i,x_{k-2}^i,...,x_{k-m}^i).$

4、重复3直到机械故障预示完成。

评价指标:

为了评价本发明的预测性能和比较仿真实验的结果，我们使用以下指标来说明。

均方根误差(RMSE)：

M指数据点的总数量，y_i和分别表示i个实际值和预测值，RMSE的值越小，表明预测精度越高。

仿真条件：仿真实验中的滚动轴承数据来自美国宇航局(NASA)的数据库。该滚动轴承的数据是四个轴承由正常状态到故障状态的连续35天运转的测试值，其中最后3天被认为是故障状态。转速恒定为2000转，在轴和轴承上附加了6000磅的径向载荷。每个轴承的垂直方向和水平方向上各安装了一个振动加速度传感器，采样频率为20kHz。在本实验中，用该数据中第三个轴承的内圈故障数据做故障预示算法的验证。

仿真结果：如表1、表2所示。该表中的p是指向前预测p步,‘RNN’是递归神经网络p步向前预测的均方根误差,‘LS-SVR’是最小二乘支持向量机的p步向前预测的均方根误差,‘LS-SVR 1阶粒子滤波’是最小二乘支持向量机与一阶粒子滤波结合的p步向前预测的均方根误差。‘本发明’是本发明的方法的p步向前预测的均方根误差。

表1.轴承故障早期(前32天)状态预示实验数据表

表2.轴承故障状态(后3天)预示实验数据表

表1～表2中列出了递归神经网络、最小二乘支持向量回归和本发明得到的p步前向预测的均方根误差，同样实验条件下最高小的均方根误差加粗。由实验数据可知，本发明所采用基于LS-SVR建模的高阶粒子滤波方法进行机械故障预示比使用传统的预示方法具有更高的精度。特别是在轴承故障状态的5步前向预示时，均方根误差降低了0.102。

Claims (2)

1.一种基于LS-SVR建模的高阶粒子滤波器的故障预示方法，其特征在于：包括以下

步骤：

步骤1：对一组从正常状态演化到故障状态的信号通过滑动时间窗截取信号，截取的每一段信号作为一个样本，计算每一个样本的6种时域统计特征，时域统计特征包括峰值、平均幅值、均方根幅值、偏斜度指标和峭度指标；然后建立样本特征空间，并用局部线性嵌入方法对构建的样本特征空间进行低维特征提取，得到低维特征集；

步骤2：将步骤1中得到的低维特征集作为样本训练和测试预测器，其中低维特征集中的奇数点作为训练集，偶数点作为测试集，训练集是用来得到预测模型，测试集是用来测试预测模型的预测的效果；

步骤3：采用训练集训练最小二乘支持向量机，得到最小二乘支持向量机预测模型，将最小二乘支持向量机预测模型加上模型噪声组合成高阶隐马尔科夫模型：

$x_k={\hat{x}}_k+v_{k-1}$

${\hat{x}}_k=f(x_{k-1},x_{k-2},......,x_{k-m})$

其中f()为最小二乘支持向量机预测模型，m为阶数，x_k为在k时刻的状态值，为在k时刻的状态估计值，x_k-1,x_k-2,......,x_k-m为k时刻之前1时刻到前m时刻的状态值，v_k-1为k-1时刻的模型噪声；

步骤4：根据步骤3得到的高阶隐马尔科夫模型产生粒子集粒子集中的粒子个数为N_s，第i粒子表示0到k时刻的状态集，粒子集中各粒子对应权值为权值由得到，其中由k时刻目标状态的后验概率密度得到，z_k表示k时刻的量测值；并将权值归一化为 $\underset{i}{\mathrm{\Σ}}{\mathrm{\ω}}_k^i=1;$

步骤5：根据步骤4得到的权值判定粒子退化程度，若粒子退化到达设定条件，则进行粒子重采样，并计算重采样后粒子权值；

步骤6：由粒子权值和粒子值得到状态估计值：

${\hat{x}}_k=\underset{i=1}{\overset{\mathrm{Ns}}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^i{x}_k^i$

再将状态估计值代入预测模型得到状态值；

步骤7：根据步骤4至步骤6进行下一时刻运算，并循环进行直至机械故障预示完成。

2.根据权利要求1所述一种基于LS-SVR建模的高阶粒子滤波器的故障预示方法，其特征在于：步骤6中，由粒子权值和粒子值通过以下多步状态预测得到状态估计值：

${\hat{x}}_k=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^i{x}_k^i$

$x_k^i=f(x_{k-1}^i,x_{k-2}^i,...,x_{k-m}^i).$