CN110726653A - PM based on heaven and earth integration information2.5Concentration monitoring method - Google Patents

PM based on heaven and earth integration information2.5Concentration monitoring method Download PDF

Info

Publication number
CN110726653A
CN110726653A CN201910910189.0A CN201910910189A CN110726653A CN 110726653 A CN110726653 A CN 110726653A CN 201910910189 A CN201910910189 A CN 201910910189A CN 110726653 A CN110726653 A CN 110726653A
Authority
CN
China
Prior art keywords
data
atmospheric
satellite
reflectivity
concentration
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
Application number
CN201910910189.0A
Other languages
Chinese (zh)
Other versions
CN110726653B (en
Inventor
秦臻
张明
赵永辉
李春晓
连芳
徐金晓
段云龙
郭新望
刘颐芳
冯继锋
阚俊峰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Henan Fangda Space Information Technology Co ltd
CETC 27 Research Institute
Original Assignee
CETC 27 Research Institute
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by CETC 27 Research Institute filed Critical CETC 27 Research Institute
Priority to CN201910910189.0A priority Critical patent/CN110726653B/en
Publication of CN110726653A publication Critical patent/CN110726653A/en
Application granted granted Critical
Publication of CN110726653B publication Critical patent/CN110726653B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N15/00Investigating characteristics of particles; Investigating permeability, pore-volume, or surface-area of porous materials
    • G01N15/06Investigating concentration of particle suspensions
    • G01N15/075

Abstract

The invention discloses a PM based on heaven and earth integration information2.5The concentration monitoring method comprises the following steps: a: receiving and acquiring remote sensing satellite image data from a polar orbit satellite receiving station in real time; b: the preprocessing module converts the remote sensing satellite image data into image products of all levels, and finally the brightness temperature value and the earth surface reflectance value are obtained through calculation; c: carrying out inversion on the optical thickness of the atmospheric aerosol to obtain the optical thickness AOD of the atmospheric aerosol in the determined aerosol mode; d: construction of satellite remote sensing monitoring PM2.5A regression model of spatial variation coefficients of concentration; obtaining PM at points to be estimated2.5A concentration calculation formula; e: calculating to obtain PM at all pixel points in the image of the region to be monitored2.5Concentration yi. The invention can realize PM in the designated monitoring area efficiently and accurately2.5And (5) monitoring the concentration.

Description

PM based on heaven and earth integration information2.5Concentration monitoring method
Technical Field
The invention relates to a PM2.5A concentration monitoring method, in particular to a PM based on heaven and earth integration information2.5And (3) a concentration monitoring method.
Background
Atmospheric fine Particulate Matter (PM)2.5) The particles with the aerodynamic diameter less than or equal to 2.5 mu m are mainly composed of organic carbon, element carbon, nitrate, sulfate and the like, are one of important components of air pollution and are also the main reason of dust haze. The atmospheric fine particles are extremely easy to attach toxic and harmful substances, and can be breathed into respiratory tracts, bronchus and lungs by human bodies after being suspended in the space for a long time, thereby causing great influence on the health of the human bodies.
In recent years, various measures are taken by countries and environmental protection departments at all levels to improve the air quality, and the monitoring system is also particularly emphasized. Compared with the traditional ground monitoring means, polar orbit satellite remote sensing is near-ground PM due to the advantages of large-scale spatial scale, continuous earth observation imaging, high time resolution transit and the like2.5Concentration information acquisition provides a number of possible means.
Currently, the PM is systematically and fully automatically monitored by means of remote sensing2.5The method of the concentration of the particulate matter is not popularized yet, and the relevant environmental protection departments lack a PM based on the integrated information of the heaven and the earth2.5The automatic process monitoring method of the particulate matter concentration supports the haze treatment work.
Disclosure of Invention
The invention aims to provide a PM based on heaven and earth integration information2.5The concentration monitoring method can efficiently and accurately realize PM in the designated monitoring area2.5And (5) monitoring the concentration.
The invention adopts the following technical scheme:
PM based on world integration information2.5The concentration monitoring method comprises the following steps:
a: receiving and acquiring remote sensing satellite image data from a polar orbit satellite receiving station in real time;
b: the preprocessing module is used for preprocessing the acquired remote sensing satellite image data, converting the remote sensing satellite image data into image products of each level, and finally obtaining a brightness temperature value and a surface reflectance value through geometric correction and radiance calculation;
c: performing inversion of the optical thickness of the atmospheric aerosol, namely performing cloud identification and eliminating interference of interference pixels, determining an aerosol mode, constructing a lookup table by using a 6S atmospheric radiation transmission model, and performing inversion to obtain the optical thickness of the atmospheric aerosol in the determined aerosol mode by combining all levels of image products obtained in the step B with the lookup table;
d: construction of satellite remote sensing monitoring PM2.5A regression model of spatial variation coefficients of concentration; obtaining PM at point i to be estimated2.5Concentration yiComprises the following steps:
Figure BDA0002214461750000021
i=1,2,…,n;
wherein (u)i,vi) The longitude and latitude coordinates of the point i to be estimated; beta is a0(ui,vi) Is a regression constant; beta is ad(ui,vi) Expressed as observation points (u)i,vi) The d-th regression coefficient of (1), 2, 3 …; beta is ad(ui,vi) Is a function of geographical location, xidD and PM representing point i to be estimated2.5Relevant factor, xidIncluding near-surface extinction coefficient, boundary layer height, relativeHumidity, near-ground temperature, and wind speed; random error term obeys normal distribution epsiloni~N(0,σ2) Random errors are independent of each other, and no sequence-dependent Cov (epsilon)ij)=0(i≠j,(i,j=1,2,…,n));
E: utilizing the PM at the point i to be estimated obtained in the step D2.5Concentration yiThe calculation formula (i) calculates and obtains the PM of all pixel points i in the image of the area to be monitored2.5Concentration yiI.e. the area PM to be monitored2.5And (5) inversion results of the concentration distribution.
And in the step A, receiving and acquiring MODIS load multi-remote sensing satellite image data of Terra and Aqua in real time from a polar orbit satellite receiving station.
The step A comprises the following specific steps:
a1: the receiving and processing system downloads two lines of data of the satellite orbit, and after formatting and decomposition of the two lines of data are completed, the two lines of data are resolved to generate an orbit prediction file;
a2: according to the satellite orbit data in the orbit forecast file, the receiving processing system carries out satellite orbit association and conflict presetting resolution work, and finally generates a tracking receiving plan;
a3: the receiving processing system carries out time coincidence according to the generated tracking receiving plan, analyzes the operation scheduling control instruction after the time coincidence is successful, and sends the operation scheduling control instruction to the antenna equipment to finish the acquisition and tracking of the satellite; meanwhile, the receiving and processing system completes satellite data locking and receiving according to preset satellite data receiving equipment parameters;
a4: the data transmission frequency band of the satellite is an X band, and a receiving processing system adopts 26.25Mbs-1And 15.0Mbs-1Respectively receiving X frequency band signals of Terra and Aqua satellites at code rates, processing and outputting the signals by a down converter to obtain intermediate frequency signals, and outputting digital baseband signals after demodulation by a receiver;
a5: according to the requirements of different satellite data transmission data formats, the receiving and processing system searches a frame synchronization field in a digital baseband signal to complete byte alignment, fault tolerance and filling of a data frame;
a6: the receiving processing system generates a PN code according to a polynomial and carries out XOR operation with image data behind a data frame synchronization field for descrambling according to a scrambling mode of satellite data transmission data; and then according to the satellite data transmission data coding mode, RS decoding and CRC checking are carried out on the image data after descrambling, and finally the original code stream data of the image, namely the remote sensing satellite image data, is obtained.
The step B comprises the following specific steps:
b1: the preprocessing module monitors an original code stream data storage path in real time, and starts to preprocess a new original code stream data file when monitoring that a new original code stream data file is generated and the condition that the new original code stream data file is not occupied and the size of the new original code stream data file is larger than 0KB is met;
b2: the preprocessing module automatically loads xml configuration files corresponding to different satellite loads according to different satellite load names by analyzing the file names of new original code stream data files;
b3: the preprocessing module analyzes the frame format of the high-level on-orbit system of the new original code stream data file according to the data frame formats of different satellite loads, namely automatically calls a standard frame of the new original code stream data file and decodes the standard frame, the decoding comprises the separation of virtual channel data units, the extraction of different information source packets from different virtual channel data units, the output of the virtual channel data units after the reformatting is carried out, and an L0-level data product is generated;
wherein, the L0 level data product refers to original data which is reconstructed and is not processed; removing all communication information in the L0-level data product, wherein the communication information comprises a data frame synchronization field, a communication header and repeated data;
b4: the preprocessing module unpacks and restores a scanning data set from a CCSDS (consultative committee for space data system) package in the L0-level data product by combining the auxiliary data information, namely the preprocessing module corrects the image by extracting ephemeris, attitude, skip second and world standard time information in the auxiliary data information to complete the production of the L1A-level data product;
the auxiliary data information comprises ephemeris, attitude, skip second, world standard time, radiometric calibration parameters, reflectivity parameters and quality parameter information; an L1A-level product, which refers to raw data without any processing, reconstructed from data, and having time reference and geographic coordinate parameter information;
b5: on the basis of an L1A-grade data product, a preprocessing module performs positioning and calibration operation according to different satellite loads by using radiometric calibration parameters, reflectivity parameters and quality parameter information in auxiliary data information and combining geographic coordinate parameter information in the L1A-grade product to generate a data set comprising longitude and latitude information, detector height and atmospheric apparent reflectivity and radiance of different channels, namely the L1B-grade data product;
the L1B-level data product refers to a data set which is subjected to positioning and calibration processing on the basis of the L1A-level data product;
b6: the preprocessing module automatically calls an L1B-grade data product containing the atmospheric apparent reflectivity, radiance and longitude and latitude information, and finally obtains a brightness temperature value and an earth surface reflectivity value through geometric correction and radiance calculation.
In the step C, the aerosol mode is a continental aerosol mode.
The step C comprises the following specific steps:
c1: the cloud identification module sets corresponding threshold values from the atmospheric apparent reflectance values and the brightness temperature values of the different channels obtained in the step B6 to respectively identify interference pixels and eliminate the interference of the interference pixels, wherein the interference pixels comprise cloud, fog and high-reflectivity ground objects;
the channels of the MODIS used are 1, 2, 3, 18, 20, 22, 26, 31 and 35, ρRApparent reflectivity of the atmosphere, rho, for the visible red band of the channel 1BApparent reflectivity of the atmosphere, ρ, for the channel 3 visible blue bandNR2Apparent reflectivity of atmosphere in near infrared band of channel 2, pNR18Apparent reflectivity, p, of the atmosphere in the near infrared band of the channel 18NR26Is the atmospheric apparent reflectivity, BT, of the near infrared band of the channel 2620、BT22、BT31And BT35Are respectively a throughLuminance temperatures in the thermal infrared bands of lanes 20, 22, 31 and 35;
apparent reflectance ρ of the atmosphere in the visible red band of the channel 1RIf the area is more than 0.25, judging that the area is a cloud, fog and high-reflectivity ground object; apparent reflectance ρ of the atmosphere in the visible red band of the channel 1RWhen the reflection coefficient is less than or equal to 0.25, and the apparent reflectivity rho of the atmosphere in the visible blue wave band of the channel 3BIf the area is larger than 0.28, judging that the area is a cloud, fog and high-reflectivity ground object;
if the luminance temperature BT of the thermal infrared band of the channel 3535If the cloud area is not less than 230 and not more than 250, determining that the area is cloud;
if the ratio rho of the atmospheric apparent reflectivity of the near infrared band of the channel 2 to the visible red band of the channel 1NR2RIs [0.9, 1.6 ]]The ratio ρ of the apparent reflectivity of the atmosphere of the near infrared band of the channel 18 to the near infrared band of the channel 2NR18NR2Is [0.28, 0.5 ]]If so, judging the area to be cloud;
if the channel 26 has an apparent atmospheric reflectance ρNR26If the area is larger than 0.05, judging the area to be cloud;
if the luminance temperature BT of the thermal infrared band of the channel 2020Brightness temperature BT in thermal infrared band with channel 2222Has a difference of [2.8, 7.6 ]]Or luminance temperature BT in thermal infrared band of channel 2020Brightness temperature BT in thermal infrared band with channel 3131Has a difference of [10, 19 ]]Judging that the area is cloud and fog;
c2: determining that the aerosol mode is continental and constructing a lookup table;
the determined continental aerosol mode, the determined aerosol content parameter, the determined satellite observation geometric parameter, the determined central wavelength of the waveband and the determined detector height are substituted into a 6S atmospheric radiation transmission model, and atmospheric parameters under different wavebands and different aerosol optical thicknesses are obtained through calculation, so that a lookup table is constructed and obtained; selecting a received data atmospheric mode, spectral parameters and earth surface reflectivity types from 6S radiation transmission model parameter setting options;
wherein the aerosol mode is a continental aerosol mode, and the aerosol content parameter is set to 6 atmospheric aerosol optical thickness values at a wavelength of 0.55 μm: namely 0, 0.25, 0.5, 1.0, 1.5 and 1.95; setting and selecting 9 sun zenith angles of 0 degree, 6 degrees, 12 degrees, 24 degrees, 35.2 degrees, 48 degrees, 54 degrees, 60 degrees and 66 degrees according to the set geometric parameters of the satellite, observing 12 zenith angles at intervals of 6 degrees in the range of [0 degrees and 66 degrees ], and observing 16 opposite azimuth angles at intervals of 12 degrees in the range of [0 degrees and 180 degrees ]; the central wavelength of the wave band is 0.47 μm and 0.66 μm; the probe height was obtained from the L1B data product obtained in step B5;
c3: according to the atmospheric radiation transmission principle, when a standard method is adopted for land aerosol inversion, the atmospheric apparent reflectivity rho is under the assumption that the surface Lambert and the atmospheric level are uniformTOAExpressed as:
Figure BDA0002214461750000051
in the formula (I), the compound is shown in the specification,
Figure BDA0002214461750000055
θsand thetavRespectively a sun zenith angle and an observation zenith angle, phi is a relative azimuth angle, and T is an atmospheric transmittance; s is the hemispherical reflectivity of the lower atmospheric boundary, ρ0Is the path radiation term equivalent reflectivity of the atmosphere, psIs the surface reflectivity; the lower corner symbol s represents the sun, v represents the observation;
selecting the atmospheric parameter closest to the satellite observation geometry from the lookup table obtained in the step C2, and performing linear interpolation according to the L1B-level data product to obtain the hemispherical reflectivity S of the atmospheric lower boundary and the equivalent reflectivity rho of the atmospheric path radiation term in the atmospheric parameters under different wave bands and different aerosol optical thicknesses0And the numerical value of the atmospheric transmittance T, searching the earth surface reflectivity data under the clear weather state from the earth surface reflectivity data of MODIS multi-day history, constructing a clear weather earth surface reflectivity library, and obtaining the earth surface reflectivity rho under the clear weather state in a set time periodsAnd simulating the real atmospheric apparent reflectivity rho by adopting least square fitting through the formula (1)TOALinear interpolation is carried out on the real atmospheric appearance to finally obtain the atmospheric gas solutionThe optical thickness of the glue.
In the step D, the particulate matter inversion model is adopted based on the PM of the geographical weighted regression2.5An inverse model, the construction principle of which is based on near-surface PM2.5The prior relationship existing between concentration and atmospheric aerosol optical thickness AOD, i.e. first according to the near-surface PM2.5PM is obtained through derivation of prior relation between concentration and atmospheric aerosol optical thickness AOD2.5The logarithmic form of the multiple linear regression model of (1);
near-surface PM2.5The prior relationship between concentration and atmospheric aerosol optical thickness AOD is:
Figure BDA0002214461750000052
where ρ is PM2.5R is the effective radius, Q is the extinction efficiency factor, AOD is the atmospheric aerosol optical thickness, HPBL is the boundary layer height, f (RH) is the aerosol scattering moisture absorption growth factor, generally expressed as a function of relative humidity f (RH) 1/(1.0-RH/100); taking natural logarithm from both sides of the above formula to obtain PM2.5The logarithmic form of the multiple linear regression model of (a) is:
wherein, beta1,β2,β3Is a regression coefficient, beta1=1,β2=-1,β3=-1;β0Is a regression constant;
embedding the geographic location of data into PM2.5Obtaining the PM of the satellite remote sensing monitoring from the regression coefficient in the logarithmic form of the multivariate linear regression model2.5A regression model of spatial variation coefficients of concentration;
satellite remote sensing monitoring PM2.5The regression model of the spatial variation coefficient of concentration is:
Figure BDA0002214461750000054
PM at point i to be estimated2.5Concentration yiComprises the following steps:
wherein (u)i,vi) The longitude and latitude coordinates of the point i to be estimated; beta is a0(ui,vi) Is a regression constant; beta is ad(ui,vi) Expressed as observation points (u)i,vi) The d-th regression coefficient of (1), 2, 3 …; beta is ad(ui,vi) Is a function of geographical location, xidD and PM representing point i to be estimated2.5Relevant factor, xidIncluding near-ground extinction coefficient, boundary layer height, relative humidity, near-ground temperature and wind speed; random error term obeys normal distribution epsiloni~N(0,σ2) Random errors are independent of each other, and no sequence-dependent Cov (epsilon)ij)=0(i≠j,(i,j=1,2,…,n))。
Further comprising step F;
f: utilizing the PM of the area to be monitored obtained in the step E2.5The inversion result of the concentration distribution is used for manufacturing PM of the area to be monitored2.5Concentration distribution thematic map.
The method comprises the steps of obtaining remote sensing satellite image data, converting the remote sensing satellite image data into image products of various levels, obtaining a brightness temperature value and a surface reflectance value through geometric correction and radiance calculation, and then performing inversion on the optical thickness of the atmospheric aerosol by using the obtained brightness temperature value and surface reflectance value to obtain the optical thickness AOD of the atmospheric aerosol; then the invention monitors PM by constructing satellite remote sensing2.5The concentration space change coefficient regression model calculates to obtain PM at the point i to be estimated2.5Concentration yiAnd finally obtaining the PM of the area to be monitored2.5The PM based on the integrated information of the heaven and earth is realized by the inversion result of the concentration distribution2.5And (5) monitoring the concentration. In the process, different thresholds are set based on different reflection characteristics of image spectrum wave bands, and cloud identification is passedThe other module detects the interference pixels in the image, effectively reduces the influence of cloud, fog and high-reflectivity ground objects on the inversion accuracy of the optical thickness of the aerosol, and monitors PM in satellite remote sensing2.5Embedding the geographic location of the data into the PM in a regression model of the spatial variation coefficient of concentration2.5Among the regression coefficients in the logarithmic form of the multiple linear regression model, the regression coefficients are changed along with the geographic space position, the inversion result is prevented from generating larger errors due to single regression coefficient, the defects that meteorological parameters and earth surface types in different areas cannot be introduced into the existing linear regression model are also avoided, and the PM is further improved2.5Accuracy of concentration monitoring.
Drawings
FIG. 1 is a schematic flow diagram of the present invention;
fig. 2 is a schematic view of an identification process of the cloud identification module according to the present invention.
Detailed Description
The invention is described in detail below with reference to the following figures and examples:
as shown in FIG. 1, the PM based on the integration of heaven and earth information of the invention2.5The concentration monitoring method comprises the following steps:
a: receiving and acquiring remote sensing satellite image data from a polar orbit satellite receiving station in real time;
in the invention, MODIS load multi-remote sensing satellite image data of Terra and Aqua are received and acquired in real time from a polar orbit satellite receiving station. Terra and Aqua are two important satellites in the United states earth observation system plan, Terra is the morning star, Aqua is the afternoon star, and MODIS is a medium resolution imaging spectrometer carried on Terra and Aqua satellites.
In this embodiment, step a includes the following specific steps:
a1: and the receiving and processing system downloads two lines of data of the satellite orbit, and after formatting and decomposition of the two lines of data are completed, the two lines of data are resolved to generate an orbit prediction file.
A2: and according to the satellite orbit data in the orbit prediction file, the receiving and processing system performs satellite orbit association and conflict presetting resolution work, and finally generates a tracking receiving plan.
A3: the receiving processing system carries out time coincidence according to the generated tracking receiving plan, analyzes the operation scheduling control instruction after the time coincidence is successful, and sends the operation scheduling control instruction to the antenna equipment to finish the acquisition and tracking of the satellite; meanwhile, the receiving processing system completes satellite data locking and receiving according to preset satellite data receiving equipment parameters.
A4: the data transmission frequency band of the satellite is an X band, and a receiving processing system adopts 26.25Mbs-1And 15.0Mbs-1The code rate of the satellite receiving device respectively receives X frequency band signals of Terra and Aqua satellites, the signals are processed and output by a down converter to obtain intermediate frequency signals, and then the intermediate frequency signals are demodulated by a receiver to output digital baseband signals.
A5: according to the requirements of different satellite data transmission data formats, the receiving and processing system completes byte alignment, fault tolerance and filling of data frames by searching frame synchronization fields in the digital baseband signals.
A6: the receiving processing system generates a PN code (pseudo-Noise pseudo-random code) and carries out exclusive-or operation on image data behind a data frame synchronization field according to a satellite data transmission data scrambling mode and a polynomial to descramble; and then according to the satellite data transmission data coding mode, RS (Reed-Solomon) decoding and CRC (cyclic redundancy check) are carried out on the image data after descrambling, and finally, the original code stream data of the image, namely the remote sensing satellite image data, is obtained.
In this embodiment, the operator may monitor the working state of the generation and the issue of the entire tracking reception plan and the control response condition of the satellite tracking and receiving device in real time.
B: the preprocessing module is used for preprocessing the acquired remote sensing satellite image data, converting the remote sensing satellite image data into image products of each level, and finally obtaining a brightness temperature value and a surface reflectance value through geometric correction and radiance calculation;
in the invention, the step B comprises the following specific steps:
b1: the preprocessing module monitors an original code stream data storage path in real time, and starts to preprocess a new original code stream data file when monitoring that a new original code stream data file is generated and the condition that the new original code stream data file is not occupied and the size of the new original code stream data file is larger than 0KB is met;
b2: the preprocessing module automatically loads xml configuration files corresponding to different satellite loads according to different satellite load names by analyzing the file names of new original code stream data files;
b3: the preprocessing module analyzes the frame format of the high-level on-orbit system of the new original code stream data file according to the data frame formats of different satellite loads, namely automatically calls a standard frame of the new original code stream data file and decodes the standard frame, the decoding comprises the separation of virtual channel data units, the extraction of different information source packets from different virtual channel data units, the output of the virtual channel data units after the reformatting is carried out, and an L0-level data product is generated; an L0-level data product, which refers to raw data that has been reconstructed without any processing; all communication information in the L0 level data product is removed, including data frame sync fields, communication headers, and duplicate data.
B4: the preprocessing module unpacks and restores a scanning data set from a CCSDS (consultative committee for space data system) package in the L0-level data product by combining the auxiliary data information, namely the preprocessing module corrects the image by extracting ephemeris, attitude, skip second and world standard time information in the auxiliary data information to complete the production of the L1A-level data product; the auxiliary data information comprises ephemeris, attitude, skip second, world standard time, radiometric calibration parameters, reflectivity parameters and quality parameter information; level L1A product, refers to raw data without any processing, reconstructed from data, with time reference and geographical coordinate parameter information.
B5: on the basis of an L1A-level data product, a preprocessing module performs positioning and calibration operation according to different satellite loads by using radiometric calibration parameters, reflectivity parameters and quality parameter information in auxiliary data information and combining geographic coordinate parameter information in the L1A-level product to generate a data set comprising longitude and latitude information, detector height and apparent reflectivity and radiance of different channels, namely the L1B-level data product. The L1B level data product refers to a data set after positioning and scaling processing is carried out on the basis of the L1A level data product.
B6: the preprocessing module automatically calls an L1B-grade data product containing the atmospheric apparent reflectivity, radiance and longitude and latitude information, and finally obtains a brightness temperature value and an earth surface reflectivity value through geometric correction and radiance calculation.
And B, inverting the optical thickness of the aerosol in a certain aerosol mode by an L1B-grade data product obtained after preprocessing in the step B through identifying the cloud pixel, eliminating the interference of cloud, fog and ground objects with high reflectivity and constructing a lookup table, wherein the ground objects with high reflectivity refer to non-selective reflectors with reflectivity not less than 50%. Thus, in the present application, step C is as follows:
c: performing inversion of the optical thickness of the atmospheric aerosol, namely performing cloud identification and eliminating interference of interference pixels, determining an aerosol mode, constructing a lookup table by using a 6S atmospheric radiation transmission model, and performing inversion to obtain the optical thickness of the atmospheric aerosol in the determined aerosol mode by combining all levels of image products obtained in the step B with the lookup table;
atmospheric aerosol Optical thickness, known in english as aod (aerosol Optical depth) or aot (aerosol Optical thickness), is defined as the integral of the extinction coefficient of the medium in the vertical direction and describes the attenuation of light by the aerosol. The optical thickness of the atmospheric aerosol is one of the most important parameters of the aerosol, and the key physical quantity for representing the turbidity degree of the atmosphere is also an important factor for determining the climate effect of the aerosol.
In the application, the cloud identification module sets different thresholds to detect the cloud in the image through different reflection characteristics based on the spectral wave band of the image, and eliminates the influence of the cloud, the fog and high-reflectivity ground objects on the optical thickness inversion of the aerosol.
Aerosol mode refers to a set of well-defined aerosol types, such as continental, marine and urban in standard aerosol models, each aerosol mode consisting of four different particle components combined in certain ratios, as detailed in table 1. The so-called aerosol types are different, i.e., the ratio of aerosol components or combinations thereof is different. In practical inversion, a look-up table is constructed by using a radiation transmission model according to an assumed aerosol mode, and different influences are brought to a final inversion result due to different particle components under different atmospheric conditions, so that the aerosol mode is an important factor influencing the inversion accuracy. Referring to the research results of aerosol moisture absorption growth and closing experiments in the research area, the aerosol mode of the area to be monitored is set to be a continental aerosol mode.
TABLE 1 comparison of compositions of different particles for urban, continental and marine aerosol models
Figure BDA0002214461750000091
After determining the aerosol pattern, a look-up table is first constructed. The lookup table is carried out by setting different satellite observation geometric parameters, atmospheric modes, aerosol content parameters, detector height, spectral parameters, earth surface reflectivity types and the central wavelength of a wave band where observation data are located and adopting a certain radiation transmission equation. The invention uses the 6S radiation transmission model to carry out radiation transmission calculation to construct the lookup table, and the 6S radiation transmission model adopts the latest scattering calculation method, so that the scattering calculation precision of the solar spectrum wave band is obviously improved compared with 5S.
In the invention, the step C comprises the following specific steps:
c1: and B6, setting corresponding thresholds by the cloud identification module according to the atmospheric apparent reflectance values and brightness temperature values of different channels obtained in the step B6 to respectively identify interference pixels and eliminate the interference of the interference pixels, wherein the interference pixels comprise cloud, fog and high-reflectivity ground objects. According to the characteristics of MODIS load and multi-image analysis in North China, threshold values are adjusted and optimized through the atmospheric apparent reflectivity, brightness temperature and spatial variation characteristics of a plurality of visible light bands, near infrared bands and infrared bands, and finally the identification process of the cloud identification module used in the application is shown in FIG. 2.
Among them, the channels of MODIS used are 1, 2, 3, 18, 20, 22, 26, 31 and 35. RhoRApparent reflectivity of the atmosphere, rho, for the visible red band of the channel 1BApparent reflectivity of the atmosphere, ρ, for the channel 3 visible blue bandNR2Apparent reflectivity of atmosphere in near infrared band of channel 2, pNR18Apparent reflectivity, p, of the atmosphere in the near infrared band of the channel 18NR26Is the atmospheric apparent reflectivity, BT, of the near infrared band of the channel 2620、BT22、BT31And BT35Brightness temperatures in the thermal infrared bands of channels 20, 22, 31 and 35, respectively.
In order to further improve the accuracy of cloud identification module identification and reduce the influence of cloud, fog and high-reflectivity ground objects on the inversion accuracy of the optical thickness of the aerosol as far as possible, the method and the device set different atmospheric apparent reflectivity thresholds for the cloud, fog and high-reflectivity ground objects according to the characteristics of visible light, near infrared and infrared of different wave bands under different channels for judgment.
The application utilizes MODIS multi-temporal image analysis, if the atmospheric apparent reflectivity rho of the channel 1 visible light red light wave bandRIf the area is more than 0.25, judging that the area is a cloud, fog and high-reflectivity ground object; apparent reflectance ρ of the atmosphere in the visible red band of the channel 1RWhen the reflection coefficient is less than or equal to 0.25, and the apparent reflectivity rho of the atmosphere in the visible blue wave band of the channel 3BIf the area is larger than 0.28, judging that the area is a cloud, fog and high-reflectivity ground object; the threshold value can effectively judge the conditions of thick clouds, dense fog and high-reflectivity ground objects in the interference pixels;
if the luminance temperature BT of the thermal infrared band of the channel 3535If the cloud area is not less than 230 and not more than 250, determining that the area is cloud; in the threshold setting, in order to eliminate the interference of haze on cloud in the cloud identification process, the analysis shows that the luminance temperature value of the cloud region in the air on the land is lower, and the luminance temperature value of the haze region is relatively higher than that of the cloud region, so that the luminance temperature BT of the thermal infrared band of the channel 35 is higher than that of the cloud region35Is set to [230, 250 ]]So as to accurately distinguish the haze area from the cloud area;
if channel 2 is near-infraredRatio rho of atmospheric apparent reflectivity of wave band and visible red light wave band of channel 1NR2RIs [0.9, 1.6 ]]The ratio ρ of the apparent reflectivity of the atmosphere of the near infrared band of the channel 18 to the near infrared band of the channel 2NR18NR2Is [0.28, 0.5 ]]If so, judging the area to be cloud; it was analyzed that channel 1, channel 2, and channel 18 have significant advantages for determining thin, ragged, and laminated clouds, and therefore ρ is set in this applicationNR2RHas a threshold value of [0.9, 1.6 ]],ρNR18NR2Has a threshold value of [0.28, 0.5 ]]Effectively judging thin clouds, broken clouds and laminated clouds in the interference pixels;
if the channel 26 has an apparent atmospheric reflectance ρNR26If the area is larger than 0.05, judging the area to be cloud; analysis has shown that the channel 26 has a significant advantage for determining the cloud, and therefore ρ is set in this applicationNR26The threshold value of (2) is 0.05, so as to effectively judge the cirrus cloud existing in the interference pixel;
if the luminance temperature BT of the thermal infrared band of the channel 2020Brightness temperature BT in thermal infrared band with channel 2222Has a difference of [2.8, 7.6 ]]Or luminance temperature BT in thermal infrared band of channel 2020Brightness temperature BT in thermal infrared band with channel 3131Has a difference of [10, 19 ]]Judging that the area is cloud and fog; the analysis shows that the identification by using the difference of the brightness and the temperature of the thermal infrared bands of the channel 20 and the channel 22 and the channel 20 and the channel 31 has obvious advantages for determining the low-layer cloud and fog particles, so the BT setting is carried out by the BT self20And BT22Has a difference threshold of [2.8, 7.6 ]]Or BT20And BT31Is [10, 19 ]]And effectively judging the low-layer cloud and fog particles existing in the interference pixel.
In the above process, thick cloud, dense fog, thin cloud, broken cloud, laminated cloud, rolling cloud, low-layer cloud and fog particles are all terms in the field of meteorology.
C2: determining that the aerosol mode is continental and constructing a lookup table;
the determined continental aerosol mode, the determined aerosol content parameter, the determined satellite observation geometric parameter, the determined central wavelength of the waveband and the determined detector height are substituted into a 6S atmospheric radiation transmission model, and atmospheric parameters under different wavebands and different aerosol optical thicknesses are obtained through calculation, so that a lookup table is constructed and obtained; the received data atmospheric mode, spectral parameters, and surface reflectivity type are selected from the 6S radiation transmission model parameter set options.
Wherein the aerosol mode is a continental aerosol mode, and the aerosol content parameter is set to 6 atmospheric aerosol optical thickness values at a wavelength of 0.55 μm: namely 0, 0.25, 0.5, 1.0, 1.5 and 1.95; setting and selecting 9 sun zenith angles of 0 degree, 6 degrees, 12 degrees, 24 degrees, 35.2 degrees, 48 degrees, 54 degrees, 60 degrees and 66 degrees according to the set geometric parameters of the satellite, observing 12 zenith angles at intervals of 6 degrees in the range of [0 degrees and 66 degrees ], and observing 16 opposite azimuth angles at intervals of 12 degrees in the range of [0 degrees and 180 degrees ]; the central wavelength of the wave band is 0.47 μm and 0.66 μm; the probe height was obtained from the L1B data product obtained in step B5.
C3: according to the atmospheric radiation transmission principle, when a standard method is adopted for land aerosol inversion, the atmospheric apparent reflectivity rho is under the assumption that the surface Lambert and the atmospheric level are uniformTOAExpressed as:
in the formula (I), the compound is shown in the specification,
Figure BDA0002214461750000112
θsand thetavRespectively a sun zenith angle and an observation zenith angle, phi is a relative azimuth angle, and T is an atmospheric transmittance; s is the hemispherical reflectivity of the lower atmospheric boundary, ρ0Is the path radiation term equivalent reflectivity of the atmosphere, psIs the surface reflectivity; the lower corner symbol s represents the sun, v represents the observation;
selecting the atmospheric parameter closest to the satellite observation geometry from the lookup table obtained in the step C2, and performing linear interpolation according to the L1B-level data product to obtain the hemispherical reflectivity S of the atmospheric lower boundary and the equivalent reflectivity rho of the atmospheric path radiation term in the atmospheric parameters under different wave bands and different aerosol optical thicknesses0Is permeable to the atmosphereThe value of the excess T. Searching the earth surface reflectivity data under the clear weather state from the earth surface reflectivity data of MODIS multi-day history, constructing a clear weather earth surface reflectivity library, and obtaining the earth surface reflectivity rho under the clear weather state in a set time periodsAnd simulating the real atmospheric apparent reflectivity rho by adopting least square fitting through the formula (1)TOAAnd performing linear interpolation on the real atmospheric appearance to finally obtain the atmospheric aerosol optical thickness AOD.
D: construction of satellite remote sensing monitoring PM2.5A regression model of spatial variation coefficients of concentration;
the particulate matter inversion model adopted by the invention is PM based on geographical weighted regression2.5Inverting the model; the construction principle is based on near-ground PM2.5There is a certain a priori relationship between concentration and atmospheric aerosol optical thickness AOD, namely first according to near-surface PM2.5PM is obtained through derivation of prior relation between concentration and atmospheric aerosol optical thickness AOD2.5The logarithmic form of the multiple linear regression model of (1);
near-surface PM2.5The prior relationship between concentration and atmospheric aerosol optical thickness AOD is:
where ρ is PM2.5R is the effective radius, Q is the extinction efficiency factor, AOD is the atmospheric aerosol optical thickness, HPBL is the boundary layer height, f (RH) is the aerosol scattering moisture absorption growth factor, generally expressed as a function of relative humidity f (RH) 1/(1.0-RH/100); taking natural logarithm from both sides of the above formula to obtain PM2.5The logarithmic form of the multiple linear regression model of (a) is:
Figure BDA0002214461750000122
wherein, beta1,β2,β3Is a regression coefficient, beta1=1,β2=-1,β3=-1;β0Is a regression constant;
PM2.5the logarithmic form of the multiple linear regression model of (1) is a global space regression model, beta1,β2,β3The regression coefficient varies with the geospatial position. However, most of the existing region-based AOD concentration remote sensing inversion does not consider the change characteristics of local parameters and is based on general experience beta1,β2,β31, -1 and-1 are respectively taken. Therefore, the existing linear regression model with a single regression coefficient can cause a large error of an inversion result, and therefore, the method comprises the following steps:
embedding the geographic location of data into PM2.5Obtaining the PM of the satellite remote sensing monitoring from the regression coefficient in the logarithmic form of the multivariate linear regression model2.5A regression model of spatial variation coefficients of concentration;
satellite remote sensing monitoring PM2.5The regression model of the spatial variation coefficient of concentration is:
Figure BDA0002214461750000123
PM at point i to be estimated2.5Concentration yiComprises the following steps:
Figure BDA0002214461750000124
wherein (u)i,vi) The longitude and latitude coordinates of the point i to be estimated; beta is a0(ui,vi) Is a regression constant; beta is ad(ui,vi) Expressed as observation points (u)i,vi) The d-th regression coefficient of (1), 2, 3 …; beta is ad(ui,vi) The inversion method is a function of the geographic position, the coefficient can be different along with the difference of variable geographic positions, and larger errors of inversion results can be effectively avoided; x is the number ofidD and PM representing point i to be estimated2.5Relevant factor, xidThe near-ground extinction coefficient, the boundary layer height, the relative humidity, the near-ground temperature and the wind speed are obtained from ground observation; random errorThe difference terms obey a normal distribution epsiloni~N(0,σ2) Random errors are independent of each other, and no sequence-dependent Cov (epsilon)i,εj)=0(i≠j,(i,j=1,2,…,n));
E: utilizing the PM at the point i to be estimated obtained in the step D2.5Concentration yiThe calculation formula (i) calculates and obtains the PM of all pixel points i in the image of the area to be monitored2.5Concentration yiI.e. the area PM to be monitored2.5The inversion result of the concentration distribution;
the invention also comprises a step F of obtaining the PM of the area to be monitored2.5The inversion result of the concentration distribution is used for manufacturing PM of the area to be monitored2.5A concentration distribution thematic map;
f: utilizing the PM of the area to be monitored obtained in the step E2.5The inversion result of the concentration distribution is used for manufacturing PM of the area to be monitored2.5Concentration distribution thematic map.
In the invention, when the system receives PM of the area to be monitored2.5When the inversion result of the concentration distribution is obtained, the production of thematic maps and thematic reports can be automatically carried out.
The files required by thematic map making comprise: PM of area to be monitored2.5The inversion result of the concentration distribution, the vector file of the display range in the thematic map, the rendering mode configuration file, the scale, the compass and the legend.
The files required by the thematic report production comprise: the method comprises the steps of obtaining an inversion result of PM2.5 concentration distribution of an area to be monitored, a thematic map and a thematic report template.
The specific steps for making the thematic map are as follows:
1. loading PM of area to be monitored2.5Loading the boundary data of the administrative region of the region to be monitored according to the inversion result of the concentration distribution, and reading a rendering configuration file; adopting shp format vector files for administrative region boundary data of a region to be monitored;
2. according to administrative region boundary data of the region to be monitored, PM of the region to be monitored2.5Clipping the inversion result of the concentration distribution;
3. according to the rendering mode configuration file, the PM of the cut area to be monitored is processed2.5Rendering an inversion result of the concentration distribution;
4. adding a scale, a compass, a thematic map name and a text label;
5. extracting metadata of the thematic map, registering the metadata in a database, and outputting the thematic map to a specified path.
The special report is made by the following steps:
1. reading PM of area to be monitored2.5Inversion result of concentration distribution according to PM2.5The pixel values of the pollution classification method are divided into 6 pollution grades, namely, excellent (the pixel value is more than 0 and less than or equal to 35), good (the pixel value is more than 35 and less than or equal to 75), light pollution (the pixel value is more than 75 and less than or equal to 115), moderate pollution (the pixel value is more than 115 and less than or equal to 150), heavy pollution (the pixel value is more than 150 and less than or equal to 250) and severe pollution (the pixel value is more than 250);
2. calculating basic data: effective monitoring area, PM, of each administrative district of the area to be monitored2.5The concentration mean value and the areas occupied by 6 pollution grades of excellent, good, light pollution, moderate pollution, severe pollution and severe pollution in the area to be monitored;
3. calculating statistical data: and (4) carrying out statistical analysis according to the basic data obtained in the step (2), counting the areas of the highest pollution levels and the distribution conditions of the administrative areas, and ranking the administrative areas according to ranking rules.
4. Reading the thematic report template and thematic map, inserting the thematic map, making a statistical table according to statistical data, and generating a text paragraph. And finally, outputting the thematic report to a specified path, extracting metadata of the thematic report, and registering the metadata to a database.
The invention can also read PM of different time period ranges2.5Performing pixel value weighted average on the inversion result of the concentration distribution to obtain PM2.5The data of day, week, month and year of concentration distribution is used to generate PM with different time frequencies of day, week, month and year according to the thematic map and newspaper production steps2.5Concentration distribution and special report.
Furthermore, the invention can also be applied to inverted PMs2.5And analyzing the precision.
The invention is in the production of PM2.5Simultaneous concentration distribution monitoring thematic map and thematic reportAnd carrying out precision analysis on the inversion result registered in the database. Utilization of provincial/municipal/regional scale range internal control foundation site PM2.5The observation value is a reference value, and the arithmetic mean S of all province/city/regional foundation sites in provinces in half hour range of the transit time of the satellite is calculatedn. Meanwhile, according to the longitude and latitude of the national control station, the effective PM measured by the method within the range of 10 kilometers around is calculated2.5The arithmetic mean value is used as the average value P of satellite observation of province/city/region where the ground site is locatedm. Calculate province/city/region PM2.5The accuracy formula of the concentration product is as follows:
Figure BDA0002214461750000141
through the formula calculation analysis, the province/city/region scale PM measured by the invention can be measured2.5The concentration precision result is directly verified, which is beneficial for related departments to intuitively master the PM of related province/city/region2.5A contamination situation.

Claims (8)

1. PM based on world integration information2.5The concentration monitoring method is characterized by comprising the following steps:
a: receiving and acquiring remote sensing satellite image data from a polar orbit satellite receiving station in real time;
b: the preprocessing module is used for preprocessing the acquired remote sensing satellite image data, converting the remote sensing satellite image data into image products of each level, and finally obtaining a brightness temperature value and a surface reflectance value through geometric correction and radiance calculation;
c: performing inversion of the optical thickness of the atmospheric aerosol, namely performing cloud identification and eliminating interference of interference pixels, determining an aerosol mode, constructing a lookup table by using a 6S atmospheric radiation transmission model, and performing inversion to obtain the optical thickness of the atmospheric aerosol in the determined aerosol mode by combining all levels of image products obtained in the step B with the lookup table;
d: construction of satellite remote sensing monitoring PM2.5A regression model of spatial variation coefficients of concentration; is obtained to standEstimating PM at point i2.5Concentration yiComprises the following steps:
Figure FDA0002214461740000011
wherein (u)i,vi) The longitude and latitude coordinates of the point i to be estimated; beta is a0(ui,vi) Is a regression constant; beta is ad(ui,vi) Expressed as observation points (u)i,vi) The d-th regression coefficient of (1), 2, 3 …; beta is ad(ui,vi) Is a function of geographical location, xidD and PM representing point i to be estimated2.5Relevant factor, xidIncluding near-ground extinction coefficient, boundary layer height, relative humidity, near-ground temperature and wind speed; random error term obeys normal distribution epsiloni~N(0,σ2) Random errors are independent of each other, and no sequence-dependent Cov (epsilon)i,εj)=0(i≠j,(i,j=1,2,…,n));
E: utilizing the PM at the point i to be estimated obtained in the step D2.5Concentration yiThe calculation formula (i) calculates and obtains the PM of all pixel points i in the image of the area to be monitored2.5Concentration yiI.e. the area PM to be monitored2.5And (5) inversion results of the concentration distribution.
2. The PM based on heaven and earth integration information according to claim 12.5The concentration monitoring method is characterized in that: and in the step A, receiving and acquiring MODIS load multi-remote sensing satellite image data of Terra and Aqua in real time from a polar orbit satellite receiving station.
3. PM based on heaven and earth integration information according to claim 22.5The concentration monitoring method is characterized in that the step A comprises the following specific steps:
a1: the receiving and processing system downloads two lines of data of the satellite orbit, and after formatting and decomposition of the two lines of data are completed, the two lines of data are resolved to generate an orbit prediction file;
a2: according to the satellite orbit data in the orbit forecast file, the receiving processing system carries out satellite orbit association and conflict presetting resolution work, and finally generates a tracking receiving plan;
a3: the receiving processing system carries out time coincidence according to the generated tracking receiving plan, analyzes the operation scheduling control instruction after the time coincidence is successful, and sends the operation scheduling control instruction to the antenna equipment to finish the acquisition and tracking of the satellite; meanwhile, the receiving and processing system completes satellite data locking and receiving according to preset satellite data receiving equipment parameters;
a4: the data transmission frequency band of the satellite is an X band, and a receiving processing system adopts 26.25Mbs-1And 15.0Mbs-1Respectively receiving X frequency band signals of Terra and Aqua satellites at code rates, processing and outputting the signals by a down converter to obtain intermediate frequency signals, and outputting digital baseband signals after demodulation by a receiver;
a5: according to the requirements of different satellite data transmission data formats, the receiving and processing system searches a frame synchronization field in a digital baseband signal to complete byte alignment, fault tolerance and filling of a data frame;
a6: the receiving processing system generates a PN code according to a polynomial and carries out XOR operation with image data behind a data frame synchronization field for descrambling according to a scrambling mode of satellite data transmission data; and then according to the satellite data transmission data coding mode, RS decoding and CRC checking are carried out on the image data after descrambling, and finally the original code stream data of the image, namely the remote sensing satellite image data, is obtained.
4. PM based on heaven and earth integration information according to claim 32.5The concentration monitoring method is characterized in that the step B comprises the following specific steps:
b1: the preprocessing module monitors an original code stream data storage path in real time, and starts to preprocess a new original code stream data file when monitoring that a new original code stream data file is generated and the condition that the new original code stream data file is not occupied and the size of the new original code stream data file is larger than 0KB is met;
b2: the preprocessing module automatically loads xml configuration files corresponding to different satellite loads according to different satellite load names by analyzing the file names of new original code stream data files;
b3: the preprocessing module analyzes the frame format of the high-level on-orbit system of the new original code stream data file according to the data frame formats of different satellite loads, namely automatically calls a standard frame of the new original code stream data file and decodes the standard frame, the decoding comprises the separation of virtual channel data units, the extraction of different information source packets from different virtual channel data units, the output of the virtual channel data units after the reformatting is carried out, and an L0-level data product is generated;
wherein, the L0 level data product refers to original data which is reconstructed and is not processed; removing all communication information in the L0-level data product, wherein the communication information comprises a data frame synchronization field, a communication header and repeated data;
b4: the preprocessing module unpacks and restores a scanning data set from a CCSDS (consultative committee for space data system) package in the L0-level data product by combining the auxiliary data information, namely the preprocessing module corrects the image by extracting ephemeris, attitude, skip second and world standard time information in the auxiliary data information to complete the production of the L1A-level data product;
the auxiliary data information comprises ephemeris, attitude, skip second, world standard time, radiometric calibration parameters, reflectivity parameters and quality parameter information; an L1A-level product, which refers to raw data without any processing, reconstructed from data, and having time reference and geographic coordinate parameter information;
b5: on the basis of an L1A-grade data product, a preprocessing module performs positioning and calibration operation according to different satellite loads by using radiometric calibration parameters, reflectivity parameters and quality parameter information in auxiliary data information and combining geographic coordinate parameter information in the L1A-grade product to generate a data set comprising longitude and latitude information, detector height and atmospheric apparent reflectivity and radiance of different channels, namely the L1B-grade data product;
the L1B-level data product refers to a data set which is subjected to positioning and calibration processing on the basis of the L1A-level data product;
b6: the preprocessing module automatically calls an L1B-grade data product containing the atmospheric apparent reflectivity, radiance and longitude and latitude information, and finally obtains a brightness temperature value and an earth surface reflectivity value through geometric correction and radiance calculation.
5. The PM based on heaven and earth integration information according to claim 12.5The concentration monitoring method is characterized in that: in the step C, the aerosol mode is a continental aerosol mode.
6. PM based on heaven and earth integration information according to claim 42.5The concentration monitoring method is characterized in that: the step C comprises the following specific steps:
c1: the cloud identification module sets corresponding threshold values from the atmospheric apparent reflectance values and the brightness temperature values of the different channels obtained in the step B6 to respectively identify interference pixels and eliminate the interference of the interference pixels, wherein the interference pixels comprise cloud, fog and high-reflectivity ground objects;
the channels of the MODIS used are 1, 2, 3, 18, 20, 22, 26, 31 and 35, ρRApparent reflectivity of the atmosphere, rho, for the visible red band of the channel 1BApparent reflectivity of the atmosphere, ρ, for the channel 3 visible blue bandNR2Apparent reflectivity of atmosphere in near infrared band of channel 2, pNR18Apparent reflectivity, p, of the atmosphere in the near infrared band of the channel 18NR26Is the atmospheric apparent reflectivity, BT, of the near infrared band of the channel 2620、BT22、BT31And BT35Luminance temperatures in the thermal infrared bands for channels 20, 22, 31 and 35, respectively;
apparent reflectance ρ of the atmosphere in the visible red band of the channel 1RIf the area is more than 0.25, judging that the area is a cloud, fog and high-reflectivity ground object; apparent reflectance ρ of the atmosphere in the visible red band of the channel 1RWhen the reflection coefficient is less than or equal to 0.25, and the apparent reflectivity rho of the atmosphere in the visible blue wave band of the channel 3BIf the area is larger than 0.28, judging that the area is a cloud, fog and high-reflectivity ground object;
if the channel 35 is thermally infraredBrightness temperature BT of wave band35If the cloud area is not less than 230 and not more than 250, determining that the area is cloud;
if the ratio rho of the atmospheric apparent reflectivity of the near infrared band of the channel 2 to the visible red band of the channel 1NR2RIs [0.9, 1.6 ]]The ratio ρ of the apparent reflectivity of the atmosphere of the near infrared band of the channel 18 to the near infrared band of the channel 2NR18NR2Is [0.28, 0.5 ]]If so, judging the area to be cloud;
if the channel 26 has an apparent atmospheric reflectance ρNR26If the area is larger than 0.05, judging the area to be cloud;
if the luminance temperature BT of the thermal infrared band of the channel 2020Brightness temperature BT in thermal infrared band with channel 2222Has a difference of [2.8, 7.6 ]]Or luminance temperature BT in thermal infrared band of channel 2020Brightness temperature BT in thermal infrared band with channel 3131Has a difference of [10, 19 ]]Judging that the area is cloud and fog;
c2: determining that the aerosol mode is continental and constructing a lookup table;
the determined continental aerosol mode, the determined aerosol content parameter, the determined satellite observation geometric parameter, the determined central wavelength of the waveband and the determined detector height are substituted into a 6S atmospheric radiation transmission model, and atmospheric parameters under different wavebands and different aerosol optical thicknesses are obtained through calculation, so that a lookup table is constructed and obtained; selecting a received data atmospheric mode, spectral parameters and earth surface reflectivity types from 6S radiation transmission model parameter setting options;
wherein the aerosol mode is a continental aerosol mode, and the aerosol content parameter is set to 6 atmospheric aerosol optical thickness values at a wavelength of 0.55 μm: namely 0, 0.25, 0.5, 1.0, 1.5 and 1.95; setting and selecting 9 sun zenith angles of 0 degree, 6 degrees, 12 degrees, 24 degrees, 35.2 degrees, 48 degrees, 54 degrees, 60 degrees and 66 degrees according to the set geometric parameters of the satellite, observing 12 zenith angles at intervals of 6 degrees in the range of [0 degrees and 66 degrees ], and observing 16 opposite azimuth angles at intervals of 12 degrees in the range of [0 degrees and 180 degrees ]; the central wavelength of the wave band is 0.47 μm and 0.66 μm; the probe height was obtained from the L1B data product obtained in step B5;
c3: according to the principle of atmospheric radiation transmissionWhen the land aerosol inversion is carried out by adopting a scalar method, the atmospheric apparent reflectivity rho is under the assumption that the surface lambert and the atmospheric level are uniformTOAExpressed as:
Figure FDA0002214461740000041
in the formula (I), the compound is shown in the specification,
Figure FDA0002214461740000042
θsand thetavRespectively a sun zenith angle and an observation zenith angle, phi is a relative azimuth angle, and T is an atmospheric transmittance; s is the hemispherical reflectivity of the lower atmospheric boundary, ρ0Is the path radiation term equivalent reflectivity of the atmosphere, psIs the surface reflectivity; the lower corner symbol s represents the sun, v represents the observation;
selecting the atmospheric parameter closest to the satellite observation geometry from the lookup table obtained in the step C2, and performing linear interpolation according to the L1B-level data product to obtain the hemispherical reflectivity S of the atmospheric lower boundary and the equivalent reflectivity rho of the atmospheric path radiation term in the atmospheric parameters under different wave bands and different aerosol optical thicknesses0And the numerical value of the atmospheric transmittance T, searching the earth surface reflectivity data under the clear weather state from the earth surface reflectivity data of MODIS multi-day history, constructing a clear weather earth surface reflectivity library, and obtaining the earth surface reflectivity rho under the clear weather state in a set time periodsAnd simulating the real atmospheric apparent reflectivity rho by adopting least square fitting through the formula (1)TOAAnd performing linear interpolation on the real atmospheric appearance to finally obtain the optical thickness of the atmospheric aerosol.
7. PM based on heaven and earth integration information according to claim 62.5The concentration monitoring method is characterized in that: in the step D, the particulate matter inversion model is adopted based on the PM of the geographical weighted regression2.5An inverse model, the construction principle of which is based on near-surface PM2.5The prior relationship existing between concentration and atmospheric aerosol optical thickness AOD, i.e. first according to the near-surface PM2.5Concentration and atmosphere gasPM is obtained through deduction of prior relation between optical thickness of adhesive AOD2.5The logarithmic form of the multiple linear regression model of (1);
near-surface PM2.5The prior relationship between concentration and atmospheric aerosol optical thickness AOD is:
Figure FDA0002214461740000051
where ρ is PM2.5R is the effective radius, Q is the extinction efficiency factor, AOD is the atmospheric aerosol optical thickness, HPBL is the boundary layer height, f (RH) is the aerosol scattering moisture absorption growth factor, generally expressed as a function of relative humidity f (RH) 1/(1.0-RH/100); taking natural logarithm from both sides of the above formula to obtain PM2.5The logarithmic form of the multiple linear regression model of (a) is:
wherein, beta1,β2,β3Is a regression coefficient, beta1=1,β2=-1,β3=-1;β0Is a regression constant;
embedding the geographic location of data into PM2.5Obtaining the PM of the satellite remote sensing monitoring from the regression coefficient in the logarithmic form of the multivariate linear regression model2.5A regression model of spatial variation coefficients of concentration;
satellite remote sensing monitoring PM2.5The regression model of the spatial variation coefficient of concentration is:
Figure FDA0002214461740000053
PM at point i to be estimated2.5Concentration yiComprises the following steps:
Figure FDA0002214461740000054
wherein (u)i,vi) The longitude and latitude coordinates of the point i to be estimated; beta is a0(ui,vi) Is a regression constant; beta is ad(ui,vi) Expressed as observation points (u)i,vi) The d-th regression coefficient of (1), 2, 3 …; beta is ad(ui,vi) Is a function of geographical location, xidD and PM representing point i to be estimated2.5Relevant factor, xidIncluding near-ground extinction coefficient, boundary layer height, relative humidity, near-ground temperature and wind speed; random error term obeys normal distribution epsiloni~N(0,σ2) Random errors are independent of each other, and no sequence-dependent Cov (epsilon)i,εj)=0(i≠j,(i,j=1,2,…,n))。
8. PM based on heaven and earth integration information according to claim 72.5The concentration monitoring method is characterized in that: further comprising step F;
f: utilizing the PM of the area to be monitored obtained in the step E2.5The inversion result of the concentration distribution is used for manufacturing PM of the area to be monitored2.5Concentration distribution thematic map.
CN201910910189.0A 2019-09-25 2019-09-25 PM based on heaven and earth integration information2.5Concentration monitoring method Active CN110726653B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910910189.0A CN110726653B (en) 2019-09-25 2019-09-25 PM based on heaven and earth integration information2.5Concentration monitoring method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910910189.0A CN110726653B (en) 2019-09-25 2019-09-25 PM based on heaven and earth integration information2.5Concentration monitoring method

Publications (2)

Publication Number Publication Date
CN110726653A true CN110726653A (en) 2020-01-24
CN110726653B CN110726653B (en) 2022-03-04

Family

ID=69219371

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910910189.0A Active CN110726653B (en) 2019-09-25 2019-09-25 PM based on heaven and earth integration information2.5Concentration monitoring method

Country Status (1)

Country Link
CN (1) CN110726653B (en)

Cited By (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111650128A (en) * 2020-06-08 2020-09-11 电子科技大学 High-resolution atmospheric aerosol inversion method based on surface reflectivity library
CN112378828A (en) * 2020-12-11 2021-02-19 中科三清科技有限公司 Method and device for inverting concentration of atmospheric fine particulate matters based on satellite remote sensing data
CN112504144A (en) * 2020-12-04 2021-03-16 南京大学 Remote sensing estimation method for accumulated snow thickness on sea ice surface
CN112798482A (en) * 2020-12-30 2021-05-14 联雨科技(天津)有限公司 PM2.5 and PM10 estimation method based on satellite remote sensing
CN113065098A (en) * 2021-03-23 2021-07-02 北京航天创智科技有限公司 PM2.5 concentration detection method and device based on satellite remote sensing and ground gas system data
CN114324074A (en) * 2021-12-22 2022-04-12 陕西智星空间科技有限公司 Atmosphere monitoring sailing method based on backward trajectory model and remote sensing image
CN114943142A (en) * 2022-04-29 2022-08-26 中国科学院空天信息创新研究院 Hyperspectral earth surface reflectivity and atmospheric parameter integrated inversion method and device
CN115859026A (en) * 2022-11-18 2023-03-28 二十一世纪空间技术应用股份有限公司 High-resolution near-surface PM2.5 concentration remote sensing inversion method and device
CN116297288A (en) * 2023-05-24 2023-06-23 上海航天空间技术有限公司 Rapid remote sensing inversion method and system for atmospheric methane dry air mixing ratio
CN116466368A (en) * 2023-06-16 2023-07-21 成都远望科技有限责任公司 Dust extinction coefficient profile estimation method based on laser radar and satellite data
CN116698688A (en) * 2023-04-20 2023-09-05 兰州大学 Method for estimating concentration of atmospheric particulates based on double-star of cloud number 4
CN117390971A (en) * 2023-12-08 2024-01-12 北京英视睿达科技股份有限公司 Method and system for estimating NO2 concentration based on heating season
CN117576570A (en) * 2024-01-16 2024-02-20 山东大学 All-weather water vapor inversion elasticity method, system, equipment and medium

Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101419282A (en) * 2008-12-05 2009-04-29 航天恒星科技有限公司 Integration high speed remote sensing data receiving and processing equipment
CN102103204A (en) * 2011-01-26 2011-06-22 环境保护部卫星环境应用中心 Inversion method for land aerosols optical thickness based on environment satellite 1
KR101382507B1 (en) * 2012-10-19 2014-04-10 사단법인대기환경모델링센터 Air quality forecast and management system
CN104535979A (en) * 2014-12-23 2015-04-22 中国科学院遥感与数字地球研究所 Remote sensing inversion method and system for land cloud optical thickness
CN105092575A (en) * 2014-04-22 2015-11-25 株式会社日立制作所 Method and apparatus for evaluating sand duststorm intensity
CN105678085A (en) * 2016-01-12 2016-06-15 环境保护部卫星环境应用中心 PM2.5 concentration estimation method and system
CN109001091A (en) * 2018-07-18 2018-12-14 北京航天宏图信息技术股份有限公司 Satellite remote-sensing monitoring method, device and the computer-readable medium of atmosphere pollution
CN109030301A (en) * 2018-06-05 2018-12-18 中南林业科技大学 Aerosol optical depth inversion method based on remotely-sensed data
CN110163035A (en) * 2018-02-11 2019-08-23 青岛星科瑞升信息科技有限公司 A kind of cloud Shadow recognition method that priori data is supported
CN110186823A (en) * 2019-06-26 2019-08-30 中国科学院遥感与数字地球研究所 A kind of aerosol optical depth inversion method

Patent Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101419282A (en) * 2008-12-05 2009-04-29 航天恒星科技有限公司 Integration high speed remote sensing data receiving and processing equipment
CN102103204A (en) * 2011-01-26 2011-06-22 环境保护部卫星环境应用中心 Inversion method for land aerosols optical thickness based on environment satellite 1
KR101382507B1 (en) * 2012-10-19 2014-04-10 사단법인대기환경모델링센터 Air quality forecast and management system
CN105092575A (en) * 2014-04-22 2015-11-25 株式会社日立制作所 Method and apparatus for evaluating sand duststorm intensity
CN104535979A (en) * 2014-12-23 2015-04-22 中国科学院遥感与数字地球研究所 Remote sensing inversion method and system for land cloud optical thickness
CN105678085A (en) * 2016-01-12 2016-06-15 环境保护部卫星环境应用中心 PM2.5 concentration estimation method and system
CN110163035A (en) * 2018-02-11 2019-08-23 青岛星科瑞升信息科技有限公司 A kind of cloud Shadow recognition method that priori data is supported
CN109030301A (en) * 2018-06-05 2018-12-18 中南林业科技大学 Aerosol optical depth inversion method based on remotely-sensed data
CN109001091A (en) * 2018-07-18 2018-12-14 北京航天宏图信息技术股份有限公司 Satellite remote-sensing monitoring method, device and the computer-readable medium of atmosphere pollution
CN110186823A (en) * 2019-06-26 2019-08-30 中国科学院遥感与数字地球研究所 A kind of aerosol optical depth inversion method

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
吴黎: "《MODIS遥感信息处理方法及应用》", 30 June 2017, 哈尔滨工程大学出版社 *
夏安全: "基于MODIS影像的济南市气溶胶反演研究与系统开发", 《中国优秀博硕士学位论文全文数据库 工程科技Ⅰ辑》 *
邵彦川 等: "基于地理加权回归克里金的中国PM2.5浓度空间制图方法", 《遥感技术与应用》 *
陈辉 等: "基于地理加权模型的我国冬季PM2.5遥感估算方法研究", 《环境科学学报》 *

Cited By (23)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111650128A (en) * 2020-06-08 2020-09-11 电子科技大学 High-resolution atmospheric aerosol inversion method based on surface reflectivity library
CN111650128B (en) * 2020-06-08 2021-07-13 电子科技大学 High-resolution atmospheric aerosol inversion method based on surface reflectivity library
CN112504144A (en) * 2020-12-04 2021-03-16 南京大学 Remote sensing estimation method for accumulated snow thickness on sea ice surface
CN112504144B (en) * 2020-12-04 2021-10-29 南京大学 Remote sensing estimation method for accumulated snow thickness on sea ice surface
CN112378828B (en) * 2020-12-11 2021-09-17 中科三清科技有限公司 Method and device for inverting concentration of atmospheric fine particulate matters based on satellite remote sensing data
CN112378828A (en) * 2020-12-11 2021-02-19 中科三清科技有限公司 Method and device for inverting concentration of atmospheric fine particulate matters based on satellite remote sensing data
CN112798482B (en) * 2020-12-30 2022-09-06 联雨科技(天津)有限公司 PM2.5 and PM10 estimation method based on satellite remote sensing
CN112798482A (en) * 2020-12-30 2021-05-14 联雨科技(天津)有限公司 PM2.5 and PM10 estimation method based on satellite remote sensing
CN113065098A (en) * 2021-03-23 2021-07-02 北京航天创智科技有限公司 PM2.5 concentration detection method and device based on satellite remote sensing and ground gas system data
CN114324074A (en) * 2021-12-22 2022-04-12 陕西智星空间科技有限公司 Atmosphere monitoring sailing method based on backward trajectory model and remote sensing image
CN114943142B (en) * 2022-04-29 2023-11-28 中国科学院空天信息创新研究院 Integrated inversion method and device for hyperspectral earth surface reflectivity and atmospheric parameters
CN114943142A (en) * 2022-04-29 2022-08-26 中国科学院空天信息创新研究院 Hyperspectral earth surface reflectivity and atmospheric parameter integrated inversion method and device
CN115859026A (en) * 2022-11-18 2023-03-28 二十一世纪空间技术应用股份有限公司 High-resolution near-surface PM2.5 concentration remote sensing inversion method and device
CN115859026B (en) * 2022-11-18 2023-12-05 二十一世纪空间技术应用股份有限公司 High-resolution near-ground PM2.5 concentration remote sensing inversion method and device
CN116698688A (en) * 2023-04-20 2023-09-05 兰州大学 Method for estimating concentration of atmospheric particulates based on double-star of cloud number 4
CN116297288A (en) * 2023-05-24 2023-06-23 上海航天空间技术有限公司 Rapid remote sensing inversion method and system for atmospheric methane dry air mixing ratio
CN116297288B (en) * 2023-05-24 2023-08-04 上海航天空间技术有限公司 Rapid remote sensing inversion method and system for atmospheric methane dry air mixing ratio
CN116466368B (en) * 2023-06-16 2023-08-22 成都远望科技有限责任公司 Dust extinction coefficient profile estimation method based on laser radar and satellite data
CN116466368A (en) * 2023-06-16 2023-07-21 成都远望科技有限责任公司 Dust extinction coefficient profile estimation method based on laser radar and satellite data
CN117390971A (en) * 2023-12-08 2024-01-12 北京英视睿达科技股份有限公司 Method and system for estimating NO2 concentration based on heating season
CN117390971B (en) * 2023-12-08 2024-04-12 北京英视睿达科技股份有限公司 Method and system for estimating NO2 concentration based on heating season
CN117576570A (en) * 2024-01-16 2024-02-20 山东大学 All-weather water vapor inversion elasticity method, system, equipment and medium
CN117576570B (en) * 2024-01-16 2024-04-12 山东大学 All-weather water vapor inversion elasticity method, system, equipment and medium

Also Published As

Publication number Publication date
CN110726653B (en) 2022-03-04

Similar Documents

Publication Publication Date Title
CN110726653B (en) PM based on heaven and earth integration information2.5Concentration monitoring method
Martens et al. Evaluating the land-surface energy partitioning in ERA5
Kumar et al. Wintertime characteristics of aerosols at middle Indo-Gangetic Plain: Impacts of regional meteorology and long range transport
Jackson et al. Suomi‐NPP VIIRS aerosol algorithms and data products
Jeong et al. Quality, compatibility, and synergy analyses of global aerosol products derived from the advanced very high resolution radiometer and Total Ozone Mapping Spectrometer
Levy et al. Algorithm for remote sensing of tropospheric aerosol over dark targets from MODIS: Collections 005 and 051: Revision 2; Feb 2009
CN102539336B (en) Method and system for estimating inhalable particles based on HJ-1 satellite
CN102103204B (en) Inversion method for land aerosols optical thickness based on environment satellite 1
US7242803B2 (en) System and method for significant dust detection and enhancement of dust images over land and ocean
US20050012035A1 (en) System and method for significant dust detection and enhancement of dust images over land and ocean
Gupta et al. Applying the Dark Target aerosol algorithm with Advanced Himawari Imager observations during the KORUS-AQ field campaign
Lu et al. An algorithm for estimating downward shortwave radiation from GMS 5 visible imagery and its evaluation over China
Thieuleux et al. Remote sensing of aerosols over the oceans using MSG/SEVIRI imagery
Zavody et al. The ATSR data processing scheme developed for the EODC
Glantz et al. Trends in MODIS and AERONET derived aerosol optical thickness over Northern Europe
Mishra Retrieval of aerosol optical depth from INSAT‐3D imager over Asian landmass and adjoining ocean: Retrieval uncertainty and validation
Barnes et al. An assessment of diurnal and seasonal cloud cover changes over the Hawaiian Islands using Terra and Aqua MODIS
Verstraete et al. Generating 275-m resolution land surface products from the Multi-Angle Imaging Spectroradiometer data
Kim et al. On air pollutant variations in the cases of long-range transport of dust particles observed in central Korea in the leeside of China in 2010
Whitlock et al. WCRP surface radiation budget shortwave data product description, version 1.1
He et al. Modis aerosol optical thickness product algorithm verification and analysis
Ye et al. Improving atmospheric CO2 retrieval based on the collaborative use of Greenhouse gases Monitoring Instrument and Directional Polarimetric Camera sensors on Chinese hyperspectral satellite GF5-02
CN107966747A (en) High-precision area weather forecasting system
CN111914396B (en) Sub-grid terrain three-dimensional earth surface solar radiation forced effect rapid parameterization method based on high-resolution DEM data
Gupta et al. Retrieval of aerosols over Asia from the Advanced Himawari Imager: Expansion of temporal coverage of the global Dark Target aerosol product

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
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20230828

Address after: 450000 Boxue Road, Zhengdong New District, Zhengzhou City, Henan Province, 36

Patentee after: NO.27 RESEARCH INSTITUTE, CHINA ELECTRONICS TECHNOLOGY Group Corp.

Patentee after: Henan Fangda Space Information Technology Co.,Ltd.

Address before: 450000 Boxue Road, Zhengdong New District, Zhengzhou City, Henan Province, 36

Patentee before: NO.27 RESEARCH INSTITUTE, CHINA ELECTRONICS TECHNOLOGY Group Corp.