CN112798482B - PM2.5 and PM10 estimation method based on satellite remote sensing - Google Patents
PM2.5 and PM10 estimation method based on satellite remote sensing Download PDFInfo
- Publication number
- CN112798482B CN112798482B CN202011613388.4A CN202011613388A CN112798482B CN 112798482 B CN112798482 B CN 112798482B CN 202011613388 A CN202011613388 A CN 202011613388A CN 112798482 B CN112798482 B CN 112798482B
- Authority
- CN
- China
- Prior art keywords
- central wavelength
- remote sensing
- satellite
- optical thickness
- red
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 38
- 238000012937 correction Methods 0.000 claims abstract description 55
- 230000003287 optical effect Effects 0.000 claims abstract description 51
- 239000005427 atmospheric aerosol Substances 0.000 claims abstract description 43
- 238000011160 research Methods 0.000 claims abstract description 20
- 238000001914 filtration Methods 0.000 claims abstract description 14
- 238000001228 spectrum Methods 0.000 claims abstract description 11
- 239000000443 aerosol Substances 0.000 claims abstract description 10
- 230000005540 biological transmission Effects 0.000 claims abstract description 9
- 238000004364 calculation method Methods 0.000 claims description 34
- 230000005855 radiation Effects 0.000 claims description 12
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 8
- 239000000428 dust Substances 0.000 claims description 7
- 239000000203 mixture Substances 0.000 claims description 7
- 239000004576 sand Substances 0.000 claims description 7
- 238000007781 pre-processing Methods 0.000 claims description 5
- 238000004422 calculation algorithm Methods 0.000 claims description 3
- 239000003897 fog Substances 0.000 claims description 3
- 238000012546 transfer Methods 0.000 claims description 2
- 238000009825 accumulation Methods 0.000 claims 1
- 238000012544 monitoring process Methods 0.000 abstract description 8
- 238000010521 absorption reaction Methods 0.000 abstract description 3
- 230000007613 environmental effect Effects 0.000 abstract description 2
- 239000003344 environmental pollutant Substances 0.000 description 2
- 238000011835 investigation Methods 0.000 description 2
- 239000013618 particulate matter Substances 0.000 description 2
- 231100000719 pollutant Toxicity 0.000 description 2
- 238000002310 reflectometry Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000000295 complement effect Effects 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000003203 everyday effect Effects 0.000 description 1
- 239000010419 fine particle Substances 0.000 description 1
- 230000004907 flux Effects 0.000 description 1
- 229910001385 heavy metal Inorganic materials 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 244000005700 microbiome Species 0.000 description 1
- 238000012821 model calculation Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000002245 particle Substances 0.000 description 1
- 210000002345 respiratory system Anatomy 0.000 description 1
- 238000012216 screening Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 231100000331 toxic Toxicity 0.000 description 1
- 230000002588 toxic effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N15/00—Investigating characteristics of particles; Investigating permeability, pore-volume or surface-area of porous materials
- G01N15/06—Investigating concentration of particle suspensions
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N15/00—Investigating characteristics of particles; Investigating permeability, pore-volume or surface-area of porous materials
- G01N15/06—Investigating concentration of particle suspensions
- G01N15/075—Investigating concentration of particle suspensions by optical means
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Chemical & Material Sciences (AREA)
- Dispersion Chemistry (AREA)
- Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Analytical Chemistry (AREA)
- Biochemistry (AREA)
- General Health & Medical Sciences (AREA)
- General Physics & Mathematics (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
Abstract
The invention relates to the technical field of environmental remote sensing monitoring, and discloses a PM2.5 and PM10 estimation method based on satellite remote sensing, which comprises the following steps of S1: acquiring relative positions of a remote sensing satellite, the sun and a ground surface research site; s2: filtering interference pixels; s3: calculating the optical thickness of the atmospheric aerosol according to the related parameters obtained in the step S1; s4: correcting the optical thickness of the atmospheric aerosol to obtain a correction factor; s5: and (4) estimating the concentrations of PM2.5 and PM10 according to the inversion of the corrected optical thickness of the aerosol. The method establishes a relational equation based on the absorption, transmission and scattering characteristics of the specific waveband spectrum to the concentrations of aerosol, PM2.5 and PM10, so that the PM2.5 and PM10 are estimated, the characteristics of large-scale survey difficulty of atmospheric pollution and high timeliness requirement are solved, and the method has the advantages of simplicity in operation, quickness in estimation and high accuracy.
Description
Technical Field
The invention relates to the technical field of environmental remote sensing monitoring, in particular to a PM2.5 and PM10 estimation method based on satellite remote sensing.
Background
PM2.5 is the air equivalent diameter less than 2.5 m fine particulate matter. The fine particles meeting the condition can float in the atmosphere for a long time, have small particles, large specific surface area and strong activity, are easy to attach toxic and harmful substances such as heavy metal, harmful microorganisms and the like, and have great influence on human bodies and air.
PM10 is particulate matter with equivalent diameter less than 10 μm in air, also called inhalable particulate; after being inhaled, a human body can accumulate in a respiratory system, which is harmful to the health of people and also influences the visibility of the atmosphere, and is one of the key elements influencing the quality of the atmosphere.
At present, two monitoring methods for PM2.5 and PM10 mainly include ground station monitoring and satellite remote sensing monitoring. The ground station monitoring has the characteristics of accurate monitoring result and complete time sequence, but due to the limitation of factors such as manpower and capital investment, the ground monitoring stations are relatively insufficient and are unevenly distributed, the spatial distribution of PM2.5 and PM10 in a large area cannot be reflected, and the macroscopic analysis on the pollutant source and the variation trend is difficult to perform. Remote sensing has the advantages of wide investigation range and strong timeliness, and estimation of PM2.5 and PM10 concentrations through a remote sensing technology is particularly important for timely assessment of atmospheric environment quality, reasonable arrangement of life and production and proposal of a control scheme for solving regional pollution.
Disclosure of Invention
Aiming at the defects in the prior art, the invention aims to provide a PM2.5 and PM10 estimation method based on satellite remote sensing, which has the advantages of simple operation, quick estimation and high accuracy.
In order to achieve the above purpose, the invention provides the following technical scheme:
a PM2.5 and PM10 estimation method based on satellite remote sensing comprises the following steps:
s1: acquiring relative positions of a remote sensing satellite, the sun and a ground surface research site;
s2: filtering interference pixels;
s3: calculating the optical thickness of the atmospheric aerosol according to the related parameters acquired in the step S1;
s4: correcting the optical thickness of the atmospheric aerosol to obtain a correction factor;
s5: and (4) estimating the concentrations of PM2.5 and PM10 according to the inversion of the corrected optical thickness of the aerosol.
In the present invention, further, the step S3 includes calculating the optical thickness of the atmospheric aerosol by calculating the weighting factor of the spectrum, and the calculation formula is:
AOD＝w1·AOD blue +w2·AOD red
wherein AOD is the optical thickness of the atmospheric aerosol, AOD blue For the optical thickness blue component, AOD, of atmospheric aerosols red The weighting coefficients of W1, W2 to the spectrum are the atmospheric aerosol optical thickness red light component.
In the present invention, preferably, the step S4 includes:
altitude correction, namely obtaining the AOD of the optical thickness of the atmospheric aerosol after altitude correction H ；
Height fastener for mixed layerPositive, obtain the mixed layer height correction factor k mix ；
Correcting air humidity to obtain air humidity correction factor k rh 。
In the present invention, preferably, the step S5 includes estimating the concentrations of PM2.5 and PM10 according to the aerosol optical thickness inversion after altitude correction, mixed layer height correction and air humidity correction, where the estimation equation is as follows:
wherein PM is PM2.5 or PM10, AOD H Optical thickness, k, of atmospheric aerosol after correction of altitude rh Correction factor for air humidity, k mix Correction factors k1, k2 and k3 of the mixed layer are fixed coefficients, and values are different when PM2.5 or PM10 is calculated.
In the invention, further, the altitude-corrected optical thickness AOD of the atmospheric aerosol H The calculation formula of (2) is as follows:
wherein f is H And H 0 H is the altitude of the place to be corrected, and AOD is the optical thickness of the atmospheric aerosol.
In the present invention, further, the interference image elements filtered in step S2 at least include cloud, water, ice, snow, sand, dust, and fog, and the filtering algorithm is: and calculating parameters at least comprising snow cover indexes, haze indexes, central wavelength 6.2 mu m thermal infrared channel radiation brightness temperature values and central wavelength 6.2 mu m and central wavelength 6.9 mu m thermal infrared channel radiation brightness temperature difference values, comparing the parameters with related thresholds, and rejecting interference pixels in a filtering mode.
Normalizing the accumulated snow index calculation formula to be:
NDSI＝(Green-SWIR 1.600 )/(Green+SWIR 1.600 )
the NDSI is an asiatic index, the Green is a Green light channel with the central wavelength of 0.510 mu m, and the SWIR1.600 is a short-wave infrared channel with the central wavelength of 2.300 mu m;
the haze index calculation formula is normalized as follows:
NDHI＝(NIR-SWIR 1.600 )/(NIR+SWIR 1.600 )
wherein NDHI is haze index, NIR is infrared channel with central wavelength of 0.860 μm, and SWIR1.600 is short wave infrared channel with central wavelength of 2.300 μm;
the calculation formula of the thermal infrared bright temperature difference between the central wavelength of 6.2 μm and the central wavelength of 6.9 μm is as follows:
T1＝(TIR 6.200 -TIR 6.900 )
wherein, T1 is the thermal infrared bright temperature difference between the central wavelength of 6.2 μm and the central wavelength of 6.9 μm, TIR6.200 is the thermal infrared channel with the central wavelength of 6.200 μm, and TIR6.900 is the thermal infrared channel with the central wavelength of 6.900 μm.
In the present invention, further, the method for acquiring the relative positions of the remote sensing satellite, the sun and the ground surface research site in step S1 includes: according to satellite remote sensing observation time, solar declination, longitude and latitude of a ground surface research site, satellite orbit height and coordinates of satellite down-satellite points, the solar zenith angle, the solar azimuth angle, the satellite zenith angle and the satellite azimuth angle of the research site at the observation moment of the sensor are calculated according to a spherical trigonometric formula, and then the relative position relation of the sun, the satellite and the ground surface research site is calculated.
In the present invention, preferably, the method for calculating the weighting factor of the spectrum is:
w1＝|ρ red-sim -ρ red |/(|ρ bluse-sim -ρ bluse |+|ρ red-sim -ρ red |)
w2＝|ρ bluse-sim -ρ bluse |/(|ρ bluse-sim -ρ bluse |+|ρ red-sim -ρ red |)
wherein ρ bluse-sim Simulated apparent reflectance of blue light obtained for the 6S radiative transfer mode, ρ blue being the apparent reflectance of blue light, ρ blue red-sim Red light simulation acquired for 6S radiation transmission modeApparent reflectance, ρ red Is the red apparent reflectance.
In the present invention, further, step S0 is provided before step S1,
s0: and preprocessing the remote sensing data.
Compared with the prior art, the invention has the beneficial effects that:
the invention discloses a PM2.5 and PM10 concentration estimation method based on satellite remote sensing, which is characterized in that a relational equation is established based on absorption, transmission and scattering characteristics of specific waveband spectrums to aerosol, PM2.5 and PM10 concentrations, specifically, satellite remote sensing data is preprocessed in advance to ensure that data input into a model meet calculation requirements, the relative positions of the sun, a remote sensing satellite and an earth surface research site are calculated to lay a foundation for atmospheric aerosol optical thickness inversion, interference pixels such as cloud, water, ice and snow, sand dust, fog and the like are filtered, influence on PM2.5 and PM10 during estimation is prevented, and estimation accuracy is improved.
The invention carries out altitude correction, mixed layer correction and humidity correction to improve the accuracy of data, corrects correction factors of obtaining different influence estimation values and finally realizes the estimation of PM2.5 and PM10 concentration. According to the method, the concentration estimation of PM2.5 and PM10 is realized by using satellite remote sensing data, the characteristics of large difficulty in large-scale atmospheric pollution survey and high timeliness requirement are solved, and the method has the advantages of simplicity in operation, rapidness in estimation and high accuracy.
Drawings
The accompanying drawings, which are included to provide a further understanding of the invention and are incorporated in and constitute a part of this specification, illustrate embodiments of the invention and together with the description serve to explain the principles of the invention and not to limit the invention. In the drawings:
FIG. 1 is a flow chart of a method for estimating PM2.5 and PM10 based on satellite remote sensing according to the present invention;
FIG. 2 is a calculation block diagram of a PM2.5 and PM10 estimation method based on satellite remote sensing according to the invention;
FIG. 3 is a graph of the output result of PM2.5 obtained by the PM2.5 and PM10 estimation method based on satellite remote sensing according to the present invention;
fig. 4 is a graph of an output result of PM10 obtained by a method for estimating PM2.5 and PM10 based on satellite remote sensing according to the present invention.
Detailed Description
The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the 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.
It will be understood that when an element is referred to as being "secured to" another element, it can be directly on the other element or intervening elements may also be present. When a component is referred to as being "connected" to another component, it can be directly connected to the other component or intervening components may also be present. When a component is referred to as being "disposed on" another component, it can be directly on the other component or intervening components may also be present. The terms "vertical," "horizontal," "left," "right," and the like as used herein are for illustrative purposes only.
Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs. The terminology used in the description of the invention herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. As used herein, the term "and/or" includes any and all combinations of one or more of the associated listed items.
Referring to fig. 1 to 2, a preferred embodiment of the present invention provides a method for estimating PM2.5 and PM10 based on satellite remote sensing, comprising the steps of:
s0: preprocessing the remote sensing data;
s1: acquiring relative positions of a remote sensing satellite, the sun and a ground surface research site;
s2: filtering out interference pixels;
s3: calculating the optical thickness of the atmospheric aerosol according to the related parameters acquired in the step S1;
s4: correcting the optical thickness of the atmospheric aerosol to obtain a correction factor;
s5: and (4) estimating the concentration of PM2.5 and PM10 according to the inversion of the optical thickness of the aerosol after correction.
The invention relates to a PM2.5 and PM10 concentration estimation method based on satellite remote sensing, which is characterized in that a relational equation is established based on absorption, transmission and scattering characteristics of specific waveband spectrums on aerosol, PM2.5 and PM10 concentrations, the specific estimation method comprises the steps of satellite remote sensing data preprocessing, calculation of relative positions of a sun, a remote sensing satellite and an earth surface research site, calculation of atmospheric aerosol optical thickness related correction coefficients, estimation of PM2.5 and PM10 concentrations and the like, and finally estimation of PM2.5 and PM10 concentrations is achieved. According to the method, the concentration estimation of PM2.5 and PM10 is realized by using satellite remote sensing data, the problems of large atmospheric pollution survey difficulty and high timeliness requirement in the prior art are solved, and the method has the advantages of simplicity in operation, quickness in estimation and high accuracy.
In the present invention, further, step S0 is to preprocess the remote sensing data, that is, to perform radiometric calibration, geometric correction and radiometric correction, image registration and atmospheric correction on the obtained original remote sensing data, specifically, according to the gain, bias and conversion parameters provided by the selected remote sensing satellite, convert the gray value recorded by the sensor into the corresponding physical quantity (reflectivity, radiant flux density or brightness temperature), further complete the geometric correction, and resample the selected data according to the bilinear method to ensure that the data of the input model meets the requirements, thereby laying the foundation for model calculation.
In the present invention, step S1 further obtains the relative positions of the remote sensing satellite, the sun, and the earth surface research site. The relative position calculation is to calculate the sun zenith angle, the sun azimuth angle, the satellite zenith angle and the satellite azimuth angle of the researched site at the observation moment of the sensor by utilizing a spherical trigonometric formula according to the satellite remote sensing observation time, the sun declination, the longitude and latitude of the ground surface research site, the satellite orbit height and the satellite subsatellite point coordinate, and acquire the azimuth angles of the sun and the satellite.
Specifically, the solar zenith angle, the solar azimuth angle, the satellite zenith angle and the satellite azimuth angle are obtained, and the calculation method for obtaining the azimuth angles of the sun and the satellite is as follows:
the solar altitude θ is:
wherein, the first and the second end of the pipe are connected with each other,
the latitude of the ground research site is shown, delta is solar declination, omega s The solar hour angle.
The solar declination delta is calculated by the formula:
δ＝＝0.006918-0.399912·cos(ω se )+0.070257·sin(ω se )-0.006758·cos(2·ω se )+0.000907·sin(2·ω se )-0.002697·cos(3·ω se )+0.00148·sin(3·ω se )
wherein, ω is se An orbital foot for the earth to move around the sun, defined as follows:
ω se ＝2×π×(J-1)/365
wherein J is the sequence, 1 is the sequence of 1 month and 1 day every year, 1 is added to the sequence of every day, and the rest are analogized in turn.
The sun's zenith leg sza is complementary to the sun's altitude theta, i.e.
sza-pi/2-theta solar time angle omega s The calculation formula is as follows:
ω s ＝15×(ST-12)×π/180
wherein, when ST is true sun, 24 hours timing method, the year, month, day, hour, minute and second of satellite remote sensing observation time are respectively year, month, day, hour, minute and second world, when the satellite remote sensing observation time is converted to the place, the longitude is calculated according to the longitude of the satellite remote sensing observation time, and the longitude is increased to east every 15 degrees for 1 hour (or decreased by 1 hour).
Calculating the solar azimuth angle omega according to the parameters solar The concrete formula is as follows:
ω solar the formula for calculating ar csin (cos (δ) × sin (ω)/sin (sza)) + pi satellite zenith angle zva is:
wherein phi sat Is the longitude of the sub-satellite point and phi is the longitude of the underground study site.
Satellite observation azimuth angle omega sat The calculation formula is as follows:
finally, the relative azimuth raa of the sun and the satellite is obtained based on the variables, and the specific calculation formula is as follows:
raa＝|ω solar -ω sat |
therefore, the relative position relation between the sun, the satellite and the ground surface research site is obtained through the parameters, and a foundation is laid for the inversion of the optical thickness of the atmospheric aerosol.
In the present invention, further, the interference pixels filtered in step S2 at least include cloud, water, ice and snow, sand and dust, and fog, and the filtering algorithm is: and calculating parameters at least comprising snow cover indexes, haze indexes, central wavelength 6.2 mu m thermal infrared channel radiation brightness temperature values and central wavelength 6.2 mu m and central wavelength 6.9 mu m thermal infrared channel radiation brightness temperature difference values, comparing the parameters with related thresholds, and rejecting interference pixels in a filtering mode.
Normalizing the accumulated snow index calculation formula to be:
NDSI＝(Green-SWIR 1.600 )/(Green+SWIR 1.600 )
the NDSI is an asiatic index, the Green is a Green light channel with the central wavelength of 0.510 mu m, and the SWIR1.600 is a short-wave infrared channel with the central wavelength of 2.300 mu m;
normalizing the haze index calculation formula as follows:
NDHI＝(NIR-SWIR 1.600 )/(NIR+SWIR 1.600 )
wherein NDHI is haze index, NIR is a near infrared channel with a central wavelength of 0.860 mu m, and SWIR1.600 is a short wave infrared channel with a central wavelength of 2.300 mu m;
the calculation formula of the thermal infrared bright temperature difference between the central wavelength of 6.2 μm and the central wavelength of 6.9 μm is as follows:
T1＝(TIR 6.200 -TIR 6.900 )
wherein, T1 is the thermal infrared bright temperature difference between the central wavelength of 6.2 μm and the central wavelength of 6.9 μm, TIR6.200 is the thermal infrared channel with the central wavelength of 6.200 μm, and TIR6.900 is the thermal infrared channel with the central wavelength of 6.900 μm.
Specifically, the filtering logic of cloud, water, ice, snow, sand, fog and the like is that the value of a wave band with the central wavelength of 0.47 microns is more than 0.01, or the value of the central wavelength of 0.47 microns is more than 0.35, and the value of a thermal infrared channel with the central wavelength of 6.9 microns is more than 285, or the normalized snow index is more than 0.13, or the temperature difference between the thermal infrared light with the central wavelength of 6.2 microns and the thermal infrared light with the central wavelength of 6.9 microns is less than 0, and the pixel point with the haze index of more than 0.02 is filtered, so that the phenomenon that the characteristic is mixed into PM2.5 and PM10 during estimation to cause larger errors of data results is prevented. Therefore, the influence of the characteristics on the estimated values of PM2.5 and PM10 is effectively prevented by filtering clouds, water bodies, ice and snow, sand dust, fog and the like, and the estimation accuracy is improved.
In the present invention, further, step S3 is to calculate the optical thickness of the atmospheric aerosol, and to perform TIR on the channel with the central waveband of 2.3 μm 2.300 And (4) screening (dark pixels), and selecting the grid points with the numerical value less than 0.2. The optical thickness calculation is carried out by adopting a 6S radiation transmission mode, specifically, the optical thickness of the atmospheric aerosol is calculated by calculating the weight coefficient of the spectrum, and the calculation formula is as follows:
AOD＝w1·AOD blue +w2·AOD red
wherein AOD is the optical thickness of the atmospheric aerosol, AOD blue For the optical thickness blue component, AOD, of atmospheric aerosols red Is the red light component of the optical thickness of the atmospheric aerosol,w1, W2 weight coefficient for spectra.
In the present invention, preferably, the calculation of the relevant parameters is performed by the 6S radiation transmission mode of the weight coefficient of the spectrum. The specific calculation method comprises the following steps:
w1＝|ρ red-sim -ρ red |/(|ρ bluse-sim -ρ bluse |+|ρ red-sim -ρ red |)
w2＝|ρ bluse-sim -ρ bluse |/(|ρ bluse-sim -ρ bluse |+|ρ red-sim -ρ red |)
where ρ is bluse-sim Blue light acquired for the 6S radiation transmission mode simulates an apparent reflectance, ρ blue being the blue light apparent reflectance, ρ red-sim Simulated apparent reflectance, ρ, of red light acquired for 6S radiation transmission mode red Is the red apparent reflectance.
The simulated apparent reflectivity is mainly calculated through a wave band with a central wave band of 2.3 mu m, parameters such as a solar zenith angle, a satellite zenith angle, a relative azimuth angle of the sun and the satellite and the like are required to be introduced during calculation, and logic operation is carried out through a method of searching related parameters and calling a determine function.
Therefore, the optical thickness of the atmospheric aerosol is calculated, the pollutant concentration distribution can be reflected through the data, and the obtained optical thickness of the atmospheric aerosol is corrected to enable the data to reflect the concentrations of PM2.5 and PM10 more truly in order to further improve the estimation accuracy of PM2.5 and PM10 due to different factors influencing the optical thickness of the atmospheric aerosol in different regions.
In the present invention, preferably, the step S4 of correcting the optical thickness of the atmospheric aerosol, and obtaining the correction factor includes:
altitude correction, namely obtaining the AOD of the optical thickness of the atmospheric aerosol after altitude correction H ；
The height of the mixed layer is corrected to obtain a height correction factor k of the mixed layer mix ；
Correcting air humidity to obtain air humidity correction factor k rh 。
Specifically, aiming at the influence factors of the optical thickness of the atmospheric aerosol, three factors of larger altitude, mixed layer height and air humidity are preferably selected, and the three factors can be adjusted according to different regions in different seasons and different influence factors, and in a specific embodiment provided by the invention, the Tianjin region is selected as a test point, and a specific calculation method is as follows:
calculating an atmospheric aerosol optical thickness altitude correction factor, setting the altitude of a ground research site as H for correction, and adopting a calculation formula as follows:
wherein AOD H Is the optical thickness f of the atmospheric aerosol after altitude correction H And H 0 To correct the coefficient, take f specifically H ＝0.09874，H 0 0.11765, H is the altitude of the site to be corrected and AOD is the atmospheric aerosol optical thickness.
The height correction factor of the mixed layer is calculated by the height of the convection layer, and the calculation formula is as follows:
k mix ＝1/η
where η is the mixed layer height (km).
The air humidity correction factor is a function related to the relative humidity of air, and the calculation formula is as follows:
k rh ＝α+β·k+γ·k 2
wherein k is rh The correction factor for the air humidity is defined, alpha, beta and gamma are coefficients related to seasons, k is a function related to the air humidity, and the calculation formula is as follows:
k＝(1-rh) -1
wherein rh is the relative humidity of air.
In the specific embodiment of the present invention, the humidity rh is 0.2, k is 1.25 calculated according to the above formula, and k is calculated rh ＝α+β·k+γ·k 2 And alpha, beta and gamma are coefficients dependent on the season
In the present invention, preferably, step S5 includes estimating the concentrations of PM2.5 and PM10 according to the aerosol optical thickness inversion after altitude correction, mixed layer height correction and air humidity correction, where the estimation equation is:
wherein PM is PM2.5 or PM10, AOD H Optical thickness, k, of atmospheric aerosol after correction of altitude rh Correction factor for air humidity, k mix Correction factor, k, for the mixing layer 1 、k 2 、k 3 For a fixed coefficient, the values are different when calculating PM2.5 or PM 10.
Wherein when PM is PM2.5, k is 1 ＝-298、k 2 -0.59172、k 3 The PM2.5 concentration was calculated 320 and the result output is shown in fig. 3.
When PM is PM10, k is 1 ＝-527、k 2 ＝-0.60606、k 3 The PM10 concentration was calculated 603 and the result output is shown in fig. 4.
In the present embodiment, it is preferred that,
the working principle is as follows: the method comprises the steps of preprocessing satellite remote sensing data to ensure that data input into a model meet calculation requirements, laying a foundation for atmospheric aerosol optical thickness inversion by calculating the relative positions of the sun, the remote sensing satellite and a ground surface research site, filtering interference pixels such as cloud, water, ice and snow, sand dust, fog and the like, preventing influences on PM2.5 and PM10 during estimation, and improving estimation accuracy.
The invention carries out altitude correction, mixed layer correction and humidity correction to improve the accuracy of data, corrects correction factors of obtaining different influence estimation values and finally realizes the estimation of PM2.5 and PM10 concentration. According to the method, the concentration estimation of PM2.5 and PM10 is realized by using satellite remote sensing data, the characteristics of large atmospheric pollution range investigation difficulty and high timeliness requirement are solved, and the method has the advantages of simplicity in operation, quickness in estimation and high accuracy.
The above description is intended to describe in detail the preferred embodiments of the present invention, but the embodiments are not intended to limit the scope of the claims of the present invention, and all equivalent changes and modifications made within the technical spirit of the present invention should fall within the scope of the claims of the present invention.
Claims (8)
1. A PM2.5 and PM10 estimation method based on satellite remote sensing is characterized by comprising the following steps:
s1: obtaining relative positions of a remote sensing satellite, the sun and a ground surface research site;
s2: filtering out interference pixels;
s3: calculating the optical thickness of the atmospheric aerosol according to the related parameters acquired in the step S1;
s4: correcting the optical thickness of the atmospheric aerosol to obtain a correction factor;
s5: estimating the concentrations of PM2.5 and PM10 according to the inversion of the corrected optical thickness of the aerosol;
the step S3 includes calculating the optical thickness of the atmospheric aerosol by calculating the weighting factor of the spectrum, where the calculation formula is:
AOD＝w1·AOD blue +w2·AOD red
wherein AOD is the optical thickness of the atmospheric aerosol, AOD blue For the optical thickness blue component, AOD, of atmospheric aerosols red For the red light component of the optical thickness of the atmospheric aerosol, the weighting coefficients of W1 and W2 to the spectrum, the calculation formula of W1 and W2 is:
w1＝|ρ red-sim -ρ red |/(|ρ bluse-sim -ρ bluse |+|ρ red-sim -ρ red |)
w2＝|ρ bluse-sim -ρ bluse |/(|ρ bluse-sim -ρ bluse |+|ρ red-sim -ρ red |)
wherein ρ bluse-sim Simulated apparent reflectance, ρ, of blue light acquired for 6S radiation transmission mode blue Apparent reflectance for blue light, ρ red-sim Simulated apparent reflectance, ρ, of red light acquired for 6S radiative transfer mode red Is the red apparent reflectance.
2. The method for estimating PM2.5 and PM10 based on satellite remote sensing according to claim 1, wherein the step S4 includes:
altitude correction, namely obtaining the AOD of the optical thickness of the atmospheric aerosol after altitude correction H ；
The height of the mixed layer is corrected to obtain a height correction factor k of the mixed layer mix ；
Correcting air humidity to obtain air humidity correction factor k rh 。
3. The method for estimating PM2.5 and PM10 based on satellite remote sensing according to claim 2, wherein the step S5 comprises estimating the concentrations of PM2.5 and PM10 according to the aerosol optical thickness inversion after altitude correction, mixed layer height correction and air humidity correction, and the estimation equation is as follows:
wherein PM is PM2.5 or PM10, AOD H Optical thickness, k, of atmospheric aerosol after correction of altitude rh Correction factor for air humidity, k mix For the correction factor of the mixed layer, k1, k2 and k3 are fixed coefficients, and values are different when PM2.5 or PM10 is calculated.
4. The method for estimating PM2.5 and PM10 based on satellite remote sensing as claimed in claim 3, wherein the altitude-corrected AOD is an optical thickness of the atmospheric aerosol H The calculation formula of (2) is as follows:
wherein, f H And H 0 H is the altitude of the place to be corrected, and AOD is the optical thickness of the atmospheric aerosol.
5. The method for estimating PM2.5 and PM10 based on satellite remote sensing according to claim 1, wherein the interference pixels filtered in step S2 at least include cloud, water, ice and snow, sand and dust, and fog, and the filtering algorithm is as follows: and calculating parameters at least comprising snow cover indexes, haze indexes, central wavelength 6.2 mu m thermal infrared channel radiation brightness temperature values and central wavelength 6.2 mu m and central wavelength 6.9 mu m thermal infrared channel radiation brightness temperature difference values, comparing the parameters with related thresholds, and rejecting interference pixels in a filtering mode.
6. The method for estimating PM2.5 and PM10 based on satellite remote sensing according to claim 5, wherein the calculation formula for normalizing the snow cover index is as follows:
NDSI＝(Green-SWI R 1.600 )/(Creen+SWI R 1.600 )
the NDSI is an asiatic index, the Green is a Green light channel with the central wavelength of 0.510 mu m, and the SWIR1.600 is a short-wave infrared channel with the central wavelength of 2.300 mu m;
normalizing the haze index calculation formula as follows:
NDHI＝(NI R-SWI R 1.600 )/(NI R+SWI R 1.600 )
wherein NDHI is haze index, NIR is an infrared channel with the central wavelength of 0.860 mu m, and SWIR1.600 is a short-wave infrared channel with the central wavelength of 2.300 mu m;
the calculation formula of the thermal infrared bright temperature difference between the central wavelength of 6.2 μm and the central wavelength of 6.9 μm is as follows:
T1＝(TIR 6.200 -TIR 6.900 )
wherein, T1 is the thermal infrared bright temperature difference between the central wavelength of 6.2 μm and the central wavelength of 6.9 μm, TIR6.200 is the thermal infrared channel with the central wavelength of 6.200 μm, TIR6.900 is the thermal infrared channel with the central wavelength of 6.900 μm;
the filtering logic of cloud, water, ice, snow, dust and fog is that the value of a wave band with a central wavelength of 0.47 mu m is more than 0.01, or the value of the central wavelength of 0.47 mu m is more than 0.35 and the value of a thermal infrared channel with the central wavelength of 6.9 mu m is more than 285, or the normalized snow accumulation index is more than 0.13, or the thermal infrared bright temperature difference between the central wavelength of 6.2 mu m and the central wavelength of 6.9 mu m is less than 0, and the pixel point with the haze index more than 0.02 is filtered.
7. The method for estimating PM2.5 and PM10 based on satellite remote sensing according to claim 1, wherein the method for acquiring the relative positions of the remote sensing satellite, the sun and the ground surface research site in the step S1 comprises: according to satellite remote sensing observation time, solar declination, longitude and latitude of a ground surface research site, satellite orbit height and satellite sub-satellite point coordinates, the solar zenith angle, the solar azimuth angle, the satellite zenith angle and the satellite azimuth angle of the research site at the observation moment of the sensor are calculated according to a spherical trigonometric formula, and then the relative position relation of the sun, the satellite and the ground surface research site is calculated.
8. The method for estimating PM2.5 and PM10 based on satellite remote sensing according to claim 1, wherein step S0 is further provided before step S1,
s0: and (5) preprocessing remote sensing data.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011613388.4A CN112798482B (en) | 2020-12-30 | 2020-12-30 | PM2.5 and PM10 estimation method based on satellite remote sensing |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011613388.4A CN112798482B (en) | 2020-12-30 | 2020-12-30 | PM2.5 and PM10 estimation method based on satellite remote sensing |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112798482A CN112798482A (en) | 2021-05-14 |
| CN112798482B true CN112798482B (en) | 2022-09-06 |
Family
ID=75805858
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011613388.4A Active CN112798482B (en) | 2020-12-30 | 2020-12-30 | PM2.5 and PM10 estimation method based on satellite remote sensing |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112798482B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116363047B (en) * | 2022-08-23 | 2024-01-30 | 生态环境部卫星环境应用中心 | Early warning method for pollution risk of straw incineration atmosphere |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102103204A (en) * | 2011-01-26 | 2011-06-22 | 环境保护部卫星环境应用中心 | Inversion method for land aerosols optical thickness based on environment satellite 1 |
| CN105678085A (en) * | 2016-01-12 | 2016-06-15 | 环境保护部卫星环境应用中心 | PM2.5 concentration estimation method and system |
| CN110411918A (en) * | 2019-08-02 | 2019-11-05 | 中国科学院遥感与数字地球研究所 | A kind of PM2.5 concentration remote-sensing evaluation method based on satellite polarization technology |
| CN110726653A (en) * | 2019-09-25 | 2020-01-24 | 中国电子科技集团公司第二十七研究所 | PM based on heaven and earth integration information2.5Concentration monitoring method |
| CN111487171A (en) * | 2020-04-28 | 2020-08-04 | 武汉拓宝科技股份有限公司 | Forward and backward combined dual-wavelength dispersed fire smoke detector and method |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7379592B2 (en) * | 2003-01-21 | 2008-05-27 | United States Of America As Represented By The Secretary Of The Navy | System and method for significant dust detection and enhancement of dust images over land and ocean |
| CN103018736B (en) * | 2012-12-03 | 2014-11-26 | 北京航空航天大学 | Satellite-borne remote sensor radiation calibration method based on atmospheric parameter remote sensing retrieval |
| CN103389494B (en) * | 2013-07-24 | 2015-04-22 | 中国科学院南海海洋研究所 | Novel atmospheric correction method of water color remote sensing data of case II water body |
-
2020
- 2020-12-30 CN CN202011613388.4A patent/CN112798482B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102103204A (en) * | 2011-01-26 | 2011-06-22 | 环境保护部卫星环境应用中心 | Inversion method for land aerosols optical thickness based on environment satellite 1 |
| CN105678085A (en) * | 2016-01-12 | 2016-06-15 | 环境保护部卫星环境应用中心 | PM2.5 concentration estimation method and system |
| CN110411918A (en) * | 2019-08-02 | 2019-11-05 | 中国科学院遥感与数字地球研究所 | A kind of PM2.5 concentration remote-sensing evaluation method based on satellite polarization technology |
| CN110726653A (en) * | 2019-09-25 | 2020-01-24 | 中国电子科技集团公司第二十七研究所 | PM based on heaven and earth integration information2.5Concentration monitoring method |
| CN111487171A (en) * | 2020-04-28 | 2020-08-04 | 武汉拓宝科技股份有限公司 | Forward and backward combined dual-wavelength dispersed fire smoke detector and method |
Non-Patent Citations (1)
| Title |
|---|
| 曹红业等.基于辐射传输模型的高分二号影像大气校正方法研究.《红外技术》.2020,第534-541页. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112798482A (en) | 2021-05-14 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111487216B (en) | Carbon dioxide flux inversion method and system | |
| CN110726653B (en) | PM based on heaven and earth integration information2.5Concentration monitoring method | |
| CN102628940B (en) | Remote sensing image atmospheric correction method | |
| CN102539336B (en) | Method and system for estimating inhalable particles based on HJ-1 satellite | |
| CN106446564A (en) | Method for remote sensing estimation of net primary productivity of plants | |
| CN110411918B (en) | PM2.5 concentration remote sensing estimation method based on satellite polarization technology | |
| CN113804829A (en) | Atmospheric pollution space-air-ground integrated real-time monitoring system and method | |
| CN106446307B (en) | Aerosol foundation data-based AOD (automated optical inspection) vertical correction effect evaluation method and system | |
| CN101866385B (en) | Target parcel ground surface temperature simulation and optimization method | |
| CN103901420A (en) | Method for dynamic threshold method remote sensing data cloud identification supported by prior surface reflectance | |
| CN111723524A (en) | PM2.5 satellite remote sensing inversion method based on daily variation constraint | |
| CN112798482B (en) | PM2.5 and PM10 estimation method based on satellite remote sensing | |
| CN105261026B (en) | A kind of atmospheric correction processing method of satellite-borne multispectral camera | |
| Yang et al. | Optical sky brightness and transparency during the winter season at dome a antarctica from the gattini-all-sky camera | |
| CN110030934A (en) | The acquisition methods of aerosol optical depth based on MODIS satellite sensor | |
| CN103954974A (en) | Particulate matter optical thickness remote sensing monitoring method used in urban area | |
| CN116822141A (en) | Method for inverting optical thickness of night atmospheric aerosol by utilizing satellite micro-optic remote sensing | |
| CN116558652A (en) | Meteorological satellite thermal infrared channel on-orbit radiation calibration method | |
| CN116519557A (en) | Aerosol optical thickness inversion method | |
| CN106501454A (en) | A kind of satellite remote-sensing monitoring method of jujube tree canopy content of chlorophyll b | |
| CN116185616A (en) | FY-3D MERSI L1B data automatic reprocessing method | |
| CN105678091B (en) | A kind of method of fast inversion muskeg primary productivity | |
| CN111914396B (en) | Sub-grid terrain three-dimensional earth surface solar radiation forced effect rapid parameterization method based on high-resolution DEM data | |
| CN113987778A (en) | Water and soil loss analog value space-time weighting correction method based on field station | |
| CN104217128B (en) | Satellite side-sway imaging air kindred effect analogy method under a kind of rolling topography |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |