CN112798482A - 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
- CN112798482A CN112798482A CN202011613388.4A CN202011613388A CN112798482A CN 112798482 A CN112798482 A CN 112798482A CN 202011613388 A CN202011613388 A CN 202011613388A CN 112798482 A CN112798482 A CN 112798482A
- Authority
- CN
- China
- Prior art keywords
- remote sensing
- satellite
- optical thickness
- red
- central wavelength
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 41
- 238000012937 correction Methods 0.000 claims abstract description 54
- 230000003287 optical effect Effects 0.000 claims abstract description 48
- 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 13
- 238000001228 spectrum Methods 0.000 claims abstract description 12
- 230000005540 biological transmission Effects 0.000 claims abstract description 11
- 239000000443 aerosol Substances 0.000 claims abstract description 10
- 238000004364 calculation method Methods 0.000 claims description 33
- 230000005855 radiation Effects 0.000 claims description 14
- 239000000428 dust Substances 0.000 claims description 9
- 239000000203 mixture Substances 0.000 claims description 7
- 239000004576 sand Substances 0.000 claims description 7
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 7
- 238000007781 pre-processing Methods 0.000 claims description 5
- 239000003897 fog Substances 0.000 claims description 4
- 238000004422 calculation algorithm Methods 0.000 claims description 3
- 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
- 238000011835 investigation Methods 0.000 description 3
- 239000003344 environmental pollutant Substances 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
- 239000000126 substance Substances 0.000 description 2
- 238000009825 accumulation Methods 0.000 description 1
- 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
- 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 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. 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 fine particulate matter in air with an equivalent diameter of less than 2.5 μm. 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 known as 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 site 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·AODblue+w2·AODred
wherein AOD is the optical thickness of the atmospheric aerosol, AODblueFor the optical thickness blue component, AOD, of atmospheric aerosolsredThe 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 correctionH；
The height of the mixed layer is corrected to obtain a height correction factor k of the mixed layermix；
Correcting air humidity to obtain air humidity correction factor krh。
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, AODHOptical thickness, k, of atmospheric aerosol after correction of altituderhCorrection factor for air humidity, kmixCorrection 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 aerosolHThe calculation formula of (2) is as follows:
wherein f isHAnd H0To make a correctionAnd 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-SWIR1.600)/(Green+SWIR1.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 ash haze index calculation formula as follows:
NDHI＝(NIR-SWIR1.600)/(NIR+SWIR1.600)
wherein NDHI is a haze index, NIR is an 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＝(TIR6.200-TIR6.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-sim|+|ρred-sim-ρred|)
w2＝|ρbluse-sim-ρbluse|/(|ρbluse-sim-ρbluse|+|ρred-sim-ρred|)
where ρ isbluse-simSimulated apparent reflectance, ρ, of blue light acquired for 6S radiation transmission modeblueApparent reflectance for blue light, ρred-simSimulated apparent reflectance, ρ, of red light acquired for 6S radiation transmission moderedIs 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 provides 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 concentrations of aerosol, PM2.5 and PM10, specifically, satellite remote sensing data is preprocessed to ensure that data of an input model meets calculation requirements, a foundation is laid for atmospheric aerosol optical thickness inversion by calculating relative positions of a sun, a remote sensing satellite and an earth surface research site, interference pixels such as cloud, water, ice and snow, sand dust, fog and the like are filtered, influence on estimation of PM2.5 and PM10 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 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.
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 block diagram of the calculation of the PM2.5 and PM10 estimation method based on satellite remote sensing according to the present 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 the output result of PM10 obtained by the method for estimating PM2.5 and PM10 based on satellite remote sensing.
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 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.
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 content of the first and second substances,
the latitude of the ground research site is shown, delta is solar declination, omegasThe 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, ω isseAn 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＝π/2-θ
Solar time angle omegasThe 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 parameterssolarThe concrete formula is as follows:
ωsolar＝arcsin(cos(δ)*sin(ω)/sin(sza))+π
the satellite zenith angle zva is calculated as:
wherein phi issatIs the longitude of the sub-satellite point and phi is the longitude of the underground study site.
Satellite observation azimuth angle omegasatThe 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-SWIR1.600)/(Green+SWIR1.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 ash haze index calculation formula as follows:
NDHI＝(NIR-SWIR1.600)/(NIR+SWIR1.600)
wherein NDHI is a dust 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＝(TIR6.200-TIR6.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, dust, 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 accumulation index is more than 0.13, or the thermal infrared bright temperature difference between the central wavelength of 6.2 microns and the central wavelength of 6.9 microns is less than 0, and the pixel point with the dust haze index of more than 0.02 is filtered, so that the situation that the characteristic is mixed into PM2.5 and PM10 during estimation to cause a large error of a data result 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 μm2300And (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·AODblue+w2·AODred
wherein AOD is the optical thickness of the atmospheric aerosol, AODblueFor the optical thickness blue component, AOD, of atmospheric aerosolsredThe weighting coefficients of W1, W2 to the spectrum are the atmospheric aerosol optical thickness red light component.
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-sim|+|ρred-sim-ρred|)
w2＝|ρbluse-sim-ρbluse|/(|ρbluse-sim-ρbluse|+|ρred-sim-ρred|)
where ρ isbluse-simSimulated apparent reflectance, ρ, of blue light acquired for 6S radiation transmission modeblueApparent reflectance for blue light, ρred-simSimulated apparent reflectance, ρ, of red light acquired for 6S radiation transmission moderedIs 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 the obtaining the correction factor includes:
altitude correction, namely obtaining the AOD of the optical thickness of the atmospheric aerosol after altitude correctionH；
The height of the mixed layer is corrected to obtain a height correction factor k of the mixed layermix；
Correcting air humidity to obtain air humidity correction factor krh。
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 AODHIs the optical thickness f of the atmospheric aerosol after altitude correctionHAnd H0To correct the coefficient, take f specificallyH＝0.09874，H00.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:
kmix＝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:
krh＝α+β·k+γ·k2
wherein k isrhThe 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 provided by the invention, the humidity rh is 0.2, k is 1.25 according to the formula, and k is used asrh＝α+β·k+γ·k2And 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, AODHOptical thickness, k, of atmospheric aerosol after correction of altituderhCorrection factor for air humidity, kmixCorrection factor, k, for the mixed layer1、k2、k3For a fixed coefficient, the values are different when calculating PM2.5 or PM 10.
Wherein when PM is PM2.5, k is1＝-298、k2-0.59172、k3The PM2.5 concentration was calculated 320 and the result output is shown in fig. 3.
When PM is PM10, k is1＝-527、k2＝-0.60606、k3The 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 (10)
1. A PM2.5 and PM10 estimation method based on satellite remote sensing is characterized by comprising 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.
2. The method for estimating PM2.5 and PM10 based on satellite remote sensing according to claim 1, wherein the step S3 includes calculating the optical thickness of the atmospheric aerosol by calculating the weighting coefficient of the spectrum, and the calculation formula is as follows:
AOD＝w1·AODblue+w2·AODred
wherein AOD is the optical thickness of the atmospheric aerosol, AODblueFor the optical thickness blue component, AOD, of atmospheric aerosolsredThe weighting coefficients of W1, W2 to the spectrum are the atmospheric aerosol optical thickness red light component.
3. The method for estimating PM2.5 and PM10 based on satellite remote sensing according to claim 2, wherein the step S4 includes:
altitude correction, namely obtaining the AOD of the optical thickness of the atmospheric aerosol after altitude correctionH；
The height of the mixed layer is corrected to obtain a height correction factor k of the mixed layermix；
Correcting air humidity to obtain air humidity correction factor krh。
4. The method for estimating PM2.5 and PM10 based on satellite remote sensing according to claim 3, 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, AODHOptical thickness, k, of atmospheric aerosol after correction of altituderhCorrection factor for air humidity, kmixFor 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.
5. 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 aerosolHThe calculation formula of (2) is as follows:
wherein f isHAnd H0H is the altitude of the place to be corrected, and AOD is the optical thickness of the atmospheric aerosol.
6. 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.
7. The method for estimating PM2.5 and PM10 based on satellite remote sensing according to claim 6, wherein the calculation formula for normalizing the snow cover index is as follows:
NDSI＝(Green-SWIR1.600)/(Green+SWIR1.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 ash haze index calculation formula as follows:
NDHI＝(NIR-SWIR1.600)/(NIR+SWIR1.600)
wherein NDHI is a haze index, NIR is an 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＝(TIR6.200-TIR6.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.
8. The method for estimating PM2.5 and PM10 based on satellite remote sensing according to claim 2, 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 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.
9. The method for estimating PM2.5 and PM10 based on satellite remote sensing according to claim 8, wherein the method for calculating the weighting coefficients of the spectrum is as follows:
w1＝|ρred-sim-ρred|/(|ρbluse-sim-ρbluser-sim|+|ρred-sim-ρred|)
w2＝|ρbluse-sim-ρbluse|/(|ρbluse-sim-ρbluse|+|ρred-sim-ρred|)
where ρ isbluse-simSimulated apparent reflectance, ρ, of blue light acquired for 6S radiation transmission modeblueApparent reflectance for blue light, ρred-simSimulated apparent reflectance, ρ, of red light acquired for 6S radiation transmission moderedIs the red apparent reflectance.
10. 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 true CN112798482A (en) | 2021-05-14 |
| CN112798482B 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) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116363047A (en) * | 2022-08-23 | 2023-06-30 | 生态环境部卫星环境应用中心 | Early warning method for pollution risk of straw incineration atmosphere |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20050012035A1 (en) * | 2003-01-21 | 2005-01-20 | Miller Steven D. | System and method for significant dust detection and enhancement of dust images over land and ocean |
| CN102103204A (en) * | 2011-01-26 | 2011-06-22 | 环境保护部卫星环境应用中心 | Inversion method for land aerosols optical thickness based on environment satellite 1 |
| CN103018736A (en) * | 2012-12-03 | 2013-04-03 | 北京航空航天大学 | Satellite-borne remote sensor radiation calibration method based on atmospheric parameter remote sensing retrieval |
| CN103389494A (en) * | 2013-07-24 | 2013-11-13 | 中国科学院南海海洋研究所 | Novel atmospheric correction method of water color remote sensing data of case II water body |
| 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 |
-
2020
- 2020-12-30 CN CN202011613388.4A patent/CN112798482B/en active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20050012035A1 (en) * | 2003-01-21 | 2005-01-20 | Miller Steven D. | System and method for significant dust detection and enhancement of dust images over land and ocean |
| CN102103204A (en) * | 2011-01-26 | 2011-06-22 | 环境保护部卫星环境应用中心 | Inversion method for land aerosols optical thickness based on environment satellite 1 |
| CN103018736A (en) * | 2012-12-03 | 2013-04-03 | 北京航空航天大学 | Satellite-borne remote sensor radiation calibration method based on atmospheric parameter remote sensing retrieval |
| CN103389494A (en) * | 2013-07-24 | 2013-11-13 | 中国科学院南海海洋研究所 | Novel atmospheric correction method of water color remote sensing data of case II water body |
| 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 |
|---|
| 曹红业等: "基于辐射传输模型的高分二号影像大气校正方法研究", 《红外技术》, 30 June 2020 (2020-06-30), pages 534 - 540 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116363047A (en) * | 2022-08-23 | 2023-06-30 | 生态环境部卫星环境应用中心 | Early warning method for pollution risk of straw incineration atmosphere |
| CN116363047B (en) * | 2022-08-23 | 2024-01-30 | 生态环境部卫星环境应用中心 | Early warning method for pollution risk of straw incineration atmosphere |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112798482B (en) | 2022-09-06 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| 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 | |
| CN100362318C (en) | Atmosphere correction method of airosol optical thickness of aeronautical high-spectrum remote-sensing inversion boundary layer | |
| CN110411918B (en) | PM2.5 concentration remote sensing estimation method based on satellite polarization technology | |
| CN106446564A (en) | Method for remote sensing estimation of net primary productivity of plants | |
| CN113804829A (en) | Atmospheric pollution space-air-ground integrated real-time monitoring system and method | |
| Li et al. | Retrieval of the haze optical thickness in North China Plain using MODIS data | |
| CN106407656A (en) | Retrieval method for aerosol optical thickness based on high resolution satellite image data | |
| CN106446307B (en) | Aerosol foundation data-based AOD (automated optical inspection) vertical correction effect evaluation method and system | |
| CN103901420A (en) | Method for dynamic threshold method remote sensing data cloud identification supported by prior surface reflectance | |
| CN101866385B (en) | Target parcel ground surface temperature simulation and optimization method | |
| CN106979911A (en) | The method that PM 2.5 and PM 10 is estimated is carried out using satellite multispectral image data | |
| CN111487216A (en) | Carbon dioxide flux inversion method and system | |
| CN102636143A (en) | Aerosol optical depth remote sensing retrieval method | |
| CN112798482B (en) | PM2.5 and PM10 estimation method based on satellite remote sensing | |
| CN107389617A (en) | The inversion method and equipment of aerosol optical depth based on No. four satellites of high score | |
| CN105261026B (en) | A kind of atmospheric correction processing method of satellite-borne multispectral camera | |
| CN110030934A (en) | The acquisition methods of aerosol optical depth based on MODIS satellite sensor | |
| CN116337701A (en) | Urban high-resolution aerosol optical thickness inversion method based on double-star networking | |
| CN116519557A (en) | Aerosol optical thickness inversion method | |
| CN116822141A (en) | Method for inverting optical thickness of night atmospheric aerosol by utilizing satellite micro-optic remote sensing | |
| Wang et al. | Retrieval of aerosol optical depth for Chongqing using the HJ-1 satellite data | |
| CN116558652A (en) | Meteorological satellite thermal infrared channel on-orbit radiation calibration method | |
| CN114925997B (en) | Method for screening effective data of multispectral sensor of farmland Internet of things |
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 |