CN112964666A - Atmospheric carbon dioxide content calculation method based on earth surface bidirectional reflection - Google Patents
Atmospheric carbon dioxide content calculation method based on earth surface bidirectional reflection Download PDFInfo
- Publication number
- CN112964666A CN112964666A CN202110148447.3A CN202110148447A CN112964666A CN 112964666 A CN112964666 A CN 112964666A CN 202110148447 A CN202110148447 A CN 202110148447A CN 112964666 A CN112964666 A CN 112964666A
- Authority
- CN
- China
- Prior art keywords
- profile
- inversion
- atmospheric
- earth surface
- model
- 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.)
- Pending
Links
- 229910002092 carbon dioxide Inorganic materials 0.000 title claims abstract description 45
- 238000004364 calculation method Methods 0.000 title claims abstract description 24
- 230000002457 bidirectional effect Effects 0.000 title claims abstract description 22
- CURLTUGMZLYLDI-UHFFFAOYSA-N Carbon dioxide Chemical compound O=C=O CURLTUGMZLYLDI-UHFFFAOYSA-N 0.000 title claims abstract description 20
- 239000001569 carbon dioxide Substances 0.000 title claims abstract description 10
- 239000000443 aerosol Substances 0.000 claims abstract description 17
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 6
- 230000003287 optical effect Effects 0.000 claims description 24
- 238000002310 reflectometry Methods 0.000 claims description 21
- 239000011159 matrix material Substances 0.000 claims description 20
- 238000000034 method Methods 0.000 claims description 19
- 230000005855 radiation Effects 0.000 claims description 10
- 230000005540 biological transmission Effects 0.000 claims description 8
- 102000006463 Talin Human genes 0.000 claims description 6
- 108010083809 Talin Proteins 0.000 claims description 6
- 238000004088 simulation Methods 0.000 claims description 4
- OKTJSMMVPCPJKN-UHFFFAOYSA-N Carbon Chemical compound [C] OKTJSMMVPCPJKN-UHFFFAOYSA-N 0.000 claims description 3
- 229910052799 carbon Inorganic materials 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 238000003384 imaging method Methods 0.000 claims description 3
- 239000000203 mixture Substances 0.000 claims description 3
- 238000001228 spectrum Methods 0.000 claims description 3
- 238000013517 stratification Methods 0.000 claims description 3
- 230000001131 transforming effect Effects 0.000 claims description 3
- 238000012937 correction Methods 0.000 claims description 2
- 238000010276 construction Methods 0.000 claims 1
- 238000012545 processing Methods 0.000 abstract description 2
- 230000008878 coupling Effects 0.000 abstract 1
- 238000010168 coupling process Methods 0.000 abstract 1
- 238000005859 coupling reaction Methods 0.000 abstract 1
- 230000000694 effects Effects 0.000 description 3
- 230000008859 change Effects 0.000 description 2
- 238000005034 decoration Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 239000000047 product Substances 0.000 description 2
- 230000001808 coupling effect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 239000005431 greenhouse gas Substances 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 239000013589 supplement Substances 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Chemical compound O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
- G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
- G01N21/25—Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
- G01N21/31—Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry
- G01N21/35—Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry using infrared light
- G01N21/359—Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry using infrared light using near infrared light
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Theoretical Computer Science (AREA)
- Chemical & Material Sciences (AREA)
- Biochemistry (AREA)
- General Health & Medical Sciences (AREA)
- Analytical Chemistry (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Life Sciences & Earth Sciences (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
Abstract
The invention provides an atmospheric carbon dioxide content calculation method based on earth surface bidirectional reflection. Experimental results show that the inversion algorithm provided by the invention can effectively correct inversion errors caused by ground surface two-way reflection and ground-gas coupling, and improves atmospheric CO2And (4) content inversion accuracy. In addition, the invention also aims to realize high-precision atmospheric CO under the complex surface reflection and the environment with larger aerosol concentration in cities and the like2Inversion provides the possibility. It is the current and next generation to monitor urban atmospheric CO2The satellite load data processing mainly based on the emission source provides a high-precision inversion algorithm, and has a wide application prospect.
Description
Technical Field
The invention belongs to the technical field of satellite near-infrared remote sensing data processing, and particularly relates to an atmospheric carbon dioxide content calculation method based on earth surface bidirectional reflection.
Background
CO2Is the main greenhouse gas in the atmosphere, the atmospheric CO2The content of (a) will have a significant influence on global climate change. At present, satellite observation becomes monitoring atmospheric CO2Content prevailing method, but WeiweiThe application of the satellite remote sensing data has extremely high requirements on the observation precision, and the error of a detection result is generally required to be controlled within 1%. Multiple studies have shown that CO is present in high-precision atmospheres2In the inversion process, the surface reflection is not sufficiently evaluated, which causes large inversion errors. The current inversion method treats the earth surface as a uniform lambertian body, and has the following defects:
(1) the earth surface has obvious two-way reflection characteristics, the common inversion method at present assumes the earth surface as an isotropic lambertian body, and although the characteristics of partial earth surface can be approximately expressed, the direction characteristics of the earth surface reflection along with the change of the angle are lost, and the inversion accuracy cannot be further improved.
(2) The atmospheric environment of a part of inversion regions is poor, the probability of obtaining an ideal observation condition with the optical thickness of aerosol less than 0.3 is low, the atmospheric and earth surface radiation coupling effect is more obvious under the condition of a larger optical thickness value of aerosol, and errors are easy to generate. For example, the inversion method in the industry is low in inversion accuracy and unreliable in result for the pixel with the aerosol optical thickness being more than 0.3.
Disclosure of Invention
The purpose of the invention is as follows: the invention aims to solve the technical problem of providing an atmospheric carbon dioxide content calculation method based on earth surface bidirectional reflection aiming at the defects of the prior art so as to obtain high-precision atmospheric CO2And (5) detecting the content.
The invention specifically comprises the following steps:
step 1, inputting observation spectrums of L1 grade 0.76 mu m, 1.58 mu m and 2.0 mu m three-band;
step 2, judging the observation position, if the observation position is in the urban area, executing step 3, otherwise executing step 5;
step 3, counting three-parameter changes of the Ross-Li earth surface reflection model within the year range of an observation point N (generally taking the value of 5) according to the time sequence;
step 4, rewriting three parameters of the Ross-Li earth surface reflection model into a single parameter form;
step 5, constructing a Rahman-pinti-Wirstreet (Rahman-Pinty-Verstraete) RPV surface reflection model, and setting a kernel function parameter;
step 6, according to a data set of a Moderate-resolution Imaging spectrometer MODIS (model-resolution Imaging Spectrophotometer), giving surface reflectance values of three wave bands of an observation point of 0.76 μm, 1.58 μm and 2.0 μm as initial values of inversion;
step 7, constructing temperature, humidity and pressure profiles in the atmosphere and an optical thickness and wavelength index database of the aerosol as inversion initial values;
step 8, constructing atmospheric CO2An initial value profile and a prior covariance matrix;
step 9, performing inversion calculation by using a Gaussian iteration model, judging whether a convergence condition is reached, if the convergence condition is not reached, executing step 10, otherwise executing step 11;
step 11, calculating to obtain atmospheric CO2Column average dry air mix ratio.
In step 2, whether the observation position is in the urban area is judged by the following method:
calculating a Normalized Difference creation Index (NDBI) based on MODIS satellite data:
NDBI＝(Band6-Band2)/(Band6+Band2)
wherein Band2 and Band6 represent the 2 nd wave Band and the 6 th wave Band data of MODIS satellite data respectively;
and if the NDBI is more than 0.1, judging that the observation position is positioned in the urban area.
The step 3 comprises the following steps:
step 3-1, Ross-Li earth surface bidirectional reflection model
The expression of (a) is:
wherein theta issAt the zenith angle of the sun, thetavIn order to observe the zenith angle,
is a relative azimuth angle, λ is a wavelength; k is a radical ofiso、kgeoAnd kvolRespectively isotropic, geometric and bulk scattering kernel functions, fiso、fgeoAnd fvolRespectively an isotropic nuclear coefficient, a geometric optical scattering nuclear coefficient and a volume scattering nuclear coefficient;
step 3-2, mixing fisoExtracting the formula, transforming the earth surface bidirectional model into
Step 3-3, to kernel function coefficient fiso、fgeoAnd fvolAccording to the time sequence, statistics is carried out in a month unit in N years, and the following are respectively calculated:
f′geo＝fgeo/fiso，
f′vol＝fvol/fiso，
obtaining a geometric optical scattering nuclear coefficient f 'after deformation'geoAnd a back-body scattering nuclear coefficient f 'of deformation'volAverage value of (a).
Step 4 comprises the following steps: constructing a Rahmann-pinto-Wirsteret (Rahman-Pinty-Verstraete) RPV surface reflection model:
where ρ is0The parameter is a hot spot parameter of a surface bidirectional reflection model, K is a kernel function, K is an anisotropic parameter of the surface bidirectional reflection model, and theta is an empirical parameter related to a scattering phase function. Will rho0The initial value is set to 0.05, k is 0.75, and theta is-0.1.
The step 5 comprises the following steps:
step 5-1, selecting MODIS afternoon star Aqua surface reflectivity daily product MYD09GA, and selecting the second reflectivity of the wave band as O2Band A with surface reflectivity and band six reflectivity as CO21 the reflectance of the earth surface and the reflectance of the wave band seven as CO2-2 zone surface reflectivity;
step 5-2, counting the data from 2006 to 2016 in 10 years according to longitude, latitude and season (divided into four seasons: spring 2, spring 3 and spring 4, summer 5, spring 6 and spring 7, autumn 8, spring 9 and spring 10, winter 11, spring 12 and spring 1) and making a lookup table; and (3) taking the longitude and latitude of the center of the observation point, searching surrounding 10 multiplied by 10 grid points, and taking the average of all effective values as the prior value of the earth surface reflectivity of the point.
The step 6 comprises the following steps:
step 6-1, selecting an ERA-Interim data set of a European middle-term Weather forecast center (European Centre for Medium-Range Weather strategies) ECMWF, setting the central longitude and latitude coordinates of a projection point to be (x, y) by adopting a bilinear interpolation algorithm, setting the coordinates of four surrounding nearest grid points to be (x1, y1), (x2, y2), (x3, y3) and (x4, y4), and setting the final profile to be:
P＝(x1-x)(y1-y)p1+(x2-x)(y2-y)p2+(x3-x)(y3-y)p3+(x4-x)(y4-y)p4
wherein P is a matrix with 3 rows and 20 rows, 3 rows are respectively the temperature, humidity and pressure profiles of observation points, and P is1，p2，p3，p4Four corner point profiles;
6-2, according to the optical thickness tau of the aerosol at 550nm550And wavelength index α, converted to optical thickness calculations for inversion of three channels:
τ760＝τ550·1.382-α，
τ1580＝τ550·2.873-α，
τ2050＝τ550·3.727-α。
the step 7 comprises the following steps:
step 7-1, CO2The profile prior values are from the Carbon tracking plan CT2017(Carbon tracker) dataset, CO2The profile only contains CO below 50km of atmosphere2Profiles, for use in radiation calculations, using standard CO2Profile and standard pressure profile, for CO2The profile is expanded to 100km of the top of the atmospheric layer;
step 7-2, covariance matrix XcovIs calculated as follows:
wherein, the elements of the p row and the p column of the covariance matrix are shown, n represents the number of profile samples, represents the concentration value of the i layer of the k profile, represents the concentration mean value of the i layer, the subscripts i and j represent the number of the profile layers, and p is the total number 46 of the profile layers; the parameters k, i and j have a value in the range of [1,46 ].
The step 8 comprises the following steps:
step 8-1, calculating radiation transmission y:
y＝F(x,b)+ε
wherein F (x) is the simulation calculation result of the forward radiation transmission model; x is the state vector to be inverted, including carbon dioxide CO2Content, aerosol optical thickness, and surface parameters; b is an auxiliary parameter vector; epsilon is an error term between the forward model and the measured data;
step 8-2, calculating a cost function J (x):
wherein F (x) is the simulation calculation result of the forward radiation transmission model, SεTo measure the error covariance matrix, xaAs a priori information, SaIs a prior error covariance matrix; cost functionJ (X) is less than X (generally 0.001), convergence is judged, and step 10 is executed;
if not, according to the Gaussian iterative formula, after updating the state vector x, repeatedly executing the step 8:
in the formula, x(i+1)For the state vector, K, to continue to participate in the calculation after one iterationiIs the weight matrix and gamma is the iterative correction factor. And when the state vector x is calculated, the minimum value of J (x) is ensured.
In step 9, the CO is iteratively updated2Profile, surface reflection parameters, aerosol optical thickness parameters. CO22The initial value of the profile is the value in step 7, and the iteration mode is that the last profile is multiplied by a scaling factor, such as 390ppm (1+ 0.01); the initial value of the surface reflection parameter is the input of the step 4, and the iteration mode is the result of the previous round plus increment, such as 0.05+ 0.01; the initial value of the aerosol optical thickness comes from step 6, iteratively by adding an increment, such as 0.1+0.02, to the previous round of results.
In step 10, the CO to be exported2And O2Conversion of the profile into atmospheric CO2Column average dry air mixing ratio XCO2The formula is as follows:
in the formula NCO2(z) is CO2Number concentration in atmosphere with high degree of stratification, NO2(z) is O2The number concentrations are highly stratified in the atmosphere. dz is the sign of the integral, integrating z.
Compared with the prior art, the invention has the following remarkable advantages: (1) the inversion error caused by the earth surface bidirectional reflection effect during inversion can be effectively corrected. (2) The model which can represent the earth surface bidirectional reflection only by a single parameter is designed for the urban area, and the inversion in the urban area is possible. (3) The effective inversion can also be carried out for the observation condition that the optical thickness of the aerosol is more than 0.3. (4) The initial values of inversion are all from the satellite historical data set, and are used as initial inversion guess values after preprocessing, so that the true values are reasonably solved, and the convergence speed of the inversion algorithm is accelerated.
Drawings
The foregoing and/or other advantages of the invention will become further apparent from the following detailed description of the invention when taken in conjunction with the accompanying drawings.
FIG. 1 is an atmospheric CO of the present invention taking into account bi-directional reflection from the earth's surface2Flow chart of the inversion method.
Fig. 2 is a diagram of a calculation result of judging the position of the observation point.
Fig. 3 is a time series chart of statistical five year period urban surface reflection parameters.
FIG. 4 is a map of surface reflection parameters modified to a single parameter.
Fig. 5 is a table for looking up the initial values of the reflectivity.
FIG. 6 CO acquisition from CT20172And calculating the concentration profile.
Fig. 7 is a diagram of the covariance matrix calculated.
Fig. 8 is a comparison graph of the effect of the bi-directional reflection of the earth's surface with and without consideration.
Detailed Description
With reference to fig. 1, the invention provides an atmospheric carbon dioxide content calculation method based on earth surface bidirectional reflection, which comprises the following steps:
step 1, inputting satellite L1 level spectrum data which can be GMI, GOSAT and other satellite observation data;
step 2, detecting whether the observation position is in the urban area comprises the following specific steps:
firstly, calculating a normalized building index NDBI based on MODIS satellite data, wherein the calculation formula is as follows:
NDBI＝(Band6-Band2)/(Band6+Band2)
wherein Band2 and Band6 represent the 2 nd wave Band and the 6 th wave Band data of MODIS satellite data respectively;
secondly, referring to fig. 2, if NDBI is greater than 0.1, it is determined that it is a city area, and step 3 is performed, otherwise, it is determined that it is a non-city area, and step 5 is performed.
Step 3, rewriting the earth surface bidirectional reflection model into a single-parameter form, and specifically realizing the method by the following steps:
first, the expression of the earth surface bidirectional reflection model is:
wherein theta issAt the zenith angle of the sun, thetavIn order to observe the zenith angle,
λ is the wavelength, relative azimuth. k is a radical ofiso、kgeoAnd kvolRespectively isotropic, geometric and bulk scattering kernel functions, fiso、fgeoAnd fvolRespectively an isotropic nuclear coefficient, a volume scattering nuclear coefficient and a geometric optical scattering nuclear coefficient;
second, f in the first stepisoExtracting a formula, and transforming the earth surface bidirectional reflection model into:
step 3, counting three parameters of the earth surface reflectivity according to the time sequence; for kernel function coefficient fiso、fgeoAnd fvolAccording to the time sequence, statistics are carried out in units of months in five-year time, and the following are respectively calculated by combining the graph 3:
f′geo＝fgeo/fiso
f′vol＝fvol/fiso
and calculating f'geoAnd f'volAs shown in fig. 4.
And 5, constructing an RPV earth surface reflection model, comprising the following steps:
will rho0The initial value is set to 0.05, k is 0.75, and theta is-0.1, and the calculation formula is as follows:
wherein theta issAt the zenith angle of the sun, thetavIn order to observe the zenith angle,
λ is the wavelength, relative azimuth. Rho0For the hot spot parameters, k is the anisotropy parameter and Θ is the empirical parameter for the scattering phase function.
And6, giving reflectance values of 0.76 μm, 1.58 μm and 2.0 μm three-band observation points according to the MODIS data set by the specific steps of:
firstly, selecting a MODIS afternoon star Aqua surface reflectivity daily product MYD09GA, and selecting a waveband secondary reflectivity thereof as O2Band A with surface reflectivity and band six reflectivity as CO21 the reflectance of the earth surface and the reflectance of the wave band seven as CO2-2 zone surface reflectivity;
secondly, data from 2006 to 2016 (year 10) are counted according to longitude, latitude and season (four seasons: spring, summer, 5, 6 and 7, autumn, 11, 12 and 1, winter) to make a look-up table, and the look-up table is combined with the graph shown in FIG. 5, wherein (a) is O2-zone A surface reflectance, (b) is CO2-1 zone surface reflectance, (c) is CO2-2 zone surface reflectivity. And (3) taking the longitude and latitude of the center of the observation point, searching surrounding 10 multiplied by 10 grid points, and taking the average of all effective values as the prior value of the earth surface reflectivity of the point.
Step 7, constructing temperature, humidity and pressure profiles in the atmosphere, and an optical thickness and wavelength index database of the aerosol, and inputting initial values:
firstly, selecting an ECMWF data set, adopting a bilinear interpolation algorithm, assuming that the longitude and latitude coordinates of the center of a projection point are (x, y), the coordinates of the four nearest grid points around are respectively (x1, y1), (x2, y2), (x3, y3) and (x4, y4), and the final profile P is:
P＝(x1-x)(y1-y)p1+(x2-x)(y2-y)p2+(x3-x)(y3-y)p3+(x4-x)(y4-y)p4
wherein P comprises the temperature, humidity, pressure profile of the observation point, P1，p2，p3，p4Four corner point profiles.
Second, according to optical thickness tau of aerosol 550nm550And wavelength index α, converted to optical thickness calculations for inversion of three channels:
τ760＝τ550·1.382-α
τ1580＝τ550·2.873-α
τ2050＝τ550·3.727-α
step 8, constructing atmospheric CO2The specific steps of selecting an inversion initial value and a prior covariance matrix by a profile library are as follows:
first, in conjunction with FIG. 6, profile data and standard atmospheric CO were selected for a CT2017 dataset of 13:30 per day2Concentration Profile data, CO2The concentration interpolation is 0km,1km,2km,3km,4km,5km,6km,7km,8km,9km,10km,11km,12km,13km,14km,15km,16km,17km,18km,19km,20km,25km,30km,40km and 50km, and a profile below 50km is formed;
second, utilizing standard CO2The profile and the standard pressure profile supplement missing parts in the height of more than 50km, and the specific calculation formula is as follows:
wherein xco2' is the adjusted profile, x50kmFor profiles above 50km for the CT2017 data set,
concentration values at 50km of the standard profile.
Third, the covariance matrix XcovIs calculated as follows:
wherein x represents the covariance matrix element, n represents the number of profile samples, and xk,iRepresenting the concentration value of the ith layer of the kth profile,
the concentration mean of the ith layer is represented, the subscripts i and j each represent the number of profile layers, and p is the total number of profile layers, as shown in fig. 7.
Step 9, performing inversion calculation by using a Gaussian iteration model, judging whether a convergence condition is reached, and if the convergence condition is not reached, executing the step 10:
firstly, calculating radiation transmission:
y＝F(x,b)+ε
wherein F is a forward radiation transmission model, and x is a state vector to be inverted including CO2Content, aerosol optical thickness and surface parameters, and b is an auxiliary parameter vector comprising temperature, water vapor, pressure and elevation. ε is the error term between the forward model and the measured data.
Second, calculate the cost function j (x):
in the formula SεFor measuring the covariance matrix of errors, take values of [1.0, 1.0%]，xaAs a priori information, SaTaking the value of the prior error covariance matrix as X in step 8cov. A cost function j (x) less than 0.001 is judged to be convergent,step 10 is performed.
Thirdly, if the state vector x is not converged, updating the state vector x according to a Gaussian iteration formula, and then executing the step 9:
in the formula NCO2(z) is CO2Number concentration in atmosphere with high degree of stratification, NO2(z) is O2The number concentrations are highly stratified in the atmosphere. In combination with fig. 8, after the earth surface bidirectional reflection effect is considered, the correlation between the inversion result and the standard value is improved from 0.56 to 0.85, the average error is reduced from 1.75ppm to 1.13ppm, and the inversion accuracy is obviously improved.
The invention provides a method for calculating the content of atmospheric carbon dioxide based on earth surface bidirectional reflection, and a plurality of methods and ways for implementing the technical scheme, and the above description is only a preferred embodiment of the invention, and it should be noted that, for those skilled in the art, a plurality of improvements and decorations can be made without departing from the principle of the invention, and these improvements and decorations should also be regarded as the protection scope of the invention. All the components not specified in the present embodiment can be realized by the prior art.
Claims (10)
1. An atmospheric carbon dioxide content calculation method based on earth surface bidirectional reflection is characterized by comprising the following steps:
step 1, inputting observation spectrums of L1 grade 0.76 mu m, 1.58 mu m and 2.0 mu m three-band;
step 2, judging the observation position, if the observation position is in the urban area, executing step 3, otherwise executing step 5;
step 3, counting three-parameter changes of the Ross-Li earth surface reflection model within the N-year range of the observation points according to the time sequence; rewriting three parameters of the Ross-Li earth surface reflection model into a single parameter form;
step 4, constructing an RPV (resilient packet virtual) surface reflection model and setting kernel function parameters;
step 5, surface reflectance values of 0.76 mu m, 1.58 mu m and 2.0 mu m observation points are given according to the MODIS data set of the medium-resolution imaging spectrometer and are used as initial values of inversion;
step 6, constructing temperature, humidity and pressure profiles in the atmosphere and an optical thickness and wavelength index database of the aerosol as inversion initial values;
step 7, construction of atmospheric CO2An initial value profile and a prior covariance matrix;
step 8, performing inversion calculation by using a Gaussian iteration model, judging whether a convergence condition is reached, if the convergence condition is not reached, executing step 9, otherwise executing step 10;
step 9, iteratively updating the parameters, and returning to the step 8;
step 10, calculating to obtain atmospheric CO2Column average dry air mix ratio.
2. The method according to claim 1, wherein in step 2, it is determined whether the observed location is in an urban area by:
calculating a normalized building index NDBI based on MODIS satellite data:
NDBI＝(Band6-Band2)/(Band6+Band2)
wherein Band2 and Band6 represent the 2 nd wave Band and the 6 th wave Band data of MODIS satellite data respectively;
and if the NDBI is more than 0.1, judging that the observation position is positioned in the urban area.
3. The method of claim 2, wherein step 3 comprises:
step 3-1, Ross-Li earth surface bidirectional reflection model
The expression of (a) is:
wherein theta issAt the zenith angle of the sun, thetavIn order to observe the zenith angle,
is a relative azimuth angle, λ is a wavelength; k is a radical ofiso、kgeoAnd kvolRespectively isotropic, geometric and bulk scattering kernel functions, fiso、fgeoAnd fvolRespectively an isotropic nuclear coefficient, a geometric optical scattering nuclear coefficient and a volume scattering nuclear coefficient;
step 3-2, mixing fisoExtracting the formula, transforming the earth surface bidirectional model into
Step 3-3, to kernel function coefficient fiso、fgeoAnd fvolAccording to the time sequence, statistics is carried out in a month unit in N years, and the following are respectively calculated:
f′geo＝fgeo/fiso，
f′vol＝fvol/fiso，
obtaining a geometric optical scattering nuclear coefficient f 'after deformation'geoAnd a back-body scattering nuclear coefficient f 'of deformation'volAverage value of (a).
4. The method of claim 3, wherein step 4 comprises: constructing an RPV (resilient packet virtual) surface reflection model:
where ρ is0The parameter is a hot spot parameter of a surface bidirectional reflection model, K is a kernel function, K is an anisotropic parameter of the surface bidirectional reflection model, and theta is an empirical parameter related to a scattering phase function.
5. The method of claim 4, wherein step 5 comprises:
step 5-1, selecting MODIS afternoon star Aqua surface reflectivity daily product MYD09GA, and selecting the second reflectivity of the wave band as O2Band A with surface reflectivity and band six reflectivity as CO21 the reflectance of the earth surface and the reflectance of the wave band seven as CO2-2 zone surface reflectivity;
step 5-2, counting the data from 2006 to 2016 for 10 years according to longitude, latitude and season and making a lookup table; and (3) taking the longitude and latitude of the center of the observation point, searching surrounding 10 multiplied by 10 grid points, and taking the average of all effective values as the prior value of the earth surface reflectivity of the point.
6. The method of claim 5, wherein step 6 comprises:
step 6-1, selecting an ERA-Interim data set of the ECMWF, setting the coordinates of the center longitude and latitude of a projection point as (x, y) by adopting a bilinear interpolation algorithm, setting the coordinates of the nearest four grid points around as (x1, y1), (x2, y2), (x3, y3) and (x4, y4), and setting the final profile as:
P＝(x1-x)(y1-y)p1+(x2-x)(y2-y)p2+(x3-x)(y3-y)p3+(x4-x)(y4-y)p4
wherein P is a matrix of 3 columns and 20 rows, and 3 columns are respectively an observationMeasuring point temperature, humidity, pressure profile, p1，p2，p3，p4Four corner point profiles;
6-2, according to the optical thickness tau of the aerosol at 550nm550And wavelength index α, converted to optical thickness calculations for inversion of three channels:
τ760＝τ550·1.382-α，
τ1580＝τ550·2.873-α，
τ2050＝τ550·3.727-α。
7. the method of claim 6, wherein step 7 comprises:
step 7-1, CO2The profile prior values are from the CT2017 dataset, CO, of the carbon tracking plan2The profile only contains CO below 50km of atmosphere2Profile, using standard CO2Profile and standard pressure profile, for CO2The profile is expanded to the top 100km of the atmosphere layer;
step 7-2, covariance matrix XcovIs calculated as follows:
in the formula, xp,pElements representing the p row and p column of the covariance matrix, n representing the number of profile samples, xk,iRepresenting the concentration value of the ith layer of the kth profile,
representing the concentration mean value of the ith layer, wherein subscripts i and j represent the number of the profile layers, and p is the total number of the profile layers; the parameters k, i and j have a value in the range of [1,46]]。
8. The method of claim 7, wherein step 8 comprises:
step 8-1, calculating radiation transmission y:
y＝F(x,b)+ε
wherein F (x) is the simulation calculation result of the forward radiation transmission model; x is a state vector to be inverted, and comprises carbon dioxide content, aerosol optical thickness and surface parameters; b is an auxiliary parameter vector; epsilon is an error term between the forward model and the measured data;
step 8-2, calculating a cost function J (x):
wherein F (x) is the simulation calculation result of the forward radiation transmission model, SεTo measure the error covariance matrix, xaAs a priori information, SaIs a prior error covariance matrix; if the cost function J (X) is less than X, the convergence is judged, and the step 10 is executed;
if not, according to the Gaussian iterative formula, after updating the state vector x, repeatedly executing the step 8:
in the formula, x(i+1)For the state vector, K, to continue to participate in the calculation after one iterationiIs the weight matrix and gamma is the iterative correction factor.
9. The method of claim 8, wherein in step 9, the CO is updated iteratively2Profile, surface reflection parameters, aerosol optical thickness parameters.
10. The method of claim 9, wherein in step 10, the exported CO is used2And O2Conversion of the profile into atmospheric CO2Column average dry air mixing ratio XCO2The formula is as follows:
in the formula NCO2(z) is CO2Number concentration in atmosphere with high degree of stratification, NO2(z) is O2The number concentrations are highly stratified in the atmosphere.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110148447.3A CN112964666A (en) | 2021-02-03 | 2021-02-03 | Atmospheric carbon dioxide content calculation method based on earth surface bidirectional reflection |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110148447.3A CN112964666A (en) | 2021-02-03 | 2021-02-03 | Atmospheric carbon dioxide content calculation method based on earth surface bidirectional reflection |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN112964666A true CN112964666A (en) | 2021-06-15 |
Family
ID=76273953
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110148447.3A Pending CN112964666A (en) | 2021-02-03 | 2021-02-03 | Atmospheric carbon dioxide content calculation method based on earth surface bidirectional reflection |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112964666A (en) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113533241A (en) * | 2021-07-19 | 2021-10-22 | 中国科学技术大学 | High-precision inversion system for atmospheric carbon dioxide concentration based on satellite infrared hyperspectral |
| CN113836731A (en) * | 2021-09-28 | 2021-12-24 | 中国科学院空天信息创新研究院 | Method and device for constructing land surface stabilized target atmospheric layer top reflectivity model |
| CN115201074A (en) * | 2022-06-30 | 2022-10-18 | 中国科学院大气物理研究所 | Method, system, device and computer-readable storage medium for remote sensing inversion of aerosol component distribution |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106407656A (en) * | 2016-08-29 | 2017-02-15 | 中国科学院遥感与数字地球研究所 | Retrieval method for aerosol optical thickness based on high resolution satellite image data |
| CN106483050A (en) * | 2015-09-02 | 2017-03-08 | 中国科学院遥感与数字地球研究所 | The inversion method of aerosol optical depth and system |
| CN111257241A (en) * | 2020-01-20 | 2020-06-09 | 中国科学院地理科学与资源研究所 | Atmospheric carbon dioxide concentration inversion algorithm based on DEEI (DeEI) |
| CN111413296A (en) * | 2019-09-05 | 2020-07-14 | 中国科学院烟台海岸带研究所 | Aerosol optical thickness remote sensing inversion method considering surface non-Lambert characteristics |
-
2021
- 2021-02-03 CN CN202110148447.3A patent/CN112964666A/en active Pending
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106483050A (en) * | 2015-09-02 | 2017-03-08 | 中国科学院遥感与数字地球研究所 | The inversion method of aerosol optical depth and system |
| CN106407656A (en) * | 2016-08-29 | 2017-02-15 | 中国科学院遥感与数字地球研究所 | Retrieval method for aerosol optical thickness based on high resolution satellite image data |
| CN111413296A (en) * | 2019-09-05 | 2020-07-14 | 中国科学院烟台海岸带研究所 | Aerosol optical thickness remote sensing inversion method considering surface non-Lambert characteristics |
| CN111257241A (en) * | 2020-01-20 | 2020-06-09 | 中国科学院地理科学与资源研究所 | Atmospheric carbon dioxide concentration inversion algorithm based on DEEI (DeEI) |
Non-Patent Citations (1)
| Title |
|---|
| 吴浩: "大气CO2反演中云和气溶胶的影响及其校正方法", 《中国博士学位论文全文数据库 基础科学辑》 * |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113533241A (en) * | 2021-07-19 | 2021-10-22 | 中国科学技术大学 | High-precision inversion system for atmospheric carbon dioxide concentration based on satellite infrared hyperspectral |
| CN113533241B (en) * | 2021-07-19 | 2022-12-20 | 中国科学技术大学 | High-precision inversion system for atmospheric carbon dioxide concentration based on satellite infrared hyperspectral |
| CN113836731A (en) * | 2021-09-28 | 2021-12-24 | 中国科学院空天信息创新研究院 | Method and device for constructing land surface stabilized target atmospheric layer top reflectivity model |
| CN113836731B (en) * | 2021-09-28 | 2023-06-20 | 中国科学院空天信息创新研究院 | Construction method and device of land surface stable target atmosphere top reflectivity model |
| CN115201074A (en) * | 2022-06-30 | 2022-10-18 | 中国科学院大气物理研究所 | Method, system, device and computer-readable storage medium for remote sensing inversion of aerosol component distribution |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112964666A (en) | Atmospheric carbon dioxide content calculation method based on earth surface bidirectional reflection | |
| Basu et al. | The impact of transport model differences on CO 2 surface flux estimates from OCO-2 retrievals of column average CO 2 | |
| CN110750904B (en) | Regional carbon reserve space pattern monitoring system and method based on remote sensing data | |
| Kulawik et al. | Characterization of Tropospheric Emission Spectrometer (TES) CO 2 for carbon cycle science | |
| Manney et al. | Solar occultation satellite data and derived meteorological products: Sampling issues and comparisons with Aura Microwave Limb Sounder | |
| CN106776481B (en) | Downscaling correction method for satellite precipitation data | |
| Shephard et al. | TES ammonia retrieval strategy and global observations of the spatial and seasonal variability of ammonia | |
| Santer et al. | Atmospheric correction over land for MERIS | |
| CN106407656A (en) | Retrieval method for aerosol optical thickness based on high resolution satellite image data | |
| CN101403790B (en) | Accurate one-point positioning method for single-frequency GPS receiver | |
| CN100362318C (en) | Atmosphere correction method of airosol optical thickness of aeronautical high-spectrum remote-sensing inversion boundary layer | |
| Walker et al. | An effective method for the detection of trace species demonstrated using the MetOp Infrared Atmospheric Sounding Interferometer | |
| Yin et al. | Accelerating methane growth rate from 2010 to 2017: leading contributions from the tropics and East Asia | |
| Byrne et al. | Improved constraints on northern extratropical CO2 fluxes obtained by combining surface‐based and space‐based atmospheric CO2 measurements | |
| CN104181515A (en) | Shallow sea water depth inversion method based on high-spectrum data of blue-yellow wave band | |
| Xie et al. | Calculating NDVI for Landsat7-ETM data after atmospheric correction using 6S model: A case study in Zhangye city, China | |
| Li et al. | A method for estimating hourly photosynthetically active radiation (PAR) in China by combining geostationary and polar-orbiting satellite data | |
| CN106324620A (en) | Tropospheric zenith delay method based not on real-time measurement of surface meteorological data | |
| CN106680273B (en) | High-spatial-resolution satellite earth surface reflectivity inversion method | |
| Li et al. | Satellite remote sensing for estimating PM 2.5 and its components | |
| Xu et al. | A global long-term (1981–2019) daily land surface radiation budget product from AVHRR satellite data using a residual convolutional neural network | |
| CN113534194B (en) | Troposphere temperature and humidity profile inversion method combining GNSS and wind lidar | |
| Li et al. | Generating daily high-resolution and full-coverage XCO2 across China from 2015 to 2020 based on OCO-2 and CAMS data | |
| CN114415208A (en) | Foundation GNSS convection layer top detection method with additional external data set information | |
| Kulawik et al. | TES atmospheric profile retrieval characterization: An orbit of simulated observations |
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 | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20210615 |