A method for determining priority supervision areas of aviation pollution sources based on multi-level evaluation

Abstract

The application discloses a kind of based on multistage evaluation's aviation pollution source priority supervision area determination method, belong to airport environmental supervision technical field, including: first acquisition multi-source data, construct the multistage evaluation index system of distinguishing aircraft type and flight phase of aviation pollution source and complete quantization;Using time-space-level three-dimensional dynamic weight distribution method gives index dynamic weight and correction, by simulating diffusion process, obtain pollutant concentration spatial distribution;The target area is divided into adaptive spatial unit grid, by adaptive evaluation algorithm coupling concentration data and dynamic weight, obtain each grid comprehensive evaluation score;According to regional environmental standard and bearing capacity, set different threshold, generate pollution evaluation spatial distribution map by spatial interpolation, demarcate priority supervision area boundary and complete internal priority ranking, realize aviation pollution source supervision area accurate determination, improve airport pollution control efficiency and pertinence.

Description

Technical Field

[0001] This invention belongs to the field of airport environmental supervision technology, specifically a method for determining priority supervision areas for aviation pollution sources based on multi-level evaluation. Background Technology

[0002] Existing technologies cannot achieve dynamic assessment and precise spatial identification of pollution risks, making it difficult to prioritize and classify high-pollution areas for control. This results in weak regulatory targeting and low efficiency, failing to meet the needs of refined and dynamic environmental regulation in airport areas. Therefore, there is an urgent need for a method to determine priority regulatory areas for aviation pollution sources that integrates multi-source data, dynamic weights, adaptive evaluation, and differentiated judgment, in order to solve the problems of inaccurate simulation, imprecise evaluation, and extensive regulation in existing technologies. Summary of the Invention

[0003] To address the shortcomings of existing technologies, this invention proposes a method for determining priority monitoring areas for aviation pollution sources based on multi-level evaluation. First, multi-source data is acquired, and a multi-level evaluation index system for aviation pollution sources, distinguishing between aircraft type and flight phase, is constructed and quantified. A three-dimensional dynamic weight allocation method based on time, space, and hierarchy is used to assign dynamic weights to the indicators and perform collaborative correction. By simulating the diffusion process, the spatial distribution of pollutant concentrations is obtained. The target area is divided into adaptive spatial unit grids, and an adaptive evaluation algorithm couples concentration data with dynamic weights to obtain the comprehensive evaluation score for each grid. Differential thresholds are set according to regional environmental standards and carrying capacity. A pollution evaluation spatial distribution map is generated through spatial interpolation, delineating the boundaries of priority monitoring areas and completing internal priority ranking. This achieves accurate determination of aviation pollution source monitoring areas, improving the efficiency and targeting of airport pollution control.

[0004] To achieve the above objectives, the present invention provides the following technical solution:

[0005] A method for determining priority regulatory areas for aviation pollution sources based on multi-level evaluation includes:

[0006] Acquire multi-source data for the target airport area, construct a multi-level evaluation index system for aviation pollution sources that distinguishes between aircraft type and flight phase based on the multi-source data, and complete the quantification; the multi-source data includes at least aircraft operation data, pollution monitoring data, environmental sensitive point data, meteorological data, and topographic elevation data;

[0007] A time-space-hierarchical three-dimensional dynamic weight allocation method is adopted to assign dynamic weights to each level of the multi-level evaluation index system and perform collaborative correction. At the same time, QAR data and ADS-B trajectory data in the aircraft operation data are analyzed, and meteorological data and terrain elevation data are combined to simulate the dynamic diffusion process of aviation pollutants and output the spatial distribution data of pollutant concentration in the target airport area.

[0008] The target airport area is divided into spatial cell grids. An adaptive evaluation algorithm is used to couple the spatial distribution data of pollutant concentration with dynamic weights to obtain the comprehensive evaluation score of each spatial cell grid.

[0009] Based on the environmental standards and pollution carrying capacity of different sub-regions within the target airport area, differentiated pollution judgment thresholds are set. Using a spatial interpolation algorithm, a pollution assessment spatial distribution map is generated based on the comprehensive evaluation score. By comparing the comprehensive evaluation score of each location in the pollution assessment spatial distribution map with the pollution judgment threshold of its respective sub-region, the boundary of the priority supervision area is delineated. Based on the comprehensive evaluation score, the sub-regions within the boundary are internally prioritized, and the priority supervision area judgment result is output.

[0010] Specifically, the construction and quantification of a multi-level evaluation index system for aviation pollution sources that distinguishes between aircraft type and flight phase includes:

[0011] Based on aircraft operation data, the ICAO baseline emissions database is queried to obtain the baseline emission factors of various pollutants per unit time for aircraft engines under standard sea-level static conditions. The aircraft's flight trajectory is then divided into flight phases, and the baseline emission factors are corrected by incorporating temperature and pressure information from meteorological data to obtain the total emissions per flight for each aircraft type in each flight phase, thus constructing an initial pollution source intensity index. The baseline emission factors at least include the emission mass of hydrocarbons, carbon monoxide, and nitrogen oxides per unit time under standard sea-level static conditions. The flight phases include taxiing, takeoff, climb, approach, and landing.

[0012] Based on environmentally sensitive point data, a buffer analysis of environmentally sensitive points is performed using a geographic information system. According to the type of sensitive point and the distribution density of sensitive points within the buffer zone, an environmental sensitivity index is assigned to each spatial unit grid in the assessment area.

[0013] Specifically, the construction and quantification of a multi-level evaluation index system for aviation pollution sources that distinguishes between aircraft type and flight phase also includes:

[0014] Based on meteorological and topographic elevation data, a meteorological field and elevation distribution covering the target airport area are generated through spatial interpolation. Based on the influence of wind speed, atmospheric stability and topography on pollutant diffusion, diffusion condition indicators are constructed for each spatial unit grid.

[0015] The initial pollution source intensity index, environmental sensitivity index, and diffusion condition index are used as the underlying quantitative indexes. The weight of each underlying quantitative index is determined by the analytic hierarchy process (AHP). A multi-level evaluation index system for aviation pollution sources that distinguishes between aircraft type and flight phase is constructed, and the quantitative values ​​of all indicators are assigned.

[0016] Specifically, the method of employing a time-space-hierarchical three-dimensional dynamic weight allocation to assign dynamic weights to each level of indicators in the multi-level evaluation index system and performing collaborative correction includes:

[0017] Based on the timestamps of the underlying quantitative indicators in the aircraft operation data, the time weight of each indicator data is calculated using a time decay function.

[0018] Based on the spatial cell grid division results, combined with the spatial distribution of terrain elevation and environmentally sensitive point data, the local Moran index calculation method is used to calculate the spatial autocorrelation between each grid and its neighboring grids, and generate spatial weights.

[0019] Based on the hierarchical structure of the multi-level evaluation index system, the entropy weight method is used to calculate the entropy value of the index data at each level and generate the hierarchical weights.

[0020] The time weight, spatial weight, and hierarchical weight are multiplied to synthesize the initial three-dimensional dynamic weight. Based on the historical measured concentrations in the pollution monitoring data, the objective function is constructed using the least squares method. The initial three-dimensional dynamic weight is adjusted through an iterative optimization algorithm to minimize the error between the simulated pollutant concentration distribution obtained based on the three-dimensional dynamic weight and the historical measured data, and the optimized three-dimensional dynamic weight is output.

[0021] Specifically, the simulated dynamic diffusion process of aviation pollutants includes:

[0022] By analyzing QAR data and ADS-B trajectory data from aircraft operation data, the instantaneous emission rates of each pollutant at different times and locations of each aircraft are calculated, and a sequence of emission intensity of mobile point sources is generated.

[0023] The emission intensity sequence of mobile point sources, meteorological data, and terrain elevation data are input into a pre-constructed Lagrange particle diffusion model. In the Lagrange particle diffusion model, each aircraft is simulated as a point source releasing particles representing pollutants at each time step, and the particle movement is driven by the input meteorological data and terrain elevation data.

[0024] The motion trajectories of all simulated particles output by the Lagrange particle diffusion model are statistically analyzed. The pollutant mass contribution of each particle through the spatial cell grid during the simulation period is calculated. The contributions of all particles in the same grid are integrated over time and accumulated spatially to generate spatial distribution data of pollutant concentration covering the target airport area and distinguishing pollutant types.

[0025] Specifically, the calculation process for the comprehensive evaluation score includes:

[0026] Based on terrain elevation data and environmental sensitive point data, the subdivision demand index of each predetermined initial grid is calculated. The quadtree partitioning algorithm is used, with the target airport area as the root node, to recursively determine whether the grid subdivision demand index in the current area is higher than the preset first index threshold. If it is higher than the first index threshold, the area is divided into four sub-areas until the subdivision demand index of all sub-areas is lower than the first index threshold.

[0027] The generated spatial distribution data of pollutant concentration is mapped to each spatial cell grid using a bilinear interpolation method to obtain the representative value of pollutant concentration for each grid. The global pollutant concentration variation coefficient of the entire target airport area and the local pollutant concentration variation coefficient of each grid local area are calculated. The global pollutant concentration variation coefficient and the local pollutant concentration variation coefficient are compared with the preset global threshold and local threshold, respectively.

[0028] If the global pollutant concentration variation coefficient is greater than a preset global threshold, or the local pollutant concentration variation coefficient is greater than a preset local threshold, then an evaluation function based on extreme value theory is selected, with the extreme value of pollutant concentration within the grid as the evaluation benchmark. Otherwise, an evaluation function based on mean is selected, with the average value of pollutant concentration within the grid as the evaluation benchmark. Then, the representative value of pollutant concentration of each spatial unit grid and the optimized three-dimensional dynamic weight are input into the selected evaluation function to obtain the comprehensive evaluation score of the grid.

[0029] Specifically, the calculation of the subdivision demand index for each predetermined initial grid includes:

[0030] Based on terrain elevation data, the elevation standard deviation within the grid is calculated and normalized to a first exponential component. Based on environmentally sensitive point data, the weighted density of various environmentally sensitive points within the grid and adjacent grids is calculated and normalized to a second exponential component. The first exponential component and the second exponential component are weighted and summed to obtain the subdivision demand index of the grid.

[0031] Specifically, the pollution determination threshold is set as follows:

[0032] Based on environmentally sensitive point data and the functional zoning map of the target airport area, the target airport area is divided into a core protection zone, a buffer control zone, and a general impact zone. Using the environmental carrying capacity assessment method, combined with historical pollution monitoring data and meteorological statistics, the pollution carrying capacity of each sub-region is calculated. Based on the ratio of environmental quality standards to pollution carrying capacity of each sub-region, pollution judgment thresholds are set from low to high for the core protection zone, buffer control zone, and general impact zone.

[0033] Specifically, the process of generating the pollution assessment spatial distribution map includes:

[0034] Using the comprehensive evaluation score of all spatial unit grids as input variables and terrain elevation as a co-variant, a variogram model is constructed using the co-kriging spatial interpolation method. The nugget value, sill value, and range parameters in the variogram model are optimized using cross-validation. Using the optimized variogram model, interpolation calculations are performed on the continuous geographic space of the entire target airport area to generate a continuously covering evaluation score surface, which is then used as a spatial distribution map of pollution evaluation.

[0035] Specifically, the process of delineating the boundaries of priority monitoring areas by comparing the comprehensive evaluation scores of each location on the pollution assessment spatial distribution map with the pollution determination threshold of their respective sub-regions, and then ranking the internal priorities of the sub-regions within the boundaries based on the comprehensive evaluation scores, includes:

[0036] The spatial distribution map of pollution assessment is overlaid with the pollution judgment threshold in a raster overlay analysis. The comprehensive evaluation score is compared with the pollution judgment threshold of the sub-region on a raster-by-raster basis, and raster units whose comprehensive evaluation score exceeds the pollution judgment threshold are selected.

[0037] The selected grid cells are spatially clustered using the eight-neighbor connectivity analysis method. Spatially adjacent grid cells are aggregated into candidate regions. Each candidate region is filtered according to a preset minimum area threshold, and candidate regions with an area greater than or equal to the minimum area threshold are retained. The minimum bounding polygon of each retained region is extracted as the boundary of the priority monitoring region.

[0038] For each designated priority regulatory area, the comprehensive evaluation score of all spatial unit grids within its boundary is extracted. At least one statistical characteristic value among the maximum, average, or median comprehensive evaluation scores of each area is calculated. The selected statistical characteristic value is used as the sorting basis, and all priority regulatory areas are sorted in descending order to generate an internal priority ranking list. The internal priority ranking list and the geographical boundary information of the priority regulatory areas are used together as the output of the priority regulatory area determination result.

[0039] Compared with the prior art, the beneficial effects of the present invention are:

[0040] 1. This invention constructs a multi-level evaluation index system that distinguishes between aircraft type and flight phase, and combines multi-source data such as aircraft operation, pollution monitoring, environmental sensitive points, meteorology and topography to achieve accurate quantification. It breaks through the limitations of traditional evaluation that ignores differences in aircraft type, characteristics of flight phase and environmental sensitivity, and improves the comprehensiveness and scientificity of aviation pollution source assessment.

[0041] 2. By adopting a three-dimensional dynamic weight allocation method based on time, space, and hierarchy, combined with measured data for collaborative correction, and in conjunction with a multi-factor coupled diffusion model of moving point sources, the dynamic diffusion process of aviation pollutants can be realistically reproduced. This solves the problems of fixed weights and insufficient diffusion simulation accuracy in traditional methods, improves the accuracy and timeliness of spatial distribution simulation of pollutant concentrations, and enables dynamic characterization of pollution status in airport areas.

[0042] 3. By coupling weights and concentration data with adaptive grid division and adaptive evaluation algorithm, the grid fineness can be automatically adjusted according to terrain undulation and distribution of sensitive points, taking into account both computational efficiency and evaluation accuracy. At the same time, an appropriate evaluation function is selected according to regional differences, so that the comprehensive evaluation score of each spatial unit can objectively reflect the actual pollution risk and avoid misjudgment caused by uniform evaluation standards.

[0043] 4. This invention sets differentiated pollution judgment thresholds according to environmental functional zoning, and accurately delineates priority supervision areas by combining spatial interpolation and grid analysis, and completes internal priority ranking. This enables the transformation of airport pollution supervision from overall control to precise control of key areas, effectively improving the efficiency of supervision resource utilization, reducing environmental management costs, and enhancing the pertinence and effectiveness of airport area pollution prevention and control. Attached Figure Description

[0044] Figure 1 This is a schematic diagram of a method for determining priority regulatory areas for aviation pollution sources based on multi-level evaluation, according to the present invention.

[0045] Figure 2 This is a flowchart illustrating the working principle of a method for determining priority regulatory areas for aviation pollution sources based on multi-level evaluation, as described in this invention. Detailed Implementation

[0046] The technical solutions in the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments.

[0047] Please see Figure 1 and Figure 2 The present invention provides an embodiment of a method for determining priority regulatory areas for aviation pollution sources based on multi-level evaluation, comprising steps S1 to S4, wherein:

[0048] S1: Acquire multi-source data for the target airport area, construct a multi-level evaluation index system for aviation pollution sources that distinguishes between aircraft types and flight phases based on the multi-source data, and complete the quantification; the multi-source data includes at least aircraft operation data, pollution monitoring data, environmental sensitive point data, meteorological data, and topographic elevation data;

[0049] Furthermore, the acquisition of multi-source data for the target airport area specifically includes: taking a domestic trunk passenger airport as the target airport for this implementation. This airport has two parallel runways and three terminals. It is surrounded by environmentally sensitive areas such as residential areas, schools, and hospitals. It is located in a subtropical monsoon climate zone, with significant differences in wind speed, wind direction, and atmospheric stability throughout the four seasons. The terrain is mainly plains with some low hills, which meets the typical scenario of aviation pollution source diffusion simulation. The data acquisition team completed the comprehensive collection and standardized processing of multi-source data through four major data sources: the airport operation management center, local ecological environment monitoring stations, meteorological departments, and geographic information surveying and mapping units. All data adopt a unified geographic coordinate system and timestamp to ensure the spatiotemporal consistency of the analysis.

[0050] Furthermore, the aircraft operation data covers core information such as flight takeoff and landing plans, actual operating trajectories, aircraft type information, takeoff and landing numbers, taxi routes, and engine models for the airport over the past year. It covers mainstream civil aircraft types including Boeing 737, Airbus A320, Boeing 787, and Airbus A330, encompassing narrow-body, wide-body, and regional aircraft. It comprehensively records the entire operational process of each aircraft from approach taxiing, takeoff, climb, cruise, approach, landing, to departure taxiing. Pollution monitoring data comes from real-time monitoring data from five fixed environmental monitoring stations and three mobile monitoring vehicles around the airport. Monitoring indicators include typical aviation pollutants such as hydrocarbons, carbon monoxide, nitrogen oxides, and fine particulate matter. Concentration data is recorded hourly, and historical monitoring data from the past three years is retained for model calibration. The data on environmentally sensitive points was collected through a geographic information system, including coordinates, names, protection levels, population density, and area of ​​sensitive locations within a 10-kilometer radius of the airport, such as residential areas, primary and secondary schools, kindergartens, hospitals, nature reserves, and drinking water source protection areas. A total of 126 sensitive points were marked. Meteorological data consisted of parameters such as wind speed, wind direction, temperature, humidity, air pressure, precipitation, atmospheric stability, and cloud cover recorded minute by minute by the airport's meteorological station. Supplementary data from three surrounding regional meteorological stations were also integrated to form a meteorological spatiotemporal sequence covering the entire airport area. Topographic elevation data was acquired through surveying drones and satellite remote sensing imagery, with a resolution of 1 meter, accurately reconstructing the altitude, topographic relief, and land cover type of the airport and its surrounding areas, providing high-precision topographic support for the diffusion model.

[0051] This embodiment takes the west runway of the trunk airport and the surrounding 5-kilometer radius as the core evaluation area. This area is the core airspace for aircraft take-off, climb, approach and landing, as well as the main passage for ground taxiing. The pollution source emission intensity is high, the diffusion conditions are complex, and the environmentally sensitive points are dense, making it a key area for the supervision of aviation pollution sources.

[0052] The construction and quantification of a multi-level evaluation index system for aviation pollution sources that distinguishes between aircraft type and flight phase includes:

[0053] A1: Based on aircraft operation data, the ICAO baseline emissions database is queried to obtain the baseline emission factors of each pollutant per unit time of the aircraft engine under standard sea-level static conditions. The aircraft's flight trajectory is divided into flight phases, and the baseline emission factors are corrected by combining temperature and pressure information from meteorological data to obtain the total emissions per flight for each aircraft type in each flight phase, thus constructing an initial pollution source intensity index. The baseline emission factors include at least the mass of hydrocarbons, carbon monoxide, and nitrogen oxides emitted per unit time under standard sea-level static conditions. The flight phases include taxiing, takeoff, climb, approach, and landing.

[0054] In this embodiment, the data processing team, based on the collected aircraft operation data, logged into the ICAO's publicly available baseline emissions database. According to the aircraft engine model, they queried the baseline emission factors for the three core pollutants—hydrocarbons, carbon monoxide, and nitrogen oxides—per unit time under standard sea-level static conditions, clarifying the basic emission levels of different aircraft models under conditions without external environmental interference. Secondly, combining ADS-B trajectory data and QAR data, the flight trajectory of each aircraft was standardized into flight phases, precisely breaking down the entire flight process into five key phases: ground taxiing, takeoff roll, initial climb, approach, and landing. The duration, speed, altitude, and engine operating conditions of each phase were accurately extracted from the raw data. Subsequently, the baseline emission factors were dynamically corrected according to actual environmental conditions, eliminating the impact of environmental factors such as high altitude, high temperature, and low pressure on engine emissions, obtaining the actual emission factors per unit time for each aircraft model under the corresponding flight phase and meteorological conditions. Finally, by combining the actual duration of each flight phase, the total emissions of each pollutant of a single aircraft during each flight phase are calculated. Then, by classifying and statistically analyzing by aircraft type and flight phase, the total emission intensity of different aircraft types in each flight phase within the core area of ​​the target airport is obtained, forming an initial pollution source intensity index. This index can accurately distinguish the differences in pollution emissions of different aircraft types and different flight phases, solving the problem of calculating emission intensity in a one-size-fits-all manner in traditional evaluation.

[0055] For example, taking the Airbus A320 as an example, this aircraft type has the most takeoffs and landings at the target airport. During the ground taxiing phase, the engines operate at low power, resulting in a lower corrected NOx emission factor. During takeoff and initial climb, the engines operate at full load, increasing NOx and carbon monoxide emission factors. During approach and landing, engine power decreases, and emission intensity decreases accordingly. Through the above calculations, the total emission differences of this aircraft type in the five flight phases can be accurately obtained, providing accurate basic data on pollution source intensity for subsequent evaluations.

[0056] Furthermore, the specific steps of A1 include:

[0057] First, acquire aircraft operation data, ICAO baseline emissions database, and meteorological data. The aircraft operation data includes Fast Access Recorder (FAR) data and Automatic Dependent Surveillance-Broadcast (ADS-B) trajectory data. Extract the baseline emission factors of the corresponding aircraft type's engine under standard sea-level static conditions from the baseline emissions database. The baseline emission factors include the emission mass of hydrocarbons, carbon monoxide, and nitrogen oxides per unit time under standard sea-level static conditions. Perform unified spatiotemporal matching of aircraft operation data, baseline emission factors, and meteorological data to form the basic input data for trajectory analysis and emission calculation.

[0058] Second, based on the fusion of data from the Fast Access Recorder (HARCR) and Automatic Dependent Surveillance-Broadcast (ADS-B) trajectory, and according to the aircraft's actual motion state and engine operating characteristics within the airport area, the aircraft's trajectory is divided into five interconnected and non-overlapping flight phases, using the aircraft's takeoff time, touchdown time, altitude threshold, speed threshold, and thrust status as the criteria: ground taxiing phase, takeoff roll phase, initial climb phase, approach phase, and landing phase. The ground taxiing phase is defined as an aircraft altitude below two meters, a speed below 10 kilometers per hour, and engine thrust at idle. The takeoff roll phase is defined as an aircraft altitude increasing from two meters to 10 meters, and a speed increasing from 10 kilometers per hour to takeoff speed. The engine thrust is at maximum takeoff thrust. The initial climb phase is determined when the aircraft's altitude increases from 10 meters to 300 meters, the speed is maintained at the climb rate, and the engine thrust is at the climb thrust level. The approach phase is determined when the aircraft's altitude decreases from 300 meters to 10 meters, the speed gradually decreases to the landing approach speed, and the engine thrust is between idle and approach thrust. The landing phase is determined when the aircraft's altitude is below 10 meters, the speed is below the landing approach speed until it comes to a complete stop, and the engine thrust is at idle. All division conditions use fixed numerical standards to ensure that the phase division results are consistent for the same aircraft type and the same operating scenario. After the division is completed, the timestamp interval, altitude sequence, speed sequence, and thrust status information for each flight phase are output.

[0059] Third, temperature and pressure information that precisely match the timestamp intervals of each flight phase are extracted from meteorological data. The temperature information uses the actual atmospheric temperature at the location of the aircraft, and the pressure information uses the actual atmospheric static pressure at the location of the aircraft. The acquisition step size for both temperature and pressure data is set to one second. The average value of the continuous temperature and pressure sequences corresponding to each flight phase is calculated to obtain the average temperature and average static pressure values ​​for each flight phase, which are used as engine environmental parameters.

[0060] Fourth, using the standard sea level temperature of 15 degrees Celsius and the standard sea level hydrostatic pressure of 1013.25 hPa as benchmarks, the temperature difference between the average temperature value and the standard sea level temperature value is calculated. Temperature correction factors of 0.002, 0.0015, and 0.003 are applied to hydrocarbons, carbon monoxide, and nitrogen oxides, respectively. The temperature difference is multiplied by the temperature correction factor of the corresponding pollutant and then added to the benchmark emission factor to obtain the temperature-corrected intermediate emission factor.

[0061] Fifth, calculate the pressure ratio between the average static pressure value and the standard sea level static pressure value, and use a pressure correction index of 0.7. Then, calculate the intermediate emission factor after temperature correction and the pressure ratio using the pressure correction index to obtain the intermediate emission factor after temperature and pressure correction.

[0062] Sixth, a fixed operating condition weighting coefficient is configured based on the engine thrust percentage for each flight phase. Specifically, the thrust percentage for the ground taxiing and landing phases is 5%, and the operating condition weighting coefficient is 0.1; the thrust percentage for the takeoff roll phase is 100%, and the operating condition weighting coefficient is 1.0; the thrust percentage for the initial climb phase is 95%, and the operating condition weighting coefficient is 0.95; and the thrust percentage for the approach phase is 30%, and the operating condition weighting coefficient is 0.3. The intermediate emission factor after temperature and pressure correction is multiplied by the operating condition weighting coefficient for the corresponding flight phase to obtain the actual emission factor per unit time for the corresponding aircraft type and flight phase.

[0063] Seventh, calculate the duration based on the start and end timestamps of each flight phase, convert it to hourly units, and multiply the actual emission factor per unit time by the duration of the flight phase to obtain the total hydrocarbon emissions, total carbon monoxide emissions, and total nitrogen oxide emissions of the corresponding aircraft type in a single flight phase.

[0064] Eighth, fixed integration weights are configured according to the degree of environmental impact of aviation pollution sources, with hydrocarbons having a weight of 0.2, carbon monoxide having a weight of 0.3, and nitrogen oxides having a weight of 0.5. The total emissions of hydrocarbons, carbon monoxide, and nitrogen oxides are multiplied by their respective weights and then summed to obtain the total comprehensive emissions of the aircraft type during the flight phase, which is the quantitative value of the initial pollution source intensity index.

[0065] A2: Based on environmental sensitive point data, a buffer analysis of environmental sensitive points is performed using a geographic information system. According to the type of sensitive point and the distribution density of sensitive points within the buffer zone, an environmental sensitivity index is assigned to each spatial unit grid in the assessment area.

[0066] Furthermore, the specific steps of A2 include:

[0067] First, obtain all environmentally sensitive point data within the airport assessment area from the ecological and environmental authorities and the geographic information public platform. The environmentally sensitive point data includes the spatial coordinates of the sensitive points, the type of sensitive points, the area they occupy, the number of people inside, and the population density of the surrounding area. Unify all sensitive point data under the same geographic coordinate system, complete the standardization of data format, remove duplicate records and invalid information, and form a standardized environmentally sensitive point dataset that can be used for spatial analysis.

[0068] Second, based on the technical specifications for environmental impact assessment and the requirements for environmental management around the airport, standardized environmental sensitive points are divided into three fixed levels. The first level is core sensitive points, including nature reserves, drinking water source protection areas, hospitals, kindergartens, and primary and secondary schools. The second level is important sensitive points, including concentrated residential areas, large office buildings, and cultural and sports venues. The third level is general sensitive points, including small shops, scattered residential areas, and ordinary green areas. After classification, classified environmental sensitive point data with level labels are output.

[0069] Third, activate the spatial analysis tool of the geographic information system, import the classified environmental sensitive point data with grade labels, and set differentiated buffer ranges for sensitive points of different grades. For first-level sensitive points, set a 500-meter core buffer and a 1000-meter impact buffer; for second-level sensitive points, set a 300-meter core buffer and an 800-meter impact buffer; and for third-level sensitive points, set a 200-meter core buffer and a 500-meter impact buffer. Generate a planar buffer with the spatial coordinates of each sensitive point as the center. The core buffer and the outer impact buffer of the same sensitive point are connected to each other and do not overlap. Spatially overlay the buffers of all sensitive points to form a multi-level environmental sensitive buffer layer covering the entire assessment area.

[0070] Fourth, based on the complete boundary of the airport assessment area, a regular grid of uniform size is divided. The length and width of a single spatial unit grid are set to 100 meters × 100 meters. All grids are square and there are no gaps or overlaps between them. After the grid is divided, a unique spatial number is assigned to each spatial unit grid, and the vertex coordinates and center coordinates of each grid are recorded to form a spatial unit grid layer.

[0071] Fifth, the multi-level environmental sensitive buffer layer and the spatial unit grid layer are spatially overlaid and matched. Each spatial unit grid is traversed one by one to determine the buffer type and sensitive point level of the grid. The number of sensitive points of each level within the grid is counted. A fixed number of weighting values ​​are set for sensitive points of different levels. The weighting value of the first level core sensitive point is set to 3, the weighting value of the second level important sensitive point is set to 2, and the weighting value of the third level general sensitive point is set to 1. The number of sensitive points of each level in the grid is multiplied by the corresponding weighting value and summed. Then, the result is divided by the grid area to obtain the weighted distribution density of sensitive points in the grid. The grid density calculation of the entire area is completed.

[0072] Sixth, based on the sensitivity level and buffer range of the spatial unit grid, a fixed initial sensitivity score is assigned. Among them, the initial sensitivity score is set to 90 points for grids that are completely within the core buffer of the first-level core sensitivity point; 80 points for grids that are within the outer influence buffer of the first-level core sensitivity point or the core buffer of the second-level important sensitivity point; 60 points for grids that are within the outer influence buffer of the second-level important sensitivity point or the core buffer of the third-level general sensitivity point; 40 points for grids that are within the outer influence buffer of the third-level general sensitivity point; and 20 points for grids that do not fall into any sensitivity buffer.

[0073] Seventh, set fixed density correction segmentation parameters and adjust the initial sensitivity score. For grids with a weighted distribution density of 5 sensitive points per hectare or more, add 10 points to the initial score. For grids with a weighted distribution density of 2 sensitive points per hectare or more but less than 5 sensitive points per hectare, add 5 points to the initial score. For grids with a weighted distribution density of 0 sensitive points but less than 2 sensitive points per hectare, keep the initial score unchanged. For grids with no sensitive points, deduct 10 points from the initial score. After correction, the corrected sensitivity score of each spatial unit grid is obtained.

[0074] Eighth, the modified sensitivity scores of all spatial unit grids are mapped to the range of 0 to 100 using a fixed linear transformation rule to complete the standardization process. The standardized scores are the environmental sensitivity index of the corresponding spatial unit grid. The higher the index value, the higher the environmental sensitivity of the area.

[0075] A3: Based on meteorological and topographic elevation data, a meteorological field and elevation distribution covering the target airport area are generated through spatial interpolation. Based on the influence of wind speed, atmospheric stability and topography on pollutant diffusion, diffusion condition indicators are constructed for each spatial unit grid.

[0076] Furthermore, the specific steps in A3 include:

[0077] First, measured meteorological data of the target airport area are obtained from airport meteorological stations and local meteorological departments. The meteorological data includes the location of the observation point, wind speed, wind direction, atmospheric stability level, temperature and air pressure information. The data collection time interval is uniformly set to one hour. Topographic elevation data of the target airport area is obtained from surveying and mapping departments and satellite remote sensing data. The topographic elevation data includes ground elevation, topographic relief and distribution information of surface obstacles. The spatial resolution of the data is set to 10 meters. The meteorological data and topographic elevation data are unified under the same geographic coordinate system, and the data format is standardized. Abnormal records and missing data are removed to form a standardized meteorological dataset and a standardized topographic elevation dataset.

[0078] Second, the standardized meteorological dataset is processed using the inverse distance weighted interpolation method. The interpolation search radius is set to 500 meters, and the number of nearest neighbor meteorological observation points participating in the interpolation calculation is set to 6. First, the wind speed data is subjected to global spatial interpolation calculation to generate continuous wind speed spatial distribution data covering the entire target airport area. Then, the atmospheric stability data is subjected to global spatial interpolation calculation. The atmospheric stability is numerically expressed according to the six-level classification standard, and each level uses a fixed numerical label to generate continuous atmospheric stability spatial distribution data. The wind speed spatial distribution data and the atmospheric stability spatial distribution data are integrated to form a meteorological field that completely covers the target airport area. The inverse distance weighted interpolation method is existing technology in this field and is not an inventive solution of this application, so it will not be described in detail here.

[0079] Third, the standardized terrain elevation dataset is processed using the Kriging interpolation method. The range parameter of the interpolation model is set to 1000 meters, the nugget value parameter is set to 0.5, and the pedestal value parameter is set to 2.5. Based on the original terrain elevation data with a resolution of 10 meters, continuous interpolation calculations are performed throughout the target airport area to generate smooth and complete elevation distribution data, which is used to reflect the elevation changes and terrain undulations in the area. The Kriging interpolation method is existing technology in this field and is not an inventive solution of this application, so it will not be described in detail here.

[0080] Fourth, based on the completed spatial unit grid layer, the spatial distribution data of wind speed, atmospheric stability and elevation in the meteorological field are matched with each spatial unit grid. The wind speed, atmospheric stability and elevation values ​​corresponding to the center of each spatial unit grid are extracted as the meteorological and topographic basic parameters of the grid, ensuring that each grid has unique basic parameters.

[0081] Fifth, wind speed diffusion contribution scores are determined for each spatial unit grid according to a fixed segmented assignment standard. Specifically, the wind speed diffusion contribution score is set to 20 points for grids with wind speed values ​​less than or equal to 1 meter per second, 40 points for grids with wind speed values ​​greater than 1 meter per second and less than or equal to 3 meters per second, 60 points for grids with wind speed values ​​greater than 3 meters per second and less than or equal to 5 meters per second, 80 points for grids with wind speed values ​​greater than 5 meters per second and less than or equal to 8 meters per second, and 100 points for grids with wind speed values ​​greater than 8 meters per second.

[0082] Sixth, a stability diffusion contribution score is determined for each spatial cell grid according to a fixed segmented assignment standard. Specifically, the stability diffusion contribution score is set to 20 points for the grid with the worst atmospheric stability level and the middle ...

[0083] Seventh, according to a fixed segmented assignment standard, determine the terrain diffusion contribution score for each spatial unit grid. For grids with a terrain elevation standard deviation greater than or equal to 50 meters and containing closed depressions or tall obstacles, the terrain diffusion contribution score is set to 20 points; for grids with a terrain elevation standard deviation greater than or equal to 30 meters and less than 50 meters, the terrain diffusion contribution score is set to 40 points; for grids with a terrain elevation standard deviation greater than or equal to 10 meters and less than 30 meters, the terrain diffusion contribution score is set to 60 points; for grids with a terrain elevation standard deviation greater than or equal to 5 meters and less than 10 meters, the terrain diffusion contribution score is set to 80 points; and for grids with a terrain elevation standard deviation less than 5 meters and flat and open terrain, the terrain diffusion contribution score is set to 100 points.

[0084] Eighth, set fixed weight parameters: the weight of wind speed diffusion contribution score is set to 0.4, the weight of atmospheric stability diffusion contribution score is set to 0.3, and the weight of topographic diffusion contribution score is set to 0.3. Multiply the three types of diffusion contribution scores of each spatial cell grid by their corresponding weights, and then sum the results to obtain the initial diffusion condition score of the grid.

[0085] Ninth, the initial diffusion condition scores of all spatial unit grids are uniformly mapped to a numerical range of zero to one hundred using a fixed linear transformation rule to complete the standardization process. The final score after standardization is the diffusion condition index of the current spatial unit grid. The higher the index value, the stronger the regional pollutant diffusion capacity; the lower the index value, the easier it is for pollutants to accumulate and remain.

[0086] A4: The initial pollution source intensity index, environmental sensitivity index, and diffusion condition index are used as the underlying quantitative indexes. The weight of each underlying quantitative index is determined by the analytic hierarchy process (AHP). A multi-level evaluation index system for aviation pollution sources that distinguishes between aircraft type and flight phase is constructed, and the quantitative values ​​of all indicators are assigned. The AHP is a prior art in this field and is not an inventive solution of this application. It will not be elaborated here.

[0087] Furthermore, the specific steps for A4 include:

[0088] First, obtain initial pollution source intensity indicators, environmental sensitivity indicators, and diffusion condition indicators;

[0089] Second, a hierarchical framework of the analytic hierarchy process is constructed. In the three-level multi-level evaluation index system for aviation pollution sources, the target layer is the comprehensive pollution risk of aviation pollution sources, the criterion layer includes three categories: pollution source emission characteristics, environmental sensitivity, and diffusion and transmission conditions, and the index layer includes three specific bottom-level quantitative indicators: initial pollution source intensity index, environmental sensitivity index, and diffusion condition index. The entire hierarchical structure is unidirectionally related from top to bottom, with no overlap or loop.

[0090] Third, a judgment matrix is ​​constructed using a scale method from one to nine. Scale value one represents that the two indicators are equally important, scale value three represents that the former indicator is slightly more important than the latter, scale value five represents that the former indicator is significantly more important than the latter, scale value seven represents that the former indicator is strongly more important than the latter, scale value nine represents that the former indicator is extremely more important than the latter, and scale values ​​two, four, and six are the median values ​​between adjacent scales. An evaluation expert group is formed by professionals in aviation environmental protection, environmental monitoring, meteorological analysis, and airport operation management. The number of experts in the expert group is fixed at ten. All experts independently score the importance of the indicators in the criterion layer and the indicator layer using a unified scale standard.

[0091] Fourth, collect all the judgment matrix data of ten experts, calculate the arithmetic mean of the scale values ​​between each group of indicators, remove abnormal scoring data that deviates too much from the overall mean, retain the valid scoring results and calculate the mean, and construct a comprehensive judgment matrix for the underlying quantitative indicators based on the mean results. The matrix has 3 rows and 3 columns, corresponding to the initial pollution source intensity index, environmental sensitivity index, and diffusion condition index, respectively. The value of each position in the matrix is ​​a fixed scale value after averaging the expert scores.

[0092] Fifth, perform row-by-row numerical calculations on the constructed comprehensive judgment matrix. First, calculate the product of all elements in each row of the matrix, and then take the cube root of the product of each row to obtain the geometric mean of the corresponding indicator. Normalize the geometric mean of the three indicators so that the sum of the three values ​​equals one. The result after processing is the initial weight vector of the three underlying quantitative indicators. The initial weight calculation result of the initial pollution source intensity indicator is fixed at 0.5, the initial weight calculation result of the environmental sensitivity indicator is fixed at 0.3, and the initial weight calculation result of the diffusion condition indicator is fixed at 0.2. All weight values ​​are retained to two decimal places.

[0093] Sixth, calculate the largest eigenvalue of the comprehensive judgment matrix, and use the largest eigenvalue to calculate the consistency index value. The calculation formula adopts the analytic hierarchy process (AHP). The calculated consistency index value is then compared with the average random consistency index value. The average random consistency index value is fixed at 0.58 under the condition of a third-order matrix. The ratio result is the consistency ratio. The test threshold for the consistency ratio is set to 0.1. When the consistency ratio is less than or equal to 0.1, the judgment matrix meets the consistency requirements, and the weight result is valid. When the consistency ratio is greater than 0.1, the process returns to the expert scoring stage to readjust the scaling value until the consistency test is passed. In this embodiment, the consistency ratio calculation result is fixed at 0.06, which meets the test requirements. The AHP is existing technology in this field and is not an inventive solution of this application, so it will not be described in detail here.

[0094] Seventh, the initial weights that pass the consistency test are determined as the final weights of the three underlying quantitative indicators. The initial weight of the pollution source intensity indicator is 0.5, the weight of the environmental sensitivity indicator is 0.3, and the weight of the diffusion condition indicator is 0.2. The target layer, criterion layer, and indicator layer are fully assembled. Combined with the indicator calculation rules that differentiate between aircraft type and flight stage in the early stage, a multi-level evaluation indicator system for aviation pollution sources that differentiates between aircraft type and flight stage is constructed. The corresponding underlying indicator data can be called according to different aircraft types and different flight stages to achieve multi-dimensional differentiated evaluation.

[0095] Eighth, the original scores of the initial pollution source intensity index, environmental sensitivity index, and diffusion condition index are standardized and calibrated. All indicators adopt the same linear transformation rule to uniformly lock the numerical range between zero and one hundred. The standardization process does not change the relative size relationship of the original values ​​of the indicators, but only unifies the units and intervals to ensure that the three indicators have equal participation conditions when weighted.

[0096] Ninth, according to the determined weight allocation scheme, the standardized initial pollution source intensity index value is multiplied by 0.5, the environmental sensitivity index value is multiplied by 0.3, and the diffusion condition index value is multiplied by 0.2. The three product results are added together to obtain the comprehensive evaluation score for each spatial unit grid corresponding to different aircraft types and different flight stages. The comprehensive evaluation score, the three underlying index values, index weights, spatial grid information, aircraft type information, and flight stage information are all integrated to complete the quantitative assignment of all indicators in the multi-level evaluation index system.

[0097] S2: The time-space-hierarchical three-dimensional dynamic weight allocation method is used to assign dynamic weights to each level of the multi-level evaluation index system and perform collaborative correction. At the same time, the QAR data and ADS-B trajectory data in the aircraft operation data are analyzed, and combined with meteorological data and terrain elevation data, the dynamic diffusion process of aviation pollutants is simulated, and the spatial distribution data of pollutant concentration in the target airport area is output.

[0098] The method of employing a three-dimensional dynamic weight allocation approach based on time, space, and hierarchy to assign dynamic weights to each level of indicators in the multi-level evaluation index system and to perform collaborative correction includes:

[0099] B1: Based on the timestamps of the underlying quantitative indicators in the aircraft operation data, the time weight of each indicator data is calculated using a time decay function. The time decay function is existing technology in this field and is not an inventive solution of this application, so it will not be described in detail here.

[0100] Furthermore, the core logic for calculating the time weight of each indicator data using the time decay function is as follows: the closer the data collection time is to the current moment, the stronger the representativeness of the current pollution status, and the higher the time weight; the older the historical data, the lower the time weight. In this way, the uneven impact of different time data on the evaluation results is eliminated, and the time dynamic change characteristics of aviation pollution source emissions are reflected. For example, the time weight of emission data during the morning peak flight period is higher than that during the early morning flight period when there are few flights.

[0101] B2: Based on the spatial cell grid division results, combined with the spatial distribution of terrain elevation and environmentally sensitive point data, the local Moran index calculation method is used to calculate the spatial autocorrelation between each grid and its neighboring grids, and generate spatial weights.

[0102] Furthermore, the specific steps of B2 include:

[0103] First, retrieve the spatial unit grid data of the target airport area, including the unique number, center coordinates, vertex coordinates and boundary range information of each grid. Simultaneously retrieve the pre-processed terrain elevation data and environmental sensitive point data, unify all data to the same geographic coordinate reference system, and ensure that the grid data, terrain elevation data and environmental sensitive point data are completely matched in spatial location to form a basic spatial dataset that can be used for spatial analysis.

[0104] Second, a fixed vehicle-type neighborhood connection rule is used to determine the adjacency relationship between grids. Specifically, the direct adjacency range of each spatial unit grid is defined as the grids that are directly adjacent to it in the four directions of up, down, left, and right. Grids that are diagonally adjacent are not included. The neighborhood search distance for spatial analysis is set to 100 meters. This distance is consistent with the side length of a single spatial unit grid to ensure that only directly adjacent grids are included in the neighborhood calculation range. All grids use a unified neighborhood determination standard.

[0105] Third, taking all spatial cell grids as the analysis object, a two-dimensional spatial adjacency matrix is ​​constructed according to the determined neighborhood determination rules. The rows and columns of the matrix correspond to the unique numbers of the spatial cell grids. For each grid cell, grid cells that are determined to be directly adjacent are assigned a value of one, and grid cells that are determined to be non-adjacent are assigned a value of zero. The matrix construction process adopts fixed numerical rules to ensure that the matrix structure is stable and has spatial consistency, and finally forms an initial spatial adjacency matrix covering all grids.

[0106] Fourth, the initial spatial adjacency matrix is ​​corrected for the first time based on terrain elevation data. The elevation difference between adjacent grids is calculated, and the elevation difference threshold is set to 30 meters. When the elevation difference between adjacent grids exceeds 30 meters, the value in the corresponding adjacency matrix is ​​adjusted from 1 to 0. The second correction is based on environmental sensitive point data. If there are core environmental sensitive points blocking adjacent grids, the value in the corresponding adjacency matrix is ​​adjusted from 1 to 0. Both corrections use fixed parameters. After correction, a final spatial adjacency matrix that integrates terrain and sensitive point features is formed.

[0107] Fifth, the initial pollution source intensity index, environmental sensitivity index and diffusion condition index that have been quantified and assigned are matched to the corresponding spatial unit grids, and the comprehensive evaluation score is selected as the core attribute value for local Moran index calculation. This ensures that each spatial unit grid has a unique corresponding attribute calculation value, and all attribute values ​​are standardized values ​​in the range of 0 to 100 to ensure a unified calculation basis.

[0108] Sixth, a fixed-parameter local Moran index calculation method is adopted. Based on the corrected final spatial adjacency matrix, calculations are performed on each spatial cell grid one by one. During the calculation process, the difference between the attribute value of a single grid and the average value of all grid attributes is calculated first, and then the weighted sum of the differences between the attribute values ​​of adjacent grids and the overall average value is calculated. The two types of differences are correlated and calculated to obtain the local Moran index value of a single grid. All grids use the same calculation process and parameters.

[0109] Seventh, a significance threshold for local Moran's index values ​​is set, and the calculation results are divided into four fixed types: grids with local Moran's index values ​​greater than zero and passing the significance test are judged as high-value clustering or low-value clustering; grids with local Moran's index values ​​less than zero and passing the significance test are judged as high-value and low-value discrete type; and grids with local Moran's index values ​​close to zero are judged as spatially non-autocorrelated type. The confidence level of the significance test is fixed at 95%, and all judgment criteria remain consistent throughout the process.

[0110] Eighth, based on the local Moran index calculation results and spatial clustering type of each spatial unit grid, a corresponding spatial weight is assigned to each pair of adjacent grid relationships. The adjacent weight corresponding to high-value clustered grids is set to 0.9, the adjacent weight corresponding to low-value clustered grids is set to 0.7, the adjacent weight corresponding to spatially discrete grids is set to 0.3, and the adjacent weight corresponding to grids without autocorrelation is set to 0.1. The weight assignment adopts a fixed segmentation standard, and all weight values ​​are within the range of zero to one, forming an initial spatial weight result that matches the spatial autocorrelation characteristics.

[0111] Ninth, the generated initial spatial weight values ​​are standardized by adjusting the sum of the weight values ​​of all adjacent grids in each row to 1, ensuring that the sum of the neighborhood weights of each grid remains consistent, and finally forming a standardized spatial weight matrix.

[0112] B3: Based on the hierarchical structure of the multi-level evaluation index system, the entropy value of the index data at each level is calculated using the entropy weight method to generate hierarchical weights. The smaller the entropy value, the greater the dispersion of the index data, the more information it provides, and the higher the weight; the larger the entropy value, the smaller the dispersion of the index data, the less information it provides, and the lower the weight. The calculation formula of the entropy weight method is existing technology in this field and is not an inventive solution of this application, so it will not be elaborated here. In this embodiment, the hierarchical weights corresponding to the initial pollution source intensity index, environmental sensitivity index, and diffusion condition index are determined to be 0.52, 0.31, and 0.17, respectively.

[0113] B4: Multiply the time weight, spatial weight, and hierarchical weight to obtain the initial three-dimensional dynamic weight. Based on the historical measured concentrations in the pollution monitoring data, construct the objective function using the least squares method. Adjust the initial three-dimensional dynamic weight through an iterative optimization algorithm to minimize the error between the simulated pollutant concentration distribution obtained based on the three-dimensional dynamic weight and the historical measured data, and output the optimized three-dimensional dynamic weight.

[0114] Furthermore, the specific steps of B4 include:

[0115] First, time weights, spatial weights, and hierarchical weights are matched one-to-one according to the same spatial unit grid, the same aircraft type, the same flight stage, and the same monitoring time. The completeness of all matching results is checked to ensure that each combination of space, time, aircraft type, and flight stage has complete weight values ​​of the three types and that all weight values ​​are within the valid range of 0 to 1, forming a standardized weight dataset that can be used for synthetic calculation.

[0116] Second, based on the standardized weight dataset, the time weight, spatial weight and hierarchical weight under the same spatial unit grid, the same aircraft type, the same flight phase and the same monitoring time are multiplied sequentially. No additional adjustment coefficients are added during the calculation process, and the initial three-dimensional dynamic weights of the corresponding grid and the corresponding time period are directly obtained. All spatial unit grids and all time periods adopt completely consistent synthesis rules to ensure that the initial three-dimensional dynamic weight results are consistent and comparable.

[0117] Third, obtain continuous historical measured concentration data of pollutants from environmental monitoring stations within the airport area. The data includes the location of the monitoring station, the monitoring time, and the measured pollutant concentration value. The original sampling interval of the data is set to 1 hour. The historical measured concentration data is converted into a geographic coordinate system consistent with the spatial unit grid. Abnormal values ​​and missing records are removed. The data format is standardized to form a historical measured concentration dataset that can be used for error calculation.

[0118] Fourth, based on the spatial coordinates of the monitoring stations, the historical measured concentration values ​​are mapped to the corresponding spatial unit grids. Each grid covered by a monitoring station is assigned a unique measured concentration value. For spatial unit grids without monitoring stations, spatial interpolation is used to supplement the concentration values. The interpolation search radius is set to 500 meters to ensure that all spatial unit grids have corresponding historical measured concentration values.

[0119] Fifth, the initial three-dimensional dynamic weights are correlated with the values ​​of the multi-level evaluation index system for aviation pollution sources. The calculation is performed grid by grid and time period by time according to the fixed weight superposition rules to generate the initial simulated pollutant concentration distribution results covering the entire area. The correspondence between the weights and the indicators remains unchanged during the calculation process. All grids use the same calculation process to ensure that the simulated concentration distribution results truly reflect the effect of the initial weights.

[0120] Sixth, with the goal of minimizing the error between simulated pollutant concentration values ​​and historical measured concentration values, a least squares objective function is constructed. The difference between the simulated concentration value and the measured concentration value for each spatial cell grid is calculated one by one. All difference results are squared, and the squared results of all grids are summed. The total error value is used as the core output of the objective function. The specific calculation process of the least squares method is the existing technology in this field and is not an inventive solution of this application. It will not be described in detail here.

[0121] Seventh, the gradient descent iterative optimization algorithm is used to adjust the initial three-dimensional dynamic weights. All running parameters of the algorithm are set to fixed values, the maximum number of iterations is set to 1000, the iteration step size is set to 0.001, the error convergence threshold is set to 0.0001, and the upper limit of the weight adjustment range in a single iteration is set to 0.05. All parameters are not dynamically modified during the optimization process. The gradient descent iterative optimization algorithm is the prior art in this field and is not an inventive solution of this application, so it will not be described in detail here.

[0122] Eighth, start the gradient descent iterative optimization algorithm according to the set parameters. In each iteration, the initial three-dimensional dynamic weights are adjusted according to a fixed step size, and the weight adjustment range is controlled within the preset upper limit. After each weight adjustment is completed, the global simulated pollutant concentration distribution is recalculated, and the total error value of the objective function is updated synchronously. The weight value and the corresponding total error value after each iteration are recorded in real time.

[0123] Ninth, continuously monitor the change of total error during the iteration process. When the total error is less than the set error convergence threshold for 50 consecutive iterations, it is determined that the weight optimization result has reached the convergence standard, and the system automatically terminates the iteration calculation. If the number of iterations reaches the set upper limit of 1000 and the convergence condition is still not met, the weight value corresponding to the minimum total error is retained as the optimal result to ensure that the optimization process ends stably.

[0124] Tenth, extract the weight value corresponding to the minimum total error value from all iteration records, formally determine this value as the optimized three-dimensional dynamic weight, check the value range of the optimized three-dimensional dynamic weight, ensure that all values ​​are within the range of zero to one, and maintain spatial and temporal continuity.

[0125] The simulated dynamic diffusion process of aviation pollutants includes:

[0126] C1: Analyze the QAR data and ADS-B trajectory data in the aircraft operation data, calculate the instantaneous emission rate of each pollutant at different times and locations of each aircraft, and generate a mobile point source emission intensity sequence;

[0127] Furthermore, the specific steps of C1 include:

[0128] First, retrieve and preprocess aircraft operation data, including data from the Rapid Access Recorder and Automatic Dependent Surveillance-Broadcast (ADS-B) trajectory, from the airport operation monitoring system. At the same time, download the baseline emission data of the corresponding aircraft engine from the official database of the International Civil Aviation Organization (ICAO). The baseline emission data includes emission reference values ​​for hydrocarbons, carbon monoxide, and nitrogen oxides under standard operating conditions. Unify all data to the same geographic coordinate system and timestamp system to complete the format standardization process.

[0129] Second, the preprocessed fast access recorder data is parsed frame by frame, and the engine fuel flow, engine high-pressure rotor speed, engine low-pressure rotor speed, engine operating mode, and thrust setting value are extracted at fixed one-second time intervals. The extracted parameters are arranged in chronological order to form a continuous and complete sequence of engine operating parameters.

[0130] Third, the preprocessed broadcast automatic dependent surveillance trajectory data is analyzed point by point, and the aircraft's longitude position, latitude position, flight altitude, ground speed, vacuum speed, heading angle and other spatial parameters are extracted according to the one-second time interval consistent with the data of the fast access recorder. The extracted spatial parameters are aligned with the engine operating condition parameter sequence by timestamp to form a set of aircraft operating parameters after spatiotemporal synchronization.

[0131] Fourth, based on the aircraft type information in the aircraft operating parameter set, the baseline emission values ​​for the corresponding aircraft type are retrieved from the ICAO engine emission database as the basis for calculating the instantaneous emission rate. Then, using the engine fuel flow rate in the aircraft operating parameter set, a fixed conversion factor of 0.002 is applied to correct the baseline emission value for fuel flow rate, resulting in the fuel flow rate-corrected emission value. Next, using the engine high-pressure rotor speed and low-pressure rotor speed in the aircraft operating parameter set, a fixed correction factor of 0.0015 is applied to correct the fuel flow rate-corrected emission value for speed, resulting in the speed-corrected emission value. Finally, using the flight altitude and ambient temperature at the location in the aircraft operating parameter set, a fixed environmental correction factor of 0.003 is applied to correct the speed-corrected emission value for altitude and temperature, resulting in the instantaneous emission rates of hydrocarbons, carbon monoxide, and nitrogen oxides for each time and location.

[0132] Fifth, the instantaneous emission rate at each moment is bound with the corresponding timestamp, longitude, latitude, altitude, and speed information to form a complete spatiotemporal emission data item. All emission rate values ​​are uniformly standardized and the unit is unified to grams per second.

[0133] Sixth, taking the complete flight trajectory of the aircraft as the main time and space line, the standardized instantaneous emission rates at each moment are arranged in chronological order to generate a sequence of mobile point source emission intensity covering the entire flight phase, including ground taxiing, takeoff run, initial climb, approach, and landing.

[0134] Seventh, the generated mobile point source emission intensity sequence is checked item by item to confirm that there are no time gaps, no spatial jumps, and no abnormal emission rate values, ensuring the continuity and validity of the sequence. The verified mobile point source emission intensity sequence is then solidified into a standard data structure containing fixed fields such as timestamp, longitude, latitude, altitude, and instantaneous emission rates of various pollutants.

[0135] C2: Input the emission intensity sequence of mobile point source, meteorological data and terrain elevation data into a pre-constructed Lagrange particle diffusion model. In the Lagrange particle diffusion model, each aircraft is simulated as a point source releasing particles representing pollutants at each time step, and the particle movement is driven according to the input meteorological data and terrain elevation data. The Lagrange particle diffusion model is the prior art in this field and is not an inventive solution of this application, so it will not be described in detail here.

[0136] Furthermore, the implementation process of C2 includes:

[0137] First, retrieve the emission intensity sequence of mobile point sources, and simultaneously retrieve standardized meteorological data and topographic elevation data. Unify all data to the same geographic coordinate system and timestamp system, complete the format normalization process, and remove outliers and missing items to form a standardized input dataset that can be directly used for simulation.

[0138] Second, a forward computation Lagrange particle diffusion model based on particle tracking is adopted. The model runs with fixed physical parameters and does not involve machine learning training. The model time step is fixed at one second, the total number of particles released by a single aircraft in a single time step is 100, the effective lifespan of the particles is 300 seconds, the simulation space range is consistent with the target airport evaluation area, the horizontal simulation resolution is 100 meters, and the vertical simulation resolution is 10 meters.

[0139] Third, the model is initialized based on the standardized input dataset to ensure that the particle release rules, weather-driven rules, and terrain constraint rules are fully matched with the input data.

[0140] Fourth, based on the emission intensity sequence of mobile point sources, the current aircraft is simulated as a particle point source releasing pollutants in each time step. The particle release location is consistent with the longitude, latitude, and altitude of the aircraft. The number of particles released is positively correlated with the instantaneous emission rate of pollutants at the current moment. Each particle carries unique identification information including the aircraft number, release time, release location, and pollutant type.

[0141] Fifth, read the meteorological data at the corresponding time and location, extract the wind speed and wind direction parameters, and calculate the horizontal velocity of the particles using a fixed scaling factor of 0.95, with the horizontal velocity of the particles consistent with the wind direction; at the same time, calculate the vertical velocity of the particles based on the atmospheric stability parameters, where the vertical velocity is 0.1 meters per second under unstable stratification conditions, 0.05 meters per second under neutral stratification conditions, and 0.01 meters per second under stable stratification conditions;

[0142] Sixth, read the terrain elevation data corresponding to the current position of the particle, and determine whether there is terrain obstruction or closed area at the particle's position. If the vertical height of the particle is lower than the corresponding terrain elevation, the particle height is forcibly corrected to 10 meters above the terrain elevation. If the particle enters a closed terrain area, the particle's horizontal movement speed is reduced by 50%. All constraint rules use fixed thresholds and are not dynamically adjusted.

[0143] Seventh, based on the horizontal and vertical motion parameters obtained by meteorological driving and the terrain constraints, the particles are driven to perform spatial displacement of the current time step. The particles synchronously complete horizontal and vertical motion in each time step, and the motion trajectory is continuous and conforms to the physical laws of atmospheric diffusion. After each time step of motion is completed, the particle's current longitude, latitude, and altitude information are updated, and it is checked whether the particle's life cycle has expired. Particles that have expired are automatically removed from the simulation system.

[0144] Eighth, continuously track each surviving particle and record its location information from release to extinction. The recorded information includes particle identifier, release time, current time, release location, current location, and type of pollutant. All trajectory information is stored continuously in chronological order.

[0145] Ninth, the particle release, weather-driven, terrain-constrained, particle motion, and trajectory recording operations are performed in a loop with a fixed time step of one second. The loop process fully covers the entire flight phase of the aircraft from ground taxiing to landing, until the simulation of the mobile point source emission intensity sequence is completed and the life cycle of all released particles ends. The simulation is then stopped. After the simulation stops, the complete motion trajectories of all particles are integrated to form a set of pollutant particle motion trajectories covering the target airport area.

[0146] C3: Statistically analyze the motion trajectories of all simulated particles output by the Lagrange particle diffusion model, calculate the pollutant mass contribution of each particle through the spatial cell grid during the simulation period, perform time integration and spatial accumulation of the contributions of all particles in the same grid, and generate spatial distribution data of pollutant concentration covering the target airport area and distinguishing pollutant types.

[0147] Furthermore, the specific steps of C3 include:

[0148] First, retrieve the complete set of particle motion trajectories output by the Lagrange particle diffusion model. The trajectory set includes particle number, release time, disappearance time, continuous motion position, pollutant type and corresponding emission source strength information. Simultaneously retrieve the spatial unit grid data of the target airport area and verify the integrity of the particle motion trajectory data to ensure that there are no missing positions, no time gaps and no abnormal coordinates, thus forming a standardized trajectory dataset.

[0149] Second, based on the boundary coordinates of the spatial unit grid, the position information of the particle at each time step is precisely matched with the corresponding spatial unit grid to determine the unique grid where the particle is located at each moment. All trajectory points of all particles are traversed to complete the global mapping of all particle positions to the spatial unit grid and generate a particle trajectory grid mapping table.

[0150] Third, based on the emission source strength sequence value corresponding to the particle and the particle lifetime duration fixed at 300 seconds, the total mass of pollutants represented by a single particle is calculated. The total mass of pollutants of a single particle is then evenly distributed according to the dwell time of the particle in each spatial cell grid. The dwell time is consistent with the model time step and is fixed at one second. Thus, the instantaneous pollutant mass contribution value of the particle in a single grid is obtained.

[0151] Fourth, traverse the particle trajectory grid mapping table, and for each particle's dwell time record in each spatial cell grid, calculate the particle's contribution to the pollutant mass of the current grid according to a fixed allocation rule. The calculation process is simultaneously bound to the pollutant type, independently distinguishing between hydrocarbons, carbon monoxide, and nitrogen oxides, ensuring that the mass contributions of different types of pollutants are calculated separately and do not interfere with each other.

[0152] Fifth, using spatial cell grids as statistical units, the pollutant mass contribution of all particles falling into the current grid is summarized by type. The contribution values ​​of all particles of the same pollutant type in the same grid are accumulated to obtain the total pollutant mass contribution of all types of pollutants in the grid. The entire spatial cell grid is traversed to complete the global spatial accumulation calculation and generate a grid-level pollutant total mass dataset.

[0153] Sixth, using the total simulation duration as the integration interval, a fixed one-second step accumulation method consistent with the model time step is adopted to perform time integration on the total mass contribution of various pollutants in each spatial cell grid. The pollutant mass contributions at all times during the simulation period are continuously superimposed to obtain the cumulative mass value of various pollutants in each grid during the complete simulation cycle, ensuring the integrity and continuity of the time dimension.

[0154] Seventh, retrieve the fixed physical dimensions of the spatial unit grid, with a horizontal side length of 100 meters and a vertical height of 10 meters, calculate the standard volume of a single grid, divide the cumulative mass value of various pollutants in each grid by the corresponding grid volume to obtain the pollutant concentration value of that grid, with the concentration unit uniformly set to milligrams per cubic meter;

[0155] Eighth, the concentration values ​​of the entire grid are classified and organized according to the pollutant type, and independent concentration datasets for three types of pollutants, namely hydrocarbons, carbon monoxide and nitrogen oxides, are generated. Each dataset includes spatial cell grid number, grid center latitude and longitude, elevation, concentration value and coverage information, so as to clearly distinguish the pollutant type.

[0156] Ninth, conduct full-domain verification of the spatial distribution data of various pollutant concentrations. The verification includes the non-negativity of concentration values, spatial continuity and numerical rationality, removing outlier extreme values, and supplementing missing values ​​in edge grids. The verification standard adopts a fixed threshold, and the concentration values ​​are controlled between 0 mg and 100 mg per cubic meter to ensure that the concentration distribution conforms to the physical laws of pollutant diffusion.

[0157] Tenth, the validated classification concentration dataset will be formatted and solidified to form spatial distribution data of pollutant concentrations covering the entire target airport area and distinguished by pollutant type.

[0158] S3: Divide the target airport area into spatial cell grids, and use an adaptive evaluation algorithm to couple the spatial distribution data of pollutant concentration with dynamic weights to obtain the comprehensive evaluation score of each spatial cell grid;

[0159] The calculation process for the comprehensive evaluation score includes:

[0160] D1: Based on terrain elevation data and environmental sensitive point data, calculate the subdivision demand index of each predetermined initial grid. Using the quadtree partitioning algorithm, with the target airport area as the root node, recursively determine whether the grid subdivision demand index in the current area is higher than the preset first index threshold. If it is higher than the first index threshold, divide the area into four sub-areas until the subdivision demand index of all sub-areas is lower than the first index threshold.

[0161] The calculation of the subdivision demand index for each predetermined initial grid includes:

[0162] D1.1: Based on the terrain elevation data, calculate the elevation standard deviation within the grid. The larger the standard deviation, the greater the terrain undulation. Normalize the standard deviation into the first exponential component. The standard deviation calculation formula is existing technology in this field and is not an inventive solution of this application. It will not be elaborated here.

[0163] D1.2: Based on environmental sensitive point data, calculate the weighted density of various environmental sensitive points in the grid and adjacent grids. The higher the density of sensitive points, the higher the importance of the grid, and normalize it into the second exponential component.

[0164] In this embodiment, environmental sensitive points are divided into three core types: residential area sensitive points, ecological protection area sensitive points, and public facility sensitive points. Each type of sensitive point has a fixed weighting coefficient, with the weighting coefficient for residential area sensitive points set to 0.4, the weighting coefficient for ecological protection area sensitive points set to 0.3, and the weighting coefficient for public facility sensitive points set to 0.3.

[0165] Furthermore, when calculating the weighted density of environmental sensitive points in a single grid, the total spatial area of ​​the current grid and four adjacent grids is retrieved. The area of ​​a single grid is 100×100 meters, and the total spatial area of ​​the five grids is fixed at 5000 square meters. The weighted sum of environmental sensitive points in the grid is divided by the total spatial area of ​​the five grids to obtain the weighted density of environmental sensitive points in the grid. The weighted sum of environmental sensitive points is achieved by accumulating the weighted quantities of various types of sensitive points in the current grid and adjacent grids, and the weighted quantity of each type of sensitive point is the product of the quantity of each type of sensitive point and the corresponding weighting coefficient.

[0166] D1.3: The first index component and the second index component are weighted and summed to obtain the subdivision demand index of the grid.

[0167] Furthermore, a quadtree segmentation algorithm is adopted, taking the entire target airport area as the root node, and recursively determining whether the grid subdivision demand index in the current area is higher than a preset first index threshold. If it is higher than the first index threshold, the area is divided into four sub-areas, until the subdivision demand index of all sub-areas is lower than the first index threshold, including:

[0168] First, retrieve the target airport's entire spatial boundary data, spatial cell grid data, and pollutant concentration spatial distribution data, using a grid with a horizontal resolution of 100 meters. Simultaneously initialize the algorithm parameters, setting the preset first index threshold to 0.7, the maximum grid subdivision depth to 8 levels, and the minimum subdivision cell size to 10 meters.

[0169] Second, the entire spatial area of ​​the target airport is used as the initial root node of the quadtree. This root node completely covers all evaluation areas of the airport and contains all the original spatial unit grids within the area. The spatial distribution data of pollutant concentration is bound to the root node, and the spatial and data information of the root node is initialized, serving as the starting node for recursive segmentation.

[0170] Third, for the spatial region covered by the current quadtree node, traverse the pollutant concentration values ​​of all spatial unit grids within the region, calculate the concentration gradient change value and the grid density weighted value, and obtain the grid subdivision demand index of the current node through a fixed weighted calculation method. The value range of the subdivision demand index is standardized between zero and one.

[0171] Fourth, compare the current node's grid subdivision demand index with the first index threshold, and check whether the current node's subdivision depth has reached the maximum subdivision depth and whether the current node's region size is smaller than the minimum subdivision unit size. Only when the subdivision demand index is higher than the first index threshold, has not reached the maximum subdivision depth, and the region size is greater than or equal to the minimum subdivision unit size, is the current node determined to meet the subdivision conditions.

[0172] Fifth, the spatial region covered by the current node that meets the segmentation conditions is evenly divided into four rectangular sub-regions of equal area. The sub-regions are divided in order of the top left, top right, bottom left, and bottom right directions to form four quadtree child nodes of the same level.

[0173] Sixth, the four newly created sub-nodes are matched with their corresponding spatial ranges, and each sub-node is associated with the spatial unit grid and pollutant concentration spatial distribution data within its range, ensuring that each sub-node has independent and complete calculation data and can directly carry out the next level of processing;

[0174] Seventh, treat each newly generated child node as the new current node, and repeat steps three through six in a loop, recursively executing them level by level according to the depth-first principle;

[0175] Eighth, if the current node's subdivision demand index is lower than or equal to the first index threshold, or reaches the maximum subdivision depth, or the region size is smaller than the minimum subdivision unit size, immediately terminate the recursive segmentation process of the node, mark the node as a quadtree leaf node, and use it as the final adaptive subdivision unit.

[0176] Ninth, continue to execute the recursive segmentation and node marking process until all nodes in the quadtree have completed condition judgment and processing, all nodes are marked as leaf nodes, and there are no intermediate nodes that can be further segmented, thus achieving adaptive quadtree segmentation of the entire target airport area.

[0177] Tenth, summarize the spatial region, grid information and pollutant concentration data corresponding to all leaf nodes, construct a complete quadtree segmentation structure, and generate adaptive subdivided grid data covering the entire target airport area.

[0178] D2: The generated spatial distribution data of pollutant concentration is mapped to each spatial cell grid through bilinear interpolation to obtain the representative value of pollutant concentration for each grid. The global pollutant concentration variation coefficient of the entire target airport area and the local pollutant concentration variation coefficient of each grid local area are calculated. The global pollutant concentration variation coefficient and the local pollutant concentration variation coefficient are compared with the preset global threshold and local threshold, respectively.

[0179] Furthermore, the calculation process for the global pollutant concentration variation coefficient and the local pollutant concentration variation coefficient includes:

[0180] First, retrieve the spatial distribution data of pollutant concentration in the target airport area generated in the early stage, and simultaneously retrieve the spatial unit grid data of the target airport area. The horizontal resolution of the grid is fixed at 100 meters, and all grids have unique numbers and latitude and longitude boundary information. Perform unified coordinate matching and integrity verification on the two types of data.

[0181] Second, for each spatial cell grid to be calculated, four directly adjacent surrounding grids (top left, top right, bottom left, and bottom right) are selected as data sources. The neighborhood search radius is fixed at 100 meters. Based on the concentration values ​​and spatial coordinates of the four neighboring grids, the spatial distribution data of pollutant concentration is mapped to the center of the current grid using a bilinear interpolation method. After traversing all grids in the entire domain, a unique representative value of pollutant concentration for each grid is obtained. The bilinear interpolation method is existing technology in this field and is not an inventive solution of this application, so it will not be described in detail here.

[0182] Third, collect the representative values ​​of pollutant concentrations of all grids in the entire region and sum them up. Divide the sum by the total number of grids to obtain the global average pollutant concentration. Using the global average as a benchmark, iterate through the representative values ​​of concentrations of all grids to calculate the sum of squared differences. Divide the result by the total number of grids and take the square root to obtain the global standard deviation of pollutant concentration.

[0183] Fourth, the standard deviation of global pollutant concentration is divided by the mean of global pollutant concentration to obtain the coefficient of variation of pollutant concentration in the entire target airport area, which is used to reflect the overall spatial dispersion of pollutant concentration in the entire airport area.

[0184] Fifth, taking each spatial unit grid as the center, a fixed local calculation area with a side length of 500 meters is defined. The sum of the representative values ​​of pollutant concentration of all grids in the area is calculated. The sum is divided by the number of grids in the local area to obtain the local average pollutant concentration of the current central grid.

[0185] Sixth, using the average local pollutant concentration as a benchmark, calculate the sum of squared differences of the concentration representative values ​​of all grids in the local area, divide the result by the number of grids in the local area, and take the square root to obtain the standard deviation of the pollutant concentration in the current local area.

[0186] Seventh, divide the standard deviation of the concentration in each local area by the corresponding local mean to obtain the local pollutant concentration variation coefficient of the current central grid. Then, traverse all grids in the entire region to obtain the local pollutant concentration variation coefficient of each grid.

[0187] In this embodiment, the global threshold for the global pollutant concentration variation coefficient is set to 0.5, and the local threshold for the local pollutant concentration variation coefficient is set to 0.8.

[0188] D3: If the global pollutant concentration variation coefficient is greater than the preset global threshold, or the local pollutant concentration variation coefficient is greater than the preset local threshold, then an evaluation function based on extreme value theory is selected, with the extreme value of pollutant concentration within the grid as the evaluation benchmark. Otherwise, an evaluation function based on mean is selected, with the average value of pollutant concentration within the grid as the evaluation benchmark. Then, the representative value of pollutant concentration of each spatial unit grid and the optimized three-dimensional dynamic weight are input into the selected evaluation function to obtain the comprehensive evaluation score of the grid.

[0189] Furthermore, the specific steps of D3 include:

[0190] First, retrieve the global pollutant concentration variation coefficient, the local pollutant concentration variation coefficient, and the preset global and local thresholds. Simultaneously retrieve the representative value of pollutant concentration for each spatial cell grid and the optimized three-dimensional dynamic weight data to ensure that all data corresponds one-to-one with the unique number of the spatial cell grid and complete the data integrity verification.

[0191] Second, two evaluation functions are preset for calculating the comprehensive evaluation score: an evaluation function based on extreme value theory and an evaluation function based on the mean. Both functions adopt a fixed structure and do not involve dynamic adjustment. For the evaluation function based on extreme value theory, the core parameters are set to an extreme value weight coefficient of 0.7 and a concentration correction coefficient of 0.3; for the evaluation function based on the mean, the core parameters are set to a mean weight coefficient of 0.6 and a concentration correction coefficient of 0.4. All parameters remain fixed throughout the calculation process and are not dynamically modified.

[0192] Third, using a single spatial cell grid as the computational unit, all grids are traversed one by one. The previously calculated global pollutant concentration variation coefficient is compared with the preset global threshold, and the local pollutant concentration variation coefficient of the current grid is compared with the preset local threshold. If either condition is met, i.e., the global pollutant concentration variation coefficient is greater than the preset global threshold, or the local pollutant concentration variation coefficient of the current grid is greater than the preset local threshold, the grid is determined to use the evaluation function based on extreme value theory. If neither condition is met, i.e., the global pollutant concentration variation coefficient is less than or equal to the preset global threshold, and the local pollutant concentration variation coefficient of the current grid is less than or equal to the preset local threshold, the grid is determined to use the evaluation function based on the mean.

[0193] Fourth, for grids that are selected based on the extreme value theory evaluation function, the maximum value among all representative pollutant concentrations within the grid is used as the evaluation benchmark. During the screening process, only the maximum value among the representative concentration values ​​is retained, and the minimum value and other intermediate values ​​are not considered to ensure the uniqueness of the evaluation benchmark. For grids that are selected based on the mean evaluation function, the average value of all representative pollutant concentrations within the grid is used as the evaluation benchmark. When calculating the average value, the representative concentration values ​​of all pollutants within the grid are summed and then divided by the total number of representative concentration values ​​to obtain the average concentration value of the grid, which is used as the evaluation benchmark.

[0194] Fifth, for the evaluation function based on extreme value theory, the calculation rule is fixed as follows: multiply the maximum value of the representative pollutant concentration within the grid by the extreme value weight coefficient, and add the representative pollutant concentration of the grid by the concentration correction coefficient. The sum of the two is the preliminary evaluation value of the grid. For the evaluation function based on the mean, the calculation rule is fixed as follows: multiply the average value of the pollutant concentration within the grid by the mean weight coefficient, and add the representative pollutant concentration of the grid by the concentration correction coefficient. The sum of the two is the preliminary evaluation value of the grid. The calculation rules for both functions are fixed throughout the process, without adding any additional adjustment terms.

[0195] Sixth, the optimized three-dimensional dynamic weights are precisely matched with each spatial unit grid to ensure that each grid has a unique corresponding three-dimensional dynamic weight value, and the weight value is between zero and one. The matched three-dimensional dynamic weights, together with the representative value of the pollutant concentration of the grid, are substituted into the previously determined evaluation function. According to the fixed calculation rules of the corresponding evaluation function, the specific calculation operation is performed to obtain the preliminary evaluation score of the grid.

[0196] Seventh, collect the preliminary evaluation scores of all spatial unit grids within the target airport area, and statistically obtain the minimum and maximum values ​​of the preliminary evaluation scores for the entire region. The standardization range is fixed at zero to one hundred points. Subtract the minimum value of the preliminary evaluation score for the entire region from the preliminary evaluation score of a single grid, and then divide by the difference between the maximum and minimum values ​​of the preliminary evaluation scores for the entire region to obtain the standardized evaluation score of that grid.

[0197] Eighth, the standardized evaluation score is retained as an integer, not a decimal, and is used as the final comprehensive evaluation score for that grid to ensure that the scores are standardized and consistent.

[0198] S4: Based on the environmental standards and pollution carrying capacity of different sub-regions within the target airport area, set differentiated pollution judgment thresholds. Using a spatial interpolation algorithm, generate a pollution assessment spatial distribution map based on the comprehensive evaluation score. By comparing the comprehensive evaluation score of each location in the pollution assessment spatial distribution map with the pollution judgment threshold of its respective sub-region, delineate the boundary of the priority supervision area. Then, based on the comprehensive evaluation score, sort the sub-regions within the boundary by internal priority and output the priority supervision area judgment result.

[0199] The pollution determination threshold is set as follows:

[0200] E1: Based on environmental sensitive point data and the functional zoning map of the target airport area, the target airport area is divided into core protection zone, buffer control zone and general impact zone;

[0201] Furthermore, based on the spatial coordinates of the primary sensitive points, a circular protection zone with a radius of one kilometer is delineated as the core protection zone. Based on the outer boundary of the core protection zone, a fixed 500-meter extension is extended outward as the buffer control zone. In conjunction with the protection zone of the secondary sensitive points, all areas within 500 meters of the core protection distance of the secondary sensitive points that are not included in the core protection zone are included in the buffer control zone. The remaining areas within the entire target airport area that are not included in the core protection zone and the buffer control zone are all designated as general impact zones.

[0202] E2: An environmental carrying capacity assessment method is adopted, which combines historical pollution monitoring data and meteorological statistics to calculate the pollution carrying capacity of each sub-region. Among them, the core protection area has a dense population and fragile ecology, and has the lowest pollution carrying capacity. The buffer control area has a moderate population density and a medium pollution carrying capacity. The general impact area is industrial and transportation land, and has the highest pollution carrying capacity. The environmental carrying capacity assessment method is the existing technology in this field and is not an inventive solution of this application, so it will not be described in detail here.

[0203] Furthermore, the pollution carrying capacity calculation process includes:

[0204] First, we retrieved the boundary data of the three sub-regions: the core protection zone, the buffer control zone, and the general impact zone. Simultaneously, we retrieved the historical pollution monitoring data and meteorological statistics data of the target airport area for the past five years. The historical pollution monitoring data included the concentration monitoring values ​​of hydrocarbons, carbon monoxide, and nitrogen oxides at different monitoring points and at different times in each sub-region. The monitoring frequency was fixed at four times a day: 2:00 AM, 8:00 AM, 2:00 PM, and 8:00 PM. The meteorological statistics data included the daily average wind speed, wind direction, atmospheric stability, precipitation, and temperature data for each sub-region. The meteorological data was recorded on the hourly basis for 24 hours a day. All data were verified for integrity, and all types of data were unified to the same geographic coordinate system and timestamp system to form a standardized basic dataset.

[0205] Second, in the environmental carrying capacity assessment method, fixed parameters and fixed calculation rules are used throughout the process, without involving machine learning training. The parameter settings include: pollution monitoring data calibration coefficient of 0.95, meteorological parameter correction coefficient of 0.9, carrying capacity safety threshold of 50 milligrams per cubic meter, and pollutant degradation coefficient of 0.01. All parameters remain fixed throughout the process.

[0206] Third, the historical pollution monitoring data in the standardized basic dataset is calibrated using the pollution monitoring data calibration coefficient to obtain calibrated pollutant concentration values. For each sub-region, the calibrated pollutant concentration data for the past five years is statistically analyzed, and the annual average concentration, daily maximum concentration, and daily minimum concentration values ​​of various pollutants are calculated. The frequency and duration of pollutant concentration exceedances are also statistically analyzed to form a pollution monitoring statistical data set for each sub-region. At the same time, meteorological statistical data are classified by sub-region, and the daily average wind speed, daily average wind direction, daily average atmospheric stability, daily average precipitation, and daily average temperature for each sub-region over the past five years are extracted and statistically analyzed. Atmospheric stability levels are divided according to the proportion of unstable stratification (30%), neutral stratification (40%), and stable stratification (30%). The average values ​​of meteorological parameters corresponding to each stability level are extracted to form a meteorological parameter statistical data set for each sub-region.

[0207] Fourth, the annual average concentration of pollutants in each sub-region is corrected by combining meteorological parameter correction coefficients, and further corrected by atmospheric stability parameters. In the case of unstable stratification, an additional correction coefficient of 0.8 is applied to obtain the corrected baseline concentration of pollutants, so that the pollution data can be accurately matched with meteorological conditions.

[0208] Fifth, calculate the pollution carrying capacity according to fixed rules. First, multiply the annual average wind speed of the sub-region by the pollutant degradation coefficient to obtain the pollutant diffusion and degradation rate. Then, subtract the corrected pollutant baseline concentration value from the carrying capacity safety threshold to obtain the remaining concentration capacity. Finally, multiply the remaining concentration capacity by the area of ​​the sub-region and divide by the pollutant diffusion and degradation rate to obtain the pollution carrying capacity of the sub-region. The unit is uniformly kilograms.

[0209] Sixth, calculations are performed separately for the core protection zone, buffer control zone, and general impact zone. The corrected pollutant baseline concentration, annual average wind speed, and area of ​​each sub-region are substituted into the above fixed calculation rules to obtain the pollutant diffusion and degradation rate, residual concentration capacity, and pollution carrying capacity. The individual carrying capacity of hydrocarbons, carbon monoxide, and nitrogen oxides is calculated for each sub-region, and the individual carrying capacities of the three pollutants are summed to obtain the comprehensive carrying capacity of the corresponding sub-region, so that each sub-region corresponds to a unique comprehensive carrying capacity and individual carrying capacity values ​​for each type of pollutant.

[0210] E3: Based on the ratio of environmental quality standards to pollution carrying capacity of each sub-region, pollution judgment thresholds are set from low to high for the core protection zone, buffer control zone and general impact zone respectively.

[0211] Furthermore, the specific steps of E3 include:

[0212] First, pollution carrying capacity data for three types of sub-regions—core protection zone, buffer control zone, and general impact zone—were retrieved, and the corresponding environmental quality standards for each sub-region were retrieved simultaneously. The environmental quality standard for the core protection zone was fixed at 10 mg / m³, the environmental quality standard for the buffer control zone was fixed at 20 mg / m³, and the environmental quality standard for the general impact zone was fixed at 30 mg / m³. Correlation verification was performed to ensure that the environmental quality standard for each sub-region corresponds one-to-one with its pollution carrying capacity data, thus forming a standardized threshold setting base dataset.

[0213] Second, a rule for calculating the ratio of environmental quality standards to pollution carrying capacity is set. The specific logic is to divide the environmental quality standard value of each sub-region by the comprehensive pollution carrying capacity value corresponding to that sub-region. The calculated ratio is kept to three decimal places. At the same time, a ratio correction coefficient is set and fixed at 1.05 to fine-tune the calculated ratio to ensure that the ratio value is reasonable and meets the control requirements of the sub-region.

[0214] Third, taking the core protection zone, buffer control zone, and general impact zone as calculation units, each sub-region is traversed one by one. The environmental quality standard value of the sub-region is substituted into the ratio calculation rule and divided by the comprehensive pollution carrying capacity value of the sub-region to obtain the preliminary ratio result. Then, the preliminary ratio result is multiplied by the correction coefficient 1.05 to obtain the final environmental quality standard to pollution carrying capacity ratio. Each sub-region corresponds to a unique ratio value. The ratio results of all sub-regions are recorded to form the ratio dataset of each sub-region.

[0215] Fourth, the principle of setting pollution judgment thresholds is strictly followed, which is an increasing rule from low to high. That is, the pollution judgment threshold of the core protection zone is the lowest, the pollution judgment threshold of the buffer control zone is higher than that of the core protection zone, and the pollution judgment threshold of the general impact zone is the highest. At the same time, threshold increment gradient parameters are set. The threshold increment gradient between the core protection zone and the buffer control zone is fixed at 0.2, and the threshold increment gradient between the buffer control zone and the general impact zone is fixed at 0.35. The gradient parameters are fixed throughout to ensure that the threshold increment pattern is consistent and reproducible.

[0216] Fifth, based on the results of all core protection zone ratios, the average of all core protection zone ratios is taken as the benchmark value for the pollution determination threshold of the core protection zone. The benchmark value is then multiplied by a fixed coefficient of 0.9 to obtain the final pollution determination threshold of the core protection zone.

[0217] Sixth, based on the final pollution determination threshold of the core protection zone, and combined with the set incremental gradient parameter, the threshold of the core protection zone (0.15) is added to the incremental gradient (0.2) to obtain the preliminary threshold of the buffer control zone; then, the preliminary threshold of the buffer control zone is multiplied by the set correction coefficient (1.05) to obtain the final pollution determination threshold of the buffer control zone.

[0218] Seventh, based on the final pollution determination threshold of the buffer control zone, and combined with the set incremental gradient parameter, the threshold of the buffer control zone is added to the incremental gradient to obtain the preliminary threshold of the general impact zone; then the preliminary threshold of the general impact zone is multiplied by the set correction coefficient to obtain the final pollution determination threshold of the general impact zone.

[0219] The process of generating the spatial distribution map of the pollution assessment includes:

[0220] G1: Using the comprehensive evaluation score of all spatial unit grids as input variables and terrain elevation as a covariance variable, a variogram model is constructed using the co-kriging spatial interpolation method.

[0221] Furthermore, the specific steps of G1 include:

[0222] First, retrieve the comprehensive evaluation scores of all spatial unit grids in the entire target airport area, and simultaneously retrieve the terrain elevation data of the entire target airport area. After completing the data format normalization, a standardized model input dataset is formed.

[0223] Second, a co-kriging spatial interpolation method with bivariate co-interpolation is adopted, using the comprehensive evaluation score as the interpolation variable and the terrain elevation data as the co-variable to carry out interpolation calculation. The co-kriging spatial interpolation method is the prior art in this field and is not an inventive solution of this application, so it will not be described in detail here.

[0224] Third, the comprehensive evaluation scores and terrain elevation data in the input dataset of the standardized model are scaled separately to unify the two types of variables into a numerical range of 0 to 1, ensuring that the two types of variables are on the same scale, thus forming a standardized variable dataset. The specific processing method is to subtract the minimum value of the comprehensive evaluation score of the entire region from the comprehensive evaluation score of a single grid, and then divide by the difference between the maximum and minimum values ​​of the comprehensive evaluation score of the entire region. The terrain elevation data is standardized using the same rules to ensure that the two types of variables are on the same scale. After processing, a standardized variable dataset is formed.

[0225] Fourth, a spherical model is used to construct the variogram model. The core parameters of the model are fixed and not dynamically adjusted. The sill value is fixed at 0.8 to reflect the total amplitude of the spatial variation of the variables; the range is fixed at 500 meters to define the influence range of the spatial correlation of the variables. The spatial correlation of the comprehensive evaluation scores of two grids with a distance of more than 500 meters can be ignored; the nugget value is fixed at 0.1 to reflect the amplitude of variation caused by random factors; the cross sill value of the covariogram is fixed at 0.6 to reflect the amplitude of covariance between the input variables and covariates. The spherical model is prior art in this field and is not an inventive solution of this application, so it will not be described in detail here.

[0226] Fifth, based on a dataset of fixed parameters and standardized variables, spatial relationship analysis is performed on all spatial cell grids. The spatial distance between any two grids is calculated using the Euclidean distance formula. The spatial covariance and variogram values ​​of the comprehensive evaluation score, the spatial covariance and variogram values ​​of the terrain elevation data, and the cross covariance and cross variogram values ​​between the two types of variables are calculated respectively, forming a complete dataset of covariance and variogram values. The Euclidean distance, covariance, and variogram calculation formulas are existing technologies in this field and are not the inventive solutions of this application, and will not be elaborated here.

[0227] Sixth, based on the spherical model, the least squares method is used to fit the covariance and variogram value dataset. The number of iterations is set to 50, the iteration convergence threshold is 0.001, and the fitting stops when the iteration error meets the convergence requirement. During the fitting process, the core parameters are kept constant to ensure the stability of the model structure.

[0228] Seventh, select 20% of the grid data from the standardized variable dataset as validation samples and the remaining 80% as fitting samples. Substitute the input variables and covariance values ​​of the validation samples into the fitted variogram model, calculate the model's predicted values, and then compare the predicted values ​​with the actual values ​​of the validation samples. Use the mean absolute error to calculate the fitting error, with the allowable error range fixed within 0.05. If the fitting error exceeds this range, do not adjust the core parameters of the model, return to step five, recalculate the spatial covariance and variogram values, and perform model fitting again until the fitting error is within the allowable range, thus obtaining the trained variogram model.

[0229] G2: The nugget value, sill value and range parameter in the variogram model are optimized by cross-validation. The cross-validation method is existing technology in the field and is not an inventive solution of this application, so it will not be described in detail here.

[0230] G3: Using the optimized variogram model, interpolation calculations are performed on the continuous geographic space of the entire target airport area to generate a continuous evaluation score surface, and the evaluation score surface is used as a spatial distribution map of pollution evaluation.

[0231] Furthermore, the specific steps of G3 include:

[0232] First, retrieve the completed variogram model, and simultaneously retrieve the comprehensive evaluation scores, terrain elevation data, and geographic boundary data of all spatial unit grids in the target airport area. Unify all data into the same geographic coordinate system, complete the data integrity and correlation verification, and form a standardized interpolation dataset.

[0233] Second, a bivariate collaborative interpolation algorithm consistent with the structure of the variogram model is adopted, with all parameters set fixed throughout the process. Specifically, the interpolation search radius is fixed at 500 meters, consistent with the range of the variogram model, ensuring spatial relevance adaptation of the interpolation calculation; the interpolation point density is fixed at one interpolation point every 100 meters, and is only deployed in the horizontal continuous geographic space, without setting additional interpolation levels in the vertical direction; the interpolation weight coefficient is fixed at 0.95 to balance the influence weight of the input variables. All parameters are fixed throughout the process without any dynamic adjustment.

[0234] Third, the standardized interpolation dataset is substituted into the bivariate collaborative interpolation algorithm to verify the consistency between the interpolation algorithm parameters and the variogram model parameters, ensuring that the spatial correlation range of the interpolation calculation is consistent with the spatial range of the model fitting, thus preparing for continuous spatial interpolation across the entire domain.

[0235] Fourth, taking the entire geographic boundary of the target airport as the interpolation range, interpolation points are set up on the continuous geographic space of the entire region according to a set density, and each interpolation point corresponds to a unique geographic coordinate. For each interpolation point, the grid comprehensive evaluation score and terrain elevation data within a specified search radius are extracted, substituted into the variogram model and calculated according to the co-kriging interpolation rule to obtain the predicted comprehensive evaluation score of the interpolation point. The predicted score is bound to the corresponding geographic coordinate. After completing the calculation point by point, a dataset of predicted scores and geographic coordinates of the entire region interpolation points is formed.

[0236] Fifth, the predicted scores of all interpolation points across the entire domain are checked one by one to confirm that the scores are within the range of zero to one. At the same time, 15% of the interpolation points are randomly selected, and their predicted scores are compared with the comprehensive evaluation scores of the surrounding grids to confirm that the difference between the two is within 0.05. If there are interpolation points that exceed the error range, the model and algorithm parameters are kept unchanged, the grid data around the interpolation point is retrieved again, and the interpolation is recalculated until the error meets the requirements.

[0237] Sixth, the predicted scores of the verified interpolation points are linearly fitted with their corresponding geographic coordinates using a fixed smoothing coefficient of 0.9, transforming the discrete interpolation points into a continuous surface. This forms an evaluation score surface covering the entire continuous geographic space of the target airport. This surface coincides with the geographic boundary of the entire airport, and each geographic coordinate point corresponds to a unique comprehensive evaluation score, achieving full coverage of the evaluation score in the continuous geographic space. Each score point on the evaluation score surface corresponds to a pollution evaluation marker with the same geographic coordinates on the pollution evaluation spatial distribution map. The scores on the evaluation score surface are divided into intervals, with the division criteria fixed as follows: 0 to 0.3 is a low-pollution area, 0.3-0.7 is a medium-pollution area, and 0.7-1.0 is a high-pollution area.

[0238] The process involves comparing the comprehensive evaluation score of each location on the pollution assessment spatial distribution map with the pollution judgment threshold of its respective sub-region to delineate the boundary of the priority monitoring area, and then ranking the sub-regions within the boundary according to their comprehensive evaluation scores. This includes:

[0239] H1: Perform raster overlay analysis on the spatial distribution map of pollution assessment and the pollution judgment threshold, compare the comprehensive evaluation score with the pollution judgment threshold of the sub-region on a raster-by-raster basis, and screen out the raster units whose comprehensive evaluation score exceeds the pollution judgment threshold, i.e. the dataset of raster units exceeding the standard.

[0240] Furthermore, the specific steps of H1 include:

[0241] First, retrieve the pollution assessment spatial distribution map and pollution judgment threshold data for each sub-region. The pollution judgment threshold for the core protection zone is 0.15, for the buffer control zone it is 0.35, and for the general impact zone it is 0.70. Bind the thresholds of the three types of sub-regions to their corresponding boundaries one by one to form a standardized threshold raster dataset. At the same time, retrieve the raster boundary data of the core protection zone, buffer control zone, and general impact zone to ensure that the boundary data is completely consistent with the raster resolution and geographic coordinate system of the pollution assessment spatial distribution map. Perform integrity and consistency checks on all data and unify the format to form a standardized overlay analysis dataset.

[0242] Second, based on the standardized overlay analysis dataset, a grid-by-grid overlay comparison method is used for grid overlay analysis. The grid matching accuracy is fixed at 100% to ensure that each grid cell in the pollution assessment spatial distribution map can accurately match the boundary grid and threshold data of the corresponding sub-region. The overlay operation type is fixed as an overlay comparison operation, which only realizes the association matching between the grid score and the corresponding threshold, without performing additional calculations. The score comparison accuracy is fixed at three decimal places to ensure that the comparison results are accurate. All parameters are fixed throughout the process and are not dynamically adjusted.

[0243] Third, the grid cells of the pollution assessment spatial distribution map are matched grid by grid with the boundary grids of each sub-region, so that each grid cell uniquely corresponds to one of the sub-regions of the core protection zone, buffer control zone, and general impact zone, without cross-regional or unmatched situations. Then, each grid cell is bound to the pollution judgment threshold of its sub-region, with the grid cells of the core protection zone bound to the threshold of 0.15, the grid cells of the buffer control zone bound to the threshold of 0.35, and the grid cells of the general impact zone bound to the threshold of 0.70, forming a grid-threshold association dataset.

[0244] Fourth, the comprehensive evaluation score of each grid cell is associated and superimposed with the corresponding pollution judgment threshold. While keeping the grid coordinates, scores and threshold values ​​unchanged, a one-to-one correspondence is established to ensure that the scores and thresholds of each grid cell can be extracted and compared synchronously. After the superposition operation is completed, a complete superimposed dataset containing grid coordinates, comprehensive evaluation scores, corresponding sub-regions and pollution judgment thresholds is formed.

[0245] Fifth, compare the comprehensive evaluation score with the pollution judgment threshold of the sub-region grid by grid. The core protection zone grid is compared with 0.15, the buffer control zone grid is compared with 0.35, and the general impact zone grid is compared with 0.70. Record the results of each grid score being greater than or less than or equal to the threshold to form a grid comparison result dataset.

[0246] Sixth, based on the raster comparison result dataset, raster units with comprehensive evaluation scores exceeding the corresponding pollution judgment thresholds are selected. Specifically, raster units with scores greater than 0.15 are selected for the core protection zone, raster units with scores greater than 0.35 are selected for the buffer control zone, and raster units with scores greater than 0.70 are selected for the general impact zone. At the same time, detailed information of each raster unit exceeding the standard is recorded, including raster coordinates, comprehensive evaluation score, corresponding sub-region, corresponding pollution judgment threshold, and the difference between the score and the threshold, forming a dataset of raster units exceeding the standard.

[0247] Seventh, check each of the out-of-standard raster cells in the dataset to confirm that the comprehensive evaluation score of each out-of-standard raster cell does indeed exceed the pollution judgment threshold of the corresponding sub-region. At the same time, randomly select 20% of the out-of-standard raster cells and compare them again with the corresponding threshold to confirm that the comparison result is consistent with the screening result. If a screening error is found, do not adjust the algorithm parameters, but simply return to step five to perform raster-by-raster comparison again until the screening result is accurate.

[0248] H2: Spatial clustering is performed on the selected grid cells using the eight-neighbor connectivity analysis method. Spatially adjacent grid cells are aggregated into candidate regions. Each candidate region is filtered according to a preset minimum area threshold. Candidate regions with an area greater than or equal to the minimum area threshold are retained. The minimum bounding polygon of each retained region is extracted as the boundary of the priority monitoring region.

[0249] Furthermore, the specific steps for H2 include:

[0250] First, retrieve the dataset of out-of-standard raster units, the geographic boundary data of the target airport, and the raster boundary data of each sub-region. Unify all the data into the same geographic coordinate system, perform integrity verification on the dataset of out-of-standard raster units, remove abnormal raster units with missing raster coordinates or missing comprehensive evaluation scores, and form a standardized cluster analysis dataset.

[0251] Second, based on the standardized clustering analysis dataset, an eight-neighborhood connectivity analysis method is adopted. The eight adjacent grid cells in each grid cell in the eight directions of up, down, left, right, upper left, upper right, lower left, and lower right are used as the connectivity determination range. The connection is determined only when all adjacent grid cells are out-of-standard grid cells, and an initial candidate region can be formed only when the number of interconnected out-of-standard grid cells is not less than four. The connectivity analysis method is the prior art in this field and is not an inventive solution of this application, so it will not be described in detail here.

[0252] Third, perform a neighborhood search for each out-of-standard raster cell. For each out-of-standard raster cell, search for its eight neighboring raster cells in eight directions according to the set eight-neighbor search range. Determine whether each neighboring raster cell belongs to the out-of-standard raster cell. If the neighboring raster cell is an out-of-standard raster cell, the two raster cells are determined to be connected and the connection relationship is recorded. If the neighboring raster cell is not an out-of-standard raster cell, it is determined to be disconnected and the connection relationship is not recorded. The neighborhood search and connectivity determination are completed for each raster cell to form a list of connected raster cells for each out-of-standard raster cell.

[0253] Fourth, based on the list of connected raster cells, interconnected out-of-standard raster cells are aggregated into initial candidate regions. That is, all out-of-standard raster cells that are connected through eight-neighborhood relationships are grouped into the same initial candidate region. Individual out-of-standard raster cells that do not reach the minimum number of connected raster cells are removed, and multiple initial candidate regions are finally formed. Each initial candidate region contains multiple interconnected out-of-standard raster cells. At the same time, basic information such as the number of raster cells and the range of raster coordinates of each initial candidate region are recorded to form an initial candidate region dataset.

[0254] Fifth, the area of ​​a single grid cell is determined based on a grid resolution of 100 meters, and a minimum area threshold for candidate regions is set. In this embodiment, it is set to 400 square meters. Only initial candidate regions with an area greater than or equal to the minimum area threshold are retained, and regions with insufficient area are eliminated. The screening process maintains the original grid information and connectivity unchanged. The area is calculated as follows: the area of ​​each initial candidate region is equal to the number of grid cells exceeding the standard contained in the region, multiplied by the area of ​​a single grid cell. The area of ​​a single grid cell is calculated based on the determined grid resolution of 100 meters and is fixed at 10,000 square meters.

[0255] Sixth, calculate the area of ​​each initial candidate region one by one, compare the calculation results with the minimum area threshold, retain candidate regions with an area greater than or equal to 400 square meters, remove candidate regions with an area less than 400 square meters, record the number of grid cells, total area, coordinate range and sub-region information of qualified candidate regions, and form a qualified candidate region dataset.

[0256] Seventh, for each qualified candidate region, fit the smallest bounding polygon based on the coordinates of the outermost grid within the region, so that the polygon completely encloses all grid cells of the corresponding candidate region. Use the smallest bounding polygon as the boundary of the priority monitoring region, and record the vertex coordinates, boundary length and corresponding area of ​​each boundary to form the priority monitoring region boundary dataset.

[0257] H3: For each designated priority regulatory area, extract the comprehensive evaluation score of all spatial unit grids within its boundary, calculate at least one statistical characteristic value among the maximum, average, or median of the comprehensive evaluation score within each area, use the selected statistical characteristic value as the sorting basis, and sort all priority regulatory areas in descending order to generate an internal priority ranking list. The internal priority ranking list and the geographical boundary information of the priority regulatory areas are used together as the output of the priority regulatory area determination result. The calculation formulas for the maximum, average, or median are existing technical content in this field and are not the inventive solution of this application, and will not be elaborated here.

[0258] The embodiments of the present invention have been described above with reference to the accompanying drawings. However, the present invention is not limited to the specific embodiments described above. The specific embodiments described above are merely illustrative and not restrictive. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments under the guidance of the present invention without departing from the spirit and scope of the present invention. All of these variations are within the protection scope of the present invention.

Claims

1. A method for determining priority regulatory areas for aviation pollution sources based on multi-level evaluation, characterized in that, include: Acquire multi-source data for the target airport area, construct a multi-level evaluation index system for aviation pollution sources that distinguishes between aircraft type and flight phase based on the multi-source data, and complete the quantification; the multi-source data includes at least aircraft operation data, pollution monitoring data, environmental sensitive point data, meteorological data, and topographic elevation data; A time-space-hierarchical three-dimensional dynamic weight allocation method is adopted to assign dynamic weights to each level of the multi-level evaluation index system and perform collaborative correction. At the same time, QAR data and ADS-B trajectory data in the aircraft operation data are analyzed, and meteorological data and terrain elevation data are combined to simulate the dynamic diffusion process of aviation pollutants and output the spatial distribution data of pollutant concentration in the target airport area. The target airport area is divided into spatial cell grids. An adaptive evaluation algorithm is used to couple the spatial distribution data of pollutant concentration with dynamic weights to obtain the comprehensive evaluation score of each spatial cell grid. Based on the environmental standards and pollution carrying capacity of different sub-regions within the target airport area, differentiated pollution judgment thresholds are set. Using a spatial interpolation algorithm, a pollution assessment spatial distribution map is generated based on the comprehensive evaluation score. By comparing the comprehensive evaluation score of each location in the pollution assessment spatial distribution map with the pollution judgment threshold of its respective sub-region, the boundary of the priority supervision area is delineated. Based on the comprehensive evaluation score, the sub-regions within the boundary are internally prioritized, and the priority supervision area judgment result is output.

2. The method for determining priority regulatory areas for aviation pollution sources based on multi-level evaluation as described in claim 1, characterized in that, The construction and quantification of a multi-level evaluation index system for aviation pollution sources that distinguishes between aircraft type and flight phase includes: Based on aircraft operation data, the ICAO baseline emissions database is queried to obtain the baseline emission factors of various pollutants per unit time for aircraft engines under standard sea-level static conditions. The aircraft's flight trajectory is then divided into flight phases, and the baseline emission factors are corrected by incorporating temperature and pressure information from meteorological data to obtain the total emissions per flight for each aircraft type in each flight phase, thus constructing an initial pollution source intensity index. The baseline emission factors at least include the emission mass of hydrocarbons, carbon monoxide, and nitrogen oxides per unit time under standard sea-level static conditions. The flight phases include taxiing, takeoff, climb, approach, and landing. Based on environmentally sensitive point data, a buffer analysis of environmentally sensitive points is performed using a geographic information system. According to the type of sensitive point and the distribution density of sensitive points within the buffer zone, an environmental sensitivity index is assigned to each spatial unit grid in the assessment area.

3. The method for determining priority regulatory areas for aviation pollution sources based on multi-level evaluation as described in claim 2, characterized in that, The construction and quantification of a multi-level evaluation index system for aviation pollution sources that distinguishes between aircraft type and flight phase also includes: Based on meteorological and topographic elevation data, a meteorological field and elevation distribution covering the target airport area are generated through spatial interpolation. Based on the influence of wind speed, atmospheric stability and topography on pollutant diffusion, diffusion condition indicators are constructed for each spatial unit grid. The initial pollution source intensity index, environmental sensitivity index, and diffusion condition index are used as the underlying quantitative indexes. The weight of each underlying quantitative index is determined by the analytic hierarchy process (AHP). A multi-level evaluation index system for aviation pollution sources that distinguishes between aircraft type and flight phase is constructed, and the quantitative values ​​of all indicators are assigned.

4. The method for determining priority regulatory areas for aviation pollution sources based on multi-level evaluation as described in claim 1, characterized in that, The method of employing a three-dimensional dynamic weight allocation approach based on time, space, and hierarchy to assign dynamic weights to each level of indicators in the multi-level evaluation index system and to perform collaborative correction includes: Based on the timestamps of the underlying quantitative indicators in the aircraft operation data, the time weight of each indicator data is calculated using a time decay function. Based on the spatial cell grid division results, combined with the spatial distribution of terrain elevation and environmentally sensitive point data, the local Moran index calculation method is used to calculate the spatial autocorrelation between each grid and its neighboring grids, and generate spatial weights. Based on the hierarchical structure of the multi-level evaluation index system, the entropy weight method is used to calculate the entropy value of the index data at each level and generate the hierarchical weights. The time weight, spatial weight, and hierarchical weight are multiplied to synthesize the initial three-dimensional dynamic weight. Based on the historical measured concentrations in the pollution monitoring data, the objective function is constructed using the least squares method. The initial three-dimensional dynamic weight is adjusted through an iterative optimization algorithm to minimize the error between the simulated pollutant concentration distribution obtained based on the three-dimensional dynamic weight and the historical measured data, and the optimized three-dimensional dynamic weight is output.

5. The method for determining priority regulatory areas for aviation pollution sources based on multi-level evaluation as described in claim 1, characterized in that, The simulated dynamic diffusion process of aviation pollutants includes: By analyzing QAR data and ADS-B trajectory data from aircraft operation data, the instantaneous emission rates of each pollutant at different times and locations of each aircraft are calculated, and a sequence of emission intensity of mobile point sources is generated. The emission intensity sequence of mobile point sources, meteorological data, and terrain elevation data are input into a pre-constructed Lagrange particle diffusion model. In the Lagrange particle diffusion model, each aircraft is simulated as a point source releasing particles representing pollutants at each time step, and the particle movement is driven by the input meteorological data and terrain elevation data. The motion trajectories of all simulated particles output by the Lagrange particle diffusion model are statistically analyzed. The pollutant mass contribution of each particle through the spatial cell grid during the simulation period is calculated. The contributions of all particles in the same grid are integrated over time and accumulated spatially to generate spatial distribution data of pollutant concentration covering the target airport area and distinguishing pollutant types.

6. The method for determining priority regulatory areas for aviation pollution sources based on multi-level evaluation as described in claim 1, characterized in that, The calculation process for the comprehensive evaluation score includes: Based on terrain elevation data and environmental sensitive point data, the subdivision demand index of each predetermined initial grid is calculated. The quadtree partitioning algorithm is used, with the target airport area as the root node, to recursively determine whether the grid subdivision demand index in the current area is higher than the preset first index threshold. If it is higher than the first index threshold, the area is divided into four sub-areas until the subdivision demand index of all sub-areas is lower than the first index threshold. The generated spatial distribution data of pollutant concentration is mapped to each spatial cell grid using a bilinear interpolation method to obtain the representative value of pollutant concentration for each grid. The global pollutant concentration variation coefficient of the entire target airport area and the local pollutant concentration variation coefficient of each grid local area are calculated. The global pollutant concentration variation coefficient and the local pollutant concentration variation coefficient are compared with the preset global threshold and local threshold, respectively. If the global pollutant concentration variation coefficient is greater than a preset global threshold, or the local pollutant concentration variation coefficient is greater than a preset local threshold, then an evaluation function based on extreme value theory is selected, with the extreme value of pollutant concentration within the grid as the evaluation benchmark. Otherwise, an evaluation function based on mean is selected, with the average value of pollutant concentration within the grid as the evaluation benchmark. Then, the representative value of pollutant concentration of each spatial unit grid and the optimized three-dimensional dynamic weight are input into the selected evaluation function to obtain the comprehensive evaluation score of the grid.

7. The method for determining priority regulatory areas for aviation pollution sources based on multi-level evaluation as described in claim 6, characterized in that, The calculation of the subdivision demand index for each predetermined initial grid includes: Based on terrain elevation data, the elevation standard deviation within the grid is calculated and normalized to a first exponential component. Based on environmentally sensitive point data, the weighted density of various environmentally sensitive points within the grid and adjacent grids is calculated and normalized to a second exponential component. The first exponential component and the second exponential component are weighted and summed to obtain the subdivision demand index of the grid.

8. The method for determining priority regulatory areas for aviation pollution sources based on multi-level evaluation as described in claim 1, characterized in that, The pollution determination threshold is set as follows: Based on environmentally sensitive point data and the functional zoning map of the target airport area, the target airport area is divided into a core protection zone, a buffer control zone, and a general impact zone. Using the environmental carrying capacity assessment method, combined with historical pollution monitoring data and meteorological statistics, the pollution carrying capacity of each sub-region is calculated. Based on the ratio of environmental quality standards to pollution carrying capacity of each sub-region, pollution judgment thresholds are set from low to high for the core protection zone, buffer control zone, and general impact zone.

9. The method for determining priority regulatory areas for aviation pollution sources based on multi-level evaluation as described in claim 1, characterized in that, The process of generating the spatial distribution map of the pollution assessment includes: Using the comprehensive evaluation score of all spatial unit grids as input variables and terrain elevation as a co-variant, a variogram model is constructed using the co-kriging spatial interpolation method. The nugget value, sill value, and range parameters in the variogram model are optimized using cross-validation. Using the optimized variogram model, interpolation calculations are performed on the continuous geographic space of the entire target airport area to generate a continuously covering evaluation score surface, which is then used as a spatial distribution map of pollution evaluation.

10. The method for determining priority regulatory areas for aviation pollution sources based on multi-level evaluation as described in claim 1, characterized in that, The process involves comparing the comprehensive evaluation score of each location on the pollution assessment spatial distribution map with the pollution judgment threshold of its respective sub-region to delineate the boundary of the priority monitoring area, and then ranking the sub-regions within the boundary according to their comprehensive evaluation scores. This includes: The spatial distribution map of pollution assessment is overlaid with the pollution judgment threshold in a raster overlay analysis. The comprehensive evaluation score is compared with the pollution judgment threshold of the sub-region on a raster-by-raster basis, and raster units whose comprehensive evaluation score exceeds the pollution judgment threshold are selected. The selected grid cells are spatially clustered using the eight-neighbor connectivity analysis method. Spatially adjacent grid cells are aggregated into candidate regions. Each candidate region is filtered according to a preset minimum area threshold, and candidate regions with an area greater than or equal to the minimum area threshold are retained. The minimum bounding polygon of each retained region is extracted as the boundary of the priority monitoring region. For each designated priority regulatory area, the comprehensive evaluation score of all spatial unit grids within its boundary is extracted. At least one statistical characteristic value among the maximum, average, or median comprehensive evaluation scores of each area is calculated. The selected statistical characteristic value is used as the sorting basis, and all priority regulatory areas are sorted in descending order to generate an internal priority ranking list. The internal priority ranking list and the geographical boundary information of the priority regulatory areas are used together as the output of the priority regulatory area determination result.

Images