杨飞+乔铁柱+庞宇松+阎高伟

摘 要： 针对球磨机系统多模态复杂过程中的料位不确定性，球磨机振动信号存在非线性、噪声和外界干扰等问题，采用一种基于改进的高斯混合回归（GMR）的球磨机料位软测量方法，解决传统高斯混合模型初始化含有噪声和异常值的数据难以聚类的问题。首先，利用改进的K?medoids聚类算法与EM算法分别初始化和优化高斯混合模型（GMM）的最佳高斯分量个数、最优模型参数，然后采用GMR预测输出球磨机料位。最后实验验证了改进GMR模型得到的预测料位可以很好地跟踪真实料位，并且通过实验结果的对比分析，验证了改进模型的有效性和实用性以及较好的预测精度。

关键词： 球磨机料位； 多模态； 振动信号； GMM； 聚类； 软测量； GMR

中图分类号： TN98?34； TP181 文献标识码： A 文章编号： 1004?373X（2018）05?0153?06

Abstract： Since the fill level of the ball mill system in multimode complicated process is uncertain， and the vibration signal of ball mill has the characteristics of nonlinearity， noise and outside interference， a soft measurement method for ball mill fill level based on improved Gaussian mixture regression （GMR） is proposed to solve the problem that it is difficult to cluster the data embedding noise and abnormal value of the traditional Gaussian mixture model （GMM） initialization. The improved K?medoids clustering algorithm and EM algorithm are used respectively to initialize and optimize the optimal Gaussian component quantity and optimal model parameters. The GMR is used to predict the output level of the ball mill. The experimental results verify that the predicted fill level obtained by improved GMR model can track the real fill level accurately. The comparative analysis of experimental results verifies that the improved model is feasible and practical， and has high prediction accuracy.

Keywords： ball mill fill level； multimode； vibration signal； Gaussian mixture model； clustering； soft measurement； Gaussian mixture regression

0 引 言

球磨機是磨矿工业中广泛使用的基础设备，准确测量球磨机料位是实现优化控制和节能降耗的关键技术之一。近年来，科研人员提出了许多测量球磨机料位的方法，其一般步骤是：先对采集的信号计算频谱特征，然后进行降维和特征提取，最后建立频谱特征与料位之间的模型。文献[1]利用快速傅里叶变换方法求取信号的功率谱密度（Power Spectrum Density，PSD），然后采用主元分析（Principal Component Analysis，PCA）方法提取特征和减少冗余信息，最后用支持向量回归（Support Vector Regression，SVR）建立回归模型。文献[2]利用PCA方法进行振动信号的特征提取，然后采用偏最小二乘回归方法（Partial Least Square Regression，PLSR）建立球磨机料位软测量模型。文献[3]采用偏最小二乘法（Partial Least Square，PLS）对球磨机振动信号进行特征提取和降维，然后采用极限学习机（Extreme Learning Machine，ELM）建立料位软测量模型。上述方法均取得了较好的结果。

文献[4]通过分析球磨机研磨过程的机理和筒体振动的加速度频域信号，将球磨机运行过程分解成若干个不同的特征模态，然后采用KPLS算法建立预测子模型，最后利用信息熵值的集成加权融合方法预测输出。

在此，对球磨机在不同料位条件下多工况的复杂多模态过程，由于不同料位采集数据的均值和协方差不断变化，基本的单高斯分布建立的模型较差，无法准确预测料位。高斯混合模型（GMM）是单高斯概率密度函数的延伸，具有平滑逼近任意形状密度分布的特性[5]，因此可以对复杂的球磨机数据的多模态性进行很好的描述。针对球磨机系统多模态复杂过程中料位的不确定性，球磨机振动信号存在非线性、噪声和外界干扰等问题，以及GMM在多模态过程监督和软测量方面的成功应用[6?7]，采用一种基于改进的高斯混合回归（GMR）在球磨机多模态过程料位软测量的方法。首先，利用改进的K?medoids算法[8]和EM[9]算法优化高斯混合模型的最佳分量个数和最优模型参数，然后利用GMR预测料位的输出。实验结果表明该方法具有较高的测量精度。endprint

1 高斯混合模型及EM算法

1.1 高斯混合模型

高斯混合模型（Gaussian Mixture Model，GMM）是描述混合密度分布的模型，用有限多个单高斯概率密度函数的加权形式逼近任意的连续分布。高斯混合模型的概率密度函数[10]定义为：

式中：是GMM中单高斯分量的个数；是第个单高斯分量的概率密度函数；是第个单高斯分量的权重，且满足：是第个单高斯分量的参数集，和分别为第个单高斯的均值向量和协方差矩阵；表示个单高斯模型参数组成的全局参数集。因此，含有个高斯分量的高斯混合模型的所有参数表示为：。

给定个独立同分布的训练数据样本的似然函数和对数似然函数分别为：

1.2 改进的GMM初始化聚类K?medoids方法

传统GMM利用K?means聚类算法对模型参数进行初始化，然后通过期望最大化（Expectation?Maximization，EM）算法估计参数，由于EM算法对初始值要求较高，初值不恰当可能导致算法收敛到局部最大值点，在此对传统的GMM进行改进，利用改进的K?medoids算法优化GMM模型的初始化参数。

与K?means中心点选取不同，K?medoids中心点选取的规则是从当前cluster中選取的中心点到其他所有cluster点的距离之和最小。K?means聚类算法产生类的大小相差不会很大，对于含有噪音或者异常值的数据很敏感，而K?medoids不容易受到那些由于误差之类的原因产生的噪音数据的影响。由于球磨机运行过程中，料位测量容易受到各种因素的干扰而出现噪音和异常值，K?medoids[11]在处理这样的数据方面鲁棒性较强，具有一定的优势，可以提高GMM模型的性能。

K?medoids初始化GMM参数过程如下：

1） 从训练数据样本中随机选择个点作为初始的聚类中心点

2） 将训练数据样本中其他的数据点根据最近原则划分到个聚类中；

3） 根据下列公式更新每个聚类的中心点；

4） 返回步骤2），当各个聚类中心点不再变化时，执行步骤5）；

5） 步骤4）得到样本的个聚类中心点和其他样本所属类的标签属性，根据每个聚类中的样本个数可以计算出每个高斯分量的权重比值即将中心点作为GMM模型的均值，即根据每个聚类中的样本，利用协方差公式可以计算出每个高斯分量的协方差

通过以上步骤可得到GMM模型的初始化参数：

1.3 EM算法

EM算法是一种从“不完全数据”中求解模型分布参数的极大似然估计方法。EM算法通过不断重复E步骤（E?step）和M步骤（M?step），直到对数似然函数收敛到一定的阈值，最终获得GMM模型中的未知参数和权值

EM算法求解未知参数和权值的迭代步骤[12]如下：

2 基于高斯混合回归的软测量

2.1 高斯混合回归（GMR）

高斯混合回归（Gaussian Mixture Regression，GMR）[12]主要是基于高斯条件和高斯分布的线性组合特性。假设数据向量是由两部分组成：输入和输出如果服从含有个分量的高斯混合模型分布，同时，每个高斯分量的均值向量和协方差矩阵可以划分为以下输入和输出部分：

对于第个高斯分量，给定输入时，相应的输出也服从高斯分布，关于的条件概率定义为：

其中，均值和协方差的计算公式为：

由于输入是由混合模型生成的，输出的分布式由部分组成的，就整个混合模型而言，关于输入的输出的期望条件分布也是一个高斯混合模型：

式中：是关于输入的第个高斯分量的后验概率，由贝叶斯公式可得：

最后，给定一个输入输出的条件期望可以由高斯分布估计出，基于高斯分布的线性转换特性，输出的均值向量和协方差矩阵的计算公式为：

2.2 拟合混合模型

在预处理混合模型时，无论是GMM模型还是GMR模型，都需要给定高斯分量个数，此外，EM算法的参数估计过程也需要预先定义混合分量的个数在本文中，采用最小梯度准则（Minimum Gradient Citerion，MGC）选择最优的高斯分量个数，MGC定义如下：

式中：分别代表的最小值和最大值；代表含有个高斯分量的GMM的概率密度的对数似然函数；grad是对所有的对数似然函数值求梯度值；是最佳高斯分量个数。在高斯混合模型中，模型性能一般都随着值的增加而提高，但是过大的值会增加模型复杂度，并且模型的泛化能力会变差，当模型的对数似然函数不再增大或者增大较缓慢时，取此时的模型为最优模型，所以在此采用最小梯度准则来选择恰当的值。

在上述过程中，当变化时需要反复执行EM算法来估计参数，EM算法的初始化参数由K?medoids聚类算法得到。拟合混合模型的算法步骤如下：

预处理：预处理训练集，给定和选取准则函数MGC

初始化：令

使用改进K?medoids初始化个分量的高斯混合模型的参数，记为

主循环：当重复：

采用EM算法估计模型的新参数并且计算相应的对数似然函数；

删除掉最不可能的高斯分量，并且与它最接近的一个高斯分量合并为一个新的高斯分量。这样得到分量个数为的高斯混合模型，并且具有参数；

令

最佳高斯分量：通过最小梯度准则函数选取

最终混合模型的参数估计为

上述步骤中关键的一步是删除掉最不可能的高斯分量，并且与它最接近的一个高斯分量合并为一个新的高斯分量。首先，最不可能的高斯分量可以通过寻找权重最小得到，即：

然后，与第个高斯分量最接近的高斯分量可以得到，对称Kullvack?Leibler（KL）[13]散度是度量概率相似性的常见准则。对于高斯密度，利用对称KL散度准则可以求出GMM中其他的高斯分量与第个高斯分量之间的相似度，如下：endprint

与第个分量最相似的另一个高斯分量是

这样，第个和第个高斯分量将会合并为一个新的分量，根据文献[14]，新高斯分量的权重均值和协方差通过下式计算：

2.3 基于GMR的软测量

假定过程输入变量是输出变量是将输入变量和输出变量合并为一个新的向量那么和的联合概率密度就是的概率密度，用GMM分布表达为：

如果新的输入变量是相应的输出变量是为了预测输出，需要构建基于GMR的软测量模型。首先，均值和协方差可以划分为输入、输出部分，如式（7）。然后，对第个高斯分量，关于的后验概率和条件概率估计为：

式中：和是第个条件高斯分布的均值和协方差参数，可以通过式（9）计算。

最后，预测输出结果是各个高斯分量的加权和：

基于GMR的软测量过程主要分为两个阶段：离线建模阶段和在线测量阶段，其流程图如图1所示。

3 基于改进的GMR球磨机料位软测量

3.1 实验数据处理

为验证本文所采用方法的有效性，针对实验室小型球磨机进行料位软测量实验研究。料位样本分别为1 L，2 L，…，20 L，实验球磨机筒体振动信号由现场同步采集，对每个样本下的料位进行振动信号采集，最终获得20组振動信号数据。然后，采用Welch方法[15]计算振动信号的功率谱密度（Power Spectrum Density，PSD），图2为振动信号的功率谱图。将每组信号的PSD平均分成22个样本，则样本总数为440，并把每组的22个样本随机分成15个训练样本和7个测试样本，因此最终得到训练集有300个样本，测试集有140个样本。

由图2可以看出，球磨机研磨过程振动信号的能量主要集中在600～6 000 Hz频率范围内，即振动信号的有效频谱。综合考虑频谱的波动范围及计算效率，以100 Hz为单位频段对有效频谱进行分割并求均值。因此，最终得到的训练集数据维数大小为300×54，测试集数据维数大小为140×54。

3.2 测量结果

为了评价所提供方法的建模能力，采用均方根误差（Root Mean Square Error，RMSE）和平均绝对误差（Mean Absolute Error，MAE）作为测量精度的评价指标，其中将RMSE作为主要的评价指标。其计算公式为：

式中：和分别表示第个样本的实际值和估计值；为测试样本个数。

根据拟合高斯混合模型，采用训练集数据进行离线建模，设置初始化参数最后通过拟合混合模型方法优化得到的最佳高斯分量个数是然后采用GMR方法对球磨机料位进行在线软测量，测试集数据预测结果如图3所示，图中实线代表球磨机料位的真实值，星号代表料位的预测值，从图中可以看出球磨机料位预测值可以很好地跟踪期望的真实值。本实验软测量结果的测量精度均方根误差RMSE为0.358 8，平均绝对误差MAE为0.261 2。

3.3 结果分析与对比

为验证所采用改进的GMR模型的软测量建模能力和有效性，将其与传统GMR方法、主元回归（PCR）、偏最小二乘回归（PLSR）和支持向量回归（SVR）方法进行比较。采用RMSE作为主要评价指标。表1给出了上述五种方法的料位预测评价指标对比结果。图4为所提方法与对照方法测试的结果比较。

根据表1中的评价指标RMSE和图4可以看出，改进GMR模型得到的测试集预测结果最好，且要明显好于其他四种方法，改进GMR模型得到的预测料位可以很好地跟踪真实料位，达到较好的预测精度。分析其原因，主要是改进GMR在处理含有噪声和外界干扰的球磨机数据时，K?medoids初始化GMM得到较优的参数，由于球磨机系统运行时存在料位的不确定性，改进GMR可以建立表征不同工况的球磨机状态模型，并且可以对球磨机数据的多模态性进行很好的描述。除此之外，当球磨机运行过程中工况改变或者受到其他干扰时，球磨机料位测量会出现异常值，改进GMR模型受到的影响较小，并可以很好地跟踪其真实料位，实验结果说明改进GMR预测模型的准确性要优于传统GMR，PCR，PLSR和SVR方法。

4 结 语

针对球磨机系统复杂过程中的料位不确定性，球磨机振动信号存在非线性特性等问题，本文采用一种基于改进的高斯混合回归的球磨机料位软测量方法。首先，利用改进的K?medoids算法和EM算法优化GMM的最佳高斯分量个数和最优模型参数，然后采用GMR预测输出球磨机料位。实验结果表明，改进GMR模型在处理球磨机系统多模态复杂过程中具有较好的优势，显示了GMR在多模态过程中软测量的有效性和可行性。下一步将研究球磨机在实际的磨矿过程中的模态分析方法，对球磨机特有的多模态过程具体划分模态工况，从而针对不同模态具体分析球磨机料位，进一步提高球磨机料位的预测精度。

参考文献

[1] TANG Jian， ZHAO Lijie， YU Wen， et al. Soft sensor modeling of ball mill load via principal component analysis and support vector machines [J]. Lecture notes in electrical engineering， 2009（67）： 803?810.

[2] 汤健，郑秀萍，赵立杰，等.基于频域特征提取与信息融合的磨机负荷软测量[J].仪器仪表学报，2010，31（10）：2161?2167.

TANG Jian， ZHENG Xiuping， ZHAO Lijie， et al. Soft sensing of ball mill based on frequency domain feature extraction and information fusion [J]. Chinese journal of scientific instrument， 2010， 31（10）： 2161?2167.endprint

[3] TANG Jian， WANG Dianhui， CHAI Tianyou. Predicting mill load using partial least squares and extreme learning machines [J]. Soft computing， 2012， 16（9）： 1585?1594.

[4] 汤健，柴天佑，赵立杰，等.基于振动频谱的磨矿过程球磨机负荷参数集成建模方法[J].控制理论与应用，2012（2）：184?186.

TANG Jian， CHAI Tianyou， ZHAO Lijie， et al. Ensemble modeling for parameters of ball?mill load in griding process based on frequency spectrum of shell vibration [J]. Control theory and applications， 2012（2）： 184?186.

[5] BISHOP C M. Pattern recognition and machine learning [M]. New York： Springer， 2007： 432?443.

[6] YU J. Multiway Gaussian mixture model based adaptive kernel partial least squares regression method for soft sensor estimation and reliable quality prediction of nonlinear multiphase batch processes [J]. Industrial and engineering chemistry research， 2012， 51（40）： 13227?13237.

[7] GRBI? R， SLI?KOVI? D， KADLEC P. Adaptive soft sensor for online prediction and process monitoring based on a mixture of Gaussian process models [J]. Computers and chemical engineering， 2013， 58（22）： 84?97.

[8] YUAN Xiaofeng， GE Zhiqiang， SONG Zhihuan. Soft sensor model development in multiphase/multimode processes based on Gaussian mixture regression [J]. Chemometrics and intelligent laboratory systems， 2014（138）： 99?103.

[9] GHAHRAMANI Z， JORDAN M I. Supervised learning from incomplete data via an EM approach [J]. Advances in neural information processing systems， 1993（6）： 120?127.

[10] KAUFMAN L， ROUSSEEUW P J. Finding groups in data： anintroduction to cluster analysis [M]. New York： John Wiley & Sons， 1990： 108?123.

[11] MCLACHLAN G， KRISHNAN T. The EM algorithm and extensions [M]. New York： John Wiley & Sons， 2007： 105?132.

[12] KULLBACK S. Information theory and statistics [M]. US： Courier Dover Publications， 2013： 189?207.

[13] 周越，司刚全，曹晖，等.功率谱分析在筒式钢球磨煤机内存煤量测量中的应用研究[J].工业仪表与自动化装置，2006，35（6）：21?24.

ZHOU Yue， SI Gangquan， CAO Hui， et al. A study of the power spectrum analysis application to ball mill load measurement [J]. Industrial instrumentation & automation， 2006， 35（6）： 21?24.

[14] YU J， QIN S J. Multimode process monitoring with Bayesian inference?based finite Gaussian mixture models [J]. Aiche journal， 2008， 54（7）： 1811?1829.

[15] CHEN Xinquan， PENG Hong， HU Jingsong. K?medoids su?batitution clustering method and a new clustering validity index method [C]// Proceedings of the 6th World Congress on Intelligent Control and Automation. [S.l.]： IEEE， 2006： 5896?5900.endprint