Thermal infrared hyperspectral surface temperature emissivity spectrum inversion method and device

Abstract

The present invention provides a method and device for inverting a thermal infrared hyperspectral surface temperature emissivity spectrum, and relates to the field of thermal infrared remote sensing technology. The method comprises: obtaining satellite remote sensing data of a thermal infrared hyperspectral spectrum, performing cloud detection on the satellite remote sensing data, screening out valid data without cloud coverage, and extracting the entrance pupil radiance and entrance pupil brightness temperature from the valid data; using a thermal infrared hyperspectral atmospheric correction model to perform atmospheric correction on the entrance pupil radiance and entrance pupil brightness temperature, and inverting to obtain the surface exit radiance and surface exit brightness temperature; based on the surface exit radiance and surface exit brightness temperature, using a pre-constructed emissivity sparse dictionary to reconstruct the emissivity spectrum; and using a pre-established minimum value empirical constraint to limit the value of the reconstructed emissivity spectrum, and simultaneously inverting to obtain the surface temperature and the surface emissivity spectrum. The present invention can realize the simultaneous inversion of the surface temperature and the surface emissivity spectrum in the 8-14μm thermal infrared spectrum.

Description

Technical Field

[0001] The present invention relates to the field of thermal infrared remote sensing technology, and in particular to a method and device for simultaneously inverting thermal infrared hyperspectral surface temperature and surface emissivity spectrum based on a sparse dictionary. Background Art

[0002] As an important branch of the earth observation field, thermal infrared remote sensing technology plays an indispensable role in the research of global earth-atmosphere energy, global climate change, urban heat island effect, environmental change, etc., because of its advantages of all-day and all-region observation. At the same time, with the rapid development of hyperspectral remote sensing technology, some substances that cannot be captured by wide-band remote sensors can be sensed. Compared with visible light and short-wave infrared, thermal infrared hyperspectral remote sensing has its own unique advantages. It is not only an important supporting technology for thermal infrared quantitative remote sensing, but also can be widely used in mineral exploration, atmospheric trace gas inversion and other fields.

[0003] Land Surface Temperature (LST) and Land Surface Emissivity (LSE), as inherent properties of the surface, are key parameters in the field of thermal infrared quantitative remote sensing. Thermal infrared (8-14μm) hyperspectral remote sensing surface temperature / emissivity inversion has always been one of the hot spots and difficulties in thermal infrared quantitative remote sensing surface parameter inversion. On the one hand, in the signals received by the satellite thermal infrared hyperspectral sensor, the outgoing radiation of the surface, the downward thermal radiation of the atmosphere reflected by the surface, and the upward thermal radiation of the atmosphere are strongly coupled and difficult to separate; on the other hand, in the radiation brightness / brightness temperature (hereinafter referred to as "brightness temperature") observed by the N channels of the thermal infrared hyperspectral sensor, assuming that the influence of atmospheric effects has been eliminated, there will always be N+1 unknowns to be solved (emissivity of N channels and 1 temperature).

[0004] At present, the thermal infrared remote sensing surface temperature / emissivity simultaneous inversion algorithms mainly include the dual-temperature method, MODIS day-night algorithm, and ASTER temperature / emissivity separation algorithm (Temperature-Emissivity Separation, TES). However, the dual-temperature method needs to assume that the surface emissivity does not change with time, and uses multiple day / night satellite observation data to increase the number of known conditions, thereby inverting the surface temperature and surface emissivity; the MODIS day-night algorithm uses the observation data of the seven MODIS medium / thermal infrared channels during the day and at night to form 14 nonlinear equations as known conditions, and inverts the surface temperature and emissivity through statistical regression and least squares fitting and other solving methods, but it is necessary to assume that some parameters are the same during the day and at night. The parameters mainly include the surface emissivity, the atmospheric water vapor content and temperature profile provided by the MODIS sounding channel and the corresponding algorithm; the TES algorithm uses the data of the five thermal infrared channels of ASTER, and uses the relationship between the maximum and minimum differences in relative emissivity and the minimum emissivity as a constraint. Under the condition of accurate atmospheric correction, the emissivity and temperature are inverted simultaneously with the help of the emissivity normalization method. However, this method not only depends on the accuracy of atmospheric parameters, but also has relatively poor inversion accuracy for land types with relatively flat emissivity spectra such as vegetation / water bodies.

[0005] Therefore, the existing simultaneous inversion method of surface temperature / emissivity still has the above-mentioned technical defects and needs further improvement. Summary of the invention

[0006] In view of the above problems, the present invention provides a thermal infrared hyperspectral surface temperature emissivity spectrum inversion method and device.

[0007] On the one hand, the present invention provides a method for inverting a thermal infrared hyperspectral surface temperature emissivity spectrum, comprising: acquiring satellite remote sensing data of a thermal infrared hyperspectral spectrum, performing cloud detection on the satellite remote sensing data, screening out valid data without cloud cover, and extracting the entrance pupil radiance and entrance pupil brightness temperature from the valid data; performing atmospheric correction on the entrance pupil radiance and entrance pupil brightness temperature using a thermal infrared hyperspectral atmospheric correction model, and inverting to obtain the surface exit radiance and surface exit brightness temperature; reconstructing the emissivity spectrum based on the surface exit radiance and surface exit brightness temperature using a pre-constructed emissivity sparse dictionary; and limiting the value of the reconstructed emissivity spectrum using a pre-established minimum value empirical constraint, and simultaneously inverting to obtain the surface temperature and the surface emissivity spectrum.

[0008] Furthermore, cloud products or threshold methods are used to perform cloud detection on satellite remote sensing data.

[0009] Furthermore, the valid data include pixel brightness values ​​without cloud coverage, and the entrance pupil radiation brightness and entrance pupil brightness temperature in the valid data are extracted, including: using the absolute radiation calibration coefficient of the thermal infrared hyperspectral load to perform absolute radiation correction on the pixel brightness value to obtain the entrance pupil radiation brightness; according to the entrance pupil radiation brightness, using the Planck function and the thermal infrared radiation transfer equation to calculate the entrance pupil brightness temperature.

[0010] Furthermore, the formula of the thermal infrared hyperspectral atmospheric correction model is as follows:

[0011]

[0012] Where, T gi is the surface brightness temperature of channel i; T i is the entrance pupil brightness temperature of channel i; T j is the pupil brightness temperature of other inversion channels j; N is the total number of channels involved in the inversion; θ v is the observation zenith angle; k0, k1, k 2j , k 3j , k4 and k5 are the empirical coefficients of fitting.

[0013] Furthermore, the emissivity sparse dictionary is constructed in the following way: the thermal infrared surface emissivity spectrum is selected from the public ground object emissivity spectrum library, and the channel response function of the thermal infrared hyperspectral sensor is used for convolution to obtain multiple channel equivalent emissivities; the orthogonal matching pursuit algorithm and K singular value decomposition algorithm are used to perform sparse coding and dictionary training on multiple channel equivalent emissivities to establish an emissivity sparse dictionary.

[0014] Furthermore, the emissivity sparse dictionary is constructed according to the following formula:

[0015] y h =Φ OMP ε n

[0016] y h =D ε x t

[0017] D ε =K-SVD update (D ε , x t )

[0018] In the formula, ε n is the thermal infrared emissivity spectrum in the ground object emissivity spectrum library, and the vector length is N; Φ OMP is the compressed observation matrix of the OMP algorithm in sparse coding; h is n The compressed observation value of length H; D εis the emission rate sparse dictionary; vector x t for y h In D ε sparse representation in the transform domain, and x t is T-sparse, that is, x t The number of non-zero atoms in is T, T<H<N; K-SVD update Indicates update D ε The process of atoms.

[0019] Furthermore, a pre-built emissivity sparse dictionary is used to reconstruct the emissivity spectrum, including: sparse and empirical reconstruction of emissivity unknowns using an orthogonal matching pursuit algorithm and an emissivity sparse dictionary based on the surface emergent brightness temperature.

[0020] Furthermore, the emissivity spectrum is reconstructed according to the following formula:

[0021]

[0022] In the formula, ε i is the i-th unknown number, i=1, 2, ..., N; is Φ OMP The inverse function of B(T s ) is the surface temperature T s Radiant brightness under gi ) is the surface outgoing radiation brightness of channel i after reconstruction.

[0023] Furthermore, the minimum value of the reconstructed emissivity spectrum is calculated according to the following formula:

[0024]

[0025]

[0026]

[0027]

[0028] In the formula, is an unbiased estimate of the emissivity spectrum; ε i is the emissivity spectrum value; is the average value of the emissivity spectra of all channels; MMD is the difference between the maximum and minimum values ​​of the unbiased estimate of the emissivity; STD is the standard deviation of the unbiased estimate of the emissivity; ε min is the minimum value of the emissivity spectrum; u, v1, v2, w1 and w2 are all empirical coefficients of fitting.

[0029] Furthermore, the minimum empirical constraint is established according to the following formula:

[0030]

[0031] In the formula, is the minimum value of the unbiased estimate spectrum of the emissivity; ε′ i is the reconstructed emissivity spectrum.

[0032] Furthermore, before simultaneously inverting the surface temperature and the surface emissivity spectrum, the method further includes: constructing a cost function and iterating based on the emissivity normalization method to simultaneously invert the surface temperature and the surface emissivity spectrum, wherein the cost function is constructed according to the following formula:

[0033]

[0034]

[0035] Where t is the number of iterations; is the reconstructed emissivity spectrum at the tth iteration, T0 is the initial temperature, which is also the average of the surface brightness temperatures of all channels; B(T0) is the radiation brightness at the initial temperature; is the surface outgoing radiation brightness after the t-th iteration reconstruction; F cost To construct a good cost function;

[0036] When the value of the cost function obtained by each iterative reconstruction is less than the preset cost threshold, or the number of iterations is greater than the preset threshold, the iteration is stopped, and the surface temperature and surface emissivity spectrum are inverted at the same time.

[0037] On the other hand, the present invention provides a thermal infrared hyperspectral surface temperature emissivity spectrum inversion device, comprising: a satellite entrance pupil data extraction module, used to obtain satellite remote sensing data of the thermal infrared hyperspectral, perform cloud detection on the satellite remote sensing data, screen out valid data without cloud cover, and extract the entrance pupil radiance and entrance pupil brightness temperature in the valid data; a surface exit data inversion module, used to use a thermal infrared hyperspectral atmospheric correction model to perform atmospheric correction on the entrance pupil radiance and entrance pupil brightness temperature, and invert the surface exit radiance and surface exit brightness temperature; a sparse dictionary reconstruction module, used to reconstruct the emissivity spectrum according to the surface exit radiance and surface exit brightness temperature using a pre-constructed emissivity sparse dictionary; and a temperature emissivity simultaneous inversion module, used to use a pre-established minimum value empirical constraint to limit the value of the reconstructed emissivity spectrum, and simultaneously invert the surface temperature and surface emissivity spectrum.

[0038] Compared with the prior art, the thermal infrared hyperspectral surface temperature emissivity spectrum inversion method and device provided by the present invention has at least the following beneficial effects:

[0039] (1) This method uses thermal infrared hyperspectral remote sensing data as input and realizes atmospheric correction of sensor entrance pupil radiance / brightness temperature based only on the characteristic information of hyperspectral thermal infrared data. This atmospheric correction model is free from dependence on external atmospheric parameters.

[0040] (2) This method uses the emissivity sparse dictionary and the minimum emissivity constraint as prior knowledge, and uses the sparsity of the sparse dictionary to sparse the unknowns of the surface emissivity, greatly reducing the number of unknowns in the inversion process and eliminating the pathological nature of the surface parameter inversion of the traditional method. The sparse dictionary is used to empirically reconstruct the emissivity spectrum characteristics, eliminating the influence of the atmospheric correction residual on the emissivity spectrum, ensuring the inversion accuracy of the surface temperature / emissivity in the atmospheric window and fully restoring the emissivity information in the strong absorption band of the atmosphere. The minimum emissivity constraint is used to further limit the emissivity value to a reasonable range, so that the iterative solution process can converge quickly.

[0041] (3) This method can achieve the simultaneous inversion of surface temperature and surface emissivity spectrum in the 8-14μm thermal infrared spectrum. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] The above and other objects, features and advantages of the present invention will become more apparent through the following description of the embodiments of the present invention with reference to the accompanying drawings, in which:

[0043] Figure 1 A flow chart of a thermal infrared hyperspectral surface temperature emissivity spectrum inversion method according to an embodiment of the present invention is schematically shown;

[0044] Figure 2 A flowchart of a process for extracting entrance pupil radiance and entrance pupil brightness temperature according to an embodiment of the present invention is schematically shown;

[0045] Figure 3 Schematically shows a flow chart of a process for constructing an emissivity sparse dictionary according to an embodiment of the present invention;

[0046] Figure 4 Schematically shows a flow chart of a process of reconstructing an emissivity spectrum using an emissivity sparse dictionary according to an embodiment of the present invention;

[0047] Figure 5 Schematically shows a flowchart of a cost function construction and iterative reconstruction process according to an embodiment of the present invention;

[0048] Figure 6 The block diagram of the thermal infrared hyperspectral ground surface temperature emissivity spectrum inversion device according to an embodiment of the present invention is schematically shown. DETAILED DESCRIPTION

[0049] In order to make the purpose, technical scheme and advantages of the present invention more clearly understood, the present invention is further described in detail below in combination with specific embodiments and with reference to the accompanying drawings. Obviously, the described embodiments are 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 ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0050] The terms used herein are only for describing specific embodiments and are not intended to limit the present invention. The terms "comprise", "include", etc. used herein indicate the existence of the features, steps, operations and / or components, but do not exclude the existence or addition of one or more other features, steps, operations or components.

[0051] All terms (including technical and scientific terms) used herein have the meanings commonly understood by those skilled in the art, unless otherwise defined. It should be noted that the terms used herein should be interpreted as having a meaning consistent with the context of this specification, and should not be interpreted in an idealized or overly rigid manner.

[0052] Figure 1 The flowchart of the thermal infrared hyperspectral surface temperature emissivity spectrum inversion method according to an embodiment of the present invention is schematically shown.

[0053] like Figure 1 As shown, the thermal infrared hyperspectral land surface temperature emissivity spectrum inversion method according to this embodiment may include operations S110 to S160.

[0054] In operation S110, satellite remote sensing data of thermal infrared hyperspectral is acquired, cloud detection is performed on the satellite remote sensing data, valid data without cloud cover is screened out, and entrance pupil radiance and entrance pupil brightness temperature are extracted from the valid data.

[0055] In operation S120, the entrance pupil radiance and the entrance pupil brightness temperature are atmospherically corrected using a thermal infrared hyperspectral atmospheric correction model, and the surface exit radiance and the surface exit brightness temperature are inverted to obtain.

[0056] In operation S130 , an emissivity spectrum is reconstructed using a pre-constructed emissivity sparse dictionary according to the surface emergent radiance and the surface emergent brightness temperature.

[0057] In operation S140, a pre-established minimum value empirical constraint is used to restrict the value of the reconstructed emissivity spectrum, and the surface temperature and the surface emissivity spectrum are simultaneously inverted.

[0058] Through the embodiments of the present invention, the effective data eliminates the influence of cloud cover, iteratively reconstructs the emissivity spectrum based on the sparseness of unknowns in the emissivity sparse dictionary and the recovery of the emissivity spectrum information, eliminates the atmospheric correction residual error in the emissivity spectrum, and limits the emissivity value with the help of minimum constraint, so as to realize the simultaneous inversion of surface temperature and surface emissivity spectrum in the 8-14μm thermal infrared spectrum.

[0059] The following will be combined Figure 2 to Figure 5 ,right Figure 1 The disclosed thermal infrared hyperspectral surface temperature emissivity spectrum inversion method is described in detail.

[0060] In this embodiment, cloud products or threshold methods are used to detect clouds in satellite remote sensing data. Specifically, cloud products produced on the same satellite platform as the thermal infrared hyperspectral payload are preferred. If there is no cloud product assistance, a cloud discrimination threshold can be set based on the representation of clouds in thermal infrared hyperspectral data to perform cloud detection and shield or remove cloud-containing pixels.

[0061] Figure 2 The flowchart of the process of extracting entrance pupil radiance and entrance pupil brightness temperature according to an embodiment of the present invention is schematically shown.

[0062] like Figure 2 As shown, in this embodiment, the valid data in the above operation S110 includes the brightness value of the pixel without cloud coverage, and extracting the entrance pupil radiation brightness and entrance pupil brightness temperature in the valid data may include operations S1101 to S1102.

[0063] In operation S1101, the absolute radiation calibration coefficient of the thermal infrared hyperspectral load is used to perform absolute radiation correction on the pixel brightness value to obtain the entrance pupil radiation brightness.

[0064] The pixel brightness value (Digital Number, DN), that is, the digital measurement value obtained by the sensor, is related to the sensor's radiation resolution, ground object emission frequency, atmospheric transmittance, and scattering rate, and can reflect the ground object's emissivity.

[0065] In this embodiment, the absolute radiation correction formula is:

[0066] L i =B(T i )=gain i ·DN i +bias i

[0067] Where, DN i Indicates the pixel brightness value of channel i; gain i Indicates the gain value of the sensor; bias i Indicates the deviation value of the sensor; Li Represents the entrance pupil radiance.

[0068] In operation S1102, the entrance pupil brightness temperature is calculated based on the entrance pupil radiance using the Planck function and the thermal infrared radiation transfer equation.

[0069] Usually, the Planck function is used to realize the conversion between radiation brightness and brightness temperature. The Planck function can be expressed as:

[0070]

[0071] In the formula, M(λ i , T i ) indicates that the wavelength is λ i And the temperature is T i The radiant brightness at that time is usually expressed in W / (m 2 ·sr·μm); constant c1=1.191×10 8 W / (m 2 ·sr·μm -4 ); constant c2=1.439×10 4 K·μm.

[0072] Then, based on the thermal infrared radiation transfer equation, for channel i, the entrance pupil brightness temperature of the thermal infrared hyperspectral remote sensing satellite can be expressed as:

[0073]

[0074] Where B represents the Planck function; B -1 represents the inverse Planck function; τ i represents the atmospheric transmittance; ε i represents the surface emissivity; T s represents the surface temperature; It represents the downward thermal radiation of the atmosphere; Represents the upward thermal radiation of the atmosphere; T i is the satellite entrance pupil brightness temperature of channel i.

[0075] The surface brightness temperature T of channel i gi Defined as:

[0076] T gi =B -1 [ε i B(T s )]

[0077] In the above operation S120, the obtained entrance pupil brightness temperature is atmospherically corrected using the thermal infrared hyperspectral atmospheric correction model and the fitted empirical coefficients, and the surface exit brightness temperature is inverted. The formula of the thermal infrared hyperspectral atmospheric correction model is as follows:

[0078]

[0079] Where, T i is the entrance pupil brightness temperature of channel i; T j is the pupil brightness temperature of other inversion channels j; N is the total number of channels involved in the inversion; θ v is the observation zenith angle; k0, k1, k 2j , k 3j , k4 and k5 are the empirical coefficients of fitting.

[0080] For the formula of the above thermal infrared hyperspectral atmospheric correction model, the method for obtaining the fitting empirical coefficients is as follows: based on the atmospheric radiation transfer model MODTRAN (MODerate Resolution Atmospheric TRANsmission), with the support of the TIGR atmospheric profile database, by reasonably setting the surface temperature, surface emissivity and observation angle, combined with the satellite sensor channel response function, the satellite entrance pupil radiation brightness / brightness temperature and the surface exit radiation brightness / brightness temperature are simulated. The simulated data covers the thermal infrared radiation transmission under different surface, atmospheric and observation conditions to ensure the universality of the atmospheric correction algorithm; then, a large amount of entrance pupil brightness temperature simulation data, surface exit brightness temperature simulation data and the corresponding observation angles are used to determine the fitting empirical coefficients through multivariate regression fitting.

[0081] Among them, the TIGR atmospheric profile database contains data on atmospheric temperature profiles, water vapor profiles, ozone profiles, etc. in polar, tropical and mid-latitude regions. The surface temperature is set with reference to the bottom temperature of the atmospheric profile; the surface emissivity is selected from the ASTRE2.0 land object emissivity spectral library, which includes emissivity spectra of various typical land objects such as water bodies, vegetation, soil, rocks and man-made objects.

[0082] Through the embodiments of the present invention, thermal infrared hyperspectral remote sensing data is used as input, and atmospheric correction of sensor entrance pupil radiance / brightness temperature is achieved only based on the characteristic information of the hyperspectral thermal infrared data itself. The atmospheric correction model gets rid of the dependence on external atmospheric parameters.

[0083] Due to the strong absorption of certain internal particles in the atmosphere (such as water vapor, carbon dioxide, etc.) in the thermal infrared spectrum, there are residual errors in the atmospheric correction process, which makes the surface emission brightness temperature obtained by atmospheric correction have certain noise and information loss, and further destroys the emissivity spectrum information in the surface emission brightness temperature. In order to suppress the influence of noise on the inversion accuracy of surface parameters as much as possible, and at the same time restore the emissivity spectrum characteristic information destroyed by atmospheric effects, the embodiment of the present invention also needs to use sparse dictionary learning to perform high-precision empirical reconstruction of the emissivity spectrum to achieve simultaneous inversion of surface temperature / emissivity spectrum.

[0084] Figure 3 The flowchart of the emissivity sparse dictionary construction process according to an embodiment of the present invention is schematically shown.

[0085] like Figure 3 As shown, in this embodiment, the emissivity sparse dictionary in the above operation S130 can be constructed through the following operations S1301 to S1302.

[0086] In operation S1301, a thermal infrared surface emissivity spectrum is selected from a public surface object emissivity spectrum library, and convolution is performed using a channel response function of a thermal infrared hyperspectral sensor to obtain multiple channel equivalent emissivities.

[0087] Preferably, the disclosed ground object emissivity spectrum library may be the ASTRE 2.0 ground object emissivity spectrum library, which includes a variety of thermal infrared ground surface emissivity spectra.

[0088] In operation S1302, an orthogonal matching pursuit algorithm (Orthogonal Matching Pursuit, OMP) and a K singular value decomposition algorithm (K-SVD) are used to perform sparse coding and dictionary training on equivalent emissivity of multiple channels to establish an emissivity sparse dictionary.

[0089] Specifically, the above operations S1301 to S1302 may refer to the following formula:

[0090] y h =Φ OMP ε n

[0091] y h =D ε x t

[0092] D ε =K-SVD update (D ε , x t )

[0093] In the formula, ε n is the thermal infrared emissivity spectrum in the ground object emissivity spectrum library, and the vector length is N; Φ OMP is the compressed observation matrix of the OMP algorithm in sparse coding; h is n The compressed observation value of length H; D ε is the emission rate sparse dictionary; vector x t for y h In D εsparse representation in the transform domain, and x t is T-sparse, that is, x t The number of non-zero atoms in is T. In general, T<H <N; K-SVD update Indicates update D ε The process of atoms.

[0094] Figure 4 The flowchart of the process of reconstructing the emissivity spectrum using the emissivity sparse dictionary according to an embodiment of the present invention is schematically shown.

[0095] like Figure 4 As shown, next, the above operation S130 of using the pre-constructed emissivity sparse dictionary to reconstruct the emissivity spectrum may further include operation S1303.

[0096] In operation S1303, the emissivity unknowns are sparsely and empirically reconstructed using an orthogonal matching pursuit algorithm and an emissivity sparse dictionary according to the surface emission brightness temperature.

[0097] Specifically, the surface brightness temperature T of channel i is known after atmospheric correction. gi When using the OMP algorithm and the emissivity sparse dictionary D ε , the unknown emission rate is sparsely and empirically reconstructed, that is:

[0098]

[0099] In the formula, ε i is the i-th unknown number, i=1, 2, ..., N; is Φ OMP The inverse function of B(T s ) is the surface temperature T s Radiant brightness under gi ) is the surface outgoing radiation brightness of channel i after reconstruction.

[0100] At this time, N conditional equations can be established based on the surface brightness temperatures of N channels, and the number of unknowns (that is, T x t The number of non-zero values ​​in + 1 surface temperature T s ) is less than the number of equations, the formula is solved. Since it is impossible to directly obtain x t T non-zero values ​​and the surface temperature T s Therefore, the iterative calculation method is used to solve the unknown variables.

[0101] In order to ensure the uniqueness of the solution during iterative solving and improve the speed of iterative convergence, the embodiment of the present invention establishes an emissivity minimum constraint, which mainly uses the statistical characteristics of the emissivity spectrum itself.

[0102] Specifically, the minimum value of the reconstructed emissivity spectrum is calculated according to the following formula:

[0103]

[0104]

[0105]

[0106]

[0107] In the formula, is an unbiased estimate of the emissivity spectrum; ε i is the emissivity spectrum value; is the average value of the emissivity spectra of all channels; MMD is the difference between the maximum and minimum values ​​of the unbiased estimate of the emissivity; STD is the standard deviation of the unbiased estimate of the emissivity; ε min is the minimum value of the emissivity spectrum; u, v1, v2, w1 and w2 are all empirical coefficients of fitting.

[0108] Among them, the emissivity in the spectral library is convolved by the response function of the thermal infrared hyperspectral sensor to obtain the channel equivalent emissivity, and the empirical coefficient can be fitted through the above statistical relationship.

[0109] Furthermore, the minimum empirical constraint is established according to the following formula:

[0110]

[0111] In the formula, is the minimum value of the unbiased estimate spectrum of the emissivity; ε′ i is the reconstructed emissivity spectrum.

[0112] Therefore, the minimum empirical constraint is derived based on the statistical relationship of the emissivity spectrum. Based on the minimum empirical constraint, the value of the emissivity spectrum after sparse dictionary reconstruction is restricted to ensure the uniqueness of the solution, so that the iterative solution process can converge quickly.

[0113] Figure 5 The flowchart of the cost function construction and iterative reconstruction process according to an embodiment of the present invention is schematically shown.

[0114] like Figure 5 As shown, in this embodiment, before the simultaneous inversion to obtain the surface temperature and the surface emissivity spectrum in the above operation S140, operations S1401 to S1402 may also be included.

[0115] In operation S1401 , a cost function is constructed and iterated based on an emissivity normalization method to simultaneously invert the surface temperature and the surface emissivity spectrum.

[0116] Among them, the cost function is constructed according to the following formula:

[0117]

[0118]

[0119] Where t is the number of iterations; is the reconstructed emissivity spectrum at the tth iteration, T0 is the initial temperature, which is also the average of the surface brightness temperatures of all channels; B(T0) is the radiation brightness at the initial temperature; is the surface outgoing radiation brightness after the t-th iteration reconstruction; F cost is the constructed cost function.

[0120] In the iterative solution process, the average value of the surface brightness temperature needs to be taken as the current initial temperature T t (t represents the number of iterations), the current initial temperature T t Substitute it into the emissivity sparse dictionary to reconstruct the emissivity, remove the atmospheric correction residual error contained in the emissivity spectrum, and use the minimum emissivity constraint to limit the emissivity value to obtain the reconstructed emissivity The surface brightness temperature is recalculated by each iterative reconstruction Construct the cost function F cost .

[0121] In operation S1402, when the value of the cost function obtained by each iterative reconstruction is less than a preset cost threshold, or the number of iterations is greater than a preset number threshold, the iteration is stopped, and the surface temperature and surface emissivity spectrum are inverted at the same time.

[0122] For example, the preset cost threshold may be 0.01, and the preset number threshold may be 12. Thus, when the cost function F is reconstructed each time, cost <0.01, or the current iteration number t>12, the iteration stops and the T t and They are respectively used as the final simultaneously inverted surface temperature and surface emissivity spectra.

[0123] Through the embodiments of the present invention, the emissivity can be sparsely realized through the above-mentioned sparse dictionary prior knowledge, which greatly reduces the unknowns in the inversion. In addition, the empirical reconstruction of the emissivity spectrum information through the sparse dictionary eliminates the noise introduced in the emissivity spectrum due to the residual error of atmospheric correction. With the support of the emissivity minimum constraint, the simultaneous inversion of thermal infrared surface temperature and emissivity spectrum is achieved.

[0124] Based on the above disclosed method, the present invention also provides a thermal infrared hyperspectral surface temperature emissivity spectrum inversion device. Figure 6 The device is described in detail.

[0125] Figure 6 The block diagram of the thermal infrared hyperspectral ground surface temperature emissivity spectrum inversion device according to an embodiment of the present invention is schematically shown.

[0126] like Figure 6 As shown, according to this embodiment, the thermal infrared hyperspectral surface temperature emissivity spectrum inversion device 600 includes a satellite entrance pupil data extraction module 610, a surface exit data inversion module 620, a sparse dictionary reconstruction module 630 and a temperature emissivity simultaneous inversion module 640.

[0127] The satellite entrance pupil data extraction module 610 is used to obtain satellite remote sensing data of thermal infrared hyperspectral, perform cloud detection on the satellite remote sensing data, filter out valid data without cloud cover, and extract entrance pupil radiation brightness and entrance pupil brightness temperature from the valid data.

[0128] The surface emission data inversion module 620 is used to perform atmospheric correction on the entrance pupil radiance and entrance pupil brightness temperature using a thermal infrared hyperspectral atmospheric correction model, and invert the surface emission radiance and surface emission brightness temperature.

[0129] The sparse dictionary reconstruction module 630 is used to reconstruct the emissivity spectrum using a pre-constructed emissivity sparse dictionary according to the surface emergent radiance and the surface emergent brightness temperature.

[0130] The temperature-emissivity simultaneous inversion module 640 is used to use a pre-established minimum value empirical constraint to restrict the value of the reconstructed emissivity spectrum, and simultaneously invert the surface temperature and the surface emissivity spectrum.

[0131] It should be noted that the implementation method of the device part is similar to the implementation method of the aforementioned method part, and the technical effects achieved are also similar. Please refer to the above-mentioned method implementation method part for specific details, which will not be repeated here.

[0132] According to an embodiment of the present invention, any multiple of the satellite entrance pupil data extraction module 610, the surface exit data inversion module 620, the sparse dictionary reconstruction module 630 and the temperature emissivity simultaneous inversion module 640 can be combined in one module for implementation, or any one of the modules can be split into multiple modules. Alternatively, at least part of the functions of one or more of these modules can be combined with at least part of the functions of other modules and implemented in one module. According to an embodiment of the present invention, at least one of the satellite entrance pupil data extraction module 610, the surface exit data inversion module 620, the sparse dictionary reconstruction module 630 and the temperature emissivity simultaneous inversion module 640 can be at least partially implemented as a hardware circuit, such as a field programmable gate array (FPGA), a programmable logic array (PLA), a system on a chip, a system on a substrate, a system on a package, an application specific integrated circuit (ASIC), or can be implemented by hardware or firmware such as any other reasonable way of integrating or packaging the circuit, or implemented in any one of the three implementation methods of software, hardware and firmware or in a suitable combination of any of them. Alternatively, at least one of the satellite entrance pupil data extraction module 610, the surface exit data inversion module 620, the sparse dictionary reconstruction module 630 and the temperature emissivity simultaneous inversion module 640 can be at least partially implemented as a computer program module, which can perform the corresponding function when the computer program module is executed.

[0133] Some block diagrams and / or flow charts are shown in the accompanying drawings. It should be understood that some blocks or combinations thereof in the block diagrams and / or flow charts may be implemented by computer program instructions. These computer program instructions may be provided to a processor of a general-purpose computer, a special-purpose computer, or other programmable data processing device, so that these instructions, when executed by the processor, may create a device for implementing the functions / operations described in these block diagrams and / or flow charts.

[0134] In addition, the terms "first" and "second" are used only for descriptive purposes and should not be understood as indicating or implying relative importance or implicitly indicating the number of the indicated technical features. Therefore, a feature defined as "first" or "second" may explicitly or implicitly include one or more of the feature. In the description of the present invention, the meaning of "plurality" is at least two, such as two, three, etc., unless otherwise clearly and specifically defined. In addition, the word "one" or "an" preceding an element does not exclude the presence of multiple such elements.

[0135] The specific embodiments described above further illustrate the objectives, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A thermal infrared hyperspectral surface temperature emissivity spectrum inversion method, characterized in that: include: Acquire satellite remote sensing data of thermal infrared hyperspectral, perform cloud detection on the satellite remote sensing data, screen out valid data without cloud cover, and extract entrance pupil radiation brightness and entrance pupil brightness temperature from the valid data; Using a thermal infrared hyperspectral atmospheric correction model, the entrance pupil radiance and entrance pupil brightness temperature are atmospherically corrected to obtain the surface exit radiance and surface exit brightness temperature; According to the surface emergent radiance and the surface emergent brightness temperature, the emissivity spectrum is reconstructed using a pre-constructed emissivity sparse dictionary; wherein the emissivity sparse dictionary is constructed in the following manner: a thermal infrared surface emissivity spectrum in a public surface object emissivity spectrum library is selected, and a channel response function of a thermal infrared hyperspectral sensor is used for convolution to obtain a plurality of channel equivalent emissivities; an orthogonal matching pursuit algorithm and a K singular value decomposition algorithm are used to perform sparse coding and dictionary training on the plurality of channel equivalent emissivities to establish the emissivity sparse dictionary; the method of reconstructing the emissivity spectrum using the pre-constructed emissivity sparse dictionary includes: according to the surface emergent brightness temperature, an orthogonal matching pursuit algorithm and the emissivity sparse dictionary are used to perform sparse and empirical reconstruction on emissivity unknowns; and The pre-established minimum empirical constraint is used to limit the value of the reconstructed emissivity spectrum, and the surface temperature and the surface emissivity spectrum are inverted at the same time, wherein the minimum empirical constraint is established according to the following formula: In the formula, is an unbiased estimate of the emissivity spectrum; is the emissivity spectrum value; is the average of the emissivity spectra of all channels; is the minimum of the spectrum for the unbiased estimate of the emissivity; is the minimum value of the emissivity spectrum; is the reconstructed emissivity spectrum.

2. The method according to claim 1, characterized in that Cloud detection is performed on the satellite remote sensing data using cloud products or a threshold method.

3. The method according to claim 1, characterized in that The valid data includes the brightness value of the pixel without cloud coverage, and the entrance pupil radiation brightness and the entrance pupil brightness temperature in the valid data are extracted, including: Using the absolute radiation calibration coefficient of the thermal infrared hyperspectral load, the pixel brightness value is subjected to absolute radiation correction to obtain the entrance pupil radiation brightness; The entrance pupil brightness temperature is calculated based on the entrance pupil radiation brightness using the Planck function and the thermal infrared radiation transfer equation.

4. The method according to claim 1, characterized in that: The formula of the thermal infrared hyperspectral atmospheric correction model is as follows: Where, T gi is the surface brightness temperature of channel i; T i is the entrance pupil brightness temperature of channel i; T j is the pupil brightness temperature of other inversion channels j; N is the total number of channels involved in the inversion; θ v is the observation zenith angle; k0, k1, k 2j , k 3j , k4 and k5 are the empirical coefficients of fitting.

5. The method according to claim 1, characterized in that The emissivity sparse dictionary is constructed according to the following formula: In the formula, is the thermal infrared emissivity spectrum in the ground object emissivity spectrum library, and the vector length is N; It is the compressed observation matrix of OMP algorithm in sparse coding; for The compressed observation value of is of length H; is the emissivity sparse dictionary; vector for exist sparse representation in the transform domain, and is T-sparse, that is The number of non-zero atoms in is T, T<H<N; Indicates update The process of atoms.

6. The method according to claim 5, characterized in that The emissivity spectrum is reconstructed according to the following formula: In the formula, is the i-th unknown number, i=1,2,…,N; for The inverse function of is the surface temperature Radiance under is the surface outgoing radiation brightness of channel i after reconstruction.

7. The method according to claim 1, characterized in that The minimum value of the reconstructed emissivity spectrum is calculated according to the following formula: Where MMD is the difference between the maximum and minimum values ​​of the unbiased estimate of the emissivity; STD is the standard deviation of the unbiased estimate of the emissivity; u, v1, v2, w1 and w2 are all empirical coefficients of the fitting.

8. The method according to claim 1, characterized in that Before the simultaneous inversion to obtain the surface temperature and surface emissivity spectrum, the method further includes: A cost function is constructed and iterated based on the emissivity normalization method to simultaneously invert the surface temperature and the surface emissivity spectrum, wherein the cost function is constructed according to the following formula: Where t is the number of iterations; is the reconstructed emissivity spectrum at the tth iteration, ; T0 is the initial temperature, which is also the average of the surface brightness temperatures of all channels; is the radiance at the initial temperature; is the surface outgoing radiation brightness after the t-th iteration reconstruction; is the constructed cost function; When the iterative cost function is less than a preset cost threshold, or the number of iterations is greater than a preset number threshold, the iteration is stopped to obtain the simultaneous inversion results of the surface temperature and the surface emissivity spectrum.

9. A thermal infrared hyperspectral surface temperature emissivity spectrum inversion device, characterized in that: include: A satellite entrance pupil data extraction module is used to obtain satellite remote sensing data of thermal infrared hyperspectral, perform cloud detection on the satellite remote sensing data, filter out valid data without cloud cover, and extract entrance pupil radiation brightness and entrance pupil brightness temperature from the valid data; A surface emission data inversion module is used to perform atmospheric correction on the entrance pupil radiance and entrance pupil brightness temperature using a thermal infrared hyperspectral atmospheric correction model, and invert the surface emission radiance and surface emission brightness temperature; A sparse dictionary reconstruction module is used to reconstruct the emissivity spectrum according to the surface emergent radiance and the surface emergent brightness temperature using a pre-constructed emissivity sparse dictionary; wherein the emissivity sparse dictionary is constructed in the following manner: selecting a thermal infrared surface emissivity spectrum from a public surface object emissivity spectrum library, and convolving it using a channel response function of a thermal infrared hyperspectral sensor to obtain a plurality of channel equivalent emissivities; using an orthogonal matching pursuit algorithm and a K singular value decomposition algorithm, sparsely encoding and dictionary training the plurality of channel equivalent emissivities to establish the emissivity sparse dictionary; the method of reconstructing the emissivity spectrum using a pre-constructed emissivity sparse dictionary includes: sparsely and empirically reconstructing emissivity unknowns using an orthogonal matching pursuit algorithm and the emissivity sparse dictionary according to the surface emergent brightness temperature; and The temperature emissivity simultaneous inversion module is used to use a pre-established minimum value empirical constraint to limit the value of the reconstructed emissivity spectrum, and simultaneously invert the surface temperature and the surface emissivity spectrum, wherein the minimum value empirical constraint is established according to the following formula: In the formula, is an unbiased estimate of the emissivity spectrum; is the emissivity spectrum value; is the average of the emissivity spectra of all channels; is the minimum of the spectrum for the unbiased estimate of the emissivity; is the minimum value of the emissivity spectrum; is the reconstructed emissivity spectrum.

Images