Urban forest SO2 purification benefit remote sensing inversion method based on coupling model

Abstract

The invention relates to the technical field of ecological environment remote sensing monitoring and evaluation, and discloses an urban forest SO2 purification benefit remote sensing inversion method based on a coupling model, and the method comprises the steps: firstly generating a normalized digital surface model and a leaf area index based on multi-source data; secondly, urban canopy three-dimensional aerodynamic resistance modeling is executed, and the friction speed and the microenvironment ventilation coefficient are solved; further constructing a particulate matter stomatal co-stress model, calculating an actual blade dust accumulation load and calculating corrected canopy stomatal resistance; meanwhile, a sulfur dioxide effective exposure concentration field is reconstructed based on the source-sink neighbor degree; and finally, synthesizing the dry settling velocity by integrating the aerodynamic resistance, the quasi-laminar flow boundary layer resistance and the corrected canopy pore resistance, and inverting the purification flux and the total amount. By coupling microenvironment ventilation, particulate matter physical blockage and non-uniform concentration field characteristics, the physical authenticity and accuracy of forest ecological benefit evaluation in a complex urban environment are improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of ecological environment remote sensing monitoring and evaluation, in particular to a remote sensing inversion method for urban forest SO2 purification benefit based on a coupling model. BACKGROUND

[0002] As an important part of urban ecosystems, urban forests play a key ecological service function in improving air quality by absorbing gaseous pollutants such as sulfur dioxide (SO2) through the process of dry deposition. Currently, the mainstream method for quantifying the SO2 purification benefit of urban forests is mostly based on the Big-Leaf Model or the Multi-layer Canopy Model, which estimates the purification flux by calculating the product of the dry deposition velocity and the concentration of atmospheric pollutants. The dry deposition velocity is usually calculated using the resistance model of an analog circuit, which involves the combined effects of aerodynamic resistance, quasi-laminar boundary layer resistance, and canopy stomatal resistance.

[0003] Although existing dry deposition models have been relatively mature in homogeneous underlying surfaces such as farmland and contiguous forests, they still have significant technical limitations when applied to highly heterogeneous urban environments. In terms of aerodynamic parameter acquisition, existing techniques rely on looking up empirical coefficient tables for land use types to determine the zero plane displacement and roughness length. However, the complex three-dimensional structure of buildings and vegetation in urban underlying surfaces produces strong drag and disturbance effects on airflow, and simple lookup tables cannot represent the turbulence characteristics in the microenvironment, resulting in significant deviations in the calculation of friction velocity and aerodynamic resistance, which further affects the estimation accuracy of quasi-laminar boundary layer resistance.

[0004] In terms of plant physiological response mechanisms, traditional models mainly focus on the physiological limiting effects of meteorological factors such as light, temperature, and vapor pressure deficit on stomatal conductance, but generally ignore the physical attachment effects of high concentrations of particulate matter (such as PM2.5, PM10, and dust) on leaves in urban environments. In fact, urban forest leaves are often under dust stress, and the non-uniform accumulation of particulate matter on leaf surfaces can physically block stomatal channels, significantly increasing gas exchange resistance. Existing models assume clean leaves or only make simple general corrections, without considering the dynamic regulation mechanisms of microenvironment ventilation conditions on dust accumulation and resuspension, thus often underestimating stomatal resistance and overestimating SO2 absorption capacity in typical urban habitats with combined pollution of particulate matter and SO2.

[0005] Furthermore, in constructing pollutant concentration fields, due to the sparse distribution of ground-based air quality monitoring stations, existing technologies often employ Kriging interpolation or inverse distance weighting to generate regional background concentration fields. This smoothing approach flattens out the localized high-concentration characteristics around major traffic arteries and industrial emission points. Since the amount of pollutants absorbed by plants is directly related to environmental exposure concentration, ignoring high concentration gradients near emission sources leads to a severe underestimation of the vegetation's contribution to purification in these key areas (hotspots). In summary, existing technologies lack a comprehensive inversion model that can organically couple the complex urban canopy aerodynamic environment, particulate matter physical blockage stress mechanisms, and non-uniform pollutant exposure characteristics, making it difficult to meet the precision requirements of refined urban ecosystem management. Summary of the Invention

[0006] To address the shortcomings of existing technologies, this invention provides a remote sensing inversion method for assessing the SO2 purification benefits of urban forests based on a coupled model. Existing technologies for assessing the SO2 purification benefits of urban forests suffer from inaccurate assessment results because they neglect the complex aerodynamic characteristics of urban canopies, the physical blocking effect of particulate matter accumulation on pores, and the uneven spatial distribution of pollutants.

[0007] To achieve the above objectives, the present invention provides a remote sensing inversion method for SO2 purification benefits in urban forests based on a coupled model, which aims to solve the problem that the existing technology fails to fully consider the complex aerodynamic environment of urban underlying surfaces, the physical blocking effect of particulate matter on pores, and the spatial heterogeneity of pollution source distribution, thus leading to biases in the evaluation of purification benefits.

[0008] This method includes multi-source data acquisition and basic parameter preprocessing steps. By acquiring multispectral remote sensing images, digital surface models, digital elevation models, and regional meteorological and atmospheric environmental monitoring data of the target area, normalized digital surface model data representing the absolute height of surface objects is constructed, and vegetation leaf area index data is generated through inversion. In this process, surface reflectance is obtained through radiometric calibration and atmospheric correction, and a nonlinear mapping relationship between the normalized vegetation index and the leaf area index is established; at the same time, difference calculation and extreme value filtering are used to eliminate the influence of topographic relief and extract pure surface cover height information.

[0009] After acquiring basic data, this method performs three-dimensional aerodynamic drag modeling of the urban canopy and calculates the microenvironment ventilation factor. This step utilizes a moving window algorithm to traverse a normalized digital surface model, extracting morphological parameters such as average obstacle height, windward area density, and sky openness. Based on the morphological density model, the zero-plane displacement and aerodynamic roughness length are inverted, and then the friction velocity and aerodynamic drag are iteratively calculated using the Moning-Obukhov similarity theory. Furthermore, this method calculates the microenvironment ventilation coefficient, which is positively correlated with the reference wind speed and sky openness, and positively correlated with the complement of the windward area density, used to quantify the local canopy's airflow exchange capacity and pollutant diffusion potential.

[0010] To address the impact of particulate matter adhesion on plant leaves in urban environments, this method constructs a stomatal stress model based on microenvironment ventilation efficiency. First, the microenvironment ventilation coefficient is used as a nonlinear gain regulator to map aerosol optical thickness data. In areas with low ventilation coefficients, the exponential term is increased to reflect localized pollutant enrichment, calculating the potential dust load. Then, wind-induced resuspension correction based on friction velocity is performed: by comparing the friction velocity with the critical friction velocity threshold for particulate resuspension pixel by pixel, static retention or dynamic cleaning mode is determined. If dynamic cleaning mode is in effect, the amount of particulate matter stripping due to excess shear force is calculated to reduce the potential dust load, yielding the actual leaf dust load. Furthermore, a leaf surface roughness sensitivity coefficient is introduced, and a stomatal blockage stress factor is constructed using a Langmuir adsorption model to quantify the degree of physical obstruction of gas exchange channels by particulate matter coverage. Finally, by combining the response functions of light, temperature and water, the stomatal blockage stress factor was introduced as a physical correction term into the multiplicative resistance model to correct the minimum stomatal resistance of vegetation, and the corrected canopy stomatal resistance, which reflects the dual physical and physiological mechanisms, was obtained.

[0011] To address the uneven distribution of pollutants within cities, this method reconstructs the effective exposure concentration field of sulfur dioxide based on source-sink proximity. This step identifies sulfur dioxide emission sources within the study area and assigns them intensity coefficients. Based on the inverse distance law, the source-nearest neighbor weights of the pixels to be evaluated relative to each emission source are calculated. These weights are then used to perform pixel-level reallocation corrections on the smoothed baseline concentration field generated by interpolation from ground monitoring stations: concentration values ​​are increased in areas where the source-nearest neighbor weights are higher than the overall average, and vice versa, thereby generating effective exposure concentration data of sulfur dioxide that reflects local pollution hotspots.

[0012] Finally, the method performs purification efficiency flux calculation and total amount inversion. Quasi-laminar boundary layer resistance is calculated based on friction velocity and gas molecule characteristic parameters. A series resistance model is used to sum and reciprocate the aerodynamic resistance, quasi-laminar boundary layer resistance, and corrected canopy porosity resistance to synthesize the dry deposition velocity. Instantaneous purification flux is obtained by performing pixel-by-pixel multiplication of the dry deposition velocity and effective sulfur dioxide exposure concentration data. Integration over the spatiotemporal dimensions yields the total sulfur dioxide purification amount for urban forests.

[0013] This invention provides a remote sensing inversion method for SO2 purification benefits in urban forests based on a coupled model. It has the following beneficial effects:

[0014] 1. This invention achieves dynamic extrapolation of the actual dust load on blades by constructing a nonlinear mapping mechanism based on the microenvironment ventilation coefficient and a wind-induced resuspension correction model based on friction velocity. Furthermore, it introduces a stomatal blockage stress factor to physically correct the stomatal resistance of the canopy. This mechanism overcomes the defect of traditional models that ignore the shielding effect of particulate matter physical retention on the gas exchange channels of stomata, and significantly improves the accuracy of dry deposition velocity calculation in complex urban atmospheric environments.

[0015] 2. This invention employs a concentration field reconstruction technology based on source-sink nearest neighbor. It utilizes the spatial distribution characteristics of emission sources to construct a source-nearest neighbor weight field and performs pixel-level non-uniform redistribution of background concentration. This feature can effectively restore the local high-concentration exposure characteristics around traffic arteries and industrial areas, ensuring that when calculating purification flux, the vegetation in high-exposure areas can match the corresponding high-concentration boundary conditions, thus eliminating the peak-shaving and valley-filling errors caused by traditional spatial interpolation methods.

[0016] 3. This invention utilizes a normalized digital surface model to extract three-dimensional morphological parameters such as average obstacle height and windward area density, and iteratively calculates friction velocity and aerodynamic drag based on the Moning-Obukhov similarity theory. Compared to traditional methods that rely solely on empirical constants assigned by land use type, this invention accurately reflects the drag and disturbance effects of mixed urban building and vegetation surfaces on airflow, providing more precise dynamic input parameters for calculating quasi-laminar boundary layer resistance and final purification flux. Attached Figure Description

[0017] Figure 1 This is a schematic diagram of the method flow of the present invention;

[0018] Figure 2 This is a schematic diagram of the wind-induced resuspension correction logic based on friction speed according to the present invention. Detailed Implementation

[0019] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0020] See attached document Figure 1 This invention provides a remote sensing inversion method for SO2 purification benefits in urban forests based on a coupled model. The method includes the following steps:

[0021] Step S100 involves multi-source data acquisition and basic parameter preprocessing. This step specifically includes: acquiring multispectral remote sensing image data, digital surface model data, digital elevation model data, regional meteorological data, and atmospheric environmental monitoring data for the target area. Radiometric calibration and atmospheric correction are performed on the multispectral remote sensing image data, the normalized vegetation index (NVI) is calculated, and the leaf area index is inverted based on the NVI. Difference calculations are performed between the digital surface model data and the digital elevation model data to generate normalized digital surface model data, which is used to characterize the absolute height information of surface objects. The regional meteorological data includes wind speed, wind direction, temperature, humidity, and photosynthetically active radiation; the atmospheric environmental monitoring data includes aerosol optical thickness and sulfur dioxide background concentration.

[0022] Step S200 involves performing three-dimensional aerodynamic drag modeling of the urban canopy and calculating the microenvironment ventilation factor. This step specifically includes: calculating pixel-level morphological parameters based on the normalized digital surface model data using a moving window algorithm. These morphological parameters include average obstacle height, windward area density, and sky openness. Based on these morphological parameters, aerodynamic roughness parameters are inverted using a morphological density model. These aerodynamic roughness parameters include zero-plane displacement and roughness length. Based on the Moning-Obukhov similarity theory, friction velocity and aerodynamic drag are calculated using the aerodynamic roughness parameters and the regional meteorological data. Based on the sky openness, the windward area density, and the reference wind speed in the regional meteorological data, the microenvironment ventilation coefficient is calculated. This microenvironment ventilation coefficient characterizes the airflow exchange capacity of the local canopy.

[0023] Step S300: Construct a particulate matter stomatal synergistic stress model based on microenvironment ventilation efficiency. This step specifically includes: using the microenvironment ventilation coefficient to nonlinearly map the aerosol optical thickness, calculating the potential dust load on the blade surface; using the friction velocity as a dynamic criterion to perform wind-induced resuspension correction on the potential dust load, calculating the actual blade dust load; constructing a stomatal blockage stress factor based on the actual blade dust load; introducing the stomatal blockage stress factor into the stomatal resistance calculation model, and combining it with the physiological response functions of light, temperature, and saturated vapor pressure difference to calculate the corrected canopy stomatal resistance.

[0024] Step S400: Reconstruct the effective SO2 exposure concentration field based on source-sink proximity. This step specifically includes: identifying sulfur dioxide emission sources based on land use classification data, including line sources and point sources; calculating the Euclidean distance from the target pixel to the emission source and constructing an emission source weight field based on the Euclidean distance; and using the emission source weight field to weight and correct the background sulfur dioxide concentration at ground monitoring stations, generating spatially continuous effective sulfur dioxide exposure concentration data.

[0025] Step S500 involves calculating the purification efficiency flux and performing total inversion. This step specifically includes: calculating the quasi-laminar boundary layer resistance based on the friction velocity; synthesizing the aerodynamic resistance, the quasi-laminar boundary layer resistance, and the corrected canopy stomatal resistance to calculate the dry deposition velocity; calculating the purification flux per pixel based on the dry deposition velocity and the effective sulfur dioxide exposure concentration data; and performing regional and temporal integration on the purification flux to output the total sulfur dioxide purification amount for the urban forest.

[0026] Step S100 aims to construct a unified spatiotemporal reference set of urban canopy environmental parameters, providing standardized input data for subsequent aerodynamic and physiological coupling calculations. This step specifically includes the following technical execution processes:

[0027] Phase 1: Acquisition and Screening of Multi-Source Heterogeneous Data. The system first establishes a multi-source data index based on the geographical extent of the area to be evaluated and the evaluation period. For remote sensing imagery data, optical satellite imagery with high revisit cycles (e.g., Sentinel-2 or Landsat-8 series) is preferred, while historical archived data with cloud cover below a preset threshold (e.g., 10%) is screened. The acquired imagery bands must include at least the visible red band and the near-infrared band (NIR) for vegetation spectral feature extraction. For three-dimensional geometric data, airborne LiDAR point cloud data of the target area or high-resolution digital surface models (DSM) and digital elevation models (DEM) generated based on stereo satellite imagery are acquired. The digital surface model represents the actual surface elevation including building tops, tree canopies, and bare ground, while the digital elevation model represents the topographic elevation after removing surface cover. For meteorological and environmental data, hourly observation records from meteorological stations within the area are accessed to obtain wind speed, wind direction (expressed as angle), dry-bulb temperature, relative humidity, and photosynthetically active radiation intensity during the evaluation period. Simultaneously, atmospheric environmental monitoring data for the corresponding time period is acquired, including the hourly average concentration of sulfur dioxide (SO2) at ground stations and aerosol optical thickness (AOD) products retrieved from satellites.

[0028] Phase Two: Unifying the Spatiotemporal Reference and Calculating Geometric Parameters. Since multi-source data have different spatial resolutions and coordinate systems, the system performs unified resampling and registration operations. A unified spatial projection coordinate system (e.g., WGS84 / UTM projection) and a standard grid resolution (e.g., 10m × 10m) are set. For point data from meteorological stations and ground environmental monitoring stations, spatial interpolation algorithms are used to map them to the standard grid, generating a meteorological parameter raster layer. For aerosol optical thickness products with lower resolution, bilinear interpolation is used to resample to the standard grid resolution. Based on this, a normalized digital surface model (nDSM) is constructed. The system reads the digital surface model values ​​and digital elevation model values ​​under the same geographic coordinates and performs a difference operation. Specifically, the pixel values ​​of the digital surface model are subtracted from the corresponding pixel values ​​of the digital elevation model to obtain the absolute height of surface objects. The system further performs extreme value filtering on the nDSM data to remove abnormally high and negative values ​​caused by sensor noise, outputting normalized height field data reflecting the true height distribution of urban buildings and forest vegetation.

[0029] The third stage: Radiometric calibration and vegetation physiological parameter inversion. For the original multispectral remote sensing image, the system performs radiometric calibration and atmospheric correction. Using sensor calibration parameters, the digital quantization (DN) values ​​are converted to top-atmospheric radiance, and atmospheric scattering and absorption effects are eliminated based on a radiative transfer model (such as the 6S model or the SEN2COR algorithm) to obtain the true surface reflectance data. Based on the corrected surface reflectance, the Normalized Difference Vegetation Index (NDVI) is calculated. The calculation logic is: take the difference between the near-infrared band reflectance and the red band reflectance, and divide by the sum of the two. The system sets a vegetation mask threshold to remove interfering pixels from water bodies and non-vegetated areas. Subsequently, the Leaf Area Index (LAI) is inverted based on the NDVI. This embodiment uses an empirical regression model or a physical model lookup table method to establish the mapping relationship between NDVI and LAI. The specific mapping logic follows the vegetation growth law, that is, within a certain range, the Leaf Area Index increases exponentially or non-linearly with the increase of the NDVI. The generated LAI raster data serves as a key physiological parameter for subsequent calculations of canopy stomatal resistance and pollutant adsorption surface area.

[0030] Through the above processing, the system completed the transformation from raw multi-source data to standardized geometric height field (nDSM), physiological parameter field (LAI), and environmental meteorological field (meteorological grid), realizing accurate registration of data in spatial dimensions and quantitative inversion of physical attributes.

[0031] Step S200 aims to address the distortion in wind speed estimation caused by using a fixed roughness constant in traditional models on non-uniform urban surfaces, and to provide key dynamic parameters for subsequent particulate matter co-correction. This step specifically includes the following technical execution processes:

[0032] Phase 1: Pixel-level extraction of urban canopy morphological parameters

[0033] Based on the normalized digital surface model (nDSM) generated in step S100, the system employs a moving window traversal algorithm to extract refined morphological features. First, a physically meaningful computational window (e.g., 3×3 or 5×5 pixels) is defined, its size adapted to the average block scale of the evaluation area. Within each computational window, the system first calculates the average obstacle height (…). ), which is the arithmetic mean of all non-zero elevation pixel values ​​within the window.

[0034] Subsequently, the system calculates the windward area density. This parameter is not a fixed value, but is dynamically related to the prevailing wind direction. The system reads the prevailing wind direction from the regional meteorological data, calculates the total projected area of ​​all land features (buildings and vegetation) within the window on the plane perpendicular to the wind direction, and divides this total projected area by the bottom area of ​​the calculation window to obtain the dimensionless windward area density. This parameter characterizes the degree of physical obstruction to the airflow in the horizontal direction.

[0035] Simultaneously, the system calculates the sky openness (SVF). Centered on the target pixel, the system emits virtual rays to various azimuth and elevation angles in the hemispherical space, and counts the proportion of rays not obstructed by surrounding obstacles (such as tall buildings or dense tree canopies). The sky openness value ranges from 0 to 1, with values ​​closer to 0 indicating the bottom of a deep street canyon, and values ​​closer to 1 indicating an open area or the top of a canopy.

[0036] Phase 2: Dynamic Inversion of Aerodynamic Roughness Parameters

[0037] Traditional methods often assign a single roughness constant to urban areas, while this embodiment uses a morphological density model for pixel-level dynamic inversion. The system establishes a mapping relationship from geometric shape to aerodynamic parameters. First, the zero-plane displacement is calculated ( The physical meaning of this parameter is the average height of the momentum sink, that is, the height at which the airflow feels the ground being lifted. The calculation logic follows a non-linear growth law: as the windward area density increases, the zero-plane displacement shows a trend of first rising rapidly and then leveling off, indicating that the airflow is gradually lifted to the canopy or the top of the building.

[0038] Secondly, calculate the aerodynamic roughness length. This parameter characterizes the frictional drag capability of the underlying surface against airflow. The calculation logic considers the standard deviation and distribution density of obstacle height. In low-density areas, The density increases; in high-density areas (such as closely packed buildings), due to the skimming flow effect, the fluid slides over the top and no longer penetrates into the gaps, resulting in... Instead, it decreases. The system precisely characterizes this hydrodynamic feature through parametric formulas.

[0039] Phase 3: Spatial calculation of friction velocity and aerodynamic drag

[0040] Based on the Monin-Obukhov similarity theory, the system combines regional meteorological observations of wind speed with the aforementioned inverted aerodynamic parameters to calculate the frictional velocity of each pixel. Friction velocity is a key velocity scale characterizing turbulent shear stress in the near-surface layer. In the calculation, the system uses a logarithmic wind profile equation to extrapolate the macroscopic wind speed from the reference height to the canopy height. During this process, the unique zero-plane displacement and roughness length of each pixel correct for the vertical attenuation rate of the wind speed.

[0041] Based on the friction velocity, the system further calculates the aerodynamic drag ( This drag characterizes the turbulent transport resistance of SO2 molecules from a reference height to the top of the quasi-laminar boundary layer. Due to the introduction of pixel-level aerodynamic parameters, this drag field exhibits significant spatial heterogeneity: the drag is low in open parklands and significantly increases in densely built-up leeward areas.

[0042] Phase 4: Construction of the Microenvironment Ventilation Coefficient (VC)

[0043] This is the key step in achieving coordinated correction of aerodynamics and dust accumulation in this invention. The system does not directly use macroscopic wind speed, but instead constructs a comprehensive index that reflects the cleaning capability of the microenvironment—the microenvironment ventilation coefficient (VC).

[0044] The calculation logic of this coefficient comprehensively considers both macroscopic driving forces and local constraints. The system uses wind speed at a reference height as the basic driving term, sky openness (SVF) as the weighting term for vertical exchange capacity, and windward area density (…). The complement of the ventilation coefficient (i.e., 1 minus the obstruction rate) is used as the weighting term for horizontal airflow capacity. Specifically, the value of the ventilation coefficient is directly proportional to the reference wind speed, positively correlated with the openness of the sky, and negatively correlated with the windward area density.

[0045] The calculated ventilation coefficient field (VCField) is a spatially continuous raster map. Regions with high VC values ​​represent areas with good airflow permeability and strong turbulent exchange, which are conducive to the diffusion of pollutants and the resuspension of particulate matter. Regions with low VC values ​​represent dead zones in the flow field (such as the bottom of a street canyon), where pollutants are prone to accumulation and retention. This parameter will be directly used as an input variable and passed to the subsequent particulate matter aeroporotic stress model.

[0046] As the first stage in constructing a particulate matter stomatal synergistic stress model, step S301 aims to address the technical challenge that traditional methods relying solely on atmospheric optical parameters cannot characterize the actual physical cover state of near-ground blades. This step specifically includes the following technical execution processes:

[0047] Phase 1: The system constructs a logic for converting atmospheric column concentration into near-surface deposition potential energy. It reads the aerosol optical thickness (AOD) raster data obtained in step S100 and the microenvironment ventilation coefficient (VC) raster data calculated in step S200. This embodiment constructs a mapping model based on the following physical mechanism: Aerosol optical thickness characterizes the total extinction effect within the atmospheric vertical column. Within the urban canopy, the vertical distribution of particulate matter is strongly modulated by the microenvironmental airflow structure. In open areas with high VC, airflow exchange is smooth, pollutants are easily diffused and diluted, and the near-surface blade deposition probability is relatively close to the background level. However, in areas with low VC (such as deep street canyons or the leeward side of dense building complexes), airflow stagnation and vortex backflow occur, leading to particulate matter enrichment in the near-surface layer, thus significantly increasing the contact deposition probability of the blades. Therefore, the system not only uses AOD as input but also as a baseline background quantity, utilizing the microenvironmental ventilation coefficient as a nonlinear gain regulator to calculate the potential dust load for each pixel.

[0048] The second stage: The solution system for the nonlinear mapping model performs nonlinear mapping calculations for aerosol blade dust accumulation for each spatial cell. This calculation uses an exponential decay or enrichment function to characterize the nonlinear control effect of ventilation capacity on dust accumulation. The core calculation formula is as follows:

[0049] ;

[0050] in, This represents the potential dust load, with physical units normalized to the particulate matter cover mass per unit leaf area (e.g., ...). This parameter represents the maximum theoretical dust accumulation without considering strong wind cleaning. This indicates the aerosol optical thickness value corresponding to that pixel. This represents the microenvironment ventilation coefficient output in step S200. This represents the conversion factor used to convert dimensionless optical thickness to a mass concentration standard. This factor can be set as an empirical constant based on the physicochemical properties of particles in the study area (such as particle size distribution and hygroscopic growth factor), or determined based on regression analysis of PM2.5 concentrations and concurrent AOD values ​​from ground-based monitoring stations. The numerical stability constant (e.g., with a value of 10). -6 This is used to prevent the denominator from becoming invalid due to the ventilation coefficient approaching zero in extremely calm weather or completely enclosed spaces.

[0051] Phase 3: Generation of Spatial Heterogeneity Features

[0052] Based on the above calculations, the system generates a potential dust load distribution map.

[0053] This distribution map exhibits significant spatial heterogeneity: under the same AOD background value, vegetation pixels located at the bottom of street canyons, on the leeward side of buildings, and in other ventilation dead zones, show different spatial heterogeneity. The term value is large, and the calculated value is... Significantly higher than background values, reflecting the hotspot accumulation effect of pollutants; while vegetation pixels located in open parks or on windward slopes, their... The value is relatively large, and the exponential term is close to 1, reflecting the normal background settlement level.

[0054] The output of this step The data will be used as an intermediate variable and directly input into the subsequent resuspension correction step (S302) to ultimately determine the total amount of particulate matter actually remaining on the blade surface.

[0055] See attached document Figure 2 Step S302 aims to introduce boundary layer shear theory from fluid mechanics to dynamically correct the potential dust load calculated in step S301, eliminating the particulate resuspension component caused by strong wind turbulence. This step specifically includes the following technical execution processes:

[0056] Phase 1: Setting Dynamic Parameters and Thresholds

[0057] The system first reads the global friction velocity calculated in step S203. ) Raster data and potential dust load output from step 5301 Raster data.

[0058] The system defines the critical frictional velocity threshold for particulate resuspension. This threshold represents the minimum turbulent shear velocity required to overcome the van der Waals forces and capillary adhesion forces between particles and the leaf surface. This threshold, as a system preset parameter, is set according to the main vegetation type (e.g., broadleaf, coniferous) and particle size characteristics of the target area. For example, a lower threshold is set for broadleaf tree species with relatively smooth surfaces, and a higher threshold is set for leaves with rough surfaces or hairy surfaces.

[0059] Phase Two: Pixel-level Re-floating Detection and Correction Calculation

[0060] The system iterates through each pixel and calculates the friction velocity value of that pixel. ) and the preset critical friction speed threshold ( Numerical comparisons are performed, and piecewise function calculations are executed based on the comparison results to output the actual blade dust accumulation load. .

[0061] The specific calculation logic follows the following physical rules:

[0062] Static dwell mode: When the friction velocity of a pixel is less than or equal to a critical threshold When the system determines that the turbulent shear force in the current microenvironment is insufficient to overcome the adhesion of particulate matter, the particulate matter on the blade surface remains in an accumulated state without significant resuspension or loss. The system directly assigns the potential dust load calculated in step S301 to the actual blade dust load.

[0063] Dynamic cleaning mode: When the friction speed of a pixel exceeds a critical threshold ( When the turbulent shear force strips away some of the settled particles, the system determines that wind-induced resuspension has occurred. At this point, the system uses an exponential decay model to calculate the remaining actual dust accumulation; the degree of decay is positively correlated with the extent to which the friction velocity exceeds a threshold.

[0064] The core formula is described as follows:

[0065] ;

[0066] in, This indicates the corrected actual dust accumulation load on the blades. This indicates the potential dust load input in step S301. This represents the pixel friction velocity calculated in step S203. This represents the critical frictional velocity threshold for particulate resuspension. This represents the resuspension attenuation coefficient, used to quantify the sensitivity of dust load to excess shear force.

[0067] Phase 3: Output of the actual dust accumulation load field

[0068] Based on the above calculations, the system generates the actual blade dust load ( Spatial distribution data of the wind. In this dataset, areas initially identified as highly polluted due to high AOD in step S301 will have their dust load values ​​significantly reduced if they are also located in highly turbulent areas with high friction velocities (such as vents and windward slopes); while areas with low friction velocities (such as leeward sides and vegetation interiors) will retain their dust load values. This result truly reflects the dual role of wind in particulate matter deposition: acting as both a transport carrier and, under specific conditions, a cleaning force. This data will serve as the final physical stress input, passed to the subsequent stomatal blockage factor calculation step.

[0069] Step S303 aims to base the actual blade dust load output in step S302 on the actual blade dust accumulation load. To construct a dimensionless stress index—the stomatal blockage stress factor ( This factor is used to quantify the degree to which particulate matter coverage physically obstructs the gas exchange channels of leaf stomata. This step specifically includes the following technical procedures:

[0070] Phase 1: Setting the characteristic parameters of particulate matter-stomatal interaction

[0071] The system first reads the actual dust accumulation load on the blades ( Data. To accurately describe the differentiated impact of dust accumulation on different vegetation types, the system introduces a leaf surface roughness sensitivity coefficient (FGSC). Sensitivity coefficient ( Stomatal density (SHD) is an empirical parameter related to the microstructure of plant leaves. The value of this parameter depends on the stomatal density, stomatal size, and microtexture of the leaf surface (such as the presence of pubescence and the thickness of the waxy layer) of the dominant tree species in the target area.

[0072] Specifically, for tree species with small stomata that are easily embedded by fine particles, or tree species with abundant leaf hairs that easily capture particles, a higher setting is adopted. The value represents the extreme sensitivity of stomatal conductance to dust accumulation; for tree species with large stomata or smooth leaf surfaces where particles easily slip off, a lower value is set. The system obtains this coefficient by consulting a plant physiology database or a preset vegetation type lookup table.

[0073] Phase 2: Nonlinear calculation of the stomatal blockage stress factor

[0074] The system is constructed based on Langmuir adsorption-type or saturated growth-type nonlinear response models to calculate the stomatal blockage stress factor.

[0075] This model reflects the physical process that as the amount of dust on the leaf surface increases, the proportion of stomata that are blocked gradually increases but eventually tends to saturate.

[0076] The core calculation formula is as follows:

[0077] ;

[0078] in, (Stomatal Blocking Stress Factor) represents the stomatal blockage stress factor, with its value strictly limited to the interval [0,1). When the value approaches 0, it indicates that the blade surface is clean and there is no pore blockage; when... When the value approaches 1, it indicates that the blade surface is completely covered by particulate matter, and the stomatal channels are in a state of maximum resistance. The actual blade dust load calculated in step S302, For the blade surface roughness sensitivity coefficient defined above, the fractional term in the formula... The physical meaning is the effective pore conductivity, which is the proportion of pores that are not blocked by particulate matter.

[0079] Phase 3: Spatialized Output of Stress Factors

[0080] Based on the above calculations, the system generates raster data of stomatal blockage stress factors for the entire region. This data visually represents the mass concentration ( Converted to physiological resistance gain ratio ( .

[0081] This factor does not directly represent the final resistance value, but rather acts as an independent external forcing term. In subsequent steps, it will act as a multiplicative factor on the minimum stomatal resistance, which was originally controlled only by light, temperature, and water. This mathematically embeds the physical blockage effect into the plant physiological and ecological model, providing the necessary correction coefficient for the final accurate inversion of SO₂ uptake flux under pollution stress.

[0082] Step S304 aims to comprehensively consider the physiological regulatory effects of meteorological environmental factors on stomatal opening and closing, as well as the physical blocking effect of particulate matter accumulation on stomata, to calculate the final stomatal resistance of the canopy surface. This step specifically includes the following technical execution processes:

[0083] Phase 1: Calculation of Basic Physiological Response Function

[0084] The system reads the regional meteorological data (photosynthetically active radiation, air temperature, relative humidity) and vegetation parameters (leaf area index) obtained in step S100.

[0085] Based on the Jarvis multiplicative model framework, the system first calculates the three main environmental stress functions that limit stomatal conductance. These functions all have a range of [0,1], representing the opening ratio of the stomata relative to their maximum opening state under specific environmental conditions.

[0086] Light response function ( The calculations were performed using a rectangular hyperbolic model. As photosynthetically active radiation (PAR) increases, stomatal conductance increases rapidly and tends to saturate.

[0087] Temperature response function ): Calculated using a bell-shaped curve model. Porous conductance at the optimum temperature ( It reaches its peak at a certain time and decreases rapidly under low or high temperature stress.

[0088] Moisture response function ( The calculations were performed using linear or nonlinear decay models. As the atmospheric saturated vapor pressure difference (VPD) increases, plants will actively close their stomata to prevent excessive water loss, leading to a decrease in the function value.

[0089] Phase Two: Construction of the Physical-Physiological Coupling Modified Model

[0090] The traditional Jarvis model only considers the aforementioned physiological responses and assumes the leaf surface is clean. This embodiment, building upon this, introduces the stomatal blockage stress factor calculated in step S303 (…). ) as a physical correction term.

[0091] The system's coupling logic is defined as follows: the physical blocking effect of particulate matter is equivalent to increasing the basic diffusion resistance on the blade surface. Therefore, the system... Minimum stomatal resistance of plants ( Gain correction is performed.

[0092] The corrected logic indicates that: when When =0 (clean), resistance is only physiologically controlled; when When the concentration increases, even under suitable environmental conditions (such as sufficient light and suitable temperature), the effective pore resistance will be forced to increase, thereby reducing the SO2 absorption rate.

[0093] Phase 3: The final solution system for canopy stomatal resistance will incorporate physiological parameters (leaf area index). ), environmental stress function ( ) and physical stress factors ( Substituting the improved resistance formula, the final canopy stomatal resistance is calculated. ).

[0094] The core calculation formula is as follows:

[0095] ;

[0096] in, The corrected stomatal resistance of the canopy. This is the minimum stomatal resistance (i.e., the reciprocal of the maximum stomatal conductance) for a specific vegetation type. This parameter is obtained through a vegetation type lookup table. Leaf area index, used to extend single-leaf resistance to the canopy scale. , , These are the normalized physiological regulatory functions for light intensity, temperature, and water vapor pressure difference, respectively. The physical blocking correction term introduced in this invention, wherein This is the blocking factor output in step S303.

[0097] Fourth stage: Generation of spatiotemporal dynamic resistance field

[0098] Through the above calculations, the system generated the canopy stomatal drag field during the evaluation period. The results exhibit significant spatiotemporal coupling characteristics: in the temporal dimension, the drag shows a U-shaped or dynamic change with sunrise / sunset and weather variations; in the spatial dimension, the drag is affected not only by vegetation density (… The impact is compounded by the dust accumulation and blockage effect caused by the urban form. For example, in busy and poorly ventilated street valleys, even with lush vegetation, the calculated... Also because The significantly increased value accurately reflects the reduced purification capacity of vegetation in the area due to dust accumulation. This parameter is the core denominator term in subsequent calculations of dry deposition velocity.

[0099] Step S400 aims to construct a concentration field that reflects the uneven distribution of pollution within the city. This concentration field will serve as a direct input variable (SinkConcentration) for subsequent calculations of the purification flux. This step specifically includes the following technical procedures:

[0100] Phase 1: Spatial Feature Extraction and Classification of Emission Sources

[0101] The system first calls upon the auxiliary geographic information data obtained in step S100 or performs land use classification on the high-resolution remote sensing imagery. Based on spectral and texture features, the system automatically identifies and extracts the vector boundaries of the main sulfur dioxide emission sources within the study area.

[0102] The system classifies emission sources into two categories for treatment:

[0103] Linear emission sources: mainly include urban arterial roads, expressways, and elevated bridges. The system extracts their centerlines as vector skeletons.

[0104] Area-based or point-based emission sources: mainly including industrial parks, thermal power plants, and large heating stations. The system extracts their polygonal outlines or centroid coordinates. The system assigns relative intensity coefficients to different types of emission sources. This coefficient is set based on the emissions inventory database, or based on an empirical weight value set according to the scale of the emissions source (such as road grade, factory area), and is used to distinguish the differences in the potential contribution of different sources to the surrounding environment.

[0105] Phase Two: Construction of the Source-SinkWeightField

[0106] The system employs a distance-attenuation-based geographic weighting algorithm to calculate the degree to which each pixel in the entire region is affected by emission sources.

[0107] For each pixel to be evaluated (i.e., sink) within the study area, the system iterates through all emission sources within its defined radius. The system calculates the Euclidean distance from the center of the pixel to each emission source.

[0108] A decay function is constructed based on the inverse distance law, and the source nearest neighbor weight value at this pixel is calculated. The calculation logic is as follows: the closer a pixel is to the emission source, the higher its weight value; the farther away it is, the lower the weight value becomes, decreasing exponentially. The system linearly superimposes the multi-source influences on the pixel to obtain the comprehensive weight.

[0109] The core computational logic is described as follows:

[0110] ;

[0111] in, For the first The source nearest neighbor weights of each pixel, This represents the total number of emission sources within the search radius. For the first Intensity coefficient of each emission source For pixels To the The distance between each emission source It serves as a smoothing factor to prevent computational overflow when the distance approaches zero.

[0112] Third stage: Interpolation generation of background concentration field

[0113] The system reads hourly sulfur dioxide concentration observations from ground environmental monitoring stations within and around the region. Using either Ordinary Kriging or Inverse Distance Weighted (IDW) interpolation, the discrete station data is mapped to a continuous grid surface covering the entire area. The interpolation result generated in this step is recorded as the baseline concentration field. This concentration field only reflects the regional-scale pollution background gradient, and is relatively smooth in space, failing to reflect the local high-value characteristics near roads or industrial areas.

[0114] Phase 4: Spatial Reconstruction and Correction of the Effective Exposure Concentration Field

[0115] The system utilizes the source nearest neighbor weight field generated in the second stage ( The basic concentration field generated in the third stage ( Pixel-level spatial non-uniformity correction is performed. The correction logic is based on the principle of relative stripping and redistribution: under the premise that the base concentration remains constant or trend is consistent, the concentration value of pixels closer to the emission source is increased, and the concentration value of pixels farther away from the emission source is decreased.

[0116] The system first calculates the average weight value of the entire domain ( Then execute the following reconstruction formula:

[0117] ;

[0118] in, The reconstructed effective exposure concentration of sulfur dioxide is the actual environmental concentration present in the plant leaves. The baseline concentration obtained from Kriging interpolation. This is an adjustment coefficient (typically between 0.1 and 0.5) used to control the magnitude of the source nearest neighbor effect's correction to the background field.

[0119] Phase 5: Output and Verification of Concentration Field

[0120] Through the above steps, the system outputs the final result. Raster data. In this dataset, distinct concentration ridges or high points are observed along major traffic arteries and around industrial areas, while concentration depressions are observed within urban parks far from emission sources. This distribution characteristic better reflects the actual physical dispersion patterns of urban air pollutants, resolving the peak-shaving and valley-filling errors caused by traditional single interpolation methods. This ensures that subsequent calculations of plant uptake fluxes can be performed with matching high-concentration boundary conditions for highly exposed vegetation.

[0121] Step S500 aims to quantify the actual removal capacity of urban forest vegetation for atmospheric sulfur dioxide (SO2), outputting purification flux data with high spatiotemporal resolution and regional statistical totals. This step specifically includes the following technical execution processes:

[0122] Phase 1: Quasi-Laminar Boundary Layer Resistance The physical calculation system first calculates the quasi-laminar boundary layer drag on the blade surface. This drag characterizes the resistance required for SO2 molecules to pass through a thin layer of stationary or laminar air close to the blade surface, and its magnitude is mainly controlled by the molecular diffusion rate and the intensity of micro-turbulence.

[0123] The friction velocity calculated by the system call step S203 is ( The dynamic input is used as the calculation logic. The calculation logic is based on the fluid dynamics boundary layer theory: the greater the friction velocity, the stronger the turbulent exchange near the blade surface, the thinner the quasi-laminar boundary layer, and thus the smaller the drag.

[0124] Simultaneously, the system incorporates molecular characteristic parameters of SO2 gas, including the Schmitt number ( The ratio of momentum diffusion to molecular diffusion (Pr) and the Zenont number (Pr) are used to calculate the system using empirical formulas. The formula shows It is proportional to the reciprocal of the friction velocity and is subject to nonlinear correction for the molecular diffusion coefficient. This step ensures that the model can distinguish the physical differences between different gases (such as SO2 and PM2.5) as they cross the microboundary layer of the leaf surface.

[0125] Second stage: Dry settlement velocity ( Multi-resistance synthesis

[0126] The system adopts a series resistance model architecture to synthesize the dry deposition rate of sulfur dioxide (S₂O₂). This model decomposes the atmospheric transport process to the vegetation surface into three consecutive resistance stages.

[0127] The system performs the following resistance superposition calculation:

[0128] Input aerodynamic drag ( Read the calculation results from step S203, which represent the turbulent resistance to the transport of pollutants from the reference height to the top of the canopy.

[0129] Input quasi-laminar boundary layer resistance ( : Read the calculation results of the first stage of this step, which represents the diffusion resistance of pollutants through the thin leaf surface layer.

[0130] Input the corrected canopy stomatal resistance ( ): Read the calculation results from step S304, representing the physiological and physical mixed resistance to the absorption of pollutants through stomata into mesophyll cells. Note that here... This already includes the physical blockage effect caused by particulate matter accumulation. ).

[0131] The system calculates the total settlement velocity based on the principle of the inverse of resistance. Because... , , The transport path is a series relationship (assuming cuticle deposition is negligible or treated as a parallel high-resistance term), and the dry deposition velocity is calculated using the following formula:

[0132] ;

[0133] The physical meaning of this formula is that the settling velocity is determined by the maximum resistance term along the transport path. In urban environments, this invention introduces particulate stress correction. Typically increases significantly, thus leading to the calculated This is lower than that of traditional models, which is more in line with the objective fact that plant stomatal efficiency is reduced under the background of high particulate matter pollution.

[0134] Phase 3: Instantaneous calculation of single-pixel purification flux

[0135] The system combines the sink's absorption capacity with the source's exposure level to calculate the instantaneous purification flux per pixel. .

[0136] The system reads the SO2 effective exposure concentration field reconstructed in step S400. ) and the dry settling velocity calculated in this step ( .

[0137] The system performs a multiplication operation on each raster cell. The calculation logic is: the purification flux equals the product of the settling velocity and the local effective concentration.

[0138] ;

[0139] The calculation results generated a map showing the distribution of purification efficiency. This map not only reflects the growth status of the vegetation itself (through...) It is manifested in the fact that it is also coupled with the proximity effect of pollution sources (through...). (Embody). For example, located next to a major traffic artery (high... And it grows well and has good ventilation (high) Protective forest belts will show extremely high purification flux values; while vegetation located in clean areas or areas with severe dust accumulation will have correspondingly lower flux values.

[0140] Phase 4: Spatiotemporal Integral and Regional Total Inversion

[0141] Finally, the system performs spatiotemporal integration to convert the instantaneous flux into the total purification volume within the evaluation period. In the spatial dimension, the system multiplies the flux value of each pixel by the actual physical area of ​​that pixel (e.g., 100 square meters at 10m resolution).

[0142] In terms of time, the system evaluates all time steps within the assessment period (such as one month or one year). (Accumulate)

[0143] ;

[0144] in, This represents the total number of vegetation pixels within the region. For time step, This represents the number of vegetation pixels within the region. This represents the total number of time steps. The pixel area.

[0145] The system ultimately outputs the total SO2 retention of urban forests in the region during the assessment period, and can be classified and statistically output according to administrative divisions, functional zones or vegetation types, providing quantitative data support for urban ecological benefit assessment and green space planning and management.

Claims

1. A method for remote sensing inversion of SO2 purification benefit of urban forest based on coupling model, characterized in that, The method comprises the following steps: Step S100: performing multi-source data acquisition and basic parameter preprocessing, acquiring multi-spectral remote sensing image data, digital surface model data, digital elevation model data, and regional meteorological and atmospheric environment monitoring data of a target area, and generating normalized digital surface model data and vegetation leaf area index data based on the data; Step S200: performing urban canopy three-dimensional aerodynamic resistance modeling and micro-environment ventilation factor calculation, extracting morphological parameters based on the normalized digital surface model data, inverting aerodynamic roughness parameters, and calculating friction velocity, aerodynamic resistance, and micro-environment ventilation coefficient in combination with the regional meteorological data; Step S300: constructing a particle gas pore collaborative stress model based on micro-environment ventilation efficiency, mapping the micro-environment ventilation coefficient to the atmospheric environment monitoring data to obtain actual leaf dust load, constructing a pore blockage stress factor based on the actual leaf dust load, and calculating a corrected canopy pore resistance in combination with meteorological environmental elements; Step S400: performing source-sink neighborhood degree-based reconstruction of a sulfur dioxide effective exposure concentration field, constructing a source neighborhood weight field based on the spatial distribution characteristics of emission sources, performing spatial non-uniformity correction on the background concentration of a ground monitoring station by using the source neighborhood weight field, and generating sulfur dioxide effective exposure concentration data; Step S500: performing purification benefit flux calculation and total amount inversion, calculating quasi-laminar boundary layer resistance by using the friction velocity calculated in step S200, calculating dry deposition velocity by comprehensively considering the aerodynamic resistance, the quasi-laminar boundary layer resistance, and the corrected canopy pore resistance, and calculating purification flux and purification total amount in combination with the sulfur dioxide effective exposure concentration data. The step of generating normalized digital surface model data and vegetation leaf area index data in step S100 specifically comprises:

2. The method according to claim 1, wherein, performing radiation calibration and atmospheric correction on the multi-spectral remote sensing image data to obtain ground reflectance, calculating normalized vegetation index based on the ground reflectance, establishing a nonlinear mapping relationship between normalized vegetation index and leaf area index, and generating leaf area index grid data; performing difference operation on the digital surface model data and the digital elevation model data under the same geographic coordinates, subtracting the pixel value of the digital elevation model at the corresponding position from the pixel value of the digital surface model, and performing extreme value filtering to generate normalized digital surface model data representing the absolute height of the ground object. The step of inverting aerodynamic roughness parameters and calculating friction velocity in step S200 specifically comprises:

3. The method according to claim 1, wherein, using a moving window algorithm to traverse the normalized digital surface model data, calculating the average obstacle height, windward area density, and sky openness in the window; establishing a morphological density model based on the windward area density, calculating zero plane displacement and aerodynamic roughness length; the zero plane displacement represents the height of the airflow being lifted, and the aerodynamic roughness length represents the frictional drag ability of the underlying surface to the airflow; ​ The friction velocity and the aerodynamic resistance are iteratively calculated by using the Monin-Obukhov similarity theory, combining the reference wind speed in the regional meteorological data, the zero plane displacement, and the aerodynamic roughness length.

4. The method according to claim 3, wherein, The micro-environment ventilation coefficient is calculated in the step S200, which specifically follows the following logic: The micro-environment ventilation coefficient is used to represent the air exchange capacity of the local canopy. The numerical value of the micro-environment ventilation coefficient is positively correlated with the reference wind speed, the sky openness, and the supplement of the windward area density. In the calculation, the reference wind speed is taken as the basic driving term, the sky openness is taken as the weight term of the vertical exchange capacity, and the supplement of the windward area density is taken as the weight term of the horizontal flow capacity, so as to generate the spatially continuous micro-environment ventilation coefficient grid data.

5. The method according to claim 1, wherein, The actual leaf dust load is calculated in the step S300, including the process of calculating the potential dust load: The aerosol optical depth data in the atmospheric environment monitoring data are obtained as the atmospheric background quantity; The micro-environment ventilation coefficient is used as a nonlinear gain adjuster to map the aerosol optical depth data; The mapping logic is that the relationship between the micro-environment ventilation coefficient and the dust accumulation effect is constructed by using an exponential decay function, in the area with low micro-environment ventilation coefficient, the local enrichment effect of pollutants is reflected by increasing the numerical value of the exponential term, and the potential dust load without considering the wind cleaning is calculated.

6. The method according to claim 5, wherein, The actual leaf dust load is calculated in the step S300, which also includes the step of performing the wind-induced resuspension correction process based on the friction velocity: The particle resuspension critical friction velocity threshold is set; the friction velocity calculated in the step S200 is compared with the particle resuspension critical friction velocity threshold on a pixel-by-pixel basis; When the friction velocity is less than or equal to the particle resuspension critical friction velocity threshold, it is determined to be in the static retention mode, and the potential dust load is directly taken as the actual leaf dust load; When the friction velocity is greater than the particle resuspension critical friction velocity threshold, it is determined to be in the dynamic cleaning mode, the amount of particle peeling caused by excessive shear force is calculated by using an exponential decay model, the potential dust load is reduced, and the actual leaf dust load is obtained.

7. The method according to claim 6, wherein, The step S300 of calculating the corrected canopy stomatal resistance specifically includes: The leaf surface roughness sensitivity coefficient is introduced, the actual leaf dust load is combined, and the stomatal blockage stress factor is calculated by using the Langmuir adsorption model; the stomatal blockage stress factor is used to quantify the physical shielding degree of the particle coverage on the leaf stomatal gas exchange channel; The photosynthetically active radiation, the air temperature, and the relative humidity data in the regional meteorological data are used to calculate the light response function, the temperature response function, and the water response function, respectively; The stomatal blockage stress factor is introduced into the multiplication resistance model as a physical correction term, the minimum stomatal resistance of the vegetation is gain-corrected by the stomatal blockage stress factor, and the corrected canopy stomatal resistance is calculated by combining the leaf area index and the response functions.

8. The method according to claim 1, wherein, The step S400 of reconstructing the effective exposure concentration field of sulfur dioxide based on the source-sink nearest neighbor degree specifically comprises: Identifying the sulfur dioxide emission sources in the research area, and assigning the intensity coefficient according to the type of the emission source; searching for all the emission sources within a preset radius for the pixel to be evaluated, calculating the Euclidean distance from the pixel to each emission source, and calculating the source nearest neighbor weight of the pixel based on the inverse distance law; Generating the basic concentration field by spatial interpolation using the observation data of the ground monitoring sites; Performing pixel-level redistribution correction on the basic concentration field using the source nearest neighbor weight, adjusting the concentration value of the pixel with the source nearest neighbor weight higher than the global average level, and adjusting the concentration value of the pixel with the source nearest neighbor weight lower than the global average level, to generate the effective exposure concentration data of sulfur dioxide.

9. The method according to claim 1, wherein, The step S500 of calculating the dry deposition velocity specifically comprises: Calculating the quasi-laminar boundary layer resistance based on the friction velocity and the characteristic parameters of the sulfur dioxide gas molecules, the quasi-laminar boundary layer resistance representing the diffusion resistance of the gas passing through the laminar air thin layer on the surface of the blade; Synthesizing the total deposition velocity by using the series resistance model, and the calculation logic is: calculating the dry deposition velocity by calculating the reciprocal of the sum of the aerodynamic resistance, the quasi-laminar boundary layer resistance and the corrected crown layer stomatal resistance.

10. The method according to claim 9, wherein, The step S500 of calculating the purification flux and the total purification amount specifically comprises: Performing pixel-by-pixel multiplication operation on the dry deposition velocity and the effective exposure concentration data of sulfur dioxide to obtain the instantaneous purification flux of a single pixel; Based on the physical area of the pixel and the evaluation time step, performing spatial region and time dimension integral operation on the instantaneous purification flux to obtain the total sulfur dioxide purification amount of the urban forest.

Images