Online monitoring method and system for unorganized emission based on edge calculation

Abstract

The invention provides an unorganized emission online monitoring method and system based on edge calculation, and relates to the technical field of environment monitoring, and the method comprises the steps: obtaining emission feature data of distributed sensing nodes, building an inter-node adjacency relation, carrying out the space-time dimension feature extraction, determining an abnormal score, and carrying out the detection of the abnormal score; and an emission source candidate position is searched based on actually measured concentration distribution and an emission state evaluation result, high and low frequency components are extracted through wavelet decomposition to analyze emission characteristics, an influence range boundary is determined in combination with meteorological parameters, and an emission early warning signal and a traceability positioning result are generated. Accurate monitoring, traceability positioning and early warning of unorganized emission are realized, and the environmental monitoring efficiency is improved.

Description

Technical Field

[0001] This invention relates to the field of environmental monitoring technology, and in particular to an online monitoring method and system for fugitive emissions based on edge computing. Background Technology

[0002] With the acceleration of industrialization, fugitive emissions have become a significant source of air pollution. Fugitive emissions are characterized by numerous and widespread locations, high randomness, and difficulty in monitoring, posing a significant challenge to environmental regulation. Traditional methods for monitoring fugitive emissions mainly rely on periodic manual sampling and analysis, which cannot achieve continuous real-time monitoring; while fixed monitoring equipment has limited deployment points, making it difficult to comprehensively capture the spatial distribution characteristics of emission behavior.

[0003] With the development of IoT and edge computing technologies, distributed sensing networks are gradually being applied to the field of environmental monitoring, providing new technical means for monitoring fugitive emissions. By deploying multiple low-cost sensor nodes, high-density coverage of the monitoring area can be achieved, and the data can be transmitted to the edge computing platform in real time for processing and analysis, which greatly improves the spatiotemporal resolution of monitoring. However, existing technologies for monitoring fugitive emissions still have problems such as lack of an effective fusion mechanism for multi-node data, difficulty in making full use of the spatial correlation between distributed nodes, reliance on simple threshold judgment, inability to effectively distinguish between environmental background fluctuations and real emission signals, especially under complex meteorological conditions, which can easily lead to false alarms or missed alarms, and reliance on complex numerical simulations or simplified Gaussian diffusion models, which have large computational load and poor real-time performance, making it difficult to meet the needs of rapid response. Summary of the Invention

[0004] This invention provides a method and system for online monitoring of fugitive emissions based on edge computing, which can at least solve some of the problems existing in the prior art.

[0005] A first aspect of this invention provides a method for online monitoring of fugitive emissions based on edge computing, comprising: Emission feature data from multiple distributed sensing nodes are acquired, and the adjacency relationships and propagation connectivity between the nodes are established by combining the spatial coordinates of the distributed sensing nodes to obtain the emission feature dataset. The emission feature dataset is encoded in the time dimension to obtain a temporal evolution feature vector. In the spatial dimension, the emission feature dataset is subjected to graph convolution and the feature information of adjacent nodes is aggregated based on the adjacency relationship to obtain aggregated spatial features. Based on the temporal evolution feature vector and aggregated spatial features, the emission anomaly score and anomaly type label of each node are determined to obtain the emission status assessment result. Based on the measured concentration distribution and emission status assessment results in the emission feature dataset, the candidate locations of emission sources are obtained by iterative search using the particle swarm optimization algorithm. The concentration time series data within a preset distance of the candidate locations of emission sources are extracted and wavelet decomposition is performed to obtain high-frequency fluctuation components and low-frequency trend components. The pulse characteristics of emissions are determined based on the high-frequency fluctuation components, the emission intensity variation pattern is identified based on the low-frequency trend components, and the theoretical diffusion distance of pollutants in each wind sector is calculated and the boundary of the influence range is determined by combining the candidate emission source locations and real-time meteorological parameters. Based on the emission intensity variation pattern and the influence range boundary, and combined with preset emission threshold constraints and duration determination criteria, an emission early warning signal and source location result are generated.

[0006] In one alternative implementation, The emission feature dataset is obtained by acquiring emission feature data from multiple distributed sensing nodes and establishing adjacency relationships and propagation connectivity between the nodes using their spatial coordinates. The dataset includes: The emission characteristic data reported by each distributed sensing node is obtained and the timestamp information and spatial coordinate information are extracted. Based on the timestamp information, the emission characteristic data of different distributed sensing nodes are aligned with the time reference to obtain time-aligned data. Based on the spatial coordinate information, the Euclidean distance between each distributed sensing node is calculated and neighboring node pairs are filtered according to a preset distance threshold. Connecting edges are established between neighboring node pairs that meet the distance conditions and an undirected graph is constructed to obtain the adjacency relationship. Acquire the prevailing wind direction data of the monitoring area and calculate the azimuth difference between each pair of neighboring nodes. Use the azimuth difference and the prevailing wind direction data to determine the neighboring node pairs in the downwind propagation link and mark the propagation connectivity. Combine the time alignment data and the adjacency relationship to obtain the emission characteristic dataset through correlation integration.

[0007] In one alternative implementation, Encoding the emission feature dataset in the time dimension yields a temporal evolution feature vector. Then, graph convolution is performed on the emission feature dataset in the spatial dimension, and feature information from adjacent nodes is aggregated based on the adjacency relationships to obtain aggregated spatial features, including: The observation values ​​of each distributed sensing node at multiple moments within a continuous time window are extracted from the emission feature dataset to form node time series data. The node time series data is sinusoidally and cosinely encoded to generate time position vectors and concatenated to obtain a time-labeled sequence. The query vector, key vector, and value vector corresponding to the time-labeled sequence are constructed along the time dimension. The dot product of the query vector and the key vector is calculated to obtain the time attention weight. The time attention weight and the value vector are weighted and summed, and residual connections and layer normalization are performed to obtain the intermediate result of time encoding. The deep semantic representation corresponding to the intermediate result of time encoding is extracted through a preset feedforward neural network and global average pooling is performed along the time dimension to obtain the time series evolution feature vector. An adjacency matrix is ​​constructed based on the adjacency relationship. The degree matrix corresponding to the adjacency matrix is ​​calculated and symmetric normalized to obtain a normalized adjacency matrix. The current feature values ​​of each node in the emission feature dataset are used to construct a node feature matrix. The node feature matrix is ​​linearly transformed with a preset initial weight matrix to obtain a transformed feature matrix. The normalized adjacency matrix and the transformed feature matrix are multiplied to obtain neighborhood aggregation features. A nonlinear activation function is applied to the neighborhood aggregation features and layer normalization is performed to obtain aggregation space features.

[0008] In one alternative implementation, Based on the temporal evolution feature vector and aggregated spatial features, the emission anomaly score and anomaly type label of each node are determined to obtain the emission status assessment results, including: The temporal evolution feature vector is concatenated with the aggregated spatial features to obtain a fused feature vector. The fused feature vector is then transformed in terms of feature dimension and mapped through a fully connected layer to obtain an intermediate representation vector. The intermediate representation vector is then weighted by channel attention along the feature dimension and key feature components are extracted. The key feature components are then input into the anomaly score prediction branch to obtain continuous anomaly score values ​​for each node. The continuous anomaly score values ​​are then normalized and mapped to a preset score interval to obtain a standardized anomaly score. The intermediate representation vector is input into the multi-class prediction branch and nonlinearly transformed through a multi-layer fully connected network to obtain the unnormalized classification vector corresponding to each anomaly type. The unnormalized classification vector is then subjected to softmax normalization to obtain the anomaly type probability distribution. The category index with the highest probability is extracted from the anomaly type probability distribution as the anomaly type label, and the highest probability value is extracted as the category confidence. The category confidence and the standardized anomaly score are weighted and fused to obtain a comprehensive anomaly index. The standardized anomaly score is combined with the anomaly type label, and the comprehensive anomaly index is associated with the node identifier and timestamp information of each node to obtain the emission status assessment result.

[0009] In one alternative implementation, Based on the measured concentration distribution and emission status assessment results in the emission feature dataset, candidate emission source locations are obtained through iterative search using a particle swarm optimization algorithm. Concentration time-series data within a preset distance of the candidate emission source locations are extracted and decomposed using wavelet decomposition to obtain high-frequency fluctuation components and low-frequency trend components, including: The measured concentration values ​​are extracted from the emission feature dataset and combined with the corresponding spatial coordinates to construct a concentration spatial distribution field. The set of nodes whose standardized anomaly scores exceed the preset anomaly threshold is extracted from the emission status assessment results, and the initial search center position is determined. The initial distribution of the particle swarm is initialized based on the initial search center position and the preset spatial constraints, and the position vector and velocity vector of each particle are set. Based on the current position of each particle, the concentration contribution value is determined in combination with the preset distance decay relationship and compared with the corresponding measured concentration values ​​in the concentration spatial distribution field to obtain the fitness value. The individual optimal position and global optimal position of the particle are updated according to the fitness value. The velocity update amount of the particle is determined in combination with the preset inertia weight. Boundary constraints are applied to the velocity update amount and iterative search is performed until the preset maximum number of iterations is reached. The global optimal position is used as the candidate position of the emission source. The concentration time series data is obtained by determining the set of neighboring nodes within a preset distance range centered on the candidate emission source location and extracting the concentration observation sequence within a continuous time window. The concentration time series data is then separated into an initial high-frequency component and an initial low-frequency component through multi-scale wavelet decomposition. The initial high-frequency component is subjected to threshold denoising processing to obtain a high-frequency fluctuation component, and the initial low-frequency component is subjected to trend fitting to obtain a low-frequency trend component.

[0010] In one alternative implementation, The pulse characteristics of emissions are determined based on the high-frequency fluctuation components, the emission intensity variation patterns are identified based on the low-frequency trend components, and the theoretical diffusion distance of pollutants in each wind sector is calculated and the boundary of the influence range is determined by combining the candidate emission source locations and real-time meteorological parameters. Peak detection is performed on the high-frequency fluctuation component to extract the peak position and peak intensity. The time interval between the peak positions is counted to obtain the pulse frequency. The rising slope and falling slope corresponding to the peak intensity are extracted and combined to obtain the response rate. The pulse frequency and the response rate are combined to obtain the emission pulse characteristics. The low-frequency trend component is fitted with a polynomial to obtain a baseline curve, and the derivative is used to obtain a rate of change curve. The rate of change curve is transformed in the frequency domain and the dominant frequency component is extracted to obtain the dominant period. The deviation amplitude and deviation peak between the low-frequency trend component and the baseline curve are determined to obtain the fluctuation characteristics. The dominant period is compared with the pulse frequency to obtain the period matching degree. Based on the period matching degree and the fluctuation characteristics, the intermittent characteristics are determined. The emission intensity change law is obtained by combining the evolution trend of the baseline curve. Real-time meteorological data is acquired and the wind direction angle corresponding to the candidate emission source location is determined. Based on the wind direction angle, wind sector is divided and the average wind speed in each wind sector is calculated to obtain the sector wind speed. The atmospheric stability level is determined based on the real-time meteorological data and the diffusion coefficient is determined in combination with the sector wind speed. The emission source intensity correction coefficient is calculated based on the pulse characteristics and the emission intensity variation law. The theoretical diffusion distance of pollutants in each wind sector is calculated in combination with the diffusion coefficient and a fan-shaped region is generated. The boundary coordinate points of the fan-shaped region are connected to obtain the boundary of the influence range.

[0011] In one alternative implementation, Based on the emission intensity variation pattern and the boundary of the impact range, and combined with preset emission threshold constraints and duration determination criteria, a system is generated that includes emission early warning signals and source location results, including: The current emission intensity value and intensity change trend are extracted from the emission intensity change pattern. The current emission intensity value is compared with the preset emission threshold constraint to obtain the exceeding state. The duration of continuous exceeding is counted and the persistence judgment result is determined according to the preset duration judgment criterion. The duration is predicted based on the intensity change trend to obtain the expected duration. The warning level is determined according to the expected duration and the persistence judgment result and combined with the exceeding state to obtain the emission warning signal. The locations of monitoring points within the boundary of the influence range are extracted, the spatial distance between each monitoring point and the candidate emission source is calculated, and the wind sector to which each monitoring point belongs is determined. The measured concentration values ​​of each monitoring point are spatially interpolated with the current emission intensity value to obtain an interpolated concentration field. The gradient corresponding to the interpolated concentration field is calculated, and the direction of maximum gradient is extracted. The source coordinates are obtained by tracing back along the direction of maximum gradient to the center of the boundary of the influence range. The source coordinates are fused with the warning level to obtain the source location result. The emission warning signal and the source location result are encapsulated and output.

[0012] A second aspect of this invention provides an online monitoring system for fugitive emissions based on edge computing, comprising: The data acquisition unit is used to acquire emission characteristic data from multiple distributed sensing nodes and combine the spatial coordinates of the distributed sensing nodes to establish the adjacency relationship and propagation connectivity between the nodes to obtain the emission characteristic dataset. The feature extraction unit is used to encode the emission feature dataset in the time dimension to obtain a temporal evolution feature vector, perform graph convolution on the emission feature dataset in the spatial dimension and aggregate the feature information of adjacent nodes based on the adjacency relationship to obtain aggregated spatial features, and determine the emission anomaly score and anomaly type label of each node based on the temporal evolution feature vector and aggregated spatial features to obtain the emission status assessment result. The source location unit is used to obtain candidate emission source locations through iterative search using a particle swarm optimization algorithm based on the measured concentration distribution and emission status assessment results in the emission feature dataset, extract the concentration time series data within a preset distance of the candidate emission source locations and perform wavelet decomposition to obtain high-frequency fluctuation components and low-frequency trend components. The characteristic analysis unit is used to determine the pulse characteristics of emissions based on the high-frequency fluctuation component, identify the emission intensity change pattern based on the low-frequency trend component, and calculate the theoretical diffusion distance of pollutants in each wind sector and determine the boundary of the influence range by combining the candidate emission source locations and real-time meteorological parameters. The result output unit is used to generate emission warning signals and source location results based on the emission intensity change pattern and the boundary of the influence range, combined with preset emission threshold constraints and duration determination criteria.

[0013] A third aspect of the present invention provides an electronic device, comprising: A processor and a memory for storing processor-executable instructions, wherein the processor is configured to invoke instructions stored in the memory to perform the aforementioned method.

[0014] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0015] In this invention, by establishing adjacency relationships and propagation connectivity between nodes, efficient integration of distributed sensing data is achieved, overcoming the problem that traditional single-point monitoring is unable to cover large-scale unorganized emissions, and improving the spatial coverage capability of the monitoring network. Feature extraction is performed by combining temporal dimension encoding and spatial dimension graph convolution, which can simultaneously capture the temporal evolution and spatial diffusion characteristics of emissions, significantly improving the accuracy of abnormal emission identification. Particle swarm optimization algorithm is used to iteratively search for candidate emission source locations, and wavelet decomposition is combined with concentration time-series data analysis to achieve precise source localization and comprehensive characterization of emission characteristics. Through the separation analysis of high-frequency fluctuation components and low-frequency trend components, the pulse characteristics and intensity variation patterns of emissions can be accurately identified. Combined with real-time meteorological parameters to calculate the theoretical diffusion distance, the dynamic boundary of the pollutant's influence range is determined, providing a scientific basis for emergency response. Attached Figure Description

[0016] Figure 1 This is a schematic flowchart of the online monitoring method for fugitive emissions based on edge computing, according to an embodiment of the present invention. Figure 2 This is a flowchart illustrating the pollution source location analysis process of the online monitoring method for fugitive emissions based on edge computing, as described in this embodiment of the invention. Detailed Implementation

[0017] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, 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.

[0018] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0019] Figure 1 This is a schematic flowchart of the online monitoring method for fugitive emissions based on edge computing according to an embodiment of the present invention. Figure 1 As shown, the method includes: Emission feature data from multiple distributed sensing nodes are acquired, and the adjacency relationships and propagation connectivity between the nodes are established by combining the spatial coordinates of the distributed sensing nodes to obtain the emission feature dataset. The emission feature dataset is encoded in the time dimension to obtain a temporal evolution feature vector. In the spatial dimension, the emission feature dataset is subjected to graph convolution and the feature information of adjacent nodes is aggregated based on the adjacency relationship to obtain aggregated spatial features. Based on the temporal evolution feature vector and aggregated spatial features, the emission anomaly score and anomaly type label of each node are determined to obtain the emission status assessment result. Based on the measured concentration distribution and emission status assessment results in the emission feature dataset, the candidate locations of emission sources are obtained by iterative search using the particle swarm optimization algorithm. The concentration time series data within a preset distance of the candidate locations of emission sources are extracted and wavelet decomposition is performed to obtain high-frequency fluctuation components and low-frequency trend components. The pulse characteristics of emissions are determined based on the high-frequency fluctuation components, the emission intensity variation pattern is identified based on the low-frequency trend components, and the theoretical diffusion distance of pollutants in each wind sector is calculated and the boundary of the influence range is determined by combining the candidate emission source locations and real-time meteorological parameters. Based on the emission intensity variation pattern and the influence range boundary, and combined with preset emission threshold constraints and duration determination criteria, an emission early warning signal and source location result are generated.

[0020] In one alternative implementation, The emission feature dataset is obtained by acquiring emission feature data from multiple distributed sensing nodes and establishing adjacency relationships and propagation connectivity between the nodes using their spatial coordinates. The dataset includes: The emission characteristic data reported by each distributed sensing node is obtained and the timestamp information and spatial coordinate information are extracted. Based on the timestamp information, the emission characteristic data of different distributed sensing nodes are aligned with the time reference to obtain time-aligned data. Based on the spatial coordinate information, the Euclidean distance between each distributed sensing node is calculated and neighboring node pairs are filtered according to a preset distance threshold. Connecting edges are established between neighboring node pairs that meet the distance conditions and an undirected graph is constructed to obtain the adjacency relationship. Acquire the prevailing wind direction data of the monitoring area and calculate the azimuth difference between each pair of neighboring nodes. Use the azimuth difference and the prevailing wind direction data to determine the neighboring node pairs in the downwind propagation link and mark the propagation connectivity. Combine the time alignment data and the adjacency relationship to obtain the emission characteristic dataset through correlation integration.

[0021] The edge computing unit receives emission characteristic data reported by each distributed sensing node. Each data point includes timestamp information, spatial coordinate information, and emission characteristic parameters such as the concentration of various gases. Emission characteristic parameters may include environmental parameters such as particulate matter concentration, volatile organic compound concentration, sulfur dioxide concentration, and nitrogen oxide concentration, which are used to characterize the features of fugitive emissions.

[0022] Edge computing units extract timestamp and spatial coordinate information from the raw data. The timestamp information is recorded in the format of year, month, day, hour, minute, and second, such as "2024-12-31 15:30:45", accurate to the second. The spatial coordinate information includes three-dimensional information: longitude, latitude, and altitude. Latitude and longitude are expressed in decimal units, and altitude is expressed in meters relative to the ground, such as "12.5".

[0023] To ensure temporal consistency of data across nodes, the edge computing unit aligns the emission characteristic data of different distributed sensing nodes to a time reference based on timestamp information. The time alignment process sets a reference time point, such as the hour of each minute, and then transforms the data from each node within a similar time window (e.g., ±15 seconds) to the reference time point using linear interpolation. For example, if the first node measures a volatile organic compound (VOC) concentration of 120 ppb at 15:00:10 and a VOC concentration of 150 ppb at 15:00:40, then the concentration at the reference time point of 15:00:00, calculated through linear interpolation, is approximately 110 ppb. This aligns the data from all nodes to the same time series, forming a time-aligned dataset.

[0024] The edge computing unit calculates the Euclidean distance between distributed sensing nodes based on spatial coordinate information. For example, for any two nodes, designated as the first node and the second node, with spatial coordinates (A, B, C) and (D, E, F) respectively, the three-dimensional Euclidean distance between the two nodes is calculated. After distance conversion, the actual distance is approximately 142 meters. Based on a preset distance threshold, such as 200 meters, neighboring node pairs are selected. Nodes with a distance less than the threshold are considered neighboring nodes. In this embodiment, the distance between the first node and the second node is less than the preset threshold, therefore they are identified as a neighboring node pair.

[0025] Connecting edges are established between neighboring node pairs that meet the distance condition, and an undirected graph is constructed to represent the adjacency relationship between nodes. The undirected graph is stored in the form of an adjacency matrix or an adjacency list. The element values ​​in the adjacency matrix indicate whether there is a connecting edge between nodes, with 1 indicating a connection and 0 indicating no connection. For example, in a monitoring network containing 10 distributed sensing nodes, approximately 25 pairs of neighboring node relationships may be formed, and the constructed adjacency matrix has approximately 50 non-zero elements.

[0026] Edge computing units acquire prevailing wind direction data for the monitoring area. Wind direction data is represented by angles, with true north as 0 degrees and increasing clockwise. For example, a prevailing northeast wind is recorded as 45 degrees, and a southeast wind as 135 degrees. Wind direction data can come from weather stations within the monitoring area or wind direction sensors integrated into distributed sensing nodes. The azimuth angle between each pair of neighboring nodes is calculated. The azimuth angle is defined as the angle between the line connecting the first node to the second node and true north. For example, if the first node is southwest of the second node, the azimuth angle from the first node to the second node is approximately 45 degrees.

[0027] The edge computing unit calculates the angle difference between the azimuth and the prevailing wind direction. When the angle difference is less than a preset angle threshold (e.g., 30 degrees), it is determined that the two nodes are on the wind propagation path. For example, when the prevailing wind direction is 45 degrees (northeast), and the azimuth of the first node to the second node is 40 degrees, the angle difference is 5 degrees, which is less than the preset threshold, and it is determined that the second node is downwind of the first node. Neighboring node pairs that meet the angle condition are marked as node pairs with propagation connectivity and are specially marked in the adjacency matrix, such as using 2 to represent a connection edge with propagation connectivity.

[0028] Edge computing units integrate time-aligned data and adjacency relationships to generate an emission feature dataset. This dataset includes not only the original emission feature parameters but also the spatiotemporal relationships between nodes and wind propagation link information. Specifically, for each pair of nodes with propagation connectivity, the corresponding time-aligned emission feature data sequence is extracted, and the time lag response pattern of the second node relative to the first node is calculated. For example, under a wind speed of 3 m / s, the theoretical propagation delay between two nodes 150 meters apart is approximately 50 seconds. This propagation delay can be verified and refined by analyzing concentration change patterns in actual data.

[0029] In this embodiment, by introducing a time reference alignment mechanism, emission characteristics from different distributed sensing nodes are compared and fused under a unified time reference, improving the accuracy and stability of emission change analysis. By calculating the distance between nodes based on spatial coordinates and screening neighboring node pairs, an adjacency relationship that conforms to the actual physical distribution characteristics is constructed, effectively reducing noise interference introduced by irrelevant or distant nodes. The directional relationship of neighboring node pairs is analyzed and downwind propagation connectivity is marked, so that the constructed node association not only reflects geometric proximity but also reflects the directionality of pollutant propagation in the actual environment, effectively enhancing the ability to characterize the fugitive emission diffusion process and improving the reliability and discrimination accuracy of emission monitoring results under complex environmental conditions.

[0030] In one alternative implementation, Encoding the emission feature dataset in the time dimension yields a temporal evolution feature vector. Then, graph convolution is performed on the emission feature dataset in the spatial dimension, and feature information from adjacent nodes is aggregated based on the adjacency relationships to obtain aggregated spatial features, including: The observation values ​​of each distributed sensing node at multiple moments within a continuous time window are extracted from the emission feature dataset to form node time series data. The node time series data is sinusoidally and cosinely encoded to generate time position vectors and concatenated to obtain a time-labeled sequence. The query vector, key vector, and value vector corresponding to the time-labeled sequence are constructed along the time dimension. The dot product of the query vector and the key vector is calculated to obtain the time attention weight. The time attention weight and the value vector are weighted and summed, and residual connections and layer normalization are performed to obtain the intermediate result of time encoding. The deep semantic representation corresponding to the intermediate result of time encoding is extracted through a preset feedforward neural network and global average pooling is performed along the time dimension to obtain the time series evolution feature vector. An adjacency matrix is ​​constructed based on the adjacency relationship. The degree matrix corresponding to the adjacency matrix is ​​calculated and symmetric normalized to obtain a normalized adjacency matrix. The current feature values ​​of each node in the emission feature dataset are used to construct a node feature matrix. The node feature matrix is ​​linearly transformed with a preset initial weight matrix to obtain a transformed feature matrix. The normalized adjacency matrix and the transformed feature matrix are multiplied to obtain neighborhood aggregation features. A nonlinear activation function is applied to the neighborhood aggregation features and layer normalization is performed to obtain aggregation space features.

[0031] The edge computing unit extracts observations from multiple moments within a continuous time window from the previously generated emission characteristic dataset to form node time-series data. For each distributed sensing node, monitoring data for 24 consecutive hours is selected, with a sampling interval of 5 minutes, forming an observation sequence containing 288 time points. The observations at each time point include multi-dimensional emission characteristic data such as particulate matter concentration, volatile organic compound (VOC) concentration, and sulfur dioxide concentration. For example, the average particulate matter concentration of the first node between 10:00 and 10:05 on a certain day is 35 micrograms per cubic meter, and the VOC concentration is 75 ppb.

[0032] The edge computing unit performs sine and cosine encoding on the node time-series data to generate time-location vectors, mapping the location index of each time point to a vector in a high-dimensional space. For each location index value, a series of sine and cosine function values ​​at different frequencies are generated. For example, for the 60th time point in the sequence, i.e., 5 hours after the start of sampling, sine and cosine values ​​at different frequencies are calculated, resulting in a 64-dimensional time-location vector. Even-numbered dimensions are encoded using sine functions, and odd-numbered dimensions are encoded using cosine functions, with the frequencies decreasing exponentially. The location vectors of all time points are concatenated to form a time-labeled sequence used to capture temporal patterns in the data.

[0033] The edge computing unit constructs query vectors, key vectors, and value vectors corresponding to the time-labeled sequence along the time dimension, and maps the time-labeled sequence into query vectors, key vectors, and value vectors through three different linear transformation layers. Each linear transformation layer consists of a set of trainable parameter matrices, the dimensions of which match the dimensions of the time-labeled sequence. For example, if the dimension of the time-labeled sequence is 288×64, i.e., 288 time points, each with 64-dimensional features, then the parameter matrix dimension of the linear transformation layer is 64×64, and the generated query vector, key vector, and value vector are all 288×64 in dimension.

[0034] The edge computing unit calculates the dot product of the query vector and the key vector to obtain the temporal attention weights. It then performs matrix multiplication on the query vector and key vector to calculate the relevance score between each pair of time points. For example, for the 60th and 63rd time points in the sequence (15 minutes apart), the dot product of the corresponding query vector and key vector is calculated. If the result is 0.75, it indicates that these two time points have a high relevance. The calculated original attention score matrix has a dimension of 288×288, representing the relevance between each pair of time points in the sequence. To prevent gradient vanishing, the original attention scores are scaled by dividing by 8. Subsequently, a softmax normalization function is applied to normalize the scores, ensuring that the sum of the attention weights of each time point to other time points is 1.

[0035] The edge computing unit performs a weighted summation of the temporal attention weights and the value vectors, followed by residual concatenation and layer normalization to obtain intermediate temporal encoding results. For each time point in the sequence, the value vectors are weighted and summed based on their attention weights relative to all other time points. For example, the output representation at the 60th time point is the weighted sum of the value vectors of all time points based on the attention weights in the 60th row. After the weighted summation, the result is residually concatenated with the original time-labeled sequence, i.e., element-wise added, forming a new feature representation with a dimension of 288×64. Layer normalization is then applied, calculating the mean and standard deviation of each feature dimension, and performing normalization operations to stabilize the training process.

[0036] The edge computing unit extracts the deep semantic representation corresponding to the intermediate time-encoded result through a pre-defined feedforward neural network and performs global average pooling along the time dimension to obtain a temporal evolution feature vector. The feedforward neural network contains two layers of linear transformation and a ReLU activation function. The first layer of linear transformation expands the feature dimension from 64 to 128. After applying the ReLU activation function, the second layer of linear transformation remaps the dimension back to 64. The network output and the intermediate time-encoded result are then subjected to residual connections and layer normalization to obtain a deep semantic representation with a dimension of 288×64. Global average pooling is performed along the time dimension, i.e., the average value of each feature dimension at all time points is calculated, resulting in a 64-dimensional temporal evolution feature vector. This temporal evolution feature vector encodes the overall pattern of node emission features changing over time.

[0037] The edge computing unit constructs an adjacency matrix based on the aforementioned adjacency relationships. This adjacency matrix is ​​a two-dimensional matrix where rows and columns represent nodes in the monitoring network. If a connection exists between the first node and the second node, the corresponding element in the adjacency matrix has a value of 1; otherwise, it has a value of 0. For example, in a monitoring network containing 30 distributed sensing nodes, the adjacency matrix has a dimension of 30×30.

[0038] The edge computing unit calculates the degree matrix corresponding to the adjacency matrix and performs symmetric normalization to obtain a normalized adjacency matrix. The degree matrix is ​​a diagonal matrix, and the elements on the diagonal are the number of connections (degrees) of the corresponding nodes. For example, if the first node is connected to 5 other nodes, the diagonal element corresponding to the first node in the degree matrix is ​​5. The symmetric normalization process multiplies the adjacency matrix by the negative square root of the degree matrix on both sides to obtain the normalized adjacency matrix.

[0039] The edge computing unit constructs a node feature matrix from the current feature values ​​of each node in the emission feature dataset. The number of rows in the node feature matrix equals the number of nodes, and the number of columns equals the feature dimension of each node. For example, for 30 nodes, each node has 8 features (such as the concentrations of different pollutants such as particulate matter, volatile organic compounds, and sulfur dioxide), and the dimension of the node feature matrix is ​​30×8.

[0040] The edge computing unit linearly transforms the node feature matrix with a preset initial weight matrix to obtain a transformed feature matrix. The initial weight matrix has a dimension of 8×16, and the transformed feature matrix after linear transformation has a dimension of 30×16. The edge computing unit performs matrix multiplication between the normalized adjacency matrix and the transformed feature matrix to obtain neighborhood aggregation features. For example, the aggregated features of the first node are a weighted combination of its own features and the features of all connected nodes, with the weights determined by the corresponding elements in the normalized adjacency matrix.

[0041] The edge computing unit applies a nonlinear activation function to the neighborhood aggregation features and performs layer normalization to obtain the aggregation spatial features. The nonlinear activation function uses the ReLU function; for each element in the input features, it remains unchanged when the element value is greater than 0, otherwise it is set to 0. Layer normalization calculates the mean and standard deviation of each dimension of the feature for each node, performing standardization processing to make the feature distribution of different nodes more consistent. The dimension of the aggregation spatial features is the same as that of the neighborhood aggregation features, which is 30×16.

[0042] In this embodiment, by uniformly modeling node observations within a continuous time window and introducing time location encoding and attention mechanisms, it is possible to adaptively focus on time segments with key impacts during emission changes, thereby improving the sensitivity and stability of emission time-series evolution characteristics. By further mapping the time encoding results into deep semantic representations and performing global aggregation, the resulting time-series evolution characteristics can comprehensively reflect emission intensity changes, fluctuation trends, and stage characteristics, effectively enhancing the expressive power for complex and non-stationary emission processes. By normalizing the node adjacency relationship model and aggregating neighborhood features, the spatial features of each node not only contain its own observation information but also integrate the emission status of its neighboring nodes, improving the robustness and consistency of spatial features in complex monitoring scenarios.

[0043] In one alternative implementation, Based on the temporal evolution feature vector and aggregated spatial features, the emission anomaly score and anomaly type label of each node are determined to obtain the emission status assessment results, including: The temporal evolution feature vector is concatenated with the aggregated spatial features to obtain a fused feature vector. The fused feature vector is then transformed in terms of feature dimension and mapped through a fully connected layer to obtain an intermediate representation vector. The intermediate representation vector is then weighted by channel attention along the feature dimension and key feature components are extracted. The key feature components are then input into the anomaly score prediction branch to obtain continuous anomaly score values ​​for each node. The continuous anomaly score values ​​are then normalized and mapped to a preset score interval to obtain a standardized anomaly score. The intermediate representation vector is input into the multi-class prediction branch and nonlinearly transformed through a multi-layer fully connected network to obtain the unnormalized classification vector corresponding to each anomaly type. The unnormalized classification vector is then subjected to softmax normalization to obtain the anomaly type probability distribution. The category index with the highest probability is extracted from the anomaly type probability distribution as the anomaly type label, and the highest probability value is extracted as the category confidence. The category confidence and the standardized anomaly score are weighted and fused to obtain a comprehensive anomaly index. The standardized anomaly score is combined with the anomaly type label, and the comprehensive anomaly index is associated with the node identifier and timestamp information of each node to obtain the emission status assessment result.

[0044] The aggregated spatial features are globally averaged along the node dimension to obtain a 16-dimensional spatial feature vector, which is then concatenated with the temporal evolution feature vector along the feature dimension to obtain an 80-dimensional fused feature vector.

[0045] The edge computing unit performs feature dimension transformation on the fused feature vector and maps it through a fully connected layer to obtain an intermediate representation vector. The feature dimension transformation is achieved through linear mapping, using a weight matrix to map the 80-dimensional fused feature vector to a 120-dimensional intermediate representation. The fully connected layer contains weight parameters and bias parameters; the weight parameter matrix has a dimension of 80×120, and the bias parameter vector has a dimension of 120. The mapping process can be described as the product of the input vector and the weight matrix plus the bias vector. For example, if the first element of the fused feature vector is 0.75, the first element of the first row and first column of the weight matrix is ​​0.2, and the first element of the bias vector is 0.1, then the first element of the intermediate representation vector is calculated as 0.75 multiplied by 0.2 plus 0.1 and the contributions of other elements. After mapping, a ReLU activation function is applied to introduce a nonlinear transformation to obtain the intermediate representation vector.

[0046] The edge computing unit applies channel attention weights to the intermediate representation vector along its feature dimensions and extracts key feature components. The channel attention mechanism is implemented through a two-layer fully connected network. The first layer compresses the 120-dimensional intermediate representation vector to 30 dimensions, and the second layer expands the 30-dimensional features back to 120 dimensions, forming the channel attention weights. A ReLU activation function is applied between the two layers to enhance non-linear expressiveness. The generated channel attention weights are mapped to a value between 0 and 1 using a sigmoid function, representing the importance of each channel. For example, if the attention weight of a channel is 0.85, it indicates that the channel contains relatively important information. The channel attention weights are then multiplied element-wise with the intermediate representation vector to obtain the weighted feature vector.

[0047] The edge computing unit inputs the key feature components into the anomaly scoring prediction branch to obtain continuous anomaly scores for each node. The anomaly scoring prediction branch consists of a fully connected layer, taking 120-dimensional key feature components as input and outputting a 1-dimensional anomaly score. The weight matrix of the fully connected layer has a dimension of 120×1, and the bias parameter is a scalar. A linear transformation maps the high-dimensional features to a single anomaly score. For example, for the first distributed sensing node, if the anomaly score calculated from its key feature components is 0.78, it indicates that the node has anomaly characteristics.

[0048] The edge computing unit normalizes the continuous anomaly scores and maps them to a preset scoring interval to obtain a standardized anomaly score. The normalization process uses a maximum-minimum normalization method, mapping the anomaly scores of all nodes to a range of 0 to 1. The maximum and minimum values ​​of the anomaly scores are calculated; for example, if the maximum value is 0.95 and the minimum value is 0.15, the normalization operation subtracts 0.15 from the original score and divides by 0.8. The normalized score is then mapped to a preset scoring interval, such as 1 to 100, to obtain the standardized anomaly score. For instance, if the normalized score is 0.75, the standardized anomaly score mapped to the 1 to 100 interval is 75.

[0049] The edge computing unit inputs the intermediate representation vector into the multi-class prediction branch and performs a nonlinear transformation through a multi-layer fully connected network to obtain the unnormalized classification vector corresponding to each anomaly type. The multi-class prediction branch contains two fully connected network layers. The first layer transforms the 120-dimensional intermediate representation vector into a 64-dimensional hidden layer representation, and the second layer maps the 64-dimensional hidden layer representation into an unnormalized classification vector representing the number of anomaly types. Six emission anomaly types are defined, such as equipment leakage, emission exceeding standards, abnormal emission patterns, and sensor failure, resulting in an unnormalized classification vector with a dimension of 6. The ReLU activation function is applied between the two layers to enhance the nonlinear expressive power of the model.

[0050] The edge computing unit performs softmax normalization on the unnormalized classification vector to obtain the anomaly type probability distribution. Softmax normalization transforms the unnormalized classification vector into a probability distribution, ensuring that the sum of the probabilities of each category is 1. For example, for the first distributed sensing node, if its unnormalized classification vector is [2.5, 0.8, 1.3, 0.2, 0.5, 0.1], the probability distribution obtained after softmax normalization is approximately [0.65, 0.12, 0.19, 0.01, 0.02, 0.01], indicating that the probability of the current node belonging to the first type of anomaly is 0.65, and the probability of belonging to other types of anomalies is low.

[0051] The edge computing unit extracts the category index with the highest probability from the anomaly type probability distribution as the anomaly type label, and simultaneously extracts the highest probability value as the category confidence score. In the example above, the anomaly type label is 0, and the category confidence score is 0.65. When the maximum probability value is lower than a preset threshold of 0.5, the anomaly type label can be set to either "undetermined" or "normal."

[0052] The edge computing unit weights and fuses the category confidence score and the standardized anomaly score to obtain a comprehensive anomaly index. The weighted fusion uses a linear combination method; for example, the standardized anomaly score is multiplied by a weight of 0.7, and the category confidence score is multiplied by a weight of 0.3. The two are then added together to obtain the comprehensive anomaly index. If the standardized anomaly score is 75 and the category confidence score is 0.65, then the comprehensive anomaly index is 75 × 0.7 + 0.65 × 100 × 0.3 ≈ 71.95.

[0053] The edge computing unit combines the standardized anomaly score with the anomaly type label, and associates the comprehensive anomaly index with the node identifier and timestamp information of each node to obtain the emission status assessment result. The emission status assessment result includes information such as node identifier, timestamp, standardized anomaly score, anomaly type label, and comprehensive anomaly index. For example, for the distributed sensing node with node ID "Node-001", the assessment result at "2025-12-30 15:30:00" is: standardized anomaly score 75, anomaly type is equipment leakage, and comprehensive anomaly index is 71.95.

[0054] In this embodiment, by unifying the representation of temporal and spatial features and reconstructing the feature dimensions, it is possible to simultaneously perceive the dynamic changes and spatial correlation features of emission behavior. This avoids the problem of insufficient discrimination ability caused by the separate modeling of temporal and spatial information. By introducing a channel attention mechanism to weight the fused features, the feature components that are highly correlated with abnormal states are enhanced, thereby improving the sensitivity and stability of the anomaly score in response to real abnormal emissions. By setting an independent multi-classification prediction branch, the anomaly type is probabilistically modeled, which can not only identify whether anomalies exist, but also distinguish different anomaly forms, thereby improving the guiding value of emission anomaly analysis in actual regulatory and operational scenarios.

[0055] In one alternative implementation, Based on the measured concentration distribution and emission status assessment results in the emission feature dataset, candidate emission source locations are obtained through iterative search using a particle swarm optimization algorithm. Concentration time-series data within a preset distance of the candidate emission source locations are extracted and decomposed using wavelet decomposition to obtain high-frequency fluctuation components and low-frequency trend components, including: The measured concentration values ​​are extracted from the emission feature dataset and combined with the corresponding spatial coordinates to construct a concentration spatial distribution field. The set of nodes whose standardized anomaly scores exceed the preset anomaly threshold is extracted from the emission status assessment results, and the initial search center position is determined. The initial distribution of the particle swarm is initialized based on the initial search center position and the preset spatial constraints, and the position vector and velocity vector of each particle are set. Based on the current position of each particle, the concentration contribution value is determined in combination with the preset distance decay relationship and compared with the corresponding measured concentration values ​​in the concentration spatial distribution field to obtain the fitness value. The individual optimal position and global optimal position of the particle are updated according to the fitness value. The velocity update amount of the particle is determined in combination with the preset inertia weight. Boundary constraints are applied to the velocity update amount and iterative search is performed until the preset maximum number of iterations is reached. The global optimal position is used as the candidate position of the emission source. The concentration time series data is obtained by determining the set of neighboring nodes within a preset distance range centered on the candidate emission source location and extracting the concentration observation sequence within a continuous time window. The concentration time series data is then separated into an initial high-frequency component and an initial low-frequency component through multi-scale wavelet decomposition. The initial high-frequency component is subjected to threshold denoising processing to obtain a high-frequency fluctuation component, and the initial low-frequency component is subjected to trend fitting to obtain a low-frequency trend component.

[0056] Concentration values ​​of the target pollutant are extracted from the concentration data of various pollutants collected by distributed sensing nodes, such as the measured concentration value of sulfur dioxide, in ppb. Simultaneously, the spatial coordinate information of each node is acquired; for example, the spatial coordinates of node A are (120.5, 30.7, 5.2), representing longitude, latitude, and altitude information in degrees and meters, respectively. The spatial coordinates of the nodes are associated with the corresponding measured concentration values ​​to form a set of sampling points, such as {(120.5, 30.7, 5.2, 35.6), (120.6, 30.7, 4.8, 42.3), ...}, where the last value represents the pollutant concentration. A continuous spatial distribution field of concentration is generated within the monitoring area using Kriging interpolation, with a spatial resolution of 5 meters × 5 meters, covering the entire monitoring area.

[0057] The edge computing unit extracts a set of nodes with standardized anomaly scores exceeding a preset anomaly threshold from the emission status assessment results and determines the initial search center location. The preset anomaly threshold is set to 70. The emission status assessment results are traversed, and nodes with standardized anomaly scores greater than 70 are selected to form an anomaly node set. For example, in a certain assessment, three nodes with anomaly scores of 83, 75, and 72 are detected, constituting an anomaly node set. The initial search center location is determined by calculating the weighted average of the coordinates of each node in the anomaly node set, with the weight proportional to the node's standardized anomaly score. For the aforementioned three nodes, assuming their spatial coordinates are (120.5, 30.7, 5.2), (120.6, 30.7, 4.8), and (120.55, 30.75, 5.0), the calculated initial search center location is (120.54, 30.71, 5.05).

[0058] The edge computing unit initializes the initial distribution of the particle swarm based on the initial search center position and preset spatial constraints, and sets the position and velocity vectors of each particle. During the initialization phase, the particle swarm size is set to 50 particles. Within a spherical space with a radius of 100 meters, 50 initial positions of the particles are randomly generated, centered on the initial search center position. Each particle's position vector contains three components, representing longitude, latitude, and altitude. For example, the initial position vector of the first particle is (120.53, 30.72, 4.9), which has a certain offset from the initial search center position. The particle's velocity vector also contains three components, with initial values ​​randomly generated. The numerical range is controlled between [-0.01, 0.01] degrees / iteration and [-5, 5] meters / iteration, representing the amount of movement of the particle in the longitude, latitude, and altitude directions during each iteration update.

[0059] The edge computing unit determines the concentration contribution value based on the current position of each particle, combined with a preset distance decay relationship, and compares it with the corresponding measured concentration value in the concentration spatial distribution field to obtain the fitness value. The preset distance decay relationship adopts a Gaussian diffusion model, considering the emission source intensity and the distance decay factor. For each particle, assuming that the location is a potential emission source, its concentration contribution value to each monitoring node is calculated. For example, for a particle at position (120.53, 30.72, 4.9), assuming an emission source intensity of 100 ppb and a distance decay factor of 0.05, the calculated concentration contribution value of this particle to node H, i.e., node (120.5, 30.7, 5.2), is 32.4 ppb. The calculated concentration contribution value is compared with the measured concentration value at the corresponding location in the concentration spatial distribution field, and the root mean square error is calculated as the fitness value. The smaller the fitness value, the closer the particle position is to the actual emission source position. For the example above, if the measured concentration at node H is 35.6 ppb, then the prediction error for this particle at that node is 3.2 ppb. Taking into account the prediction errors of all monitored nodes, the calculated fitness value for this particle is 5.8.

[0060] The edge computing unit updates the individual optimal position and global optimal position of particles based on their fitness values. For each particle, if the fitness value in the current iteration is better than the particle's historical best fitness value, then the individual optimal position is updated to the current position. The position with the smallest fitness value among all the individual optimal positions of particles is then updated as the global optimal position. For example, if after the 10th iteration, the first particle's individual optimal position is (120.52, 30.73, 5.1) with a fitness value of 4.2, and its global optimal position is (120.51, 30.74, 5.0) with a fitness value of 3.8.

[0061] The edge computing unit determines the particle's velocity update amount by combining a preset inertia weight. The velocity update formula consists of three parts: an inertia term, a cognitive term, and a social term. The inertia term is related to the particle's current velocity, with an inertia weight set to 0.7; the cognitive term reflects the particle's use of its own experience, with a cognitive factor set to 1.5; and the social term represents the particle's learning from group experience, with a social factor set to 2.0. Taking the first particle as an example, if its current velocity is (0.005, -0.008, 2.0), its current position is (120.53, 30.72, 4.9), its individual optimal position is (120.52, 30.73, 5.1), and its global optimal position is (120.51, 30.74, 5.0), then the calculated velocity update amount is (0.0028, 0.0075, 1.4).

[0062] The edge computing unit applies boundary constraints to the velocity update and iteratively searches until a preset maximum number of iterations is reached, using the globally optimal location as a candidate emission source location. To prevent excessively high particle velocities from causing search instability, a velocity upper limit is set at ±0.02 degrees in longitude and latitude, and ±10 meters in altitude per iteration cycle. When the calculated velocity exceeds the upper limit, it is truncated to the boundary value. The particle position is updated based on the updated velocity, ensuring that the particle position does not exceed the preset search space range. The maximum number of iterations is set to 100, and the iteration stops when the maximum number of iterations is reached or the change in the globally optimal fitness value is less than 0.01 after 10 consecutive iterations. After the iteration ends, the globally optimal location (120.51, 30.74, 5.0) is used as a candidate emission source location, with a corresponding fitness value of 3.6.

[0063] The edge computing unit determines a set of neighboring nodes within a preset distance range centered on the candidate emission source location and extracts the concentration observation sequence within a continuous time window to obtain concentration time series data. The preset distance range is set to 200 meters, and all monitoring nodes within 200 meters of the candidate emission source location (120.51, 30.74, 5.0) are selected to form a set of neighboring nodes. For each node in the set of neighboring nodes, the concentration observation sequence for the most recent 24 hours is extracted, with a sampling interval of 5 minutes, resulting in concentration time series data containing 288 time points. For example, the concentration time series data for node H is [32.5, 33.8, 35.6, 34.2, ...].

[0064] The edge computing unit separates the concentration time-series data into initial high-frequency and initial low-frequency components through multi-scale wavelet decomposition. Discrete wavelet transform is employed, and the db4 wavelet basis is selected to decompose the concentration time-series data to a third level. Taking the concentration time-series data of node H as an example, the third-level approximation coefficients obtained through wavelet decomposition correspond to the initial low-frequency components, while the detailed coefficients correspond to the initial high-frequency components. The initial low-frequency components reflect the overall trend of concentration changes, while the initial high-frequency components contain rapid fluctuations and noise in the concentration data.

[0065] The edge computing unit performs threshold denoising on the initial high-frequency components to obtain high-frequency fluctuation components. A soft threshold denoising method is used, with the threshold value set to 1.5 times the standard deviation of the initial high-frequency components. For each wavelet coefficient, if the absolute value of the wavelet coefficient is less than the preset threshold, it is set to zero; if the absolute value is greater than the threshold, the absolute value is subtracted from the threshold, while maintaining the original sign. After denoising, the obtained high-frequency fluctuation components retain significant concentration fluctuation information and remove the influence of random noise.

[0066] Edge computing units perform trend fitting on the initial low-frequency components to obtain low-frequency trend components. A polynomial fitting method, with an order of 3, is used to fit the initial low-frequency components, resulting in smooth low-frequency trend components. These low-frequency trend components reflect the long-term variation trend of pollutant concentrations, such as diurnal variation patterns or variation patterns influenced by meteorological conditions.

[0067] In this embodiment, a continuous spatial distribution field of concentration is constructed using measured concentration values, and nodes with anomaly scores exceeding a threshold are used as initial search centers. This allows the emission source search process to focus on high-risk areas, improving source tracing efficiency and convergence reliability. By introducing a particle swarm optimization mechanism and comprehensively considering the distance decay relationship and the degree of matching with measured concentrations during the search process, misjudgments caused by local noise or occasional anomalies are effectively suppressed, improving the robustness of emission source location results under complex environmental conditions. By performing multi-scale wavelet decomposition on the concentration time series data of nodes surrounding candidate emission sources, short-term fluctuation components are separated from long-term trends, significantly reducing the impact of environmental noise and random disturbances on source tracing judgment and enhancing the ability to distinguish between stable and intermittent emission sources.

[0068] Figure 2 This is a flowchart illustrating the pollution source location analysis process of the online monitoring method for fugitive emissions based on edge computing, as described in this embodiment of the invention.

[0069] In one alternative implementation, The pulse characteristics of emissions are determined based on the high-frequency fluctuation components, the emission intensity variation patterns are identified based on the low-frequency trend components, and the theoretical diffusion distance of pollutants in each wind sector is calculated and the boundary of the influence range is determined by combining the candidate emission source locations and real-time meteorological parameters. Peak detection is performed on the high-frequency fluctuation component to extract the peak position and peak intensity. The time interval between the peak positions is counted to obtain the pulse frequency. The rising slope and falling slope corresponding to the peak intensity are extracted and combined to obtain the response rate. The pulse frequency and the response rate are combined to obtain the emission pulse characteristics. The low-frequency trend component is fitted with a polynomial to obtain a baseline curve, and the derivative is used to obtain a rate of change curve. The rate of change curve is transformed in the frequency domain and the dominant frequency component is extracted to obtain the dominant period. The deviation amplitude and deviation peak between the low-frequency trend component and the baseline curve are determined to obtain the fluctuation characteristics. The dominant period is compared with the pulse frequency to obtain the period matching degree. Based on the period matching degree and the fluctuation characteristics, the intermittent characteristics are determined. The emission intensity change law is obtained by combining the evolution trend of the baseline curve. Real-time meteorological data is acquired and the wind direction angle corresponding to the candidate emission source location is determined. Based on the wind direction angle, wind sector is divided and the average wind speed in each wind sector is calculated to obtain the sector wind speed. The atmospheric stability level is determined based on the real-time meteorological data and the diffusion coefficient is determined in combination with the sector wind speed. The emission source intensity correction coefficient is calculated based on the pulse characteristics and the emission intensity variation law. The theoretical diffusion distance of pollutants in each wind sector is calculated in combination with the diffusion coefficient and a fan-shaped region is generated. The boundary coordinate points of the fan-shaped region are connected to obtain the boundary of the influence range.

[0070] The edge computing unit performs peak detection on high-frequency fluctuation components to extract peak positions and intensities. Peak detection employs a sliding window method with a window width of 11 time points. Within each sliding window, if the value at the center point is greater than the values ​​of all other points within the window and exceeds a preset threshold of 5.0 ppb, it is marked as a peak point. For each detected peak point, the corresponding time position and peak intensity are recorded. Taking node H as an example, six peak points are detected in the high-frequency fluctuation components, with time positions of 07:15, 09:30, 11:45, 14:00, 16:15, and 18:30, and corresponding peak intensities of 8.6 ppb, 12.3 ppb, 9.5 ppb, 11.8 ppb, 7.2 ppb, and 10.5 ppb, respectively.

[0071] The edge computing unit calculates the pulse frequency by statistically analyzing the time intervals between peak positions. The time difference between adjacent peaks is calculated; for example, the time interval between the first and second peaks is 135 minutes, and the time interval between the second and third peaks is also 135 minutes. The average pulse period is obtained by averaging all time intervals. For node H, the average pulse period is 135 minutes, corresponding to a pulse frequency of 0.44 pulses per hour. Simultaneously, the standard deviation of the time intervals is calculated to assess pulse regularity. For node H, the standard deviation of the time intervals is 0 minutes, indicating a high degree of pulse regularity.

[0072] The edge computing unit extracts the rise and fall slopes corresponding to the peak intensities and combines them to obtain the response rate. For each peak, the rise slope from the previous trough to the peak and the fall slope from the peak to the next trough are calculated. The rise slope is expressed as the ratio of peak intensity to rise time, and the fall slope is expressed as the ratio of peak intensity to fall time. Taking the second peak as an example, the peak intensity is 12.3 ppb, the rise time is 45 minutes, and the fall time is 40 minutes, resulting in a calculated rise slope of 0.273 ppb / min and a fall slope of 0.308 ppb / min. The average rise and fall slopes of all peaks are then averaged to obtain an average rise slope of 0.26 ppb / min and an average fall slope of 0.29 ppb / min. The response rate is defined as a weighted average of the rise and fall slopes, with weights of 0.4 and 0.6, respectively, resulting in a calculated response rate of 0.278 ppb / min.

[0073] Edge computing units combine pulse frequency and response rate to obtain the pulse characteristics of emissions. The pulse characteristic index is calculated by multiplying the pulse frequency and response rate. For node H, the pulse characteristic index is 0.44 times / hour × 0.278 ppb / minute = 0.122 ppb / minute·times / hour. This index reflects the instantaneous response characteristics of the emission source; a higher value indicates stronger impulsivity and a faster response rate. Based on a preset threshold range, the pulse characteristic index is divided into three levels: low impulsivity, medium impulsivity, and high impulsivity. For node H, the pulse characteristic index falls into the medium impulsivity range and is therefore classified as having medium impulsivity.

[0074] The edge computing unit performs polynomial fitting on the low-frequency trend components to obtain the baseline curve and then differentiates it to obtain the rate of change curve. The polynomial fitting uses the least squares method with an order of 4, and the coefficients of the fitting equation are obtained through matrix operations. Taking the low-frequency trend components of node H as an example, the coefficients obtained through polynomial fitting are [0.000012, -0.00086, 0.019, -0.12, 36.5], representing the coefficients of the fourth, third, fourth, first, and constant terms, respectively. The baseline curve is generated based on the coefficients obtained from the polynomial fitting, covering the entire 24-hour time window. The rate of change curve is obtained by differentiating the baseline curve. The derivative is calculated using the central difference method, specifically by dividing the difference in function values ​​between two consecutive time points by the time interval. For example, at time point t=180 minutes, the function value of the baseline curve is 38.2 ppb, and at t=185 minutes, the function value is 38.3 ppb. With a time interval of 5 minutes, the rate of change at this point is calculated to be 0.02 ppb / minute.

[0075] The edge computing unit performs a frequency domain transformation on the rate of change curve and extracts the dominant frequency component to obtain the dominant period. The frequency domain transformation uses the Fast Fourier Transform algorithm to transform the rate of change curve from the time domain to the frequency domain. In the frequency domain representation, the frequency component with the largest amplitude is extracted as the dominant frequency component. The calculated dominant frequency is 0.0074 times / minute, and the corresponding dominant period is 135 minutes.

[0076] The edge computing unit determines the deviation magnitude and peak deviation between the low-frequency trend component and the baseline curve to obtain the fluctuation characteristics. The difference between the low-frequency trend component and the baseline curve at each time point is calculated to obtain the deviation sequence. The root mean square value of the deviation sequence is calculated to obtain the deviation magnitude; for node H, the deviation magnitude is 1.8 ppb. The maximum absolute value in the deviation sequence is extracted as the peak deviation; for node H, the peak deviation is 4.2 ppb, occurring at t=345 minutes. The fluctuation characteristic is defined as the weighted sum of the deviation magnitude and the peak deviation, with weights of 0.6 and 0.4, respectively; the calculated fluctuation characteristic value is 2.76 ppb.

[0077] The edge computing unit compares the dominant period with the pulse frequency to obtain the period matching degree, and determines the intermittent characteristic based on the period matching degree and fluctuation characteristics. The period matching degree is defined as the ratio of the dominant period to the pulse period. For node H, the dominant period is 135 minutes, and the pulse period is also 135 minutes, resulting in a calculated period matching degree of 1.0. When the period matching degree is between 0.9 and 1.1 and the fluctuation characteristic is greater than 2.0 ppb, it is determined to have significant intermittent characteristics. For node H, the period matching degree is 1.0, and the fluctuation characteristic is 2.76 ppb, meeting the determination criteria and thus determining that it has significant intermittent characteristics.

[0078] Edge computing units combine the evolution trend of the baseline curve to obtain the emission intensity variation pattern. The evolution trend of the baseline curve is determined by calculating the mean, maximum, minimum values ​​and their occurrence times over a 24-hour period. For node H, the mean of the baseline curve is 37.6 ppb, the maximum value of 41.8 ppb occurs at 14:30, and the minimum value of 34.2 ppb occurs at 03:15. Based on the aforementioned characteristics, a diurnal variation pattern of high emission intensity during the day and low emission intensity at night is identified, with the peak occurring in the afternoon. The emission intensity variation pattern is stored in text description form to provide a basis for subsequent emission source intensity estimation.

[0079] The edge computing unit acquires real-time meteorological data and determines the wind direction angle corresponding to the candidate emission source locations. Real-time meteorological data, including parameters such as wind speed, wind direction, temperature, humidity, and air pressure, is obtained from meteorological monitoring stations. For each candidate emission source location, wind direction data for the most recent hour is extracted, and the prevailing wind direction angle is calculated. Wind direction angles follow meteorological conventions: northerly winds are 0 degrees, increasing clockwise to 90 degrees for easterly winds, 180 degrees for southerly winds, and 270 degrees for westerly winds. For the candidate emission source location in the aforementioned example, the prevailing wind direction angle is 135 degrees, representing a southeasterly wind.

[0080] The wind sector is determined by dividing the wind direction into sectors based on the wind direction angle and calculating the average wind speed within each sector. A 360-degree wind direction circle is divided into 8 sectors, each covering a 45-degree angle range. For a prevailing wind direction of 135 degrees, its corresponding sector ranges from 112.5 degrees to 157.5 degrees. The average wind speed within this sector is calculated as the sector wind speed; for the emission source in the aforementioned example, the sector wind speed is 3.2 m / s.

[0081] The edge computing unit determines the atmospheric stability level based on real-time meteorological data and the diffusion coefficient by combining it with sector wind speed. The atmospheric stability level is determined comprehensively based on solar radiation intensity, cloud cover, and wind speed, using the Pasquill-Gifford stability classification method to divide atmospheric stability into six levels, A to F. For the meteorological conditions in this case, the atmospheric stability level is determined to be C, indicating slightly unstable atmosphere. Based on the atmospheric stability level and sector wind speed, the horizontal and vertical diffusion coefficients are determined by referring to tables or applying empirical formulas. For a stability level of C, the horizontal diffusion coefficient is 0.22, and the vertical diffusion coefficient is 0.16.

[0082] The edge computing unit calculates the emission source intensity correction coefficient based on pulse characteristics and emission intensity variation patterns. Combined with the diffusion coefficient, it calculates the theoretical diffusion distance of pollutants in each wind direction sector and generates sector-shaped regions. The emission source intensity correction coefficient considers pulse characteristic indicators and daily variation amplitude of emission intensity; for node H, the correction coefficient is calculated to be 1.25. A Gaussian diffusion model is applied to calculate the theoretical diffusion distance of pollutants in each wind direction sector, with the distance reaching a specified concentration threshold serving as the diffusion boundary. The specified concentration threshold is 0.5 times the ambient air quality standard limit; for sulfur dioxide, the threshold is 75 μg / m³. 3 For each wind sector, considering sector wind speed, atmospheric stability, and corrected source strength, the diffusion boundary distance is calculated. For the prevailing wind sector, the calculated diffusion distance is 850 meters.

[0083] Edge computing units connect the boundary coordinates of the sector regions to obtain the boundary of the influence range. Within each sector, a sector is generated with the candidate emission source location as the center and the calculated diffusion distance as the radius. The boundary coordinates of the sector regions of adjacent sectors are connected sequentially to form a closed polygon, which serves as the boundary of the pollutant's influence range. The sequence of influence range boundary coordinates is stored in a GIS-compatible format for easy visualization in the monitoring system.

[0084] In this embodiment, by extracting peak and slope features from high-frequency fluctuation components, the transient pulse behavior during the emission process can be quantitatively described, effectively distinguishing between continuous stable emissions and short-term sudden emissions, and significantly enhancing the ability to identify intermittent and pulsed emission characteristics. By performing baseline modeling and periodic feature analysis on low-frequency trend components, the long-term variation trend and periodicity of emission concentration can be clearly characterized. By matching periodic features with high-frequency pulse characteristics, the accuracy and stability of emission behavior discrimination are improved. By introducing real-time meteorological conditions, environmental corrections are made to the emission source behavior characteristics, enabling emission intensity assessment to dynamically reflect the impact of changes in wind direction, wind speed, and atmospheric stability on pollutant diffusion, overcoming the shortcomings of fixed diffusion parameters or empirical assumptions that are difficult to adapt to changes in actual operating conditions.

[0085] In one alternative implementation, Based on the emission intensity variation pattern and the boundary of the impact range, and combined with preset emission threshold constraints and duration determination criteria, a system is generated that includes emission early warning signals and source location results, including: The current emission intensity value and intensity change trend are extracted from the emission intensity change pattern. The current emission intensity value is compared with the preset emission threshold constraint to obtain the exceeding state. The duration of continuous exceeding is counted and the persistence judgment result is determined according to the preset duration judgment criterion. The duration is predicted based on the intensity change trend to obtain the expected duration. The warning level is determined according to the expected duration and the persistence judgment result and combined with the exceeding state to obtain the emission warning signal. The locations of monitoring points within the boundary of the influence range are extracted, the spatial distance between each monitoring point and the candidate emission source is calculated, and the wind sector to which each monitoring point belongs is determined. The measured concentration values ​​of each monitoring point are spatially interpolated with the current emission intensity value to obtain an interpolated concentration field. The gradient corresponding to the interpolated concentration field is calculated, and the direction of maximum gradient is extracted. The source coordinates are obtained by tracing back along the direction of maximum gradient to the center of the boundary of the influence range. The source coordinates are fused with the warning level to obtain the source location result. The emission warning signal and the source location result are encapsulated and output.

[0086] The edge computing unit extracts the current emission intensity value and its trend from the emission intensity variation pattern. The current emission intensity value is obtained by finding the function value corresponding to the current time point in the baseline curve. For example, the current time is 15:20, and the corresponding emission intensity value is 40.8 ppb. The intensity variation trend is determined by calculating the average slope of the baseline curve over the previous 30 minutes. This is done by dividing the difference between the function value at the current time point and the time point 30 minutes ago by the time interval; the resulting slope represents the rate of concentration change. For the time point 15:20 in the aforementioned example, the average slope over the previous 30 minutes is -0.03 ppb / min, indicating that the emission intensity is slowly decreasing.

[0087] The edge computing unit compares the current emission intensity value with preset emission threshold constraints to determine the exceedance status. The emission threshold constraints are set according to the pollutant type and environmental standards; for sulfur dioxide, the emission threshold is set to 35 ppb. The comparison method is to directly compare the current emission intensity value with the emission threshold. If the current value is greater than the threshold, it is determined to be an exceedance status; otherwise, it is a non-exceedance status. In the aforementioned example, the current emission intensity value of 40.8 ppb is greater than the threshold of 35 ppb, and is therefore determined to be an exceedance status. The exceedance margin is calculated as the difference between the current value and the threshold divided by the threshold, expressed as a percentage; the exceedance margin is 16.6%.

[0088] The edge computing unit counts the duration of continuous exceedances and determines the persistence result based on preset duration judgment criteria. The duration counting method traces historical data and calculates the cumulative duration from the earliest continuous exceedance time point to the current time point. In the example above, the continuous exceedance started at 14:05 and lasted until the current time 15:20, a duration of 75 minutes. The duration judgment criteria are set according to pollutant characteristics and environmental risk level. For sulfur dioxide, a continuous exceedance duration of less than 30 minutes is considered a short-term exceedance, 30 to 120 minutes is considered a moderate continuous exceedance, and more than 120 minutes is considered a long-term continuous exceedance. A duration of 75 minutes is considered a moderate continuous exceedance, and the persistence judgment result is determined to be "moderate continuous".

[0089] The edge computing unit predicts the duration based on the intensity change trend. The prediction method estimates the time required for the emission intensity to decrease below the threshold based on the current emission intensity value, the intensity change trend, and the emission threshold. A linear extrapolation algorithm is used, subtracting the emission threshold from the current emission intensity and then dividing by the absolute value of the intensity change trend to obtain the predicted duration. In the example above, the current emission intensity is 40.8 ppb, the emission threshold is 35 ppb, and the intensity change trend is -0.03 ppb / minute, resulting in a predicted duration of 193.3 minutes. Considering the uncertainty of the prediction, the result is rounded to 190 minutes.

[0090] The edge computing unit determines the warning level based on the expected duration and persistence determination result, and combines this with the exceedance status to obtain an emission warning signal. The warning level determination criteria are based on a combination of the persistence determination result and the expected duration, and are divided into four levels: blue, yellow, orange, and red, corresponding to minor, low, moderate, and severe risk levels, respectively. A short-term exceedance with an expected duration of less than 60 minutes is classified as a blue warning; a short-term exceedance with an expected duration of 60 minutes or more, or a moderately persistent exceedance with an expected duration of less than 120 minutes, is classified as a yellow warning; a moderately persistent exceedance with an expected duration of 120 minutes or more, or a long-term persistent exceedance with an expected duration of less than 180 minutes, is classified as an orange warning; and a long-term persistent exceedance with an expected duration of 180 minutes or more is classified as a red warning. In the aforementioned example, the persistence determination result is "moderately persistent," the expected duration is 190 minutes, and the warning level is determined to be orange. The emission warning signal consists of the warning level, the state of exceeding the standard, the magnitude of exceeding the standard, the duration, and the expected duration. It is formatted into a uniform data structure. In the example, the emission warning signal is {"Warning Level": "Orange", "State of Exceeding Standard": true, "Magnitude of Exceeding Standard": "16.6%", "Duration": "75 minutes", "Expected Duration": "190 minutes"}.

[0091] The edge computing unit extracts the locations of monitoring points within the boundary of the impact area. It retrieves all monitoring points located within the previously calculated boundary of the impact area from the monitoring network database, obtaining their geographic coordinates and the latest monitoring data. In this case, there are five monitoring points within the boundary of the impact area, numbered M1 to M5. The relative locations of these monitoring points are distributed to the west, south, southeast, north, and south of the emission source. The latest monitored sulfur dioxide concentrations are 38.5 ppb, 42.3 ppb, 35.2 ppb, 36.8 ppb, and 39.1 ppb, respectively.

[0092] The edge computing unit calculates the spatial distance between each monitoring point and the candidate emission source location and determines the wind sector to which each monitoring point belongs. The spatial distance calculation uses the spherical distance formula, taking into account the Earth's curvature, and calculates the actual distance based on the coordinates of the monitoring point and the candidate emission source location. The wind sector determination method calculates the azimuth angle from the candidate emission source location to the monitoring point and determines the sector to which that azimuth angle belongs. In this case, the calculated distances between monitoring points M1 to M5 and the emission source are 1.57 km, 1.12 km, 1.83 km, 2.09 km, and 2.24 km, respectively, and their corresponding wind sectors are the west sector, south sector, southeast sector, north sector, and south sector.

[0093] The edge computing unit spatially interpolates the measured concentration values ​​at each monitoring point with the current emission intensity value to obtain an interpolated concentration field. Spatial interpolation employs an inverse distance weighted interpolation algorithm, calculating the concentration values ​​at interpolation points on a two-dimensional spatial grid. The grid resolution is set to 50 meters × 50 meters, covering the entire boundary of the influence area. For each point in the grid, a weighting coefficient is calculated based on the distance between each point and all monitoring points; the closer the distance, the greater the weight. The measured concentrations at all monitoring points are then weighted and averaged to obtain the concentration value at the interpolated point. The weight calculation uses an inverse square ratio of distance to ensure that nearby monitoring points have a greater impact. The interpolation calculation considers the current emission intensity value as a constraint, fixing the concentration values ​​at candidate emission source locations to the current emission intensity value. In the aforementioned example, the generated interpolated concentration field covers an area of ​​approximately 4 square kilometers, containing 1600 grid points, with concentration values ​​ranging from 34.2 ppb to 43.5 ppb.

[0094] The edge computing unit calculates the gradient corresponding to the interpolated concentration field and extracts the direction of maximum gradient. The gradient calculation uses the finite difference method; for each grid point, the concentration change rate in the x and y directions is calculated and combined to obtain the gradient vector. The magnitude of the gradient vector represents the severity of the concentration change, and the direction represents the direction of the fastest concentration increase. The region with the largest gradient magnitude is searched throughout the entire interpolated concentration field, and the average gradient direction of this region is extracted as the direction of maximum gradient. In the aforementioned example, the region with the maximum gradient is located near monitoring point M2, with an average gradient magnitude of 0.15 ppb / m and a maximum gradient direction in the northeast direction with an azimuth angle of 42°.

[0095] The edge computing unit performs reverse tracing along the direction of maximum gradient to the center of the influence range boundary to obtain the source coordinates. The reverse tracing algorithm is based on gradient ascent, starting from the region of maximum gradient and stepping in the opposite direction of the gradient to find local maxima in the concentration field. The step distance is adaptively adjusted according to the gradient magnitude; a smaller step size for a large gradient and a larger step size for a small gradient ensures search accuracy and efficiency. Path points are recorded during the search until the influence range boundary is reached or a local maximum concentration is detected. The search results are smoothed to remove anomalous jumps, and the final source coordinates are determined. In the aforementioned example, the reverse tracing path contains 27 path points, and the source coordinates are located at the center of the influence range boundary, with a deviation of less than 50 meters from the actual emission source coordinates.

[0096] The edge computing unit fuses the source coordinates and warning level to obtain the source location result, encapsulates the emission warning signal and source location result, and outputs them. The source location result includes source coordinates, location accuracy assessment, and warning level, formatted into a unified data structure. The location accuracy assessment is calculated based on the convergence of the reverse tracing path and the gradient characteristics of the interpolated concentration field, and is divided into three levels: high, medium, and low. In this case, the source location result is {"Source Coordinates": "Northeast corner of the storage tank area in the north area of ​​the factory", "Location Accuracy": "High", "Warning Level": "Orange"}. The edge computing unit encapsulates the emission warning signal and source location result into a complete monitoring report, adds a timestamp and monitoring node identifier, and outputs it to the monitoring center and mobile terminals through a preset communication interface. The data encapsulation uses a lightweight JSON format to ensure transmission efficiency and cross-platform compatibility. The monitoring report is also stored in the local database of the edge computing unit for subsequent querying and analysis.

[0097] In this embodiment, by comparing the current emission intensity with the threshold constraint and introducing statistics and predictions on the duration of continuous exceedances, the persistence of the exceedance state and its future development trend can be determined. This effectively reduces the probability of false alarms caused by short-term fluctuations, improves the stability and practical value of the early warning results, and uses the emission intensity change trend to predict the expected duration and determine the graded early warning level accordingly. This results in higher risk identification accuracy under complex operating conditions. By constructing a spatial interpolation concentration field within the influence range and extracting the direction of maximum gradient for reverse tracing, the emission source tracing process can comprehensively utilize the spatial information of multiple monitoring points, significantly improving the accuracy and consistency of source tracing and positioning under complex diffusion conditions.

[0098] A second aspect of this invention provides an online monitoring system for fugitive emissions based on edge computing, comprising: The data acquisition unit is used to acquire emission characteristic data from multiple distributed sensing nodes and combine the spatial coordinates of the distributed sensing nodes to establish the adjacency relationship and propagation connectivity between the nodes to obtain the emission characteristic dataset. The feature extraction unit is used to encode the emission feature dataset in the time dimension to obtain a temporal evolution feature vector, perform graph convolution on the emission feature dataset in the spatial dimension and aggregate the feature information of adjacent nodes based on the adjacency relationship to obtain aggregated spatial features, and determine the emission anomaly score and anomaly type label of each node based on the temporal evolution feature vector and aggregated spatial features to obtain the emission status assessment result. The source location unit is used to obtain candidate emission source locations through iterative search using a particle swarm optimization algorithm based on the measured concentration distribution and emission status assessment results in the emission feature dataset, extract the concentration time series data within a preset distance of the candidate emission source locations and perform wavelet decomposition to obtain high-frequency fluctuation components and low-frequency trend components. The characteristic analysis unit is used to determine the pulse characteristics of emissions based on the high-frequency fluctuation component, identify the emission intensity change pattern based on the low-frequency trend component, and calculate the theoretical diffusion distance of pollutants in each wind sector and determine the boundary of the influence range by combining the candidate emission source locations and real-time meteorological parameters. The result output unit is used to generate emission warning signals and source location results based on the emission intensity change pattern and the boundary of the influence range, combined with preset emission threshold constraints and duration determination criteria.

[0099] A third aspect of the present invention provides an electronic device, comprising: A processor and a memory for storing processor-executable instructions, wherein the processor is configured to invoke instructions stored in the memory to perform the aforementioned method.

[0100] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0101] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0102] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for online monitoring of fugitive emissions based on edge computing, characterized in that, include: Emission feature data from multiple distributed sensing nodes are acquired, and the adjacency relationships and propagation connectivity between the nodes are established by combining the spatial coordinates of the distributed sensing nodes to obtain the emission feature dataset. The emission feature dataset is encoded in the time dimension to obtain a temporal evolution feature vector. In the spatial dimension, the emission feature dataset is subjected to graph convolution and the feature information of adjacent nodes is aggregated based on the adjacency relationship to obtain aggregated spatial features. Based on the temporal evolution feature vector and aggregated spatial features, the emission anomaly score and anomaly type label of each node are determined to obtain the emission status assessment result. Based on the measured concentration distribution and emission status assessment results in the emission feature dataset, the candidate locations of emission sources are obtained by iterative search using the particle swarm optimization algorithm. The concentration time series data within a preset distance of the candidate locations of emission sources are extracted and wavelet decomposition is performed to obtain high-frequency fluctuation components and low-frequency trend components. The pulse characteristics of emissions are determined based on the high-frequency fluctuation components, the emission intensity variation pattern is identified based on the low-frequency trend components, and the theoretical diffusion distance of pollutants in each wind sector is calculated and the boundary of the influence range is determined by combining the candidate emission source locations and real-time meteorological parameters. Based on the emission intensity variation pattern and the influence range boundary, and combined with preset emission threshold constraints and duration determination criteria, an emission early warning signal and source location result are generated.

2. The method according to claim 1, characterized in that, The emission feature dataset is obtained by acquiring emission feature data from multiple distributed sensing nodes and establishing adjacency relationships and propagation connectivity between the nodes using their spatial coordinates. The dataset includes: The emission characteristic data reported by each distributed sensing node is obtained and the timestamp information and spatial coordinate information are extracted. Based on the timestamp information, the emission characteristic data of different distributed sensing nodes are aligned with the time reference to obtain time-aligned data. Based on the spatial coordinate information, the Euclidean distance between each distributed sensing node is calculated and neighboring node pairs are filtered according to a preset distance threshold. Connecting edges are established between neighboring node pairs that meet the distance conditions and an undirected graph is constructed to obtain the adjacency relationship. Acquire the prevailing wind direction data of the monitoring area and calculate the azimuth difference between each pair of neighboring nodes. Use the azimuth difference and the prevailing wind direction data to determine the neighboring node pairs in the downwind propagation link and mark the propagation connectivity. Combine the time alignment data and the adjacency relationship to obtain the emission characteristic dataset through correlation integration.

3. The method according to claim 1, characterized in that, Encoding the emission feature dataset in the time dimension yields a temporal evolution feature vector. Then, graph convolution is performed on the emission feature dataset in the spatial dimension, and feature information from adjacent nodes is aggregated based on the adjacency relationships to obtain aggregated spatial features, including: The observation values ​​of each distributed sensing node at multiple moments within a continuous time window are extracted from the emission feature dataset to form node time series data. The node time series data is sinusoidally and cosinely encoded to generate time position vectors and concatenated to obtain a time-labeled sequence. The query vector, key vector, and value vector corresponding to the time-labeled sequence are constructed along the time dimension. The dot product of the query vector and the key vector is calculated to obtain the time attention weight. The time attention weight and the value vector are weighted and summed, and residual connections and layer normalization are performed to obtain the intermediate result of time encoding. The deep semantic representation corresponding to the intermediate result of time encoding is extracted through a preset feedforward neural network and global average pooling is performed along the time dimension to obtain the time series evolution feature vector. An adjacency matrix is ​​constructed based on the adjacency relationship. The degree matrix corresponding to the adjacency matrix is ​​calculated and symmetric normalized to obtain a normalized adjacency matrix. The current feature values ​​of each node in the emission feature dataset are used to construct a node feature matrix. The node feature matrix is ​​linearly transformed with a preset initial weight matrix to obtain a transformed feature matrix. The normalized adjacency matrix and the transformed feature matrix are multiplied to obtain neighborhood aggregation features. A nonlinear activation function is applied to the neighborhood aggregation features and layer normalization is performed to obtain aggregation space features.

4. The method according to claim 1, characterized in that, Based on the temporal evolution feature vector and aggregated spatial features, the emission anomaly score and anomaly type label of each node are determined to obtain the emission status assessment results, including: The temporal evolution feature vector is concatenated with the aggregated spatial features to obtain a fused feature vector. The fused feature vector is then transformed in terms of feature dimension and mapped through a fully connected layer to obtain an intermediate representation vector. The intermediate representation vector is then weighted by channel attention along the feature dimension and key feature components are extracted. The key feature components are then input into the anomaly score prediction branch to obtain continuous anomaly score values ​​for each node. The continuous anomaly score values ​​are then normalized and mapped to a preset score interval to obtain a standardized anomaly score. The intermediate representation vector is input into the multi-class prediction branch and nonlinearly transformed through a multi-layer fully connected network to obtain the unnormalized classification vector corresponding to each anomaly type. The unnormalized classification vector is then subjected to softmax normalization to obtain the anomaly type probability distribution. The category index with the highest probability is extracted from the anomaly type probability distribution as the anomaly type label, and the highest probability value is extracted as the category confidence. The category confidence and the standardized anomaly score are weighted and fused to obtain a comprehensive anomaly index. The standardized anomaly score is combined with the anomaly type label, and the comprehensive anomaly index is associated with the node identifier and timestamp information of each node to obtain the emission status assessment result.

5. The method according to claim 1, characterized in that, Based on the measured concentration distribution and emission status assessment results in the emission feature dataset, candidate emission source locations are obtained through iterative search using a particle swarm optimization algorithm. Concentration time-series data within a preset distance of the candidate emission source locations are extracted and decomposed using wavelet decomposition to obtain high-frequency fluctuation components and low-frequency trend components, including: The measured concentration values ​​are extracted from the emission feature dataset and combined with the corresponding spatial coordinates to construct a concentration spatial distribution field. The set of nodes whose standardized anomaly scores exceed the preset anomaly threshold is extracted from the emission status assessment results, and the initial search center position is determined. The initial distribution of the particle swarm is initialized based on the initial search center position and the preset spatial constraints, and the position vector and velocity vector of each particle are set. Based on the current position of each particle, the concentration contribution value is determined in combination with the preset distance decay relationship and compared with the corresponding measured concentration values ​​in the concentration spatial distribution field to obtain the fitness value. The individual optimal position and global optimal position of the particle are updated according to the fitness value. The velocity update amount of the particle is determined in combination with the preset inertia weight. Boundary constraints are applied to the velocity update amount and iterative search is performed until the preset maximum number of iterations is reached. The global optimal position is used as the candidate position of the emission source. The concentration time series data is obtained by determining the set of neighboring nodes within a preset distance range centered on the candidate emission source location and extracting the concentration observation sequence within a continuous time window. The concentration time series data is then separated into an initial high-frequency component and an initial low-frequency component through multi-scale wavelet decomposition. The initial high-frequency component is subjected to threshold denoising processing to obtain a high-frequency fluctuation component, and the initial low-frequency component is subjected to trend fitting to obtain a low-frequency trend component.

6. The method according to claim 1, characterized in that, The pulse characteristics of emissions are determined based on the high-frequency fluctuation components, the emission intensity variation patterns are identified based on the low-frequency trend components, and the theoretical diffusion distance of pollutants in each wind sector is calculated and the boundary of the influence range is determined by combining the candidate emission source locations and real-time meteorological parameters. Peak detection is performed on the high-frequency fluctuation component to extract the peak position and peak intensity. The time interval between the peak positions is counted to obtain the pulse frequency. The rising slope and falling slope corresponding to the peak intensity are extracted and combined to obtain the response rate. The pulse frequency and the response rate are combined to obtain the emission pulse characteristics. The low-frequency trend component is fitted with a polynomial to obtain a baseline curve, and the derivative is used to obtain a rate of change curve. The rate of change curve is transformed in the frequency domain and the dominant frequency component is extracted to obtain the dominant period. The deviation amplitude and deviation peak between the low-frequency trend component and the baseline curve are determined to obtain the fluctuation characteristics. The dominant period is compared with the pulse frequency to obtain the period matching degree. Based on the period matching degree and the fluctuation characteristics, the intermittent characteristics are determined. The emission intensity change law is obtained by combining the evolution trend of the baseline curve. Real-time meteorological data is acquired and the wind direction angle corresponding to the candidate emission source location is determined. Based on the wind direction angle, wind sector is divided and the average wind speed in each wind sector is calculated to obtain the sector wind speed. The atmospheric stability level is determined based on the real-time meteorological data and the diffusion coefficient is determined in combination with the sector wind speed. The emission source intensity correction coefficient is calculated based on the pulse characteristics and the emission intensity variation law. The theoretical diffusion distance of pollutants in each wind sector is calculated in combination with the diffusion coefficient and a fan-shaped region is generated. The boundary coordinate points of the fan-shaped region are connected to obtain the boundary of the influence range.

7. The method according to claim 1, characterized in that, Based on the emission intensity variation pattern and the boundary of the impact range, and combined with preset emission threshold constraints and duration determination criteria, a system is generated that includes emission early warning signals and source location results, including: The current emission intensity value and intensity change trend are extracted from the emission intensity change pattern. The current emission intensity value is compared with the preset emission threshold constraint to obtain the exceeding state. The duration of continuous exceeding is counted and the persistence judgment result is determined according to the preset duration judgment criterion. The duration is predicted based on the intensity change trend to obtain the expected duration. The warning level is determined according to the expected duration and the persistence judgment result and combined with the exceeding state to obtain the emission warning signal. The locations of monitoring points within the boundary of the influence range are extracted, the spatial distance between each monitoring point and the candidate emission source is calculated, and the wind sector to which each monitoring point belongs is determined. The measured concentration values ​​of each monitoring point are spatially interpolated with the current emission intensity value to obtain an interpolated concentration field. The gradient corresponding to the interpolated concentration field is calculated, and the direction of maximum gradient is extracted. The source coordinates are obtained by tracing back along the direction of maximum gradient to the center of the boundary of the influence range. The source coordinates are fused with the warning level to obtain the source location result. The emission warning signal and the source location result are encapsulated and output.

8. An online monitoring system for fugitive emissions based on edge computing, used to implement the method of any one of claims 1-7, characterized in that, include: The data acquisition unit is used to acquire emission characteristic data from multiple distributed sensing nodes and combine the spatial coordinates of the distributed sensing nodes to establish the adjacency relationship and propagation connectivity between the nodes to obtain the emission characteristic dataset. The feature extraction unit is used to encode the emission feature dataset in the time dimension to obtain a temporal evolution feature vector, perform graph convolution on the emission feature dataset in the spatial dimension and aggregate the feature information of adjacent nodes based on the adjacency relationship to obtain aggregated spatial features, and determine the emission anomaly score and anomaly type label of each node based on the temporal evolution feature vector and aggregated spatial features to obtain the emission status assessment result. The source location unit is used to obtain candidate emission source locations through iterative search using a particle swarm optimization algorithm based on the measured concentration distribution and emission status assessment results in the emission feature dataset, extract the concentration time series data within a preset distance of the candidate emission source locations and perform wavelet decomposition to obtain high-frequency fluctuation components and low-frequency trend components. The characteristic analysis unit is used to determine the pulse characteristics of emissions based on the high-frequency fluctuation component, identify the emission intensity change pattern based on the low-frequency trend component, and calculate the theoretical diffusion distance of pollutants in each wind sector and determine the boundary of the influence range by combining the candidate emission source locations and real-time meteorological parameters. The result output unit is used to generate emission warning signals and source location results based on the emission intensity change pattern and the boundary of the influence range, combined with preset emission threshold constraints and duration determination criteria.

9. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 7.

Images