In this paper, the time series of dissolved oxygen in aquaculture water of marine ranches in Shandong Province, China, was used as a data sample to construct an accurate water quality prediction model. In recent years, Shandong Province has vigorously developed the construction of marine ranching, and more than 150 high-quality marine ranching projects have been built thus far. The authors of this paper deployed water quality sensors in 8 marine ranches on the Shandong Peninsula and established multiple aquaculture water quality data collection networks. Starting from June 2020, time series data, such as dissolved oxygen, pH, salinity, and water temperature, for aquaculture water quality indicators were collected year round, and the data were collected every ten minutes. Figure 1 shows the approximate locations of the marine ranches where water quality data were collected. Table 1 shows the descriptive statistics of the dissolved oxygen time series collected from 8 marine ranches from June 1, 2021 to July 1, 2022.

VMD is an adaptive, completely non recursive variational mode signal decomposition method that determines the frequency center and bandwidth of each component by iteratively searching for the optimal solution of the variational model. The nonstationary time series is decomposed into k relatively stationary intrinsic mode function (IMF) subsequences u k ( t ) with different frequency scales. The center frequency of

each subsequence is v k . The constraint condition is that the sum of the IMF components is equal to the original signal. The specific construction steps are as follows:

(1) Obtain the signal of each IMF component through the Hilbert transform, calculate its unilateral spectrum, and use exponential correction to modulate the spectrum of each mode function to the corresponding baseband [13], as shown in the following equation (1).

In the above formula, u k = { u 1 , u 2 ⋯ u k } represents the k IMF components obtained by the decomposition process, and v k = { v 1 , v 2 ⋯ v k } represents the center frequency of each component, δ t is a generalized function, that is, the Dirac delta function, j is the imaginary unit, and exponential e − j v k t correction is used to modulate the spectrum of each mode function to the corresponding fundamental frequency band (Niu et al., 2020).

(2) Calculate the square of the norm of the gradient of the demodulated signal and perform Gaussian smoothing to obtain the bandwidth of each mode signal, as shown in the following equation (2)

In the above formula, ‖ • 2 ‖ indicates 2-norm processing and * is the convolution operation.

(3) When construct a variational problem, the sum of the estimated bandwidths of each mode is the smallest, the constraint condition is that the sum of all modes is equal to the original signal f ( t ) , and the corresponding constraint variational expression is shown in equation (3).

(4) To find the optimal solution of the constrained variational problem, the augmented Lagrangian function is introduced by taking advantage of the quadratic penalty term and the Lagrangian multiplier to transform the constrained variational problem into an unconstrained variational problem, as shown in Equation (4).

In the above formula, α is the quadratic penalty term, which can ensure the accuracy of signal reconstruction in a Gaussian noise environment. λ is the Lagrangian multiplication operator, which can maintain the strictness of the constraints (Dragomiretskiy and Zosso, 2014; Yang and Liu, 2022).

(5) Use the alternating direction method of multipliers (ADMM) to iteratively solve each IMF component u k ( t ) and its corresponding center frequency v k . The formulas are (5) and (6), respectively.

where w is the frequency, f ^ ( w ) is the Fourier transform of the original signal f(t), and u ^ k ( w ) and λ ^ ( w ) are the Fourier transforms of u k ( t ) and λ ( t ) , respectively.

The SMA was proposed by Li et al. in 2020. The algorithm mainly simulates the behavior and morphological changes in slime moulds during foraging foraging (Li et al., 2020). The slime mould judges the location of the food according to the concentration of the food. The higher the concentration of food is, the stronger the propagating waves generated by the biological oscillator in the body. This causes an increase in the width of the veins, forming a network of veins with different thicknesses between multiple food sources and establishing the best path for finding food. In addition, after the slime mould obtains food, there is still a certain probability of searching the unknown area (Zubaidi et al., 2020). The position update formula of a slime mould individual is shown in (7).

Among them, U b and L b represent the upper and lower boundaries of the search range, R a is a random oscillation ranging from -a to a and gradually approaches 0 as the number of iterations increases, R b oscillates from -1to 1 and finally tends to 0. Moreover, r a n d represents a random value in the [0,1] interval, z is a custom parameter that represents the proportion of randomly distributed slime mould individuals to the population and is generally 0.03, t is the current iteration number, and X b ( t ) represents the slime mould individual with the best fitness. X ( t ) represents the position of the current iteration of the slime mould individual, and X A ( t ) and X B ( t ) are the positions of two random individuals. W represents the weight coefficient.

The update formulas of the control parameter p , parameter a and parameter R b in formula (7) are formulas (8), (9) and (10), respectively.

In the above formula, n is the number of slime mould populations, and S ( i ) represents the fitness value of the i-th slime mould individual. S b is the current optimal fitness value. The t a n h function is a nonlinear activation function whose output value is in the interval [-1,1].

In the above formula, t is the current number of iterations, and T max is the maximum number of iterations. The a r c t a n h ( x ) function is the reverse hyperbolic tangent function.

The value of R b oscillates in the interval [ − 1 , 1 ] , and the synergy between R b and R a simulates the selection behavior of slime moulds. The weight coefficient W of slime mould individuals is related to their fitness value. The weight formula of individuals whose fitness value sequence is in the top 50% in the slime mould population is provided in (11).

The weight coefficient formula W O of other slime mould individuals is provided in (12).

In Equations (11) and (12), r represents a random number ranging from 0 to 1, and F b and F w represent the best and worst fitness values for the current iteration, respectively.

A DBN is a deep neural network composed of multiple layers of restricted Boltzmann machines (RBMs) and one layer of a back propagation (BP) neural network (Hinton and Salakhutdinov, 2006; Hinton et al., 2006). The RBM consists of a visible layer and a hidden layer. In the DBN structure, the hidden layer of the bottom RBM is the visible layer of the next RBM, and the output of the last RBM is used as the input of the BP neural network. The training process of the DBN is carried out layer by layer, and the RBM of the current layer can be trained only after the RBM of the previous layer is fully trained. The structure of a simple DBN model is shown in Figure 2 .

The optimization of the DBN model consists of two steps: unsupervised training and supervised fine-tuning. The DBN trains each layer of RBM separately, determines the weights and offsets of the first RBM, and uses the states of its hidden neurons as the input vector for the second RBM. After the second RBM is fully trained, the second RBM is stacked on top of the first RBM, and each layer of RMB is trained in turn, making sure to extract as much feature information as possible.

The last layer of the DBN is the BP neural network. During supervised tuning training, the BP algorithm is used to update the weights and biases of the model. Each layer of RBM can only ensure that the weights and bias values in its own layer are optimal for the feature vector mapping of this layer and not for the feature vector mapping of the entire DBN. The BP algorithm propagates the error information from top to bottom to each layer of the RBM and fine-tunes the entire DBN layer by layer. RBM is a shallow neural network with only two layers. The training time of RBM is significantly reduced compared with the deep neural network. The overall training of the DBN deep neural network is simplified to the training of multiple RBMs, which can improve the training efficiency and speed up the convergence speed. After training, the network is fine-tuned by the BP algorithm, so that the model converges to the local optimum, improves the convergence accuracy, and realizes the efficient training of the DBN neural network. The DBN prevents falling into a local optimum through unsupervised layer-by-layer training and improves the convergence speed and accuracy through supervised fine-tuning.

The accuracy of the prediction model is measured by calculating the similarity between the predicted value and the actual value. Common evaluation indicators include mean absolute error (MAE), mean absolute percentage error (MAPE), mean squared error (MSE), Willmott indexr (WI)(Willmott et al., 2012) and goodness of fit ( R 2 ). In this paper, the above indicators are used to evaluate the prediction performance of our model and other comparative models. Assuming that the predicted value is y ^ = { y 1 , y 2 ⋯ y n } , the actual value is y = { y 1 , y 2 , ⋯ y n } , and n is the number of samples, the formulas of each indicator are as follows: (13), (14), (15), (16) and (17).

where y ¯ is the mean of all samples.

Aquaculture water quality data are easily affected by the external environment, biological activities and fishermen intervention. They have the characteristics of nonlinearity and instability and are affected by the accuracy and performance of water quality sensor acquisition equipment, and the data contain much noise. Therefore, training the model directly with the original dissolved oxygen time series limits the accuracy and stability of the prediction to a certain extent. Through the mode decomposition of VMD, multiple IMF subsequences with different frequency scales and relative stability are obtained. Using each subsequence for model training can greatly improve the prediction accuracy.

When applying VMD for mode decomposition, we can set the number of mode decomposition components according to the data characteristics. When the number of modes is too small, some important information in the original signal will be filtered, which will affect the accuracy of subsequent predictions. When the number of modes is too large, the center frequencies of adjacent mode components will be closer, resulting in mode overlap and excessive decomposition. In this paper, when determining the number of VMD modes, the MAE of the original signal and the decomposed mode reconstruction sequence is used as the optimization index, and the mode decomposition number is determined according to the change trend of the MAE value. The calculation of the MAE is shown in Equation (18).

In the above formula, f ( t ) is the original signal to be decomposed, u k ( t ) is the decomposed mode component, k is the number of mode components, and N is the length of the time series. The VMD of 8640 dissolved oxygen data from Weihai Demingmarine ranching from June 1, 2021, to August 1, 2021, is taken as an example to optimize the mode decomposition number. In this paper, the number of decomposition components is set in the range of 2 to 10. After each decomposition step, the decomposed IMF sequence is reconstructed, and the MAE value of the reconstructed sequence and the original sequence are calculated. From the variation trend of the MAE in Figure 3 , it can be found that when the number of decomposition modes is 4, the value of the MAE tends to be stable, so the number of VMD modes of this dissolved oxygen time series is set to 4.

Figure 4 shows the decomposition effect of this dissolved oxygen time series.

The SMA has a simple structure and strong global search ability, but it also has shortcomings, such as a weak shrinkage mechanism and the tendency to easily fall into a local optimum. To make the SMA have better convergence accuracy and more stability, this paper uses elite opposition-based learning and a nonlinear convergence factor to improve it and proposes a multi strategy ISMA.

The opposition-based learning (OBL) model proposed by Tizhoosh calculates the current solution and the opposite solution in the process of searching for a feasible solution to a problem and selects an item that is closer to the optimal solution for the next iteration (Rahnamayan et al., 2008). OBL has been applied to improve various intelligent algorithms (Gupta et al., 2020; Shekhawat and Saxena, 2020; Tubishat et al., 2020; Wang et al., 2021; Li et al., 2022a).

Elite OBL (EOBL) forms an elite population and its opposite population based on the individual with the best fitness and then selects an individual with better fitness from the two populations to form a new population (Abualigah et al., 2021).

Assuming that at the t-th iteration, X i j e ( t ) is the elite solution of individual X i j in the j dimension, then its opposite solution is X ¯ i j E ( t ) . Its solution formula is shown in equation (19).

Among them, λ is a random value in the range of 0 to 1, and l j and u j are the minimum and maximum values of the j-th dimension individual, respectively. If the elite opposition solution X ¯ i j E ( t ) exceeds the boundary, it is reset by random generation, as shown in formula (20).

The performance of the slime mould optimization algorithm is greatly affected by the quality of the initial population. A high-quality initial population can speed up the convergence of the algorithm and help to find the global optimal solution. This paper applies the EOBL strategy to the initialization of the SMA population, calculates the fitness value, obtains the elite slime mould individuals and dynamic boundaries, and improves the quality and diversity of the initial population.

In the iterative optimization process of SMA, the change in parameter a has an important influence on balancing the global search ability and local exploration ability. Parameter a decreases rapidly in the early iterations, this decrease slows down inthe later iterations, and a smaller a in the early stage is not conducive to global exploration. Therefore, this paper proposes a new nonlinear decreasing strategy, and the updated formula for parameter a is shown in (21).

Among them, t is the current number of iterations, and T is the maximum number of iterations.

As shown in Figure 5 , the parameter a in this paper decreases slowly in the early stage, which is beneficial to global exploration. In later iterations, the speed of convergence is accelerated, which is beneficial to the local search.

In the construction of the DBN prediction model, several hyperparameter settings are involved, including the number of hidden layers, the number of nodes in each layer, the training times and learning rate in the unsupervised pretraining stage and the momentum. If the hyperparameter optimization is carried out by means of experimental verification, the process is slow, and it is difficult to obtain the best combination of hyperparameters. Therefore, this paper uses the ISMA to optimize the hyperparameters to build the prediction model.

The construction steps of an aquaculture water quality prediction model (VMD-ISMA-DBN) based on VMD and the ISMA to optimize the DBN are as follows.

The structure diagram of the VMD-ISMA-DBN prediction model is shown in Figure 6 .

Using the dissolved oxygen data of the Laizhou marine ranch as the experimental sample of the VMD-ISMA-DBN model, the four IMF subsequences obtained by mode decomposition were trained and predicted. The first 80% of the data were used as the training set, and the last 20% of the data were used as the test set. Figure 7 shows the convergence curves for each IMF sequence trained using the SMA-DBN model and the ISMA-DBN model respectively

As shown in Figure 7 , in the iterative optimization process of the four IMF subsequences, compared with SMA-DBN, ISMA-DBN has faster convergence speed, higher convergence accuracy, and stronger ability to avoid falling into local optimum. Table 2 shows the optimal values of the hyperparameters of the DBN model optimized by the ISMA for each IMF sequence, including the learning rate of the pretraining (lr_p), the number of iterations for pretraining (epo_p), the learning rate for fine-tuning (lr_f), the number of iterations for fine-tuning (epo_f), the momentum, the batch size (B_S) and the number of hidden layers (nlayer). From the data in Table 2 , it can be seen that the prediction of the two high-frequency components of IMF3 and IMF4 is relatively difficult, and the network structure is more complex than the two low-frequency components of IMF1 and IMF2.

Figure 8 shows the fitting effect of each IMF subsequence predicted by the VMD-ISMA-DBN model proposed in this paper.

Figure 8 shows that after dividing the original time series into relatively stable IMF subsequences, the complexity of the series is reduced, and the prediction accuracy is greatly improved. The low-frequency components with strong regularity and large amplitudes have an accurate prediction effect. There are some errors in the prediction of the high-frequency signal, but due to its small amplitude, it has little effect on the final prediction after being superimposed with the prediction result of the low-frequency component.

The prediction results of each IMF component are superimposed to obtain the final prediction results of the original dissolved oxygen time series. As shown in Figure 9 , the predicted value after superposition fits well with the actual value, and the model in this paper shows very accurate prediction performance.

To verify the effectiveness of each component of VMD-ISMA-DBN in improving model performance, the following five prediction models were constructed for ablation research.

The above five models are compared with the VMD-ISMA-DBN model proposed in this paper, and the experimental results are shown in Figure 10 .

The evaluation index values of each experimental model are shown in Table 3 .

As shown in Figure 10 and Table 3 , the DBN model needs to adjust many hyperparameters, and the prediction error is large if adjusted purely by experience. The training process of the DBN model is optimized by the SMA, the setting of hyperparameters is more reasonable, and the prediction error is reduced. The two strategies of EOBL and the improved convergence factor are separately used to optimize the SMA. The ISMA further optimizes the hyperparameter combination of the DBN, and the prediction accuracy is further improved. The original dissolved oxygen time series is preprocessed by VMD, and the decomposed signal is smoother, which can reduce the difficulty of prediction and significantly improve the prediction performance. The prediction curve of the VMD-ISMA-DBN model proposed in this paper has a better fitting effect than the other methods in the maximum value and minimum value, and the prediction trend is closer to the actual curve, which illustrates the more accurate prediction effect.

To further verify the prediction performance of the VMD-ISMA-DBN model, it is compared with six other intelligent prediction models, namely, the autoregressive integrated moving average model (ARIMA) (Xuan et al. (2021)), random forest (RF) regressor (Tiyasha et al. (2021)), temporal convolutional network (TCN) (Li et al. (2022b)), extreme learning machine (ELM) (An et al. (2022)), gated recurrent unit (GRU) (Cao et al. (2020) and long short-term memory (LSTM) (Barzegar et al. (2020)) The main parameters of each model are shown in Table 4 .

Figure 11 shows the prediction effect of each comparative model and VMD-ISMA-DBN on a 12-hour dissolved oxygen time series. It shows that the VMD-ISMA-DBN model proposed in this paper has a better prediction effect on dissolved oxygen water quality data than the other comparison models, and the prediction curve is more in line with the actual time series.

The ARIMA model performs stationary preprocessing on nonstationary sequences by means of differencing, and it cannot capture nonlinear relationships well. For unstable and nonlinear water quality data, the hyperparameter optimization process is complicated, and the prediction accuracy is lower than that of other models. The RF regressor model integrates multiple decision trees, and overfitting may occur when it is used for the prediction of time series with considerable noise. Due to the complexity of its structure, it takes more time to train than other models. The ELM has the characteristics of few training parameters and a fast training speed. However, for aquaculture water quality data with noise interference, the prediction effect of the ELM model is unstable. The LSTM and GRU neural networks achieve good results in time series prediction and solve the long-term dependency problem of RNNs. The vanishing gradient problem may occur during the training process. The model structure of the LSTM and GRU models is relatively complex, and it is easy to consume much memory when the training sequence is long. The TCN prediction model has the advantages of a stable gradient, low training memory requirements and accepting variable length inputs. The TCN model achieves better performance than the LSTM and GRU models in the prediction of dissolved oxygen time series. The DBN has the advantages of easy expansion, fast training, and less convergence time. In this paper, the improved SMA is used to optimize the hyperparameters of the DBN, and VMD is used to decompose the dissolved oxygen time series. The model in this paper achieves higher prediction accuracy than the other comparative models.

To verify the generalization performance of the model in this paper, the dissolved oxygen data of 6 other marine ranches were selected for experimental verification, and different time spans were selected to study the short-term and long-term prediction performance of the model. The water quality data of all marine ranches were collected every 10 minutes, and approximately 140 pieces of effective dissolved oxygen data were collected every day. During the experiment, 80% of the data are selected as the training set, and 20% of the data are used as the test set. Table 5 shows the statistics of the experimental data.

Figure 12 shows that the model in this paper achieves good prediction performance for the dissolved oxygen time series from the six marine ranches, and the prediction curve fits the actual data curve well, indicating that the prediction model has strong generalization performance and stability. Observing the prediction effect of different time spans from 15 days to 1 year, the prediction error of the long-term series increases to a certain extent, and the fitting effect of the data mutation decreases. The reason for this phenomenon may be the accumulation of prediction errors. The model in this paper eliminates the influence of data non stationarity on the prediction effect by decomposing the water quality data by VMD and optimizes the parameters of the DBN model through the ISMA, which makes the hyperparameter settings reasonable and improves the prediction accuracy. The experimental results show that the VMD-ISMA-DBN model proposed in this paper is suitable for the prediction of aquaculture water quality time series.