Ozone profile inversion method based on thermal infrared remote sensing

Abstract

The invention provides an ozone profile inversion method based on thermal infrared remote sensing, and belongs to the technical field of atmospheric environment remote sensing. The technical problem of low precision of ozone profile inversion by using thermal infrared remote sensing is solved. According to the technical scheme, the method comprises the following steps: S1, collecting and processing data; s2, establishing an ozone prior profile model; s3, establishing a temperature and humidity profile error correction model; s4, performing thermal infrared ozone profile optimization inversion; and S5, determining a scale factor. The method has the advantages that the total stratosphere column amount and the total troposphere column amount are used for constraint in the optimization inversion, so that the precision of the ozone profile is improved, especially the precision of the middle layer and the lower layer of the troposphere is improved, and the defects of a thermal infrared sensor in a nadir observation mode are overcome.

Description

Technical Field

[0001] The invention relates to the technical field of atmospheric environment remote sensing, and in particular to an ozone profile inversion method based on thermal infrared remote sensing. Background Art

[0002] In the thermal infrared channel, since the main signal source comes from temperature and there is continuous absorption of water vapor, the error of temperature and humidity profile during the inversion process will affect the weight function of radiation and ozone, and thus affect the inversion accuracy of ozone profile. Figure 2 It is based on the sensitivity analysis of the HIRAS channel. When the temperature profile is +1K, the humidity profile is Q-10%, and the O3 concentration is -10%, the radiance at the HIRAS entrance pupil changes. Figure 2 It can be seen that the ozone absorption band is interfered by the temperature and humidity profile.

[0003] Figure 2 Sensitivity analysis of the HIRAS channel located near 9.6μm. The radiance change at the HIRAS entrance pupil when the temperature profile is +1K, the humidity profile is Q-10%, and the O3 concentration is -10% (based on the ECMWF profile sample library).

[0004] Thermal infrared uses brightness temperature difference to invert gas concentration. In the lower troposphere, the difference in brightness temperature between the atmosphere and the ground is very small, so the radiation information available for inversion is very limited. Secondly, the ultraviolet scattering spectrum is affected by factors such as absorption line band dependence and surface reflectivity. The absorption of ultraviolet radiation by tropospheric aerosols often leads to an underestimate of the total amount of ozone. The error caused by absorbing dust aerosols is greater than the inversion error caused by non-absorbing aerosols, and the vertical distribution information of ultraviolet ozone is limited. Simple joint inversion cannot solve the inversion uncertainty caused by the ultraviolet surface reflectivity and absorption line band dependence and the infrared temperature error. This inversion model has a partial effect on improving the information freedom of tropospheric ozone, but it does not improve the inversion accuracy. Since 10% of ozone is concentrated in the troposphere, satellite-borne sensors that observe in a nadir mode to detect ozone are only sensitive to the entire layer of ozone and ozone above the upper troposphere, whether in the ultraviolet band or the infrared band, which leads to a very large uncertainty in the inversion of tropospheric ozone. Summary of the invention

[0005] The purpose of the present invention is to provide an ozone profile inversion method based on thermal infrared remote sensing. This part is mainly based on thermal infrared sensor data, and the optimization method is used to invert the ozone profile based on the total stratospheric column and the total tropospheric column, respectively. The optimization model is used to improve the accuracy of the ozone profile in the nadir observation mode, especially to improve the accuracy of the middle and lower troposphere.

[0006] In order to achieve the above-mentioned invention object, the technical solution adopted by the present invention is specifically: an ozone profile inversion method based on thermal infrared remote sensing, comprising the following steps:

[0007] Step (1): Data collection and processing: Obtain satellite thermal infrared spectral data (such as FY3D / E / F-HIRAS, NOAA20-CrIS), cloud product data, temperature and humidity profiles of reanalysis ERA5 data, tropopause data, total ozone column data and total ozone column data of tropospheric ozone (such as Sentinel5P-TROPOMI), remove clouds from thermal infrared spectral data using cloud product data (MERSI cloud product), and perform spatiotemporal matching of reanalysis data, tropopause data, total ozone column data and total tropospheric column data based on thermal infrared spectral data;

[0008] Thermal infrared hyperspectral data: The infrared hyperspectral detection instrument HIRAS is carried on the FY-3D (launched in 2017), FY-3E (launched in 2021) and FY-3F (planned to be launched in 2022) satellites. The HIRAS carried on FY3D is the first infrared hyperspectral detection instrument on my country's polar-orbiting meteorological satellite. It has three wavelength ranges: long-wave infrared (8.80-15.39μm), medium-wave infrared (5.71-8.26μm) and short-wave infrared (3.92-4.64μm). The spectral resolution and spatial resolution are 0.625cm-1 and 16km respectively. At present, the primary spectral data of FY3D-HIRAS can be downloaded from the National Satellite Meteorological Center (http: / / www.nsmc.org.cn / NSMC / Home / Index.html). The Moderate Resolution Spectral Imager (MERSI) achieves high-precision quantitative inversion of atmospheric, land, and ocean parameters such as clouds, water vapor, and land surface characteristics through continuous and comprehensive multi-spectral observations in the low-light and infrared bands.

[0009] Reanalysis data: ERA5 is a new generation of reanalysis data, created by the Copernicus Climate Change Service (C3S) operated by the European Centre for Medium-Range Weather Forecasts (ECMWF). It has a spatial resolution of 62 kilometers and a temporal resolution of 1 hour, and can provide users with relatively accurate real-time atmospheric information. At present, the Integrated Global Radiosonde Archive (IGRA) is distributed in various regions of the world, with the number of stations as high as 1,500. There are about 92 IGRA stations in my country. The meteorological parameters that this detection can provide mainly include temperature, dew point temperature, air pressure, geopotential height, tropopause, ozone column concentration and other data, covering different pressure ranges from near the ground to the lower stratosphere.

[0010] TROPOMI (TROPOspheric Monitoring Instrument) is a sensor carried on the Sentinel-5 Precursor (S5P). It was successfully launched on October 13, 2017 and is expected to have an on-orbit life of 7 years. The payload has global daily coverage detection capabilities in the ultraviolet-visible-near infrared band. The spectral resolution of the ultraviolet band is 0.5nm. TROPOMI can provide high spatial resolution tropospheric trace gases, such as sampling column concentration, tropospheric ozone column concentration, NO2, SO2, HCHO, and aerosol information.

[0011] Step (2): Ozone prior profile model: The ozone profile is affected by a variety of factors, so its changes are subject to great uncertainty. In this study, the LSTM model is used to generate the ozone prior profile. LSTM is known as the long short-term memory network, which is the inheritance, development and optimization of the recurrent neural network (RNN). The LSTM network improves the shortcomings of the RNN network, such as gradient vanishing and gradient explosion, and is suitable for solving problems with long time series. At the same time, the LSTM network also solves the problem of long dependency of the RNN. By adding an adaptive "forget gate" to the network, the LSTM model can learn automatically and determine whether previous information is stored. This project uses this feature of LSTM to learn and model complex nonlinear relationships, build relationships between factors, and establish a priori estimation model for ozone profiles: first, collect ERA5 meteorological data (temperature and humidity, potential vorticity, wind speed, wind direction, etc.), infrared and ultraviolet hyperspectral data near the ozone absorption band (using thermal infrared and ultraviolet hyperspectral data with long time series (such as AIRS and OMI on Aqua), ozone precursors and industrial production data, ozone data from lidar, soundings, ground stations, and ERA5 reanalysis data; secondly, unify different data sets in space and time and resample them to the same vertical grid; use the random forest method to analyze the feature importance of the collected data to obtain important features of the model input; finally, use LSTM to train and test the data, select the appropriate hidden layer and number of nodes, and obtain the optimal network structure, thereby obtaining the best generalized ozone profile prior model.

[0012] Based on ground-based site data, sounding data, meteorological data, ozone precursor data, etc., the atmospheric background profile library is updated, and the ozone prior profile and ozone prior covariance matrix are obtained using a deep learning model. The covariance matrix B represents the background field covariance matrix or the error covariance of the initial profile. The background covariance matrix is ​​obtained by statistics of the prior profile library:

[0013]

[0014] Where X is the profile in the profile library, is the statistical mean of the profile library.

[0015] Step (3): Temperature and humidity profile error correction model: Use the ground-based temperature and humidity profile data and the ERA5 data of the input radiation transfer model to correct the temperature and humidity profile to obtain the correction factor: According to the climate characteristics, the input ERA5 temperature and humidity profile products and IGRA ground-based data are statistically analyzed by region to obtain the temperature and humidity profile errors at different heights, and then based on the forward radiation transfer model RTTOV, the impact of the temperature and humidity profile error on radiation is studied, and appropriate parameters are selected for correction to reduce the uncertainty of the temperature and humidity profile error.

[0016] Step (4): Optimization inversion of thermal infrared ozone profile: The ozone absorption band of 9.6 μm in the thermal infrared absorption band is used for inversion. The total ozone column and tropospheric column data are introduced to constrain the tropospheric and stratospheric ozone inversion results during the inversion process, and the Levenberg-Marquardt (LM) model is used to invert the tropospheric ozone profile. The cost function with the addition of the total ozone column is redefined as:

[0017]

[0018] The vector X = X(O3) is the n-layer ozone concentration profile, the Y vector is the brightness temperature observed by the m inversion channel satellite, γ is the introduced Lagrange multiplier, λ and ε are proportional coefficients, TO3 and TO tropo are the total ozone column and the total tropospheric ozone column, respectively. Tropo represents the number of layers of tropospheric pressure. F(X) is the atmospheric forward radiation transmission mode. The norm here is the Mahalanobis distance, so the cost function can be written as:

[0019]

[0020] R is for Y m The error covariance matrix of , which mainly includes the errors caused by the instrument and the system; B represents the initial error covariance matrix of the desired variable O3 profile. After defining the cost function, the iterative formula is derived and the constraints on the iterative formula are strengthened according to the Levenberg-Marquardt model. The positive definite matrix and the γ factor in the Levenberg-Marquardt method are adaptively adjusted through the diagonal elements and nonlinearity of the Hessian matrix;

[0021] According to the principle of steepest descent, setting J'(X) = 0, we can get the linear solution form of this problem, namely:

[0022]

[0023] Where K is the first-order derivative of Y with respect to X, also known as the Jacobian matrix. Due to the nonlinearity of the inversion problem, the Newton nonlinear iteration can be used to improve the inversion accuracy and convergence speed of the problem, that is, the cost function is derived again to obtain:

[0024] 0.5J”(X)=γB -1 +K T R -1 K (5)

[0025] J”(X) is called the Hessian matrix. Substituting it into formula (4) yields:

[0026]

[0027] or:

[0028]

[0029] This method requires a more accurate estimate of the covariance of the background error and the setting of an appropriate Lagrange multiplier γ. The Lagrange multiplier is affected to a certain extent by human factors and the variance of the inversion error.

[0030] In order to solve the two shortcomings of the Gauss-Newton iteration method: first, the Hessian matrix in the algorithm may be singular, which will cause the successor point to be unable to be determined; second, the non-singular Hessian matrix may also be non-positive definite, resulting in the invalidity of the model. There are many methods to solve the problem of non-positive definite Hessian matrix. The basic idea is to correct the matrix to construct a symmetric positive definite matrix. According to the method proposed by Levenberg-Marquardt to solve the least squares problem, a positive definite matrix is ​​constructed. Let

[0031]

[0032] Where H is a positive diagonal matrix, and w is re-determined during each iteration so that the cost function can decrease rapidly. When w is large, it is the direction of the fastest descent, and when w is small, it is the direction of Gauss-Newton descent. Therefore, the Levenberg-Marquardt (LM) method is the fastest descent method when the iteration point is far from the minimum value; when it is close to the minimum value, it is the Gauss-Newton method. These two methods complement each other, the former converges quickly when it is far from the minimum value, and the latter converges quickly when it is close to the minimum value. In the iterative process, w is the factor that reconciles the two methods. After defining the cost function in this paper, the iterative formula can be derived through the Gauss-Newton nonlinear iteration method and the iterative formula can be constrained according to the Levenberg-Marquardt model. In addition to strengthening the constraints on the iterative formula by the LM model, this section constrains the iterative results with the total stratospheric column and the total tropospheric column, respectively.

[0033] Step (5): Determination of scaling factors: The scaling factors λ and ε are determined experimentally. The existing ground-based, sounding and lidar ozone profile data are used as a benchmark. The root mean square error and root mean square relative error of the ozone profiles inverted in different regions with different scaling factors and the benchmark values ​​are analyzed to obtain reasonable allocation of scaling factors λ and ε.

[0034] Compared with the prior art, the present invention has the following beneficial effects:

[0035] 1. The present invention utilizes the total amount of stratospheric column and the total amount of tropospheric column as constraints in the optimization inversion, thereby improving the accuracy of the ozone profile, especially the accuracy of the middle and lower troposphere, thereby compensating for the shortcomings of the thermal infrared sensor in the nadir observation mode.

[0036] 2. The present invention uses the total amount of tropospheric column and the total amount of stratospheric column as constraints in the optimization inversion, which can effectively improve the accuracy and reliability of thermal infrared ozone profile inversion. At the same time, it is insensitive to ozone in the middle and lower layers of the troposphere in the nadir observation mode, reducing the uncertainty of tropospheric ozone inversion.

[0037] 3. In the traditional thermal infrared ozone profile inversion, the information of the ozone profile in the middle and lower troposphere is often derived from prior knowledge, and the inversion accuracy is low. However, after the total amount of tropospheric column is constrained, the information freedom of the troposphere is increased, making the inversion results of the middle and lower troposphere closer to the actual situation. After the stratosphere is constrained, the accuracy of the inversion results and the degree of freedom of information of the stratosphere will be greatly improved. Therefore, by combining the total amount of stratospheric and tropospheric column constraints with thermal infrared observations, the accuracy of ozone profile inversion can be improved.

[0038] 4. In practical applications, the thermal infrared ozone profile obtained by constraining the total amount of stratospheric and tropospheric columns can provide more accurate information on the vertical distribution of ozone, providing important data support for atmospheric chemistry and climate research. Accurate ozone profiles are of vital importance for understanding atmospheric circulation, climate change, and the impact of human activities on the atmospheric environment. Secondly, they are of great significance for the horizontal and vertical transport of tropospheric ozone between regions, and provide high-precision and large-scale ozone profile information for PM2.5 and ozone pollution prevention and control. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] The accompanying drawings are used to provide further understanding of the present invention and constitute a part of the specification. They are used to explain the present invention together with the embodiments of the present invention and do not constitute a limitation of the present invention.

[0040] Figure 1 It is a schematic diagram of the detailed inversion results of the Legionowo station in the mid-latitude region of the Northern Hemisphere on July 4, 2019 in an embodiment of the present invention;

[0041] Among them, (a) is a schematic diagram of the inversion results of the Legionowo site with the total column amount constraint; (b) is a schematic diagram of the inversion results of the Legionowo site without constraints.

[0042] Figure 2 A schematic flow chart of an ozone profile inversion method based on thermal infrared remote sensing provided by the present invention.

[0043] Figure 3 It is a schematic diagram of sensitivity analysis based on HIRAS channel in the prior art recorded in the background technology of the present invention. DETAILED DESCRIPTION

[0044] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. Of course, the specific embodiments described here are only used to explain the present invention and are not used to limit the present invention.

[0045] Example 1

[0046] See also Figures 2 to 3 This embodiment provides a technical solution, which is a method for inverting ozone profiles based on thermal infrared remote sensing. This embodiment uses the results of the total amount of tropospheric and stratospheric columns constrained by the method without constraints. The comparison results are as follows: Figure 1 As shown: Figure 1 Inversion results of the Legionowo station with (a) and without (b) constraints on the total column volume; comparison of sounding data, smoothed sounding data, prior profiles and inversion results (left), relative errors of prior and inversion results with sounding smoothing results (middle), and average kernel function of inversion results at different pressures (right).

[0047] Figure 1 This is a detailed description of the inversion results for the selected days at the Legionowo station in the mid-latitudes of the Northern Hemisphere. Before the Legionowo station was unconstrained, when the pressure was >500hPa, most of the inversion results came from the prior profile, and after the column total constraint, the accuracy of the inversion results and smoothing results was improved by about 15% compared with before the constraint, but the inversion results of the 200-300hPa column total constraint were slightly lower than before the constraint, about 12%, and the maximum peak of the average kernel function of 200-390hPa after the column total was about 0.98, and the maximum peak of the average kernel function before the constraint was about 1.2, and the maximum peak of the average kernel function from 103hPa to near the ground was 200-300hPa, the maximum peak after the column total constraint was 2.5, and the maximum peak without constraint was 3, but the kernel function was improved after the column total constraint at other pressures, which is because of the proportional factor after the column total constraint, which improves the inversion accuracy of the entire tropospheric ozone profile.

[0048] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for inverting ozone profiles based on thermal infrared remote sensing, characterized in that: The following steps are involved: S1. Data collection and processing: Obtain satellite thermal infrared spectral data, FY3D / E / F-HIRAS, NOAA20-CrIS, cloud product data, temperature and humidity profiles of reanalysis ERA5 data, tropospheric top data, total ozone column data and total ozone column data of tropospheric ozone, remove clouds from thermal infrared spectral data using cloud product data, and perform spatiotemporal matching of reanalysis data, tropospheric top data, total ozone column data and total tropospheric column data based on thermal infrared spectral data; S2. Ozone prior profile model: Based on ground-based site data, sounding data, meteorological data, and ozone precursor data, the atmospheric background profile library is updated. The ozone prior profile and ozone prior covariance matrix are obtained using a deep learning model. The covariance matrix B represents the background field covariance matrix or the error covariance of the initial profile. The background covariance matrix is ​​obtained from the prior profile library statistics: Where X is the profile in the profile library, is the statistical mean of the profile library; S3. Temperature and humidity profile error correction model: Use the ground-based temperature and humidity profile data and the ERA5 data of the input radiation transfer model to correct the temperature and humidity profile to obtain the correction factor: According to the climate characteristics, the input ERA5 temperature and humidity profile products and IGRA ground-based data are statistically analyzed by region to obtain the temperature and humidity profile errors at different heights. Then, based on the forward radiation transfer model RTTOV, the influence of the temperature and humidity profile error on radiation is studied, and the radiation amount is corrected to reduce the uncertainty of the temperature and humidity profile error; S4. Optimization inversion of thermal infrared ozone profile: The ozone absorption band of 9.6 μm in the thermal infrared absorption band is used for inversion. The total ozone column and tropospheric column data are introduced to constrain the tropospheric and stratospheric ozone inversion results during the inversion process, and the Levenberg-Marquardt model is used to invert the tropospheric ozone profile. S5. Determination of proportional factors: The proportional factors λ and ε are determined by experimental methods. Using the existing ground-based, sounding and lidar ozone profile data as a benchmark, the root mean square error and root mean square relative error of the ozone profiles inverted in different regions with different proportional factors and the benchmark values ​​are analyzed to obtain reasonable allocation of proportional factors λ and ε.

2. The ozone profile inversion method based on thermal infrared remote sensing according to claim 1 is characterized in that: In step S4, the cost function of adding the total amount of ozone column is redefined as: Among them, vector X = X (O3) is the n-layer ozone concentration profile, Y vector is the brightness temperature observed by the m inversion channel satellite, γ is the introduced Lagrange multiplier, λ and ε are proportional coefficients, TO3 and TO tropo are the total ozone column and the total tropospheric ozone column, tropo represents the number of layers of tropospheric pressure, F(X) is the atmospheric forward radiation transmission mode, and the norm here is the Mahalanobis distance. The cost function is written as: R is for Y m The error covariance matrix includes the errors caused by the instrument and the system; B represents the initial error covariance matrix of the desired variable O3 profile. After defining the cost function, the iterative formula is derived and the constraints on the iterative formula are strengthened according to the Levenberg-Marquardt model. The positive definite matrix and the γ factor in the Levenberg-Marquardt method are adaptively adjusted through the diagonal elements and nonlinearity of the Hessian matrix. According to the principle of steepest descent, let J'(X) = 0 to obtain the linear solution form of this problem, namely: Among them, K is the first-order derivative of Y with respect to X, also known as the Jacobian matrix. Due to the nonlinearity of the inversion problem, the Gauss-Newton iteration method is used to improve the inversion accuracy and convergence speed of the problem, that is, the cost function is derived again to obtain: 0.5J”(X)=γB -1 +K T R -1 K (5) J”(X) is called the Hessian matrix. Substituting it into formula (4) yields: or: This method should give an accurate estimate of the covariance of the background error and set a suitable Lagrange multiplier γ, which is affected by human factors and inverts the variance of the error.

3. The ozone profile inversion method based on thermal infrared remote sensing according to claim 2 is characterized in that: In the Gauss-Newton iterative method, according to the least squares method proposed by Levenberg-Marquardt, a positive definite matrix is ​​constructed, Among them, H is a diagonal matrix, and w is re-determined in each iteration process, so that the cost function decreases rapidly.

4. The ozone profile inversion method based on thermal infrared remote sensing according to claim 1 is characterized in that: In step S1, thermal infrared hyperspectral data: the infrared hyperspectral detection instrument HIRAS is carried on the FY-3D, FY-3E and FY-3F satellites; the HIRAS carried on the FY3D is the first infrared hyperspectral detection instrument on the polar-orbiting meteorological satellite, with three wavelength ranges of long-wave infrared 8.80-15.39μm, medium-wave infrared 5.71-8.26μm and short-wave infrared 3.92-4.64μm, and the spectral resolution and spatial resolution are 0.625cm respectively -1 and 16km.

5. The ozone profile inversion method based on thermal infrared remote sensing according to claim 1 is characterized in that: In step S1, reanalysis data: ERA5 is a new generation of reanalysis data with a spatial resolution of 62 kilometers and a temporal resolution of 1 hour. It provides users with accurate real-time atmospheric information. The meteorological parameters provided include temperature, dew point temperature, air pressure, geopotential height, tropopause, and ozone column concentration data, covering different pressure ranges from the ground to the lower stratosphere.

Images