A method and system for simultaneous simulation and inversion of land surface temperature and land surface emissivity

Abstract

The present invention belongs to the field of climate and environmental monitoring technology and specifically provides a method and system for the simultaneous simulation and inversion of land surface temperature and land surface emissivity. The method comprises: simulating and generating global thermal infrared channel brightness temperature data based on observational parameters of thermal infrared channels from different geostationary meteorological satellites using the thermal infrared radiation transfer equation; adding random Gaussian noise to the brightness temperature data, using this as the observed brightness temperature, and performing simultaneous inversion of land surface temperature and emissivity. By combining observational data from geostationary meteorological satellites from different countries around the world and leveraging the multispectral and high-frequency observational characteristics of the new generation of geostationary meteorological satellites, a physical algorithm for simultaneous inversion of land surface temperature (LST) and land surface emissivity (LSE) is developed. This allows for simultaneous simulation and inversion of LST / LSE on a global scale, providing a reference for quantitative inversion of LST / LSE at high temporal and spatial resolutions worldwide. Satellite observations confirm that underlying surfaces with large diurnal temperature variations, such as deserts, are more conducive to the separation of land surface temperature and emissivity.

Description

Technical Field

[0001] The present invention relates to the field of climate and environmental monitoring technology, and more specifically, to a method and system for synchronous simulation and inversion of land surface temperature and land surface emissivity. Background Art

[0002] Land surface temperature and land surface emissivity are important underlying surface physical characteristic parameters. Land surface temperature (LST) refers to the temperature of the land surface and plays a vital role in the exchange of matter and energy between the earth and the atmosphere. It is a key indicator of the Earth's surface energy balance and the greenhouse effect, and has an important indicative role in regional and global climate change characteristics. The International Geosphere-Biosphere Program (IGBP) lists land surface temperature as one of its priority parameters for measurement. Land surface emissivity (LSE), also known as surface emissivity, is a measure of the ability of ground objects to radiate electromagnetic waves.

[0003] Remote sensing allows for large-scale, synchronous (quasi-synchronous) observations of the Earth. It is fast, low-cost, and offers high temporal, spatial, and spectral resolution, making it an indispensable data source for obtaining land surface temperature and emissivity. With the continuous development of remote sensing technology and the widespread application of thermal infrared remote sensing, remote sensing inversion of land surface temperature and emissivity has attracted increasing attention from scholars both domestically and internationally.

[0004] Currently, simultaneous retrieval of land surface temperature and emissivity from remote sensing falls into three main categories: emissivity spectral signature-based methods, physical inversion methods, and other methods. This article focuses on the widely used temperature / emissivity separation method (based on emissivity spectral signatures), the two-step physical inversion algorithm (based on physical inversion methods), and the daytime / nighttime operational method based on physical principles. The Temperature / emissivity Separation Method (TES) is a well-known algorithm for retrieving land surface temperature and emissivity from the ASTER satellite. It requires data from at least three atmospheric window channels (Gillespie et al. 1998; Abrams 2000) and can simultaneously retrieve land surface temperature and emissivity for all underlying surface types. This method primarily involves three steps: normalized emissivity method (NEM), spectral ratio (SR), and minimum-maximum difference (MMD). First, the normalized emissivity method is used to estimate initial values ​​of land surface temperature and land surface emissivity based on atmospherically corrected radiance values. Second, the spectral ratio method is used to calculate the ratio of the normalized emissivity relative to the mean. Finally, the minimum-maximum difference model is used to calculate the spectral contrast values ​​of N channels. The land surface emissivity inversion results are obtained using the empirical relationship between the minimum emissivity and the minimum-maximum difference of the N channels. A two-step physical inversion algorithm simplifies the radiative transfer equation using a special underlying surface reflectivity and a constant angular factor. Principal component analysis is used to reduce the dimensionality of the unknowns without compromising the accuracy of the inversion results. The method inverts land surface temperature and atmospheric temperature-humidity parameters in two steps. In the first step, given the initial land surface emissivity, land surface temperature, and atmospheric temperature-humidity profile, the Tikhonov regularization rule is used to establish a system of equations based on the radiative transfer equation. In the second step, Newton iteration is used to further improve the inversion results and better determine the solution to the system. This method does not require a priori atmospheric correction. The application of principal component analysis and Tikhonov regularization makes the inversion results more stable and accurate, and can simultaneously invert land surface temperature, land surface emissivity, and atmospheric profiles. However, the algorithm is relatively complex, computationally time-consuming, and inefficient. It requires a sufficient number of channel data and initial values ​​of land surface temperature, land surface emissivity, and atmospheric temperature-humidity profiles, and the results are significantly dependent on these initial values ​​(Li et al. 2007; Ma et al. 2000, 2002). The daytime / nighttime operational method based on physical principles simultaneously inverts land surface temperature and land surface emissivity without accurate land surface emissivity and atmospheric parameters. Wan and Li (1997) developed a physical inversion algorithm based on daytime / nighttime data pairs from mid-infrared and thermal infrared channels.The purpose of this method is to invert land surface temperature and land surface emissivity in arid and semi-arid regions, where land surface emissivity varies greatly in space (Wan 1999). Based on the assumptions about the optical properties of the land surface: (1) the land surface emissivity does not change significantly during the day / night period between days, unless rain or snow occurs within a short period of time; (2) in the mid-infrared band, the angular factor changes very little (<2%); (3) in the 3-14 μm range, the surface's reflection of downward solar radiation conforms to the Lambertian body approximation, and atmospheric thermal infrared radiation does not introduce significant errors. The method requires mid-infrared and thermal infrared data from multiple channels at different times, precise geometric correction, and prior atmospheric profile data. By continuously correcting the atmospheric profile data during the inversion process, the inversion accuracy of land surface temperature and land surface emissivity is greatly improved. The introduction of the mid-infrared band makes the inversion results more stable and more accurate. The inversion process is relatively complex, requires initial values, and the observation angles at different times are similar (Wan 1999; Wan and Li 1997).

[0005] However, with the need for regional and global weather, climate and ecological environment monitoring, higher requirements are placed on the temporal and spatial resolutions of LST / LSE remote sensing products. However, at present, relatively little work has been done on the inversion of satellite remote sensing LST / LSE products with high temporal and spatial resolutions on a global scale. Although polar-orbiting meteorological satellite remote sensing can provide regional or global LST / LSE products with a spatial resolution of 1 km or even 250 m, the temporal resolution is relatively low, with only four morning and afternoon satellite joint inversions per day. In terms of remote sensing LST / LSE synchronous inversion algorithms, in order to establish a set of equations that matches the number of inversion parameters as much as possible, some algorithms make full use of the multispectral characteristics of remote sensing data, some algorithms make full use of the multi-temporal characteristics of remote sensing data, and some algorithms make full use of the multi-phase characteristics of remote sensing data. Few algorithms make use of both characteristics at the same time. The widely used method based on the multispectral characteristics of remote sensing data introduces the same number of desired emissivity values ​​while introducing infrared radiation observation channels. It usually requires the use of some statistical empirical relationships or assumptions to perform the simultaneous inversion of land surface temperature and land surface emissivity. Moreover, each channel generally needs to undergo precise atmospheric correction. The established statistical empirical relationships or assumptions have a great influence on the inversion results of land surface temperature and land surface emissivity. Summary of the Invention

[0006] The present invention aims to solve the technical problem of difficult LST / LSE synchronous inversion in the prior art.

[0007] The present invention provides a method for synchronous simulation and inversion of land surface temperature and land surface emissivity, comprising the following steps:

[0008] S1, based on the thermal infrared channel observation parameters of different geostationary meteorological satellites, uses the thermal infrared radiation transfer equation to simulate and generate the brightness temperature data of the thermal infrared channel covering the global area;

[0009] S2, adding random Gaussian noise to the brightness temperature data, using it as the observed brightness temperature, and performing synchronous inversion of land surface temperature and land surface emissivity.

[0010] Preferably, the S1 specifically includes: simulating the brightness temperature of the global thermal infrared channel using the spectral response functions and related observation parameters of FY-4A AGRI, MSG SEVIRI and GOES-R ABI as inputs respectively;

[0011] The relevant observation parameters include the longitude, latitude, local zenith angle, and remote sensor noise (NeDR) matrix corresponding to the satellite remote sensing disk image.

[0012] Preferably, the simulation process in S1 adopts the precise transmittance model PFAAST, divides the atmosphere from 1100-0.05hPa into 101 layers, and calculates the contribution of the atmospheric gas absorption lines one by one in the order of the infrared spectrum.

[0013] Preferably, the S1 specifically includes:

[0014] It is assumed that within a 3-hour time range, the land surface temperature changes with time, while the land surface emissivity does not change with time;

[0015] Ignoring atmospheric scattering, the radiation from the Earth-atmosphere system in the infrared window band can be approximately expressed by the radiative transfer equation (RTE) as follows:

[0016]

[0017] Where R is the spectral radiation of the top of the atmosphere or the infrared radiation of the multi-channel scanning imaging radiometer, ε is the land surface emissivity, B(T) is the Planck function, τ(0, p) is the atmospheric transmittance from the top of the atmosphere to the atmosphere layer P, s represents the surface, is the atmospheric downward transmittance, e represents the uncertainty of the forward model, and R′ represents the reflected solar radiation.

[0018] Preferably, the S2 further includes:

[0019] Linearize RTE and ignore the influence of ozone and other trace gases. The first-order linearization of equation (1) is:

[0020]

[0021] where δR is the radiation perturbation, which is the difference between the observed radiation and the radiation calculated by the radiative transfer equation based on the initial values; K is a weighting function defined as Where x is the variable to be inverted, T S represents the land surface temperature, ε is the land surface emissivity, T is the air temperature, and Q is the water vapor, which shows the sensitivity of the top-of-atmosphere radiation to changes in the variable x; ∑ is the sum of different atmospheric layers; the logarithmic form of the mixing ratio is used instead of the direct mixing ratio because it has a better linear relationship with the radiation; e in formula (2) includes both the uncertainty of the forward model and the observation noise.

[0022] Preferably, the S2 specifically includes:

[0023] In the LST / LSE simultaneous simulation inversion algorithm of the multi-channel scanning imaging radiometer, the atmospheric contribution is represented by a variable:

[0024]

[0025] in It is a combination of temperature and humidity profiles. For each channel, the radiation error caused by the atmospheric profile error can be expressed as:

[0026]

[0027] make The new linearized equation is:

[0028]

[0029] The three radiation values ​​of the current time step T0, three hours ago T0-3, and six hours ago T0-6 are combined with formula (5) to solve LST, LSE and

[0030] Preferably, the 6-hour forecast fields with a resolution of 0.5° provided by NCEP GFS are used as initial values ​​for the temperature and humidity profiles, and each temperature and humidity profile is temporally and spatially interpolated to match the satellite data in time and position.

[0031] The present invention also provides a land surface temperature and land surface emissivity synchronous simulation and inversion system, characterized in that the system is used to implement a land surface temperature and land surface emissivity synchronous simulation and inversion method, including:

[0032] The brightness temperature simulation module is used to simulate and generate brightness temperature data of the thermal infrared channel covering the global area based on the observation parameters of the thermal infrared channel of different geostationary meteorological satellites and the thermal infrared radiation transfer equation;

[0033] The inversion module is used to add random Gaussian noise to the brightness temperature data and use it as the observed brightness temperature to perform synchronous inversion of land surface temperature and land surface emissivity.

[0034] The present invention also provides an electronic device comprising a memory and a processor, wherein the processor is configured to implement the steps of a method for synchronously simulating and inverting land surface temperature and land surface emissivity when executing a computer management program stored in the memory.

[0035] The present invention also provides a computer-readable storage medium storing a computer management program, which, when executed by a processor, implements the steps of a method for synchronous simulation and inversion of land surface temperature and land surface emissivity.

[0036] Beneficial Effects: The present invention provides a method and system for the simultaneous simulation and inversion of land surface temperature and land surface emissivity. The method comprises: simulating and generating global thermal infrared channel brightness temperature data based on observational parameters of thermal infrared channels from different geostationary meteorological satellites using the thermal infrared radiation transfer equation; adding random Gaussian noise to the brightness temperature data, using this as the observed brightness temperature, for the simultaneous inversion of land surface temperature and land surface emissivity. By combining observational data from geostationary meteorological satellites from different countries around the world and leveraging the multispectral and high-frequency observational characteristics of the new generation of geostationary meteorological satellites, a physical algorithm for simultaneous inversion of land surface temperature (LST) and land surface emissivity (LSE) is developed. This allows for simultaneous simulation and inversion of global-scale LST / LSE, providing a reference for quantitative inversion of global LST / LSE with high temporal and spatial resolution. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] Figure 1 A flow chart of a method for simultaneous simulation and inversion of land surface temperature and land surface emissivity provided by the present invention;

[0038] Figure 2 The global simulated brightness temperature maps of different thermal infrared channels of FY-4AGRI, MSG SEVIRI and GOES-R ABI provided by the present invention;

[0039] Figure 3 The global land surface emissivity map obtained by inverting the brightness temperature simulated by FY-4AGRI, MSG SEVIRI and GOES-R ABI provided by the present invention;

[0040] Figure 4 The global land surface temperature map obtained by inverting the brightness temperature simulation based on FY-4AGRI, MSG SEVIRI and GOES-R ABI provided by the present invention;

[0041] Figure 5The sensitivity map of the FY-4AAGRI and GOES-RABI synchronous simulation inversion land surface emissivity / land surface temperature results to the observation area provided by the present invention;

[0042] Figure 6 This is a sensitivity diagram of the FY-4AAGRI and GOES-RABI synchronous simulation inversion land surface emissivity / land surface temperature results to observation noise provided by the present invention;

[0043] Figure 7 The sensitivity of the FY-4AAGRI and GOES-RABI synchronous simulation inversion land surface emissivity / land surface temperature results to the observation time provided by the present invention;

[0044] Figure 8 A schematic diagram of the hardware structure of a possible electronic device provided by the present invention;

[0045] Figure 9 A schematic diagram of the hardware structure of a possible computer-readable storage medium provided by the present invention. DETAILED DESCRIPTION

[0046] The following embodiments of the present invention are described in further detail with reference to the accompanying drawings and examples. The following examples are used to illustrate the present invention but are not intended to limit the scope of the present invention.

[0047] Figure 1 The present invention provides a method for synchronously simulating and inverting land surface temperature and land surface emissivity, comprising the following steps:

[0048] S1, based on the thermal infrared channel observation parameters of different geostationary meteorological satellites, uses the thermal infrared radiation transfer equation to simulate and generate the brightness temperature data of the thermal infrared channel covering the global area;

[0049] S2, adding random Gaussian noise to the brightness temperature data, using it as the observed brightness temperature, and performing synchronous inversion of land surface temperature and land surface emissivity.

[0050] This solution leverages the significant improvements in spatial, temporal, and spectral resolution of new-generation geostationary meteorological satellites, making it possible to retrieve global land surface temperature and land surface emissivity. The European Meteosat Second Generation Spinning Enhanced Visible and Infrared Imager (MSG SEVIRI) offers a 3km spatial resolution and 15-minute temporal resolution for its thermal infrared channel, with three thermal infrared channels suitable for LST / LSE retrieval. The US Geostationary Operational Environmental Satellite-R Advanced Baseline Imager (GOES-RABI) offers a 2km spatial resolution and 15-minute temporal resolution for its thermal infrared channel, with five thermal infrared channels suitable for LST / LSE retrieval. China's FengYun-4AAGRI offers a 4km spatial resolution at the subsatellite point for its thermal infrared channel, a 15-minute temporal resolution for conventional observations, and three thermal infrared channels suitable for LST / LSE retrieval.

[0051] This embodiment of the present invention aims to combine the thermal infrared channels of MSG SEVIRI, GOES-RABI, and FY-4AGRI. Based on numerical simulation, the method first obtains simulated brightness temperatures from the three satellites, then performs a simultaneous inversion of global LST / LSE based on the simulated brightness temperatures, and conducts a sensitivity analysis of the results. The primary data sources used include information such as longitude, latitude, and local zenith angle corresponding to the satellite remote sensing disk image; the spectral response functions of the three thermal infrared channels; the noise of the remote sensing device (NeDR) matrix; Numerical Weather Prediction (NWP) forecast fields, analysis fields, or profile climate datasets; and the Earth Observing System Moderate-Resolution Imaging Spectroradiometer (EOS MODIS) monthly land surface emissivity dataset.

[0052] Specifically, a method based on simulated brightness temperature was used to simulate and invert global land surface temperature and land surface emissivity based on data from a new generation of geostationary meteorological satellites. First, the thermal infrared radiation transfer equation was used to simulate and generate global thermal infrared channel brightness temperature data based on relevant parameters from different geostationary meteorological satellites, namely FY-4AAGRI, MSG SEVIRI, and GOES-RABI. Random Gaussian noise was then added to the simulated brightness temperatures, which served as the observed brightness temperatures for the simultaneous inversion of land surface temperature and land surface emissivity. The brightness temperature simulation utilizes the Pressure-Layer Fast Algorithm for Atmospheric Transmittances (PFAAST). Based on the HITRAN2000 spectral data, PFAAST employs a line-by-line (LBL) method. This method divides the atmosphere from 1100 to 0.05 hPa into 101 layers and calculates the contributions of each atmospheric gas absorption line, line by line, in the order of the infrared spectrum (Clough and Iacono, 1995, Rothman et al., 1992). This model accounts for the satellite's zenith angle, absorption by mixed gases (including nitrogen, oxygen, and carbon dioxide), water vapor (including the water vapor continuum), and ozone. Consequently, PFAAST not only provides rapid computational speed but also achieves high brightness temperature simulation accuracy. The profile data required by the model are uniformly based on the Global Forecast System (GFS) data from the National Centers for Environmental Prediction (NCEP). The initial values ​​of land surface emissivity and land surface temperature are derived from the climate dataset (MODIS LSE monthly product) and the GFS, respectively. The specific input and output of brightness temperature simulation and its relationship with the LST / LSE inversion process are as follows: Figure 1 shown.

[0053] It should be noted that the embodiment of the present invention is based on the premise that within a certain time range (3 hours), the land surface temperature changes with time, while the land surface emissivity does not change with time.

[0054] Specifically, ignoring the scattering of the atmosphere, the radiation from the infrared window band of the Earth-atmosphere system can be approximately expressed by the radiative transfer equation (RTE) as follows:

[0055]

[0056] Where R is the spectral radiation of the top of the atmosphere or the infrared radiation of the multi-channel scanning imaging radiometer, ε is the land surface emissivity, B(T) is the Planck function, τ(0, p) is the atmospheric transmittance from the top of the atmosphere to the atmosphere layer P, s represents the surface, is the atmospheric downward transmittance, e represents the uncertainty of the forward model, and R′ represents the reflected solar radiation, which is usually ignored in the long-wave infrared window. As shown in Equation (1), the infrared radiation of the multi-channel scanning imaging radiometer has three main contributions: surface radiation, atmospheric upward radiation, and atmospheric downward radiation reflected from the surface.

[0057] The inversion problem of LST / LSE is to solve the land surface emissivity and surface temperature on the right side of Equation (1) given the observed radiation. Due to the contribution of the atmosphere to the radiation of the window channel, atmospheric correction is also necessary.

[0058] (1) Linearization of the radiation transfer equation

[0059] Due to the nonlinearity and ill-posedness of the inversion problem, there is no analytical solution for the inversion of LST / LSE and regularization is required. Therefore, the first step is to linearize RTE. Ignoring the effects of ozone and other trace gases, the first-order linearization of Equation (1) is:

[0060]

[0061] where δR is the radiation perturbation, which is the difference between the observed radiation and the radiation calculated by the radiative transfer equation based on the initial values. K is a weighting function defined as Where x is the variable to be inverted, indicating the sensitivity of top-of-atmosphere radiation to changes in x. ∑ is the sum of the different atmospheric layers. The logarithmic form of the mixing ratio is used rather than the direct mixing ratio because it provides a better linear relationship with the radiation. The e in Equation (2) includes both the uncertainty of the forward model and the observation noise. Equation (2) shows that the radiation perturbation consists of three components: LST, LSE, and the atmosphere (including temperature and humidity profiles). Perturbations in any of these components can lead to differences between the calculated and observed radiation.

[0062] (2) Atmospheric correction

[0063] Solving Equation (2) using only three thermal infrared channels is difficult. The above linear approximation places more complex demands on the atmospheric profile, especially the humidity profile. Therefore, it is necessary to simplify Equation (2) to eliminate significant errors, especially those related to the atmosphere. This simplified Equation (2) not only makes it easier to solve for LST and LSE, but also improves the accuracy of the inversion results.

[0064] The simplest way to simplify Equation (2) is to remove the atmospheric contribution (terms 3 and 4 on the right). This removal is equivalent to assuming that the atmospheric state is known and that the initial atmospheric profile represents the actual state well. The initial value can be a NWP forecast profile, satellite retrieval results, or a climate background.

[0065] In the LST / LSE algorithm of the multi-channel scanning imaging radiometer, the atmospheric contribution is represented by a variable:

[0066]

[0067] in is the combination of temperature and humidity profiles. For each channel, the radiation error caused by the atmospheric profile error can be expressed as:

[0068]

[0069] make The new linearized equation is:

[0070]

[0071] This formula will be used to solve LST, LSE and

[0072] (3) Temporal continuity

[0073] In general, assuming there are N channels, there are N+2 unknowns in equation (5): 1 LST, N LSEs, and 1 Therefore, for a single moment, the number of unknowns (N+2) is always greater than the number of equations (N). Therefore, it is difficult to obtain high inversion accuracy.

[0074] Due to the advantage of high temporal resolution, the LST / LSE algorithm of the embodiment of the present invention is based on the assumption that in a short period of time, the infrared LSE does not change with time, while the LST changes with time. Assume that the number of time steps is M, the total number of equations is M×N, and the number of unknowns is N+2M (each time step has 1 LST and 1 ). In order to obtain higher inversion accuracy, the number of equations should be equal to or greater than the number of unknowns, or:

[0075] M×N≥N+2M (6)

[0076] In the embodiment of the present invention, the number of channels is N=3, and the solution of equation (6) is M≥3; that is, at least 3 time steps are required.

[0077] Because experiments have shown that three time steps of 3 hours can better guarantee the inversion accuracy of the results, the radiation values ​​of the current time step (T0), three hours ago (T0-3), and six hours ago (T0-6) will be used together for the inversion of LST and LSE.

[0078] (4) Inversion algorithm

[0079] For the case of 3 time steps and 3 channels, there are 9 equations and 9 unknowns. Let:

[0080]

[0081] The superscript indicates the time step, and the table below indicates the channel index.

[0082] Formula (5) can be written as:

[0083] δY=KδX+e (7)

[0084] Where K is the linear or tangent model of the forward radiation transfer mode, also known as the Jacobian matrix or K matrix. The iterative solution of Equation (7) can be obtained using a simple least squares method:

[0085] δX n+1 =(K′ n E -1 K n ) -1 K′ n E -1 (δY n +K n δX n ) (8)

[0086] where δX n =X n -X0, δY n =Y m -Y(X n ), K n is the Jacobian matrix of the nth iteration, and E is the covariance matrix of the observation error, including the instrument noise and the uncertainty of the forward model. n is the vector of parameters to be inverted, X0 is the initial state or initial value, Y m is the vector of the observed radiation used in the inversion process, Y(X n ) is based on the atmospheric and surface conditions X n For a given initial value and satellite observations, if E′ n E -1 K n Reversible, the parameters can be inverted by formula (8).

[0087] However, since K′ n E -1 K n The singularity or near-singularity of δY makes the iteration unstable, so that Equation (8) has no solution or the solution does not exist. n Any noise in the inversion will be greatly amplified, and the inversion result will be non-existent. Therefore, an optimization estimation algorithm is needed to solve Equation (7). The general form of the solution is to minimize the following cost function (Rodgers 1976; Li et al., 2000).

[0088] J(X)=[Y m -Y(X)]′E -1 [Y m -Y(X)]+[X n -X0]′H[X n -X0] (9)

[0089] Where H is the prior matrix containing the solution, which can be the inverse of the initial value error covariance matrix or other types of matrices. By applying the following Newton iteration:

[0090] X n+1 =X n +J″(X n ) -1 ·J′(X n ) (10)

[0091] The following quasi-nonlinear iterative form can be obtained:

[0092] δX n+1 =(K′ n E -1 K n +H) -1 K′ n E -1 (δY n +K n δX n ) (11)

[0093] Compared with the least squares solution of equation (8), the only difference in equation (11) is the addition of an extra term H. The physical significance of this term is that it provides background information, so that the adjustment of the inversion parameters in the iterative process has a basis. The mathematical significance of this term is that in the matrix K′ n E -1 K n The addition of additional positive values ​​on the diagonal of reduces its singularity, making its inverse matrix (K′ n E -1 K n ) -1 Possibly existent and stable.

[0094] ① Initial value (X0)

[0095] For nonlinear, ill-posed inversion problems, the quality of the initial values ​​is crucial to inversion accuracy. In this simulation, the 6-hour forecast fields from the NCEP PFS with a 0.5° resolution were used as initial values ​​for the temperature and humidity profiles. Each profile was temporally and spatially interpolated to match the satellite data in time and location.

[0096] ② Initial value error covariance matrix

[0097] The error covariance matrix of the initial value must be consistent with the initial field. There is no universal covariance matrix that is suitable for all initial values. Ideally, the inverse matrix of the initial value error covariance matrix is ​​obtained by transforming the initial field error covariance matrix. Since the initial field related errors of LST and LSE are small, they can be ignored, that is, the elements on the non-diagonal line of the inverse matrix of the error covariance matrix can be set to 0. The diagonal elements are obtained by assuming that the error of LST is 10K, the errors of the three-channel LSE are 10%, 2% and 2% respectively, and The error is calculated to be 1K. Therefore, H is defined as:

[0098]

[0099] 0.01 is obtained by 1 / 10*10, and 100 is obtained by 1 / 0.1 / 0.1.

[0100] ③Observation error covariance matrix

[0101] Typically, the matrix E consists of two parts: observation noise and uncertainty of the radiation transfer model. The observation noise is estimated based on the instrument characteristics and is typically less than 0.15K for the three window channels. The uncertainty of the forward model is estimated by cross-comparison of different radiation transfer models, and the uncertainty of the three window channels is set to 0.2K. Similar to the inverse of the background field error covariance matrix, the matrix E is a diagonal matrix. This is equivalent to assuming that there is no correlated error between the observed radiation and the forward model uncertainty. The form of the matrix E is:

[0102]

[0103] where e i，j is the mixed error between the observation noise at the jth time step of the i-th channel and the forward model uncertainty.

[0104] ④Jacobian matrix

[0105] Jacobian matrix or k matrix K n(The subscript n indicates the nth iteration of the physical inversion process) represents the variation of the top-of-atmosphere radiation with the inversion parameters, which can be obtained by the difference method or the analytical method (Li et al., 2000). Considering the computational efficiency, the analytical method was used.

[0106] ⑤Convergence conditions

[0107] At each iteration, convergence is obtained by the increment of the last iteration, and the average increment is defined as:

[0108]

[0109] in is the increment of the i-th unknown quantity at the n-th iteration, and m is the number of unknowns. If the value is greater than a given threshold, the iteration continues until the maximum number of iterations is reached.

[0110] Based on the observation parameters of the three thermal infrared channels of FY-4AGRI, MSG SEVIRI and GOES-R ABI, and the PFAAST numerical model, the brightness temperature at 0:00, 03:00 and 06:00 UTC on March 28, 2018, covering the entire globe, is simulated. Figure 2 As shown in Figure 2, it can be seen that the spatial distribution of the simulation results shows the brightness temperature differences of different underlying surfaces well, and has good spatial continuity.

[0111] Figure 3 Shows the Figure 2 From the land surface emissivity of the corresponding channel obtained by the synchronous inversion of the simulated brightness temperature, it can be seen that the land surface emissivity of different underlying surfaces in a single channel is relatively reasonable. For example, in the 10.8μm channel, the Sahara Desert, the Australian desert and the sparse vegetation area all show low emissivity. By comparing different channels, it can be seen that the land surface emissivity generally increases with the increase of thermal infrared wavelength, which is also very reasonable.

[0112] Figure 4 Shows the Figure 2 From the land surface temperatures at different times obtained by the synchronous inversion of the simulated brightness temperature, it can be seen that the land surface temperatures in different regions and on different underlying surfaces at a single moment are quite reasonable. For example, at 03:00 UTC, Australia, India, and southern China all exhibited relatively high land surface temperatures. By comparing different moments, it can be seen that the overall process of the land surface temperature gradually increasing from 00:00 UTC to 06:00 UTC is also very reasonable.

[0113] Figure 5The accuracy of the simulated inversion results of land surface emissivity and land surface temperature by FY-4A AGRI and GOES-R ABI for different regions (Asia and Oceania, Europe and Africa, North America and South America, which correspond to the full disk observation areas of FY-4A AGRI, MSG SEVIRI and GOES-R ABI respectively) is shown. Figure 5 .a It can be seen that the simulation inversion accuracy of land surface emissivity of FY-4A AGRI and GOES-R ABI channels is basically the same in Asia and Oceania, North America and South America. However, in Europe and Africa, there is a slight difference in the inversion accuracy of the two channels, possibly due to the influence of the special underlying surface of the Sahara Desert. The inversion accuracy differences of FY-4A AGRI and GOES-R ABI for the 8.5μm, 10.8μm and 12.0μm channels are 0.0044, -0.0034 and -0.0018 respectively. Figure 5 As can be seen from the figure, for the land surface temperature simulation accuracy at the same time (03:00 UTC), the FY-4AAGRI result is lower than the GOES-R ABI result in all three regions, especially in North and South America, where the difference is -0.0022K, followed by Europe and Africa, where the difference is -0.0017K, and the smallest difference is in Asia and Oceania, where the difference is -0.0004K.

[0114] Three different Gaussian noises were added to the simulated brightness temperatures of the three thermal infrared channels of FY-4A AGRI and GOES-R ABI on March 28, 2018, with the standard deviations of the added noise being 0.5*NeDT+0.1 (NeDT is the instrument observation noise, and 0.1 is the model average noise), 1.0*NeDT+0.2, and 2*NeDT+0.4, respectively.

[0115] Figure 6The accuracy of the land surface emissivity and temperature (LT) derived from simulated brightness temperature inversions with varying amounts of noise added to the data is shown. The accuracy of both the FY-4A AGRI and GOES-R ABI LST inversions decreases with increasing noise, but the magnitude of the effect is small. The greatest impact on LST occurs when the noise level is 2*NeDT + 0.4, with the lowest accuracy at 8.7 μm, reaching less than 0.0006 for both the FY-4A AGRI and GOES-R ABI, respectively. Similarly, the greatest impact on LST occurs when the noise level is 2*NeDT + 0.4, with the lowest accuracy reaching 0.9904 K for the FY-4A AGRI and 1.5533 K for the GOES-R ABI. The accuracy of FY-4A AGRI is slightly higher than that of GOES-RABI, which may be because the thermal infrared channel band of FY-4A AGRI is slightly wider than that of GOES-RABI.

[0116] Finally, based on the simulated brightness temperatures of the three thermal infrared channels of FY-4A AGRI and GOES-R ABI every 3 hours on March 28, 2018, the land surface emissivity / land surface temperature were simultaneously inverted. Figure 7 The retrieval accuracy of land surface emissivity and temperature for both FY-4A AGRI and GOES-R ABI at different observation times is shown. It can be seen that the retrieval accuracy of land surface emissivity for each channel of the FY-4A AGRI and GOES-R ABI channels remains largely unchanged with observation time, while the retrieval accuracy of land surface temperature shows significant variation with observation time. This is primarily due to the distinct diurnal variation of land surface temperature. The retrieval accuracy of both FY-4A AGRI and GOES-R ABI increases gradually from 0:00 to 6:00, then decreases from 6:00 to 15:00, remains relatively stable from 15:00 to 21:00, and decreases from 21:00 to 24:00. However, the standard deviation of the retrieval results for both channels within a day is less than 0.082 K. At different times of the day, the accuracy of FY-4A AGRI is slightly higher than that of GOES-R ABI.

[0117] An embodiment of the present invention further provides a system for synchronous simulation and inversion of land surface temperature and land surface emissivity. The system is used to implement a method for synchronous simulation and inversion of land surface temperature and land surface emissivity, including:

[0118] The brightness temperature simulation module is used to simulate and generate brightness temperature data of the thermal infrared channel covering the global area based on the observation parameters of the thermal infrared channel of different geostationary meteorological satellites and the thermal infrared radiation transfer equation;

[0119] The inversion module is used to add random Gaussian noise to the brightness temperature data and use it as the observed brightness temperature to perform synchronous inversion of land surface temperature and land surface emissivity.

[0120] See also Figure 8 Schematic diagram of an embodiment of an electronic device provided by an embodiment of the present invention. Figure 8 As shown, an embodiment of the present invention provides an electronic device, including a memory 1310, a processor 1320, and a computer program 1311 stored in the memory 1310 and executable on the processor 1320. When the processor 1320 executes the computer program 1311, the following steps are implemented: S1, based on thermal infrared channel observation parameters of different geostationary meteorological satellites, using the thermal infrared radiation transfer equation, simulate and generate brightness temperature data of the thermal infrared channel covering the global area;

[0121] S2, adding random Gaussian noise to the brightness temperature data, using it as the observed brightness temperature, and performing synchronous inversion of land surface temperature and land surface emissivity.

[0122] See also Figure 9 Schematic diagram of an embodiment of a computer-readable storage medium provided by the present invention. Figure 9 As shown, this embodiment provides a computer-readable storage medium 1400 having a computer program 1411 stored thereon. When the computer program 1411 is executed by a processor, the following steps are implemented: S1, based on thermal infrared channel observation parameters of different geostationary meteorological satellites, using the thermal infrared radiation transfer equation, simulate and generate brightness temperature data of the thermal infrared channel covering the global area;

[0123] S2, adding random Gaussian noise to the brightness temperature data, using it as the observed brightness temperature, and performing synchronous inversion of land surface temperature and land surface emissivity.

[0124] Beneficial effects:

[0125] 1. The method used in this paper can be applied to the simultaneous simulation and inversion of global land surface emissivity and land surface temperature, and the results are quite reasonable.

[0126] 2. This paper's sensitivity experiments using simulated brightness temperatures to invert global land surface emissivity and land surface temperature demonstrate that the FY-4A AGRI's thermal infrared remote sensing system has comparable detection capabilities to the GOES-R ABI and slightly surpasses the GOES-R ABI in land surface temperature inversion accuracy. This provides theoretical support for the widespread application of the FY-4A AGRI's thermal infrared remote sensing technology.

[0127] 3. The characteristics of the temporal variation of the land surface temperature simulation inversion accuracy of the present invention show that it can be confirmed from the perspective of satellite observation that underlying surfaces with large daily temperature differences, such as deserts, are more conducive to the separation of land surface temperature and emissivity.

[0128] It should be noted that, in the above embodiments, the description of each embodiment has its own focus. For parts that are not described in detail in a certain embodiment, reference can be made to the relevant description of other embodiments.

[0129] Obviously, those skilled in the art may make various changes and modifications to the present invention without departing from the spirit and scope of the present invention. Thus, if such changes and modifications fall within the scope of the claims and their equivalents, the present invention is intended to include such changes and modifications.

Claims

1. A method for simultaneous simulation and inversion of land surface temperature and land surface emissivity, characterized in that: The following steps are involved: S1, based on the thermal infrared channel observation parameters of different geostationary meteorological satellites, uses the thermal infrared radiation transfer equation to simulate and generate the brightness temperature data of the thermal infrared channel covering the global area; S2, adding random Gaussian noise to the brightness temperature data, using it as the observed brightness temperature, and performing simultaneous inversion of land surface temperature and land surface emissivity; Wherein, the S1 specifically includes: It is assumed that within a 3-hour time range, the land surface temperature changes with time, while the land surface emissivity does not change with time; Ignoring atmospheric scattering, the radiation from the Earth-atmosphere system in the infrared window band can be approximately expressed by the radiative transfer equation (RTE) as follows: Where R is the spectral radiation of the top of the atmosphere or the infrared radiation of the multi-channel scanning imaging radiometer, ε is the land surface emissivity, B(T) is the Planck function, τ(0, p) is the atmospheric transmittance from the top of the atmosphere to the atmosphere layer P, s represents the surface, is the atmospheric downward transmittance, e represents the uncertainty of the forward model, R′ represents the reflected solar radiation, Ts represents the land surface temperature, B(Ts) represents the blackbody radiation value corresponding to the surface temperature Ts, τs represents the atmospheric transmittance from the surface to the top of the atmosphere, and Ps represents the atmospheric pressure on the land surface; Wherein, the S2 further includes: Linearize RTE and ignore the influence of ozone and other trace gases. The first-order linearization of equation (1) is: where δR is the radiation perturbation, which is the difference between the observed radiation and the radiation calculated by the radiative transfer equation based on the initial values; K is a weighting function defined as where x is the variable to be inverted, T S represents the land surface temperature, ε is the land surface emissivity, T is the air temperature, and Q is the water vapor, which shows the sensitivity of the top atmosphere radiation to the change of the variable x; ∑ is the sum of different atmospheric layers; the logarithmic form of the mixing ratio is used instead of the direct mixing ratio because it has a better linear relationship with the radiation; δε v represents the perturbation of the land surface emissivity ε of channel v relative to the initial value; the e in formula (2) includes both the uncertainty of the forward model and the observation noise; Wherein, the S2 specifically includes: In the LST / LSE simultaneous simulation inversion algorithm of the multi-channel scanning imaging radiometer, the atmospheric contribution is represented by a variable: in It is a combination of temperature and humidity profiles. For each channel, the radiation error caused by the atmospheric profile error can be expressed as: make The new linearized equation is: The three radiation values ​​of the current time step T0, three hours ago T0-3, and six hours ago T0-6 are combined with formula (5) to solve LST, LSE and 2. The method for simultaneous simulation and inversion of land surface temperature and land surface emissivity according to claim 1, characterized in that: Said S1 specifically comprises: simulating the brightness temperature of the global thermal infrared channel using the spectral response functions and related observation parameters of FY-4A AGRI, MSG SEVIRI and GOES-R ABI as inputs respectively; The relevant observation parameters include the longitude, latitude, local zenith angle, and remote sensor noise (NeDR) matrix corresponding to the satellite remote sensing disk image.

3. The method for simultaneous simulation and inversion of land surface temperature and land surface emissivity according to claim 1, characterized in that: The simulation process in S1 adopts the precise transmittance model PFAAST, divides the atmosphere from 1100-0.05hPa into 101 layers, and calculates the contribution of the atmospheric gas absorption lines one by one in the order of the infrared spectrum.

4. The method for simultaneous simulation and inversion of land surface temperature and land surface emissivity according to claim 1, characterized in that: The 6-hour forecast fields with a resolution of 0.5° provided by NCEP GFS are used as the initial values ​​of temperature and humidity profiles. Each temperature and humidity profile is temporally and spatially interpolated to match the satellite data in time and position.

5. A land surface temperature and land surface emissivity simultaneous simulation and inversion system, characterized in that: The system is used to implement the method for simultaneous simulation and inversion of land surface temperature and land surface emissivity according to any one of claims 1 to 4, comprising: The brightness temperature simulation module is used to simulate and generate brightness temperature data of the thermal infrared channel covering the global area based on the observation parameters of the thermal infrared channel of different geostationary meteorological satellites and the thermal infrared radiation transfer equation; The inversion module is used to add random Gaussian noise to the brightness temperature data and use it as the observed brightness temperature to perform synchronous inversion of land surface temperature and land surface emissivity.

6. An electronic device, characterized in that: The method comprises a memory and a processor, wherein the processor is used to implement the steps of the method for synchronous simulation and inversion of land surface temperature and land surface emissivity as described in any one of claims 1 to 4 when executing a computer management program stored in the memory.

7. A computer-readable storage medium, characterized in that A computer management program is stored thereon, and when the computer management program is executed by the processor, the steps of the method for synchronous simulation and inversion of land surface temperature and land surface emissivity as described in any one of claims 1 to 4 are implemented.

Images