Satellite technology-based surface temperature and surface emissivity synchronous inversion method

Abstract

The invention discloses a satellite technology-based surface temperature and surface emissivity synchronous inversion method. The method comprises the following steps of: acquiring adjacent observation data; atmospheric parameter initial values are calculated; performing estimation based on the atmospheric parameter initial value and an I-DN algorithm to obtain a first surface temperature, a first surface emissivity, a temperature profile correction value and a humidity profile correction coefficient; the initial atmospheric temperature profile and the initial humidity profile are corrected, and corrected atmospheric parameters are calculated through an RTTOV model; estimating based on the first surface temperature, the first surface emissivity, the corrected atmospheric parameters and an I-TES algorithm to obtain a second surface temperature and a second surface emissivity; calculating weighted surface emissivity; estimating based on the weighted surface emissivity and an I-SW algorithm to obtain a third surface temperature; calculating a weighted surface temperature; and completing integrated inversion of the surface temperature and the surface emissivity of the target area. According to the method, the accuracy and the stability of synchronous inversion of the surface temperature and the emissivity are improved.

Description

Technical Field

[0001] The present invention relates to the field of communication technology, and in particular to a method for synchronously inverting surface temperature and surface emissivity based on satellite technology. Background Art

[0002] In the inversion of thermal infrared remote sensing surface temperature and surface emissivity, the main technical methods currently include single-channel algorithm, temperature and emissivity separation (TES) algorithm and split window (SW) algorithm. Among them, the single-channel algorithm is more used in remote sensing data with only a single thermal infrared band because it is highly dependent on the accuracy of atmospheric correction. The TES algorithm is essentially a synchronous inversion algorithm for surface emissivity and surface temperature, which is suitable for remote sensing data with three or more thermal infrared bands. However, on the one hand, it requires high-precision atmospheric correction, and on the other hand, its algorithm implementation process relies on some empirical formulas, which makes it more sensitive to atmospheric errors, and the accuracy will be significantly reduced in low-emissivity surface areas. Compared with the single-channel and TES algorithms, although the SW algorithm does not require atmospheric correction, its algorithm accuracy depends on the prior surface emissivity accuracy, especially in arid and semi-arid surfaces, the accuracy will be significantly reduced.

[0003] In summary, the existing thermal infrared remote sensing surface temperature and emissivity inversion methods are generally based on certain assumptions or require some prior knowledge to reduce the number of unknowns in the radiation transfer equation, thereby solving the ill-conditioned solution problem in the inversion process of surface temperature and surface emissivity. However, this approach does not solve the coupling problem between parameters well, resulting in unstable inversion accuracy. Summary of the invention

[0004] The invention provides a method for synchronous inversion of surface temperature and surface emissivity based on satellite technology, which is used to improve the accuracy and stability of synchronous inversion of surface temperature and emissivity.

[0005] In a first aspect, the present invention provides a method for synchronously inverting land surface temperature and land surface emissivity based on satellite technology, comprising:

[0006] Step 1: Collect nearby observation data of the target area based on thermal infrared remote sensing acquisition equipment;

[0007] Step 2: Calculate the initial values ​​of atmospheric parameters based on the RTTOV model and nearby observation data; the initial values ​​of atmospheric parameters include the atmospheric downwelling radiance Atmospheric upward radiance and atmospheric transmittance τ t ;

[0008] Step 3: Estimate based on the initial values ​​of atmospheric parameters and the I-DN algorithm to obtain the first surface temperature First surface emissivity ε 1, Temperature profile correction value T a,t and humidity profile correction factor w t ; where t represents the observation time;

[0009] Step 4: Correct T based on temperature profile a,t and humidity profile correction factor w t For the initial atmospheric temperature profile T prof，t and the initial humidity profile Q prof，t Correction is performed and the corrected atmospheric parameters are calculated by the RTTOV model;

[0010] Step 5: Based on the first surface temperature First surface emissivity ε 1 , the corrected atmospheric parameters and the I-TES algorithm are used to estimate the second surface temperature and the second surface emissivity ε 2 ;

[0011] Step 6: Based on the first surface emissivity ε 1 Second surface emissivity ε 2 , the weighted surface emissivity ε is calculated w ;

[0012] Step 7: Based on the weighted surface emissivity ε w The third surface temperature is estimated by using the I-SW algorithm

[0013] Step 8: Based on the first surface temperature Second surface temperature The third surface temperature Calculate weighted land surface temperature

[0014] Step 9: Based on weighted surface temperature and the surface emissivity ε w , calculate the cost function and repeat steps 3-8 until the function threshold is met or the maximum number of iterations is reached;

[0015] Step 10: Repeat steps 3-9 to complete the inversion of the surface temperature and surface emissivity of the target area;

[0016] Step 11: For the uninverted areas that do not meet the observation number greater than 2, first use the interpolation method to obtain the corresponding atmospheric correction parameters, and then use the I-TES algorithm to calculate the corresponding surface temperature and emissivity to complete the regional integrated inversion.

[0017] According to the satellite-based synchronous inversion method for land surface temperature and land surface emissivity provided by the present invention, the weighted land surface temperature The formula is:

[0018]

[0019] Among them, θ j Represents the surface emissivity weight coefficient, which can be determined by statistical regression or a small neural network. The default value is: θ j =[0.5,0.3,0.2].

[0020] According to the synchronous inversion method of surface temperature and surface emissivity based on satellite technology provided by the present invention, the cost function is:

[0021]

[0022]

[0023] Where C represents the cost function, N represents the number of thermal infrared bands, i represents the i-th thermal infrared band, and L i,t represents the on-board radiance observed in the i-th thermal infrared band, L(i,t) represents the on-board radiance estimated in the i-th thermal infrared band, B represents the Planck function, represents the weighted surface temperature, represents the weighted surface emissivity of the i-th thermal infrared band, represents the corrected atmospheric downward radiation in the ith thermal infrared band, represents the corrected atmospheric upwelling radiation of the ith thermal infrared band, and τ′i,t represents the corrected atmospheric transmittance of the ith thermal infrared band.

[0024] According to the satellite-based synchronous inversion method for surface temperature and surface emissivity provided by the present invention, the first surface temperature is obtained by estimating based on the initial values ​​of atmospheric parameters and the I-DN algorithm. First surface emissivity ε 1 , Temperature profile correction value T a ,t and humidity profile correction factor w t ,include:

[0025] Step 3-1: Calculate the initialization vector x0. The expression of the initialization vector x0 is as follows:

[0026] x0=[ΔT a，t , w t , T s，t , ε i ], t∈[1,D],i∈[1,N];

[0027] t represents the tth observation; D represents the number of observations, and D≥2; i represents the i-th thermal infrared band; N represents the number of thermal infrared bands; ΔT a,trepresents the correction value of the atmospheric temperature profile of the tth observation; wt represents the correction coefficient of the atmospheric humidity profile of the tth observation; T s,t is the surface temperature observed at the tth time, and the near-surface air temperature is used as the initial value by default; ε i represents the surface emissivity of the i-th thermal infrared band, which can be obtained by the following formula:

[0028]

[0029] T b_max，t =max(B -1 (L i，t ));

[0030] Among them, B i represents the Planck function, L i，t is the directly observed on-star radiance, represents the atmospheric downwelling radiance, represents the atmospheric upward radiance, τ i，t represents the atmospheric transmittance, represents the radiance above the ground, B -1 represents the inverse function of Planck function;

[0031] Step 3-2: Substitute the surface emissivity ε i Substituting into the following formula, we get the characteristic vector e of the surface emissivity;

[0032]

[0033] Among them, e i represents the i-th feature; f k represents the prediction function, which can be determined by statistical regression, machine learning, deep learning and other algorithms; k represents the kth prediction function, and K represents the number of prediction functions;

[0034] Step 3-3: Replace the surface emissivity ε in the initialization vector x0 with the eigenvector e i , get the target vector x;

[0035] x=[ΔT a，t , w t , T s，t , e];

[0036] Step 3-4: By ΔT a，t and w t Modify the atmospheric profile and reconstruct the surface emissivity through the eigenvector e

[0037] Step 3-5: Set the modified atmospheric parameters and T s，t and Input the RTTOV model to estimate the on-board radiance, denoted as Li，t (x);

[0038] Step 3-6: Calculate the corrected value of x. The specific steps are as follows:

[0039] Construct the correction matrix Δx:

[0040]

[0041] Where Δx represents the correction matrix; I n represents the unit matrix with dimension n; n is the length of the target vector; γ is the correction factor, whose value is set to 1e -6 ;

[0042] Based on the correction matrix and RTTOV model, calculate the correction error ΔL of the on-board radiance i，t :

[0043] ΔL i，t =L i，t (x+Δx)-L i，t ;

[0044] Where i represents the i-th thermal infrared band, x represents the target vector, Δx represents the correction matrix, and L i，t (x+Δx) represents the estimated corrected on-board radiance, L i，t It represents the on-board radiance of thermal infrared remote sensing observation;

[0045] The x correction value δ is calculated by the following formula:

[0046]

[0047] Among them, pinv represents the pseudo-inverse function, γ represents the radiance scaling factor, when the wave number unit is used, the default value is 1000, when the wavelength unit is used, the default value is 1;

[0048] Step 3-7: Correct the target vector x:

[0049] x=x+α·δ·x

[0050] Among them, α is the correction weight coefficient, which is determined by the observation conditions:

[0051]

[0052] a Ta +a w +a Ts +a ε =1;

[0053] Among them, a Ta 、a w 、a Ts and a εThey represent the correction weights of temperature profile, humidity profile, surface temperature and surface emissivity, respectively; i represents the i-th thermal infrared band, c Ta,1 、c w,1 、c Ts,1 、c ε,1 Respectively represent a Ta 、a w 、a Ts 、a ε The weight coefficient can be calculated by numerical fitting or machine learning model algorithms. W is the water vapor content of the atmospheric profile, Ta represents the atmospheric temperature, and T b,i represents the equivalent brightness temperature on the satellite in the i-th band; N represents the number of thermal infrared bands;

[0054] Step 3-8: Execute steps 3-1 to 3-7 repeatedly. When the norm of δ is less than the specified threshold or reaches the predetermined number of iterations, the estimation of the surface temperature and surface emissivity is completed and recorded as the first surface temperature. and the first surface emissivity ε 1 .

[0055] According to the satellite-based synchronous inversion method for surface temperature and surface emissivity provided by the present invention, the specific steps of atmospheric profile correction in step 4 are as follows:

[0056] According to the current atmospheric temperature profile, starting from the ground height, the temperature offset is added layer by layer in the atmospheric temperature profile. The specific formula is as follows:

[0057]

[0058] Among them, T′ a,l represents the corrected temperature of the first layer of atmosphere; T a,0 represents the near-surface (0th layer) atmospheric temperature; T a,l represents the initial temperature of the first layer of atmosphere; ΔT a represents the temperature offset; nl represents the number of vertical layers of the atmospheric profile;

[0059] According to the current atmospheric humidity profile, the atmospheric humidity profile is scaled as a whole. The calculation formula is as follows:

[0060] Q′ l =Q l w,l∈[0,nl-1];

[0061] Among them, Q′ l represents the corrected humidity of the lth layer of atmosphere; Q represents the initial humidity of the lth layer of atmosphere; w represents the scaling factor of atmospheric water vapor; nl represents the number of vertical layers of the atmospheric profile.

[0062] According to the synchronous inversion method of surface temperature and surface emissivity based on satellite technology provided by the present invention, the I-TES algorithm in step 5 is an optimization algorithm based on the ASTERTES algorithm, and the implementation steps are as follows:

[0063] Step 5-1: From the surface emissivity ε 1 Extract the maximum value from

[0064] Step 5-2: Calculate the initial value of the surface temperature using the following formula:

[0065]

[0066] Where i represents the i-th thermal infrared band, B represents the Planck function, T i,t represents the equivalent surface temperature, B -1 represents the inverse Planck function, It represents the band's ground radiance. represents the corrected atmospheric downward radiation;

[0067] Step 5-3: Using the initial surface temperature Calculate the corresponding surface emissivity and its normalized spectrum The formula is as follows:

[0068]

[0069] Step 5-4: Calculate the minimum surface emissivity using the β-MMD hybrid model

[0070]

[0071] where f k represents the prediction function, which can be determined by algorithms such as statistical regression or machine learning, k represents the kth prediction function, and K represents the number of prediction functions. and represent the maximum and minimum values ​​of the normalized spectral emissivity, respectively;

[0072] Step 5-5: Reconstruct the surface emissivity spectrum using the following formula:

[0073]

[0074] Step 5-6: Exploitation Re-estimate the surface temperature Denoted as the second surface emissivity and the second surface temperature respectively;

[0075] Step 5-7: Loop through steps 5-1 to 5-6 until the cost function changes between the previous and next two times are less than the specified threshold. The cost function used is as follows:

[0076]

[0077] i represents the i-th thermal infrared band, t represents the t-th observation, L i,t is the observed on-star radiance, represents the estimated on-star radiance, N represents the number of thermal infrared bands, and D represents the number of observations.

[0078] According to the satellite-based synchronous inversion method for surface temperature and surface emissivity provided by the present invention, based on the first surface emissivity ε 1 and the second surface emissivity ε 2 Calculate and get the weighted surface emissivity ε w ,include:

[0079] The first surface emissivity ε 1 and the second surface emission ε 2 Substituting into the following formula, the weighted surface emissivity ε is calculated w ;

[0080]

[0081] Among them, φ j is the surface emissivity weight coefficient, which can be determined by statistical regression or a small neural network. The default value is: j =0.5.

[0082] According to the satellite-based synchronous inversion method for surface temperature and surface emissivity provided by the present invention, based on the weighted surface emissivity ε w The third surface temperature is estimated by using the I-SW algorithm The steps include:

[0083]

[0084] Among them, f k represents the prediction function, which can be determined by statistical regression or machine learning algorithms, k represents the kth prediction function, K represents the number of prediction functions, i represents the i-th thermal infrared band, T b,i,t is the equivalent stellar brightness temperature, Represents the surface emissivity.

[0085] In a second aspect, the present invention further provides a device for synchronously inverting land surface temperature and land surface emissivity based on satellite technology, comprising:

[0086] An acquisition module is used to acquire nearby observation data of a target area based on a thermal infrared remote sensing acquisition device;

[0087] The first calculation module is used to calculate the initial values ​​of atmospheric parameters based on the RTTOV model and nearby observation data; the initial values ​​of atmospheric parameters include the atmospheric downwelling radiance Atmospheric upward radiance and atmospheric transmittance τ t ;

[0088] The first estimation module is used to estimate the first surface temperature based on the initial value of the atmospheric parameters and the I-DN algorithm. First surface emissivity ε 1 , Temperature profile correction value T a,t and humidity profile correction factor w t ; where t represents the observation time;

[0089] Correction module, used to correct the value T based on the temperature profile a,t and humidity profile correction factor w t For the initial atmospheric temperature profile T prof，t and the initial humidity profile Q prof，t Correction is performed and the corrected atmospheric parameters are calculated by the RTTOV model;

[0090] The second estimation module is used to estimate the surface temperature based on the first First surface emissivity ε 1 , the corrected atmospheric parameters and the I-TES algorithm are used to estimate the second surface temperature and the second surface emissivity ε 2 ;

[0091] The second calculation module is used to calculate the surface emissivity ε based on the first surface emissivity ε 1 Second surface emissivity ε 2 , the weighted surface emissivity ε is calculated w ;

[0092] The third estimation module is used to estimate the surface emissivity based on the weighted surface emissivity ε w The third surface temperature is estimated by using the I-SW algorithm

[0093] The third calculation module is used to calculate the temperature of the first surface of the Second surface temperature The third surface temperature Calculate weighted land surface temperature

[0094] Synchronous inversion module for:

[0095] Based on weighted surface temperature and the surface emissivity ε w , calculate the cost function and repeat steps 3-8 until the function threshold is met or the maximum number of iterations is reached;

[0096] Repeat the above steps to complete the inversion of surface temperature and surface emissivity in the target area;

[0097] For the uninverted areas that do not meet the observation number greater than 2, the interpolation method is first used to obtain the corresponding atmospheric correction parameters, and then the I-TES algorithm is used to calculate the corresponding surface temperature and emissivity to complete the regional integrated inversion.

[0098] In a third aspect, the present invention further provides an electronic device, comprising: a memory for storing a computer software program; a processor for reading and executing the computer software program, thereby implementing the method for synchronous inversion of surface temperature and surface emissivity based on satellite technology as in the first aspect mentioned above.

[0099] In a fourth aspect, the present invention further provides a non-transitory computer-readable storage medium, in which a computer software program is stored. When the computer software program is executed by a processor, the method for synchronous inversion of surface temperature and surface emissivity based on satellite technology as described in the first aspect above is implemented.

[0100] In a fifth aspect, the present invention further provides a computer program product, including a computer program, which, when executed by a processor, implements the method for synchronous inversion of surface temperature and surface emissivity based on satellite technology as described in the first aspect above.

[0101] The satellite-based synchronous inversion method for land surface temperature and land surface emissivity provided by the present invention realizes effective decoupling in the process of solving land-atmosphere parameters based on a multi-algorithm integrated inversion model based on uncoupled error optimization. At the same time, by weighting the land surface temperature and land surface emissivity inverted by multiple algorithms, the differences in the principle error levels of different algorithms are effectively balanced, thereby improving the accuracy and stability of the synchronous inversion of land surface temperature and emissivity. BRIEF DESCRIPTION OF THE DRAWINGS

[0102] Figure 1 It is a schematic flow chart of a method for synchronously inverting land surface temperature and land surface emissivity based on satellite technology provided by the present invention;

[0103] Figure 2 It is a schematic structural diagram of a device for synchronously inverting ground surface temperature and ground surface emissivity based on satellite technology provided by the present invention;

[0104] Figure 3 A schematic diagram of an electronic device according to an embodiment of the present invention;

[0105] Figure 4A schematic diagram of an embodiment of a computer-readable storage medium provided for an embodiment of the present invention. DETAILED DESCRIPTION

[0106] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative work are within the scope of protection of the present invention.

[0107] In the description of the present invention, the terms "first" and "second" are used for descriptive purposes only and should not be understood as indicating or implying relative importance or implicitly indicating the number of technical features indicated. Therefore, the features defined as "first" and "second" may explicitly or implicitly include one or more features. In the description of the present invention, the meaning of "plurality" is two or more, unless otherwise clearly and specifically defined.

[0108] In the description of the present invention, the term "for example" is used to mean "used as an example, illustration or explanation". Any embodiment described as "for example" in the present invention is not necessarily to be construed as being more preferred or advantageous than other embodiments. The following description is given to enable any technician in the field to implement and use the present invention. In the following description, details are listed for the purpose of explanation. It should be understood that a person of ordinary skill in the art can recognize that the present invention can be implemented without using these specific details. In other examples, well-known structures and processes will not be elaborated in detail to avoid obscuring the description of the present invention with unnecessary details. Therefore, the present invention is not intended to be limited to the embodiments shown, but is consistent with the widest scope consistent with the principles and features disclosed in the present invention.

[0109] Optional, see Figure 1 , Figure 1 The present invention provides a method for synchronously inverting the surface temperature and the surface emissivity based on satellite technology, which includes the following steps:

[0110] Step 1: Collect nearby observation data of the target area based on thermal infrared remote sensing acquisition equipment.

[0111] Among them, the adjacent observation data include: ① Observed radiance, denoted as L i,p,t , i represents the i-th thermal infrared band, t represents the observation time, and p represents the regional pixel. The following parameters will contain pixel information by default, that is, p, and will not be marked separately. ② Pixel longitude and latitude; ③ Atmospheric profile: (T prof,t , Q prof,t , Pprof,t , O3 prof,t );④ Ground elevation (Elevation); ⑤ Surface type (Surface Type and Water Type).

[0112] Optionally, the radiance of adjacent observations in step (1) should have nonlinear characteristics, that is, there should be significant differences in certain characteristic parameters, such as the near-surface air temperature change of the atmospheric temperature profile > 2K or the satellite-observed radiance brightness temperature at the same location > 1K, etc., to avoid linear correlation between the two observation results and the inability to solve the equations in the I-DN algorithm.

[0113] (2) Collecting a priori surface emissivity, which can be either a priori data product or obtained through calculation, such as estimation based on the NDVI index.

[0114] Step 2: Calculate the initial values ​​of atmospheric parameters based on the RTTOV model and nearby observation data.

[0115] Among them, the RTTOV (Radiative Transfer for TOVS) model is a tool for simulating atmospheric radiation transfer. The initial values ​​of atmospheric parameters include the atmospheric downwelling radiance Atmospheric upward radiance and atmospheric transmittance τ t Therefore, it can be understood that the initial values ​​of atmospheric parameters under the current atmospheric profile, surface type and ground elevation conditions are calculated through the RTTOV model.

[0116] Step 3: Estimate based on the initial values ​​of atmospheric parameters and the I-DN algorithm to obtain the first surface temperature First surface emissivity ε 1 , Temperature profile correction value T a,t and humidity profile correction factor.

[0117] Where t represents the observation time.

[0118] It can be understood that, assuming that the surface emissivity remains unchanged in adjacent observations, the initial values ​​of atmospheric parameters and the I-DN algorithm are used to estimate the first surface temperature First surface emissivity ε 1 , Temperature profile correction value T a ,t and humidity profile correction factor w t .

[0119] In one embodiment, the first surface temperature is obtained by estimating based on the initial value of the atmospheric parameter and the I-DN algorithm in step 3. First surface emissivity ε 1 , Temperature profile correction value T a,t and humidity profile correction factor wt ,include:

[0120] Step 3-1: Calculate the initialization vector x0. The expression of the initialization vector x0 is as follows:

[0121] x0=[ΔT a，t , w t , T s,t , ε i ], t∈[1,D],i∈[1,N];

[0122] t represents the tth observation; D represents the number of observations, and D≥2; i represents the i-th thermal infrared band; N represents the number of thermal infrared bands; ΔT a,t represents the corrected value of the atmospheric temperature profile at the tth observation; w t represents the atmospheric humidity profile correction coefficient of the tth observation; T s,t is the surface temperature observed at the tth time, and the near-surface air temperature is used as the initial value by default; ε i represents the surface emissivity of the i-th thermal infrared band, which can be obtained by the following formula:

[0123]

[0124]

[0125] T b_max，t =max(B -1 (L i，t );

[0126] Where B represents the Planck function, L i，t is the directly observed on-star radiance, represents the atmospheric downwelling radiance, represents the atmospheric upward radiance, τ i，t represents the atmospheric transmittance, represents the radiance above the ground, B -1 represents the inverse function of Planck function;

[0127] Step 3-2: Substitute the surface emissivity ε i Substituting into the following formula, we get the characteristic vector e of the surface emissivity;

[0128]

[0129] Among them, e i represents the i-th feature; f k represents the prediction function, which can be determined by statistical regression, machine learning, deep learning and other algorithms; k represents the kth prediction function, and K represents the number of prediction functions;

[0130] Step 3-3: Replace the surface emissivity ε in the initialization vector x0 with the eigenvector e i , get the target vector x;

[0131] x=[ΔT a，t w t , T s，t , e];

[0132] Step 3-4: By ΔT a，t and w t Modify the atmospheric profile and reconstruct the surface emissivity through the eigenvector e

[0133] Step 3-5: Set the modified atmospheric parameters and T s，t and Input the RTTOV model to estimate the on-board radiance, denoted as L i，t (x);

[0134] Step 3-6: Calculate the corrected value of x. The specific steps are as follows:

[0135] Construct the correction matrix Δx:

[0136]

[0137] Where Δx represents the correction matrix; I n represents the unit matrix with dimension n; n is the length of the target vector; γ is the correction factor, whose value is set to 1e -6 ;

[0138] Based on the correction matrix and RTTOV model, calculate the correction error ΔL of the on-board radiance i，t :

[0139] ΔL i，t =L i，t (x+Δx)-L i，t ;

[0140] Where i represents the i-th thermal infrared band, x represents the target vector, Δx represents the correction matrix, and L i，t (x+Δx) represents the estimated corrected on-board radiance, L i ,t represents the on-board radiance of thermal infrared remote sensing observation;

[0141] The x correction value δ is calculated by the following formula:

[0142]

[0143] Among them, pinv represents the pseudo-inverse function, γ represents the radiance scaling factor, when the wave number unit is used, the default value is 1000, when the wavelength unit is used, the default value is 1;

[0144] Step 3-7: Correct the target vector x:

[0145] x=x+α·δ·x

[0146] Among them, α is the correction weight coefficient, which is determined by the observation conditions:

[0147]

[0148] a Ta +a w +a Ts +a ε =1;

[0149] Among them, a Ta 、a w 、a Ts and a ε They represent the correction weights of temperature profile, humidity profile, surface temperature and surface emissivity, respectively; i represents the i-th thermal infrared band, c Ta,1 、c w,1 、c Ts,1 、c ε,1 Respectively represent a Ta 、a w 、a Ts 、a ε The weight coefficient can be calculated by numerical fitting or machine learning model algorithms. W is the water vapor content of the atmospheric profile, Ta represents the atmospheric temperature, and T b,i represents the equivalent brightness temperature on the satellite in the i-th band; N represents the number of thermal infrared bands;

[0150] Step 3-8: Execute steps 3-1 to 3-7 repeatedly. When the norm of δ is less than the specified threshold or reaches the predetermined number of iterations, the estimation of the surface temperature and surface emissivity is completed and recorded as the first surface temperature. and the first surface emissivity ε 1 .

[0151] Step 4: Correct T based on temperature profile a ,t and humidity profile correction factor w t For the initial atmospheric temperature profile T prof，t and the initial humidity profile Q prof，t Corrections are made and the corrected atmospheric parameters are calculated using the RTTOV model.

[0152] In one embodiment, the specific steps of correcting the atmospheric profile in step 4 are as follows:

[0153] According to the current atmospheric temperature profile, starting from the ground height, the temperature offset is added layer by layer in the atmospheric temperature profile. The specific formula is as follows:

[0154]

[0155] Among them, T′ a,l represents the corrected temperature of the first layer of atmosphere; T a,0 represents the near-surface (0th layer) atmospheric temperature; T a,l represents the initial temperature of the first layer of atmosphere; ΔT a represents the temperature offset; nl represents the number of vertical layers of the atmospheric profile;

[0156] According to the current atmospheric humidity profile, the atmospheric humidity profile is scaled as a whole. The calculation formula is as follows:

[0157] Q′ l =Q l w, l∈[0,nl-1];

[0158] Among them, Q′ l represents the corrected humidity of the lth layer of atmosphere; Q represents the initial humidity of the lth layer of atmosphere; w represents the scaling factor of atmospheric water vapor; nl represents the number of vertical layers of the atmospheric profile.

[0159] Step 5: Based on the first surface temperature First surface emissivity ε 1 , the corrected atmospheric parameters and the I-TES algorithm are used to estimate the second surface temperature and the second surface emissivity ε 2 .

[0160] Step 6: Based on the first surface emissivity ε 1 Second surface emissivity ε 2 , the weighted surface emissivity ε is calculated w .

[0161] Weighted surface emissivity ε w The calculation process is as follows:

[0162] The first surface emissivity ε 1 and the second surface emission ε 2 Substituting into the following formula, the weighted surface emissivity ε is calculated w ;

[0163]

[0164] Among them, φ j is the surface emissivity weight coefficient, which can be determined by statistical regression or a small neural network. The default value is: j =0.5.

[0165] Step 7: Based on the weighted surface emissivity ε wThe third surface temperature is estimated by using the I-SW algorithm

[0166] Based on the weighted surface emissivity ε w The third surface temperature is estimated by using the I-SW algorithm The steps include:

[0167]

[0168] Among them, f k represents the prediction function, which can be determined by statistical regression or machine learning algorithms, k represents the kth prediction function, K represents the number of prediction functions, i represents the i-th thermal infrared band, T b,i,t is the equivalent stellar brightness temperature, Represents the surface emissivity.

[0169] Step 8: Based on the first surface temperature Second surface temperature The third surface temperature Calculate weighted land surface temperature

[0170] Weighted surface temperature The formula is:

[0171]

[0172] Among them, θ j Represents the surface emissivity weight coefficient, which can be determined by statistical regression or a small neural network. The default value is: θ j =[0.5,0.3,0.2].

[0173] Step 9: Based on weighted surface temperature and the surface emissivity ε w , calculate the cost function, and repeat steps 3-8 until the function threshold is met or the maximum number of iterations is reached.

[0174] The cost function is:

[0175]

[0176]

[0177] Where C represents the cost function, N represents the number of thermal infrared bands, i represents the i-th thermal infrared band, and L i,t represents the on-board radiance observed in the i-th thermal infrared band, L(i,t) represents the on-board radiance estimated in the i-th thermal infrared band, B represents the Planck function, represents the weighted surface temperature of the i-th thermal infrared band, represents the weighted surface emissivity of the i-th thermal infrared band, represents the corrected atmospheric downward radiation in the ith thermal infrared band, represents the corrected atmospheric upwelling radiation of the ith thermal infrared band, τ′ i,t represents the corrected atmospheric transmittance of the ith thermal infrared band.

[0178] Step 10: Repeat steps 3-9 to complete the inversion of the surface temperature and surface emissivity of the target area.

[0179] Specifically, the weighted surface temperature set Weighted surface emissivity ε w , Temperature profile correction value T a,t and humidity profile correction factor w t As the initial value of the I-DN algorithm, repeat steps 3 to 9 until the change in the loss function is less than the specified threshold or the maximum number of iterations is reached.

[0180] Step 11: For the uninverted areas that do not meet the observation number greater than 2, first use the interpolation method to obtain the corresponding atmospheric correction parameters, and then use the I-TES algorithm to calculate the corresponding surface temperature and emissivity to complete the regional integrated inversion.

[0181] In one embodiment, the I-TES algorithm in step 5 is an optimization algorithm based on the ASTER TES algorithm, and the implementation steps are as follows:

[0182] Step 5-1: From the surface emissivity ε 1 Extract the maximum value from

[0183] Step 5-2: Calculate the initial value of the surface temperature using the following formula:

[0184]

[0185] Where i represents the i-th thermal infrared band, B represents the Planck function, T i,t represents the equivalent surface temperature, B -1 represents the inverse Planck function, It represents the band's ground radiance. represents the corrected atmospheric downward radiation;

[0186] Step 5-3: Using the initial surface temperature Calculate the corresponding surface emissivity and its normalized spectrum The formula is as follows:

[0187]

[0188] Step 5-4: Calculate the minimum surface emissivity using the β-MMD hybrid model

[0189]

[0190] where f k represents the prediction function, which can be determined by algorithms such as statistical regression or machine learning, k represents the kth prediction function, and K represents the number of prediction functions. and represent the maximum and minimum values ​​of the normalized spectral emissivity, respectively;

[0191] Step 5-5: Reconstruct the surface emissivity spectrum using the following formula:

[0192]

[0193] Step 5-6: Exploitation Re-estimate the surface temperature Denoted as the second surface emissivity and the second surface temperature respectively;

[0194] Step 5-7: Loop through steps 5-1 to 5-6 until the cost function changes between the previous and next two times are less than the specified threshold. The cost function used is as follows:

[0195]

[0196] i represents the i-th thermal infrared band, t represents the t-th observation, L i,t is the observed on-star radiance, represents the estimated on-star radiance, N represents the number of thermal infrared bands, and D represents the number of observations.

[0197] In one embodiment, the β-mmd model in step 5-4 is a model structure based on a machine learning algorithm, and its purpose is to replace the Ratio and MMD algorithms in the original TES algorithm to improve the accuracy.

[0198] The embodiment of the present invention uses the XGBoost algorithm by default, and based on the ASTER band characteristics, the model accuracy is tested using simulated data. The experimental results show that the root mean square error (RMSE) of estimating the natural surface emissivity (≥0.9) using the β-mmd model is about 0.004, and the bias is about 0.0002. The RMSE of the original TES method is about 0.016, and the bias is about -0.0047, and the accuracy is improved by about 4 times.

[0199] The main parameters of the XGBoost model are as follows:

[0200] parameter Value Maximum tree depth (max_depth) 10 Learning rate (learning_rate) 0.001 subsample 0.9 colsample_bytree 0.9

[0201] The above-mentioned β-mmd model can also use other model frameworks such as NN and LSTM, and the specific test selection can be based on the number of bands used and the specific situation of the simulated data set.

[0202] The I-SW algorithm in step 7 above is an improved algorithm based on the principle of the SW algorithm. The specific steps are as follows:

[0203] By simulating data, we build a T based on machine learning methods. b -T s Model. This paper takes the XGboost algorithm as an example and sets the main parameters of the model as follows:

[0204] parameter Value Maximum tree depth (max_depth) 6 Learning rate (learning_rate) 0.1 subsample 0.8 colsample_bytree 0.8

[0205] The input of the model is the surface emissivity characteristic e, the atmospheric water vapor content W and the radiation equivalent brightness temperature Tb; the output of the model is the surface temperature.

[0206] The verification results based on simulated data show that the RMSE of the I-SW algorithm under natural surface conditions (ε>0.9) is about 0.15K, and its RMSE is still better than 0.7K when the surface emissivity error is 0.01, the water vapor error is 10%, and the equivalent brightness temperature error is 0.2K.

[0207] Therefore, it can be understood that the simulation data set implemented in the present invention refers to the observed radiation, surface emissivity and atmospheric parameters simulated by using the radiation transfer model after considering various typical atmosphere, surface, observation angle and other factors.

[0208] In the solution of the embodiment of the present invention, the method for constructing the simulation data set involved is as follows:

[0209] (1) The radiative transfer model uses RTTOV by default;

[0210] (2) TIGR2000 and US 6M Profiles were used as the basic profiles, and a temperature perturbation of [-6K, 6K] was added to the temperature profile with a step size of 2K; the humidity profile was scaled to [0.5, 2.0] with a step size of 0.1;

[0211] (3) The observation angle VZA was gradually increased from 0° to 75° in 5° steps;

[0212] (4) The surface temperature is set according to the observation time and near-surface air temperature. Specifically:

[0213] Under daytime clear sky conditions, the value range is [Ta, Ta+25K], changing in steps of 5K, where Ta is the atmospheric near-surface temperature;

[0214] Under clear sky conditions at night, the value range is [Ta-15, Ta+5], changing in steps of 5K.

[0215] (5) The surface emissivity is calculated using ASTER LIB and UW-IR MAP, and the band equivalent emissivity is calculated using the band response function.

[0216] (6) Finally, the above parameters are arbitrarily combined and the radiation transfer model R is driven to obtain the atmospheric uplink and downlink radiation, transmittance parameters and observed radiance under different conditions, thus completing the construction of the simulation data set.

[0217] The embodiment of the present invention is based on a multi-algorithm integrated inversion model with uncoupled error optimization, which realizes effective decoupling in the process of solving the ground-air parameters. At the same time, by weighting the surface temperature and surface emissivity inverted by multiple algorithms, the difference in the principle error level of different algorithms is effectively balanced, thereby improving the accuracy and stability of the synchronous inversion of surface temperature and emissivity. At the same time, the embodiment of the present invention also constructs a β-mmd model and improves the TES algorithm of the original ASTER. The verification results based on simulation data show that the β-mmd model can significantly improve the inversion accuracy of the surface emissivity in the TES algorithm.

[0218] See also Figure 2 , Figure 2 A schematic diagram of the structure of a synchronous inversion device for surface temperature and surface emissivity based on satellite technology provided by an embodiment of the present invention, wherein the synchronous inversion device for surface temperature and surface emissivity based on satellite technology comprises:

[0219] The acquisition module 210 is used to acquire nearby observation data of the target area based on the thermal infrared remote sensing acquisition equipment;

[0220] The first calculation module 220 is used to calculate the initial values ​​of atmospheric parameters based on the RTTOV model and nearby observation data; the initial values ​​of atmospheric parameters include the atmospheric downwelling radiance Atmospheric upward radiance and atmospheric transmittance τ t ;

[0221] The first estimation module 230 is used to estimate the first surface temperature based on the initial value of the atmospheric parameter and the I-DN algorithm. First surface emissivity ε 1 , Temperature profile correction value T a,t and humidity profile correction factor w t ; where t represents the observation time;

[0222] Correction module 240, used to correct the value T based on the temperature profile a,t and humidity profile correction factor w tFor the initial atmospheric temperature profile T prof，t and the initial humidity profile Q prof，t Correction is performed and the corrected atmospheric parameters are calculated by the RTTOV model;

[0223] The second estimation module 350 is used to estimate the temperature of the first ground surface based on the First surface emissivity ε 1 , the corrected atmospheric parameters and the I-TES algorithm are used to estimate the second surface temperature and the second surface emissivity ε 2 ;

[0224] The second calculation module 260 is used to calculate the first surface emissivity ε based on the first surface emissivity ε 1 Second surface emissivity ε 2 , the weighted surface emissivity ε is calculated w ;

[0225] The third estimation module 270 is used to estimate the weighted surface emissivity ε w The third surface temperature is estimated by using the I-SW algorithm

[0226] The third calculation module 280 is used to calculate the temperature of the first ground surface based on the Second surface temperature The third surface temperature Calculate weighted land surface temperature

[0227] The synchronous inversion module 290 is used to:

[0228] Based on weighted surface temperature and the surface emissivity ε w , calculate the cost function and repeat steps 3-8 until the function threshold is met or the maximum number of iterations is reached;

[0229] Repeat the above steps to complete the inversion of surface temperature and surface emissivity in the target area;

[0230] For the uninverted areas that do not meet the observation number greater than 2, the interpolation method is first used to obtain the corresponding atmospheric correction parameters, and then the I-TES algorithm is used to calculate the corresponding surface temperature and emissivity to complete the regional integrated inversion.

[0231] The satellite-based synchronous inversion device for land surface temperature and land surface emissivity provided by the present invention realizes effective decoupling in the process of solving land-atmosphere parameters based on a multi-algorithm integrated inversion model based on non-coupled error optimization. At the same time, by weighting the land surface temperature and land surface emissivity inverted by multiple algorithms, the differences in the principle error levels of different algorithms are effectively balanced, thereby improving the accuracy and stability of the synchronous inversion of land surface temperature and emissivity.

[0232] See also Figure 3 , Figure 3 Schematic diagram of an electronic device provided by an embodiment of the present invention. Figure 3 As shown, an embodiment of the present invention provides an electronic device 300, including a memory 310, a processor 320, and a computer program 311 stored in the memory 310 and executable on the processor 320. When the processor 320 executes the computer program 311, the following steps are implemented:

[0233] Step 1: Collect nearby observation data of the target area based on thermal infrared remote sensing acquisition equipment;

[0234] Step 2: Calculate the initial values ​​of atmospheric parameters based on the RTTOV model and nearby observation data; the initial values ​​of atmospheric parameters include the atmospheric downwelling radiance Atmospheric upward radiance and atmospheric transmittance τ t ;

[0235] Step 3: Estimate based on the initial values ​​of atmospheric parameters and the I-DN algorithm to obtain the first surface temperature First surface emissivity ε 1 , Temperature profile correction value T a,t and humidity profile correction factor w t ; where t represents the observation time;

[0236] Step 4: Correct T based on temperature profile a,t and humidity profile correction factor w t For the initial atmospheric temperature profile T prof，t and the initial humidity profile Q prof，t Correction is performed and the corrected atmospheric parameters are calculated by the RTTOV model;

[0237] Step 5: Based on the first surface temperature First surface emissivity ε 1 , the corrected atmospheric parameters and the I-TES algorithm are used to estimate the second surface temperature and the second surface emissivity ε 2 ;

[0238] Step 6: Based on the first surface emissivity ε 1Second surface emissivity ε 2 , the weighted surface emissivity ε is calculated w ;

[0239] Step 7: Based on the weighted surface emissivity ε w The third surface temperature is estimated by using the I-SW algorithm

[0240] Step 8: Based on the first surface temperature Second surface temperature The third surface temperature Calculate weighted land surface temperature

[0241] Step 9: Based on weighted surface temperature and the surface emissivity ε w , calculate the cost function and repeat steps 3-8 until the function threshold is met or the maximum number of iterations is reached;

[0242] Step 10: Repeat steps 3-9 to complete the inversion of the surface temperature and surface emissivity of the target area;

[0243] Step 11: For the uninverted areas that do not meet the observation number greater than 2, first use the interpolation method to obtain the corresponding atmospheric correction parameters, and then use the I-TES algorithm to calculate the corresponding surface temperature and emissivity to complete the regional integrated inversion.

[0244] See also Figure 4 , Figure 4 Schematic diagram of an embodiment of a computer-readable storage medium provided in an embodiment of the present invention. Figure 4 As shown, this embodiment provides a computer-readable storage medium 400, on which a computer program 311 is stored. When the computer program 311 is executed by a processor, the following steps are implemented:

[0245] Step 1: Collect nearby observation data of the target area based on thermal infrared remote sensing acquisition equipment;

[0246] Step 2: Calculate the initial values ​​of atmospheric parameters based on the RTTOV model and nearby observation data; the initial values ​​of atmospheric parameters include the atmospheric downwelling radiance Atmospheric upward radiance and atmospheric transmittance τ t ;

[0247] Step 3: Estimate based on the initial values ​​of atmospheric parameters and the I-DN algorithm to obtain the first surface temperature First surface emissivity ε 1 , Temperature profile correction value T a,t and humidity profile correction factor w t; where t represents the observation time;

[0248] Step 4: Correct T based on temperature profile a,t and humidity profile correction factor w t For the initial atmospheric temperature profile T prof，t and the initial humidity profile Q prof，t Correction is performed and the corrected atmospheric parameters are calculated by the RTTOV model;

[0249] Step 5: Based on the first surface temperature First surface emissivity ε 1 , the corrected atmospheric parameters and the I-TES algorithm are used to estimate the second surface temperature and the second surface emissivity ε 2 ;

[0250] Step 6: Based on the first surface emissivity ε 1 Second surface emissivity ε 2 , the weighted surface emissivity ε is calculated w ;

[0251] Step 7: Based on the weighted surface emissivity ε w The third surface temperature is estimated by using the I-SW algorithm

[0252] Step 8: Based on the first surface temperature Second surface temperature The third surface temperature Calculate weighted land surface temperature

[0253] Step 9: Based on weighted surface temperature and the surface emissivity ε w , calculate the cost function and repeat steps 3-8 until the function threshold is met or the maximum number of iterations is reached;

[0254] Step 10: Repeat steps 3-9 to complete the inversion of the surface temperature and surface emissivity of the target area;

[0255] Step 11: For the uninverted areas that do not meet the observation number greater than 2, first use the interpolation method to obtain the corresponding atmospheric correction parameters, and then use the I-TES algorithm to calculate the corresponding surface temperature and emissivity to complete the regional integrated inversion.

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

[0257] It will be appreciated by those skilled in the art that embodiments of the present invention may be provided as methods, devices, or computer program products. Therefore, the present invention may take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Moreover, the present invention may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0258] The present invention is described with reference to flowcharts and / or block diagrams of methods, devices (apparatus), and computer program products according to embodiments of the present invention. It should be understood that each process and / or block in the flowchart and / or block diagram, as well as the combination of processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded computer, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 Process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0259] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 Process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0260] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process. Figure 1 Process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0261] Although the preferred embodiments of the present invention have been described, those skilled in the art may make other changes and modifications to these embodiments once they have learned the basic creative concept. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments and all changes and modifications that fall within the scope of the present invention.

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

Claims

1. A method for synchronous inversion of surface temperature and surface emissivity based on satellite technology, characterized in that: include: Step 1: Collect nearby observation data of the target area based on thermal infrared remote sensing acquisition equipment; Step 2: Calculate the initial values ​​of atmospheric parameters based on the RTTOV model and nearby observation data; The initial values ​​of atmospheric parameters include the atmospheric downwelling radiance Atmospheric upward radiance and atmospheric transmittance τ t ; Step 3: Estimate based on the initial values ​​of atmospheric parameters and the I-DN algorithm to obtain the first surface temperature First surface emissivity ε 1 , Temperature profile correction value T a，t and humidity profile correction factor w t ; where t represents the observation time; Step 4: Correct T based on temperature profile a，t and humidity profile correction factor w t For the initial atmospheric temperature profile T prof，t and the initial humidity profile Q prof，t Correction is performed and the corrected atmospheric parameters are calculated by the RTTOV model; Step 5: Based on the first surface temperature First surface emissivity ε 1 , the corrected atmospheric parameters and the I-TES algorithm are used to estimate the second surface temperature and the second surface emissivity ε 2 ; Step 6: Based on the first surface emissivity ε 1 Second surface emissivity ε 2 , the weighted surface emissivity ε is calculated w ; Step 7: Based on the weighted surface emissivity ε w The third surface temperature is estimated by using the I-SW algorithm Step 8: Based on the first surface temperature Second surface temperature The third surface temperature Calculate weighted land surface temperature Step 9: Based on weighted surface temperature and the surface emissivity ε w , calculate the cost function and repeat steps 3-8 until the function threshold is met or the maximum number of iterations is reached; Step 10: Repeat steps 3-9 to complete the inversion of the surface temperature and surface emissivity of the target area; Step 11: For the uninverted areas that do not meet the observation number greater than 2, first use the interpolation method to obtain the corresponding atmospheric correction parameters, and then use the I-TES algorithm to calculate the corresponding surface temperature and emissivity to complete the regional integrated inversion.

2. The method for synchronously inverting the surface temperature and the surface emissivity based on satellite technology according to claim 1, characterized in that: Weighted surface temperature The formula is: Among them, θ j Represents the surface emissivity weight coefficient, which can be determined by statistical regression or a small neural network. The default value is: θ j =[0.5, 0.3, 0.2].

3. The method for synchronously inverting the surface temperature and the surface emissivity based on satellite technology according to claim 1, characterized in that: The cost function is: Where C represents the cost function, N represents the number of thermal infrared bands, i represents the i-th thermal infrared band, and L i，t represents the on-board radiance observed in the i-th thermal infrared band, L(i, t) represents the on-board radiance estimated in the i-th thermal infrared band, B represents the Planck function, represents the weighted surface temperature, represents the weighted surface emissivity of the i-th thermal infrared band, represents the corrected atmospheric downward radiation in the ith thermal infrared band, represents the corrected atmospheric upward radiation of the ith thermal infrared band, T′ i，t represents the corrected atmospheric transmittance of the ith thermal infrared band.

4. The method for synchronously inverting land surface temperature and land surface emissivity based on satellite technology according to claim 1, characterized in that: The first surface temperature is estimated based on the initial values ​​of atmospheric parameters and the I-DN algorithm. First surface emissivity ε 1 , Temperature profile correction value T a，t and humidity profile correction factor w t , include: Step 3-1: Calculate the initialization vector x0. The expression of the initialization vector x0 is as follows: x0=[ΔT a，t ,w t ,T s，t ,ε i ],t∈[1,D],i∈[1,N]; t represents the tth observation; D represents the number of observations, and D≥2; i represents the i-th thermal infrared band; N represents the number of thermal infrared bands; ΔT a，t represents the corrected value of the atmospheric temperature profile at the tth observation; w t represents the atmospheric humidity profile correction coefficient of the tth observation; T s，t is the surface temperature observed at the tth time, and the near-surface air temperature is used as the initial value by default; ε i represents the surface emissivity of the i-th thermal infrared band, which can be obtained by the following formula: T b_max，t =max(B -1 (L i，t )); Where B represents the Planck function, Lit represents the directly observed radiance on the star, represents the atmospheric downwelling radiance, represents the atmospheric upward radiance, τ i，t represents the atmospheric transmittance, represents the radiance above the ground, B -1 represents the inverse function of Planck function; Step 3-2: Substitute the surface emissivity ε i Substituting into the following formula, we get the characteristic vector e of the surface emissivity; Among them, e i represents the i-th feature; f k represents the prediction function, which can be determined by statistical regression, machine learning, deep learning and other algorithms; k represents the kth prediction function, and K represents the number of prediction functions; Step 3-3: Replace the surface emissivity ε in the initialization vector x0 with the eigenvector e i , get the target vector x; x=[ΔT a，t ,w t ,T s，t，e ]; Step 3-4: By ΔT a，t and w t Modify the atmospheric profile and reconstruct the surface emissivity through the eigenvector e Step 3-5: Set the modified atmospheric parameters and T s，t and Input the RTTOV model to estimate the on-board radiance, denoted as L i，t (x); Step 3-6: Calculate the corrected value of x. The specific steps are as follows: Construct the correction matrix Δx: Where Δx represents the correction matrix; I n represents the unit matrix with dimension n; n is the length of the target vector; γ is the correction factor, whose value is set to 1e -6 ; Based on the correction matrix and RTTOV model, calculate the correction error ΔL of the on-board radiance i，t : ΔL i，t =L i，t (x+Δx)-L i，t ; Where i represents the i-th thermal infrared band, x represents the target vector, Δx represents the correction matrix, Lit(x+Δx) represents the estimated corrected on-board radiance, and L i，t It represents the on-board radiance of thermal infrared remote sensing observation; The x correction value δ is calculated by the following formula: Among them, pinv represents the pseudo-inverse function, γ represents the radiance scaling factor, when the wave number unit is used, the default value is 1000, when the wavelength unit is used, the default value is 1; Step 3-7: Correct the target vector x: x=x+α·δ·x Among them, α is the correction weight coefficient, which is determined by the observation conditions: a Ta +a w +a Ts +a ε =1; Among them, a Ta 、a w 、a Ts and a ε They represent the correction weights of temperature profile, humidity profile, surface temperature and surface emissivity, respectively; i represents the i-th thermal infrared band, c Ta，1 、c w，1 、c Ts，1 、c ε，1 Respectively represent a Ta 、a w 、a Ts 、a ε The weight coefficient can be calculated by numerical fitting or machine learning model algorithms. W is the water vapor content of the atmospheric profile, Ta represents the atmospheric temperature, and T b，i represents the equivalent brightness temperature on the satellite in the i-th band; N represents the number of thermal infrared bands; Step 3-8: Execute steps 3-1 to 3-7 repeatedly. When the norm of δ is less than the specified threshold or reaches the predetermined number of iterations, the estimation of the surface temperature and surface emissivity is completed and recorded as the first surface temperature. and the first surface emissivity ε 1 .

5. The method for synchronously inverting land surface temperature and land surface emissivity based on satellite technology according to claim 1, characterized in that: The specific steps for atmospheric profile correction in step 4 are as follows: According to the current atmospheric temperature profile, starting from the ground height, the temperature offset is added layer by layer in the atmospheric temperature profile. The specific formula is as follows: Among them, T′ a，l represents the corrected temperature of the first layer of the atmosphere; T a,0 represents the near-surface (0th layer) atmospheric temperature; T a，1 Indicates the initial temperature of the first layer of atmosphere; ΔT a represents the temperature offset; n1 represents the number of vertical layers of the atmospheric profile; According to the current atmospheric humidity profile, the atmospheric humidity profile is scaled as a whole. The calculation formula is as follows: Q’ l =Q l ·w,l∈[0,nl-1]; Among them, Q′ l represents the corrected humidity of the lth layer of atmosphere; Q represents the initial humidity of the lth layer of atmosphere; w represents the scaling factor of atmospheric water vapor; n1 represents the number of vertical layers of the atmospheric profile.

6. The method for synchronously inverting land surface temperature and land surface emissivity based on satellite technology according to claim 1, characterized in that: The I-TES algorithm in step 5 is an optimization algorithm based on the ASTER TES algorithm. The implementation steps are as follows: Step 5-1: From the surface emissivity ε 1 Extract the maximum value from Step 5-2: Calculate the initial value of the surface temperature using the following formula: Where i represents the i-th thermal infrared band, B represents the Planck function, T i,t represents the equivalent surface temperature, B -1 represents the inverse Planck function, It represents the band's ground radiance. represents the corrected atmospheric downward radiation; Step 5-3: Using the initial surface temperature Calculate the corresponding surface emissivity and its normalized spectrum The formula is as follows: Step 5-4: Calculate the minimum surface emissivity using the β-MMD hybrid model where f k represents the prediction function, which can be determined by algorithms such as statistical regression or machine learning, k represents the kth prediction function, and K represents the number of prediction functions. and represent the maximum and minimum values ​​of the normalized spectral emissivity, respectively; Step 5-5: Reconstruct the surface emissivity spectrum using the following formula: Step 5-6: Exploitation Re-estimate the surface temperature Denoted as the second surface emissivity and the second surface temperature respectively; Step 5-7: Loop through steps 5-1 to 5-6 until the cost function changes between the previous and next two times are less than the specified threshold. The cost function used is as follows: Where i represents the i-th thermal infrared band, t represents the t-th observation, and L i，t is the observed on-star radiance, represents the estimated on-star radiance, N represents the number of thermal infrared bands, and D represents the number of observations.

7. The method for synchronously inverting land surface temperature and land surface emissivity based on satellite technology according to claim 2, characterized in that: Based on the first surface emissivity ε 1 and the second surface emissivity ε 2 Calculate and get the weighted surface emissivity ε w ,include: The first surface emissivity ε 1 and the second surface emission ε 2 Substituting into the following formula, the weighted surface emissivity ε is calculated w ; Among them, φ j is the surface emissivity weight coefficient, which can be determined by statistical regression or a small neural network. The default value is: j =0.

5.

8. The method for synchronously inverting land surface temperature and land surface emissivity based on satellite technology according to claim 1, characterized in that: Based on the weighted surface emissivity ε w The third surface temperature is estimated by using the I-SW algorithm The steps include: Among them, f k represents the prediction function, which can be determined by statistical regression or machine learning algorithms, k represents the kth prediction function, K represents the number of prediction functions, i represents the i-th thermal infrared band, T b，i，t represents the equivalent stellar brightness temperature. Represents the surface emissivity.

9. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the program, the method for synchronous inversion of surface temperature and surface emissivity based on satellite technology as claimed in any one of claims 1 to 8 is implemented.

10. A non-transitory computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the method for synchronous inversion of surface temperature and surface emissivity based on satellite technology as claimed in any one of claims 1 to 8 is implemented.

Images