Water quality prediction method, system, device and readable storage medium

Abstract

The application discloses a water quality prediction method, system and device and a readable storage medium, and particularly relates to the field of prediction calculation based on signal decomposition, and is used for solving the problem that the cross-scale time lag coupling relationship is difficult to be stably captured due to the mismatch between the covariant full-band information and the target variable single-band component information structure in the existing water quality prediction; the target variable sequence and the covariant sequence are both subjected to modal decomposition, the false components are removed by using the phase randomization to replace the nonlinear complexity difference test of the data, the stable components are screened and combined into prediction units by means of the double-domain coincidence judgment of the time-domain stable domain and the frequency-domain resonance region, each prediction unit is predicted by the time series prediction network embedded with the self-attention mechanism, and finally the water quality prediction result is obtained by superposition and reconstruction.

Description

Technical Field

[0001] This invention relates to the field of predictive computing technology based on signal decomposition, and more specifically, to water quality prediction methods, systems, devices, and readable storage media. Background Technology

[0002] In nearshore aquaculture, accurate prediction of water quality indicators is of significant reference value for aquaculture management decisions. Water quality parameters, such as dissolved oxygen, are influenced by a combination of environmental factors including temperature, salinity, and chlorophyll content, resulting in monitoring data exhibiting obvious non-stationarity and multi-scale fluctuations. To reduce the modeling difficulty of the original sequences, existing technologies typically employ signal decomposition methods to preprocess the target water quality variable, breaking it down into several relatively uniform frequency sub-sequence components. Then, neural network prediction models are constructed separately for each component, and these components are superimposed and reconstructed. Within this framework, the input to the prediction model, in addition to the sub-sequences of the target variable, usually includes the unprocessed sequence of related covariates.

[0003] In the existing processing methods described above, the covariate sequence directly participates in the predictive modeling of each sub-sequence component in its original form. Since the original covariates contain fluctuation information across the entire frequency band, while a single sub-sequence component of the target variable represents a local feature of a specific frequency band, their information structures do not match. This means that when the model predicts a component in a certain frequency band, it needs to manually distinguish the driving components actually relevant to that frequency band from the mixed full-frequency signal of the covariates, while fluctuation information unrelated to that frequency band constitutes interference. When there is a cross-scale time-delay coupling relationship between the covariates and the target variable—that is, when different target components in different frequency bands correspond to different key covariates and their lag times—the above-mentioned information mixing problem is further exacerbated, making it difficult to stably capture the co-evolutionary patterns between variables and limiting the accuracy of the prediction model in characterizing dynamic changes in water quality. Summary of the Invention

[0004] In order to overcome the above-mentioned defects of the prior art, the present invention provides a water quality prediction method, system, device and readable storage medium to solve the problems mentioned in the background art.

[0005] Water quality prediction methods include: S1: Obtain the target variable sequence and covariate sequence of the water area to be predicted; S2: Perform mode decomposition on the target variable sequence and the covariate sequence respectively to obtain the target variable mode components and the covariate mode components; S3: A significance test is performed by comparing the nonlinear complexity differences between the target variable modal components and covariate modal components and their respective corresponding phase randomized substitute data. Components that fail the test are removed, and stable target variable modal components and stable covariate modal components are retained. S4: Combine the stable target variable modal components and stable covariate modal components whose overlap between the time domain stable domain and the frequency domain resonance region exceeds the overlap threshold into a prediction unit. The time domain stable domain is determined by the instantaneous coupling strength of the sliding window, and the frequency domain resonance region is determined by the cross wavelet power spectrum. S5: By embedding a time series prediction network with a self-attention mechanism, the predicted values ​​of the stable target variable mode components are obtained by taking the stable covariate mode components in the prediction unit as input. S6: Superimpose and reconstruct the predicted values ​​of the modal components of all stable target variables to obtain the water quality prediction results.

[0006] Furthermore, the target variable sequence and covariate sequence of the water area to be predicted are obtained, including: Obtain the location information of monitoring buoys in the water area to be predicted, and determine the target monitoring station in the water area to be predicted based on the location information; Collect historical monitoring data recorded at the target monitoring stations. The historical monitoring data includes the original target variable sequence and the original covariate sequence. Calculate the kurtosis values ​​of the original target variable sequence and the original covariate sequence within the sliding window, and remove outliers whose kurtosis values ​​exceed the kurtosis threshold; Iterative principal component analysis was used to fill missing values ​​in the sequences after outlier removal to obtain the filled sequences. The filled sequence is standardized to obtain the target variable sequence and the covariate sequence.

[0007] Furthermore, mode decomposition is performed on the target variable sequence and the covariate sequence respectively to obtain the target variable mode components and covariate mode components, including: The target variable sequence is divided into multiple target variable sub-segments by identifying extreme points and determining the segmentation boundary based on the local density of the extreme point distribution. For each target variable segment, calculate the power spectral density of the target variable segment, and adaptively set the number of decomposition layers of the target variable segment according to the number of main peaks of the power spectral density. Adaptive noise complete set empirical mode decomposition is performed on each target variable sub-segment to obtain the sub-mode components corresponding to each target variable sub-segment; The submodal components of all target variable segments are concatenated in chronological order to obtain the target variable modal components; The covariate sequence is processed using the same method as the target variable sequence to obtain the covariate modal components.

[0008] Furthermore, a significance test is performed by comparing the nonlinear complexity differences between the target variable modal components and covariate modal components and their corresponding phase randomized substitute data. Components that fail the test are removed, and stable target variable modal components and stable covariate modal components are retained, including: Fourier transforms are performed on the target variable modal components and covariate modal components respectively, keeping the amplitude spectrum unchanged. The phase spectrum is randomly shuffled and then inverse Fourier transform is performed. This process is repeated multiple times to generate a corresponding set of phase randomization replacement data. Calculate the permutation entropy of the modal components of the target variable and the modal components of the covariate, as well as the mean and standard deviation of the permutation entropy of the corresponding phase randomized substitute data; For the modal components of the target variable, if their permutation entropy falls within the range of the permutation entropy mean plus or minus a preset multiple of the permutation entropy standard deviation of the corresponding phase randomized substitute data, then it is determined that the significance test has not been passed and the variable is removed; otherwise, it is determined that the significance test has been passed and the variable is retained. The same significance test and elimination / retention method were used for the covariate modal components as for the target variable modal components; The target variable modal components that pass the significance test are taken as stable target variable modal components, and the covariate modal components that pass the significance test are taken as stable covariate modal components.

[0009] Furthermore, the stable target variable modal components and stable covariate modal components whose overlap between the time-domain stable domain and the frequency-domain resonance region exceeds the overlap threshold are combined into a prediction unit, including: For each stable target variable mode component and each stable covariate mode component, a sliding window is used to extract the synchronization segments of the stable covariate mode component and the stable target variable mode component. The average mutual information rate of the two in each synchronization segment is calculated, and continuous synchronization segments with an average mutual information rate higher than the set mutual information threshold are marked as time-domain stable regions. Cross wavelet transform is performed on the stable target variable modal components and the stable covariate modal components. The cross wavelet power spectrum is calculated, and the frequency interval where the local maxima of the cross wavelet power spectrum are located is marked as the frequency domain resonance region. Calculate the area of ​​intersection between the time span of the time domain stable region and the frequency span of the frequency domain resonance region on the time-frequency plane, and take the ratio of the intersection area to the area corresponding to the frequency span of the frequency domain resonance region as the degree of coincidence. Stable target variable modal components with overlap exceeding the overlap threshold are combined with stable covariate modal components into a single prediction unit.

[0010] Furthermore, by embedding a self-attention mechanism into a time series prediction network, and using the stable covariate mode components in the prediction unit as input, the predicted values ​​of the stable target variable mode components are obtained, including: The stable covariate mode components and stable target variable mode components in the prediction unit are aligned by time to construct the input sequence and label sequence; The input sequence is fed into a temporal convolutional network, and multi-scale local temporal features are extracted through dilated causal convolution, outputting a feature sequence. The feature sequence is input into a bidirectional gated recurrent unit, and the context representation sequence is obtained by concatenating the forward hidden state and the backward hidden state. A linear mapping is performed on the context representation sequence to generate a query matrix, a key matrix, and a value matrix. The self-attention weight matrix is ​​calculated, and the self-attention weight matrix and the value matrix are weighted and summed to obtain the self-attention output sequence. The self-attention output sequence is mapped to the predicted values ​​of the modal components of the stable target variable through a fully connected layer.

[0011] Furthermore, by superimposing and reconstructing the predicted values ​​of all stable objective variable modal components, water quality prediction results are obtained, including: The predicted values ​​of the stable target variable modal components in all prediction units are grouped according to the sequence number of the corresponding stable target variable modal components in the original decomposition. The predicted values ​​in each group are linearly superimposed in time order to obtain the preliminary reconstruction sequence of the stable target variable modal components in each group. Calculate the residual sequence between the preliminary reconstructed sequence and the corresponding original target variable sequence. Perform empirical mode decomposition on the residual sequence to obtain residual mode components. Select the components in the residual mode components whose frequencies fall within the frequency band of the stable target variable mode components and perform secondary superposition to obtain the final reconstructed sequence of each group of stable target variable mode components. All the final reconstructed sequences are aligned by time and summed to obtain the reconstructed result sequence. The reconstructed result sequence is then de-standardized to obtain the water quality prediction result.

[0012] On the other hand, the present invention provides a water quality prediction system, comprising: The sequence acquisition module is used to acquire the target variable sequence and covariate sequence of the water area to be predicted; The mode decomposition module is used to perform mode decomposition on the target variable sequence and the covariate sequence respectively, to obtain the target variable mode components and the covariate mode components; The component filtering module is used to perform a significance test by comparing the nonlinear complexity differences between the target variable modal components and covariate modal components and their respective corresponding phase randomized substitute data, eliminating components that fail the test, and retaining stable target variable modal components and stable covariate modal components. The coupling identification module is used to combine the stable target variable modal components and stable covariate modal components whose overlap between the time domain stable domain and the frequency domain resonance region exceeds the overlap threshold into a prediction unit. The time domain stable domain is determined by the instantaneous coupling strength of the sliding window, and the frequency domain resonance region is determined by the cross wavelet power spectrum. The prediction execution module is used to obtain the predicted values ​​of the stable target variable mode components by taking the stable covariate mode components in the prediction unit as input to the time series prediction network with embedded self-attention mechanism. The result reconstruction module is used to overlay and reconstruct the predicted values ​​of the modal components of all stable target variables to obtain water quality prediction results.

[0013] On the other hand, the present invention provides a water quality prediction device, which includes a processor, a memory, and a program or instructions stored in the memory and executable on the processor. When the program or instructions are executed by the processor, a water quality prediction method is implemented.

[0014] On the other hand, the present invention provides a readable storage medium on which a program or instruction is stored, and when the program or instruction is executed by a processor, a water quality prediction method is implemented.

[0015] Compared with the prior art, the present invention has the following beneficial effects: 1. The beneficial effects of this invention are mainly reflected in the improvement of the method for capturing cross-scale coupling relationships in multivariate water quality time series. Since both the target variable sequence and the covariate sequence are simultaneously decomposed into modal components, different frequency components of the covariate can establish direct matching channels with the corresponding frequency band components of the target variable. This eliminates the problem of information structure mismatch between the full-band covariate and the single-band target component in the prior art from the data organization level. By using phase randomization to replace the nonlinear complexity significance test of the data, the modal components are screened for authenticity, eliminating random artifacts that may be mixed in during the decomposition process. This ensures that all components participating in subsequent coupling pairing are stable components with statistically significant nonlinear structures. The dual-domain overlap of the time-domain stable domain and the frequency-domain resonance region is used as the coupling criterion to replace the conventional single-dimensional correlation analysis. This allows the pairing screening to simultaneously consider the temporal persistence and frequency-domain energy synergy of the coupling relationship. The temporal mutual information captures the instantaneous coupling strength of the nonlinearity, and the frequency-domain cross-wavelet power spectrum locates the synchronous oscillation energy. The structure of the intersection of the dual-domain overlap ensures that the selected prediction units have real physical driving significance rather than accidental statistical correlations.

[0016] 2. At the prediction execution level, each prediction unit is independently modeled through a time series prediction network with an embedded self-attention mechanism. The network dynamically weights the context representation of the entire time step, so that historical information that is far from the current time but highly relevant to the current prediction can be directly extracted and utilized. This solves the problem of the decay of key information at distant times step by step in long sequence modeling. After the prediction is completed, the residual decomposition and frequency band matching are used to compensate for the loss of signal components in the same frequency band in the reconstructed sequence due to coupling screening or prediction errors. This makes the water quality prediction results obtained by superposition and reconstruction more completely restore the fluctuation characteristics of each frequency scale in the original sequence. Attached Figure Description

[0017] Figure 1 This is a flowchart of the water quality prediction method of the present invention; Figure 2This is a schematic diagram of the water quality prediction system of the present invention. Detailed Implementation

[0018] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0019] Example 1: Figure 1 The present invention provides a water quality prediction method, comprising: S1: Obtain the target variable sequence and covariate sequence of the water area to be predicted; S2: Perform mode decomposition on the target variable sequence and the covariate sequence respectively to obtain the target variable mode components and the covariate mode components; S3: A significance test is performed by comparing the nonlinear complexity differences between the target variable modal components and covariate modal components and their respective corresponding phase randomized substitute data. Components that fail the test are removed, and stable target variable modal components and stable covariate modal components are retained. S4: Combine the stable target variable modal components and stable covariate modal components whose overlap between the time domain stable domain and the frequency domain resonance region exceeds the overlap threshold into a prediction unit. The time domain stable domain is determined by the instantaneous coupling strength of the sliding window, and the frequency domain resonance region is determined by the cross wavelet power spectrum. S5: By embedding a time series prediction network with a self-attention mechanism, the predicted values ​​of the stable target variable mode components are obtained by taking the stable covariate mode components in the prediction unit as input. S6: Superimpose and reconstruct the predicted values ​​of the modal components of all stable target variables to obtain the water quality prediction results.

[0020] In a specific implementation of S1, the water area to be predicted is a near-shore aquaculture area where water quality monitoring buoys are deployed. Multiple monitoring buoys are distributed within this water area, each with unique location information expressed in latitude and longitude coordinates.

[0021] Export the latitude and longitude coordinate records of each monitoring buoy from the monitoring buoy management system to obtain the location information of the monitoring buoys in the water area to be predicted. Traverse the latitude and longitude coordinates of all monitoring buoys, and mark the monitoring buoy numbers that fall within the boundary coordinate range of the water area to be predicted as target monitoring stations, thus completing the determination of the target monitoring station locations within the water area to be predicted.

[0022] Historical monitoring data recorded at the target monitoring stations are collected. This historical data includes the original target variable sequence and the original covariate sequence. The original target variable sequence is the sequence of dissolved oxygen concentration observations at consecutive, equally spaced sampling times, for example, a time interval of 1 hour, expressed in micromoles per kilogram. The original covariate sequence is at least one of the following: water temperature observation sequence, salinity observation sequence, or chlorophyll concentration observation sequence, all within the same time dimension as the original target variable sequence.

[0023] Outlier removal is performed on the original target variable sequence and the original covariate sequence. A sliding window with a fixed step size is set, the width of which is determined based on the sampling frequency of historical monitoring data, for example, a window width of 24 sampling points. The sliding window is moved sequentially along the time direction, and the kurtosis value of the original target variable sequence within the window is calculated after each movement. The kurtosis value is calculated as the ratio of the fourth central moment of the original target variable sequence within the sliding window to the square of the variance of the original target variable sequence within the sliding window. When the kurtosis value of the sliding window containing a sampling point exceeds a preset kurtosis threshold, the sampling point with the kurtosis value exceeding the threshold is marked as an outlier and removed. The kurtosis threshold is determined based on the upper limit of the baseline range of kurtosis values ​​obtained from historical monitoring data of the target monitoring station under normal operating conditions. The baseline range is taken as the percentile interval of the kurtosis values ​​of the same historical period, with the 95th percentile of the kurtosis values ​​of the same historical period as the upper limit. The original covariate sequence is rearranged column-wise, and the same process of kurtosis value calculation and outlier removal is performed as with the original target variable sequence.

[0024] Missing value imputation is performed on the sequences after outlier removal. The original target variable sequences and original covariate sequences after outlier removal are arranged column-wise to form the imputation matrix. The missing value positions in the imputation matrix are initialized to the mean of the corresponding sequences. In each iteration, principal component analysis is performed on the imputation matrix to extract a preset number of principal components. The preset number of components to be retained is determined based on the proportion of the number of sequences in the imputation matrix, for example, 1 / 3 of the total number of sequences. An approximate matrix is ​​reconstructed using the extracted principal components, and the current value at the missing value position is replaced with the value at the corresponding position in the approximate matrix. This iteration is repeated until the reconstruction error is lower than a preset reconstruction error threshold. The reconstruction error is calculated as the Frobenius norm of the difference between the imputation matrices of two consecutive iterations. The reconstruction error threshold is set based on the magnitude of the non-missing value values ​​in the imputation matrix, for example, 1% of the standard deviation of the non-missing values. After the iteration terminates, the imputed sequence is obtained.

[0025] The filled sequences are standardized to obtain the target variable sequence and covariate sequence. The standardization process involves calculating the minimum and maximum values ​​for each sequence in the filled sequence. The minimum value is subtracted from the value of each sample point, and then divided by the difference between the maximum and minimum values, mapping each sample point to a value range of 0 to 1. The normalized sequences are the target variable sequence and the covariate sequence. The minimum and maximum values ​​of each sequence are saved during the standardization process for use in subsequent destandardization of the water quality prediction results.

[0026] In a specific implementation of S2, when performing modal decomposition on the target variable sequence, the target variable sequence is considered as a curve extending along the time axis. All extreme points are scanned point-by-point along the time axis and identified. Extreme points include maxima larger than their immediate neighbors and minima smaller than their immediate neighbors. The time interval between adjacent extreme points is calculated, and all adjacent time intervals are sorted by value and the median is taken as the local density measure. The boundaries of segments in the target variable sequence where consecutive intervals deviate from the median by a preset deviation factor are defined as segment boundaries. The preset deviation factor is set according to the non-stationarity of the target variable sequence, for example, three times the median. At the segment boundaries, the target variable sequence is divided into multiple target variable sub-segments. The difference between the time interval of extreme points within each target variable sub-segment and the median of the time intervals of extreme points within that sub-segment is less than the difference between the boundary and the median of the time intervals of extreme points in adjacent sub-segments.

[0027] For each target variable segment, the power spectral density (PSD) of that segment is calculated. The calculation involves performing a Fast Fourier Transform (FFT) on the target variable segment to obtain a spectral sequence. The amplitude value at each frequency point in the spectral sequence is squared and divided by the length of the target variable segment to obtain its PSD curve. All peak points are detected on the PSD curve. A frequency point whose power value is simultaneously greater than the power values ​​of its two adjacent frequencies is marked as a peak point. All peak points are sorted from largest to smallest power value. Peak points whose power values ​​exceed a power threshold are selected as the dominant peaks. The power threshold is a fixed percentage of the maximum power value on the PSD curve, for example, 10%. This percentage is empirically determined based on the ratio of the base power value to the peak power value in a sequence segment without significant periodic components. The number of dominant peaks is counted and directly used as the decomposition level for the adaptive noise completeness set empirical mode decomposition of that target variable segment. When the number of dominant peaks in a sub-segment of the target variable is large, the frequency components within that sub-segment are rich, requiring more decomposition levels to separate fluctuations at different scales. When the number of dominant peaks in a sub-segment of the target variable is small, the frequency components within that sub-segment are homogeneous, requiring fewer decomposition levels to achieve separation and avoiding the generation of spurious components with no physical meaning.

[0028] Adaptive noise-complete ensemble empirical mode decomposition (EMD) is performed on each target variable sub-segment. The process involves repeatedly adding paired positive and negative Gaussian white noise sequences to the target variable sub-segment. The amplitude of the added noise sequences is proportionally set based on the standard deviation of the target variable sub-segment; for example, the noise sequence amplitude is 20% of the standard deviation. EMD is then performed on the synthesized sequence after each addition of noise sequences. The first intrinsic mode function (IMF) obtained from each EMD is extracted, and the overall average of all obtained IMFs is calculated to obtain the first sub-mode component corresponding to the target variable sub-segment. The first sub-mode component is subtracted from the target variable sub-segment to obtain the first residual sequence. Paired positive and negative Gaussian white noise sequences are then added to the first residual sequence, and EMD is performed again. The overall average is then calculated to obtain the second sub-mode component. Repeat the above process of adding noise sequence, empirical mode decomposition, and calculating the overall average. After each iteration of adding noise sequence, empirical mode decomposition, and calculating the overall average, one sub-mode component is obtained and a new residual sequence is generated. Continue operating on the new residual sequence until the number of sub-mode components obtained reaches the decomposition level set for the target variable segment, and then stop to obtain all sub-mode components corresponding to the target variable segment.

[0029] The submodal components of all target variable segments are concatenated in chronological order. Following the original temporal order of the sampling points within each target variable segment, submodal components belonging to the same decomposition level from different target variable segments are concatenated end-to-end. Each concatenated sequence represents one target variable modal component. All target variable modal components are arranged in order from high frequency to low frequency.

[0030] The covariate sequence is processed to obtain covariate modal components. All covariate extreme points are identified by traversing the covariate sequence. Based on the local density of the covariate extreme point distribution, segment boundaries are determined using the same segmentation boundary determination method as the target variable sequence, dividing the covariate sequence into multiple covariate sub-segments. The power spectral density of each covariate sub-segment is calculated, and the decomposition level of the covariate sub-segment is adaptively set according to the number of main peaks in the power spectral density. Adaptive noise complete ensemble empirical mode decomposition is performed on each covariate sub-segment to obtain the corresponding sub-modal components. The sub-modal components of all covariate sub-segments are concatenated in chronological order to obtain the covariate modal components.

[0031] In a specific implementation of S3, a significance test is performed on each component of the target variable modal component and covariate modal component obtained in S2. Taking a single target variable modal component as an example, the process of generating phase randomized substitute data involves performing a Fourier transform on the target variable modal component to obtain an amplitude spectrum and a phase spectrum. Keeping the amplitude spectrum values ​​unchanged, the phase value at each frequency point in the phase spectrum is replaced with a random number uniformly generated in the interval between 0 and 2π, forming a new phase spectrum. The amplitude spectrum and the new phase spectrum are combined and then subjected to an inverse Fourier transform to obtain one phase randomized substitute data. The above phase replacement and inverse Fourier transform process is repeated multiple times, for example, 100 times, to obtain 100 phase randomized substitute data, constituting a set of phase randomized substitute data corresponding to the target variable modal component. The phase randomization process destroys the phase self-organization relationship generated by nonlinear dynamics in the original sequence, while retaining all the linear autocorrelation structure contained in the amplitude spectrum. Therefore, the phase randomized substitute data represents a null hypothesis signal with only linear structure and no nonlinear structure.

[0032] Calculate the permutation entropy of the modal components of the target variable. Permutation entropy measures the nonlinear complexity of a sequence. The calculation process involves reconstructing the phase space of the sequence according to a preset embedding dimension and time delay, resulting in a set of state vectors. Elements within each state vector are sorted by value, and the corresponding index permutation pattern is recorded. The frequency of each permutation pattern in all state vectors is counted, and the probability of each permutation pattern is calculated. The permutation entropy is calculated as the probability-weighted information entropy of all permutation patterns, expressed as: H = -Σp(π) × log(p(π)), where H represents the permutation entropy, π represents the permutation pattern, p(π) represents the probability of permutation pattern π, and Σ represents the summation over all possible permutations. The embedding dimension is set according to the sequence length; for example, when the sequence length exceeds 1000 sampling points, the embedding dimension is set to 6, and the time delay is set to 1 sampling step. The embedding dimension determines the total number of permutation patterns. The total number of permutation patterns is the factorial of the embedding dimension. If the embedding dimension is too large, the number of permutation patterns will exceed the statistical quantity that the sequence length can support. If the embedding dimension is too small, it will be unable to characterize the fine nonlinear structure of the sequence. Therefore, the embedding dimension must be matched with the sequence length.

[0033] Calculate the mean and standard deviation of the permutation entropy of a set of phase randomized substitute data corresponding to the modal component of the target variable. Calculate the permutation entropy of each phase randomized substitute data point in the set sequentially. The mean permutation entropy is calculated as the sum of the permutation entropies of all phase randomized substitute data points divided by the number of phase randomized substitute data points. The standard deviation of the permutation entropy is calculated as the square root of the sum of the squares of the differences between the permutation entropy of each phase randomized substitute data point and the mean permutation entropy, divided by the number of phase randomized substitute data points.

[0034] The permutation entropy of the target variable's modal component is compared with the mean and standard deviation of the permutation entropy of the corresponding set of phase-randomized surrogate data. The lower bound of the comparison interval is the mean permutation entropy minus a preset multiple multiplied by the standard deviation of the permutation entropy, and the upper bound is the mean permutation entropy plus the preset multiple multiplied by the standard deviation of the permutation entropy. The preset multiple is set according to the confidence level required in the significance test. The correspondence between the preset multiple and the confidence level is the quantile value under the standard normal distribution. For example, the preset multiple is 1.96 when a confidence level of 95% is required, and 2.58 when a confidence level of 99% is required. If the permutation entropy of the target variable's modal component is not less than the lower bound of the comparison interval and not greater than the upper bound of the comparison interval, then the target variable's modal component is determined to have failed the significance test and is removed. Failing the significance test means that the nonlinear complexity of the target variable's modal component is not significantly different from the null hypothesis signal, belonging to random perturbation artifacts rather than the intrinsic deterministic mode of the data. If the permutation entropy of the modal component of the target variable is less than the lower bound of the comparison interval or greater than the upper bound of the comparison interval, then the modal component of the target variable is determined to have passed the significance test and is retained. The modal components of the target variable that pass the significance test are regarded as stable modal components of the target variable. These components contain a nonlinear phase structure that is statistically significantly different from that of a linear random signal.

[0035] The same significance test and elimination / retention methods were applied to the covariate modal components as to the target variable modal components. For each covariate modal component, a corresponding set of phase randomized substitute data was generated. The permutation entropy of the covariate modal component and the mean and standard deviation of the permutation entropy of the corresponding set of phase randomized substitute data were calculated, and the same interval comparisons were performed. Covariate modal components that passed the significance test were considered stable covariate modal components, while those that failed the significance test were eliminated.

[0036] In a specific implementation of S4, for each stable target variable modal component and each stable covariate modal component retained in S3, a two-dimensional coupling determination in the time domain and frequency domain is performed to determine whether the two are combined into a prediction unit.

[0037] In the time domain, the instantaneous coupling strength is calculated and the time-domain stability region is determined using a sliding window. The stable covariate modal components and the stable target variable modal components are aligned along the time axis at the same sampling time. A fixed-length sliding window is set, the length of which is determined by the principal period of the stable target variable modal component; for example, the number of sampling points corresponding to twice the period of the maximum peak value of the power spectral density curve of the stable target variable modal component. Starting from the beginning of both sequences, the sliding window moves along the time axis point by point, moving one sampling step at a time. A segment of the stable covariate modal component and a segment of the stable target variable modal component within the coverage of the sliding window are extracted, forming a synchronization segment. For each synchronization segment, the average mutual information rate between the stable covariate modal component segment and the stable target variable modal component segment is calculated. The average mutual information rate is calculated by dividing the numerical range of each segment into several equal intervals, the number of intervals determined by the segment length; for example, the square root of the segment length is rounded up.

[0038] The joint frequency of the values ​​of two sub-segments falling into each interval combination and the marginal frequency within each interval are statistically analyzed. The mutual information is calculated as: I = Σp(xi,yj) × log[p(xi,yj) / (p(xi) × p(yj))]; where I represents the mutual information, xi represents the i-th interval into which the values ​​of the stable covariate modal component sub-segments fall, yj represents the j-th interval into which the values ​​of the stable target variable modal component sub-segments fall, p(xi,yj) represents the joint probability estimate of the values ​​of the stable covariate modal component sub-segments falling into the i-th interval and the values ​​of the stable target variable modal component sub-segments falling into the j-th interval, p(xi) represents the marginal probability estimate of the values ​​of the stable covariate modal component sub-segments falling into the i-th interval, p(yj) represents the marginal probability estimate of the values ​​of the stable target variable modal component sub-segments falling into the j-th interval, log represents the logarithm to the base 2, and Σ represents the summation over all interval combinations. The average mutual information rate is calculated as: R = I / L; where R represents the average mutual information rate, I represents the mutual information content, and L represents the length of the synchronization segment. The average mutual information rate value is obtained sequentially after each movement of the sliding window. All average mutual information rates are arranged in chronological order, and the time positions of synchronization segments with average mutual information rates higher than the set mutual information threshold are marked as effective coupling points.

[0039] The mutual information threshold is determined based on the overall distribution of the average mutual information rate between the stable target variable modal components and the stable covariate modal components. Specifically, the median and interquartile range of the average mutual information rate of all synchronization segments are calculated, and the median plus 1.5 times the interquartile range is taken as the mutual information threshold.

[0040] Consecutive adjacent effective coupling points on the time axis are merged into a single continuous coupling period, and the time span of each coupling period constitutes a time-domain stability domain. If there is an interval point between two effective coupling points that does not reach the mutual information threshold, then the two sides of the interval point belong to different time-domain stability domains. The average mutual information rate, rather than the nonlinear correlation coefficient, is used to measure the instantaneous coupling strength because mutual information can capture the nonlinear coupling relationship that may exist between two modal components, while the linear correlation coefficient is only sensitive to linear relationships and will miss the nonlinear cooperative variation patterns between driving and response variables that are common in marine environments.

[0041] In the frequency domain, the cross-wavelet power spectrum is calculated and the frequency domain resonance region is determined through cross-wavelet transform. Continuous cross-wavelet transforms are performed on the stable target variable mode components and the stable covariate mode components. One mother wavelet function is selected, such as the complex Morlet wavelet, and continuous wavelet transforms are performed on the stable target variable mode components and the stable covariate mode components respectively, yielding the wavelet coefficient time-frequency matrices of the stable target variable mode components and the stable covariate mode components. The cross-wavelet power spectrum is calculated as: P(f,t)=|Wx(f,t)| 2 ×|Wy(f,t)| 2 Where P(f,t) represents the power value of the cross wavelet power spectrum at frequency f and time t, Wx(f,t) represents the wavelet coefficients at frequency f and time t in the time-frequency matrix of the wavelet coefficients of the stable objective variable modal component, and Wy(f,t) represents the wavelet coefficients at frequency f and time t in the time-frequency matrix of the wavelet coefficients of the stable covariate modal component. 2 The modulus is squared. The cross-wavelet power spectrum is a two-dimensional matrix, with rows corresponding to frequency scales and columns corresponding to time positions. The value of each element in the matrix represents the covariant power intensity of the two components at the corresponding frequency and time position. Local maxima are detected in the cross-wavelet power spectrum. When the value of an element in the cross-wavelet power spectrum is simultaneously greater than its two adjacent elements in both the frequency dimension and the time dimension, that element is marked as a local maximum. Expanding in both upward and downward directions along the frequency dimension around each local maximum, the expansion stops when the power value of the element drops below half of the power value of the local maximum. The frequency range covered by the upward and downward expansion is the frequency range corresponding to that local maximum. The frequency ranges with the highest power values ​​of the local maximums among all the frequency ranges corresponding to all local maximums are merged to form the frequency domain resonance region. The frequency domain resonance region indicates that at these frequency points, the stable target variable mode component and the stable covariate mode component have significant frequency domain resonance energy, that is, they oscillate synchronously at that frequency and their energies reinforce each other, rather than being accidental frequency overlaps.

[0042] The intersection of the time-domain stable domain and the frequency-domain resonance region is calculated on the time-frequency plane. The time-frequency plane has time as the horizontal axis and frequency as the vertical axis. The time span occupied by the time-domain stable domain on the time-frequency plane is from its start to its end, and the frequency span is the entire frequency range of the stable target variable's modal components. The time span occupied by the frequency-domain resonance region on the time-frequency plane is the entire time range of the stable target variable's modal components, and the frequency span is from the lower frequency limit to the upper frequency limit of the frequency-domain resonance region. The overlap degree is calculated as: C = Soverlap / Stotal; where C represents the overlap degree, Soverlap represents the area of ​​the intersection region of the time span of the time-domain stable domain and the frequency span of the frequency-domain resonance region on the time-frequency plane, and Stotal represents the area corresponding to the frequency span of the frequency-domain resonance region. The area corresponding to the frequency span of the frequency-domain resonance region is calculated as the product of the frequency span of the frequency-domain resonance region and the entire time span of the stable target variable's modal components. The closer the overlap is to 1, the more it indicates that the stable time domain contains most of the resonant frequency components of the frequency domain resonance region, and the coupling relationship holds simultaneously in both the time and frequency domains. The overlap is compared to an overlap threshold. Stable target variable modal components with overlap exceeding the threshold are combined with stable covariate modal components into a single prediction unit. A prediction unit consisting of one stable target variable modal component and one stable covariate modal component indicates that the change in the stable covariate modal component has a stable drive-response coupling relationship with the change in the stable target variable modal component in both the time and frequency domains.

[0043] The overlap threshold is set according to the stringency of the prediction unit combination. For example, an overlap threshold of 0.7 means that more than 70% of the resonant frequency components in the frequency domain resonance region must be within a stable coupling period in the time domain. A higher overlap threshold results in fewer prediction units but more reliable coupling relationships, while a lower overlap threshold results in more prediction units but may include some false positive couplings. The overlap threshold is selected based on the emphasis on accuracy and recall in the actual prediction task.

[0044] In a specific implementation of S5, each prediction unit obtained from the combination in S4 is predicted one by one using a time-series prediction network with an embedded self-attention mechanism. Taking a single prediction unit as an example, the stable covariate mode components and stable target variable mode components in the prediction unit are mapped one-to-one according to the sampling times on the time axis. All sampled values ​​of the stable covariate mode components from the start time to the end time on the time axis are used as the input sequence, and all sampled values ​​of the stable target variable mode components from the start time to the end time on the same time axis are used as the label sequence. The input sequence and the label sequence correspond one-to-one in the time dimension, forming a training sample pair for supervised learning.

[0045] The input sequence is fed into a temporal convolutional network (TCNN), whose input layer receives the raw numerical values ​​of the input sequence. The TCNN consists of multiple stacked residual blocks, each containing two dilated causal convolutional layers. During convolution operations, the kernels of these layers operate only on the input elements from the current time step and previous time steps. The dilation factor increases exponentially with the stacking depth of the residual blocks; for example, the dilation factor of the first residual block is 1, the second is 2, and the third is 4. Through layer-by-layer processing of these dilated causal convolutions, the receptive field expands exponentially with each layer. Each dilated causal convolution is followed by weight normalization and an activation function using Gaussian error linear units. The output of each residual block is the element-wise sum of the output of the convolutional path within the residual block and the input of the residual block aligned with the 1x1 convolution dimension. After processing through all residual blocks, the input sequence outputs a single feature sequence. The length of the feature sequence is equal to the length of the input sequence. The feature vector at each time step in the feature sequence integrates multi-scale local temporal information within that time step and its historical dependencies. Temporal convolutional networks enable effective backpropagation of gradients in deep networks through residual connections, avoiding the gradient vanishing problem and making them suitable for handling long-term dependencies of water quality modal components.

[0046] The feature sequence is input into a bidirectional gated recurrent unit (BRN). The BRN consists of one forward-gated recurrent unit and one backward-gated recurrent unit connected in parallel. The forward-gated recurrent unit processes the feature sequence step-by-step along the forward temporal direction, outputting a forward hidden state at each time step. The backward-gated recurrent unit processes the feature sequence step-by-step along the reverse temporal direction, outputting a backward hidden state at each time step. The forward and backward hidden states at each time step are concatenated as vectors, with the vectors joined end-to-end to form a single context representation vector. The context representation vectors from all time steps are arranged chronologically to form a context representation sequence. The BRN gates within the bidirectional gated recurrent unit control the retention and forgetting of historical information through reset and update gates. The reset gate controls the proportion of the previous hidden state written into the candidate hidden state, and the update gate controls the weighted combination of the previous hidden state and the candidate hidden state when generating the current hidden state. The output of the reset gate is calculated by adding the previous hidden state and the current input after linear transformation of their respective weight matrices, and then mapping them to the 0-1 interval using the Sigmoid function. The output of the update gate is also mapped using the Sigmoid function. The candidate hidden state is jointly determined by the current input and the previous hidden state adjusted by the reset gate. This bidirectional processing ensures that the vector at each time step in the context representation sequence contains both the temporal evolution information before that time step and the temporal context information after that time step, overcoming the limitation of unidirectional recurrent networks that can only use past information for prediction.

[0047] Self-attention is computed on the context representation sequence. The context representation vector at each time step is simultaneously input into three independent linear mapping layers, each with its own independent weight matrix and bias vector. The first linear mapping layer maps the context representation vector to a query vector, the second to a key vector, and the third to a value vector. The query vectors from all time steps are arranged column-wise to form a query matrix, the key vectors from all time steps are arranged column-wise to form a key matrix, and the value vectors from all time steps are arranged column-wise to form a value matrix.

[0048] The self-attention weight matrix is ​​calculated as the product of the transpose of the query matrix and the key matrix, divided by a scaling factor. The scaling factor is the square root of the key vector dimension. Each column of the product is then normalized using the Softmax function, which maps each element in a column to a probability distribution between 0 and 1, where the sum of all elements in the column is 1. The scaling factor is introduced to prevent the gradient of the Softmax function from approaching zero when the key vector dimension is large, as this would result in an excessively large inner product of the query and key matrices. The output vector at each time step in the self-attention output sequence is calculated as the sum of the products of the corresponding columns of the value matrix and the self-attention weight matrix. The self-attention mechanism allows the network to automatically allocate attention to the value vectors at each time step based on the similarity between the query vector and the key vectors at each time step, thereby highlighting the contextual information of the time steps crucial to the current prediction task and suppressing interference from irrelevant time steps. Compared to using only the hidden state of the final time step of the bidirectional gated recurrent unit for prediction, the self-attention mechanism preserves and dynamically weights the context representation of all time steps, solving the problem that key information at the far end is gradually diluted during the time-step transmission process in long sequence modeling.

[0049] The output vector at each time step of the self-attention output sequence is input into a fully connected layer. The fully connected layer compresses the dimension of the output vector layer by layer, ultimately mapping it to a single scalar value. This scalar value is the predicted value of the stable target variable mode component corresponding to that time step. Gaussian error linear units are used as activation functions between fully connected layers. No activation function is used in the final output layer to ensure that the numerical range of the predicted values ​​is consistent with the numerical range of the stable target variable mode components. The predicted values ​​of the stable target variable mode components at all time steps are arranged in chronological order to obtain the complete prediction sequence of the stable target variable mode components corresponding to that prediction unit.

[0050] In a specific implementation of S6, the predicted values ​​of the stable target variable modal components corresponding to each prediction unit output from S5 are superimposed and reconstructed to obtain the water quality prediction result. In S2, the target variable modal components are arranged from high to low frequency and assigned a unique original decomposition index. The stable target variable modal components retain the original decomposition index of the corresponding target variable modal components in S2. The predicted values ​​of the stable target variable modal components in all prediction units are grouped according to their original decomposition index, and the predicted values ​​of stable target variable modal components with the same original decomposition index are grouped together. For stable target variable modal components that did not form a prediction unit with any stable covariate modal components in S4, their predicted values ​​are recorded as a zero-value sequence. The time length of the zero-value sequence is the same as the time length of the original target variable sequence, and the value at each time step is 0.

[0051] For each group, the predicted values ​​of all stable target variable modal components within the group are summed point-by-point at the same time step. The summation yields a preliminary reconstruction sequence corresponding to the original decomposition sequence number of that group. If a group contains only one predicted value for a stable target variable modal component, the preliminary reconstruction sequence is that predicted value itself. The length of the preliminary reconstruction sequence is equal to the length of the original target variable sequence.

[0052] For each group of stable target variable modal components, the residual sequence between the preliminary reconstructed sequence and the original target variable sequence is calculated. The original target variable sequence is the target variable sequence obtained in S1, i.e., the standardized target variable sequence. The residual sequence is calculated as: Rk(n) = Y(n) − Pk(n); where Rk(n) represents the residual value of the residual sequence corresponding to the k-th stable target variable modal component at the n-th time step, Y(n) represents the value of the original target variable sequence at the n-th time step, Pk(n) represents the value of the preliminary reconstructed sequence corresponding to the k-th stable target variable modal component at the n-th time step, k is the original decomposition index, and n is the time step index. The residual sequence reflects the deviation between the preliminary reconstructed sequence and the original target variable sequence at each time step. This deviation includes the contributions of other modal components not captured by the current prediction unit during the preliminary reconstruction process, as well as noise.

[0053] Empirical Mode Decomposition (EMD) is performed on the residual sequence to obtain multiple residual mode components and one residual trend term. The EMD process involves finding all maxima and minima of the residual sequence, fitting the upper and lower envelopes using cubic spline interpolation, calculating the mean curves of the upper and lower envelopes, subtracting the mean curves from the residual sequence to obtain a candidate component, and repeating this process until a candidate component satisfies the intrinsic mode function condition. This candidate component is then used as a residual mode component, and the residual sequence is subtracted from it to obtain the remaining sequence. This process is repeated on the remaining sequence until the number of extreme points is less than two, finally yielding all residual mode components. For each residual mode component, a Hilbert transform is used to calculate the instantaneous frequency sequence, and the time average of the instantaneous frequency sequence is taken as the average frequency of the residual mode component.

[0054] The frequency band range of the k-th stable target variable mode component is determined. The lower and upper limits of the frequency band range are determined based on the power spectral density of the corresponding target variable mode component in S2. The frequencies corresponding to the decrease in power spectral density from the main peak to half of the peak power are used as the lower and upper frequency limits. Residual mode components whose average frequency falls between these lower and upper frequency limits are selected as the residual components to be superimposed. The selected residual components to be superimposed are added step-by-step with the preliminary reconstruction sequence of the k-th stable target variable mode component to obtain the final reconstruction sequence of the k-th stable target variable mode component. The secondary superposition reintroduces the residual fluctuation signal in the same frequency band as the mode component during the prediction process, compensating for the loss of information within the frequency band caused by coupling selection or prediction errors.

[0055] The final reconstructed sequences of all stable target variable modal components are time-aligned and summed step-by-step to obtain the reconstructed result sequence. This reconstructed result sequence contains the contributions of all stable target variable modal components, and its numerical range remains within the normalized scale. The reconstructed result sequence is then de-standardized by calling the minimum and maximum values ​​of the target variable sequence saved during the normalization process in S1. The de-standardization calculation is: X(n) = Z(n) × (max − min) + min; where X(n) represents the water quality prediction result at the nth time step, Z(n) represents the value of the reconstructed result sequence at the nth time step, max represents the maximum value of the target variable sequence, and min represents the minimum value of the target variable sequence. The de-standardization process restores the prediction results from the normalized interval to the original physical numerical range and dimensions of the target variables, yielding the final water quality prediction result, with dissolved oxygen concentration in micromoles per kilogram.

[0056] It is worth noting that in the complete scheme consisting of S1 to S6, each step addresses a different aspect of the core problem of the difficulty in capturing the cross-scale time-delay coupling relationship between covariates and target variables in water quality prediction, forming a progressive solution chain. S2 simultaneously decomposes the target variable sequence and the covariate sequence, eliminating the structural mismatch between the full-band information of the covariates and the single-band components of the target variables in existing technologies from the source. S3 eliminates spurious components with no physical meaning that may be generated during the decomposition process by using phase randomization to replace the nonlinear complexity significance test of the data, ensuring the reliability of the components used in the subsequent prediction unit combination. S4 introduces a dual-domain overlap judgment of the time-domain stable domain and the frequency-domain resonance region, cross-validating the coupling relationship between modal components from two independent dimensions: time-domain persistence and frequency-domain resonance. The time-domain mutual information captures the nonlinear coupling strength, and the frequency-domain cross-wavelet power spectrum identifies the synchronous oscillation energy. The operation of taking the intersection of the dual-domain overlap ensures that the selected prediction units simultaneously meet the two conditions of time-domain stability and frequency-domain resonance, eliminating false positive pairings that may be retained by single-dimensional judgment. S5 utilizes a time-series prediction network with an embedded self-attention mechanism to independently model each prediction unit, retaining contextual information across all time steps for dynamic weighting. S6 compensates for the loss of in-band signal components during prediction through residual decomposition and secondary superposition of frequency band matching.

[0057] Example 2: Figure 2 A schematic diagram of the water quality prediction system of the present invention is provided. The water quality prediction system includes: The sequence acquisition module is used to acquire the target variable sequence and covariate sequence of the water area to be predicted; The mode decomposition module is used to perform mode decomposition on the target variable sequence and the covariate sequence respectively, to obtain the target variable mode components and the covariate mode components; The component filtering module is used to perform a significance test by comparing the nonlinear complexity differences between the target variable modal components and covariate modal components and their respective corresponding phase randomized substitute data, eliminating components that fail the test, and retaining stable target variable modal components and stable covariate modal components. The coupling identification module is used to combine the stable target variable modal components and stable covariate modal components whose overlap between the time domain stable domain and the frequency domain resonance region exceeds the overlap threshold into a prediction unit. The time domain stable domain is determined by the instantaneous coupling strength of the sliding window, and the frequency domain resonance region is determined by the cross wavelet power spectrum. The prediction execution module is used to obtain the predicted values ​​of the stable target variable mode components by taking the stable covariate mode components in the prediction unit as input to the time series prediction network with embedded self-attention mechanism. The result reconstruction module is used to overlay and reconstruct the predicted values ​​of the modal components of all stable target variables to obtain water quality prediction results.

[0058] Example 3: The present invention provides a water quality prediction device, which includes a processor, a memory, and a program or instructions stored in the memory and executable on the processor. When the program or instructions are executed by the processor, a water quality prediction method is implemented.

[0059] Example 4: The present invention provides a readable storage medium on which a program or instruction is stored, and when the program or instruction is executed by a processor, a water quality prediction method is implemented.

[0060] All calculations involved in the embodiments are performed using dimensionless numerical values, and the preset parameters and thresholds in the calculations can be set by those skilled in the art according to actual conditions.

[0061] This technical solution can be flexibly deployed, for example, as embedded software running on device hardware, or installed on personal computers or other smart terminals with user interfaces, thus adapting to various hardware environments and usage requirements.

[0062] The above solutions can be implemented in software, hardware, firmware, or a combination thereof. When implemented in software, they are presented, in whole or in part, as a computer program product, containing one or more computer instructions or programs. When these instructions or programs are loaded and executed on a computer, results are produced corresponding to the processes or functions of the embodiments of this application. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. Computer instructions can be stored in a computer-readable storage medium or transmitted from one medium to another via wired or wireless means, such as from a website, server, or data center via wired means such as fiber optic cables, twisted-pair cables, or coaxial cables, or wireless means such as infrared or microwaves to another site. A computer-readable storage medium refers to any usable medium that a computer can access or a data storage device such as a server or data center that contains one or more usable media, including magnetic media such as floppy disks, hard disks, and magnetic tapes, optical media such as DVDs, and semiconductor media such as solid-state drives.

[0063] The specific working process of the system, device and module can be found in the method embodiment, and will not be repeated here.

[0064] The disclosed systems, devices, and methods can be implemented in other ways. The device embodiments are for illustrative purposes only, and the module division is only a logical division. In practice, different divisions can be implemented, such as merging or integrating multiple modules or components, or omitting some features. Coupling, direct coupling, or communication connections between the components can be achieved through interfaces, while indirect coupling or communication connections can take electrical, mechanical, or other forms.

[0065] The modules described as separate components may or may not be physically separated. The components shown as modules can be physical hardware or software, and can be deployed centrally or distributed across multiple network nodes. Some or all of the modules can be selected to achieve the purpose of this embodiment, depending on actual needs.

[0066] The functional modules in each embodiment can be integrated into one processing module, or they can exist independently, or two or more modules can be integrated into one.

[0067] If the functionality is implemented as a software module and used as an independent product, it can be stored in a computer-readable storage medium. Under this understanding, the substantial contribution of the technical solution of this application can be embodied in the form of a software product. This computer software product is stored in a storage medium and contains instructions to cause a computer device, such as a personal computer, server, or network device, to execute all or part of the steps of the methods in the embodiments of this application. The storage medium includes any medium capable of storing program code, such as a USB flash drive, portable hard drive, read-only memory, random access memory, magnetic disk, or optical disk.

[0068] The above are merely specific embodiments of this application, and the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should fall within the scope of protection of this application.

Claims

1. A water quality prediction method, characterized in that, include: S1: Obtain the target variable sequence and covariate sequence of the water area to be predicted; S2: Perform mode decomposition on the target variable sequence and the covariate sequence respectively to obtain the target variable mode components and the covariate mode components; S3: A significance test is performed by comparing the nonlinear complexity differences between the target variable modal components and covariate modal components and their respective corresponding phase randomized substitute data. Components that fail the test are removed, and stable target variable modal components and stable covariate modal components are retained. S4: Combine the stable target variable modal components and stable covariate modal components whose overlap between the time domain stable domain and the frequency domain resonance region exceeds the overlap threshold into a prediction unit. The time domain stable domain is determined by the instantaneous coupling strength of the sliding window, and the frequency domain resonance region is determined by the cross wavelet power spectrum. S5: By embedding a time series prediction network with a self-attention mechanism, the predicted values ​​of the stable target variable mode components are obtained by taking the stable covariate mode components in the prediction unit as input. S6: Superimpose and reconstruct the predicted values ​​of the modal components of all stable target variables to obtain the water quality prediction results.

2. The water quality prediction method according to claim 1, characterized in that, Obtain the target variable sequence and covariate sequence for the water area to be predicted, including: Obtain the location information of monitoring buoys in the water area to be predicted, and determine the target monitoring station in the water area to be predicted based on the location information; Collect historical monitoring data recorded at the target monitoring stations. The historical monitoring data includes the original target variable sequence and the original covariate sequence. Calculate the kurtosis values ​​of the original target variable sequence and the original covariate sequence within the sliding window, and remove outliers whose kurtosis values ​​exceed the kurtosis threshold; Iterative principal component analysis was used to fill missing values ​​in the sequences after outlier removal to obtain the filled sequences. The filled sequence is standardized to obtain the target variable sequence and the covariate sequence.

3. The water quality prediction method according to claim 1, characterized in that, Modal decomposition is performed on the target variable sequence and the covariate sequence respectively to obtain the target variable modal components and the covariate modal components, including: The target variable sequence is divided into multiple target variable sub-segments by identifying extreme points and determining the segmentation boundary based on the local density of the extreme point distribution. For each target variable segment, calculate the power spectral density of the target variable segment, and adaptively set the number of decomposition layers of the target variable segment according to the number of main peaks of the power spectral density. Adaptive noise complete set empirical mode decomposition is performed on each target variable sub-segment to obtain the sub-mode components corresponding to each target variable sub-segment; The submodal components of all target variable segments are concatenated in chronological order to obtain the target variable modal components; The covariate sequence is processed using the same method as the target variable sequence to obtain the covariate modal components.

4. The water quality prediction method according to claim 1, characterized in that, A significance test was performed by comparing the nonlinear complexity differences between the modal components of the target variable and the modal components of the covariates and their corresponding phase randomized surrogate data. Components that failed the test were removed, and the stable modal components of the target variable and the stable modal components of the covariates were retained, including: Fourier transforms are performed on the target variable modal components and covariate modal components respectively, keeping the amplitude spectrum unchanged. The phase spectrum is randomly shuffled and then inverse Fourier transform is performed. This process is repeated multiple times to generate a corresponding set of phase randomization replacement data. Calculate the permutation entropy of the modal components of the target variable and the modal components of the covariates, as well as the mean and standard deviation of the permutation entropy of the corresponding phase randomized substitute data; For the modal components of the target variable, if their permutation entropy falls within the range of the permutation entropy mean plus or minus a preset multiple of the permutation entropy standard deviation of the corresponding phase randomized substitute data, then it is determined that the significance test has not been passed and the variable is removed; otherwise, it is determined that the significance test has been passed and the variable is retained. The same significance test and elimination / retention method were used for the covariate modal components as for the target variable modal components; The target variable modal components that pass the significance test are taken as stable target variable modal components, and the covariate modal components that pass the significance test are taken as stable covariate modal components.

5. The water quality prediction method according to claim 1, characterized in that, The stable target variable modal components and stable covariate modal components whose overlap between the time-domain stable domain and the frequency-domain resonance region exceeds the overlap threshold are combined into a prediction unit, including: For each stable target variable mode component and each stable covariate mode component, a sliding window is used to extract the synchronization segments of the stable covariate mode component and the stable target variable mode component. The average mutual information rate of the two in each synchronization segment is calculated, and continuous synchronization segments with an average mutual information rate higher than the set mutual information threshold are marked as time-domain stable regions. Cross wavelet transform is performed on the stable target variable modal components and the stable covariate modal components. The cross wavelet power spectrum is calculated, and the frequency interval where the local maxima of the cross wavelet power spectrum are located is marked as the frequency domain resonance region. Calculate the area of ​​intersection between the time span of the time domain stable region and the frequency span of the frequency domain resonance region on the time-frequency plane, and take the ratio of the intersection area to the area corresponding to the frequency span of the frequency domain resonance region as the degree of coincidence. Stable target variable modal components with overlap exceeding the overlap threshold are combined with stable covariate modal components into a single prediction unit.

6. The water quality prediction method according to claim 1, characterized in that, A time series prediction network with embedded self-attention mechanism is used to obtain the predicted values ​​of the stable target variable mode components by taking the stable covariate mode components in the prediction unit as input, including: The stable covariate mode components and stable target variable mode components in the prediction unit are aligned by time to construct the input sequence and label sequence; The input sequence is fed into a temporal convolutional network, and multi-scale local temporal features are extracted through dilated causal convolution, outputting a feature sequence. The feature sequence is input into a bidirectional gated recurrent unit, and the context representation sequence is obtained by concatenating the forward hidden state and the backward hidden state. A linear mapping is performed on the context representation sequence to generate a query matrix, a key matrix, and a value matrix. The self-attention weight matrix is ​​calculated, and the self-attention weight matrix and the value matrix are weighted and summed to obtain the self-attention output sequence. The self-attention output sequence is mapped to the predicted values ​​of the modal components of the stable target variable through a fully connected layer.

7. The water quality prediction method according to claim 1, characterized in that, By superimposing and reconstructing the predicted values ​​of all stable objective variable modal components, water quality prediction results are obtained, including: The predicted values ​​of the stable target variable modal components in all prediction units are grouped according to the sequence number of the corresponding stable target variable modal components in the original decomposition. The predicted values ​​in each group are linearly superimposed in time order to obtain the preliminary reconstruction sequence of the stable target variable modal components in each group. Calculate the residual sequence between the preliminary reconstructed sequence and the corresponding original target variable sequence. Perform empirical mode decomposition on the residual sequence to obtain residual mode components. Select the components in the residual mode components whose frequencies fall within the frequency band of the stable target variable mode components and perform secondary superposition to obtain the final reconstructed sequence of each group of stable target variable mode components. All the final reconstructed sequences are aligned by time and summed to obtain the reconstructed result sequence. The reconstructed result sequence is then de-standardized to obtain the water quality prediction result.

8. A water quality prediction system, used to implement the water quality prediction method according to any one of claims 1-7, characterized in that, include: The sequence acquisition module is used to acquire the target variable sequence and covariate sequence of the water area to be predicted; The mode decomposition module is used to perform mode decomposition on the target variable sequence and the covariate sequence respectively, to obtain the target variable mode components and the covariate mode components; The component filtering module is used to perform a significance test by comparing the nonlinear complexity differences between the target variable modal components and covariate modal components and their respective corresponding phase randomized substitute data, eliminating components that fail the test, and retaining stable target variable modal components and stable covariate modal components. The coupling identification module is used to combine the stable target variable modal components and stable covariate modal components whose overlap between the time domain stable domain and the frequency domain resonance region exceeds the overlap threshold into a prediction unit. The time domain stable domain is determined by the instantaneous coupling strength of the sliding window, and the frequency domain resonance region is determined by the cross wavelet power spectrum. The prediction execution module is used to obtain the predicted values ​​of the stable target variable mode components by taking the stable covariate mode components in the prediction unit as input to the time series prediction network with embedded self-attention mechanism. The result reconstruction module is used to overlay and reconstruct the predicted values ​​of the modal components of all stable target variables to obtain water quality prediction results.

9. A water quality prediction device, characterized in that, The water quality prediction device includes: a processor, a memory, and a program or instructions stored in the memory and executable on the processor, wherein the program or instructions, when executed by the processor, implement the water quality prediction method as described in any one of claims 1-7.

10. A readable storage medium, characterized in that, A program or instruction is stored on a readable storage medium, which, when executed by a processor, implements the water quality prediction method as described in any one of claims 1-7.

Images