1. Urban CO acquisition based on satellite-borne laser radar₂Method for gridding fluxThe method is characterized in that: calculating a reverse footprint of a specified three-dimensional space-time position by using a Lagrange particle diffusion model, wherein the specified three-dimensional space-time position comprises date, time, longitude, latitude and height; relying on the existing artificial CO₂Emission flux and biological CO₂Flux construction prior CO₂Gridding the flux field and its spatio-temporal covariance; establishing a quantitative relationship between the flux field and the concentration field, and obtaining optimized posterior CO by Bayesian inversion₂The flux product is reticulated.

2. The method of claim 1 for urban CO acquisition based on space-borne lidar₂A method of gridding flux, comprising: the implementation process comprises the following steps of,

step 1, obtaining CO as required₂Calculating satellite subsatellite point positions including longitude and latitude and transit time by using a high-precision satellite orbit simulator within the space range and the time period of the flux;

step 2, taking the space-time information of the satellite subsatellite points as input, and calculating the layered footprints of each subsatellite point by utilizing a WRF-driven STILT model;

step 3, utilizing the existing artificial CO₂Emission list and natural CO₂Flux production, construction of CO on a predetermined spatial scale₂A priori flux field, and a priori covariance is set;

step 4, utilizing atmospheric temperature, humidity and pressure vertical profiles provided by meteorological data and CO provided by HITRAN database₂Molecular absorption parameters, calculating CO of each satellite by adopting Voigt linear mode₂A weight function of the concentration product;

step 5, carrying out vector multiplication on the weight function obtained in the step 4 and the layered footprint obtained in the step 2 to obtain a column weighted footprint field;

step 6, utilizing the prior flux information obtained in the step 3, the column weighted footprint field and the satellite-borne CO obtained in the step 5₂-IPDA LIDAR observed CO₂Constructing a target equation for the column concentration;

step 7, obtaining the optimized preset equation by using the target equation given in the Bayesian inversion optimization step 6Spatial scale CO₂Flux posterior solution.

3. The method of claim 2 for urban CO acquisition based on space-borne lidar₂A method of gridding flux, comprising: urban CO₂The scale of the gridded flux was 0.1 ° or 10 km.

4. Urban CO acquisition based on satellite-borne laser radar₂Gridding flux system, its characterized in that: for realizing acquisition of urban CO based on satellite-borne lidar according to any of claims 1 to 3₂Gridding flux method.

5. The method of claim 4 for urban CO acquisition based on space-borne lidar₂Gridding flux system, its characterized in that: comprises the following modules which are used for realizing the functions of the system,

a first module for obtaining CO on demand₂Calculating satellite subsatellite point positions including longitude and latitude and transit time by using a high-precision satellite orbit simulator within the space range and the time period of the flux;

the second module is used for taking the space-time information of the satellite subsatellite points as input and calculating the layered footprints of each subsatellite point by utilizing a WRF-driven STILT model;

a third module for utilizing existing artificial CO₂Emission list and natural CO₂Flux production, construction of CO on a predetermined spatial scale₂A priori flux field, and a priori covariance is set;

a fourth module for utilizing the atmospheric temperature, humidity and pressure vertical profiles provided by meteorological data and the CO provided by the HITRAN database₂Molecular absorption parameters, calculating CO of each satellite by adopting Voigt linear mode₂A weight function of the concentration product;

the fifth module is used for carrying out vector multiplication on the weight function obtained by the fourth module and the layered footprint obtained by the second module to obtain a column weighted footprint field;

a sixth module to utilize the prior flux information, the column weighted footprint field, and the satellite-borne CO₂-IPDA LIDAR observed CO₂Constructing a target equation for the column concentration;

a seventh module, configured to optimize a target equation given by the sixth module using bayesian inversion to obtain an optimized preset spatial scale CO₂Flux posterior solution.

6. The method of claim 4 for urban CO acquisition based on space-borne lidar₂Gridding flux system, its characterized in that: comprising a processor and a memory, wherein the memory is used for storing program instructions, and the processor is used for calling the stored instructions in the memory to execute the method for acquiring the city CO based on the spaceborne lidar according to any one of claims 1 to 5₂Gridding flux method.

7. The method of claim 4 for urban CO acquisition based on space-borne lidar₂Gridding flux system, its characterized in that: comprising a readable storage medium having stored thereon a computer program which, when executed, implements a method for urban CO acquisition based on a lidar according to any of claims 1 to 5₂Gridding flux method.

Background

IPDA LIDAR is an abbreviation for integrated path differential absorption LIDAR, known in the chinese translation as path-integrated differential absorption LIDAR. Since 2010, satellite-borne CO was established successively in the United states, Europe and China₂The item is detected IPDA LIDAR. After 2021 the associated load will gradually move into the track and start to make continuous ground observation. This type of plant makes it possible to obtain CO worldwide₂Column weighted concentrations. The primary goal of these satellite programs was to obtain more accurate CO on an intercontinental scale of-1000 km₂And (4) flux production. However, in recent years, with the agreement on "Paris climate Agreement", the world is accelerating towards the goal of carbon neutralization. Therefore, accurate determination of urban area high-resolution (10km scale) grid CO₂Flux becomes an urgent technical problem to be solved. Existing CO suitable for satellite products₂The goal of the flux inversion framework design is to provide a 1000km scale of CO₂Flux products, even optimized to provide at most 100km scale CO₂And (4) flux production. Meanwhile, the inversion frameworks are CO based on passive satellite remote sensing technology₂Product designed for column concentration, not perfectly adapted for lidar-based CO₂Column concentration product. Therefore, the city CO suitable for the satellite-borne laser radar is provided₂The method for inverting the gridding flux is a great problem to be solved urgently in the field.

Disclosure of Invention

For obtaining small-scale (0.1 degree or 10km) urban CO₂The invention provides a grid-connected flux, and provides a method and a system for acquiring grid-connected flux of city CO2 based on a satellite-borne laser radar.

The technical scheme of the invention provides a method for acquiring urban CO based on a satellite-borne laser radar₂Method of screening flux, usingCalculating a reverse footprint of a specified three-dimensional space-time position by using a Lagrange particle diffusion model, wherein the specified three-dimensional space-time position comprises date, time, longitude, latitude and height; relying on the existing artificial CO₂Emission flux and biological CO₂Flux construction prior CO₂Gridding the flux field and its spatio-temporal covariance; establishing a quantitative relationship between the flux field and the concentration field, and obtaining optimized posterior CO by Bayesian inversion₂The flux product is reticulated.

Moreover, the implementation process comprises the following steps,

step 1, obtaining CO as required₂Calculating satellite subsatellite point positions including longitude and latitude and transit time by using a high-precision satellite orbit simulator within the space range and the time period of the flux;

step 2, taking the space-time information of the satellite subsatellite points as input, and calculating the layered footprints of each subsatellite point by utilizing a WRF-driven STILT model;

step 3, utilizing the existing artificial CO₂Emission list and natural CO₂Flux production, construction of CO on a predetermined spatial scale₂A priori flux field, and a priori covariance is set;

step 4, utilizing atmospheric temperature, humidity and pressure vertical profiles provided by meteorological data and CO provided by HITRAN database₂Molecular absorption parameters, calculating CO of each satellite by adopting Voigt linear mode₂A weight function of the concentration product;

step 5, carrying out vector multiplication on the weight function obtained in the step 4 and the layered footprint obtained in the step 2 to obtain a column weighted footprint field;

step 6, utilizing the prior flux information obtained in the step 3, the column weighted footprint field and the satellite-borne CO obtained in the step 5₂-IPDA LIDAR observed CO₂Constructing a target equation for the column concentration;

step 7, obtaining the optimized preset spatial scale CO by using the target equation given in the Bayesian inversion optimization step 6₂Flux posterior solution.

Furthermore, urban CO₂The scale of the gridded flux was 0.1 ° or 10 km.

The invention provides a method for acquiring urban CO based on a satellite-borne laser radar₂The grid-based flux system is used for realizing the acquisition of urban CO based on the satellite-borne laser radar₂Gridding flux method.

And, including the following modules,

a first module for obtaining CO on demand₂Calculating satellite subsatellite point positions including longitude and latitude and transit time by using a high-precision satellite orbit simulator within the space range and the time period of the flux;

the second module is used for taking the space-time information of the satellite subsatellite points as input and calculating the layered footprints of each subsatellite point by utilizing a WRF-driven STILT model;

a third module for utilizing existing artificial CO₂Emission list and natural CO₂Flux production, construction of CO on a predetermined spatial scale₂A priori flux field, and a priori covariance is set;

a fourth module for utilizing the atmospheric temperature, humidity and pressure vertical profiles provided by meteorological data and the CO provided by the HITRAN database₂Molecular absorption parameters, calculating CO of each satellite by adopting Voigt linear mode₂A weight function of the concentration product;

the fifth module is used for carrying out vector multiplication on the weight function obtained by the fourth module and the layered footprint obtained by the second module to obtain a column weighted footprint field;

a sixth module to utilize the prior flux information, the column weighted footprint field, and the satellite-borne CO₂-IPDA LIDAR observed CO₂Constructing a target equation for the column concentration;

a seventh module, configured to optimize a target equation given by the sixth module using bayesian inversion to obtain an optimized preset spatial scale CO₂Flux posterior solution.

Or, the system comprises a processor and a memory, wherein the memory is used for storing program instructions, and the processor is used for calling the stored instructions in the memory to execute the method for acquiring the city CO based on the spaceborne lidar₂Gridding flux method.

Or, includeRead storage medium, on which a computer program is stored, which, when executed, implements a method for acquiring urban CO based on a space-borne lidar as described above₂Gridding flux method.

The invention is simultaneously suitable for satellite-borne CH₄Detecting laser radar CH₄And carrying out grid flux inversion on the ground laser radar greenhouse gas concentration product.

Drawings

FIG. 1 is a flow chart of an embodiment of the present invention;

FIG. 2 is a schematic diagram of an example of the height distribution of the output atmospheric weighting function according to the embodiment of the present invention.

Detailed Description

For obtaining small-scale (0.1 degree or 10km) urban CO₂The invention provides a method for calculating reverse footprints of specified three-dimensional space-time positions (date, time, longitude, latitude and height) by utilizing a Lagrange particle diffusion model and further relying on the existing artificial CO₂Emission flux and biological CO₂Flux (NEE) construction of a priori CO₂The method comprises the steps of establishing a quantitative relation between a flux field and a concentration field by means of the two data sets and obtaining optimized posterior CO by means of Bayesian inversion₂The flux product is reticulated.

Referring to fig. 1, an embodiment of the present invention provides a method for obtaining grid-like flux of city CO2 based on a satellite-borne laser radar, including the following steps:

step S1, accurate settlement of satellite subsatellite point positions: obtaining CO on demand₂And calculating the satellite subsatellite point position including longitude and latitude and transit time by using a high-precision satellite orbit simulator within the space range and the time period of the flux.

In the embodiment, satellite orbit information is calculated by using an ephem library based on python language and two-row satellite orbit parameters defined by the north american air defense commander, and satellite intersatellite point coordinate information in a specified region and a time range is obtained by using self-defined space-time range information. For ease of reference, specific implementations are provided as follows:

first, a readtle function using an ephem library and a loaded CO₂-IPDA LIDAR, completing initialization.

Secondly, setting the starting time by using the datatime data type of python, and simultaneously limiting the spatial range to be calculated by using the longitude and latitude.

Thirdly, a computer function of the ephem library is driven by a time vector established by the starting time and the duration time to obtain all the latitude and longitude of the sub-satellite points in the specified time range, and whether the latitude and longitude of the sub-satellite points are located in the specified space range is judged.

Fourthly, outputting the satellite subsatellite point space-time information which meets the requirement, wherein the specific form of the satellite subsatellite point space-time information is a data table with the attributes including date, time, longitude and latitude.

And finally, according to the longitude and latitude information, further adding an elevation attribute to each piece of data by using a high-precision digital elevation model. And finally obtaining satellite intersatellite point information in the appointed space-time range.

Step S2, obtaining the layered footprint field.

This step preferably uses the WRF-STILT model to simulate the backward random trajectory of the particle at a given location (longitude, latitude, altitude) under the control of the weather field and to obtain the footprint field of that point. The WRF-STILT model is a Lagrange atmospheric propagation model. The embodiment takes the space-time information of the satellite subsatellite points as input, and calculates the 0-15km layered footprint of each subsatellite point by using a WRF-driven STILT model.

The aim of the invention is to obtain 0.1 ° (-10 km) resolved flux information, so the resolution of the footprint field needs to be set to 0.1 °. Thus, the meteorological field resolution needs to be better than 0.1 °. The WRF settings required the use of three-tier nesting, with either (27km, 9km, 3km) or (32km, 8km, 2km) embodiments. The time range of the WRF simulation needs to be reversed according to the time of the sub-satellite point, the recommended reversing time is 7 days, the shortest time is not less than 3 days, and the longest time is not more than 15 days. Other settings for the WRF (atmospheric scheme, radiation scheme, terrain scheme, etc.) may be determined based on the conditions and time of the designated area. And converting the format of the weather field output by the WRF, inputting the weather field into the STILT model, and calculating the footprint field of the specified position.

The STILT model is a transport model of Lagrangian random walk theory that relates source (sink) flux upstream of an observation point to concentration changes at the observation point with footprint weights. STILT can accurately simulate ultrahigh resolution in urban areas<1km) capacity of atmospheric transmission process, so the model fully satisfies urban CO₂Flux inversion requires a 0.1 ° resolution footprint. The specific principle is that a backward movement track of gas under the drive of turbulent flow and average wind direction is simulated by releasing a large number of air particles backward, and a footprint weight value f () is calculated quantitatively by calculating the number of all particles in a certain height of a boundary layer of a certain upstream area and the residence time of each particle, wherein the calculation formula is as follows:

wherein x is_rIs the position and height of the observation point, t_rFor corresponding observation times, m_airIs the molar weight of air, h is the height of the underlying influence layer, ρ () is the average density of all particles, t_mRepresents a discrete time step, (x)_i,y_j) Representing the underlying surface of a region, p represents a particle, Δ t_p,i,jRepresenting each particle corresponding to a certain region (x)_i,y_j) The underlying surface of (a) affects the time the layer resides. N is a radical of_totIs the total number of particles released.

Since the satellite observation object is the whole layer CO from the earth surface to the top of the atmospheric layer₂And (4) calculating the layered footprints of different high layers of each observation point theoretically according to the column concentration. This process is extremely computationally intensive, and this step thus gives a preferred hierarchical scheme obtained over multiple experiments. The method specifically comprises the following steps:

(1) layering is carried out within the range of 0-0.8 km at the speed of 0.2km, and 4 layered footprints are calculated;

(2) layering is carried out within the range of 0.8-2.0 km at 0.4km, and footprint fields of 3 layers are calculated;

(3) layering at 1km within the range of 2.0-5.0 km, and calculating footprint fields of 3 layers;

(4) and calculating footprint fields at 10km and 12km within the range of 5-12 km.

(5) And calculating no more than 12km, and finally obtaining 12 layers of layered footprints for each observation point.

And step S3, constructing a priori flux field and covariance.

Example utilizing existing artificial CO₂Emission list and natural CO₂Flux (NEE) product, CO construction on a spatial scale of 0.1 ° (-10 km)₂A priori flux field and a priori covariance is set.

The prior flux field is composed of artificial and natural CO₂Flux field composition. Anthropogenic CO in non-U.S. regions₂Construction of the flux field recommends the use of the EDGAR version 5.0 plus emissions list, and the U.S. region recommends the use of VULCAN. Natural CO₂The flux field recommends using a NEE product set based on SIF (solar induced fluorescence), and recommends using an open-source SMUrF code to calculate the high-resolution NEE of a specified region by itself.

The prior covariance matrix B (m × m dimension) can be represented by:

wherein D (dimension D × D) is a temporal covariance matrix and E (dimension E × E) is a spatial covariance matrix,is the Kronecker product, so n ═ d × e.

The temporal covariance matrix and the spatial covariance matrix may be represented by:

D＝V_t ^1/2M_tV_t ^1/2 3

E＝V_s ^1/2M_sV_s ^1/2 4

wherein M is_tAnd M_sRespectively a temporal correlation matrix and a spatial correlation matrix, V_tAnd V_sRespectively, a diagonal matrix of temporal and spatial variances.

Step S4, a weight function is calculated.

Embodiments utilize atmospheric temperature, humidity, and pressure vertical profiles provided by meteorological data and CO provided by a HITRAN database₂Molecular absorption parameters, calculating CO of each satellite by adopting Voigt linear mode₂Weight function of concentration product.

Equation 5 gives the inversion of atmospheric CO using IPDA (path-integrated differential absorption)₂Theoretical equation of column concentration, where P_λAnd E_λFor the echo signal intensity and emission energy at a wavelength λ, ω (p) is CO₂A concentration profile. The denominator part is the atmospheric weight function (see equation 6). In formula 6, σ is a molecular absorption cross-sectional area, M is a molecular mass, ρ is a density, p is an atmospheric pressure, and g is a gravitational acceleration. IPDA LIDAR the 2 operating wavelengths are commonly denoted as λ_on、λ_offSo that λ is λ ═ λ_onOr λ_off。

As can be seen from equation 6, the variables in WF include σ and ρ_H2O. It should also be noted that although generally expressed in terms of σ (p, λ) for atmospheric detection, σ () is actually a variable determined collectively by temperature, pressure, and wavelength after specifying a particular gas molecule. Specific expressions are given by equations 7-11, wherein, v is the wave number (reciprocal of wavelength), v_cIs the central wavenumber of the absorption spectrum, S (T) is the absorption line intensity, t is the integral variable, γ_LIs Lorentz half-width, gamma_DIs the doppler half-width. p is a certain atmospheric pressure, p₀And T₀Standard atmospheric pressure and standard temperature, 101325Pa and 296.15K, respectively. Gamma ray₀Is a standard lineWidth, n is the temperature dependent coefficient of the pressure broadening coefficient, k_BIs the boltzmann constant, E "is the low energy state, h is the planck constant, and c is the speed of light.

ν_c＝ν₀+p×(δ⁰(T₀)+δ'(T-T₀)) 10

Wherein, v₀Is atmospheric pressure p-0 Pa and temperature T-T₀Wave number of absorption peak of time, delta⁰(T₀) Is T ═ T₀Pressure translation coefficient of time, delta' (T-T)₀) Is the temperature dependent coefficient of the pressure translation coefficient, T₀Is the standard temperature; v is_cIs the central wavenumber of the absorption spectrum;

m is the molecular mass, c is the speed of light, k_BIs the Boltzmann constant, γ_DIs the Doppler half-width;

γ_Lis Lorentz half-width, gamma₀Is the standard line width, p₀Is at standard atmospheric pressure, n isA temperature dependent coefficient of pressure broadening coefficient;

s (T) is the absorption line intensity, S₀Is the absorption line intensity at a temperature of 296.13K and a pressure of 101325Pa, E' is a low energy state, h is a Planckian constant, exp () is an exponential function;

σ (p, T, v) is the molecular absorption cross-sectional area at temperature T, pressure p and wave number v, ν is the wavenumber (reciprocal of wavelength);

WF (p) is a weight function describing the absorption capacity of light by specific atmospheric molecules at different pressures, σ (p, λ)_on) Is the pressure p, the wavelength lambda_onLower molecular absorption cross section, σ (p, λ)_off) Is the pressure p, the wavelength lambda_offThe absorption cross-sectional area of the lower molecule,is the mass of water molecule, M_airIs the molecular weight of the air and is,is the water vapor density, g is the acceleration of gravity;

XCO₂is column average CO₂Dry air mixture ratio, Psurf surface pressure, ω (p) CO₂Concentration profile, P_onIs the echo signal intensity at a wavelength λ, P_offIs the echo signal intensity at a wavelength λ, E_onIs a wavelength of λ_onIntensity of time-dependent echo signal, E_offIs a wavelength of λ_offThe echo signal strength of time.

Fig. 2 shows the WF distribution calculated by S4, which is calculated according to the current atmospheric environment star parameters and can be directly used. If other similar laser CO is used in the future₂The satellites are detected and need to be recalculated according to the step of S4.

Step S5, generation of a column weighted footprint field.

And in the embodiment, the weight function obtained in the step 4 and the layered footprint obtained in the step 2 are subjected to vector multiplication to obtain a column weighted footprint field.

And multiplying the layered footprints obtained in the step S2 by the WF obtained in the step S4 according to the corresponding elevations, and summing the N products to obtain the column weighted footprint field. This process can be expressed by equation 12. z represents different layering numbers, and a 12-layer (N-12) layering scheme is recommended in step 2. In addition, other layering schemes can be designed according to requirements.

Wherein, fotopprint (z) represents the surface flux footprint at elevation z.

And step S6, constructing an objective equation of the relation among the flux, the footprint and the concentration.

In the embodiment, the prior flux information obtained in the step 3, the column weighted footprint field obtained in the step 5 and the satellite-borne CO are utilized₂-IPDA LIDAR observed CO₂Column concentrations build the target equation.

x^*Represents CO₂A priori field of flux with CO₂The relationship of concentration can be given by:

y^*＝Hx^*+ε. 12

wherein, y^*Is XCO₂Is the mean value of_bAnd the diagonal covariance matrix is the Nx 1 vector of R, which is obeyed from ε to N (μ)_bAnd R) normally distributed noise. R is determined by the nature of the observation device itself and the observation environment, and in practice can be characterized by the fluctuation level of the observation result in a specified time range. H is footprint matrix output by WRF-STILT and has physical meaning of carbon source sink at different positions to CO at specified site₂Contribution of concentration values. The mismatch errors are assumed to be uncorrelated and the average deviation is zero. These synthetic observations can be used in a bayesian inference framework to estimate the best CO₂Flux.

Step S7, Bayesian inversion is conducted on the optimal solution of the flux field.

In the embodiment, the target equation given in the step 6 is optimized by using a Bayesian inversion technology to obtain the optimized 0.1-degree resolution CO₂Flux posterior solution.

The posterior CO can be given assuming that the probability distribution of the prior error obeys a Gaussian distribution₂Closed form solution of flux:

wherein x is_pIs a priori CO₂The flux, B, is the m x m prior error covariance matrix, given by step S3.For the high resolution meshed CO finally obtained by the invention₂And (4) flux production.

In specific implementation, a person skilled in the art can implement the automatic operation process by using a computer software technology, and a system device for implementing the method, such as a computer-readable storage medium storing a corresponding computer program according to the technical solution of the present invention and a computer device including a corresponding computer program for operating the computer program, should also be within the scope of the present invention.

In some possible embodiments, a system for acquiring urban CO2 grid-based traffic based on satellite-borne lidar is provided, which includes the following modules,

a first module for obtaining CO on demand₂Calculating satellite subsatellite point positions including longitude and latitude and transit time by using a high-precision satellite orbit simulator within the space range and the time period of the flux;

the second module is used for taking the space-time information of the satellite subsatellite points as input and calculating the layered footprints of each subsatellite point by utilizing a WRF-driven STILT model;

a third module for utilizing existing artificial CO₂Emission list and natural CO₂Flux production, construction of CO on a predetermined spatial scale₂A priori flux field and a priori co-square is setA difference;

a fourth module for utilizing the atmospheric temperature, humidity and pressure vertical profiles provided by meteorological data and the CO provided by the HITRAN database₂Molecular absorption parameters, calculating CO of each satellite by adopting Voigt linear mode₂A weight function of the concentration product;

the fifth module is used for carrying out vector multiplication on the weight function obtained by the fourth module and the layered footprint obtained by the second module to obtain a column weighted footprint field;

a sixth module to utilize the prior flux information, the column weighted footprint field, and the satellite-borne CO₂-IPDA LIDAR observed CO₂Constructing a target equation for the column concentration;

a seventh module, configured to optimize a target equation given by the sixth module using bayesian inversion to obtain an optimized preset spatial scale CO₂Flux posterior solution.

In some possible embodiments, there is provided a system for acquiring urban CO2 gridded flux based on-board lidar, comprising a processor and a memory, the memory storing program instructions, the processor being configured to invoke the stored instructions in the memory to perform a method for acquiring urban CO2 gridded flux based on-board lidar as described above.

In some possible embodiments, a system for acquiring urban CO2 gridded traffic based on satellite-borne lidar is provided, which includes a readable storage medium on which is stored a computer program that, when executed, implements a method for acquiring urban CO2 gridded traffic based on satellite-borne lidar as described above.

The specific embodiments described herein are merely illustrative of the spirit of the invention. Various modifications or additions may be made to the described embodiments or alternatives may be employed by those skilled in the art without departing from the spirit or ambit of the invention as defined in the appended claims.