Precision prediction of aquaculture water quality: a spatiotemporal model integrating optimized-LSTM and radial basis function neural networks

Study area and data acquisition

The study focuses on an irregularly shaped aquaculture water area located at coordinates 32^∘18′50″N, 120^∘35′31″E, in Changzhou City, China, as depicted in Figs. 1A and 1B. This area measures approximately 150 m in length, 45 m in width, and 2 m in depth. The spatial environment of this region is well-preserved, with no surrounding pollution sources. During the deployment of the buoys, they were strategically positioned at various locations. Each buoy was equipped with multiple sensors installed at the same depth, ensuring uniformity across the sensors on each individual buoy. However, the deployment depths of the sensors varied among different buoys, ranging from 0.3 m to 1.9 m. To provide a clearer illustration of this arrangement, we have included a schematic diagram, as shown in Figs. 1C and 1D. The sensors, housed within waterproof enclosures, are detailed in Table 1. The collected water quality data includes DO concentration, water temperature, pH, ammonia nitrogen level, and ORP, with measurements recorded at 10-min intervals.

Data preparation

For the missing data in sensor collection, the adjacent data collected under the same conditions are utilized to estimate the missing values, as demonstrated in Eq. (1).

(1) $$X_{k+i}=X_k+\frac{i\cdot \left(X_{k+j}-X_k\right)}{j},0<i<j$$where $X_{k+i}$ represents the missing sensor data at time $k+i$, while $X_k$ and $X_{k+j}$ denote the original sensor data at times $k$ and $k+j$, respectively.

Given that water quality data and meteorological data exhibit continuity and temporal sequence, any instance where the variation range of continuously collected data exceeds that of the preceding and succeeding moments is identified as an outlier. These outliers are subsequently processed horizontally utilizing the mean smoothing method, as illustrated in Eq. (2).

(2) $$X_k^{\mathrm{\prime }}=\frac{X_{k-1}^{\mathrm{\prime }}+X_{k+1}^{\mathrm{\prime }}}{2}$$where $X_k^{\mathrm{\prime }}$ denotes the outlier at time $k$, while $X_{k-1}^{\mathrm{\prime }}$ and $X_{k+1}^{\mathrm{\prime }}$ represent the valid neighboring values at the adjacent time points.

Upon partitioning the collected data into training and testing sets, normalization is subsequently applied, as demonstrated in Eq. (3).

(3) $$\tilde{X}=\frac{X-X_{min}}{X_{max}-X_{min}}$$where $\tilde{X}$ represents the normalized value, X denotes the original data, and $X_{max}$ and $X_{min}$ correspond to the maximum and minimum values of the original sequence, respectively.

The research collected data samples over a 32-day period, specifically from March 15 to April 15, 2024. The dataset was divided into training and testing sets with a ratio of 70:30. Table 2 presents a representative subset of the raw data collected during the monitoring period of 2024.

Model’s performance evaluation metrics

To comprehensively evaluate the predictive performance of the model, this study utilizes three key evaluation metrics: RMSE, MAE, and the Coefficient of Determination (R-squared, R²). Superior model is indicated by lower values of RMSE and MAE, along with a higher R² value. The mathematical formulas for these metrics are presented in Eqs. (4) to (6).

(4) $$RMSE=\sqrt{\frac{1}{n}\sum _{i=1}^n{\left(y_i-{\hat{y}}_i\right)}^2}$$

(5) $$MAE=\frac{1}{n}\sum _{i=1}^n\left|{\hat{y}}_i-y_i\right|$$

(6) $$R^2=1-\frac{{\sum }_{i=1}^n{\left(y_i-{\hat{y}}_i\right)}^2}{{\sum }_{i=1}^n{\left(y_i-\overline{y_i}\right)}^2}$$where $y_i$ represents the observed value of the water quality parameter, $\widehat{y_i}$ represents the predicted value of the water quality parameter, $n$ is the sample size, and $\overline{y_i}$ is the mean value of the water quality parameter.

Descriptions of the proposed models

To improve the spatiotemporal prediction accuracy of aquaculture water quality parameters, this study proposes an innovative prediction model that integrates ESSA-SA-LSM with RBF neural networks.

(1) SA-LSTM model

The primary objective of this study is to address the challenge of time series prediction for forecasting future water quality parameters based on historical data. To achieve this, an LSTM model is employed to construct a foundational prediction framework, which effectively captures the sequential temporal dependencies inherent in water quality parameters. Furthermore, the model is enhanced by incorporating a self-attention mechanism, leading to the proposed SA-LSTM prediction model. The architecture of the proposed model is depicted in Fig. 2.

The time step of the LSTM is represented by d, with a value of 6. The processed data from time t−d+1 to t, including pH, ammonia nitrogen, water temperature, DO, and ORP ( $X_{t-d+1}^i$ to $X_t^i$), are utilized as inputs to the SA-LSTM model. Initially, feature extraction is conducted using the LSTM network. Subsequently, the data are fed into the self-attention mechanism to capture the relationships between the current time step and historical water quality parameter data, assigning higher weights to more significant time steps. An improved SSA is employed to optimize the hyperparameters of the LSTM network. Finally, the model predicts water quality parameters, outputting the spatiotemporal predictions of pH, ammonia nitrogen, water temperature, DO, and ORP at time t+1 as $X_{t+1}^1$ to $X_{t+1}^5$.

The LSTM model is a specialized architecture of RNN that incorporates gating mechanisms and hidden states to handle the problems of gradient explosion and vanishing gradients when managing long-term dependencies. In contrast to traditional RNNs, LSTM introduces gating units that regulate the flow of information, allowing the model to retain significant data and essential features while discarding irrelevant information. The cell architecture of the LSTM is depicted in Fig. 3, and the computational processes for the forget gate ( $f_t$), input gate ( $i_t$), output gate ( $o_t$), cell state ( $c_t$) and output ( $h_t$) are presented in Eqs. (7) to (12).

(7) $$f_t=\sigma \left(W_f\left[h_{t-1};x_t\right]+b_f\right)$$

(8) $$i_t=\sigma \left(W_i\left[h_{t-1};x_t\right]+b_i\right)$$

(9) $$o_t=\sigma \left(W_o\left[h_{t-1};x_t\right]+b_o\right)$$

(10) $${\tilde{c}}_t=\tanh (W_c\left[h_{t-1};x_t\right]+b_c$$

(11) $$c_t=i_t\times {\tilde{c}}_t+f_t\times c_{t-1}$$

(12) $$h_t=o_t\times \tanh \left(c_t\right)$$where $\tilde{X}$ denotes the value labeled by time; tanh represents the hyperbolic tangent function, which outputs values in the range (−1, 1); $b_o$, $b_i$ and $b_c$ are constant biases; and $W_i$, $W_o$ and $W_c$ correspond to the weight matrices associated with the input gate, output gate, and the conveyor belt mechanism, respectively.

SA represents a significant advancement in attention mechanisms, offering two primary benefits: it facilitates the efficient extraction of critical information while minimizing attention to irrelevant data, and it concurrently reduces reliance on external information sources, thereby enhancing the capacity to capture intrinsic relationships within input data. By integrating the self-attention mechanism, the LSTM can effectively address challenges related to information overload, while simultaneously improving predictive accuracy and system robustness.

The attention mechanism is fundamentally based on three parameters: Q, K, and V. The output of LSTM, denoted as $h_t$, serves as the input for the SA module. The output of the SA module, denoted as $Attention(Q,K,V{)}_i$, is calculated as illustrated in Eq. (13).

(13) $$\{\begin{array}{c}Q={\omega }_Q{h}_t\hfill \\ V={\omega }_v{h}_t\hfill \\ K={\omega }_k{h}_t\hfill \\ Attention{(Q,K,V)}_i=\mathrm{softmax}\left(\frac{Q{K}^T}{\sqrt{d_k}}\right)v_i\hfill \end{array}.$$

The predicted value $y_t^i$ of the target water quality parameter is presented in Eq. (14).

(14) $$y_t^i=\mathrm{softmax}\left({\omega }_{\mathrm{o}\mathrm{u}\mathrm{t}}\cdot Attention(Q,K,V)+b_{\mathrm{o}\mathrm{u}\mathrm{t}}\right)$$where ${\omega }_{out}$ represents the output weight matrix of the fully connected network, and $b_{out}$ denotes the output bias vector across the entire network structure.

(2) ESSA module

The SSA, a swarm intelligence-based optimization technique, is widely recognized for its robust optimization capabilities and rapid convergence. However, the standard SSA implementation exhibits certain limitations that hinder its performance. Firstly, the random initialization of sparrow population positions at the algorithm’s inception results in limited population diversity, consequently yielding suboptimal target solutions that adversely affect the algorithm’s iterative performance and error rates. Secondly, the discoverers in SSA tend to exhibit excessive aggressiveness during food source exploration. Upon locating an optimal solution, other individuals rapidly converge towards it, thereby diminishing population diversity and increasing susceptibility to local optima. Furthermore, the underutilization of the finder’s positional information may cause the algorithm to overlook potential optimal regions, thereby missing valuable exploration opportunities. To address these limitations, this study proposes the ESSA module that incorporates two key improvements: (1) composite chaos mapping for population initialization and (2) an adaptive inertial weight factor. These modifications aim to optimize the hyperparameters of LSTM neural networks, effectively resolving issues related to slow parameter convergence and inadequate global search capabilities. Specifically, the improved Tent-Logistic-Cosine composite mapping enhances the initial positioning of sparrows, as detailed in Eqs. (15)–(17), while the adaptive inertia weight factor and finder position optimization are mathematically formulated in Eqs. (18)–(19).

(15) $$y_{i+1}=\{\begin{array}{c}\cos (\pi (2r{y}_i+4(1-r)y_i(1-y_i))-0.5)),y_i<0.5\hfill \\ \cos (\pi (2r(1-y_i))+4(1-r)y_i(1-y_i))-0.5)),y_i\ge 0.5\hfill \end{array}$$where $y_i$ represents the sequence value at the i-th iteration, and $r$ denotes the control parameter that governs the system’s behavior.

When $y_i\in [0,1]$ and $r\in [0,1]$, Eq. (13) is normalized to constrain its range within the interval [0, 1], as expressed by Eq. (16).

(16) $$z_i=\mathrm{a}\mathrm{b}\mathrm{s}\left(\omega y_i\right)$$where $\mathrm{abs}(\cdot )$ represents the absolute value function, while $\omega$ serves as the normalization constant.

The final initialization of sparrow positions is given by Eq. (17).

(17) $$X_{i,j}^1=l{b}_j+z_i\left(u{b}_j-l{b}_j\right)$$where $X_{i,j}^1$ represents the initial coordinate of the $i\text{-}th$ sparrow in the $j\text{-}th$ dimension, where $l{b}_j$ and $u{b}_j$ denote the lower and upper bounds of the $j\text{-}th$ dimension, respectively. The dimension j corresponds to the parameter to be optimized, with $j=4$ in the current implementation. Specifically, the parameter ranges are defined as follows: the learning rate is bounded within [0.001, 0.1], the iteration number is constrained to [10, 1,000], and the number of hidden layer neurons is limited to [1, 100].

(18) $${\omega }_t=0.5\left(1-{\left(\frac{2t}{{\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}}_{max}}-1\right)}^3\right)$$

(19) $$X_{i,j}^{t+1}=\{\begin{array}{c}{\omega }_t\cdot X_{i,j}^t\cdot \exp \left(\frac{-i}{\alpha \cdot {\mathrm{iter}}_{max}}\right),R_2<T_s\hfill \\ {\omega }_t\cdot \left(X_{i,j}^t+Q\cdot L\right),R_2\ge T_s\hfill \end{array}$$where ${\omega }_t$ represents the adaptive inertia weight factor, where ${\omega }_{max}$ is set to 0.9 and ${\omega }_{min}$ to 0.3. The variable $t$ denotes the current iteration index, and j indicates the dimensionality of the parameters to be optimized $X_{i,j}^t$ represents the position of the $i\text{-}th$ sparrow in the $j\text{-}th$ dimension at the $t\text{-}th$ iteration. The term L denotes a unit row vector, and $\alpha$ is a random number uniformly distributed in the interval [0, 1]. The maximum iteration count is denoted by ${\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}}_{max}$ and Q represents a random number following the standard normal distribution. The warning threshold $R_2$ satisfies $R_2\in [0,1]$, and the safety threshold $T_s$ is bounded within the range [0.5, 1].

The ESSA module is employed to optimize four critical hyperparameters of the SA-LSTM model: the learning rate, iteration count, neuron count in the first hidden layer, and neuron count in the second hidden layer. The optimization procedure consists of the following specific steps:

• (a)

  Initialization: Configure the population size and the proportion of discoverers, joiners, and alarmers within the SSA module. The population is initialized using the Tent-Logistic-Cosine composite mapping to ensure a more uniform distribution of the initial population.

• (b)

  Parameter Space Definition: Define the dimensionality of the parameter space and establish the search boundaries for each parameter, thereby determining the initial positions of the sparrow population.

• (c)

  Fitness Evaluation: Calculate the fitness value for each sparrow using the error function, identifying the optimal fitness value and its corresponding position.

• (d)

  Position Update: Incorporate the adaptive inertia factor to update the positions of discoverers, joiners, and alarmers through the improved dynamic step size mechanism.

• (e)

  Termination Check: Re-evaluate the fitness values of the sparrow population. If the error falls within the target threshold, terminate the iteration; otherwise, return to step (d) for further optimization.

The workflow of the ESSA module is depicted in Fig. 4.

(3) Spatial prediction models

The key water quality parameters exhibit distinct three-dimensional distribution characteristics in large surface ponds. Significant variations are observed both across different water layers within the same vertical plane and at different positions within the same water layer. Consequently, measurements from a single point cannot adequately represent the overall water quality dynamics of the entire pond. To better characterize the spatiotemporal distribution of these critical parameters, a three-dimensional spatial analysis method is employed, enabling comprehensive monitoring of water quality parameter distributions throughout the pond.

To enhance the interpolation accuracy of water quality parameters, this study utilizes the RBF for interpolating prediction results. The coordinates of the RBF centers are determined from the training data using the k-means clustering algorithm. For each RBF center, the width parameter $\sigma$ is calculated based on the average distance from this center to its $p$ nearest neighboring centers, expressed by Eq. (20). To mitigate the risk of overfitting, $L2$ regularization is applied when solving for the output layer weights.

(20) $${\sigma }_i=\frac{1}{p}\sum _{j=1}^p\Vert c_i-c_j\Vert$$where $p$ = 3 is a predefined hyperparameter that represents the number of the nearest neighboring centers considered; $c_i$ is the $i_{th}$ RBF center; $c_j$ is the $j_{th}$ nearest neighbor of center $c_i$; $\Vert c_i-c_j\Vert$ represents the Euclidean distance between the two centers.

The spatiotemporal prediction outcomes are demonstrated using DO as a representative example. Figure 5 illustrates the structural framework of the RBF-based spatiotemporal prediction model. The model incorporates data from the three sensors nearest to the interpolation point at time $t+1$ as input, generating predictions for both water temperature and DO values at the specified interpolation point and time.

This process begins with the construction of a RBF network and the initialization of its parameters. Subsequently, the time series prediction results at time $t+1$ from the three monitoring points nearest to the target prediction location are used as inputs to the RBF network. Specifically, 15 data points, represented as $X_{t+1}^{1,1}\sim X_{t+1}^{1,5}$, $X_{t+1}^{2,1}\sim X_{t+1}^{2,5}$ and $X_{t+1}^{3,1}\sim X_{t+1}^{3,5}$, are employed to predict the target point’s pH, ammonia nitrogen, water temperature, DO, and ORP values ( $y_{t+1}^1\sim y_{t+1}^5$) at time $t+1$.

Ablation experiments

To comprehensively evaluate the predictive performance of the proposed ESSA-SA-LSTM model, systematic ablation experiments were conducted using LSTM as the baseline model. The experimental design included three key components: the baseline LSTM, SA mechanism and ESSA optimization. The comparative results of these experiments are presented in Table 5.

The proposed model demonstrates significant performance improvements across all evaluation metrics compared to alternative architectures. Specifically, when compared to the SA-LSTM model, it achieves an average reduction of 48.52% in RMSE and 54.81% in MAE, while improving the $R^2$ value by an average of 3.46%. More substantial enhancements are observed relative to the LSTM model, with average reductions of 59.11% in RMSE and 57.49% in MAE, respectively, along with an average increase of 6.42% in $R^2$.

These performance gains can be attributed to three key innovations: (1) the implementation of correlation analysis on the original data, which effectively reduces computational redundancy; (2) the integration of a self-attention mechanism that mitigates information loss during sequential processing; (3) the incorporation of composite mapping and adaptive inertia factors within the enhanced SSA algorithm, which optimizes LSTM parameters, accelerates convergence, and enhances global search capabilities, ultimately improving overall model performance. The proposed model requires less than 10 MB of memory and achieves an average inference time of under 20 ms during training, both of which are entirely acceptable. In scenarios that demand even faster execution, the inference time can be significantly reduced by strategically compromising the density of the spatial prediction points.

Comparison of experimental results on time-series sequence

The model proposed in this study is benchmarked against the SA-LSTM, and standard LSTM models, all of which utilize identical datasets. Using Point 1 (Fig. 1D) as a case study, the comparative performance across five key variables is illustrated in Figs. 6A to 6E. The performance metrics reveal substantial improvements: compared to the SA-LSTM model, the proposed model achieves reductions of 75.00% in RMSE and 73.24% in MAE, along with a 3.37% increase in R². More significant enhancements are observed relative to the LSTM model, with 80.91% and 79.12% reductions in RMSE and MAE, respectively, and a 6.40% improvement in $R^2$. These results substantiate that the proposed model not only accurately tracks the temporal dynamics of dissolved oxygen and water temperature variations but also effectively predicts their absolute values. The comprehensive evaluation confirms the superior predictive capabilities of the proposed architecture over existing models.

The proposed model demonstrates a significant enhancement in performance across all five water quality parameters compared to the other two models. Significance testing reveals that all $p$-values are substantially below the 0.05 threshold, indicating highly statistically significant differences. These results further substantiate the superiority of the model presented in this study.

Spatiotemporal distribution prediction of water quality parameters

Using DO and water temperature as representative cases, the collected dataset reveals distinct diurnal patterns with observable peaks occurring at 12:00 and 24:00. This temporal variation underscores the necessity of developing accurate spatiotemporal distribution models for these critical water quality parameters at these specific time points. To address this, the study employs the ESSA-SA-LSTM model to predict DO and water temperature levels at 12:00 and 24:00 on April 15, 2024. These prediction results subsequently serve as input data for comprehensive spatiotemporal analysis. The spatial prediction component was implemented using the Python computational platform. For enhanced visualization and analysis of the spatiotemporal distribution patterns, the model incorporates RBF interpolation to generate detailed distribution maps. Figure 7 presents the spatial distribution characteristics of DO and water temperature at both 12:00 and 24:00, with subfigures 7A through 7D providing comprehensive visual representations of these patterns.

As illustrated in Figs. 7A and 7C, the spatiotemporal distribution of dissolved oxygen exhibits distinct vertical stratification, with DO concentration progressively decreasing with increasing depth from the water surface. This phenomenon can be attributed to the primary mechanism of oxygen production in aquatic systems-photosynthesis by phytoplankton and aquatic plants, which predominantly inhabit the surface and littoral zones. As depth increases, photosynthetic activity diminishes due to light attenuation caused by water refraction and turbidity, consequently reducing oxygen production and resulting in lower dissolved oxygen levels in deeper water layers. Similarly, Figs. 7B and 7D reveal a comparable vertical temperature gradient, with water temperature decreasing as distance from the surface increases. This thermal stratification pattern results from the extensive exposure of the surface layer to solar radiation, creating a temperature differential between the warmer surface waters and cooler deeper layers.

To comprehensively evaluate the performance of the proposed RBF method, we conducted comparative experiments using ESSA-SA-LSTM integrated with two alternative interpolation approaches: IDW and linear triangular interpolation. The evaluation was performed on identical datasets to ensure fair comparison. The study employed a cross-validation methodology to assess the accuracy of each algorithm. Table 6. presents the spatiotemporal prediction accuracy metrics for DO and water temperature across all evaluated algorithms. The results demonstrate that, compared to ESSA-SA-LSTM with linear triangular interpolation and ESSA-SA-LSTM with IDW, the proposed ESSA-SA-LSTM-RBF approach achieved reductions in MSE of 38.49% and 48.08%, respectively. Similarly, RMSE decreased by 35.15% and 27.96%, while MAE showed reductions of 15.09% and 3.68%, respectively. Through a comparative analysis of interpolation accuracy for DO and water temperature, the results indicate that the proposed algorithm exhibits superior performance, particularly in scenarios with limited monitoring points. This enhanced performance is attributed to the RBF method’s ability to effectively handle sparse data distributions while maintaining interpolation accuracy.