A multi-sensor cooperative air pollution monitoring method and system

Abstract

The application relates to the technical field of environmental monitoring and discloses a multi-sensor cooperative air pollution monitoring method and system. The method comprises the following steps: fusing fixed base station and mobile device monitoring data, correcting the data in time and space to obtain a fused data set; correcting the base station reference concentration through filtering and dynamically adjusting and updating the frequency, mapping a grid to obtain a reference concentration distribution; generating a prediction field through three-dimensional space-time mapping, optimizing when the accuracy threshold is exceeded; extracting diffusion trends to construct a refined concentration field, correcting the concentration field to be smooth and continuous through boundary correction; when there is a migration trend, dynamically adjusting the grid to generate a dynamic concentration prediction field, and iteratively correcting model parameters to obtain a high-precision local concentration field; combining real-time sensor observation values, high-frequency resampling and filter reconstruction are performed when the mutation threshold is exceeded, a closed-loop correction factor is generated, and a prediction model is iteratively optimized. The method can realize multi-source data fusion and dynamic accurate prediction of air pollution monitoring, and improves the monitoring accuracy and pollution migration tracking capability.

Description

Technical Field

[0001] This invention relates to the field of environmental monitoring technology, and in particular to a multi-sensor collaborative method and system for air pollution monitoring. Background Technology

[0002] Against the backdrop of increasingly stringent requirements for refined ecological and environmental protection and public health safeguards, air pollution monitoring is a core means of understanding air quality and formulating scientific prevention and control strategies. The widespread adoption of smart sensors provides technical support for the collection of multi-source pollution data, and their monitoring accuracy and data fusion capabilities directly determine the effectiveness of pollution prevention and control decisions.

[0003] Currently, air pollution monitoring combines fixed base stations and mobile monitoring devices, both equipped with smart sensors for automated sampling. However, on the one hand, while fixed base station monitoring is stable, it struggles to cover small urban areas and cannot capture localized, sudden pollution events; mobile devices offer high spatiotemporal resolution, but the data is fragmented and transient, and sensor accuracy deviations lead to uneven data quality. On the other hand, existing methods lack a professional fusion mechanism for multi-source smart sensor data, failing to address the inconsistency in spatiotemporal scales between fixed and mobile sampling, and only perform basic error correction on sensor data, unable to achieve refined analysis by incorporating pollutant diffusion characteristics. This type of monitoring model struggles to accurately depict the dynamic spatiotemporal distribution of pollutants, lacks the ability to capture sudden pollution events and pollutant migration processes, and results in significant discrepancies between monitoring results and actual pollution conditions. Summary of the Invention

[0004] This invention provides a multi-sensor collaborative air pollution monitoring method and system to integrate scattered monitoring data such as fixed base station benchmark concentration and mobile device sampling data into a unified concentration field characterization basis, and to carry out intelligent and high-precision control throughout the entire process.

[0005] In a first aspect, to address the aforementioned technical problems, the present invention provides a multi-sensor collaborative air pollution monitoring method, comprising: The base station baseline concentration is collected by a fixed base station, and the mobile sampling concentration is obtained by a mobile device to form a raw dataset. The data is then aligned by spatiotemporal inconsistency correction to obtain a fused dataset. Based on the fused dataset, Kalman filtering is used to predict and correct the base station baseline concentration. The measurement update frequency is adjusted according to the correction result to obtain an optimized concentration sequence, which is then mapped to a grid to obtain the baseline concentration distribution. The baseline concentration distribution is mapped to three-dimensional spatiotemporal space to obtain a baseline concentration field, which is then input into a preset prediction model to generate a prediction field. The concentration error within the evaluation range is evaluated. If the evaluation result exceeds a preset accuracy threshold, the original dataset is supplemented by historical data smoothing and Kriging interpolation to obtain an optimized concentration distribution. The concentration diffusion trend is extracted from the optimized concentration distribution to generate a baseline concentration weight distribution. Sampled concentration interpolation data is obtained, and a fine concentration distribution field is constructed by applying a grid resolution adjustment method in combination with the baseline concentration weight distribution. Boundary correction is performed on the concentration distribution field, discontinuous boundaries are identified and numerical assimilation is performed to obtain a smooth continuous concentration field; If the smooth continuous concentration field shows a migration trend, the spatiotemporal changes of the migration path are tracked and a grid deformation control vector is generated to dynamically adjust the spatiotemporal grid of the prediction field, thereby obtaining a dynamic concentration prediction field. Based on the dynamic concentration prediction field, a weighted baseline concentration distribution is calculated. A gradient deviation field is constructed based on the weighted baseline concentration distribution and the sampled concentration interpolation data. Based on the gradient deviation field, the parameters of the preset prediction model are iteratively corrected to obtain a high-precision local concentration field.

[0006] Secondly, the present invention provides a multi-sensor collaborative air pollution monitoring system, comprising: The data acquisition and fusion module is used to collect base station baseline concentrations through fixed base stations and obtain mobile sampling concentrations through mobile devices to form a raw dataset. The data is then aligned using spatiotemporal inconsistency correction to obtain a fused dataset. The Kalman filter correction module is used to predict and correct the base station reference concentration based on the fused dataset using Kalman filtering, adjust the measurement update frequency according to the correction result, obtain an optimized concentration sequence, and map it to a grid to obtain the reference concentration distribution. The concentration field prediction and optimization module is used to map the baseline concentration distribution to three-dimensional spatiotemporal space to obtain the baseline concentration field and input it into the preset prediction model to generate the prediction field. It evaluates the concentration error within the range. If the evaluation result exceeds the preset accuracy threshold, it supplements the original dataset by smoothing historical data and Kriging interpolation to obtain the optimized concentration distribution. The refined construction module is used to extract the concentration diffusion trend from the optimized concentration distribution, generate a baseline concentration weight distribution, obtain sampled concentration interpolation data, and combine the baseline concentration weight distribution with a grid resolution adjustment method to construct a refined concentration distribution field. The boundary correction processing module is used to correct the boundary of the concentration distribution field, identify discontinuous boundaries and perform numerical assimilation processing to obtain a smooth continuous concentration field. The dynamic migration prediction module is used to track the spatiotemporal changes of the migration path and generate a grid deformation control vector if the smooth continuous concentration field shows a migration trend, and dynamically adjust the spatiotemporal grid of the prediction field to obtain a dynamic concentration prediction field. The final output module is used to calculate the weighted baseline concentration distribution based on the dynamic concentration prediction field, construct a gradient deviation field based on the weighted baseline concentration distribution and the sampled concentration interpolation data, and iteratively correct the parameters of the preset prediction model based on the gradient deviation field to obtain a high-precision local concentration field.

[0007] Compared with the prior art, the present invention has the following beneficial effects: (1) This invention uses a fusion strategy of co-sampling between fixed base stations and mobile devices and spatiotemporal inconsistency correction to align the static base station reference concentration and dynamic mobile sampling concentration to a unified spatiotemporal framework, effectively eliminating the problem of spatiotemporal scale misalignment of multi-source monitoring data.

[0008] (2) By constructing a full-process monitoring system, the present invention first optimizes the background concentration distribution through Kalman filtering, then captures the pollutant migration characteristics and dynamically adjusts the grid by combining optical flow method and Lagrange particle tracking method, and finally corrects the model parameters through gradient deviation field iteration, which greatly improves the spatiotemporal characterization accuracy and dynamic prediction capability of the pollution concentration field.

[0009] (3) This invention uses a linkage mechanism of dynamic error assessment and Kalman filter closed-loop feedback to capture changes in concentration prediction error in real time and trigger high-frequency resampling, reconstruct the filter observation vector and generate a closed-loop correction factor to continuously optimize the model, quickly respond to sudden pollution and continuously improve prediction accuracy. Attached Figure Description

[0010] Figure 1 This is a schematic diagram of a multi-sensor collaborative air pollution monitoring method provided in the first embodiment of the present invention; Figure 2 This is a schematic diagram of a multi-sensor collaborative air pollution monitoring system provided in the second embodiment of the present invention. Detailed Implementation

[0011] 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.

[0012] Reference Figure 1 The first embodiment of the present invention provides a multi-sensor collaborative air pollution monitoring method, including steps S11 to S17: S11: Base station reference concentration is collected by fixed base station and mobile sampling concentration is obtained by mobile device to form raw dataset. Spatiotemporal inconsistency correction is used to align the data to obtain fused dataset. S12, Based on the fused dataset, Kalman filtering is used to predict and correct the base station reference concentration. The measurement update frequency is adjusted according to the correction result to obtain an optimized concentration sequence and map it to a grid to obtain the reference concentration distribution. S13, the baseline concentration distribution is mapped to three-dimensional spatiotemporal space to obtain the baseline concentration field and input into the preset prediction model to generate the prediction field. The concentration error within the evaluation range is evaluated. If the evaluation result exceeds the preset accuracy threshold, the original dataset is supplemented by historical data smoothing and Kriging interpolation to obtain the optimized concentration distribution. S14, extract the concentration diffusion trend from the optimized concentration distribution, generate a baseline concentration weight distribution, obtain sampled concentration interpolation data, and combine the baseline concentration weight distribution with a grid resolution adjustment method to construct a refined concentration distribution field; S15, perform boundary correction on the concentration distribution field, identify discontinuous boundaries and perform numerical assimilation processing to obtain a smooth continuous concentration field; S16, If the smooth continuous concentration field shows a migration trend, then track the spatiotemporal changes of the migration path and generate a grid deformation control vector to dynamically adjust the spatiotemporal grid of the prediction field to obtain a dynamic concentration prediction field. S17. Based on the dynamic concentration prediction field, a weighted baseline concentration distribution is calculated. A gradient deviation field is constructed based on the weighted baseline concentration distribution and the sampled concentration interpolation data. Based on the gradient deviation field, the parameters of the preset prediction model are iteratively corrected to obtain a high-precision local concentration field.

[0013] In step S11, a base station reference concentration is collected using a fixed base station, and a mobile sampling concentration is obtained using a mobile device to form an original dataset. Spatiotemporal inconsistency correction is then used to align the data, resulting in a fused dataset, including: Obtain the base station baseline concentration collected by fixed base stations and the mobile sampling concentration collected by mobile devices to construct the original dataset; The original dataset is parsed to calculate the spatiotemporal difference feature values ​​between the base station baseline concentration and the mobile sampling concentration; If the spatiotemporal difference feature value has a nonlinear deviation, then the corresponding spatiotemporal inconsistency distribution vector is extracted. The spatiotemporal inconsistency distribution vector includes the spatial deviation direction, the degree of time lag, and the magnitude of concentration change. A correction mapping matrix is ​​constructed based on the spatiotemporal inconsistency distribution vector, and alignment processing is performed on the original dataset to obtain the fused dataset.

[0014] First, it should be noted that the data sources used for data collection include fixed urban air quality monitoring base stations (hereinafter referred to as "fixed base stations") and mobile monitoring devices equipped with multiple sensors (hereinafter referred to as "mobile devices"). Fixed base stations are deployed at national / provincial monitoring sites in the city and are capable of detecting six parameters: PM2.5, PM10, SO2, NO2, CO, and O3, with a data collection frequency of once per minute. Mobile devices are equipped with miniaturized optical and electrochemical sensors and are deployed on mobile carriers such as sanitation vehicles and buses, collecting concentration data of the same six parameters at a frequency of once every 30 seconds, and simultaneously recording latitude, longitude, and timestamps via GPS modules.

[0015] The data acquisition process is initiated, and the base station baseline concentration and mobile sampling concentration within the same monitoring period are acquired synchronously. These are then summarized and organized according to the dimensions of "monitored pollutant type - acquisition timestamp - device identifier" to form the original dataset. For example, the PM2.5 base station baseline concentration collected by the fixed base station numbered GD-08 is [28, 30, 29, 32, 31, ...] μg / m³, with timestamps accurate to the minute, such as 2026-03-01 08:01:00; the PM2.5 mobile sampling concentration collected by the mobile device numbered YD-12 during the same period is [25, 27, 29, 33, 30, 32, ...] μg / m³, with timestamps accurate to the second, such as 2026-03-01 08:01:20, and the synchronously recorded GPS location data is [121.502 degrees, 31.231 degrees].

[0016] Because fixed base stations and mobile devices have different data collection frequencies and dynamic differences in spatial location, and because mobile devices suffer from sensor response delays and GPS time synchronization errors, it is necessary to analyze the raw dataset and calculate the spatiotemporal difference characteristic values ​​between the two. These spatiotemporal difference characteristic values ​​include three indicators: average time difference, average spatial distance, and average absolute concentration difference. Using the fixed base station's minute-by-minute timestamp as a benchmark, data from mobile devices during the same period are filtered, and the average difference between the collection time and the benchmark, the average Euclidean distance between the collection location and the base station, and the average absolute difference between the concentration value and the benchmark concentration are calculated respectively. For example, for the base station baseline concentration of 28 μg / m³ at 08:01:00 on 2026-03-01, two data collection records from mobile devices between 08:01:00 and 08:01:59 were selected, with concentrations of 25 μg / m³ and 27 μg / m³, respectively. The average time difference was calculated to be 25 seconds, the average spatial distance was 250 meters, and the average absolute concentration difference was 2.5 μg / m³. The combination of these three factors forms the spatiotemporal difference characteristic value [25, 250, 2.5] for this time point.

[0017] To determine whether nonlinear deviations exist in the spatiotemporal difference characteristic values, three sets of spatiotemporal consistency thresholds were set. The time lag threshold was set to 60 seconds, determined by the sum of the maximum response delay of the mobile device sensor (30 seconds) and the maximum time synchronization error of the GPS module (30 seconds). The spatial influence threshold was set to 500 meters, based on the horizontal diffusion patterns of near-surface atmospheric pollutants in cities and measured data from urban building obstructions, determining 500 meters as the maximum spatial range within which fixed base station monitoring data has reference value for surrounding mobile device data. The concentration deviation threshold was set to 10 μg / m³, determined based on cross-device calibration experiments of fixed base stations and mobile devices simultaneously monitoring at the same location, establishing a normal upper limit for concentration deviation between the two as 10 μg / m³. If any indicator in the spatiotemporal difference characteristic values ​​exceeds the corresponding threshold, a nonlinear deviation is determined, and the corresponding spatiotemporal inconsistency distribution vector needs to be extracted. This vector includes the spatial deviation direction, the degree of time lag, and the magnitude of concentration abrupt changes. The spatial deviation direction is the average azimuth angle (0-360 degrees) of the mobile device data collection points with the base station as the origin; the time lag is the ratio of the average time difference to the time lag threshold (reflecting the proportion of time lag); and the concentration mutation amplitude is the ratio of the average absolute concentration difference to the base station reference concentration (eliminating the influence of concentration magnitude on the deviation). For example, the spatiotemporal difference characteristic value at a certain time node is [70, 600, 13] (all exceeding the corresponding threshold), indicating a nonlinear deviation; after calculation, the average azimuth angle is 135 degrees (southeast direction), the time lag is 1.17, and the concentration mutation amplitude is 0.37, ultimately forming a spatiotemporal inconsistency distribution vector [135 degrees, 1.17, 0.37].

[0018] The correction mapping matrix is ​​a third-order diagonal matrix, with the three diagonal elements representing the spatial correction coefficient, temporal correction coefficient, and concentration correction coefficient, respectively, and all other elements being 0. The spatial correction coefficient is calculated using a distance attenuation rule: when the average spatial distance is ≤ the spatial influence threshold of 500 meters, the spatial correction coefficient is 1; when the average spatial distance is > 500 meters, the spatial correction coefficient is 500 divided by the average spatial distance, used to reduce the weight of data from long-distance mobile devices. For example, for the aforementioned spatiotemporal inconsistency distribution vector [135 degrees, 1.17, 0.37], with a corresponding average spatial distance of 600 meters, the calculated correction coefficients are 0.833, 0.85, and 0.63, respectively, based on which the correction mapping matrix is ​​constructed.

[0019] The alignment process for the original dataset is performed in three steps: time, space, and concentration. First, based on the time correction coefficient, linear interpolation is used to interpolate the 30-second data from mobile devices to the 1-minute data, aligning it with the base station timestamps. Second, based on the spatial correction coefficient, weights are assigned to the interpolated mobile sampling concentrations, with data exceeding a 500-meter spatial threshold having reduced weights to minimize interference. Third, based on the concentration correction coefficient, bias compensation is applied to the weighted data, and spatiotemporal joint correction is completed using a correction mapping matrix. Simultaneously, the mean and standard deviation of the corrected data are calculated, and outliers exceeding ±3 times the standard deviation are removed. This threshold is set based on the normal distribution characteristics of air pollutant concentrations, and a range of 3 times the standard deviation covers 99.7% of the normal concentration data. For example, the base station baseline concentration at 08:01:00 on March 1, 2026 is 28 μg / m³. The interpolated concentrations of mobile devices are 26 μg / m³, 27 μg / m³, and 25 μg / m³, respectively. After correction with a spatial weight of 0.833 and a concentration coefficient of 0.63, the corrected concentrations are 13.7 μg / m³, 14.2 μg / m³, and 13.1 μg / m³, respectively. The base station baseline concentration and the corrected concentrations of each mobile device are integrated to form a fused dataset. A complete record is: 2026-03-01 08:01:00 | PM2.5 | 28 μg / m³ | [13.7,14.2,13.1] μg / m³.

[0020] In step S12, based on the fused dataset, Kalman filtering is used to predict and correct the base station reference concentration. The measurement update frequency is adjusted according to the correction result to obtain an optimized concentration sequence, which is then mapped to a grid to obtain the reference concentration distribution. This includes: The fused dataset was analyzed to obtain the base station baseline concentration sequence and the mobile sampling concentration subset; Calculate the short-term rate of change of the current reference concentration based on the base station reference concentration sequence, and couple the short-term rate of change with the current reference concentration to construct a reference concentration state vector. Using a preset state transition matrix, the baseline concentration state vector is predicted along the time dimension to obtain the prior state estimate and the state prediction error. Based on the data fluctuation and noise characteristics of the moving sampling concentration subset, the observation noise weight matrix is ​​calculated, and then the Kalman gain matrix is ​​calculated by combining the state prediction error. The prior state estimate is corrected using the Kalman gain matrix to obtain the posterior state estimate; The measurement update frequency is adjusted according to the convergence magnitude of the posterior state estimate, and the posterior state estimate is dynamically generated and integrated according to the adjusted frequency to obtain the optimized concentration sequence. Spatial interpolation combined with trend surface fitting is used to map the optimized concentration sequence onto the urban grid to obtain the corrected baseline concentration distribution.

[0021] First, it should be noted that the fused dataset is analyzed and split along two dimensions: "pollutant type + timestamp," resulting in a base station baseline concentration sequence and a mobile sampling concentration subset. The base station baseline concentration sequence represents stable concentration values ​​output by fixed base stations at continuous time points, reflecting the average pollution level in the region. The mobile sampling concentration subset is the corrected collection of sampling data from all mobile devices at the same time point, capturing local concentration fluctuation characteristics. For example, for PM2.5 pollutants, the base station baseline concentration sequence from 08:01:00 to 08:05:00 on March 1, 2026 is [28,30,29,32,31] μg / m³, and the corresponding mobile sampling concentration subsets are: 08:01:00 [24.3,25.1,23.8] μg / m³, 08:02:00 [26.5,27.2,25.9] μg / m³, 08:03:00 [25.7,26.3,24.9] μg / m³, 08:04:00 [28.9,29.5,27.8] μg / m³, and 08:05:00 [27.6,28.3,26.9] μg / m³.

[0022] The short-term rate of change is calculated using a sliding window method, with three consecutive time points as the window length. The window length is obtained by averaging the concentration differences between consecutive time points. This window length is set based on the temporal continuity of atmospheric pollutant concentration changes, balancing trend capture and noise suppression. The short-term rate of change is obtained by dividing the difference between the current concentration and the average concentration of the previous two time points by the time interval (minutes), with the time interval uniformly set to 1 minute. For example, if the base station baseline concentration at the current time (08:03:00) is 29 μg / m³, and the average concentration of the previous two time points is also 29 μg / m³, the calculated short-term rate of change is 0 μg / m³·min, and the baseline concentration state vector is constructed as [29,0].

[0023] The state transition matrix is ​​a 2×2 matrix based on the natural diffusion law of atmospheric pollutants. The core assumption is that concentration changes are gradual in the short term, and it is calibrated to [[0.99,1.0],[0.01,0.98]] using historical data. Diagonal elements are used to retain current state information, while off-diagonal elements are used to convey the trend of change. The state prediction error is obtained through statistical analysis of historical prediction deviations. The initial state prediction error matrix is ​​set to [[1.5,0],[0,0.5]], reflecting the uncertainty range of concentration and rate of change predictions. During prediction calculation, the baseline concentration state vector is multiplied by the state transition matrix to obtain the prior state estimate; the state prediction error is calculated by multiplying the state transition matrix, the initial state prediction error matrix, and the transpose of the transition matrix. For example, predicting the state vector [29,0] at 08:03:00 yields the prior state estimate vector [28.71,0.29]; the state prediction error is calculated to be [[1.485,0.015],[0.015,0.495]].

[0024] The construction of the observation noise weight matrix is ​​related to the reliability of the mobile sampling points. Reliability is determined by the distance of the sampling point from the base station, the sensor accuracy level, and the degree of environmental interference. The noise weight of each mobile sampling point is obtained by multiplying the sum of the distance coefficient and the environmental interference coefficient by the base noise weight. The base noise weight is set to 2.0 (μg / m³)² based on sensor calibration data. The distance coefficient is obtained by dividing the sampling point distance by 500 meters (spatial influence threshold). The environmental interference coefficient is set based on real-time meteorological data (1.0 for calm and stable weather, 1.2 for light wind, and 1.5 for strong wind, determined by actual measurements of the interference level to sensor detection under different meteorological conditions). The noise weights of all sampling points at the same time node are averaged to obtain the observation noise weight matrix for that node (a 1×1 matrix, as the observation value is a single dimension of concentration). For example, the mobile sampling concentration subset [25.7, 26.3, 24.9] μg / m³ at 08:03:00 corresponds to sampling points that are 300 meters away from the base station. The real-time weather is stable. The calculated distance coefficient for each sampling point is 0.6, the environmental interference coefficient is 1.0, the individual noise weight is 3.2 (μg / m³)², and the observation noise weight matrix is ​​[[3.2]]. The Kalman gain matrix is ​​calculated as follows: the inverse of the sum of the state prediction error and the observation noise weight matrix is ​​multiplied by the state prediction error to obtain the gain, which is used to balance the confidence levels of the predicted and observed values. For example, combining the aforementioned state prediction error matrix and observation noise weight matrix, the calculated Kalman gain matrix is ​​[[0.32, 0.01], [0.01, 0.11]].

[0025] The posterior state estimate is obtained by subtracting the predicted concentration from the prior state estimate from the mean concentration of the moving samples, multiplying by the Kalman gain, and then adding the prior state estimate. The mean concentration of the moving samples is the arithmetic mean of all corrected data at the same time point, used to reduce the impact of random errors at single sampling points. For example, the mean concentration of the moving samples at 08:03:00 is 25.63 μg / m³, with a concentration correction of -0.98 μg / m³ and a rate of change correction of -0.03 μg / m³·min; the corrected posterior state estimate vector is [27.73, 0.26].

[0026] The convergence amplitude is defined as the absolute difference between the concentration values ​​in the posterior state estimates of two consecutive time points, reflecting the stability of the corrected concentration. Two sets of update frequency adjustment thresholds are set, determined based on historical monitoring data: when the convergence amplitude is <1.5 μg / m³, it indicates that the pollution concentration changes slowly, and the measurement update frequency is reduced from the default 1 minute / time to 15 minutes / time to save computing resources; when the convergence amplitude is >4 μg / m³, it indicates that the concentration changes drastically, and the update frequency is increased to 2 minutes / time to quickly track the changing trend; when the convergence amplitude is between 1.5 and 4 μg / m³, the default frequency of 1 minute / time is maintained. For example, the corrected concentration value at 08:02:00 is 28.5 μg / m³, and the corrected concentration value at 08:03:00 is 27.73 μg / m³, with a convergence range of 0.77 μg / m³ < 1.5 μg / m³. The subsequent measurement update frequency is adjusted to 15 minutes / time. The posterior state estimates for each time point are continuously calculated at the adjusted frequency, and the optimized concentration values ​​are extracted and arranged in chronological order to obtain the optimized concentration sequence. For example, the adjusted optimized PM2.5 concentration sequence is [28.2, 28.5, 27.73, 27.9, 28.1, ...] μg / m³.

[0027] A three-dimensional spatiotemporal grid coordinate system is established for the urban monitoring area. The horizontal coordinate system is 500m x 500m (set according to the spatial resolution requirements of urban pollutant diffusion, balancing accuracy and computational load). The vertical coordinate system is layered in 50-meter intervals (covering the ground to a height of 200 meters, adapting to the vertical diffusion characteristics of near-ground pollutants). The temporal dimension is consistent with the update frequency of the optimized concentration sequence. Spatial interpolation uses the inverse distance-weighted method. The optimized concentration at the base station is multiplied by the sum of distance weights, and then divided by the sum of distance weights to obtain the grid point concentration. The distance weight is calculated as 1 divided by the distance from the grid point to the base station plus 50 meters (50 meters is a smoothing coefficient to avoid excessive weighting due to close distances).

[0028] Trend surface fitting is used to correct the spatial distribution trend after interpolation. The fitting model is set based on regional topography and pollution source distribution characteristics, and adopts a quadratic polynomial trend surface equation, which is calibrated through historical concentration data and topographic and pollution source data. For example, after interpolation and fitting of the grid in the main urban area of ​​a city, the corrected baseline concentration distribution is obtained: the concentration in the grid near the main traffic axis is 32-35 μg / m³, the concentration in the residential area is 27-30 μg / m³, and the concentration in the suburban area is 23-26 μg / m³, showing a continuous distribution trend of gradually decreasing from the main traffic axis to the suburbs. A complete baseline concentration distribution record is: 2026-03-01 08:15:00 | Grid number (121.50°E,31.23°N,50m) | PM2.5 | 33.2 μg / m³.

[0029] In step S13, the baseline concentration distribution is mapped to three-dimensional spatiotemporal space to obtain a baseline concentration field, which is then input into a preset prediction model to generate a prediction field. The concentration error within the evaluation range is assessed. If the evaluation result exceeds a preset accuracy threshold, the original dataset is supplemented by historical data smoothing and Kriging interpolation to obtain an optimized concentration distribution, including: The baseline concentration distribution is mapped to a preset three-dimensional spatiotemporal grid coordinate system to obtain the baseline concentration field, which is then input into a preset atmospheric pollutant diffusion prediction model to generate a prediction field. Calculate the difference between the predicted field and the reference concentration field within the effective coverage area to generate a grid point residual sequence; The root mean square error of the grid point residual sequence is calculated as the prediction error evaluation result; If the prediction error assessment result exceeds the preset accuracy threshold, then historical data of the base station is obtained, records matching the meteorological conditions are selected from them, and a subset of historical reference data is constructed. The historical reference data subset is subjected to exponential smoothing to generate a smoothed supplementary sequence. The smoothed supplementary sequence is then inserted into the original dataset, and kriging interpolation is performed to obtain the optimized concentration distribution.

[0030] First, it should be noted that the parameter settings of the three-dimensional spatiotemporal grid coordinate system are consistent with S12. The horizontal direction uses a 500m × 500m spatial grid, and the vertical direction is divided into 50m layers (a total of 4 layers, covering the near-surface layer from 0m to 200m). The time dimension is synchronized with the update frequency of the dynamically optimized concentration sequence. The concentration value mapping rule in the vertical dimension is set according to the characteristics of the atmospheric boundary layer structure and is calculated using an exponential decay model. That is, for every 50m increase in altitude, the concentration value is 0.9 times the concentration of the bottom layer. This decay coefficient is calibrated based on measured data of the vertical distribution of PM2.5 near the ground. For example, the baseline concentration of the bottom grid point (121.50°E, 31.23°N) is 33.2μg / m³, then the concentrations of the layers above are 29.9μg / m³, 26.9μg / m³, and 24.2μg / m³, respectively. Combined with the timestamp, a complete three-dimensional baseline concentration field record is formed.

[0031] The preset atmospheric pollutant diffusion prediction model adopts an improved Gaussian plume model, the specific form of which is as follows: ,in coordinates The concentration of pollutants, Because the pollution source is strong, The average wind speed, , These are the horizontal and vertical diffusion coefficients, respectively. For high-efficiency sources, , , Using spatial coordinates, the diffusion coefficient is adjusted to adapt to complex urban terrain. The model uses nearly one year of measured urban pollutant diffusion data as its training set, including concentration data from fixed base stations and mobile devices simultaneously monitored under different meteorological conditions and pollution source types, totaling over 100,000 samples. The training process employs gradient descent to iteratively optimize core parameters such as the diffusion coefficient and source height correction factor. The root mean square error between the model's predicted concentration and the measured concentration is used as the loss function until the loss function converges to below 0.8 μg / m³, completing model training and parameter calibration. It can simulate the transport and diffusion process of pollutants in three-dimensional space by incorporating real-time meteorological parameters such as wind speed, wind direction, temperature, and humidity. Using a baseline concentration field as the initial input field, and simultaneously inputting real-time meteorological data, the model extrapolates one hour along the time dimension to generate a prediction field. For example, the predicted concentration at the bottom layer of the aforementioned grid points at 09:15:00 on 2026-03-01 is 31.5 μg / m³. The effective coverage area is defined as the spatial grid actually traversed by the mobile device during the current monitoring period. The baseline concentration data within this area has been verified by mobile sampling data and has a higher confidence level than areas without mobile sampling coverage. The residual of a single grid point is obtained by subtracting the baseline concentration value from the predicted field concentration value. The residuals of all grid points within the effective coverage area are aggregated according to the dimension of "timestamp - grid number - height layer" to form a grid point residual sequence. For example, within the effective coverage area from 08:15:00 to 09:15:00 on March 1, 2026, there are 12 bottom-level grid points with residual values ​​of [-1.7, -2.1, 0.8, -3.2, -0.5, 1.2, -2.8, -1.3, 0.3, -4.1, -1.9, -0.9] μg / m³, which are combined to form the grid point residual sequence for this period. The root mean square error of the grid point residual sequence is calculated as the prediction error evaluation result. For example, the sum of squares of the above residual sequence is calculated to be 48.34, the mean square error is 4.03, and the root mean square error is approximately 2.01 μg / m³. This value is the error assessment result of this prediction.

[0032] Based on a statistical comparison of model predictions and measured data over the past three months, when the root mean square error (RMSE) is ≤1.5 μg / m³, the spatial distribution trend of the predicted field matches the actual concentration field by more than 90%. Therefore, the accuracy threshold is set at 1.5 μg / m³. If this threshold is exceeded, it indicates insufficient spatiotemporal coverage of the original dataset, leading to bias in the initial field of the model. The historical base station data consists of monitoring data from fixed base stations over the past year, including concentration data, real-time meteorological data, and timestamps, with a data granularity of 1 minute per instance. The screening rules for meteorological condition matching use Euclidean distance similarity calculation. Before the calculation, the four core meteorological parameters of wind speed, wind direction, temperature, and relative humidity are normalized to eliminate dimensional differences. Wind speed is linearly normalized according to the monitoring range of 0-8 m / s; wind direction is converted from a 360° circumference angle to two-dimensional coordinate values ​​and then normalized; temperature is linearly normalized according to the common urban monitoring range of 5-35℃; and relative humidity is linearly normalized according to the range of 0-100%. The similarity threshold is set at 0.8 (based on experimental comparison statistics), that is, data with Euclidean distance similarity ≥ 0.8 are selected. For example, given the current meteorological parameters as wind speed 3 m / s, wind direction southeast, temperature 18℃, and relative humidity 65%, 15 matching records are selected from historical data to construct a historical reference data subset. The PM2.5 concentration data within this subset is [30.2, 28.9, 31.5, 29.8, 32.1, 27.6, 28.3, 30.9, 29.2, 31.8, 28.7, 29.5, 30.1, 27.9, 31.2] μg / m³. A single exponential smoothing model is used, with a smoothing coefficient set to 0.7 (derived from the fluctuation characteristics analysis of historical data). The length of the smoothed supplementary sequence is consistent with the number of missing spatiotemporal grids in the current original dataset. For example, for the aforementioned historical reference data subset containing 15 data points, after a single exponential smoothing calculation, the generated smoothed supplementary sequence is [29.8, 30.1, 29.5, 30.7, 29.2] μg / m³. According to the corresponding spatiotemporal grid positions, it is inserted into the baseline concentration sequence of the original dataset to complete the spatiotemporal coverage supplementation of the dataset.

[0033] Ordinary Kriging interpolation (spherical variogram, range 1500 meters, determined by measured spatial correlation of pollutants) was used to interpolate the supplemented dataset to obtain an optimized concentration distribution. For example, after interpolating all 50 bottom-level grids in the main urban area, the PM2.5 concentration in the blank grids was accurately filled. For instance, the optimized concentration for grid (121.51°E, 31.24°N) was 29.8 μg / m³, and the optimized concentration for grid (121.49°E, 31.22°N) was 30.7 μg / m³. The final optimized concentration distribution is recorded as follows: 2026-03-01 08:15:00 | Grid number (121.51°E, 31.24°N, 50m) | PM2.5 | 29.8 μg / m³.

[0034] In step S14, the concentration diffusion trend is extracted from the optimized concentration distribution to generate a baseline concentration weight distribution. Sampled concentration interpolation data is obtained, and a refined concentration distribution field is constructed by applying a grid resolution adjustment method in conjunction with the baseline concentration weight distribution. This includes: The concentration diffusion trend of the optimized concentration distribution is extracted using optical flow. The concentration diffusion trend is projected onto the baseline concentration field to generate a baseline concentration weight distribution; The moving sampling concentration is interpolated to obtain the sampling concentration interpolation data; Using the aforementioned baseline concentration weight distribution as the density control function, quadtree mesh refinement is performed to construct a multi-scale adaptive mesh set; Resampling is performed based on the multi-scale adaptive grid set and the sampled concentration interpolation data to construct a refined concentration distribution field.

[0035] First, it should be noted that when using the optical flow method, the time window is set to 3 consecutive update cycles (based on the update frequency of S12). For each grid point, the concentration change rate scalar of that point is obtained by the ratio of the concentration difference between adjacent time frames to the time interval. The concentration diffusion trend vector of the region is determined by combining the concentration gradient direction of the adjacent grid points (composed of two dimensions: diffusion main direction and diffusion rate). The diffusion main direction is determined by calculating the direction of the maximum concentration gradient between the target grid point and the surrounding 8 neighboring grid points. The diffusion rate is the result of spatial normalization of the concentration change rate scalar of the target grid point. For example, for the PM2.5 optimized concentration distribution of three consecutive update cycles on March 1, 2026, at 08:15:00, 08:30:00, and 08:45:00, the concentration of the target grid point (121.50°E, 31.23°N) decreased from 33.2 μg / m³ to 32.1 μg / m³ and then to 31.5 μg / m³ within this time window, and the scalar of the concentration change rate was calculated to be 0.36 km / h. At the same time, the concentration gradient of the 8-neighborhood of this grid point was calculated, and it was found that the concentration of the northeast grid point (121.505°E, 31.23°N) decreased from 30.1 μg / m³ to 29.9 μg / m³ and then to 29.8 μg / m³, which is the direction of the maximum concentration gradient. Therefore, the concentration diffusion trend vector of this region was determined to be [northeast direction, 0.36 km / h].

[0036] The projection rules are set based on the influence range and intensity of the diffusion trend. The weight calculation uses a Gaussian kernel function, where the baseline weight is set to 1.0 (the default weight when there is no diffusion influence). The diffusion radius (0.27 km in this case) is obtained by multiplying the diffusion velocity by the time window duration, ensuring that the weight decreases reasonably with distance. For example, the weight of the core grid point (121.50°E, 31.23°N) is 1.2 (a 20% increase in weight for the core diffusion source region), the weight of the northeast grid point (0.1 km away) is 1.1, the weight of the grid point (0.27 km away) is 0.8, and the weight of the grid point (more than 0.5 km away) remains unchanged at 1.0, forming a gradient-decreasing baseline concentration weight distribution.

[0037] The interpolation method employed was inverse distance weighted interpolation, with a search radius set to 1 km (to suit the sampling density). The weighting index was set to 2, which, after multiple experimental verifications, effectively balanced the influence of local and surrounding data. For example, the PM2.5 concentrations at three mobile sampling points within a certain area were 24.3 μg / m³ (121.502°E, 31.231°N), 25.1 μg / m³ (121.503°E, 31.232°N), and 23.8 μg / m³ (121.504°E, 31.233°N), respectively. Interpolation was performed on the blank grid point (121.503°E, 31.2315°N), yielding an interpolated concentration of 24.7 μg / m³.

[0038] Quadtree partitioning sets partitioning level thresholds based on weight gradients: a weight ≥ 1.1 indicates Level 1 partitioning (100m × 100m resolution), a weight ≤ 0.9 < 1.1 indicates Level 2 partitioning (250m × 250m resolution), and a weight < 0.9 maintains the original resolution (500m × 500m). These thresholds are determined based on the gradient characteristics of the weight distribution to ensure that the partitioned grid accurately matches sensitive areas of concentration changes. For example, in the aforementioned baseline concentration weight distribution, the core region with a weight of 1.2 is partitioned into a fine 100m × 100m grid, and the diffusion influence region with weights of 1.0-1.1 is partitioned into a medium 250m × 250m grid, constructing a multi-scale adaptive grid set.

[0039] The resampling rules are as follows: for the fine grid of the first-level subdivision, bilinear interpolation is used to extract the concentration value at the corresponding location from the sampled concentration interpolation data; for the medium grid of the second-level subdivision, the nearest neighbor interpolation method is used to quickly obtain the concentration value; for the coarse grid, the original data of the optimized concentration distribution is directly used. At the same time, a concentration continuity verification rule is set: the concentration difference between adjacent grids shall not exceed 5 μg / m³ (set according to the spatial continuity characteristics of atmospheric pollutant concentrations). If the threshold is exceeded, the average concentration of the surrounding grids is used for correction. For example, in a 100m × 100m fine grid in the core area, the resampling concentration of a certain intersection grid is 27.3μg / m³, and the concentrations of the four surrounding adjacent grids are 26.8μg / m³, 27.5μg / m³, 26.9μg / m³, and 27.1μg / m³, respectively. The concentration transition is smooth and there is no abrupt change. The concentration of the medium grid in the diffusion influence area is gradient-distributed between 24-26μg / m³, and the concentration of the coarse grid in the suburbs is maintained at 23-25μg / m³, thus forming a refined concentration distribution field.

[0040] In step S15, the concentration distribution field is subjected to boundary correction, discontinuous boundaries are identified and numerical assimilation is performed to obtain a smooth, continuous concentration field, including: A global scan of the concentration distribution field is performed to identify discontinuous boundaries where concentration values ​​abruptly change. By combining the effective concentration data of the surrounding grid, the discontinuous boundary is interpolated and corrected to obtain the corrected concentration distribution field; A global smoothing process is performed on the modified concentration distribution field to obtain a smooth continuous concentration field.

[0041] First, it should be noted that the global scan employs a neighborhood gradient detection method, covering all grid points in the refined concentration distribution field, and calculating the concentration gradient between each grid point and its eight neighbors. Based on the physical characteristics of atmospheric turbulent diffusion, the natural gradient of pollutant concentration within 500 meters of the ground, in the absence of sudden emission sources, typically does not exceed 3.5 μg / m³ between adjacent grid points; therefore, the concentration gradient threshold is set to 4.0 μg / m³. If the absolute value of the concentration difference between a grid point and any neighboring grid point exceeds 4.0 μg / m³, and this abrupt change occurs at three or more consecutive grid points, then the set of these grid points is marked as a discontinuous boundary. For example, in the PM2.5 fine concentration distribution field at 09:00:00 on March 1, 2026, the scan found that the concentration of the grid point (121.505°E, 31.235°N) was 35.2 μg / m³, and the concentration of the grid point in its western neighboring area (121.504°E, 31.235°N) was 28.1 μg / m³, with a gradient of 7.1 μg / m³. Moreover, this abrupt change extended continuously for 4 grids in the northeast direction, indicating that there was a discontinuous boundary in this area.

[0042] The interpolation correction employs a numerical assimilation approach, using grid points at discontinuous boundaries as the assimilation core to construct a local influence window. The window radius is set to 3 grid cells (determined based on the 1500-meter range of pollutant spatial correlation). Abnormal grid points marked as discontinuous boundaries within the assimilation window are removed, retaining only normal, valid grid points as interpolation samples. A weighted inverse distance interpolation method is used to fit the concentration values ​​of the sample points to the abnormal grid points, with a weighting index of 3 to enhance the correction effect of nearest-neighbor valid data on abnormal points and quickly eliminate abrupt changes. For example, for the core abnormal point at the aforementioned discontinuous boundary, the valid sample point concentrations within its 3 grid cell range are 28.1 μg / m³, 27.9 μg / m³, and 29.2 μg / m³, respectively. After interpolation correction, the concentration of this abnormal point is adjusted from 35.2 μg / m³ to 28.5 μg / m³, and the gradient with the surrounding concentration decreases to 0.4 μg / m³. The corrected concentration distribution field is obtained after completing all boundary corrections.

[0043] Global smoothing employs a Gaussian filter convolution algorithm, setting the kernel function window to a 3×3 grid and the standard deviation σ to 1.0. Verification using experimental data shows that this setting effectively smooths minor fluctuations remaining after interpolation correction without blurring the boundary between the core contaminated area and the background area. For each grid point, the concentration values ​​of that point and its eight neighboring grid points are multiplied by a Gaussian kernel weight coefficient and summed to obtain the smoothed concentration value for that grid point. For example, if the corrected concentration of a grid point is 28.5 μg / m³, and its eight neighboring concentrations are between 27.8 μg / m³ and 29.0 μg / m³, after Gaussian filtering smoothing, the concentration is determined to be 28.4 μg / m³. After smoothing, it is ensured that the concentration values ​​of all grid points are non-negative, and the concentration gradients of adjacent grids are controlled within a threshold range of 4.0 μg / m³. After successful verification, a smoothed continuous concentration field is obtained.

[0044] In step S16, if the smooth continuous concentration field shows a migration trend, the spatiotemporal changes of the migration path are tracked and a grid deformation control vector is generated. The spatiotemporal grid of the prediction field is dynamically adjusted to obtain a dynamic concentration prediction field, including: Based on the smooth, continuous concentration field, the Lagrange particle tracking method is used to identify regions of significant concentration migration. Based on the regions of significant concentration migration, a spatiotemporal evolution manifold is constructed to deduce the spatiotemporal trajectory of the migration path; The spatial drift of the spatiotemporal change trajectory is calculated by combining the preset parameter update cycle, and a grid deformation control vector is generated. The predicted field spatiotemporal grid is stretched using the grid deformation control vector to construct an adaptive spatiotemporal grid system; The baseline concentration field, the optimized concentration distribution, and the sampled concentration interpolation data are mapped to the adaptive spatiotemporal grid system, and flux conservation interpolation calculations are performed to obtain a preliminary dynamic concentration prediction field.

[0045] First, it should be noted that the Lagrange particle tracking method releases tracer particles only in the high-concentration core area and concentration gradient zone. The threshold for determining the high-concentration core area is 1.5 times the current field average concentration. This threshold is set according to the ambient air quality index classification standard to ensure that particles are concentrated in areas with significant pollution impact. Ten tracer particles are released for each target grid, and the initial particle velocity is initialized by the optical flow velocity field of the smooth continuous concentration field. If, within three consecutive time steps, the displacement of the tracer particle swarm in a certain area exceeds 200 meters (i.e., exceeds 40% of the basic grid resolution), and the particle swarm aggregation degree remains above 0.8, then the area is marked as a region of significant concentration migration. For example, in the smooth continuous PM2.5 concentration field at 09:00:00 on March 1, 2026, the tracer particle swarm around the core grid (121.50°E, 31.23°N) shifted 280 meters northeastward within 45 minutes, with an aggregation degree of 0.85, and this area was determined to be a region of significant concentration migration.

[0046] The spatiotemporal evolution manifold is based on regions of significant concentration migration over a continuous time series. The region contour at each time step is treated as a slice of the manifold. Adjacent slices are spatiotemporally fused using thin-plate spline interpolation to form a continuous four-dimensional evolution manifold. Based on the local curvature of the manifold and historical displacement rates, the migration direction and distance within the next parameter update cycle are predicted. For example, based on the spatiotemporal evolution manifold from 08:15:00 to 09:00:00, it is deduced that the PM2.5 pollution plume will migrate in a northeast-east direction at a speed of 0.42 km / h within the next 15 minutes, with an estimated displacement of 105 meters.

[0047] The parameter update cycle is consistent with the S12 dynamic update frequency (15 minutes in this case, set according to the short-term diffusion response characteristics of atmospheric pollutants, taking into account both timeliness and deformation rationality). The displacement vector of the migration path within the update cycle is used as the core area grid drift reference, and the drift amount of the non-core area decreases linearly according to the distance from the core area (attenuation coefficient 0.5, to ensure smooth deformation). The grid deformation control vector is a two-dimensional vector field (including eastward and northward drift components). For example: the core migration area is [80 meters eastward, 65 meters northward], the attenuation of the outer area 500 meters away is [40 meters eastward, 32.5 meters northward], and the background area is [0,0].

[0048] The mesh stretching employs a dynamic mesh adjustment algorithm based on a spring-mass model. The spatiotemporal mesh nodes of the original prediction field are treated as masses, the connections between nodes as springs, and the mesh deformation control vector as an external force driving the mass motion. The mesh in the core migration region is stretched in the migration direction and compressed in the vertical direction. The stretched mesh resolution remains 100m × 100m in the core region and smoothly transitions to 500m × 500m in the transition region, avoiding mesh distortion. For example, the five 500m × 500m meshes covering the pollution plume in the original prediction field are stretched and merged into eight 100m × 100m fine meshes, extending in a northeast-east direction, constructing an adaptive spatiotemporal mesh system that follows the movement of the pollution plume.

[0049] The flux conservation interpolation method first calculates the concentration flux of multi-source concentration data in the original grid, and then distributes the flux proportionally to the new grid according to the area and volume of the adaptive spatiotemporal grid. The baseline concentration field provides regional background concentration constraints, the optimized concentration distribution provides historical supplementary data support, and the sampled concentration interpolation data provides real-time local calibration. For example, after performing flux conservation interpolation on the baseline concentration field (background concentration 28.0 μg / m³), optimized concentration distribution (core concentration 33.2 μg / m³), and sampled concentration interpolation data (local concentration 24.7 μg / m³) at 09:00:00 on March 1, 2026, the predicted concentration in the core migration grid of the adaptive grid system is 31.8 μg / m³, the transition region is 29.5 μg / m³, and the background region is 28.2 μg / m³. The concentration distribution dynamically adjusts with the migration path while maintaining total concentration conservation.

[0050] In step S17, a weighted baseline concentration distribution is calculated based on the dynamic concentration prediction field. A gradient deviation field is constructed based on the weighted baseline concentration distribution and the sampled concentration interpolation data. The parameters of the preset prediction model are iteratively corrected based on the gradient deviation field to obtain a high-precision local concentration field, including: Based on the dynamic concentration prediction field, a weighted baseline concentration distribution is generated by combining the preset historical average concentration and spatial distance weights. The concentration difference between the weighted baseline concentration distribution and the corresponding sampled concentration interpolation data is calculated for each grid point to obtain the moving sampling residual set. Spatial interpolation is then performed on the moving sampling residual set to obtain the gradient bias field. The gradient correction vector is calculated based on the spatial distribution of the gradient deviation field. The gradient correction vector is then used to adjust the local diffusion coefficient and model decay factor of the preset prediction model to generate the corrected model variable matrix. The modified model variable matrix is ​​input into the concentration evolution equation to perform iterative calculations until the results converge, ultimately yielding a high-precision local concentration field.

[0051] First, it should be noted that the preset historical average concentration is the statistical average of pollutant concentrations at fixed base stations during the same period and under the same meteorological conditions over the past year. This value is set based on the seasonal and diurnal variations in air pollutant concentrations. For example, the preset historical average PM2.5 concentration during the morning rush hour in the main urban area of ​​a city in spring is 26 μg / m³. Spatial distance weighting is calculated using an inverse distance weighting method, with the fixed base station as the center. Let the straight-line distance between the grid point and the base station be... The spatial distance weight is defined as follows: The 500-meter threshold is the spatial influence threshold set by S11, and 0.1 is a smoothing coefficient to avoid meaningless weighting when the distance is 0. Based on measured data, a dynamic prediction field weight of 0.6 ensures the tracking of real-time pollution trends, and a historical average weight of 0.4 is used to suppress local random fluctuations in the prediction field. The dynamic concentration prediction field concentration is multiplied by 0.6, and the product of the preset historical average concentration and the spatial distance weight of 0.4 is calculated separately. The weighted baseline concentration is obtained after multiplying the two products. For example, the dynamic prediction concentration of a grid point 200 meters away from the base station is 31 μg / m³, and the spatial distance weight is 500 / (200+500) = 0.714 (rounded to two decimal places). The weighted baseline concentration of this point is: 31*0.6+26*0.4*0.714=26.0 μg / m³; the dynamic prediction concentration of a grid point 600 meters away from the base station is 28 μg / m³, so the weighted baseline concentration can be calculated as 21.5 μg / m³. Calculate the weighted baseline concentration for all grid points to form the weighted baseline concentration distribution covering the monitoring area.

[0052] Subtracting the weighted baseline concentration distribution from the sampled concentration interpolation data yields the residual for each grid point. The residual values ​​of all grid points are then aggregated along the dimension of "grid number - residual value" to form a moving sample residual set. For example, this set could be [(grid 1, 3.4), (grid 2, -1.5), (grid 3, 2.8), (grid 4, -0.9), ...] μg / m³. Spatial interpolation of the moving sample residual set is performed using ordinary Kriging interpolation (set consistent with S13) to obtain the gradient deviation field. For example, the gradient deviation field near major traffic arteries shows a high positive deviation value, while suburban areas show a slight negative deviation.

[0053] The gradient correction vector is a two-dimensional vector containing the spatial gradient correction coefficient and the concentration deviation correction coefficient. The spatial gradient correction coefficient is calculated by dividing the rate of change of concentration in a certain direction of the gradient deviation field by the preset baseline rate of change, reflecting the trend of deviation with spatial location. The concentration deviation correction coefficient is calculated by dividing the average value of the gradient deviation field by the average value of the dynamic concentration prediction field, reflecting the overall degree of deviation.

[0054] The preset baseline change rate is the natural spatial change rate of air pollutants under stable weather conditions, set at 0.05 μg / m³ / m. The local diffusion coefficient of the prediction model is used to characterize the diffusion capacity of pollutants in space, and the model attenuation factor is used to characterize the degree of concentration attenuation of pollutants during transport. The corrected local diffusion coefficient is obtained by adding the spatial gradient correction coefficient to 1 and multiplying it by the original local diffusion coefficient; the corrected model attenuation factor is obtained by subtracting the concentration deviation correction coefficient from 1 and multiplying it by the original model attenuation factor. For example, if the original local diffusion coefficient is 0.8 m² / s, the original model attenuation factor is 0.015 / h, and the gradient correction vector for a certain region is [0.23, 0.12], the calculated corrected local diffusion coefficient is 0.984 m² / s, and the corrected model attenuation factor is 0.0132 / h. The corrected local diffusion coefficient, model attenuation factor, and other core parameters of the prediction model are integrated dimensionally to generate the corrected model variable matrix. For example, a row of data in the model variable matrix is ​​[grid 5, 0.984 m² / s, 0.0132 / h, 1.05, 0.98], which corresponds to the grid number, local diffusion coefficient, model attenuation factor, wind speed correction coefficient, and wind direction influence coefficient, respectively.

[0055] The concentration evolution equation adopts a convection-diffusion equation based on atmospheric dynamics, specifically in the form of: ,in For pollutant concentration, For time, , , These represent the wind speed components in three dimensions. , , These are the local diffusion coefficients in three dimensions. This is the model decay factor. The pollution source strength is set to 0 (when there are no new emissions). The initial value for the iterative calculation is the concentration value of the dynamic concentration prediction field. The corrected model variable matrix is ​​substituted into the equation, and the iterative calculation is performed with a time step of 1 minute and a spatial step consistent with the grid resolution. The convergence criterion is that after three consecutive iterations, the concentration change value of all grid points is less than 0.1 μg / m³ (set according to the accuracy requirements of environmental monitoring). For example, after 5 iterations, the calculation results meet the convergence criterion. The concentration of grid points near the main traffic artery is corrected from the original prediction of 31 μg / m³ to 32.3 μg / m³, and the deviation from the actual sampling data is reduced from 3.4 μg / m³ to 0.2 μg / m³; the concentration of grid points in the suburbs is corrected from the original prediction of 25 μg / m³ to 24.8 μg / m³, and the deviation is reduced from -1.5 μg / m³ to 0.4 μg / m³. After iterative convergence, the final concentrations of each grid point are integrated according to "timestamp - grid number - pollutant type - concentration value" to form a high-precision local concentration field.

[0056] To achieve overall closed-loop optimization of the system, after obtaining the local concentration field in step S17, the following steps are also included: The concentration observations from multiple real-time sensors are acquired, and the instantaneous prediction error tensor is generated by differential operation in combination with the local concentration field. The Holt-Winters smoothing algorithm is applied to the instantaneous prediction error tensor to obtain a smoothed error curve. If the second derivative of the smoothed error curve is greater than a preset abrupt change threshold, a high-frequency resampling instruction is generated. The state observation vector of the Kalman filter is reconstructed according to the high-frequency resampling command, and the dynamic feedback gain matrix is ​​calculated based on the state observation vector. A closed-loop correction factor is generated based on the dynamic feedback gain matrix. The closed-loop correction factor is then incorporated into the concentration evolution iteration calculation of the preset prediction model at the next time step to obtain the corrected prediction model.

[0057] First, it should be noted that the multi-source real-time sensors include the fixed base stations and mobile monitoring equipment mentioned in S11, as well as the supplementary roadside miniature monitoring terminals. The data acquisition frequency is uniformly set at the highest priority of 30 seconds per acquisition to ensure the capture of instantaneous concentration changes. The instantaneous prediction error tensor is constructed using the dimension of "time-space-pollutant type". The corresponding predicted value is extracted from the local concentration field using the nearest neighbor interpolation method according to the sensor's geographical coordinates. The error value is obtained by subtracting the model prediction value from the sensor observation value. All sensor error values ​​within five consecutive time steps are arranged in tensor form to form the instantaneous prediction error tensor. For example, at 09:32:30 on March 1, 2026, the observed value of sensor YD-12 was 35.6 μg / m³, the corresponding predicted value of the local concentration field was 32.3 μg / m³, and the error value was 3.3 μg / m³; the observed value of sensor MX-05 (roadside terminal) was 26.1 μg / m³, the predicted value was 28.4 μg / m³, and the error value was -2.3 μg / m³. Combining the data from the previous four time steps, a three-dimensional instantaneous prediction error tensor containing multiple positive and negative error values ​​is generated.

[0058] A Holt-Winters model without seasonal terms is adopted, with a horizontal smoothing coefficient of 0.8 and a trend smoothing coefficient of 0.2. This parameter combination is calibrated based on the sensor noise characteristics, enabling rapid response to the latest error changes and smoothing trend fluctuations. The preset abrupt change threshold is set to 0.05 μg / m³ / s², which is determined based on the rate of change of the sensor's maximum measurement error. When the second derivative of the smoothed error curve exceeds this value, high-frequency monitoring needs to be initiated. For example, after smoothing the above error tensor, the error trend rapidly rises from 0.5 μg / m³ to 3.3 μg / m³, and the calculated second derivative is 0.07 μg / m³ / s², which is greater than the threshold of 0.05. The system immediately generates a high-frequency resampling command, increasing the sampling frequency of the mobile device and micro-terminal from 30 seconds / time to 5 seconds / time.

[0059] The Kalman filter is reconstructed, expanding the original two-dimensional state vector [concentration, rate of change] into a three-dimensional vector [concentration, rate of change, error trend], where the error trend term is directly provided by the trend component of the Holt-Winters model. The state transition matrix is ​​synchronously updated to a 3×3 matrix, with the added dimension used to transmit the historical trend of the error. Based on the covariance matrix of the high-frequency resampled data, the observation noise covariance is updated in real time; the time-varying Kalman gain matrix is ​​calculated by combining the state prediction error covariance. For example, in the high-frequency resampling mode, the observation data density is increased by 6 times, the observation noise covariance is reduced from 3.2 to 0.8, and the calculated dynamic feedback gain matrix is ​​[[0.75,0.02,0.15],[0.02,0.60,0.10],[0.15,0.10,0.85]].

[0060] The closed-loop correction factor is divided into a state correction factor and a parameter correction factor. The state correction factor is determined by the product of the Kalman gain and the latest observation residual, and is used to correct the prior state estimate at the next time step. The parameter correction factor is the trace of the dynamic feedback gain matrix, and is used to dynamically adjust the diffusion coefficient of the concentration evolution equation in S17. Before each model iteration, the initial concentration values ​​of the grid points are first updated using the state correction factor; then, the parameter correction factor is multiplied by the original local diffusion coefficient to obtain the adaptive diffusion coefficient. For example, the trace of the dynamic feedback gain matrix is ​​2.2, and the parameter correction factor is taken from the trace itself. The original diffusion coefficient of 0.984 m² / s is adjusted to 2.165 m² / s (0.984 × 2.2 ≈ 2.165), which significantly improves the model's ability to simulate sudden diffusion. The above closed-loop correction mechanism is solidified into the computational framework of the original prediction model to form a corrected prediction model.

[0061] In summary, this invention discloses a multi-sensor collaborative air pollution monitoring method. Through a fusion strategy of coordinated sampling by fixed base stations and mobile devices and spatiotemporal inconsistency correction, the static base station baseline concentration and the dynamic mobile sampling concentration are aligned to a unified spatiotemporal framework, effectively eliminating the spatiotemporal scale misalignment problem of multi-source monitoring data. By constructing a full-process monitoring system, the background concentration distribution is first optimized using Kalman filtering, then pollutant migration characteristics are captured and the grid is dynamically adjusted using optical flow and Lagrange particle tracking methods, and finally, model parameters are iteratively corrected through gradient bias field, significantly improving the spatiotemporal characterization accuracy and dynamic prediction capability of the pollution concentration field. Through a linkage mechanism of dynamic error evaluation and Kalman filter closed-loop feedback, changes in concentration prediction error are captured in real time and high-frequency resampling is triggered. The filtered observation vector is reconstructed and a closed-loop correction factor is generated to continuously optimize the model, enabling rapid response to sudden pollution and continuous improvement of prediction accuracy.

[0062] Reference Figure 2 The second embodiment of the present invention provides a multi-sensor collaborative air pollution monitoring system, comprising: The data acquisition and fusion module is used to collect base station baseline concentrations through fixed base stations and obtain mobile sampling concentrations through mobile devices to form a raw dataset. The data is then aligned using spatiotemporal inconsistency correction to obtain a fused dataset. The Kalman filter correction module is used to predict and correct the base station reference concentration based on the fused dataset using Kalman filtering, adjust the measurement update frequency according to the correction result, obtain an optimized concentration sequence, and map it to a grid to obtain the reference concentration distribution. The concentration field prediction and optimization module is used to map the baseline concentration distribution to three-dimensional spatiotemporal space to obtain the baseline concentration field and input it into the preset prediction model to generate the prediction field. It evaluates the concentration error within the range. If the evaluation result exceeds the preset accuracy threshold, it supplements the original dataset by smoothing historical data and Kriging interpolation to obtain the optimized concentration distribution. The refined construction module is used to extract the concentration diffusion trend from the optimized concentration distribution, generate a baseline concentration weight distribution, obtain sampled concentration interpolation data, and combine the baseline concentration weight distribution with a grid resolution adjustment method to construct a refined concentration distribution field. The boundary correction processing module is used to correct the boundary of the concentration distribution field, identify discontinuous boundaries and perform numerical assimilation processing to obtain a smooth continuous concentration field. The dynamic migration prediction module is used to track the spatiotemporal changes of the migration path and generate a grid deformation control vector if the smooth continuous concentration field shows a migration trend, and dynamically adjust the spatiotemporal grid of the prediction field to obtain a dynamic concentration prediction field. The final output module is used to calculate the weighted baseline concentration distribution based on the dynamic concentration prediction field, construct a gradient deviation field based on the weighted baseline concentration distribution and the sampled concentration interpolation data, and iteratively correct the parameters of the preset prediction model based on the gradient deviation field to obtain a high-precision local concentration field.

[0063] It should be noted that the multi-sensor collaborative air pollution monitoring system provided in this embodiment of the invention is used to execute all the process steps of the multi-sensor collaborative air pollution monitoring method in the above embodiment. The working principles and beneficial effects of the two are one-to-one, so they will not be described again.

[0064] It should be noted that the system embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs. Furthermore, in the accompanying drawings of the system embodiments provided by this invention, the connection relationships between modules indicate that they have communication connections, which can be specifically implemented as one or more communication buses or signal lines. Those skilled in the art can understand and implement this without any creative effort.

[0065] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the scope of protection of the present invention. In particular, it should be noted that any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention for those skilled in the art.

Claims

1. A method of air pollution monitoring by multi-sensor coordination, characterized by, include: The base station baseline concentration is collected by a fixed base station, and the mobile sampling concentration is obtained by a mobile device to form a raw dataset. The data is then aligned by spatiotemporal inconsistency correction to obtain a fused dataset. Based on the fused dataset, Kalman filtering is used to predict and correct the base station baseline concentration. The measurement update frequency is adjusted according to the correction result to obtain an optimized concentration sequence, which is then mapped to a grid to obtain the baseline concentration distribution. The baseline concentration distribution is mapped to three-dimensional spatiotemporal space to obtain a baseline concentration field, which is then input into a preset prediction model to generate a prediction field. The concentration error within the evaluation range is evaluated. If the evaluation result exceeds a preset accuracy threshold, the original dataset is supplemented by historical data smoothing and Kriging interpolation to obtain an optimized concentration distribution. The concentration diffusion trend is extracted from the optimized concentration distribution to generate a baseline concentration weight distribution. Sampled concentration interpolation data is obtained, and a fine concentration distribution field is constructed by applying a grid resolution adjustment method in combination with the baseline concentration weight distribution. Boundary correction is performed on the concentration distribution field, discontinuous boundaries are identified and numerical assimilation is performed to obtain a smooth continuous concentration field; If the smooth continuous concentration field shows a migration trend, the spatiotemporal changes of the migration path are tracked and a grid deformation control vector is generated to dynamically adjust the spatiotemporal grid of the prediction field, thereby obtaining a dynamic concentration prediction field. Based on the dynamic concentration prediction field, a weighted baseline concentration distribution is calculated. A gradient deviation field is constructed based on the weighted baseline concentration distribution and the sampled concentration interpolation data. The parameters of the preset prediction model are iteratively corrected based on the gradient deviation field to obtain a high-precision local concentration field.

2. The multi-sensor synergistic air pollution monitoring method of claim 1, wherein, The process involves collecting base station baseline concentrations from fixed base stations and acquiring mobile sampling concentrations using mobile devices to form an original dataset. Spatiotemporal inconsistency correction is then used to align the data, resulting in a fused dataset, which includes: Obtain the base station baseline concentration collected by fixed base stations and the mobile sampling concentration collected by mobile devices to construct the original dataset; The original dataset is parsed to calculate the spatiotemporal difference feature values ​​between the base station baseline concentration and the mobile sampling concentration; If the spatiotemporal difference feature value has a nonlinear deviation, then the corresponding spatiotemporal inconsistency distribution vector is extracted. The spatiotemporal inconsistency distribution vector includes the spatial deviation direction, the degree of time lag, and the magnitude of concentration change. A correction mapping matrix is ​​constructed based on the spatiotemporal inconsistency distribution vector, and alignment processing is performed on the original dataset to obtain the fused dataset.

3. The multi-sensor coordinated air pollution monitoring method of claim 1, wherein, Based on the fused dataset, Kalman filtering is used to predict and correct the base station baseline concentration. The measurement update frequency is adjusted according to the correction result to obtain an optimized concentration sequence, which is then mapped to a grid to obtain the baseline concentration distribution. This includes: The fused dataset was analyzed to obtain the base station baseline concentration sequence and the mobile sampling concentration subset; Calculate the short-term rate of change of the current reference concentration based on the base station reference concentration sequence, and couple the short-term rate of change with the current reference concentration to construct a reference concentration state vector. Using a preset state transition matrix, the baseline concentration state vector is predicted along the time dimension to obtain the prior state estimate and the state prediction error. Based on the data fluctuation and noise characteristics of the moving sampling concentration subset, the observation noise weight matrix is ​​calculated, and then the Kalman gain matrix is ​​calculated by combining the state prediction error. The prior state estimate is corrected using the Kalman gain matrix to obtain the posterior state estimate; The measurement update frequency is adjusted according to the convergence magnitude of the posterior state estimate, and the posterior state estimate is dynamically generated and integrated according to the adjusted frequency to obtain the optimized concentration sequence. Spatial interpolation combined with trend surface fitting is used to map the optimized concentration sequence onto the urban grid to obtain the corrected baseline concentration distribution.

4. The multi-sensor collaborative air pollution monitoring method according to claim 1, characterized in that, The process involves mapping the baseline concentration distribution to three-dimensional spatiotemporal space to obtain a baseline concentration field, which is then input into a preset prediction model to generate a prediction field. The concentration error within the evaluation range is assessed. If the evaluation result exceeds a preset accuracy threshold, the original dataset is supplemented using historical data smoothing and Kriging interpolation to obtain an optimized concentration distribution. This includes: The baseline concentration distribution is mapped to a preset three-dimensional spatiotemporal grid coordinate system to obtain the baseline concentration field, which is then input into a preset atmospheric pollutant diffusion prediction model to generate a prediction field. Calculate the difference between the predicted field and the reference concentration field within the effective coverage area to generate a grid point residual sequence; The root mean square error of the grid point residual sequence is calculated as the prediction error evaluation result; If the prediction error assessment result exceeds the preset accuracy threshold, then historical data of the base station is obtained, records matching the meteorological conditions are selected from them, and a subset of historical reference data is constructed. The historical reference data subset is subjected to exponential smoothing to generate a smoothed supplementary sequence. The smoothed supplementary sequence is then inserted into the original dataset, and kriging interpolation is performed to obtain the optimized concentration distribution.

5. The multi-sensor collaborative air pollution monitoring method according to claim 1, characterized in that, The process of extracting the concentration diffusion trend from the optimized concentration distribution, generating a baseline concentration weight distribution, obtaining sampled concentration interpolation data, and constructing a refined concentration distribution field by applying a grid resolution adjustment method in conjunction with the baseline concentration weight distribution includes: The concentration diffusion trend of the optimized concentration distribution is extracted using optical flow. The concentration diffusion trend is projected onto the baseline concentration field to generate a baseline concentration weight distribution; The moving sampling concentration is interpolated to obtain the sampling concentration interpolation data; Using the aforementioned baseline concentration weight distribution as the density control function, quadtree mesh refinement is performed to construct a multi-scale adaptive mesh set; Resampling is performed based on the multi-scale adaptive grid set and the sampled concentration interpolation data to construct a refined concentration distribution field.

6. The multi-sensor collaborative air pollution monitoring method according to claim 1, characterized in that, The step of boundary correction of the concentration distribution field, identifying discontinuous boundaries and performing numerical assimilation to obtain a smooth, continuous concentration field includes: A global scan of the concentration distribution field is performed to identify discontinuous boundaries where concentration values ​​abruptly change. By combining the effective concentration data of the surrounding grid, the discontinuous boundary is interpolated and corrected to obtain the corrected concentration distribution field; A global smoothing process is performed on the modified concentration distribution field to obtain a smooth continuous concentration field.

7. The multi-sensor collaborative air pollution monitoring method according to claim 1, characterized in that, If the smooth, continuous concentration field shows a migration trend, then the spatiotemporal changes of the migration path are tracked and a grid deformation control vector is generated to dynamically adjust the spatiotemporal grid of the prediction field, resulting in a dynamic concentration prediction field, including: Based on the smooth, continuous concentration field, the Lagrange particle tracking method is used to identify regions of significant concentration migration. Based on the regions of significant concentration migration, a spatiotemporal evolution manifold is constructed to deduce the spatiotemporal trajectory of the migration path; The spatial drift of the spatiotemporal change trajectory is calculated by combining the preset parameter update cycle, and a grid deformation control vector is generated. The predicted field spatiotemporal grid is stretched using the grid deformation control vector to construct an adaptive spatiotemporal grid system; The baseline concentration field, the optimized concentration distribution, and the sampled concentration interpolation data are mapped to the adaptive spatiotemporal grid system, and flux conservation interpolation calculations are performed to obtain a preliminary dynamic concentration prediction field.

8. The multi-sensor collaborative air pollution monitoring method according to claim 1, characterized in that, The process involves calculating a weighted baseline concentration distribution based on the dynamic concentration prediction field, constructing a gradient deviation field based on the weighted baseline concentration distribution and the sampled concentration interpolation data, and iteratively correcting the parameters of the preset prediction model based on the gradient deviation field to obtain a high-precision local concentration field. This includes: Based on the dynamic concentration prediction field, a weighted baseline concentration distribution is generated by combining the preset historical average concentration and spatial distance weights. The concentration difference between the weighted baseline concentration distribution and the corresponding sampled concentration interpolation data is calculated for each grid point to obtain the moving sampling residual set. Spatial interpolation is then performed on the moving sampling residual set to obtain the gradient bias field. The gradient correction vector is calculated based on the spatial distribution of the gradient deviation field. The gradient correction vector is then used to adjust the local diffusion coefficient and model decay factor of the preset prediction model to generate the corrected model variable matrix. The modified model variable matrix is ​​input into the concentration evolution equation to perform iterative calculations until the results converge, ultimately yielding a high-precision local concentration field.

9. The multi-sensor collaborative air pollution monitoring method according to claim 1, characterized in that, After obtaining the high-precision local concentration field, the process also includes: The concentration observations from multiple real-time sensors are acquired, and the instantaneous prediction error tensor is generated by differential operation in combination with the local concentration field. The Holt-Winters smoothing algorithm is applied to the instantaneous prediction error tensor to obtain a smoothed error curve. If the second derivative of the smoothed error curve is greater than a preset abrupt change threshold, a high-frequency resampling instruction is generated. The state observation vector of the Kalman filter is reconstructed according to the high-frequency resampling command, and the dynamic feedback gain matrix is ​​calculated based on the state observation vector. A closed-loop correction factor is generated based on the dynamic feedback gain matrix. The closed-loop correction factor is then incorporated into the concentration evolution iteration calculation of the preset prediction model at the next time step to obtain the corrected prediction model.

10. A multi-sensor collaborative air pollution monitoring system, characterized in that, include: The data acquisition and fusion module is used to collect base station baseline concentrations through fixed base stations and obtain mobile sampling concentrations through mobile devices to form a raw dataset. The data is then aligned using spatiotemporal inconsistency correction to obtain a fused dataset. The Kalman filter correction module is used to predict and correct the base station reference concentration based on the fused dataset using Kalman filtering, adjust the measurement update frequency according to the correction result, obtain an optimized concentration sequence, and map it to a grid to obtain the reference concentration distribution. The concentration field prediction and optimization module is used to map the baseline concentration distribution to three-dimensional spatiotemporal space to obtain the baseline concentration field and input it into the preset prediction model to generate the prediction field. It evaluates the concentration error within the range. If the evaluation result exceeds the preset accuracy threshold, it supplements the original dataset by smoothing historical data and Kriging interpolation to obtain the optimized concentration distribution. The refined construction module is used to extract the concentration diffusion trend from the optimized concentration distribution, generate a baseline concentration weight distribution, obtain sampled concentration interpolation data, and combine the baseline concentration weight distribution with a grid resolution adjustment method to construct a refined concentration distribution field. The boundary correction processing module is used to correct the boundary of the concentration distribution field, identify discontinuous boundaries and perform numerical assimilation processing to obtain a smooth continuous concentration field. The dynamic migration prediction module is used to track the spatiotemporal changes of the migration path and generate a grid deformation control vector if the smooth continuous concentration field shows a migration trend, and dynamically adjust the spatiotemporal grid of the prediction field to obtain a dynamic concentration prediction field. The final output module is used to calculate the weighted baseline concentration distribution based on the dynamic concentration prediction field, construct a gradient deviation field based on the weighted baseline concentration distribution and the sampled concentration interpolation data, and iteratively correct the parameters of the preset prediction model based on the gradient deviation field to obtain a high-precision local concentration field.

Images