Dissolved oxygen prediction method for tidal river network based on timing and dynamic process guidance

Abstract

This invention discloses a dissolved oxygen (DO) prediction method for tidal river networks based on time-series and dynamic processes, belonging to the field of intelligent water environment prediction technology. The method acquires hydrological and meteorological observation sequences, water quality monitoring sequences, and dynamic field simulation results output from a three-dimensional hydrodynamic model of the target water area to calculate the DO deficit sequence. It constructs a dynamic transport time distribution feature library based on particle transport tracking and generates physical prior biases for future prediction moments. The logarithmically transformed physical prior biases are embedded into the attention scoring function of a time-series prediction network, trained using multi-source time-series data and dynamic features as input, and outputs predicted DO deficit values ​​for future moments. This invention can improve the prediction accuracy and extrapolation robustness of tidal river networks under extreme hydrological conditions, and can identify the dominant periods of physically transported oxygen-consuming substances and the active periods of local biochemical oxygen consumption, providing a reliable prediction and diagnostic tool for water quality management in complex water bodies like tidal river networks.

Description

Technical Field

[0001] This invention belongs to the field of intelligent prediction technology for water environment, and specifically relates to a method for predicting dissolved oxygen in tidal river networks based on time-series and dynamic processes. Background Technology

[0002] Dissolved oxygen (DO) is a core indicator for measuring the health of water bodies, and its concentration directly affects the survival and metabolic activities of aquatic organisms. In complex water systems like tidal river networks, excessively low dissolved oxygen concentrations can lead to serious problems such as organism asphyxiation, black and odorous water, and ecosystem degradation. Accurately predicting the changing trends of dissolved oxygen concentration is of significant practical importance for water quality early warning, ecological protection, and water resource management.

[0003] Currently, dissolved oxygen (DO) prediction methods mainly fall into two categories: numerical models based on physical mechanisms and machine learning models based on data-driven approaches. Numerical models based on physical mechanisms simulate pollutant migration, diffusion, and biochemical reactions by solving hydrodynamic governing equations and water quality transport equations. Three-dimensional baroclinic dynamic models can simulate the flow field structure and temperature-salinity distribution under multiple dynamic influences such as tides, runoff, and wind, and can be combined with water quality modules to predict DOL. However, these methods are highly dependent on biochemical reaction parameters, and model calibration under complex boundary conditions is labor-intensive and computationally expensive, making it difficult to meet the needs of high-frequency rolling prediction. In recent years, data-driven models, represented by Long Short-Term Memory (LSTM) networks, temporal convolutional networks, and attention mechanisms, have been widely used in DOL prediction. These methods can learn temporal correlations from historical monitoring data and have advantages such as high computational efficiency and strong fitting ability. However, purely data-driven models lack constraints on the physical processes of pollutant transport, making it difficult to accurately characterize the transport time delays of pollutants under complex dynamic influences such as tidal flow, runoff changes, and density circulation. When upstream pollutants affect target sites after being transported for days or even weeks, models struggle to establish current dissolved oxygen changes and the spatiotemporal correlation between these changes and previous pollution events in neighboring areas. Furthermore, typical hypoxic events are usually closely related to specific hydrological and meteorological conditions such as heavy rainfall and large runoff, and their frequency is low. The limited number of historical samples available for model training makes it difficult for machine learning models to fully learn the evolution of dissolved oxygen under extreme conditions. With insufficient training samples, purely data-driven models have weak generalization and extrapolation capabilities for typical hypoxic events and other extreme scenarios, leading to a decline in prediction accuracy.

[0004] Therefore, it is urgent to combine the ability of physical mechanisms to characterize transport time delays with the ability of data-driven models to fit complex nonlinear relationships. This would allow the introduction of hydrodynamic physical mechanisms that can reflect the transport patterns of pollutants into data-driven prediction models, providing the models with physical prior constraints that conform to the actual transport process, thereby improving the accuracy and reliability of dissolved oxygen prediction under complex dynamic conditions and insufficient sample size. Summary of the Invention

[0005] This invention provides a method for predicting dissolved oxygen in tidal river networks based on time-series and dynamic processes, aiming to solve the technical problems mentioned in the background art.

[0006] This invention provides a method for predicting dissolved oxygen in tidal river networks based on time-series and dynamic processes, comprising: Step S1: Obtain the hydrological and meteorological observation sequences, water quality monitoring sequences, and multi-parameter dynamic field simulation results output by the three-dimensional hydrodynamic model from the target station; Step S2: Calculate the theoretical saturated dissolved oxygen concentration and dissolved oxygen deficit sequence based on the water quality monitoring sequence and the hydrological and meteorological observation sequence; Step S3: Based on the three-dimensional velocity field tracking neutral tracer particles in the dynamic field simulation results, construct a dynamic transport time distribution feature library corresponding to tidal phase, tidal pattern, and runoff conditions; Step S4: For each future predicted time, determine the dynamic state based on hydrological and meteorological forecast data and the forecast results of the three-dimensional hydrodynamic model, and combine the transport time interval from each historical time to the future predicted time to obtain the conditional arrival probability from the dynamic transport time distribution feature library and generate the corresponding physical prior bias vector. Step S5: Embed the physical prior bias obtained by logarithmic transformation of the conditional arrival probability into the attention scoring function of the time series prediction network. Take the hydro-meteorological observation sequence, water quality monitoring sequence, dynamic field simulation results and hydro-meteorological forecast features and dynamic forecast features corresponding to the prediction period as input, and take the dissolved oxygen deficit sequence as the supervision target. Train the network parameters and the non-negative scaling coefficient used to adjust the influence intensity of the physical prior bias through backpropagation. Step S6: Input the hydrological and meteorological observation sequence, water quality monitoring sequence, dynamic field simulation results, hydrological and meteorological forecast data covering the prediction period, and three-dimensional hydrodynamic model forecast results, as well as the physical prior bias vectors corresponding to each future prediction time, into the trained time series prediction network, and output the predicted dissolved oxygen deficit value for the future time.

[0007] Further, in step S1, the hydro-meteorological observation sequence is collected by automatic weather stations and hydrological stations, including rainfall, air temperature, wind speed, wind direction, river flow, and atmospheric pressure; the water quality monitoring sequence is collected by water quality monitoring buoy equipment, including dissolved oxygen concentration, water temperature, salinity, chemical oxygen demand, ammonia nitrogen concentration, and chlorophyll concentration; the three-dimensional hydrodynamic model is a three-dimensional baroclinic dynamic model, and the dynamic field simulation results include the flow velocity vector, water level, temperature, salinity, and bottom shear stress of each grid point over time; The tidal boundary conditions, runoff boundary conditions, and tidal level change sequence of the target water area of ​​the three-dimensional hydrodynamic model are obtained to determine the release time of neutral tracer particles, the transport period of neutral tracer particles, and the dynamic state corresponding to the future prediction time. Using the timestamps of the water quality monitoring sequence as a reference, the hydrological and meteorological observation sequence and the dynamic field simulation results are time-aligned, and the dynamic field simulation results are spatially interpolated to extract the dynamic characteristics of the target station, forming training input data with a unified time interval.

[0008] Furthermore, step S2 specifically includes: Water temperature, salinity, and actual dissolved oxygen concentration at each time point are extracted from the water quality monitoring sequence, and the corresponding atmospheric pressure is extracted from the hydro-meteorological observation sequence. The theoretical saturated dissolved oxygen concentration is calculated based on the water temperature and salinity, and the theoretical saturated dissolved oxygen concentration is corrected for atmospheric pressure based on the atmospheric pressure. The difference between the corrected theoretical saturated dissolved oxygen concentration and the actual dissolved oxygen concentration is calculated to obtain the dissolved oxygen deficit value at each time point, and a dissolved oxygen deficit sequence is formed. The theoretical saturated dissolved oxygen concentration was calculated using the APHA formula under low salinity conditions and the Weiss formula under seawater salinity conditions.

[0009] Furthermore, step S3 specifically includes: The three-dimensional hydrodynamic model is driven by atmospheric forcing, tidal and runoff conditions; Neutral tracer particles are released at the upstream boundary section of the three-dimensional hydrodynamic model, and the particle positions are updated step by step based on the three-dimensional velocity field in the dynamic field simulation results. Record the timestamp of the particle's first entry into the target section of the target site and the repeated arrival timestamp of the particle's re-entry into the target section after leaving. Obtain particle transport time samples based on the difference between each arrival timestamp and the corresponding release timestamp. The tidal phase is determined based on the tidal level change at the moment of particle release, the tidal pattern is determined based on the tidal range, and the runoff conditions are determined based on the river flow during the moment of particle release and the corresponding transport period. Kernel density estimation is performed on particle transport time samples under the same dynamic state to obtain the transport time condition arrival probability density function, which is then associated with and stored with the corresponding tidal phase, tidal pattern, and runoff conditions to form the dynamic transport time distribution feature library.

[0010] Furthermore, kernel density estimation is performed on particle transport time samples under the same dynamic state to obtain the transport time conditional arrival probability density function, which is then associated and stored with the corresponding tidal phase, tidal type, and runoff conditions to form the dynamic transport time distribution feature library, specifically including: Kernel density estimation is performed on particle transport time samples under the same dynamic state to obtain the transport time conditional arrival probability density function; Based on the arrival probability density function of the transport time condition, the expected transport time, transport dispersion and reciprocating flow return rate are calculated as auxiliary transport dynamics characteristics. These are then associated with the corresponding tidal phase, tidal pattern and runoff conditions and stored in the dynamic transport time distribution feature library for feature library quality assessment and index optimization. Wherein, the expected transport time is the first moment of the transport time conditional arrival probability density function, which is used to characterize the statistical average time for the water mass to reach the target section from the release location; The degree of transport dispersion is the standard deviation of the arrival probability density function under the transport time condition, which is used to characterize the variability of transport time; The reciprocating flow return rate is the proportion of the number of particles that arrive at least once via a reciprocating flow to the total number of particles that arrive for the first time, and is used to characterize the contribution of tidal reciprocating flow to material retention. The expected transport time, transport dispersion, and reciprocating flow return rate are associated with the corresponding tidal phase, tidal pattern, and runoff conditions and stored in the dynamic transport time distribution feature library.

[0011] Furthermore, step S4 specifically includes: For each future prediction time, iterate through each historical time in the historical time window, take each historical time as the hypothetical particle release time, and determine the dynamic state corresponding to each historical time and the transport period up to the future prediction time. The dynamic state includes tidal phase, tidal pattern and runoff conditions. Based on the dynamic state, the corresponding transport time condition arrival probability density function is obtained by matching or interpolating from the dynamic transport time distribution feature library; Calculate the transport time interval between each historical moment and the predicted future moment, and obtain the conditional arrival probability corresponding to each transport time interval; The probability of each condition being reached is added to a preset positive number, and the natural logarithm is taken. The logarithm is then arranged in the order of historical time moments to generate a physical prior bias vector of the same length as the historical time window. A corresponding physical prior bias vector is generated for each future prediction time moment.

[0012] Furthermore, step S5 specifically includes: The time series prediction network includes an encoder and a decoder based on an attention mechanism. The encoder is used to extract time series features in historical time windows, and the decoder is used to combine time attention weights and input features of future time periods to generate prediction results. The auxiliary water quality parameters in the hydrological and meteorological observation sequence and water quality monitoring sequence, as well as the target station dynamic features extracted from the dynamic field simulation results, are organized into an input feature matrix according to historical time windows. The rows of the input feature matrix correspond to each moment in the historical time window, and the columns correspond to each input feature dimension. The hydrological and meteorological features and target station dynamic features corresponding to the prediction period of the training samples are used as input features for the future period. During the model training phase, the hydrological and meteorological features and target station dynamic features observed synchronously in history are used as input features for the future period, i.e., assuming perfect forecast conditions. The physical prior bias vector is used as a logarithmic prior bias, and its influence intensity is adjusted by a non-negative scaling factor obtained by non-negative mapping of learnable parameters, and then fused with the temporal attention score. Using the dissolved oxygen deficit sequence as the supervision target, the network parameters and the learnable parameters are jointly updated through backpropagation to obtain the attention posterior probability distribution and the trained dissolved oxygen prediction model. Specifically, a loss function is constructed based on the prediction error between the predicted dissolved oxygen deficit at each training prediction time and the actual value at the corresponding time in the dissolved oxygen deficit sequence described in step S2. The gradient of the loss function relative to the temporal prediction network parameters and the learnable parameters used to generate the non-negative scaling coefficients is calculated through backpropagation, and the network parameters and learnable parameters are jointly updated based on the gradient. Based on the updated network parameters and learnable parameters, the corrected attention score at each training prediction time is recalculated according to the above method and normalized by the Softmax function to obtain the temporal attention weight distribution, which is used as the attention posterior probability distribution. After training convergence, the network parameters and learnable parameters are fixed, and forward propagation is performed on the validation set, test set, or sample to be analyzed to obtain the attention posterior probability distribution at each prediction time.

[0013] Furthermore, the physical prior bias embedding temporal attention scoring function specifically includes: The hidden state of the decoder at the current future prediction time is obtained as the query feature, the hidden state of each historical time step output by the encoder is obtained as the key feature, and the initial attention score of each historical time step is calculated. Obtain the physical prior bias vector corresponding to the current future prediction time, adjust the influence intensity of the physical prior bias vector through the non-negative scaling factor, and fuse the adjusted physical prior bias with the initial attention score of the corresponding historical time step to obtain the corrected attention score; The temporal attention weights are calculated based on the corrected attention scores, and the historical hidden state sequences output by the encoder are weighted and aggregated to generate temporal context features. The time context features and the future time period input features corresponding to the current future prediction time are input into the decoder to obtain the dissolved oxygen deficit prediction value for that future prediction time.

[0014] Furthermore, step S6 specifically includes: Obtain the hydrological and meteorological observation sequence and water quality monitoring sequence up to the current moment, as well as the hydrological and meteorological forecast data and three-dimensional hydrodynamic model forecast results covering the forecast period; form a historical window input feature matrix based on the hydrological and meteorological observation sequence, water quality monitoring sequence and target station dynamic characteristics extracted from the dynamic field simulation results; and form a future forecast feature sequence based on the hydrological and meteorological forecast data and three-dimensional hydrodynamic model forecast results. For each future prediction time, a physical prior bias vector is generated according to step S4, and arranged according to the future prediction time to form a physical prior bias matrix; the rows of the physical prior bias matrix correspond to the future prediction time, the columns correspond to each historical time in the historical time window, and the matrix elements are the physical prior bias values ​​of the corresponding historical time. The historical window input feature matrix, future forecast feature sequence, and physical prior bias matrix are input into the trained time series prediction network. The physical prior bias vector corresponding to each future prediction time is called respectively, and the dissolved oxygen deficit prediction value for each future prediction time is output.

[0015] Furthermore, following step S6, the process also includes the identification of the dominant process and the time window for hypoxia risk: For each prediction time, the arrival probabilities of conditions at each historical time used when generating the corresponding physical prior bias vector are normalized to obtain the physical prior probability distribution, and the attention posterior probability distribution at the same prediction time is obtained. Calculate the Pearson correlation coefficient and peak position deviation between the physical prior probability distribution and the attention posterior probability distribution; When the correlation coefficient is greater than the first threshold and the peak position deviation is less than the second threshold, it is determined to be the physical transport-dominated period; when the correlation coefficient is less than the third threshold or the peak position deviation is greater than the fourth threshold, it is determined to be the local biochemical process-dominated period; all other cases are determined to be the transition period. An interpretable attribution method was used to analyze the contribution of each input feature to the dissolved oxygen deficit prediction value, identify the stage contributions of exogenous oxygen-consuming substance transport processes and local biochemical oxygen-consuming processes, and combine the dominant processes in continuous future prediction periods with the dissolved oxygen deficit prediction results to determine the time window of hypoxia risk under high runoff conditions.

[0016] The present invention has the following beneficial effects: First, this invention significantly improves the prediction accuracy and robustness of dissolved oxygen deficit in tidal river networks. By employing a three-dimensional baroclinic dynamic model and Lagrange particle tracking technology, this invention constructs a dynamic transport time distribution feature library characterizing the probability distribution of upstream water mass transport time. This feature library is then transformed into an attention scoring function embedded in a time-series prediction network using physical prior biases. This allows the network to focus on key historical moments when upstream pollutants arrive at target sites after being transported through complex hydrodynamic processes, guided by physical laws. The influence strength of the physical prior bias is adjusted by non-negative scaling coefficients obtained from learnable parameters through non-negative mapping, and this is fused with the initial attention score to form a collaborative attention allocation mechanism combining physical priors and data-driven approaches. This leverages the advantages of physical models in characterizing long-range transport time lags while retaining the data-driven model's ability to fit local biochemical processes, thereby effectively improving the prediction accuracy and extrapolation robustness under extreme hydrological conditions in tidal river networks.

[0017] Second, adaptive adjustment of the physical prior guidance intensity is achieved. This invention sets the scaling factor used to adjust the influence intensity of the physical prior bias as a learnable parameter jointly optimized by backpropagation. This scaling factor, through non-negative mapping, ensures the correctness of the physical guidance direction. During training, the network jointly optimizes the scaling factor based on the prediction error, enabling the model to adaptively determine the optimal contribution level of the physical prior bias according to the training data. This effectively reduces the risk of physical model error propagating to the prediction results and maintains stable prediction performance under different water areas and hydrological conditions.

[0018] Third, this invention enhances the interpretability of the prediction process and its ability to support management decisions. By comparing the Pearson correlation coefficient and peak position deviation between the physical prior probability distribution and the attention posterior probability distribution actually generated after network training convergence, this invention can quantitatively identify the dominant influencing mechanism of dissolved oxygen deficit in the current period, dividing the hypoxia process into a physical transport-dominated period, a local biochemical process-dominated period, and a transition period. Simultaneously, an interpretable attribution method is used to analyze the contribution of each input feature to the dissolved oxygen deficit prediction value, identifying the stage-specific contributions of exogenous oxygen-consuming substance transport processes and local biochemical oxygen-consuming processes. Combined with the dominant processes and dissolved oxygen deficit prediction results for consecutive future prediction periods, the invention determines the hypoxia risk time window under high runoff conditions, providing a clear predictive and diagnostic basis for differentiated water quality management. Attached Figure Description

[0019] Figure 1 A flowchart of the dissolved oxygen prediction method for tidal river networks based on time-series and dynamic processes provided by the present invention; Figure 2 Radar chart comparing the comprehensive performance of the deep fusion scheme of physical transport laws and data-driven learning provided by this invention with ordinary fusion schemes. Detailed Implementation

[0020] To further understand the content of this invention, a detailed description of the invention is provided in conjunction with the accompanying drawings and embodiments. The specific embodiments described herein are for illustrative purposes only and are not intended to limit the invention. It should also be noted that, for ease of description, only the parts relevant to the invention are shown in the accompanying drawings.

[0021] The application scenario of this invention is to provide water quality early warning for tidal river networks such as river estuaries and nearshore sea areas, provide environmental monitoring departments with rolling predictions of future dissolved oxygen depletion, and support operational applications such as pollution source tracing, ecological scheduling, and hypoxia early warning.

[0022] This embodiment provides a method for predicting dissolved oxygen in tidal river networks based on time-series and dynamic processes, such as... Figure 1 As shown, it includes: Step S1: Obtain the hydrological and meteorological observation sequence, water quality monitoring sequence, and multi-parameter dynamic field simulation results output by the three-dimensional hydrodynamic model from the target station.

[0023] In this embodiment, the target water area is a tidal river network, estuary, or nearshore sea area that is significantly affected by both tides and runoff. Multi-source time-series data of the target water area are acquired. This multi-source time-series data includes hydrological and meteorological observation sequences, water quality monitoring sequences, and dynamic field simulation sequences generated by a three-dimensional fluid dynamics model. These three data are arranged in a time series according to their respective sampling frequencies.

[0024] The hydro-meteorological observation sequence is obtained through automatic weather stations and hydrological stations deployed along the target water area. The automatic weather stations integrate temperature sensors, tipping bucket rain gauges, and wind speed and direction sensors, measuring rainfall, air temperature, wind speed, and wind direction on an hourly sampling cycle. The hydrological stations simultaneously record river flow at the boundary of the target water area; the flow data is calculated from the water level-flow curve by the hydrological stations, with a time resolution typically hourly. The observation data from the automatic weather stations and hydrological stations are transmitted in real-time to the data center's storage system via wireless communication modules. The hydro-meteorological observation sequence includes rainfall, air temperature, wind speed, wind direction, river flow, and atmospheric pressure.

[0025] Water quality monitoring data is collected by automatic monitoring buoys deployed in the target water area. These buoys are equipped with multi-parameter water quality analyzers, using thermistor method to measure water temperature, electrochemical sensor to measure dissolved oxygen concentration, electrode method to measure salinity, and ultraviolet absorption spectroscopy or chemical analysis to measure chemical oxygen demand (COD), ammonia nitrogen concentration, and chlorophyll concentration. The sensor sampling cycle is set to once per hour. Data is converted into digital signals by the built-in data processing unit and transmitted to the data center in real time via a wireless communication module. The automatic monitoring buoys operate in continuous mode, automatically collecting water quality data according to a preset cycle. The water quality monitoring sequence includes dissolved oxygen concentration, water temperature, salinity, COD, ammonia nitrogen concentration, and chlorophyll concentration.

[0026] The dynamic field simulation sequence was generated by a three-dimensional hydrodynamic model. This three-dimensional hydrodynamic model employs a three-dimensional baroclinic dynamic mode and is constructed based on numerical frameworks such as environmental hydrodynamic models or finite volume oceanographic models. The landward boundary of the model is set above the tidal boundary, and the seaward boundary extends to the shallow waters of the continental shelf. The model inputs include tidal level data synthesized from multiple tidal harmonic constants provided by the global tidal model, initial field and boundary conditions of three-dimensional temperature and salinity profiles provided by regional ocean reanalysis data, and atmospheric forcing data such as wind stress, heat flux, and net shortwave radiation provided by atmospheric reanalysis data. The model uses sigma coordinates or a hybrid coordinate system in the vertical direction, and outputs time-series data of velocity vectors, water level, temperature, salinity, and bottom shear stress at each grid point in the target water area, with an output time interval of one hour. The dynamic field simulation results include the time-varying velocity vectors, water level, temperature, salinity, and bottom shear stress at each grid point.

[0027] The tidal boundary conditions, runoff boundary conditions, and tidal level change sequence of the target water area of ​​the three-dimensional hydrodynamic model are obtained to determine the dynamic state corresponding to the release time of neutral tracer particles, the transport period of neutral tracer particles, and the future prediction time.

[0028] The three data sequences mentioned above differ in temporal resolution, spatial scale, and data format, necessitating spatiotemporal alignment and preprocessing. Using the timestamps of the water quality monitoring sequence as a benchmark, linear interpolation resampling is performed on the hydrological and meteorological observation sequence and the dynamic field simulation sequence: For daily or hourly records of rainfall and river flow in the hydrological and meteorological observation sequence, linear interpolation is used between adjacent observation times to calculate missing data values; for data with a 30-minute output interval in the dynamic field simulation sequence, the output values ​​at the hourly interval are extracted or the arithmetic mean of adjacent time steps is calculated as the data values ​​at the hourly interval, ensuring that each data source corresponds to a unique data value at the same timestamp.

[0029] Spatial interpolation is performed on the dynamic field simulation sequence. The geographical coordinates of the target monitoring station typically do not coincide with the grid center point of the 3D fluid dynamics model. Spatial interpolation is needed to calculate the variable values ​​at the target station location from the simulated values ​​of surrounding grid points. A bilinear interpolation method is used: the model grid cell containing the target station is determined, defined by four grid center points. The normalized coordinates of the target station relative to the grid cell are calculated. Linear interpolation is performed first in the horizontal direction, and then in the vertical direction to obtain the temporal characteristics of flow velocity, water level, and bottom shear stress at the target station. For vertically distributed variables, horizontal interpolation is performed first to obtain the variable values ​​for each vertical layer, and then vertical interpolation or vertical averaging is performed based on the water depth at the target station. The target station is a discretely distributed water quality monitoring point within the river network, serving as the final output point for dissolved oxygen prediction results. Multi-source time-series monitoring data can be directly collected.

[0030] Outliers in multi-source time-series data are identified and removed using pre-defined statistical rules. A sliding window-based statistical test method is employed: a sliding window is set on the time series, with a window length ranging from 24 to 72 hours. The mean and standard deviation of each data variable within the window are calculated. If a data point deviates from the mean by more than three times the standard deviation, it is identified as an outlier and removed. For isolated single-point outliers, a reasonable value is estimated from adjacent normal values ​​using linear interpolation and then replaced. For consecutive multi-point outliers, the data for the entire outlier period is deleted, and this period is not included in subsequent model training and prediction. When the dissolved oxygen concentration exceeds the conventional physical threshold based on historical data (below 0 mg / L or above 20 mg / L, where 20 mg / L is the 99.5th percentile value of the training set), the data for the corresponding time step is removed.

[0031] Through the above processing, the hydrological and meteorological observation sequence, water quality monitoring sequence, and dynamic field simulation sequence are aligned to a unified timestamp in the time dimension, and each variable corresponds to the geographical location of the target monitoring station in the spatial dimension. After the outliers in each sequence are identified and processed, a complete and continuous multi-source time series data matrix is ​​formed, which forms training input data with a unified time interval.

[0032] Step S2: Calculate the theoretical saturated dissolved oxygen concentration and dissolved oxygen deficit sequence based on the water quality monitoring sequence and the hydrological and meteorological observation sequence.

[0033] In this embodiment, water temperature, salinity, and actual dissolved oxygen concentration at each time point are extracted sequentially from the water quality monitoring sequence, and the corresponding atmospheric pressure is extracted from the hydrological and meteorological observation sequence. This water quality monitoring sequence is obtained by automatic monitoring buoys deployed in the target water area, collecting data hourly and then performing spatiotemporal alignment. Each time point record includes multiple parameter fields such as water temperature, salinity, and actual dissolved oxygen concentration. Water temperature is measured using the thermistor method, salinity using the electrode method, and actual dissolved oxygen concentration using the electrochemical sensor method. Atmospheric pressure data is obtained by a barometric pressure sensor mounted on the automatic monitoring buoy or by synchronous observation from coastal meteorological stations, and then correlated with the water quality monitoring sequence after spatiotemporal alignment. These four parameters, at the same timestamp, together constitute the input data set for calculating dissolved oxygen deficit.

[0034] The theoretical saturated dissolved oxygen concentration (DSOC) is calculated based on water temperature and salinity. The theoretical saturated DOC is calculated using either the well-known APHA formula (applicable to salinities below 5‰) or the Weiss formula (applicable to salinities not less than 5‰). The choice of either formula depends solely on the salinity value, which is provided by synchronous salinity data extracted from a water quality monitoring sequence. The theoretical saturated DOC calculated using the above formula corresponds to a value under current standard atmospheric pressure conditions.

[0035] The theoretical saturated dissolved oxygen concentration is corrected for atmospheric pressure using synchronous atmospheric pressure. The synchronous atmospheric pressure data in the water quality monitoring sequence reflects the actual atmospheric pressure at the current monitoring location, while the theoretical saturated dissolved oxygen concentration calculated using the APHA and Weiss formulas corresponds to standard atmospheric pressure conditions. When the actual atmospheric pressure deviates from standard atmospheric pressure, pressure correction is required. The correction method involves multiplying the theoretical saturated dissolved oxygen concentration by the ratio of the actual atmospheric pressure to the standard atmospheric pressure. After pressure correction, the theoretical saturated dissolved oxygen concentration at the current time point, corrected for temperature, salinity, and pressure, is obtained.

[0036] The difference between the corrected theoretical saturated dissolved oxygen concentration and the actual dissolved oxygen concentration is calculated to obtain the dissolved oxygen deficit value at each time point. The dissolved oxygen deficit is defined as the difference between the theoretical saturated dissolved oxygen concentration and the actual dissolved oxygen concentration, reflecting the degree of oxygen deficiency in the water body relative to saturation. A positive difference indicates an oxygen deficit, meaning the water body is unsaturated; a zero difference indicates saturation; and a negative difference indicates supersaturation. The above steps of extracting water temperature and salinity, calculating the theoretical saturated dissolved oxygen concentration, pressure correction, and calculating the difference are repeated for each time point to generate a dissolved oxygen deficit sequence of the same length as the water quality monitoring sequence. Each value in this dissolved oxygen deficit sequence corresponds to the dissolved oxygen deficit state at the same time point in the water quality monitoring sequence, forming continuous time-series data.

[0037] Using dissolved oxygen deficit as a predictive target variable offers several advantages. First, the numerical range of dissolved oxygen deficit is relatively stable, avoiding baseline drift caused by seasonal temperature variations. Second, the relationship between dissolved oxygen deficit and pollution loads and ecological processes is more direct, facilitating the model's capture of the impacts of biochemical reactions and physical transport on water quality. Third, the positive and negative values ​​of dissolved oxygen deficit have clear ecological significance; positive values ​​indicate hypoxia risk, while negative values ​​indicate supersaturation, aiding in subsequent water quality health assessments and early warning decisions.

[0038] Step S3: Based on the three-dimensional velocity field tracking neutral tracer particles in the dynamic field simulation results, construct a dynamic transport time distribution feature library corresponding to tidal phase, tidal pattern and runoff conditions.

[0039] In this embodiment, the three-dimensional hydrodynamic model is driven by atmospheric forcing, tidal, and runoff conditions. The landward boundary of the model is set above the tidal boundary, and the seaward boundary extends to the shallow shelf region. Setting the landward boundary above the tidal boundary ensures that the water flow at this boundary is controlled by runoff, and the influence of tidal backpropagation is attenuated to a negligible level, thus allowing flow rate to be used as the boundary condition. The seaward boundary extending to the shallow shelf region is far from the direct influence of estuarine plumes and density fronts, and the tidal amplitude and phase at the boundary can be accurately provided by a global tidal model. The model uses sigma coordinates or a hybrid coordinate system in the vertical direction to fit underwater irregular topography. This coordinate system divides the vertical direction into equal layers according to water depth, maintaining the same number of vertical layers in both shallow and deep water areas, which can adapt to the characteristics of tidal river networks and the dramatic water depth variations in estuary areas. The three-dimensional hydrodynamic model operates based on the three-dimensional velocity field data in the dynamic field simulation sequence obtained in step S1.

[0040] Neutral tracer particles are released in batches at the upstream boundary section of a three-dimensional hydrodynamic model. These neutral tracer particles have the same density as the surrounding water, are not subject to gravitational settling or buoyancy lifting, and move entirely with the water mass. The release locations cover the entire water depth and lateral range of the upstream boundary section. Vertically, particles are distributed at equal intervals according to water depth, ensuring that the transport characteristics of water masses at different depths are sampled. Laterally, particles are evenly distributed in the main channel and sidebank areas to avoid missing transport paths in low-velocity areas along the shore. The number of particles released per batch is set based on the upstream boundary section area, typically between five thousand and twenty thousand, to ensure the stability of the statistical results.

[0041] After particle release, the 3D spatial coordinates of each particle are updated time-step by time based on the 3D velocity field in the dynamic field simulation sequence. The 3D velocity field includes east-west, north-south, and vertical components. The particle position update uses a numerical integration method, progressing along the time direction in a second-order Runge-Kutta scheme, with the particle tracking time step consistent with the output time step of the 3D hydrodynamic model. During the particle's motion, tracking stops when it reaches the seaward or lateral open boundary of the model and its velocity direction points outward from the boundary; when it touches the shoreline or seabed, the particle's normal velocity component reverses, while the tangential component remains unchanged, simulating the particle's sliding behavior along the boundary. The 3D velocity components required for particle coordinate update are obtained from the velocity values ​​of surrounding mesh nodes using trilinear interpolation.

[0042] The system records the timestamp of a particle's first entry into the target section corresponding to the target site, as well as the timestamp of its re-entry into the target section after leaving. When a particle's horizontal coordinate enters the target section and its direction of movement points towards the target section, the system records the timestamp of the particle's first arrival at the target section. This timestamp is based on the particle's release time, recording the time interval from release to arrival at the target section. In tidal river networks and estuary areas, the effects of tidal currents are significant. Some particles, after initially passing the target section, will move downstream with the ebb current and then return to the vicinity of the target section during the next high tide. Therefore, after the particle's initial arrival, the system continues to track the particle and records the timestamp of its re-entry into the target section. The condition for a re-arrival is that the particle left the grid containing the target section after its initial arrival and subsequently re-entered the grid. The particle transport time sample is obtained based on the difference between each arrival timestamp and its corresponding release timestamp.

[0043] The tidal phase is determined based on the tidal level changes at the moment of particle release, the tidal pattern is determined based on the tidal range, and the runoff conditions are determined based on the river flow during the particle release time and the corresponding transport period. The tidal phase is divided according to the tidal level changes in the target water area, including six phases: early tide, mid-tidal period, high tide, early ebb tide, mid-ebb tide, and low tide. It also distinguishes between spring tides and neap tides. The determination of the period from early tide to mid-ebb tide is based on the sign of the tidal level change rate and the normalized position of the tidal level relative to the high and low tide levels of the day: a positive tidal level change rate with the normalized tidal level between 0 and 1 / 3 indicates early tide, between 1 / 3 and 2 / 3 indicates mid-tidal period, and between 2 / 3 and 1 indicates high tide; a negative tidal level change rate with the normalized tidal level between 2 / 3 and 1 indicates early ebb tide, between 1 / 3 and 2 / 3 indicates mid-ebb tide, and between 0 and 1 / 3 indicates low tide. The distinction between spring tides and neap tides is based on the ratio of the daily tidal range to the monthly average tidal range: a ratio greater than 1.2 is considered a spring tide, less than 0.8 is considered a neap tide, and a ratio between 0.8 and 1.2 is considered a transitional tide.

[0044] Kernel density estimation is performed on particle transport time samples under the same dynamic state to obtain the transport time conditional arrival probability density function. Kernel density estimation uses each arrival timestamp as a sample point, sets a Gaussian kernel function, and determines the bandwidth parameter of the kernel function through cross-validation: the samples are randomly divided into training and validation subsets in a 4:1 ratio. The log-likelihood function value of the validation subset is used as the evaluation index. The search is performed within a preset bandwidth search range [0.5h, 12h] with a step size of 0.5h, selecting the bandwidth value that maximizes the log-likelihood function value of the validation subset to obtain a probability density function that balances fitting accuracy and smoothness. After kernel density estimation, a continuous curve is obtained with transport time as the horizontal axis and probability density value as the vertical axis. This curve represents the transport time distribution of water masses released from the upstream boundary section reaching the target section via different paths.

[0045] Based on the transport time conditional arrival probability density function, the expected transport time, transport dispersion, and reciprocating flow return rate are calculated and stored as auxiliary transport dynamics features in the dynamic transport time distribution feature library for feature library quality assessment and index optimization. The expected transport time is the first moment of the transport time conditional arrival probability density function, used to characterize the statistically average time for a water mass to reach the target section from its release location; the transport dispersion is the standard deviation of the transport time conditional arrival probability density function, used to characterize the variability of transport time; and the reciprocating flow return rate is the proportion of particles that have undergone at least one reciprocating arrival to the total number of particles arriving for the first time, used to characterize the contribution of tidal reciprocating flow to material retention.

[0046] Expected transport time is the mathematical expectation of the transport time probability density function, calculated as the integral of the product of transport time and the corresponding probability density value over the entire time domain. Expected transport time represents the statistically average time for a water mass to reach the target station from upstream and is a key indicator describing transport efficiency. Under different hydrological and meteorological conditions, expected transport time exhibits regular variations. During the wet season, large runoff and stronger ebb tide velocity shorten the expected transport time. During the dry season, small runoff and relatively stronger tidal effects lead to significant reciprocating flow, extending the expected transport time. During spring tides, large tidal ranges and fast tidal currents improve transport efficiency and decrease the expected transport time. During neap tides, small tidal ranges and weak tidal currents increase the water mass's residence time, increasing the expected transport time. By analyzing the time series of expected transport time, the main controlling factors of transport forces and their seasonal variation patterns can be identified.

[0047] Transport dispersion is quantified by calculating the standard deviation of the transport time probability density function. The standard deviation is the square root of the weighted average of the squares of the differences between each transport time and the expected transport time, with the weights being the corresponding probability densities. A larger standard deviation indicates greater variability in transport time, and more significant differences in transport time caused by different paths or dynamic processes. In estuaries with complex topography and multiple dynamic processes, transport dispersion is typically large, with the standard deviation reaching more than 50% of the expected transport time. High dispersion implies strong uncertainty in pollutant arrival times; some water masses may arrive rapidly, causing short-term water quality deterioration, while others may remain for extended periods, continuously releasing pollution loads. Changes in transport dispersion are closely related to tidal type. During spring tides, the flow field structure is simple, particle transport paths converge, and dispersion is low. During neap tides, flow velocity is weak, and particles are significantly affected by secondary circulation and local topography, resulting in dispersed paths and increased dispersion.

[0048] The reciprocating flow return rate is defined as the proportion of particles that arrive at least once via a round trip to the total number of particles arriving on their first visit. This indicator directly reflects the contribution of tidal reciprocating flows to material retention. A high reciprocating flow return rate indicates that the water mass repeatedly lingers near the target area, and the actual residence time of pollutants far exceeds the expected transport time, which is conducive to local biochemical reactions but also increases the risk of pollutant accumulation. In the Pearl River Estuary case study, the reciprocating flow return rate can reach over 0.6 during neap tides, indicating that more than half of the water mass undergoes at least one round trip. However, during spring tides with large runoff conditions, the reciprocating flow return rate drops to around 0.3, showing a clear unidirectional downward discharge characteristic of the water mass. The seasonal and tidal variation patterns of the reciprocating flow return rate provide important physical constraints for understanding estuarine material retention and water quality evolution.

[0049] The arrival probability density function of the transport time conditions is associated and stored with the corresponding tidal phase, tidal pattern, and runoff conditions to form the dynamic transport time distribution feature library. The dynamic transport time distribution feature library is a three-dimensional data structure. The first dimension corresponds to the release time (tidal phase and runoff conditions), the second dimension corresponds to the transport time, and the third dimension is the probability density value. Each row in the library represents the transport time probability density function corresponding to a certain release time, and each column represents the probability density change of a specific transport time at different release times. This library comprehensively reflects the modulation effect of multiple factors such as tidal cycle, runoff variation, and wind field fluctuations on the transport process, and serves as the foundational data for the subsequent construction of physical prior biases.

[0050] In the Pearl River Estuary case study, particle tracking analysis revealed that during the dry season (neap tide), the expected transport time for particles released from the Xijiang Modaomen channel to reach the Lingdingyang monitoring station was approximately 36 hours, with a standard deviation of about 12 hours and a return flow rate of about 0.4, indicating that 40% of the particles would return after their initial passage. However, during the wet season (spring tide), with increased runoff and tidal range, the expected transport time shortened to 24 hours, the standard deviation increased to 18 hours, and the return flow rate rose to 0.6, indicating faster but more variable water mass transport and more pronounced back-and-forth retention. These temporal evolution patterns of transport dynamics provide important physical prior knowledge for subsequent water quality prediction models.

[0051] Step S4: For each future predicted time, determine the dynamic state based on hydrological and meteorological forecast data and the forecast results of the three-dimensional hydrodynamic model, and combine the transport time interval from each historical time to the future predicted time to obtain the conditional arrival probability from the dynamic transport time distribution feature library, and generate the corresponding physical prior bias vector.

[0052] In this embodiment, the rows of the dynamic transport time distribution feature library correspond to different release times (tidal phase and runoff conditions), the columns correspond to different transport times, and the library elements are probability density values. This library represents the probability distribution of water masses released upstream under specific tidal phase and runoff conditions reaching the target site via different transport times. When generating physical prior biases, it is necessary to map the physical transport patterns in this library to the historical time windows of the time series prediction network.

[0053] For each future prediction time, the historical time in the historical time window is traversed, and each historical time is taken as the hypothetical particle release time. The dynamic state corresponding to each historical time and the transport period up to the future prediction time is determined. The dynamic state includes tidal phase, tidal pattern and runoff conditions.

[0054] The tidal phase and runoff conditions corresponding to the target time to be predicted are determined. The target time is the future time at which the time-series prediction network needs to output the dissolved oxygen deficit prediction value. The tidal phase at this time is determined based on the tidal level forecast data of the target water area or the water level sequence output by the three-dimensional hydrodynamic model. The determination of the tidal phase is based on the rising or falling tide stage and the magnitude of the tidal range at the target time, and is divided into six phases: early rising tide, mid-rising tide, high tide, early falling tide, mid-falling tide, and low tide. At the same time, the conditions of spring tide and neap tide are distinguished. The criteria for determining the ebb and flow of tide are the sign of the tidal level change rate and the normalized position of the tide relative to the high and low tide levels of the day: a positive tidal level change rate with the normalized tide level between 0 and 1 / 3 indicates the ebb and flow of tide, between 1 / 3 and 2 / 3 indicates the middle of the ebb and flow of tide, and between 2 / 3 and 1 indicates high tide; a negative tidal level change rate with the normalized tide level between 2 / 3 and 1 indicates the ebb and flow of tide, between 1 / 3 and 2 / 3 indicates the middle of the ebb and flow of tide, and between 0 and 1 / 3 indicates low tide. The distinction between spring tides and neap tides is based on the ratio of the daily tidal range to the monthly average tidal range: a ratio greater than 1.2 indicates a spring tide, less than 0.8 indicates a neap tide, and between 0.8 and 1.2 indicates a transitional tide. The runoff condition sequence is determined.

[0055] The corresponding transport time condition arrival probability density function is obtained by matching or interpolating from the dynamic transport time distribution feature library based on the described dynamic state. The dynamic transport time distribution feature library has already generated corresponding probability density functions for release times at different tidal phases during its construction; therefore, the corresponding row in the library can be located based on the tidal phase at the target time. If a row in the library perfectly matches the tidal phase at the target time, the probability density function corresponding to that row is directly extracted; if no exact match exists, the two rows with the closest tidal range and tidal phase angle are selected for weighted interpolation to obtain the probability density function corresponding to the target tidal phase. The interpolation weights are determined based on the tidal phase difference; the smaller the difference, the greater the weight.

[0056] The transport time interval between each historical moment and the predicted future moment is calculated, and the conditional arrival probability corresponding to each transport time interval is obtained. The preset historical time window length for the time series prediction network is obtained. The historical time window length is set to 168 hours, covering the main transport time scales including multiple complete tidal days and typical runoff impulse response periods. Its value is determined based on the hydrodynamic transport time scale of the target water area, ensuring coverage of the main transport time distribution range of upstream water masses arriving at the target station in the tidal river network. The historical time window length determines the range of historical data that the network can backtrack when processing the current prediction task. Using the target moment as a reference, the time difference between each historical moment in the historical time window and the target moment is used as the transport time interval, thus making each historical moment in the historical time window correspond to a unique transport time interval. This process maps the historical time window from absolute time coordinates to relative time coordinates with the target moment as the endpoint. This relative coordinate has the same time dimension as the transport time dimension of the dynamic transport time distribution feature library.

[0057] For each historical moment within the historical time window, the time difference between that historical moment and the target moment is used as the transport time interval. The corresponding probability density value is then retrieved from the extracted probability density function. If the transport time interval falls within the domain of the probability density function, i.e., within the transport time coverage of the dynamic transport time distribution feature library, the corresponding probability density value is obtained. If the transport time interval exceeds the transport time coverage, it indicates that the upstream water mass cannot reach the target station within that time interval at that physical transport time scale, and the corresponding probability density value is set to zero. The upper limit of the transport time coverage is determined by the longest duration of particle tracking, and the lower limit is determined by the temporal resolution of the tracking.

[0058] The arrival probabilities of each condition are added to a preset positive number, and the natural logarithm is taken. These logarithms are then arranged in chronological order according to historical timeframes to generate a physical prior bias vector of the same length as the historical time window. A corresponding physical prior bias vector is generated for each future prediction timeframe. The dimension of the physical prior bias vector is equal to the sequence length of the historical time window. The values ​​of its elements are distributed between negative and zero (the result of taking the logarithm after adding the preset positive number). Elements with larger values ​​(i.e., smaller absolute values) correspond to historical times with higher physical transport probabilities, while elements with smaller values ​​(i.e., larger absolute values) correspond to historical times with lower physical transport probabilities. The purpose of adding the preset positive number is to avoid negative infinity when taking the logarithm of zero probability; the preset positive number is usually a very small positive number, such as 1e-8. This vector serves as the physical prior bias input to guide the network's attention mechanism during subsequent training and prediction processes of the temporal prediction network.

[0059] For example, taking a tidal estuary as an example, suppose the target time for a prediction is 10:00 AM on a certain day, at which time the tidal phase is mid-ebb tide. Based on the predicted runoff conditions, the historical time window length is 168 hours. The particle transport time distribution in the dynamic transport time distribution feature library during mid-ebb tide exhibits a bimodal structure, with the main peak at 24 hours and the secondary peak at 48 hours. After mapping this distribution inversely to the historical time window, the physical prior bias vector shows significant peaks at positions 24 hours and 48 hours before the target time.

[0060] Step S5: Embed the physical prior bias obtained by logarithmic transformation of the conditional arrival probability into the attention scoring function of the time series prediction network. Take the hydro-meteorological observation sequence, water quality monitoring sequence, dynamic field simulation results and hydro-meteorological forecast characteristics and dynamic forecast characteristics corresponding to the prediction period as input, and take the dissolved oxygen deficit sequence as the supervision target. Train the network parameters and physical prior strength parameters through backpropagation.

[0061] In this embodiment, the temporal prediction network adopts an encoder-decoder architecture, with both the encoder and decoder using a long short-term memory network as the basic recurrent unit. The encoder is used to extract temporal features from historical time windows, and the decoder is used to combine temporal attention weights and input features from future time periods to generate prediction results.

[0062] The encoder consists of two stacked long short-term memory (LSM) networks, each with a 128-dimensional hidden state. The historical time window is set to 168 hours, covering major transport timescales including multiple complete tidal days and typical runoff impulse response periods. At each time step, after receiving a multi-dimensional input feature vector, the encoder updates the cell and hidden states through LSM network units, outputting a hidden state sequence of the same length as the historical time window. The encoder includes a two-stage spatial attention weight optimization module embedded between the two LSM networks: the first stage performs a dot product operation between the initial hidden state of the first LSM network at the current time step and the input water quality spatial feature matrix, then normalizes using the Softmax function to obtain preliminary spatial weights; the second stage inputs the spatial context vector obtained by weighting and aggregating the preliminary spatial weights into the second LSM network to obtain the updated hidden states, performs another dot product operation with the water quality spatial feature matrix, and outputs the corrected spatial weights, which are used for the final spatial weighting of the input features.

[0063] The decoder consists of two stacked long short-term memory (LSM) networks, with a hidden state dimension of 128 and an output dimension of 1 (corresponding to the dissolved oxygen deficit prediction). The decoder includes a two-stage temporal attention weight optimization module embedded between the two LSM networks: the first stage matches the hidden state of the decoder's first-layer LSM network at the current prediction time step with the entire sequence of historical hidden states output by the encoder, and obtains preliminary temporal weights through Softmax normalization; the second stage concatenates the temporal context features obtained by weighting and aggregating the preliminary temporal weights with the decoder's current hidden state and inputs it into the second-layer LSM network to obtain the updated hidden state, which is then matched again with the entire sequence of historical hidden states output by the encoder to output the corrected temporal weights.

[0064] The auxiliary water quality parameters from the hydrological and meteorological observation sequences and water quality monitoring sequences, as well as the target station dynamic features extracted from the dynamic field simulation results, are organized into an input feature matrix according to historical time windows. The rows of the input feature matrix correspond to each moment in the historical time window, and the columns correspond to each input feature dimension, for a total of 15 dimensions. During the model training phase, the hydrological and meteorological features and target station dynamic features from historical synchronous observations are used as input features for future periods, i.e., assuming perfect forecast conditions. During the operational forecast phase, the hydrological and meteorological forecast features for future periods output by numerical weather prediction models or statistical forecast models, as well as the target station dynamic forecast features from the three-dimensional hydrodynamic model forecast results, are used as input features for future periods.

[0065] The physical prior bias vector generated in step S4 is used as the logarithmic prior bias. Its influence intensity is adjusted by a non-negative scaling factor obtained from the learnable parameters via Softplus non-negative mapping, and then fused with the temporal attention score. The mathematical form of the Softplus non-negative mapping is s = ln(1 + e^γ), where γ is the learnable parameter, s is the output non-negative scaling factor, and the initial value is set to 1.0. The embedding operation of the physical prior bias is implemented in the temporal attention scoring function. Traditional scoring functions only include the dot product of the query vector and the key vector; this dot product value reflects the similarity between the two vectors, with higher similarity resulting in a higher score. This invention, based on the dot product scoring, superimposes the corresponding physical prior bias value at each moment in the historical sequence. The physical bias value comes from a dynamic transport time distribution feature library. Specifically, the transport time probability density function corresponding to the current prediction moment is mapped backward to the historical time axis. Each historical moment obtains a probability density value, which, after adding a preset positive number and taking the natural logarithm, becomes the physical bias at that moment.

[0066] The superposition operation is implemented through addition. The corrected score equals the original dot product score plus the physical bias value multiplied by a learnable scaling factor. The scaling factor adjusts the contribution of the physical bias to the score, initially set to one, and automatically learned and optimized by the network during training. When the scaling factor is close to zero, the influence of the physical bias is suppressed, and the network degenerates into a purely data-driven mode. When the scaling factor is large, the physical bias significantly alters the score distribution, forcing the network to focus on key moments indicated by the physical model. The adaptive adjustment of the scaling factor allows the network to dynamically balance the weights of the data and the physical model based on their consistency. The corrected score sequence is normalized using the Softmax function to transform it into the final attention weight distribution. Normalization ensures that the sum of the weights at all historical moments is one, satisfying the probability distribution constraints.

[0067] Using dissolved oxygen deficit sequences as the supervised learning target, network parameters and the learnable parameters are jointly updated via backpropagation to obtain the attention posterior probability distribution and a trained dissolved oxygen prediction model. Training employs the backpropagation algorithm, updating network parameters by minimizing the loss function between predicted and actual values. During training, the batch size is set to 64, meaning 64 samples are processed simultaneously for gradient calculation and parameter updates, achieving a balance between computational efficiency and gradient stability. The initial learning rate is set to 0.001, and a cosine annealing decay strategy is used to gradually reduce the learning rate to 1×10⁻⁶. -5 The optimizer chosen is Adam, with a first-order moment decay coefficient β1 = 0.9, a second-order moment decay coefficient β2 = 0.999, and a numerical stability constant ε = 1 × 10⁻⁶. -8 The training iterations are set to 200, with each iteration iterating through all training samples once. Training progress is monitored using the prediction error on the validation set. An early stopping strategy is used to terminate the training iterations: after each training iteration, the loss value on the validation set is calculated. If the validation loss does not decrease for 15 consecutive iterations, training is terminated and the network parameters are restored to their minimum value to prevent overfitting.

[0068] The loss function is designed to address the stratified anoxic characteristics of tidal river networks. Dissolved oxygen deficit in tidal river networks often exhibits a significant difference between the surface and bottom layers—the surface layer experiences smaller deficits due to oxygen production from photosynthesis and atmospheric reoxygenation, while the bottom layer suffers larger deficits due to oxygen consumption by organic matter and the stratification of water hindering vertical exchange. The loss function adds a stratification-weighted term and a temporal gradient penalty term to the standard mean square error. L_total = L_MSE + α·L_layer + β·L_grad Wherein, L_MSE is the mean square error between the predicted dissolved oxygen deficit value and the actual value; L_layer is the layered weighted loss, which divides the water column into three layers according to the water depth of the target site: surface (0-2m), middle layer (2m to 2m above the bottom), and bottom layer (within 2m of the bottom). When the water depth is less than 4m, the middle layer is canceled, and only the surface layer (0-2m) and the bottom layer (2m to the bottom) are divided. At this time, the original weight of the middle layer of 1.0 is distributed to the surface layer and the bottom layer according to the weight ratio of the surface layer and the bottom layer (0.5:1.5), that is, the weight of the surface layer is adjusted to 0.75 and the weight of the bottom layer is adjusted to 1.75; each layer is assigned a different weight (surface layer weight 0.5, middle layer weight 1.0, bottom layer weight 1.5) to enhance the model's ability to fit the hypoxia event at the bottom layer; L_grad is the temporal gradient penalty term, which calculates the mean square error between the adjacent time step difference values ​​of the predicted sequence and the adjacent time step difference values ​​of the actual sequence, constraining the time change rate of the predicted sequence to be consistent with the actual sequence, and avoiding non-physical drastic fluctuations in the prediction results. α and β are balance coefficients, set to 0.3 and 0.1 respectively.

[0069] Specifically, a loss function is constructed based on the prediction error between the predicted dissolved oxygen deficit at each training prediction time and the actual value at the corresponding time in the dissolved oxygen deficit sequence described in step S2. The gradient of the loss function relative to the temporal prediction network parameters and the learnable parameters used to generate the non-negative scaling coefficients is calculated through backpropagation, and the network parameters and learnable parameters are jointly updated based on the gradient. Based on the updated network parameters and learnable parameters, the corrected attention score at each training prediction time is recalculated according to the above method and normalized by the Softmax function to obtain the temporal attention weight distribution, which is used as the attention posterior probability distribution. After training convergence, the network parameters and learnable parameters are fixed, and forward propagation is performed on the validation set, test set, or sample to be analyzed to obtain the attention posterior probability distribution at each prediction time.

[0070] During training, the scaling factor of the physical prior bias is used as part of the network parameters. The gradient of the loss function with respect to the scaling factor is calculated using the chain rule, and the scaling factor is updated synchronously according to the gradient descent update rule: a negative gradient drives the scaling factor to increase, and a positive gradient drives it to decrease. When the prediction error is large, the gradient calculated by backpropagation indicates how to adjust the scaling factors of the dot product score and the physical bias to reduce the error. If the error is mainly caused by errors guided by the physical bias, the gradient will cause the scaling factor to decrease, reducing the weight of the physical bias and allowing the dot product score to dominate the attention distribution. If the error stems from the dot product score failing to capture key historical moments, the gradient will adjust the transformation parameters of the query and key, improve the calculation of the dot product score, and may increase the scaling factor to strengthen the supplementary role of the physical bias. This adaptive adjustment mechanism enables the method to maintain robust performance under different water conditions and hydrological conditions.

[0071] After training is completed, the trained dissolved oxygen prediction model is obtained. All parameters of the encoder, decoder, two-stage attention weight optimization module, and scaling factor of physical prior bias are fixed and will no longer be updated.

[0072] Step S6: Input the hydrological and meteorological observation sequence, water quality monitoring sequence, dynamic field simulation results, hydrological and meteorological forecast data covering the prediction period, and three-dimensional hydrodynamic model forecast results, as well as the physical prior bias vectors corresponding to each future prediction time, into the trained time series prediction network, and output the predicted dissolved oxygen deficit value for the future time.

[0073] In this embodiment, the prediction phase is based on the trained time-series prediction network. The latest multi-source time-series data up to the current moment is acquired. This data includes hydrological and meteorological observation sequences, water quality monitoring sequences, and dynamic field simulation sequences, and its source is consistent with the data acquisition method described in step S1. The latest multi-source time-series data is preprocessed using the same spatiotemporal alignment, anomaly removal, and standardization methods as in the training phase. Spatiotemporal alignment uses the timestamp of the water quality monitoring sequence as a reference, and linear interpolation resampling is performed on the hydrological and meteorological observation sequences and dynamic field simulation sequences. Anomaly removal uses a sliding window-based statistical test method to remove outlier data points that deviate from the mean by more than three standard deviations. Standardization uses the Z-score method, converting the data into a standard normal distribution using the mean and standard deviation of each input feature saved in the training phase.

[0074] Obtain the hydrological and meteorological observation sequence and water quality monitoring sequence up to the current moment, as well as the hydrological and meteorological forecast data and three-dimensional hydrodynamic model forecast results covering the forecast period; form a historical window input feature matrix based on the hydrological and meteorological observation sequence, water quality monitoring sequence and target station dynamic characteristics extracted from the dynamic field simulation results; and form a future forecast feature sequence based on the hydrological and meteorological forecast data and three-dimensional hydrodynamic model forecast results.

[0075] The preprocessed data is truncated according to the historical time window length to obtain the historical window input feature matrix for the current moment. The number of rows in this matrix equals the historical time window length, and the number of columns equals the input feature dimension. Each column contains features such as rainfall, temperature, wind speed, wind direction, and river flow from hydrological and meteorological observation sequences, auxiliary water quality parameters from water quality monitoring sequences, and flow velocity, water level, and bottom shear stress from dynamic field simulation sequences. The historical time window length is consistent with that used in the training phase.

[0076] For each future prediction time, a physical prior bias vector is generated according to step S4, and arranged in order of future prediction time to form a physical prior bias matrix. The rows of the physical prior bias matrix correspond to the future prediction time, the columns correspond to each historical time in the historical time window, and the matrix elements are the physical prior bias values ​​for the corresponding historical time. The tidal phase and runoff conditions corresponding to the current prediction time are determined, and the corresponding probability density functions are extracted from the dynamic transport time distribution feature library. The current prediction time is the future time from which the dissolved oxygen deficit prediction value needs to be output. Its tidal phase is determined based on the tidal level forecast data of the target water area or the water level sequence output by the three-dimensional fluid dynamics model, and the runoff conditions are determined by the future runoff sequence. Using the current prediction time as a reference, the time difference between each historical time in the historical time window and the current prediction time is used as the transport time interval. The probability density value corresponding to each historical time is queried from the extracted probability density function. If the transport time interval exceeds the transport time coverage range of the dynamic transport time distribution feature library, the probability density value is set to zero. All probability density values ​​are arranged in chronological order to generate a physical prior bias vector of the same length as the historical time window.

[0077] Future runoff sequences are obtained from at least one of operational hydrological forecasts, publicly available hydrological ensemble forecasts, and rainfall-runoff forecasts generated from meteorological forecasts. Runoff volumes are primarily obtained from stations upstream of the tidal boundary to avoid the influence of tidal back-and-forth flow on the flow. Spatial matching and bias correction based on historical measured flow are performed on the obtained runoff forecasts. Based on historical forecasting skills for different lead times, the corrected runoff forecasts are adaptively fused with multi-year daily runoff climatology. When external runoff forecasts are missing or their historical forecasting skills are below a set threshold, future runoff sequences are generated using multi-year daily runoff climatology or high-water, normal-water, and low-water scenarios.

[0078] Specifically, the future runoff sequence can be calculated using "bias-corrected runoff forecasts, daily runoff climatology, and fusion weights that vary with the forecast period," as shown in the following formula: (t+h)=Q_clim(d)+w(h)×[Q_corr(t+h)-Q_clim(d)] In the formula: (t+h) represents the runoff forecast value at the current time t to the future time t+h; Q_clim(d) represents the multi-year daily runoff climatology corresponding to the target date d; Q_corr(t+h) represents the future runoff forecast value after correction for historical measured flow deviations at the local hydrological station; w(h) represents the runoff forecast fusion weight when the lead time is h. This weight is determined based on historical return experiments or historical forecasting techniques and gradually decreases as the lead time increases; h represents the lead time of the runoff forecast.

[0079] The fusion weights corresponding to different forecast periods can be determined as follows: (1) When the forecast period is 1 to 3 days, the reliability of the runoff forecast is relatively high, and the fusion weight w(h) can be 0.8 to 1.0, so that the future runoff sequence mainly retains the runoff forecast information after deviation correction. (2) When the forecast period is 4 to 10 days, the fusion weight is gradually reduced according to the root mean square error, Nash efficiency coefficient, correlation coefficient or other forecasting skill indicators of the runoff forecast in the historical return experiment. (3) When the forecast period is 11 to 15 days, the proportion of multi-year daily runoff climatology is appropriately increased, while retaining the abnormal change information of the runoff forecast relative to the daily runoff climatology. (4) When the forecast period exceeds 15 days, the sub-seasonal ensemble runoff forecast is adopted, or the high-water, normal-water and low-water scenarios are constructed according to the multi-year daily runoff climatology to characterize the runoff changes and their uncertainties within a longer forecast period.

[0080] Considering that abnormal runoff processes and the number of historical samples may affect the multi-year diurnal runoff climatology, the multi-year diurnal runoff climatology is preferably calculated using a sliding time window. Specifically, with the target date as the center, a sliding time window of 7 to 15 days before and after the target date is selected, the measured average daily flow within the same window over many years is extracted, and the diurnal runoff climatology for the corresponding target date is calculated based on the extracted measured average daily flow.

[0081] The historical window input feature matrix, future forecast feature sequence, and physical prior bias matrix are input into the trained temporal prediction network. The physical prior bias vector corresponding to each future prediction time is called, and the dissolved oxygen deficit prediction value for each future prediction time is output. During forward propagation, the temporal prediction network synchronously outputs the attention posterior probability distribution for each prediction time, which is used to execute the dominant process after step S6 and to identify the low-oxygen risk time window. The encoder processes the historical window input feature matrix step-by-step. At each time step, the multi-source input features are adaptively weighted through a two-stage spatial attention weight optimization module, and the historical hidden state sequence is output after being updated by two layers of long short-term memory network units. At each prediction time step, the decoder uses the encoder's hidden state from the last time step as the initial decoded hidden state, obtains the decoder's hidden state at the current prediction time step as the query vector, and obtains the entire sequence of historical hidden states output by the encoder as the key vector. The dot product of the query vector and each key vector is calculated to obtain the dot product attention score. Each element of the physical prior bias vector is multiplied by the scaling factor optimized during training, and then added to the corresponding historical time step's dot product attention score to obtain the corrected attention score. All historical time step's corrected attention scores are normalized using the Softmax function to obtain the attention weight distribution. Based on this attention weight distribution, the encoder's historical hidden state sequence is weighted and aggregated to generate temporal context features. These temporal context features are concatenated with the decoder's current hidden state and mapped through a fully connected layer before being input into a Long Short-Term Memory (LSTM) network unit to update the decoder's hidden state. Finally, the output layer is fully connected and mapped to the standardized dissolved oxygen deficit prediction value for the current prediction time step. The decoder uses the prediction value from the previous time step as the input for the next time step, and autoregressively generates standardized dissolved oxygen deficit prediction values ​​for multiple future time steps.

[0082] The standardized dissolved oxygen deficit prediction values ​​are destandardized based on the mean and standard deviation of the dissolved oxygen deficit sequence saved during the training phase, restoring them to predicted values ​​with physical dimensions, and then output. These predictions can be directly used for water quality early warning and management decisions; an early warning is triggered when the predicted dissolved oxygen deficit exceeds a preset threshold. In practical applications, a rolling prediction method is used: each time new measured data is obtained, it is added to the historical sequence, the oldest data in the time window is removed, the historical time window length remains unchanged, and the above prediction process is repeated to generate updated future time-based prediction values.

[0083] Step S6 and subsequent steps include identifying the dominant process and the time window for low oxygen risk.

[0084] In this embodiment, after training convergence, the attention weight distribution generated by the decoder at each prediction time step integrates physical prior guidance and data-driven correction. This distribution reflects the network's judgment on the importance of historical moments based on all available information and is used as the posterior probability distribution of attention. The physical prior bias vector characterizes the contribution strength of each historical moment to the current prediction under pure hydrodynamic driving. The two constitute the basic data pair for comparative analysis.

[0085] For each prediction time, the arrival probabilities of historical time conditions used when generating the corresponding physical prior bias vector are normalized to obtain the physical prior probability distribution, and the attention posterior probability distribution for the same prediction time is obtained. The Pearson correlation coefficient and peak position deviation between the physical prior probability distribution and the attention posterior probability distribution are calculated. The Pearson correlation coefficient is obtained by calculating the ratio of the product of the covariance and standard deviation of two equal-length vectors, ranging from -1 to +1. The closer the coefficient is to +1, the more consistent the two distributions are, and the higher the agreement between the key historical time the network focuses on and the arrival time predicted by the physical transport model. The peak position deviation is calculated by matching the time index corresponding to the maximum weight in the attention posterior probability distribution with the time index corresponding to the maximum weight in the physical prior probability distribution. The time difference between the two is calculated, and this deviation is expressed in hours, quantifying the degree of deviation of the network's judgment of key influencing times relative to the physical transport prediction.

[0086] When the correlation coefficient is greater than the first threshold and the peak position deviation is less than the second threshold, it is determined to be a period dominated by physical transport; when the correlation coefficient is less than the third threshold or the peak position deviation is greater than the fourth threshold, it is determined to be a period dominated by local biochemical processes; all other cases are determined to be transitional periods. Before determining the thresholds, the daily prediction results for 730 days of training and validation sets were statistically analyzed: the correlation coefficient distribution of the physical transport-dominated period samples had a mean of 0.82 and a standard deviation of 0.08, while that of the local biochemical-dominated period samples had a mean of 0.28 and a standard deviation of 0.10. Using 0.65 and 0.40 as the cutoff points for correlation coefficients, and 6 hours and 12 hours as the cutoff points for peak deviations, ROC curve analysis showed that this set of thresholds minimized the sum of the misclassification rates for the two dominant periods (approximately 12%). Quantitative determination adopted a threshold comparison mechanism: the first and third thresholds were used to divide the intervals of correlation coefficients, and the second and fourth thresholds were used to divide the intervals of peak deviations. The first threshold was greater than the third threshold, and the second threshold was less than the fourth threshold. The specific values ​​of each threshold were determined based on the statistical characteristics of historical observation data of the target water area.

[0087] The physical transport-dominated period refers to a time when dissolved oxygen depletion is primarily caused by upstream pollutants transported to the target site via water flow. During this period, water quality deterioration at the target site is closely related to upstream discharges or changes in hydrological conditions, with the time lag determined by the transport time. The local biochemical process-dominated period refers to a time when dissolved oxygen depletion is mainly driven by local biochemical processes at the target site, such as organic matter degradation, phytoplankton respiration, or sediment oxygen consumption. During this period, water quality evolution is controlled by local environmental factors such as light and temperature, and its correlation with upstream transport is weaker. Accurately identifying the dominant mechanisms is crucial for developing targeted water quality management strategies.

[0088] Interpretable attribution methods were employed to analyze the contribution of each input feature to the dissolved oxygen deficit prediction value. The stage-specific contributions of exogenous oxygen-consuming substance transport processes and local biochemical oxygen-consuming processes were identified. Combined with the dominant processes and dissolved oxygen deficit prediction results over consecutive future prediction periods, the hypoxia risk time window under high runoff conditions was determined. Interpretable attribution methods could employ gradient-based methods such as the integral gradient method, or perturbation-based methods such as the Shapley value method, to calculate the contribution of each input feature to the prediction results. Analysis of these contributions identified which features (such as upstream chemical oxygen demand, ammonia nitrogen concentration, and other exogenous oxygen-consuming substance indicators, or local water temperature, chlorophyll concentration, and other biochemical process indicators) played a dominant role in dissolved oxygen deficit prediction during specific time periods, thereby further validating and refining the identification results of the dominant mechanisms.

[0089] In the Pearl River Estuary case study, during the summer months of July to September, water temperatures reached over 30 degrees Celsius, leading to a surge in phytoplankton growth and a significant increase in dissolved oxygen depletion during the early morning hours due to nighttime respiration. Comparison revealed that the peak of the posterior distribution of attention was concentrated in historical moments six to twelve hours prior, corresponding to the evening to nighttime of the previous day. These moments coincide with the period when respiration begins to dominate after photosynthesis ceases. Furthermore, the posterior weights of transport moments 24 and 48 hours prior, indicated by physical prior bias, were significantly lower than the physical prior weights, with the correlation coefficient dropping to 0.4. This period was determined to be dominated by local biochemical processes, and water quality management should shift its focus to controlling local nutrient loading and regulating phytoplankton biomass.

[0090] Conversely, during the high-water season from March to May, strong runoff carries large amounts of pollutants from upstream and rapidly discharges them. The posterior distribution of attention closely matches the physical prior, with peak values ​​occurring 20 to 30 hours prior, showing a correlation coefficient of 0.85. This period is identified as the physical transport-dominated phase, and the outbreak of dissolved oxygen depletion is directly related to the heavy rainfall and pollutant runoff two days prior upstream. At this time, water quality management requires strengthened monitoring of upstream pollution sources and runoff regulation, mitigating downstream transport loads through source reduction and reservoir flooding measures.

[0091] By continuously monitoring and periodically comparing data, a spatiotemporal distribution map of the dominant mechanisms throughout the year can be created. This map reveals that the period dominated by physical transport is mainly concentrated during the spring tides of the high-water season, lasting approximately 30% of the year. The period dominated by local biochemical processes mainly occurs during the neap tides of the low-water season and the high-temperature periods of summer, lasting approximately 40% of the year. The remaining 30% is a transitional period, in which physical transport and local biochemical processes work together and are of equal importance. This mechanism division provides a dynamic zoning basis for water quality prediction and management, allowing for differentiated prediction strategies and management measures in different periods, avoiding a one-size-fits-all approach to extensive management, and improving the accuracy and effectiveness of water quality protection.

[0092] To verify the effectiveness of the technical solution of this invention, a comparative experiment was conducted in a typical tidal river network in a tidal estuary. The experiment used observation data from three consecutive years for model training and testing (training in the first two years and testing in the third year), with the target being the dissolved oxygen deficit value for the next 24 hours.

[0093] The comparison objects include: (1) a pure three-dimensional baroclinic mechanism model (based on the environmental fluid dynamics model EFDC, relying only on physical equations for simulation, without data-driven modules); (2) a general long short-term memory network model (standard two-layer LSTM, without attention mechanism, without physical prior); (3) a common fusion scheme (directly using the dynamic field simulation sequence as an additional input feature to the LSTM input layer, without physical prior bias embedding attention mechanism); (4) the scheme of this invention (a time-series prediction network containing a physical prior bias embedding two-stage attention mechanism).

[0094] Evaluation metrics include root mean square error (RMSE, mg / L), mean absolute error (MAE, mg / L), coefficient of determination (R²), and time taken for a single prediction calculation (seconds).

[0095] The RMSE of the proposed solution is 0.62 mg / L, which is approximately 66.5% lower than the pure three-dimensional baroclinic mechanism model (RMSE 1.85 mg / L), approximately 48.3% lower than the general long short-term memory network model (RMSE 1.20 mg / L), and approximately 34.7% lower than the ordinary fusion solution (RMSE 0.95 mg / L). 2 The result of 0.91 indicates that the model can explain more than 90% of the variation in dissolved oxygen deficit. The computation time was 2.0 seconds, which is slightly longer than that of the general LSTM (0.8 seconds) and the ordinary fusion scheme (1.2 seconds), but much shorter than that of the pure mechanistic model (3600 seconds, 1 hour), thus meeting the timeliness requirements of operational rolling forecasts.

[0096] Regarding the accuracy of stratified prediction, the RMSE of the proposed scheme for predicting dissolved oxygen deficit in the surface layer is 0.48 mg / L, in the middle layer it is 0.58 mg / L, and in the bottom layer it is 0.71 mg / L, with an overall weighted RMSE of 0.62 mg / L. In comparison, the overall RMSE of the general long short-term memory network model, after weighting each layer by sample size, is 1.20 mg / L, and that of the ordinary fusion scheme is 0.95 mg / L, consistent with the overall RMSE of the proposed scheme. During periods when the tidal flow has a significant impact (during spring tides), the RMSE of the proposed scheme is 0.58 mg / L, while the RMSE of the ordinary fusion scheme rises to 1.12 mg / L during this period, indicating that the guiding effect of physical prior bias on the transport lag of the tidal flow is particularly significant during periods of complex tidal dynamics.

[0097] Regarding the diagnosis of the dominant mechanism, the daily prediction results for the entire third test year were statistically analyzed: During the spring high tide period from March to May, the Pearson correlation coefficient between the posterior distribution of attention and the physical prior reached 0.85, with a peak position deviation of less than 4 hours, indicating a period dominated by physical transport, accounting for approximately 30% of the total days of the year; during the high-temperature period from July to September in summer, the correlation coefficient dropped to below 0.35, with a peak deviation exceeding 14 hours, indicating a period dominated by local biochemical processes, accounting for approximately 40% of the total days of the year; the remaining periods were transitional periods, accounting for approximately 30%. These diagnostic results are consistent with the concurrent hydrological, meteorological, and ecological observation data of a certain tidal estuary, verifying the effectiveness of the diagnostic mechanism of this invention.

[0098] Figure 2 A radar chart comparing the normalized overall performance of the present invention's solution and a conventional fusion solution, as provided in the embodiments of the present invention. Figure 2 Scores for each dimension were calculated using a linear normalization method, mapping each scheme's best value to 1 and worst value to 0 for that dimension. The computational efficiency dimension was inversely normalized with a baseline of 0 points for a pure mechanistic model (3600 seconds) and 1 point for a general LSTM (0.8 seconds) (shorter time, higher score). Therefore, the present invention scored 0.90 for 2.0 seconds, while the ordinary fusion scheme scored 1.0 for 1.2 seconds. The radar chart evaluated five dimensions: fitting ability, computational efficiency, tidal adaptability, interpretability, and prediction accuracy. The normalized scores for each dimension ranged from 0 to 1, with higher scores indicating better performance. Solid polygons represent the present invention's scheme, and dashed polygons represent the ordinary fusion scheme. Figure 2As can be seen, the normalized scores of the proposed solution are higher than those of the ordinary fusion solution in terms of fitting ability, tidal adaptability, interpretability, and prediction accuracy, while slightly lower in terms of computational efficiency. Specifically, the normalized scores of the proposed solution in terms of fitting ability, computational efficiency, tidal adaptability, interpretability, and prediction accuracy are approximately 0.93, 0.90, 1.00, 1.00, and 0.82, respectively; while the normalized scores of the ordinary fusion solution in the corresponding dimensions are approximately 0.70, 1.00, 0.65, 0.60, and 0.55, respectively. The difference in computational efficiency between the two solutions is relatively small, while the proposed solution shows a more significant advantage in fitting ability, tidal adaptability, interpretability, and prediction accuracy. The ordinary fusion solution uses the simulated dynamic field variables as ordinary input features and concatenates them with other monitoring variables. Although it can utilize some hydrodynamic information, its ability to represent the long time delays, repeated arrivals, and multi-peak transport time distribution caused by tidal reciprocating flow is relatively limited, and it is also difficult to directly reveal the correspondence between the historical moments of interest in the model and the water mass transport process. The present invention converts the dynamic transport time probability distribution obtained by particle tracking into a physical prior bias and embeds it into the attention scoring function of the temporal prediction network. This enables the network to prioritize historical moments with higher physical arrival probabilities, thus achieving higher evaluation in terms of tidal condition adaptability and predictive process interpretability, while also improving fitting ability and prediction accuracy. Figure 2 The comparison results show that, compared with ordinary feature splicing fusion, the physical prior-guided attention mechanism adopted in this invention has better overall performance.

[0099] In summary, this embodiment verifies the dissolved oxygen prediction method based on time-series and dynamic process guidance proposed in this invention, using continuous monitoring data of tidal river networks and three-dimensional fluid dynamics simulation data. Experimental results show that the proposed scheme can generate a physical prior bias using the transport time probability distribution obtained from dynamic field simulation and embed it into the attention scoring function of the time-series prediction network, achieving a deep integration of physical transport laws and data-driven learning. The dissolved oxygen deficit value predicted by the proposed scheme can track the overall trend of measured values ​​well. The main peaks, troughs, and periodic changes of the predicted curve and the measured curve have a good temporal correspondence, indicating that the proposed scheme can maintain good continuous prediction and trend tracking capabilities. Compared with conventional fusion schemes that use dynamic field simulation sequences as ordinary input features, the proposed solution exhibits superior overall performance in terms of fitting ability, prediction accuracy, tidal adaptability, and predictive interpretability. In particular, it is more effectively able to adapt to the complex transport delays caused by tidal reciprocating flows and can identify the dominant periods of physical transport and local biochemical processes based on the consistency between physical priors and attention posteriors. This provides more reliable technical support for water quality prediction and anoxic causes analysis in tidal river networks. The above implementation results demonstrate that the proposed solution can balance prediction accuracy, physical consistency, and predictive interpretability, making it suitable for dissolved oxygen rolling prediction, water quality early warning, and refined water environment management in complex hydrodynamic environments such as tidal river networks.

[0100] The above description is merely a specific embodiment of this specification. Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the specific working processes of the systems, modules, and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here. It should be understood that the scope of protection of this specification is not limited thereto. Any person skilled in the art can easily conceive of various equivalent modifications or substitutions within the technical scope disclosed in this specification, and these modifications or substitutions should all be covered within the scope of protection of this specification.

Claims

1. A method for predicting dissolved oxygen in tidal river networks based on time-series and dynamic processes, characterized in that, include: Step S1: Obtain the hydrological and meteorological observation sequences, water quality monitoring sequences, and multi-parameter dynamic field simulation results output by the three-dimensional hydrodynamic model from the target station; Step S2: Calculate the theoretical saturated dissolved oxygen concentration and dissolved oxygen deficit sequence based on the water quality monitoring sequence and the hydrological and meteorological observation sequence; Step S3: Based on the three-dimensional velocity field tracking neutral tracer particles in the dynamic field simulation results, construct a dynamic transport time distribution feature library corresponding to tidal phase, tidal pattern, and runoff conditions; Step S4: For each future predicted time, determine the dynamic state based on hydrological and meteorological forecast data and the forecast results of the three-dimensional hydrodynamic model, and combine the transport time interval from each historical time to the future predicted time to obtain the conditional arrival probability from the dynamic transport time distribution feature library and generate the corresponding physical prior bias vector. Step S5: Embed the physical prior bias obtained by logarithmic transformation of the conditional arrival probability into the attention scoring function of the time series prediction network. Take the hydro-meteorological observation sequence, water quality monitoring sequence, dynamic field simulation results and hydro-meteorological forecast features and dynamic forecast features corresponding to the prediction period as input, and take the dissolved oxygen deficit sequence as the supervision target. Train the network parameters and the non-negative scaling coefficient used to adjust the influence intensity of the physical prior bias through backpropagation. Step S6: Input the hydrological and meteorological observation sequence, water quality monitoring sequence, dynamic field simulation results, hydrological and meteorological forecast data covering the prediction period, and three-dimensional hydrodynamic model forecast results, as well as the physical prior bias vectors corresponding to each future prediction time, into the trained time series prediction network, and output the predicted dissolved oxygen deficit value for the future time.

2. The method for predicting dissolved oxygen in tidal river networks based on time-series and dynamic processes according to claim 1, characterized in that, In step S1, the hydro-meteorological observation sequence is collected by automatic weather stations and hydrological stations, including rainfall, air temperature, wind speed, wind direction, river flow, and atmospheric pressure; the water quality monitoring sequence is collected by water quality monitoring buoy equipment, including dissolved oxygen concentration, water temperature, salinity, chemical oxygen demand, ammonia nitrogen concentration, and chlorophyll concentration; the three-dimensional hydrodynamic model is a three-dimensional baroclinic dynamic model, and the dynamic field simulation results include the flow velocity vector, water level, temperature, salinity, and bottom shear stress of each grid point over time. The tidal boundary conditions, runoff boundary conditions, and tidal level change sequence of the target water area of ​​the three-dimensional hydrodynamic model are obtained to determine the release time of neutral tracer particles, the transport period of neutral tracer particles, and the dynamic state corresponding to the future prediction time. Using the timestamps of the water quality monitoring sequence as a reference, the hydrological and meteorological observation sequence and the dynamic field simulation results are time-aligned, and the dynamic field simulation results are spatially interpolated to extract the dynamic characteristics of the target station, forming training input data with a unified time interval.

3. The method for predicting dissolved oxygen in tidal river networks based on time-series and dynamic processes according to claim 1, characterized in that, Step S2 specifically includes: Water temperature, salinity, and actual dissolved oxygen concentration at each time point are extracted from the water quality monitoring sequence, and the corresponding atmospheric pressure is extracted from the hydro-meteorological observation sequence. The theoretical saturated dissolved oxygen concentration is calculated based on the water temperature and salinity, and the theoretical saturated dissolved oxygen concentration is corrected for atmospheric pressure based on the atmospheric pressure. The difference between the corrected theoretical saturated dissolved oxygen concentration and the actual dissolved oxygen concentration is calculated to obtain the dissolved oxygen deficit value at each time point, and a dissolved oxygen deficit sequence is formed. The theoretical saturated dissolved oxygen concentration was calculated using the APHA formula under low salinity conditions and the Weiss formula under seawater salinity conditions.

4. The method for predicting dissolved oxygen in tidal river networks based on time-series and dynamic processes according to claim 1, characterized in that, Step S3 specifically includes: The three-dimensional hydrodynamic model is driven by atmospheric forcing, tidal and runoff conditions; Neutral tracer particles are released at the upstream boundary section of the three-dimensional hydrodynamic model, and the particle positions are updated step by step based on the three-dimensional velocity field in the dynamic field simulation results. Record the timestamp of the particle's first entry into the target section of the target site and the repeated arrival timestamp of the particle's re-entry into the target section after leaving. Obtain particle transport time samples based on the difference between each arrival timestamp and the corresponding release timestamp. The tidal phase is determined based on the tidal level change at the moment of particle release, the tidal pattern is determined based on the tidal range, and the runoff conditions are determined based on the river flow during the moment of particle release and the corresponding transport period. Kernel density estimation is performed on particle transport time samples under the same dynamic state to obtain the transport time condition arrival probability density function, which is then associated with and stored with the corresponding tidal phase, tidal pattern, and runoff conditions to form the dynamic transport time distribution feature library.

5. The method for predicting dissolved oxygen in tidal river networks based on time-series and dynamic processes according to claim 4, characterized in that, Kernel density estimation is performed on particle transport time samples under the same dynamic state to obtain the transport time conditional arrival probability density function, which is then associated with and stored in relation to the corresponding tidal phase, tidal type, and runoff conditions to form the dynamic transport time distribution feature library, specifically including: Kernel density estimation is performed on particle transport time samples under the same dynamic state to obtain the transport time conditional arrival probability density function; Based on the arrival probability density function of the transport time condition, the expected transport time, transport dispersion and reciprocating flow return rate are calculated as auxiliary transport dynamics characteristics. These are then associated with the corresponding tidal phase, tidal pattern and runoff conditions and stored in the dynamic transport time distribution feature library for feature library quality assessment and index optimization. Wherein, the expected transport time is the first moment of the transport time conditional arrival probability density function, which is used to characterize the statistical average time for the water mass to reach the target section from the release location; The degree of transport dispersion is the standard deviation of the arrival probability density function under the transport time condition, which is used to characterize the variability of transport time; The reciprocating flow return rate is the proportion of the number of particles that arrive at least once via a reciprocating flow to the total number of particles that arrive for the first time, and is used to characterize the contribution of tidal reciprocating flow to material retention. The expected transport time, transport dispersion, and reciprocating flow return rate are associated with the corresponding tidal phase, tidal pattern, and runoff conditions and stored in the dynamic transport time distribution feature library.

6. The method for predicting dissolved oxygen in tidal river networks based on time-series and dynamic processes according to claim 1, characterized in that, Step S4 specifically includes: For each future prediction time, iterate through each historical time in the historical time window, take each historical time as the hypothetical particle release time, and determine the dynamic state corresponding to each historical time and the transport period up to the future prediction time. The dynamic state includes tidal phase, tidal pattern and runoff conditions. Based on the dynamic state, the corresponding transport time condition arrival probability density function is obtained by matching or interpolating from the dynamic transport time distribution feature library; Calculate the transport time interval between each historical moment and the predicted future moment, and obtain the conditional arrival probability corresponding to each transport time interval; The probability of each condition being reached is added to a preset positive number, and the natural logarithm is taken. The logarithm is then arranged in the order of historical time moments to generate a physical prior bias vector of the same length as the historical time window. A corresponding physical prior bias vector is generated for each future prediction time moment.

7. The method for predicting dissolved oxygen in tidal river networks based on time-series and dynamic processes according to claim 1, characterized in that, Step S5 specifically includes: The time series prediction network includes an encoder and a decoder based on an attention mechanism. The encoder is used to extract time series features in historical time windows, and the decoder is used to combine time attention weights and input features of future time periods to generate prediction results. The auxiliary water quality parameters in the hydrological and meteorological observation sequence and water quality monitoring sequence, as well as the target station dynamic features extracted from the dynamic field simulation results, are organized into an input feature matrix according to historical time windows. The rows of the input feature matrix correspond to each moment in the historical time window, and the columns correspond to each input feature dimension. The hydrological and meteorological features and target station dynamic features corresponding to the prediction period of the training samples are used as input features for the future period. During the model training phase, the hydrological and meteorological features and target station dynamic features observed synchronously in history are used as input features for the future period, i.e., assuming perfect forecast conditions. The physical prior bias vector is used as a logarithmic prior bias, and its influence intensity is adjusted by a non-negative scaling factor obtained by non-negative mapping of learnable parameters, and then fused with the temporal attention score. Using the dissolved oxygen deficit sequence as the supervision target, the network parameters and the learnable parameters are jointly updated through backpropagation to obtain the attention posterior probability distribution and the trained dissolved oxygen prediction model. Specifically, a loss function is constructed based on the prediction error between the predicted dissolved oxygen deficit at each training prediction time and the actual value at the corresponding time in the dissolved oxygen deficit sequence described in step S2. The gradient of the loss function relative to the temporal prediction network parameters and the learnable parameters used to generate the non-negative scaling coefficients is calculated through backpropagation, and the network parameters and learnable parameters are jointly updated based on the gradient. Based on the updated network parameters and learnable parameters, the corrected attention score at each training prediction time is recalculated according to the above method and normalized by the Softmax function to obtain the temporal attention weight distribution, which is used as the attention posterior probability distribution. After training convergence, the network parameters and learnable parameters are fixed, and forward propagation is performed on the validation set, test set, or sample to be analyzed to obtain the attention posterior probability distribution at each prediction time.

8. The method for predicting dissolved oxygen in tidal river networks based on time-series and dynamic processes according to claim 7, characterized in that, The physical prior bias embedding time attention scoring function specifically includes: Obtain the hidden state of the decoder at the current future prediction time as the query feature, obtain the hidden state of each historical time step output by the encoder as the key feature, and calculate the initial attention score for each historical time step. Obtain the physical prior bias vector corresponding to the current future prediction time, adjust the influence intensity of the physical prior bias vector through the non-negative scaling factor, and fuse the adjusted physical prior bias with the initial attention score of the corresponding historical time step to obtain the corrected attention score; The temporal attention weights are calculated based on the corrected attention scores, and the historical hidden state sequences output by the encoder are weighted and aggregated to generate temporal context features. The time context features and the future time period input features corresponding to the current future prediction time are input into the decoder to obtain the dissolved oxygen deficit prediction value for that future prediction time.

9. The method for predicting dissolved oxygen in tidal river networks based on time-series and dynamic processes according to claim 1, characterized in that, Step S6 specifically includes: Obtain the hydrological and meteorological observation sequence and water quality monitoring sequence up to the current moment, as well as the hydrological and meteorological forecast data and three-dimensional hydrodynamic model forecast results covering the forecast period; form a historical window input feature matrix based on the hydrological and meteorological observation sequence, water quality monitoring sequence and target station dynamic characteristics extracted from the dynamic field simulation results; and form a future forecast feature sequence based on the hydrological and meteorological forecast data and three-dimensional hydrodynamic model forecast results. For each future prediction time, a physical prior bias vector is generated according to step S4, and arranged according to the future prediction time to form a physical prior bias matrix; the rows of the physical prior bias matrix correspond to the future prediction time, the columns correspond to each historical time in the historical time window, and the matrix elements are the physical prior bias values ​​of the corresponding historical time. The historical window input feature matrix, future forecast feature sequence, and physical prior bias matrix are input into the trained time series prediction network. The physical prior bias vector corresponding to each future prediction time is called respectively, and the dissolved oxygen deficit prediction value for each future prediction time is output.

10. The method for predicting dissolved oxygen in tidal river networks based on time-series and dynamic processes according to claim 1, characterized in that, Following step S6, the process also includes the identification of the dominant process and the time window for hypoxia risk: For each prediction time, the arrival probabilities of conditions at each historical time used when generating the corresponding physical prior bias vector are normalized to obtain the physical prior probability distribution, and the attention posterior probability distribution at the same prediction time is obtained. Calculate the Pearson correlation coefficient and peak position deviation between the physical prior probability distribution and the attention posterior probability distribution; When the correlation coefficient is greater than the first threshold and the peak position deviation is less than the second threshold, it is determined to be the physical transport-dominated period; when the correlation coefficient is less than the third threshold or the peak position deviation is greater than the fourth threshold, it is determined to be the local biochemical process-dominated period; all other cases are determined to be the transition period. An interpretable attribution method was used to analyze the contribution of each input feature to the dissolved oxygen deficit prediction value, identify the stage contributions of exogenous oxygen-consuming substance transport processes and local biochemical oxygen-consuming processes, and combine the dominant processes in continuous future prediction periods with the dissolved oxygen deficit prediction results to determine the time window of hypoxia risk under high runoff conditions.

Images