A road traffic pollution diffusion fast simulation method based on a physical structured agent model

Abstract

The application belongs to the technical field of near-road side pollution diffusion simulation, and particularly relates to a road traffic pollution diffusion fast simulation method based on a physical structured proxy model. The method comprises the following steps: S1. meteorological feature engineering and sampling; S2. road network standardization decomposition, which discretizes and decomposes a complex urban road network topology into standardization unit line sources that can be independently calculated; S3. spatial geometric feature reconstruction; S4. physical mechanism double partition, which constructs an orthogonal physical mechanism partition system according to the atmospheric stability level and the relative position of the receptor point; S5. proxy model training, which trains an XGBoost regression model for each physical partition; S6. regional field fast reconstruction, which receives real-time boundary conditions for parallel inference and reconstructs the regional pollutant concentration field through linear superposition. The application significantly improves the calculation efficiency under the premise of ensuring the simulation accuracy consistent with the full physical model, and is suitable for real-time air quality simulation of large-scale urban road networks.

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:

[0007] A rapid simulation method for road traffic pollution diffusion based on a physically structured surrogate model includes the following steps:

[0008] 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;

[0009] 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;

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

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

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

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

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

[0015] Beneficial effects

[0016] Compared with the prior art, the present invention has the following significant advantages:

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

[0018] 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".

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

[0020] 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

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

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

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

[0024] 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).

[0025] 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;

[0026] 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;

[0027] Figure 7This 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

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

[0029] 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:

[0030] 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;

[0031] S2. Road Network Standardization Decomposition: Discretizes and decomposes the complex urban road network topology into independently computable standardized unit line sources;

[0032] 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;

[0033] 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;

[0034] S5. Proxy Model Training: Train XGBoost regression models for each physical partition separately to learn the nonlinear response surface of the physical model;

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

[0036] 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:

[0037] Step 1.1 To meet the meteorological input requirements of the RLINE model, construct a differentiated feature system based on atmospheric stability.

[0038] 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. MThe 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.

[0039] 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).

[0040] The formula for calculating the Pearson correlation coefficient and its explanation are as follows:

[0041]

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

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

[0044] The formula and explanation for calculating VIF are as follows:

[0045]

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

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

[0048] 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:

[0049] 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).

[0050] 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).

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

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

[0053] The SHAP framework is used to quantify the marginal contributions of each variable and identify the key physical factors driving diffusion. Features iThe formula for calculating the SHAP value is as follows:

[0054]

[0055] 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 S The predicted value is below.

[0056] 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:

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

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

[0059] 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:

[0060]

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

[0062] 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:

[0063] 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 to be divided 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:

[0064] The recursive formula for the cumulative arc length sequence S of the road is:

[0065]

[0066] No. The center position of each infinitesimal segment d m for:

[0067]

[0068] Determine the index j Make The center coordinates are calculated as follows:

[0069]

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

[0071] 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:

[0072] Vector direction angle calculate:

[0073]

[0074] Geographic azimuth conversion and normalization:

[0075]

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

[0077] Step 2.3 Based on all-time traffic emission data, 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.

[0078] 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:

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

[0080] 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:

[0081]

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

[0083] 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:

[0084]

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

[0086] 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:

[0087] 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:

[0088]

[0089] 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:

[0090]

[0091] 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).

[0092] 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:

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

[0094] 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:

[0095]

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

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

[0098] 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:

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

[0100] 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:

[0101]

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

[0103] Example

[0104] 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:

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

[0106] 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:

[0107] (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:

[0108]

[0109] 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 variablesx 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.

[0110] The formula and explanation for calculating VIF are as follows:

[0111]

[0112] 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. :

[0113] 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).

[0114] 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).

[0115] 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. M A 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.

[0116] Table 1 List of candidate meteorological features

[0117]

[0118] Table 2 VIF Analysis under Stable Meteorological Conditions

[0119]

[0120] Table 3. VIF Analysis under Unstable Meteorological Conditions

[0121]

[0122] (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:

[0123]

[0124] in, f For pre-trained models; N For the input feature set; SA 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 S The predicted value is below.

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

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

[0127] (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:

[0128]

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

[0130] 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:

[0131] (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:

[0132] The recursive formula for the cumulative arc length sequence S of the road is:

[0133]

[0134] The center position of the m-th infinitesimal segment d m for:

[0135]

[0136] Determine the index j Make The center coordinates are calculated as follows:

[0137]

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

[0139] (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. kThe 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:

[0140] Vector direction angle calculate:

[0141]

[0142] Geographic azimuth conversion and normalization:

[0143]

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

[0145] (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.

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

[0147] Table 4. Standardized decomposition of road network unit information (partial list)

[0148]

[0149] 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:

[0150] (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.

[0151] (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:

[0152]

[0153] 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).

[0154] (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:

[0155]

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

[0157] 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:

[0158] (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:

[0159]

[0160] (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:

[0161]

[0162] (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).

[0163] 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:

[0164] (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.

[0165] (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:

[0166]

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

[0168] (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.

[0169] 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 5 As 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).

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

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

[0172] 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:

[0173]

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

[0175] 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 7As 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.

[0176] Table 5 Overall Predictive Performance Statistics of the Model

[0177]

[0178] Table 6. Details of concentration errors between the receptor-by-receptor surrogate model and RLINE in typical plains scenarios (partial list only)

[0179]

[0180] Table 7. Concentration error details between the receptor-per-point surrogate model and RLINE in typical suburban town scenarios (partial list only).

[0181]

[0182] 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).

[0183]

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

[0185] Table 9 Comparison of computational efficiency between the proxy model and the baseline model RLINE in the same scenario

[0186]

[0187] 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 as equivalent substitutions and included within the protection scope of the present invention.

Claims

1. A rapid simulation method for road traffic pollution diffusion based on a physically structured surrogate model, characterized in that, Includes the following steps: S1. Meteorological Feature Engineering and Sampling: Multicollinearity diagnosis of meteorological parameters is performed based on Pearson correlation coefficient and variance inflation factor (VIF), and key features driving diffusion are screened out by combining SHAP attribution analysis; a clustering algorithm is used to perform stratified sampling of historical meteorological data 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 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, a wind-rotating coordinate system is constructed to convert the geographic absolute coordinates of the receiver point into the wind direction projection distance and crosswind vertical distance relative to the line source. A non-uniform spatial sampling strategy is used to generate a receiver point sample set containing the near-source high gradient region and the far-end background region. The rotated actual coordinates are further converted to the wind direction coordinate system for feature engineering expression. S4. Physical mechanism dual zoning: Based on the Moning-Obukhov length, the atmospheric turbulence state is divided into different stability levels, and based on the relative position of the receptor point, it is divided into plume-dominated region and meander-dominated region, 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.

2. The method for rapid simulation of road traffic pollution diffusion based on a physically structured proxy model as described in claim 1, characterized in that, In step S1, the specific method of meteorological feature engineering includes: The specific steps are as follows: Step 1.1 For the meteorological input parameters required by the RLINE model, sampling was performed under both stable and unstable conditions. Statistical diagnosis was conducted using both Pearson correlation coefficient and variance inflation factor (VIF) to initially eliminate redundant variables; the resulting preliminary meteorological feature set was then obtained. ; Step 1.2 Construct standardized line source scenarios under stable and unstable conditions, use the SHAP framework to quantify the marginal contribution of each variable, and identify the key physical factors driving diffusion; Step 1.3 Perform stratified representative sampling on the original meteorological dataset; under each meteorological condition, use the K-means clustering method to group the hourly observation data, determine the optimal cluster center by minimizing the sum of squared errors within the cluster, and select the sample closest to the cluster center as a representative meteorological instance to ensure the efficiency and generalization ability of the surrogate model during training.

3. The method for rapid simulation of road traffic pollution diffusion based on a physically structured surrogate model as described in claim 2, characterized in that, In step 1.1, Select the basic variable set from the meteorological input parameters required for the RLINE model: under steady-state conditions, select wind direction WDIR, sensible heat flux H, surface friction velocity USTAR, and mechanical boundary layer mixing height MixHGT. 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. Under unstable conditions, the convective velocity scale WSTAR, potential temperature gradient ThetaGrad, and convective mixing height MixHGT are added to the basic variable set for stable conditions. C Three parameters to reflect the atmospheric turbulence structure when convection is dominant; Sampling was performed under both stable and unstable conditions, and statistical diagnosis was conducted using a dual index of Pearson correlation coefficient and variance inflation factor (VIF). The formula for calculating the Pearson correlation coefficient is 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; by As a reference for determining whether there is a strong linear correlation between variables, it is used to help identify specific variable pairs that lead to high VIF values; 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. The screening process employs a VIF iterative elimination strategy with physical meaning constraints: 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 overlaps in physical meaning or has functional redundancy with other retained variables, it is eliminated. Conversely, if the variable represents an independent physical dimension, 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 the 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. The initial meteorological feature set obtained after preliminary screening based on this strategy 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. 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 are: Morning-Obukhov length (L), Bowen-B0, wind speed (WSPD), precipitation (PrecipAmt), relative humidity (RH), and cloud cover (CloudCover).

4. The method for rapid simulation of road traffic pollution diffusion based on a physically structured surrogate model as described in claim 2, characterized in that, Step 1.2 specifically involves: Standardized line source scenarios under stable and unstable weather conditions are constructed, a single receptor point is fixed, and the concentration at the receptor point is simulated using the benchmark model RLINE; Constructing an XGBoost model to replace the RLINE model for simulating receptor concentrations: using a preliminary screening of meteorological feature sets. Using the receptor concentration simulated by the baseline model RLINE as the label, the XGBoost model is trained to obtain a pre-trained XGBoost model in the normalized line source scenario. The SHAP framework is used to quantify the marginal contribution of each variable and identify the key physical factors driving diffusion. feature 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; 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 S The predicted value below; The global feature importance ranking is determined based on the absolute value of the SHAP value of each feature, and key features are selected and feature vectors are constructed based on their magnitude; the final selected key features are as follows: Under stable conditions, the key feature vector 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.

5. The method for rapid simulation of road traffic pollution diffusion based on a physically structured proxy model as described in claim 2, characterized in that, Step 1.3 specifically involves: 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, The centroid of the cluster is then selected. Subsequently, the sample closest to the cluster center is selected 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.

6. The method for rapid simulation of road traffic pollution diffusion based on a physically structured proxy model as described in claim 1, characterized in that, Step S2 is 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... Calculate the number of micro-segments required for the road 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; 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 polar coordinates is converted to a geographic azimuth angle with true south as 0 degrees. 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 The vector direction angle in the mathematical polar coordinate system. It is the geographical azimuth. This represents the floor function, used to remove redundant half-cycles to ensure... ; Step 2.3 Based on all-time traffic emission data, 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.

7. The method for rapid simulation of road traffic pollution diffusion based on a physically structured proxy model as described in claim 1, characterized in that, Step S3 is as follows: Step 3.1 Construct an initial set of receptor points with a non-uniform distribution under the reference wind direction; with the midpoint of the unit line source as the origin, cover a continuous section [-100m, 2000m] from the upwind background area to the downwind dilution area on the prevailing wind axis, and construct the initial set of receptor points by adopting a variable step size distribution strategy of dense near source and sparse far end. 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, in order to avoid instability in numerical calculations; 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; Step 3.3: Transform the actual receptor point coordinates to a wind-rotating coordinate system for feature representation; orthogonally decouple the receptor point position into a wind-direction projected distance. Vertical distance from crosswind direction These represent the advection transport attenuation of pollutants along the streamline and the turbulent diffusion structure in the lateral direction, respectively; the calculation formulas are as follows: in, This represents the flow field vector angle 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.

8. The method for rapid simulation of road traffic pollution diffusion based on a physically structured proxy model as described in claim 1, characterized in that, Step S4 is as follows: Step 4.1 Construct meteorological stability partitions based on the Moning-Obukhov length L; divide 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-meander" dual-component physical mechanism in the RLINE model, and based on the windward projection distance of the receptor point in the wind direction coordinate system... Constructing spatial geometric partitions :when When, it is determined to be a meandering dominant area; when At that time, it was determined to be the area dominated by the plume; Step 4.3 By orthogonally combining meteorological stability and spatial geometric location, construct sub-surrogate models that correspond one-to-one with the RLINE physics solution module, and fit the physical response surface of each sub-model using the XGBoost model.

9. The method for rapid simulation of road traffic pollution diffusion based on a physically structured proxy model as described in claim 1, characterized in that, Step S5 is as follows: Step 5.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; At the feature construction level, 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 to construct the meteorological characteristics of the entire process; simultaneously, the projected distance along the wind direction is also considered. Vertical distance from crosswind direction This forms a complete model input space; Step 5.2 Define the objective function that includes training loss and structural risk regularization. Assuming that the training set contains n samples, XGBoost approximates the real physical field by iteratively accumulating K decision trees. Its objective function The differencer convex loss function measures the prediction error and the regularization term penalizes the complexity of the tree model. 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; Let L2 be the squared value of the leaf node score; γ and λ are the corresponding penalty coefficients, respectively. Step 5.3 determines the optimal hyperparameter combination and makes targeted settings for the key hyperparameters of the XGBoost model: the maximum tree depth max_depth of the decision tree is limited to 5; the learning rate is set to 0.15 and the number of iterations n_estimators is set to 1000; and the subsample ratio is set to 0.8 to enhance robustness.

10. The method for rapid simulation of road traffic pollution diffusion based on a physically structured surrogate model as described in claim 1, characterized in that, Step S6 is as follows: Receive the real-time meteorological boundary conditions and dynamic road network emission inventory for the period to be simulated. For each receptor point in the target area grid, traverse all unit line sources after standardization in step S2 and construct "source-receptor" feature pairs. Based on real-time wind direction data, the wind direction projection distance of the receiver 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. Each physical sub-proxy model performs parallel inference on the aforementioned feature batches to obtain the predicted diffusion concentration of each unit line source at the acceptor point under 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 moment. The total 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.

Images