CN101776531B - 一种球磨机负荷参数软测量方法 - Google Patents

Abstract

一种球磨机负荷参数软测量方法，该方法通过硬件支撑平台获得的磨机筒体振动信号、振声信号及电流信号对表征磨机负荷的磨机内部参数(料球比、矿浆浓度、充填率)进行软测量。该方法包括如下步骤：采集磨机筒体振动、振声和电流数据及时域滤波，对振动及振声数据进行时频转换，对振动及振声数据进行频域内的分频段的基于核主元分析的非线性特征提取，对时域电流数据进行非线性特征提取，对融合后非线性特征数据进行特征选择，建立基于最小二乘-支持向量机的软测量模型。本发明的软测量方法灵敏度高，测量结果准确，具有很好实用价值和推广前景，有助于实现磨矿生产过程的稳定控制、优化控制和节能降耗。

Description

一种球磨机负荷参数软测量方法

技术领域：

本发明属于自动化技术测量领域，特别是涉及磨矿生产过程的一种湿式球磨机负荷参数即球磨机内部料球比、矿浆浓度、充填率的软测量方法。 

背景技术：

筒式球磨机在金属矿及非金属选矿厂、建材、耐火材料、水泥、煤炭、化工、电力和冶金等行业的粉磨作业中应用广泛。粉磨作业的电耗占全世界总发电量的3～4％，在选矿、电力、水泥等行业分别占各自工业过程能耗的30％～70％、15％、60％～70％，磨矿介质和衬板的消耗达0.4-3.0kg/t。根据试验研究，球磨机至少有10％以上的节能潜力和9％以上的节约钢材原材料的潜力。由于球磨机具有非线性、大时滞、强耦合、随机干扰大等综合复杂特性，旋转筒体内部运行环境复杂，对磨机运行状态及其内部参数难以进行有效测量，生产中经常发生空磨、饱磨、堵磨等现象，直接影响到产品的质量、产量以及技术经济指标。 

目前，实际生产中一般是由操作员通过振声(俗称电耳)、振动或功率等间接方法，凭经验对球磨机负荷状态进行人工判断与处理。为克服人工操作的主观性和随意性，结合领域专家知识、规则推理和统计过程控制，出现了融合多种信号的智能球磨机负荷状态检测方法。但针对能够反映球磨机负荷的球磨机负荷参数(料球比、矿浆浓度、充填率)的研究较少，所建模型的输入多为轴承振动、振声、轴承压力、功率等信号的时域或频域能量之和。国外针对球磨机的轴承振动、振声的研究，采用频谱平均分段求和的方式提取频谱变量的特征，但其检测的球磨机内部参数仅限于矿浆浓度。 

发明内容：

针对现有的球磨机负荷的确定主要是依靠现场操作人员的长期工作经验，难以测量表征球磨机负荷的内部参数，而且现有检测方法难以充分提取信号的非线性，导致工业过程中为保证设备的安全性，球磨机常运行在欠负荷状态，影响球磨机处理量和产品质量，无法根据市场需求进行有针对性地生产而获得最大的经济效益的问题，本发明提供一种可测量现场球磨机负荷参数的基于多源数据特征融合的球磨机负荷软测量方法。通过对球磨机筒体振动及振声信号在频域内的低、中、高三个频段采用核主元分析(KPCA)分别进行非线性特征提取，并融合时域电流信号的非线性特征，结合最小二乘-支持向量机(LSSVM)方法进行特征选择，以达到准确检测球磨机负荷参数的目的。 

本发明的软测量方法由球磨机硬件支撑平台和软测量软件组成，硬件平台提供软测量方法所需数据，软测量软件负责实施本发明所提出的软测量方法。软测量软件即可以安装在分 布式计算机控制***的监控计算机上，也可以运行于独立的计算机上。该方法通过与硬件支撑平台进行通讯，获得实时的球磨机筒体振动、振声、电流信号的数据，检测球磨机负荷参数(料球比、矿浆浓度、充填率)。 

本发明包括步骤如下：如图1和2所示 

步骤一：采集球磨机振动、振声及电流数据 

通过硬件支撑平台采集如下信号：球磨机筒体的振动加速度信号X_V；球磨机研磨区域下方的振声信号X_A；球磨机驱动电动机的电流信号X_I； 

对以上信号进行信号处理，将提取和选择的非线性特征做为软测量模型的输入，在线估计球磨机负荷内部参数即料球比、矿浆浓度、充填率，球磨机负荷参数软测量方法按照如下的结构组成样本，并收集样本数据。样本表达式为{x_k，y_k}，其中，x_k为样本的输入即采集的信号——球磨机筒体的振动信号X_V、振声信号X_A、电流信号X_I，样本的输出y_k为待估的主导变量——球磨机负荷参数：料球比y₁、矿浆浓度y₂、充填率y₃。样本的记录机构如表1，时间为样本获取的时间。 

表1 样本数据结构 

考虑到样本数据应该具有代表性，并且尽可能的覆盖范围较款，应该包括球磨机研磨过程中的高矿浆浓度、低矿浆浓度；高充填率、低充填率；高料球比、低料球比等工况。将样本分为两组，一组用于用于模型误差最小的训练来选择模型参数(包括模型训练样本和误差训练样本)，另一组用于模型的验证。 

步骤二：对采集的球磨机振动、振声及电流数据进行滤波 

采用有限冲激响应滤波器分别对振动及振声信号进行带通滤波和低通滤波。滤波器公式见下式： 

$y\left(m\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}h\left(n\right)x(m-n)---\left(1\right)$

式中，h(n)——滤波器的系数； 

N——滤波器的长度； 

y(m)——滤波器在第m点的输出； 

x(m-n)——滤波器的第(m-n)个输入； 

n——采样时刻n＝0，1，...N-1； 

m——当前采样时刻； 

采用均值滤波法对电流信号进行滤波，公式见下式： 

$y\left(m\right)=\frac{1}{N}\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_n---\left(2\right)$

式中，y(m)——滤波器的输出，即N次采样的算术平均值； 

x_n——第n次的采样值； 

N——采样点数；。 

步骤三：对滤波后的振动及振声数据进行时频转换 

球磨机机械研磨过程中产生的平稳周期性振动和振声信号中包含着大量与研磨过程相关的有用信息。但是时域内，该信息却被随机噪声信号“白噪声”所覆盖。球磨机筒体振动和振声信号可在频域内分解为叠加的正弦信号，而这些不同频段的正弦信号的幅值包含着直接与球磨机的操作状态相关信息。为保证对功率谱的无偏估计，采用改进的平均周期图法求球磨机旋转整周期的球磨机筒体振动及振声信号的功率谱密度(PSD)，按如下公式： 

$\tilde{P}\left(k\right)=\frac{1}{R}\underset{r=1}{\overset{R}{\mathrm{\Σ}}}{P}^r\left(k\right),r=1,\·\·\·,R---\left(3\right)$

$P^r\left(k\right)=\frac{1}{N}{\left|X_N\left(k\right)\right|}^2---\left(4\right)$

式中， 

——总功率谱； 

P^r(k)——第r段的功率谱； 

R——分段的段数； 

N——X数字序列的数据点数； 

k——为求得的频率点数。 

步骤四：对振动及振声信号的频谱数据和时域电流数据进行非线性特征提取 

4.1对振动及振声信号的频谱数据数据进行非线性特征提取 

分析球磨机不同工况下球磨机频域信号的特征，可知空转时频域信号主要集中在低频段、干磨时信号扩展至中频段；而湿磨时其中频段特征明显，而且在球磨机内的浓度较低，而充 填率较高时，存在明显的高频段；水磨时则高频段特征明显。因此，可将湿式球磨机频域信号划分为低频段，中频段及高频段，其分别对应着球磨机的固有频率段、主要冲击频率段、次要冲击频率段。将球磨机筒体振动信号在低、中、高频的数据分别记做X_{V_LF_org}、X_{V_MF_org}、X_{V_HF_org}；振声信号在低、中、高频的非线性特征分别记做X_{A_LF_org}、X_{A_MF_org}、X_{A_HF_org}，这些频谱数据存在高维和共线性的问题，采用核主元分析(KPCA)分别对X_{V_LF_org}、X_{V_MF_org}、X_{V_HF_org}、X_{A_LF_org}、X_{A_MF_org}、X_{A_HF_org}进行降维和非线性特征提取。 

基于核主元分析的非线性特征提取过程描述如下： 

核主成分分析的基本思想是先通过一个非线性变换Φ把输入数据X(X∈R^N)映射到一个高维特征空间F上，F＝{Φ(X)：X∈R^N}，然后在特征空间F上执行经典的线性主成分分析，以便实现输入空间无法实现的线性分类。假设映射后的样本矢量已经被中心化，即满足其均值为0，可表示为 $\underset{k=1}{\overset{M}{\mathrm{\Σ}}}\mathrm{\Φ}\left(X_k\right)=0,$ 则映射后在F空间上样本集的协方差矩阵可定义为 

$C=\frac{1}{M}\underset{k=1}{\overset{M}{\mathrm{\Σ}}}\mathrm{\Φ}\left(X_k\right)\mathrm{\Φ}{\left(X_k\right)}^T---\left(5\right)$

与该矩阵对应的特征值方程为 

λv＝Cv                                            (6) 

根据再生核(theory of reproducing kernel)理论，由于在特征空间F中，对于任一特征值λ≠0所对应的特征向量v一定位于由{Φ(X₁)，...，Φ(X_M)}张成的空间内。考虑到特征空间F中的所有向量都可由Φ(X_i)线性组合，所以有： 

$v=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_k\mathrm{\Φ}\left(X_k\right)---\left(7\right)$

M×M阶核矩阵K定义为： 

K_i，j＝K(X_i，X_j)＝Φ(X_i)^TΦ(X_j)                    (8) 

式中：i、j均表示具体的第几个样本； 

将式(5)、式(7)、式(8)代入式(6)后，可得以下全部用核函数表示的特征值方程 

$\mathrm{K\α}=\mathrm{M\λ\α}=\widehat{\mathrm{\λ}}\mathrm{\α}---\left(9\right)$

这样就将求解式(6)特征向量v的问题转换成求解式(9)特征值方程的特征向量α的问题。设该特征值方程的解由大到小为：λ₁≥λ₁≥...≥λ_M≥0，对应的特征向量为α₁，α₂，...，α_M，对应于式(6)的特征矢量为v₁，v₂，...，v_M，并假设λ_M为最后一个非零特征值，则从特征矢量的正交规一条件 

$v_g^T{v}_r={\mathrm{\δ}}_{g,r},\∀g,r=\mathrm{1,2},\·\·\·,M---\left(10\right)$

表示对所有的g，r＝1，2，...，M都有 $v_g^T{v}_r=1$

就可以推导出关于特征向量α_k的正交规一条件，即 

$\widehat{\mathrm{\λ}}{\mathrm{\α}}_g^T{\mathrm{\α}}_r={\mathrm{\δ}}_{g,r},\∀g,r=\mathrm{1,2},\·\·\·,M---\left(11\right)$

进而由式(9)与条件式(11)就可以得出特征向量的唯一解{α_k，k＝1，2，...M}； 

于是对输入空间的任一样本X，可按下式求得其在高维空间F中Φ(X)的第h个主成分t_h，即它在第h个主成分特征矢量v_h方向上的投影： 

$t_h=(v^h\·\mathrm{\Φ}\left(X\right))=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_k^{\left(h\right)}(\mathrm{\Φ}\left(X_k\right)\·\mathrm{\Φ}\left(X\right))=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_k^{\left(h\right)}K(X_k,X),k=1,\·\·\·M---\left(12\right)$

式中：t_h——第h个主元； 

α_k——协方差阵的特征向量； 

λ——协方差阵的特征值； 

Φ——非线性变换函数； 

M——样本数据的个数； 

v^h——高维空间中，第h个主成份的特征矢量； 

K——M×M阶核矩阵； 

h——主元的个数。 

之前假设高维特征空间的样本矢量已经中心化，故需要对由原始数据空间映射后的核矩阵按下式进行中心化处理： 

$\tilde{K}=K-I_{M}K-{\mathrm{KI}}_M+I_M{\mathrm{KI}}_M---\left(13\right)$

式中： 

——中心化处理后的核矩阵； 

I_M——系数为1/M的单位阵； 

K——中心化处理之前的核矩阵。 

从而，对振动及振声信号共六个频段提取非线性特征均按如下公式进行提取： 

$t_h=(v^h\·\mathrm{\Φ}\left(X\right))=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_k^{\left(h\right)}\tilde{K}(X_k,X)---\left(14\right)$

将球磨机筒体振动信号在低、中、高频的频谱变量X_{V_LF_org}、X_{V_MF_org}、X_{V_HF_org}提取的非线性特征分别记做X_{V_LF}、X_{V_MF}、X_{V_HF}；球磨机振声信号在低、中、高频的频谱变量X_{A_LF_org}、X_{A_MF_org}、X_{A_HF_org}提取的非线性特征分别记做X_{A_LF}、X_{A_MF}、X_{A_HF}，其值分别用下式表示： 

$X_{V\_\mathrm{LF}}=[t_{V\_\mathrm{LF}\_1},t_{V\_\mathrm{LF}\_2},\·\·\·,t_{V\_\mathrm{LF}\_h_{{\mathrm{VLF}}_{\max }}}];$

$X_{V\_\mathrm{MF}}=[t_{V\_\mathrm{LF}\_1},t_{V\_\mathrm{MF}\_2},\·\·\·,t_{V\_\mathrm{MF}\_h_{{\mathrm{VMF}}_{\max }}}];$

$X_{V\_\mathrm{HF}}=[t_{V\_\mathrm{HF}\_1},t_{V\_\mathrm{HF}\_2},\·\·\·,t_{V\_\mathrm{HF}\_h_{{\mathrm{VHF}}_{\max }}}];$

$X_{A\_\mathrm{LF}}=[t_{A\_\mathrm{LF}\_1},t_{A\_\mathrm{LF}\_2},\·\·\·,t_{A\_\mathrm{LF}\_h_{{\mathrm{ALF}}_{\max }}}];---\left(15\right)$

$X_{A\_\mathrm{MF}}=[t_{A\_\mathrm{MF}\_1},t_{A\_\mathrm{MF}\_2},\·\·\·,t_{A\_\mathrm{MF}\_h_{{\mathrm{AMF}}_{\max }}}];$

$X_{A\_\mathrm{HF}}=[t_{A\_\mathrm{HF}\_1},t_{A\_\mathrm{HF}\_2},\·\·\·,t_{A\_\mathrm{HF}\_h_{{\mathrm{AHF}}_{\max }}}];$

式中， 

——表示各个频段提取的最大主元个数。其值的采用如下公式： 

$h_{{\mathrm{MNF}}_{\max }}={\arg }_{h_M}\min \{{\mathrm{CPV}}_{\mathrm{hM}}\≥99\%\};$

式中，CPV_hM——表示振动或振声信号的前 

个主元的方差累积之和； 

M＝[V，A]——表示振动(V)和振声信号(A)； 

N＝[L，M，H]——本别表示频域信号的低(L)、中(M)、高(H)频段； 

对新样本数据，按下式生成核矩阵： 

$K_{i^{\′},j^{\′}}^{\mathrm{test}}=K(X_{i^{\′}}^{\mathrm{test}},X_{j^{\′}}),i^{\′},j^{\′}=1,\·\·\·,L$

式中，X^test——表示新样本数据； 

K^test——表示新样本数据的核矩阵； 

L——是新样本数据的样本数量即行数。 

I_N′是一个系数为1/L L×M的矩阵，L是新数据的行数。 

新样本数据核矩阵按下式进行中心化： 

${\tilde{K}}_{i,j}^{\mathrm{test}}=K^{\mathrm{test}}-I_M^{\′}K-K^{\mathrm{test}}{I}_M+I_M^{\′}{\mathrm{KI}}_M---\left(17\right)$

式中， 

——中心化后的新样本核矩阵； 

I_M′——系数为1/L L×M的矩阵； 

4.2对电流数据进行非线性特征提取 

诸多理论研究及实践应用表明，电流信号与球磨机负荷间存在明显的非线性关系。采用与步骤4.1中相同的方法即KPCA提取时域球磨机电流信号中的非线性特征。提取后的非线性特征记做X_I，其值可用下式表示： 

$X_I=[t_{I\_1},t_{I\_2},\·\·\·,t_{I\_h_{I\max }}]---\left(18\right)$

式中，t——主元，即提取的非线性特征；h_Imax——表示针对时域电流信号进行非线性特征提取后的最大主元个数。其值采用如下公式确定： 

$h_{I\max }={\arg }_{h_I}\min \{{\mathrm{CPV}}_{h_I}\≥99\%\}.---\left(19\right)$

式中， ——表示电流信号前h_Imax个主元的方差累积之和。 

步骤五：对振动、振声、电流非线性特征提取后的融合数据进行特征选择 

主元分析(PCA)基于线性相关和高斯统计假设，当应用于非线性***时，较小的主元虽然代表不重要的方差，但可能包含重要的***信息。虽然核主元分析(KPCA)方法是在高维特征空间中确定主元，具有很好的非线性逼近能力和特征抽取能力，但其本质仍然是描述响应变量X数据空间内的特征信息，而X中的特征信息的大小并不一定与预测变量Y相关。因此，需要合理确定KPCA的主元个数。同时，本文所要预测的磨机负荷参数，料球比、矿浆浓度及充填率与磨机振动、振声信号中的不同频段相关，需要为不同的预测模型选择不同数量的主元个数。模型输入变量的合理选择还可以提高磨机负荷参数模型的预测精度，防止过拟和。本发明将非线性特征输入变量中每类的的最佳变量个数(KPCA的主元个数)的集合称为特征参数集，并表示为： ${\mathrm{FP}}_{\mathrm{sel}}=\left\{(h_{{\mathrm{VLF}}_{\mathrm{sel}}},h_{{\mathrm{VMF}}_{\mathrm{sel}}},h_{{\mathrm{VLF}}_{\mathrm{sel}}})(h_{{\mathrm{ALF}}_{\mathrm{sel}}},h_{{\mathrm{AMF}}_{\mathrm{sel}}},h_{{\mathrm{ALF}}_{\mathrm{sel}}})\left(h_{I_{\mathrm{sel}}}\right)\right\}.$ 结合最小二乘-支持向量机模型，料球比、矿浆浓度及充填率软测量模型的特征参数集的确定均采 用如下的特征参数优化模型： 

$s.t.\left\{\begin{array}{c}1\≤h_{{\mathrm{VLF}}_{\mathrm{sel}}}\≤h_{{\mathrm{VLF}}_{\max }}\\ 1\≤h_{{\mathrm{VMF}}_{\mathrm{sel}}}\≤h_{{\mathrm{VMF}}_{\max }}\\ 1\≤h_{{\mathrm{VHF}}_{\mathrm{sel}}}\≤h_{{\mathrm{VHF}}_{\max }}\\ 1\≤h_{{\mathrm{ALF}}_{\mathrm{sel}}}\≤h_{{\mathrm{ALF}}_{\max }}\\ 1\≤h_{{\mathrm{AMF}}_{\mathrm{sel}}}\≤h_{{\mathrm{AMF}}_{\max }}\\ 1\≤h_{{\mathrm{AHF}}_{\mathrm{sel}}}\≤h_{{\mathrm{AHF}}_{\max }}\\ 1\≤h_{I_{\mathrm{sel}}}\≤h_{I_{\max }}\end{array}\right.---\left(20\right)$

式中，RMSE_val——最小二乘-支持向量机模型的最小评价误差； 

$h_{{\mathrm{MNF}}_{\max }}={\arg }_{h_{\mathrm{MN}}}\min \{{\mathrm{CPV}}_{\mathrm{hMN}}\≥99\%\}$ ——振动及振声信号在频域内的低中高 

三个频段进行KPCA降维和特征提取后能够选择的最大主元个数； 

CPV_hMN——表示振动或振声信号的前 

个主元的方差累积之和； 

M＝[V，A]——表示振动(V)和振声信号(A)； 

N＝[L，M，H]——本别表示频域信号的低(L)、中(M)、高(H)频段； 

$h_{I\max }={\arg }_{h_I}\min \{{\mathrm{CPV}}_{h_I}\≥99\%\}$ ——是时域电流信号进行KPCA非线性特征提取后的最大主元个数，一般取1； 

n——输入变量的行数即样本数； 

y——磨机负荷参数的真值； 

——磨机负荷参数的估计值即模型的计算值。 

按以下步骤对料球比模型、浓度模型、充填率模型进行特征选择，其选择过程首先是采用全部特征变量作为最小二乘-支持向量机模型的输入得到模型的参数，其次进行特征参数的选择。其流程如附图3所示。具体如下： 

(5.1)开始：进行初始化操作，按下式组合特征变量并记为X_all，作为最小二乘-支持向量机模型的输入变量：X_all＝[X_{V_LF}，X_{V_MF}，X_{V_HF}，X_{A_LF}，X_{A_MF}，X_{A_HF}，X_I]，此时的特征参数集 ${\mathrm{FP}}_{\mathrm{sel}\_\mathrm{all}}=\left\{(h_{{\mathrm{VLF}}_{\max }},h_{{\mathrm{VMF}}_{\max }},h_{{\mathrm{VLF}}_{\max }})(h_{{\mathrm{ALF}}_{\max }},h_{{\mathrm{AMF}}_{\max }},h_{{\mathrm{ALF}}_{\max }})\left(h_{I_{\max }}\right)\right\};$

(5.2)对全部样本数据进行标准化处理：将样本数据标准化为0均值1方差的数据； 

(5.3)选择误差惩罚参数和核参数的调整范围：依据经验确定对模型训练使用的误差惩罚函数和核函数的区间范围，便于选择最佳的模型参数； 

(5.4)调整误差惩罚参数和核参数：从各自参数的范围设置的下限开始，每次每个参数循环一个步长，作为调整后的参数，用于建立相应的模型并对该组参数进行误差评价； 

(5.5)建立以选择全部特征X_all为模型输入变量的模型：建立本模型的主要目的是初步选择模型参数：误差惩罚参数和核参数； 

基于最小二乘-支持向量机的模型的建立过程描述如下： 

对于给定训练集为[x(k)，y(k)]，k∈[1，M]，其中x(k)∈R^d，y(k)∈R，d为辅助变量的个数，支持向量机的基本建模思想是首先通过非线性映射Φ(X)将输入向量从原空间R^d映射到一个高维特征空间F，在这个高维特征空间中采用结构风险最小化原则构造最优决策函数，并利用原空间的核函数取代高维特征空间的点积运算以避免复杂运算，从而将非线性函数估计问题转化为高维特征空间的线性函数问题，其构造的最优决策函数具有如下形式： 

f(x)＝W^TΦ(x)+b                                     (21) 

式中，W——权向量； 

b——偏置量； 

Φ(x)——低维到高维空间的非线性映射； 

求解的目的就是利用结构风险最小化原则，寻找参数W^T和b，使得对于样本外的输入x，有|y_k-(W^TΦ(x_k)-b|≤ε，寻找W^T和b等价于求解下面的优化问题； 

$\min J=\frac{1}{2}{W}^{T}W+c\·R_{\mathrm{emp}}$

其中，W控制模型的复杂度；J是需要优化的性能指标；c＞0是误差惩罚函数，表示函数的平滑度和允许误差大于ε的数值之间的折衷；R_emp为经验风险，即ε不敏感损失函数，其定义为： 

$R_{\mathrm{emp}}\left(\mathrm{\θ}\right)=\frac{1}{M}\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{|y_k-(W^T\mathrm{\Φ}\left(x_k\right)+b)|}_{\mathrm{\ϵ}},k=1,\·\·\·,M;---\left(22\right)$

式中，M——样本数量； 

通过定义不同的损失函数，就可以构造不同的支持向量机，最小二乘支持向量机就是选择了误差的二次项ξ_k ²代替标准支持向量机优化目标中损失函数为ξ_k即允许错分的松弛变量，从而，最小二乘支持向量机的优化问题为： 

$\underset{W\·b}{\min }{J}_p=\frac{1}{2}{W}^{T}W+\frac{1}{2}c\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{{\mathrm{\ξ}}_k}^2---\left(23\right)$

s.t：y_k＝W^TΦ(x_k)+b+ξ_k

式中J_p是新的需要优化的性能指标，用拉格朗日方法求解上述优化问题，定义拉格朗日函数如下： 

$L(W,b,\mathrm{\ξ},\mathrm{\α})=\frac{1}{2}{W}^{T}W+\frac{1}{2}c\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{{\mathrm{\ξ}}_k}^2-\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_k\left[(y_k-(W^T\mathrm{\Φ}\left(x_k\right)+b+{\mathrm{\ξ}}_k)\right]---\left(24\right)$

根据优化条件 

$\frac{\∂L}{\∂W}=0,$ $\frac{\∂L}{\∂b}=0,$ $\frac{\∂L}{\∂\mathrm{\ξ}}=0,$ $\frac{\∂L}{\∂\mathrm{\α}}=0,$

可得 $W-\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_k\mathrm{\Φ}\left(x_k\right)=0,$ $\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_k=0,$ α_k＝cξ_k，W^TΦ(x_k)+b+ξ_k-y_k＝0，      (25) 

定义核函数k(x，x_k)＝(Φ(x)·Φ(x_k))代替非线性映射，根据上式可得求解的优化问题转化为求解线性方程： 

$\left[\begin{array}{cccc}0& 1& \·\·\·& 1\\ 1& k(x_1,x_1)+\frac{1}{c}& \·\·\·& k(x_1,x_M)\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ 1& k(x_M,x_1)& \·\·\·& k(x_M,x_M)+\frac{1}{c}\end{array}\right]\·\left[\begin{array}{c}b\\ {\mathrm{\α}}_1\\ \·\\ \·\\ \·\\ {\mathrm{\α}}_M\end{array}\right]=\left[\begin{array}{c}1\\ y_1\\ \·\\ \·\\ \·\\ y_M\end{array}\right]---\left(26\right)$

通过上式确定系数b和α＝[α₁，...，α_M]，最后可得软测量模型为 

$f\left(x\right)=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_{k}k(x,x_k)+b,k=1,\·\·\·,M---\left(27\right)$

其中核函数是满足默瑟即Mercer条件的任意对称函数，本方法采用径向基核函数代替非线性映射，建立软测量模型，该核函数形式为： 

$k(x_k,x_l)=\exp (-\frac{{\left|\right|x_k-x_l\left|\right|}^2}{2{\mathrm{\δ}}^2}),k,l=1,\·\·\·,M$ 式中：δ——核参数； 

(5.6)读取误差评价训练样本数据：读入准备用于误差评价的一组样本数据； 

(5.7)记录误差评价结果和参数：取误差训练样本数据集，定义误差函数为： 

$e=\sqrt{\frac{1}{L}\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{e_k}^2}=\sqrt{\frac{1}{L}\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{(y_k-(\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\α}}_{k}k(x,x_k)+b))}^2},k=1,\·\·\·,L---\left(28\right)$

式中，L——误差评价训练样本的样本数量即行数； 

最终的误差评价函数为： 

E(c，δ)＝min(e)； 

式中，c——为误差惩罚参数； 

在给定的参数区间内利用软测量模型得到误差评价指标，并记录对应的参数； 

(5.8)参数调整是否已经达到范围上限：定义l为调整步长，c为误差惩罚参数，c_up为惩罚参数范围的上限，δ为核参数，δ_up为核参数范围的上限。若c+l≥c_up和δ+l≥δ_up同时满足，则说明所有参数组合的误差评价工作完成，否则，重复步骤(5.4)～步骤(5.8)的工作； 

(5.9)记录误差评价最好的模型参数：将步骤(5.7)中记录的误差评价指标，寻找其中的极小值，选择对应的参数为模型参数，该模型参数将用做以下步骤中的LSSVM模型参数； 

(5.10)重新读入标准化后的训练样本：重新读入步骤(5.2)特征变量标准化的训练样本； 

(5.11)选择的特征个数是否达到设定值：对输入矩阵X_all的维数进行判断，如何维数小于设定值，则认为特征选择的计算量较小，采用网格算法进行特征选择，转到步骤(5.24)；否则采用遗传算法，转到步骤(5.12)，加快搜索速度； 

(5.12)特征变量编码：采用基于遗传算法的特征选择，直接采用二进制编码对特征变量进行编码，染色体的长度设为特征变量的个数； 

(5.13)种群初始化：随机初始化种群； 

(5.14)特征变量编码进行译码：将特征变量的二进制编码译码为实际的特征变量，得到此时的特征参数集FP_{sel_GA_ini}； 

(5.15)建立LSSVM模型：按步骤(5.5)的方法建立LSSVM模型，采用步骤(5.9)存储的模型参数和特征参数集； 

(5.16)读取误差评价训练样本数据； 

(5.17)计算个体适应度值：采用LSSVM校正模型中的均方根误差作为种群个体的适应度函数； 

(5.18)记录特征参数集：记录适应度最好的特征参数集为FP_{sel_GA_opt}； 

(5.19)达到设定繁殖代数？判断GA是否达到设定的繁殖代数，若达到，转步骤(5.30)；否则，转步骤(5.20)； 

(5.20)种群复制：根据个体适应度值的大小对个体进行排序，选择最佳个体，随机排序后替换适应度较弱的个体； 

(5.21)种***叉：随机选择交配的父体和交叉点进行单点交叉； 

(5.22)种群变异：采用基本变异算子进行变异操作； 

(5.23)种群更新：采用更新后的种群替换步骤(5.13)产生的初始化种群，并转步骤(5.14)； 

(5.24)读取一组特征参数：读取网格内的一组特征参数FP_{sel_grid_ini}； 

(5.25)LSSVM校正模型：按步骤(5.5)的方法建立LSSVM模型，采用步骤(5.9)存储的模型参数和步骤(5.24)选择的特征参数； 

(5.26)读取误差评价训练样本； 

(5.27)记录均方根误差最小的特征参数集FP_{sel_grid_opt}； 

(5.29)全部网格搜索完毕？判断网格算法是否结束，是则转步骤(5.30)；否则，转步骤(5.20)，存储优化后的特征参数集步骤(5.24)； 

(5.30)存储优化后的特征参数集：选择由遗传算法选择的FP_{sel_GA_opt}或网格算法选择的FP_{sel_GA_opt}的特征参数集，并改写为FP_sel； 

(5.31)结束：完成LSSVM模型的最优特征选择； 

对料球比、矿浆浓度及充填率模型分别执行上述过程，得到各自的优选特征参数集FP_{sel_mbr}，FP_{sel_den}，FP_{sel_bmwr}可统一表示为FP_seli，详见下式： 

${\mathrm{FP}}_{\mathrm{sel}}i=\left\{(h_{{\mathrm{VLF}}_{\mathrm{sel}}i},h_{{\mathrm{VMF}}_{\mathrm{sel}}i},h_{{\mathrm{VLF}}_{\mathrm{sel}}i})(h_{{\mathrm{ALF}}_{\mathrm{sel}}i},h_{{\mathrm{AMF}}_{\mathrm{sel}}i},h_{{\mathrm{ALF}}_{\mathrm{sel}}i})\left(h_{I_{\mathrm{sel}}i}\right)\right\}---\left(29\right)$

式中，i＝1，2，3分别对应料球比、矿浆浓度及充填率的特征选择参数。 

步骤六：确定球磨机的负荷参数料球比、矿浆浓度、充填率 

球磨机负荷参数预测模型的建立流程见附图4，其详细步骤如下。 

(6.1)开始：采用三个球磨机负荷参数的特征参数集FP_{sel_i}得到料球比、矿浆浓度及充填率模型的输入变量，统一表示为X_i，详见下式： 

$X_i=[X_{V\_{\mathrm{LF}}_{\mathrm{sel}}i},X_{V\_{\mathrm{MF}}_{\mathrm{sel}}i},X_{V\_{\mathrm{HF}}_{\mathrm{sel}}i},X_{A\_{\mathrm{LF}}_{\mathrm{sel}}i},X_{A\_{\mathrm{MF}}_{\mathrm{sel}}i},X_{A\_{\mathrm{HF}}_{\mathrm{sel}}i},X_{I_{\mathrm{sel}}i}]$

$X_{V\_{\mathrm{LF}}_{\mathrm{sel}}i}=[t_{V\_\mathrm{LF}\_1},\·\·\·,t_{V\_\mathrm{LF}\_h_{{\mathrm{VLF}}_{\mathrm{sel}}i}}]$

$X_{V\_{\mathrm{MF}}_{\mathrm{sel}}i}=[t_{V\_\mathrm{MF}\_1},\·\·\·,t_{V\_\mathrm{MF}\_h_{{\mathrm{VMF}}_{\mathrm{sel}}i}}]$

$X_{V\_{\mathrm{HF}}_{\mathrm{sel}}i}=[t_{V\_\mathrm{HF}\_1},\·\·\·,t_{V\_\mathrm{HF}\_h_{{\mathrm{VHF}}_{\mathrm{sel}}i}}]$

$X_{A\_{\mathrm{LF}}_{\mathrm{sel}}i}=[t_{A\_\mathrm{LF}\_1},\·\·\·,t_{A\_\mathrm{LF}\_h_{{\mathrm{ALF}}_{\mathrm{sel}}i}}]---\left(30\right)$

$X_{A\_{\mathrm{MF}}_{\mathrm{sel}}i}=[t_{A\_\mathrm{MF}\_1},\·\·\·,t_{A\_\mathrm{MF}\_h_{{\mathrm{AMF}}_{\mathrm{sel}}i}}]$

$X_{A\_{\mathrm{HF}}_{\mathrm{sel}}i}=[t_{A\_\mathrm{HF}\_1},\·\·\·,t_{A\_\mathrm{HF}\_h_{{\mathrm{AHF}}_{\mathrm{sel}}i}}]$

$X_{I_{\mathrm{sel}}i}=[t_{I\_1},\·\·\·,t_{I\_h_{I_{\mathrm{sel}}i}}]$

(6.2)选择误差惩罚参数和核参数的调整范围：方法同步骤(5.3)； 

(6.3)调整误差惩罚参数和核参数：方法同步骤(5.4)； 

(6.4)建立球磨机负荷参数的模型：方法同步骤(5.5)，建立的料球比、矿浆浓度及充填率模型如下： 

$y_i=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{ki}}{k}_i(x_i,x_{\mathrm{ki}})+b_i,k=1,\·\·\·,M---\left(31\right)$

式中，i＝1，2，3分别对应料球比、矿浆浓度及充填率； 

(6.5)读取误差评价训练样本数据：方法同步骤(5.6)； 

(6.6)记录误差评价结果和参数：方法同步骤(5.7)； 

(6.7)参数调整是否已经达到范围上限：方法同步骤(5.8)； 

(6.8)记录误差评价最好的模型参数：方法同步骤(5.9)； 

(6.9)确定模型：根据步骤(6.8)中选择的模型参数，确定模型训练结果，确定软测量模型； 

(6.10)读入验证样本数据：读入准备用于模型验证的一组样本数据； 

(6.11)模型验证：采用步骤(6.9)确定的模型按照步骤(5.7)中的均方根误差； 

(6.12)模型是否满足精度？精度满足要求，则转步骤(6.14)结束模型的建立，采用建立的模型对料球比、矿浆浓度、充填率进行软测量；否则，转步骤(6.13)； 

(6.13)重新收集构造样本数据：验证精度不能满足软测量的需要，需要增加试验次数，重新构造训练样本，转到步骤二； 

(6.14)结束。 

本发明的优点：利用计算机***和检测仪表提供的球磨机筒体振动、振声及电流信号，实现了球磨机负荷参数(料球比、矿浆浓度、充填率)软测量。本发明通过测量灵敏度高、抗干扰性强的球磨机筒体振动信号，增加了模型的灵敏度和预测精度。同时，本发明提出的 对筒体振动及振声信号的低、中、高三个频段分别进行非线性特征提取的方法，结合了球磨机的研球磨机理及振动与振声信号的产生机理，充分提取了信号中的球磨机负荷参数信息，克服了振动及振声信号的时域特征难以提取、频域变量超高维及常用特征频段求和方法的不足。而且，本发明提出的特征选择方法，同时解决了非线性特征数量(主元个数)难以合理确定和最小二乘-支持向量机模型的输入变量维数过多造成的训练速度慢和过拟合的问题，提高了模型的鲁棒性和预测性能。该方法有助于实现粉磨过程的优化控制和稳定运行，降低球磨机的电耗和钢耗，提高球磨机的台时处理量，提高了企业的经济效益。 

附图说明：

图1本发明一种球磨机负荷参数软测量方法的结构图； 

图2本发明一种球磨机负荷参数软测量方法的总的程序流程图； 

图3本发明一种球磨机负荷参数软测量方法的非线性特征选择流程图； 

图4本发明一种球磨机负荷参数软测量方法的建模过程流程图； 

图5(a)本发明一种球磨机负荷参数软测量方法的振动信号频域图； 

图5(b)本发明一种球磨机负荷参数软测量方法的振声信号频域图。 

具体实施方式：

以实验室内小型球磨机***为例。该球磨机的筒体尺寸为Φ460mm×460mm，内壁装有锰钢衬板，球磨机中部开口，可以随时停机、改变球磨机内矿物量和钢球量。该球磨机所装钢球的最大重量为80kg，矿物处理量为10kg/h。球磨机转速为53转/分钟。实验过程中使用Φ30mm、Φ20mm和Φ15mm三种锰钢钢球。实验中处理的矿石为经过破碎的铜矿石，颗粒尺寸小于6mm。 

如图1和2所示，一种球磨机负荷参数软测量方法基于KPCA-LSSVM方法，包括如下步骤： 

步骤一：数据采集 

按实验设计方案在球磨机内添加一定质量的钢球、矿石和水，待均匀混合后，启动球磨机，记录研磨过程的振动、振声及电流信号，达到设定的采集时间后，停止球磨机并清洗，准备下一样本的采集。 

步骤二：时域滤波 

筒体振动信号采用频率范围为100～12000Hz的带通滤波器；振声信号则采用4000Hz的低通滤波器；电流信号则采用均值法滤波。除采集湿磨条件下的数据外，还采集了球磨机空转，空砸，水磨和干磨状态下的数据。 

步骤三：振动及振声信号的时频转换： 

采用改进的平均周期图法(Welch)求得球磨机旋转每周的振动及振声信号的功率谱密度(PSD)，取全部周期的均值为最后结果。某一样本振动及振声信号如图5(a)和图5(b)所示。 

步骤四：对振动及振声信号的频谱数据和时域电流数据进行非线性特征提取： 

4.1对振动及振声信号的频谱数据数据进行非线性特征提取： 

分析球磨机空转、干磨、水磨及湿磨时的频域信号特征，可得球磨机筒体振动信号的低频率段为100～1800Hz、中频段为1800～4000Hz，高频段为4000～11000Hz。对训练样本的三个频段在100～1800，1800～4000，4000～7500三个频段采用以RBF核的KPCA做非线性特征提取，其结果如表2所示。其第一个主元在低、中、高频段的方差贡献率分别为87.42％，96.89％及98.92％。 

表2筒体振动信号不同频段进行KPCA分析的结果 

球磨机振声信号可以分为1～800，800～2500，2500～4000三个频段。对训练样本的三个频段以RBF核的KPCA做非线性特征提取，其结果如表3所示。其第一个主元在低、中、高频段的方差贡献率分别为46.72％，42.31％及65.5％，可见振声信号的灵敏度低于筒体振动信号。 

表3 振声信号不同频段进行KPCA分析的结果 

4.2对电流数据进行非线性特征提取； 

对电流信号进行非线性特征提取后，其第一主元的方差贡献率达到97.5％。用提取后的信号建立与充填率的模型，发先精度比提取前增加了一倍，表明了之间非线性的存在。 

步骤五：对振动、振声、电流非线性特征的组合数据进行特征选择； 

结合模型输入变量选择模型，球磨机负荷参数料球比、矿浆浓度、充填率模型的特征参数分别为mbr_FP_{sel_opt}＝{(4，2，1)(3，1，1)(1)}，density_FP_{sel_opt}＝{(4，2，1)(1，1，6)(1)}，bmwr_FP_{sel_opt}＝{(4，2，1)(2，1，6)(1)}。 

步骤六：确定球磨机负荷参料球比、矿浆浓度、充填率 

按上步选定的模型输入的特征参数，分别建立料球比、矿浆浓度、充填率针最小二乘-支持向量机软测量模型，误差惩罚参数和核参数选为(28510.7，15)、(19700，447)、(30000，83)。在该参数下，三个模型的参数α和b的值分别为： 

料球比模型： 

α＝[-0.811 0.631 -1.71 -1.29 -11.8 6.49 -3.04 3.14 -2.68 0.0436 1.82 3.875.39]；b＝0.966； 

矿浆浓度模型： 

α＝[0.514 -0.685 2.94 -41.2 0.517 18.7 2.81 -11.1 40.3 16.8 2.63 -54.89 23.8-0.719]；b＝0.686。 

充填率模型： 

α＝[0.693 -6.99 25.3 3.11 -11.4 -0.516 1.55 -26.8 -10.1 -3.66 36.2 -10.43.25]；b＝0.289。 

采用均方根误差作为模型精度的统计结果，其计算按下式进行： 

$\mathrm{RMSE}=\sqrt{\frac{1}{N}\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{e_k}^2}=\sqrt{\frac{1}{N}\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{({\hat{y}}_k-y_k)}^2}$

统计结果见表4。该表同时给出了建模的特征参数集及基于PCA-LSSVM及选择全部主元的建模结果。 

表4 不同建模方法的球磨机负荷参数软测量建模结果 

本实例是基于实验室球磨机***的球磨机负荷参数软测量，其料球比、矿浆浓度及充填率的变化范围大于工业现场球磨机，但是其测量精度仍然较高。因此，本发明即基于多源数据特征融合的球磨机负荷参数软测量方法能够有效的提取振动、振声及电流信号中的非线性特征并对特征参数进行有效的选择。本发明可以根据球磨机运行过程中产生的筒体振动信号、振声信号及电流信号实时估计球磨机负荷参数(料球比、矿浆浓度及充填率)，精度较一般的主元分析方法平均提高一倍，其均方根误差分别为0.0889，0.1006和0.09553。因此，本发明是一个具有很好实用价值和推广前景的球磨机负荷参数(料球比、矿浆浓度及充填率)测量手段。 

Claims (1)

1.一种球磨机负荷参数软测量方法，其特征在于：包括步骤如下

步骤一：采集球磨机筒体振动、振声及电流数据；

通过硬件支撑平台采集如下信号：球磨机筒体的振动加速度信号X_V；球磨机研磨区域下方的振声信号X_A；球磨机驱动电动机的电流信号X_I；

步骤二：对采集的球磨机振动、振声及电流数据进行滤波；

采用有限冲激响应滤波器分别对振动及振声信号进行带通滤波和低通滤波，滤波器公式见下式：

$y_1\left(m\right)=\underset{n_1=0}{\overset{N_1-1}{\mathrm{\Σ}}}h\left(n_1\right)x(m-n_1)$

式中， h(n₁)——滤波器的系数；

N₁——滤波器的长度；

y₁(m)——滤波器在第m点的输出；

x(m-n₁)——滤波器的第(m-n₁)个输入；

n₁——采样时刻n₁=0,1,…N₁-1；

m ——当前采样时刻；

采用均值滤波法对电流信号进行滤波，公式见下式：

$y_2\left(m\right)=\frac{1}{N^2}\underset{n_2=0}{\overset{N_2-1}{\mathrm{\Σ}}}{x}_{n_2}$

式中，y₂(m)——滤波器的输出，即N₂次采样的算术平均值；

——第n₂次的采样值；

N₂——采样点数；

步骤三：对滤波后的振动及振声数据进行时频转换；

采用改进的平均周期图法求球磨机旋转整周期的球磨机筒体振动及振声信号的功率谱密度；按如下公式：

$\tilde{P}\left(k\right)=\frac{1}{R}\underset{r=1}{\overset{R}{\mathrm{\Σ}}}{P}^r\left(k\right)r=1,...,R$

$P^r\left(k\right)=\frac{1}{N_3}{\left|X_{N_3}\left(k\right)\right|}^2$

式中，

——总功率谱；

P^r(k)——第r段的功率谱；

R——分段的段数；

N₃——X数字序列的数据点数；

k——为求得的频率点数；

步骤四：对振动及振声信号的频谱数据和时域电流数据进行非线性特征提取；

4.1对振动及振声信号的频谱数据进行非线性特征提取；

将湿式球磨机频域信号划分为低频段，中频段及高频段，其分别对应着球磨机的固有频率段、主要冲击频率段、次要冲击频率段，将球磨机筒体振动信号在低、中、高频的数据分别记做X_{V_LF_org}、X_{V_MF_org}、X_{V_HF_org}；振声信号在低、中、高频的非线性特征分别记做X_{A_LF_org}、X_{A_MF_org}、X_{A_HF_org}，这些频谱数据存在高维和共线性的问题，采用核主元分析分别对X_{V_LF_org}、X_{V_MF_org}、X_{V_HF_org}、X_{A_LF_org}、X_{A_MF_org}、X_{A_HF_org}进行降维和非线性特征提取；

对振动及振声信号共六个频段提取非线性特征均按如下公式进行提取：

$t_h=(v^h\·\mathrm{\Φ}\left(X\right))=\underset{k=1}{\overset{M_1}{\mathrm{\Σ}}}{\mathrm{\α}}_k^{\left(h\right)}\tilde{K}(X_k,X)$

式中：k表示第k个样本，范围从1到M₁；t_h——第h个主元；α_k——协方差阵的特征向量；Φ——非线性变换函数；Φ(X)表示输入样本数据X在高维特征空间的非线性映射；v^h——高维空间中，第h个主成份的特征矢量；K——M₁×M₁阶核矩阵；h——主元的个数；

——中心化处理后的核矩阵；

将球磨机筒体振动信号在低、中、高频的频谱变量X_{V_LF_org}、X_{V_MF_org}、X_{V_HF_org}提取的非线性特征分别记做X_{V_LF}、X_{V_MF}、X_{V_HF}；球磨机振声信号在低、中、高频的频谱变量X_{A_LF_org}、X_{A_MF_org}、X_{A_HF_org}提取的非线性特征分别记做X_{A_LF}、X_{A_MF}、X_{A_HF}，其值分别用下式表示：

$X_{V\_\mathrm{LF}}=[t_{V\_\mathrm{LF}\_1},t_{V\_\mathrm{LF}\_2},...,t_{V\_\mathrm{LF}\_h_{\mathrm{VL}{F}_{\max }}}]$ ;

$X_{V\_\mathrm{MF}}=[t_{V\_\mathrm{MF}\_1},t_{V\_\mathrm{MF}\_2},...,t_{V\_\mathrm{MF}\_h_{\mathrm{VM}{F}_{\max }}}]$ ；

$X_{V\_\mathrm{HF}}=[t_{V\_\mathrm{HF}\_1},t_{V\_\mathrm{HF}\_2},...,t_{V\_\mathrm{HF}\_h_{\mathrm{VH}{F}_{\max }}}]$ ；

$X_{A\_\mathrm{LF}}=[t_{A\_\mathrm{LF}\_1},t_{A\_\mathrm{LF}\_2},...,t_{A\_\mathrm{LF}\_h_{\mathrm{AL}{F}_{\max }}}]$ ；

$X_{A\_\mathrm{MF}}=[t_{A\_\mathrm{MF}\_1},t_{A\_\mathrm{MF}\_2},...,t_{A\_\mathrm{MF}\_h_{\mathrm{AM}{F}_{\max }}}]$ ；

$X_{A\_\mathrm{HF}}=[t_{A\_\mathrm{HF}\_1},t_{A\_\mathrm{HF}\_2},...,t_{A\_\mathrm{HF}\_h_{\mathrm{AH}{F}_{\max }}}]$ ；

式中，——振动及振声信号在频域内的低中高三个频段进行KPCA降维和特征提取后能够选择的最大主元个数；

——表示振动或振声信号的前

个主元的方差累积之和；

M₂=[V,A]——表示振动（V）和振声信号（A）；

N₄=[L,M,H]——分别表示频域信号的低(L)、中(M)、高(H)频段。

4.2对电流数据进行非线性特征提取；

采用与步骤4.1中相同的方法即核主元分析提取时域球磨机电流信号中的非线性特征，提取后的非线性特征记做X_I，其值可用下式表示：

$X_I=[t_{I\_1},t_{I\_2},...,t_{I\_h_{I\max }}]$

式中，t——主元，即提取的非线性特征；h_Imax——是时域电流信号进行KPCA非线性特征提取后的最大主元个数，其值采用如下公式确定：

$h_{I\max }={\arg }_{h_I}\min \{{\mathrm{CPV}}_{h_I}\≥99\%\}$ ；

式中，

——表示电流信号前

个主元的方差累积之和；

步骤五：对振动、振声、电流非线性特征提取后的融合数据进行特征选择；

将非线性特征输入变量中每类的最佳变量个数即KPCA的主元个数的集合称为特征参数集，并统一表示为： ${\mathrm{FP}}_{\mathrm{sel}}=\left\{(h_{{\mathrm{VLF}}_{\mathrm{sel}}},h_{{\mathrm{VMF}}_{\mathrm{sel}}},h_{{\mathrm{VLF}}_{\mathrm{sel}}})(h_{{\mathrm{ALF}}_{\mathrm{sel}}},h_{\mathrm{AM}{F}_{\mathrm{sel}}},h_{{\mathrm{ALF}}_{\mathrm{sel}}})\left(h_{I_{\mathrm{sel}}}\right)\right\}$ ，采用LSSVM模型选择误差惩罚参数和核参数的调整范围，获取存储优化后的特征参数集，最后确定料球比优选特征参数集、矿浆浓度优选特征参数集及充填率优选特征参数集，

特征参数集的确定采用如下特征参数优化模型：

式中，RMAE_val——最小二乘-支持向量机模型的最小评价误差；

——振动及振声信号在频域内的低中高三个频段进行KPCA降维和特征提取后能够选择的最大主元个数；

——表示振动或振声信号的前

个主元的方差累积之和；

M₂=[V,A]——表示振动（V）和振声信号（A）；

N₄=[L,M,H]——分别表示频域信号的低(L)、中(M)、高(H)频段；

——是时域电流信号进行KPCA非线性特征提取后的最大主元个数，一般取1；

n₃———输入变量的行数即样本数；

y———磨机负荷参数的真值；

———磨机负荷参数的估计值即模型的计算值；

步骤六：确定球磨机的负荷参数料球比、矿浆浓度、充填率；

通过料球比、矿浆浓度及充填率模型的输入变量x_i，采用LSSVM模型，确定的料球比、矿浆浓度及充填率：

$Y_i=\underset{k=1}{\overset{M_1}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{ki}}{k}_i(x_i,x_{\mathrm{ki}})+b_{i}k=1,...,M_1$ ，i=1,2,3分别对应料球比、矿浆浓度及充填率，α_ki、b_i为模型的系数；

得到球磨机的负荷参数料球比、矿浆浓度和充填率。

Images