基于改进多项式响应面的VPMCD方法及其在滚动轴承故障诊断中的应用

杨 宇，潘海洋，李 杰，程军圣

（湖南大学 汽车车身先进设计制造国家重点实验室，长沙 410082）

滚动轴承的故障诊断本质上是一个模式识别的过程。在滚动轴承的故障模式识别中，神经网络、支持向量机等模式识别方法得到了广泛的应用。但是它们都存在着一定的局限性，如人工神经网络具有收敛速度慢、网络结构难以确定［1］等问题，支持向量机具有核函数及核参数难以确定［2］等问题。除此之外，神经网络和支持向量机在进行模式识别时都忽略了从原始数据中所提取的特征值之间的相互内在关系。

然而，在机械故障诊断中，所有或部分特征值之间大都具有一定的内在关系，而且这种内在关系在不同的系统或类别（相同的系统在不同的工作状态下）间具有明显的不同。因此，可以对各个特征值之间的相互内在关系建立数学模型，对于不同的类别可以得到不同的数学模型，从而可以采用这些数学模型对被测试样本的特征值进行预测，把预测结果作为分类的依据，进一步进行模式识别。基于此，Raghuraj与Lakshminarayanan提出了一种新的模式识别方法——基于变量预测模型的模式识别（Variable Predictive Model Based Class Discriminate，简称VPMCD）方法，同时还将该方法与神经网络、支持向量机等其它模式识别方法进行了对比，结果验证了VPMCD方法的有效性和优越性［3］。

在VPMCD方法中，模型的训练过程也就是变量预测模型的建立过程，在建模过程中，变量预测模型的拟合精度影响着整个模式识别的效果。因此，获得合理的变量预测模型至关重要。VPMCD是通过特征值之间存在相互内在关系来建立变量预测模型的，但是这种相互关系的具体情况却难以确定，而且特征值之间相互内在关系的实际预测模型也无法得到。因此，该方法选取了线性（L）模型、线性交互（LI）模型、二次（Q）模型、二次交互（QI）模型四种模型作为实际模型的代理模型，并从中选取最佳代理模型，这四种模型都属于多项式响应面（Polynomial Response Surface，简称PRS）模型［4］。多项式响应面模型是一种常用的响应面（Response Surface Method，简称 RSM）［5］模型，也是应用最广泛的近似模型，该方法计算简单，但其高阶计算量却很大，且提高的精度有限。另外，响应面法不能够随着样本容量的增大而有效提高其近似精度，这两点严重限制了该方法的使用［6］。究其原因，主要是PRS方法采用最小二乘拟合参数，而最小二乘拟合与真实数据往往存在偏差，并且放弃这些残差项，从而导致模型近似过程中忽略了残差项，而大多数情况下残差项包含了很多重要信息，将其忽略容易使模型拟合误差偏大。针对这一缺陷，本文将原方法中的多项式响应面法进行了改进，改进的多项式响应面法对PRS法忽略掉的残差重新进行近似，保留了模型拟合的重要信息，因而比原方法的模型拟合效果更理想。基于此，本文提出了基于改进多项式响应面（Improved Polynomial Response Surface，简称IPRS）的VPMCD方法，并将其运用于滚动轴承故障诊断。

1 基于改进多项式响应面的VPMCD方法

1.1 VPMCD方法

以机械故障诊断问题为例，采用p个不同的特征值 X＝［X1，X2，…，Xp］来描述一个故障类别，对于其中的特征值Xi来说，当故障类别不同时，其他的一个或者多个特征值对其影响也会发生变化。因此，特征值Xi与其余的一个或者多个特征值之间存在着一定的函数关系，而这种关系可以是线性的，也可以是非线性的。为了识别滚动轴承的故障模式，需要有能够描述这些函数关系的数学模型，以便对测试样本的特征值进行预测，进一步对测试样本进行分类，这种模型称为变量预测模型（Variable predictive model，简称 VPM）。

为特征值Xi定义的变量预测模型是一个线性或非线性的回归模型，可以选择以下四种模型之一：

① 线性模型（L）：

② 线性交互模型（LI）：

③ 二次交互模型（QI）：

④ 二次模型（Q）：

式中：r≤p－1为模型阶数。这四种模型都属于多项式响应面（PRS）模型，以p个特征值为例，选取四种模型中任意一个模型，用特征值Xj（j≠i）对Xi进行预测，都可以得到：

式（5）称为变量Xi的变量预测模型VPMi。其中，特征值Xi称为被预测变量；Xj（j≠i）称为预测变量；e为预测误差；b0，bj，bjj，bjk为模型参数，可以通过所有训练样本的特征值对它们进行参数估计。

对于g类故障分类问题，提取p个特征值X＝［X1，X2，…，Xp］．在分别采用不同故障类别下的训练样本数据对预测模型进行训练后，不同故障类别下的不同特征值就可以分别建立g×p个预测模型VPMki，其中k＝1，2，…，g代表不同的类别，i＝1，2，…，p代表不同的特征值。然后针对测试样本，提取其特征值，并用g×p个模型VPMki分别对它们进行预测，得到g×p个预测值X～i，以同一类别下所有特征值的预测误差平方和最小为判别函数，对测试样本的故障类型和工作状态进行分类，该方法称为基于变量预测模型的模式识别（VPMCD）方法。

1.2 多项式响应面法

VPMCD在模型训练过程中是从四种多项式响应面模型中选择最佳模型为变量预测模型，在进行模型近似时主要是利用最小二乘拟合的方法，通过对数据的分析，找出相关因素与响应量之间的函数关系，并用这种函数代替原来的函数关系。PRS法通常是高阶多项式，基本形式如下：

式中：xi是 m维自变量 x的第 i个分量；β0，βi和 βij是未知参数，将它们按顺序构成列向量β，求解PRS模型最重要的就是求解向量β。将样本点的值代入式（6）中，再利用最小二乘法可以求得其估计值，即 β＝

PRS方法采用的是泰勒展开式的思想，对于函数f（x），当采用常用的二次多项式时，其自变量x在x0处展开后得到：

式中：Δf（x0）为函数 f在 x0点的梯度与 Δx的点积；Hf（x0）为函数f在 x0点的 Hessian矩阵。PRS法根据样本数据利用最小二乘法计算函数f在某个定义域区间内的梯度及Hessian矩阵，并将项o（‖x‖3）作为残差而将其忽略。但是，该项往往会比较大，而且包含很多重要信息，忽略掉该项会使近似值和实际值之间存在较大差异。另外，由式（6）可知，如果一个PRS模型具有n个变量，阶数为 m，则相应的未知参数为未知参数的数目会随着阶数的增加而增加，所需的计算量也会相应增加，而模型的计算精度却没有明显提高。因此，本文在分析PRS模型缺陷的基础上，提出了改进的多项式响应面（Improved Polynomial Response Surface，简称 IPRS）模型。

1．3 改进的多项式响应面法

相对于PRS方法，IPRS方法在近似过程中并没有忽略残差项，而是对残差重新进行近似。在PRS方法中，设生成的近似曲面为f1，近似曲面在采样点处和已知结果f存在一定的差别，记为残差R。IPRS方法则是进一步对残差进行近似估计，并把残差的估计结果加入到近似曲面中去。在对残差进行插值估计的过程中，选择合适的插值估计方法至关重要，在插值估计方法中，径向基神经网络是一种性能良好的前向网络，具有最佳逼近，能克服局部极小值问题的性能［7－8］，用径向基神经网络对残差进行插值处理可以达到较好的插值效果。

基于此，改进的多项式响应面法运用径向基神经网络对残差进行插值处理，具体的步骤如下：

（1）对初始采样点生成PRS近似曲面f1；

（2）计算近似曲面f1和实际曲面之间的残差，设实际曲面为Y，则残差R＝Y－f1；

（3）利用残差矩阵R进行径向基神经网络插值估计，得到插值函数f2；

（4）利用曲面f1和曲面f2分别进行插值计算，得到插值结果y1和y2，并叠加y1和y2，作为最终的插值结果；

计算流程可以参考图1，这样，可以得到IPRS方法在n维向量x点处的插值结果fa（x）：

1．4 仿真对比分析

为了验证IPRS法的拟合效果，利用如下的测试函数对其进行分析：

其中 0．5≤x1，x2≤3．5。

图1 改进响应面法计算流程图Fig．1 The calculation flow chart of improved polynomial response surface method

试验在定义域内随机对自变量进行取值，为了验证训练样本容量对拟合结果的影响，分别随机产生30，100，200组训练样本，同时在区间内按照线性分布进行插值，插值点数选择30个，分别采用PRS法和IPRS法对该函数在定义域区间内进行拟合。PRS方法中，与式（8）最接近的模型为二次交互（QI）模型，因此选择QI模型。在用IPRS法拟合时，需利用径向基神经网络对残差进行处理，本文选择创建的径向基网络，其径向基神经元数目等于输入样本数目，创建该网络只需要进行一次运算，因而能够达到较快的训练速度。然而，创建网络时需要合理选择散布常数spread，散布常数控制着网络输出的光滑情况，其值越大输出结果越光滑，但是太大的结果又会造成计算上的困难，经过多次试验，选择散布常数spread＝1。用两种方法拟合后，可以得到各自的近似曲面，图2～图4分别是实际图形以及两种方法分别在训练样本为30、100、200组情况下的插值拟合三维图。

从三图中可以明显看出，在不同训练样本容量的情况下，采用IPRS法拟合得到的近似模型相对于采用PRS法拟合得到的近似模型要更加接近于实际模型，随着采样样本数量的增加，PRS法的拟合效果没有显著提高，而IPRS法的拟合效果有了明显的提高。

对于两种方法的计算结果，可以分别通过以下三种方法来验证它们的拟合精度：

（1） 预测点的均方根差：MSE＝E［（y＾－y）2］．

（2）为了评价整个模型的精度，可以选择经验积累方差标准：其中 m为误差取样点总数，在这里m＝30×30；

图2 训练样本为30组情况下的拟合图形和原始图形Fig．2 The fitting graphics and original graphics at the training sample group of 30

图3 训练样本为100组情况下的拟合图形和原始图形Fig．3 The fitting graphics and original graphics at the training sample group of 100

图4 训练样本为200组情况下的拟合图形和原始图形Fig．4 The fitting graphics and original graphics at the training sample group of 200

两种方法在训练样本数目不同的情况下的拟合精度对比如表1所示。从表中可以看出，PRS方法在各验证方法下的拟合精度都比IPRS方法低，而且随着训练样本数目的增加，其精度提高的也不明显；相反，IPRS的拟合精度随着训练样本数量的增加有着显著地提升，这说明IPRS方法在训练样本数目较多的情况下比PRS方法的拟合能力更好。因此，将IPRS法代替PRS法而引入VPMCD算法可以达到更理想的模式分类效果。

表1 两种方法在训练样本容量不同的情况下的拟合精度对比Tab．1The contrast of fitting precision of the two methods under different capacity of training sample

1．5 基于改进多项式响应面的VPMCD方法

本文将改进的多项式响应面法代替原始VPMCD中的多项式响应面法，提出了基于改进多项式响应面的VPMCD方法（以下简称IPRS－VPMCD），具体的步骤如下所示：

（1）模型训练过程：

① 对于g类状态分类问题，共收集n个训练样本，每一类状态样本数分别为 n1，n2，…，ng。

② 对所有训练样本提取特征值X＝［X1，X2，…，Xp］。

③ 对任意被预测变量 Xi，选择代理模型（L、LI、QI、Q四种PRS模型之一）、预测变量和模型阶数。对于不同的特征值，其预测模型类型、预测变量和模型阶数都有可能不同。

④ 令k＝1，对于nk个第k类训练样本中的任意一个样本，分别对每一个特征值Xi建立PRS模型，则对每一个特征值可以建立nk个方程；然后用PRS法获得Xi的插值结果，再以预测误差平方和为判别函数，得到预测误差平方和最小的PRS模型；再利用此模型，得到被预测变量的拟合值和实际值之间的残差Ri，利用径向基神经网络对其进行插值估计，得到插值结果和对应的插值模型，并将其与之前的PRS模型相叠加，所得到的模型即为特征值Xi的预测模型

⑤ 令k＝k＋1，循环步骤④，至k＝g结束。

⑥至此，对所有模型类别下的所有特征值都分别建立了预测模型VPMik，其中 k＝1，2，…，g代表不同类别，i＝1，2，…，p代表不同特征值。

（2）模型分类过程：

① 选择测试样本，并提取其特征值 X＝［X1，X2，…，Xp］。

② 对于测试样本的所有特征值 Xi（i＝1，2，…，p），分别采用 VPMki（k＝1，2，…，g）对其进行预测，得到测试值，其中 k＝1，2，…，g代表不同类别，i＝1，2，…，p代表不同特征值。

③ 计算同一类别下所有特征值的预测误差平方最小为判别函数对测试样本进行分类，其中k＝1，2，…，g代表不同类别。当在g个预测误差平方和值中小时，将测试样本识别为第k类。

2 基于LCD近似熵和IPRS-VPMCD的滚动轴承故障诊断方法

在滚动轴承故障振动信号特征提取算法中，近似熵在描述信号的复杂性时具有较好的抗噪、抗干扰能力，而且利用较短数据即可以较稳健地估计出信号的近似熵。同时，近似熵能用于随机过程和确定性过程，其取值大小会随着随机过程和确定过程的混合比例不同而不同［9］。因此，近似熵能表征信号的复杂程度和产生新模式的概率，可将其应用于故障诊断领域。而且，文献［10］指出，在实际计算中，数据长度n为有限值，当近似熵的嵌入维数m＝2，相似容量r＝0．1SDx～0．2SDx（SDx为原始数据 x（i）的标准差）时，熵值对 N的依赖程度最小，具有较合理的统计特性。因此，可以通过近似熵算法对滚动轴承振动信号进行特征提取。

然而滚动轴承振动信号往往表现出非平稳性，若直接进行近似熵计算会影响诊断精度，因此，必须先对原始振动信号进行处理。在非平稳信号的分析中，时频分析方法能同时提供非平稳信号在时域和频域的局部化信息，因而得到了广泛的应用。在时频分析方法中，局部特征尺度分解（Local Characteristic Scale Decomposition，简称LCD）算法［10］是一种新的基于极值点的局部特征尺度参数的自适应、非平稳信号处理方法，它能将信号自适应地分解为一系列瞬时频率具有物理意义的内禀尺度分量（Intrinsic Scale Component，简称ISC）。LCD在分解速度、信号的还原度方面都能够达到较满意的效果，因而具有非常大的应用前景。本文将近似熵算法和LCD算法相结合应用于滚动轴承故障诊断，通过LCD方法将滚动轴承振动信号分解为若干个平稳的ISC分量，计算每个ISC分量的近似熵，再利用不同ISC分量中提取的近似熵之间存在相互关系这一特点，采用IPRS-VPMCD方法建立预测模型，从而进行模式分类。

基于LCD近似熵和IPRS-VPMCD的滚动轴承故障诊断方法的具体步骤为：

（1）在一定转速下以采样率fs对滚动轴承正常、内圈故障、外圈故障和滚动体故障四种状态进行采样，每种状态采集N组样本。

（2）分别对各类别状态的原始信号进行LCD分解，每个信号得到若干个ISC分量和余量。选择合适的i个ISC分量，对所选的ISC分量分别提取近似熵作为特征值，组成特征值向量，每种状态下得到N×i阶的特征值矩阵。

（3）对于不同状态下的N组特征值，选择n组样本作为训练样本，采用IPRS-VPMCD方法进行训练得到预测模型；

（4）将剩余的N－n组样本作为测试样本，用训练得到的预测模型对其进行预测并分类，根据 IPRSVPMCD分类器的输出结果来确定滚动轴承的工作状态和故障类型。

3 应 用

本文将IPRS-VPMCD方法应用于滚动轴承故障诊断，采用美国西储大学电气工程实验室的滚动轴承实验数据来对该方法的有效性和优越性进行验证，所采用的轴承型号、参数和实验装置见文献［11］。采样频率为48 kHz，电机负载为 0．746 kW，转速为 1 772 r／min，状态类型分别为：正常状态、外圈故障、内圈故障、滚动体故障。故障点的直径为0．177 8 mm，故障深度为0．279 4 mm，每种状态各得到200个样本。对各样本的原始信号进行LCD分解，选择标准偏差法［12］作为终止判据，选择镜像对称延拓方法［13］减少边界效应，图5所示的是某一滚动轴承滚动体故障状态下的信号及其LCD分解后的分量。

由于滚动轴承的故障信息主要集中在高频段，因此，可选取前四个ISC分量，并对各分量求取近似熵值（算法中选择 m＝2，r＝0．2SDx），分别标记为 X1，X2，X3，X4。将所得的近似熵值组成特征向量，以此作为分类器的输入进行模式识别。各类数据可以提取出200组近似熵，每组近似熵为4个。各类别状态选择50组数据作为测试样本，训练样本分别选择30、90、150组，用IPRS-VPMCD方法进行训练和测试，得到不同训练样本容量下的分类结果。同时，将IPRS-VPMCD方法与VPMCD方法在相同情况下的分类结果进行对比分析，以此验证两种模式识别方法在不同训练样本容量下的模式分类能力。

图5 滚动体故障状态下的振动信号及其LCD分解后的分量Fig．5Vibration signal of rolling element bearing fault state and its components after LCD decomposition

在IPRS-VPMCD方法的残差估计中，选择径向基神经网络时，网络输入中散布常数spread的取值不同会对分类结果有一定的影响，需合理选择，本实验参照交叉验证［14］的方法，并通过多次试验，选择 spread＝1．1。另外，为了验证IPRS-VPMCD方法的适用性及优越性，分别与RBF神经网络、支持向量机和VPMCD比较，首先经过优化选择，设置RBF神经网络训练误差的平方和为0．05，设置支持向量机的折衷系数为10，核函数为RBF核函数，然后用四种分类器分别进行训练分类，最后得到不同训练样本数目情况下的分类结果。具体情况如表2所示。

表2 四种分类方法的分类时间和识别率对比Tab．2 The comparison of four kinds of classification methods in the classification time and recognition rate

从表中可以看出，随着样本容量的增大，四种方法的分类精度都明显有所提高，提高了3%左右，但是无论样本容量多少，RBF神经网络和VPMCD的分类精度都明显低于支持向量机和IPRS-VPMCD。另外，再从分类时间上看，VPMCD和IPRS-VPMCD的分类时间明显快于另外两种，尤其是支持向量机，综合分析，IPRSVPMCD方法在训练过程中除了获得最佳PRS模型之外，对其残差还进行了插值处理，保留了部分重要信息，因此，IPRS-VPMCD方法在保证分类效率的同时，有效地提高了计算精度，并且分类精度能够随训练样本容量增加而明显提高。

4 结 论

改进的多项式响应面法将多项式相应面法中忽略掉的残差项重新进行了插值近似，保留了模型拟合过程中的重要信息，因而能够获得更加合理的近似模型。通过仿真分析，将多项式响应面法和改进的多项式响应面法进行拟合效果的比较，结果证明了改进多项式响应面法具有更好的拟合效果，而且能够随着训练样本容量的增加较明显地提高拟合精度。

将改进的多项式响应面法代替原始VPMCD方法中的多项式响应面法，提出了IPRS-VPMCD方法，并将其应用于滚动轴承故障诊断中。通过实验，将IPRSVPMCD方法和VPMCD方法在训练样本容量不同情况下的分类结果进行对比分析，结果证明了它的优越性。

［1］Wang C C，Kang Y，Shen P C，et al．Applications of fault diagnosis in rotating machinery by using time series analysis with neural network［J］．Expert System with Application，2010，37（2）：1696－1702．

［2］Fei Sheng-wei，Zhang Xiao-bin．Fault diagnosis of power transformer based on support vector machine with genetic algorithm［J］．Expert Systems with Application，2009，36（8）：11352－11357．

［3］Raghuraj R， Lakshminarayanan S． Variable predictive models-A new multivariate classification approach for pattern recognition applications［J］．Pattern recognition，2009，42（1）：7－16．

［4］穆雪峰，姚卫星，余雄庆，等．多学科设计优化中常用代理模型的研究［J］．计算力学学报，2005，36（22）：5－8．MU Xue-feng，YAO Wei-xing，YU Xiong-qing，et al．A survey of surrogate models used in MDO ［J］．Chinese Journal of Computational Mechanics，2005，36（22）：5－8．

［5］Joyce A P，Leung SS．Use of response surface methods and path of steepest ascent to optimize ligand-binding assay sensitivity［J］．Journal of Immunological Methods，2013，392（1－2）：12－23．

［6］潘雷，谷良贤，阎代维．改进响应面法及其近似性能研究［J］．宇航学报，2009，30（2）：806－810．PAN Lei， GU Liang-xian， YAN Dai-wei． Research on improved response surface method and its performance［J］．Journal of Astronautics，2009，30（2）：806－810．

［7］彭磊，刘莉，龙腾．基于动态径向基函数代理模型的优化策略［J］．机械工程学报，2011，47（7）：164－170．PENG Lei，LIU Li，LONG Teng．Optimization strategy using dynamic radial basis function metamodel［J］．Journal of Mechanical Engineering，2011，47（7）：164－170．

［8］吴月明，王益群，李莉．BP神经网络与广义RBF神经网络在产品寿命分布模型识别中的应用研究［J］．中国机械工程，2006，17（20）：2140－2143．WU Yue-ming，WANG Yi-qun，LI Li．Application study on BP network and generalized RBF network in estimating distribution model of mechanical products［J］． China mechanical engineering，2006，17（20）：2140－2143．

［9］Lee JH，Kim J，Kim H J．Development of enhanced wignerville distribution function［J］．Mechanical Systems and Signal Processing，2001，13（2）：367－398．

［10］Classen T，Mecklenbrauker W．The aliasing problem in discrete-time wigner distribution［J］．IEEE Transactions on Acoustics，Speech，and Signal Processing，1983，31（5）：1067－1072．

［11］Case Western Reserve University Bearing Data Center．Bearing Data Center Fault Test Data．［EB／OL］．［2009－10－01］． http： ／／www．eecs． case． edu／laboratory／bearing．

［12］Huang N E，Shen Z，Long S R．The empirical mode decomposition and the Hilbert spectrum for nonlinear nonstationary time Series Analysis［J］． Proc． Roy． Soc．London，1998，454：903－995．

［13］Rilling G，Flandrin P，Goncalves P．On empirical mode decomposition and its algorithms［J］．In IEEE-EURASIP Workshop on Nonlinear Signaland Image Processing，Grado，Italy，2003：8－11．

［14］Breiman L．Heuristics of instability and stabilization in model selection［J］．The Annals of Statistics，1996，24（4）：2350－2383．