A method for calculating pollutant flux based on periodic and atlas wavelet combination

By combining periodicity and graph wavelet transform, a graph convolutional network is constructed, and a multilayer perceptron network is used to predict pollutant concentrations. This solves the problem of insufficient differentiation of pollutant flux frequency characteristics in existing technologies, and achieves more accurate pollutant flux prediction and environmental early warning.

CN122245494APending Publication Date: 2026-06-19MEIZHOU HYDROLOGY BRANCH OF GUANGDONG PROVINCIAL HYDROLOGY BUREAU

Patent Information

Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
MEIZHOU HYDROLOGY BRANCH OF GUANGDONG PROVINCIAL HYDROLOGY BUREAU
Filing Date
2026-03-09
Publication Date
2026-06-19

AI Technical Summary

Technical Problem

Existing technologies cannot effectively distinguish the frequency characteristics of flow rate and concentration when calculating pollutant flux, leading to distorted predictions under extreme weather conditions, making it impossible to accurately infer pollutant concentrations, and affecting the accuracy of environmental early warnings.

Method used

A method combining periodicity and graph wavelet analysis was adopted. By collecting historical monitoring data of the reservoir, the periodicity of timestamps was identified, a graph convolutional network was constructed, and a multilayer perceptron network was used to predict pollutant concentrations and calculate pollutant fluxes.

Benefits of technology

It enables detailed analysis of different frequency components, improves the accuracy of pollutant flux prediction and the effectiveness of environmental early warning, can distinguish between low-frequency and high-frequency changes in system dynamics, and provides frequency domain analysis to support unified modeling of complex systems.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN122245494A_ABST
    Figure CN122245494A_ABST
Patent Text Reader

Abstract

This invention relates to the field of environmental monitoring and water environment data analysis technology, and particularly to a method for calculating pollutant flux based on a combination of periodicity and spectral wavelet analysis. The method includes the following steps: collecting historical monitoring data from a reservoir, including sampling timestamps, pollutant concentrations, and flow rates; identifying the periodicity characteristics of the timestamps in the historical monitoring data and extracting periodic feature sequences using these characteristics; abstracting the periodic feature sequences into network nodes; calculating a similarity score matrix based on the network nodes; and constructing an adjacency matrix based on the similarity score matrix. This invention achieves multi-scale characterization and fusion modeling of the spatiotemporal evolution characteristics of pollutant concentrations based on environmental monitoring and water environment data analysis technology; and performs concentration prediction based on the fused global feature vector, thereby improving the dynamic characterization and prediction accuracy of pollutant flux calculation.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This invention relates to the field of environmental monitoring and water environment data analysis technology, and in particular to a method for calculating pollutant flux based on a combination of periodicity and spectral wavelet. Background Technology

[0002] Pollutant flux, referring to the total amount of pollutants passing through a cross-section of a river per unit time, is a core indicator for assessing water pollution load and environmental health. Its calculation formula is typically: Pollutant flux = Pollutant concentration (C) × Flow rate (Q) × Conversion factor. In practical water environment management and early warning work, a key challenge lies in prediction under asymmetric information: Flow rate (Q) data can usually be easily and continuously obtained through hydrological models or online monitoring stations; however, the determination of pollutant concentration (C) often requires manual sampling and laboratory analysis, resulting in data that is typically low-frequency, sparse, and expensive. Therefore, how to accurately infer the daily pollutant concentration (C) and thus calculate the pollutant flux when only the predicted daily flow rate (Q) is available has become a pressing technical problem in this field. Existing technologies completely ignore the frequency characteristics when analyzing the CQ relationship. They cannot distinguish systematic seasonal variations; for example, during the high-water season, flow rate and background concentration rise slowly and synchronously, which is a low-frequency coupling, while sudden flood events, such as heavy rain causing a surge in flow rate, result in a sharper and shorter initial scouring pulse in concentration, which is a high-frequency coupling. This frequency blind spot often causes models to be severely distorted when predicting pollution peaks caused by extreme weather, which is precisely the part that environmental early warning systems need to pay the most attention to. Summary of the Invention

[0003] Therefore, it is necessary for the present invention to provide a method for calculating pollutant flux based on the combination of periodicity and spectral wavelet to solve at least one of the above-mentioned technical problems.

[0004] To achieve the above objectives, a method for calculating pollutant flux based on a combination of periodicity and spectral wavelet analysis is proposed, comprising the following steps: Step S1: Collect historical monitoring data of the reservoir, including sampling timestamps, pollutant concentrations, and flow rates; identify the periodicity of timestamps in the historical monitoring data of the reservoir, and extract periodic feature sequences using the periodicity. Step S2: Abstract the periodic feature sequence into network nodes; calculate the similarity score matrix based on the network nodes; construct the adjacency matrix based on the similarity score matrix, and calculate the degree matrix using the adjacency matrix; calculate the graph Laplacian matrix based on the degree matrix; perform eigenvalue decomposition based on the graph Laplacian matrix to obtain the decomposition results; Step S3: Perform graph convolution based on the decomposition results and calculate wavelet coefficients; perform frequency band feature aggregation based on the wavelet coefficients to determine the global feature vector; calculate the weighted fusion feature value based on the global feature vector; Step S4: Input the weighted fusion feature values ​​into the multilayer perceptron network and output the pollutant concentration prediction value; calculate the pollutant flux based on the pollutant concentration prediction value.

[0005] The beneficial effects of this invention are as follows: (1) This invention introduces graph wavelet transform, which elevates the prediction model from analyzing only in the spatial domain of the graph to performing joint analysis in the frequency domain of the graph. This is a deep application of graph information processing theory, enabling the model to distinguish and utilize different frequency components in the system dynamics, thus achieving a new depth of analysis.

[0006] (2) The model constructed in this invention can distinguish between two essentially different but superficially similar dynamics: multiple cycles rising synchronously and smoothly (low frequency dominance); and multiple cycles oscillating violently and canceling each other out (high frequency dominance). By modeling different frequency bands separately, the model can capture more refined and essential physical laws, thereby improving prediction accuracy.

[0007] (3) The output of this invention is not only the predicted value, but also the frequency domain analysis of the system dynamics. By assigning weights to different frequency band features in the analysis model, it can be revealed whether the system is currently in a stable evolution state (high weight for low frequency) or a violent turbulent state (high weight for high frequency). This invention retains the structured decomposition of the periodicity of the time series, and at the same time captures the dynamic and multi-frequency interaction between these structured units through dynamic graphs and graph wavelet transforms, thus realizing unified modeling of complex systems. Attached Figure Description

[0008] Other features, objects, and advantages of the invention will become more apparent from the following detailed description of non-limiting embodiments with reference to the accompanying drawings: Figure 1 This is a schematic diagram of the steps in the method for calculating pollutant flux based on the combination of periodicity and spectral wavelet in this invention. Figure 2 This is a schematic diagram of the process for extracting periodic feature sequences based on periodic features in this invention; Figure 3 This is a schematic diagram of the structure of the multilayer perceptron network in this invention; The realization of the objective, functional features and advantages of the present invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

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

[0010] Furthermore, the accompanying drawings are merely illustrative of the invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and therefore repeated descriptions of them will be omitted. Some block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. These functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor methods and / or microcontroller methods.

[0011] It should be understood that although the terms "first," "second," etc., may be used herein to describe various units, these units should not be limited by these terms. These terms are used merely to distinguish one unit from another. For example, without departing from the scope of the exemplary embodiments, a first unit may be referred to as a second unit, and similarly, a second unit may be referred to as a first unit. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0012] To achieve the above objectives, please refer to Figures 1 to 3 This invention provides a method for calculating pollutant flux based on a combination of periodicity and spectral wavelet, the method comprising the following steps: Step S1: Collect historical monitoring data of the reservoir, including sampling timestamps, pollutant concentrations, and flow rates; identify the periodicity of timestamps in the historical monitoring data of the reservoir, and extract periodic feature sequences using the periodicity. In one embodiment, historical monitoring data of a reservoir for the past 5 years were continuously collected at 1-hour intervals. Data fields included timestamps, total phosphorus concentration (mg / L), ammonia nitrogen concentration (mg / L), and inflow rate (m³ / s). First, Fourier transform (FFT) analysis was performed on the timestamp sequence to identify significant periodic peaks (such as 24-hour daily cycles, 168-hour weekly cycles, and 8760-hour annual cycles). Major periodic components were screened based on spectral energy thresholds, and corresponding periodic feature subsequences were extracted using a sliding window (window length equal to one major period) to form the periodic feature sequence P(t).

[0013] In another embodiment, it is assumed that the data collection spans 3 years, totaling approximately 26,000 time sampling points, with a sampling interval of 2 hours. Periodic spectral analysis identifies the main periods as 12 hours, 24 hours, and 365 days. Assuming the main period energy accounts for 72% of the total spectral energy, three periodic feature subsequences are constructed, each with lengths of 12, 24, and 4380 points, respectively. The original concentration sequence is then rearranged and aligned according to the period to obtain the periodic feature matrix (assuming a dimension of 300×24).

[0014] Step S2: Abstract the periodic feature sequence into network nodes; calculate the similarity score matrix based on the network nodes; construct the adjacency matrix based on the similarity score matrix, and calculate the degree matrix using the adjacency matrix; calculate the graph Laplacian matrix based on the degree matrix; perform eigenvalue decomposition based on the graph Laplacian matrix to obtain the decomposition results; In one embodiment, each periodic feature subsequence is considered a network node, forming a total of N nodes. Similarity scores between nodes are calculated based on the Pearson correlation coefficient, resulting in an N×N similarity score matrix S. A threshold is then set. ,when At the node With nodes Establish edges between them and construct an adjacency matrix. Then the degree matrix is ​​calculated. ,in According to the formula Constructing the Laplacian matrix of an undirected graph and to Perform eigenvalue decomposition to obtain eigenvalues. and corresponding feature vectors .

[0015] In another embodiment, assuming there are 60 periodic nodes, similarity scores are calculated using Dynamic Time Warping (DTW) distance and normalized to the [0,1] interval. A threshold is set. This ultimately forms a graph structure with 60 nodes and 420 edges. An adjacency matrix A (60×60) is constructed, and the degree matrix D is calculated, with the maximum degree assumed to be 18. The graph Laplacian matrix L is then obtained and subjected to eigenvalue decomposition, yielding 60 eigenvalues. The first five eigenvalues ​​account for 65% of the total spectral energy, and their corresponding eigenvectors serve as the basis for graph frequency domain analysis.

[0016] Step S3: Perform graph convolution based on the decomposition results and calculate wavelet coefficients; perform frequency band feature aggregation based on the wavelet coefficients to determine the global feature vector; calculate the weighted fusion feature value based on the global feature vector; In one embodiment, a graph convolution kernel is constructed based on the eigenvalues ​​and eigenvectors of the graph Laplacian matrix. Graph convolution is performed on the periodic node features to obtain the convolution response matrix. The spectral wavelet transform function is then used... Construct a multi-scale graph wavelet basis, perform multi-scale decomposition on the graph convolution result, and calculate the corresponding wavelet coefficients. The wavelet coefficients are grouped and aggregated according to the low-frequency, mid-frequency, and high-frequency bands. The mean and variance of each band are calculated separately, and then concatenated to form a global feature vector. Weights are set based on the energy proportion of each frequency band. , , ,right Perform weighted fusion to obtain fusion feature values. .

[0017] In another embodiment, assuming scale parameters t=1, 2, and 4 are selected for three-scale wavelet decomposition, three sets of wavelet coefficient matrices (60×3) are obtained for each node. The average value of the low-frequency coefficients is 0.62, the mid-frequency coefficients are 0.35, and the high-frequency coefficients are 0.18. The weights are set to 0.5, 0.3, and 0.2, respectively, and a global fused feature vector (assuming a dimension of 1×128) is obtained after weighted fusion.

[0018] Step S4: Input the weighted fusion feature values ​​into the multilayer perceptron network and output the pollutant concentration prediction value; calculate the pollutant flux based on the pollutant concentration prediction value.

[0019] In one embodiment, the fused feature value Z is input into a multilayer perceptron (MLP) containing three hidden layers, with the following structure: input layer → 128 neurons → ReLU → 64 neurons → ReLU → 32 neurons → ReLU → output layer (1 neuron). The output is a predicted pollutant concentration. According to the flux formula (in The pollutant flux is calculated based on the flow rate, enabling the prediction of pollutant emission loads for future periods.

[0020] In another embodiment, the input feature dimension is assumed to be 128, and the number of training samples is 2000. The MLP network uses the Adam optimizer with a learning rate of 0.001 and converges after 100 training epochs. The mean square error (MSE) of the model's predicted concentration is 0.015 (mg / L)². Assuming a predicted concentration of 0.45 mg / L at a certain moment corresponds to a flow rate of 120 m³ / s, the instantaneous pollutant flux is calculated to be 54 mg / s. Further integration of the predicted flux over 24 consecutive hours yields the total daily pollution load.

[0021] It should be noted that you should refer to [link / reference]. Figure 2The figure illustrates the multi-period feature extraction process: starting from the original concentration history sequence, frequency domain analysis is first performed to identify the main periodic components, and then the periods at different scales (such as year, month, week) are divided into independent branches for processing; each branch reconstructs and transforms the corresponding periodic sequence, extracts fixed-dimensional feature representations, generates their own periodic feature vectors through convolutional or linear layers, and finally fuses the multi-scale periodic features to form a unified multi-period feature expression result.

[0022] It should be noted that you should refer to [link / reference]. Figure 3 The figure shows a schematic diagram of a multilayer perceptron network structure, which consists of an input layer, three hidden layers, and an output layer. Each layer contains several neuron nodes, and feature information is passed between nodes layer by layer through a fully connected manner. The input layer receives the original feature vector, which is then abstracted and mapped step by step after being weighted and nonlinearly transformed by the hidden layers. Finally, the output layer provides the prediction result or regression value, which reflects a typical feedforward fully connected neural network structure.

[0023] Preferably, step S1, which involves identifying the periodicity of timestamps in historical reservoir monitoring data and extracting periodic feature sequences using these periodic features, includes: Spectral analysis was performed based on historical reservoir monitoring data to identify periodic characteristics of timestamps, including annual, monthly, and weekly periodic characteristics; periodic feature sequences were then extracted based on these periodic characteristics.

[0024] In one embodiment, historical monitoring data of a reservoir for five consecutive years were collected, with a sampling interval of one hour. The data fields included timestamps, total nitrogen concentration, total phosphorus concentration, and inflow rate. The pollutant concentration sequences corresponding to the timestamps were detrended and normalized. Then, Fast Fourier Transform (FFT) was used to perform spectral analysis on the concentration time series to obtain a frequency-amplitude spectrum. Significant periodic components, such as corresponding frequencies, were identified based on the peak positions of the spectra. (Annual cycle) (Monthly cycle) (Periodic cycle). Further calculate the spectral energy percentage of each periodic component, and select periods with an energy percentage higher than a set threshold (e.g., 10%) as valid periodic features. Based on the identified annual, monthly, and weekly cycles, decompose the original concentration sequence into annual, monthly, and weekly periodic feature subsequences. The corresponding frequency band signals can be extracted using bandpass filtering or frequency domain reconstruction methods, and the data can be rearranged according to the period length to form a periodic feature matrix. For example, align the data from 8760 hours of data per year by year to form a matrix. A matrix is ​​formed with a 720-hour cycle. A matrix is ​​formed with 168 hours as a cycle. , as a periodic feature sequence.

[0025] Of particular importance is the performance of spectral analysis based on historical reservoir monitoring data to identify periodic characteristics of timestamps, including: Historical reservoir monitoring data were used to construct a continuous time series in chronological order of timestamps. Missing value imputation, outlier removal, and mean normalization were performed on the time series to form a regular sampling sequence. In one embodiment, multi-dimensional historical data, including water level, inflow, outflow, rainfall, and reservoir temperature, are read from the reservoir monitoring system and sorted from earliest to latest according to timestamps to construct a continuous time series under a unified time axis. For data with inconsistent sampling intervals, time alignment is performed using linear interpolation or spline interpolation to ensure that all variables correspond to the same sampling interval, for example, a fixed sampling period of 1 hour. For missing values ​​in the time series, different strategies are adopted based on the length of the missing interval: when the consecutive missing length is less than 3 sampling points, adjacent point linear interpolation is used; when the consecutive missing length is greater than 3 sampling points, imputation based on the historical mean of the same period is used. Outlier removal can be performed using a sliding window statistical method, calculating the mean and standard deviation within a 24-hour window. When a data point deviates from the window mean by more than 3 times the standard deviation, it is replaced or smoothed. After data cleaning, mean normalization is performed on each variable to make the mean of the processed regular sampling series close to 0 and the numerical range stable at a unified scale, providing a consistent data input structure for subsequent frequency analysis.

[0026] Identify the frequency distribution of the regular sampling sequence to record the frequency distribution results; identify the annual frequency component in the frequency distribution results to determine the annual cycle characteristics; identify the monthly frequency component in the frequency distribution results to determine the monthly cycle characteristics; In one embodiment, assuming that frequency domain analysis is performed on 1825 daily data points to obtain the frequency distribution results, a significant peak with an amplitude of 0.85 appears near frequency 0.0027 (approximately 1 / 365), which is identified as the annual cycle component; a secondary peak with an amplitude of 0.42 appears near frequency 0.033 (approximately 1 / 30), which is identified as the monthly cycle component. The annual cycle amplitude accounts for approximately 38% of the total spectral energy, and the monthly cycle amplitude accounts for approximately 17%. Phase angle calculations indicate that the annual cycle peak occurs in mid-July each year, and the monthly cycle peak occurs between the 5th and 8th of each month.

[0027] Identify the weekly frequency component in the frequency distribution results and determine the weekly periodic characteristics; mark the annual periodic characteristics, monthly periodic characteristics, and weekly periodic characteristics as the periodic characteristics of the timestamp.

[0028] In one embodiment, the frequency distribution results are further searched for frequency components corresponding to a weekly scale. When the sampling interval is 1 day, the frequency is 1 / 7 ( The components near this frequency correspond to periodic variations. By setting a preset bandwidth around this frequency, the local maximum amplitude points are filtered out, and their amplitude and phase information are extracted to determine the periodic characteristic parameters. The annual, monthly, and weekly periodic characteristics are mapped back to the time domain, and corresponding periodic marker parameters are generated for each timestamp. For example, a three-dimensional periodic feature vector is constructed, including the annual phase, monthly phase, and weekly phase.

[0029] In another embodiment, it is assumed that a weekly frequency peak with an amplitude of 0.26 is detected near a frequency of 0.142 (approximately 1 / 7), accounting for about 9% of the total spectral energy. This identifies three main periodic components: an annual periodic amplitude of 0.85, a monthly periodic amplitude of 0.42, and a weekly periodic amplitude of 0.26. The phase angles corresponding to these three types of periodic components are denoted as follows: , This is then converted into a time-domain periodic label, and a periodic feature vector corresponding to each timestamp is constructed, such as [annual phase value 0.65, monthly phase value 0.42, weekly phase value 0.88] etc.

[0030] Preferably, the periodic feature extraction periodic feature sequence includes: To address the annual cycle characteristics, historical reservoir monitoring data is divided into two-dimensional images with a 365-day cycle. A pre-defined convolutional neural network is then used to perform convolution operations on the two-dimensional images to extract annual variation feature maps. Based on the annual variation feature maps, the annual cycle feature sequence is determined. In one embodiment, daily average pollutant concentration data and corresponding flow data for a reservoir over five consecutive years were acquired, totaling 1825 sampling points. Using 365 days as a complete cycle, the historical monitoring data was segmented and rearranged, constructing a 365×M two-dimensional matrix for each year, where M represents the number of monitoring indicators (e.g., if there are three indicators: total nitrogen, total phosphorus, and flow, then M=3). This two-dimensional matrix was then treated as a single-channel or multi-channel image and input into a pre-defined convolutional neural network (CNN). The CNN includes two convolutional layers: the first layer has a 3×3 convolutional kernel, 32 channels, and a stride of 1; the second layer has a 3×3 convolutional kernel, 64 channels, and features are compressed using the ReLU activation function and a max-pooling layer. After convolution operations, an annual variation feature map is output. Its size is, for example, 365×64. Global average pooling is performed on the feature map along the channel dimension to obtain a one-dimensional feature vector of length 64, which serves as the annual cycle feature sequence.

[0031] In another embodiment, assume that daily average data for a reservoir over three years is obtained, totaling 1095 sampling points. After dividing the data into 365 days, three 365×2 two-dimensional matrices are obtained (assuming there are two types of monitoring indicators). These matrices are input into a CNN network containing three convolutional layers (with 16, 32, and 64 kernels respectively), outputting an annual variation feature map of size 365×64. The average response values ​​of the top 10 main channels in this feature map are calculated as [0.82, 0.75, 0.69, 0.64, 0.58, 0.55, 0.51, 0.48, 0.44, 0.40]. After average pooling along the 365-dimensional time axis, an annual periodic feature vector of length 64 is obtained, assuming the largest eigenvalue is 1.25 and the smallest is 0.12.

[0032] To address the monthly cycle characteristics, historical reservoir monitoring data were filled and reconstructed according to a 30-day cycle to form a two-dimensional structure tensor for the monthly cycle; the monthly cycle characteristic sequence was then determined based on the two-dimensional structure tensor for the monthly cycle. In one embodiment, historical daily average monitoring data is divided into 30-day periods. Periods shorter than 30 days are padded with zeros or the mean, resulting in a 30×M two-dimensional matrix for each monthly period. Multiple periods are then stacked to form a three-dimensional tensor. The structure has dimensions K×30×M, where K represents the number of monthly cycles. One-dimensional convolution or average pooling operations are performed on the two-dimensional structure tensor of the monthly cycles along the time dimension (30 days) to extract local variation features for each cycle, and a fixed-length feature vector is output for each cycle. All cycle feature vectors are concatenated in chronological order to form a monthly cycle feature sequence. .

[0033] In another embodiment, assuming 730 daily average sampling points are obtained over two years, this data is padded and divided into 25 monthly periods, each period being reshaped into a 30×3 two-dimensional matrix. The monthly period tensor is constructed with dimensions of 25×30×3. A 1×3 convolution is performed on each period to extract inter-channel coupling features, and average pooling is performed over the 30-day dimension to obtain a feature vector of length 32 for each month. Assuming the mean of the feature vectors for months 1 to 5 are [0.42, 0.55, 0.61, 0.58, 0.49] and the standard deviation ranges from [0.08, 0.12], a monthly periodic feature sequence of length 25×32 is formed.

[0034] To address the periodic characteristics, historical monitoring data of the reservoir were filled and reconstructed into a periodic two-dimensional image tensor with a 7-day cycle, and the periodic characteristic sequence was determined from the periodic two-dimensional image tensor.

[0035] In one embodiment, historical daily average monitoring data is segmented and rearranged according to a 7-day cycle. For segments shorter than 7 days, padding with previous values ​​is used to construct a 7×M periodic two-dimensional matrix. Multiple consecutive periodic data are stacked to form a periodic image tensor. The dimensions are L×7×M. A two-dimensional convolution operation is performed on the periodic two-dimensional image tensor with a kernel size of 2×2 and 16 channels, and the ReLU activation function is used to extract local time-index coupling features. Pooling is then performed on the convolution result along the time axis to output a periodic feature vector of length 16, which is then concatenated according to the period sequence to form a periodic feature sequence. .

[0036] In another embodiment, assuming 365 sampling points from one year of data are acquired, after padding, 53 weekly cycles are formed, each cycle being reshaped into a 7×2 two-dimensional matrix. The weekly cycle tensor is constructed with dimensions of 53×7×2. Two convolutional layers are performed on each cycle (8 and 16 channels respectively), outputting a feature map of size 7×16. Max pooling is then applied to the time dimension to obtain a feature vector of length 16. Assuming the maximum response value for a particular week is 0.93 and the minimum is 0.15, the overall weekly cycle feature sequence has a mean of 0.57 and a standard deviation of 0.11.

[0037] Preferably, determining the annual cycle characteristic sequence based on the annual variation characteristic map includes: Normalize and perform nonlinear activation on the annual variation feature map to record the enhanced annual feature response map; perform spatial dimension compression on the enhanced annual feature response map to convert it into a one-dimensional feature vector; In one embodiment, after obtaining the annual variation feature map After (with a size of 365×C, where C is the number of channels, e.g., C=64), the feature values ​​of each channel are normalized. Specifically, min-max normalization or batch normalization is performed on each channel to map the channel feature values ​​to the [0,1] interval or the standard normal distribution interval. The normalization result is then input into a non-linear activation function (e.g., ReLU or LeakyReLU) to enhance high-response regions and suppress low-response noise, resulting in an enhanced annual feature response map. To achieve spatial compression, global average pooling or max pooling is performed on the enhanced annual feature response map in the spatial dimension (365 days). Assuming global average pooling is used, averaging over each channel at 365 time steps yields a one-dimensional vector of length C. (For example, the dimension is 1×64).

[0038] In another embodiment, assume the annual variation feature map size is 365×128. Calculate the mean for each channel. with standard deviation For example, the mean of the first 5 channels is [0.42, 0.55, 0.61, 0.48, 0.37], and the standard deviation is [0.12, 0.15, 0.10, 0.09, 0.11]. After standardization, the signal is input into the Sigmoid activation function to obtain an enhanced annual feature response map, with an assumed overall mean of 0.63, a maximum response value of 0.94, and a minimum of 0.08. Subsequently, global average pooling is performed on the 365-dimensional time axis to obtain a one-dimensional feature vector of length 128, where the values ​​of the first 10 dimensions are assumed to be [0.71, 0.68, 0.65, 0.62, 0.59, 0.57, 0.54, 0.52, 0.49, 0.47].

[0039] Of particular importance is the normalization and nonlinear activation processing of the annual variation feature map, which records the enhanced annual feature response map, including: Calculate the channel mean in the annual change feature map, and use the channel mean to scale and shift the feature response values ​​in the annual change feature map to form a feature map; In one embodiment, the annual variation feature map is assumed to be 256×256×8 pixels in size, containing 8 channels. The average values ​​of the 8 channels are statistically obtained as [0.42, 0.37, 0.51, 0.46, 0.39, 0.44, 0.48, 0.41]. For the first channel, 0.42 is subtracted from each of the 65536 pixels, and the image is scaled by a ratio of 1.2, then shifted by 0.05. For the second channel, 0.37 is subtracted, and the image is scaled by a ratio of 1.1 and shifted by 0.04. The remaining channels are scaled and shifted according to their respective average values. After processing, the numerical distribution of the data in each channel is adjusted from the original range [0.10, 0.92] to approximately [−0.35, 0.78], and the amplitude differences between channels converge significantly, resulting in a normalized feature map.

[0040] The feature map is input into a preset nonlinear activation function to perform numerical mapping, and the mapping result is recorded. The mapping result is then recombined according to the channel structure of the feature map to form an enhanced annual feature response map.

[0041] In one embodiment, the output feature map size is assumed to remain 256×256×8, with a data range of approximately [−0.35, 0.78]. All channel data are input into a nonlinear activation function for point-by-point mapping, resulting in a numerical range between [0, 1]. Taking the third channel as an example, its mean value before mapping is 0.12, and after mapping, the mean value increases to 0.56, with the extreme values ​​changing from −0.28 and 0.73 to 0.21 and 0.92. After mapping all eight channels, eight 256×256 matrices are generated, which are then stacked and combined according to the original channel index order to form a new 256×256×8 tensor structure. This enhanced annual feature response map, while maintaining the same spatial resolution, exhibits a more concentrated channel response distribution and a more prominent representation of changing regions.

[0042] Generate annual feature representation vectors from one-dimensional feature vectors; arrange the annual feature representation vectors in chronological order to form an annual cyclical feature sequence.

[0043] In one embodiment, it is assumed that annual feature vectors for four years are obtained, each vector having a dimension of 128. These vectors are mapped to 64 dimensions through a fully connected layer, resulting in four annual feature representation vectors. It is assumed that the means of these four vectors are [0.58, 0.62, 0.67, 0.71], and the standard deviations are [0.09, 0.11, 0.10, 0.12], showing an increasing trend year by year. These vectors are arranged in chronological order to form a 4×64 annual periodic feature sequence matrix. The value of a certain feature dimension over the four years is assumed to be [0.45, 0.49, 0.53, 0.60], reflecting a year-on-year increase in the pollution response corresponding to that feature.

[0044] Preferably, step S2, calculating the similarity score matrix based on network nodes, includes: Network nodes are mapped to query vectors, key vectors, and value vectors, respectively; matrix multiplication is performed on the query vectors and key vectors to calculate the dot product, and an initial similarity score matrix is ​​generated based on the dot product. In one embodiment, assuming the number of nodes N=128 and the original feature dimension d=64, the following settings are made: , The dimensions of the weight matrices are then respectively... , , After mapping, we get , , Perform matrix multiplication. Later obtained The initial similarity score matrix. For example, the value of the element in row 10, column 25 might be 6.84, the value in row 10, column 26 might be 2.17, and the value in row 10, column 27 might be -1.35.

[0045] The initial similarity score matrix is ​​scaled by the square root of the value vector to obtain the scaled similarity score matrix; the scaled similarity score matrix is ​​then normalized to obtain the final similarity score matrix.

[0046] In one embodiment, for the initial similarity score matrix According to the dimension of the value vector The square root is scaled, i.e., a scaling matrix is ​​constructed. .when When the value is large, scaling using a square root factor can suppress excessively large amplitudes in the dot product, making the numerical distribution more stable in subsequent normalization processes. Then, normalization is performed on each row vector, and a similarity score matrix in probability distribution form is generated using exponential mapping and row-level summation. Specifically, for the elements of the i-th row... Perform an exponential transformation and divide by the sum of all exponential values ​​in that row, so that the sum of the elements in each row is 1, thus forming the final similarity score matrix.

[0047] Preferably, step S2, which involves constructing an adjacency matrix based on the similarity score matrix and calculating the degree matrix using the adjacency matrix, includes: A threshold filtering process is performed on the similarity score matrix, retaining similarity scores greater than a preset threshold as valid connection weights to generate an adjacency matrix; the degree matrix of the corresponding node is obtained by summing the elements of each row of the adjacency matrix.

[0048] In one embodiment, assuming there are 6 nodes, the weights between some nodes in the similarity score matrix are 0.82, 0.34, 0.67, 0.75, 0.58, 0.31, etc. A threshold of 0.50 is set, retaining weights greater than 0.50 such as 0.82, 0.67, 0.75, and 0.58, while setting the rest to zero, such as 0.34 and 0.31. After filtering, the first node has valid connections with the second and fifth nodes; the fourth node has valid connections with the second, third, and fifth nodes; the sixth node is an isolated node because all its connection weights are less than 0.50. The sums are calculated for each row; for example, the sum of valid connection weights for the first node is 1.38, for the second node it is 1.45, for the third node it is 0.73, for the fourth node it is 2.02, for the fifth node it is 1.13, and for the sixth node it is 0. Fill the above values ​​into the diagonal positions of the degree matrix to form the corresponding degree matrix structure.

[0049] Preferably, the degree matrix of the corresponding node is obtained by summing the elements of each row of the adjacency matrix, including: The connection weights corresponding to the elements in each row of the adjacency matrix are accumulated sequentially and recorded as node degree values; according to the node order of the node degree values, the node degree values ​​are filled into the adjacency matrix sequentially, and the diagonal positions are determined. In one embodiment, assume the adjacency matrix has a 6×6 structure, where the cumulative non-zero weights of each row are: 1.38 for the first node, 1.45 for the second node, 0.73 for the third node, 2.02 for the fourth node, 1.13 for the fifth node, and 0 for the sixth node. Following the node order, these values ​​are filled into positions (1,1), (2,2), (3,3), (4,4), (5,5), and (6,6) of the matrix. At this point, the diagonal elements of the matrix are 1.38, 1.45, 0.73, 2.02, 1.13, and 0, respectively, while the values ​​in the remaining positions remain unchanged from the original adjacency structure.

[0050] All positions in the adjacency matrix except the diagonal positions are assigned zero to form a degree matrix.

[0051] In one embodiment, after filling the diagonal elements, a masking process is performed on the matrix, retaining only the main diagonal elements and setting all off-diagonal elements to zero. This can be achieved by checking each element's index to determine if the rows and columns are equal, thus forming a matrix structure containing only diagonal elements; this matrix is ​​the degree matrix.

[0052] In another embodiment, continuing with the above six nodes as an example, after filling the diagonal with 1.38, 1.45, 0.73, 2.02, 1.13, and 0, all off-diagonal weights in the original adjacency matrix, such as 0.82, 0.67, 0.75, and 0.58, are cleared to zero. The final matrix retains only the above six values ​​in the diagonal positions, with the remaining 30 off-diagonal elements all being 0, forming a 6×6 degree matrix structure, which is used in subsequent graph structure normalization or Laplacian matrix construction processes.

[0053] Preferably, in step S2, eigenvalue decomposition is performed based on the graph Laplacian matrix to obtain the decomposition results, including: The lowest and highest frequency DC components are determined based on the Laplace matrix. Three Gaussian wavelet kernels are then set within the range of these components: a low-frequency wavelet kernel, a mid-frequency wavelet kernel, and a high-frequency wavelet kernel. The low-frequency wavelet kernel is set with a center parameter of [0, 0.5] and a bandwidth parameter of (0, 0.5]. The mid-frequency wavelet kernel is set with a center parameter of (0.5, 1.5) and a bandwidth parameter of (0, 0.7]. The high-frequency wavelet kernel is set with a center parameter of [1.5,2] and a bandwidth parameter of (0,0.5]. In one embodiment, the normalized graph Laplacian matrix is ​​calculated for the constructed graph structure, and eigenvalue decomposition is performed on the matrix to obtain an ascending sequence of eigenvalues ​​λ0≤λ1≤…≤λ n Where λ0 corresponds to the lowest frequency DC component, which is usually close to 0, λ n Corresponding to the highest frequency components. Using λ0 and λ n Based on the determined spectral range, the entire spectral domain is divided into three sub-ranges, and three Gaussian wavelet kernel functions are constructed within each corresponding range. The center parameter of the low-frequency wavelet kernel is set in the range [0, 0.5], and its bandwidth is limited to (0, 0.5], giving it a high response near λ0. The center parameter of the mid-frequency wavelet kernel is set in the range (0.5, 1.5) and its bandwidth is limited to (0, 0.7], covering the main structural changes in the middle of the spectral domain. The center parameter of the high-frequency wavelet kernel is set in the range [1.5, 2] and its bandwidth is limited to (0, 0.5], used to emphasize the areas of drastic local changes in the graph structure. The three Gaussian wavelet kernels are distributed around their respective center frequencies, and their response curves exhibit continuous decay characteristics in the spectral domain, thus forming a multi-scale filtering structure covering the entire frequency range.

[0054] In another embodiment, assuming a graph contains 256 nodes, the calculated eigenvalue range of the graph's Laplacian matrix is ​​approximately [0.002, 1.98]. Using this range as the range for the lowest frequency DC component and the highest frequency component, the center of the low-frequency wavelet kernel is set to 0.25, and the bandwidth to 0.4; the center of the mid-frequency wavelet kernel is set to 1.0, and the bandwidth to 0.6; and the center of the high-frequency wavelet kernel is set to 1.75, and the bandwidth to 0.35. The three Gaussian wavelet kernels constructed accordingly cover approximately [0, 0.65], [0.4, 1.6], and [1.4, 2] intervals on the spectral axis, respectively. By statistically analyzing the eigenvalue distribution density, if the proportion of low-frequency eigenvalues ​​is approximately 30%, mid-frequency approximately 45%, and high-frequency approximately 25%, then the above parameter configuration can form a relatively balanced spectral response distribution across different frequency bands.

[0055] The low-frequency wavelet kernel, mid-frequency wavelet kernel, and high-frequency wavelet kernel are mapped to the preset spectral domain space, and the frequency band filtering results are recorded. Eigenvalue decomposition is performed based on the frequency band filtering results to obtain the decomposition results.

[0056] In one embodiment, assuming the input image signal has a dimension of 256×1, a 256-dimensional spectral coefficient vector is obtained after spectral projection. After filtering in the low-frequency band, approximately 80 principal non-zero response coefficients are retained, approximately 120 in the mid-frequency band, and approximately 60 in the high-frequency band. A 256×256 covariance matrix is ​​constructed for the recovered vertex domain signal of each frequency band, and the top 10 principal eigenvalues ​​and their corresponding eigenvectors are extracted. It is assumed that the cumulative percentage of the top 3 eigenvalues ​​in the low-frequency band reaches 72%, the top 5 in the mid-frequency band reaches 68%, and the top 2 in the high-frequency band reaches 65%. By comparing the decomposition results of different frequency bands, a multi-scale spectral decomposition structure can be formed.

[0057] Preferably, the low-frequency wavelet kernel, mid-frequency wavelet kernel, and high-frequency wavelet kernel are mapped to a preset spectral domain space, and the frequency band filtering results are recorded, including: The low-frequency wavelet kernel is loaded into the low-frequency segment of the spectral space, and the low-frequency filtering result is recorded; the mid-frequency wavelet kernel is loaded into the mid-frequency segment of the spectral space, and the mid-frequency filtering result is recorded. In one embodiment, assuming a graph contains 300 nodes, 300 eigenvalues ​​are obtained after eigenvalue decomposition, ranging from approximately [0.01, 2.05]. The components corresponding to the first 90 eigenvalues ​​are classified as low-frequency segments, and the middle 120 are classified as mid-frequency segments. After loading a wavelet kernel in the low-frequency segment, approximately 85 components in the low-spectrum coefficients show significant responses. After recovery to the vertex domain, a low-frequency filtered signal is obtained, with an overall amplitude range of approximately [-0.8, 0.9]. After loading a wavelet kernel in the mid-frequency segment, approximately 110 components in the mid-spectrum coefficients produce responses, and the amplitude range of the recovered mid-frequency filtered signal is approximately [-0.5, 0.6]. By statistically analyzing the energy proportions, the energy of the low-frequency filtered result accounts for approximately 48% of the total energy of the original signal, and the energy of the mid-frequency filtered result accounts for approximately 34%.

[0058] The high-frequency wavelet kernel is loaded into the high-frequency region of the spectral domain, and the high-frequency filtering results are recorded. The low-frequency filtering results, mid-frequency filtering results, and high-frequency filtering results are encapsulated into a frequency band filtering result.

[0059] In one embodiment, after low-frequency and mid-frequency filtering, a high-frequency wavelet kernel is loaded into the high-frequency segment of the spectral domain. The high-frequency segment typically corresponds to the latter part of the eigenvalue sequence, and its spectral components reflect abrupt changes, edges, or local anomalies in the graph structure. High-frequency filtered spectral coefficients are obtained by applying high-frequency wavelet weighting to the spectral coefficients of this segment, and then inversely mapped to the vertex domain using the eigenvector matrix to form the high-frequency filtering result. The low-frequency, mid-frequency, and high-frequency filtering results are integrated using a unified data structure, such as stacking them in a three-channel format or indexing and encapsulating them using frequency band labels, to form a complete set of frequency band filtering results.

[0060] In another embodiment, it is assumed that the high-frequency band corresponds to the last 90 eigenvalues. After loading the high-frequency wavelet kernel, approximately 70 spectral components produce significant responses. After recovery to the vertex domain, the amplitude range of the high-frequency filtered signal is approximately [-0.3, 0.4], accounting for 18% of the total energy of the original signal. At this point, the energy proportions of the three types of filtering results are 48% for low frequency, 34% for mid frequency, and 18% for high frequency. The three types of results are constructed into a 300×3 matrix structure according to the node order, with each column corresponding to the low-frequency, mid-frequency, and high-frequency filtered components, respectively; or encapsulated as a band filtering data object containing three sub-signal arrays, along with the energy statistics of each band and the number of effective response nodes, thereby forming a complete band filtering result output structure.

[0061] Preferably, in step S4, the weighted fusion feature values ​​are input into the multilayer perceptron network, and the predicted pollutant concentration values ​​are output, including: The weighted fused feature values ​​are input into the input layer of the multilayer perceptron network; the feature vectors of the weighted fused feature values ​​are then passed to each neuron node in the first hidden layer through the input layer. In one embodiment, the weighted fusion feature values ​​obtained in the preceding steps are vectorized to construct a fixed-dimensional input feature vector. For example, multi-source monitoring parameters, frequency domain components, or time-series statistics are concatenated to form a one-dimensional feature vector X. This feature vector is then input to the input layer of a multilayer perceptron network (MLP). The input layer does not perform nonlinear operations; it only aligns and batch normalizes the feature dimensions and, according to the network structure connections, transmits each feature component to the respective neuron node of the first hidden layer. During data transmission, a fully connected mapping relationship is established between each input feature component and the first hidden layer neuron, forming an initial weight matrix structure. The input layer outputs a standardized feature representation, which is used for nonlinear transformations in subsequent layers.

[0062] In another embodiment, it is assumed that the fused feature vector has a dimension of 64, corresponding to information such as temperature, humidity, wind speed, historical concentration sequence statistics, and frequency band characteristics in air quality monitoring. This 64-dimensional vector is input into a first hidden layer containing 128 neurons to construct a 64×128 weight matrix. The input layer output values ​​are standardized and distributed in the interval [-1, 1], serving as the initial input data for the first hidden layer.

[0063] In the first hidden layer, each neuron node performs weighted summation and non-linear activation on the feature vector to obtain the first hidden feature layer; the first hidden feature layer is then passed to the second hidden layer, and the weighted mapping and non-linear activation operations are repeated to obtain the second hidden feature layer. In one embodiment, each neuron node in the first hidden layer receives all feature components from the input layer and performs a weighted summation based on the corresponding weight parameters, while simultaneously adding a bias term to form a linear combination result. A non-linear activation function, such as ReLU or LeakyReLU, is then applied to the linear combination result for mapping, to enhance the non-linearity of feature representation. This processing yields the first hidden feature vector H1. The first hidden feature vector is then passed as input to the second hidden layer. The second hidden layer also employs a fully connected structure. The weighted summation and bias stacking operations are performed, and then mapped using a non-linear activation function to form the second layer of hidden features. The features in this layer can be compressed or expanded in dimensionality compared to the first layer to adjust the network's representational power and complexity.

[0064] In another embodiment, assume the first hidden layer has 128 neurons and the second hidden layer has 64 neurons. After weighted mapping, the first hidden layer outputs 128-dimensional hidden features, where approximately 60% of the neurons output positive values ​​upon activation. The second hidden layer maps these 128-dimensional features to construct a 128×64 weight matrix, outputting a 64-dimensional hidden feature vector. Assume that after training, the mean of the second layer's hidden features is approximately 0.35 and the variance is approximately 0.12, representing the network's intermediate abstract representation of the input features.

[0065] The second hidden feature is passed to the third hidden layer, and the output is input to the output layer to determine the predicted pollutant concentration.

[0066] In one embodiment, the second layer of hidden features As input to the third hidden layer, the third hidden layer performs a further weighted mapping and non-linear activation process on it to form the third layer hidden features. The third layer of hidden features is then passed to the output layer. The output layer typically employs a linear mapping structure to generate continuous predictions. The output layer can be configured with a single node or multiple nodes, corresponding to single pollutant prediction or multi-pollutant joint prediction, respectively.

[0067] In another embodiment, assuming the third hidden layer contains 32 neurons, the 64-dimensional hidden features of the second layer are mapped to a 32-dimensional feature representation; the output layer has one neuron to output the predicted PM2.5 concentration value. After training, the predicted concentration value in the test samples ranges from approximately [value missing]. The average absolute error between the actual and actual monitoring values ​​is approximately 4.7 μg / m³. If multiple output nodes are set, for example, three nodes corresponding to PM2.5, PM10, and NO2 concentrations respectively, the output vector dimension is 3-dimensional, and the predicted values ​​can fall within the intervals of [20,140], [30,180], and [10,90] respectively, which can be used to construct multi-pollutant concentration prediction results.

[0068] Therefore, the embodiments should be considered as exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of the equivalents of the application are intended to be included within the invention.

[0069] The above description is merely a specific embodiment of the present invention, enabling those skilled in the art to understand or implement the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the present invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features of the invention herein.

Claims

1. A method for calculating pollutant flux based on a combination of periodicity and spectral wavelet, characterized in that, Includes the following steps: Step S1: Collect historical monitoring data of the reservoir, including sampling timestamps, pollutant concentrations, and flow rates; identify the periodicity of timestamps in the historical monitoring data of the reservoir, and extract periodic feature sequences using the periodicity. Step S2: Abstract the periodic feature sequence into network nodes; Calculate the similarity score matrix based on the network nodes; Construct an adjacency matrix based on the similarity score matrix, and then use the adjacency matrix to calculate the degree matrix; Calculate the graph Laplacian matrix based on the degree matrix; Perform eigenvalue decomposition based on the graph Laplacian matrix to obtain the decomposition results; Step S3: Perform graph convolution based on the decomposition results and calculate wavelet coefficients; Frequency band feature aggregation is performed based on wavelet coefficients to determine the global feature vector; Calculate the weighted fusion feature value based on the global feature vector; Step S4: Input the weighted fusion feature values ​​into the multilayer perceptron network and output the pollutant concentration prediction value; calculate the pollutant flux based on the pollutant concentration prediction value.

2. The method for calculating pollutant flux based on the combination of periodicity and spectral wavelet as described in claim 1, characterized in that, Step S1 involves identifying the periodicity of timestamps in historical reservoir monitoring data and extracting periodic feature sequences using these periodic features, including: Spectral analysis was performed based on historical reservoir monitoring data to identify periodic characteristics of timestamps, including annual, monthly, and weekly periodic characteristics; periodic feature sequences were then extracted based on these periodic characteristics.

3. The method for calculating pollutant flux based on the combination of periodicity and spectral wavelet as described in claim 2, characterized in that, Extracting periodic feature sequences based on periodic characteristics includes: To address the annual cycle characteristics, historical reservoir monitoring data is divided into two-dimensional images with a 365-day cycle. A pre-defined convolutional neural network is then used to perform convolution operations on the two-dimensional images to extract annual variation feature maps. Based on the annual variation feature maps, the annual cycle feature sequence is determined. To address the monthly cycle characteristics, historical reservoir monitoring data were filled and reconstructed according to a 30-day cycle to form a two-dimensional structure tensor for the monthly cycle; the monthly cycle characteristic sequence was then determined based on the two-dimensional structure tensor for the monthly cycle. To address the periodic characteristics, historical monitoring data of the reservoir were filled and reconstructed into a periodic two-dimensional image tensor with a 7-day cycle, and the periodic characteristic sequence was determined from the periodic two-dimensional image tensor.

4. The method for calculating pollutant flux based on the combination of periodicity and spectral wavelet as described in claim 3, characterized in that, The annual cycle characteristic sequence determined based on the annual variation characteristic map includes: Normalize and perform nonlinear activation on the annual variation feature map to record the enhanced annual feature response map; perform spatial dimension compression on the enhanced annual feature response map to convert it into a one-dimensional feature vector; Generate annual feature representation vectors from one-dimensional feature vectors; arrange the annual feature representation vectors in chronological order to form an annual cyclical feature sequence.

5. The method for calculating pollutant flux based on the combination of periodicity and spectral wavelet as described in claim 1, characterized in that, Step S2, which calculates the similarity score matrix based on network nodes, includes: Network nodes are mapped to query vectors, key vectors, and value vectors, respectively; matrix multiplication is performed on the query vectors and key vectors to calculate the dot product, and an initial similarity score matrix is ​​generated based on the dot product. The initial similarity score matrix is ​​scaled by the square root of the value vector to obtain the scaled similarity score matrix; the scaled similarity score matrix is ​​then normalized to obtain the final similarity score matrix.

6. The method for calculating pollutant flux based on the combination of periodicity and spectral wavelet as described in claim 1, characterized in that, Step S2 involves constructing an adjacency matrix based on the similarity score matrix and calculating the degree matrix using the adjacency matrix, including: A threshold filtering process is performed on the similarity score matrix, retaining similarity scores greater than a preset threshold as valid connection weights to generate an adjacency matrix; the degree matrix of the corresponding node is obtained by summing the elements of each row of the adjacency matrix.

7. The method for calculating pollutant flux based on the combination of periodicity and spectral wavelet as described in claim 6, characterized in that, Summing the elements of each row of the adjacency matrix yields the degree matrix of the corresponding node, including: The connection weights corresponding to the elements in each row of the adjacency matrix are accumulated sequentially and recorded as node degree values; according to the node order of the node degree values, the node degree values ​​are filled into the adjacency matrix sequentially, and the diagonal positions are determined. All positions in the adjacency matrix except the diagonal positions are assigned zero to form a degree matrix.

8. The method for calculating pollutant flux based on the combination of periodicity and spectral wavelet as described in claim 1, characterized in that, In step S2, eigenvalue decomposition is performed based on the graph Laplacian matrix to obtain the decomposition results, including: The lowest and highest frequency DC components are determined based on the Laplace matrix. Three Gaussian wavelet kernels are then set within the range of these components: a low-frequency wavelet kernel, a mid-frequency wavelet kernel, and a high-frequency wavelet kernel. The low-frequency wavelet kernel is set with a center parameter of [0, 0.5] and a bandwidth parameter of (0, 0.5]. The mid-frequency wavelet kernel is set with a center parameter of (0.5, 1.5) and a bandwidth parameter of (0, 0.7]. The high-frequency wavelet kernel is set with a center parameter of [1.5,2] and a bandwidth parameter of (0,0.5]. The low-frequency wavelet kernel, mid-frequency wavelet kernel, and high-frequency wavelet kernel are mapped to the preset spectral domain space, and the frequency band filtering results are recorded. Eigenvalue decomposition is performed based on the frequency band filtering results to obtain the decomposition results.

9. The method for calculating pollutant flux based on the combination of periodicity and spectral wavelet as described in claim 8, characterized in that, The low-frequency wavelet kernel, mid-frequency wavelet kernel, and high-frequency wavelet kernel are mapped to a preset spectral domain space, and the frequency band filtering results are recorded as follows: The low-frequency wavelet kernel is loaded into the low-frequency segment of the spectral space, and the low-frequency filtering result is recorded; the mid-frequency wavelet kernel is loaded into the mid-frequency segment of the spectral space, and the mid-frequency filtering result is recorded. The high-frequency wavelet kernel is loaded into the high-frequency region of the spectral domain, and the high-frequency filtering results are recorded. The low-frequency filtering results, mid-frequency filtering results, and high-frequency filtering results are encapsulated into a frequency band filtering result.

10. The method for calculating pollutant flux based on the combination of periodicity and spectral wavelet as described in claim 1, characterized in that, In step S4, the weighted fused feature values ​​are input into the multilayer perceptron network, and the predicted pollutant concentration values ​​are output, including: The weighted fused feature values ​​are input into the input layer of the multilayer perceptron network; the feature vectors of the weighted fused feature values ​​are then passed to each neuron node in the first hidden layer through the input layer. In the first hidden layer, each neuron node performs weighted summation and non-linear activation on the feature vector to obtain the first hidden feature layer; the first hidden feature layer is then passed to the second hidden layer, and the weighted mapping and non-linear activation operations are repeated to obtain the second hidden feature layer. The second hidden feature is passed to the third hidden layer, and the output is input to the output layer to determine the predicted pollutant concentration.