Road traffic pollution diffusion rapid simulation method based on physical structured proxy model

Abstract

The invention belongs to the technical field of near-roadside pollution diffusion simulation, and particularly relates to a road traffic pollution diffusion rapid simulation method based on a physical structured proxy model. Comprising the following steps: S1, meteorological characteristic engineering and sampling; s2, road network standardization decomposition: carrying out topology discretization decomposition on the complex urban road network into standardized unit line sources which can be independently calculated; s3, reconstructing spatial geometric features; s4, physical mechanism dual partitioning: constructing an orthogonal physical mechanism partitioning system according to the atmospheric stability level and the relative position of the receptor point; s5, proxy model training: respectively training an XGBoost regression model for each physical partition; and S6, area field rapid reconstruction: receiving real-time boundary conditions for parallel inference, and rapidly reconstructing an area pollutant concentration field through linear superposition. On the premise of ensuring the simulation precision highly consistent with that of a full-physical model, the calculation efficiency is remarkably improved, and the method is suitable for real-time air quality simulation of a large-scale urban road network.

Description

Technical Field

[0001] This invention belongs to the field of near-roadside pollution diffusion simulation technology, and relates to a rapid simulation method for road traffic pollution diffusion based on a physically structured surrogate model. Background Technology

[0002] Road traffic emissions are a major source of near-surface air pollution in modern cities. These emissions are low-lying and located close to areas of high population density, posing a significant threat to public health. Unlike elevated industrial sources, the concentration field of pollutants generated by road traffic exhibits strong spatial heterogeneity, typically displaying a non-linear, rapidly decaying gradient within a 150-200 meter radius around the road. To accurately assess exposure risks at such fine scales, high-resolution diffusion models are essential to capture these microscale concentration changes.

[0003] Currently, the AERMOD diffusion model system recommended by the U.S. Environmental Protection Agency (EPA) and its built-in RLINE line source algorithm are the mainstream regulatory tools for traffic environmental impact assessments and compliance checks, and are also widely used in academic research. The RLINE algorithm employs a rigorous adaptive numerical integration scheme, dynamically discretizing continuous line sources into point sources to ensure the physical consistency of near-road concentration gradient simulations under different wind directions and atmospheric stability conditions. However, existing RLINE algorithms face insurmountable computational efficiency bottlenecks when applied to regional-scale or long-term simulations. In practical city-level applications, involving thousands of kilometers of road networks, tens of thousands of receiving points, and 8760 hours of meteorological conditions annually, the combined computational workload of source-receiving point-hourly conditions explodes (typically reaching 10). 11 Up to 10 12 (On a scale of magnitude). Because RLINE requires complex iterative integration calculations for each combination, its computational cost increases superlinearly with the size of the road network and the density of receiving points. For example, calculating the annual PM2.5 concentration for a medium-sized metropolitan area... 2.5 Hourly simulations often require months or even longer of computation time. This enormous computational overhead makes it difficult to apply physical mechanism-based regulatory models to large-scale regional exposure assessments, real-time traffic policy analysis, or emergency responses to sudden pollution events. While existing simplified empirical models are computationally fast, they often sacrifice fidelity to physical mechanisms under complex meteorological conditions, making them unsuitable for high-precision assessments.

[0004] Therefore, given the contradiction between physical fidelity and computational speed in existing regulatory models, it is crucial to find a method that can both preserve the rigorous physical mechanism of RLINE and significantly reduce computational complexity through surrogate model technology, thereby achieving rapid and high-resolution simulation of road traffic pollution at the urban scale. Summary of the Invention

[0005] The technical problem this invention aims to solve is to address the issue that existing regulatory-grade air quality models (AERMOD-RLINE) suffer from excessive computational load and long processing time due to complex numerical integration mechanisms when simulating regional-scale, long-term road traffic pollution, thus failing to meet the requirements for high-resolution exposure assessment and real-time emergency response. This invention provides a rapid simulation method for road traffic pollution diffusion based on a physically structured proxy model.

[0006] Technical solution of the present invention: A rapid simulation method for road traffic pollution diffusion based on a physically structured surrogate model includes the following steps: S1. Meteorological Feature Engineering and Sampling: Multicollinearity diagnosis of meteorological parameters was performed based on Pearson correlation coefficient and variance expansion factor (VIF), and key features driving diffusion (including wind speed, wind direction, Moning-Obukhov length, sensible heat flux, and mixing layer height) were screened using SHAP attribution analysis; historical meteorological data were stratified and sampled using K-means clustering algorithm to construct a representative meteorological sample set covering multiple boundary layer features; S2. Road Network Standardization Decomposition: Based on the superposition property of the Gaussian line source model, the complex urban road network topology is decomposed into standardized unit line sources, transforming the multi-source road network simulation problem into a single-source independent computation problem; S3. Spatial Geometric Feature Reconstruction: For each unit line source, construct a wind-rotating coordinate system to convert the absolute geographic coordinates of the receiver point into a wind direction projection distance (x) relative to the line source. ^ ) and vertical distance from the crosswind direction (y ^ A non-uniform spatial sampling strategy was used to generate a receptor point sample set containing a near-source high gradient region and a far-end background region. The rotated actual coordinates were further transformed to the wind direction coordinate system for feature engineering expression. S4. Physical mechanism dual zoning: Based on the Moning-Obukhov length (L), the atmospheric turbulence state is divided into different stability levels (such as very stable, stable, neutral, unstable, and very unstable), and based on the relative position of the receptor point, it is divided into the plume-dominant region (downwind) and the meandering-dominant region (upwind / near source), forming an orthogonal physical zoning system. S5. Proxy Model Training: High-precision true concentration values ​​are generated using the RLINE full physical model, which combines the representative meteorological sample set and the receptor point sample set. For each physical partition, an XGBoost nonlinear regression model is constructed, with meteorological and spatial geometric features as input, and the model is trained to learn the response surface of physical diffusion. S6. Rapid Reconstruction of Regional Field: Receives road network, emission, and meteorological data of the area to be simulated, performs parallel inference through a group of trained surrogate models, maps the standardized prediction results back to the actual geographic space, and overlays them to generate a regional concentration field.

[0007] This invention achieves rapid and high-precision simulation of road traffic pollution through the following logical process: First, in the data preparation stage (S1), this invention scientifically selects key physical variables driving atmospheric diffusion through Pearson correlation analysis, VIF statistical diagnosis, and SHAP physical attribution analysis. Cluster sampling is used to ensure that the training data fully covers the entire spectrum of meteorological conditions, from extreme stability to strong convection, laying a solid foundation of physical data. Second, through the "line-source decomposition" strategy (S2), the real road network with infinite geometric possibilities is transformed into standardized feature inputs with finite dimensions, solving the problem that traditional machine learning models struggle to handle complex topological road networks and achieving geometric generalization for different urban road network layouts. In the modeling stage (S3-S5), this invention differs from traditional black-box models by adopting a "physically structured" architecture. Geographic coordinates are converted into a relative coordinate system that rotates with the wind, and a layered sub-model system is established based on atmospheric physical laws (such as wind direction and convective instability, neutrality, and stability stratification). This design ensures that the surrogate model focuses on learning the universal geometric law of "pollutant decay with distance" and strictly follows the laws of physical diffusion, maintaining the physical plausibility of predictions, especially under extreme weather conditions. In the application phase (S6), the efficient reasoning capability of the tree model is used to replace the time-consuming adaptive numerical integration calculation in the original physical model, thereby achieving an order-of-magnitude improvement in computational efficiency while maintaining physical fidelity.

[0008] Beneficial effects Compared with the prior art, the present invention has the following significant advantages: Order-of-magnitude improvement in computational efficiency: Compared to the traditional AERMOD-RLINE regulatory model, this invention achieves more than 15 times the computational speedup when dealing with regional road networks (such as those larger than 1.5 square kilometers) by replacing complex numerical integration with efficient tree model inference (e.g., controlling the simulation time per hour to within 0.6 seconds), making high-resolution hourly simulations possible throughout the year.

[0009] Excellent physical fidelity and robustness: This invention employs a physical partitioning-based modeling strategy, avoiding the failure of traditional statistical models under extreme conditions. Validation results show that this method achieves a high degree of accuracy (R0) in predicting near-path pollutant concentrations, consistent with the original physical model. 2 >0.98, MAE <1.0 μg / m 3 This effectively resolves the contradiction between "computation speed" and "physical accuracy".

[0010] Powerful geometric generalization ability: Through the standardized line source decomposition technique (S1), this model possesses the ability to transfer across regions. The model only needs to be trained once on a specific dataset (such as Texas data) and can be directly applied to completely different geographical regions (such as California or Shanghai, China), without the need for large-scale physical modeling and training for each new city, greatly reducing the application threshold.

[0011] Support for large-scale, refined exposure assessment: This invention breaks through the bottleneck of computing power, making it possible to conduct long-term, high spatial resolution traffic-related air pollution exposure assessments at the scale of large cities, providing refined data support for epidemiological research, urban planning, and traffic policy formulation. Attached Figure Description

[0012] Figure 1 This is the main flow of a rapid simulation method for road traffic pollution diffusion based on a physically structured proxy model, as described in this invention. Figure 2 This is a heatmap of the Pearson correlation coefficient matrix of meteorological input parameters under stable and unstable conditions according to embodiments of the present invention. Figure 3 This is a schematic diagram showing the marginal contribution results of various meteorological variables under stable and unstable conditions based on the XGBoost-SHAP framework in an embodiment of the present invention. Figure 4 The following are examples of receiver point layout near a unit line source and corresponding receiver point layout schemes under special wind directions in embodiments of the present invention ((a) receiver point layout scheme for a 10 m unit line source, (b) receiver point layout scheme with a wind direction of 60° east of north as an example). Figure 5 This is a heatmap of the RMSE test results of each sub-model in the embodiments of the present invention on the validation set; Figure 6 This is a schematic diagram of the road network and receiver point layout scheme for a test scenario according to an embodiment of the present invention; Figure 7 This is a schematic diagram of the overall prediction performance of the surrogate model under three typical surface roughness conditions in an embodiment of the present invention. Detailed Implementation

[0013] To provide a clearer understanding of the technical features, objectives, and effects of the present invention, specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0014] like Figure 1 As shown, a rapid simulation method for road traffic pollution diffusion based on a physically structured surrogate model includes the following steps: S1. Meteorological Feature Engineering and Sampling: Key meteorological factors are screened based on statistical diagnosis and SHAP attribution, and a representative meteorological sample set is constructed using clustering algorithms; S2. Road Network Standardization Decomposition: Discretizes and decomposes the complex urban road network topology into independently computable standardized unit line sources; S3. Spatial geometric feature reconstruction: Construct a wind-rotating coordinate system to generate a non-uniform receptor point space containing near-source high gradient features, and use the features in this coordinate system as input for model training; S4. Dual physical mechanism zoning: An orthogonal physical mechanism zoning system is constructed based on atmospheric stability levels and the relative positions of receptor points; S5. Proxy Model Training: Train XGBoost regression models for each physical partition separately to learn the nonlinear response surface of the physical model; S6. Rapid reconstruction of regional field: Receives real-time boundary conditions for parallel inference and rapidly reconstructs the regional pollutant concentration field through linear superposition.

[0015] Furthermore, in step S1, the meteorological feature engineering and sampling aim to eliminate redundant variable interference and construct a physically representative training set. The specific steps are as follows: Step 1.1 To meet the meteorological input requirements of the RLINE model, construct a differentiated feature system based on atmospheric stability.

[0016] In this system, wind direction (WDIR), as a key variable coupling the meteorological field and the spatial geometry, is set as a mandatory retention feature and is not involved in subsequent elimination screening. Other meteorological parameters are configured differently based on the stability state: under stable conditions, wind direction (WDIR), sensible heat flux (H), surface friction velocity (USTAR), and mechanical boundary layer mixing height (MixHGT) are selected. M The basic variable set consists of the Moning-Obukhov length (L), Bowen ratio (B0), albedo, wind speed (WSPD), temperature (Temp), precipitation (PrecipAmt), relative humidity (RH), local air pressure (Pressure), and cloud cover (CloudCover). Under unstable conditions, the convective velocity scale (WSTAR), potential temperature gradient (ThetaGrad), and convective mixing height (MixHGT) are further introduced on this basis. C Three variables are used to reflect the atmospheric turbulence structure when convection is dominant.

[0017] Sampling was conducted under both stable and unstable conditions, and statistical diagnosis was performed using a dual approach: the Pearson correlation coefficient (reference indicator) and the variance inflation factor (VIF) (primary criterion).

[0018] The formula for calculating the Pearson correlation coefficient and its explanation are as follows: in, Representing meteorological variables x With meteorological variables y The Pearson correlation coefficient between them; n The total length of the time series of meteorological data samples; and They are respectively t = i The observed values ​​of two meteorological variables at a given time; and meteorological variables x and y Arithmetic mean over a time series.

[0019] This invention is based on It serves as a reference for determining whether there is a strong linear correlation between variables and is used to help identify specific variable pairs that lead to high VIF values.

[0020] The formula and explanation for calculating VIF are as follows: in, It is the first k The determination coefficients were obtained by performing linear regression analysis with one meteorological variable as the dependent variable and all other meteorological variables as independent variables. This invention uses... A warning threshold is set to determine if there is strong collinearity among variables.

[0021] Specifically, the screening process employs a "physically constrained VIF iterative elimination strategy": In each round, the VIF values ​​of all variables in the current feature set are first calculated, and candidate variables with the largest values ​​exceeding the warning threshold (VIF>10) are identified and their physical attributes are determined. If the variable has physical overlap and functional redundancy with other retained variables (e.g., surface friction velocity and wind speed both represent momentum magnitude, and the former can be deduced from the latter), it is eliminated. Conversely, if the variable represents an independent physical dimension (e.g., wind speed represents momentum transport, and sensible heat flux represents heat exchange), even if the VIF value is high, it is still determined that it cannot be eliminated and is forcibly retained to ensure the integrity of the model's description of advection and diffusion physical processes, and the second highest value variable is examined instead. After a variable is eliminated, the remaining variables are used to update and calculate the VIF values ​​for a new round, and the above process is repeated until the VIF values ​​of all non-forced retained variables meet the threshold requirements.

[0022] Based on this strategy, the initial meteorological feature set obtained after preliminary screening of the input variables of the aforementioned RLINE meteorology was obtained. as follows: The characteristics after initial screening under stable conditions are sensible heat flux (H), Moning-Obukhov length (L), Bowen ratio (B0), albedo, wind speed (WSPD), precipitation (PrecipAmt), relative humidity (RH), and cloud cover (CloudCover). The characteristics of the initial screening under unstable conditions are sensible heat flux (H), potential temperature gradient (ThetaGrad), and convective mixing height (MixHGT). C The parameters include: Morning-Obukhov length (L), Bowen ratio (B0), wind speed (WSPD), precipitation (PrecipAmt), relative humidity (RH), and cloud cover (CloudCover).

[0023] Step 1.2 Construct standardized line source scenarios under stable and unstable weather conditions, fix a single receptor point (10m directly in front of the line source), and use the baseline model RLINE to simulate the concentration at the receptor point.

[0024] An XGBoost model was constructed to replace the RLINE model for calculating the receptor concentration at the sampling point: using the initial meteorological feature set obtained in step 1.1. Using the receptor concentration simulated by the baseline model RLINE as the input, the XGBoost model is trained to obtain a pre-trained XGBoost model under the normalized line source scenario.

[0025] The SHAP framework is used to quantify the marginal contributions of each variable and identify the key physical factors driving diffusion. Features i The formula for calculating the SHAP value is as follows: in, f For pre-trained standardized XGBoost models in line source scenarios; N For the input feature set (i.e. ); S A subset of features; The set difference operation refers to the operation of subtracting elements from the entire set. N Mid-removal features i The feature set formed afterward; Representing a subset relation, the summation symbol traverses... All possible feature subsets S ; The cardinality of a set is the number of features contained in the set within the double vertical lines. This represents the set union operation. Indicates target features i Add to subset S New feature combinations formed in the process; For the model in feature subset SThe predicted value is below.

[0026] The global feature importance ranking is determined based on the absolute value of the SHAP values ​​of each feature, and key features are selected according to their contribution to construct feature vectors. The results are as follows: Under stable conditions, the feature vector formed by the key features ultimately selected in this invention is: = [ WSPD, WDIR, L,H ], corresponding to wind speed, wind direction, Moning-Obukhov length, and sensible heat flux, respectively; the key feature vector under unstable conditions is = [ WSPD,WDIR, L,H, MixHGT C ],in MixHGT C This is the height for convective mixing.

[0027] Step 1.3 To further ensure the efficiency and generalization ability of the surrogate model during training, the original meteorological dataset composed of the above-selected features is subjected to stratified representative sampling.

[0028] Under each meteorological condition, the hourly observation data were grouped using the K-means clustering method, and the sum of squared errors within each cluster was minimized. The optimal cluster center is determined using the following formula: Where K is the number of clusters, For the k-th meteorological cluster, Let be the centroid of the cluster. Then, select the sample closest to the cluster center from each cluster as a representative meteorological instance to construct the final representative meteorological set, which serves as the meteorological condition input for subsequent surrogate model training.

[0029] Furthermore, in step S2, the road network standardization decomposition aims to transform the complex urban road geometry and topology into independently computable standardized units based on the linear superposition principle of the Gaussian model. The specific steps are as follows: Step 2.1 First, based on the node coordinate sequence of the road network topology data, calculate the Euclidean distance between adjacent nodes and generate the cumulative arc length sequence S of the road; then, according to the set standard discrete interval... (In this embodiment, the value is set to 10 m) Calculate the number of micro-segments required for the division, and determine the distance from the center of each micro-segment along the road direction. d m Finally, a linear interpolation algorithm is used to map the one-dimensional center distance back to the two-dimensional plane, obtaining the geometric center coordinates of each standardized micro-element. x m , y m The calculation formula is as follows: The recursive formula for the cumulative arc length sequence S of the road is: No. The center position of each infinitesimal segment d m for: Determine the index j Make The center coordinates are calculated as follows: in, The total length of the road; n The total number of road nodes; x j ， y j ) is the first j The coordinates of the original road nodes.

[0030] Step 2.2 For each standardized micro-segment, first locate its original road segment index by retrieving the cumulative arc length sequence. k The vector direction angle in radians is calculated based on the difference in coordinates between the start and end points of the road segment. Then, the radian value in the mathematical polar coordinate system is converted into a geographical azimuth angle with due south as 0 degrees (clockwise direction). Finally, considering the axisymmetric characteristics of line source diffusion, the angle is normalized to a modulus of 180 degrees to obtain a standardized road orientation angle. The calculation formula is as follows: Vector direction angle calculate: Geographic azimuth conversion and normalization: in,( , ) is the index k The starting coordinates of the original road segment It is the vector direction angle (in radians) in the mathematical polar coordinate system. It is the geographical azimuth. This represents the floor function, used to remove redundant half-cycles to ensure... .

[0031] Step 2.3 Based on all-time traffic emission data, firstly, using the road ID as the association key, traverse each time step. tThe static geometric micro-element set containing spatial coordinates and azimuth angles generated by discretization is matched with the dynamic emission intensity at that moment to construct a standardized spatiotemporal input tensor containing timestamp, micro-element ID, center coordinates, direction angle and instantaneous emission intensity.

[0032] Furthermore, in step S3, the spatial geometric feature reconstruction aims to capture concentration gradient features through non-uniform sampling and construct a physically consistent model input space using coordinate transformation. The specific steps are as follows: Step 3.1 Based on the physical attenuation law of pollutant concentration field with large gradient in the near-source region and gradual gradient in the far-source region, a non-uniformly distributed initial set of receptor points is constructed under the reference wind direction (set as 270° westerly wind in this invention); with the midpoint of the unit line source as the origin, a continuous segment [-100m, 2000m] from the upwind background area to the downwind dilution area is covered on the dominant wind direction axis. A variable step size distribution strategy of dense near-source and sparse far-end is adopted to fully analyze the concentration characteristics of the high gradient area while ensuring computational efficiency.

[0033] Step 3.2 Based on the meteorological wind direction data at each moment, construct a two-dimensional rotation transformation matrix to map the initial sampling points under the reference wind direction to the actual geographic space, thereby obtaining the coordinates of the actual receiver point. (x r , y r During this process, a threshold for the geometric buffer zone on both sides of the line source is set to eliminate singular points falling within the buffer zone, thus avoiding instability in the numerical calculation. The rotation transformation calculation formula is as follows: in, The initial coordinates are under the reference wind direction; The angle difference between the wind direction and the reference wind direction at this moment; The coordinates of the source center of the unit line are given.

[0034] Step 3.3 To eliminate the geometric differences caused by different wind directions and improve the generalization ability of the model, the actual receiver point coordinates are transformed into a wind-rotating coordinate system for feature representation; the receiver point position is orthogonally decoupled into the projected distance along the wind direction. Vertical distance from crosswind direction These characteristics represent the advectional transport attenuation of pollutants along the streamline and the turbulent diffusion structure in the lateral direction, respectively. This allows the model input features to be independent of absolute geographic coordinates, focusing on learning universal physical diffusion laws. The calculation formula is as follows: in, The flow field vector angle (in radians) corresponding to the measured wind direction; This represents the longitudinal distance in the wind direction coordinate system; This represents the lateral distance in the wind direction coordinate system.

[0035] Furthermore, in step S4, the dual partitioning of the physical mechanism aims to decouple the global simulation problem into several sub-problems with consistent analytical forms based on the switching boundary of the computational mechanism within the RLINE model. The specific steps are as follows: Step 4.1 Based on the decisive influence of atmospheric turbulence on the vertical diffusion parameterization scheme, meteorological stability zones are constructed based on the Moning-Obukhov length (L); in the RLINE model, different vertical diffusion coefficients are used under stable (L>0) and unstable (L<0) conditions. The calculation formula and probability density function (PDF) are used. Therefore, this invention sets a physical threshold to classify meteorological samples into 5 meteorological stability categories (very stable, stable, neutral, unstable, and very unstable). The discrimination rules are as follows: Step 4.2 Based on the dominant region differences of the "plume-meandering" dual-component physical mechanism in the RLINE model, spatial geometric partitioning is constructed. The RLINE basic framework decouples the diffusion process into two fundamentally different physical processes: the "plume component," describing the conventional downwind Gaussian diffusion structure, and the "meandering component," describing the irregular oscillation of wind direction caused by low wind speeds or strong turbulence. Physically, the concentration in the upwind and near-source regions originates entirely from the lateral symmetrical diffusion caused by the meandering effect, while the downwind region is mainly dominated by the plume effect stretched along the streamline direction. Given the inherent physical heterogeneity of the concentration formation mechanism when the receiver point is upwind or downwind, this invention divides the space into an "upwind / meandering-dominated region" influenced only by random oscillations and a "downwind / plume-dominated region" controlled by advection transport, using the normal position perpendicular to the wind direction as the boundary. Partitioning Indicator Variables The definition is as follows: Step 4.3 Construct XGBoost sub-proxy models corresponding one-to-one with the RLINE physics solution module by orthogonally combining meteorological stability and spatial geometric location; for each subdomain (e.g., "stable-downwind"), due to its relationship with the underlying RLINE analytical expression (including diffusion parameters) , By maintaining consistency with physical response mechanisms (such as calculations), the XGBoost model can more efficiently capture the physical response patterns, thereby avoiding gradient abrupt changes and fitting errors that can occur when a single model crosses physical mechanisms (such as from unstable to stable, or from upwind to downwind).

[0036] Furthermore, in step S5, the surrogate model training aims to fit the input-output relationship of the physical model using machine learning algorithms. By learning the nonlinear response surface between physical parameters and the concentration field, a group of surrogate models capable of replacing RLINE for fast inference is constructed. The specific steps are as follows: Step 5.1 Constructs a hybrid feature vector containing physical meteorological attributes and relative geometric attributes, and selects Extreme Gradient Boosting Tree (XGBoost) as the core regression engine. Compared with deep neural networks such as Multilayer Perceptron (MLP), XGBoost utilizes histogram splitting and parallel computing mechanisms, exhibiting superior computational efficiency and interpretability when processing hundreds of millions of structured tabular data. At the feature construction level, to eliminate the negative impact of discontinuous wind direction angle values, the wind direction is decomposed into orthogonal wind direction sine and cosine components, and combined with wind speed, Moning-Obukhov length, sensible heat flux, and convective mixing layer height (the results after meteorological feature screening in Step S1) to constitute the full-process meteorological features. For the meteorological features in the neutral case, this invention sets when L... At 1000, its input characteristics are the same as in the steady state, while L At -1000, the situation is the same as the unstable case; simultaneously, combined with the windward projected distance generated in step S3... Vertical distance from crosswind direction This forms a complete model input space.

[0037] Step 5.2 Define an objective function that includes training loss and structural risk regularization to constrain model complexity while ensuring fitting accuracy; assuming the training set contains n samples, XGBoost approximates the real physical field by iteratively accumulating K decision trees. Its objective function... The differentiable convex loss function (mean squared error, MSE) measures the prediction error, and the regularization term for the complexity of the penalty tree model is used. The composition and calculation formula are as follows: in To predict concentration; The truth value calculated for RLINE; Let the mean squared error loss function be used. This represents the number of leaf nodes; γ is the L2 modulus square of the leaf node score; γ and λ are the corresponding penalty coefficients, and in this embodiment, γ and λ are 0.1 and 1 respectively; by introducing L1 / L2 regularization, the generalization ability of the model under unseen weather conditions is effectively improved.

[0038] Step 5.3 Finally, the optimal hyperparameter combination is determined using a strategy of "initial selection based on experience combined with local fine-tuning," maximizing the real-time inference speed for engineering deployment while ensuring fitting accuracy. To meet the high-frequency inference requirements of millions of receptor points, this invention specifically sets the key hyperparameters of the XGBoost model: the maximum tree depth (max_depth) of the decision tree is limited to 5 to significantly shorten the computation path of a single prediction and reduce latency; the learning rate (learning_rate) is set to 0.15 and the number of iterations (n_estimators) is set to 1000 to ensure sufficient residual fitting; and the subsample ratio (subsample) is set to 0.8 to enhance robustness. Through the above configuration, the trained surrogate model achieves the best balance between computational cost and prediction performance while maintaining physical consistency.

[0039] Furthermore, in step S6, the rapid reconstruction of the regional field aims to apply the offline-trained surrogate model to online inference, generating a high-resolution concentration distribution through parallel computing and superposition principles. The specific steps are as follows: The system receives real-time meteorological boundary conditions and dynamic road network emission inventories for the simulated period. For each receptor point within the target area grid, it iterates through all unit line sources standardized in step S2 to construct "source-receptor" feature pairs. Based on real-time wind direction data, it batch calculates the wind direction projection distance of the receptor point relative to each line source. Vertical distance from crosswind direction The stability level is determined by combining the atmospheric Moning-Obukhov length at the current moment, and each feature vector is matched and input into the corresponding physical sub-proxy model interface trained in step S5.

[0040] Each physical sub-proxy model performs parallel inference on the above feature batches to obtain the normalized source strength of each unit line source ( = 1) The predicted diffusion concentration at the acceptor point; finally, based on the linear additivity principle of the Gaussian model, the predicted value is multiplied by the actual emission intensity of each linear source at that moment. The total regional pollutant concentration field with high spatiotemporal resolution is reconstructed by summing the contributions from all line sources. The calculation formula is as follows: in, The total number of unit line sources affecting the receptor point within the region; For the first Each line source is in The actual emission source strength at any given time; In order to be in accordance with physical conditions k The corresponding sub-proxy model; For the first Each line source to receptor point The relative geometric and meteorological characteristic vectors.

[0041] Example Taking a 1.5×1.5 km square urban road network centered on the Yan'an Elevated Road Station in Jing'an District, Shanghai, as the research object, hourly meteorological observation data and hourly NO data of the road network from the Xuhui District Meteorological Station in Shanghai throughout 2023 were selected. X Using emission data as input, and relying on the rapid simulation method for road traffic pollution diffusion based on a physically structured surrogate model proposed in this invention, a regional concentration distribution field with high spatiotemporal resolution is generated. The specific implementation steps are as follows: S1. Using AERMOD-RLINE meteorological data regenerated from 256 towns in Texas, USA for five consecutive years (2016-2021) as input (containing meteorological data under three typical surface roughness conditions: open plains, general suburban towns, and dense urban centers), multicollinearity diagnosis and physical driving factor screening were performed, and stratified representative sampling was conducted.

[0042] Furthermore, in step S1, the meteorological feature engineering and sampling aim to eliminate redundant variable interference and construct a physically representative training set. The specific steps are as follows: (1) For the meteorological input parameters required by the RLINE model, a differentiated feature system is constructed based on atmospheric stability. In this system, wind direction (WDIR), as a key variable coupling the meteorological field and the spatial geometric field, is set as a mandatory feature and is not involved in subsequent elimination screening. The remaining meteorological parameters are configured differently according to the stability state, as shown in Table 1. Statistical diagnosis is performed using both Pearson correlation coefficient and variance inflation factor (VIF) under stable and unstable conditions. The calculation formula and explanation of the Pearson correlation coefficient are as follows: in, Representing meteorological variables x With meteorological variables y The Pearson correlation coefficient between them; n The total length of the time series of meteorological data samples; and They are respectively t = i The observed values ​​of two meteorological variables at a given time; and meteorological variables x and y The arithmetic mean over a time series. This invention specifies that... It serves as a reference for determining whether there is a strong linear correlation between variables and is used to help identify specific variable pairs that lead to high VIF values.

[0043] The formula and explanation for calculating VIF are as follows: in, It is the first k The determination coefficients were obtained by performing linear regression analysis with one meteorological variable as the dependent variable and all other meteorological variables as independent variables. This invention uses... To determine the warning threshold for strong collinearity among variables, the screening process employs a "physically constrained VIF iterative elimination strategy": In each round, the VIF values ​​of all variables in the current feature set are first calculated, identifying candidate variables with the highest values ​​exceeding the warning threshold (VIF>10), and their physical properties are then assessed. If this variable overlaps in physical meaning or has functional redundancy with other retained variables (e.g., surface friction velocity and wind speed both represent momentum magnitude, and the former can be deduced from the latter), it is eliminated. Conversely, if the variable represents an independent physical dimension (e.g., wind speed represents momentum transport, and sensible heat flux represents heat exchange), even with a high VIF value, it is deemed unremovable and forcibly retained to ensure the model's complete description of advection and diffusion processes, and the next highest-valued variable is considered. After a variable is eliminated, the remaining variables are used to update and calculate the VIF values ​​for a new round, repeating the above process until the VIF values ​​of all non-forced-retained variables meet the threshold requirements. Based on this strategy, using VIF as the primary criterion and Pearson correlation coefficient as a reference, the input variables for the aforementioned RLINE meteorological data are initially screened to obtain a preliminary set of meteorological features. : The characteristics after initial screening under stable conditions are sensible heat flux (H), Moning-Obukhov length (L), Bowen ratio (B0), albedo, wind speed (WSPD), precipitation (PrecipAmt), relative humidity (RH), and cloud cover (CloudCover). The characteristics of the initial screening under unstable conditions are sensible heat flux (H), potential temperature gradient (ThetaGrad), convective mixing height (MixHGTC), Moning-Obukhov length (L), Bowen ratio (B0), wind speed (WSPD), precipitation (PrecipAmt), relative humidity (RH), and cloud cover (CloudCover).

[0044] The Pearson correlation coefficient matrices under stable and unstable conditions are as follows: Figure 2 As shown, the calculation results of VIF under the two conditions are presented in Tables 2 and 3. Initial statistical diagnostics under steady-state conditions indicate that the sensible heat flux (H), surface friction velocity (USTAR), and mechanical mixing layer height (MixHGT) are significantly affected. MA significant linear correlation exists between sensible heat flux (H) and wind speed (WSPD). Although the variance expansion factor (VIF) values ​​of sensible heat flux (H) and wind speed (WSPD) are still relatively high after initial screening, given that they independently characterize the thermal and dynamic driving mechanisms of the atmospheric boundary layer, respectively, and their physical significance is irreplaceable, they are forcibly retained in this embodiment. Conversely, variables with redundant physical information are removed, and a variable set containing 8 features is finally constructed, as shown in Table 2. The VIF values ​​of all retained variables have fallen back to an acceptable range. Under unstable conditions, a strong negative correlation is observed between sensible heat flux (H) and the convective velocity scale (WSTAR). After screening and removing functionally redundant variables, a variable set including the convective mixing layer height (MixHGT) is obtained. C Nine final characteristics, including ( ). As shown in Table 3, the VIF values ​​of all retained variables are below 6, indicating that the multicollinearity problem has been effectively controlled.

[0045] Table 1 List of candidate meteorological features Table 2 VIF Analysis under Stable Meteorological Conditions Table 3. VIF Analysis under Unstable Meteorological Conditions (2) Construct standardized line source scenarios under stable and unstable conditions, fix a single receiver point (10m directly in front of the line source), and use the initial screening meteorological feature set obtained in step (1.1) Using the receptor concentration simulated by the baseline model RLINE as input and the concentration at the target site as output, the XGBoost model is trained. The SHAP framework is used to quantify the marginal contributions of each variable and identify the key physical factors driving diffusion. The formula for calculating the SHAP value of feature i is as follows: in, f For pre-trained models; N For the input feature set; S A subset of features; The set difference operation refers to the operation of subtracting elements from the entire set. N Mid-removal features i The feature set formed afterward; Representing a subset relation, the summation symbol traverses... All possible feature subsets S ; The cardinality of a set is the number of features contained in the set within the double vertical lines. This represents the set union operation. Indicates target features i Add to subsetS New feature combinations formed in the process; For the model in feature subset S The predicted value is below.

[0046] The global feature importance ranking is determined based on the absolute value of the SHAP value of each feature, and key features are selected according to their contribution to construct a feature vector. The key meteorological input features derived from a large number of meteorological data samples in this embodiment are as follows: Under stable conditions, the feature vector constructed from the key features selected in this invention is... = [ WSPD,WDIR, L, H ], corresponding to wind speed, wind direction, Moning-Obukhov length, and sensible heat flux, respectively; the key feature vector under unstable conditions is = [ WSPD,WDIR, L,H, MixHGT C ],in MixHGT C This is the height for convective mixing.

[0047] The marginal contributions of each variable under stable and unstable conditions based on the XGBoost-SHAP framework are as follows: Figure 3 As shown.

[0048] (3) To further ensure the efficiency and generalization ability of the surrogate model during training, this invention performs stratified representative sampling on the original meteorological dataset. Under each meteorological condition, the hourly observation data is grouped using the K-means clustering method, and the sum of squared errors within each cluster is minimized. The optimal cluster center is determined using the following formula: Where K is the number of clusters, and CK is the k-th meteorological cluster. Let K be the centroid of the cluster. Then, the sample closest to the cluster center is selected from each cluster as a representative meteorological instance. In this invention, K takes the value of 8760 under each type of meteorological condition. Finally, a representative meteorological set is constructed as the meteorological training set for subsequent surrogate model training.

[0049] S2. Standardized decomposition of the road network, taking a square urban road network with a radius of 1.5 × 1.5 km centered on the Yan'an Elevated Road Station in Jing'an District, Shanghai as an example, specifically includes the following steps: (1) First, based on the node coordinate sequence of the road network topology data, calculate the Euclidean distance between adjacent nodes and generate the cumulative arc length sequence S of the road; then, according to the set standard discrete interval (In this embodiment, the value is set to 10 m) Calculate the number of micro-segments required for the division, and determine the distance from the center of each micro-segment along the road direction. d mFinally, a linear interpolation algorithm is used to map the one-dimensional center distance back to the two-dimensional plane, obtaining the geometric center coordinates of each standardized micro-element. x m , y m The calculation formula is as follows: The recursive formula for the cumulative arc length sequence S of the road is: The center position of the m-th infinitesimal segment d m for: Determine the index j Make The center coordinates are calculated as follows: in, The total length of the road; n The total number of road nodes; x j ， y j ) is the first j The coordinates of the original road nodes.

[0050] (2) For each standardized micro-segment, the original road segment index to which it belongs is first located by retrieving the cumulative arc length sequence. k The vector direction angle in radians is calculated based on the difference in coordinates between the start and end points of the road segment. Then, the radian value in the mathematical polar coordinate system is converted into a geographic azimuth angle with true north as 0 degrees (clockwise). Finally, considering the axisymmetric characteristics of line source diffusion, the angle is normalized to a modulus of 180 degrees to obtain a standardized road orientation angle. The calculation formula is as follows: Vector direction angle calculate: Geographic azimuth conversion and normalization: in( , ) is the index k The starting coordinates of the original road segment It is the vector direction angle (in radians) in the mathematical polar coordinate system. It is the geographical azimuth. This represents the floor function, used to remove redundant half-cycles to ensure... .

[0051] (3) Based on the traffic emission data for all time periods, firstly, using the road ID as the association key, traverse each time step. t The static geometric micro-element set containing spatial coordinates and azimuth angles generated by discretization is matched with the dynamic emission intensity at that moment to construct a standardized spatiotemporal input tensor containing timestamp, micro-element ID, center coordinates, direction angle and instantaneous emission intensity.

[0052] The final standard road network unit information is shown in Table 4.

[0053] Table 4. Standardized decomposition of road network unit information (partial list) S3. With a length of 10 m, and y Taking a unit line source segment with parallel axes and the origin as its midpoint as an example, this step captures the concentration gradient features of the unit line source through non-uniform sampling and uses coordinate transformation to construct a physically consistent model input space. The specific steps are as follows: (1) Based on the physical attenuation law of pollutant concentration field with large gradient in the near-source region and gradual flattening in the far-source region, a non-uniformly distributed initial set of receptor points is constructed under the reference wind direction (set as 270° westerly wind in this invention). Taking the midpoint of the unit line source as the origin, a continuous segment [-100m, 2000m] from the upwind background area to the downwind dilution area is covered on the dominant wind axis. A variable step size distribution strategy with dense near-source and sparse far-source points is adopted. Along the wind direction dimension, a fine step size of 0.5~1.0m is used in the near-source region [-20m, 20m], and the step size gradually increases to 50~100m in the far-source region. In the vertical wind direction dimension, a small spacing of 0.5m is set for the high gradient region (±10m range) at the center of the plume, and increases to 10m towards the outer edge. This strategy fully analyzes the concentration characteristics of the high gradient region while ensuring computational efficiency. A schematic diagram of the receptor point layout is shown below. Figure 4 In (a), the arrangement of the receiver points for a 10 m unit line source is shown.

[0054] (2) Based on the meteorological wind direction data at each moment, a two-dimensional rotation transformation matrix is ​​constructed to map the initial sampling points under the reference wind direction to the actual geographic space, thereby obtaining the coordinates of the actual receiver point. (x r , y r During this process, a threshold for the geometric buffer zone on both sides of the line source is set to eliminate singular points falling within the buffer zone, thus avoiding instability in the numerical calculation. The rotation transformation calculation formula is as follows: in, The initial coordinates are under the reference wind direction; The angle difference between the wind direction and the reference wind direction at this moment; The coordinates of the source center are shown in the diagram. The receiver point arrangement scheme, taking a wind direction of 60° east of north as an example, is illustrated below. Figure 4 As shown in (b).

[0055] (3) In order to eliminate the geometric differences caused by different wind directions and improve the generalization ability of the model, the present invention transforms the actual receiver point coordinates to a wind-rotating coordinate system for feature engineering during subsequent model training. Specifically, the receiver point position is orthogonally decoupled into the projected distance along the wind direction. Vertical distance from crosswind direction These characteristics represent the advectional transport attenuation of pollutants along the streamline and the turbulent diffusion structure in the lateral direction, respectively. This allows the model input features to be independent of absolute geographic coordinates, focusing on learning universal physical diffusion laws. The calculation formula is as follows: in, The flow field vector angle (in radians) corresponding to the measured wind direction; This represents the longitudinal distance in the wind direction coordinate system; This represents the lateral distance in the wind direction coordinate system.

[0056] S4. After physically representing the spatial characteristics of the unit line source, this step further utilizes a physical mechanism for dual partitioning. Based on the switching boundary of the computational mechanism within the RLINE model, the simulation problem of the global unit line source is decoupled into several sub-problems with consistent analytical forms. The specific steps are as follows: (1) Based on the decisive influence of atmospheric turbulence on the parameterization scheme of vertical diffusion, meteorological stability zones are constructed based on the Monin-Obukhov length (L); in the RLINE model, drastically different vertical diffusion coefficients are used under stable (L>0) and unstable (L<0) conditions. The calculation formula and probability density function (PDF) are used. Therefore, this invention sets a physical threshold to classify meteorological samples into 5 meteorological stability categories (very stable, stable, neutral, unstable, and very unstable). The discrimination rules are as follows: (2) Based on the difference in the dominant regions of the "plume-meandering" dual-component physical mechanism in the RLINE model, spatial geometric partitioning is constructed; the RLINE basic framework decouples the diffusion process into two fundamentally different physical processes: the "plume component" describing the conventional downwind Gaussian diffusion structure, and the "meandering component" describing the irregular oscillation of wind direction caused by low wind speed or strong turbulence. Physically, the concentration in the upwind and near-source regions is entirely due to the lateral symmetrical diffusion caused by the meandering effect, while the downwind region is mainly dominated by the plume effect stretched along the streamline direction. Given that the concentration formation mechanism of the receptor point is inherently physically heterogeneous when it is upwind or downwind, this invention divides the space into the "upwind / meandering-dominated region" which is only affected by random oscillation and the "downwind / plume-dominated region" which is controlled by advection transport, using the normal position perpendicular to the wind direction as the boundary. Partitioning indicator variables The definition is as follows: (3) By orthogonally combining meteorological stability and spatial geometric location, a sub-proxy model corresponding one-to-one with the RLINE physical solution module is constructed; for each subdomain (e.g., "stable-downwind"), due to its relationship with the underlying RLINE analytical expression (including diffusion parameters) , To maintain consistency with the calculations, etc., the XGBoost model in subsequent steps can more efficiently capture the physical response patterns during training, thereby avoiding gradient abrupt changes and fitting errors caused by a single model when crossing physical mechanisms (such as from unstable to stable, or from upwind to downwind).

[0057] S5. This step selects the XGBoost model as the core algorithm model, fits the input-output relationship of the RLINE model in each physical partition through appropriate hyperparameter configuration, and constructs a surrogate model group that can replace RLINE for fast inference by learning the nonlinear response surface between meteorological parameters, spatial parameters and concentration fields. The specific steps are as follows: (1) Construct a hybrid feature vector containing physical meteorological attributes and relative geometric attributes, and select Extreme Gradient Boosting Tree (XGBoost) as the core regression engine; compared with deep neural networks such as Multilayer Perceptron (MLP), XGBoost utilizes histogram splitting and parallel computing mechanisms, which has better computational efficiency and interpretability when processing hundreds of millions of structured tabular data. At the feature construction level, in order to eliminate the negative impact of discontinuous wind direction angle values, the wind direction is decomposed into orthogonal wind direction sine and wind direction cosine components, and combined with wind speed, Moning-Obukhov length, sensible heat flux and convective mixing layer height (the result after meteorological feature screening in step S1) to form the meteorological features of the whole process. Among them, for the meteorological features of the neutral case, this invention sets when L At 1000, its input characteristics are the same as in the steady state, while L At -1000, the situation is the same as the unstable case; simultaneously, combined with the windward projected distance generated in step S3... Vertical distance from crosswind direction This forms a complete model input space.

[0058] (2) Define an objective function that includes training loss and structural risk regularization to constrain model complexity while ensuring fitting accuracy; assuming the training set contains n samples, XGBoost approximates the real physical field by iteratively accumulating K decision trees. Its objective function... The differentiable convex loss function (mean squared error, MSE) measures the prediction error, and the regularization term for the complexity of the penalty tree model is used. The composition and calculation formula are as follows: in To predict concentration; The truth value calculated for RLINE; The mean squared error loss function during training; This represents the number of leaf nodes; γ is the L2 modulus square of the leaf node score; γ and λ are the corresponding penalty coefficients, and in the embodiments of the present invention, the values ​​of γ and λ are 0.1 and 1, respectively; by introducing L1 / L2 regularization, the generalization ability of the model under unseen meteorological conditions is effectively improved.

[0059] (3) The optimal hyperparameter combination is determined by adopting a strategy of "initial selection based on experience combined with local fine-tuning" to maximize the real-time inference speed of engineering deployment while ensuring fitting accuracy. In order to meet the high-frequency inference requirements of millions of receptor points, the key hyperparameters of the XGBoost model are set in a targeted manner: the maximum tree depth (max_depth) of the decision tree is limited to 5 to significantly shorten the computation path of a single prediction and reduce latency; the learning rate (learning_rate) is set to 0.15 and the number of iterations (n_estimators) is set to 1000 to ensure sufficient residual fitting; and the subsample ratio (subsample) is set to 0.8 to enhance robustness. Through the above configuration, the surrogate model generated by training achieves the best balance between computational cost and prediction performance while maintaining physical consistency.

[0060] Finally, the RMSE test results of the sub-models on the validation set under three typical surface roughness conditions (representing open plains, general suburban towns, and dense urban centers, respectively) are as follows: Figure 5As shown, surface roughness Z0=0.05 represents an open plain area, Z0=0.5 represents a typical suburban town area, and Z0=1 represents a densely populated urban center area; VS, S, N1, N2, U, and VU represent very stable (…). ),Stablize( ), Neutral case 1 (L) 1000), Neutral Case 2 (L) -1000), unstable ( ), very unstable ( (meteorological conditions).

[0061] S6. Apply the sub-model group trained in step S5 to the actual road network and meteorological scenario. Since the meteorological data used for training the sub-model comes from Texas, USA, and the test scenario of this application example is the urban road network around Yan'an West Road in Jing'an District, Shanghai, and uses the meteorological data of Shanghai in 2023, it can also be regarded as a stress test of the model's generalization performance.

[0062] Based on the test road network and meteorological scenario in S1, uniformly spaced receptor points with a 50 m background interval were constructed (e.g., Figure 6 (as shown by the green dot in the middle) and the receptor points densified on both sides of the road (such as...) Figure 6 (As shown by the blue dot in the middle), the specific scene layout is as follows: Figure 6 As shown, the red area represents the isolation zone on both sides of the road, where no receptor points are deployed. For each receptor point within the target area grid, all unit line sources standardized in step S2 are traversed to construct a massive number of "source-receptor" feature pairs. Based on real-time wind direction data, the wind direction projection distance of the receptor point relative to each line source is calculated in batches. Vertical distance from crosswind direction The stability level is determined by combining the atmospheric Moning-Obukhov length at the current moment, and each feature vector is matched and input into the corresponding physical sub-proxy model interface trained in step S5.

[0063] Parallel inference is performed on the above feature batches using computer matrix operations to obtain the normalized source strength of each unit line source. = 1) The predicted diffusion concentration at the acceptor point; finally, based on the linear additivity principle of the Gaussian model, the predicted value is multiplied by the actual emission intensity of each linear source at that moment. The total regional pollutant concentration field with high spatiotemporal resolution is reconstructed by summing the contributions from all line sources. The calculation formula is as follows: in, The total number of unit line sources affecting the receptor point within the region; For the first Each line source is in The actual emission source strength at any given time; In order to be in accordance with physical conditions k The corresponding sub-proxy model; For the first Each line source to receptor point The relative geometric and meteorological characteristic vectors.

[0064] The overall predictive performance statistics of the model under three typical surface roughness conditions are shown in Table 5 and Figure 7 As shown, macroscopic statistical results indicate that the surrogate model of this invention exhibits excellent fitting ability under different landforms, and its coefficient of determination R with the benchmark RLINE model is [value missing]. 2 All exceeded 0.98, and the mean absolute error (MAE) was strictly controlled to 1.0 μg / m 3 Within. For example Figure 7 As shown, the surrogate model accurately captures the spatial distribution characteristics of the annual average hourly concentration, confirming the model's high accuracy performance. Tables 6 to 8 further list hourly prediction error samples under three typical scenarios. Comparative verification shows that the point-to-point prediction error of the surrogate model remains at an extremely low level. Combining the statistical data in Table 5 with the sample verification results in Tables 6-8, it can be seen that the model successfully balances global statistical robustness with local physical consistency, achieving high-precision reproduction of the pollutant diffusion process.

[0065] Table 5 Overall Predictive Performance Statistics of the Model Table 6. Details of concentration errors between the receptor-by-receptor surrogate model and RLINE in typical plains scenarios (partial list only) Table 7. Concentration error details between the receptor-per-point surrogate model and RLINE in typical suburban town scenarios (partial list only). Table 8. Details of concentration errors between the receptor-by-receptor surrogate model and RLINE in high-density central urban area scenarios (partial list only). Furthermore, to quantitatively evaluate the efficiency of practical engineering applications, this embodiment conducted a 24-hour performance comparison test based on the refined road network of Yan'an West Road in Shanghai (1,045 original roads discretized into 2,911 computational units). The test results are shown in Table 8. Under the single-process environment of the unified hardware platform of Intel Core i5-14600KF @ 5.30 GHz, compared with the 213 seconds (8.875 seconds / simulation hour) consumed by the AERMOD-RLINE benchmark model due to complex numerical integration, the PS-XGB-RLINE proxy model proposed in this invention only took 13.5 seconds (0.563 seconds / simulation hour), achieving a computational speedup of approximately 15.7 times. This significant improvement is attributed to the model's transformation of high-cost physics solutions into efficient vectorized matrix operations and tree model inference. Even with the source list refinement leading to a nearly three-fold increase in computational load, it still effectively breaks through the traditional computing power bottleneck. Theoretically, it can complete the full-process simulation of 8,760 hours throughout the year in just 1.37 hours. While maintaining high-precision prediction results, it fully meets the stringent timeliness requirements of large-scale long-term exposure assessment and real-time environmental emergency response.

[0066] Table 9 Comparison of computational efficiency between the proxy model and the baseline model RLINE in the same scenario The above embodiments are intended to illustrate the present invention in detail; however, the present invention is not limited to this single example. Any modifications, alterations, substitutions, combinations, simplifications, etc., made to the embodiments without departing from the spirit and principle of the present invention should be considered equivalent substitutions and included within the protection scope of the present invention.

Claims

1. A method for fast simulation of road traffic pollution dispersion based on physically structured agent models, characterized in that, The method comprises the following steps: S1. Meteorological feature engineering and sampling: based on the Pearson correlation coefficient and the variance inflation factor VIF, the meteorological parameters are subjected to multiple collinearity diagnosis, and key features driving diffusion are screened out by combining SHAP attribution analysis; hierarchical sampling is performed on historical meteorological data by using a clustering algorithm to construct a representative meteorological sample set covering multiple boundary layer features; S2. Road network standardization decomposition: based on the superposition of the Gaussian line source model, the city road network topology is decomposed into standardized unit line sources, and the multi-source road network simulation problem is converted into a single-source independent calculation problem; S3. Spatial geometric feature reconstruction: for each unit line source, a wind-rotating coordinate system is constructed, the geographical absolute coordinates of the receptor points are converted into wind direction projection distances and crosswind vertical distances relative to the line source, a non-uniform spatial sampling strategy is used to generate a receptor point sample set containing a near-source high-gradient area and a far-end background area, and the rotated actual coordinates are further converted to wind direction coordinates for feature engineering expression; S4. Dual partition of physical mechanism: the atmospheric turbulence state is divided into different stability levels according to the Monin-Obukhov length, and the receptor points are divided into plume-dominated areas and meandering-dominated areas according to their relative positions, forming an orthogonal physical partition system; S5. Training of proxy model: high-precision concentration true values are generated by using the RLINE full-physical model under the combination of the representative meteorological sample set and the receptor point sample set, and XGBoost nonlinear regression models are constructed for each physical partition, inputting meteorological and spatial geometric features, and training the model to learn the response surface of physical diffusion; S6. Rapid reconstruction of regional field: receiving the road network, emission and meteorological data of the region to be simulated, performing parallel inference by using the trained proxy model group, mapping the standardized prediction results back to the actual geographical space, and superimposing to generate a regional concentration field.

2. The method of claim 1, wherein the method comprises: In step S1, the specific method of the meteorological feature engineering comprises: The specific steps are as follows: Step 1.1 For the meteorological input parameters required by the RLINE model, the Pearson correlation coefficient and the variance inflation factor VIF are used to statistically diagnose the correlation between the stable and unstable cases, and the redundant variables are preliminarily removed; the obtained preliminary screening meteorological feature set ; Step 1.2: Construct a standardized line source scenario under stable and unstable conditions, quantify the marginal contribution of each variable using the SHAP framework, and identify key physical factors driving diffusion; Step 1.3: Hierarchical representative sampling of the original meteorological data set; in each type of meteorological condition, K-means clustering method is used to group the hourly observation data, the optimal cluster center is determined by minimizing the sum of squared errors within the cluster, and the sample closest to the cluster center is selected as the representative meteorological instance to ensure the efficiency and generalization ability of the proxy model training.

3. The method of claim 2, wherein the method comprises: In step 1.1, Select the basic variable set among the meteorological input parameters required by the RLINE model: under stable conditions, select wind direction WDIR, sensible heat flux H, surface friction velocity USTAR, mechanical boundary layer mixing height MixHGT M , Monin-Obukhov length L, Bowen ratio B0, albedo Albedo, wind speed WSPD, temperature Temp, rainfall PrecipAmt, relative humidity RH, station pressure Pressure and cloud cover CloudCover to form the basic variable set; under unstable conditions, add convective velocity scale WSTAR, potential temperature gradient ThetaGrad and convective mixing height MixHGT C to the basic variable set under stable conditions to reflect the atmospheric turbulence structure when convection dominates; In stable and unstable conditions, the Pearson correlation coefficient and the variance inflation factor VIF are used for statistical diagnosis; The calculation formula of the Pearson correlation coefficient is as follows: wherein, denotes the Pearson correlation coefficient between the weather variables x and y ; n is the total length of the time series of weather data samples; and are the mean values of the weather variables t = i and are the two weather variable observations at time instant and x are the mean values of the weather variables y over the time series. With As a reference basis for determining the existence of strong linear correlation between variables, it is used to assist in identifying specific variable pairs that lead to high VIF values; The calculation formula and explanation of VIF are as follows: wherein, is the coefficient of determination obtained from a linear regression analysis taking the first k meteorological variable as the dependent variable and all remaining meteorological variables as the independent variables; The screening process adopts the VIF iterative elimination strategy with physical meaning constraints: in each round, first calculate the VIF values of all variables in the current feature set, identify the candidate variables with the largest value and exceeding the warning threshold VIF>10, and determine the physical properties; if the variable and other retained variables have physical meaning overlap and functional redundancy, it will be removed; otherwise, if the variable represents an independent physical dimension, even if the VIF value is high, in order to ensure the integrity of the model in describing the physical processes of advection and diffusion, it is still determined that it cannot be removed and is forced to be retained, and the next highest value variable is investigated; when variable elimination occurs, the remaining variables are used to update the VIF values in the new round, and the above process is repeated until the VIF values of all non-forced retention variables meet the threshold requirements; Based on this strategy, the preliminary screening of the weather feature set obtained after the preliminary screening As follows: The features after preliminary screening under stable conditions are heat flux H, Monin-Obukhov length L, Bowen ratio B0, albedo Albedo, wind speed WSPD, rainfall PrecipAmt, relative humidity RH, and cloud cover CloudCover; The characteristics after the initial screening under unstable conditions are heat flux H, theta gradient ThetaGrad, convective mixing height MixHGT C , Monin-Obukhov length L, Bowen ratio B0, wind speed WSPD, precipitation PrecipAmt, relative humidity RH, and cloud cover CloudCover.

4. The method of claim 2, wherein the physical structured agent-based model is a physical structured agent-based model of road traffic pollution dispersion. Step 1.2 is specifically: A standardized line source scenario under stable and unstable weather conditions is constructed, a single receptor point is fixed, and the RLINE benchmark model is used to simulate the concentration at the receptor point; Constructing an XGBoost model to simulate concentrations at receptor points instead of RLINE model: with the initial set of meteorological features As input, the RLINE-simulated concentrations at receptor points as labels, train the XGBoost model to get the pre-trained XGBoost model under the standard line source scenario; The SHAP framework is used to quantify the marginal contribution of each variable and identify the key physical factors driving diffusion; Features i The formula for calculating the SHAP value of Xj is as follows: wherein, f is the XGBoost model under the pre-trained normalized line source scenario; N is the input feature set; S is the feature subset; represents the set difference operation, which means the feature set formed after removing the feature N from the universal set i ; represents the subset relationship, which means the summation symbol traverses all possible feature subsets ; S ; represents the cardinality of the set, that is, the number of features contained in the set within the double vertical lines; represents the set union operation, represents the new feature combination formed by adding the target feature i to the subset S ; is the prediction value of the model under the feature subset S ; The absolute value of the SHAP value of each feature is used to determine the global feature importance ranking, and the key features are selected according to their size to construct a feature vector; the selected key features are as follows: In stable conditions, the key feature vector is = [ WSPD, WDIR, L, H ] corresponding to wind speed, wind direction, Monin-Obukhov length, sensible heat flux, respectively. In unstable conditions, the key feature vector is = [ WSPD, WDIR, L, H, MixHGT C ] where MixHGT C is the convective mixing height.

5. The method of claim 2, wherein the physical structured agent-based model is a physical structured agent-based model of road traffic pollution dispersion. Step 1.3 is specifically: Under each weather condition, the K-means clustering method is used to group the hourly observation data by minimizing the sum of squared errors within clusters The optimal cluster center is determined, and the calculation formula is as follows: wherein K is the number of clustering clusters, is the kth weather cluster, is the centroid of the cluster; subsequently, from each cluster, the sample closest to the cluster center is selected as a representative weather instance, and a final representative weather set is constructed as the weather condition input for the subsequent agent model training.

6. The method of claim 1, wherein the method is characterized by, Step S2 is specifically: Step 2.1 Firstly, based on the node coordinate sequence of the road network topology data, the Euclidean distance between adjacent nodes is calculated and the road cumulative arc length sequence S is generated; then, according to the set standard discrete interval The number of microelement segments required for division is calculated, and the center position distance of each microelement segment along the road direction is determined d m ; finally, the linear interpolation algorithm is used to map the one-dimensional center position distance back to the two-dimensional plane to obtain the geometric center coordinates of each standardized microelement segment x m , y m ); the calculation formula is as follows: The recursive formula of the road cumulative arc length sequence S is The center position of the first micro-element section The center position of the second micro-element section d m is: determining the index j such that The center coordinates are then calculated as follows: wherein, is the total length of the road; n is the total number of road nodes; ( x j ， y j ) is the coordinate of the j th original road node; Step 2.2 For each normalized infinitesimal segment, first locate its belonging original road segment index by searching the cumulative arc length sequence k , and then calculate the radian value of the vector direction angle according to the coordinate difference of the road segment endpoints; further, convert the radian value in the mathematical polar coordinate system to the geographic azimuth angle with 0 degree as the south; finally, combine the axial symmetry characteristic of the line source diffusion to normalize the angle by modulo 180 degrees to obtain the standardized road alignment angle ; the calculation formula is as follows: Vector direction angle Calculation: Geographical azimuth conversion and normalization: wherein, , is the start coordinate of the original road segment with index k , is the vector direction angle in the mathematical polar coordinate system, is the geographic azimuth angle, denotes the floor function for eliminating the redundant half periods, ensuring ; Step 2.3 Based on the full-period traffic emission data, first traverse each time step with road ID as the association key t The discrete generated static geometric microelement set containing spatial coordinates and azimuth angle is matched with the dynamic emission intensity at this moment to construct a standardized space-time input tensor containing time stamp, microelement segment ID, center coordinates, strike angle, and instantaneous emission intensity.

7. The method of claim 1, wherein the method is characterized by, Step S3 is specifically: Step 3.1 constructs a non-uniformly distributed initial receptor point set under the benchmark wind direction; taking the midpoint of the unit line source as the origin, covering the continuous segment [-100m, 2000m] from the upwind background area to the downwind dilution area on the dominant wind direction axis, using a variable step size distribution strategy to construct an initial receptor point set with dense near-source and sparse far-end; Step 3.

2. According to the meteorological wind direction data at each moment, a two-dimensional rotation transformation matrix is constructed to map the initial sampling point under the reference wind direction to the actual geographical space to obtain the actual receptor point coordinates (x r , y r ); in this process, the geometric buffer zone threshold on both sides of the line source is set, and the singular points falling into the buffer zone are removed to avoid numerical calculation instability; the rotation transformation calculation formula is as follows: wherein, is the initial coordinate for the reference wind direction; is the angle difference between the wind direction at this moment and the reference wind direction; is the unit line source center coordinate; Step 3.3 Convert the actual receptor point coordinates into the wind-rotated coordinate system for feature representation; decouple the receptor point position orthogonally into the along-wind projection distance Distance perpendicular to the crosswind direction , respectively representing the structure of the pollutant along the streamline advection decay and the transverse turbulent diffusion; the calculation formula is as follows: wherein, is the measured flow field vector angle (in radian) corresponding to the wind direction; is the longitudinal distance in the wind direction coordinate system; is the lateral distance in the wind direction coordinate system.

8. The method of claim 1, wherein the method is characterized by, Step S4 is specifically: Step 4.1 Construct the meteorological stability partition based on the Monin-Obukhov length L; divide the meteorological samples into 5 meteorological stability categories: very stable, stable, neutral, unstable, very unstable, meteorological stability category The discrimination rule is as follows: Step 4.2: According to the dominant area difference of the "plume-wavy" double-component physical mechanism in the RLINE model, the along-wind projection distance of the receptor point in the wind direction coordinate system is taken as the criterion Constructing spatial geometric partition : when , it is determined as the wavy dominant area; when , it is determined as the plume dominant area; Step 4.3 constructs a sub-agent model corresponding to the RLINE physical solving module through the orthogonal combination of meteorological stability and spatial geometric position, and uses the XGBoost model to fit the physical response surface of each sub-model.

9. The method of claim 1, wherein the method is characterized by, Step S5 is specifically: Step 5.1 constructs a mixed feature vector containing physical meteorological properties and relative geometric properties, and selects the extremely gradient boosting tree XGBoost as the core regression engine; At the feature construction level, the wind direction is decomposed into orthogonal wind direction sine and cosine components, and combined with wind speed, Monin-Obukhov length, sensible heat flux and convective mixed layer height to form the whole-process meteorological features; at the same time, combined with the projection distance along the wind direction and the distance perpendicular to the side wind direction , a complete model input space is formed; Step 5.2 defines the objective function containing the training loss and structural risk regularization, assuming that the training set contains n samples, XGBoost approximates the true physical field by iteratively accumulating K decision trees; Objective function a differentiable convex loss function that measures prediction error and a regularization term that penalizes tree model complexity The composition is calculated as follows: wherein is the predicted concentration; is the true value for RLINE; is the mean squared error loss function; is the number of leaf nodes; is the L2 norm square of the leaf node scores; and γ and λ are respective penalty coefficients. Step 5.3 Determine the optimal hyperparameter combination, the key hyperparameters of the XGBoost model are set as follows: limit the maximum tree depth max_depth of the decision tree to 5; set the learning rate learning_rate to 0.15 and the number of iterations n_estimators to 1000; and set the subsampling ratio subsample to 0.8 to enhance robustness.

10. The method of claim 1, wherein the method is characterized by, Step S6 is specifically: Receiving real-time meteorological boundary conditions and dynamic road network emission inventory of the to-be-simulated period, for each receptor point in the target regional grid, traversing all unit line sources standardized in step S2, a "source-receptor" feature pair is constructed. According to real-time wind direction data, batch calculation of wind direction projection distance of receptor point relative to each line source Distance perpendicular to crosswind direction And combined with the atmospheric Monin-Obukhov length at the current time to determine the stability level, each feature vector is matched and input to the corresponding physical sub-agent model interface trained in step S5. Each physical sub-agent model performs parallel inference on the above feature batch to obtain the predicted value of the diffusion concentration of each unit line source at the receptor point under the normalized source strength; finally, based on the linear additivity principle of the Gaussian model, the predicted value is multiplied by the actual emission intensity of each line source at that time , and all line source contributions are reconstructed to obtain the total concentration field of regional pollutants with high spatiotemporal resolution; the calculation formula is as follows: wherein, is the total number of unit line sources affecting the receptor point in the area; is the actual emission source strength of the th line source at the moment in time; is the sub-agent model corresponding to the physical condition k ; is the th line source relative geometry and meteorology feature vector for the receptor point .

Images