# Development of Hybrid Methods for Prediction of Principal Mineral Resources

Authors: Maria Qurban; Xiang Zhang; Hafiza Mamona Nazir; Ijaz Hussain; Muhammad Faisal; Elsayed Elsherbini Elashkar; Jameel Ahmad Khader; Sadaf Shamshoddin Soudagar; Alaa Mohamd Shoukry; Fares Fawzi Al-Deek
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Accurate estimation of the mining process is vital for the optimal allocation of mineral resources. The development of any country is precisely connected with the management of mineral resources. Therefore, the forecasting of mineral resources contributed much to management, planning, and a maximum allocation of mineral resources. However, it is challenging because of its multiscale variability, nonlinearity, nonstationarity, and high irregularity. In this paper, we proposed two revised hybrid methods to address these issues to predict mineral resources. Our methods are based on denoising, decomposition, prediction, and ensemble principles that are applied to the production of mineral resource time-series data. The performance of the proposed methods is compared with the existing traditional one-stage model (without denoised and decomposition strategies) and two-stage hybrid models (based on denoised strategy), and three-stage hybrid models (with denoised and decomposition strategies). The performance of these methods is evaluated using mean relative error (MRE), mean absolute error (MAE), and mean square error (MSE) as evaluation measures for the production of four principle mineral resources of Pakistan. It is concluded that the proposed framework for the prediction of mineral resources indicated better performance as compared to other existing one-stage, two-stage, and three-stage models. Furthermore, the prediction accuracy of the revised hybrid model is improved by reducing the complexity of the production of mineral resource time-series data.


---

## Body

## 1. Introduction

In the development of the national economy of any country, there is a significant contribution of minerals. The development of effective policies is necessary for mining industries. The income obtained by the mineral sector can be optimized, and the industry can alleviate the shortage of metals through accurate prediction of the mining process, as minerals have a massive impact on the economy of Pakistan. Since the mineral production data are nonstationary and have multiscale stochastic attributes, the mining process is affected by climate change and other projects related to socioeconomic developments that result in a challenging task for predicting the production of minerals resources. Two models are commonly used to predict such types of data: a process-based model and a data-driven model. The process-based model needs an extensive calibration and validation dataset [1]. The data-driven model considered the physical mechanism and scientific knowledge of stochastic geologic processes [2]. Huang et al. [2] used the process-based models for the prediction of hydrological data that were also nonlinear and having complex characteristics. They concluded that the lack of scientific knowledge and unavailability of sufficient data make it difficult to predict. Data-driven models were efficiently used to overcome the drawback of process-based models [3].The data-driven models are further categorized as traditional statistical and machine learning (ML) models. The traditional statistical methods, i.e., autoregressive integrated moving average (ARIMA), only consider the stationary and linear data. The ARIMA model is successfully applied to predict the production of mineral resources [4]. However, the problem with this traditional model is that it needed data to be stationary. Furthermore, since the production data of mineral resources are nonstationary and have nonlinearity in them, they also contain complex time-varying characteristics. Therefore traditional statistical models are not enough to capture the nonlinearity and time-varying characteristics of the nonlinear data. However, ML techniques can be used to deal with the drawbacks and instability of the traditional statistical models, as ML methods are used to deal with nonstationary and nonlinear time-series data [5]. Therefore, ML methods also are suitable for capturing the nonstationary and nonlinear behavior of mineral production data.Commonly used ML techniques are artificial neural network (ANN) and support vector machine (SVM), which also suffered from the problem of parameter selection and overfitting, and do not consider time-varying characteristics of nonstationary time-series data. Moreover, the production data of mineral resources contain noises that forbid the researcher to predict them accurately. Therefore, hybrid techniques are needed to model the noise corrupted and time-varying characteristics of the production of mineral resources data. Bokde et al. [6] proposed a hybrid model (EMD-ANN) for the nonstationary wind speed prediction and concluded that its performance was highly satisfactory and robust for jumping samplings than ANN and ARIMA models.To deal with the limitations of simple existing models, some of the preprocessing strategies are utilized with different data-driven models to increase the prediction accuracy of a different kind of nonlinear and nonstationary time series. The hybris models are developed to get the time-varying characteristics through noise reduction. The different hybrid models with preprocessing strategies have already been applied for complex and nonlinear data, having time-varying characteristics [7]. The strategy of the hybrid technique is based on decomposition and denoising, prediction, and ensemble stages [1, 8, 9]. Several algorithms are developed, such as spectral analysis, wavelet analysis (WA), Fourier analysis, and empirical mode decomposition (EMD), to reduce these noises or stochastic volatiles from the data [10, 11]. Fourier analysis and spectral analysis are used for those kinds of data that are stationary or linear. However, EMD and WA are the most commonly used preprocessing algorithms for nonlinear or nonstationary data and provide better results. The algorithms of WA decompose the nonlinear and nonstationary data of mineral resources into multiscale components [12]. These components are used as inputs at the prediction stage, and then these predicted components are ensemble for final prediction. The authors in [5] proposed a hybrid model for forecasting streamflow by using wavelet transform to decompose streamflows and use ANN for forecasting purposes. They concluded that their hybrid model was efficient than simple existing models. In the present paper, EMD- and WA-based thresholds are used to reduce noises from the mineral production data. The WA is considered as a powerful tool for converting a signal into a stationary signal with specific effectiveness. In the literature, different hybrid models with wavelet decomposition are used for the prediction of different kinds of nonlinear and nonstationary time-series data [13]. Wu et al. [14] proposed a hybrid forecasting model in combination with the particle swarm optimization algorithm and wavelet neural network (PSO-WNN) to predict China’s natural gas consumption. They used the PSO algorithm to optimize the initial weights, and the parameters of wavelet are updated through dynamic learning to have improvement in forecasting precision and reduce fluctuation of the WNN. They concluded that the proposed model is superior as compared to ANN- and WNN-based models. Some improvements have been made in hybrid modes consisting WA decomposition to get more accurate results of the prediction of daily flows which is also nonlinear time-series data [13]. However, the performance of the WA depends upon the selection of the type of mother wavelet. Prior knowledge about the signal, which is to be analyzed, and prior knowledge about its frequency content are needed for a suitable choice of the mother wavelet. Wu et al. [15] proposed an EMD method to overcome the shortcomings of WA for scrutinizing the nonlinear data and nonstationary datasets. Complex time-series data can be decomposed into a small and finite number of IMFs by using EMD. The EMD strategy has the advantage of converting the nonstationary series into stationary series. There exist different studies in the literature that used EMD with different data-driven models such as EMD-artificial neural network (EMD-ANN), EMD-radial basis function (EMD-RBF), EMD-support vector machines (EMD-SVMs), EMD-relevant vector machine (EMD-RVM), and EMD-ARIMA, and these hybrid models improve the prediction accuracy [16–18]. The EMD is combined with the ANN in many past studies, especially in hydrology [19], and also a novel model based on EMD and deep learning is used by Mi et al. [20] to reduce the noises and extract the information of trend of the original data of wind speed. However, EMD suffers from the problem of mode mixing between the IMFs [13]. To overcome the mode mixing problem, the authors in [13] proposed a new technique of decomposition called ensemble empirical mode decomposition (EEMD), in which they used white Gauss noise. The EEMD method can separate the signals without inappropriate mode mixing. It uses white noise that helps to establish the dyadic reference frame on time-scale space. They concluded that removing the problem of mode mixing EEMD produced a set of IMFs that carry the full physical meanings. Many hybrid techniques based on EEMD are used for streamflow and wind speed prediction and in hydrology [21, 22]. Di et al. [1] proposed a four-stage hybrid model based on EEMD decomposition to bring improvement in the prediction accuracy by minimizing the noises. They concluded that EEMD in combination with data-driven models could improve the prediction accuracy compared to EMD-based hybrid models. Liu [23] proposed a hybrid model by combining EEMD with the grey SVM model (GSVM) to forecast high-speed rail passenger flow. They concluded that their two-stage hybrid model is more suitable for short-term predictions than other single and hybrid models.Ghumman et al. [24] proposed two hybrid models based on preprocessing techniques EEMD, called EEMD pattern sequence-based forecasting (EEMD-PSF) and EEMD difference pattern sequence-based forecasting (EEMD-DPSF) models. The first model decomposes the series into number of IMFs and improves the prediction performance using the PSF model, and the second model reduced the effects of seasonality, trend, and irregularity from the wind speed data. They compared their performances with simple PSF, ARIMA, and LSSVM models for wind speed and concluded that EEMD-PSF and EEMD-DPSF models outperformed than simple models. Jaitly and Hinton [25] presented a comprehensive review of decomposition-based methods to improve the wind forecasting accuracy. They used and discussed wavelet analysis, EMD, seasonal adjusted, intrinsic time-scale decomposition, variational mode decomposition, and Bernaola-Galvan algorithm for decomposition of wind speed time-series data as an alternative forecasting algorithm. Through comparative analysis of various decomposition-based models, they concluded that the decomposition-based model provided the accuracy in prediction as compared to other forecasting algorithms.Although EEMD proved helpful in solving mode mixing issues, it created one other problem: some residual noise during signal reconstruction. As a result, there may be different modes because of the various realizations of signal and noise. To solve this problem, Mehrsai et al. [26] proposed the complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) technique. A specific noise is added in their proposed model at each decomposition stage, and a residue is calculated to get each mode. The final complete results of decomposition are obtained numerically with the negligible error term.This paper proposed a revised hybrid model to decompose the time-varying characteristics of the production data of mineral resources built upon CEEMDAN based on a 48-step ahead of the direct recursive prediction strategy. The study focuses on the preprocessing techniques (denoising and decomposing) and their effect on the prediction of mineral resources. An important point to note from the past studies is that they used a one-step-ahead prediction strategy rather than a multistep-ahead prediction [27]. In one-step-ahead prediction, all the observations are used by the predictor to estimate the parameter for the time step, and we used 48-month-ahead prediction to predict the mineral resources. For improved and better performance of the revised hybrid model, the long-term behavior of mineral resource production is observed by considering the 48-step-ahead direct recursive prediction strategy.In the current study, we aimed to make a revised hybrid model using the preprocessing techniques to decompose the production data of mineral resources using the WA-CEEMDAN technique, which follows the steps of EEMD by adding white noise at each level of decomposition. Therefore, the purpose of using the CEEMDAN algorithm in the current study is to find an efficient way to decompose the time-series data of mineral resources, which increased its prediction accuracy. Furthermore, the current study explores its prediction performance using 48-step-ahead prediction, especially for nonlinear data, by considering this emerging hybrid modeling technique.Based on the details outlined above, this paper is focused on developing a modified model to predict the production of mineral resources, improving the EMD/WA-CEEMDAN-based multimodels (EMD/WA-CEEMDAN). In this paper, the prediction performance of CEEMDAN-based hybrid models used for long-term to review the prediction of mineral resources is observed. Section2 is focused on the motivations behind the production data of mineral resource prediction. A short review of models used for prediction of mineral’s data and an introduction to EMD, EEMD, CEEMDAN, and their modified versions are discussed in Section 3. Additionally, a short review of various approaches used in hybrid CEEMDAN models to select the appropriate prediction methods with the consideration of the characteristics of respective IMFs is also discussed in Section 3. Section 4 describes details of the study area and data. Section 5 presents and discusses the results of the case study, while conclusions are made for this research in Section 6.


## 2. Motivation behind the Prediction of the Mineral Resources

In recent years, prediction of mineral resource prediction has become a challenging task for researchers. Although there is abundance of coal reserves in Pakistan, still share of gas and oil is about 65% in the energy mix. Pakistan, despite being a mineral-enriched country, is facing an alarming situation as its power generation is based on foreign exchange. The mineral sector of Pakistan is dominated by four principal minerals which are gas, oil, gypsum, and coal. There is a need to analyze and accurate prediction of the production of these four major minerals to deal with emerging challenges. The main purpose towards accurate prediction of mineral production data is to get the efficient and optimum utilization of natural resources in the development of the economy.Prediction of the mineral resources is considered to be an important study because of nonlinear, nonstationary, and having time-varying characteristics of datasets. Many researchers are working towards accurate prediction of nonlinear datasets, having complex time-varying characteristics. Although these researchers have the same goal, most of them have varieties in their motivations.Saadat et al. [28] used the CEEMDAN algorithm in their proposed model to predict the daily peak load. The proposed three-stage hybrid model comprised the CEEMDAN technique, which showed a robust decomposed ability of reliable prediction. Nazir et al. [29] developed a robust hybrid model by utilizing the CEEMDAN algorithm to decompose hydrological time-series data, which showed the robust decomposition ability to forecast the nonlinear and complex time-series data. The purpose of including such article is that those inferred the practicability and serviceability of the hybrid models for the decomposition of nonlinear datasets.


## 3. Proposed Methodology

In the current study, we proposed two hybrid methods to enhance the prediction accuracy of the time-series data. Both methods have the same map and formation except in the denoising stage, where the approaches for eliminating the noises from the data are different. At the decomposition stage of both methods, CEEMDAN (an improved form of EEMD) [25] is used to identify the oscillation. At the prediction stage, the multimodels are used after considering the nature of IMFs to predict the obtained IMFs accurately. There are two main objectives of using the multimodels: predicting the IMF by looking up their nature accurately. The second objective is to determine the attainments of the complex and straightforward models after the reduction of the complexity of mineral resource time-series data by using decomposition strategies. Finally, all the IMFs obtained through prediction are combined to predict the time-series data. The proposed hybrid models consist of denoising, denoted as D-step, decomposition, and prediction, denoted as Decompose-step and P-step, respectively [29]. A short description of these steps is given as follows:(1)
In D-step, the WA and EMD methods are presented to remove the noise from the production of mineral resource time-series data.(2)
In Decompose-step, CEEMDAN is used to decompose the denoised series into thek IMF components and one residual term.(3)
In P-step, the series obtained after denoising is decomposed into IMF components and one residual term, which is predicted by using linear stochastic and nonlinear machine learning methods. The model, which has a lower rate of error for prediction than other models, is utilized for further consideration based on three performance evaluation criteria. At last, the predicted results are combined to get the final prediction.Compared with the earlier review papers, the main contribution of the current study is to review the prediction literature from the perspective of decomposition-based hybrid methods. More precisely, we compile the decomposition-based models and discuss their algorithms and structures. Furthermore, we examine that how decomposition-based hybrid models improve prediction accuracy. Finally, existing techniques adopted in recent years to improve the performance of decomposition-based models are also compiled and discussed.To better understand and ease the reader, by combining all the steps mentioned above, the proposed strategy can be named as EMD (denoising)-CEEMDAN (decomposing)-MM (multimodels), and WA (denoising)-CEEMDAN (decomposing)-MM (multimodels) shortly named as EMD-CEEMDA-MM and WA-CEEMDAN-MM models. Both models are exhibited in Figure1.Figure 1
The suggested structure of EMD/WA-CEEMDAN-MM for the prediction of production data of mineral resources.


### 3.1. Denoising Step (D-Step)

In many time-series data, noises are certain components that sometimes reduce the prediction accuracy. Many algorithms exist in the literature, such as spectral analysis, WA, Fourier analysis, and EMD, to reduce these noises or stochastic volatiles from the data [30]. Fourier analysis and spectral analysis are used for those kinds of data that are stationary or linear. However, EMD and WA can also deal with nonlinear or nonstationary data and provide better results. In this study, EMD- and WA-based thresholds are used to reduce noises from the mineral production data.


#### 3.1.1. Wavelet Analysis (WA)

WA has developed as an efficient tool for converting a signal into a stationary signal with specific effectiveness. The wavelet with finite energy (small wave) is established in time, or frequency for analysis of temporal phenomena serves as a basis function. One of the favorite uses of the discrete wavelet transformation (DWT) is to eliminate the noises from the signal. This idea contains a strategy of applying the threshold wavelet on the coefficients of the signal contains the noises called the wavelet technique for denoising. The wavelet threshold is recognized as a powerful method for removing the noises from the signal. For using the wavelet basis soft and hard threshold, it is necessary to decide about the choice of the mother wavelet function. Using the Symlet 8 mother wavelet [31], the mineral production data are decomposed into the following coefficients [26, 30]:(1)bk,l=∑l=02N−k−12−k/2ϕ2−kp−l,(2)ck,l=∑k=1K∑l=02N−k−12−k/2θ2−kp−l,where bk,l is the approximation coefficient and ck,l is the detail coefficient.Signal’s energy can be divided over some of the wavelet coefficients with high magnitude by using this technique. On the other side, noise energy can be divided into many wavelet coefficients that have a low magnitude. The denoised signal can be obtained by removing wavelet coefficient noise and by reconstructing the coefficients. Wavelet denoising has three steps: (1) decomposition of the signal into coefficients of wavelet, (2) applying thresholds on the coefficients obtained in step 1 [32], and (3) reconstruction of the wavelet coefficients.


#### 3.1.2. Hard and Soft Thresholds

The assignment of the actual threshold is necessary for the elimination of the noises during denoising. As in wavelet transformation, the efficiency of the denoising technique depends on the value of the threshold. Several selection procedures are most likely for deciding about the thresholding rule, but two symbolic ones are hard thresholding and soft thresholding [33]. In the current study, for calculating the wavelet coefficients, soft and hard thresholding are applied to implement the wavelet denoising technique, which is listed below [26]. The soft thresholding is(3)ck,l′=ck,l,ck,l≥Thk,0,ck,l<Thk,and the hard thresholding is(4)ck,l′=sgnck,lck,l−Thk,ck,l≥Thk,0,ck,l<Thk,where Thk is the threshold which is calculated as Thk=a2EklnN,k=1,2,3,…,K, a is the constant which considers the values in the interval of 0.4 and 1.4 with the jump of 0.1, and Md^k is the median deviation, i.e.,(5)Md^k=medianck−1,l|k=1,2,3,…,2k−1−10.6745.By utilizing the following equation, the decomposed signal is reconstructed using the approximations and noise-free details:(6)xp=∑k=02N−K−1ak,l′2−k/2ϕ2−kp−l+∑k=1k∑k=02N−k−1ck,l′2−k/2θ2−kp−l,where ak,l′ and ck,l′ are the threshold approximation and detailed coefficients, respectively.


#### 3.1.3. Empirical Mode Decomposition (EMD)

The EMD is a flexible kind of decomposition technique and a purely data-driven technique to decompose the sophisticated signals into the series of components of the nonlinear and nonstationary data. After decomposition, by combining all the components for reconstructing, the original signal without any loss of the information can be obtained. The primary purpose of EMD is to get the IMFs from the complex signals. The following two conditions should be satisfied by the extracted IMFs [34]: (a) from all the data, the total of zero-crossings and the extrema should either be equal or differ by at most one; (b) at any point, the value of the envelope should be zero.The main steps of the EMD for an original time seriesxp,p=1,2,3,…,N are as follows [33]:(1)
Identification of all the local extrema of original time seriesxp.(2)
Create the upper and lower envelope asUp and Gp by using a cubic spline.(3)
Estimate the mean value of the upper envelope and the lower envelopemp=Up+Gp/2.(4)
Find the difference in the mean of the envelope from the original seriesxp. The difference dp is calculated as dp=xp−mp. Then, examine the properties of dp.(5)
Repeat 1–4 steps before the number of extrema is less than or equal to one so that no more IMF can be extracted or residueep becomes a monotonic function.At last, the signal can be shown as the sum of all the IMFs and residueep, where n is the number of IMFs, Cip,i=1,2,3,…,n is the ith IMF, and epis the residue. The way to denoise the IMF is the same as described in steps (1)–(4), except for the last two because of the low-frequency IMFs, which are used without denoising them [26, 35]:(7)xp=∑indi+ep,where n shows the number of sifted IMFs, as i=1,2,3,…,nep is the trend of the signal and dip is the ith IMF. Except for the last two IMFs, which are used without applying denoising techniques on them because of the low magnitude in their frequencies, the rest of the procedure of denoising process is the same as in wavelet-based denoising by using (3), (4), and (7) equations. In equations (3), (4), and (7), according to the number of IMFs, the subscript is replaced by i. Before reconstructing the signal, a smooth signal for input can be obtained by applying the thresholds on IMFs. The rebuilding of the denoised signal is generalized as(8)xp=∑i=1n1−2hip+∑i=n1−2nhip+ep,where the parameter n1 is the number of IMFs which provides us with the easiness of the elimination of the low-order IMF which are noisy and also of higher-order IMFs, which are a little bit noisy in Gaussian noise conditions as i=1,2,3,…,nhip, and ep is the ith IMF and trend of the signal, respectively.


### 3.2. Decomposition Step (Decompose-Step)

In the decomposition step, we need to apply a helpful tool to reduce the noises, providing accurate results. EEMD technique is applied in the decomposition step to handle the trouble of mode mixing.


#### 3.2.1. Ensemble Empirical Mode Decomposition (EEMD)

To improve the EMD and mitigate the mode mixing, EEMD is developed. In this technique, the white noise added by EEMD is distributed equally among all the time-frequency space, which helps in the separation of the frequency scales and decreases the occurrence of mode mixing. The procedure is presented as follows [36, 37]:(i)
Initialization of the ensemble numberQ.(ii)
Set the amplitude of the added white noisei=1.(iii)
Add the random white noise signalwnip in the original signal xp:(9)xip=xp+wnip,wherewnip is the ith included series of white noise and xip denotes the ith included noise signal i=1∼Q,Q>1.(iv)
By using EMD, decompose the noise signalxip into N IMFsCj,ip,j=1,2,3,…,P, where Cj,ip shows the jth IMF of the ith noise signal and P is the total number of IMFs.(a)
To obtain the pro-IMFs1p, subtract the mean envelope m1p from the original x1p, i.e., s1p=x1p−m1p, where m1p=Up+Gp/2.(b)
Considers1p as a new signal if the average of the lower and upper envelope becomes zero and if the number of zero-crossing and the number of extrema are equal or almost one.(c)
Consider the points (a)–(c) of4th step as far as the resulting signal as proper IMF C1p.(d)
From the original signal,x1p, subtract the resulting IMF C1p. Consider the residue e1p as the new data and go back to step 1:(10)x1p=C1p+e1p.(e)
In the substep, if the residue becomes the monotonic function, then complete the algorithm, and ifi<Q, then go back to step 3. The last residual is treated as the trend:(v)
Estimate the ensemble meanIMFjp of all trials of each IMF:(12)IMFjp=1P∑j=1PCj,ip,j=1,2,3,…,P,i=1,2,3,…,Q.(vi)
Consider the meanIMFjp as the final mean of all the P IMFs.


#### 3.2.2. Complete Ensemble Empirical Mode Decomposition with Added Noise (CEEMDAN)

Although EEMD can bring down the problem of mode mixing to a certain degree with included white noise sequence, the error cannot be eliminated after the computation of the averaging to a finite number. It affects the sequence of reconstruction. For the elimination of the mode mixing, CEEMDAN adds the adaptive white noise smoothing pulse interference in decomposition, and for making the decomposition of the data more complete, it uses the properties of the mean Gaussian white noise whose mean is zero. The detailed procedure of the CEEMDAN is as follows.Persistent with EEMD, in the calculation of CEEMDAN,P times decompose the original signal xp, i.e., xp+riwnip, where ri is the parameter which deals the signal to noise ratio. The first component of the IMF is(13)IMF1p=1P∑i=1Pcj,1p.The residual of the signal is(14)e1p=xp−IMF1p.Thedp is defined as the lth IMF component obtained from EMD. The sequence v1p+r1d1pnjp is decomposed as follows to get the second IMF component:(15)IMF2p=1P∑j=1Pd1e1p+r1d1njp.The second residual signal is(16)e2p=e1p−IMF2p.Similarly, by following the above procedure, the expression of thelth residual signal will be as follows:(17)elp=el−1p−IMFlp.The expression of thel+1th residual signal is(18)IMFl+1p=1P∑j=1Pdlelp+rldlnjp.The above procedure is repeated until the criterion is met. The expression of the original sequence, if the number of IMF components isM, is as follows:(19)xp=∑i=1MIMFip+ep,where IMFip is the ith IMF, ep is the overall residual signal, and xp is the signal obtained after decomposition. The decomposed IMFs are further used in the prediction stage.


### 3.3. Prediction Step (P-Step)

Multistep-ahead prediction is used to predict a sequence of values in time-series data. This approach is applied to predict a model step-by-step and use the predicted value of the current time step to determine its value in the next time step. We split the data into training and testing datasets. For the training dataset, 70% of the observations are used, and 30% are used for the testing dataset.In the prediction stage, by using machine learning time series and stochastic methods, the denoised IMFs are used as input to predict the production data of mineral resources. For this, we used multistep-ahead forecasting strategy. The reason for using two different types of models in prediction is that the IMFs with high frequencies are predicted accurately through ML methods and do not provide accurate results for the IMFs having low frequencies. Stochastic models provide better outcomes for the prediction of the IMFs with low frequencies. These two types of models are used for the direct forecast of IMFs having high and low frequencies.The models used for the prediction purpose are briefly described as follows.


#### 3.3.1. Autoregressive Integrated Moving Average (ARIMA) Model

For the prediction of the IMF, the autoregressive moving average model is used as follows:(20)IMFki=ψ1IMFk−1i+⋯+ψpIMFk−pi+εki+φ1εk−1i+⋯+φqεk−qi,where IMFki shows the ith IMF and εpi shows ith residual obtained through CEEMDAN, and p and q are the lag values of autoregressive and moving average term. Sometimes the time-series data are nonstationary; in such a situation, the series can be made stationary by differencing to an appropriate degree. Then, in the case of differencing the series, the model is called ARIMAp,d,q, where d is a degree of difference that makes the series stationary.


#### 3.3.2. Group Method of Data Handling (GMDH)-Type NN Model

GMDH is a type of unexplored neural network. The GMDH-NN models are established by considering the evolutionary method of modeling (GEvoM), which is a program that generates a polynomial type neural network for modeling the data. The input variables, hidden layers containing neurons, best model structure, and the number of layers are determined automatically in these networks. By considering the evaluation criteria, some of the neurons are chosen then the output of these selected neurons turns into the input of the next layer. For the selection of the neurons, the prediction means square criterion is considered using some transfer functions shown in Table1.Table 1
Transfer functions for GMDH-NN algorithms.
Transfer functionsSigmoid functionv=1/1+exp−uTangent functionv=tanuRadial basis functionv=exp−u2Polynomial functionv=uThe procedure is repeated until the final layer. In the last segment, only one predicted neuron is considered. However, GMDH-NN selects the relation of only two variables and ignores the effect of an individual variable. The relationship between input and output variables is generally expressed through the Volterra functional series called Kolmogorov–Gabor polynomial [38]:(21)v=b0+∑k=1nbkuk+∑k=1n∑l=1nbklukul+∑k=1n∑l=1n∑m=1nbklmukulum+⋯.A refined form of GMDN-NN is the architecture group method of data handling (RGMDH-NN), which not only considers two variables but also considers them individually; the remaining procedure of RGMDH is the same as GMDH.


#### 3.3.3. Radial Basis Function Neural Network (RBFNN)

As an ANN technique, RBFNN is used to predict the decomposed IMFs and components of residual. The reason for the selection of RBFNN is the simplicity of its structure and flexibility in selecting the number of neurons [39]. On the contrary to other feedforward neural networks, that RBFNN consists of one layer. Moreover, the RBFNN has a sound capability of approximation and fast convergence speed [40].The structure of the RBFNN consists of two sections, one section is the nonlinear conversion of input (the first layer in Figure2) to the hidden layer (the second layer in Figure 2), and the other part is the linear conversion of the hidden layer to the output layer which is the third layer described in Figure 2. Some of the expressions of RBF are shown in Table 2.Figure 2
Structure of the radial basis function neural network.Table 2
Transfer functions for RBFN algorithms.
Radial basis functions for RBFN algorithmsPower functionψw=wc,c=odd,w∈RGaussian functionψw=exp−w2/2c2,c>0,w∈RSquare root functionψw=w2+c2,c>0,w∈RHyperbolic tangent functionψw=1−e−2w/1+e−2w,w∈RSigmoid functionψw=1+e−w−1,w∈RThin plate spline functionψw=w2logw,w∈RReciprocal square root functionψw=w2+c2−1/2,c>0,w∈R


#### 3.3.4. Mutistep-Ahead Prediction

Time-series prediction can be used for both single (one-step-ahead prediction) and multiple periods (multistep-ahead prediction). Multistep-ahead prediction has to deal with problems, such as accumulation error, uncertainty, and accuracy, unlike one-step-ahead prediction. However, accurate time-series prediction for long horizon has become challenging. A multistep-ahead time-series prediction consists of predicting the nextH values of a time series consists of N observations, where the forecasting horizon is denoted by H > 1.The DirREC strategy [41] combines the steps and the principles of direct and iterated strategies. The recursive strategy is used for the prediction of all IMFs using three different models, and at last, all the predicted IMFs are ensembled and predicted by direct approach. DirREC makes the prediction with different models for each horizon, and also, as the recursive approach, it increases the set of inputs by adding more variables related to the predictions of the previous step. However, the embedding size is different for all horizons.


## 3.1. Denoising Step (D-Step)

In many time-series data, noises are certain components that sometimes reduce the prediction accuracy. Many algorithms exist in the literature, such as spectral analysis, WA, Fourier analysis, and EMD, to reduce these noises or stochastic volatiles from the data [30]. Fourier analysis and spectral analysis are used for those kinds of data that are stationary or linear. However, EMD and WA can also deal with nonlinear or nonstationary data and provide better results. In this study, EMD- and WA-based thresholds are used to reduce noises from the mineral production data.


### 3.1.1. Wavelet Analysis (WA)

WA has developed as an efficient tool for converting a signal into a stationary signal with specific effectiveness. The wavelet with finite energy (small wave) is established in time, or frequency for analysis of temporal phenomena serves as a basis function. One of the favorite uses of the discrete wavelet transformation (DWT) is to eliminate the noises from the signal. This idea contains a strategy of applying the threshold wavelet on the coefficients of the signal contains the noises called the wavelet technique for denoising. The wavelet threshold is recognized as a powerful method for removing the noises from the signal. For using the wavelet basis soft and hard threshold, it is necessary to decide about the choice of the mother wavelet function. Using the Symlet 8 mother wavelet [31], the mineral production data are decomposed into the following coefficients [26, 30]:(1)bk,l=∑l=02N−k−12−k/2ϕ2−kp−l,(2)ck,l=∑k=1K∑l=02N−k−12−k/2θ2−kp−l,where bk,l is the approximation coefficient and ck,l is the detail coefficient.Signal’s energy can be divided over some of the wavelet coefficients with high magnitude by using this technique. On the other side, noise energy can be divided into many wavelet coefficients that have a low magnitude. The denoised signal can be obtained by removing wavelet coefficient noise and by reconstructing the coefficients. Wavelet denoising has three steps: (1) decomposition of the signal into coefficients of wavelet, (2) applying thresholds on the coefficients obtained in step 1 [32], and (3) reconstruction of the wavelet coefficients.


### 3.1.2. Hard and Soft Thresholds

The assignment of the actual threshold is necessary for the elimination of the noises during denoising. As in wavelet transformation, the efficiency of the denoising technique depends on the value of the threshold. Several selection procedures are most likely for deciding about the thresholding rule, but two symbolic ones are hard thresholding and soft thresholding [33]. In the current study, for calculating the wavelet coefficients, soft and hard thresholding are applied to implement the wavelet denoising technique, which is listed below [26]. The soft thresholding is(3)ck,l′=ck,l,ck,l≥Thk,0,ck,l<Thk,and the hard thresholding is(4)ck,l′=sgnck,lck,l−Thk,ck,l≥Thk,0,ck,l<Thk,where Thk is the threshold which is calculated as Thk=a2EklnN,k=1,2,3,…,K, a is the constant which considers the values in the interval of 0.4 and 1.4 with the jump of 0.1, and Md^k is the median deviation, i.e.,(5)Md^k=medianck−1,l|k=1,2,3,…,2k−1−10.6745.By utilizing the following equation, the decomposed signal is reconstructed using the approximations and noise-free details:(6)xp=∑k=02N−K−1ak,l′2−k/2ϕ2−kp−l+∑k=1k∑k=02N−k−1ck,l′2−k/2θ2−kp−l,where ak,l′ and ck,l′ are the threshold approximation and detailed coefficients, respectively.


### 3.1.3. Empirical Mode Decomposition (EMD)

The EMD is a flexible kind of decomposition technique and a purely data-driven technique to decompose the sophisticated signals into the series of components of the nonlinear and nonstationary data. After decomposition, by combining all the components for reconstructing, the original signal without any loss of the information can be obtained. The primary purpose of EMD is to get the IMFs from the complex signals. The following two conditions should be satisfied by the extracted IMFs [34]: (a) from all the data, the total of zero-crossings and the extrema should either be equal or differ by at most one; (b) at any point, the value of the envelope should be zero.The main steps of the EMD for an original time seriesxp,p=1,2,3,…,N are as follows [33]:(1)
Identification of all the local extrema of original time seriesxp.(2)
Create the upper and lower envelope asUp and Gp by using a cubic spline.(3)
Estimate the mean value of the upper envelope and the lower envelopemp=Up+Gp/2.(4)
Find the difference in the mean of the envelope from the original seriesxp. The difference dp is calculated as dp=xp−mp. Then, examine the properties of dp.(5)
Repeat 1–4 steps before the number of extrema is less than or equal to one so that no more IMF can be extracted or residueep becomes a monotonic function.At last, the signal can be shown as the sum of all the IMFs and residueep, where n is the number of IMFs, Cip,i=1,2,3,…,n is the ith IMF, and epis the residue. The way to denoise the IMF is the same as described in steps (1)–(4), except for the last two because of the low-frequency IMFs, which are used without denoising them [26, 35]:(7)xp=∑indi+ep,where n shows the number of sifted IMFs, as i=1,2,3,…,nep is the trend of the signal and dip is the ith IMF. Except for the last two IMFs, which are used without applying denoising techniques on them because of the low magnitude in their frequencies, the rest of the procedure of denoising process is the same as in wavelet-based denoising by using (3), (4), and (7) equations. In equations (3), (4), and (7), according to the number of IMFs, the subscript is replaced by i. Before reconstructing the signal, a smooth signal for input can be obtained by applying the thresholds on IMFs. The rebuilding of the denoised signal is generalized as(8)xp=∑i=1n1−2hip+∑i=n1−2nhip+ep,where the parameter n1 is the number of IMFs which provides us with the easiness of the elimination of the low-order IMF which are noisy and also of higher-order IMFs, which are a little bit noisy in Gaussian noise conditions as i=1,2,3,…,nhip, and ep is the ith IMF and trend of the signal, respectively.


## 3.1.1. Wavelet Analysis (WA)

WA has developed as an efficient tool for converting a signal into a stationary signal with specific effectiveness. The wavelet with finite energy (small wave) is established in time, or frequency for analysis of temporal phenomena serves as a basis function. One of the favorite uses of the discrete wavelet transformation (DWT) is to eliminate the noises from the signal. This idea contains a strategy of applying the threshold wavelet on the coefficients of the signal contains the noises called the wavelet technique for denoising. The wavelet threshold is recognized as a powerful method for removing the noises from the signal. For using the wavelet basis soft and hard threshold, it is necessary to decide about the choice of the mother wavelet function. Using the Symlet 8 mother wavelet [31], the mineral production data are decomposed into the following coefficients [26, 30]:(1)bk,l=∑l=02N−k−12−k/2ϕ2−kp−l,(2)ck,l=∑k=1K∑l=02N−k−12−k/2θ2−kp−l,where bk,l is the approximation coefficient and ck,l is the detail coefficient.Signal’s energy can be divided over some of the wavelet coefficients with high magnitude by using this technique. On the other side, noise energy can be divided into many wavelet coefficients that have a low magnitude. The denoised signal can be obtained by removing wavelet coefficient noise and by reconstructing the coefficients. Wavelet denoising has three steps: (1) decomposition of the signal into coefficients of wavelet, (2) applying thresholds on the coefficients obtained in step 1 [32], and (3) reconstruction of the wavelet coefficients.


## 3.1.2. Hard and Soft Thresholds

The assignment of the actual threshold is necessary for the elimination of the noises during denoising. As in wavelet transformation, the efficiency of the denoising technique depends on the value of the threshold. Several selection procedures are most likely for deciding about the thresholding rule, but two symbolic ones are hard thresholding and soft thresholding [33]. In the current study, for calculating the wavelet coefficients, soft and hard thresholding are applied to implement the wavelet denoising technique, which is listed below [26]. The soft thresholding is(3)ck,l′=ck,l,ck,l≥Thk,0,ck,l<Thk,and the hard thresholding is(4)ck,l′=sgnck,lck,l−Thk,ck,l≥Thk,0,ck,l<Thk,where Thk is the threshold which is calculated as Thk=a2EklnN,k=1,2,3,…,K, a is the constant which considers the values in the interval of 0.4 and 1.4 with the jump of 0.1, and Md^k is the median deviation, i.e.,(5)Md^k=medianck−1,l|k=1,2,3,…,2k−1−10.6745.By utilizing the following equation, the decomposed signal is reconstructed using the approximations and noise-free details:(6)xp=∑k=02N−K−1ak,l′2−k/2ϕ2−kp−l+∑k=1k∑k=02N−k−1ck,l′2−k/2θ2−kp−l,where ak,l′ and ck,l′ are the threshold approximation and detailed coefficients, respectively.


## 3.1.3. Empirical Mode Decomposition (EMD)

The EMD is a flexible kind of decomposition technique and a purely data-driven technique to decompose the sophisticated signals into the series of components of the nonlinear and nonstationary data. After decomposition, by combining all the components for reconstructing, the original signal without any loss of the information can be obtained. The primary purpose of EMD is to get the IMFs from the complex signals. The following two conditions should be satisfied by the extracted IMFs [34]: (a) from all the data, the total of zero-crossings and the extrema should either be equal or differ by at most one; (b) at any point, the value of the envelope should be zero.The main steps of the EMD for an original time seriesxp,p=1,2,3,…,N are as follows [33]:(1)
Identification of all the local extrema of original time seriesxp.(2)
Create the upper and lower envelope asUp and Gp by using a cubic spline.(3)
Estimate the mean value of the upper envelope and the lower envelopemp=Up+Gp/2.(4)
Find the difference in the mean of the envelope from the original seriesxp. The difference dp is calculated as dp=xp−mp. Then, examine the properties of dp.(5)
Repeat 1–4 steps before the number of extrema is less than or equal to one so that no more IMF can be extracted or residueep becomes a monotonic function.At last, the signal can be shown as the sum of all the IMFs and residueep, where n is the number of IMFs, Cip,i=1,2,3,…,n is the ith IMF, and epis the residue. The way to denoise the IMF is the same as described in steps (1)–(4), except for the last two because of the low-frequency IMFs, which are used without denoising them [26, 35]:(7)xp=∑indi+ep,where n shows the number of sifted IMFs, as i=1,2,3,…,nep is the trend of the signal and dip is the ith IMF. Except for the last two IMFs, which are used without applying denoising techniques on them because of the low magnitude in their frequencies, the rest of the procedure of denoising process is the same as in wavelet-based denoising by using (3), (4), and (7) equations. In equations (3), (4), and (7), according to the number of IMFs, the subscript is replaced by i. Before reconstructing the signal, a smooth signal for input can be obtained by applying the thresholds on IMFs. The rebuilding of the denoised signal is generalized as(8)xp=∑i=1n1−2hip+∑i=n1−2nhip+ep,where the parameter n1 is the number of IMFs which provides us with the easiness of the elimination of the low-order IMF which are noisy and also of higher-order IMFs, which are a little bit noisy in Gaussian noise conditions as i=1,2,3,…,nhip, and ep is the ith IMF and trend of the signal, respectively.


## 3.2. Decomposition Step (Decompose-Step)

In the decomposition step, we need to apply a helpful tool to reduce the noises, providing accurate results. EEMD technique is applied in the decomposition step to handle the trouble of mode mixing.


### 3.2.1. Ensemble Empirical Mode Decomposition (EEMD)

To improve the EMD and mitigate the mode mixing, EEMD is developed. In this technique, the white noise added by EEMD is distributed equally among all the time-frequency space, which helps in the separation of the frequency scales and decreases the occurrence of mode mixing. The procedure is presented as follows [36, 37]:(i)
Initialization of the ensemble numberQ.(ii)
Set the amplitude of the added white noisei=1.(iii)
Add the random white noise signalwnip in the original signal xp:(9)xip=xp+wnip,wherewnip is the ith included series of white noise and xip denotes the ith included noise signal i=1∼Q,Q>1.(iv)
By using EMD, decompose the noise signalxip into N IMFsCj,ip,j=1,2,3,…,P, where Cj,ip shows the jth IMF of the ith noise signal and P is the total number of IMFs.(a)
To obtain the pro-IMFs1p, subtract the mean envelope m1p from the original x1p, i.e., s1p=x1p−m1p, where m1p=Up+Gp/2.(b)
Considers1p as a new signal if the average of the lower and upper envelope becomes zero and if the number of zero-crossing and the number of extrema are equal or almost one.(c)
Consider the points (a)–(c) of4th step as far as the resulting signal as proper IMF C1p.(d)
From the original signal,x1p, subtract the resulting IMF C1p. Consider the residue e1p as the new data and go back to step 1:(10)x1p=C1p+e1p.(e)
In the substep, if the residue becomes the monotonic function, then complete the algorithm, and ifi<Q, then go back to step 3. The last residual is treated as the trend:(v)
Estimate the ensemble meanIMFjp of all trials of each IMF:(12)IMFjp=1P∑j=1PCj,ip,j=1,2,3,…,P,i=1,2,3,…,Q.(vi)
Consider the meanIMFjp as the final mean of all the P IMFs.


### 3.2.2. Complete Ensemble Empirical Mode Decomposition with Added Noise (CEEMDAN)

Although EEMD can bring down the problem of mode mixing to a certain degree with included white noise sequence, the error cannot be eliminated after the computation of the averaging to a finite number. It affects the sequence of reconstruction. For the elimination of the mode mixing, CEEMDAN adds the adaptive white noise smoothing pulse interference in decomposition, and for making the decomposition of the data more complete, it uses the properties of the mean Gaussian white noise whose mean is zero. The detailed procedure of the CEEMDAN is as follows.Persistent with EEMD, in the calculation of CEEMDAN,P times decompose the original signal xp, i.e., xp+riwnip, where ri is the parameter which deals the signal to noise ratio. The first component of the IMF is(13)IMF1p=1P∑i=1Pcj,1p.The residual of the signal is(14)e1p=xp−IMF1p.Thedp is defined as the lth IMF component obtained from EMD. The sequence v1p+r1d1pnjp is decomposed as follows to get the second IMF component:(15)IMF2p=1P∑j=1Pd1e1p+r1d1njp.The second residual signal is(16)e2p=e1p−IMF2p.Similarly, by following the above procedure, the expression of thelth residual signal will be as follows:(17)elp=el−1p−IMFlp.The expression of thel+1th residual signal is(18)IMFl+1p=1P∑j=1Pdlelp+rldlnjp.The above procedure is repeated until the criterion is met. The expression of the original sequence, if the number of IMF components isM, is as follows:(19)xp=∑i=1MIMFip+ep,where IMFip is the ith IMF, ep is the overall residual signal, and xp is the signal obtained after decomposition. The decomposed IMFs are further used in the prediction stage.


## 3.2.1. Ensemble Empirical Mode Decomposition (EEMD)

To improve the EMD and mitigate the mode mixing, EEMD is developed. In this technique, the white noise added by EEMD is distributed equally among all the time-frequency space, which helps in the separation of the frequency scales and decreases the occurrence of mode mixing. The procedure is presented as follows [36, 37]:(i)
Initialization of the ensemble numberQ.(ii)
Set the amplitude of the added white noisei=1.(iii)
Add the random white noise signalwnip in the original signal xp:(9)xip=xp+wnip,wherewnip is the ith included series of white noise and xip denotes the ith included noise signal i=1∼Q,Q>1.(iv)
By using EMD, decompose the noise signalxip into N IMFsCj,ip,j=1,2,3,…,P, where Cj,ip shows the jth IMF of the ith noise signal and P is the total number of IMFs.(a)
To obtain the pro-IMFs1p, subtract the mean envelope m1p from the original x1p, i.e., s1p=x1p−m1p, where m1p=Up+Gp/2.(b)
Considers1p as a new signal if the average of the lower and upper envelope becomes zero and if the number of zero-crossing and the number of extrema are equal or almost one.(c)
Consider the points (a)–(c) of4th step as far as the resulting signal as proper IMF C1p.(d)
From the original signal,x1p, subtract the resulting IMF C1p. Consider the residue e1p as the new data and go back to step 1:(10)x1p=C1p+e1p.(e)
In the substep, if the residue becomes the monotonic function, then complete the algorithm, and ifi<Q, then go back to step 3. The last residual is treated as the trend:(v)
Estimate the ensemble meanIMFjp of all trials of each IMF:(12)IMFjp=1P∑j=1PCj,ip,j=1,2,3,…,P,i=1,2,3,…,Q.(vi)
Consider the meanIMFjp as the final mean of all the P IMFs.


## 3.2.2. Complete Ensemble Empirical Mode Decomposition with Added Noise (CEEMDAN)

Although EEMD can bring down the problem of mode mixing to a certain degree with included white noise sequence, the error cannot be eliminated after the computation of the averaging to a finite number. It affects the sequence of reconstruction. For the elimination of the mode mixing, CEEMDAN adds the adaptive white noise smoothing pulse interference in decomposition, and for making the decomposition of the data more complete, it uses the properties of the mean Gaussian white noise whose mean is zero. The detailed procedure of the CEEMDAN is as follows.Persistent with EEMD, in the calculation of CEEMDAN,P times decompose the original signal xp, i.e., xp+riwnip, where ri is the parameter which deals the signal to noise ratio. The first component of the IMF is(13)IMF1p=1P∑i=1Pcj,1p.The residual of the signal is(14)e1p=xp−IMF1p.Thedp is defined as the lth IMF component obtained from EMD. The sequence v1p+r1d1pnjp is decomposed as follows to get the second IMF component:(15)IMF2p=1P∑j=1Pd1e1p+r1d1njp.The second residual signal is(16)e2p=e1p−IMF2p.Similarly, by following the above procedure, the expression of thelth residual signal will be as follows:(17)elp=el−1p−IMFlp.The expression of thel+1th residual signal is(18)IMFl+1p=1P∑j=1Pdlelp+rldlnjp.The above procedure is repeated until the criterion is met. The expression of the original sequence, if the number of IMF components isM, is as follows:(19)xp=∑i=1MIMFip+ep,where IMFip is the ith IMF, ep is the overall residual signal, and xp is the signal obtained after decomposition. The decomposed IMFs are further used in the prediction stage.


## 3.3. Prediction Step (P-Step)

Multistep-ahead prediction is used to predict a sequence of values in time-series data. This approach is applied to predict a model step-by-step and use the predicted value of the current time step to determine its value in the next time step. We split the data into training and testing datasets. For the training dataset, 70% of the observations are used, and 30% are used for the testing dataset.In the prediction stage, by using machine learning time series and stochastic methods, the denoised IMFs are used as input to predict the production data of mineral resources. For this, we used multistep-ahead forecasting strategy. The reason for using two different types of models in prediction is that the IMFs with high frequencies are predicted accurately through ML methods and do not provide accurate results for the IMFs having low frequencies. Stochastic models provide better outcomes for the prediction of the IMFs with low frequencies. These two types of models are used for the direct forecast of IMFs having high and low frequencies.The models used for the prediction purpose are briefly described as follows.


### 3.3.1. Autoregressive Integrated Moving Average (ARIMA) Model

For the prediction of the IMF, the autoregressive moving average model is used as follows:(20)IMFki=ψ1IMFk−1i+⋯+ψpIMFk−pi+εki+φ1εk−1i+⋯+φqεk−qi,where IMFki shows the ith IMF and εpi shows ith residual obtained through CEEMDAN, and p and q are the lag values of autoregressive and moving average term. Sometimes the time-series data are nonstationary; in such a situation, the series can be made stationary by differencing to an appropriate degree. Then, in the case of differencing the series, the model is called ARIMAp,d,q, where d is a degree of difference that makes the series stationary.


### 3.3.2. Group Method of Data Handling (GMDH)-Type NN Model

GMDH is a type of unexplored neural network. The GMDH-NN models are established by considering the evolutionary method of modeling (GEvoM), which is a program that generates a polynomial type neural network for modeling the data. The input variables, hidden layers containing neurons, best model structure, and the number of layers are determined automatically in these networks. By considering the evaluation criteria, some of the neurons are chosen then the output of these selected neurons turns into the input of the next layer. For the selection of the neurons, the prediction means square criterion is considered using some transfer functions shown in Table1.Table 1
Transfer functions for GMDH-NN algorithms.
Transfer functionsSigmoid functionv=1/1+exp−uTangent functionv=tanuRadial basis functionv=exp−u2Polynomial functionv=uThe procedure is repeated until the final layer. In the last segment, only one predicted neuron is considered. However, GMDH-NN selects the relation of only two variables and ignores the effect of an individual variable. The relationship between input and output variables is generally expressed through the Volterra functional series called Kolmogorov–Gabor polynomial [38]:(21)v=b0+∑k=1nbkuk+∑k=1n∑l=1nbklukul+∑k=1n∑l=1n∑m=1nbklmukulum+⋯.A refined form of GMDN-NN is the architecture group method of data handling (RGMDH-NN), which not only considers two variables but also considers them individually; the remaining procedure of RGMDH is the same as GMDH.


### 3.3.3. Radial Basis Function Neural Network (RBFNN)

As an ANN technique, RBFNN is used to predict the decomposed IMFs and components of residual. The reason for the selection of RBFNN is the simplicity of its structure and flexibility in selecting the number of neurons [39]. On the contrary to other feedforward neural networks, that RBFNN consists of one layer. Moreover, the RBFNN has a sound capability of approximation and fast convergence speed [40].The structure of the RBFNN consists of two sections, one section is the nonlinear conversion of input (the first layer in Figure2) to the hidden layer (the second layer in Figure 2), and the other part is the linear conversion of the hidden layer to the output layer which is the third layer described in Figure 2. Some of the expressions of RBF are shown in Table 2.Figure 2
Structure of the radial basis function neural network.Table 2
Transfer functions for RBFN algorithms.
Radial basis functions for RBFN algorithmsPower functionψw=wc,c=odd,w∈RGaussian functionψw=exp−w2/2c2,c>0,w∈RSquare root functionψw=w2+c2,c>0,w∈RHyperbolic tangent functionψw=1−e−2w/1+e−2w,w∈RSigmoid functionψw=1+e−w−1,w∈RThin plate spline functionψw=w2logw,w∈RReciprocal square root functionψw=w2+c2−1/2,c>0,w∈R


### 3.3.4. Mutistep-Ahead Prediction

Time-series prediction can be used for both single (one-step-ahead prediction) and multiple periods (multistep-ahead prediction). Multistep-ahead prediction has to deal with problems, such as accumulation error, uncertainty, and accuracy, unlike one-step-ahead prediction. However, accurate time-series prediction for long horizon has become challenging. A multistep-ahead time-series prediction consists of predicting the nextH values of a time series consists of N observations, where the forecasting horizon is denoted by H > 1.The DirREC strategy [41] combines the steps and the principles of direct and iterated strategies. The recursive strategy is used for the prediction of all IMFs using three different models, and at last, all the predicted IMFs are ensembled and predicted by direct approach. DirREC makes the prediction with different models for each horizon, and also, as the recursive approach, it increases the set of inputs by adding more variables related to the predictions of the previous step. However, the embedding size is different for all horizons.


## 3.3.1. Autoregressive Integrated Moving Average (ARIMA) Model

For the prediction of the IMF, the autoregressive moving average model is used as follows:(20)IMFki=ψ1IMFk−1i+⋯+ψpIMFk−pi+εki+φ1εk−1i+⋯+φqεk−qi,where IMFki shows the ith IMF and εpi shows ith residual obtained through CEEMDAN, and p and q are the lag values of autoregressive and moving average term. Sometimes the time-series data are nonstationary; in such a situation, the series can be made stationary by differencing to an appropriate degree. Then, in the case of differencing the series, the model is called ARIMAp,d,q, where d is a degree of difference that makes the series stationary.


## 3.3.2. Group Method of Data Handling (GMDH)-Type NN Model

GMDH is a type of unexplored neural network. The GMDH-NN models are established by considering the evolutionary method of modeling (GEvoM), which is a program that generates a polynomial type neural network for modeling the data. The input variables, hidden layers containing neurons, best model structure, and the number of layers are determined automatically in these networks. By considering the evaluation criteria, some of the neurons are chosen then the output of these selected neurons turns into the input of the next layer. For the selection of the neurons, the prediction means square criterion is considered using some transfer functions shown in Table1.Table 1
Transfer functions for GMDH-NN algorithms.
Transfer functionsSigmoid functionv=1/1+exp−uTangent functionv=tanuRadial basis functionv=exp−u2Polynomial functionv=uThe procedure is repeated until the final layer. In the last segment, only one predicted neuron is considered. However, GMDH-NN selects the relation of only two variables and ignores the effect of an individual variable. The relationship between input and output variables is generally expressed through the Volterra functional series called Kolmogorov–Gabor polynomial [38]:(21)v=b0+∑k=1nbkuk+∑k=1n∑l=1nbklukul+∑k=1n∑l=1n∑m=1nbklmukulum+⋯.A refined form of GMDN-NN is the architecture group method of data handling (RGMDH-NN), which not only considers two variables but also considers them individually; the remaining procedure of RGMDH is the same as GMDH.


## 3.3.3. Radial Basis Function Neural Network (RBFNN)

As an ANN technique, RBFNN is used to predict the decomposed IMFs and components of residual. The reason for the selection of RBFNN is the simplicity of its structure and flexibility in selecting the number of neurons [39]. On the contrary to other feedforward neural networks, that RBFNN consists of one layer. Moreover, the RBFNN has a sound capability of approximation and fast convergence speed [40].The structure of the RBFNN consists of two sections, one section is the nonlinear conversion of input (the first layer in Figure2) to the hidden layer (the second layer in Figure 2), and the other part is the linear conversion of the hidden layer to the output layer which is the third layer described in Figure 2. Some of the expressions of RBF are shown in Table 2.Figure 2
Structure of the radial basis function neural network.Table 2
Transfer functions for RBFN algorithms.
Radial basis functions for RBFN algorithmsPower functionψw=wc,c=odd,w∈RGaussian functionψw=exp−w2/2c2,c>0,w∈RSquare root functionψw=w2+c2,c>0,w∈RHyperbolic tangent functionψw=1−e−2w/1+e−2w,w∈RSigmoid functionψw=1+e−w−1,w∈RThin plate spline functionψw=w2logw,w∈RReciprocal square root functionψw=w2+c2−1/2,c>0,w∈R


## 3.3.4. Mutistep-Ahead Prediction

Time-series prediction can be used for both single (one-step-ahead prediction) and multiple periods (multistep-ahead prediction). Multistep-ahead prediction has to deal with problems, such as accumulation error, uncertainty, and accuracy, unlike one-step-ahead prediction. However, accurate time-series prediction for long horizon has become challenging. A multistep-ahead time-series prediction consists of predicting the nextH values of a time series consists of N observations, where the forecasting horizon is denoted by H > 1.The DirREC strategy [41] combines the steps and the principles of direct and iterated strategies. The recursive strategy is used for the prediction of all IMFs using three different models, and at last, all the predicted IMFs are ensembled and predicted by direct approach. DirREC makes the prediction with different models for each horizon, and also, as the recursive approach, it increases the set of inputs by adding more variables related to the predictions of the previous step. However, the embedding size is different for all horizons.


## 4. Study Area and Experimental Design

### 4.1. Selection of Area for Study

Pakistan is blessed with numerous geological potential and has numerous reserves of minerals such as coal, copper, gold, and limestone, which are much useful for industrial development. However, we have not yet promoted growth and eliminated poverty in the country by utilizing our natural resources to the maximum level. The wealthiest province of Pakistan is Baluchistan, with approximately 80 to 85% minerals. The rest of the minerals, 10 to 15%, are present in KPK, Sindh, and Punjab. Despite Pakistan’s precious mineral resources and two continuous mineral policies, this sector contributes poorly to the country’s GDP. The reason may be insufficient assessments and monitoring, political instability, problems related to weather, shortage of foreign investments, and insecurity in the mineral-rich areas (https://www.pc.gov.pk/uploads/pub/FIRST_05_PAGES_STRATEGY_FOR_MINERAL_SECTOR_DEVELOPMENT_IN_PAKISTAN.pdf).


### 4.2. Description of Data

The observed data consist of the production of principle minerals resources of Pakistan, which are named natural gas, oil, coal, and gypsum measured in terms of metric tons. It consists of 168 monthly observations recorded from July 2005 to June 2019. The data are divided into the training dataset and testing data for observing the model performance. The dataset of data contains 118 observations from the month of July 2005 to April 2015, and the testing dataset contains 50 observations from May 2015 to June 2019. The training dataset consists of 80% observations of the observed series, and the testing dataset includes 20% observations of the observed series.


### 4.3. Comparison of the Proposed Hybrid Model with Other Models

Both suggested models are compared with other models used for prediction, with and without decomposed and denoised techniques. We named them 1-stage, 2-stage, and 3-stage models for our convenience, which we used for comparison purposes:(i)
1-stage model: models without having denoising and decomposition techniques are selected in this stage, i.e., ARIMA. We called them 1-stage as used in [29].(ii)
2-stage model: in 2-stage models, denoised techniques (EMD/WA) are selected for comparison, having noise removal capacity, i.e., EMD/WA-ARIMA, EMD/WA-RGMDH, and EMD/WA-RBFNN. For prediction purposes, different models are selected to compare the statistical model with the models based on artificial intelligence, i.e., RGMDH and RBFNN models. These 2-stage models are selected from [42] for comparison purpose.(iii)
3-stage model: in these models, both denoised and decomposed strategies are accessed; that is, EMD-EEMD-MM is selected from [15] for the purpose of comparison. Multiple models are selected under 3-stage models for the prediction by keeping the same strategy as in the proposed model. A direct 48-step-ahead forecasting strategy is used in the prediction step. We used a multistep-ahead forecasting methodology. Three methods for prediction are used: one traditional statistical model, i.e., ARIMA (p, d, and q), and two machine learning methods, i.e., GMDH and RBFNN.


### 4.4. Accuracy Measure Techniques

The performance of the model can be achieved by measuring the closeness of the predicted values and the observed values for the test dataset. By using three evaluation measures, i.e., signal-to-noise ratio (SNR), mean relative error (MRE), mean absolute error (MAE), mean square error (MSE), and mean absolute percentage error (MAPE) [43], the prediction accuracy of the selected and proposed models is obtained. The following are their mathematical expressions:(22)MRE=1N∑p=1ngp0−gppredgp0,MAE=1N∑p=1ngp0−gppred,MSE=1N∑p=1ngp0−gppred2,MAPE=gpm0−gpmpredgpm0∗100,where gppred and gp0 are the predicted and real data, respectively, and n is the data size; also, gpm0 is the mean of the observed data, and gpmpred is the mean of the predicted values. MRE, MAE, and MSE measures the departure between the original values and predicted values. The neurons of neural networks (GMDH and RGMDH) are selected according to their MSEs.


## 4.1. Selection of Area for Study

Pakistan is blessed with numerous geological potential and has numerous reserves of minerals such as coal, copper, gold, and limestone, which are much useful for industrial development. However, we have not yet promoted growth and eliminated poverty in the country by utilizing our natural resources to the maximum level. The wealthiest province of Pakistan is Baluchistan, with approximately 80 to 85% minerals. The rest of the minerals, 10 to 15%, are present in KPK, Sindh, and Punjab. Despite Pakistan’s precious mineral resources and two continuous mineral policies, this sector contributes poorly to the country’s GDP. The reason may be insufficient assessments and monitoring, political instability, problems related to weather, shortage of foreign investments, and insecurity in the mineral-rich areas (https://www.pc.gov.pk/uploads/pub/FIRST_05_PAGES_STRATEGY_FOR_MINERAL_SECTOR_DEVELOPMENT_IN_PAKISTAN.pdf).


## 4.2. Description of Data

The observed data consist of the production of principle minerals resources of Pakistan, which are named natural gas, oil, coal, and gypsum measured in terms of metric tons. It consists of 168 monthly observations recorded from July 2005 to June 2019. The data are divided into the training dataset and testing data for observing the model performance. The dataset of data contains 118 observations from the month of July 2005 to April 2015, and the testing dataset contains 50 observations from May 2015 to June 2019. The training dataset consists of 80% observations of the observed series, and the testing dataset includes 20% observations of the observed series.


## 4.3. Comparison of the Proposed Hybrid Model with Other Models

Both suggested models are compared with other models used for prediction, with and without decomposed and denoised techniques. We named them 1-stage, 2-stage, and 3-stage models for our convenience, which we used for comparison purposes:(i)
1-stage model: models without having denoising and decomposition techniques are selected in this stage, i.e., ARIMA. We called them 1-stage as used in [29].(ii)
2-stage model: in 2-stage models, denoised techniques (EMD/WA) are selected for comparison, having noise removal capacity, i.e., EMD/WA-ARIMA, EMD/WA-RGMDH, and EMD/WA-RBFNN. For prediction purposes, different models are selected to compare the statistical model with the models based on artificial intelligence, i.e., RGMDH and RBFNN models. These 2-stage models are selected from [42] for comparison purpose.(iii)
3-stage model: in these models, both denoised and decomposed strategies are accessed; that is, EMD-EEMD-MM is selected from [15] for the purpose of comparison. Multiple models are selected under 3-stage models for the prediction by keeping the same strategy as in the proposed model. A direct 48-step-ahead forecasting strategy is used in the prediction step. We used a multistep-ahead forecasting methodology. Three methods for prediction are used: one traditional statistical model, i.e., ARIMA (p, d, and q), and two machine learning methods, i.e., GMDH and RBFNN.


## 4.4. Accuracy Measure Techniques

The performance of the model can be achieved by measuring the closeness of the predicted values and the observed values for the test dataset. By using three evaluation measures, i.e., signal-to-noise ratio (SNR), mean relative error (MRE), mean absolute error (MAE), mean square error (MSE), and mean absolute percentage error (MAPE) [43], the prediction accuracy of the selected and proposed models is obtained. The following are their mathematical expressions:(22)MRE=1N∑p=1ngp0−gppredgp0,MAE=1N∑p=1ngp0−gppred,MSE=1N∑p=1ngp0−gppred2,MAPE=gpm0−gpmpredgpm0∗100,where gppred and gp0 are the predicted and real data, respectively, and n is the data size; also, gpm0 is the mean of the observed data, and gpmpred is the mean of the predicted values. MRE, MAE, and MSE measures the departure between the original values and predicted values. The neurons of neural networks (GMDH and RGMDH) are selected according to their MSEs.


## 5. Results and Discussion

In Figure1, step-by-step procedure of the proposed method is described. In which firstly, in D-step, the original series is decomposed using EMD or WA. In the second step (Decompose-step), the denoised series is decomposed using the CEEMDAN technique. In the third step of Figure 1, the IMFs which are obtained through CEEMDAN are predicted using ARIMA, RGMDH, and RBFNN models. At the final stage, all the predicted IMFs are ensembled :D-stage results: according to Figure1, the original series of data is decomposed in the first step by using two noise removal filters; the results of denoising are described as follows:Wavelet-based denoising: by following the ELT process shown in Figure1 and using equations (1) and (2), the approximations are calculated on the data of mineral resources of Pakistan, and soft and hard thresholding are used to remove the noises from the coefficients of mineral production time-series data. Rules of soft and hard thresholding are calculated by using equations (3) and (4), respectively. Then, based on the lower value of MSE, denoised series, hard thresholds are reconstructed for wavelet analysis.EMD-based denoising: for removing the noises from the production of mineral resource data using EMD, IMFs are calculated using equation (7). For denoising, these calculated IMFs except for the last two IMFs, soft and hard thresholding rules, are used as in WA by using equations (3) and (4), respectively. Since the last two IMFs are smooth and have low-frequency characteristics, that is why there is no need for denoising the last two IMFs. The denoised IMFs based on hard thresholding reconstruct noise-free mineral resource time-series data from equation (8). The WA- and EMD-based denoising are combined in Figure 3 for gas and oil. The statistical measures of the original series and denoised series of all four minerals, i.e., mean, standard deviation, MRE, MAE, and MSE for EMD and wavelet noise removal techniques, are shown in Table 3. According to these findings, it can be observed that EMD behaves better than WA.According to statistical measures, the results show that WA and EMD behave differently for all four minerals. The mean and standard deviation of the original series, denoised series by EMD, and denoised by WA are almost the same in all minerals. However, for gas and coal production, the standard deviation becomes less by using WA. According to other statistical measures, i.e., MRE, MAE, and MSE, EMD performs better than WA, as, in all minerals, these measures have lower values than WA. Therefore, it is concluded that EMD and WA performed equally to denoise the mineral resources in the long run. In the decomposition stage, EMD and WA, both denoising techniques are used separately as input to get those characteristics that change in terms of varying frequencies, i.e., high frequencies and low frequencies.Decompose-stage results: to get the local changing features for the time from denoised time-series data of mineral resources by WA/EMD, which are decomposed further into six IMF components and one residual term. The EMD and CEEMDAN decomposition methods are used for the extraction of IMFs from four minerals. The decomposition results of EMD-CEEMDAN techniques of gas and oil are presented in Figure4, where all four minerals are decomposed into six IMF components, and one residual term is presented. The drawn-out IMFs represent mineral production time-series data characteristics where the starting IMFs show the higher frequency. In contrast, the last three IMFs represent lower frequency and residual shown as a trend.The results of WA-CEEMDAN decomposition are shown in Figure5. The extracted IMFs represent the attributes of production of mineral resource time-series data, where at the beginning, some IMFs have a higher frequency and then slowly frequency reduces till sixth IMF and residual shown as the trend in Figures 4 and 5. The amplitude’s value of the white noise is chosen as 0.2, and a maximum number of the ensemble members (1000) are selected. Almost all IMFs and residuals for all minerals show identical characteristics for EMD-CEEMDAN and WA-CEEMDAN decomposition methods.P-step results: three methods are adopted to predict all extracted IMFs and for residuals to get precise results. For this purpose, one traditional and two other nonlinear methods are used, i.e., ARIMA (p, d, and q) as traditional statistical method and GMDH-NN and RBFNN as nonlinear methods, which are used for the prediction of IMFs and residual for all four-mineral production. The mineral resource data of four minerals are partitioned into 70% and 30% for the training dataset and testing dataset, respectively. The parameters of the model and its structure are estimated using 118 observations of mineral resources. Then, the suggested model and other models used for comparison purposes are tested in terms of their validity using 30% mineral resource data. After estimating multimodels for every IMF component and residual, the model with the least value of MRE, MAE, and MSE is considered as the most appropriate and selected for the prediction of each IMF. The findings for the training dataset of the suggested model and all other models in comparison to gas, oil, coal, and gypsum production are given in Table 4. The prediction results of suggested models EMD-CEEMDAN-MM and WA-CEEMDAN-MM illustrated the effectiveness of all four minerals with a minimum value of MRE, MAE, and MSE compared with 1-stage, 2-stage, and 3-stage evaluation models. However, the suggested model WA-CEEMDAN-MM acquires the lowest value of MSE than another proposed EMD-CEEMDAN-MM model. The model with the worst prediction is a 1-stage model, as shown in Table 4, with the maximum value of MRE, MAE, and MSE.Here, the 1-stage model, i.e., the ARIMA model, attains the maximum value of MSE without applying the techniques of denoising and decomposition on the mineral production time-series data. The predicted graph of the suggested model EMD-CEEMDAN-MM in comparison with 2-stage, i.e., EMD-ARIMA and EMD-RGMDH models, is presented in Figure6, and also, the predicted graph of WA-CEEMDAN-MM in comparison with 2-stage, i.e., WA-ARIMA and WA-RGMDH models, is presented in Figure 7. Here, in Table 4, the MSE of the RBFNN model is greater than all the MSEs of other models, so for comparison, we skip RBFNN in the graphical presentation of predicted results. Denoising, decomposition, and ensemble principles can be used to predict the mineral resources’ production. In Table 4, 2-stage, i.e., ARIMA and RGMDH models, performs better than the single-stage model, and from the 2-stage model, the WA-based model performs better than EMD-based models.However, it is concluded from Table4 that the suggested models perform well to predict the production of mineral resources by decreasing its intricacy and increasing the performance of prediction over 1-stage, 2-stage, and 3-stage models.The prediction errors for the testing dataset are presented in Table5. It can be concluded that our suggested model performed better than all other benchmark models, i.e., the error obtained through prediction using the testing dataset in comparison with the error of prediction using the training dataset. From Table 5, it can be examined that the performance of the WA-CEEMDAN-MM model with minimum values of MAD, MAPE, and MSE is better than all other models. The graphical presentation of suggested models by considering the testing dataset for only gas production is exhibited in Figure 8, based on the EMD technique. Figure 9 is based on the WA technique. From these figures, it can be observed that our suggested model performed better than all other existing models.Overall comparison of the proposed model with denoised and decomposed models: in general, from Tables4 and 5, it can be observed that removing the noises from the production data of mineral resources by using EMD and WA techniques provides better results as compared to single or model without denoising and decomposing. It can be observed that the value of MAE, MRE, and MSE from Tables 4 and 5 of all four minerals performs well for 2-stage models as compared to 1-stage models. It is noticed that the 1-stage ARIMA model predicted some of the IMFs having low frequencies, precisely but not for the IMFs having high frequencies as they contained more time-varying characteristics. Moreover, different statistical and machine learning models are used with these denoised series to predict mineral production and explore the performances of simple and complex models. These statistical and machine learning models can be seen from Tables 4 and 5, where RBFNN with EMD and WA performs worst than EMD- and WA-based ARIMA and RGMDH. This shows that one can use simple models for predicting mineral production compared to complex models such as RBFNN.It is observed that the two-stage model performed on average 6.7% better than the 1-stage model, and the 3-stage existing model performed 76.3% better than the 1-stage model on average and 87.3% on average better than the two-stage model.Our proposed model 3-stage EMD/WA-CEEMDAN-MM attained 29.4% on average less value of MSE compared to the 1-stage model, 49.2% less value compared to 2-stage, and 68.75% less value of MSE as compared to the existing 3-stage model, similarly, on average 10.87% less value of MRE as compared to 1-stage, 2.0% on average less value of MRE than the 2-stage model, and 75.5% on average than the 3-stage model. Also, on average, there is a 22.49%, 52.38%, and −46.63% decrease in MAE by using the proposed model as compared with 1-stage, 2-stage, and existing 3-stage models.However, the performance of 3-stage models, i.e., EMD–CEEMDAN-MM and WA-CEEMDAN-MM, is better than 1-stage and 2-stage models as it reduces the complexity of mineral production data in many ways by combining denoising and decomposing techniques. The integrated features of denoising and decomposing of 3-stage models enhance the prediction accuracy of mineral production, as shown in Tables4 and 5. Moreover, the WA-based hybrid models provided better accuracy of prediction as compared to EMD-based models.Figure 3
The denoised series for the production of gas and oil. The denoised series are obtained by EMD-threshold (in blue color) and wavelet threshold (in red color).Figure 4
The EMD-CEEMDAN decomposition of gas (a) and oil (b) production. The four series are decomposed into six IMF components and one residual.
(a)(b)Figure 5
The WA-CEEMDAN decomposition of gas (a) and oil (b) production. The four series are decomposed into six IMF components and one residual.
(a)(b)Figure 6
Prediction results of gas and oil production using suggested EMD-CEEMDAN-MM in comparison with EMD-ARIMA and EMD-RGMDH predicted models.
(a)(b)Figure 7
Prediction results of oil production using suggested WA-CEEMDAN-MM in association with WA-ARIMA and WA-RGMDH models.
(a)(b)Figure 8
The prediction of gas production using the hybrid model (EMD-CEEMDAN-MM) in comparison with EMD-ARIMA and EMD-RGMDH models.Figure 9
The prediction of gas production using the hybrid model (WA-CEEMDAN-MM) in comparison with WA-ARIMA and WA-RGMDH models.Table 3
Statistical measures of WA- and EMD-based denoised production of four-mineral time-series datasets.
Mineral resourcesModeμσMREMAEMSEGas productionOriginal series1223555477.39EMD122352.45452.54564.560.005849.80WA1223553757.3833258.050.0314792170Oil productionOriginal series2341219353450.2EMD2341179353338.3441.920.00266404.30WA2341219343459.162569.990.036181270027Gypsum productionOriginal series109677.656078.37EMD109673.256031.31178.880.0043509.52WA109675.954138.909578.650.10143169479Coal productionOriginal series297596.673584.54EMD297615.973382.25315.160.00144581.50WA297596.658222.6730897.950.121552598008Table 4
The evaluation of the prediction error of the suggested model (EMD-CEEMDAN-MM and WA-CEEMDAN-MM) in comparison with other existing models for gas, oil, coal, and gypsum production.
Mineral productionModel nameModelsMREMAEMSEGas production1-SARIMA5562.080.05346886722-SWA-ARIMA564.790.00500307.2WA-RGMDH1041.510.011621799WA-RBFN22919.260.19545273283EMD-ARIMA3750.390.0323240715EMD-RGMDH3753.870.0323340680EMD-RBFN22916.160.195569133243-SEMD-EEMD-MM2020.750.027766789WA-CEEMDAN-MM158.460.0038472.38EMD-CEEMDAN-MM1195.000.012336537Oil production1-SARIMA86924.620.04129330644622-SWA-ARIMA32680.160.011778236730WA-RGMDH24693.970.01981106280WA-RBFN441208.800.17313742347809EMD-ARIMA86826.350.0412907664665EMD-RGMDH84082.760.0411935015062EMD-RBFN445083.300.173187040000003-SEMD-EEMD-MM34349.940.012014443387WA-CEEMDAN-MM5676.630.0054755189EMD-CEEMDAN-MM26569.210.011213733733Coal production1-SARIMA47665.830.1935541834582-SWA-ARIMA4919.990.0263677954WA-RGMDH7587.040.03109980128WA-RBFN65501.760.206500657737EMD-ARIMA47592.770.193543414179EMD-RGMDH46715.480.193480224927EMD-RBFN74976.310.2584614746513-SEMD-EEMD-MM14999.790.06365101582WA-CEEMDAN-MM2384.120.0110215110EMD-CEEMDAN-MM13761.370.05328725693Gypsum production1-SARIMA16095.920.164945536742-SWA-ARIMA10218.460.09239738995WA-RGMDH9318.040.08217996433WA-RBFN44669.820.423399773324EMD-ARIMA16040.730.16492041414EMD-RGMDH16495.230.16501503414EMD-RBFN45538.210.4335493289503-SEMD-EEMD-MM6876.550.0877558497WA-CEEMDAN-MM2569.180.0225152289EMD-CEEMDAN-MM21073.190.22773453648Table 5
The evaluation of the prediction error of the suggested model (EMD-CEEMDAN-MM and WA-CEEMDAN-MM) in comparison with other models for all four minerals having the testing dataset.
Mineral productionModel nameModelsMREMAEMSEGas production1-SARIMA3253.460.03165057152-SWA-ARIMA3357.080.0314297822WA-RGMDH85.120.0010700.46WA-RBFN85178.710.707484485277EMD-ARIMA3243.620.0316313393EMD-RGMDH3564.160.0317443358EMD-RBFN85223.530.7075076980543-SEMD-EEMD-MM120692.201.0014573026239WA-CEEMDAN-MM55.780.004528.72EMD-CEEMDAN-MM1121.050.011835955Oil production1-SARIMA104425.900.04169126347272-SWA-ARIMA92291.300.0311294925595WA-RGMDH5925.510.0064583522WA-RBFN1919080.000.703802014000000EMD-ARIMA104445.200.0416895725118EMD-RGMDH100369.300.0415464333281EMD-RBFN1919040.000.7038146790000003-SEMD-EEMD-MM2722506.001.007416098000000WA-CEEMDAN-MM7888.740.00147703705EMD-CEEMDAN-MM34261.480.011730054305Coal production1-SARIMA52706.140.1739864559452-SWA-ARIMA34790.810.111834667964WA-RGMDH3399.370.0117296854WA-RBFN234617.900.6961083423522EMD-ARIMA52689.360.173984927942EMD-RGMDH47981.980.163413398302EMD-RBFN234083.700.68633182689863-SEMD-EEMD-MM332627.701.06111602246835WA-CEEMDAN-MM2207.270.019144544EMD-CEEMDAN-MM20414.310.0720414.31Gypsum production1-SARIMA26046.030.1411649472752-SWA-ARIMA29409.920.161482365241WA-RGMDH1007.130.001651773WA-RBFN135000.600.7019562632811EMD-ARIMA26005.650.1411618477EMD-RGMDH27240.400.141228486486EMD-RBFN134601.000.69201659773853-SEMD-EEMD-MM193557.501.0238076758501WA-CEEMDAN-MM1023.040.012727448EMD-CEEMDAN-MM11205.510.06213535838


## 6. Conclusion

For the optimal mineral supply and the purposes of mineral resources, an accurate prediction of mineral resources is necessary. Here, some of the data processing methods are utilized to increase the prediction accuracy of such stochastic type data by efficiently using decomposition techniques. By using the three strategies, denoising, decomposing, and ensemble, two hybrid models with 3-stages are suggested, EMD-CEEMDAN-MM and WA-CEEMDAN-MM, from which WA-CEEMDAN-MM performed better than the EMD-CEEMDAN-MM model for the decomposition of the nonstationary and nonlinear data of minerals. For evaluating the performance of both models, the production data of four minerals are used. In general, our suggested model performed well for all four minerals than the other 1-stage, 2-stage, and existing 3-stage models. For evaluation measures, three techniques are utilized, i.e., MAE, MSE, and MAPE. These 3-stage hybrid models can be used for the decomposition of any nonlinear and nonstationary data for prediction [44–47].



---

*Source: 6362660-2021-08-09.xml*