基于关联及多阶拟合的螺杆泵故障诊断

Fault Diagnosis of Screw Pump Based on Correlation and Multi-Order Fitting

摘要: 螺杆泵在不同中心频率下的振动加速度偏离正常范围会导致螺杆泵产生异常振动，因此，螺杆泵在出厂前需要设置不同测点检测其振动加速度。为此，本文基于螺杆泵在正常工作下不同测点的振动加速度应用Apriori算法得到螺杆泵不同测点间振动加速度并无显著性差异且呈强正相关的结论；然后应用多阶拟合法建立螺杆泵中心频率和振动加速度的函数关系式，以正常运行参数所对应的系数向量作为模板，计算出待诊断设备的测点运行参数所构造的系数向量与模板之间的系数距离，并应用统计方法设计适当的阈值与之作比较，当且仅当系数距离大于阈值时，将待诊断样本判定为故障。经实际应用，本文方法能够及时发现螺杆泵的异常振动，从而达到故障诊断的目的。

Abstract: The vibration acceleration of the screw pump at different center frequencies deviates from the normal range, which will cause abnormal vibration of the screw pump. Therefore, the screw pump needs to set different measurement points to detect its vibration acceleration before leaving the factory. For this reason, based on the vibration acceleration of different measurement points of the screw pump under normal operation, this paper applies the Apriori algorithm to obtain the conclu-sion that there is no significant difference in vibration acceleration between different measure-ment points of the screw pump and that there is a strong positive correlation. The multi-order fit-ting method is used to establish the functional relationship between the center frequency and the vibration acceleration of the screw pump. The coefficient vector corresponding to the normal op-erating parameters is used as a template, and calculate the coefficient distance between the coeffi-cient vector constructed by the operating parameters of the measurement points of the equipment to be diagnosed and the template, then apply statistical methods to design an appropriate thresh-old to compare with it, and determine the sample to be diagnosed as a fault if and only if the coeffi-cient distance is bigger than the threshold. After practical application, the method in this paper can find the abnormal vibration of the screw pump in time, so as to achieve the purpose of fault diagno-sis.

1. 引言

螺杆泵作为一种排水及动力装置，因其具有平衡性好、启动时间短和振动噪声小等优点，常常作为辅机，装备在船舶、航母及潜艇等大型军事设备中。因此，螺杆泵能正常、安全、高效地运行不仅关系到泵器件本身，更关乎舰船的运行安全 [1]。不仅如此，诸如潜艇此类军事设备，为了保证作战时的隐蔽性，必须要求舰艇内部设备及机械装置的噪声、振动、声波严格控制在安全范围中。螺杆泵的异常振动会极大地影响舰艇的隐蔽性，从而对整个舰艇和舰艇上的人员的安全造成威胁。因此，维持螺杆泵的正常运行，检测螺杆泵存在的异常振动具有实际意义，基于螺杆泵的振动指标对螺杆泵进行异常检测成为近些年的研究热点。

之前已经有许多学者提出了异常检测所存在的问题和解决方法。由于采集设备故障数据的成本较高且难度较大，因此大部分的实验数据都是在设备正常工作下检测得到的，Zaher等人 [2] 将缺乏故障记录描述为早期故障检测的主要挑战，从而强调了正常行为建模的重要性。Shen Zhang等人 [3] 提出了一种使用基于变分自动编码器(VAE)的深度生成模型进行轴承异常检测的半监督学习方法，当只有一小部分数据带有标签时，该方法可以有效利用数据集。Sreelekha Guggilam等人 [4] 提供了一种流式聚类和异常检测算法，该算法在执行概率异常检测和聚类时不需要对异常分数或簇数有严格的阈值，这确保了群集的形成不会受到异常数据的存在的影响，从而能够更可靠地定义“正常与异常”行为。

针对上述背景，本文基于螺杆泵在正常工作下的特征参数(不同中心频率下的振动加速度)，深入研究螺杆泵的振动异常检测方法。首先根据企业方提供的螺杆泵运行数据应用K-Means进行离散化处理，然后利用Apriori算法深入挖掘正常工作下不同测点振动加速度间的潜在关系。并在此基础上，基于螺杆泵中心频率和振动加速度的函数关系应用多阶拟合法，求解多项式拟合系数，以正常运行参数所对应的系数向量作为模板，计算出待诊断设备的测点运行参数所构造的系数向量与模板之间的系数距离，根据最小距离原则，距离越小表明待诊断样本与模板越相似。为进一步量化诊断规则，本文通过统计方法，设计并计算适当的阈值，将其与系数距离作比较，从而得出故障诊断结果。最后，为了验证方法的有效性和可行性，本文通过模拟仿真试验，收集故障数据并构建测试数据集，对测试数据集中的样本应用多阶拟合法并计算相应的系数距离，根据阈值得出诊断结果。另外，计算相应的评估指标(正确率、精度及召回率等)，从而对实验结果进行分析。

2. 螺杆泵实验数据

本文分析的数据源是企业对正常工作的泵进行检测所得到的实验数据。针对同一型号的泵，按GJB4058-2000《J船设备噪声、振动测量方法》中的相关规定，对机脚振动加速度进行振动测量，泵组机脚振动加速度测点位置如下图1所示。

Figure 1. Location of vibration acceleration measurement point of pump set feet

图1. 泵组机脚振动加速度测点位置

企业方已经应用傅里叶变换将测量数据从时域变换到频域，因此，本文的数据源字段主要包括中心频率和各个测点的振动加速度，描述的是不同中心频率所对应的各个测点的振动加速度，展示部分数据参见表1。

Table 1. Some data sources

表1. 部分数据源

3. 基于Apriori算法的关联规则分析

由表1可以看出，测点之间的差异并不显著，因此本文进一步挖掘各测点之间是否存在互相影响的关系。首先应用K-Means对实验数据进行离散化处理，然后基于Apriori算法 [5] 和相关系数探究测点之间的潜在关系。一方面能够避免相似测点对诊断结果的影响，另一方面能够简化后续的模型求解过程。

3.1. 基于K-Means的数据离散化

由于Apriori算法要求离散的数据格式，而实验数据是连续数值型数据，故需要先对数据进行离散化处理。本文采用基于K-means的一维聚类离散化方法，先确定最佳K值，再将连续属性值划分为K个区间，让簇内的点差值小，而让簇间点差值大。

本文采用肘部法则(Elbow Method) [6] 确定K值。此方法适用于K值相对较小的情况，将每个簇的质点与簇内样本点的平方距离误差和称为畸变程度(Distortion)，当选择的k值小于真正的K时，k每增加1，Distortion值就会大幅的减小；当选择的k值大于真正的K时，k每增加1，Distortion值的变化就不会那么明显。这样，正确的K值就会在这个转折点，类似手肘的地方。

如图2所示，通过绘制k与Distortion的关系曲线图确定最佳K值，肘部的值(Distortion开始时下降很快，在肘部开始平缓了)作为K值，可以看出K应取3。对每个测点都应用肘部法求K值，得到各个测点的最优K值均为3。

Figure 2. Curve: applying measurement point 1 to the elbow method

图2. 测点1应用肘部法

根据每个测点的最优K值，应用K-means对各个测点数据进行聚类，从而达到离散化的目的。K-means的基本思想是先从样本集中随机选取K个样本作为簇中心，并计算所有样本与这K个“簇中心”的距离，对于每一个样本，将其划分到与其距离最近的“簇中心”所在的簇中，对于新的簇计算各个簇的新的“簇中心”。直到“簇中心”不再变化或达到300次迭代为止。

将离散化结果中的每个测点数据标记为K个振动加速度级，例如测点1的K值为3，则将其离散化结果分别标记为“一级”，“二级”和“三级”。

3.2. 关联规则分析

为挖掘不同测点之间的关系，基于离散化后的数据，本文采用Apriori算法分析不同测点的振动加速度级之间的关联规则。该算法的基本思想是：首先利用递归的方法 [7] 找出所有的频繁项集(诸如{测点1 = “二级”，测点2 = “二级”}等)，然后由频繁项集产生强关联规则(诸如测点2 = “二级” à 测点1 = “二级”等)，将大于给定最小置信度的规则保留下来，具体的算法流程如图3所示。

Figure 3. Flow chart of association rule analysis based on Apriori

图3. 基于Apriori的关联规则分析流程图

基于Apriori的关联规则分析流程图如图3所示，将最小支持度和最小置信度分别设置为0.3和0.75，总共得到48条关联规则，单属性对的关联关系共12条，多属性对之间的关联关系共36条。观察得知多属性组之间的关联关系均依托于单属性对之间的关联关系，数量较多，故本文仅就部分单属性关联规则结果进行展示，参见表2。

Table 2. Association rules

表2. 关联规则

由表2可以发现，测点之间的关联规则能够形成一个循环，即：测点3 = “二级” à 测点2 = “二级” à 测点4 = “二级”测点1 = “二级” à 测点3 = “二级”。因此，对于正常工作下的螺杆泵，在不同的中心频率下，不同测点的相同振动加速度级总是同时出现的，即不同测点间的振动加速度并无显著性差异。

3.3. 基于相关系数的不同测点间的相关性分析

为验证不同测点间的相关性，本文考虑基于相关系数进行相关性分析。相关系数是用以反映变量之间相关关系密切程度的统计指标 [8]。由于研究对象的不同，相关系数有多种定义方式，当研究数据呈正态分布时较多采用皮尔逊相关系数，不满足正态分布时较多使用斯皮尔曼相关系数。本研究中的数据不满足正态分布，故使用斯皮尔曼相关系数。

斯皮尔曼相关系数 [9] 定义如下：对于样本容量为n的样本，n个原始数据被转换成等级数据，斯皮尔曼相关系数ρ为

$\rho =\frac{{\displaystyle {\sum }_i\left(x_i-\bar{x}\right)\left(y_i-\bar{y}\right)}}{\sqrt{{{\displaystyle {\sum }_i\left(x_i-\bar{x}\right)}}^2{{\displaystyle {\sum }_i\left(y_i-\bar{y}\right)}}^2}}$

x_i，y_i代表降序后的等级数据， $\bar{x}$， $\bar{y}$ 分别代表x和y的均值。

计算测点间的相关系数，结果如表3所示。可以发现，测点间的相关系数均在0.92以上，呈强正相关，印证了由关联分析得到的不同测点间并无显著性差异的结论。

Table 3. Correlation coefficient between each measurement point

表3. 各测点间的相关系数

4. 多阶拟合法

由上文的分析可知，测点间并无显著性差异，因此单测点分析和多测点分析并不对诊断结果造成影响。在此基础上，应用多阶拟合法建立中心频率和振动加速度的函数关系式，以系数向量间的欧式距离作为诊断依据，将其与适当的阈值作比较从而得出故障诊断结果。

4.1. 模型的建立

在螺杆泵正常运行的前提下，测点间无显著性差异，因此可用任一测点代表其余测点。本文以中心频率为自变量，任一测点作为因变量求解数学模型，函数关系式形式如下，其中x表示中心频率， ${\omega }_i$ 表示各阶的权重。

$f\left(x\right)={\omega }_0+{\omega }_{1}x+\cdots +{\omega }_n{x}^n={\displaystyle \underset{i=0}{\overset{n}{\sum }}{\omega }_i{x}^i}$

多项式拟合是曲线拟合方法之一，目的是使求得的拟合曲线能够反应数据的总体分布，从而使得数据能够分布在拟合曲线的附近。为了评价拟合函数的优劣，需要建立损失函数，测量每个样本点目标值与预测值之间的误差，拟合的目标是让误差最小。

假设有数据点 $\left(\left(x_1,y_1\right),\left(x_2,y_2\right),\cdots ,\left(x_m,y_m\right)\right)$，其中x_i表示中心频率，y_i表示振动加速度。本文应用最小二乘法 [10]，求解拟合函数解析式，均方误差函数为

$\delta \left({\omega }_0,\cdots ,{\omega }_n\right)={\displaystyle \underset{i=1}{\overset{m}{\sum }}{\left[f\left(x\right)-y_i\right]}^2}$

求解出最小 $\delta$ 值，即可解出f(x)的函数系数 ${\omega }_1,{\omega }_2,\cdots ,{\omega }_n$。为求得误差函数的最小值，对函数两边求导。

$\begin{array}{l}\frac{\partial \delta }{\partial {\omega }_0}=2{\displaystyle \underset{i=1}{\overset{m}{\sum }}\left[f\left(x\right)-y_i\right]},\\ \frac{\partial \delta }{\partial {\omega }_1}=2{\displaystyle \underset{i=1}{\overset{m}{\sum }}\left[f\left(x\right)-y_i\right]x},\\ ⋮\\ \frac{\partial \delta }{\partial {\omega }_n}=2{\displaystyle \underset{i=1}{\overset{m}{\sum }}\left[f\left(x\right)-y_i\right]x^n}\end{array}$

令上述等式为零，得到如下正规方程组：

$\{\begin{array}{l}m{\omega }_0+{\displaystyle \underset{i=1}{\overset{m}{\sum }}{\omega }_1{x}_i+\cdots +}{\displaystyle \underset{i=1}{\overset{m}{\sum }}{\omega }_n{x}_i^n}={\displaystyle \underset{i=1}{\overset{m}{\sum }}{y}_i}\\ {\displaystyle \underset{i=1}{\overset{m}{\sum }}{\omega }_0{x}_i}+{\displaystyle \underset{i=1}{\overset{m}{\sum }}{\omega }_1{x}_i^2+\cdots +}{\displaystyle \underset{i=1}{\overset{m}{\sum }}{\omega }_n{x}_i^{n+1}}={\displaystyle \underset{i=1}{\overset{m}{\sum }}{x}_i{y}_i}\\ ⋮\\ {\displaystyle \underset{i=1}{\overset{m}{\sum }}{\omega }_0{x}_i^n}+{\displaystyle \underset{i=1}{\overset{m}{\sum }}{\omega }_1{x}_i^{n+1}+\cdots +}{\displaystyle \underset{i=1}{\overset{m}{\sum }}{\omega }_n{x}_i^{2n}}={\displaystyle \underset{i=1}{\overset{m}{\sum }}{x}_i^n{y}_i}\end{array}$

将上式改写成矩阵形式：

$\left[\begin{array}{cccc}m& {\displaystyle \underset{i=1}{\overset{m}{\sum }}{x}_i}& \cdots & {\displaystyle \underset{i=1}{\overset{m}{\sum }}{x}_i^n}\\ {\displaystyle \underset{i=1}{\overset{m}{\sum }}{x}_i}& {\displaystyle \underset{i=1}{\overset{m}{\sum }}{x}_i^2}& \cdots & {\displaystyle \underset{i=1}{\overset{m}{\sum }}{x}_i^{n+1}}\\ ⋮& ⋮& \ddots & ⋮\\ {\displaystyle \underset{i=1}{\overset{m}{\sum }}{x}_i^n}& {\displaystyle \underset{i=1}{\overset{m}{\sum }}{x}_i^{n+1}}& \cdots & {\displaystyle \underset{i=1}{\overset{m}{\sum }}{x}_i^{2n}}\end{array}\right]\left[\begin{array}{c}{\omega }_0\\ {\omega }_1\\ ⋮\\ {\omega }_n\end{array}\right]=\left[\begin{array}{c}{\displaystyle \underset{i=1}{\overset{m}{\sum }}{y}_i}\\ {\displaystyle \underset{i=1}{\overset{m}{\sum }}{x}_i{y}_i}\\ ⋮\\ {\displaystyle \underset{i=1}{\overset{m}{\sum }}{x}_i^n{y}_i}\end{array}\right]$

可以证明，由于方程组的系数矩阵是一个对称正定矩阵，因此存在唯一解 ${\omega }_1,{\omega }_2,\cdots ,{\omega }_n$。

为探究在本实验中函数拟合的最佳阶数n，本文利用测点1的数据，拟合1至9阶的多阶函数进行对比分析。不同阶数所对应的均方误差及R²如表4所示，均方误差回归损失越小且R²越接近1，则说明拟合效果越好，从而确定最优的“最大阶数”。

Table 4. Mean square error and R2 of different orders

表4. 不同阶数的均方误差及R²

由表4可知，阶数为4所对应的RMSE最小，R²非负且接近于1，因此所对应的模型最佳。为了更直观的展示拟合效果，图4是在不同阶数下，拟合中心频率和测点1的关系曲线示意图，可以发现，四阶及以上拟合程度较高。当设置阶数越高，震荡越明显，出现过度拟合现象。因此，阶数为4的拟合函数的确能够得到最佳模型。在阶数为4时，可得到如下函数模型：

$y=-2.532\times {10}^{-13}\cdot x^4+4.304\times {10}^{-9}\cdot x^3-2.37\times {10}^{-5}\cdot x^2+0.04688\cdot x+77.41$

Figure 4. Fitting curves of center frequency and measurement Point 1 under different orders

图4. 中心频率和测点1在不同阶数下的拟合曲线

4.2. 阈值的确定

由4.1节可知，阶数为4的拟合函数最优。针对正常工作下的数据集，建立多个4阶一元回归模型，以中心频率为自变量，每个测点的振动加速度为因变量，相应的得到多个拟合系数向量。计算每个拟合系数向量的阶数相同的各个分量的标准差 ${\sigma }_i,i=1,2,3,4$。记X为测点拟合系数的数字矩阵，则

$X=\left[\begin{array}{cccc}{x}_{11}& x_{12}& x_{13}& x_{14}\\ x_{21}& x_{22}& x_{23}& x_{24}\\ ⋮& ⋮& ⋮& ⋮\\ x_{n1}& x_{n2}& x_{n3}& x_{n4}\end{array}\right]$

其中n代表测点数，列数为4对应模型的最佳拟合阶数，也就是说为每个测点都建立了一个一元4次回归模型，每一行对应每个测点的拟合系数向量。各个测点所对应的拟合系数向量参见表5。

Table 5. Fitting coefficient vector of each measurement point

表5. 各个测点的拟合系数向量

根据拟合系数向量求解 ${\sigma }_i,i=1,2,3,4$ 的计算公式如下，其中 ${\bar{x}}_i$ 为第i列的平均值。

${\sigma }_i=\frac{{\displaystyle \underset{k=1}{\overset{n}{\sum }}{\left(x_{ki}-{\bar{x}}_i\right)}^2}}{n-1}$

最后根据 $3\sigma$ 原则，阈值大小由一个正常分布在99.73%概率的置信区间的 $3\sigma$ 正态分布值来决定，具体计算公式如下。

$\theta ={\Vert 3{\sigma }_1,3{\sigma }_2,3{\sigma }_3,3{\sigma }_4\Vert }_2$

根据公式计算得到阈值大小为0.0051。为利用上述阈值进行故障诊断，应用多项式拟合法得到任一正常工作下测点的拟合系数向量 $x$ 与待诊断设备的测点所对应的拟合系数向量 $y$，二者之间的欧式距离Distance与阈值作比较，大于阈值即为故障。

$Distance=\sqrt{{\left(x-y\right)}^{\text{T}}\left(x-y\right)}$

因此，判断故障就可以通过评估Distance的大小是否超过一个代表正常工作下测点所对应的拟合系数向量的可变化程度的阈值。

5. 基于多阶拟合法的故障诊断

为了验证多阶拟合法行之有效，本文通过仿真实验，模拟泵在故障时的运行状态，收集故障测点数据。将正常数据和故障数据构成测试数据集，并应用上述方法进行实验。由于不同测点近乎等价，因此可以选用正常工作下的任一测点作为模板数据，模板数据所对应的拟合系数向量为模板系数向量，分别计算模板系数向量与测试数据集中所有数据样本的系数向量的欧氏距离，部分结果如表6所示。

Table 6. Calculation results of template data and test set data

表6. 模板数据与测试集数据的计算结果

当且仅当系数距离不超过4.2节中得到的阈值0.0051时，才将测点判定为正常，反之为故障。由上表可以看出，多阶拟合法能够诊断出所有的故障测点，但是会出现将正常测点划分为故障测点的情况，这很有可能意味着该测点所检测的螺杆泵存在潜在的故障因素，如设备运行时间过长从而导致内部零件受损等。由于将故障样本判定出错的代价远高于将正常样本判定错误的代价，因此，本文所提出的多阶拟合法及阈值的确定方法能够较为准确的对螺杆泵进行异常检测，将判定结果以混淆矩阵的形式展现，参见表7。

Table 7. Confusion matrix of fault diagnosis by multi-order fitting method

表7. 应用多阶拟合法进行故障诊断的混淆矩阵

基于混淆矩阵计算的评价指标参见表8，基于多阶拟合法的故障诊断的准确度高达0.9，就F值(0.89)而言，查全率和精确度之间达到了良好的平衡。因此，本文所提出的多阶拟合法具有良好的应用前景，具有实际意义。

Table 8. Evaluation index of multi-order fitting

表8. 多阶拟合法的评价指标

6. 结束语

本文的数据源是基于螺杆泵在正常工作下得到的，为挖掘数据之间的关系，本文通过Apriori算法发现螺杆泵正常工作下不同测点振动加速度间并无显著性差异，甚至近乎等价，且呈现出强正相关。因此，尽管数据项多，但是单属性分析和多属性分析并没有显著区别。

本文根据上述的数据分析结果，将中心频率作为自变量，测点数据作为因变量，建立多阶拟合模型。由于正常测点间并无显著差异，因此可以任选一个测点作为模板，分别建立模板数据及待诊断样本的多阶拟合模型，计算两个系数向量之间的欧氏距离，并通过统计方法设计适当的阈值，当系数距离大于指定阈值时，说明该测点所检测的设备很有可能出现异常，反之设备正常运行。将测试数据集应用到上述方法，故障诊断正确率高达0.9，F值为0.89，因此，本文所提出的多阶拟合法具有实际意义。

本文创新性的提出多阶拟合法，为螺杆泵的异常检测提供理论依据。但是仍存在一定的缺陷，比如数据量有限，在将来若能收集更多的数据就能够帮助我们得到统计上更显著的结论，从而不断拓展实验思路。此外，本文提出的多阶拟合法具有很强的可拓展性，其研究思想和实验思路除了可以应用于螺杆泵的异常检测之外，还能够应用于其他设备的故障诊断，应用前景十分广阔。

基金项目

天津市智能制造专项资金项目201810602，201907206，201907210；天津市互联网先进制造专项资金项目18ZXRHGX00110。