Carbon emission anomaly data detection method based on spatio-temporal graph neural network

Abstract

The present application relates to the technical field of carbon emission detection, in particular to a carbon emission anomaly data detection method based on a space-time graph neural network. The method comprises collecting data through a sensor network and a micro-meteorological station array to obtain a network input carbon emission tensor and an adjacency matrix, generating a hidden layer state tensor through information aggregation operation, and generating a prediction output tensor through time evolution operation; a composite objective function is constructed based on the network input carbon emission tensor and the prediction output tensor; the composite objective function contains a reconstruction error term of the Euclidean distance between the network input carbon emission tensor and the prediction output tensor, and a physical constraint penalty term converted from the gas mass conservation equation; the gradient vector of the composite objective function with respect to the network input carbon emission tensor is calculated, and a back propagation update operation is performed to output an optimized input tensor. The present application realizes seamless connection from environmental perception to physical constraint to data correction, and improves the overall repair capability for carbon emission anomaly data.

Description

Technical Field

[0001] This invention relates to the field of carbon emission detection technology, specifically a method for detecting abnormal carbon emission data based on spatiotemporal graph neural networks. Background Technology

[0002] In complex environments, continuous monitoring of greenhouse gases (such as carbon dioxide and methane) with high spatiotemporal resolution is fundamental for achieving quantitative assessment of carbon emissions and environmental compliance management. Current technologies generally employ an architecture combining distributed sensor networks with a central monitoring system. This involves real-time concentration data collection via monitoring terminals deployed in different geographical locations, followed by screening of uploaded data using preset statistical thresholds. When carbon emission monitoring values ​​exceed safety thresholds, alarms are triggered and data anomalies are flagged. This is currently the mainstream method in the industry for controlling the quality of environmental monitoring data.

[0003] However, in practical applications, mainstream technical solutions face problems such as rigid resource allocation and passive response, specifically in the following aspects: existing anomaly judgment logic is mostly based on the geographical coordinates of monitoring points to build a fixed distance model, assuming that the data of nearby nodes should be similar; and data sampling and transmission are usually carried out at a fixed frequency, lacking the ability to adaptively adjust resource input according to data fluctuations; in addition, there is the problem of poor physical consistency, ignoring the nonlinear modulation effect of sensor data on gas transmission paths, resulting in the inability to establish a connection relationship that conforms to the diffusion law, causing the monitoring system to suffer from severe distortion in the extraction of spatiotemporal correlation features when facing asymmetric diffusion dominated by wind direction, making it difficult to distinguish abnormal numerical anomalies.

[0004] To address this, a method for detecting abnormal carbon emission data based on spatiotemporal graph neural networks is proposed. Summary of the Invention

[0005] The purpose of this invention is to provide a method for detecting abnormal carbon emission data based on spatiotemporal graph neural networks, so as to solve the problems mentioned in the background art.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a method for detecting anomaly data in carbon emissions based on spatiotemporal graph neural networks, comprising: Greenhouse gas concentration and spatial elevation topological coordinates of a three-dimensional sensor network within an industrial park, as well as near-surface turbulent field data from a micro-weather station array, are collected. Through feature high-dimensional mapping, a network input carbon emission tensor is generated. Based on the near-surface turbulent field data, a nonlinear diffusion trajectory is calculated, and an adjacency matrix is ​​constructed by combining the spatial elevation topological coordinates. The carbon emission tensor and adjacency matrix of the network input are loaded into the spatiotemporal graph neural network model. The graph convolution operator is used to perform information aggregation on the carbon emission tensor of the network input to generate the hidden state tensor. Then, the gated recurrent unit is used to perform temporal evolution on the hidden state tensor to generate the predicted output tensor. A composite objective function is constructed based on the network input carbon emission tensor and the predicted output tensor. The composite objective function includes a reconstruction error term for the Euclidean distance between the network input carbon emission tensor and the predicted output tensor, and transforms the gas mass conservation equation into a physical constraint penalty term for physical constraints. With the weight parameters of the fixed spatiotemporal graph neural network model, the gradient vector of the composite objective function with respect to the network input carbon emission tensor is calculated; based on the gradient vector, backpropagation is performed to update the network input carbon emission tensor, and the optimized input tensor is output.

[0007] Preferably, the specific generation process of the network input carbon emission tensor and adjacency matrix includes: constructing a three-dimensional sensor network at ground monitoring stations, high-altitude towers, and the top of equipment in the industrial park to obtain greenhouse gas concentrations and simultaneously analyze the spatial elevation topological coordinates of each sensor node; deploying a micro-weather station array around the wind duct and emission source to obtain near-surface turbulent field data including horizontal wind speed vectors, vertical wind shear index, and atmospheric turbulence intensity; extracting the time-domain trend component and frequency-domain fluctuation component of the greenhouse gas concentration; mapping the near-surface turbulent field data into an environmental state embedding vector using an environmental encoder network; and... Tensor concatenation operations are performed on the temporal trend component, frequency fluctuation component, and environmental state embedding vector to generate the network input carbon emission tensor. Based on the horizontal wind speed vector and atmospheric turbulence intensity in the near-surface turbulence field data, the Gaussian plume model is driven to solve the nonlinear diffusion trajectory of greenhouse gases. Combining spatial elevation topological coordinates, the projected distance and lateral offset distance of the target node relative to the source node on the nonlinear diffusion trajectory are calculated. The coverage probability of the target node is calculated based on the projected distance and lateral offset distance, and the vertical connection weight is corrected by introducing a gravity settling factor to generate an adjacency matrix.

[0008] Preferably, the specific operation process of the graph convolution operator includes: reading a preset node activity state vector, applying the node activity state vector as a mask to the adjacency matrix to generate a carbon emission topology matrix; accumulating the element values ​​of the carbon emission topology matrix row by row to generate a carbon emission degree matrix, performing a normalization operation on the adjacency matrix based on the carbon emission degree matrix to generate a normalized Laplacian matrix; initializing and configuring a trainable weight matrix and a bias vector; performing a matrix left multiplication operation on the network input carbon emission tensor using the normalized Laplacian matrix to generate a carbon emission spatial aggregation feature tensor; performing a linear mapping operation on the carbon emission spatial aggregation feature tensor using the trainable weight matrix and superimposing the bias vector to generate a linear transformation feature tensor; inputting the linear transformation feature tensor into a nonlinear activation function layer to perform an element-wise activation operation, and outputting a hidden layer state tensor.

[0009] Preferably, the specific computation process of the gated recurrent unit includes: introducing the carbon emission memory state tensor from the previous time step, performing concatenation and linear transformation with the hidden layer state tensor, and generating a reset gate tensor and an update gate tensor via a Sigmoid activation function; performing a Hadamard product operation on the carbon emission memory state tensor from the previous time step and the reset gate tensor to generate a retained memory tensor; fusing the retained memory tensor with the hidden layer state tensor and inputting it into a hyperbolic tangent activation function layer to generate a candidate state tensor; performing a linear weighted combination operation on the carbon emission memory state tensor from the previous time step and the candidate state tensor according to the weight ratio defined by the update gate tensor, outputting the carbon emission memory state tensor from the current time step; inputting the carbon emission memory state tensor from the current time step into a fully connected distribution mapping layer to perform a linear regression operation, outputting a carbon emission prediction mean tensor and a confidence interval tensor in parallel, and combining the carbon emission prediction mean tensor and the confidence interval tensor into a prediction output tensor.

[0010] Preferably, the specific generation process of the reconstruction error term includes: using the network input carbon emission tensor as the minuend and the predicted mean tensor in the predicted output tensor as the subtrahend, performing element-wise subtraction to generate the original carbon emission residual tensor; extracting the numerical boundary of the confidence interval tensor to calculate the interval span, and performing a reciprocal operation on the interval span to generate an inverse proportional weight tensor; performing a square operation on the original carbon emission residual tensor, and performing a Hadamard product operation on the result and the inverse proportional weight tensor to generate a confidence error tensor; and performing a full-dimensional cumulative operation on the confidence error tensor to generate the reconstruction error term.

[0011] Preferably, the specific generation process of the physical constraint penalty term includes: based on the node connection relationship defined by the adjacency matrix, extracting the in-degree neighbor set flowing into the target node and the out-degree neighbor set flowing out of the target node; using the carbon emission prediction mean tensor and combining it with the near-surface turbulent field data, calculating the carbon flux values ​​of the target node and the in-degree neighbor set and the out-degree neighbor set, and constructing a flux balance equation; calculating the flux divergence residual between the sum of inflow flux and the sum of outflow flux in the flux balance equation; defining the absolute value of the flux divergence residual as the physical violation amount, and multiplying the physical violation amount by a preset Lagrange multiplier coefficient used as the penalty weight coefficient to generate the physical constraint penalty term.

[0012] Preferably, the specific process of the backpropagation update operation includes: initializing the iteration counter and step size hyperparameter; calculating the partial derivative tensor of the composite objective function with respect to the network input carbon emission tensor, performing a weighted operation on the partial derivative tensor using the step size hyperparameter to generate an updated gradient tensor; performing a tensor subtraction operation between the network input carbon emission tensor and the updated gradient tensor to generate an intermediate update tensor; constructing a hard threshold truncation function defining the effective numerical domain, performing a numerical projection operation on the intermediate update tensor using the hard threshold truncation function to forcibly map elements exceeding the effective numerical domain to boundary values, generating an updated network input carbon emission tensor; calculating the numerical variation of the composite objective function, and generating a carbon emission iteration loop trigger signal when the numerical variation is greater than a preset convergence threshold; responding to the carbon emission iteration loop trigger signal, accumulating the iteration counter and calculating the partial derivative tensor until the numerical variation is no greater than the convergence threshold, and outputting the optimized input tensor.

[0013] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. Carbon emission data of the target area is collected through a three-dimensional sensor network and a micro-weather station array. A high-dimensional feature mapping is used to obtain the network input carbon emission tensor. An adjacency matrix defining node connectivity is constructed based on near-surface turbulent field data, effectively integrating the environmental wind field's influence on gas transport into the graph structure. This allows the network input carbon emission tensor and adjacency matrix to be loaded into a spatiotemporal graph neural network model, enabling the generation of hidden state tensors through graph convolution operators by performing neighborhood information aggregation operations on the network input carbon emission tensor. Since the adjacency matrix directly reflects the dynamic connectivity defined by the wind field vector, the graph convolution operator no longer relies solely on static geographical distance when aggregating neighborhood information, but fully considers the impact of wind speed and direction on carbon emission diffusion, improving the accuracy of spatiotemporal feature extraction and adaptability to dynamic environments.

[0014] 2. A composite objective function is constructed based on the network input carbon emission tensor and the predicted output tensor. This composite objective function includes a reconstruction error term for the Euclidean distance between the network input carbon emission tensor and the predicted output tensor, and the gas mass conservation equation is transformed into a physical constraint penalty term, which helps the prediction results conform to objective physical laws. During the process of generating the predicted output tensor by performing time-series evolution operations on the hidden state tensor using gated recurrent units, the physical constraint penalty term forces the model output to follow the gas mass conservation law, avoiding prediction biases that might arise from purely data-driven models that violate physical laws. By introducing the physical constraint penalty term as a regularization factor into the optimization objective, while minimizing the Euclidean distance reconstruction error, the flux balance constraint of gas diffusion must be satisfied. This ensures that the generated predicted output tensor not only numerically approximates the actual observation value but also eliminates prediction biases that might arise from purely data-driven models that violate physical laws.

[0015] 3. By fixing the weight parameters of the spatiotemporal graph neural network model, the gradient vector of the composite objective function with respect to the network input carbon emission tensor is calculated. Then, backpropagation is performed on the network input carbon emission tensor based on the gradient vector until the composite objective function converges numerically. The updated optimized input tensor is then output, achieving high-precision reconstruction and repair of anomalous data. Through backpropagation updates based on the gradient vector, anomalous values ​​in the network input carbon emission tensor can be adaptively adjusted, gradually approximating a normal distribution state that conforms to spatiotemporal evolution and physical conservation constraints. When the composite objective function converges numerically, the updated optimized input tensor output is the normal value verified by physical and spatiotemporal consistency checks. This effectively distinguishes between normal nonlinear fluctuations caused by atmospheric turbulence and numerical abrupt changes caused by sensor malfunctions, avoiding false anomaly detections caused by drastic changes in the environmental flow field.

[0016] 4. By calculating the nonlinear diffusion trajectory based on near-surface turbulent field data and constructing an adjacency matrix using spatial elevation topology, a dynamic spatial topology conforming to the three-dimensional Gaussian plume diffusion law of industrial parks was first established. Based on this, a composite objective function was constructed using the network input carbon emission tensor and the predicted output tensor, and the gas mass conservation equation was transformed into a physical constraint penalty term, thus superimposing a carbon flux balance rule on the three-dimensional topology. This design allows the model to calculate the gradient vector of the composite objective function with respect to the network input carbon emission tensor when the weight parameters of the spatiotemporal graph neural network model are fixed. This gradient vector not only contains numerical reconstruction error information but also incorporates the physical violation signal propagated backward along the nonlinear diffusion trajectory. Furthermore, based on the gradient vector, a backpropagation update operation is performed on the network input carbon emission tensor, driving the data to self-repair within the multi-dimensional constraint space until the composite objective function converges numerically, outputting the updated optimized input tensor. This mechanism achieves a seamless connection from three-dimensional environmental perception to physical flux constraints to adaptive data correction, effectively improving the overall repair capability for abnormal carbon emission data under complex flow fields in industrial parks. Attached Figure Description

[0017] Figure 1 The flowchart shows the carbon emission anomaly data detection method based on spatiotemporal graph neural network proposed in the embodiments of this invention application; Figure 2 This is a flowchart illustrating the multi-source data tensor quantization and dynamic wind field map construction process proposed in an embodiment of this invention. Figure 3 This is a flowchart of the physical constraint-driven iterative repair process for abnormal data proposed in an embodiment of this 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] Please see Figures 1-3 The carbon emission anomaly data detection method based on spatiotemporal graph neural network provided by this invention has the following specific steps: Greenhouse gas concentration and spatial elevation topological coordinates of a three-dimensional sensor network within an industrial park, as well as near-surface turbulent field data from a micro-weather station array, are collected. Through feature high-dimensional mapping, a network input carbon emission tensor is generated. Based on the near-surface turbulent field data, a nonlinear diffusion trajectory is calculated, and an adjacency matrix is ​​constructed by combining the spatial elevation topological coordinates. The carbon emission tensor and adjacency matrix of the network input are loaded into the spatiotemporal graph neural network model. The graph convolution operator is used to perform information aggregation on the carbon emission tensor of the network input to generate the hidden state tensor. Then, the gated recurrent unit is used to perform temporal evolution on the hidden state tensor to generate the predicted output tensor. A composite objective function is constructed based on the network input carbon emission tensor and the predicted output tensor. The composite objective function includes a reconstruction error term for the Euclidean distance between the network input carbon emission tensor and the predicted output tensor, and transforms the gas mass conservation equation into a physical constraint penalty term for physical constraints. With the weight parameters of the fixed spatiotemporal graph neural network model, the gradient vector of the composite objective function with respect to the network input carbon emission tensor is calculated; based on the gradient vector, backpropagation is performed to update the network input carbon emission tensor, and the optimized input tensor is output.

[0020] The technical solution of the present invention will be further described in detail below with reference to specific embodiments.

[0021] Example 1 This application discloses a method for detecting abnormal carbon emission data based on a spatiotemporal graph neural network. See also... Figure 1 The specific steps proposed in this invention include: S1, collecting greenhouse gas concentration and spatial elevation topological coordinates from a three-dimensional sensor network within an industrial park, and near-surface turbulent field data from a micro-weather station array; generating a network input carbon emission tensor through feature high-dimensional mapping; calculating nonlinear diffusion trajectories based on near-surface turbulent field data, and constructing an adjacency matrix by combining spatial elevation topological coordinates; S2, loading the network input carbon emission tensor and adjacency matrix into a spatiotemporal graph neural network model, performing information aggregation on the network input carbon emission tensor through graph convolution operators to generate a hidden state tensor, and utilizing gated loops... Unit S3 performs temporal evolution on the hidden state tensor to generate the predicted output tensor; S4 constructs a composite objective function based on the network input carbon emission tensor and the predicted output tensor; the composite objective function includes a reconstruction error term of the Euclidean distance between the network input carbon emission tensor and the predicted output tensor, and transforms the gas mass conservation equation into a physical constraint penalty term of physical constraints; S5 fixes the weight parameters of the spatiotemporal graph neural network model and calculates the gradient vector of the composite objective function with respect to the network input carbon emission tensor; backpropagation is performed on the network input carbon emission tensor based on the gradient vector to output the optimized input tensor.

[0022] Further, greenhouse gas concentration and spatial elevation topological coordinates from a three-dimensional sensor network within the industrial park, as well as near-surface turbulent field data from a micro-weather station array, are collected. Through feature high-dimensional mapping, a network input carbon emission tensor is generated. Based on the near-surface turbulent field data, a nonlinear diffusion trajectory is calculated, and combined with the spatial elevation topological coordinates, an adjacency matrix is ​​constructed; this corresponds to step S1 above; see [reference needed]. Figure 2 The specific implementation process includes: A three-dimensional sensor network was constructed at the gas ground monitoring station, high-altitude tower, and top of the equipment in the industrial park to collect greenhouse gas concentrations and simultaneously analyze the spatial elevation topological coordinates of each sensor node. Micro-weather station arrays were deployed around wind ducts and emission sources to collect near-surface turbulent field data including horizontal wind speed vectors, vertical wind shear indices, and atmospheric turbulence intensity. Noise reduction was performed on the greenhouse gas concentration data, extracting the temporal trend component and frequency fluctuation component. An environmental encoder was used to map the near-surface turbulent field data into an environmental state embedding vector. Tensor fusion operations were performed on the temporal trend component, frequency fluctuation component, and environmental state embedding vector to generate the network input carbon emission tensor. Based on the horizontal wind speed vector and atmospheric turbulence intensity in the near-surface turbulent field data, a Gaussian plume model was driven to solve the nonlinear diffusion trajectory of greenhouse gases. Combining the spatial elevation topological coordinates, the projection distance and lateral offset distance of the target node relative to the source node on the nonlinear diffusion trajectory were calculated. The coverage probability of the target node was calculated based on the projection distance and lateral offset distance, and a gravity settling factor was introduced to correct the vertical connection weights, generating an adjacency matrix.

[0023] Specifically, the process of generating the network input carbon emission tensor is as follows: In practical industrial park monitoring scenarios, such as an industrial park covering an area of ​​25 square kilometers with 50 monitoring nodes, a fully covered three-dimensional sensor network is first constructed. This network not only includes ground-based gas monitoring stations but also utilizes existing infrastructure to deploy sensor nodes on high-altitude towers (e.g., 30 to 50 meters high) and on top of key emission devices (e.g., 10 to 20 meters high). It collects minute-level greenhouse gas concentrations at each node and uses differential positioning technology to simultaneously resolve the spatial elevation topological coordinates of each sensor node. Simultaneously, micro-weather station arrays are deployed around wind ducts and emission sources to collect near-surface turbulence field data at high frequencies. To address the multi-source heterogeneity of the collected data, signal preprocessing based on empirical mode decomposition is first performed. The original greenhouse gas concentration time series is used as input and decomposed into several intrinsic mode function components and one residual component. Based on the frequency characteristics, the low-frequency, high-energy intrinsic mode functions and the residual component are superimposed to reconstruct a time-domain trend component, which is used to characterize the long-term emission baseline level. The high-frequency intrinsic mode functions are superimposed to reconstruct a frequency-domain fluctuation component, which is used to characterize the instantaneous concentration fluctuations affected by atmospheric turbulence.

[0024] In this embodiment, after generating the frequency domain fluctuation components, a signal discrimination mechanism based on spectral kurtosis is introduced. The spectral kurtosis index of the fluctuation components is calculated, and based on their non-Gaussian characteristics, normal high-frequency fluctuations caused by gas turbulence are distinguished from impulse noise caused by sensor circuit failures, retaining only fluctuation components that conform to the statistical laws of turbulence. Specifically, the extracted frequency domain fluctuation components undergo spectral kurtosis calculation based on fourth-order cumulants. Since atmospheric turbulence approximately follows a Gaussian distribution at the microscopic scale, while impulse noise generated by sensor circuit aging or poor contact exhibits strong non-Gaussianity, the sliding time window length is set to 64 sampling points, with a step size of 16. For each fluctuation sequence within the window, its spectral kurtosis value is calculated. This value is calculated by dividing the fourth-order central moment of the data within the window by the square of the second-order central moment, and then subtracting 3 from the quotient. The kurtosis threshold is set to 1.5 (this threshold is set based on empirical values ​​for the electromagnetic environment noise floor of industrial parks). If the calculated spectral kurtosis value is greater than the kurtosis threshold, it is determined that the fluctuation within this time window is mainly dominated by non-natural impulse noise. Soft thresholding denoising based on wavelet packets is then performed, shrinking the wavelet coefficients to 10% of their original values ​​to suppress noise. If the spectral kurtosis value is less than or equal to the kurtosis threshold, it is determined to be a valid atmospheric turbulence fluctuation, and the original value is retained. Finally, the processed sequence is reassembled and used as the corrected frequency domain fluctuation component in subsequent tensor fusion. This step effectively eliminates interference from non-Gaussian impulse noise from the sensor, which helps the fluctuation characteristics of the input model to truly reflect the intensity of atmospheric turbulence.

[0025] Subsequently, an environment encoder based on a multilayer perceptron was constructed and trained to map low-dimensional near-surface turbulence field data into a high-dimensional environmental state embedding vector that matches the carbon emission feature dimension. The specific model architecture and parameter settings of this environment encoder are as follows: The input layer receives a 5-dimensional near-surface turbulence field feature vector (including horizontal wind speed, horizontal wind direction sine, horizontal wind direction cosine, vertical wind shear index, and atmospheric turbulence intensity); this is followed by two fully connected hidden layers, the first with 64 neurons and the second with 32 neurons. During computation, each fully connected layer performs a linear weighted summation on the input data and superimposes a bias parameter vector. Then, a nonlinear transformation is performed by modifying the linear unit activation function, truncating negative values ​​to zero to increase sparsity; the output layer outputs a 16-dimensional environmental state embedding vector. Finally, tensor fusion operations on the multidimensional features are performed. The one-dimensional time-domain trend component and one-dimensional frequency-domain fluctuation component extracted from each node are concatenated with the 16-dimensional environmental state embedding vector output by the environmental encoder, performing a concatenation operation on the feature channel dimension to generate a composite feature vector with a feature dimension of 18. For example, with a batch size of 64, a time step of 12 (corresponding to data from the past hour), and a number of nodes of 50, the final generated network input carbon emission tensor has a dimension of "64 x 12 x 50 x 18". This tensor numerically preserves the complex three-dimensional spatiotemporal flow field information and gas evolution characteristics within the monitoring area, providing a high-quality input rich in physical semantics for the subsequent spatiotemporal graph neural network.

[0026] Specifically, the process of constructing the adjacency matrix is ​​as follows: First, the spatial structure of the three-dimensional sensor network is parameterized. The longitude, latitude, and altitude values ​​of each sensor node are converted into three-dimensional coordinate vectors in a local Cartesian coordinate system using Gaussian Srkrür projection, where the vertical coordinate represents the vertical height relative to the ground (e.g., setting the vertical height of the gas ground monitoring station to 2 meters and the vertical height of the tower node to 30 meters). To construct an adjacency matrix conforming to atmospheric dynamics, diffusion parameter calculations based on a Gaussian plume model are performed, with the input being near-surface turbulence field data collected by micro-weather stations, including horizontal wind speed scalar, wind direction angle, and atmospheric turbulence intensity. The specific algorithm flow is as follows: First, using the Pasquale-Gifford classification method, atmospheric stability is divided into six levels from A to F based on horizontal wind speed scalar, daytime solar radiation level, and nighttime cloud cover (where A is extremely unstable, B is unstable, C is weakly unstable, D is neutral, E is weakly stable, and F is stable). For example, when strong solar radiation is detected and the horizontal wind speed scalar is less than 2.5 meters per second, it is classified as level A. Second, the preset Briggs diffusion coefficient formula is used to calculate the lateral and vertical diffusion coefficients at a distance downwind. The lateral diffusion coefficient is calculated by multiplying the downwind distance by the regression coefficient 'a', then by a correction term, which is 1 plus 0.0001 times the sum of downwind distances raised to the power of -0.5. The vertical diffusion coefficient is calculated by multiplying the downwind distance regression coefficient 'd' raised to the power of the regression coefficient 'c'. The proportionality coefficient and exponential parameters are read from a pre-defined lookup table based on the determined stability level. Typical lookup table data is as follows: when atmospheric stability is level A, coefficient 'a' is 0.22, coefficient 'c' is 0.20, and coefficient 'd' is 1.0. For Grade B, coefficient a is 0.16, coefficient c is 0.12, and coefficient d is 1.0; for Grade C, coefficient a is 0.11, coefficient c is 0.08, and coefficient d is 1.0; for Grade D, coefficient a is 0.08, coefficient c is 0.06, and coefficient d is 1.0; for Grade E, coefficient a is 0.06, coefficient c is 0.03, and coefficient d is 1.0; for Grade F, coefficient a is 0.04, coefficient c is 0.016, and coefficient d is 1.0.

[0027] Next, combining the spatial coordinates of each node with the calculated diffusion parameters, the coverage probability between nodes is calculated. For any source node and target node in the network, a relative position vector pointing from the source node to the target node is first constructed. This relative position vector is then used to construct a two-dimensional rotation matrix using the current horizontal wind direction angle. The relative position vector is left-multiplied by the rotation matrix, projecting it into a local wind field coordinate system with the horizontal wind direction as the X-axis (i.e., downwind) and the direction perpendicular to the wind as the Y-axis (i.e., crosswind). In this coordinate system, three key geometric features are generated: the X-axis component is the projection distance (i.e., downwind distance); the absolute value of the Y-axis component is the horizontal offset distance; and the absolute value of the altitude difference between the target node and the source node is directly calculated as the vertical offset distance. If the projection distance is less than or equal to 0, it indicates that the target node is upwind, and the connection weight is reset to 0. If the projection distance is greater than 0, the above geometric features are substituted into the Gaussian concentration distribution equation for normalization to generate the coverage probability of the source node over the target node. The specific calculation logic is as follows: The horizontal and vertical components are calculated separately. The vertical component is calculated by multiplying the natural constant by the negative square of the vertical offset distance by twice the square of the vertical diffusion coefficient. The final coverage probability is the product of the horizontal and vertical components. Example data: Assume the source node is located at the top of the tower (coordinate height 30 meters), the target node is located on the ground downwind (coordinate height 2 meters), the current horizontal wind speed is 4 meters per second, and the atmospheric stability is determined to be neutral. First, the projected distance is calculated to be 200 meters. Looking up the table, the vertical diffusion coefficient at this distance is approximately 8 meters. The vertical offset distance is calculated to be 28 meters. Substituting these values ​​into the formula, the vertical contribution component is calculated by multiplying the natural constant by the negative square of 28 by twice the square of 8. The result is approximately 0.002. This value quantifies that even downwind, due to the large height difference and insufficient vertical turbulent mixing, the actual physical transmission impact of the source node on the target node is negligible, thus assigning a very low connection weight.

[0028] To further improve the model's ability to fit the vertical transport characteristics of macromolecular greenhouse gases and associated particulate matter, an asymmetric correction operation based on the gravity settling factor is performed. The inputs include the elevation difference between the source and target nodes, the horizontal wind speed scalar, and a preset settling rate parameter. The settling rate parameter typically ranges from 0.1 to 1.0, depending on the molecular weight and aerosol adhesion characteristics of the monitored gases (e.g., 0.3 to 0.6 for carbon dioxide; in this embodiment, it is set to 0.5). The specific correction logic is as follows: First, calculate the difference between the elevation of the target node and the elevation of the source node. If this difference is less than 0, it indicates that the target node is below the source node, and the greenhouse gases exhibit a natural settling tendency under gravity, requiring a weight enhancement operation. The enhancement coefficient is calculated by multiplying the settling rate parameter by the absolute value of the elevation difference, and then dividing by the product of the horizontal wind speed scalar and the projected distance (to prevent the denominator from being zero, a preset minimum stability constant, such as 10, is added to the product term). -6 The quotient is then incremented by 1 to generate the settlement enhancement coefficient. If the difference is greater than 0, it indicates that the target node is located above and requires work to overcome gravity to reach it. In this case, a weight attenuation operation is performed, with the attenuation coefficient set to the reciprocal of the settlement enhancement coefficient. Finally, the coverage probability is multiplied by either the settlement enhancement coefficient or the attenuation coefficient to generate the corrected connection weight. For example, if the source node is 30 meters above the ground and the target node is 2 meters below the ground, the altitude difference is -28 meters, placing it in a settlement zone. The settlement rate parameter is set to 0.5, the horizontal wind speed scalar to 4 meters per second, and the projection distance to 200 meters. First, the settlement term is calculated: 0.5 multiplied by 28 equals 14, 4 multiplied by 200 equals 800, and 14 divided by 800 equals 0.0175. Adding 1 gives the settlement enhancement coefficient of 1.0175. If the coverage probability is 0.002, the final corrected weight is 0.002 multiplied by 1.0175, approximately equal to 0.002035.

[0029] After calculating and correcting the weights of all node pairs, the final adjacency matrix assembly and normalization operation is performed. First, a square matrix with dimensions equal to the product of the total number of sensor nodes and the total number of sensor nodes is constructed. The coverage probability values ​​calculated in the previous steps and corrected for settlement are filled into the corresponding positions of this square matrix. For elements on the diagonal of the square matrix, a preset self-loop weight value (usually set to 1) is assigned to ensure that nodes retain their own feature information when aggregating neighborhood information. To eliminate computational noise caused by extremely small weights and improve computational efficiency, a sparsity threshold is set, forcing elements in the square matrix with values ​​less than this threshold to zero. Subsequently, to prevent feature value explosion or gradient vanishing problems during subsequent multi-layer graph neural network operations, row normalization must be performed on this sparse square matrix. The specific normalization operation logic is as follows: traverse the square matrix row by row, calculate the sum of the values ​​of all non-zero elements in each row, i.e., the total out-degree weight of that node; then, divide the value of each element in that row by the sum of the values ​​in that row to generate the normalized connection weights. The matrix output after this series of operations is the final adjacency matrix. For example, assuming the source node is tower node number one, after calculating the weights of all downwind targets, it is found that the sum of the corrected weights of all its neighboring nodes in the first row of the matrix (including itself and the ground target node with a corrected weight of 0.002035) is 0.5. Then, a normalization operation is performed: 0.002035 is divided by 0.5, resulting in 0.00407, which is the element value of the corresponding column in the first row of the final adjacency matrix. This represents that in one graph convolution aggregation operation, the feature information of node number one will be transferred to the ground target node at a ratio of 0.407%.

[0030] A three-dimensional sensor network was constructed, and a dynamic weighting mechanism based on a Gaussian plume model was introduced. The vertical connection weights were corrected by solving the nonlinear diffusion trajectory and introducing a gravity settling factor, enabling the generated adjacency matrix to characterize the three-dimensional transport of greenhouse gases under the combined influence of non-uniform turbulent wind fields and their own gravity. This topology provides a high-quality prior knowledge graph for subsequent spatiotemporal graph neural networks to accurately extract anomalous spatiotemporal patterns, improving the detection accuracy of carbon emission anomalies under complex and variable meteorological conditions.

[0031] Furthermore, the input carbon emission tensor and adjacency matrix are loaded into the spatiotemporal graph neural network model, and information aggregation is performed on the input carbon emission tensor using a graph convolution operator to generate the hidden state tensor; this corresponds to step S2 above; the specific implementation process includes: Read the preset node activity state vector, apply the node activity state vector as a mask to the adjacency matrix to generate a carbon emission topology matrix; accumulate the element values ​​of the carbon emission topology matrix row by row to generate a carbon emission degree matrix; perform a normalization operation on the adjacency matrix based on the carbon emission degree matrix to generate a normalized Laplacian matrix; initialize and configure a trainable weight matrix and a bias vector; perform a matrix left multiplication operation on the network input carbon emission tensor using the normalized Laplacian matrix to generate a carbon emission spatial aggregation feature tensor; perform a linear mapping operation on the carbon emission spatial aggregation feature tensor using the trainable weight matrix and superimpose the bias vector to generate a linear transformation feature tensor; input the linear transformation feature tensor into a nonlinear activation function layer to perform element-wise activation operation, and output a hidden layer state tensor.

[0032] Specifically, the process of generating the normalized Laplacian matrix is ​​as follows: In actual long-term series monitoring, some sensor nodes may be temporarily inactive due to equipment failure, network outages, or maintenance / calibration. Directly incorporating these nodes into graph convolution operations would result in invalid zero values ​​or noise data propagating to neighboring nodes, contaminating the overall feature distribution. Therefore, a preset node activity state vector is first read. This vector has the same dimension as the total number of sensor nodes (e.g., 50 dimensions), with active nodes marked as 1 and faulty nodes marked as 0. This vector is then used as a mask on the adjacency matrix using a broadcast mechanism, performing element-wise multiplication to force all connection weights associated with inactive nodes to zero, thereby generating a cleaned carbon emission topology matrix. Subsequently, a row-wise summation operation is performed on this carbon emission topology matrix to calculate the sum of effective connection strengths for each node at the current time, generating a diagonal carbon emission degree matrix. The carbon emission topology matrix is ​​then normalized. The specific calculation process is as follows: Calculate the arithmetic square root of the reciprocal of the diagonal elements in the carbon emission index matrix to construct a negative half-power matrix of the carbon emission index matrix; multiply this power matrix by the carbon emission topology matrix on the left and right respectively, perform symmetric normalization operations, and finally generate a normalized Laplace matrix. The spectral radius of this matrix is ​​restricted to within 1.

[0033] In this embodiment, a soft-threshold smoothing process based on the signal-to-noise ratio is performed on the node activity state vector. The ratio of the signal fluctuation variance of each sensor's local time-series slice to the historical baseline variance is calculated, and a continuous reliability score is generated by combining this with the data missing rate. First, the proportion of non-empty data for each node within the current time window is statistically analyzed. Second, the variance of the data within this window is calculated, and the historical average variance of the node over the past 24 hours is retrieved to calculate the volatility ratio. This ratio equals the current variance divided by the sum of the historical average variance and the minimum stable value. Then, a soft-threshold state vector is generated using a smoothing function, where the reliability score of any node equals the current proportion of non-empty data multiplied by the exponent of the natural constant. This exponent is a negative attenuation coefficient multiplied by the corrected volatility deviation, where the corrected volatility deviation is the larger of zero and the volatility ratio minus the upper limit of volatility tolerance. In this embodiment, the upper limit of volatility tolerance is set to 3.0, and the attenuation coefficient is set to 0.5. If a node has complete data but the noise is four times the historical average, its score will drop from 1.0 to approximately 0.6. This continuous numerical vector is used as a soft mask and broadcast multiplication is performed with the corresponding row of the adjacency matrix. This allows graph convolution to automatically reduce the weight contribution of noisy or unstable nodes when aggregating neighborhood information, preserving potential trend information, effectively suppressing noise propagation from unhealthy nodes, and avoiding excessive information loss due to misjudgment of faults.

[0034] Specifically, the process for setting model parameters is as follows: To map low-dimensional raw monitoring data to a high-dimensional feature space to capture complex nonlinear patterns, linear transformation parameters need to be defined in the graph convolutional layer. The raw physical data collected by the sensor network may only have 5 to 7 dimensions (e.g., including carbon dioxide concentration, meridional wind speed, zonal wind speed, temperature, air pressure, etc.), but the feature dimension of the carbon emission tensor input to the spatiotemporal graph neural network has been expanded to N dimensions (in this embodiment, N is 18, i.e., including a 1-dimensional trend component, a 1-dimensional fluctuation component, and a 16-dimensional environmental state embedding vector). The feature dimension of the hidden layer state tensor is set to 64. Therefore, the trainable weight matrix is ​​initialized as a two-dimensional tensor with a dimension of 18 multiplied by 64, typically using a Kemming uniform distribution initialization strategy to maintain the consistency of input and output variance; the bias vector is initialized as a one-dimensional zero vector with a dimension of 64. In terms of hyperparameter configuration for model training, an adaptive moment estimation optimizer is selected, with an initial learning rate of 0.001, and the exponential decay rate of the first-order moment estimation is configured to be 0.9, and the exponential decay rate of the second-order moment estimation is configured to be 0.999. To prevent overfitting, an L2 regularization term was also introduced during the training process, with its weight decay coefficient set to 10. -4 .

[0035] Specifically, the generation process of the hidden layer state tensor is as follows: The dimension of the network input carbon emission tensor is "batch size multiplied by time step multiplied by number of nodes multiplied by input feature dimension" (e.g., 64×12×50×5). During aggregation, the normalized Laplacian matrix (50×50 dimension) is treated as a spatial propagation operator and multiplied left by the spatial dimension of the network input carbon emission tensor. Physically, this is equivalent to weighted summation of the feature data of each sensor node's neighbors, with the weights determined by the normalized adjacency relationships. For example, for node A, its aggregated features not only include its current state but also incorporate the high-concentration plume features transmitted from upwind node B and the background concentration features of its spatially adjacent node C. Through this operation, the generated spatially aggregated carbon emission feature tensor (64×12×50×5 dimension) aggregates the environmental context information within the local receptive field of each spatial node. After spatial aggregation, a linear mapping is performed on the spatially aggregated feature tensor using an initialized trainable weight matrix, projecting the feature dimension from 5 to 64 dimensions and superimposing a bias vector to generate a linearly transformed feature tensor. To enable the model to handle nonlinear relationships (such as the complexity of turbulent diffusion), the tensor is input into a nonlinear activation function layer. This embodiment preferably uses ReLU as the activation function. The specific operational logic is as follows: for each element in the tensor, if its value is greater than 0, the value itself is output; if it is less than or equal to 0, 0 is output. The hidden state tensor output after activation processing not only preserves the spatiotemporal topology of the input data but also extracts high-order abstract semantic features.

[0036] By reading the preset node activity state vector and applying it as a mask to the adjacency matrix, a carbon emission topology matrix is ​​generated. This effectively shields the impact of faulty or offline nodes on neighborhood aggregation, enhancing the model's fault tolerance. Normalizing the adjacency matrix based on the carbon emission degree matrix generates a normalized Laplacian matrix, effectively preventing numerical explosion or gradient vanishing problems in multi-layer graph convolution operations. Using the normalized Laplacian matrix, a matrix left multiplication operation is performed on the network input carbon emission tensor to generate a carbon emission spatial aggregation feature tensor. This tensor, combined with the output hidden state tensor of the nonlinear activation function layer, improves the numerical stability and nonlinear expressive power of spatial feature extraction.

[0037] Furthermore, the input carbon emission tensor and adjacency matrix are loaded into the spatiotemporal graph neural network model, and a gated recurrent unit is used to perform temporal evolution on the hidden state tensor to generate the predicted output tensor; this corresponds to step S2 above; the specific implementation process includes: The carbon emission memory state tensor from the previous time step is introduced and concatenated with the hidden layer state tensor using a linear transformation. A reset gate tensor and an update gate tensor are then generated via a Sigmoid activation function. A Hadamard product operation is performed on the carbon emission memory state tensor from the previous time step and the reset gate tensor to generate a retained memory tensor. The retained memory tensor is fused with the hidden layer state tensor and input into a hyperbolic tangent activation function layer to generate a candidate state tensor. Based on the weight ratio defined in the update gate tensor, a linear weighted combination operation is performed on the carbon emission memory state tensor from the previous time step and the candidate state tensor to output the carbon emission memory state tensor from the current time step. The carbon emission memory state tensor from the current time step is input into a fully connected distribution mapping layer to perform a linear regression operation, outputting a carbon emission prediction mean tensor and a confidence interval tensor in parallel. Finally, the carbon emission prediction mean tensor and the confidence interval tensor are combined into a prediction output tensor.

[0038] Specifically, the generation process of the carbon emission prediction mean tensor and confidence interval tensor is as follows: First, a carbon emission memory state tensor from the previous time step is introduced. This tensor stores all historical accumulated information up to the previous time step, and its dimension is set to "batch size multiplied by number of nodes multiplied by hidden layer dimension" (e.g., 64 multiplied by 50 multiplied by 64). At the current time step, the hidden layer state tensor (i.e., the spatial aggregation feature at the current time step) is concatenated with the memory state tensor from the previous time step along the feature dimension to generate a combined input vector. Subsequently, a linear transformation is performed on this combined vector using two independent trainable weight matrices (typically 128 multiplied by 64), and bias vectors are superimposed on them respectively. The transformed result is input into a sigmoid activation function, which compresses the values ​​to the range of 0 to 1, generating a reset gate tensor that controls the degree of information forgetting and an update gate tensor that controls the fusion ratio of old and new information.

[0039] The Hadamard product operation (element-wise multiplication) is performed on the carbon emission memory state tensor from the previous time step and the reset gate tensor, retaining only historical information useful for the current prediction (such as persistent emission source features) to generate a retained memory tensor. Next, this retained memory tensor is concatenated or fused with the hidden state tensor from the current time step, linearly mapped using a third set of trainable weight matrices, and then input into a hyperbolic tangent activation function. The hyperbolic tangent function maps the output value to the interval -1 to 1, generating a candidate state tensor containing new features from the current time step.

[0040] Based on the weight ratio defined by the update gate tensor, a linear weighted combination operation is performed on the carbon emission memory state tensor of the previous time step and the candidate state tensor of the current time step. The specific calculation logic is: the current memory equals (1 minus the update gate value) multiplied by the memory of the previous time step plus the update gate value multiplied by the candidate state. To quantify the uncertainty of the prediction, the carbon emission memory state tensor of the current time step is input into a fully connected distributed mapping layer. This layer contains two sets of parallel linear regression weights: one set is used to perform a linear mapping on the carbon emission memory state tensor and to perform a nonlinear transformation on the linear output using the Softplus activation function to smooth and constrain the predicted value to the non-negative interval, generating a carbon emission prediction mean tensor that conforms to the physical concentration; the other set is used to fit the logarithmic variance of the prediction error, and after exponential transformation, outputs a confidence interval tensor (usually defined as ±1.96 times the standard deviation, covering 95% confidence level). The prediction mean tensor and the confidence interval tensor are combined into the prediction output tensor. During the model training phase, to simultaneously optimize mean accuracy and interval coverage, a loss function based on Gaussian negative log-likelihood was constructed. The specific calculation process is as follows: for each sample, the square of the difference between the true value and the predicted mean is divided by twice the prediction variance, and then half the natural logarithm of the prediction variance is added. Minimizing this loss function not only drives the predicted mean to approximate the actual monitoring data but also forces the model to expand the confidence interval in areas with large prediction errors (such as areas with turbulent wind fields).

[0041] By inputting the carbon emission memory state tensor at the current moment into a fully connected distributed mapping layer to perform linear regression, and outputting a carbon emission prediction mean tensor and a confidence interval tensor in parallel, a probabilistic prediction containing uncertainty information is provided. Combining the prediction mean tensor and the confidence interval tensor into a prediction output tensor can quantify the reliability of the prediction results, providing rich statistical evidence for subsequent anomaly detection.

[0042] Furthermore, a composite objective function is constructed based on the network input carbon emission tensor and the predicted output tensor. This composite objective function includes a reconstruction error term representing the Euclidean distance between the network input carbon emission tensor and the predicted output tensor, and transforms the gas mass conservation equation into a physical constraint penalty term, corresponding to step S3 above. The specific implementation process includes: Using the network input carbon emission tensor as the minuend and the predicted mean tensor in the predicted output tensor as the subtrahend, perform element-wise subtraction to generate the original carbon emission residual tensor; extract the numerical boundaries of the confidence interval tensor to calculate the interval span, and perform a reciprocal operation on the interval span to generate an inverse proportional weight tensor; perform a square operation on the original carbon emission residual tensor, and perform a Hadamard product operation on the result and the inverse proportional weight tensor to generate a confidence error tensor; perform a full-dimensional summation operation on the confidence error tensor to generate a reconstruction error term.

[0043] Specifically, the specific process for generating the reconstruction error term is as follows: To capture deviations between monitoring data and normal patterns, an element-wise subtraction operation is performed using the normalized network input carbon emission tensor (dimensions 64 x 12 x 50 x 5) as the minuend and the predicted mean tensor as the subtrahend. For example, if the normalized input value of a sensor node is 1.2 (corresponding to an actual concentration of 438 PPM), while the model predicts a mean of 1.0 (corresponding to 435 PPM) based on historical patterns, the generated original carbon emission residual tensor at that point will have a value of 0.2. This tensor directly reflects the degree to which the current monitoring value deviates from the normal diffusion state.

[0044] Due to the complexity of the atmospheric environment, the reliability of predictions varies significantly under different meteorological conditions. The confidence interval tensor accompanying the prediction output tensor is extracted, defining the upper and lower boundaries of the predicted values ​​(e.g., ±1.96 times the standard deviation). The difference between the upper and lower boundaries is calculated to obtain the interval span, and the reciprocal operation is performed on this span (to prevent division by zero errors, a very small smoothing constant, such as 10, is usually added). -6 This generates an inversely proportional weight tensor. In high-uncertainty regions with turbulent wind fields or missing data, the output interval span increases (e.g., from 5 PPM to 50 PPM), causing the corresponding weight values ​​to decrease sharply.

[0045] An element-wise squaring operation is performed on the original carbon emission residual tensor to quantify the magnitude energy of the bias. This squared result is then multiplied by the inverse proportional weight tensor using a Hadamard product (i.e., multiplying corresponding elements) to generate a confidence error tensor. The confidence error tensor is then summed across four dimensions: batch, time, node, and feature, to generate the final scalar reconstruction error term. This error term serves as the dominant gradient source for backpropagation, driving the neural network to adjust its weight parameters, enabling the model to approximate the actual observations as closely as possible under relatively deterministic, stable weather conditions.

[0046] By performing a Hadamard product operation on the squared result of the original carbon emission residual tensor and the inverse proportional weight tensor to generate a confidence error tensor, an adaptive weighted reconstruction error term is constructed. This gives greater weight to prediction results with narrower confidence intervals and reduces the error weight in areas with higher uncertainty. This effectively suppresses the interference of data fluctuations caused by environmental noise on anomaly detection and improves the sensitivity of the reconstruction error term to real anomalies.

[0047] Based on the node connection relationships defined by the adjacency matrix, the in-degree neighbor set flowing into the target node and the out-degree neighbor set flowing out of the target node are extracted. Using the carbon emission prediction mean tensor and combined with the near-surface turbulent field data, the carbon flux values ​​of the target node, the in-degree neighbor set, and the out-degree neighbor set are calculated to construct a flux balance equation. The flux divergence residual between the sum of inflow flux and the sum of outflow flux in the flux balance equation is calculated. The absolute value of the flux divergence residual is defined as the physical violation quantity, and the physical violation quantity is multiplied by a preset Lagrange multiplier coefficient to generate the physical constraint penalty term.

[0048] Specifically, the generation process of the physical constraint penalty term is as follows: In each forward propagation iteration of the spatiotemporal graph neural network, the adjacency matrix dynamically generated by the wind field is first analyzed. For any target sensor node in the network (e.g., the central monitoring station numbered 35), the adjacency matrix is ​​traversed column-wise to extract the indices of all source nodes with non-zero connection weights to the target node, forming an in-degree neighbor set representing the upwind pollutant transport sources; simultaneously, the matrix is ​​traversed row-wise to extract the indices of all nodes pointed to by the target node, forming an out-degree neighbor set representing the downwind diffusion area of ​​pollutants. The carbon emission prediction mean tensor at the current moment (in parts per million, i.e., volume ratio) is used as a scalar field, combined with the horizontal wind speed vector (in meters per second) from the collected near-surface turbulent field data as a vector field. For each pair of connected edges, the convective transport flux is calculated, which represents the mass of greenhouse gases passing through the boundary per unit time. The specific numerical calculation process is as follows: The carbon concentration value of the upstream node, the effective wind speed component projected along the edge direction, and the preset grid cross-sectional area constant (e.g., set to 100 square meters) are multiplied together, and the result is multiplied by a preset unit conversion factor (this factor is determined by the conversion factor of the standard air density constant and the PPM dimension, which is 10). -6 The product of these factors (e.g., kilograms per second) is used to generate a mass flux value in a uniform unit. When constructing the flux balance equation, based on the discretization rules of the finite volume method, the change in gas mass at the target node per unit time step is defined as equal to the total flux input to all in-degree neighbors minus the total flux flowing to all out-degree neighbors, plus the estimated source emissions of the node itself (if the node is marked as an emission source). For example, if the total input flux of the in-degree set is 500 units and the total output flux of the out-degree set is 480 units, then theoretically, the node's own concentration should increase by the mass equivalent to 20 units in the current time step.

[0049] A first-order difference operation is performed on the predicted output tensor along the time axis to calculate the predicted concentration change rate of the target node at the current time step. This rate is then multiplied by the unit conversion factor and the control volume constant (i.e., the effective spatial volume covered by the node) to convert it into a mass change rate. This rate is then compared with the net flux difference calculated in the flux balance equation (i.e., inflow flux minus outflow flux plus source term emissions). The algebraic difference between the two is the flux divergence residual, which directly reflects the degree to which the model prediction results violate the law of mass conservation. The absolute values ​​of the flux divergence residuals for all nodes in the network are taken and summed to generate the total physical violation amount.

[0050] To transform this hard constraint into an optimizable soft constraint, Lagrange multipliers are introduced as penalty weights. In the early stages of training, these coefficients are initialized to 0 or a minimum value (0.01 in this example) to allow the model to preferentially learn the data distribution. As the training epochs increase, a linear growth strategy is employed to gradually increase these coefficients (e.g., increasing by 0.05 every 10 epochs until an upper limit of 10.0), forcing the model parameters to converge to a physically consistent manifold during gradient descent. The resulting physical constraint penalty term is directly added to the total loss function, and the weight matrix is ​​corrected using the backpropagation algorithm. This ensures that the optimized input tensor satisfies data fitting requirements while adhering to the physical conservation laws of atmospheric transport.

[0051] In this embodiment, an incremental PID feedback mechanism is used to dynamically adjust the Lagrange multiplier coefficients. Based on the difference between the current physical violation and the violation in the previous iteration, the proportional and derivative terms are adjusted to adaptively update the multiplier coefficients, preventing training oscillations caused by linear growth. To address the optimization overshoot or slow convergence issues that may result from a fixed linear growth strategy, PID control logic is introduced to update the Lagrange multiplier coefficients. The proportional gain is initialized to 0.05, the derivative gain to 0.02, and the integral gain to 0.01. At the end of each backpropagation iteration, the current physical violation is obtained, which is the sum of the absolute values ​​of the flux divergence residuals. The rate of change of the violation is calculated, which is the difference between the current violation and the previous violation, and the historical cumulative sum of the violations is also calculated. The Lagrange multiplier coefficients for the next round are equal to the coefficients of the current round plus the sum of the proportional, differential, and integral term increments. The proportional term increment is the proportional gain multiplied by the current physical violation, the differential term increment is the differential gain multiplied by the rate of change of the violation, and the integral term increment is the integral gain multiplied by the historical cumulative sum of the violations. The upper and lower bounds of the multiplier coefficients are set to 0.1 and 20.0, respectively. The specific adjustment logic is as follows: when the physical violation is large and decreases slowly, the integral term intervenes, accelerating the increase of the coefficients and forcing the network to prioritize physical constraints; when the physical violation decreases sharply (with a large negative rate of change), the differential term generates negative feedback, suppressing the excessively rapid growth of the coefficients. This dynamic adjustment mechanism allows the model to prioritize fitting the data distribution in the early stages of training and automatically lock onto the physically conserved manifold in the later stages, avoiding gradient explosion caused by applying excessive penalties too early and improving convergence speed.

[0052] By calculating the flux divergence residual between the sum of inflow flux and the sum of outflow flux, and defining the absolute value of the flux divergence residual as the physical violation quantity, the physical conservation law is quantified into a computable physical constraint penalty term. This penalty term is then multiplied by the Lagrange multiplier coefficient and incorporated into the optimization objective, ensuring the physical interpretability of the anomaly detection results.

[0053] Further, with the weight parameters of the spatiotemporal graph neural network model fixed, the gradient vector of the composite objective function with respect to the network input carbon emission tensor is calculated; based on the gradient vector, backpropagation is performed to update the network input carbon emission tensor, and the optimized input tensor is output; this corresponds to step S4 above; see [link to relevant documentation]. Figure 3 The specific implementation process includes: Initialize the iteration counter and step size hyperparameter; calculate the partial derivative tensor of the composite objective function with respect to the network input carbon emission tensor, and perform a weighted operation on the partial derivative tensor using the step size hyperparameter to generate an updated gradient tensor; perform a tensor subtraction operation between the network input carbon emission tensor and the updated gradient tensor to generate an intermediate updated tensor; construct a hard threshold truncation function defining the effective numerical domain, and perform a numerical projection operation on the intermediate updated tensor using the hard threshold truncation function to forcibly map elements outside the effective numerical domain to boundary values, generating an updated network input carbon emission tensor; calculate the numerical variation of the composite objective function, and generate a carbon emission iteration loop trigger signal when the numerical variation is greater than a preset convergence threshold; in response to the carbon emission iteration loop trigger signal, increment the iteration counter and calculate the partial derivative tensor until the numerical variation is no greater than the convergence threshold, and output the optimized input tensor.

[0054] Specifically, the process of generating the optimized input tensor is as follows: After the spatiotemporal graph neural network model completes pre-training and locks its weight parameters, in order to clean or reconstruct potentially anomalous input data, the iteration counter is initialized to 0, and the step size hyperparameter (i.e., the learning rate of the input space) is set, for example, to 0.1, which can be adaptively set according to the data variance. The network input carbon emission tensor is marked as a differentiable leaf node, while the weight parameters of all convolutional layers, gating units, and fully connected layers of the model are set to an untrainable state. Subsequently, a complete forward propagation is performed to calculate the value of the composite objective function (including reconstruction error and physical constraint penalty). Using the automatic differentiation engine (in this embodiment, PyTorch's Autograd mechanism), the chain rule is applied to the composite objective function with respect to the network input carbon emission tensor to calculate the partial derivative tensor corresponding to each input element (such as the carbon dioxide reading of a sensor at a certain moment).

[0055] After obtaining the partial derivative tensor, a weighted operation is performed on the partial derivative tensor using a preset step size hyperparameter to generate an updated gradient tensor. This process is similar to parameter updates in neural network training, but here the object of the update is the data itself. An intermediate updated tensor is generated by performing tensor subtraction, i.e., subtracting the updated gradient tensor (if the goal is to minimize the error) from the current network input carbon emission tensor. For example, if the monitoring value of a node is abnormally high (e.g., 2000 PPM), and the physical constraints indicate a huge unexplained increase in mass residuals at that point, the calculated gradient is usually positive, indicating that the input data should be adjusted downwards. Through this tensor subtraction step, a preliminary version of the data, after removing some noise or correcting anomalies, is generated—the intermediate updated tensor.

[0056] Since gradient descent at the purely mathematical level may result in physically meaningless values ​​in the input data (such as negative carbon concentrations or extreme values ​​exceeding the sensor's range), a hard threshold cutoff function defining the effective domain must be constructed. The lower bound of this effective domain is typically set to 0 (representing no emissions), and the upper bound is set to the upper limit of the sensor's physical range (e.g., 5000 PPM) or the maximum possible boundary of historical observations (e.g., 6000 PPM). This function is used to perform numerical projection operations on the intermediate update tensor: for any element less than zero, it is forced to map to zero to conform to the physical axiom of non-negativity of concentration; for any element greater than the upper limit, it is forced to map to the upper limit value to prevent numerical explosion or deviation from the sensor's linear response region due to gradient accumulation during backpropagation; elements within the effective domain remain unchanged, generating the updated network input carbon emission tensor.

[0057] Calculate the numerical change of the composite objective function before and after this iteration (i.e., the decrease in the loss value). When the numerical change is greater than a preset convergence threshold (e.g., 10), the convergence threshold is determined. -4 If the iteration counter has not reached the maximum preset number of rounds (e.g., 100 rounds), a carbon emission iteration loop trigger signal is generated, indicating that the optimization is not yet complete. In response to this signal, the iteration counter is incremented, and the updated network input carbon emission tensor after projection correction is used as the new starting point. It is then input into the spatiotemporal graph neural network with frozen weights, and the forward propagation and backward gradients are recalculated to enter the next round of the loop. This process is repeated until the decrease in the objective function is negligible (not greater than the convergence threshold), indicating that the input data has been adjusted to the optimal state under the current physical constraints and statistical laws. At this point, the loop terminates, and the final optimized input tensor is output.

[0058] By constructing a hard threshold truncation function that defines the effective numerical domain, and using this function to perform numerical projection operations on the intermediate update tensor, elements exceeding the effective numerical domain are forcibly mapped to boundary values. This ensures that the updated network input carbon emission tensor always remains within a reasonable physical range. Combined with an iterative loop trigger signal based on numerical changes, until numerical convergence is achieved, the optimized input tensor is output. This process improves the physical validity and numerical stability of the reconstructed data, resolving the data divergence problem that may result from unconstrained optimization.

[0059] This invention provides a method for detecting anomaly carbon emission data based on a spatiotemporal graph neural network. It collects greenhouse gas concentration and spatial elevation topological coordinates from a three-dimensional sensor network within an industrial park, as well as near-surface turbulent field data from a micro-weather station array. Through high-dimensional feature mapping, a carbon emission tensor is generated as the network input. Nonlinear diffusion trajectories are calculated based on the near-surface turbulent field data, and an adjacency matrix is ​​constructed using the spatial elevation topological coordinates. The network input carbon emission tensor and the adjacency matrix are loaded into a spatiotemporal graph neural network model. Graph convolution operators and gated recurrent units are used to extract spatial and temporal features, capturing the spatiotemporal evolution of carbon emissions. A composite objective function is constructed based on the network input carbon emission tensor and the predicted output tensor, including a physical constraint penalty term. This transforms the gas mass conservation equation into an intrinsic constraint of the model, ensuring the prediction results conform to physical laws. By fixing the weight parameters of the spatiotemporal graph neural network model, the gradient vector of the composite objective function with respect to the network input carbon emission tensor is calculated, and backpropagation is performed to update the function. This allows for iterative repair of anomaly data through gradient optimization until the composite objective function converges, outputting an updated optimized input tensor and improving the detection accuracy of anomaly carbon emission detection.

[0060] Example 2 This embodiment illustrates a practical application of a carbon emission anomaly detection method based on spatiotemporal graph neural networks in a low-carbon demonstration zone of a coastal port logistics park. The specific implementation process is as follows: The area covers land-based storage yards and near-shore berths, and is equipped with 80 high-precision non-dispersive infrared gas sensors and ultrasonic anemometers. Some sensors are installed on the top of quay cranes tens of meters high, forming a typical three-dimensional monitoring network. First, mode decomposition was performed on the collected minute-level carbon dioxide concentration data to extract the time-domain trend component reflecting the port's tidal operation cycle and the frequency-domain fluctuation component affected by strong sea winds. Simultaneously, near-surface turbulence field data, including horizontal wind speed, vertical shear, and turbulence intensity, were input into a pre-trained environmental encoder and mapped into a high-dimensional environmental state embedding vector. Then, the time-domain trend component, frequency-domain fluctuation component, and environmental state embedding vector were concatenated along the feature dimensions to generate a network input carbon emission tensor containing 80 nodes, 6 time steps, and 18 feature dimensions. This tensor fully records the complex three-dimensional emission characteristics of the port area.

[0061] Considering the variable wind direction and significant three-dimensional hierarchical differences in port areas (such as between towering quay cranes and ground-level storage yards), this approach moves away from static geographical distances and instead constructs real-time connectivity based on three-dimensional turbulent diffusion. First, a Gaussian plume model is driven by real-time acquired turbulent field data to calculate the nonlinear diffusion trajectory of emission sources. Next, combining the spatial elevation topological coordinates of the quay cranes and storage yards, the projected distance and lateral offset distance between each node pair are calculated, quantifying the coverage probability of source nodes over target nodes. Based on this, a gravity settlement factor is introduced to correct the vertical connection weights, particularly for the path of diffusion from the top of the towering quay cranes downwind to the ground-level storage yards, where weight enhancement calculations are performed. Finally, the coverage probability and settlement correction coefficient are multiplied element-wise to generate an adjacency matrix that conforms to the three-dimensional flow field characteristics of sea and land winds.

[0062] After loading the constructed network input carbon emission tensor and adjacency matrix into the model, the spatial dimension matrix left multiplication operation is first performed on the input tensor using the normalized Laplacian matrix. For example, berth nodes downwind will highly aggregate emission features from upwind storage yard nodes. The aggregated spatial feature tensor is input into a gated recurrent unit, incorporating the memory state from the previous time step. After concatenation with the current features, the reset and update gate values ​​are calculated through a fully connected layer and a sigmoid activation function. The reset gate determines how much old historical accumulation (such as nighttime background concentration) needs to be forgotten, while the update gate controls the proportion of new observations at the current time step (such as high emissions at the moment a ship docks) written into the new memory state. After nonlinear activation using a hyperbolic tangent function, the hidden state at the current time step is output, and the carbon emission prediction mean tensor and the confidence interval tensor representing uncertainty are predicted in parallel through a linear regression layer.

[0063] To ensure that the model predictions do not violate the axiom of mass conservation, a composite objective function incorporating physical constraints was constructed. Based on the predicted concentration field gradient and wind field vector, the carbon mass flux passing through the boundary of each sensor control volume was calculated. Specifically, for any node, the sum of the mass transfers from all its in-degree neighbors was calculated, and the sum of the mass transfers to all its out-degree neighbors was subtracted to obtain the net mass flux change rate. Simultaneously, the first-order difference of the predicted concentration over time at that node was calculated and multiplied by the volume constant of the control volume and an air density conversion factor including a PPM conversion factor to convert it into the time-varying rate of internal mass. According to the atmospheric convection-diffusion equation, this time-varying rate of internal mass should equal the net mass flux change rate plus the source term emission rate. The difference between these two was calculated, and its absolute value was taken as the physical violation residual. This residual was multiplied by a Lagrange multiplier that increases linearly with the training epochs and added to the mean squared error term of the data fitting to form the total loss function.

[0064] During the model inference phase, if the reconstruction error of a node exceeds a dynamic threshold, the input optimization process is triggered. First, all weight parameters of the neural network are locked, making them unupdateable. Then, the partial derivative of the composite objective function with respect to the network input carbon emission tensor is calculated. Using the direction indicated by this partial derivative, a reverse correction is performed on the original input tensor with a preset step size (e.g., 0.05). Specifically, the intermediate correction value is generated by subtracting the product of the gradient value and the step size from the current input value. To ensure that the corrected data conforms to physical reality, a hard threshold truncation is applied to the intermediate correction value: concentration values ​​less than 0 are forcibly corrected to zero, and values ​​exceeding the sensor's range are truncated to the upper limit of the range. After dozens of rounds of "calculate gradient, update input, truncate values," the objective function value no longer decreases significantly. At this point, the output optimized input tensor is the emission data after removing sensor noise and non-physical interference.

[0065] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for detecting carbon emission anomalies based on spatiotemporal graph neural networks, characterized in that, include: The greenhouse gas concentration and spatial elevation topological coordinates of the three-dimensional sensor network in the industrial park, as well as the near-surface turbulent field data of the micro-weather station array, are collected. Through feature high-dimensional mapping, the network input carbon emission tensor is generated. The nonlinear diffusion trajectory is calculated based on near-surface turbulent field data, and an adjacency matrix is ​​constructed by combining spatial elevation topological coordinates. The carbon emission tensor and adjacency matrix of the network input are loaded into the spatiotemporal graph neural network model. The graph convolution operator is used to perform information aggregation on the carbon emission tensor of the network input to generate the hidden state tensor. Then, the gated recurrent unit is used to perform temporal evolution on the hidden state tensor to generate the predicted output tensor. A composite objective function is constructed based on the network input carbon emission tensor and the predicted output tensor. The composite objective function includes a reconstruction error term for the Euclidean distance between the network input carbon emission tensor and the predicted output tensor, and transforms the gas mass conservation equation into a physical constraint penalty term for physical constraints. With the weight parameters of the fixed spatiotemporal graph neural network model, the gradient vector of the composite objective function with respect to the network input carbon emission tensor is calculated; based on the gradient vector, backpropagation is performed to update the network input carbon emission tensor, and the optimized input tensor is output.

2. The carbon emission anomaly data detection method based on spatiotemporal graph neural network according to claim 1, characterized in that, The specific generation process of the network input carbon emission tensor and adjacency matrix includes: constructing a three-dimensional sensor network at ground monitoring stations, high-altitude towers, and the top of equipment in the industrial park to obtain greenhouse gas concentrations and simultaneously analyze the spatial elevation topological coordinates of each sensor node; deploying micro-weather station arrays around wind ducts and emission sources to obtain near-surface turbulent field data including horizontal wind speed vectors, vertical wind shear indices, and atmospheric turbulence intensity; extracting time-domain trend components and frequency-domain fluctuation components from the greenhouse gas concentration; mapping the near-surface turbulent field data into an environmental state embedding vector using an environmental encoder network; and converting the time-domain data into an environmental state embedding vector. Tensor concatenation operations are performed on the trend component, frequency domain fluctuation component, and environmental state embedding vector to generate the network input carbon emission tensor. Based on the horizontal wind speed vector and atmospheric turbulence intensity in the near-surface turbulence field data, the Gaussian plume model is driven to solve the nonlinear diffusion trajectory of greenhouse gases. Combining spatial elevation topological coordinates, the projected distance and lateral offset distance of the target node relative to the source node on the nonlinear diffusion trajectory are calculated. The coverage probability of the target node is calculated based on the projected distance and lateral offset distance, and the vertical connection weight is corrected by introducing a gravity settling factor to generate an adjacency matrix.

3. The carbon emission anomaly data detection method based on spatiotemporal graph neural network according to claim 1, characterized in that, The specific computational process of the graph convolution operator includes: reading a preset node activity state vector, applying the node activity state vector as a mask to the adjacency matrix to generate a carbon emission topology matrix; accumulating the element values ​​of the carbon emission topology matrix row by row to generate a carbon emission degree matrix; performing a normalization operation on the adjacency matrix based on the carbon emission degree matrix to generate a normalized Laplacian matrix; initializing and configuring a trainable weight matrix and a bias vector; performing a matrix left multiplication operation on the network input carbon emission tensor using the normalized Laplacian matrix to generate a carbon emission spatial aggregation feature tensor; performing a linear mapping operation on the carbon emission spatial aggregation feature tensor using the trainable weight matrix and superimposing the bias vector to generate a linear transformation feature tensor; inputting the linear transformation feature tensor into a nonlinear activation function layer to perform an element-wise activation operation, and outputting a hidden layer state tensor.

4. The carbon emission anomaly data detection method based on spatiotemporal graph neural network according to claim 1, characterized in that, The specific computation process of the gated recurrent unit includes: introducing the carbon emission memory state tensor from the previous time step, performing concatenation and linear transformation with the hidden layer state tensor, and generating a reset gate tensor and an update gate tensor via a Sigmoid activation function; performing a Hadamard product operation on the carbon emission memory state tensor from the previous time step and the reset gate tensor to generate a retained memory tensor; fusing the retained memory tensor with the hidden layer state tensor and inputting it into a hyperbolic tangent activation function layer to generate a candidate state tensor; performing a linear weighted combination operation on the carbon emission memory state tensor from the previous time step and the candidate state tensor according to the weight ratio defined by the update gate tensor, outputting the carbon emission memory state tensor from the current time step; inputting the carbon emission memory state tensor from the current time step into a fully connected distribution mapping layer to perform a linear regression operation, outputting a carbon emission prediction mean tensor and a confidence interval tensor in parallel, and combining the carbon emission prediction mean tensor and the confidence interval tensor into a prediction output tensor.

5. The carbon emission anomaly data detection method based on spatiotemporal graph neural network according to claim 1, characterized in that, The specific generation process of the reconstruction error term includes: using the network input carbon emission tensor as the minuend and the predicted mean tensor in the predicted output tensor as the subtrahend, performing element-wise subtraction to generate the original carbon emission residual tensor; extracting the numerical boundary of the confidence interval tensor to calculate the interval span, and performing a reciprocal operation on the interval span to generate an inverse proportional weight tensor; performing a square operation on the original carbon emission residual tensor, and performing a Hadamard product operation on the result and the inverse proportional weight tensor to generate a confidence error tensor; and performing a full-dimensional summation operation on the confidence error tensor to generate the reconstruction error term.

6. The carbon emission anomaly data detection method based on spatiotemporal graph neural network according to claim 1, characterized in that, The specific generation process of the physical constraint penalty term includes: based on the node connection relationship defined by the adjacency matrix, extracting the in-degree neighbor set flowing into the target node and the out-degree neighbor set flowing out of the target node; using the carbon emission prediction mean tensor and combining it with the near-surface turbulent field data, calculating the carbon flux values ​​of the target node and the in-degree neighbor set and the out-degree neighbor set, and constructing a flux balance equation; calculating the flux divergence residual between the sum of inflow flux and the sum of outflow flux in the flux balance equation; defining the absolute value of the flux divergence residual as the physical violation amount, and multiplying the physical violation amount by a preset Lagrange multiplier coefficient used as the penalty weight coefficient to generate the physical constraint penalty term.

7. The carbon emission anomaly data detection method based on spatiotemporal graph neural network according to claim 1, characterized in that, The specific process of the backpropagation update operation includes: initializing the iteration counter and step size hyperparameter; calculating the partial derivative tensor of the composite objective function with respect to the network input carbon emission tensor, performing a weighted operation on the partial derivative tensor using the step size hyperparameter to generate an updated gradient tensor; performing a tensor subtraction operation between the network input carbon emission tensor and the updated gradient tensor to generate an intermediate update tensor; constructing a hard threshold truncation function defining the effective numerical domain, performing a numerical projection operation on the intermediate update tensor using the hard threshold truncation function to forcibly map elements exceeding the effective numerical domain to boundary values, generating an updated network input carbon emission tensor; calculating the numerical variation of the composite objective function, and generating a carbon emission iteration loop trigger signal when the numerical variation is greater than a preset convergence threshold; responding to the carbon emission iteration loop trigger signal, accumulating the iteration counter and calculating the partial derivative tensor until the numerical variation is no greater than the convergence threshold, and outputting the optimized input tensor.

Images