阅读说明：本技术 一种基于热点网格技术的大气污染物扩散轨迹追踪方法 (Atmospheric pollutant diffusion track tracking method based on hotspot grid technology ) 是由 刘旭东 毛璐 苏琪 李江昊 于 2019-10-23 设计创作，主要内容包括：本发明公开了一种基于热点网格技术的大气污染物扩散轨迹追踪方法,包括以下步骤：开始反演陆地气溶胶,对查算表中不同的气溶胶类型做循环,对查算表中的550nm的AOD做循环,计算改气溶胶类型时,对应550nmAOD和2.25μm观测值条件下的2.25μm的地表反射率,验证地表反射率是否合理,通过确定的0.47μm的地表反射率,观测到的0.47μm的表现反射率是否在这次和上次计算的表现反射率之间,线性插值得到550nm时的AOD,反演得到的550nm的AOD值是否合理,通过外插的方法得到550nm的AOD以及地表反射率的值,选出残差最小时对应的气溶胶类型以及对应的气溶胶光学厚度的值,反演结束。本发明能够实现局部大气污染物成因及扩散分析,最终以产品图的形式生成,可以为环保部门提供执法依据。(The invention discloses an atmospheric pollutant diffusion track tracking method based on a hotspot grid technology, which comprises the following steps of: beginning to invert the terrestrial aerosol, circulating different aerosol types in a lookup table, circulating an AOD (automated optical inspection) value of 550nm in the lookup table, when calculating the aerosol change type, corresponding to the surface reflectivity of 2.25 microns under the observation value conditions of 550nmAOD and 2.25 microns, verifying whether the surface reflectivity is reasonable, determining whether the observed surface reflectivity of 0.47 microns is between the current time and the last calculated surface reflectivity through the determined surface reflectivity of 0.47 microns, linearly interpolating to obtain the AOD value of 550nm, inverting to obtain the AOD value of 550nm and the surface reflectivity value through an extrapolation method, selecting the corresponding aerosol type and the corresponding aerosol optical thickness value when the residual error is minimum, and finishing the inversion. The invention can realize the cause and diffusion analysis of local atmospheric pollutants, and finally generate the results in the form of a product diagram, thereby providing law enforcement basis for environmental protection departments.)

1. An atmospheric pollutant diffusion track tracking method based on a hotspot grid technology comprises the following steps:

step S1: starting to invert the terrestrial aerosol;

step S2: the reflection rates of 0.47 and 0.64/0.25 μm observed by a satellite;

step S3: circulating different aerosol types in the lookup table, if the minimum residue is not found, performing step S4, and if the minimum residue is found, directly performing operation in step S14;

step S4: circulating the AOD of 550nm in the lookup table;

step S5: when the aerosol type is changed, the surface reflectivity of 2.25 mu m under the observation value conditions of 550nmAOD and 2.25 mu m is calculated;

step S6: determining the surface reflectivity of the red and blue channels according to the surface reflectivity of 0.47 and 0.64 mu m under the condition of 'dark background' and the empirical relation of 2.25 mu m;

step S7: verifying whether the surface reflectivity is reasonable (0.47, 0.64 and 2.25 mu m through), if so, directly performing the step S8, and if not, directly jumping to the step S11;

step S8: calculating the atmospheric layer top performance reflectivity at 0.47 mu m through the determined surface reflectivity at 0.47 mu m, the aerosol type under the circulation condition and the AOD at 550 nm;

step S9: whether the observed representative reflectance of 0.47 μm is between the representative reflectances calculated this time and the previous time, and if so, directly proceeds to step S10, and if not, returns to step S4;

step S10: performing linear interpolation to obtain AOD at 550nm and surface reflectivity values of three channels of 0.47, 0.64 and 2.25;

step S11: whether the AOD value of 550nm obtained by inversion is reasonable or not is judged, if not, the step S12 is carried out, and if not, the step S13 is carried out;

step S12: obtaining an AOD of 550nm and a surface reflectivity value by an extrapolation method;

step S13: calculating the residual value of the top expression reflectivity of the atmospheric layer of the channel with the diameter of 0.64 mu m, and returning to the step S3 after the calculation is finished;

step S14: selecting the corresponding aerosol type and the corresponding value of the optical thickness of the aerosol when the residual error is minimum;

step S15: and finishing the inversion.

2. The atmospheric pollutant diffusion trajectory tracking method based on the hotspot grid technology as set forth in claim 1, is characterized in that: FY-4 geostationary satellite data transmitted by the China meteorological office, sunflower-8 geostationary satellite data transmitted by the Japan meteorological hall, American MODIS polar orbit WeChat data and European NPP polar orbit satellite data are selected, and the atmospheric AOD inversion is carried out through satellite data assimilation.

3. The atmospheric pollutant diffusion trajectory tracking method based on the hotspot grid technology as set forth in claim 1, is characterized in that: classifying the aerosol types of the Chinese areas by utilizing the physical parameter information of the aerosols observed by the foundation, and applying the information to a radiation transmission mode to calculate a lookup table containing the equivalent atmospheric transmittance and the equivalent atmospheric apparent emissivity; analyzing the relation between the earth surface reflectivity of the visible light channel in the wave band of 0.47 mu m and 0.66 mu m and the earth surface reflectivity of the short wave infrared wave band, and analyzing the application condition and the influence factor of the correlation relation; developing a cloud detection inversion algorithm to obtain a clear sky area range; carrying out atmospheric absorption correction on the observed quantity of the satellite by using the equivalent quantity of water vapor and ozone in the data of the analysis field; calculating atmospheric molecular Rayleigh scattering contributions by using surface elevation data; after the molecular scattering and atmospheric absorption are deducted, the multichannel satellite observation value under the dark background condition of clear air is combined with a look-up table, and aerosol optical parameters such as the optical thickness of atmospheric aerosol, the mass concentration of a suspended particulate matter column, the Angstrom index and the like are obtained based on an optimization method.

4. The atmospheric pollutant diffusion trajectory tracking method based on the hotspot grid technology as set forth in claim 1, is characterized in that: the optical thickness of the aerosol remotely sensed by the satellite represents the total extinction of the atmosphere in the vertical direction, while PM2.5 represents the concentration of inhalable particles in the ground, so that in order to obtain the concentration of PM2.5 from the AOD remotely sensed by the satellite, the value of the AOD needs to be converted to the ground, the vertical distribution structure of the aerosol needs to be known, the aerosol profile can be determined by adopting a method of an assumed aerosol model, mode simulation (such as GEO-CHEM) and ground observation, and the laser radar can accurately invert the vertical distribution of the aerosol extinction coefficient, so that the method which is the most reliable method for obtaining the aerosol profile at present is provided. After the profile of the aerosol extinction coefficient is inverted by the laser radar, the height correction from the aerosol optical thickness to the ground extinction coefficient mapping can be performed. However, the lidar data can only represent a single-point profile, so that the ground visibility data with more stations are adopted to perform high correction of ground particulate matter inversion in the actual ground particulate matter concentration inversion process. The atmospheric aerosol optical thickness τ can be expressed as the integral of the aerosol extinction coefficient over the entire atmospheric vertical direction:

wherein beta is_zRepresents the extinction coefficient in the vertical direction, m-1. It is generally assumed that the aerosol extinction coefficient is exponentially distributed in the vertical direction:

wherein beta is₀Is the ground extinction coefficient, m-1; h is the aerosol elevation, m. Generally, aerosol density decreases exponentially, and if the spectral distribution of the particles does not change with height, the height at which the extinction coefficient decreases to e-1, the ground extinction coefficient, is referred to as the aerosol level H. Substituting the formula (2) into the formula (1) as follows:

with known site visibility, the horizontal extinction coefficient of the ground can be expressed as:

β₀＝3.912/L (4)

and L represents the horizontal visibility observed by the ground station, and the elevation of the corresponding station can be obtained by combining the optical thickness of the atmospheric aerosol remotely sensed by the satellite of the corresponding station. The elevations of these discrete sites are used to interpolate into high resolution elevation data H (x, y) that match the satellites. And (4) substituting the obtained elevation into the formula (3) to obtain the spatial distribution of the ground extinction coefficient. Since the mass concentration of the dry aerosol is measured by PM2.5, the humidity correction is performed according to the equation (4) after the ground extinction coefficient is obtained.

β_dry＝β₀/f(RH) (4)

Wherein the growth factor

According to the mie scattering theory, the aerosol extinction coefficient is a function of the particle spectral distribution and the extinction characteristics:

wherein Q_ext(r, λ, m) is the extinction efficiency of the particle, and m is the complex refractive index, which is related to the chemical composition of the aerosol and directly affects the scattering absorption characteristics of the particle. x/2 corresponds to the radius of the particle at PM10 or PM 2.5. For an aerosol of a certain particle spectral distribution, the mass concentration of PM can be expressed as:

combining the two formulas (5) and (6) and having the following structure:

from the equation (9), the extinction coefficient of the dry aerosol is directly related to PM2.5, where MEE is mass extinction efficiency, which represents the extinction capability of the aerosol per mass concentration. Beta can thus be established by observing PM2.5 for several typical sites_dryAnd PM2.5, and the obtained correlation is related to beta_dryThe PM2.5 product can be obtained by combining the spatial distribution of the components. According to the principle that the ground visibility is inversely related to the extinction coefficient and the principle that the optical thickness of the aerosol is equal to the integral of the extinction coefficient in the vertical direction, the extinction coefficient grid point data is obtained by calculation, and according to the extinction coefficient grid point data, the ground visibility scattering point observation data and the extinction coefficient grid point data are obtained by calculationAnd calculating the optical coefficient to obtain the horizontal visibility of the ground. The relationship among the optical thickness of the aerosol, the extinction coefficient and the elevation is as follows:

1) calculating the extinction coefficient according to the scattered observation data of the ground visibility, wherein the relation between the ground atmospheric horizontal visibility (meteorological visibility range) R and the atmospheric extinction coefficient is as follows:

wherein epsilon is the contrast threshold of human eyes, and epsilon is 0.02, which is the formula of the visibility of Koschmieder.

2) Calculating the height (elevation) of the boundary layer of the scatter point

And substituting the extinction coefficient and the AOD data of the optical thickness of the aerosol into a formula (3) to obtain the elevation of the discrete point.

3) Performing interpolation operation on the elevation data of the scattered points to obtain elevation grid point data;

and performing Kriging interpolation operation by using the elevation data of the scattered points to obtain lattice point data of the boundary layer height (elevation).

4) Calculating wet extinction coefficient

Substituting lattice point data of boundary layer height (elevation) and AOD data of aerosol optical thickness into formula (10) to obtain wet extinction coefficient

According to the wet extinction coefficient, the formula (11) is replaced, and the visibility can be calculated.

5. The atmospheric pollutant diffusion trajectory tracking method based on the hotspot grid technology as set forth in claim 1, is characterized in that: and adopting an online air quality prediction mode WRF-CHEM to carry out air quality prediction.

6. The atmospheric pollutant diffusion trajectory tracking method based on the hotspot grid technology as set forth in claim 5, wherein: a reasonable parameterization scheme in an air quality prediction mode is selected by utilizing satellite-ground integrated monitoring data, and particularly a boundary layer scheme (PBL) with large influence on a particulate matter prediction result is selected.

7. The atmospheric pollutant diffusion trajectory tracking method based on the hotspot grid technology as set forth in claim 5, wherein: haze numerical prediction is mainly a primary value problem, and due to the lag of an emission source list and uncertainty caused by human factors, the primary value of a haze numerical prediction system is very inaccurate when the haze numerical prediction system is used for solving, the requirement on assimilation is very strong, the mode primary value can be effectively improved through assimilation, and the prediction efficiency is improved. Assimilation application work of satellite data in an atmospheric composition numerical prediction system is carried out earlier abroad, and a certain application result is obtained, but the work in the aspect of atmospheric composition data assimilation in China is carried out very limitedly, and aerosol observation products of geostationary satellites are not assimilated in a mode, so that high-frequency geostationary satellite products are assimilated by a three-dimensional variational assimilation (3DVar) method, namely satellite observation data are continuously utilized for correction in the prediction process, the new initial field is superior to the original mode prediction field (an analysis field is superior to a background field), and the effect of improving the prediction result at the next moment can be achieved.

8. The atmospheric pollutant diffusion trajectory tracking method based on the hotspot grid technology as set forth in claim 5, wherein: since uncertainty exists in both the meteorological field and the pattern scheme, great uncertainty is brought to the final pattern prediction result, and therefore, the uncertainty can be reduced by ensemble prediction means. An estimation of field error distribution is given in advance, at each lattice point in the mode area, an atmospheric state quantity is given a random disturbance quantity based on Gaussian distribution according to a forecast error guess value of a corresponding physical quantity at the position, and a ensemble (ensemble) member is generated by adding the random disturbance. And forecasting by each set member, and finally, providing an optimal air quality forecasting conclusion by using the multi-mode output integration technology for the multi-mode results. Differences between set members can be used to model uncertainty in the initial field; the simulation results of multiple physical schemes can also be used to evaluate uncertainty of the pattern scheme.

9. The atmospheric pollutant diffusion trajectory tracking method based on the hotspot grid technology as set forth in claim 5, wherein: the forward tracing track and backward diffusion track method combines horizontal and vertical movements of the atmosphere, and has complete reference value in the aspects of pollutant transmission, diffusion, sedimentation and the like. The forward tracing trajectory and the backward diffusion trajectory have various statistical methods, such as a potential source area contribution method (PSCF), a concentration weight contribution method (CWT), a retention time weight concentration method (RTWC) and quantitative transport deviation analysis (QTBA). The results produced by the forward tracing trajectory and backward diffusion trajectory method become larger as dead center and time increase. Therefore, the area covered by the track is subjected to hierarchical gridding processing, and the grid standard adopts 2KM by 2 KM. And calculating the PM2.5 mass concentration average value of the time corresponding to each type of track based on the clustering analysis results of the forward tracing track and the backward diffusion track to represent the influence of the air mass in the track source direction on the concentration level of the monitoring area. Because the difference between the diffusion speed of nitrogen oxides and the diffusion speed of particulate matters is more than 3 hours. Whether the pollutants are generated locally is obtained through the diffusion value of the nitrogen oxides and the particulate matters in the same grid. Through data visualization, the particulate matter monitoring results are plotted on the generated PM2.5 product at a frequency of one every 10 minutes. Thereby obtaining a continuous diffusion condition analysis. Thereby realizing the functions of tracing and diffusion.

Technical Field

The invention relates to the technical field of satellite remote sensing, in particular to an atmospheric pollutant diffusion track tracking method based on a hot spot grid technology, which is based on the premise that a product graph is generated by satellite remote sensing hot spot grid, utilizes a forward diffusion track and backward tracing track technology of atmospheric pollutant diffusion conditions, has good functionality particularly on the forecasting of atmospheric pollutant short-time diffusion and the tracing of pollutant tracing and diffusion tracks.

Background

In recent years, along with the increasing severity of the problem of atmospheric pollution, various pollutants are mutually coupled and superposed, the phenomenon of atmospheric pollution appears compressibility, regionality and composite type, the urban and regional atmospheric pollution condition is serious, the daily occurrence frequency of heavy pollution is high, most basic remote sensing products of atmospheric and ecological elements are automatically produced by utilizing an automatic production platform of an environmental remote sensing product (few basic remote sensing products cannot be automatically obtained by utilizing the existing man-machine interaction mode), the automatic production platform of the environmental remote sensing product is butted with data acquired by a distributed satellite data acquisition platform in real time, the data is processed in real time, basic remote sensing products required by subsequent application are obtained, an observation area is divided into grids of 3 kilometers multiplied by 3 kilometers, a grid with relatively high annual average concentration of PM2.5 is comprehensively selected by utilizing a satellite remote sensing technology to serve as an atmospheric pollution hotspot grid, and observing and screening pollutant emission list units, monitoring pollution change in real time by utilizing satellite remote sensing and a mode in a hotspot grid with relatively concentrated gas pollution sources and combining with a ground monitoring micro-station, identifying and selecting an area with abnormal PM2.5 concentration, reducing the range to 2000 m multiplied by 2000 m, calculating and analyzing the conveying and diffusion tracks of atmospheric pollutants, and guiding the development of supervision and law enforcement.

At present, a polluted area is only divided into 2 KM-2 KM grids by satellite remote sensing monitoring, and the pollution concentration of the area is represented by different color codes, but the tracing (local generation or external area diffusion) and diffusion conditions of pollutants in the area cannot be known, so that the functions of tracing the source of the pollutants in the area and forecasting the short-time diffusion can be realized and the tracing of the pollutant diffusion track can be realized by the forward diffusion track and backward diffusion track technology.

Disclosure of Invention

The invention aims to provide an atmospheric pollutant diffusion track tracking method based on a hot spot grid technology, and aims to solve the problems that the monitoring of the satellite remote sensing atmospheric pollutants in the background technology only divides a polluted area into 2 KM-2 KM grids, and different color codes are used for expressing the pollution concentration of the area, but the tracing and diffusion conditions of the pollutants in the area cannot be known, the tracing and diffusion functions of the pollutants in the area can be realized through the forward diffusion track and backward tracing track technology, and the tracing of the pollutant diffusion track is realized.

In order to achieve the purpose, the invention provides the following technical scheme: an atmospheric pollutant diffusion track tracking method based on a hotspot grid technology comprises the following steps:

step S1: starting to invert the terrestrial aerosol;

step S2: the reflection rates of 0.47 and 0.64/0.25 μm observed by a satellite;

step S3: circulating different aerosol types in the lookup table, if the minimum residue is not found, performing step S4, and if the minimum residue is found, directly performing operation in step S14;

step S4: circulating the AOD of 550nm in the lookup table;

step S5: when the aerosol type is changed, the surface reflectivity of 2.25 mu m under the conditions of 550nmAOD and 2.25 mu m observation value is calculated;

step S6: determining the surface reflectivity of the red and blue channels according to the surface reflectivity of 0.47 and 0.64 mu m under the condition of 'dark background' and the empirical relation of 2.25 mu m;

step S7: verifying whether the surface reflectivity is reasonable (0.47, 0.64 and 2.25 mu m through), if so, directly performing the step S8, and if not, directly jumping to the step S11;

step S8: calculating the atmospheric layer top performance reflectivity at 0.47 mu m through the determined surface reflectivity at 0.47 mu m, the aerosol type under the circulation condition and the AOD at 550 nm;

step S9: whether the observed apparent reflectance of 0.47 μm is between the apparent reflectance calculated this time and the apparent reflectance calculated last time, if it is consistent, go directly to step S10, if it is not consistent, go back to step S4;

step S10: performing linear interpolation to obtain AOD at 550nm and surface reflectivity values of three channels of 0.47, 0.64 and 2.25;

step S11: whether the AOD value of 550nm obtained by inversion is reasonable or not is judged, if not, the step S12 is carried out, and if not, the step S13 is carried out;

step S12: obtaining an AOD of 550nm and a surface reflectivity value by an extrapolation method;

step S13: calculating the residual value of the top expression reflectivity of the atmospheric layer of the channel with the diameter of 0.64 mu m, and returning to the step S3 after the calculation is finished;

step S14: selecting the corresponding aerosol type and the corresponding value of the optical thickness of the aerosol when the residual error is minimum;

step S15: and finishing the inversion.

Further, FY-4 geostationary satellite data transmitted by the China meteorological office, sunflower-8 geostationary satellite data transmitted by the Japan meteorological hall, American MODIS polar orbit WeChat data and European NPP polar orbit satellite data are selected, and the atmospheric AOD inversion is carried out through satellite data assimilation.

Further, classifying the aerosol types of the Chinese area by utilizing the aerosol physical parameter information observed by the foundation, and applying the aerosol types to a radiation transmission mode to calculate a lookup table containing the equivalent atmospheric transmittance and the equivalent atmospheric apparent emissivity; analyzing the relation between the earth surface reflectivity of the visible light channel in the wave band of 0.47 mu m and 0.66 mu m and the earth surface reflectivity of the short wave infrared wave band, and analyzing the application condition and the influence factor of the correlation relation; developing a cloud detection inversion algorithm to obtain a clear sky area range; carrying out atmospheric absorption correction on the observed quantity of the satellite by using the equivalent quantity of water vapor and ozone in the data of the analysis field; calculating atmospheric molecular Rayleigh scattering contributions by using surface elevation data; after the molecular scattering and atmospheric absorption are deducted, the multichannel satellite observation value under the dark background condition of clear air is combined with a look-up table, and aerosol optical parameters such as the optical thickness of atmospheric aerosol, the mass concentration of a suspended particulate matter column, the Angstrom index and the like are obtained based on an optimization method.

Further, the optical thickness of the aerosol remotely sensed by the satellite represents the total extinction of the atmosphere in the vertical direction, and PM2.5 represents the concentration of inhalable particles in the ground, so to obtain the concentration of PM2.5 from the AOD remotely sensed by the satellite, the value of the AOD must be converted to the ground, which requires to know the vertical distribution structure of the aerosol, and the aerosol profile can be determined by adopting the methods of assumed aerosol model, mode simulation (such as GEO-CHEM) and ground observation. After the profile of the extinction coefficient of the aerosol is inverted by the laser radar, the height correction from the optical thickness of the aerosol to the mapping of the extinction coefficient of the ground can be carried out. However, the lidar data can only represent a single-point profile, so that the ground visibility data with more stations are adopted to perform high correction of ground particulate matter inversion in the actual ground particulate matter concentration inversion process. The atmospheric aerosol optical thickness τ can be expressed as the integral of the aerosol extinction coefficient over the entire atmospheric vertical direction:

wherein beta is_zRepresents the extinction coefficient in the vertical direction, m-1. It is generally assumed that the aerosol extinction coefficient is exponentially distributed in the vertical direction:

wherein beta is₀Is the ground extinction coefficient, m-1; h is the aerosol elevation, m. Typically, aerosol density decreases exponentially, and if the spectral distribution of the particles does not change with height, the height at which the extinction coefficient decreases to e-1, the ground extinction coefficient, is referred to as the aerosol level H. Substituting the formula (2) into the formula (1) as follows:

with known site visibility, the horizontal extinction coefficient of the ground can be expressed as:

β₀＝3.912/L (4)

and L represents the horizontal visibility observed by the ground station, and the elevation of the corresponding station can be obtained by combining the optical thickness of the atmospheric aerosol remotely sensed by the satellite of the corresponding station. The elevations of these discrete sites are used to interpolate into high resolution elevation data H (x, y) that match the satellites. And substituting the obtained standard height into the formula (3) to obtain the spatial distribution of the ground extinction coefficient. Since the mass concentration of the dry aerosol is measured by PM2.5, the humidity correction is performed according to the equation (4) after the ground extinction coefficient is obtained.

β_dry＝β₀/f(RH) (4)

Wherein the growth factor

Wherein RH0 is 0.4.

According to the mie scattering theory, the aerosol extinction coefficient is a function of the particle spectral distribution and the extinction characteristics:

wherein Q_ext(r, λ, m) is the extinction efficiency of the particle, and m is the complex refractive index, which is related to the chemical composition of the aerosol and directly affects the scattering absorption characteristics of the particle. x/2 corresponds to the radius of the particle at PM10 or PM 2.5. For an aerosol of a certain particle spectral distribution, the mass concentration of PM can be expressed as:

combining the two formulas (5) and (6) and having the following structure:

from the equation (9), the extinction coefficient of the dry aerosol is directly related to PM2.5, where MEE is mass extinction efficiency, which represents the extinction capability of the aerosol per mass concentration. Beta can thus be established by observing PM2.5 for several typical sites_dryAnd PM2.5, and the obtained correlation is related to beta_dryThe PM2.5 product can be obtained by combining the spatial distribution of the components. According to the principle that the ground visibility is inversely related to the extinction coefficient and the principle that the optical thickness of the aerosol is equal to the integral of the extinction coefficient in the vertical direction, the grid point data of the extinction coefficient is obtained through calculation, and the ground horizontal visibility is obtained through calculation according to the extinction coefficient. The relationship among the optical thickness of the aerosol, the extinction coefficient and the elevation is as follows:

1) calculating extinction coefficients according to scattered observation data of ground visibility under the condition of uniform level, wherein the relationship between the ground atmospheric horizontal visibility (meteorological visibility range) R and the atmospheric extinction coefficients is as follows:

wherein epsilon is the contrast sensitivity threshold of human eyes, and epsilon is 0.02, which is the Koschmieder visibility formula.

2) Calculating the height (elevation) of the boundary layer of the scatter point

And substituting the extinction coefficient and the AOD data of the aerosol optical thickness into the formula (3) to obtain the elevation of the discrete point.

3) Performing interpolation operation on the elevation data of the scattered points to obtain elevation grid point data;

and performing Kriging interpolation operation by using the elevation data of the scattered points to obtain grid point data of the height (elevation) of the boundary layer.

4) Calculating wet extinction coefficient

Substituting lattice point data of boundary layer height (elevation) and AOD data of aerosol optical thickness into formula (10) to obtain wet extinction coefficient

According to the wet extinction coefficient, the formula (11) is replaced, and the visibility can be calculated.

Further, an air quality prediction mode WRF-CHEM of online is adopted to carry out air quality prediction.

Furthermore, a reasonable parameterization scheme in an air quality forecasting mode is selected by utilizing satellite-ground integrated monitoring data, and particularly a boundary layer scheme (PBL) which has large influence on particulate matter forecasting results is selected.

Furthermore, the haze numerical value prediction is mainly the problem of initial values, the initial values of the haze numerical value prediction system are very inaccurate when the haze numerical value prediction system is used for solving due to the lag of an emission source list and the uncertainty caused by human factors, the requirement on assimilation is very strong, the mode initial values can be effectively improved through assimilation, and the prediction efficiency is improved. Assimilation application work of satellite data in an atmospheric composition numerical prediction system is carried out earlier abroad, and a certain application result is obtained, but the work in the aspect of atmospheric composition data assimilation in China is carried out very limitedly, and aerosol observation products of geostationary satellites are not assimilated in a mode, so that high-frequency geostationary satellite products are assimilated by a three-dimensional variational assimilation (3DVar) method, namely satellite observation data are continuously utilized for correction in the prediction process, the new initial field is superior to the original mode prediction field (an analysis field is superior to a background field), and the prediction result at the next moment can be improved.

Further, since uncertainty exists in both meteorological field and pattern scheme, great uncertainty is brought to the final pattern prediction result, so that it is considered that uncertainty can be reduced by assembling prediction means. Given and predicted field error distribution estimation in advance, at each lattice point in the mode region, the atmospheric state quantity is given a random disturbance quantity based on Gaussian distribution according to the predicted error guess value of the corresponding physical quantity at the position, and the ensemble (ensemble) members are generated by adding the random disturbance. And forecasting each set member, and finally giving an optimal air quality forecasting conclusion by using the multi-mode output integration technology according to the multi-mode results. Differences between set members can be used to model uncertainty in the initial field; the simulation results of multiple physical schemes can also be used to evaluate uncertainty of the pattern scheme.

Furthermore, the forward tracing track and backward diffusion track method combines horizontal and vertical movements of the atmosphere, and has complete reference value in the aspects of pollutant transmission, diffusion, sedimentation and the like. The forward tracing trajectory and the backward diffusion trajectory have various statistical methods, such as a potential source area contribution method (PSCF), a concentration weight contribution method (CWT), a retention time weight concentration method (RTWC) and quantitative transport deviation analysis (QTBA). The results produced by the forward tracing trajectory and backward diffusion trajectory method become larger as dead center and time increase. Therefore, the area covered by the track is subjected to hierarchical gridding processing, and the grid standard adopts 2KM by 2 KM. And calculating the PM2.5 mass concentration average value of each type of track corresponding time based on the clustering analysis results of the forward tracing track and the backward diffusion track to represent the influence of the air mass in the track source direction on the concentration level of the monitoring area. Since the difference between the diffusion rate of nitrogen oxides and the diffusion rate of particulate matter is more than 3 hours. Whether the pollutants are generated locally is obtained through the diffusion value of the nitrogen oxides and the particulate matters in the same grid. Through data visualization, the particulate matter monitoring results are plotted on the generated PM2.5 product at a frequency of one every 10 minutes. Thereby obtaining a continuous diffusion condition analysis. Thereby realizing the functions of tracing and diffusion.

The invention provides early warning release and emergency response aiming at the heavy pollution condition, and simultaneously applies the source analysis technology to realize the analysis of the regional and urban pollution causes, thereby providing support for making scientific pollution emission reduction measures. And (3) establishing a transportation and diffusion model of the atmospheric pollutants in the air, and forecasting the dynamic distribution and variation trend of the concentration of the atmospheric pollutants in the air. The method can realize local atmospheric pollutant cause and diffusion analysis, and finally generate the result in a product diagram form, thereby providing law enforcement basis for environmental protection departments.

Drawings

FIG. 1 is a geostationary satellite terrestrial aerosol inversion flow chart of an atmospheric pollutant diffusion trajectory tracking method based on a hot spot grid technology;

FIG. 2 is a flow chart of a satellite remote sensing ground particulate matter concentration algorithm of an atmospheric pollutant diffusion track tracking method based on a hot spot grid technology;

fig. 3 is a hotspot grid 2KM by 2KM product diagram of an atmospheric pollutant diffusion trace tracking method based on a hotspot grid technology;

FIG. 4 is a backward diffusion trace diagram of an atmospheric pollutant diffusion trace tracking method based on a hot spot grid technology;

fig. 5 is a forward trace tracing diagram of the atmospheric pollutant diffusion trace tracing method based on the hotspot grid technology.

Detailed Description

The technical solutions in the embodiments of the present invention will be described clearly and completely with reference to the accompanying drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only some embodiments of the present invention, not all embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.

Referring to fig. 1-5, the present invention provides a technical solution: an atmospheric pollutant diffusion track tracking method based on a hotspot grid technology comprises the following steps:

step S1: starting to invert the terrestrial aerosol;

step S2: the reflection rates of 0.47 and 0.64/0.25 μm observed by a satellite;

step S3: circulating different aerosol types in the lookup table, if the minimum residue is not found, performing step S4, and if the minimum residue is found, directly performing operation in step S14;

step S4: circulating the AOD of 550nm in the lookup table;

step S5: when the aerosol type is changed, the surface reflectivity of 2.25 mu m under the conditions of 550nmAOD and 2.25 mu m observation value is calculated;

step S6: determining the surface reflectivity of the red and blue channels according to the surface reflectivity of 0.47 and 0.64 mu m under the condition of 'dark background' and the empirical relation of 2.25 mu m;

step S7: verifying whether the surface reflectivity is reasonable (0.47, 0.64 and 2.25 mu m through), if so, directly performing the step S8, and if not, directly jumping to the step S11;

step S8: calculating the atmospheric layer top performance reflectivity at 0.47 mu m through the determined surface reflectivity at 0.47 mu m, the aerosol type under the circulation condition and the AOD at 550 nm;

step S9: whether the observed apparent reflectance of 0.47 μm is between the apparent reflectance calculated this time and the apparent reflectance calculated last time, if it is consistent, go directly to step S10, if it is not consistent, go back to step S4;

step S10: performing linear interpolation to obtain AOD at 550nm and surface reflectivity values of three channels of 0.47, 0.64 and 2.25;

step S11: whether the AOD value of 550nm obtained by inversion is reasonable or not is judged, if not, the step S12 is carried out, and if not, the step S13 is carried out;

step S12: obtaining an AOD of 550nm and a surface reflectivity value by an extrapolation method;

step S13: calculating the residual value of the top expression reflectivity of the atmospheric layer of the channel with the diameter of 0.64 mu m, and returning to the step S3 after the calculation is finished;

step S14: selecting the corresponding aerosol type and the corresponding value of the optical thickness of the aerosol when the residual error is minimum;

step S15: and finishing the inversion.

In the invention, FY-4 geostationary satellite data transmitted by the China meteorological office, sunflower-8 geostationary satellite data transmitted by the Japan meteorological hall, American MODIS polar orbit WeChat data and European NPP polar orbit satellite data are selected, and the atmospheric AOD inversion is carried out through satellite data assimilation.

In the invention, aerosol types in China are classified by utilizing aerosol physical parameter information observed by a foundation, and the classified aerosol types are applied to a radiation transmission mode to calculate a lookup table containing the equivalent atmospheric transmittance and atmospheric apparent emissivity; analyzing the relation between the earth surface reflectivity of the visible light channel in the wave band of 0.47 mu m and 0.66 mu m and the earth surface reflectivity of the short wave infrared wave band, and analyzing the application condition and the influence factor of the correlation relation; developing a cloud detection inversion algorithm to obtain a clear sky area range; carrying out atmospheric absorption correction on the observed quantity of the satellite by using the equivalent quantity of water vapor and ozone in the data of the analysis field; calculating atmospheric molecular Rayleigh scattering contributions by using surface elevation data; after the molecular scattering and atmospheric absorption are deducted, the multichannel satellite observation value under the dark background condition of clear air is combined with a look-up table, and aerosol optical parameters such as the optical thickness of atmospheric aerosol, the mass concentration of a suspended particulate matter column, the Angstrom index and the like are obtained based on an optimization method.

In the invention, the optical thickness of the aerosol remotely sensed by the satellite represents the total extinction of the atmosphere in the vertical direction, and PM2.5 represents the concentration of inhalable particles on the ground, so that in order to obtain the concentration of PM2.5 from the AOD remotely sensed by the satellite, the value of the AOD needs to be converted to the ground, the vertical distribution structure of the aerosol needs to be known, an assumed aerosol model, a mode simulation (such as GEO-CHEM) and a ground observation method can be adopted for determining the aerosol profile, and the laser radar can accurately invert the vertical distribution of the aerosol extinction coefficient, so that the method is the most reliable method for obtaining the aerosol profile at present. After the profile of the extinction coefficient of the aerosol is inverted by the laser radar, the height correction from the optical thickness of the aerosol to the mapping of the extinction coefficient of the ground can be carried out. However, the lidar data can only represent a single-point profile, so that the ground visibility data with more stations are adopted to perform high correction of ground particulate matter inversion in the actual ground particulate matter concentration inversion process. The atmospheric aerosol optical thickness τ can be expressed as the integral of the aerosol extinction coefficient over the entire atmospheric vertical direction:

wherein beta is_zRepresents the extinction coefficient in the vertical direction, m-1. It is generally assumed that the aerosol extinction coefficient is exponentially distributed in the vertical direction:

wherein beta is₀Is the ground extinction coefficient, m-1; h is the aerosol elevation, m. Typically, aerosol density decreases exponentially, and if the spectral distribution of the particles does not change with height, the height at which the extinction coefficient decreases to e-1, the ground extinction coefficient, is referred to as the aerosol level H. Substituting the formula (2) into the formula (1) as follows:

with known site visibility, the horizontal extinction coefficient of the ground can be expressed as:

β₀＝3.912/L (4)

and L represents the horizontal visibility observed by the ground station, and the elevation of the corresponding station can be obtained by combining the optical thickness of the atmospheric aerosol remotely sensed by the satellite of the corresponding station. The elevations of these discrete sites are used to interpolate into high resolution elevation data H (x, y) that match the satellites. And substituting the obtained standard height into the formula (3) to obtain the spatial distribution of the ground extinction coefficient. Since the mass concentration of the dry aerosol is measured by PM2.5, the humidity correction is performed according to the equation (4) after the ground extinction coefficient is obtained.

β_dry＝β₀/f(RH) (4)

Wherein the growth factor

Wherein RH0 is 0.4.

According to the mie scattering theory, the aerosol extinction coefficient is a function of the particle spectral distribution and the extinction characteristics:

wherein Q_ext(r, λ, m) is the extinction efficiency of the particle, and m is the complex refractive index, which is related to the chemical composition of the aerosol and directly affects the scattering absorption characteristics of the particle. x/2 corresponds to the radius of the particle at PM10 or PM 2.5. For an aerosol of a certain particle spectral distribution, the mass concentration of PM can be expressed as:

combining the two formulas (5) and (6) and having the following structure:

from the equation (9), the extinction coefficient of the dry aerosol is directly related to PM2.5, where MEE is mass extinction efficiency, which represents the extinction capability of the aerosol per mass concentration. Beta can thus be established by observing PM2.5 for several typical sites_dryAnd PM2.5, and the obtained correlation is related to beta_dryThe PM2.5 product can be obtained by combining the spatial distribution of the components. Calculating the optical thickness lattice point data of the existing atmospheric aerosol and the observation data of the ground visibility scattering points according to the principle that the ground visibility is inversely related to the extinction coefficient and the integral of the optical thickness of the aerosol equal to the extinction coefficient in the vertical direction to obtain the optical thickness lattice point data of the existing atmospheric aerosol and the observation data of the ground visibility scattering pointsAnd calculating the extinction coefficient grid data according to the extinction coefficient to obtain the horizontal visibility of the ground. The relationship among the optical thickness of the aerosol, the extinction coefficient and the elevation is as follows:

2) calculating extinction coefficients according to scattered observation data of ground visibility under the condition of uniform level, wherein the relationship between the ground atmospheric horizontal visibility (meteorological visibility range) R and the atmospheric extinction coefficients is as follows:

wherein epsilon is the contrast sensitivity threshold of human eyes, and epsilon is 0.02, which is the Koschmieder visibility formula.

2) Calculating the height (elevation) of the boundary layer of the scatter point

And substituting the extinction coefficient and the AOD data of the aerosol optical thickness into the formula (3) to obtain the elevation of the discrete point.

3) Performing interpolation operation on the elevation data of the scattered points to obtain elevation grid point data;

and performing Kriging interpolation operation by using the elevation data of the scattered points to obtain grid point data of the height (elevation) of the boundary layer.

4) Calculating wet extinction coefficient

Substituting lattice point data of boundary layer height (elevation) and AOD data of aerosol optical thickness into formula (10) to obtain wet extinction coefficient

According to the wet extinction coefficient, the formula (11) is replaced, and the visibility can be calculated.

In the invention, an online air quality prediction mode WRF-CHEM is adopted to carry out air quality prediction.

In the invention, a reasonable parameterization scheme in an air quality prediction mode is selected by utilizing satellite-ground integrated monitoring data, in particular to a boundary layer scheme (PBL) which has larger influence on a particulate matter prediction result.

In the invention, the haze numerical value prediction is mainly the problem of initial values, the initial values of the haze numerical value prediction system are very inaccurate when the haze numerical value prediction system is used for solving due to the lag of an emission source list and the uncertainty caused by human factors, the assimilation requirement is very strong, the mode initial values can be effectively improved through assimilation, and the prediction efficiency is improved. Assimilation application work of satellite data in an atmospheric composition numerical prediction system is carried out earlier abroad, and a certain application result is obtained, but the work in the aspect of atmospheric composition data assimilation in China is carried out very limitedly, and aerosol observation products of geostationary satellites are not assimilated in a mode, so that high-frequency geostationary satellite products are assimilated by a three-dimensional variational assimilation (3DVar) method, namely satellite observation data are continuously utilized for correction in the prediction process, the new initial field is superior to the original mode prediction field (an analysis field is superior to a background field), and the prediction result at the next moment can be improved.

In the invention, because the meteorological field and the mode scheme have uncertainty, great uncertainty is brought to the final mode forecasting result, and the uncertainty can be reduced by integrating forecasting means. Given and predicted field error distribution estimation in advance, at each lattice point in the mode region, the atmospheric state quantity is given a random disturbance quantity based on Gaussian distribution according to the predicted error guess value of the corresponding physical quantity at the position, and the ensemble (ensemble) members are generated by adding the random disturbance. And forecasting each set member, and finally giving an optimal air quality forecasting conclusion by using the multi-mode output integration technology according to the multi-mode results. Differences between set members can be used to model uncertainty in the initial field; the simulation results of multiple physical schemes can also be used to evaluate uncertainty of the pattern scheme.

In the invention, the forward tracing track and backward diffusion track method combines the horizontal and vertical movements of the atmosphere, and has complete reference value in the aspects of pollutant transmission, diffusion, sedimentation and the like. The forward tracing trajectory and the backward diffusion trajectory have various statistical methods, such as a potential source area contribution method (PSCF), a concentration weight contribution method (CWT), a retention time weight concentration method (RTWC) and quantitative transport deviation analysis (QTBA). The results produced by the forward tracing trajectory and backward diffusion trajectory method become larger as dead center and time increase. Therefore, the area covered by the track is subjected to hierarchical gridding processing, and the grid standard adopts 2KM by 2 KM. And calculating the PM2.5 mass concentration average value of each type of track corresponding time based on the clustering analysis results of the forward tracing track and the backward diffusion track to represent the influence of the air mass in the track source direction on the concentration level of the monitoring area. Since the difference between the diffusion rate of nitrogen oxides and the diffusion rate of particulate matter is more than 3 hours. Whether the pollutants are generated locally is obtained through the diffusion value of the nitrogen oxides and the particulate matters in the same grid. Through data visualization, the particulate matter monitoring results are plotted on the generated PM2.5 product at a frequency of one every 10 minutes. Thereby obtaining a continuous diffusion condition analysis and realizing the functions of tracing and diffusion.

And a professional model which is applied to the abnormal atmospheric pollution grid region and used for calculating and analyzing the atmospheric pollutant conveying and diffusing tracks is applied on the basis of the forward diffusion track and the backward tracing track of the hot spot grid. The model has a complete conveying, diffusing and settling mode for processing functions of various meteorological element input fields, various physical processes and different types of pollutant emission sources, and comprises diffusion simulation of sand dust, convection layer ozone, EC, sulfur dioxide, benzene, volcanic eruption, forest fire, mercury and PM2.5, wherein the mode is an Eulerian Lagrange mixed calculation mode, and Lagrange methods are adopted for advection and diffusion calculation. The HYSPLIT track model is generally used for tracking the moving direction of particles or gas carried by airflow, forecasting the wind field situation in real time, analyzing precipitation and researching the track, and is based on a track equation of hydromechanics; the forward trajectory analysis develops the path tracking, which can use trajectory prediction, the idea of simulating and analyzing the air mass path by an airflow three-dimensional trajectory mode is that if the particles in the air drift with the wind, the moving trajectory is the integral of the position vector in time and space, wherein the meteorological field adopts the total-sphere reanalysis data of the American national center for forecasting and the American national center for atmospheric research, 4 times are taken every day, the vertical direction is 17 layers, from the ground to 10hPa, the element field comprises temperature, humidity, pressure, height, horizontal and vertical wind speed fields and the like, the trajectory mode calculates 4 times every day, the height of the starting point is the ground, the range covers the European continent, the trajectory method can qualitatively analyze the direction and the source of the cluster influenced by a specific observation point, understand the distribution and the change of specific aerosol concentration, can give reasonable explanation, can solve the direction influencing the airflow in a larger range, the source area, the moving track and the moving speed of the affected point air mass are obtained, meanwhile, the primary and secondary paths of the air mass moving can be given out through calculation of the proportion of the tracks one by one, the guiding effect is achieved for extrapolating the pollutant development, and the Air Quality Index (AQI) range of the proximity early warning and forecasting within 6 hours is provided by utilizing a machine learning method; providing a short-term air quality prediction result of 6-72 hours by using an air quality prediction mode, wherein product information comprises air quality grade, range, primary pollutants and content, classifying aerosol types in a Chinese area by using aerosol physical parameter information observed by a foundation, and applying the aerosol types to a radiation transmission mode to calculate a lookup table containing the equivalent atmospheric transmittance and the equivalent atmospheric apparent emissivity; analyzing the relation between the earth surface reflectivity of the visible light channel in a wave band of 0.47 mu m and 0.66 mu m and the earth surface reflectivity of the short wave infrared wave band, and analyzing the application conditions and the influence factors of the correlation relation; developing a cloud detection inversion algorithm to obtain a clear sky area range; atmospheric absorption correction is carried out on the observed quantity of the satellite by utilizing the quantities of water vapor, ozone and the like in the analysis field data; calculating atmospheric molecular Rayleigh scattering contributions by using surface elevation data; after molecular scattering and atmospheric absorption are deducted, a multichannel satellite observation value under the condition of clear air and dark background is combined with a lookup table, aerosol optical parameters such as atmospheric aerosol optical thickness, suspended particulate matter column mass concentration and Angstrom index are obtained based on an optimization method, and a product which realizes environmental air quality numerical value prediction operation and display and provides environmental air quality numerical value prediction service is realized by using relevant monitoring products of a foundation and a satellite and other auxiliary data based on meteorological data and a pollutant emission list, application integration, a weather prediction mode (WRF), an emission source processing mode (Smoke), an air quality numerical value prediction mode system (WRF-Chem and the like) and an air quality statistical prediction mode system. And the forecasting precision is continuously improved by the methods of three-dimensional variation and assimilation and ensemble forecasting of observation data such as satellites and the ground. Early warning issuing and emergency response are provided aiming at the heavy pollution condition, and meanwhile, the source analysis technology is applied to realize the cause and analysis of regional and urban pollution, so that support is provided for making scientific pollution emission reduction measures. The method comprises the steps of establishing a transport and diffusion model of the atmospheric pollutants in the air, forecasting the dynamic distribution and change trend of the concentration of the atmospheric pollutants in the air, implementing air quality forecasting of the same city and the periphery, and providing early warning information when meeting and meeting relevant standards. The method is mainly divided into two parts according to forecast timeliness: the method is characterized in that a statistical method-based approach forecast adopts a Convolutional LSTM Network to forecast air quality, factors for distinguishing local pollution sources from future pollution sources are considered to improve forecast accuracy, and monitoring data outside satellite-ground integration is utilized to divide the pollution sources into external pollution and local pollution. For the external pollution sources, dye flow field velocity vectors in a future period are tracked through Taylor freezing hypothesis and key point matching, and a future pollutant source area above the target station is deduced by extrapolating the motion trail of each coordinate point in the future. And then, taking the external pollution source as the input of a convolutional neural network, taking the local pollutant as the input of a convolutional LSTM for the local pollution source, combining the external pollution source, and taking the historical data of the local and the historical data of the external pollutant as a data set input model training network, so that the prediction precision is improved, and the air quality prediction is carried out by adopting an online air quality prediction mode WRF-CHEM based on the short-term prediction of the air quality prediction mode. And the forecast precision is promoted by considering the following three aspects: the selection of the parameterization scheme utilizes satellite-ground integrated monitoring data to select a reasonable parameterization scheme in an air quality prediction mode, particularly a boundary layer scheme (PBL) which has a large influence on a particulate matter prediction result, assimilation of high-frequency secondary satellite observation data and haze numerical value prediction mainly solve the problem of an initial value, the initial value of a haze numerical value prediction system is very inaccurate when the haze numerical value prediction system is used for solving due to delay of a discharge source list and uncertainty caused by human factors, the assimilation requirement is very strong, the mode initial value can be effectively improved through assimilation, and the prediction efficiency is improved. The assimilation application work of satellite data in the atmospheric composition numerical forecasting system is developed earlier abroad, and a certain application result is obtained, the work of assimilating atmospheric composition data in China is very limited, and the aerosol observation product of a geostationary satellite is not assimilated in a mode, therefore, the high-frequency geostationary satellite products are assimilated by a three-dimensional variational assimilation (3DVar) method, namely, the satellite observation data is continuously utilized to correct in the forecasting process, the new initial field is superior to the original mode forecasting field (the analysis field is superior to the background field), further, the forecasting result at the next moment can be improved, and the WRF-CHEM-based ensemble forecasting has uncertainty due to weather fields and mode schemes, therefore, the final pattern prediction result is subjected to a great uncertainty, and therefore, the uncertainty can be reduced by ensemble prediction means. An estimation of field error distribution is given in advance, and at each lattice point in the mode area, an atmospheric state quantity is given a random disturbance quantity based on Gaussian distribution according to a forecast error guess value of a corresponding physical quantity at the position, and a ensemble (ensemble) member is generated by adding the random disturbance. Forecasting each set member, and finally giving an optimal air quality forecasting conclusion according to the mass mode result by using a multi-mode output integration technology, wherein the difference between the set members can be used for simulating the uncertainty of an initial value field; the simulation results of multiple physical schemes can also be used to assess uncertainty of the pattern scheme.

The invention provides early warning release and emergency response aiming at the heavy pollution condition, and simultaneously applies the source analysis technology to realize the analysis of the regional and urban pollution causes, thereby providing support for making scientific pollution emission reduction measures. And (3) establishing a transportation and diffusion model of the atmospheric pollutants in the air, and forecasting the dynamic distribution and variation trend of the concentration of the atmospheric pollutants in the air. The method can realize local atmospheric pollutant cause and diffusion analysis, and finally generate the result in a product diagram form, thereby providing law enforcement basis for environmental protection departments.

It is noted that, in this document, relational terms such as first and second, and the like are used solely to distinguish one entity or action from another entity or action without necessarily requiring or implying any actual such relationship or order between such entities or actions. Also, the terms "comprises," "comprising," or any other variation thereof, are intended to cover a non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements does not include only those elements but may include other elements not expressly listed or inherent to such process, method, article, or apparatus.

Although embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that changes, modifications, substitutions and alterations can be made in these embodiments without departing from the principles and spirit of the invention, the scope of which is defined in the appended claims and their equivalents.