# Hybrid Intelligence Model Based on Image Features for the Prediction of Flotation Concentrate Grade

Authors: YaLin Wang; XiaoFang Chen; XiaoLing Zhou; WeiHua Gui; Louis Caccetta; Honglei Xu
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

In flotation processes, concentrate grade is the key production index but is difficult to be measured online. The mechanism models reflect the basic tendency of concentrate grade changes but cannot provide adequate prediction precision. The data-driven models based on froth image features provide accurate prediction within well-sampled space but rely heavily on sample data with less generalization capability. So, a hybrid intelligent model combining the two kinds of model is proposed in this paper. Since the information of image features is enormous, and the relationship between image features and concentrate grade is nonlinear, a B-spline partial least squares (BS-PLS) method is adopted to construct the data-driven model for concentrate grade prediction. In order to gain better generalization capability and prediction accuracy, information entropy is introduced to integrate the mechanism model and the BS-PLS model together and modify the model output online through an output deviation compensation strategy. Moreover, a slide window scheme is employed to update the hybrid model in order to improve its adaptability. The industrial practical data testing results show that the performance of the hybrid model is better than either of the two single models and it satisfies the accuracy and stability requirements in industrial applications.


---

## Body

## 1. Introduction

Concentrate grade is one of the most important production indices for froth flotation processes, which is generally measured offline through artificially timed testing of the content of valuable minerals. However, the time interval between each artificial test is too long to reflect the current state of flotation processes. Although online concentrate grade measurement devices are available (e.g., X-ray-fluorescence analysis), they usually require a significant amount of costly maintenance and calibration. Most of the flotation industries cannot afford such an expensive investment, which motivates the approaches for online grade prediction [1].Froth flotation is a kind of ore dressing process based on the different hydrophobic nature of its composed particles, involving complicated physicochemical reactions. The complexity of the process mainly arises from the inherently chaotic nature of the underlying microscopic phenomena [2]. Thus, it is difficult to build accurate mechanism models of the flotation process, not to mention predicting concentrate grade precisely based on flotation mechanism. So far, there are scarcely any papers related to mechanism modeling; however, there is a lot of research on flotation recovery mechanism [3–8] and (the relationship between recovery and concentrate grade [9, 10]). In this paper, the mechanism-based concentrate grade predicting model is built through substituting the recovery mechanism model into the relational model between recovery and concentrate grade proposed by [9]. Since the mechanism of flotation process is very intricate, all of the conclusions mentioned above are simplified based on rational hypothesis. Thus, the proposed model can only reflect the variance tendency of the process state but with low prediction accuracy.In order to gain better prediction accuracy of concentrate grade, it is feasible to exploit the froth image information. Experienced operators usually judge the flotation condition by the visual appearance of flotation froth, and earlier work has shown that the surface features of froth could reflect the production performance of froth flotation [2, 11, 12]. With the development of machine vision and digital image processing techniques, images of froth can be acquired and processed online. Combined with the extracted image information with enormous process variables (such as PH, concentration, flow rate, and froth height), more and more data-driven models have been developed to predict concentrate grade [13]. Data-driven models can predict the process indices through building a black box based on input and output data; related works are introduced as follows.NNs (neural networks) have been used in mining engineering applications and gained popularity since the 1990s. Since it can give a better performance than linear models [14], many works have been done to improve the prediction accuracy based on NNs model. Hargrave and Hall [11] used a feedforward network with one hidden layer consisting of four neurons to construct the relationship between the froth visual parameters and flotation performance characteristics. In this article, the inputs of the NN model only involve color and bubble size with a lot of other available information excluded; the prediction accuracy is not good enough. Nakhaei et al. [15] compare the prediction accuracy of the NN and multivariate nonlinear regression approaches, and the results showed that NN approaches offer a superior predictive capability over statistical method and meanwhile provide a smaller prediction error. Many approaches have also been found in the improvement of NNs. Wang et al. [16] proposed an improved BP network algorithm based on PCA to predict concentrate grade, which has been validated that it can solve the overmatching problem and improve the generalization ability of the BP network. Tu et al. [17] adopted Jordan network to establish concentrate grade prediction model which can predict the concentrate grade directly or in steps.Another modeling method as popular as NNs is the support vector machine (SVM), which has been proved that it is superior for applications with insufficient sample data. Geng and Chai [18] adopted least square support vector machine (LS-SVM) to estimate concentrate grade and tailing grade. Yang and Zhao [19] combined the adaptive principal component analysis (APCA) and composite kernel support vector regression (CK-SVR) to predict concentrate grade, and this method has a high estimation precision.However, traditional neural network and SVM are difficult to deal with the correlation problem when the input dimensions increase. In order to utilize the information obtained from froth image analysis, Hatonen et al. [2] proposed partial least squares (PLS) to predict the concentrate grade, which sufficiently include the useful image textures such as speed, stability, and bubble size. For concentrate grade prediction, Gonzalez et al. [20] used the root mean square error and the correlation coefficient as criteria to compare several classes of local dynamic models: nonlinear ARMAX, Takagi and Sugeno, fuzzy combination, projection on latent states (PLS), and wavelet based models. It is concluded that the PLS provides a better dynamic performance. Though PLS modeling could deal with the high input dimension problem, it is not sufficient enough for nonlinear application.The data-driven models mentioned above rely heavily on the sampled data. Therefore, when some fluctuation occurs during the sampling process, the prediction performance of data-driven models may be heavily impacted. Although the mechanism model can comparatively indicate the varying tendency of flotation state with the reflection of the inner mechanism of flotation, its prediction accuracy is low for the complexity of the flotation. It is obvious that the two types of models complement each other. Thus, to improve the concentrate prediction performance, we can combine the data-driven model and mechanism model together. The most popular approaches to integrate the two models are fuzzy classification and information entropy. Gui et al. [21] propose an integrated model for agglomerate compositions prediction, which combines the expertise-and-mechanism-based model and supervised distributed neural network model together, to help improve the blending ratios.In this paper, considering that PLS modeling is not suitable for the high dimensional and nonlinear problems in flotation, we propose BS-PLS modeling method combining B-spline interpolation with PLS to predict concentrate grade. To synthesize the advantage of data-driven model and mechanism model, information entropy is used to integrate BS-PLS model and mechanism model mentioned above. In order to improve the prediction accuracy of the hybrid model, an online output deviation compensation strategy is adopted to modify the final output of the model. Moreover, sliding window scheme is applied to update its structure so as to improve the robustness of the hybrid model.This paper is organized as follows. Section1 introduces research background of flotation concentrate grade prediction model. Section 2 briefly analyzes the flotation process and froth image features. Section 3 describes the hybrid model for concentrate grade prediction. Section 4 gives the implementation and result of the proposed hybrid model, and a conclusion of this paper is given in Section 5.


## 2. Flotation Process Description and Froth Image Features

### 2.1. Description of Flotation Process

Usually, the valuable minerals in the raw ore exploited are bonded together with gangue. To separate the valuable minerals from the gangue, many kinds of separation methods are developed, such as froth flotation, gravity, and electromagnetic, among which the froth flotation is the most widely used separation technique in the mineral industry.Froth flotation utilizes the physical chemistry discrepancy of different mineral particles to extract certain minerals from raw ore. With the reaction of reagents, minerals can be divided into two different categories: hydrophobic and hydrophilic species. Hydrophobic particles can attach to bubbles more easily, while most of the hydrophilic species will remain in the liquid; thus, it is feasible to separate valuable minerals by froth flotation. It is obvious that froth flotation is a complex physical and chemical process, composed of three phases: solid, liquid, and gas. The flotation process mechanism is shown in Figure1.Figure 1
Schematic diagram of flotation mechanism.Finely ground raw ores are entrained to the stirred vessel, where the ores are sufficiently mixed with the reagent, then fed into a rougher cell. With the effect of reagents, valuable minerals are made hydrophobic and the others hydrophilic. Since hydrophobic particles can attach to the bubbles more easily, we should generate a lot of bubbles in the bottom of the cell. Then, as a result of buoyancy, bubbles will rise to the surface of the cell, during which hydrophobic particles will attach to the bubbles, and will be carried to the surface together with bubbles and form a frothy layer. Finally, upper froth layer flows out of the cell forming concentrate product. Hydrophilic particles remain in the liquid and flow out of the bottom of the cell as tailings.For the flotation industry, a typical flotation flow sheet is composed of roughing, scavenging, and cleaning. A typical industry flotation flow sheet is shown in Figure2. The slurry from the grinding process is fed to roughing cell, which is the first operation of the process, in order to separate valuable minerals from the well mixed ores of stirred vessel. Cleaning is the next procedure of roughing, refining the product of rougher cell. For the purpose of improving final concentrate grade, flotation industry generally will set two or more cleaning operations. In order to ultimately recycle valuable ores, the tailings of rougher cell usually will be fed to a scavenging procedure. The product of this operation combined with the tailings of cleaning will be fed to rougher cell, and the tailings of scavenging are the final tailings of the whole flotation process. In some processes, a regrinding part may be included.Figure 2
Flotation process flow sheet.In the froth flotation process, improving the final concentrate grade is the key objective, while its measurement is time-consuming, as long as one or two hours. The offline measurement of concentrate grade deteriorates the flotation performance. Thus, online prediction of concentrate grade is important. For flotation processes, optimal control of roughing cell can maximize the recovery of valuable minerals from raw ores, with a certain concentrate grade maintained and the subsequent procedures benefited. It indicates that the roughing procedure is the most important part of the overall flotation flow sheet. So, a roughing cell is selected as the object.For a single cell, inspired by the froth flotation mechanism showed in Figure1, it can be separated into two parts: froth and slurry, as shown in Figure 3. This figure indicates that froth features are the most important indices for concentrate grade. Therefore, in this paper, we construct the concentrate grade prediction model for rougher cell based on froth image features.Figure 3
Single flotation system.


### 2.2. Froth Image Features

As mentioned before, froth image features can reflect concentrate grade to a great extent; thus, an image acquisition system is important for the froth flotation process. Figure4 shows a typical image acquisition system. A high-definition camera is installed 110 cm above the froth layer at a certain distance to collect froth images. Usually, an industrial flotation process runs all day long to ensure its profit, which means that the lighting condition varies with natural environment. So, a sealed box is needed to avoid the disturbance of the illumination variation, where an LED light with constant power is installed to provide stable illumination. The image information is transferred to industrial computer through fiber. Then, image features are extracted by digital processing technique.Figure 4
Image acquisition system.With the development of digital image processing technology, the extracted image features become more and more abundant, such as color, texture, speed, and structure. In this paper, four categories of fourteen groups features of bubble image are extracted; the extracting process is shown in Figure5. The specific extraction methods and detailed steps will be studied later.Figure 5
Features of extraction process.


#### 2.2.1. Color Features

It is commonly accepted that the color is the most widely used visual feature. Compared with other visual features, color is more robust since it is not sensitive to size and direction. For froth flotation, color features are especially useful for the prediction of concentrate grade.Generally, we use vector space forming a corresponding color space to represent color feature, and there exist various mathematical approaches to express the color space. The most extensively used method is the RGB (red, green, and blue) color space. RGB color space is calculated based on the trichromatic theory, which means that any color can be expressed by the combination of base colors of red, green, and blue. Consider(1)F=r[R]+g[G]+b[B],
where r, g, and b represent the functions of red, green, and blue, respectively.Equation (1) shows that F would vary with any of the RGB components, which means that F can express every kind of color. Moreover, relative redness is needed when adopting RGB color space to reflect color features of froth images, since it can weaken the effects of illumination variation. Consider
(2)Rrelative=RmeanGraymean,
where Rmean and Graymean represent the average values of red and grey levels, respectively.RGB color space exhibits obvious collinearity. Therefore, to express the color information more accurately, HSI color space is adopted in this paper. The three components of HSI are hue, saturation, and Intensity. HSI and RGB are the different expressions of the same physical quantity; there exists a transformational relationship between them. Consider(3)I=R+G+B3,S=I-3R+G+B[min⁡⁡(R,G,B)],H=cos⁡-1⁡[(1/2)((R-G)+(R-B))(R-G)2+(R-B)(G-B)].


#### 2.2.2. Velocity

Generally, based on the digital image, the velocity is estimated through tracking specific objects in the images with matching rules. In the flotation process, however, some deformation of froth bubbles usually occurs, such as rationing and scaling, which are the results of ceaseless impeller stirring in the flotation cells. The SIFT (scale invariant feature transform) method can be adopted to transform an image into a large collection of feature vectors, which are invariant to image rotation and illumination changes. Mu et al. [22] extracted froth velocity in various flotation states based on SIFT. The SIFT features extraction procedures are as follows.Step 1.
Scale-space peak extraction. In this procedure, a difference of Gaussians (DoG) is constructed to search the local peaks as interesting points. Consider
(4)D(x,y,σ)=(G(x,y,kσ)-G(x,y,σ))*I(x,y)=L(x,y,kσ)-L(x,y,σ),
where * represents convolution, G(x,y,σ) is the Gaussian convolution kernel, and k is a constant.Step 2.
Key points localization. At this stage, unstable points which are sensitive to noise are eliminated to achieve tolerance to noise.Step 3.
Orientation assignment. The dominant orientations for each key point can be fixed at the peak in a histogram of local image gradient orientation. The orientation of point (x,y) can be calculated using the following equation:
(5)θ(x,y)=tan-1⁡L(x,y+1)-L(x,y-1)L(x+1,y)-L(x-1,y).Step 4.
Key point description.


#### 2.2.3. Textural Features

Texture is an innate property of bubble surface, which is homogeneous and independent of color or illumination. It contains important information about the structural arrangement of surface and its relationship to the surrounding environment. There exist many approaches to represent texture features, among which grey level cooccurrence matrix (GLCM) is the most suitable approach for statistical calculations, which is first proposed by Haralick et al. [23].Since the cooccurring grey levels is defined for grey level correlations between pixels within a certain translation distance at a certain angle, it is feasible to construct GLCM to represent image texture. GLCM is a two-dimensional matrix, the size of which is the same as image grey levels [24]. The (i,j)th element of the cooccurrence matrix C is the sum of the times when pixel values i and j are offset from pixel intensity j, defined over an image X with n×m pixels. Thus, a grey level cooccurrence matrix is a matrix of frequencies. A cooccurrence matrix is constructed by repeating the calculations for each (i,j)th combination:
(6)C(i,j)=∑p=1n∑q=1mI(p,q).P(i,j)=C(i,j)∑i=1G∑j=1GC(i,j),A variety of image informative features can be obtained from the cooccurrence matrices, such as the following.Entropy:(7)Entr=∑i=0G-1∑j=0G-1Pd(i,j)log⁡2⁡Pd(i,j).Energy:(8)Ener=∑i=0G-1∑j=0G-1Pd2(i,j).Inertia:(9)Ine=∑i=0G-1∑j=0G-1[(i-j)2Pd(i,j)].All of these definitions can represent some textural aspect of the image.


#### 2.2.4. Structural Features

Here, the froth structural features refer to bubble size, which varies with the flotation operation conditions and directly affects the concentrate grade.A methodology called watershed is proposed to fully segment the froth image based on the idea of simulating water flows in a topographic representation of image intensity, as shown in Figure6.Figure 6
Schematic diagram of watershed algorithm.We can regard the grayscale image as a three-dimensional landscape. Each pixel of the image represents an altitude of the corresponding landscape, and the minima are called valleys. Suppose that each regional minimum has punched a hole. Put the model into a lake. Then, water is immersed from the bottom through the hole at a uniform rate, forming various regional lakes until the entire landscape has been fully immersed. When two lakes meet, the watershed can be identified. Lin et al. [25] emphasized that when applying watershed segmentation to highly detailed images, oversegmentation may arise.As shown in Figure7, white spots will appear when the grayscale values are high, which will considerably affect segmentation performance. On the contrary, the location where the grayscale value is the local minimum presents weak boundaries between bubbles. Based on this perspective, Yang et al. [26] proposed an improved valley edge detection algorithm (VED) to avoid edges around the white spots.Figure 7
Grey intensity distribution of bubble’s cross-section.The flow sheet of bubble size statistics is shown in Figure8; detailed information can be found in [27].Figure 8
Flow sheet of bubble size statistic.


### 2.3. Data Processing

The image features applied in this paper are 14 groups in total, including bubble size, velocity, stability, brightness, loading, energy, entropy, inertia, tone,R-means, G-means, B-means, relative redness, and saturation. In industrial froth flotation, the sampled original images may often be deteriorated because of sudden changes of operating conditions and the industrial environment. So, we need to preprocess the image features extracted from origin images.As we all know, generally, errors can be classified into gross error and stochastic error. The gross error can use a data outlier exclusion method to eliminate it. Here, we used the3σ criteria to distinguish the gross error; for a Gaussian distributed sample data, the probability of satisfying the mission requirements defined by ±3-σ limits is 0.9973%. In engineering terms, this probability is deemed acceptable.For a set of measurement data sequencexi(i=1,2,…,n), define its average value and standard deviation σ as follows:
(10)x-=1n∑i=1nxi,σ=1n-1∑i=1n(xi-x-)2.Ifxi-x->3σ, it is identified as an outlier.While, for the stochastic error, it is subject to stochastic laws, we can weaken or even eliminate it through signal denoising. There are different digital filtering methods for signal denoising such as mean filtering, moving average filtering, low-pass filtering, and wavelet filtering. The most commonly used method is the wavelet filtering, which is based on the theory that signal and noise present different properties under wavelet transformation. The wavelet filtering data denoising process is as follows. Firstly, the signals are decomposed into components distributed at different frequency range and time interval. Since the frequency range within which real signals are located is different from noise, we determine the one-scale wavelet coefficients corresponding to the noise’s frequency range and set a threshold to eliminate it. At last, through signal reconstruction, we can get the filtered signals.


## 2.1. Description of Flotation Process

Usually, the valuable minerals in the raw ore exploited are bonded together with gangue. To separate the valuable minerals from the gangue, many kinds of separation methods are developed, such as froth flotation, gravity, and electromagnetic, among which the froth flotation is the most widely used separation technique in the mineral industry.Froth flotation utilizes the physical chemistry discrepancy of different mineral particles to extract certain minerals from raw ore. With the reaction of reagents, minerals can be divided into two different categories: hydrophobic and hydrophilic species. Hydrophobic particles can attach to bubbles more easily, while most of the hydrophilic species will remain in the liquid; thus, it is feasible to separate valuable minerals by froth flotation. It is obvious that froth flotation is a complex physical and chemical process, composed of three phases: solid, liquid, and gas. The flotation process mechanism is shown in Figure1.Figure 1
Schematic diagram of flotation mechanism.Finely ground raw ores are entrained to the stirred vessel, where the ores are sufficiently mixed with the reagent, then fed into a rougher cell. With the effect of reagents, valuable minerals are made hydrophobic and the others hydrophilic. Since hydrophobic particles can attach to the bubbles more easily, we should generate a lot of bubbles in the bottom of the cell. Then, as a result of buoyancy, bubbles will rise to the surface of the cell, during which hydrophobic particles will attach to the bubbles, and will be carried to the surface together with bubbles and form a frothy layer. Finally, upper froth layer flows out of the cell forming concentrate product. Hydrophilic particles remain in the liquid and flow out of the bottom of the cell as tailings.For the flotation industry, a typical flotation flow sheet is composed of roughing, scavenging, and cleaning. A typical industry flotation flow sheet is shown in Figure2. The slurry from the grinding process is fed to roughing cell, which is the first operation of the process, in order to separate valuable minerals from the well mixed ores of stirred vessel. Cleaning is the next procedure of roughing, refining the product of rougher cell. For the purpose of improving final concentrate grade, flotation industry generally will set two or more cleaning operations. In order to ultimately recycle valuable ores, the tailings of rougher cell usually will be fed to a scavenging procedure. The product of this operation combined with the tailings of cleaning will be fed to rougher cell, and the tailings of scavenging are the final tailings of the whole flotation process. In some processes, a regrinding part may be included.Figure 2
Flotation process flow sheet.In the froth flotation process, improving the final concentrate grade is the key objective, while its measurement is time-consuming, as long as one or two hours. The offline measurement of concentrate grade deteriorates the flotation performance. Thus, online prediction of concentrate grade is important. For flotation processes, optimal control of roughing cell can maximize the recovery of valuable minerals from raw ores, with a certain concentrate grade maintained and the subsequent procedures benefited. It indicates that the roughing procedure is the most important part of the overall flotation flow sheet. So, a roughing cell is selected as the object.For a single cell, inspired by the froth flotation mechanism showed in Figure1, it can be separated into two parts: froth and slurry, as shown in Figure 3. This figure indicates that froth features are the most important indices for concentrate grade. Therefore, in this paper, we construct the concentrate grade prediction model for rougher cell based on froth image features.Figure 3
Single flotation system.


## 2.2. Froth Image Features

As mentioned before, froth image features can reflect concentrate grade to a great extent; thus, an image acquisition system is important for the froth flotation process. Figure4 shows a typical image acquisition system. A high-definition camera is installed 110 cm above the froth layer at a certain distance to collect froth images. Usually, an industrial flotation process runs all day long to ensure its profit, which means that the lighting condition varies with natural environment. So, a sealed box is needed to avoid the disturbance of the illumination variation, where an LED light with constant power is installed to provide stable illumination. The image information is transferred to industrial computer through fiber. Then, image features are extracted by digital processing technique.Figure 4
Image acquisition system.With the development of digital image processing technology, the extracted image features become more and more abundant, such as color, texture, speed, and structure. In this paper, four categories of fourteen groups features of bubble image are extracted; the extracting process is shown in Figure5. The specific extraction methods and detailed steps will be studied later.Figure 5
Features of extraction process.


### 2.2.1. Color Features

It is commonly accepted that the color is the most widely used visual feature. Compared with other visual features, color is more robust since it is not sensitive to size and direction. For froth flotation, color features are especially useful for the prediction of concentrate grade.Generally, we use vector space forming a corresponding color space to represent color feature, and there exist various mathematical approaches to express the color space. The most extensively used method is the RGB (red, green, and blue) color space. RGB color space is calculated based on the trichromatic theory, which means that any color can be expressed by the combination of base colors of red, green, and blue. Consider(1)F=r[R]+g[G]+b[B],
where r, g, and b represent the functions of red, green, and blue, respectively.Equation (1) shows that F would vary with any of the RGB components, which means that F can express every kind of color. Moreover, relative redness is needed when adopting RGB color space to reflect color features of froth images, since it can weaken the effects of illumination variation. Consider
(2)Rrelative=RmeanGraymean,
where Rmean and Graymean represent the average values of red and grey levels, respectively.RGB color space exhibits obvious collinearity. Therefore, to express the color information more accurately, HSI color space is adopted in this paper. The three components of HSI are hue, saturation, and Intensity. HSI and RGB are the different expressions of the same physical quantity; there exists a transformational relationship between them. Consider(3)I=R+G+B3,S=I-3R+G+B[min⁡⁡(R,G,B)],H=cos⁡-1⁡[(1/2)((R-G)+(R-B))(R-G)2+(R-B)(G-B)].


### 2.2.2. Velocity

Generally, based on the digital image, the velocity is estimated through tracking specific objects in the images with matching rules. In the flotation process, however, some deformation of froth bubbles usually occurs, such as rationing and scaling, which are the results of ceaseless impeller stirring in the flotation cells. The SIFT (scale invariant feature transform) method can be adopted to transform an image into a large collection of feature vectors, which are invariant to image rotation and illumination changes. Mu et al. [22] extracted froth velocity in various flotation states based on SIFT. The SIFT features extraction procedures are as follows.Step 1.
Scale-space peak extraction. In this procedure, a difference of Gaussians (DoG) is constructed to search the local peaks as interesting points. Consider
(4)D(x,y,σ)=(G(x,y,kσ)-G(x,y,σ))*I(x,y)=L(x,y,kσ)-L(x,y,σ),
where * represents convolution, G(x,y,σ) is the Gaussian convolution kernel, and k is a constant.Step 2.
Key points localization. At this stage, unstable points which are sensitive to noise are eliminated to achieve tolerance to noise.Step 3.
Orientation assignment. The dominant orientations for each key point can be fixed at the peak in a histogram of local image gradient orientation. The orientation of point (x,y) can be calculated using the following equation:
(5)θ(x,y)=tan-1⁡L(x,y+1)-L(x,y-1)L(x+1,y)-L(x-1,y).Step 4.
Key point description.


### 2.2.3. Textural Features

Texture is an innate property of bubble surface, which is homogeneous and independent of color or illumination. It contains important information about the structural arrangement of surface and its relationship to the surrounding environment. There exist many approaches to represent texture features, among which grey level cooccurrence matrix (GLCM) is the most suitable approach for statistical calculations, which is first proposed by Haralick et al. [23].Since the cooccurring grey levels is defined for grey level correlations between pixels within a certain translation distance at a certain angle, it is feasible to construct GLCM to represent image texture. GLCM is a two-dimensional matrix, the size of which is the same as image grey levels [24]. The (i,j)th element of the cooccurrence matrix C is the sum of the times when pixel values i and j are offset from pixel intensity j, defined over an image X with n×m pixels. Thus, a grey level cooccurrence matrix is a matrix of frequencies. A cooccurrence matrix is constructed by repeating the calculations for each (i,j)th combination:
(6)C(i,j)=∑p=1n∑q=1mI(p,q).P(i,j)=C(i,j)∑i=1G∑j=1GC(i,j),A variety of image informative features can be obtained from the cooccurrence matrices, such as the following.Entropy:(7)Entr=∑i=0G-1∑j=0G-1Pd(i,j)log⁡2⁡Pd(i,j).Energy:(8)Ener=∑i=0G-1∑j=0G-1Pd2(i,j).Inertia:(9)Ine=∑i=0G-1∑j=0G-1[(i-j)2Pd(i,j)].All of these definitions can represent some textural aspect of the image.


### 2.2.4. Structural Features

Here, the froth structural features refer to bubble size, which varies with the flotation operation conditions and directly affects the concentrate grade.A methodology called watershed is proposed to fully segment the froth image based on the idea of simulating water flows in a topographic representation of image intensity, as shown in Figure6.Figure 6
Schematic diagram of watershed algorithm.We can regard the grayscale image as a three-dimensional landscape. Each pixel of the image represents an altitude of the corresponding landscape, and the minima are called valleys. Suppose that each regional minimum has punched a hole. Put the model into a lake. Then, water is immersed from the bottom through the hole at a uniform rate, forming various regional lakes until the entire landscape has been fully immersed. When two lakes meet, the watershed can be identified. Lin et al. [25] emphasized that when applying watershed segmentation to highly detailed images, oversegmentation may arise.As shown in Figure7, white spots will appear when the grayscale values are high, which will considerably affect segmentation performance. On the contrary, the location where the grayscale value is the local minimum presents weak boundaries between bubbles. Based on this perspective, Yang et al. [26] proposed an improved valley edge detection algorithm (VED) to avoid edges around the white spots.Figure 7
Grey intensity distribution of bubble’s cross-section.The flow sheet of bubble size statistics is shown in Figure8; detailed information can be found in [27].Figure 8
Flow sheet of bubble size statistic.


## 2.2.1. Color Features

It is commonly accepted that the color is the most widely used visual feature. Compared with other visual features, color is more robust since it is not sensitive to size and direction. For froth flotation, color features are especially useful for the prediction of concentrate grade.Generally, we use vector space forming a corresponding color space to represent color feature, and there exist various mathematical approaches to express the color space. The most extensively used method is the RGB (red, green, and blue) color space. RGB color space is calculated based on the trichromatic theory, which means that any color can be expressed by the combination of base colors of red, green, and blue. Consider(1)F=r[R]+g[G]+b[B],
where r, g, and b represent the functions of red, green, and blue, respectively.Equation (1) shows that F would vary with any of the RGB components, which means that F can express every kind of color. Moreover, relative redness is needed when adopting RGB color space to reflect color features of froth images, since it can weaken the effects of illumination variation. Consider
(2)Rrelative=RmeanGraymean,
where Rmean and Graymean represent the average values of red and grey levels, respectively.RGB color space exhibits obvious collinearity. Therefore, to express the color information more accurately, HSI color space is adopted in this paper. The three components of HSI are hue, saturation, and Intensity. HSI and RGB are the different expressions of the same physical quantity; there exists a transformational relationship between them. Consider(3)I=R+G+B3,S=I-3R+G+B[min⁡⁡(R,G,B)],H=cos⁡-1⁡[(1/2)((R-G)+(R-B))(R-G)2+(R-B)(G-B)].


## 2.2.2. Velocity

Generally, based on the digital image, the velocity is estimated through tracking specific objects in the images with matching rules. In the flotation process, however, some deformation of froth bubbles usually occurs, such as rationing and scaling, which are the results of ceaseless impeller stirring in the flotation cells. The SIFT (scale invariant feature transform) method can be adopted to transform an image into a large collection of feature vectors, which are invariant to image rotation and illumination changes. Mu et al. [22] extracted froth velocity in various flotation states based on SIFT. The SIFT features extraction procedures are as follows.Step 1.
Scale-space peak extraction. In this procedure, a difference of Gaussians (DoG) is constructed to search the local peaks as interesting points. Consider
(4)D(x,y,σ)=(G(x,y,kσ)-G(x,y,σ))*I(x,y)=L(x,y,kσ)-L(x,y,σ),
where * represents convolution, G(x,y,σ) is the Gaussian convolution kernel, and k is a constant.Step 2.
Key points localization. At this stage, unstable points which are sensitive to noise are eliminated to achieve tolerance to noise.Step 3.
Orientation assignment. The dominant orientations for each key point can be fixed at the peak in a histogram of local image gradient orientation. The orientation of point (x,y) can be calculated using the following equation:
(5)θ(x,y)=tan-1⁡L(x,y+1)-L(x,y-1)L(x+1,y)-L(x-1,y).Step 4.
Key point description.


## 2.2.3. Textural Features

Texture is an innate property of bubble surface, which is homogeneous and independent of color or illumination. It contains important information about the structural arrangement of surface and its relationship to the surrounding environment. There exist many approaches to represent texture features, among which grey level cooccurrence matrix (GLCM) is the most suitable approach for statistical calculations, which is first proposed by Haralick et al. [23].Since the cooccurring grey levels is defined for grey level correlations between pixels within a certain translation distance at a certain angle, it is feasible to construct GLCM to represent image texture. GLCM is a two-dimensional matrix, the size of which is the same as image grey levels [24]. The (i,j)th element of the cooccurrence matrix C is the sum of the times when pixel values i and j are offset from pixel intensity j, defined over an image X with n×m pixels. Thus, a grey level cooccurrence matrix is a matrix of frequencies. A cooccurrence matrix is constructed by repeating the calculations for each (i,j)th combination:
(6)C(i,j)=∑p=1n∑q=1mI(p,q).P(i,j)=C(i,j)∑i=1G∑j=1GC(i,j),A variety of image informative features can be obtained from the cooccurrence matrices, such as the following.Entropy:(7)Entr=∑i=0G-1∑j=0G-1Pd(i,j)log⁡2⁡Pd(i,j).Energy:(8)Ener=∑i=0G-1∑j=0G-1Pd2(i,j).Inertia:(9)Ine=∑i=0G-1∑j=0G-1[(i-j)2Pd(i,j)].All of these definitions can represent some textural aspect of the image.


## 2.2.4. Structural Features

Here, the froth structural features refer to bubble size, which varies with the flotation operation conditions and directly affects the concentrate grade.A methodology called watershed is proposed to fully segment the froth image based on the idea of simulating water flows in a topographic representation of image intensity, as shown in Figure6.Figure 6
Schematic diagram of watershed algorithm.We can regard the grayscale image as a three-dimensional landscape. Each pixel of the image represents an altitude of the corresponding landscape, and the minima are called valleys. Suppose that each regional minimum has punched a hole. Put the model into a lake. Then, water is immersed from the bottom through the hole at a uniform rate, forming various regional lakes until the entire landscape has been fully immersed. When two lakes meet, the watershed can be identified. Lin et al. [25] emphasized that when applying watershed segmentation to highly detailed images, oversegmentation may arise.As shown in Figure7, white spots will appear when the grayscale values are high, which will considerably affect segmentation performance. On the contrary, the location where the grayscale value is the local minimum presents weak boundaries between bubbles. Based on this perspective, Yang et al. [26] proposed an improved valley edge detection algorithm (VED) to avoid edges around the white spots.Figure 7
Grey intensity distribution of bubble’s cross-section.The flow sheet of bubble size statistics is shown in Figure8; detailed information can be found in [27].Figure 8
Flow sheet of bubble size statistic.


## 2.3. Data Processing

The image features applied in this paper are 14 groups in total, including bubble size, velocity, stability, brightness, loading, energy, entropy, inertia, tone,R-means, G-means, B-means, relative redness, and saturation. In industrial froth flotation, the sampled original images may often be deteriorated because of sudden changes of operating conditions and the industrial environment. So, we need to preprocess the image features extracted from origin images.As we all know, generally, errors can be classified into gross error and stochastic error. The gross error can use a data outlier exclusion method to eliminate it. Here, we used the3σ criteria to distinguish the gross error; for a Gaussian distributed sample data, the probability of satisfying the mission requirements defined by ±3-σ limits is 0.9973%. In engineering terms, this probability is deemed acceptable.For a set of measurement data sequencexi(i=1,2,…,n), define its average value and standard deviation σ as follows:
(10)x-=1n∑i=1nxi,σ=1n-1∑i=1n(xi-x-)2.Ifxi-x->3σ, it is identified as an outlier.While, for the stochastic error, it is subject to stochastic laws, we can weaken or even eliminate it through signal denoising. There are different digital filtering methods for signal denoising such as mean filtering, moving average filtering, low-pass filtering, and wavelet filtering. The most commonly used method is the wavelet filtering, which is based on the theory that signal and noise present different properties under wavelet transformation. The wavelet filtering data denoising process is as follows. Firstly, the signals are decomposed into components distributed at different frequency range and time interval. Since the frequency range within which real signals are located is different from noise, we determine the one-scale wavelet coefficients corresponding to the noise’s frequency range and set a threshold to eliminate it. At last, through signal reconstruction, we can get the filtered signals.


## 3. Hybrid Intelligent Model

The structure of hybrid modeling is shown in Figure9. It consists of three modules: hybrid modeling, slide window update strategy, and output deviation compensation. Since the recovery of a single flotation cell can be estimated based on mechanical analysis, we can construct a mechanism model through relating concentrate grade and recovery to predict concentrate grade. While the inner reactions of flotation are complex, inevitably the mechanical analysis of recovery should be based on some hypothesis; thus, the prediction accuracy of this method will be affected. Statistical analysis of the image information could compensate the disadvantage of mechanism modeling. BS-PLS modeling is appropriate to gain the statistical relationship between the 14 image features and the corresponding concentrate grade. Therefore, we adopted information entropy (IE) to integrate the mechanism model and BS-PLS model, which can distribute the weight of the two models based on their predicting error. To correct the output of hybrid model online, an output deviation compensation scheme is used. In the flotation process, the feed conditions and other disturbances will gradually vary with time; to improve the adaptability of the proposed hybrid model, an update strategy is needed. Detailed instruction of each module shown in Figure 9 will be discussed subsequently.Figure 9
The structure of hybrid modeling.Using information entropy, the output of BS-PLS modeling and mechanism modelingy^1, y^2 is integrated, gaining the output of hybrid model y^. To ensure that y^ will not deviate from laboratory values, adopt online output deviation compensation method to modify the value of y^ at time point k with errors between model output and laboratory value E(k). The value λ in Figure 9 represents the weight of error E(k) at present time point k and the previous time point k-1. The variable p in the slide window update strategy module refers to the fixed time interval when the model structure should be updated.


### 3.1. Mechanism Modeling

In the flotation process, froth recovery and concentrate grade are the most important indices for its performance. The procedures to build a mechanism model for grade prediction are as follows.Firstly, give an estimation of the overall recovery. Then, find out the relationship between grade and recovery, with which we can predict the concentrate grade according to the estimation value of recovery.


#### 3.1.1. Flotation Overall Recovery Estimation

Flotation process can be represented by a two-stage model consisting of the pulp zone and froth layer. The pulp zone is related to the collection process, and the froth layer represents the separation process, as shown in Figure10. The mass flow rate in the mineral transport streams is described as F: feed, C: concentrate, T: tailings, B: bubble-particle aggregate, E: entrainment, and D: drop-back.Figure 10
Mineral transport streams in a flotation cell.Figure10 shows that the overall flotation process can be divided into the following parts: selective transfer of material from slurry to the froth by particle-bubble attachment; nonselective transfer of material from the slurry to the froth by entrainment; drop-back (both selective and nonselective) of material from the froth to the slurry; and mechanical or hydraulic transfer of material from the froth into the concentrate.In the opinion of mass balance, the overall recovery can be divided into collection zone recovery and froth zone recovery, as shown in Figure11. The overall recovery, for true flotation (minerals are carried to the froth layer by bubbles) or entrained minerals (minerals are carried to the froth layer by pulp), can be described by the following relationship:
(11)R=RCRF1-RC(1-RF),
where R is the overall flotation recovery, RC is the collection zone recovery, and RF is the froth zone recovery. Since the major part of valuable minerals collected in the concentrate is from true flotation, the discussion about the flotation process in the following text is all about true flotation.Figure 11
Pulp and froth recovery model.Neethling [28] proposed a simple approximation for estimating froth recovery RF, which mainly considers the bubble loading and the fraction of particles attached to the film that become detached when it ruptures. According to this article, the froth recovery could be represented as follows:
(12)RF={(αvgvset)f/2(rinrout)fα<12(vg2vset)f/2(rinrout)fα≥12,
where α is the air recovery, vg is the gas velocity, vset is the particle settling velocity, rin is the initial average bubble radius into the froth zone, rout is flowing output average bubble radius, and f is the detach rate of particles, which can be estimated from surface bubble loading:
(13)f=1-Lsurf.Air recovery and particle settling velocity can be estimated by the following equations:(14)α=vghlQ,vset=2gr2(ρp-ρf)9η,
where h is the froth height, l is the notch length, Q is the air flux, g is the gravitational acceleration, ρp and ρf represent the density of solid particles and the medium, respectively, and η is the viscosity of medium.For most flotation systems, the air recovery is lower than 50%. Thus, generally, we applyα<1/2 for the industrial field. Substituting (13), (14) into (12), we can get the final froth recovery form:
(15)RF=(vgroutβ)1-Lsurf,(16)β=9ηhlrin22g(ρp-ρf)r2Q.For (16), we usually regard β as a constant when the working conditions are stable.Polat and Chander [29] proposed that the collection zone recovery at time t will be equal to
(17)RC=R∞(1-∫0∞F(k)exp⁡⁡(-kt)dk).
Here, R∞ is the ultimate recovery at long times, k is the flotation rate constant, and F(k) is the distribution function of flotation float rate k. Generally, we can simplify it with a rectangular distribution. Then, the corresponding recovery can be obtained:
(18)RC=R∞(1-1tK+e-tKtK),
where K is the max flotation rate constant; it is decided by a specific flotation process.From (18), we can see that the collection zone recovery relies on flotation time t and flotation float rate constant k. In this paper, according to the industrial field data, the flotation time t=16.5 min. The flotation float rate constant heavily depends on the properties of particles. The distribution of mineral particles in flotation cell is complex; thus, to simplify the calculation, we set k=1.21 where we can fit the collection zone recovery well. According to a longtime sampling of industrial flotation process, R∞ is approximate to 0.91. Substitute (15) and (18) into (11); then,
(19)R=0.885((vg/rout)β)1-Lsurf0.885((vg/rout)β)1-Lsurf+0.115,
where vg, rout, and Lsurf are the froth image features that can be extracted by digital image processing technique. Meanwhile, β can be identified based on vg, rout, Lsurf,  and the corresponding recovery values; detail procedure will be introduced subsequently.


#### 3.1.2. The Grade-Recovery Relationship

Neethling and Cilliers [9] proposed an approach to estimate the grade-recovery relationship. Their methodology is based on a single assumption: the separation occurring in the froth and the rate at which solids enter the froth (the pulp zone kinetics) are uniquely defined by the conditions prevalent in the pulp, which implies that the concentrate grade produced by any cell (Gins) is a function just relying on operating conditions and the various concentrations of minerals in the pulp. Further, in well-mixed flotation cells, the concentration of each particle type in the pulp is the same as its concentrations in the tailings.The concentrate grade can then be expressed as a function of cell operating conditions and tailing concentrations:(20)Gins=f(CellOperatingConditions;mmmPulpConcentrations)≈f(CellOperatingConditions;mmmTailingConcentrations).For a certain condition, it can be simplified as(21)Gins=f(PulpConcentrations)=f(TailingConcentrations).Equation (21) indicates that instantaneous concentrate grade closely connects with tailing concentration. Based on the mass balance theory, we can calculate the tailing grade Gtail from the feed grade Gfeed, the cumulative recovery R, and the cumulative grade G:
(22)Gtail=Gfeed(1-R)(1-(Gfeed/G)R).Then, theGins can be expressed as
(23)Gins=f(Gfeed(1-R)1-(Gfeed/G)R).For the batch flotation test, the concentrate grade can be defined as follows:(24)G=VM,
where V is the total mass of valuable mineral recovered and M is the total mass recovered in the concentrate.The instantaneous concentrate grade in batch flotationGins is the ratio between the rate of recovery of valuable mineral and the rate of total mass recovery:
(25)Gins=dV/dtdM/dt=(dV/dM)(dM/dt)dM/dt=dVdM.Combining (24) and (25), we can obtain
(26)Gins=G+MdGdM=G+MdRdMdGdR.IfSinitial is the initial mass of solids in the batch flotation cell, then
(27)G=SinitialRGfeedM.Combining (26) and (27) yields
(28)dGdR=SinitialGfeedM+RdSinitialGfeeddM=GR-RSinitialGfeed1M2dMdR=GR(Gins-GGins)=GR(1-GGins).As mentioned before,Gins can be expressed as a function of Gtail, which can be obtained through a statistical fit based on sampled data. Since the solution to (28) is implicit, it can be solved readily using the Runge-Kutta method with initial values t=0, R=0, and G=Gins. Finally, we can get the relationship between concentrate grade and recovery.When we get a frame of froth picture, we can extract the image features to estimate the flotation recovery. Then, based on the relationship between concentrate grade and recovery, we can predict the current concentrate grade.


### 3.2. BS-PLS Modeling

In most of the flotation industry, experienced operators usually adjust the operating conditions by observing the state of the bubble on the surface. The state of surface bubble reflects significant information about the current operating conditions. Thus, we can get the statistical relationship between froth image features and flotation indices, with which the concentrate grade is predicted. With the development of modern digital processing technology, enormous image features extracted lead to the increase of input dimension. Models with high input dimension may run into colinearity problem. Research has shown that PLS can sufficiently utilize the large input image information; meanwhile, it can eliminate the colinearity problem. PLS is a good choice for flotation concentrate prediction. However, the relationship between froth image features and concentrate grade is nonlinear, while PLS is developed for linear modeling. To solve this problem, we combine B-spline with PLS (BS-PLS) to deal with nonlinear modeling. With the help of BS-PLS, we can form a nonlinear equation between the image inputs and concentrate grade of the flotation process.


#### 3.2.1. B-Spline Theory

The term spline was first proposed by Schoenberg in 1946, and, with the development of data fitting and function approximation in the industrial field, spline theory gradually became mature. As the spline function is too complex to be used for most applications, therefore, a relatively simpler and more practical function is needed. In this paper, B-spline function can provide a satisfactory performance.Divide an intervalx∈[a,b] into n subintervals a≤x0<x1<⋯<xn-1<xn≤b, and, for each subinterval [xi-1,xi], i=0,…,n, there exists continuous function ΩK(x) of K order. Here, ΩK(x) is the B-spline function; it is defined as follows:
(29)ΩK(x)=1K!∑k=0K+1(-1)k(K+1k)(x+K+12-k)+K,
where
(30)(x+K+12-k)+K={(x+K+12-k)K(x≥k-K+12)0(x<k-K+12).
We call εk=k-(K+1)/2, k=0,1,…,K the internal node. k refers to the kth component of function ΩK(x); K represents the order of function ΩK(x).


#### 3.2.2. Partial Least Squares (PLS) Regression Theory

Advances in digital image processing improve the quality of froth image feature extraction, which leads to a problem of high input dimensionality when using image features as model input, directly affecting the model performance. Partial least squares (PLS) regression is effectively used in process modeling and monitoring to deal with this kind of problem. It firstly orthogonally projects inputs onto one-dimensional latent variables and then gets the regression. A cross-validation strategy is frequently used to choose the number of latent variables or factors applied in the model so that the prediction variance is minimized. The calculation procedure is as follows.Step 1.
Standardized data[X,Y]=[E0,F0]; then, extract (t1,u1). Here, t1 and u1 represent the first component of E0 and F0, respectively,  t1=E0w1, ||w1||=1, u1=F0c1, and ||c1=1||. In order to illustrate the information of X and Y, we should maximize the covariance of vectors t1 and u1; thus, we can transform this problem as follows:
(31)max⁡〈E0w1,F0c1〉s.t.{w1Tw1=1c1Tc1=1.To solve this optimization problem, use a Lagrange algorithm; that is, introduce an equation:(32)s=w1TE0TF0c1-λ1(w1Tw1-1)-λ2(c1Tc1-1).By partial deviation operations of (32), we obtain
(33)∂s∂w1=E0TF0c1-2λ1w1=0,(34)∂s∂c1=F0TE0w1-2λ2c1=0,(35)∂s∂λ1=-(w1Tw1-1)=0,(36)∂s∂λ2=-(c1Tc1-1)=0.Then, we can obtain the following formula:(37)θ1=2λ1=2λ2=w1TE0TF0c1=〈w1E0,F0c1〉.Equation (37) indicates that the optimization problem can be transformed into maximized parameter θ1. Combining (33) and (34), we obtain the following equations:
(38)E0TF0F0TE0w1=θ12w1,F0TE0E0TF0c1=θ12c1.The vectorsw1 and c1 are the eigenvectors of E0TF0F0TE0 and F0TE0E0TF0, respectively; θ12 is the largest corresponding eigenvalue.Then, the components can be obtained:(39)t1=E0w1,u1=F0c1.The regression equations are then obtained(40)E0=t1p1T+E1,p1=E0Tt1∥t1∥2F0=t1r1T+F1,r1=F0Tt1∥t1∥2,
where E1 and F1 are the residual matrixes of regression equations.Step 2.
SubstituteE0, F0 with E1 and F1, and repeat Step 1. Finally, we have
(41)E0=t1p1T+t2p2T+⋯+tApAT,F0=t1r1T+t2r2T+⋯+tArAT+FA,
where A=r(X) and F is the residual matrix.During the PLS modeling, the number of principle components can generally be determined by cross-validation.


#### 3.2.3. BS-PLS Modeling

Firstly, use B-spline function to transform the model input vectorsxj. Consider
(42)f^j(xj)=β0j+∑l=0Mj+2βj,lΩ3(xj-εj,l-1hj),Ω3(xj-εj,l-1hj)=13!hj3∑p=04(-1)p(4p)(xj-εj,l-3+p)+3,
where β0j and βj,l are the unknown coefficients and εj,l-1, hj, and Mj represent the division interval, interval length, and number of intervals for variable xj, respectively:
(43)εj,l-1=min⁡⁡(xj)+(l-1)hj,hj=max⁡⁡(xj)-min⁡⁡(xj)Mj,
where l=0,1,…,Mj+2.Then, the nonlinear relationship between input and output can be expressed as follows:(44)y=β0+∑j=1pf^j(xj)=β0+∑j=1p∑l=0Mj+2βj,lΩ3(xj-εj,l-1hj)+ε.Finally, use PLS regression to identify the unknown coefficients.The specific steps for BS-PLS modeling are as follows.Step 1.
Transformxj→Zj for each column of the input matrix. First, decide the number of intervals Mj, and then calculate εj,l-1 according to (43); then, carry out B-spline transforming for [xj]n×1:(45)Zj={zj,0,zj,0,…,zj,Mj+2}={Ω3(xj-εj,l-1hj),l=0,1,…,Mj+2}.Step 2.
Normalize the input matrix (transformed by B-spline function) and the output matrix (not transformed). Consider(46)z~j,li=zj,li-z-j,lsj,l,y~i=(yi-y-)sy,
where l=0,1,…,Mj+2; j=1,2,…,p; i=1,2,…,n, z-j,l, y- represent the mean value of zj,l and y, respectively, and sj,l and sy represent the variance of zj,l and y, respectively. Then, it becomes a normal linear PLS problem as follows:
(47)y~=∑j=1p∑l=0Mj+2αj,lz~j,l,Step 3.
Use PLS regression to identify the parametersαj,l(j=1,2,…,p; l=0,1,…,Mj+2).Step 4.
Substitute (46) into (47); then,
(48)y-y-sy=∑j=1p∑l=0Mj+2αj,lzj,l-z-j,lsj,l.And the linear formula can be expressed as follows:(49)y=β0+∑j=1p∑l=0Mj+2βj,lzj,l,
where
(50)βj,l=syαj,lsj,l,β0=y--∑j=1p∑l=0Mj+2βj,lz-j,l.Step 5.
Substitute (45) and the regression parameters into (49); then, the final nonlinear model can be expressed as follows:
(51)y^=β0+∑j=1pf^j(xj)=β0+∑j=1p∑l=0Mj+2βj,lΩ3(xj-εj,l-1hj).


### 3.3. Hybrid Modeling Method

To integrate the mechanism model and BS-PLS model, the reasonable weight distribution of each one in the hybrid model is very important. The information theory indicates that the amount of information is related to its disorder degree. As clarified by Shannon’s information entropy theory on the quantified disorder degree, a larger information entropy interprets a higher disorder degree and in this case produces a lower weight in the hybrid model. Thereby, we can apply information entropy to quantify the disorder degree of each model based on its sequential predicting error and calculate their weights. The specific procedures are as follows.Step 1.
Normalization. First, we should define the relative prediction error of models:
(52)eij={|yi-y^ijyi|,0≤|yi-y^ijyi|<11,|yi-y^ijyi|≥1.It represents the relative prediction error of thejth model at the ith sampling time, and then uses the following formula to normalize the errors:
(53)pij=eij∑j=1neij(i=1,…,n,j=1,2).Step 2.
Compute entropyEj. Consider
(54)Ej=-k∑i=1npijln⁡pijj=1,2,
where the constant k equals (ln⁡n)−1.Step 3.
Setdj=1-Ej(j=1,2) as the degree of diversification; then, calculate the degree of importance for each model:
(55)wj=1m-1(1-dj∑j=1mdj)j=1,2.Step 4.
Calculate the results of the hybrid model:(56)y^i=∑j=12wjy^iji=1,…,n.


### 3.4. Online Updating Scheme for the Hybrid Model

Froth flotation is a long process, and the operating conditions are always changing. As mentioned before, the parameterβ in the mechanism model needs to be identified based on vg, rout, Lsurf, and the corresponding recovery value. The nonlinear structure of BS-PLS modeling is determined by sample data. Therefore, the proposed hybrid model relies on history sample data, which will cause the model to be invalid for future prediction. Thus, we need to update the hybrid model and substitute old data with new data with a certain strategy.


#### 3.4.1. Sliding Window Update Strategy

A sliding window is a kind of update strategy where only the most “recent” elements remain active and the rest are discarded. There are two equally important types of sliding window update strategy. One is updating with fixed-size windows, by which a fixed amount of the recent elements being active is defined (e.g., items arrive one at a time, and only the most recentn items remain active for some fixed parameter n). The other is updating with timestamp-based windows, by which the validity of an element is defined by an additional parameter such as a timestamp (e.g., many items can arrive in “bursts” at a single step and where only items from the last t steps remain active again for some fixed parameter t).Fixed-size windows are important for applications where the arrival rate of the data is fixed (but still extremely fast), such as sensors or stock market measurements. Timestamp-based windows are important for applications with asynchronous data arrivals, such as networking or database applications.In this paper, the hybrid model is built on the image texture information, and the images are processed at a fixed rate; thus, the arrival rate of the data for the model is fixed, which means that the fixed-size windows are suitable for this work.Since the hybrid model is composed of two parts, mechanism model and BS-PLS model, they both need to be updated at the same time.The parameter that needs to be updated for mechanism model is the constantβ according to (19). Therefore, for mechanism model, when the window slides one step forward, recalculate the Rave, rave, Lsav, and vgav within the current window; then, get the renewed value of β. Thus far, the mechanism model update is completed.For BS-PLS model, as it is built upon the sampled data, the most important problem is how to deal with these data. According to Section3, the structure of the BS-PLS model is determined by the zj (which is the B-spline transform of xj) and the variance S. In order to consider the history data (i.e., the data in the previous window) adequately, the update strategy for the average and variance of the current window data is as follows:
(57)z-k+1=N-PNz-k+PNz~k+1,Sk+12=N-P-1N-1Sk2+PN-1(z~k+1-z-k+1)2,
where z-k, Sk are the average and variance of the previous window data, respectively, and z~k+1 is the average of the current window data. Then, the update strategy is applied to reidentify the BS-PLS structure.Since the two composing models are changed, their entropy weight should also be updated along with them, which means that every time the window slides, after the two composed models are updated, we should recalculate the entropy weight according to Section3.


#### 3.4.2. Online Correction Compensation Strategy

In a real process operation, the operating conditions may fluctuate frequently. In order to enable the hybrid model to have the adaptive ability to overcome the system disturbances and the drifting of the operation point, a practical online correction algorithm is devised.The correction algorithm takes the errors between the estimation of the hybrid model and the laboratory testing data as the correction drive and calculates a suitable compensation value to append to hybrid model prediction. The correction strategy is as follows:(58)E(k)=λ×E0(k)+(1-λ)×E(k-1),E0(k)=ylab(k)-ymod⁡(k),ycor(k)=ymod⁡(k)+E(k),
where E(k) is the modified value for errors between the hybrid model and the laboratory value. E0(k) and E(k-1) are the error of the current and one-step previous hybrid model prediction and the testing data, respectively. λ is the weight constant of the two errors. ylab is the laboratory testing value of concentrate grade, ymod⁡ is the model prediction for concentrate grade, and ycor is the final output of hybrid model compensated by the corrected errors.Through the correction compensation strategy proposed above, the output of the hybrid model can be modified online, and as a result the accuracy of the model prediction can be improved.The flow chart of overall modified scheme is shown in Figure12. The variable p refers to the fixed time interval when the model structure should be updated. Variable k represents the kth sampling data.Figure 12
Flow chart of overall modified scheme.


## 3.1. Mechanism Modeling

In the flotation process, froth recovery and concentrate grade are the most important indices for its performance. The procedures to build a mechanism model for grade prediction are as follows.Firstly, give an estimation of the overall recovery. Then, find out the relationship between grade and recovery, with which we can predict the concentrate grade according to the estimation value of recovery.


### 3.1.1. Flotation Overall Recovery Estimation

Flotation process can be represented by a two-stage model consisting of the pulp zone and froth layer. The pulp zone is related to the collection process, and the froth layer represents the separation process, as shown in Figure10. The mass flow rate in the mineral transport streams is described as F: feed, C: concentrate, T: tailings, B: bubble-particle aggregate, E: entrainment, and D: drop-back.Figure 10
Mineral transport streams in a flotation cell.Figure10 shows that the overall flotation process can be divided into the following parts: selective transfer of material from slurry to the froth by particle-bubble attachment; nonselective transfer of material from the slurry to the froth by entrainment; drop-back (both selective and nonselective) of material from the froth to the slurry; and mechanical or hydraulic transfer of material from the froth into the concentrate.In the opinion of mass balance, the overall recovery can be divided into collection zone recovery and froth zone recovery, as shown in Figure11. The overall recovery, for true flotation (minerals are carried to the froth layer by bubbles) or entrained minerals (minerals are carried to the froth layer by pulp), can be described by the following relationship:
(11)R=RCRF1-RC(1-RF),
where R is the overall flotation recovery, RC is the collection zone recovery, and RF is the froth zone recovery. Since the major part of valuable minerals collected in the concentrate is from true flotation, the discussion about the flotation process in the following text is all about true flotation.Figure 11
Pulp and froth recovery model.Neethling [28] proposed a simple approximation for estimating froth recovery RF, which mainly considers the bubble loading and the fraction of particles attached to the film that become detached when it ruptures. According to this article, the froth recovery could be represented as follows:
(12)RF={(αvgvset)f/2(rinrout)fα<12(vg2vset)f/2(rinrout)fα≥12,
where α is the air recovery, vg is the gas velocity, vset is the particle settling velocity, rin is the initial average bubble radius into the froth zone, rout is flowing output average bubble radius, and f is the detach rate of particles, which can be estimated from surface bubble loading:
(13)f=1-Lsurf.Air recovery and particle settling velocity can be estimated by the following equations:(14)α=vghlQ,vset=2gr2(ρp-ρf)9η,
where h is the froth height, l is the notch length, Q is the air flux, g is the gravitational acceleration, ρp and ρf represent the density of solid particles and the medium, respectively, and η is the viscosity of medium.For most flotation systems, the air recovery is lower than 50%. Thus, generally, we applyα<1/2 for the industrial field. Substituting (13), (14) into (12), we can get the final froth recovery form:
(15)RF=(vgroutβ)1-Lsurf,(16)β=9ηhlrin22g(ρp-ρf)r2Q.For (16), we usually regard β as a constant when the working conditions are stable.Polat and Chander [29] proposed that the collection zone recovery at time t will be equal to
(17)RC=R∞(1-∫0∞F(k)exp⁡⁡(-kt)dk).
Here, R∞ is the ultimate recovery at long times, k is the flotation rate constant, and F(k) is the distribution function of flotation float rate k. Generally, we can simplify it with a rectangular distribution. Then, the corresponding recovery can be obtained:
(18)RC=R∞(1-1tK+e-tKtK),
where K is the max flotation rate constant; it is decided by a specific flotation process.From (18), we can see that the collection zone recovery relies on flotation time t and flotation float rate constant k. In this paper, according to the industrial field data, the flotation time t=16.5 min. The flotation float rate constant heavily depends on the properties of particles. The distribution of mineral particles in flotation cell is complex; thus, to simplify the calculation, we set k=1.21 where we can fit the collection zone recovery well. According to a longtime sampling of industrial flotation process, R∞ is approximate to 0.91. Substitute (15) and (18) into (11); then,
(19)R=0.885((vg/rout)β)1-Lsurf0.885((vg/rout)β)1-Lsurf+0.115,
where vg, rout, and Lsurf are the froth image features that can be extracted by digital image processing technique. Meanwhile, β can be identified based on vg, rout, Lsurf,  and the corresponding recovery values; detail procedure will be introduced subsequently.


### 3.1.2. The Grade-Recovery Relationship

Neethling and Cilliers [9] proposed an approach to estimate the grade-recovery relationship. Their methodology is based on a single assumption: the separation occurring in the froth and the rate at which solids enter the froth (the pulp zone kinetics) are uniquely defined by the conditions prevalent in the pulp, which implies that the concentrate grade produced by any cell (Gins) is a function just relying on operating conditions and the various concentrations of minerals in the pulp. Further, in well-mixed flotation cells, the concentration of each particle type in the pulp is the same as its concentrations in the tailings.The concentrate grade can then be expressed as a function of cell operating conditions and tailing concentrations:(20)Gins=f(CellOperatingConditions;mmmPulpConcentrations)≈f(CellOperatingConditions;mmmTailingConcentrations).For a certain condition, it can be simplified as(21)Gins=f(PulpConcentrations)=f(TailingConcentrations).Equation (21) indicates that instantaneous concentrate grade closely connects with tailing concentration. Based on the mass balance theory, we can calculate the tailing grade Gtail from the feed grade Gfeed, the cumulative recovery R, and the cumulative grade G:
(22)Gtail=Gfeed(1-R)(1-(Gfeed/G)R).Then, theGins can be expressed as
(23)Gins=f(Gfeed(1-R)1-(Gfeed/G)R).For the batch flotation test, the concentrate grade can be defined as follows:(24)G=VM,
where V is the total mass of valuable mineral recovered and M is the total mass recovered in the concentrate.The instantaneous concentrate grade in batch flotationGins is the ratio between the rate of recovery of valuable mineral and the rate of total mass recovery:
(25)Gins=dV/dtdM/dt=(dV/dM)(dM/dt)dM/dt=dVdM.Combining (24) and (25), we can obtain
(26)Gins=G+MdGdM=G+MdRdMdGdR.IfSinitial is the initial mass of solids in the batch flotation cell, then
(27)G=SinitialRGfeedM.Combining (26) and (27) yields
(28)dGdR=SinitialGfeedM+RdSinitialGfeeddM=GR-RSinitialGfeed1M2dMdR=GR(Gins-GGins)=GR(1-GGins).As mentioned before,Gins can be expressed as a function of Gtail, which can be obtained through a statistical fit based on sampled data. Since the solution to (28) is implicit, it can be solved readily using the Runge-Kutta method with initial values t=0, R=0, and G=Gins. Finally, we can get the relationship between concentrate grade and recovery.When we get a frame of froth picture, we can extract the image features to estimate the flotation recovery. Then, based on the relationship between concentrate grade and recovery, we can predict the current concentrate grade.


## 3.1.1. Flotation Overall Recovery Estimation

Flotation process can be represented by a two-stage model consisting of the pulp zone and froth layer. The pulp zone is related to the collection process, and the froth layer represents the separation process, as shown in Figure10. The mass flow rate in the mineral transport streams is described as F: feed, C: concentrate, T: tailings, B: bubble-particle aggregate, E: entrainment, and D: drop-back.Figure 10
Mineral transport streams in a flotation cell.Figure10 shows that the overall flotation process can be divided into the following parts: selective transfer of material from slurry to the froth by particle-bubble attachment; nonselective transfer of material from the slurry to the froth by entrainment; drop-back (both selective and nonselective) of material from the froth to the slurry; and mechanical or hydraulic transfer of material from the froth into the concentrate.In the opinion of mass balance, the overall recovery can be divided into collection zone recovery and froth zone recovery, as shown in Figure11. The overall recovery, for true flotation (minerals are carried to the froth layer by bubbles) or entrained minerals (minerals are carried to the froth layer by pulp), can be described by the following relationship:
(11)R=RCRF1-RC(1-RF),
where R is the overall flotation recovery, RC is the collection zone recovery, and RF is the froth zone recovery. Since the major part of valuable minerals collected in the concentrate is from true flotation, the discussion about the flotation process in the following text is all about true flotation.Figure 11
Pulp and froth recovery model.Neethling [28] proposed a simple approximation for estimating froth recovery RF, which mainly considers the bubble loading and the fraction of particles attached to the film that become detached when it ruptures. According to this article, the froth recovery could be represented as follows:
(12)RF={(αvgvset)f/2(rinrout)fα<12(vg2vset)f/2(rinrout)fα≥12,
where α is the air recovery, vg is the gas velocity, vset is the particle settling velocity, rin is the initial average bubble radius into the froth zone, rout is flowing output average bubble radius, and f is the detach rate of particles, which can be estimated from surface bubble loading:
(13)f=1-Lsurf.Air recovery and particle settling velocity can be estimated by the following equations:(14)α=vghlQ,vset=2gr2(ρp-ρf)9η,
where h is the froth height, l is the notch length, Q is the air flux, g is the gravitational acceleration, ρp and ρf represent the density of solid particles and the medium, respectively, and η is the viscosity of medium.For most flotation systems, the air recovery is lower than 50%. Thus, generally, we applyα<1/2 for the industrial field. Substituting (13), (14) into (12), we can get the final froth recovery form:
(15)RF=(vgroutβ)1-Lsurf,(16)β=9ηhlrin22g(ρp-ρf)r2Q.For (16), we usually regard β as a constant when the working conditions are stable.Polat and Chander [29] proposed that the collection zone recovery at time t will be equal to
(17)RC=R∞(1-∫0∞F(k)exp⁡⁡(-kt)dk).
Here, R∞ is the ultimate recovery at long times, k is the flotation rate constant, and F(k) is the distribution function of flotation float rate k. Generally, we can simplify it with a rectangular distribution. Then, the corresponding recovery can be obtained:
(18)RC=R∞(1-1tK+e-tKtK),
where K is the max flotation rate constant; it is decided by a specific flotation process.From (18), we can see that the collection zone recovery relies on flotation time t and flotation float rate constant k. In this paper, according to the industrial field data, the flotation time t=16.5 min. The flotation float rate constant heavily depends on the properties of particles. The distribution of mineral particles in flotation cell is complex; thus, to simplify the calculation, we set k=1.21 where we can fit the collection zone recovery well. According to a longtime sampling of industrial flotation process, R∞ is approximate to 0.91. Substitute (15) and (18) into (11); then,
(19)R=0.885((vg/rout)β)1-Lsurf0.885((vg/rout)β)1-Lsurf+0.115,
where vg, rout, and Lsurf are the froth image features that can be extracted by digital image processing technique. Meanwhile, β can be identified based on vg, rout, Lsurf,  and the corresponding recovery values; detail procedure will be introduced subsequently.


## 3.1.2. The Grade-Recovery Relationship

Neethling and Cilliers [9] proposed an approach to estimate the grade-recovery relationship. Their methodology is based on a single assumption: the separation occurring in the froth and the rate at which solids enter the froth (the pulp zone kinetics) are uniquely defined by the conditions prevalent in the pulp, which implies that the concentrate grade produced by any cell (Gins) is a function just relying on operating conditions and the various concentrations of minerals in the pulp. Further, in well-mixed flotation cells, the concentration of each particle type in the pulp is the same as its concentrations in the tailings.The concentrate grade can then be expressed as a function of cell operating conditions and tailing concentrations:(20)Gins=f(CellOperatingConditions;mmmPulpConcentrations)≈f(CellOperatingConditions;mmmTailingConcentrations).For a certain condition, it can be simplified as(21)Gins=f(PulpConcentrations)=f(TailingConcentrations).Equation (21) indicates that instantaneous concentrate grade closely connects with tailing concentration. Based on the mass balance theory, we can calculate the tailing grade Gtail from the feed grade Gfeed, the cumulative recovery R, and the cumulative grade G:
(22)Gtail=Gfeed(1-R)(1-(Gfeed/G)R).Then, theGins can be expressed as
(23)Gins=f(Gfeed(1-R)1-(Gfeed/G)R).For the batch flotation test, the concentrate grade can be defined as follows:(24)G=VM,
where V is the total mass of valuable mineral recovered and M is the total mass recovered in the concentrate.The instantaneous concentrate grade in batch flotationGins is the ratio between the rate of recovery of valuable mineral and the rate of total mass recovery:
(25)Gins=dV/dtdM/dt=(dV/dM)(dM/dt)dM/dt=dVdM.Combining (24) and (25), we can obtain
(26)Gins=G+MdGdM=G+MdRdMdGdR.IfSinitial is the initial mass of solids in the batch flotation cell, then
(27)G=SinitialRGfeedM.Combining (26) and (27) yields
(28)dGdR=SinitialGfeedM+RdSinitialGfeeddM=GR-RSinitialGfeed1M2dMdR=GR(Gins-GGins)=GR(1-GGins).As mentioned before,Gins can be expressed as a function of Gtail, which can be obtained through a statistical fit based on sampled data. Since the solution to (28) is implicit, it can be solved readily using the Runge-Kutta method with initial values t=0, R=0, and G=Gins. Finally, we can get the relationship between concentrate grade and recovery.When we get a frame of froth picture, we can extract the image features to estimate the flotation recovery. Then, based on the relationship between concentrate grade and recovery, we can predict the current concentrate grade.


## 3.2. BS-PLS Modeling

In most of the flotation industry, experienced operators usually adjust the operating conditions by observing the state of the bubble on the surface. The state of surface bubble reflects significant information about the current operating conditions. Thus, we can get the statistical relationship between froth image features and flotation indices, with which the concentrate grade is predicted. With the development of modern digital processing technology, enormous image features extracted lead to the increase of input dimension. Models with high input dimension may run into colinearity problem. Research has shown that PLS can sufficiently utilize the large input image information; meanwhile, it can eliminate the colinearity problem. PLS is a good choice for flotation concentrate prediction. However, the relationship between froth image features and concentrate grade is nonlinear, while PLS is developed for linear modeling. To solve this problem, we combine B-spline with PLS (BS-PLS) to deal with nonlinear modeling. With the help of BS-PLS, we can form a nonlinear equation between the image inputs and concentrate grade of the flotation process.


### 3.2.1. B-Spline Theory

The term spline was first proposed by Schoenberg in 1946, and, with the development of data fitting and function approximation in the industrial field, spline theory gradually became mature. As the spline function is too complex to be used for most applications, therefore, a relatively simpler and more practical function is needed. In this paper, B-spline function can provide a satisfactory performance.Divide an intervalx∈[a,b] into n subintervals a≤x0<x1<⋯<xn-1<xn≤b, and, for each subinterval [xi-1,xi], i=0,…,n, there exists continuous function ΩK(x) of K order. Here, ΩK(x) is the B-spline function; it is defined as follows:
(29)ΩK(x)=1K!∑k=0K+1(-1)k(K+1k)(x+K+12-k)+K,
where
(30)(x+K+12-k)+K={(x+K+12-k)K(x≥k-K+12)0(x<k-K+12).
We call εk=k-(K+1)/2, k=0,1,…,K the internal node. k refers to the kth component of function ΩK(x); K represents the order of function ΩK(x).


### 3.2.2. Partial Least Squares (PLS) Regression Theory

Advances in digital image processing improve the quality of froth image feature extraction, which leads to a problem of high input dimensionality when using image features as model input, directly affecting the model performance. Partial least squares (PLS) regression is effectively used in process modeling and monitoring to deal with this kind of problem. It firstly orthogonally projects inputs onto one-dimensional latent variables and then gets the regression. A cross-validation strategy is frequently used to choose the number of latent variables or factors applied in the model so that the prediction variance is minimized. The calculation procedure is as follows.Step 1.
Standardized data[X,Y]=[E0,F0]; then, extract (t1,u1). Here, t1 and u1 represent the first component of E0 and F0, respectively,  t1=E0w1, ||w1||=1, u1=F0c1, and ||c1=1||. In order to illustrate the information of X and Y, we should maximize the covariance of vectors t1 and u1; thus, we can transform this problem as follows:
(31)max⁡〈E0w1,F0c1〉s.t.{w1Tw1=1c1Tc1=1.To solve this optimization problem, use a Lagrange algorithm; that is, introduce an equation:(32)s=w1TE0TF0c1-λ1(w1Tw1-1)-λ2(c1Tc1-1).By partial deviation operations of (32), we obtain
(33)∂s∂w1=E0TF0c1-2λ1w1=0,(34)∂s∂c1=F0TE0w1-2λ2c1=0,(35)∂s∂λ1=-(w1Tw1-1)=0,(36)∂s∂λ2=-(c1Tc1-1)=0.Then, we can obtain the following formula:(37)θ1=2λ1=2λ2=w1TE0TF0c1=〈w1E0,F0c1〉.Equation (37) indicates that the optimization problem can be transformed into maximized parameter θ1. Combining (33) and (34), we obtain the following equations:
(38)E0TF0F0TE0w1=θ12w1,F0TE0E0TF0c1=θ12c1.The vectorsw1 and c1 are the eigenvectors of E0TF0F0TE0 and F0TE0E0TF0, respectively; θ12 is the largest corresponding eigenvalue.Then, the components can be obtained:(39)t1=E0w1,u1=F0c1.The regression equations are then obtained(40)E0=t1p1T+E1,p1=E0Tt1∥t1∥2F0=t1r1T+F1,r1=F0Tt1∥t1∥2,
where E1 and F1 are the residual matrixes of regression equations.Step 2.
SubstituteE0, F0 with E1 and F1, and repeat Step 1. Finally, we have
(41)E0=t1p1T+t2p2T+⋯+tApAT,F0=t1r1T+t2r2T+⋯+tArAT+FA,
where A=r(X) and F is the residual matrix.During the PLS modeling, the number of principle components can generally be determined by cross-validation.


### 3.2.3. BS-PLS Modeling

Firstly, use B-spline function to transform the model input vectorsxj. Consider
(42)f^j(xj)=β0j+∑l=0Mj+2βj,lΩ3(xj-εj,l-1hj),Ω3(xj-εj,l-1hj)=13!hj3∑p=04(-1)p(4p)(xj-εj,l-3+p)+3,
where β0j and βj,l are the unknown coefficients and εj,l-1, hj, and Mj represent the division interval, interval length, and number of intervals for variable xj, respectively:
(43)εj,l-1=min⁡⁡(xj)+(l-1)hj,hj=max⁡⁡(xj)-min⁡⁡(xj)Mj,
where l=0,1,…,Mj+2.Then, the nonlinear relationship between input and output can be expressed as follows:(44)y=β0+∑j=1pf^j(xj)=β0+∑j=1p∑l=0Mj+2βj,lΩ3(xj-εj,l-1hj)+ε.Finally, use PLS regression to identify the unknown coefficients.The specific steps for BS-PLS modeling are as follows.Step 1.
Transformxj→Zj for each column of the input matrix. First, decide the number of intervals Mj, and then calculate εj,l-1 according to (43); then, carry out B-spline transforming for [xj]n×1:(45)Zj={zj,0,zj,0,…,zj,Mj+2}={Ω3(xj-εj,l-1hj),l=0,1,…,Mj+2}.Step 2.
Normalize the input matrix (transformed by B-spline function) and the output matrix (not transformed). Consider(46)z~j,li=zj,li-z-j,lsj,l,y~i=(yi-y-)sy,
where l=0,1,…,Mj+2; j=1,2,…,p; i=1,2,…,n, z-j,l, y- represent the mean value of zj,l and y, respectively, and sj,l and sy represent the variance of zj,l and y, respectively. Then, it becomes a normal linear PLS problem as follows:
(47)y~=∑j=1p∑l=0Mj+2αj,lz~j,l,Step 3.
Use PLS regression to identify the parametersαj,l(j=1,2,…,p; l=0,1,…,Mj+2).Step 4.
Substitute (46) into (47); then,
(48)y-y-sy=∑j=1p∑l=0Mj+2αj,lzj,l-z-j,lsj,l.And the linear formula can be expressed as follows:(49)y=β0+∑j=1p∑l=0Mj+2βj,lzj,l,
where
(50)βj,l=syαj,lsj,l,β0=y--∑j=1p∑l=0Mj+2βj,lz-j,l.Step 5.
Substitute (45) and the regression parameters into (49); then, the final nonlinear model can be expressed as follows:
(51)y^=β0+∑j=1pf^j(xj)=β0+∑j=1p∑l=0Mj+2βj,lΩ3(xj-εj,l-1hj).


## 3.2.1. B-Spline Theory

The term spline was first proposed by Schoenberg in 1946, and, with the development of data fitting and function approximation in the industrial field, spline theory gradually became mature. As the spline function is too complex to be used for most applications, therefore, a relatively simpler and more practical function is needed. In this paper, B-spline function can provide a satisfactory performance.Divide an intervalx∈[a,b] into n subintervals a≤x0<x1<⋯<xn-1<xn≤b, and, for each subinterval [xi-1,xi], i=0,…,n, there exists continuous function ΩK(x) of K order. Here, ΩK(x) is the B-spline function; it is defined as follows:
(29)ΩK(x)=1K!∑k=0K+1(-1)k(K+1k)(x+K+12-k)+K,
where
(30)(x+K+12-k)+K={(x+K+12-k)K(x≥k-K+12)0(x<k-K+12).
We call εk=k-(K+1)/2, k=0,1,…,K the internal node. k refers to the kth component of function ΩK(x); K represents the order of function ΩK(x).


## 3.2.2. Partial Least Squares (PLS) Regression Theory

Advances in digital image processing improve the quality of froth image feature extraction, which leads to a problem of high input dimensionality when using image features as model input, directly affecting the model performance. Partial least squares (PLS) regression is effectively used in process modeling and monitoring to deal with this kind of problem. It firstly orthogonally projects inputs onto one-dimensional latent variables and then gets the regression. A cross-validation strategy is frequently used to choose the number of latent variables or factors applied in the model so that the prediction variance is minimized. The calculation procedure is as follows.Step 1.
Standardized data[X,Y]=[E0,F0]; then, extract (t1,u1). Here, t1 and u1 represent the first component of E0 and F0, respectively,  t1=E0w1, ||w1||=1, u1=F0c1, and ||c1=1||. In order to illustrate the information of X and Y, we should maximize the covariance of vectors t1 and u1; thus, we can transform this problem as follows:
(31)max⁡〈E0w1,F0c1〉s.t.{w1Tw1=1c1Tc1=1.To solve this optimization problem, use a Lagrange algorithm; that is, introduce an equation:(32)s=w1TE0TF0c1-λ1(w1Tw1-1)-λ2(c1Tc1-1).By partial deviation operations of (32), we obtain
(33)∂s∂w1=E0TF0c1-2λ1w1=0,(34)∂s∂c1=F0TE0w1-2λ2c1=0,(35)∂s∂λ1=-(w1Tw1-1)=0,(36)∂s∂λ2=-(c1Tc1-1)=0.Then, we can obtain the following formula:(37)θ1=2λ1=2λ2=w1TE0TF0c1=〈w1E0,F0c1〉.Equation (37) indicates that the optimization problem can be transformed into maximized parameter θ1. Combining (33) and (34), we obtain the following equations:
(38)E0TF0F0TE0w1=θ12w1,F0TE0E0TF0c1=θ12c1.The vectorsw1 and c1 are the eigenvectors of E0TF0F0TE0 and F0TE0E0TF0, respectively; θ12 is the largest corresponding eigenvalue.Then, the components can be obtained:(39)t1=E0w1,u1=F0c1.The regression equations are then obtained(40)E0=t1p1T+E1,p1=E0Tt1∥t1∥2F0=t1r1T+F1,r1=F0Tt1∥t1∥2,
where E1 and F1 are the residual matrixes of regression equations.Step 2.
SubstituteE0, F0 with E1 and F1, and repeat Step 1. Finally, we have
(41)E0=t1p1T+t2p2T+⋯+tApAT,F0=t1r1T+t2r2T+⋯+tArAT+FA,
where A=r(X) and F is the residual matrix.During the PLS modeling, the number of principle components can generally be determined by cross-validation.


## 3.2.3. BS-PLS Modeling

Firstly, use B-spline function to transform the model input vectorsxj. Consider
(42)f^j(xj)=β0j+∑l=0Mj+2βj,lΩ3(xj-εj,l-1hj),Ω3(xj-εj,l-1hj)=13!hj3∑p=04(-1)p(4p)(xj-εj,l-3+p)+3,
where β0j and βj,l are the unknown coefficients and εj,l-1, hj, and Mj represent the division interval, interval length, and number of intervals for variable xj, respectively:
(43)εj,l-1=min⁡⁡(xj)+(l-1)hj,hj=max⁡⁡(xj)-min⁡⁡(xj)Mj,
where l=0,1,…,Mj+2.Then, the nonlinear relationship between input and output can be expressed as follows:(44)y=β0+∑j=1pf^j(xj)=β0+∑j=1p∑l=0Mj+2βj,lΩ3(xj-εj,l-1hj)+ε.Finally, use PLS regression to identify the unknown coefficients.The specific steps for BS-PLS modeling are as follows.Step 1.
Transformxj→Zj for each column of the input matrix. First, decide the number of intervals Mj, and then calculate εj,l-1 according to (43); then, carry out B-spline transforming for [xj]n×1:(45)Zj={zj,0,zj,0,…,zj,Mj+2}={Ω3(xj-εj,l-1hj),l=0,1,…,Mj+2}.Step 2.
Normalize the input matrix (transformed by B-spline function) and the output matrix (not transformed). Consider(46)z~j,li=zj,li-z-j,lsj,l,y~i=(yi-y-)sy,
where l=0,1,…,Mj+2; j=1,2,…,p; i=1,2,…,n, z-j,l, y- represent the mean value of zj,l and y, respectively, and sj,l and sy represent the variance of zj,l and y, respectively. Then, it becomes a normal linear PLS problem as follows:
(47)y~=∑j=1p∑l=0Mj+2αj,lz~j,l,Step 3.
Use PLS regression to identify the parametersαj,l(j=1,2,…,p; l=0,1,…,Mj+2).Step 4.
Substitute (46) into (47); then,
(48)y-y-sy=∑j=1p∑l=0Mj+2αj,lzj,l-z-j,lsj,l.And the linear formula can be expressed as follows:(49)y=β0+∑j=1p∑l=0Mj+2βj,lzj,l,
where
(50)βj,l=syαj,lsj,l,β0=y--∑j=1p∑l=0Mj+2βj,lz-j,l.Step 5.
Substitute (45) and the regression parameters into (49); then, the final nonlinear model can be expressed as follows:
(51)y^=β0+∑j=1pf^j(xj)=β0+∑j=1p∑l=0Mj+2βj,lΩ3(xj-εj,l-1hj).


## 3.3. Hybrid Modeling Method

To integrate the mechanism model and BS-PLS model, the reasonable weight distribution of each one in the hybrid model is very important. The information theory indicates that the amount of information is related to its disorder degree. As clarified by Shannon’s information entropy theory on the quantified disorder degree, a larger information entropy interprets a higher disorder degree and in this case produces a lower weight in the hybrid model. Thereby, we can apply information entropy to quantify the disorder degree of each model based on its sequential predicting error and calculate their weights. The specific procedures are as follows.Step 1.
Normalization. First, we should define the relative prediction error of models:
(52)eij={|yi-y^ijyi|,0≤|yi-y^ijyi|<11,|yi-y^ijyi|≥1.It represents the relative prediction error of thejth model at the ith sampling time, and then uses the following formula to normalize the errors:
(53)pij=eij∑j=1neij(i=1,…,n,j=1,2).Step 2.
Compute entropyEj. Consider
(54)Ej=-k∑i=1npijln⁡pijj=1,2,
where the constant k equals (ln⁡n)−1.Step 3.
Setdj=1-Ej(j=1,2) as the degree of diversification; then, calculate the degree of importance for each model:
(55)wj=1m-1(1-dj∑j=1mdj)j=1,2.Step 4.
Calculate the results of the hybrid model:(56)y^i=∑j=12wjy^iji=1,…,n.


## 3.4. Online Updating Scheme for the Hybrid Model

Froth flotation is a long process, and the operating conditions are always changing. As mentioned before, the parameterβ in the mechanism model needs to be identified based on vg, rout, Lsurf, and the corresponding recovery value. The nonlinear structure of BS-PLS modeling is determined by sample data. Therefore, the proposed hybrid model relies on history sample data, which will cause the model to be invalid for future prediction. Thus, we need to update the hybrid model and substitute old data with new data with a certain strategy.


### 3.4.1. Sliding Window Update Strategy

A sliding window is a kind of update strategy where only the most “recent” elements remain active and the rest are discarded. There are two equally important types of sliding window update strategy. One is updating with fixed-size windows, by which a fixed amount of the recent elements being active is defined (e.g., items arrive one at a time, and only the most recentn items remain active for some fixed parameter n). The other is updating with timestamp-based windows, by which the validity of an element is defined by an additional parameter such as a timestamp (e.g., many items can arrive in “bursts” at a single step and where only items from the last t steps remain active again for some fixed parameter t).Fixed-size windows are important for applications where the arrival rate of the data is fixed (but still extremely fast), such as sensors or stock market measurements. Timestamp-based windows are important for applications with asynchronous data arrivals, such as networking or database applications.In this paper, the hybrid model is built on the image texture information, and the images are processed at a fixed rate; thus, the arrival rate of the data for the model is fixed, which means that the fixed-size windows are suitable for this work.Since the hybrid model is composed of two parts, mechanism model and BS-PLS model, they both need to be updated at the same time.The parameter that needs to be updated for mechanism model is the constantβ according to (19). Therefore, for mechanism model, when the window slides one step forward, recalculate the Rave, rave, Lsav, and vgav within the current window; then, get the renewed value of β. Thus far, the mechanism model update is completed.For BS-PLS model, as it is built upon the sampled data, the most important problem is how to deal with these data. According to Section3, the structure of the BS-PLS model is determined by the zj (which is the B-spline transform of xj) and the variance S. In order to consider the history data (i.e., the data in the previous window) adequately, the update strategy for the average and variance of the current window data is as follows:
(57)z-k+1=N-PNz-k+PNz~k+1,Sk+12=N-P-1N-1Sk2+PN-1(z~k+1-z-k+1)2,
where z-k, Sk are the average and variance of the previous window data, respectively, and z~k+1 is the average of the current window data. Then, the update strategy is applied to reidentify the BS-PLS structure.Since the two composing models are changed, their entropy weight should also be updated along with them, which means that every time the window slides, after the two composed models are updated, we should recalculate the entropy weight according to Section3.


### 3.4.2. Online Correction Compensation Strategy

In a real process operation, the operating conditions may fluctuate frequently. In order to enable the hybrid model to have the adaptive ability to overcome the system disturbances and the drifting of the operation point, a practical online correction algorithm is devised.The correction algorithm takes the errors between the estimation of the hybrid model and the laboratory testing data as the correction drive and calculates a suitable compensation value to append to hybrid model prediction. The correction strategy is as follows:(58)E(k)=λ×E0(k)+(1-λ)×E(k-1),E0(k)=ylab(k)-ymod⁡(k),ycor(k)=ymod⁡(k)+E(k),
where E(k) is the modified value for errors between the hybrid model and the laboratory value. E0(k) and E(k-1) are the error of the current and one-step previous hybrid model prediction and the testing data, respectively. λ is the weight constant of the two errors. ylab is the laboratory testing value of concentrate grade, ymod⁡ is the model prediction for concentrate grade, and ycor is the final output of hybrid model compensated by the corrected errors.Through the correction compensation strategy proposed above, the output of the hybrid model can be modified online, and as a result the accuracy of the model prediction can be improved.The flow chart of overall modified scheme is shown in Figure12. The variable p refers to the fixed time interval when the model structure should be updated. Variable k represents the kth sampling data.Figure 12
Flow chart of overall modified scheme.


## 3.4.1. Sliding Window Update Strategy

A sliding window is a kind of update strategy where only the most “recent” elements remain active and the rest are discarded. There are two equally important types of sliding window update strategy. One is updating with fixed-size windows, by which a fixed amount of the recent elements being active is defined (e.g., items arrive one at a time, and only the most recentn items remain active for some fixed parameter n). The other is updating with timestamp-based windows, by which the validity of an element is defined by an additional parameter such as a timestamp (e.g., many items can arrive in “bursts” at a single step and where only items from the last t steps remain active again for some fixed parameter t).Fixed-size windows are important for applications where the arrival rate of the data is fixed (but still extremely fast), such as sensors or stock market measurements. Timestamp-based windows are important for applications with asynchronous data arrivals, such as networking or database applications.In this paper, the hybrid model is built on the image texture information, and the images are processed at a fixed rate; thus, the arrival rate of the data for the model is fixed, which means that the fixed-size windows are suitable for this work.Since the hybrid model is composed of two parts, mechanism model and BS-PLS model, they both need to be updated at the same time.The parameter that needs to be updated for mechanism model is the constantβ according to (19). Therefore, for mechanism model, when the window slides one step forward, recalculate the Rave, rave, Lsav, and vgav within the current window; then, get the renewed value of β. Thus far, the mechanism model update is completed.For BS-PLS model, as it is built upon the sampled data, the most important problem is how to deal with these data. According to Section3, the structure of the BS-PLS model is determined by the zj (which is the B-spline transform of xj) and the variance S. In order to consider the history data (i.e., the data in the previous window) adequately, the update strategy for the average and variance of the current window data is as follows:
(57)z-k+1=N-PNz-k+PNz~k+1,Sk+12=N-P-1N-1Sk2+PN-1(z~k+1-z-k+1)2,
where z-k, Sk are the average and variance of the previous window data, respectively, and z~k+1 is the average of the current window data. Then, the update strategy is applied to reidentify the BS-PLS structure.Since the two composing models are changed, their entropy weight should also be updated along with them, which means that every time the window slides, after the two composed models are updated, we should recalculate the entropy weight according to Section3.


## 3.4.2. Online Correction Compensation Strategy

In a real process operation, the operating conditions may fluctuate frequently. In order to enable the hybrid model to have the adaptive ability to overcome the system disturbances and the drifting of the operation point, a practical online correction algorithm is devised.The correction algorithm takes the errors between the estimation of the hybrid model and the laboratory testing data as the correction drive and calculates a suitable compensation value to append to hybrid model prediction. The correction strategy is as follows:(58)E(k)=λ×E0(k)+(1-λ)×E(k-1),E0(k)=ylab(k)-ymod⁡(k),ycor(k)=ymod⁡(k)+E(k),
where E(k) is the modified value for errors between the hybrid model and the laboratory value. E0(k) and E(k-1) are the error of the current and one-step previous hybrid model prediction and the testing data, respectively. λ is the weight constant of the two errors. ylab is the laboratory testing value of concentrate grade, ymod⁡ is the model prediction for concentrate grade, and ycor is the final output of hybrid model compensated by the corrected errors.Through the correction compensation strategy proposed above, the output of the hybrid model can be modified online, and as a result the accuracy of the model prediction can be improved.The flow chart of overall modified scheme is shown in Figure12. The variable p refers to the fixed time interval when the model structure should be updated. Variable k represents the kth sampling data.Figure 12
Flow chart of overall modified scheme.


## 4. Implementation and Results

To validate the feasibility, the proposed hybrid model is tested on a set of industrial data from the flotation industry, and the results are compared with the output of mechanism model and BS-PLS model, respectively.


### 4.1. Results of Mechanism Model

According to (28), to get the relationship between concentrate grade G and recovery R, we need to obtain the function relating the experimentally determined values of Gins and Gtail. A dataset of 200 groups of industrial data is used, 150 groups out of which are used to identify the unknown parameters of mechanism model, and the rest are used as testing data. The 150 training data groups are denoised by means of wavelet analysis. To verify the generalization ability of the model, the testing data groups only with the input image features denoised and the concentrate keep the original value.First of all, the statistical formula forGins and Gtail can be gained through fitting the 150 sets of training data:
(59)Gins=2.5Gtail+4.7.Assume that the feed grade is constant; here, according to the industrial data, we set it equal to 5.Then, the final formula can be deduced for the relationship betweenG and R according to (22), (28), and (59):
(60)dGdR=17.2G2-17.5G2R-23.5GR-G317.2GR-12.5GR2-23.5R2.
After obtaining the relationship between recovery and concentrate grade, to predict the concentrate grade, we need to estimate the froth recovery first. According to (19), to estimate the froth recovery, the gas velocity vg, average froth size rout, loading Lsurf, and parameter β are needed. Among them, parameter β is difficult to calculate; thus, we could estimate it by derivation from (19):
(61)β*=(0.115Rave0.885-Rave1-Lsavravevgav)2,
where Rave is the average of recovery for n sampled data groups. rave, Lsav, and vgav are the average value of rout, Lsurf, and vg for n sampled data groups (denoised by the wavelet transform), respectively.For the 50 groups of testing data, the recovery can be estimated with the image featuresrout, Lsurf, and vg as inputs:
(62)R=0.885((vg/rout)β*)1-Lsurf0.885((vg/rout)β*)1-Lsurf+0.115.Finally, based on the relationship between concentrate grade and recovery, predict the concentrate grade with the estimated recovery; results are shown in Figure13.Figure 13
Mechanism model prediction.According to the simulation result, we can see that the prediction accuracy of mechanism model is not satisfying enough; however, it is capable of tracking the variance tendency of assay value, which is the advantage of mechanism model.


### 4.2. Results of BS-PLS Model

Since the higher the input dimension is, the more complex the calculations will be, to reduce the input dimension, it is necessary to analyze the relevance between the 14 input data sets and the concentrate grade and select the most relevant variables as model inputs. The results are shown in Figure14.Figure 14
Relevance between image textures and the concentrate grade. 1. Bubble size; 2. velocity; 3. stability; 4. lightness; 5. loading; 6. energy; 7. entropy; 8. inertial; 9. tune; 10.R-means; 11. G-means; 12. B-means; 13. relative redness; 14. saturation.Considering that the relevance of entropy andR-means is relatively low, select the remaining 12 features as model inputs. Similar to mechanism modeling, 150 out of 200 data groups are used to identify the model parameters; the rest are used to validate the model prediction. The 150 training data groups are denoised by means of wavelet analysis. To verify the generalization ability of the model, the 50 testing data groups only with the input image features are denoised and the concentrate data keeps original value.Referring to the BS-PLS modeling procedure mentioned before, to identify the model parameters, firstly, we need to transform the 150 groups of model inputs with B-spline function. In this paper, we selected the number of subsectionsMj as 3; then, based on the specific procedures described in Section 3, we construct the BS-PLS concentrate grade predicting model.To predict the concentrate grade of the 50 groups of testing data, likewise, we need to transform the input variables with B-spline function first, withMj=3, then predict concentrate grade. The result is shown in Figure 15.Figure 15
BS-PLS prediction.Simulation results indicate that the prediction accuracy of the BS-PLS model is better than mechanism model. Since the structural parameters of BS-PLS are determined by training sample data, the predicting performance is not steady enough.


### 4.3. Results of Hybrid Model

Here, like the case of BS-PLS modeling, 150 out of 200 data groups are used to train the proposed model; the rest are used to validate the model prediction. The prediction performance of the hybrid model is shown in Figure16.Figure 16
Hybrid model prediction.Figure16 indicates that the hybrid model could efficiently tracks the assay value of concentrate grade, and its prediction accuracy is satisfying. It synthesizes the advantage of mechanism model and BS-PLS model.


### 4.4. Comparison of the Three Models

To compare the three models, mechanism model, BS-PLS model, and hybrid model, we can apply root mean square error (RMSE) as the criterion:(63)RMSE=1N∑i=1N(yi-y^i)2.The comparison results are shown in Table1. We can see that the hybrid model gives a much better performance considering the RMSE.Table 1
Prediction error of three models.
RMSE
Mechanism model
1.0578
BS-PLS model
0.8878
Hybrid model
0.7025


### 4.5. Results of Online Updated and Corrected Hybrid Model

As mentioned before, flotation is a complex industrial process, and an updating and modification strategy is needed to improve the prediction accuracy of the proposed hybrid model. The specific procedures are as follows.Step 1.
Set the length of windowN=150 and slide step size p=25; initialize the error of the hybrid model output E(0)=0, index k=0.Step 2.
Adding a new set of sampling data,k=k+1, and executing the online correction compensation strategy mentioned before, we get the modified output of hybrid model.Step 3.
If the indexk<p, repeat Step 2, until k=p; then, execute the slide window update strategy; that is, identify the structure of BS-PLS model, the coefficient β of mechanism model, and the entropy weight of the two models again; then, reset E(t)=0, k=0 and repeat Step 2. The result is shown in Figure 17.Figure 17
The prediction result of online updated and corrected hybrid model.Figure17 indicates that the online updated and corrected hybrid model can predict the concentrate grade well; the RMSE can reach 0.5702, which means that it will perform more stably and accurately for industrial application.


## 4.1. Results of Mechanism Model

According to (28), to get the relationship between concentrate grade G and recovery R, we need to obtain the function relating the experimentally determined values of Gins and Gtail. A dataset of 200 groups of industrial data is used, 150 groups out of which are used to identify the unknown parameters of mechanism model, and the rest are used as testing data. The 150 training data groups are denoised by means of wavelet analysis. To verify the generalization ability of the model, the testing data groups only with the input image features denoised and the concentrate keep the original value.First of all, the statistical formula forGins and Gtail can be gained through fitting the 150 sets of training data:
(59)Gins=2.5Gtail+4.7.Assume that the feed grade is constant; here, according to the industrial data, we set it equal to 5.Then, the final formula can be deduced for the relationship betweenG and R according to (22), (28), and (59):
(60)dGdR=17.2G2-17.5G2R-23.5GR-G317.2GR-12.5GR2-23.5R2.
After obtaining the relationship between recovery and concentrate grade, to predict the concentrate grade, we need to estimate the froth recovery first. According to (19), to estimate the froth recovery, the gas velocity vg, average froth size rout, loading Lsurf, and parameter β are needed. Among them, parameter β is difficult to calculate; thus, we could estimate it by derivation from (19):
(61)β*=(0.115Rave0.885-Rave1-Lsavravevgav)2,
where Rave is the average of recovery for n sampled data groups. rave, Lsav, and vgav are the average value of rout, Lsurf, and vg for n sampled data groups (denoised by the wavelet transform), respectively.For the 50 groups of testing data, the recovery can be estimated with the image featuresrout, Lsurf, and vg as inputs:
(62)R=0.885((vg/rout)β*)1-Lsurf0.885((vg/rout)β*)1-Lsurf+0.115.Finally, based on the relationship between concentrate grade and recovery, predict the concentrate grade with the estimated recovery; results are shown in Figure13.Figure 13
Mechanism model prediction.According to the simulation result, we can see that the prediction accuracy of mechanism model is not satisfying enough; however, it is capable of tracking the variance tendency of assay value, which is the advantage of mechanism model.


## 4.2. Results of BS-PLS Model

Since the higher the input dimension is, the more complex the calculations will be, to reduce the input dimension, it is necessary to analyze the relevance between the 14 input data sets and the concentrate grade and select the most relevant variables as model inputs. The results are shown in Figure14.Figure 14
Relevance between image textures and the concentrate grade. 1. Bubble size; 2. velocity; 3. stability; 4. lightness; 5. loading; 6. energy; 7. entropy; 8. inertial; 9. tune; 10.R-means; 11. G-means; 12. B-means; 13. relative redness; 14. saturation.Considering that the relevance of entropy andR-means is relatively low, select the remaining 12 features as model inputs. Similar to mechanism modeling, 150 out of 200 data groups are used to identify the model parameters; the rest are used to validate the model prediction. The 150 training data groups are denoised by means of wavelet analysis. To verify the generalization ability of the model, the 50 testing data groups only with the input image features are denoised and the concentrate data keeps original value.Referring to the BS-PLS modeling procedure mentioned before, to identify the model parameters, firstly, we need to transform the 150 groups of model inputs with B-spline function. In this paper, we selected the number of subsectionsMj as 3; then, based on the specific procedures described in Section 3, we construct the BS-PLS concentrate grade predicting model.To predict the concentrate grade of the 50 groups of testing data, likewise, we need to transform the input variables with B-spline function first, withMj=3, then predict concentrate grade. The result is shown in Figure 15.Figure 15
BS-PLS prediction.Simulation results indicate that the prediction accuracy of the BS-PLS model is better than mechanism model. Since the structural parameters of BS-PLS are determined by training sample data, the predicting performance is not steady enough.


## 4.3. Results of Hybrid Model

Here, like the case of BS-PLS modeling, 150 out of 200 data groups are used to train the proposed model; the rest are used to validate the model prediction. The prediction performance of the hybrid model is shown in Figure16.Figure 16
Hybrid model prediction.Figure16 indicates that the hybrid model could efficiently tracks the assay value of concentrate grade, and its prediction accuracy is satisfying. It synthesizes the advantage of mechanism model and BS-PLS model.


## 4.4. Comparison of the Three Models

To compare the three models, mechanism model, BS-PLS model, and hybrid model, we can apply root mean square error (RMSE) as the criterion:(63)RMSE=1N∑i=1N(yi-y^i)2.The comparison results are shown in Table1. We can see that the hybrid model gives a much better performance considering the RMSE.Table 1
Prediction error of three models.
RMSE
Mechanism model
1.0578
BS-PLS model
0.8878
Hybrid model
0.7025


## 4.5. Results of Online Updated and Corrected Hybrid Model

As mentioned before, flotation is a complex industrial process, and an updating and modification strategy is needed to improve the prediction accuracy of the proposed hybrid model. The specific procedures are as follows.Step 1.
Set the length of windowN=150 and slide step size p=25; initialize the error of the hybrid model output E(0)=0, index k=0.Step 2.
Adding a new set of sampling data,k=k+1, and executing the online correction compensation strategy mentioned before, we get the modified output of hybrid model.Step 3.
If the indexk<p, repeat Step 2, until k=p; then, execute the slide window update strategy; that is, identify the structure of BS-PLS model, the coefficient β of mechanism model, and the entropy weight of the two models again; then, reset E(t)=0, k=0 and repeat Step 2. The result is shown in Figure 17.Figure 17
The prediction result of online updated and corrected hybrid model.Figure17 indicates that the online updated and corrected hybrid model can predict the concentrate grade well; the RMSE can reach 0.5702, which means that it will perform more stably and accurately for industrial application.


## 5. Conclusion

In this paper, a practical model for flotation concentrate grade prediction is proposed based on hybrid intelligent method with a deliberate combination of mechanism analysis, image feature processing, and updating strategy. Through sufficiently utilizing the relationship between concentrate grade and recovery, a mechanism model is built to predict concentrate grade upon mechanical estimation of recovery. Since the flotation mechanism is complicated, the mechanism model constructed based on some hypothesis inevitably restricts the predicting performance. Since the mechanism model is merely able to reflect the varying tendency of concentrate grade with relatively low prediction accuracy, a hybrid model integrated mechanism model and BS-PLS model based on froth image features are proposed to improve the concentrate grade predicting performance. Information entropy is adopted to distribute the weight of each model. To modify the output of hybrid model, output deviation compensation strategy is used. Moreover, the sliding window is employed to update the hybrid model structure, aiming at improving the adaptability of it. The model is verified by practical data from industrial field and the result shows its effectiveness for production guidance. The proposed hybrid modeling method is a specially designed effort to overcome the complexity of relationships in flotation process and a referential approach for similar problems.



---

*Source: 401380-2014-06-30.xml*