Cooperative retrieval method for static satellite land surface emissivity based on background field reconstruction

Abstract

The application discloses a static satellite ground surface emissivity cooperative inversion method based on background field reconstruction and belongs to the technical field of satellite remote sensing information processing, which comprises the following steps: constructing a background field of spatial structure cooperation constraint of ground surface temperature and atmospheric temperature, adopting a three-dimensional variation assimilation cost function to assimilate the measured ground surface temperature of a ground meteorological observation station, obtaining a ground surface temperature analysis field as a ground surface temperature initial guess field of a one-dimensional variation inversion cost function to invert the ground surface emissivity; generating a ground surface emissivity background field from the ground surface emissivity in a ground surface element historical data set, constructing a ground surface item background error covariance matrix capable of representing the coupling relationship between the ground surface emissivity and the ground surface temperature based on the ground surface element historical data set, and further constructing a one-dimensional variation inversion cost function; and combining the satellite observation brightness temperature of a clear sky observation pixel observed by an infrared window channel of a geostationary orbit radiometer, so that the inversion of the ground surface emissivity can be realized, and thus the ground surface emissivity inversion level is improved.

Description

Technical Field

[0001] This invention relates to a method for co-inversion of geostationary satellite surface emissivity based on background field reconstruction, belonging to the field of satellite remote sensing information processing technology. Background Technology

[0002] The Fengyun-4 Advanced Geostationary Radiometric Imager (FY-4A / AGRI) features high spatiotemporal resolution. Its infrared detection channel, with a horizontal spatial resolution of 4 km, a full-disk scan time of 15 minutes, and advantageous channel design, provides rich atmospheric and surface observational information for numerical weather prediction (NWP), enhancing forecast accuracy. To effectively assimilate infrared channel observational data in NWP models, atmospheric and surface parameters need to be used as input information for the radiative transfer model (observation operator) to simulate the brightness temperature of the corresponding channel.

[0003] Under clear-sky atmospheric conditions, radiative transfer models require atmospheric temperature, humidity, and pressure parameters, as well as surface information including surface temperature and emissivity. Additionally, geometric parameters from satellite observations, such as the observed zenith angle, are necessary inputs. The amount of atmospheric and surface radiation received by satellite instruments is calculated using radiative transfer models and converted into satellite-observed brightness temperature. The center wavelengths of channels 11, 12, and 13 of FY-4A / AGRI are 8.5 µm, 10.7 µm, and 12 µm, respectively. Their weighting functions are near the ground and sensitive to surface information. Under clear-sky conditions, atmospheric radiation has a relatively small impact, and the accuracy of land surface brightness temperature simulation mainly depends on the accuracy of surface emissivity and surface temperature. Surface temperature can be obtained by assimilating ground observation station data. High-precision surface temperature assimilated from surface observations can be obtained by correcting numerical weather prediction model data using surface temperature observations from ground meteorological stations. Surface emissivity is an important physical parameter describing the ability of the Earth's surface to radiate energy outwards and directly affects the radiative brightness temperature value obtained from satellite infrared remote sensing observations. Satellite data assimilation fuses instrumental observation data with numerical model data, while accurate inversion of surface emissivity can reduce systematic biases between observational and model data, thus improving assimilation efficiency. Therefore, improving the accuracy of surface emissivity is one of the key factors in enhancing brightness temperature simulation.

[0004] The inversion of infrared surface emissivity (abbreviated as surface emissivity) mainly employs the split-window algorithm and the surface temperature-emissivity separation algorithm. The split-window algorithm separates surface temperature and emissivity based on the difference in atmospheric absorption between two adjacent thermal infrared channels (such as channels 31 and 32 in MODIS) using linear or nonlinear models. In recent years, the split-window algorithm has significantly improved inversion accuracy by introducing dynamic atmospheric correction parameters (such as water vapor content) and surface type classification (such as vegetation cover). The generalized split-window algorithm has become a core method in MODIS standard products. The temperature-emissivity separation algorithm utilizes multi-channel thermal infrared data and decouples temperature and emissivity through iterative optimization or statistical constraints (such as the minimum emissivity variation method). The TES algorithm significantly reduces inversion uncertainty by normalizing emissivity assumptions (such as 0.98 emissivity for the 13 µm channel).

[0005] However, the acquisition of surface emissivity mainly relies on satellite remote sensing inversion and ground observations. Existing methods still have the following shortcomings: 1) Insufficient consideration of temporal variations: Surface emissivity is affected by seasonal, weather, and land cover changes, and existing methods do not adequately characterize these temporal variations. 2) Insufficient integration with surface observation data: There are few studies on the integration of surface temperature observation data with satellite remote sensing data, leading to significant uncertainties in the surface emissivity inversion results. 3) According to the Planck function, surface emissivity and surface temperature are coupled. Surface temperature errors have a linear impact on the brightness temperature (BT) of lower atmospheric channels; a 6 K error in land temperature can amplify the emissivity inversion error by more than 30%. The nonlinear impact of surface temperature errors on emissivity inversion is particularly problematic in complex surface regions. Therefore, the inversion of infrared emissivity is highly dependent on the accuracy of surface temperature. Summary of the Invention

[0006] The technical objective of this invention is to provide a geostationary satellite surface emissivity co-inversion method based on background field reconstruction. On one hand, it considers the coupling relationship between surface temperature and surface emissivity to optimize the existing 1DVAR (one-dimensional variational inversion algorithm, integrated into a one-dimensional variational inversion system) surface emissivity retrieval, solving the surface temperature-surface emissivity coupling optimization problem and improving the simulation accuracy of satellite data. On the other hand, this invention also introduces a surface temperature term to reconstruct the background field with coordinated constraints on the spatial structure of surface temperature and atmospheric temperature. Then, it utilizes three-dimensional variational assimilation technology (integrated into a three-dimensional variational assimilation system) to assimilate more accurate measured surface temperatures from ground meteorological stations, obtaining a more accurate surface temperature analysis field as the initial guess field for the surface temperature in the one-dimensional variational inversion system. This improves the 1DVAR surface emissivity retrieval scheme and further enhances the simulation accuracy of satellite data.

[0007] To achieve the above-mentioned technical objectives, the present invention will adopt the following technical solution:

[0008] A method for co-inversion of geostationary satellite surface emissivity based on background field reconstruction includes the following steps:

[0009] The measured land surface temperature and the traditional land surface temperature background field and atmospheric temperature and humidity profiles output by the numerical weather prediction model were obtained respectively.

[0010] For the traditional land surface temperature background field output by numerical weather prediction models, a land surface temperature term is introduced to construct a background field that is jointly constrained by the spatial structure of atmospheric temperature and land surface temperature, thereby reconstructing the background field and obtaining an optimized land surface temperature background field.

[0011] Based on the measured surface temperature and the optimized surface temperature background field, an optimized three-dimensional variational assimilation cost function is constructed to assimilate the measured surface temperature and output the surface temperature analysis field.

[0012] Obtain surface features within a specified historical period to construct a historical dataset of surface features; surface features include surface temperature and surface emissivity;

[0013] Based on historical datasets of surface elements, a background error covariance matrix of surface terms is constructed to characterize the cooperative constraint relationship between surface emissivity and surface temperature. Based on the obtained background error covariance matrix of surface terms, a one-dimensional variational inversion cost function is constructed.

[0014] Obtain the satellite observation brightness temperature of clear-sky observation pixels observed in the infrared window channel of the geostationary orbit radiometric imager;

[0015] The surface emissivity from the obtained historical dataset of surface elements is used to generate the surface emissivity background field required as input to the one-dimensional variational inversion cost function. The surface temperature analysis field output by the optimized three-dimensional variational assimilation cost function is used as the initial surface temperature guess field required as input to the one-dimensional variational inversion cost function. By combining the atmospheric temperature and humidity profiles and the satellite observation brightness temperature of the observed clear sky pixels, the independent inversion of surface emissivity or the joint inversion of surface emissivity and surface temperature can be achieved, and then the optimized value of surface emissivity can be output.

[0016] As a further optimization of the present invention, a one-dimensional variational inversion cost function Represented as:

[0017] ;

[0018] In the formula: It is the surface state vector; It is the background field vector of the Earth's surface term; The background error covariance matrix for the surface term; It is a satellite observation of brightness temperature; It is the background brightness temperature output by the observation operator H in each iteration; E is the observation error covariance matrix of the brightness temperature term;

[0019] Surface term background error covariance matrix Represented as:

[0020] ;

[0021] In the formula: This is the correlation coefficient between surface emissivity and surface temperature, obtained through sample statistical calculations based on historical datasets of surface elements. The error variance of surface emissivity is obtained through sample statistical calculations based on historical datasets of surface elements. It is the error variance of surface temperature obtained by statistical calculation of samples based on historical datasets of surface elements;

[0022] Background brightness temperature Represented as:

[0023] ;

[0024] In the formula: It is the surface term state vector, including the surface temperature analysis field. Atmospheric temperature and humidity profiles and surface emissivity at the current iteration number .

[0025] As a further optimization of this invention, when iterating the one-dimensional variational inversion cost function, the Levenberg-Marquardt algorithm is used to minimize the one-dimensional variational inversion cost function. The sensitivity of surface emissivity and surface temperature to satellite-observed brightness temperature is calculated using the Jacobian matrix, and the surface emissivity is iteratively updated. Until it converges.

[0026] As a further optimization of the present invention, a one-dimensional variational inversion cost function The iterative solution involves the following steps:

[0027] Using a given surface emissivity background field and the initial guess of the surface emissivity field The initial surface temperature field was estimated by simulating the initial background brightness temperature using observation operators. Initially, define The surface emissivity error at the initial moment ;

[0028] Calculate the surface emissivity error in the m-th iteration within the iterative loop. Then, with the surface emissivity background field Add them together to calculate the surface emissivity in the m-th iteration. The background brightness temperature was calculated using the observation operator. If the surface emissivity error is at this time If the convergence condition has been met, the iteration ends, and the surface emissivity of the m-th iteration is output. If the surface emissivity is used as the optimal value, otherwise continue to the (m+1)th iteration;

[0029] In the iterative solution of the one-dimensional variational inversion cost function, the convergence degree is used when determining the convergence condition. test:

[0030] ;

[0031] Where N represents the number of infrared window channels used in the iterative calculation; It is a satellite observation of brightness temperature; is the background brightness temperature calculated by the observation operator in the m-th iteration; E is the observation error covariance matrix of the brightness temperature term; the superscript T indicates matrix transpose.

[0032] As a further optimization of the present invention, the extended background error covariance matrix... It is expressed as follows:

[0033] ;

[0034] In the formula: B represents the atmospheric state vector flow function. Velocity potential function Atmospheric temperature relative humidity air pressure and surface temperature The background error covariance matrix of cocorrelation; It is the identity matrix; Represents the velocity potential function and atmospheric state vector stream function The relationship between them Indicates atmospheric temperature With atmospheric state vector stream function The relationship between them Represents air pressure With atmospheric state vector stream function The relationship between them surface temperature Atmospheric temperature of the non-equilibrium part The relationship between them; , and Represent the velocity potential function, atmospheric temperature, and air pressure of the non-equilibrium component, respectively. This represents the non-equilibrium portion of the Earth's surface temperature.

[0035] Surface temperature Atmospheric temperature of the non-equilibrium part The relationship between Represented as:

[0036] ;

[0037] ;

[0038] In the above formula: Indicates coordinate location The non-equilibrium portion of the Earth's surface temperature; Indicates coordinate location The equilibrium portion of the Earth's surface temperature; Represents the horizontal coordinate; Represents the vertical coordinate; Indicates coordinate location The corresponding region is divided into partition numbers based on preset latitude intervals; Indicates coordinate location In partition number index The following partition; Indicates the total number of partitions; Indicates surface temperature Atmospheric temperature of the non-equilibrium part The linear correlation coefficient is calculated using the following formula:

[0039] ;

[0040] In the above formula: Indicated based on surface temperature Atmospheric temperature of the non-equilibrium part The covariance matrix obtained from historical statistics; Represents atmospheric temperature based on the non-equilibrium component. The variance matrix obtained from historical statistics.

[0041] As a further optimization of the present invention, the optimized three-dimensional variational assimilation cost function Represented as:

[0042] ;

[0043] ;

[0044] ;

[0045] ;

[0046] In the formula: Background penalty item; This is a penalty for surface temperature observation. This is a penalty for observations other than surface temperature;

[0047] x represents the value including surface temperature analysis. The state vector; x b Includes surface temperature forecast values The background field vector; For the extended background error covariance matrix;

[0048] y is the value excluding the measured surface temperature. Other observation vectors besides H; H is the observation operator; It is the observation error covariance matrix, and the superscript T indicates matrix transpose;

[0049] It is the measured surface temperature.

[0050] As a further optimization of the present invention, atmospheric temperature based on the non-equilibrium component is achieved through Cholesky decomposition. variance matrix obtained from historical statistics Find the inverse;

[0051] The gradient of the three-dimensional variational assimilation cost function is obtained using the conjugate gradient descent method. Find the gradient of the gradient of the 3D variational assimilation cost function during gradient descent. The numerical value and the gradient of the initial three-dimensional variational assimilation cost function The algorithm converges when the ratio is less than a preset value; the corresponding state vector at convergence. This is the optimal surface temperature analysis field after minimizing the three-dimensional variational assimilation cost function, which includes the analytical values ​​of atmospheric elements and surface temperature. Atmospheric elements include air temperature, humidity, and air pressure.

[0052] As a further optimization of the present invention, based on the optimized surface emissivity value obtained by inversion of the one-dimensional variational inversion cost function, the simulated brightness temperature of the satellite observation point is first calculated by simulating different infrared window channels of the geostationary orbit radiation imager using the radiative transfer mode. Then, by calculating the root mean square error between the simulated brightness temperature of the corresponding infrared window channel and the satellite observed brightness temperature, the simulation capability of the optimized surface emissivity value for the brightness temperature observation data of the infrared window channel can be verified.

[0053] As a further optimization of the present invention, when calculating the simulated brightness temperature of different infrared window channels of the geostationary radiometric imager through the radiative transfer mode, the conventional atmospheric elements and surface elements of the spatial grid points output by the numerical forecast model and the optimized surface emissivity value obtained by inversion are interpolated to the satellite observation points of the geostationary radiometric imager as input to the radiative transfer mode, thereby obtaining the simulated brightness temperature of the infrared window channels of the geostationary radiometric imager.

[0054] As a further optimization of the present invention, the infrared window channel of the geostationary orbital radiation imager includes channels 11, 12 and 13.

[0055] Based on the above-mentioned technical objectives, the present invention has the following advantages compared with the prior art:

[0056] The geostationary satellite surface emissivity collaborative inversion method based on background field reconstruction described in this invention combines a three-dimensional variational assimilation method (implemented based on a three-dimensional variational assimilation system) and a 1DVAR inversion algorithm (i.e., a one-dimensional variational inversion algorithm implemented based on a one-dimensional variational inversion system). It also incorporates the satellite observation brightness temperature of clear-sky pixels observed in the infrared window channel of a geostationary orbit radiometric imager (this invention uses the satellite observation brightness temperature of channels 11-13 of the FY-4A AGRI). This enables the inversion of surface emissivity and surface temperature, resolving the coupling problem between surface emissivity and surface temperature, and improving the poor applicability of traditional methods in complex surface regions (traditional methods mostly perform inversion by decoupling surface emissivity and surface temperature, but surface temperature errors have a nonlinear effect on surface emissivity inversion, thus traditional methods often result in distorted surface emissivity inversion due to surface temperature errors).

[0057] Specifically, this invention reconstructs the traditional surface temperature background field output by numerical weather prediction models by co-constraining the spatial structure of surface temperature and atmospheric temperature, thereby achieving background field reconstruction. Then, the measured surface temperature (obtained through observations at meteorological ground stations) and numerical model data are assimilated using three-dimensional variational assimilation technology to obtain an accurate surface temperature analysis field. The obtained surface temperature analysis field is used as the initial guess field of surface temperature in a one-dimensional variational inversion system, establishing a surface emissivity inversion scheme based on a one-dimensional variational inversion algorithm. This reduces the impact of surface temperature on surface emissivity, solves the problem of surface emissivity inversion distortion caused by surface temperature errors in traditional methods, and improves the accuracy of satellite data assimilation and surface emissivity inversion. Attached Figure Description

[0058] Figure 1 This is a flowchart of the geostationary satellite surface emissivity co-inversion method based on background field reconstruction and infrared brightness temperature simulation verification method described in this invention.

[0059] Figure 2This is a graph showing the root mean square error (RMSE) (in K) of satellite-observed brightness temperatures of different infrared window channels of FY-4A AGRI under clear sky conditions and the simulated brightness temperatures of the corresponding infrared window channels obtained by using different surface emissivity datasets as inputs for radiative transfer modes; in the graph: (a) represents channel 11, (b) represents channel 12, and (c) represents channel 13. Detailed Implementation

[0060] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. The following description of at least one exemplary embodiment is merely illustrative and is in no way intended to limit the present invention or its application or use. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention. Unless otherwise specifically stated, the relative arrangement, expressions, and values ​​of components and steps set forth in these embodiments do not limit the scope of the present invention. Techniques, methods, and devices known to those skilled in the art may not be discussed in detail, but where appropriate, such techniques, methods, and devices should be considered part of the specification. In all examples shown and discussed herein, any specific values ​​should be interpreted as merely exemplary and not as limitations. Therefore, other examples of exemplary embodiments may have different values.

[0061] like Figure 1 As shown, the geostationary satellite surface emissivity co-inversion method based on background field reconstruction disclosed in this invention introduces a surface temperature term into the traditional surface temperature background field output by numerical weather prediction models to construct a background field with co-constrained spatial structure of surface temperature and atmospheric temperature, thereby obtaining an optimized surface temperature background field. Then, an optimized three-dimensional variational assimilation cost function (integrated in a three-dimensional variational assimilation system) is constructed using the measured surface temperature and the optimized surface temperature background field. This optimized three-dimensional variational assimilation cost function is then used to assimilate the measured surface temperature, outputting a surface temperature analysis field. Finally, the obtained surface temperature analysis field is used as a one-dimensional variational inversion cost function (integrated in a one-dimensional variational inversion system) to invert the initial surface temperature guess field for surface emissivity, realizing surface emissivity inversion considering the co-constraints of surface temperature and surface emissivity. Specifically, it includes the following steps:

[0062] Step 1: Construction and optimization of the initial temperature prediction field from ground observation stations:

[0063] This step is based on the traditional land surface temperature background field derived from measured land surface temperatures (observed from ground meteorological station data) and model outputs (i.e., exponential forecast model outputs). A land surface temperature term is introduced into this traditional background field to construct a background field with co-constrained spatial structure of land surface and atmospheric temperatures. This allows for optimization of the three-dimensional variational assimilation cost function, enabling the assimilation of measured land surface temperatures to obtain a high-precision land surface temperature analysis field. This provides an optimized initial land surface temperature guess as input for subsequent separate inversions of land surface emissivity or joint inversions of land surface emissivity and land surface temperature. The specific process is as follows:

[0064] Step 1.1, Data Processing and Fusion:

[0065] Obtain measured surface temperature (Source: National-level surface meteorological stations) and traditional surface temperature background field (From the WRF numerical weather prediction model, whose output includes land surface temperature forecast values) ).

[0066] measured surface temperature Perform rapid quality control, including threshold testing (-50). o C to 70 o C) Background field outlier test (removal) , The tests include: ground meteorological station observation error, spatial consistency test (removing data with a difference of >10°C from the mean of the surrounding 220 km), temporal consistency test (removing data with a difference of >5°C from the mean of the preceding and following hour), and topographic height consistency test (removing data when the difference between the model and the measured height is >50 m).

[0067] Step 1.2, Construction of the background field for spatial structural collaborative constraints:

[0068] For the traditional land surface temperature background field output by numerical weather prediction models, a land surface temperature term is introduced to construct a background field that is synergistically constrained by the spatial structure of atmospheric temperature and land surface temperature, thereby reconstructing the background field and obtaining an optimized land surface temperature background field.

[0069] Specifically, in the atmospheric state vector stream function Velocity potential function Atmospheric temperature relative humidity and air pressure Based on the correlation, increase the surface temperature With atmospheric temperature The cocorrelated background error covariance matrix (denoted as the extended background error covariance matrix) Extended background error covariance matrix It is expressed as follows:

[0070] (1)

[0071] in, , and The corresponding values ​​are the velocity potential function, atmospheric temperature, and air pressure in the non-equilibrium part. Represents the velocity potential function and atmospheric state vector stream function The relationship between them Indicates atmospheric temperature With atmospheric state vector stream function The relationship between them Represents air pressure With atmospheric state vector stream function The relationship between them It is an identity matrix. This represents the non-equilibrium portion of the Earth's surface temperature. surface temperature Atmospheric temperature of the non-equilibrium part The relationship between them can be expressed in functional form as follows:

[0072] (2)

[0073] (3)

[0074] in, Indicates coordinate location The non-equilibrium portion of the Earth's surface temperature; Indicates coordinate location The equilibrium portion of the Earth's surface temperature; Represents the horizontal coordinate; Represents the vertical coordinate; Indicates coordinate location The corresponding region is divided into partition numbers based on preset latitude intervals; Indicates coordinate location In partition number index The following partition; Indicates the total number of partitions; Indicates surface temperature Atmospheric temperature of the non-equilibrium part The linear correlation coefficient.

[0075] Therefore, in this invention, the coordinate position The equilibrium part of the ground surface temperature By temperature The linear regression represents this part, reflecting the surface temperature. Atmospheric temperature of the non-equilibrium part The correlation between the two is determined by using a statistical sample to calculate the linear correlation coefficient. Furthermore, considering the regional differences in surface temperature, different zones are established at preset latitudinal intervals (in this invention, the latitudinal interval is 5°). It describes the characteristics under different latitude conditions.

[0076] Surface temperature Atmospheric temperature of the non-equilibrium part linear correlation coefficient Calculated using the following formula:

[0077] (4)

[0078] In the above formula: Indicated based on surface temperature Atmospheric temperature of the non-equilibrium part The covariance matrix obtained from historical statistics; Represents atmospheric temperature based on the non-equilibrium component. The variance matrix obtained from historical statistics.

[0079] Step 1.3: Three-dimensional variational assimilation optimization of the initial guess of surface temperature field.

[0080] Based on the optimized surface temperature background field, an optimized three-dimensional variational assimilation cost function is constructed. (Integrated into the three-dimensional variational assimilation system), and employing this optimized three-dimensional variational assimilation cost function. The measured surface temperature is assimilated to output a surface temperature analysis field. The resulting surface temperature analysis field is the initial guess field for the surface temperature required for subsequent inversion.

[0081] Optimized 3D variational assimilation cost function Compared to the traditional three-dimensional variational assimilation cost function, in the background penalty term... And observation penalties other than surface temperature Based on this, a new penalty item for surface temperature observation has been added. That is, the three-dimensional variational assimilation cost function used in this invention. Represented as:

[0082] (5)

[0083] In the formula, the background field penalty term is:

[0084] (6)

[0085] Observational penalties other than surface temperature:

[0086] (7)

[0087] Newly added penalty for surface temperature observation:

[0088] (8)

[0089] In the above formula: x b Includes surface temperature forecast values The background field vector, x, contains the surface temperature analysis value. The state vector, The extended background error covariance matrix (which, compared to the traditional background error covariance matrix, incorporates surface temperature) Atmospheric temperature of the non-equilibrium part The relationship between Therefore, it can be denoted as the extended background error covariance matrix, where y is the observation vector excluding the measured surface temperature, and H is the observation operator; It is the measured surface temperature. It is the observation error covariance matrix, and the superscript T indicates matrix transpose.

[0090] The three-dimensional variational assimilation cost function used in this invention is solved iteratively using the conjugate gradient method. The minimum value of the gradient is reached when the gradient converges to 1 / 1000 of the initial value. The high-precision surface temperature analysis field after fusion is then output as the initial guess field for the surface temperature of the subsequent one-dimensional variational inversion cost function inversion of surface emissivity.

[0091] Step 2, one-dimensional variational inversion of surface emissivity:

[0092] Based on the background field reconstructed by the joint constraints of the spatial structure of surface temperature and atmospheric temperature, a three-dimensional variational assimilation cost function is used to assimilate the measured surface temperature from ground meteorological observation stations, resulting in the surface temperature analysis field (which includes surface temperature analysis values). ), combined with satellite-observed brightness temperature of clear-sky observation pixels (Obtained from infrared channel observations using FY-4A / AGRI; clear-sky observation pixels were obtained from clear-sky markers corresponding to cloud mask products). Surface emissivity was jointly retrieved using a one-dimensional variational inversion algorithm. The specific steps are as follows:

[0093] Step 2.1: Select clear sky observation pixels from FY-4A AGRI.

[0094] Using FY-4A AGRI's secondary cloud mask product CLM (which categorizes pixels into four types: cloud, possible cloud, possible clear sky, and clear sky) as a reference, when a pixel is a clear sky observation and the surrounding four pixels are also clear sky observations, the pixel is determined to be a clear sky observation pixel and is retained.

[0095] Step 2.2: Construct the covariance matrix of the surface emissivity background field and the surface term background error.

[0096] (1) Construct the surface emissivity background field.

[0097] Land surface features within a specified historical period are acquired to construct a historical dataset of land surface features, including land surface temperature and land surface emissivity. Therefore, the historical dataset comprises two subsets: a historical dataset of land surface temperature and a historical dataset of land surface emissivity. The land surface temperature data in the historical dataset is derived from the model output of a numerical weather prediction model. The historical data of land surface emissivity is derived from the historical background field dataset CAMEL Climate. Simultaneously, atmospheric temperature and humidity profiles from the ERA5 reanalysis data are acquired.

[0098] The surface emissivity background field is generated based on the surface emissivity in the historical dataset of surface features, i.e., based on the historical surface emissivity dataset. This is used as the initial guess of the surface emissivity of the one-dimensional variational inversion cost function, while the surface temperature analysis field output by the three-dimensional variational assimilation cost function in step one is used as the initial guess of the surface temperature field of the one-dimensional variational inversion cost function.

[0099] (2) Construct a surface term background error covariance matrix that can characterize the coupling relationship between surface emissivity and surface temperature. .

[0100] By analyzing the error distribution of each element in the historical dataset of surface elements, a background error covariance matrix of surface terms is constructed to characterize the coupling relationship between surface emissivity and surface temperature. Surface term background error covariance matrix The format is as follows:

[0101] (9)

[0102] In the formula: This is the correlation coefficient between surface emissivity and surface temperature, obtained through sample statistical calculations based on historical datasets of surface elements. The error variance of surface emissivity is obtained through sample statistical calculations based on historical surface emissivity datasets. It is the error variance of surface temperature obtained by statistical calculation of samples based on historical surface temperature datasets.

[0103] Step 2.3, implementation of the one-dimensional variational inversion algorithm.

[0104] Step 2.3.1: Establish the one-dimensional variational inversion cost function :

[0105] (10)

[0106] in: It is the surface term state vector, which includes the surface temperature analysis field, the atmospheric temperature and humidity profile (from Global Numerical Weather Prediction Model Analysis Data GFS), and the surface emissivity at the current iteration number. ; It is the background field vector of the Earth's surface term; The background error covariance matrix for the surface term; It is a satellite-observed brightness temperature, obtained based on the infrared window channel observation of a geostationary orbit radiation imager; is the background brightness temperature output by the observation operator H (i.e., the radiative transfer mode) at each iteration; E is the observation error covariance matrix of the brightness temperature term.

[0107] Background brightness temperature Represented as:

[0108] ;

[0109] In the formula: It is the surface term state vector, which includes the surface temperature analysis field, the atmospheric temperature and humidity profile, and the surface emissivity at the current iteration number. .

[0110] Step 2.3.2, Iterative Optimization: The Levenberg-Marquardt algorithm is used to minimize the one-dimensional variational inversion cost function. The sensitivity of surface emissivity and surface temperature to satellite-observed brightness temperature is calculated using the Jacobian matrix, and the surface emissivity is iteratively updated. and surface temperature This process continues until convergence. Specifically, it includes the following steps:

[0111] Using a given surface emissivity background field and the initial guess of the surface emissivity field (definition Then the initial time The initial background brightness temperature is simulated using observation operators. If the error at this point meets the convergence condition, no further iteration is performed; otherwise, the iteration loop is entered. In the iteration loop, the surface emissivity error of the m-th iteration is calculated. Then, with the surface emissivity background field Add them together to calculate the surface emissivity in the m-th iteration. The background brightness temperature was calculated using the observation operator. If the surface emissivity error is at this time If the convergence condition has been met, the iteration ends, and the surface emissivity of the m-th iteration is output. If the surface emissivity is used as the optimal value, then continue to the (m+1)th iteration.

[0112] In the iterative solution of the one-dimensional variational inversion cost function, the convergence degree is used when determining the convergence condition. test:

[0113] (11)

[0114] Where N represents the number of infrared window channels used in the iterative calculation (here N=3). It is a satellite observation of brightness temperature; is the background brightness temperature calculated by the observation operator in the m-th iteration; E is the observation error covariance matrix of the brightness temperature term; the superscript T denotes matrix transpose. When convergence... The convergence condition is considered to have been met.

[0115] Step 3: Establish an infrared window channel brightness temperature simulation scheme to verify the ability of the optimized surface emissivity value to simulate the satellite observation brightness temperature of the infrared window channel.

[0116] The conventional atmospheric elements such as temperature, water vapor, and air pressure from the spatial grid points output by the numerical weather prediction model, along with the optimized values ​​of surface temperature and surface emissivity output by the one-dimensional variational inversion system, are used to calculate the surface emissivity. Interpolated to the satellite observation points of FY-4A AGRI, and used as input to the radiative transfer model, the satellite-observed brightness temperatures of channels 11, 12, and 13 of FY-4A AGRI are simulated to output the corresponding simulated brightness temperatures. The root mean square error (RMSE) of the satellite-observed brightness temperature and the simulated brightness temperature is calculated. It can be seen that the improvement of this simulation scheme lies in that the surface emissivity input to the radiative transfer model is derived from the constraint on the surface emissivity background field considered in step two, which takes into account the assimilation of measured surface temperature and satellite-observed brightness temperature, rather than the historical surface emissivity dataset.

[0117] This method integrates measured land surface temperatures from ground meteorological stations using an optimized three-dimensional variational assimilation cost function to obtain the initial land surface temperature field required for inversion by a one-dimensional variational inversion system. This significantly reduces the nonlinear impact of land surface temperature errors on land surface emissivity inversion. By combining the 1DVAR inversion algorithm (i.e., a one-dimensional variational inversion system) with the inverted land surface temperature and emissivity, the distortion problem of traditional decoupling methods is solved. This reduces the root mean square error (RMSE) between satellite-observed brightness temperature and simulated brightness temperature in the FY-4A / AGRI infrared window channel by 12%, making it suitable for numerical weather prediction and radiation characteristic studies in complex land surface areas.

[0118] Figure 1 The flowchart of the geostationary satellite surface emissivity co-inversion method based on background field reconstruction described in this invention is presented, mainly including three parts:

[0119] (1) The measured surface temperature from the ground meteorological station is assimilated by the three-dimensional variational assimilation cost function to obtain a more realistic surface temperature analysis field and an optimized surface temperature analysis field; the grid surface temperature in the surface temperature analysis field is interpolated to the satellite observation point.

[0120] (2) The optimized surface temperature analysis field is used as the initial guess field of surface temperature in the one-dimensional variational inversion system. The surface emissivity background field and the surface term background error covariance matrix that can characterize the coupling relationship between surface emissivity and surface temperature are constructed using the historical dataset of surface elements. The one-dimensional variational inversion cost function is established to independently invert surface emissivity or couple invert surface emissivity and surface temperature.

[0121] (3) Interpolate the gridded atmospheric temperature, humidity, and air pressure elements and the surface temperature and surface air pressure elements output by the numerical forecast model to the satellite observation points and input them into the radiative transfer model. Then input the surface emissivity obtained by the one-dimensional variational inversion cost function into the radiative transfer model to simulate and calculate the satellite observation brightness temperature of the infrared window channel and output the simulated brightness temperature of the corresponding satellite observation points.

[0122] Based on clear-sky observations by FY-4A AGRI at 12:00 UTC on October 15, 2023, the spatial distributions of surface temperatures at satellite observation points were obtained by interpolation between surface temperature Ts_Anl (the surface temperature after assimilation processing by the three-dimensional variational assimilation system described in this invention) and surface temperature TS_Bkg (the surface temperature from the Global Forecast System (GFS) analysis field). Comparison shows that surface temperatures Ts_Anl and TS_Bkg are spatially consistent, but there are differences of 2-4 K between the two surface temperatures under different surface types. Specifically, surface temperature Ts_Anl is higher in most parts of South China and lower in parts of North China. Furthermore, the surface temperature Ts_Anl, after evaluation, is closer to the measured surface temperature, indicating that Ts_Anl is a more reasonable initial guess for the surface temperature in the one-dimensional variational inversion system.

[0123] Convergence based on different inversion schemes Spatial distribution comparison shows that for the inversion scheme (inversion scheme 1) that uses the optimized surface temperature (the optimized surface temperature refers to the surface temperature analysis field obtained after assimilation processing by the three-dimensional variational assimilation system described in this invention, denoted as surface temperature Ts_Anl) as the initial surface temperature guess field of the one-dimensional variational inversion system and only inverts the surface emissivity, the convergence obtained is... The pixel ratio is 83.1%. For the inversion scheme (inversion scheme 2) that uses the surface temperature (denoted as surface temperature TS_Bkg) from the Global Forecast System (GFS) analysis field as the initial guess of the surface temperature in the one-dimensional variational inversion system and only inverts the surface emissivity, the convergence obtained is... The pixel ratio is 75.2%. Therefore, inversion scheme 1 is superior to inversion scheme 2. Furthermore, we also found that for convergence... The regions where surface temperatures differ significantly between Ts_Anl and TS_Bkg are mainly located. This indicates that the accuracy of surface temperature has a significant impact on surface emissivity inversion, and the inversion results of surface temperature Ts_Anl as the initial guess field of the one-dimensional variational inversion system are better than those of surface temperature TS_Bkg as the initial guess field of the one-dimensional variational inversion system.

[0124] Considering the coupling relationship between surface emissivity and surface temperature, and simultaneously retrieving both surface emissivity and surface temperature, using surface temperature Ts_Anl as the initial guess field for the one-dimensional variational inversion system, the convergence rate can reach 99%. This indicates that surface temperature, as a physical quantity mutually constrained with surface emissivity, can be improved by using a variational inversion algorithm that considers the coupled physical constraints between surface temperature and surface emissivity, allowing the one-dimensional variational inversion system to simultaneously invert and output surface emissivity and surface temperature, thus enhancing the inversion accuracy of surface emissivity.

[0125] Based on the comparison of the spatial distribution of the surface emissivity and the surface emissivity dataset CAMEL (The Combined ASTER and MODISEmissivity database over Land, a monthly average surface emissivity product developed by integrating MODIS and the US ASTER surface emissivity dataset, with a spatial resolution of 5 k) of different infrared window channels (including channels 11, 12 and 13) of the FY-4A AGRI obtained by the inversion method of the present invention, it can be seen that the surface emissivity and the surface emissivity dataset CAMEL obtained by the inversion method of the present invention have a high degree of spatial consistency.

[0126] Figure 2The figure presents a root mean square error (RMSE) (in K) test chart comparing the satellite-observed brightness temperatures of different infrared window channels (including channels 11, 12, and 13) of FY-4A AGRI under clear-sky conditions with the simulated brightness temperatures of the corresponding infrared window channels obtained by using different surface emissivity datasets as inputs to the radiative transfer mode. The different surface emissivity datasets used in the figure include: surface emissivity obtained based on the inversion method of this invention (corresponding to 1DVAR retrieval), surface emissivity released by the UW team (corresponding to UW), the surface emissivity dataset CAMEL (corresponding to CAMEL), and surface emissivity from the background field dataset CAMEL Climate (corresponding to CAMEL Climate). Figure 2 (a) shows the root mean square error (RMSE) test plot of the satellite observation brightness temperature of channel 11 of FY-4A AGRI under clear sky conditions and the simulated brightness temperature of the corresponding infrared window channel obtained by using different surface emissivity datasets as input for the radiative transfer mode. Figure 2 (b) shows the root mean square error (RMSE) test plot of the satellite observation brightness temperature of channel 12 of FY-4A AGRI under clear sky conditions and the simulated brightness temperature of the corresponding infrared window channel obtained by using different surface emissivity datasets as input for the radiative transfer mode. Figure 2 (c) shows the root mean square error (RMSE) (in K) comparison between the satellite-observed brightness temperature of channel 13 of FY-4A AGRI under clear sky conditions and the simulated brightness temperature of the corresponding infrared window channel obtained by using different surface emissivity datasets as inputs to the radiative transfer mode. It can be seen that the simulated brightness temperature of the corresponding infrared window channel obtained by using the surface emissivity (1DVAR retrieval) retrieved by this method as input to the radiative transfer mode has a better RMSE than the simulation results of other surface emissivity datasets. Compared with the simulation results of other surface emissivity datasets, the RMSE results of channels 12 and 13 are reduced by 12.2%. This indicates that using the optimized surface temperature value (output of the three-dimensional variational assimilation system) as the initial guess field of the surface temperature in the one-dimensional variational inversion system improves the inversion accuracy of surface emissivity, thereby improving the simulation accuracy of the satellite-observed brightness temperature of the infrared window channel of the geostationary satellite.

[0127] 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 and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A static satellite surface emissivity collaborative inversion method based on background field reconstruction, characterized in that, The method comprises the following steps: respectively acquiring measured surface temperature, and conventional surface temperature background field and atmospheric temperature and humidity profile output by a numerical prediction model; for the conventional surface temperature background field output by the numerical prediction model, introducing a surface temperature term to construct a background field with spatial structure coordination constraint of atmospheric temperature and surface temperature, to realize background field reconstruction and obtain an optimized surface temperature background field; based on the measured surface temperature and the optimized surface temperature background field, constructing an optimized three-dimensional variational assimilation cost function to assimilate the measured surface temperature and output a surface temperature analysis field; acquiring surface element in a specified historical period to construct a surface element historical data set; the surface element includes surface temperature and surface emissivity; based on the surface element historical data set, constructing a surface term background error covariance matrix representing the coordination constraint relationship between surface emissivity and surface temperature, and based on the obtained surface term background error covariance matrix, constructing a one-dimensional variational inversion cost function; acquiring satellite observation brightness temperature of a clear sky observation pixel observed by an infrared window channel of a geostationary radiation imager; generating a surface emissivity background field required by the one-dimensional variational inversion cost function from the surface emissivity in the surface element historical data set, and taking the surface temperature analysis field output by the optimized three-dimensional variational assimilation cost function as a surface temperature initial guess field required by the one-dimensional variational inversion cost function, combining the atmospheric temperature and humidity profile and the satellite observation brightness temperature of the clear sky observation pixel, to realize independent inversion of the surface emissivity or joint inversion of the surface emissivity and the surface temperature, and further output an optimized value of the surface emissivity; One-dimensional variational retrieval cost function is represented as: ； where: is the surface state vector; is the surface background field vector; is the surface background error covariance matrix; is the satellite observed brightness temperature; is the background brightness temperature output by the observation operator H at each iteration; E is a brightness temperature term observation error covariance matrix; surface element background error covariance matrix is represented as: ； In the formula: is a correlation coefficient of the surface emissivity and the surface temperature obtained by sample statistical calculation based on the historical data set of the surface element, is an error variance of the surface emissivity obtained by sample statistical calculation based on the historical data set of the surface element, is an error variance of the surface temperature obtained by sample statistical calculation based on the historical data set of the surface element. Background luminance is expressed as: ； wherein: is the surface state vector, including the surface temperature analysis field , the atmospheric temperature and humidity profiles, and the surface emissivity at the current iteration number ; In the iteration of the one-dimensional variational retrieval cost function, the Levenberg-Marquardt algorithm is used to minimize the one-dimensional variational retrieval cost function, and the sensitivity of the surface emissivity and the surface temperature to the satellite observed brightness temperature is calculated by the Jacobian matrix to iteratively update the surface emissivity until convergence. , until convergence.

2. The method of claim 1, wherein the one-dimensional. Variational inverse cost function iterative solution of the cost function, comprising the following steps: Using a given surface emissivity background field and a surface emissivity guess field , a surface temperature guess field, the initial background brightness temperature is simulated by an observation operator ; at the initial time, define , the surface emissivity error at the initial time ; The surface emissivity error for the mth iteration is calculated in an iterative loop and added to the surface emissivity background field to calculate the surface emissivity for the mth iteration ; The background brightness temperature at this time is calculated by an observation operator ; If the surface emissivity error at this time If the convergence condition has been satisfied, the iteration is exited and the surface emissivity at the mth iteration is output If the convergence condition has been satisfied, the iteration is exited and the surface emissivity at the mth iteration is output In the process of solving the one-dimensional variational inversion cost function iteratively, the convergence condition is determined using the convergence degree Test: ； Where N represents the number of infrared window channels used in the iterative calculation; is the satellite observed brightness temperature; is the background brightness temperature calculated by the observation operator in the mth iteration; E is the observation error covariance matrix of the brightness temperature term; and the superscript T represents matrix transposition.

3. The method of claim 2, wherein, Extended background error covariance matrix is represented as follows: ； where B represents the atmospheric state vector flux function , the velocity potential function , the atmospheric temperature , the relative humidity , the air pressure , and the surface temperature correlated background error covariance matrix; is the unit matrix; represents the correlation between the velocity potential function and the atmospheric state vector flux function , represents the correlation between the atmospheric temperature and the atmospheric state vector flux function , represents the correlation between the air pressure and the atmospheric state vector flux function , is the correlation between the surface temperature and the non-equilibrium part of the atmospheric temperature ; , and represent the non-equilibrium part of the velocity potential function, the atmospheric temperature, and the air pressure, respectively, is the non-equilibrium part of the surface temperature; Surface temperature Correlation between the temperature of the non-equilibrium portion of the atmosphere and the surface temperature is expressed as: ； ； In the above formula: Indicates coordinate location The non-equilibrium portion of the Earth's surface temperature; Indicates coordinate location The equilibrium portion of the Earth's surface temperature; Represents the horizontal coordinate; Represents the vertical coordinate; Indicates coordinate location The corresponding region is divided into partition numbers based on preset latitude intervals; Indicates coordinate location In partition number index The following partition; Indicates the total number of partitions; Indicates surface temperature Atmospheric temperature of the non-equilibrium part The linear correlation coefficient is calculated using the following formula: ； In the above formulae: represents the surface temperature and the temperature of the atmosphere in the non-equilibrium part covariance matrix obtained from historical statistics; represents the surface temperature variance matrix obtained from historical statistics.

4. The method of claim 3, wherein, Optimized three-dimensional variational assimilation cost function is represented as: ； ； ； ； wherein: is a background field penalty term; is a surface temperature observation penalty term; is an observation penalty term other than surface temperature. x is a state vector containing surface temperature analysis values x is a background vector containing surface temperature forecast values b x is a background vector containing surface temperature forecast values x is a background vector containing surface temperature forecast values x is a background vector containing surface temperature forecast values y is the other observation vector excluding the measured surface temperature H is the observation operator; is the observation error covariance matrix, and the superscript T denotes matrix transposition. is the measured surface temperature.

5. The method of claim 4, wherein, Inversion of a variance matrix derived from historical statistics by Cholesky decomposition of the non-equilibrium part of the atmosphere The gradient of the three-dimensional variational assimilation cost function is obtained using the conjugate gradient descent method. Find the gradient of the gradient of the 3D variational assimilation cost function during gradient descent. The numerical value and the gradient of the initial three-dimensional variational assimilation cost function The algorithm converges when the ratio is less than a preset value; the corresponding state vector at convergence. This is the optimal surface temperature analysis field after minimizing the three-dimensional variational assimilation cost function, which includes the analysis values ​​of atmospheric elements and surface temperature. Atmospheric elements include air temperature, humidity, and air pressure.

6. The method of claim 1, wherein, based on the optimized value of the surface emissivity obtained by the one-dimensional variational inversion cost function, first simulating the satellite observation point by a radiation transfer model to calculate the simulated brightness temperature of different infrared window channels of the geostationary radiation imager, and then calculating the root mean square error of the simulated brightness temperature of the corresponding infrared window channel and the satellite observation brightness temperature, to verify the simulation capability of the optimized value of the surface emissivity on the brightness temperature observation data of the infrared window channel.

7. The method according to claim 6, wherein, When the simulated brightness temperature of different infrared window channels of the geostationary radiation imager is calculated by the radiation transfer model, the conventional atmospheric element and surface element of the spatial grid output by the numerical prediction model and the optimized value of the surface emissivity obtained by inversion are interpolated to the satellite observation point of the geostationary radiation imager as the input of the radiation transfer model, to obtain the simulated brightness temperature of the infrared window channel of the geostationary radiation imager.

8. The method according to claim 7, wherein, The infrared window channel of the geostationary radiation imager includes channels 11, 12 and 13.

Images