CN104932425A - Mill load parameter soft measurement method - Google Patents

Abstract

The invention discloses a mill load parameter soft measurement method, which comprises the steps of firstly decomposing mill cylinder vibrations and vibration sound signals into sub-signals (intrinsic mode functions (IMFs)) with different time scales and different physical meanings by adopting an ensemble empirical mode decomposition (EEMD) technology; then selecting of three kinds of features of the multi-scale IMF by adopting mutual information (MI) based adaptive feature selection approach, wherein the three kinds of features are the frequency spectrum, the marginal spectrum, and the mean and the standard deviation of the instantaneous amplitude and the frequency of Hilbert transformation; and finally, constructing a soft measurement model based on selective integration kernel partial least squares (KPLS) method on the basis of selected spectrum features and a training sample. A small ball mill based simulation experiment result shows that the method disclosed by the invention can effectively detect load parameters.

Description

Soft measurement method for load parameters of mill

Technical Field

The invention relates to the field of soft measurement, in particular to a soft measurement method for a load parameter of a mill.

Background

Optimal Operational control of the grinding process requires accurate detection of the load parameters within the mill (see document [1] p. zhou, t.y. chai, h.wang, "Intelligent operation-setting control for grinding circuits of the processing," IEEE Transactions on Automation Science and Engineering,6(2009) 730-. Tens of thousands of steel balls in the mill are arranged in layers, and the impact force of the steel balls in different layers to materials in the mill and a mill cylinder has different strength and period. Typically the measured cylinder vibration signal is a mixture of multiple sub-signals having different time scales. Barrel vibration is the primary source of mill vibro-acoustic signals. Thus, these mechanical vibration and vibro-acoustic signals have unsteady and multi-component characteristics. The outstanding field experts can effectively monitor the load state of the mill and the internal load parameters of partial mills by simultaneously considering various operating conditions and various source information. Research shows that the human ear can distinguish valuable information from the mill vibration sound signal. In fact, the human ear is a set of adaptive band-pass filters and the human brain has a multi-layered cognitive structure. The domain expert can extract valuable information from multiple source characteristics and multiple operating conditions to make decisions. Differences in field expertise and limited effort make it difficult to ensure that the mill is operating in an optimum load condition for a long period of time. For these situations, it is necessary to simulate the cognitive process of a field expert to build a soft measurement model of the mill load parameters.

In the time domain, valuable information in mill barrel vibrations and rattles is implicit in broadband random noise (see document [3] Y., Zeng, E.Forssberg, "Monitoring and marking parameters by noise signaling-a priority application," Minerals Engineering,1994,7(4): 495-501.). The modeling of mill load parameters based on mechanical vibration and vibro-acoustic signals requires attention to 3 sub-problems: self-adaptive decomposition of multi-component signals, self-adaptive selection of multi-source spectral features, and construction of a soft measurement model based on selection of multiple operating conditions.

Studies have shown that signal processing can simplify the selection and extraction process of features (see document [4] S.Shukla, S.Mishra, and B.Singh, "Power Quality Event Classification Under noise Conditioning EMD-Based De-noise Techniques," IEEE Transaction on Industrial information, 10(2014) 1044.). The mill load parameters are closely related to the Power Spectral Density (PSD) of the barrel vibration and vibro-acoustic signals (see document [5] J.Tang, L.J.ZHao, J.W.ZHou, H.Yue, T.Y.Chai, "empirical analysis of wet load based on vibration signals of laboratory-scale ball mill shell," Minerals Engineering,23(2010)720-730 "), but these spectral data typically contain thousands of features. Many dimension reduction algorithms are used to process data with this characteristic (see document [6] J.Tang, T.Y.Chai, W.Yu, L.J.ZHao, "Modeling load parameters of ball in grinding process based on selective sensitive multi-sensor information," IEEE Transactions on Automation Science and engineering,10(2013) 726-. Algorithms based on Mutual Information (MI) and Partial Least Squares (PLS) can efficiently identify these features (see document [6 ]). For efficient fusion of these spectral features, Soft measurement model methods based on integrated PLS, selective integration (SEN) and kernel PLS (kpls) have been reported (see document [7] j.tang, t.y.chai, l.j.zhao, w.yu, h.yue, "Soft sensor for parameters of small loaded based on porous segments and on-line adaptive weighted fusion," neuro-sampling, 78(2012)38-47. document [8] j.ta, t.y.chai, w.yu, l.j.zhao, "Featurext and selection based on porous fusion with parameter of P.n., (999999999999991) of synthetic samples. However, Fast Fourier Transform (FFT) is not suitable for processing of Mechanical vibration and vibro-acoustic signals with unsteady characteristics (see document [9] Y.G.Lei, Z.J.He, Y.Y.Zi, "Application of the EEMDmethod to rotor fault diagnosis of rotation machinery," Mechanical Systems and Signal processing,23(2009) 1327-1338.). Discrete wavelet transformation, Continuous Wavelet Transformation (CWT), wavelet packet transformation, and the like, have been widely used for Fault diagnosis of rotating machinery (see document [10] G.K. Singh, S.A.S.AlKazzaz, "Isolation and identification of dry finding faults in indexing machine using Fault transform," Triblogy International 42(2009) "849-, "simulation Analysis Based Interactive Diagnosis in indication mechanisms," transfer on Industrial information, 10(2014)340 "350. document [14] P.K.Kankar, S.C.Sharma, S.P.Harsha," Rolling element bearing surface using auto-correlation and contact wave transform, "Journal of simulation and Control,17(2011) 2081-. However, these methods are not capable of adaptively decomposing the multi-component signals faced herein, as a suitable mother wavelet must be selected for the CWT in the face of any particular practical problem. Empirical Mode Decomposition (EMD) technique obtains intrinsic mode functions (IMFs, also called subsignals) With Different time scales by adaptive decomposition (see document [15] N.E. Huang, Z.Shen, S.R. Long, "The empirical mode decomposition and The Hilbert specification for non-linear and non-static time Analysis," Proc.Royal Soc.London A,454(1998) 903) 995), and has been widely applied to rotating equipment fault diagnosis (see document [16] J.Faiz, V.G. vertical, and B.M.Ebranchi, "EMD-Based Analysis of Motors th Brookn Rouk for simulation of Analysis, IEEE model of Analysis, 16. Analysis of simulation of Analysis, 10. Analysis of Analysis, 12. Analysis of metals," IEEE 19. 12. simulation of Analysis of metals, "10. B.S. 12. simulation of Analysis of 12, IEEE 10. 16, 10(2014)1044-1054. document [18] R.Y.Li, D.He, "Rotational machine health monitoring and fault detection using EMD-based access establishment," IEEE Transaction on Instrumentation and Measurement,61(2012)990-1001 ]. Document [19] (V.K.rai, A.R.Mohanty, "Bearing fault diagnosis FFT of internal model functions in Hilbert-Huang transform," Mechanical Systems and Signal Processing,21(2007) 2607-; document [21] (L.J.Zoha, J.Tang, W.R.Zong, "Embedded based on empirical mode decomposition and partial least squares," Journal of Theoretical and Applied Information Technology,45(2012)179 and 191.) and document [22] (J.Tang, T.Y.Chai, Q.M.Cong, B.C.Yuan, etc., "Soft sensor application for modeling of small diameter based on empirical mode and selected sensitive algorithm," Acta Automatic site, 40(2014)1-14.) A model based on EMD is constructed, but these methods have a low prediction accuracy. The methods based on MI and PLS do not allow efficient selection of the spectral characteristics of IMF. The newly proposed integrated EMD (EEMD) technology overcomes the mode mixing problem of the EMD method (see document [23] Z. Wu, N.E. Huang, "Embedded empirical mode composition: a noise-associated Data Analysis method," Advances in Adaptive Data Analysis,1(2009) 1-41.). Studies have shown that the number of valuable IMFs is limited. Therefore, how to select valuable IMFs and their spectral characteristics at the same time is a second problem to be solved. Many of the SEN modeling approaches mentioned above employ an integrated construction approach of "manipulating input features". The methods can effectively select different frequency spectrum characteristic subsets from the perspective of multi-source information fusion, but are difficult to reflect the contribution of training samples representing different operating conditions.

In summary, attention needs to be paid to selective simultaneous fusion of multi-source multi-scale spectral features and multi-condition training samples. Selective integration (GASEN) based on genetic algorithm (see [24] Zhou ZH, Wu J and Tang W, "engineering core bean tholl," engineering intellectual integration, vol.137, No.1-2, pp.239-263,2002.) adopts a "sampling training sample" method to construct integration, adopts a Back Propagation Neural Network (BPNN) to construct candidate submodel and adopts a simple average method to merge the integrated submodel. The partial least squares (KPLS) method can overcome the defects of long training time, overfitting, difficulty in small sample modeling and the like of BPNN. The optimal observed values can be obtained when Adaptive Weighted Fusion (AWF) is oriented to a multi-sensor system (see document [25] L.Xu, J.Q.Zhang, Y.Yan, "A wave-based multisensor data fusion algorithm," IEEE transactions on instrumentation and Measurement,53(2004) 1539-.

Disclosure of Invention

In view of the above, the present invention provides a new soft measurement method for mill load parameters, which comprises three steps: constructing a SEN soft measurement model based on EEMD multi-scale signal adaptive decomposition, MI and PLS IMF and multi-type spectral feature adaptive selection, and AWF improved GASEN-KPLS (IGASEN-KPLS); and is constructed based on a soft measurement model.

The soft measurement method for the load parameters of the mill comprises the following steps:

s100, detecting and acquiring a sample vibration signal and a sample vibration sound signal of a cylinder of the grinding machine operating under known load parameters;

s200, decomposing the sample vibration signal and the sample vibration sound signal into J based on an Ensemble Empirical Mode Decomposition (EEMD) method_IMFSample sub-signals, each sample sub-signal representing a single vibrational mode having a physical meaning, J_IMFIs a predetermined value;

s300, extracting the spectral features of all the sub-signals, and acquiring a training sample set and a verification sample set based on the spectral features and corresponding load parameters, wherein the spectral features are obtained by selecting from a Marginal Spectrum (MSHT) based on Hilbert transform, a mean value and variance (MVHT) of instantaneous amplitude and frequency of the Hilbert transform and a Power Spectral Density (PSD) based on fast Fourier transform of the sub-signals;

s400, generating J training subsamples from a training sample set based on a Bootstrap algorithm;

s500, training J training sub-samples to obtain J candidate sub-models by adopting the same number of nuclear latent variables and nuclear parameters based on a Kernel Partial Least Squares (KPLS) algorithm;

s600, selecting all candidate submodels with corresponding model selection weight parameters larger than a model selection threshold value from the J candidate submodels according to the verification sample set to form an integrated submodel set (also called an integrated model), wherein the optimized weight parameters are obtained by optimizing a randomly generated initial weight parameter through a genetic algorithm with a minimized prediction error as a target;

s700, calculating and acquiring the integration weight of each sub-model in the integration sub-model set according to an AWF algorithm;

s800, obtaining the spectrum characteristics of the test data of the mill needing soft measurement;

and S900, calculating load parameters corresponding to the spectrum characteristics of the test data according to the integrated sub-model set and the corresponding integrated weights.

The soft measurement method provided by the embodiment of the invention has the advantages that: (1) valuable spectral features and training samples can be selected simultaneously; (2) the method is effectively suitable for a selective integrated modeling method of high-dimensional small sample data; (3) the valuable IMF frequency spectrum and the spectrum characteristics thereof are selected in a self-adaptive manner; (4) after the multi-component signal is adaptively decomposed into multi-scale sub-signals with different physical meanings, three different types of spectral features are extracted simultaneously. Based on the points, the soft measurement method for the load parameters of the mill according to the embodiment of the invention can have higher accuracy.

Drawings

The above and other objects, features and advantages of the present invention will become more apparent from the following description of the embodiments of the present invention with reference to the accompanying drawings, in which:

FIG. 1 is a schematic hardware configuration of the mill system of the present invention and the associated soft measurement system;

FIG. 2 is a flow chart of a method of soft measurement of mill load parameters in an embodiment of the present invention;

FIG. 3 is a flowchart of the EEMD-based sub-signal decomposition step of step S200 according to an embodiment of the present invention;

FIG. 4 is a flowchart of the spectral feature extraction step of step S300 according to an embodiment of the present invention;

FIG. 5 is a graph of MI threshold versus Root Mean Square Relative Error (RMSREs) for different mill loading parameters;

FIGS. 6a-6c are plots of maximum and minimum MI values and feature numbers, respectively, of selected MSHT, MVHT and PSD spectral features for an MBVR;

FIGS. 7a-7c are plots of maximum and minimum MI values and feature numbers, respectively, of selected MSHT, MVHT and PSD spectral features for a PD;

FIGS. 8a-8c are plots of maximum and minimum MI values and feature numbers, respectively, for selected MSHT, MVHT and PSD spectral features for a CVR;

FIG. 9 is a graph of the integrated dimensions of three different load parameters;

FIG. 10 is a graph of prediction error for three different load parameters.

Detailed Description

The present invention will be described below based on examples, but the present invention is not limited to only these examples. In the following detailed description of the present invention, certain specific details are set forth. It will be apparent to one skilled in the art that the present invention may be practiced without these specific details. Well-known methods, procedures, components and circuits have not been described in detail so as not to obscure the present invention.

Further, those of ordinary skill in the art will appreciate that the drawings provided herein are for illustrative purposes and are not necessarily drawn to scale.

Unless the context clearly requires otherwise, throughout the description and the claims, the words "comprise", "comprising", and the like are to be construed in an inclusive sense as opposed to an exclusive or exhaustive sense; that is, what is meant is "including, but not limited to".

In the description of the present invention, it is to be understood that the terms "first," "second," and the like are used for descriptive purposes only and are not to be construed as indicating or implying relative importance. In addition, in the description of the present invention, "a plurality" means two or more unless otherwise specified.

FIG. 1 is a schematic diagram of the hardware configuration of the mill system of the present invention and the associated soft measurement system. As shown in fig. 1, a two-stage Grinding Circuit (GC) is widely used in mineral separation processes and generally includes, at a first end of the grinding circuit, a silo 1, a feeder 2, a wet preselector 3, a mill 4 and a pump sump 5 connected in series. Hydrocyclone 6 is connected between pump sump 5 and wet preselector 3 so that the coarser fraction is returned as underflow to the mill for regrinding. Fresh feed, fresh feed water, and periodic addition of steel balls, along with the underflow of the hydrocyclone, enter the mill 4 (typically a ball mill). The ore is impacted and ground into fine particles in the mill 4 by the steel balls, and the ore pulp obtained by mixing the ore with water in the mill 4 continuously flows out of the mill and enters the pump pool 5. The pulp is diluted by injecting fresh water into the pump basin 5 and this diluted pulp is injected into the hydrocyclone 6 at a certain pressure, and the pulp pumped into the hydrocyclone is separated into two fractions: the fraction containing the coarser particle size enters the mill as underflow for regrinding; the rest part enters the second-stage grinding (GC II).

Meanwhile, in order to perform soft measurement of the load parameters, a vibration signal acquisition device 7 and a vibration sound signal acquisition device 8 are respectively arranged in combination with the mill 4 to acquire a vibration signal and a vibration sound signal, and a data processing device 9 performs data processing soft measurement according to the vibration signal and the vibration sound signal obtained by detection to acquire the load parameters.

The grinding productivity (i.e. grinding yield) is usually obtained by maximizing the optimum cyclic load, which is often determined by the load of the GCI. Overload of the mill can cause material discharge of the mill, coarsening of granularity of materials at the outlet of the mill, blockage of the mill and even production stop in the ore grinding process. The underload of the mill can lead to empty smashing of the mill, energy waste is caused, the loss of the steel ball is increased, and even the mill is damaged. Thus, mill load is a very important parameter. The accurate measurement of the internal load parameters of the ball mill is closely related to the guarantee of the product quality and the production efficiency in the ore grinding process and the safety of the production process. In industrial fields, field experts mostly rely on multi-source information and self experience to monitor the load state of the mill. A data driving soft measurement method based on vibration signals and vibration sound signals of a mill cylinder is commonly used for overcoming subjectivity and instability caused by expert inference mill load. However, the existing method is difficult to selectively fuse valuable information hidden in multi-scale spectral data and multi-operation working conditions (training samples) of mechanical vibration signals and vibro-acoustic signals at the same time.

The ball mill is heavy rotating machinery equipment, and mainly realizes grinding by means of impact and grinding and peeling of steel balls on materials. The impact force of the outermost steel ball on any point on the mill during one cycle of the mill rotation can be expressed as:

$\begin{array}{c}{\left(F_{\mathrm{bmw}}^{1\mathrm{stlayer}}\right)}_{\mathrm{period}}=[{\left(F_{\mathrm{bmw}}\right)}_{i1}^{1\mathrm{stlayer}},...,{\left(F_{\mathrm{bmw}}\right)}_{\mathrm{in}}^{1\mathrm{stlayer}},...,\\ {\left(F_{\mathrm{bmw}}\right)}_{g1}^{1\mathrm{stlayer}},...,{\left(F_{\mathrm{bmw}}\right)}_{\mathrm{gm}}^{1\mathrm{stlayer}},...,{\left(F_{\mathrm{bmw}}\right)}_{s1}^{1\mathrm{stlayer}},...,{\left(F_{\mathrm{bmw}}\right)}_{s1}^{1\mathrm{stlayer}},...]\end{array}---\left(1\right)$

wherein,andrepresenting the impact force of the impact, rub-off and slip phases. Obviously, these forces have different impact amplitudes and impact frequencies, just divided intoThe impact force of the vibration of the refiner cylinder at a certain moment is not reasonable. Thus, the length of the cylinder vibration signal analysis is at least one revolution of the mill cylinder.

ByThe induced barrel vibration was marked as:

${\left(x_V^t\right)}_{\mathrm{period}}^{1\mathrm{stlayer}}=[{\left(x_V^t\right)}_{i1}^{1\mathrm{stlayer}},...,{\left(x_V^t\right)}_{\mathrm{in}}^{1\mathrm{stlayer}},...,{\left(x_V^t\right)}_{g1}^{1\mathrm{stlayer}},...,{\left(x_V^t\right)}_{\mathrm{gm}}^{1\mathrm{stlayer}},...,{\left(x_V^t\right)}_{s1}^{1\mathrm{stlayer}},...,{\left(x_V^t\right)}_{s1}^{1\mathrm{stlayer}},...,]---\left(2\right)$

wherein,andrepresenting barrel vibration induced during the impact, rub and slip phases.

In practice, there are tens of thousands of steel balls in a mill. The steel balls are arranged in layers and fall simultaneously with different impact forces. These different frequencies and amplitudes of impact-induced vibrations are superimposed on each other. The mill's own mass imbalance and offset mounting of the mill can also cause the mill barrel to vibrate. These vibration signals couple to each other, and the resulting measurable cylinder vibration signal can be expressed as:

<math> <mrow> <mfenced open='' close=''> <mtable> <mtr> <mtd> <msubsup> <mi>x</mi> <mi>V</mi> <mi>t</mi> </msubsup> <mo>=</mo> <msubsup> <mrow> <mo>(</mo> <msubsup> <mi>x</mi> <mi>V</mi> <mi>t</mi> </msubsup> <mo>)</mo> </mrow> <mi>period</mi> <mrow> <mn>1</mn> <mi>stlayer</mi> </mrow> </msubsup> <mo>+</mo> <msubsup> <mrow> <mo>(</mo> <msubsup> <mi>x</mi> <mi>V</mi> <mi>t</mi> </msubsup> <mo>)</mo> </mrow> <mi>period</mi> <mrow> <mn>2</mn> <mi>ndlayer</mi> </mrow> </msubsup> <mo>+</mo> <msubsup> <mrow> <mo>(</mo> <msubsup> <mi>x</mi> <mi>V</mi> <mi>t</mi> </msubsup> <mo>)</mo> </mrow> <mi>period</mi> <mrow> <mn>3</mn> <mi>rdlayer</mi> </mrow> </msubsup> <mo>+</mo> <mo>.</mo> <mo>.</mo> <mo>.</mo> <mo>,</mo> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msubsup> <mi>x</mi> <mi>V</mi> <mi>t</mi> </msubsup> <mo>)</mo> </mrow> <mi>millself</mi> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msubsup> <mi>x</mi> <mi>V</mi> <mi>t</mi> </msubsup> <mo>)</mo> </mrow> <mi>install</mi> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msubsup> <mi>x</mi> <mi>V</mi> <mi>t</mi> </msubsup> <mo>)</mo> </mrow> <mi>others</mi> </msup> </mtd> </mtr> <mtr> <mtd> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <msub> <mi>j</mi> <mi>V</mi> </msub> <mo>=</mo> <mn>1</mn> </mrow> <msubsup> <mi>J</mi> <mi>V</mi> <mi>all</mi> </msubsup> </munderover> <msubsup> <mi>x</mi> <msub> <mi>j</mi> <mi>V</mi> </msub> <mi>t</mi> </msubsup> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow> </math>

wherein,andrespectively represent the j th_Vth component and the number of barrel vibrator signals; andrespectively representing the vibrator signals caused by the steel balls at the 1 st layer, the 2 nd layer and the 3 rd layer, the unbalance of the mass of the mill, the installation bias and other reasons.

The mill loading parameters include ball-to-ball ratio (MBVR), pulp concentration (PD) and filling rate (CVR), which are related to mill loading, mill loading conditions.

The acoustic radiation of the cylinder vibration, i.e. the structural noise, is the main component of the vibro-acoustic signal. Since the mill cylinder is a strong reflection surface in acoustics, the noise inside the mill is continuously reflected to form a mixed sound field, and the parts transmitted to the outside of the mill through the mill cylinder and the mill bolts are called air noise. The measured vibro-acoustic signals outside the mill grinding zone also contain noise adjacent the mill and other equipment. Thus, the composition of the vibro-acoustic signal can be represented by:

<math> <mrow> <msubsup> <mi>x</mi> <mi>A</mi> <mi>t</mi> </msubsup> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <msub> <mi>J</mi> <mi>A</mi> </msub> <mo>=</mo> <mn>1</mn> </mrow> <msubsup> <mi>J</mi> <mi>A</mi> <mi>all</mi> </msubsup> </munderover> <msubsup> <mi>x</mi> <msub> <mi>j</mi> <mi>A</mi> </msub> <mi>t</mi> </msubsup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow> </math>

wherein,is the jth_AThe sub-signals are then transmitted to the receiver,is the number of vibronic signals.

In conclusion, the vibration and vibration sound signals of the mill cylinder have multi-component and multi-scale characteristics, and the self-adaptive decomposition of the signals is necessary.

FIG. 2 is a flow chart of a method of soft measurement of mill load parameters in an embodiment of the present invention. As shown in FIG. 2, the soft measurement method of the load parameter of the mill comprises the following steps:

s100, detecting and acquiring a sample vibration signal and a sample vibration sound signal of a cylinder of the grinding machine operating under known load parameters.

As described above, the sample vibration signal and the sample local acoustic signal may be detected by operating the mill under a known load, and the corresponding load parameter may be converted according to the load.

S200, decomposing the sample vibration signal and the sample vibration sound signal into J based on the EEMD method_IMFSample sub-signals, each sample sub-signal representing a single vibrational mode having a physical meaning, J_IMFIs a predetermined value.

The EMD method is widely used for analyzing vibration signals of rotary Mechanical equipment (see the documents [26] Y.G.Lei, J.Lin, Z.J.He, M.J.Zuo, "A review on empirical mode composition in fault diagnosis of rotating machinery," Mechanical Systems and Signal Processing,35(2013) 108. 126.). The sub-signals (i.e., IMFs) obtained by the decomposition are arranged in order of frequency from top to bottom, and represent natural vibration modes in the original signal. Theoretically, each IMF represents a single vibrational mode with physical meaning, and a detailed mathematical description can be found in document [15], the flow of which is shown in fig. 3. EEMD can solve the modal mixing problem caused by EMD, and the relationship between the two is shown in FIG. 3, wherein M represents the number of sets.

Specifically, step S200 includes:

and S201, adding predetermined noise into the sample vibration signal and the sample vibration sound signal.

In the present embodiment, the predetermined noise added is white noise

S202, acquiring M groups of initial sub-signals based on an EMD method.

Specifically, as shown in fig. 3, the EMD method includes finding a signal extreme point, obtaining an upper envelope and a lower envelope of a signal, calculating a mean value of the upper envelope and the lower envelope, calculating a difference value between the envelope mean value and the signal and using the difference value as a new component, then checking whether the new component meets an IMF criterion, if not, returning the new component as an original signal to a first step of the EMD method, otherwise, using the new component as an IMF, further calculating a difference value between the new component and the original signal as a residual signal, and then further checking whether the residual signal is a monotonous signal, if so, ending a flow of the EMD method, otherwise, returning the residual signal to the first step of the EMD method instead of the original signal.

S203, calculating the mean value of M groups of initial sub-signals to obtain J_IMFA sub-signal.

By the EEMD method of step S200, the original sample vibration signal and vibration-sound signal can be decomposed into:

<math> <mrow> <msubsup> <mi>x</mi> <mi>V</mi> <mi>t</mi> </msubsup> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <msub> <mi>j</mi> <mi>V</mi> </msub> <mo>=</mo> <mn>1</mn> </mrow> <msubsup> <mi>J</mi> <mi>V</mi> <mi>all</mi> </msubsup> </munderover> <msub> <mrow> <mo>(</mo> <msubsup> <mover> <mi>c</mi> <mo>&OverBar;</mo> </mover> <mi>V</mi> <mi>t</mi> </msubsup> <mo>)</mo> </mrow> <msub> <mi>j</mi> <mi>V</mi> </msub> </msub> <mo>+</mo> <msub> <mrow> <mo>(</mo> <msub> <mover> <mi>r</mi> <mo>&OverBar;</mo> </mover> <mi>V</mi> </msub> <mo>)</mo> </mrow> <msubsup> <mi>J</mi> <mi>V</mi> <mi>all</mi> </msubsup> </msub> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow> </math>

<math> <mrow> <msubsup> <mi>x</mi> <mi>A</mi> <mi>t</mi> </msubsup> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <msub> <mi>j</mi> <mi>A</mi> </msub> <mo>=</mo> <mn>1</mn> </mrow> <msubsup> <mi>J</mi> <mi>A</mi> <mi>all</mi> </msubsup> </munderover> <msub> <mrow> <mo>(</mo> <msubsup> <mover> <mi>c</mi> <mo>&OverBar;</mo> </mover> <mi>A</mi> <mi>t</mi> </msubsup> <mo>)</mo> </mrow> <msub> <mi>j</mi> <mi>A</mi> </msub> </msub> <mo>+</mo> <msub> <mrow> <mo>(</mo> <msub> <mover> <mi>r</mi> <mo>&OverBar;</mo> </mover> <mi>A</mi> </msub> <mo>)</mo> </mrow> <msubsup> <mi>J</mi> <mi>A</mi> <mi>all</mi> </msubsup> </msub> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow> </math>

wherein,andis the jth_VAnd j_AThe number of the IMFs is one,andis the residual signal.

The maximum number of IMFs used for feature selection needs to be determined in advance. Since different mills have different characteristics, the number can be determined empirically from the result of the decomposition of the cylinder vibration signal at zero load. In this embodiment, the numbers of the IMFs of the cylinder vibration signal and the vibration acoustic signal are respectively denoted by J_AAnd J_VAnd all IMFs are relabeled as:

<math> <mrow> <msubsup> <mrow> <mo>{</mo> <msubsup> <mover> <mi>c</mi> <mo>&OverBar;</mo> </mover> <msub> <mi>j</mi> <mi>IMF</mi> </msub> <mi>f</mi> </msubsup> <mo>}</mo> </mrow> <mrow> <msub> <mi>j</mi> <mi>IMF</mi> </msub> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>J</mi> <mi>IMF</mi> </msub> </msubsup> <mo>=</mo> <mo>[</mo> <msub> <mrow> <mo>(</mo> <msubsup> <mover> <mi>c</mi> <mo>&OverBar;</mo> </mover> <mi>V</mi> <mi>f</mi> </msubsup> <mo>)</mo> </mrow> <mn>1</mn> </msub> <mo>,</mo> <mo>.</mo> <mo>.</mo> <mo>.</mo> <mo>,</mo> <msub> <mrow> <mo>(</mo> <msubsup> <mover> <mi>c</mi> <mo>&OverBar;</mo> </mover> <mi>V</mi> <mi>f</mi> </msubsup> <mo>)</mo> </mrow> <msub> <mi>j</mi> <mi>V</mi> </msub> </msub> <mo>,</mo> <mo>.</mo> <mo>.</mo> <mo>.</mo> <mo>,</mo> <msub> <mrow> <mo>(</mo> <msubsup> <mover> <mi>c</mi> <mo>&OverBar;</mo> </mover> <mi>V</mi> <mi>f</mi> </msubsup> <mo>)</mo> </mrow> <msub> <mi>J</mi> <mi>V</mi> </msub> </msub> <mo>,</mo> <msub> <mrow> <mo>(</mo> <msubsup> <mover> <mi>c</mi> <mo>&OverBar;</mo> </mover> <mi>A</mi> <mi>f</mi> </msubsup> <mo>)</mo> </mrow> <mn>1</mn> </msub> <mo>,</mo> <mo>.</mo> <mo>.</mo> <mo>.</mo> <mo>,</mo> <msub> <mrow> <mo>(</mo> <msubsup> <mover> <mi>c</mi> <mo>&OverBar;</mo> </mover> <mi>A</mi> <mi>f</mi> </msubsup> <mo>)</mo> </mrow> <msub> <mi>j</mi> <mi>A</mi> </msub> </msub> <mo>,</mo> <mo>.</mo> <mo>.</mo> <mo>.</mo> <mo>,</mo> <msub> <mrow> <mo>(</mo> <msubsup> <mover> <mi>c</mi> <mo>&OverBar;</mo> </mover> <mi>A</mi> <mi>f</mi> </msubsup> <mo>)</mo> </mrow> <msub> <mi>J</mi> <mi>A</mi> </msub> </msub> <mo>]</mo> <mo>=</mo> <mo>[</mo> <msubsup> <mover> <mi>c</mi> <mo>&OverBar;</mo> </mover> <mn>1</mn> <mi>f</mi> </msubsup> <mo>,</mo> <mo>.</mo> <mo>.</mo> <mo>,</mo> <msubsup> <mover> <mi>c</mi> <mo>&OverBar;</mo> </mover> <msub> <mi>j</mi> <mi>IMF</mi> </msub> <mi>f</mi> </msubsup> <mo>,</mo> <mo>.</mo> <mo>.</mo> <mo>.</mo> <mo>,</mo> <msubsup> <mover> <mi>c</mi> <mo>&OverBar;</mo> </mover> <msub> <mi>J</mi> <mi>IMF</mi> </msub> <mi>f</mi> </msubsup> <mo>]</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow> </math>

wherein, J_IMF＝J_A+J_VIs the number of total IMFs, due to J as described above_AAnd J_VAre all predetermined values, therefore J_IMFAlso a predetermined value.

S300, extracting the spectral features of all the sub-signals, and acquiring a training sample set and a verification sample set based on the spectral features and corresponding load parameters, wherein the spectral features are calculated according to a Marginal Spectrum (MSHT) based on Hilbert Transform (HT) of the sub-signals, a mean value and a variance (MVHT) of instantaneous amplitude and frequency of the Hilbert transform, and a Power Spectral Density (PSD) based on fast Fourier transform.

In the present invention, three types of features are extracted from the IMF of the vibration signal and the vibration-sound signal of the cylinder: a hilbert transform based marginal spectrum (abbreviated MSHT), a hilbert transform instantaneous magnitude and frequency mean and variance (abbreviated MVHT), and an FFT based power spectral density (abbreviated PSD). These three types of features have their own characteristics and have been successfully used in various fields (see documents [26], [22], [17 ]). Thus, these spectral features can serve as multi-source information from different perspectives. The present embodiment performs selection and combination based on the above extracted features to calculate spectral features for extracting all sub-signals.

Specifically, as shown in fig. 4, step S300 includes:

and S301, calculating MSHT, MVHT and PSD of the subsignals.

The same type of features of different IMFs are combined and denoted in this embodiment $Z_{\mathrm{MVFT}}^{\mathrm{ori}}={\left\{z_{\mathrm{MVFT}}^{\mathrm{ori}}\right\}}_{l=1}^k$ And $Z_{\mathrm{PSD}}^{\mathrm{ori}}={\left\{z_{\mathrm{PSD}}^{\mathrm{ori}}\right\}}_{l=1}^k.$

and S302, calculating mutual information values between all MSHT spectrum variables and mill load parameters of the MSHT respectively.

For MSHT, a density estimation method can be used (see document [27 ]]H.C.Peng, F.H.Long, C.Ding, Feature selection based on statistical information criterion of max-dependency, and min-redundancy. IEEE Transactions on Pattern Analysis and machinery Analysis, 27(2005) 1226-1238) calculation of the p-th spectral variableAnd the mill load parameter y:

<math> <mrow> <mfenced open='' close=''> <mtable> <mtr> <mtd> <mi>Mi</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>;</mo> <msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mi>MSHT</mi> <mi>ori</mi> </msubsup> <mo>)</mo> </mrow> <mi>p</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mo>&Integral;</mo> <mo>&Integral;</mo> <mi>p</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>,</mo> <msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mi>MSHT</mi> <mi>ori</mi> </msubsup> <mo>)</mo> </mrow> <mi>p</mi> </msub> <mo>)</mo> </mrow> <mi>log</mi> <mfrac> <mrow> <mi>p</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>,</mo> <msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mi>MSHT</mi> <mi>ori</mi> </msubsup> <mo>)</mo> </mrow> <mi>p</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <mi>p</mi> <mrow> <mo>(</mo> <msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mi>MSHT</mi> <mi>ori</mi> </msubsup> <mo>)</mo> </mrow> <mi>p</mi> </msub> <mo>)</mo> </mrow> <mi>p</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mi>d</mi> <mrow> <mo>(</mo> <msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mi>MSHT</mi> <mi>ori</mi> </msubsup> <mo>)</mo> </mrow> <mi>p</mi> </msub> <mo>)</mo> </mrow> <mi>dy</mi> </mtd> </mtr> <mtr> <mtd> <mo>=</mo> <mi>H</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <mi>H</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>|</mo> <msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mi>MSHT</mi> <mi>ori</mi> </msubsup> <mo>)</mo> </mrow> <mi>p</mi> </msub> <mo>)</mo> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow> </math>

wherein,and p (y) isAnd the edge probability density of y,is the combined probability density of the two or more,is the condition entropy of the sample,it is the information moisture.

S303, obtaining the maximum value of mutual information values of MSHT spectrum variablesAnd minimum value

S304, obtaining a plurality of spectrum feature selection threshold values which change between the minimum value and the maximum value of mutual information values of MSHT spectrum variables in a preset step length, selecting all the MSHT spectrum variables of which the corresponding mutual information values are larger than the spectrum feature selection threshold value based on each spectrum feature selection threshold value to construct a PLS model, and calculating the prediction error of the PLS model corresponding to each spectrum feature selection threshold value.

Specifically, the spectral features are selected based on the following formula:

<math> <mrow> <msub> <mi>&zeta;</mi> <msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mi>MSFT</mi> <mi>ori</mi> </msubsup> <mo>)</mo> </mrow> <mi>p</mi> </msub> </msub> <mo>=</mo> <mfenced open='{' close=''> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <mi>if</mi> </mtd> <mtd> <mi>Muin</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>;</mo> <msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mi>MSHT</mi> <mi>ori</mi> </msubsup> <mo>)</mo> </mrow> <mi>p</mi> </msub> <mo>)</mo> </mrow> <mo>&GreaterEqual;</mo> <msub> <mi>&theta;</mi> <mi>Muin</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mi>else</mi> </mtd> <mtd> <mi>Muin</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>;</mo> <msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mi>MSHT</mi> <mi>ori</mi> </msubsup> <mo>)</mo> </mrow> <mi>p</mi> </msub> <mo>)</mo> </mrow> <mo>&lt;</mo> <msub> <mi>&theta;</mi> <mi>Muin</mi> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow> </math>

wherein,θ_Muinis the threshold for the selection of the spectral features,andis thatMinimum and maximum values of. The step size theta between threshold increments is calculated using the formula_step：

<math> <mrow> <msub> <mi>&theta;</mi> <mi>step</mi> </msub> <mo>=</mo> <mrow> <mo>(</mo> <msubsup> <mi>&theta;</mi> <mi>Muin</mi> <mi>Max</mi> </msubsup> <mo>-</mo> <msubsup> <mi>&theta;</mi> <mi>Muin</mi> <mi>MIn</mi> </msubsup> <mo>)</mo> </mrow> <mo>/</mo> <mn>10</mn> <mo>.</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo> </mrow> </mrow> </math>

SelectingTo construct a partial least squares regression (PLS) model. Then by step size theta_stepThe threshold for MI is increased and the process is repeated. Finally, the threshold with the smallest PLS model prediction error is selected as the final threshold, and the selected spectral features are scored as

In the above process, the IMFs used to select spectral features do not need to be selected in advance before the feature selection process is performed, and these valuable IMFs and their features are adaptively and selectively determined.

S305, the MSHT spectrum variable set corresponding to the PLS model with the minimum prediction error is combined into MSHT spectrum characteristics.

In parallel with the above steps S302-305, steps S306-309 and steps S310-313 are respectively alignedAndand respectively obtaining the MVHT spectral characteristics corresponding to the MVHT and the PSD spectral characteristics corresponding to the PSD by adopting the same process.

S306, calculating mutual information values between all MVHT spectrum variables of the MVHT and the mill load parameters respectively;

s307, acquiring the maximum value and the minimum value of mutual information values of the MVHT spectrum variables.

S308, obtaining a plurality of spectrum feature selection thresholds which change between the minimum value and the maximum value of mutual information values of MVHT spectrum variables by a preset step length, selecting all the MVHT spectrum variables of which the corresponding mutual information values are larger than the spectrum feature selection thresholds based on each spectrum feature selection threshold to construct a PLS model, and calculating the PLS model prediction error corresponding to each spectrum feature selection threshold.

S309, the MVHT spectrum variable set corresponding to the PLS model with the minimum prediction error is combined into the MVHT spectrum characteristic.

And S310, calculating mutual information values between all PSD spectrum variables of the PSD and the mill load parameters respectively.

And S311, acquiring the maximum value and the minimum value of mutual information values of the PSD spectrum variables.

S312, a plurality of spectrum feature selection threshold values which change between the minimum value and the maximum value of mutual information values of PSD spectrum variables in a preset step length are obtained, all PSD spectrum variables of which the corresponding mutual information values are larger than the spectrum feature selection threshold values are selected to construct a PLS model based on each spectrum feature selection threshold value, and a PLS model prediction error corresponding to each spectrum feature selection threshold value is calculated.

S313, the PSD spectrum variable set corresponding to the PLS model with the minimum prediction error is combined into PSD spectrum characteristics.

And S314, combining the MSHT spectral characteristics, the MVHT spectral characteristics and the PSD spectral characteristics to obtain the spectral characteristics of all the sub-signals.

The final spectral features of all sub-signals are labeled:

$Z=[Z_{\mathrm{MSHT}},Z_{\mathrm{MVFT}},Z_{\mathrm{PSD}}]={\left\{z\right\}}_{l=1}^k,---\left(11\right)$

wherein Z is_MSHT,Z_MVFTAnd Z_PSDRespectively, features selected from MSHT, MVHT and PSD.

Thus, samples for training and validation can be obtained. The subsequent step is based on the sample to carry out 3 parts of operations, namely candidate submodel construction, integrated submodel selection and integrated submodel combination to obtain a model for soft measurement.

S400, generating J training subsamples from the training sample set based on Bootstrap algorithm.

The first problem of SEN modeling is integrated construction. Method for sampling training sample based on Bootstrap algorithmMiddle schoolRaw training child sampleWhere J is the number of subsamples, i.e. the number of candidate submodels and the population number in the GA algorithm.

S500, training and obtaining J candidate sub-models from J training sub-samples by adopting the same number of nuclear latent variables and nuclear parameters based on a Kernel Partial Least Squares (KPLS) algorithm.

Based on KPLS algorithm, the sub-training samples are adopted to construct candidate sub-models, and the jth sub-training sample is usedMapping to a high-dimensional feature space:

K^j＝Φ((z^j)_l)^TΦ((z^j)_m),l,m＝1,2,…k. (12)

wherein, K^jCentralization was performed using the formula:

${\tilde{K}}^j=(I-\frac{1}{k}{1}_k{1}_k^T)K^j(I-\frac{1}{k}{1}_k{1}_k^T),---\left(13\right)$

wherein I is a k-dimensional unit matrix, 1_kIs a vector of length k with a value of 1.

Selecting the same modeling parameter, i.e. kernel parameter K, for all candidate submodels_paraAnd nuclear latent variable (KLV) h_KLVAnd marking the candidate submodels as

S600, selecting all candidate submodels with corresponding model selection weight parameters larger than a model selection threshold value from the J candidate submodels according to the verification sample set to form an integrated submodel set, wherein the optimized weight parameters are obtained by optimizing the randomly generated initial weight parameters through a genetic algorithm with the minimized prediction error as a target.

Validating a sampleThe predicted output based on the jth candidate submodel is:

wherein k is^validIs the number of validation samples.

The prediction error is calculated by the following formula:

$e_{\mathrm{valid}}^j={\hat{Y}}_{\mathrm{valid}}^j-Y_{\mathrm{valid}}.---\left(15\right)$

wherein, Y_validThe actual load parameter in the sample can be obtained by load conversion corresponding to experiments.

The correlation coefficient of the jth and the s candidate submodels is calculated using the following formula:

<math> <mrow> <msubsup> <mi>c</mi> <mi>js</mi> <mi>valid</mi> </msubsup> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <msup> <mi>k</mi> <mi>valid</mi> </msup> </munderover> <msubsup> <mi>e</mi> <mi>valid</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>j</mi> <mo>,</mo> <msup> <mi>k</mi> <mi>valid</mi> </msup> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <msubsup> <mi>e</mi> <mi>valid</mi> <mi>s</mi> </msubsup> <mrow> <mo>(</mo> <mi>s</mi> <mo>,</mo> <msup> <mi>k</mi> <mi>valid</mi> </msup> <mo>)</mo> </mrow> <mo>/</mo> <msup> <mi>k</mi> <mi>valid</mi> </msup> <mo>.</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>16</mn> <mo>)</mo> </mrow> </mrow> </math>

thus, the correlation coefficient matrix can be written as:

<math> <mrow> <msup> <mi>C</mi> <mi>valid</mi> </msup> <mo>=</mo> <msub> <mfenced open='[' close=']'> <mtable> <mtr> <mtd> <msubsup> <mi>c</mi> <mn>11</mn> <mi>valid</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>c</mi> <mn>12</mn> <mi>valid</mi> </msubsup> </mtd> <mtd> <mo>.</mo> <mo>.</mo> <mo>.</mo> </mtd> <mtd> <msubsup> <mi>c</mi> <mrow> <mn>1</mn> <mi>J</mi> </mrow> <mi>valid</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>c</mi> <mn>21</mn> <mi>valid</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>c</mi> <mn>22</mn> <mi>valid</mi> </msubsup> </mtd> <mtd> <mo>.</mo> <mo>.</mo> <mo>.</mo> </mtd> <mtd> <msubsup> <mi>c</mi> <mrow> <mn>2</mn> <mi>J</mi> </mrow> <mi>valid</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> <mo>.</mo> <mo>.</mo> </mtd> <mtd> <mo>.</mo> <mo>.</mo> <mo>.</mo> </mtd> <mtd> <msubsup> <mi>c</mi> <mi>js</mi> <mi>valid</mi> </msubsup> </mtd> <mtd> <mo>.</mo> <mo>.</mo> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>c</mi> <mrow> <mi>J</mi> <mn>1</mn> </mrow> <mi>valid</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>c</mi> <mrow> <mi>J</mi> <mn>2</mn> </mrow> <mi>valid</mi> </msubsup> </mtd> <mtd> <mo>.</mo> <mo>.</mo> <mo>.</mo> </mtd> <mtd> <msubsup> <mi>c</mi> <mi>JJ</mi> <mi>valid</mi> </msubsup> </mtd> </mtr> </mtable> </mfenced> <mrow> <mi>J</mi> <mo>&times;</mo> <mi>J</mi> </mrow> </msub> <mo>.</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>17</mn> <mo>)</mo> </mrow> </mrow> </math>

generating a set of random vectors for candidate sub-modelsAs initial weight parameter, based on C^validThese random vectors were optimized with the goal of minimizing prediction error using standard genetic algorithms and the results were scored asAs an optimization weight parameter.

According to the optimized weight parameters, the integration submodel can be selected by adopting the following criteria:

<math> <mrow> <msub> <mi>&xi;</mi> <mi>j</mi> </msub> <mo>=</mo> <mfenced open='{' close=''> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <mi>if</mi> </mtd> <mtd> <msubsup> <mi>w</mi> <mi>j</mi> <mo>*</mo> </msubsup> <mo>&GreaterEqual;</mo> <mi>&lambda;</mi> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mi>else</mi> </mtd> <mtd> <msubsup> <mi>w</mi> <mi>j</mi> <mo>*</mo> </msubsup> <mo>&lt;</mo> <mi>&lambda;</mi> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>18</mn> <mo>)</mo> </mrow> </mrow> </math>

where λ is a selection threshold for the submodel, which is a predetermined value. Xi is reduced_jThe submodel 1 is selected as the integrated submodel, and the number thereof is denoted as J^*I.e. the integration size of the SEN model.

Will j to^*Individual integrated submodel taggingThe output is:

all integration submodels are noted

That is, step S600 includes:

s601, calculating prediction output based on the jth candidate sub-model for each verification sample.

S602, calculating a prediction error based on the prediction output and the actual load parameter

S603, calculating correlation coefficients among all candidate submodels based on the following formula to form a correlation coefficient matrix;

<math> <mrow> <msubsup> <mi>c</mi> <mi>js</mi> <mi>valid</mi> </msubsup> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <msup> <mi>k</mi> <mi>valid</mi> </msup> </munderover> <msubsup> <mi>e</mi> <mi>valid</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>j</mi> <mo>,</mo> <msup> <mi>k</mi> <mi>valid</mi> </msup> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <msubsup> <mi>e</mi> <mi>valid</mi> <mi>s</mi> </msubsup> <mrow> <mo>(</mo> <mi>s</mi> <mo>,</mo> <msup> <mi>k</mi> <mi>valid</mi> </msup> <mo>)</mo> </mrow> <mo>/</mo> <msup> <mi>k</mi> <mi>valid</mi> </msup> </mrow> </math>

wherein,is the correlation coefficient, k, between the jth candidate submodel and the S candidate submodel^validTo verify the number of samples, the correlation coefficient matrix is:

<math> <mrow> <msup> <mi>C</mi> <mi>valid</mi> </msup> <mo>=</mo> <msub> <mfenced open='[' close=']'> <mtable> <mtr> <mtd> <msubsup> <mi>c</mi> <mn>11</mn> <mi>valid</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>c</mi> <mn>12</mn> <mi>valid</mi> </msubsup> </mtd> <mtd> <mo>.</mo> <mo>.</mo> <mo>.</mo> </mtd> <mtd> <msubsup> <mi>c</mi> <mrow> <mn>1</mn> <mi>J</mi> </mrow> <mi>valid</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>c</mi> <mn>21</mn> <mi>valid</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>c</mi> <mn>22</mn> <mi>valid</mi> </msubsup> </mtd> <mtd> <mo>.</mo> <mo>.</mo> <mo>.</mo> </mtd> <mtd> <msubsup> <mi>c</mi> <mrow> <mn>2</mn> <mi>J</mi> </mrow> <mi>valid</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> <mo>.</mo> <mo>.</mo> </mtd> <mtd> <mo>.</mo> <mo>.</mo> <mo>.</mo> </mtd> <mtd> <msubsup> <mi>c</mi> <mi>js</mi> <mi>valid</mi> </msubsup> </mtd> <mtd> <mo>.</mo> <mo>.</mo> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>c</mi> <mrow> <mi>J</mi> <mn>1</mn> </mrow> <mi>valid</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>c</mi> <mrow> <mi>J</mi> <mn>2</mn> </mrow> <mi>valid</mi> </msubsup> </mtd> <mtd> <mo>.</mo> <mo>.</mo> <mo>.</mo> </mtd> <mtd> <msubsup> <mi>c</mi> <mi>JJ</mi> <mi>valid</mi> </msubsup> </mtd> </mtr> </mtable> </mfenced> <mrow> <mi>J</mi> <mo>&times;</mo> <mi>J</mi> </mrow> </msub> </mrow> </math>

s604, generating a random vector for the candidate submodel

S605, using random vectorCalculating an optimized weight parameter which minimizes a prediction error by a genetic algorithm according to the correlation coefficient matrix as an initial weight parameter;

s606, selecting the candidate submodels with the corresponding optimization weight parameters larger than or equal to the preset selection threshold value as the integration submodels to obtain an integration submodel set.

S700, calculating and acquiring integration weights corresponding to all the integration sub-model sets according to an AWF algorithm.

The different integration submodels contribute differently to the SEN model. The embodiment of the invention adopts an Adaptive Weighted Fusion (AWF) algorithm to calculate the optimal weight of the integrated submodel, namely solving the following optimization problem:

<math> <mrow> <mfenced open='' close=''> <mtable> <mtr> <mtd> <mi>min</mi> </mtd> <mtd> <msup> <mi>&sigma;</mi> <mn>2</mn> </msup> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <msup> <mi>j</mi> <mo>*</mo> </msup> <mo>=</mo> <mn>1</mn> </mrow> <msup> <mi>J</mi> <mo>*</mo> </msup> </munderover> <msup> <mrow> <mo>(</mo> <msubsup> <mi>W</mi> <msup> <mi>j</mi> <mo>*</mo> </msup> <mi>AWF</mi> </msubsup> <mo>)</mo> </mrow> <mn>2</mn> </msup> <msubsup> <mi>&sigma;</mi> <msup> <mi>j</mi> <mo>*</mo> </msup> <mn>2</mn> </msubsup> </mtd> </mtr> <mtr> <mtd> <mi>s</mi> <mo>.</mo> <mi>t</mi> <mo>.</mo> </mtd> <mtd> <munderover> <mi>&Sigma;</mi> <mrow> <msup> <mi>j</mi> <mo>*</mo> </msup> <mo>=</mo> <mn>1</mn> </mrow> <msup> <mi>J</mi> <mo>*</mo> </msup> </munderover> <msubsup> <mi>W</mi> <msup> <mi>j</mi> <mo>*</mo> </msup> <mi>AWF</mi> </msubsup> <mo>=</mo> <mn>1,0</mn> <mo>&le;</mo> <msubsup> <mi>W</mi> <msup> <mi>j</mi> <mo>*</mo> </msup> <mi>AWF</mi> </msubsup> <mo>&le;</mo> <mn>1</mn> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>20</mn> <mo>)</mo> </mrow> </mrow> </math>

where σ is the combined predictorThe variance of (a) is determined,is the integrated submodel prediction isThe variance of (c). The optimal weight is calculated using the following formula:

<math> <mrow> <msubsup> <mi>W</mi> <msup> <mi>j</mi> <mo>*</mo> </msup> <mi>AWF</mi> </msubsup> <mo>=</mo> <mn>1</mn> <mo>/</mo> <mrow> <mo>(</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>&sigma;</mi> <msup> <mi>j</mi> <mo>*</mo> </msup> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <munderover> <mi>&Sigma;</mi> <mrow> <msup> <mi>j</mi> <mo>*</mo> </msup> <mo>=</mo> <mn>1</mn> </mrow> <msup> <mi>J</mi> <mo>*</mo> </msup> </munderover> <mfrac> <mn>1</mn> <msup> <mrow> <mo>(</mo> <msub> <mi>&sigma;</mi> <msup> <mi>j</mi> <mo>*</mo> </msup> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> <mo>)</mo> </mrow> <mo>.</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>21</mn> <mo>)</mo> </mrow> </mrow> </math>

minimum variance of prediction values of SEN modelThe method comprises the following steps:

<math> <mrow> <msubsup> <mi>&sigma;</mi> <mi>min</mi> <mn>2</mn> </msubsup> <mo>=</mo> <mn>1</mn> <mo>/</mo> <munderover> <mi>&Sigma;</mi> <mrow> <msup> <mi>j</mi> <mo>*</mo> </msup> <mo>=</mo> <mn>1</mn> </mrow> <msup> <mi>J</mi> <mo>*</mo> </msup> </munderover> <mfrac> <mn>1</mn> <msubsup> <mi>&sigma;</mi> <msup> <mi>j</mi> <mo>*</mo> </msup> <mn>2</mn> </msubsup> </mfrac> <mo>.</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>22</mn> <mo>)</mo> </mrow> </mrow> </math>

therefore, the establishment of a prediction model for soft measurement is successful, and soft measurement operation can be carried out based on the extracted parameters of the mill to be measured.

And S800, acquiring the spectral characteristics of the test data of the mill needing soft measurement.

The method for acquiring the spectral feature of the test data in step S800 may be performed in a manner referring to S100-S300, and will not be described herein again.

Specifically, S800 includes:

s801, detecting vibration signals to be measured and vibration sound signals to be measured of a cylinder body of the mill to be measured.

S802, decomposing the vibration signal to be measured and the vibration sound signal to be measured into JIMF sample sub-signals based on the EEMD method, wherein each sub-signal represents a single vibration mode with physical meaning, and the JIMF is a preset value.

And S803, extracting the spectral characteristics of all the sub-signals, wherein the spectral characteristics are obtained by calculating the Marginal Spectrum (MSHT) based on the Hilbert transform, the mean value and the variance (MVHT) of the instantaneous amplitude and frequency of the Hilbert transform and the Power Spectral Density (PSD) based on the fast Fourier transform of the sub-signals.

And S900, calculating load parameters corresponding to the spectrum characteristics of the test data according to the integrated sub-model set and the corresponding integrated weights.

The predicted output of the test sample based on the SEN model is:

<math> <mrow> <msub> <mover> <mi>y</mi> <mo>^</mo> </mover> <mi>test</mi> </msub> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <msup> <mi>j</mi> <mo>*</mo> </msup> <mo>=</mo> <mn>1</mn> </mrow> <msup> <mi>J</mi> <mo>*</mo> </msup> </munderover> <msubsup> <mi>W</mi> <msup> <mi>j</mi> <mo>*</mo> </msup> <mi>AWF</mi> </msubsup> <mo>&CenterDot;</mo> <msubsup> <mover> <mi>y</mi> <mo>^</mo> </mover> <mi>test</mi> <msup> <mi>j</mi> <mo>*</mo> </msup> </msubsup> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>23</mn> <mo>)</mo> </mrow> </mrow> </math>

wherein,is based on the j^*The predicted output of the individual integration submodels is calculated as follows:

${\hat{y}}_{\mathrm{test}}^{j^{*}}==f^{j*}\left({\left\{{\left(z_{\mathrm{test}}\right)}_l\right\}}_{l=1}^{k^{\mathrm{test}}}\right),---\left(24\right)$

wherein,representing the spectral characteristics of the test data.

The actual ball mill was tested based on the soft testing method of the present invention and compared with other methods. The experiment was carried out on an XMQL 420X 450 lattice ball mill with roller outer diameter and length of 460 mm. The mill is driven by a three-phase motor with the power of 2.12kw, the maximum steel ball loading capacity is 80kg, the designed milling capacity is 10 kg/h, and the rotating speed is 57 revolutions per minute. The middle part of the mill is opened for adding steel balls, materials and water load. The materials adopted in the experiment are copper ores, the diameters of the materials are less than 6mm, and the density is 4.2t/m < 3 >. Steel balls with the diameters of 30 mm, 20 mm and 15mm are used as grinding media, and the ratio is 3:4: 3. The load in the mill was a total of 4 groups of experiments including ball load, material load, and water load: first group: the material load is 10kg, and the water load is increased from 5kg to 40 kg; the second group: the material load is 20kg, and the water load is increased from 2kg to 20 kg; third group: the water load is 2kg, and the material load is increased from 10kg to 20 kg; fourth group: the water load was 10kg and the material load increased from 22kg to 50 kg.

The IMFs for the barrel vibration and vibro-acoustic signals are labeled VIMF and AIMF. Selection A_noiseThe vibration signal and the vibration-sound signal of 4 periods of mill rotation are decomposed into 17 IMFs, respectively, 0.1 and M10. The four different grinding conditions were idle (zero load of mill), impact (only 20kg ball load in mill), dry (40 kg ball load in mill, 10kg material load), wet (40 kg ball load in mill, 10kg material load, 5kg water load).

The relationship between MI threshold and Root Mean Square Relative Error (RMSREs) when different mill loading parameters were selected is shown in FIG. 5. As can be seen from fig. 5, the optimal MI thresholds for the MSHT, MVHT and PSD characteristics selected for MBVR, PD and CVR are 0.7, 0.7 and 0.6, 0.2, 0.7 and 1.0, 0.5, 0.6 and 0.3, respectively.

FIG. 6a is a plot of maximum and minimum MI values and feature numbers of MSHT spectral features selected for MBVR; FIG. 6b is a plot of maximum and minimum MI values and the number of features of the selected MVHT spectral features for MBVR; FIG. 6c is a plot of the maximum and minimum MI values of the PSD spectral features selected for MBVR versus the number of features. FIGS. 6a-6c show: MSHT from AIMF1 for MBVR maximum MI values; the MI value of the PSD characteristic of AIMF1 was slightly lower than MSHT; the MI threshold of MVHT is less than 1 and is about 50 lower than the MI values of MSHT and PSD. These conclusions are similar to previous studies. In addition, more of the MSHT features stem from VIMF, the main reason being that vibration is the main source of vibro-acoustic signals.

FIG. 7a is a graph of maximum and minimum MI values and feature numbers of MSHT spectral features selected for PD; FIG. 7b is a graph of maximum and minimum MI values and the number of features of the selected MVHT spectral features for a PD; FIG. 7c is a graph of the maximum and minimum MI values of the PSD spectral features selected for the PD as well as the number of features. FIGS. 7a-7c show: for PD, the maximum MI value is derived from the MSHT signature of VIMF 1; the maximum MI value of the PSD features is derived from VIMF 8; the number of features selected from VIMF was higher than that of AIMF, consistent with previous studies.

FIG. 8a is a plot of maximum and minimum MI values and feature numbers of MSHT spectral features selected for a CVR; FIG. 8b is a plot of maximum and minimum MI values and the number of features of the selected MVHT spectral features for CVR; FIG. 8c is a plot of the maximum and minimum MI values of the PSD spectral features selected for CVR, as well as the number of features. FIGS. 8a-8c show that for CVR, the maximum MI values are derived from the MSHT characteristic of AIMF 2; the maximum MI value of PSD is from AIMF 1; the number of features selected from the AIMF is much smaller than the VIMF, possibly because the drum vibration signal is used at a much higher frequency than the acoustic signal. Therefore, more additional experiments are required for validation in future studies.

Thus, most of the spectral features selected were from the first 6 VIMFs and AIMFs, a conclusion that is not consistent with the MVHT selected for MBVR and PD; in addition, the IMF numbers of the selected features and the selected feature order are not continuous, indicating that different spectral features have different mappings with different mill loading parameters, and also indicating that these valuable IMFs and features are adaptively selected. Therefore, the multi-scale spectrum feature adaptive selection method provided by the invention is effective.

The statistical results of the multi-scale spectral feature selection are shown in table 1.

TABLE 1 multiscale spectral feature selection statistics

These selected features are serially combined to construct a soft measurement model based on IGASEN-KPLS. The learning parameters selected for the MBVR, PD and CVR models were: the same Radial Basis Function (RBF) was chosen for the kernel type, the number of KLVs 3, 8 and 7, respectively, the kernel radius 30, 8 and 20, respectively, the population number 20, and the sub-model selection threshold 0.05.

To overcome the random initialization property of the Genetic Algorithm (GA), the soft measurement model was run repeatedly 40 times. The curves of the integrated dimensions and prediction error are shown in fig. 9 and 10, and the statistical results of the predicted performance are shown in table 2.

TABLE 2 statistical results of soft measurement model based on IGASEN-KPLS

Figure 9 shows that the median values of the integrated sizes of MBVR, PD and CVR models are 3, 4 and 3, respectively. Fig. 10 shows that the mean values of these model RMSREs are 0.1114, 0.2137, and 0.1983, respectively. Among the three soft measurement models, the PD model has the largest integration size and the largest prediction error, and the stability is also the weakest.

Tables 1 and 2 show that:

(1) the MBVR model has the minimum prediction error and the minimum input feature quantity, and the method can effectively model the MBVR.

(2) The PD and CVR models have prediction accuracy of region MBVR; and the number of selected spectral features is much higher than that of the MBVR model, indicating that the number of features associated with these process parameters is high, and also indicating that a secondary selection of features is necessary.

The EEMD and IGASEN-KPLS based multi-scale SEN model proposed by the present invention and other methods in previous studies, such as FFT based single-scale single model [8]Single-scale SEN model based on FFT [6]]EMD and FFT-based multiscale SEN model [22]A comparison was made. In table 3, the selection of different integration numbers M and noise amplitudes a is given herein_noiseEEMD decomposition is carried out to establish the statistical result of the SEN model.

TABLE 3 statistical results of different modeling methods

Table 3 shows that:

(1) for MBVR models, the methods presented herein have the best prediction accuracy in all single-scale and multi-scale models. Document [6] adopts a SEN modeling method based on single-scale spectral feature branch-bound and AWF; the document [22] adopts MI, branch and bound and SEN modeling methods based on EMD multi-scale spectral features, and fails to realize the self-adaptive selection of IMF and features thereof. Neither of these approaches effectively utilizes valuable training samples. The methods presented herein are effective for measuring MBVR.

(2) The prediction accuracy of the method provided by the invention is lower than that of the literature [6] aiming at PD and CVR models. The multi-scale nature of barrel vibration and rattle is not fully reflected in PD and CVR, since only four mill rotation cycles of data are used. The FFT-based method [6] only utilizes single-scale spectral features, and the grinding mechanism of the mill is difficult to effectively show. There is a need to develop more efficient modeling methods.

(3) The average prediction accuracy of the EEMD-based multi-scale SEN model is 2 times that of the EMDf-based method, one reason is that the EEMD can more effectively obtain the EMD, the IMF and the characteristics thereof can be automatically selected, and the multi-scale spectrum characteristics and the multi-working-condition training samples are effectively and selectively fused.

(4) The EEMD algorithm resolves the change in the parameter, i.e. the noise amplitude A_noise(when M is 10, A_noiseFrom 0.05 to 0.5) and the number of sets M (at A)_noiseWhen M is 0.1, M is from 10 to 100), the prediction accuracy of the soft measurement changes less. One reason may be that the background interference of the experimental ball mill is relatively small and the next study will be validated on an industrial mill.

The invention realizes soft measurement of the load parameter of the mill by three steps. Firstly, decomposing a vibration signal and a vibration sound signal of a mill cylinder into Intrinsic Mode Functions (IMFs) with different time scales and physical meanings by adopting an Ensemble Empirical Mode Decomposition (EEMD) technology; then, three types of characteristics (frequency spectrum, marginal spectrum, Hilbert transformation instantaneous amplitude and mean and variance of frequency) of the multi-scale IMF are selected by adopting a self-adaptive characteristic selection method based on Mutual Information (MI); finally, a soft measurement model based on a selective set Kernel Partial Least Squares (KPLS) method is constructed based on the selected spectral features and the training samples. The simulation experiment result based on the small ball mill shows that the method can effectively detect the load parameters.

The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention, and various modifications and changes may be made by those skilled in the art. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims (8)

1. A soft measurement method for load parameters of a grinding machine comprises the following steps:

s100, detecting and acquiring a sample vibration signal and a sample vibration sound signal of a cylinder of the grinding machine operating under known load parameters;

s200, decomposing the sample vibration signal and the sample vibration sound signal into J based on an Ensemble Empirical Mode Decomposition (EEMD) method_IMFSample sub-signals, each sample sub-signal representing a single vibrational mode having a physical meaning, J_IMFIs a predetermined value;

s300, extracting spectral features of all the sub-signals, and acquiring a training sample set and a verification sample set based on the spectral features and corresponding load parameters, wherein the spectral features are acquired from a Marginal Spectrum (MSHT) based on Hilbert Transform (HT) of the sub-signals, a mean value and a variance (MVHT) of instantaneous amplitude and frequency of the Hilbert transform and a Power Spectral Density (PSD) based on fast Fourier transform;

s400, generating J training subsamples from a training sample set based on a Bootstrap algorithm;

s500, training J training sub-samples to obtain J candidate sub-models by adopting the same number of nuclear latent variables and nuclear parameters based on a Kernel Partial Least Squares (KPLS) algorithm;

s600, selecting all candidate submodels with corresponding model selection weight parameters larger than a model selection threshold value from the J candidate submodels according to a verification sample set to form an integrated submodel set, wherein the optimized weight parameters are obtained by optimizing a randomly generated initial weight parameter through a genetic algorithm with a minimized prediction error as a target;

s700, calculating and acquiring the integration weight of each sub-model in the integration sub-model set according to an AWF algorithm;

s800, obtaining the spectrum characteristics of the test data of the mill needing soft measurement;

and S900, calculating load parameters corresponding to the spectrum characteristics of the test data according to the integration submodels and the corresponding integration weights.

2. The method of claim 1, wherein S200 comprises:

s201, adding preset noise into the sample vibration signal and the sample vibration sound signal;

s202, acquiring M groups of initial sub-signals based on an EMD method;

s203, calculating the mean value of M groups of initial sub-signals to obtain J_IMFA sub-signal.

3. The method of claim 1, wherein the extracting the spectral features of all the sub-signals in S300 comprises:

s301, calculating MSHT, MVHT and PSD of the sub-signals;

s302, mutual information values between all MSHT spectrum variables of the MSHT and mill load parameters are calculated respectively;

s303, acquiring the maximum value and the minimum value of mutual information values of MSHT spectrum variables;

s304, obtaining a plurality of spectrum feature selection thresholds which change between the minimum value and the maximum value of mutual information values of MSHT spectrum variables in a preset step length, selecting all MSHT spectrum variables of which the corresponding mutual information values are larger than the spectrum feature selection thresholds based on each spectrum feature selection threshold to construct a partial least squares regression (PLS) model, and calculating the prediction error of the PLS model corresponding to each spectrum feature selection threshold;

s305, combining the MSHT spectrum variable set corresponding to the PLS model with the minimum prediction error into MSHT spectrum characteristics;

s306, calculating mutual information values between all MVHT spectrum variables of the MVHT and the mill load parameters respectively;

s307, acquiring the maximum value and the minimum value of mutual information values of MVHT spectrum variables;

s308, obtaining a plurality of spectrum feature selection thresholds which change between the minimum value and the maximum value of mutual information values of MVHT spectrum variables by a preset step length, selecting all the MVHT spectrum variables of which the corresponding mutual information values are larger than the spectrum feature selection thresholds based on each spectrum feature selection threshold to construct a PLS model, and calculating a PLS model prediction error corresponding to each spectrum feature selection threshold;

s309, combining the MVHT spectrum variable set corresponding to the PLS model with the minimum prediction error into an MVHT spectrum characteristic;

s310, mutual information values between all PSD spectrum variables of the PSD and the mill load parameters are calculated respectively;

s311, obtaining the maximum value and the minimum value of mutual information values of PSD spectrum variables;

s312, obtaining a plurality of spectrum feature selection threshold values which change between the minimum value and the maximum value of mutual information values of PSD spectrum variables in a preset step length, selecting all PSD spectrum variables of which the corresponding mutual information values are larger than the spectrum feature selection threshold value based on each spectrum feature selection threshold value to construct a PLS model, and calculating a PLS model prediction error corresponding to each spectrum feature selection threshold value;

s313, combining the PSD spectrum variable set corresponding to the PLS model with the minimum prediction error into PSD spectrum characteristics;

and S314, combining the MSHT spectral characteristics, the MVHT spectral characteristics and the PSD spectral characteristics to obtain the spectral characteristics of all the sub-signals.

4. The method of claim 1, wherein S700 comprises calculating weights for each sub-model in the integrated model according to the following formula:

<math> <mrow> <msubsup> <mi>W</mi> <msup> <mi>j</mi> <mo>*</mo> </msup> <mi>AWF</mi> </msubsup> <mo>=</mo> <mn>1</mn> <mo>/</mo> <mrow> <mo>(</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>&sigma;</mi> <msup> <mi>j</mi> <mo>*</mo> </msup> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <munderover> <mi>&Sigma;</mi> <mrow> <msup> <mi>j</mi> <mo>*</mo> </msup> <mo>=</mo> <mn>1</mn> </mrow> <msup> <mi>J</mi> <mo>*</mo> </msup> </munderover> <mfrac> <mn>1</mn> <msup> <mrow> <mo>(</mo> <msub> <mi>&sigma;</mi> <msup> <mi>j</mi> <mo>*</mo> </msup> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> <mo>)</mo> </mrow> </mrow> </math>

wherein,is the predicted value of the integrated sub-model to the verification sampleVariance of J^*In order to integrate the number of sub-models in the model,is the weight.

5. The method of claim 1, wherein S800 comprises:

s801, detecting a vibration signal to be measured and a vibration sound signal to be measured of a cylinder body of the mill to be measured;

s802, decomposing the vibration signal to be measured and the vibration sound signal to be measured into J based on EEMD method_IMFSample sub-signals, each sub-signal representing a single vibration mode having a physical meaning, J_IMFIs a predetermined value;

and S803, extracting the spectral characteristics of all the sub-signals, wherein the spectral characteristics are obtained by selecting the Marginal Spectrum (MSHT) based on the Hilbert transform, the mean value and the variance (MVHT) of the instantaneous amplitude and frequency of the Hilbert transform and the Power Spectral Density (PSD) based on the fast Fourier transform of the sub-signals.

6. The method of claim 1, wherein S600 comprises:

s601, calculating prediction output based on the jth candidate sub-model for each verification sample;

s602, calculating a prediction error based on the prediction output and the actual load parameter

S603, calculating correlation coefficients among all candidate submodels based on the following formula to form a correlation coefficient matrix;

<math> <mrow> <msubsup> <mi>c</mi> <mi>js</mi> <mi>valid</mi> </msubsup> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <msup> <mi>k</mi> <mi>valid</mi> </msup> </munderover> <msubsup> <mi>e</mi> <mi>valid</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>j</mi> <mo>,</mo> <msup> <mi>k</mi> <mi>valid</mi> </msup> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <msubsup> <mi>e</mi> <mi>valid</mi> <mi>s</mi> </msubsup> <mrow> <mo>(</mo> <mi>s</mi> <mo>,</mo> <msup> <mi>k</mi> <mi>valid</mi> </msup> <mo>)</mo> </mrow> <mo>/</mo> <msup> <mi>k</mi> <mi>valid</mi> </msup> </mrow> </math>

wherein,is the correlation coefficient, k, between the jth candidate submodel and the s candidate submodel^validTo verify the number of samples, the correlation coefficient matrix is:

<math> <mrow> <msup> <mi>C</mi> <mi>valid</mi> </msup> <mo>=</mo> <msub> <mfenced open='[' close=']'> <mtable> <mtr> <mtd> <msubsup> <mi>c</mi> <mn>11</mn> <mi>valid</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>c</mi> <mn>12</mn> <mi>valid</mi> </msubsup> </mtd> <mtd> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> </mtd> <mtd> <msubsup> <mi>c</mi> <mrow> <mn>1</mn> <mi>J</mi> </mrow> <mi>valid</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>c</mi> <mn>21</mn> <mi>valid</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>c</mi> <mn>22</mn> <mi>valid</mi> </msubsup> </mtd> <mtd> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> </mtd> <mtd> <msubsup> <mi>c</mi> <mrow> <mn>2</mn> <mi>J</mi> </mrow> <mi>valid</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> </mtd> <mtd> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> </mtd> <mtd> <msubsup> <mi>c</mi> <mi>js</mi> <mi>valid</mi> </msubsup> </mtd> <mtd> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>c</mi> <mrow> <mi>J</mi> <mn>1</mn> </mrow> <mi>valid</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>c</mi> <mrow> <mi>J</mi> <mn>2</mn> </mrow> <mi>valid</mi> </msubsup> </mtd> <mtd> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> </mtd> <mtd> <msubsup> <mi>c</mi> <mi>JJ</mi> <mi>valid</mi> </msubsup> </mtd> </mtr> </mtable> </mfenced> <mrow> <mi>J</mi> <mo>&times;</mo> <mi>J</mi> </mrow> </msub> </mrow> </math>

s604, generating a random vector for the candidate submodel

S605, using random vectorCalculating an optimized weight parameter which minimizes a prediction error by a genetic algorithm according to the correlation coefficient matrix as an initial weight parameter;

s606, selecting the candidate submodels with the corresponding optimization weight parameters larger than or equal to the preset selection threshold value as the integration submodels to obtain an integration submodel set.

7. The method of claim 1, wherein S900 includes predicting the load parameter according to the following formula:

<math> <mrow> <msub> <mover> <mi>y</mi> <mo>^</mo> </mover> <mi>test</mi> </msub> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <msup> <mi>j</mi> <mo>*</mo> </msup> <mo>=</mo> <mn>1</mn> </mrow> <msup> <mi>J</mi> <mo>*</mo> </msup> </munderover> <msubsup> <mi>W</mi> <msup> <mi>j</mi> <mo>*</mo> </msup> <mi>AWF</mi> </msubsup> <mo>&CenterDot;</mo> <msubsup> <mover> <mi>y</mi> <mo>^</mo> </mover> <mi>test</mi> <msup> <mi>j</mi> <mo>*</mo> </msup> </msubsup> </mrow> </math>

wherein,is based on the j^*th the predicted output of the sub-model,for the integration weight, J^*In order to integrate the number of sub-models in the model,load parameters obtained for the prediction.

8. The method of claim 1, wherein the mill is a ball mill and the loading parameters include one or more of a feed-to-ball ratio, a slurry concentration, and a packing fraction.