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

Abstract

本发明公开了一种磨机负荷参数软测量方法，本发明的方法面向筒体振动和振声信号的多尺度和非稳态特性，基于不同视角，采用多种不同的信号分解技术将原始筒体振动信号和振声信号分解为系列子信号。将选择的子信号频谱和原始信号频谱作为多源多尺度信息构建磨机负荷参数软测量模型。采用基于自适应遗传算法(AGA)和分支定界(BB)算法的全局优化选择性集成核偏最小二乘(GOSENKPLS)优化选择候选子模型和选择性集成模型(SEN模型)的结构参数和学习参数，实现对多源多尺度信号的有效选择性融合。本发明可以提高磨机负荷参数软测量的精确性。

Description

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

技术领域

本发明涉及软测量领域，具体涉及一种磨机负荷参数软测量方法。

背景技术

球磨机是一类重型封闭连续工作的旋转机械设备。它们广泛的应用于磨矿过程，并且是这些生产过程的瓶颈设备(参见文献T.Y.Chai,“Operationaloptimizationandfeedbackcontrolforcomplexindustrialprocesses,”ActaAutomaticaSinica,Vol.39,No.11,pp.1744-1757,May.2013.)。球磨机内部负荷参数(如料球比(MBVR)、磨矿浓度(PD)、充填率(CVR))的过高或过低会导致磨机过负荷或欠负荷(参见文献P.Zhou,T.Y.Chai,H.Wang,“Intelligentoptimal-settingcontrolforgrindingcircuitsofmineralprocessing,”IEEETransactionsonAutomationScienceandEngineering,Vol.6,No.4,pp.730-743,April.2009.)。磨机过负荷会导致“磨机吐料”、“磨机出口粒度变粗”、“堵磨”，甚至“磨矿过程停产”；磨机欠负荷会导致磨机“空砸”，导致能源浪费、增加钢耗和研磨设备损坏。可见，磨机负荷参数与产品质量和产量相关，过高或过低甚至会导致研磨设备损坏。因此，它们是实现磨矿过程的优化运行和控制的关键过程参数。磨机负荷参数难以直接检测，其主要原因是球磨机封闭连续旋转的工作特点导致检测仪表不能被安装或无损害的使用。特殊设计的Sensomag检测装置(参见文献P.Keshav,B.D.Haas,B.Clermont,A.Mainza,M.Moys,“Optimisationofthesecondaryballmillusinganon-lineballandpulploadsensor-TheSensomag,”MineralsEngineering,Vol.24,No.3-4,pp.325-334,Feb.2011.)直接安装在磨机内部，其冲击防护和电力供给问题很难实现。磨机负荷参数也难以通过机理模型分析求得，主要原因：一是磨机入口给料存在的难以检测的频繁波动；二是磨矿过程的复杂物理、物理化学和力学过程。尽管磨机的研磨机理不是很清晰，但是物料破碎和研磨的实现主要依赖于钢球负荷的冲击行为。这个过程产生强烈的振动和振声信号。因此，这些信号被领域专家广泛应用于估计磨机负荷参数。

通常，磨机内部装载数以万计的钢球。理论上，这些钢球表面被矿浆覆盖并且分层排列。它们以不同的幅值和周期直接或间接的冲击磨机筒体。由这些冲击力导致的磨机筒体振动信号必定是由具有不同时间尺度的子信号组成。筒体振动信号具有多组分和非稳态特性。由于振动是振声信号的主要来源之一，振声信号的组成也许更为复杂。在工业实践中，领域专家基于多源信息通过他们的“人脑模型”估计他们所熟悉的磨机负荷及负荷参数来确保生产过程的安全。最常用的信号是磨机研磨区域附近的振声信号。优秀的领域专家能够“听音估计”磨机负荷参数。研究表明人耳在本质上是一组带通滤波器。因此，在某种角度上，“听音估计”过程可以被认为是一个认知建模过程，其主要包括：基于人耳的带通滤波、基于人脑的特征提取、基于专家经验的估计模型构建。本质上，这是个基于选择性融合多源多组分子信号的建模过程。

在时域内，有价值信息被隐含在宽带随机噪声信号“白噪声”内(参见文献Y.Zeng,E.Forssberg,“Monitoringgrindingparametersbyvibrationsignalmeasurement-aprimaryapplication,”MineralsEngineering,Vol.7,No.4,pp.495-501,April.1994.)，构建磨机负荷参数软测量模型的首要工作是对原始筒体振动和振声信号进行快速傅里叶变换(FFT)。这些频谱在本申请中称之为单尺度频谱。但是，传统的FFT是建立在被处理信号是稳态和线性的基础上的(参见文献Y.G.Lei,Z.J.He,Y.Y.Zi,“ApplicationoftheEEMDmethodtorotorfaultdiagnosisofrotatingmachinery,”MechanicalSystemsandSignalProcessing,Vol.23,No.4,pp.1327-1338,May.2009.)。短时傅里叶变换、小波分析、Wigner-Ville分布、进化谱等基于FFT的非稳态信号分析方法用于改进FFT分析的全局表达能力，但这些方法都不能对原始信号进行自适应分解。Huang带领的研究团队提出采用经验模态分解(EMD)技术处理非线性非稳态信号，该方法具有完全、正交、局部和自适应的优点(参见文献N.E.Huang,“TheempiricalmodedecompositionandtheHilbertspectrumfornonlinearandnon-stationarytimeseriesanalysis,”ProceedingsoftheRoyalSocietyofLondonAMathematicalPhysical&EngineeringSciences,Vol.454,No.1971,pp.903-995,Mar.1998.)。但是，在分解过程中可能会产生一些虚假的子信号模态，即模态混叠现象。借助噪声辅助分析技术，集合EMD(EEMD)可以有效克服这一问题(参见文献Z.H.Wu,N.E.,Huang,“EnsembleempiricalmodedecompositionforhighfrequencyECGnoisereduction,”AdvancesinAdaptiveDataAnalysis,vol.55,No.4,pp.193-201,May.2009.)。EEMD已经广泛应用在轴承故障诊断、结构模态辨识等领域。EMD和EEMD能够将原始信号按照频率由高到底分解为系列内禀模态函数(IMFs)，即将子信号按频率分布的角度进行分解。另外一种方法，Hilbert振动分解(HVD)能够从子信号能量分布的角度对原始信号进行自适应分解(参见文献M.Feldman,“Time-varyingvibrationdecompositionandanalysisbasedonHilberttransform,”JournalofSoundandVibration,2006,Vol.295,No.3-5,pp.518-530,Aug.2006)，即HVD分解后的子信号按能量从高到底进行排列。因此，这些不同分解算法获得的子信号可看作是源于多个视图的多源信息。但是，这些分解算法都只能够获得有限数量的有价值子信号。研究表明，基于IMFs的FFT变换能够较精确的确定IMFs的频率分布(参见文献V.K.Rai,A.R.Mohanty,“BearingfaultdiagnosisusingFFTofintrinsicmodefunctionsinHilbert-Huangtransform,”MechanicalSystemsandSignalProcessing,Vol.21,No.6,pp.2607-2165,Aug.2007.)。对应于单尺度频谱，这些子信号的频谱在本文中称之为多尺度频谱。

基于单尺度频谱的软测量方法包括：基于最小二乘-支持向量机(LS-SVM)的单一模型，基于偏最小二乘(PLS)和核偏最小二乘KPLS(KPLS)的集成模型(参见文献J.Tang,T.Y.Chai,L.J.Zhao,H.Yue,X.P.Zheng,“Ensemblemodelingforparametersofball-millloadingrindingprocessbasedonfrequencysepctrumofshellvibration,”ChineseControlTheory&Applications(Chinese),Vol.29,No.2,pp.183-201,May.2012.以及文献J.Tang,T.Y.Chai,L.J.Zhao,W.Yu,H.Yue,“Softsensorforparametersofmillloadbasedonmulti-spectralsegmentsPLSsub-modelsandon-lineadaptiveweightedfusionalgorithm,”Neurocomputing,Vol.78,no.1,pp.38-47,Feb.2012.)，基于KPLS和自适应加权融合算法(AWF)的选择性集成(SEN)模型(SENKPLS)(参见文献J.Tang,T.Y.Chai,W.Yu,andL.J.Zhao,“Modelingloadparametersofballmillingrindingprocessbasedonselectiveensemblemulti-sensorinformation,”IEEETransactionsonAutomationScienceandEngineering,Vol.10,no.3,pp.726-740,Jul.2013.)，以及基于KPLS的在线集成(OLENKPLS)模型(参见文献J.Tang,T.Y.Chai,W.Yu,L.J.Zhao,“On-lineKPLSalgorithmwithapplicationtoensemblemodelingparametersofmillload,”ActaAutomaticaSinica,Vol.39,No.5,pp.471-486,May.2013.)。这些方法均基于选择或提取的单尺度频谱特征或特征子集构建软测量模型。这些方法在提高模型的可解释性、解释筒体振动和振声信号的组成成分，以及分析磨机内部研磨机理等方面较为困难。汤等人首先提出了基于多尺度频谱的筒体振动信号分析方法(参见文献J.Tang,L.J.Zhao,H.Yue,W.Yu,T.Y.Chai,“Vibrationanalysisbasedonempiricalmodedecompositionandpartialleastsquares,”ProcediaEngineering,Vol.16,pp.646-652,May.2011.)。多种基于EMD的单一模型和多模型软测量方法用于检测磨机负荷参数，但未取得较好的预测性能。现有模型采用的SENKPLS算法都是采用“操纵输入特征”的集成构造方法产生候选子模型。采用不同的单尺度频谱特征子集或不同的多尺度频谱特征单独的构建候选子模型。也就是说，每个候选子模型从其自身具有最佳预测精度的角度选择结构参数和学习参数。采用分支定界(BB)和自适应加权(AWF)算法进行集成子模型的选择和组合。类似的，一种新提出的SEN神经网络模型能够分别构建优化的候选子模型和选择优化的集成子模型及其权重(参见文献S.Soares,C.H.Antunes,R.Araújo,“Comparisonofageneticalgorithmandsimulatedannealingforautomaticneuralnetworkensembledevelopment,”Neurocomputing,Vol.21,No.9,pp.498-511,May.2013.)。但是，候选子模型的良好预测性能并不意味着最终的SEN模型具有较好的泛化能力。不同集成子模型间的差异性对提高SEN模型的泛化能力具有重要的贡献。最合理的建模策略也许是同时全局优化SEN模型和候选子模型的结构参数和学习参数。

由此，亟需一种更加准确的磨机负荷参数软测量方法。

发明内容

有鉴于此，本发明提供一种磨机负荷参数软测量方法，以更加精确地进行磨机负荷参数的软测量。

所述方法包括：

S100、在预定的磨机负荷下运行磨机，获取表征磨机振动和/或振声的训练样本特征信号；

S200、将所述训练样本特征信号分解为具有不同时间或能量尺度的多种不同训练样本子信号，每个训练样本子信号代表具有物理含义的单一振动模态；

S300、对所有训练样本子信号进行相关性分析，并对与原始信号的相关性度量高于预定阈值的子信号进行时域频域转换以获得元素相互独立的特征频谱向量；

S400、以核变量数量、核参数以及集成模型尺寸进行编码获取遗传算法的染色体，并以编码后的染色体为个体根据预定种群规模生成初始种群；

S500、对当前种群的染色体进行解码以获取每个个体对应的核变量数量、核参数以及集成模型尺寸；

S600、基于每个个体对应的核变量数量和核参数对所述特征频谱向量中的每一个频谱生成对应的候选子模型；

S700、根据每个个体对应的集成模型尺寸以及对应的所有候选子模型获取每个个体对应的候选选择性集成预测模型；

S800、计算当前种群对应的所有候选选择性集成预测模型的预测输出的均方根相对误差(RMSRE)作为适应度，并根据适应度对当前种群的个体进行排序，判断遗传代数是否达到预定换代数，如果未达到，则根据交叉概率和变异概率进行选择、交叉和变异操作以获得新一代种群的染色体，转向S500，否则执行S900；

S900、输出当前种群中适应度最小的个体对应的候选选择性集成预测模型作为磨机负荷软测量模型；

S1000、获取需要进行软测量的磨机的测试数据的特征频谱向量；

S1100、根据所述磨机负荷参数软测量模型计算所述测试数据的特征频谱向量对应的负荷参数。

本发明面向筒体振动和振声信号的多尺度和非稳态特性，基于不同视角，采用多种不同的信号分解技术将原始筒体振动和振声信号分解为系列子信号。将选择的子信号频谱和原始信号频谱作为多源多尺度信息构建磨机负荷参数软测量模型。采用基于自适应遗传算法(AGA)和分子定界(BB)优化算法的全局优化选择性集成核偏最小二乘(GOSENKPLS)优化选择候选子模型和SEN模型的结构参数和学习参数，实现对多源多尺度信号的有效选择性融合。本发明可以提高磨机负荷参数软测量的精确性。

附图说明

通过以下参照附图对本发明实施例的描述，本发明的上述以及其它目的、特征和优点将更为清楚，在附图中：

图1是本发明实施例的磨机系统以及配套的软测量系统的硬件结构示意图；

图2是磨机筒体中振动信号和振声信号的产生示意图；

图3是本发明实施例的磨机负荷参数软测量方法的流程图；

图4是本发明实施例的磨机负荷参数软测量方法的信号流图；

图5A-图5C分别是本发明实施例的磨机负荷参数软测量方法中对于不同的参数预测性能前5的候选子模型和被选择次数前5的集成子模型统计结果；

图6是本发明实施例的磨机负荷参数软测量方法中不同候选子模型被选择的总次数的统计结果示意图。

具体实施方式

以下基于实施例对本发明进行描述，但是本发明并不仅仅限于这些实施例。在下文对本发明的细节描述中，详尽描述了一些特定的细节部分。对本领域技术人员来说没有这些细节部分的描述也可以完全理解本发明。为了避免混淆本发明的实质，公知的方法、过程、流程、元件和电路并没有详细叙述。

图1是本发明的磨机系统以及配套的软测量系统的硬件结构示意图。如图1所示，两段式磨矿回路(GC)广泛的应用于选矿过程，在磨矿回路的第一端，通常包括顺序连接的料仓1、给矿机2、湿式预选机3、磨机4和泵池5。水力旋流器6连接在泵池5和湿式预选机3之间，使得较粗粒度的部分作为底流被返回磨机再磨。新给矿、新给水以及周期性的添加钢球，与水力旋流器的底流一起进入磨机4(通常为球磨机)。矿石在磨机4内被钢球冲击和磨剥为较细的颗粒，与磨机4内水混合得到的矿浆连续的流出磨机，进入泵池5。通过向泵池5内注入新水对矿浆进行稀释，并将这稀释的矿浆以一定的压力注入水力旋流器6，进而这些泵入水力旋流器的矿浆被分离为两部分：包含较粗粒度的部分作为底流进入磨机再磨；其余部分则进入二段磨矿(GCII)。

同时，为了进行负荷参数的软测量，振动信号获取装置7和振声信号获取装置8分别与磨机4结合设置以获取振动信号和振声信号，数据处理装置9根据检测获得的振动信号和振声信号进行数据处理软测量获取负荷参数。在本说明书的以下描述中，以筒体振动信号和振声信号为例来对本发明的技术方案进行说明。但是，本领域技术人员容易理解，也可以筒体振动信号、轴承振动信号、磨机研磨区域振声信号中的至少一种。只要获取的样本信号和测试信号可以用于表征磨机振动和/或振声即可。

球磨机筒体振动系统可被看作是受到外界的持续扰动并对外界的扰动无反馈作用的一类机械系统。磨机负荷的检测与识别问题可以归为一类特殊的动载荷识别问题，但识别的不是冲击力，而是与此冲击力直接相关的磨机负荷参数及钢球、物料和水负荷(参见文“磨矿过程磨机负荷软测量方法研究”[D],汤健，东北大学,中国,2012.)。

工业实际中，磨机筒体内的物料、钢球和水负荷连续运动且相对比较稳定(尤其是安装有自动加球机的磨矿过程)，故磨机负荷可看作是磨机筒体结构的一部分。因此，加入磨机负荷后的球磨机筒体振动系统是由磨机筒体及其内部的物料、钢球和水负荷组成的一个受交变应力载荷作用的新机械结构体。该新机械结构体的物理参数(质量、刚度及阻尼等)中包含磨机负荷信息。磨机负荷不同时，该结构体的各种物理参数也同时发生变化。因此，磨机负荷的识别问题就包含了机械振动中的系统识别和动载荷识别两个方面的问题。常规的动载荷识别方法很难用于识别磨机负荷。

图2是磨机筒体中振动信号和振声信号的产生示意图。磨机筒体可近似看作板壳结构体，理想情况下，磨机内部钢球的运动可以分为上升阶段和下降阶段。在上升阶段，钢球随着磨机筒体的旋转从落回点B到达脱离点A；在下降阶段，钢球沿着抛物运行从脱离点A到达落回点B。因此，钢球运动可以分为四个阶段：下落、冲击、研磨和滑动。通过判断是否对物料产生破碎作用，磨机内部空间可以分为五个区域：抛落区(I)、冲击破碎区(II),研磨区(III)、滑动区(IV)、和死角区(V)。实际上，磨机内部存在数以万计的钢球。理论上，这些钢球分层排列，并且不同层的钢球以不同的冲击力同时抛落。矿石硬度、粒度分布和矿石其它属性也影响这些冲击力的大小。这些具有不同幅值和频率的冲击力相关叠加。而且，磨机筒体自身的质量不平衡，磨机安装或使用过程中导致的偏置，以及其它原因均可能激起筒体振动。这些信号相互耦合和叠加形成筒体振动信号，可用下式表示：

$x_V^t=\underset{j_V=1}{\overset{J_V}{\Σ}}{x}_{j_V}^t---\left(1\right)$

其中，和J_V分别表示第j_Vth子信号和筒体振动子信号的数量。

筒体振动的声辐射，即结构噪声，是振声信号的主要组分部分。由于磨机筒体是声学上的硬壁面即强反射面，大部分噪声在筒体内连续反射，形成混合响声场。该混合声场噪声的一部分被矿浆及筒体吸收，另外一部分通过筒体和固定磨机衬板的螺栓传输至磨机外部，构成空气噪声。而且，我们通常测量到的振声信号还包括邻近磨机和其它研磨设备产生的环境噪声。显然，其灵敏度要弱于筒体振动信号。振声信号至少由结构噪声、空气噪声和环境噪声组成，其组成成分比筒体振动信号更复杂，可用下式表示：

$x_A^t=\underset{j_A=1}{\overset{J_A}{\Σ}}{x}_{j_A}^t---\left(2\right)$

其中，和J_A分别表示第j_A个子信号和振声子信号的数量。现有技术中，工业现场经常使用的磨机振声信号是由振动辐射噪声即筒体结构噪声、磨机内部混合声场传输至磨机外部的空气噪声、与磨机负荷无关的环境噪声等部分组成，其同样具有类似筒体振动信号的多组分特性。优秀领域专家可以凭借自身经验“听音识别”所熟悉的特定磨机的负荷及其内部参数状态。但是，专家仅能对振声信号进行辨识，而空气噪声和环境噪声会对此构成较大的干扰。而振动信号并不能为专家辨识。

图3是本发明实施例的磨机负荷参数软测量方法的流程图。如图3所示，所述方法包括：

S100、在预定的磨机负荷参数下运行磨机，获取表征磨机振动和/或振声的训练样本特征信号。

S200、将所述训练样本特征信号分解为具有不同时间或能量尺度的多种不同训练样本子信号，每个训练样本子信号代表具有物理含义的单一模态。

S300、对所有训练样本子信号进行相关性分析，将与原始信号相关性较强的子信号进行时域频域转换以获得相互独立的特征频谱向量。

S400、以核变量数量、核参数以及集成模型尺寸进行编码获取遗传算法的染色体，并以编码后的染色体为个体根据预定种群规模生成初始种群。

S500、对当前种群的染色体进行解码以获取每个个体对应的核变量数量、核参数以及集成模型尺寸。

S600、基于每个个体对应的核变量数量和核参数对所述特征频谱向量中的每一个频谱生成对应的候选子模型。

S700、以每个个体对应的集成模型尺寸以及对应的所有候选子模型获取每个个体对应的候选选择性集成预测模型。

S800、计算当前种群对应的每个候选选择性集成预测模型的预测输出的均方根相对误差(RMSRE)作为每个个体的适应度，并根据适应度对当前种群的个体进行排序，判断遗传代数是否达到预定换代数，如果未达到，则根据交叉概率和变异概率进行选择、交叉和变异操作以获得新一代种群的染色体，然后转向S500，否则执行S900。

S900、输出当前种群中适应度最低的个体对应的候选选择性集成预测模型作为磨机负荷软测量模型。

S1000、获取需要进行软测量的磨机的测试数据的特征频谱向量。

S1100、根据所述磨机负荷参数软测量模型计算所述测试数据的特征频谱向量对应的负荷参数。

图4是本发明实施例的磨机负荷参数软测量方法的信号流图。以下结合图3和图4对本发明实施例的一个优选实施方式做进一步地说明。

其中，在图4中，和表示长度为N的磨机筒体振动及振声信号；和表示第和个磨机筒体振动子信号；和表示第和个振声子信号；其中， 和表示不同分解算法获得的筒体振动和振声子信号的数量； 和是经相关性分析后选择的子信号； 和是选择的筒体振动和振声子信号的数量； 是选择的筒体振动和振声子信号的特征频谱子向量；{z₁,...,z_j,...,z_J}表示重新编号后的多源多尺度的特征频谱向量；表示基于z_j构建的候选子模型；表示第j个候选子模型的预测输出；和分别表示第j_sel个集成子模型及其权重系数；表示集成子模型的数量，即集成模型尺寸；和表示候选子模型的核参数和核变量(KLVs)的数量；表示磨机负荷软测量模型的输出。

在步骤S200，收集到的磨机的振动信号和振声信号被输入到图4所示的频谱转换部分进行频谱转换。首先，以并行的方式分别通过经验模态分解(EMD)、集合经验模态分解(EEMD)、希尔伯特振动分解(HVD)、局部均值分解(LMD)中的一种或几种多组分信号分解方法将所述样本特征信号分解为具有不同时间尺度的样本子信号。在图4中，分别对振动信号进行EMD分解、EEMD分解和HVD分解，获得对应的训练样本子信号和同时分别对振声信号进行EMD分解、EEMD分解和HVD分解，获得对应的训练样本子信号和

其中，EMD和EEMD将原始信号按时间尺度大小依次分解。EEMD能够克服EMD的模态混叠效应，因需要人工设定分解参数丧失了自适应分解特性。HVD按照能量分布由强到弱对信号进行自适应分解。理论上，这些子信号都具有不同的物理含义。因此，它们被看作基于不同视图的多源多尺度信息。原始筒体振动和振声信号的分解过程可用如下公式表示：

${\left\{x_V^t\right\}}_{n=1}^N=\underset{j_{VDECOM}=1}{\overset{J_{VDECOM}}{\Σ}}{\left\{x_{j_{VDECOM}^{all}}^t\right\}}_{n=1}^N+{\left(r_{J_{VDECOM}}\right)}_{VDECOM}---\left(3\right)$

${\left\{x_A^t\right\}}_{n=1}^N=\underset{j_{ADECOM}=1}{\overset{J_{ADECOM}}{\Σ}}{\left\{x_{j_{ADECO}^{all}}^t\right\}}_{n=1}^N+{\left(r_{J_{ADECOM}}\right)}_{ADECOM}---\left(4\right)$

其中，DECOM＝{EMDorEEMDorHVD}表示所采用的分解算法；和分别表示第个和第个子信号；和表示相应的残差信号。

在步骤S300，对上述六组训练样本子信号进行子信号相关性分析以选取获得与原始信号具有较强相关性的子信号，实现对应于训练样本子信号空间的降维操作。

具体地，可以选择与原始信号的相关性度量高于预定阈值的子信号。选择与原始筒体振动和振声信号相关的子信号是构建可靠的软测量模型的基础。在本实施例中，采用如下公式计算这些子信号和原始信号间的相关系数：

${\mathrm{\λ}}_{j_{VDECOM}}=\frac{\mathrm{\Σ}(x_V^t-{\stackrel{\‾}{x}}_V^t)(x_{j_{VDECOM}^{all}}^t-{\stackrel{\‾}{x}}_{j_{VDECOM}^{all}}^t)}{\mathrm{\Σ}{(x_V^t-{\stackrel{\‾}{x}}_V^t)}^2\mathrm{\Σ}{(x_{j_{VDECOM}^{all}}^t-{\stackrel{\‾}{x}}_{j_{VDECOM}^{all}}^t)}^2}---\left(5\right)$

${\mathrm{\λ}}_{j_{ADECOM}}=\frac{\mathrm{\Σ}(x_A^t-{\stackrel{\‾}{x}}_A^t)(x_{j_{ADECOM}^{all}}^t-{\stackrel{\‾}{x}}_{j_{ADECOM}^{all}}^t)}{\mathrm{\Σ}{(x_A^t-{\stackrel{\‾}{x}}_A^t)}^2\mathrm{\Σ}{(x_{j_{ADECOM}^{all}}^t-{\stackrel{\‾}{x}}_{j_{ADECOM}^{all}}^t)}^2}---\left(6\right)$

其中， $x_{j_{VDECOM}^{all}}^t={\left\{x_{j_{VDECOM}^{all}}^t\right\}}_{n=1}^N,x_{j_{ADECOM}^{all}}^t={\left\{x_{j_{ADECOM}^{all}}^t\right\}}_{n=1}^N;$ 和表示信号均值。

当真实相关系数为0时，不相关值假设检验值 具有与随机抽取得到的观测值相同的概率(优选地，可以采用预定的相关性计算函数获得所述概率)。通常，如果和都小于0.05，相关系数和是显著的。采用如下准则选择子信号：

${\mathrm{\&xi;}}_{j_{DECOM}^{sel}}=\left\{\begin{array}{ccccc}1& if& p_{j_{VDECOM}}& or& p_{j_{ADECOM}}\&le;p_{threshold}\\ 0& others& p_{j_{VDECOM}}& or& p_{j_{ADECOM}}<p_{threshold}\end{array}\right.---\left(7\right)$

其中，p_threshold是预设定阈值。

选择的子信号，并将选择的筒体振动和振声子信号分别标记为和其中和表示被选择的子信号的数量。

容易理解，虽然上述实施例中采用相关性分析来进行子信号的选择和降维，但是，也可以采用其它相关性分析和互信息分析中的一种或几种方法的组合对所述训练样本子信号进行降维。

同时，在步骤S300中，在降维后，对降维后选择获得的振动和振声子信号进行时域频域转换，将位于时域的振动和振声子信号转换为频谱信号获得训练样本子信号的特征频谱向量。

具体地，时/频转换过程可表示为：

${\left\{x_{j_{VDECOM}^{sel}}^t\right\}}_{n=1}^N\stackrel{FFT}{\&RightArrow;}{z}_{j_{VDECOM}^{sel}}---\left(8\right)$

${\left\{x_{j_{ADECOM}^{sel}}^t\right\}}_{n=1}^N\stackrel{FFT}{\&RightArrow;}{z}_{j_{ADECOM}^{sel}}---\left(9\right)$

为简洁表示，将这些频谱进行改写和重新编号：

$\begin{array}{c}z=\{z_{1_{VEMD}^{sel}}\mathrm{...},z_{j_{VEMD}^{sel}},\mathrm{...},z_{J_{VEMD}^{sel}},z_{1_{VEEMD}^{sel}}\mathrm{...},z_{j_{VEEMD}^{sel}},\mathrm{...},z_{J_{VEEMD}^{sel}},\\ z_{1_{VHVD}^{sel}}\mathrm{...},z_{j_{AEEMD}^{sel}},\mathrm{...},z_{J_{AEEMD}^{sel}},z_{l_{AHVD}^{sel}}\mathrm{...},z_{j_{AHVD}^{sel}},\mathrm{...},z_{J_{AHVD}^{sel}},\\ z_{1_{AEEMD}^{sel}}\mathrm{...},z_{j_{AEEMD}^{sel}},\mathrm{...},z_{J_{AEEMD}^{sel}},z_{1_{AHVD}^{sel}}\mathrm{...},z_{j_{AHVD}^{sel}},\mathrm{...},z_{J_{AHVD}^{sel}}\}\\ =\{z_1,\mathrm{...},z_j,\mathrm{...},z_J\}\end{array}---\left(10\right)$

其中，

$J=J_{VEMD}^{sel}+J_{VEEMD}^{sel}+J_{VHVD}^{sel}+J_{AEMD}^{sel}+J_{AEEMD}^{sel}+J_{AHVD}^{sel}---\left(11\right)$

也即，J是全部频谱的数量，也即特征频谱向量的维度。

获得的训练样本子信号的特征频谱向量被输入到全局优化建模部分以计算训练磨机负荷参数软测量模型。

从步骤S400至步骤S900，本发明实施例的方法基于提取获得的特征频谱向量进行全局优化建模。

具体地，将全部建模样本表示为{Z,y},并将第j个频谱改写为则有如下关系式：

$Z={\left\{Z_j\right\}}_{j=1}^J={\left\{{\left\{{\left(z_j\right)}_l\right\}}_{l=1}^k\right\}}_{j=1}^J---\left(12\right)$

这些频谱可以看作为来自不同视图的多源信息，基于这些频谱的选择性集成(SEN)模型构建本质上是个多源信息融合过程，其目的是模拟基于领域专家的人脑估计模型，可描述为如下优化问题：

${\mathrm{maxJ}}_{RMSRE}={\mathrm{\&theta;}}_{th}-\sqrt{\frac{1}{k}\underset{l=1}{\overset{k}{\&Sigma;}}{\left(\frac{y^l-\underset{j_{sel}=1}{\overset{J_{sel}^{sel}}{\&Sigma;}}{w}_{j_{sel}}^{sel}{\hat{y}}_{j_{sel}}^l}{y^l}\right)}^2}---\left(13\right)$

$s.t.\left\{\begin{array}{c}2\&le;J_{sel}^{sel}\&le;J,1\&le;j_{sel}\&le;J_{sel}^{sel}\\ \underset{j_{sel}=1}{\overset{J_{sel}^{sel}}{\&Sigma;}}{w}_{j_{sel}}^{sel}=1,0<w_{j_{sel}}^{sel}<1\\ {\hat{y}}_{j_{sel}}^l=f_{j_{sel}}^{sel}(z_{j_{sel}}^{sel},{\mathrm{\&gamma;}}_{KPLS}^{sel},h_{KPLS}^{sel}),\\ {\left\{f_{j_{sel}}^{sel}(\&CenterDot;)\right\}}_{j_{sel}=1}^{J_{sel}^{sel}}\&Element;{\left\{f_j^{can}(\&CenterDot;)\right\}}_{j=1}^J\\ {\hat{y}}_{j_{sel}}={\left\{{\hat{y}}_{j_{sel}}^l\right\}}_{l=1}^k,\\ {\hat{y}}_j={\left\{{\hat{y}}_j^l\right\}}_{l=1}^k\\ {\left\{{\hat{y}}_{j_{sel}}\right\}}_{j_{esl}=1}^{J_{sel}^{sel}}\&Element;{\left\{{\hat{y}}_j\right\}}_{j=1}^J\end{array}\right.$

其中，θ_th是预先设定的阈值。

本文中，和表示集成尺寸和集成子模型的权重系数，也可看作为SEN模型的结构参数和学习参数；和分别是核变量(KLVs)数量和核参数，也可以看作为基于核偏最小二乘(KPLS)的候选子模型的结构参数和学习参数。为了得到式(13)描述的优化目标，至少需要以下三个步骤：

(1)、选择和构建具有合适预测精度和差异性的候选子模型。

(2)、从候选子模型集合中选择集成子模型

(3)、采用权重系数合并这些集成子模型。

为了实现全局优化，本发明实施例采用双层优化的建模策略：其中外层优化采用自适应遗传算法(AGA)选择和内层优化采用分支定界法(BB)和自适应加权融合算法(AWF)来选择和通过为全部候选子模型优化选择和可以在子模型的预测精度和多样性间取得均衡。通过AGA给定的固定的内层优化策略的搜索时间比在“2～J”之间搜索最优局部SEN模型的效率明显提高。因此，内层和外层优化实现了集成，本发明将上述方法称之为GOSENKPLS。

具体地，在步骤S400，采用二进制编码系统表征遗传算法(GA)中的染色体。共有三个参数需要编码，将其长度分别记为 和则编码染色体的全部长度为：

$Len=Len\left(h_{KPLS}^{sel}\right)+Len\left({\mathrm{\&gamma;}}_{KPLS}^{sel}\right)+Len\left(J_{sel}^{sel}\right)---\left(14\right)$

其中，和是的最大值和最小值； h_max是KLV数量的最大值，采用下式获得：

$h_{max}=\min \left\{{\left\{rank\left(K_j\right)\right\}}_{j=1}^J\right\}---\left(15\right)$

其中，K_j是基于Z_j的核矩阵。在本文中，SEN模型是基于高维小样本数据构造，故取h_max＝k。

对于预定的AGA的种群规模J_GA，则可用和表示种群内的个体的建模数。由此，通过在上述范围内以随机方式选择和组合可以获得初始种群，种群中有J_GA个不同的染色体作为种群个体。

在步骤S500，对于当前的种群，对每个个体对应的染色体进行解码，也即，将二进制编码形式的向量转换为对应的核变量数量、核参数以及集成模型尺寸。

具体地，以种群内的第j_GA个体为例，采用如下公式解码建模参数：

${\mathrm{\&gamma;}}_{KPLS}^{j_{GA}}={\mathrm{\&gamma;}}_{KPLS}^{\min }+\frac{{\mathrm{\&gamma;}}_{KPLS}^{\max }-{\mathrm{\&gamma;}}_{KPLS}^{\min }}{2^{Len\left({\mathrm{\&gamma;}}_{KPLS}\right)-1}}\&times;\underset{i_{bit}=1}{\overset{Len\left({\mathrm{\&gamma;}}_{KPLS}\right)}{\&Sigma;}}(2^{i_{bit}}\&CenterDot;B_{i_{bit}})---\left(16\right)$

$h_{KPLS}^{j_{GA}}=1+\frac{h_{\max }-1}{2^{Len\left(h_{KPLS}\right)-1}}\&times;\underset{i_{bit}=1}{\overset{Len\left(h_{KPLS}\right)}{\&Sigma;}}(2^{i_{bit}}\&CenterDot;B_{i_{bit}})---\left(17\right)$

$J_{sel}^{j_{GA}}=2+\frac{J-2}{2^{Len\left(J_{sel}\right)-1}}\&times;\underset{i_{bit}=1}{\overset{Len\left(J_{sel}\right)}{\&Sigma;}}(2^{i_{bit}}\&CenterDot;B_{i_{bit}})---\left(18\right)$

其中，B∈{0,1}，也即，为染色体第i_bit位的二级制编码。

在步骤S600，本实施例以KPLS方法基于当前种群内每一个个体对应的核变量数量、核参数来针对特征频谱向量中的J个不同的频谱分别生成J个不同候选子模型。由此，生成的候选子模型数量总数为J_GA*J。

以当前种群中第j_GA个体为例，构建基于第j个频谱的KPLS候选子模型。首先采用核技巧基于核参数将频谱映射到高维特征空间：

$K_j(Z_j,{\mathrm{\&gamma;}}_{KPLS}^{j_{GA}})=\mathrm{\&Phi;}{\left({\left(z_j\right)}_l\right)}^T\mathrm{\&Phi;}\left({\left(z_j\right)}_m\right),l,m=1,2,\mathrm{...},k---\left(19\right)$

然后，对核矩阵进行中心化：

${\tilde{K}}_j(Z_j,{\mathrm{\&gamma;}}_{KPLS}^{j_{GA}})=(I-\frac{1}{k}{1}_k{1}_k^T)\&CenterDot;K_j(Z_j,{\mathrm{\&gamma;}}_{KPLS}^{j_{GA}})\&CenterDot;(I-\frac{1}{k}{1}_k{1}_k^T),---\left(20\right)$

其中，I是k维的单位阵，1_k是长度为k值为1的向量。接着，在该高维空间执行线性偏最小二乘(PLS)算法。最终构建完成KPLS模型，并将基于频谱Z_j构建的候选子模型的输出记为：

${\left({\hat{y}}_j^{can}\right)}_{j_{GA}}={\left(f_j^{can}\right({\tilde{K}}_j\left(Z_j,{\mathrm{\&gamma;}}_{KPLS}^{j_{GA}}\right),h_{KPLS}^{j_{GA}}\left)\right)}_{j_{GA}}---\left(21\right)$

对当前种群的第j_GA个个体，其所有候选子模型选择相同的核参数和核变量数量，全部候选子模型可表示为：

$F_{j_{GA}}={\left\{{\left(f_j^{can}\right({\tilde{K}}_j\left(Z_j,{\mathrm{\&gamma;}}_{KPLS}^{j_{GA}}\right),h_{KPLS}^{j_{GA}}\left)\right)}_{j_{GA}}\right\}}_{j=1}^J.---\left(22\right)$

进一步，对于种群内的全部J_GA个个体，其候选子模型的集合可以表示为可见，候选子模型的全部数量为J_GA*J。

在步骤S700，所述建模方法进入到内循环中。在本实施例，基于分支定界法(BB)和自适应加权融合算法(AWF)来选择和由此，通过选择和加权集成，可以对于种群内的J_GA个个体获取到J_GA个候选的选择性集成预测模型，以供后续步骤进行进一步处理。

具体地，步骤S700包括：

S710、对于每个个体，根据对应的集成模型尺寸J_sel从对应的J个候选子模型中选取J_sel个模型选择权重参数最佳的候选子模型作为集成子模型。

S720、根据如下公式计算集成子模型对应的权重

$w_{j_{sel}}^{j_{GA}}=1/\left({\left({\mathrm{\&sigma;}}_{j_{sel}}^{j_{GA}}\right)}^2\underset{j_{sel}=1}{\overset{J_{sel}^{j_{GA}}}{\&Sigma;}}\frac{1}{{\left({\mathrm{\&sigma;}}_{j_{sel}}^{j_{GA}}\right)}^2}\right)---\left(23\right)$

其中，为集成子模型输出值的方差。

以种群第j_GA个个体为例说明集成子模型的选择和合并流程。基于预定的权重参数计算方法，可以将优化选择集成子模型当作从原始的J个子模型中选择个最佳特征子模型的问题。

这样，将最优特征子集选择算法(此处采用BB)和最佳的多传感器系统融合算法(此处采用AWF)组合后用于解决这一问题。集成子模型的权重采用上述的(23)式计算。

集成子模型和候选子模型的预测输出间的关系可用下式表示：

${\left\{{\hat{y}}_{j_{sel}}^{j_{GA}}\right\}}_{j_{sel}=1}^{J_{sel}}\&Element;{\left\{{\hat{y}}_j^{j_{GA}}\right\}}_{j=1}^J.---\left(24\right)$

可见，前者为后者的子集。种群中第j_GA个个体的局部SEN模型的输出可用下式计算：

${\hat{y}}_{j_{GA}}=\underset{j_{sel}=1}{\overset{J_{sel}^{j_{GA}}}{\&Sigma;}}{w}_{j_{sel}}^{j_{GA}}{\hat{y}}_{j_{sel}}^{j_{GA}},---\left(25\right)$

进一步，可将J_GA个种群的全部局部SEN型的输出标记为：

${\left\{{\hat{y}}_{j_{GA}}\right\}}_{j_{GA}=1}^{J_{GA}}={\left\{{\left\{{\left({\hat{y}}_l\right)}_{j_{GA}}\right\}}_{l=1}^k\right\}}_{j_{GA}=1}^{J_{GA}}---\left(26\right)$

由此，可以为种群中的每一个个体获得局部最优的选择性集成预测模型。

应理解，本实施例的选择性集成操作还可以采用人工神经网络(ANN)、支持向量机(SVM)、最小二乘-支持向量机(LS-SVM)、岭回归及其核版本算法、单隐层前馈神经网络(SLFN)、随机向量泛函连接网络(RVFL)等一种或几种候选子模型构建算法，基于遗传算法、分支定界等一种或几种优化算法的集成子模型选择算法，以及基于自适应加权融合(AWF)、基于误差信息熵加权、基于线性或非线性模型等一种或几种集成子模型加权系数计算方法。

在步骤S800，对获得候选选择性集成预测模型进行排序，并通过选择、变异和交叉操作来获得新一代的种群。

具体地，将AGA的最大换代数记为G_max，并将局部SEN的预测输出的均方根相对误差(RMSRE)作为适应度函数。以种群的第j_GA个个体为例，采用如下公式计算适应度：

$F_{fit}^{j_{GA}}=\sqrt{\left\{\underset{l=1}{\overset{k}{\&Sigma;}}{\left(\frac{y_l-{\left({\hat{y}}_l\right)}_{j_{GA}}}{y_l}\right)}^2\right\}/k}---\left(27\right)$

其中，y_l是真值，是第j_GA个个体对应的候选选择性集成预测模型对第l个样本的预测值。

将全部J_GA个个体的适应度记为并将其从低到高进行排列。若AGA的遗传代数未到达G_max，通过种群选择、交叉和变异过程获得新的种群；然后重复进行上述S500至S700的“解码”、“候选子模型构建”、“局部SEN模型构建”过程。若遗传代数达到G_max，则进入步骤S900，也即，将具有最小适应度的局部SEN模型作为磨机负荷参数软测量模型。

交叉概率p_c和变异概率p_m采用如下公式自适应计算得到:

$\left\{\begin{array}{cc}{p}_c=p_{c1}-\frac{(p_{c1}-p_{c2})(f_{max}-f_{l\arg er})}{f_{max}-f_{ave}},& f_{l\arg er}\&GreaterEqual;f_{ave}\\ p_c=p_{c1},& f_{l\arg er}<f_{ave}\end{array}\right.---\left(28\right)$

$\left\{\begin{array}{cc}{p}_m=p_{m1}-\frac{(p_{m1}-p_{m2})(f_{max}-f_{l\arg er})}{f_{max}-f_{ave}},& f_{l\arg er}\&GreaterEqual;f_{ave}\\ p_m=p_{m1},& f_{l\arg er}<f_{ave}\end{array}\right.---\left(29\right)$

其中，p_c1＝0.9,p_c2＝0.6,p_m1＝0.1,p_m2＝0.001；和是全部种群中的平均和最大适应度；f_larger是两个种群中较大的适应度。应理解，p_c1,p_c2,p_m1,p_m2也可以根据需要设定为其它的预定值。

由此，可以快速获得全局最优的磨机负荷参数软测量模型。

在测量阶段，也即步骤S1000和S1100中，根据建模过程中的降维方式进行降维以及时域频域转换获得测试数据的特征频谱向量，并进而根据建模获取的磨机负荷参数软测量模型计算负荷参数。

基于本发明实施例的方法进行了实验。实验在小型实验球磨机上进行，实验设备的参数如表1所示。

表1实验设备参数

实验所用矿石负荷在使用前被粉碎至直径小于6mm。采用直径为30、20和15mm的钢球作为研磨介质，总量40Kg，配比为3:4:3，。磨机实验共进行了4组：①第一组：料负荷为10kg，水负荷从5kg增加到40kg；②第二组：料负荷为20kg，水负荷从2kg增加到20kg；③第三组：水负荷为2kg，料负荷从10kg增加到20kg；④第四组：水负荷为10kg，料负荷从22kg增加到50kg。

实验中，对振动信号和振声信号进行分解后，子信号与原始信号间的相关性分析如表2所示。

表2分解结果的相关性分析

在实验中，将p_threshold取为0.05。通常，若可以认为这些子信号与原始信号是强相关的。表2表明只有部分筒体振动子信号的相关系数高于0.05，分别是EMD的Vsub6-Vsub10，EEMD的Vsub6-7和Vsub9，以及HVD的Vsub8；同时，表明：EEMD在一定程度上解决了EMD的虚假模态问题，提高了分解结果的相关性；若从相关性的角度考虑，HVD是最佳方法。表2同时表明：所有的振声信号的子信号均与原始信号强相关；所有的振声信号的值高于筒体振动信号。这也从另外一个角度表明筒体振动信号的组分较为复杂。实际上，筒体振动信号记录的是磨机旋转一周的振动模态，而振声信号记录的仅仅是研磨区域的噪声信号。

同时，实验中所关心的不只是子信号与原始信号间的相关性，主要是从这些相关的子信号中选择与磨机负荷参数更相关的子信号。从另外的视角出发，本研究的目的是模拟领域专家如何从这些子信号中选择有价值信息。因此，这些的子信号的频谱用于构建SEN模型的候选子模型。由于原始信号代表着全局信息，故也将其用构建候选子模型。

用于构建SEN模型的信号包括：原始筒体振动信号，HVD的Vsub1-7和Vsub9-13，EMD的Vsub1-5，EEMD的Vsub1-4和Vsub7-8，原始振声信号，HVD的Asub1-11，EMD和EEMD的Asub1-12，总共有61个不同来源和不同尺度的信号。本文中，将它们从1到61顺序编号，可见前24和后37个分别源于筒体振动和振声信号；因此，候选子模型的数量是61，即J＝61。候选子模型的编号与选择的信号之间的对应关系如表3所示：

表3候选子模型的编号与选择的信号之间的对应关系

软测量模型的构建采用如下参数：建模参数的编码长度均为20；解码范围：为1-100,为1-12,为2-12；种群规模J_GA为20；最大遗传代数G_max为10；种群收敛的百分比为98％；种群初始化时的百分比为50％。为了克服AGA算法的随机初始化特性，建模算法重复运行20次。

基于本发明实施例的SEN模型的学习参数和预测性能如表4所示。

表4SEN模型的学习参数和预测性能

表4表明：(1)、对于RMSREs的均值，在3个磨机负荷参数中：PD最小，而MBVR最大，这与之前的研究结果相符。

(2)、对于集成尺寸的均值，同RMSREs的均值类似，对于PD，CVR和MBVR分别是3.10、4.0和4.75。

(3)、对于核变量数量的均值，按从小到达的顺序对于PD，MBVR和CVR分别是5.70、8.9和9.80。

(4)、对于核参数的均值，MBVR最小，只有PD和CVR的三分之一。因此，对于不同的SEN模型，需要选择不同的学习参数。

基于这些不同的学习参数，在20次重复运行中，预测性能前5的候选子模型和被选择次数前5的集成子模型的统计结果如图5A-5C所示，不同候选子模型被选择的总次数如图6所示，由图5A-5C以及图6可知：

(1)对于MBVR：图5A表明53#候选子模型在20次重复运行中具有最佳的预测精度。通过对比，表明53#是EEMD的Asub4，其频带宽度为500-1000Hz；表明信号分解在提高单一模型的预测性能是有效的。其它的被选择次数较多的候选子模型是15#，25#，50#，52#和3#，其中只有15#和3#源于筒体振动信号，故与MBVR相关的信息更多蕴藏于振声信号。这也解释了为什么振声信号可以较为准确的检测干式磨机负荷(即MBVR)。从信号的产生机理上上，振声信号的主要来源是筒体振动。图6表明选择次数最多的集成子模型也是53#；排在第1位的被选择次数最多的集成子模型是第25#，即原始振声信号，原因是其蕴含着全局信息。基于图6可知，集成子模型15#，52#和53#具有相同的被选择次数，即EMD的VSub2、EEMD的Asub3和Asub4对预测MBVR具有最大贡献；3#和15#集成子模型的选择次数是15和14，故HVD的VSub2也具有较大的贡献。综上，本文所提方法有效选择性融合了有价值信号对MBVR进行预测。

(2)对于PD：图5B表明14#候选子模型总是具有最小的预测误差，但是与15#、23#和3#子模型的差异并不大，只有1-2次的差距；对比表4可知，这些集成子模型来自于不同的分解算法，如3#源于HVD、14#和15#源于EMD、23#源于EEMD；而且，全部子模型均源于筒体振动信号。图5B中，集成子模型的选择结果也是全部源于筒体振动信号。图6表明，15#子模型被选择了20次，即EMD的VSub2对预测PD具有最大贡献；我们可以推测该子信号主要是由球负荷的冲击引起的；因钢球表面覆盖着矿浆，故PD与这个子信号直接相关；同时14#和23#子模型被选择了14和9次，1#,4#和11#子模型也被选择了3次。可见，所提方法有效的选择了与PD相关的有价值子信号。

(3)对于CVR：图5C表明23#候选子模型具有最小的预测误差，即EEMD的VSub6在预测CVR时远胜于其它子信号；与PD只选择与筒体振动相关的信号不同，51#(EEMD的ASub2)候选子模型也具有较好的预测性能。图5C表明，集成子模型只选择了筒体振动相关的信号用于构建SEN模型；排在第1位被选择的集成子模型是15#，但是该子模型在图5C中却没有出现，事实上，15#子模型的预测性能经常排在第6位。这个现象表明，同时考虑精度和多样性，选择适当的集成子模型构建SEN模型是很重要的。图6表明，3#和23#子模型被选择了20次，4#和15#子模型被选择了17次；这些子模型分别来自于筒体振动信号的EMD、EEMD和HVD的分解结果。上述结果验证了所提方法的有效性。

从不同信息来源的视角出发，集成子模型被选择次数的统计结果如表5所示。

表5基于不同信息来源视角的集成子模型被选择的次数统计

表5表明：对于3个不同的磨机负荷参数，共选择了237个不同来源的集成子模型；基于HVD、EMD和EEMD分解结果被选择的子模型所占的比例为26％、33％和31％，只有10％基于未分解的原始信号，故对筒体振动及振声信号进行分解是合理的；并且，不同的磨机负荷参数对不同的分解算法的喜好不同，如对于MBVR、PD和CVR，44％,65％和46％的集成子模型源于EEMD,EMD和HVD，尤其是PD选择的基于EMD的子模型达到了65％。因此，本发明实施例的方法能够构建有效的SEN模型，并且在一定程度上模拟了领域专家的磨机负荷估计过程。

本发明实施例从多个不同的视角，采用不同的多组分非稳态信号分解算法获得具有不同物理含义的子信号。例如，EMD和EEMD可以获得具有不同时间尺度的子信号，前者可以完全自适应分解原始信号，后者通过预设定参数可以克服EMD的模态混叠问题；HVD可以将原始信号从能量分布的角度依次分解；原始筒体振动和振声信号包含全局信息。虽然领域专家的磨机负荷估计行为过于复杂难以描述，这些多源多尺度信号也许有效地模拟了人耳的带通滤波和大脑的信号处理过程。本发明比现有的方法具有更佳有效的和更多的信息来源。同时，本发明实施例的方法能够选择优化的SEN模型和候选子模型的结构参数和学习参数。每个候选子模型均是从构建最终有效的SEN模型的视角出发的。工业现场的领域专家也是从全局视角出发估计磨机负荷。更够更好地模拟专家的判断逻辑，提供更加精确的磨机负荷参数软测量。

以上所述仅为本发明的优选实施例，并不用于限制本发明，对于本领域技术人员而言，本发明可以有各种改动和变化。凡在本发明的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种磨机负荷参数软测量方法，包括：

S100、在预定的磨机负荷下运行磨机，获取表征磨机振动和/或振声的训练样本特征信号；

S200、将所述训练样本特征信号分解为具有不同时间或能量尺度的多种不同训练样本子信号，每个训练样本子信号代表具有物理含义的单一振动模态；

S300、对所有训练样本子信号进行相关性分析，并对与原始信号相关性度量高于预定阈值的子信号进行时域频域转换以获得元素相互独立的特征频谱向量；

S400、以核变量数量、核参数以及集成模型尺寸进行编码获取遗传算法的染色体，并以编码后的染色体为个体根据预定种群规模生成初始种群；

S500、对当前种群的染色体进行解码以获取每个个体对应的核变量数量、核参数以及集成模型尺寸；

S600、基于每个个体对应的核变量数量和核参数对所述特征频谱向量中的每一个频谱生成对应的候选子模型；

S700、根据每个个体对应的集成模型尺寸以及对应的所有候选子模型获取每个个体对应的候选选择性集成预测模型；

S800、计算当前种群对应的所有候选选择性集成预测模型的预测输出的均方根相对误差(RMSRE)作为适应度，并根据适应度对当前种群的个体进行排序，判断遗传代数是否达到预定换代数，如果未达到，则根据交叉概率和变异概率进行选择、交叉和变异操作以获得新一代种群的染色体，转向S500，否则执行S900；

S900、输出当前种群中适应度最小的个体对应的候选选择性集成预测模型作为磨机负荷软测量模型；

S1000、获取需要进行软测量的磨机的测试数据的特征频谱向量；

S1100、根据所述磨机负荷参数软测量模型计算所述测试数据的特征频谱向量对应的负荷参数。

2.根据权利要求1所述的磨机负荷参数软测量方法，其特征在于，所述步骤S200包括：

通过经验模态分解(EMD)、集合经验模态分解(EEMD)、希尔伯特振动分解(HVD)、局部均值分解(LMD)中的一种或几种多组分信号分解方法将所述训练样本特征信号分解为具有不同时间或能量尺度的多种不同训练样本子信号。

3.根据权利要求1所述的磨机负荷参数软测量方法，其特征在于，所述步骤S500包括基于如下公式解码获取对应的核变量数量、核参数以及集成模型尺寸：

${\mathrm{\&gamma;}}_{KPLS}^{j_{GA}}={\mathrm{\&gamma;}}_{KPLS}^{\min }+\frac{{\mathrm{\&gamma;}}_{KPLS}^{\max }-{\mathrm{\&gamma;}}_{KPLS}^{\min }}{2^{Len\left({\mathrm{\&gamma;}}_{KPLS}\right)-1}}\&times;\underset{i_{bit}=1}{\overset{Len\left({\mathrm{\&gamma;}}_{KPLS}\right)}{\&Sigma;}}(2^{i_{bit}}\&CenterDot;B_{i_{bit}})$

$h_{KPLS}^{j_{GA}}=1+\frac{h_{max}-1}{2^{Len\left(h_{KPLS}\right)-1}}\&times;\underset{i_{bit}=1}{\overset{Len\left(h_{KPLS}\right)}{\&Sigma;}}(2^{i_{bit}}\&CenterDot;B_{i_{bit}})$

$J_{sel}^{j_{GA}}=2+\frac{J-2}{2^{Len\left(J_{sel}\right)-1}}\&times;\underset{i_{bit}=1}{\overset{Len\left(J_{sel}\right)}{\mathrm{\&Sigma;}}}(2^{i_{bit}}\&CenterDot;B_{i_{bit}})$

其中，和分别是第j_GA个个体对应的核参数、核变量数量以及集成模型尺寸，和分别是的最大值和最小值，Len(·)为括号内参数的编码长度，K_j是基于Z_j的核矩阵，Z_j为第j个频谱，J为特征频谱向量的维度，为染色体第i_bit位的二级制编码。

4.根据权利要求1所述的磨机负荷参数软测量方法，其特征在于，步骤S600包括：

对于当前种群中的每个个体，采用该个体对应的核变量数量和核参数基于核偏最小二乘(KPLS)算法对J个特征频谱分别训练获取J个候选子模型。

5.根据权利要求1所述的磨机负荷参数软测量方法，其特征在于，步骤S700包括：

S710、对于每个个体，根据对应的集成模型尺寸J_sel从对应的J个候选子模型中选取J_sel个模型选择权重参数最佳的候选子模型作为集成子模型；

S720、根据如下公式计算集成子模型对应的权重

$w_{j_{sel}}^{j_{GA}}=1/\left({\left({\mathrm{\&sigma;}}_{j_{sel}}^{j_{GA}}\right)}^2\underset{j_{sel}=1}{\overset{J_{sel}^{j_{GA}}}{\&Sigma;}}\frac{1}{{\left({\mathrm{\&sigma;}}_{j_{sel}}^{j_{GA}}\right)}^2}\right)$

其中，为集成子模型输出值的方差。

6.根据权利要求1所述的磨机负荷参数软测量方法，其特征在于，所述交叉概率和所述变异概率根据如下公式计算：

$\left\{\begin{array}{cc}{p}_c=p_{c1}-\frac{\left(p_{c1}-p_{c2}\right)\left(f_{\max }-f_{l\arg er}\right)}{f_{\max }-f_{ave}},& f_{l\arg er}\&GreaterEqual;f_{ave}\\ p_c=p_{c1},& f_{l\arg er}<f_{ave}\end{array}\right.$

$\left\{\begin{array}{cc}{p}_m=p_{m1}-\frac{\left(p_{m1}-p_{m2}\right)\left(f_{\max }-f_{l\arg er}\right)}{f_{\max }-f_{ave}},& f_{l\arg er}\&GreaterEqual;f_{ave}\\ p_m=p_{m1},& f_{l\arg er}<f_{ave}\end{array}\right.$

其中，p_c1,p_c2,p_m1,p_m2为预定值，和是种群的平均适应度和最大适应度；f_larger是两个个体中较大的适应度。

7.根据权利要求6所述的磨机负荷参数软测量方法，其特征在于，p_c1＝0.9,p_c2＝0.6,p_m1＝0.1,p_m2＝0.001。

8.根据权利要求1所述的磨机负荷参数软测量方法，其特征在于，所述表征磨机振动和/或振声的样本特征信号包括筒体振动信号、轴承振动信号、磨机研磨区域振声信号中的至少一种。

9.根据权利要求1所述的磨机负荷参数软测量方法，其特征在于，所述磨机为球磨机，所述负荷参数包括料球比、矿浆浓度和充填率中的一项或多项。