A thermal infrared hyperspectral temperature and emissivity inversion method and system based on a Gaussian prior model

Abstract

A method and system for inverting temperature and emissivity in thermal infrared hyperspectral data based on a Gaussian prior model, relating to the field of electronic digital data processing systems, is disclosed. The method includes: acquiring hyperspectral observation data and a pre-constructed emissivity sample library; calculating the mean vector and covariance matrix of the emissivity samples in the emissivity sample library to construct a Gaussian prior model of emissivity; constructing a Bayesian inversion model of the hyperspectral observation data based on the thermal infrared radiative transfer equation, observation noise distribution, and the Gaussian prior model; iteratively optimizing the temperature and emissivity in the Bayesian inversion model until convergence to obtain the optimal temperature and optimal emissivity; calculating the approximate covariance of the Bayesian inversion model under the optimal temperature and optimal emissivity conditions to obtain an uncertainty measure; and filtering the optimal temperature and optimal emissivity based on the uncertainty measure to obtain the target temperature and target emissivity. Implementing this application can reduce numerical calculation bias.

Description

Technical Field

[0001] This application relates to the field of electronic digital data processing systems, and in particular to a method and system for inverting thermal infrared hyperspectral temperature and emissivity based on a Gaussian prior model. Background Technology

[0002] Thermal infrared hyperspectral remote sensing technology refers to the technique of capturing the Earth's surface thermal radiation energy in a continuous and narrow range of bands within the thermal infrared spectral region using sensors mounted on airborne or spacecraft platforms. Within this technology, the mathematical calculation process of deriving the true physical parameters of the Earth's surface (such as surface temperature and emissivity) from the radiance observation data acquired by the sensors is called inversion. Since a sensor with N bands can only provide N observation equations, while the parameters to be determined include emissivity in N bands and one surface temperature, the total number of unknowns is always greater than the number of observation equations. This separation of temperature and emissivity constitutes an ill-conditioned inversion problem.

[0003] In related technologies, temperature and emissivity separation schemes typically employ algorithmic models based on empirical relationships. This approach primarily involves statistical regression analysis of typical ground object spectral data measured in laboratories to establish an empirical equation relating the minimum emissivity to the difference between the maximum and minimum emissivity. During inversion calculations, this empirical equation is used as a supplementary constraint in the thermal infrared radiation transfer equations to reduce the number of unknowns, thereby closing the equation set. Finally, algebraic operations are used to derive the surface temperature and emissivity values ​​for each spectral band.

[0004] However, with the improvement of hyperspectral sensor resolution and the aggravation of atmospheric noise in actual remote sensing environments, the simple empirical regression equations established by related technologies based on limited laboratory data will produce significant numerical calculation biases when processing actual hyperspectral remote sensing data containing complex noise. Summary of the Invention

[0005] This application provides a method and system for inverting thermal infrared hyperspectral temperature and emissivity based on a Gaussian prior model, which can reduce numerical calculation bias.

[0006] Firstly, this application provides a method for inverting temperature and emissivity in thermal infrared hyperspectral data based on a Gaussian prior model, applicable to a data processing system. The method includes: acquiring hyperspectral observation data and a pre-constructed emissivity sample library; calculating the mean vector and covariance matrix of the emissivity samples in the emissivity sample library to construct a Gaussian prior model of emissivity; constructing a Bayesian inversion model of the hyperspectral observation data based on the thermal infrared radiative transfer equation, observation noise distribution, and the Gaussian prior model; iteratively optimizing the temperature and emissivity in the Bayesian inversion model until convergence to obtain the optimal temperature and optimal emissivity; calculating the approximate covariance of the Bayesian inversion model under the optimal temperature and optimal emissivity conditions to obtain an uncertainty measure; and filtering the optimal temperature and optimal emissivity based on the uncertainty measure to obtain the target temperature and target emissivity.

[0007] In the above embodiments, the data processing system constructs a Gaussian prior model by calculating the mean vector and covariance matrix of the emissivity sample library. This Gaussian prior model is then combined with the thermal infrared radiation transfer equation and the observation noise distribution to construct a Bayesian inversion model. The target temperature and target emissivity are obtained through alternating iterative optimization and uncertainty metric screening. This allows the inversion process to integrate observation information and prior knowledge within a probabilistic statistical framework, reducing numerical calculation bias in ill-conditioned inversion problems.

[0008] In conjunction with some embodiments of the first aspect, in some embodiments, the step of constructing a Gaussian prior model of emissivity by calculating the mean vector and covariance matrix of emissivity samples in the emissivity sample library specifically includes: dividing the emissivity samples in the emissivity sample library according to a clustering algorithm to obtain multiple land cover category subsets; calculating the mean vector and initial covariance matrix of the emissivity samples within the land cover category subsets; performing regularization processing on the initial covariance matrix by combining a preset positive number and an identity matrix to obtain a regularized covariance matrix; and combining the mean vector and the regularized covariance matrix to determine the Gaussian prior model of emissivity.

[0009] In the above embodiments, the data processing system divides the emissivity samples into multiple land cover category subsets using a clustering algorithm, calculates the mean vector and initial covariance matrix of each subset, and performs regularization processing on the initial covariance matrix by combining preset positive numbers and identity matrices. This enables the constructed Gaussian prior model to characterize the spectral statistical characteristics of different land cover categories, while ensuring the numerical invertibility of the covariance matrix, thus improving the precision and computational stability of the prior constraints.

[0010] In conjunction with some embodiments of the first aspect, in some embodiments, the step of constructing a Bayesian inversion model of hyperspectral observation data based on the thermal infrared radiation transfer equation, observation noise distribution, and Gaussian prior model specifically includes: calculating the observation fitting term based on the hyperspectral observation data, the forward simulation data corresponding to the thermal infrared radiation transfer equation, and the covariance matrix of the observation noise distribution; calculating the prior constraint term based on the emissivity, the mean vector in the Gaussian prior model, and the covariance matrix in the Gaussian prior model; fusing the observation fitting term and the prior constraint term to obtain the objective function; and determining the probability estimation model containing the objective function as the Bayesian inversion model of the hyperspectral observation data.

[0011] In the above embodiments, the data processing system calculates the observation fitting term and the prior constraint term separately, and merges them into an objective function to construct a probability estimation model containing the objective function as a Bayesian inversion model. This allows the degree to which the observation data fitting residual and the emissivity deviate from the prior distribution to be jointly measured in a unified objective function, providing an optimization objective with clear probabilistic significance for subsequent alternating iterative optimization.

[0012] In some embodiments of the first aspect, before the step of iteratively optimizing the temperature and emissivity in the Bayesian inversion model until convergence to obtain the optimal temperature and emissivity, the method further includes: extracting reference band data with minimal atmospheric absorption from hyperspectral observation data; calculating the initial temperature by reversing the thermal infrared radiation transfer equation based on the reference band data and a preset initial emissivity constant; substituting the initial temperature into the thermal infrared radiation transfer equation and performing forward derivation based on the hyperspectral observation data to calculate the initial emissivity vector; and using the initial temperature and the initial emissivity vector as the initial values ​​for iterative optimization to constrain the optimization space of the Bayesian inversion model.

[0013] In the above embodiments, the data processing system extracts reference band data that is least affected by atmospheric absorption from hyperspectral observation data, obtains the initial temperature based on the reference band data and the initial emissivity constant inverted thermal infrared radiative transfer equation, and then derives the initial emissivity vector from the initial temperature in a forward manner. Using the initial temperature and the initial emissivity vector as the initial values ​​for iteration, the optimization space of the Bayesian inversion model is constrained within a physically reasonable range, which accelerates the iteration convergence and reduces the risk of getting trapped in local extrema.

[0014] In conjunction with some embodiments of the first aspect, in some embodiments, the step of calculating the initial temperature by reversing the thermal infrared radiation transfer equation based on reference band data and a preset initial emissivity constant specifically includes: obtaining the radiance, atmospheric upward radiation, and atmospheric downward radiation corresponding to the reference band data; constructing a surface thermal radiation term by combining the radiance, atmospheric upward radiation, atmospheric downward radiation, and the initial emissivity constant; and solving the surface thermal radiation term based on the inverse function of the Planck function to obtain the initial temperature.

[0015] In the above embodiments, the data processing system acquires the radiance, atmospheric up-going radiation, and atmospheric down-going radiation corresponding to the reference band data, constructs the surface thermal radiation term by combining it with the initial emissivity constant, and obtains the initial temperature by solving the surface thermal radiation term based on the inverse function of the Planck function. This ensures that the calculation of the initial temperature is completed in the band with the least atmospheric influence, thus guaranteeing the physical reliability of the initial temperature estimate.

[0016] In conjunction with some embodiments of the first aspect, in some embodiments, after the step of filtering the optimal temperature and optimal emissivity based on uncertainty metric to obtain the target temperature and target emissivity, the method further includes: obtaining the spatial neighborhood pixel set corresponding to the target temperature and target emissivity; calculating the spatial similarity weight between the target temperature and the temperatures of each neighboring pixel in the spatial neighborhood pixel set; and performing spatial joint correction on the target temperature and target emissivity based on the spatial similarity weight to obtain the final temperature and final emissivity.

[0017] In the above embodiments, the data processing system obtains the set of spatial neighboring pixels corresponding to the target temperature and the target emissivity, calculates the spatial similarity weight between the target temperature and the temperature of each neighboring pixel, and performs spatial joint correction based on the spatial similarity weight. It uses the spatial correlation between adjacent pixels to smooth the pixel-by-pixel inversion results, thereby suppressing the abnormal deviation of isolated pixels caused by noise.

[0018] In conjunction with some embodiments of the first aspect, in some embodiments, the step of performing spatial joint correction of the target temperature and target emissivity based on spatial similarity weights to obtain the final temperature and final emissivity specifically includes: constructing an adaptive bilateral filter based on spatial similarity weights, performing local weighted averaging on the target temperature to obtain the final temperature; substituting the final temperature back into the thermal infrared radiative transfer equation, and performing residual calculation based on the hyperspectral observation data of the corresponding pixel to obtain the radiance residual; and performing inverse compensation update on the target emissivity based on the radiance residual to obtain the final emissivity.

[0019] In the above embodiments, the data processing system constructs an adaptive bilateral filter based on spatial similarity weights to perform local weighted averaging of the target temperature to obtain the final temperature. Then, it substitutes the final temperature into the thermal infrared radiation transfer equation to calculate the radiance residual, and performs inverse compensation update of the target emissivity based on the radiance residual to obtain the final emissivity. This ensures that the radiant energy change after temperature correction is transmitted to the compensation update of emissivity, thus guaranteeing the physical consistency between temperature and emissivity.

[0020] In a second aspect, embodiments of this application provide a data processing system comprising: one or more processors and a memory; the memory is coupled to the one or more processors and is used to store computer program code, the computer program code including computer instructions, wherein the one or more processors invoke the computer instructions to cause the data processing system to perform the method described in the first aspect and any possible implementation thereof.

[0021] Thirdly, embodiments of this application provide a computer-readable storage medium including instructions that, when executed on a data processing system, cause the data processing system to perform the method described in the first aspect and any possible implementation thereof.

[0022] Fourthly, embodiments of this application provide a computer program product containing instructions that, when the computer program product is run on a data processing system, cause the data processing system to perform the method described in the first aspect and any possible implementation thereof.

[0023] Understandably, the data processing system provided in the second aspect, the computer storage medium provided in the third aspect, and the computer program product provided in the fourth aspect are all used to execute the methods provided in the embodiments of this application. Therefore, the beneficial effects they can achieve can be referred to the beneficial effects in the corresponding methods, and will not be repeated here.

[0024] One or more technical solutions provided in the embodiments of this application have at least the following technical effects or advantages:

[0025] 1. By employing a technique of constructing a Gaussian prior model based on an emissivity sample library and combining this Gaussian prior model with the thermal infrared radiative transfer equation and observation noise distribution to construct a Bayesian inversion model, the spectral shape constraint and inter-band covariance information of emissivity are incorporated into the inversion objective function and jointly optimized with the fitting residuals of the observation data. This effectively solves the problem of numerical calculation bias caused by prior constraints based on simple empirical regression equations in existing technologies when processing hyperspectral remote sensing data with complex noise. It realizes high-precision collaborative inversion of temperature and emissivity and endogenous quantification of uncertainty in inversion results within a probabilistic statistical framework.

[0026] 2. By employing the technique of extracting reference band data with minimal atmospheric absorption from hyperspectral observation data and using this reference band data to invert the thermal infrared radiation transfer equation to calculate the initial temperature and initial emissivity vector as the initial values ​​for iteration, the optimization space of the Bayesian inversion model is constrained within a physically reasonable neighborhood. This effectively solves the problem in existing technologies where high-dimensional nonlinear optimization problems lead to iteration divergence or convergence to local extrema due to initial values ​​deviating from the true solution. It achieves rapid convergence of the iterative solution process and stable acquisition of the global optimal solution.

[0027] 3. Because the technique of calculating spatial similarity weights based on the set of spatial neighboring pixels and performing joint spatial correction on the target temperature and target emissivity is adopted, the spatial correlation between adjacent pixels is introduced into the joint processing of temperature smoothing and emissivity compensation. This effectively solves the problem of spatial discontinuity jumps in the inversion results of adjacent pixels caused by observation noise when inverting pixel by pixel independently in the existing technology, and realizes the continuity constraint of the inversion results in the spatial dimension and the maintenance of the physical consistency of temperature-emissivity. Attached Figure Description

[0028] Figure 1 This is a schematic flowchart of a method for inverting thermal infrared hyperspectral temperature and emissivity based on a Gaussian prior model in an embodiment of this application.

[0029] Figure 2 This is another schematic diagram of the process for retrieving thermal infrared hyperspectral temperature and emissivity based on a Gaussian prior model in this application embodiment;

[0030] Figure 3 This is a schematic diagram of the physical device structure of a data processing system in an embodiment of this application. Detailed Implementation

[0031] The terminology used in the following embodiments of this application is for the purpose of describing particular embodiments only and is not intended to be limiting of this application. As used in the specification of this application, the singular expressions “a,” “an,” “the,” “the,” and “this” are intended to include the plural expressions as well, unless the context clearly indicates otherwise. It should also be understood that the term “and / or” as used in this application refers to any or all possible combinations including one or more of the listed items.

[0032] Hereinafter, the terms "first" and "second" are used for descriptive purposes only and should not be construed as implying or suggesting relative importance or implicitly indicating the number of indicated technical features. Thus, a feature defined as "first" or "second" may explicitly or implicitly include one or more of that feature, and in the description of the embodiments of this application, unless otherwise stated, "multiple" means two or more.

[0033] Thermal infrared hyperspectral remote sensing inversion involves multiple technical steps, forming a complete data processing chain. The following explains the technical concepts involved in each step according to the execution order of this chain.

[0034] The input for the inversion consists of two parts.

[0035] The first part consists of thermal infrared hyperspectral observation data, which is the surface thermal radiation data collected by thermal infrared hyperspectral sensors in continuous narrow bands (32 to 128 bands, band width 0.05-0.2μm) within an 8-14μm atmospheric window. Each pixel corresponds to a spectral vector composed of a set of multi-band radiance values, which is the observation to be interpreted by the inversion.

[0036] The second part is the emissivity sample library, which is a standard data set of laboratory measurements of thermal infrared emissivity spectra of various land cover types collected in advance (typical sources include the USGS spectral library, JHU spectral library, and ASU thermal emission spectral library). Before use, it needs to be convolved to the band settings of the target sensor through spectral resampling to construct prior constraints.

[0037] Since the joint inversion of temperature and emissivity is an ill-conditioned problem (with more unknowns than observed equations), prior constraints need to be introduced to narrow the solution space. The prior constraints in this application are derived from statistical modeling of an emissivity sample library: the mean vector and covariance matrix of the emissivity spectra in the sample library are calculated, describing the emissivity distribution as a multivariate normal distribution, thus forming a Gaussian prior model. This model encodes the prior knowledge of "what a reasonable emissivity spectrum should look like," providing regularization constraints for subsequent inversion.

[0038] Combining prior constraints with a physical model requires a forward model describing "how observational data is generated from surface parameters." This forward model is the thermal infrared radiative transfer equation: for each band, the radiance received by the sensor is the superposition of three components: the component of surface-emitted radiation attenuated by the atmosphere, the component of atmospheric upward radiation, and the component of atmospheric downward radiation attenuated by the atmosphere after reflection from the surface. This equation relates surface temperature T and emissivity ε to atmospheric transmittance τ, atmospheric upward radiation, and atmospheric downward radiation, which can be predetermined by the atmospheric model, providing a deterministic mapping from surface parameters to sensor observations.

[0039] Building upon the prior constraints provided by the Gaussian prior model and the positive physical mapping provided by the radiative transfer equation, this application unifies the two within the framework of Bayes' theorem, constructing a Bayesian inversion model: the likelihood function measures the degree of fit of the current parameters to the observed data (determined by the observation residuals and noise covariance matrix), the prior probability density is given by the Gaussian prior model, and the product of the two is proportional to the posterior probability density. Taking the negative logarithm of the posterior probability yields the objective function to be minimized, which consists of two parts: an observation fit term and a prior constraint term, balancing data fit and prior consistency.

[0040] The joint optimization of the objective function with respect to temperature and emissivity is a nonlinear problem. This application employs an alternating iterative optimization strategy: in each iteration, one variable is fixed while the other is optimized. When the temperature is fixed, the objective function with respect to emissivity is a quadratic convex function, and an analytical solution can be obtained by solving a system of linear equations. When the emissivity is fixed, the objective function with respect to temperature is a one-dimensional nonlinear function, and the optimal value can be searched within the physical temperature range using the golden section method or Newton's method. This process is repeated alternately until convergence, yielding the optimal temperature and optimal emissivity.

[0041] After obtaining the optimal solution through inversion, it is necessary to evaluate the reliability of the result. This application achieves this evaluation through uncertainty measurement: based on the Laplace approximation, the Hessian matrix of the objective function is calculated at the optimal solution, and its inverse matrix is ​​the posterior approximate covariance matrix; the diagonal elements correspond to the posterior variance of each parameter, thereby obtaining the temperature inversion standard deviation and the emissivity inversion variance of each band, which serve as pixel-level quality indicators.

[0042] The aforementioned pixel-by-pixel inversion does not utilize the spatial correlation between adjacent pixels, and lacks spatial constraints when noise causes abnormal results for individual pixels. To address this, this application introduces a spatial joint correction: a spatial neighborhood pixel set is formed by selecting surrounding pixels centered on the current pixel according to a preset window size, and spatial similarity weights are calculated based on the temperature differences between adjacent pixels using a Gaussian kernel function (the smaller the temperature difference, the greater the weight, which has the same mathematical form as the range kernel of the bilateral filter). This weight is used to perform a local weighted average of the temperature and a physically consistent inverse compensation for the emissivity, thus suppressing spatial jumps caused by noise while maintaining the boundaries of ground features.

[0043] This processing link is suitable for applications requiring high-precision surface temperature and emissivity parameters, such as geological and mineral identification, urban heat island monitoring, and agricultural growth assessment.

[0044] The method provided in this embodiment is described in detail below. Please refer to [link / reference]. Figure 1 This is a schematic diagram of a method for inverting thermal infrared hyperspectral temperature and emissivity based on a Gaussian prior model in this application. The entire inversion chain can be divided into the following six logical stages:

[0045] Phase 1: Data preparation, namely: S101, acquiring hyperspectral observation data and a pre-constructed emissivity sample library.

[0046] Before performing the inversion calculation, the data processing system needs to acquire two types of basic input data:

[0047] 1. Hyperspectral observation data: Pixel-by-pixel radiance spectral vectors acquired by a thermal infrared hyperspectral sensor within the thermal infrared atmospheric window band. For a pixel with spatial location (r, c), its observation data is represented as an N-dimensional radiance vector: ;in Let be the observed radiance value of the i-th band (i=1,2,…,N), in units of W / (m²·sr·μm), and N typically ranges from 32 to 128.

[0048] At the same time, each pixel corresponds to a set of atmospheric parameter data, including:

[0049] Atmospheric transmittance vector: ;

[0050] Atmospheric upward radiation vector: ;

[0051] Downward atmospheric radiation vector: Atmospheric parameter data can be obtained by combining atmospheric radiative transfer models (such as MODTRAN or 6S) with synchronous sounding data or atmospheric reanalysis data.

[0052] 2. Emissivity Sample Library: A pre-collected set of emissivity spectral data for various land cover types in the thermal infrared band. Each emissivity spectrum in this dataset undergoes spectral resampling to match the band response function of the target sensor. The emissivity sample library is stored in matrix form. ,in , where is the m-th emissivity spectral vector (m=1,2,…,M), M is the number of samples, and N is the number of bands.

[0053] In step S101, the data processing system can simultaneously acquire two types of input data. The first type is hyperspectral observation data, which is obtained by scanning and imaging the target area using a thermal infrared hyperspectral sensor. The second type is an emissivity sample library. The data processing system extracts emissivity spectra covering various land cover types, such as rocks, minerals, soils, vegetation, water bodies, and artificial materials, from standard spectral databases (typical sources include the USGS spectral library, the JHU spectral library, and the ASU thermal emission spectral library). Each spectrum is then convolved and resampled using the sensor's spectral response function to ensure its band settings are consistent with the hyperspectral observation data.

[0054] In some embodiments, some pixels in the hyperspectral observation data may have invalid radiance values ​​due to cloud cover or sensor saturation. To address this, the data processing system performs a pixel-by-pixel data validity determination after reading the hyperspectral observation data. The validity determination function V(r, c) is defined as follows:

[0055] when (r,c)∈[ , The condition holds true for all i=1,2,…,N, and the number of consecutive zero bands is less than a preset threshold. When the value is 1, V(r,c) = 1 (valid pixel); otherwise, V(r,c) = 0 (invalid pixel).

[0056] in and This represents the boundary value of the physically reasonable radiance range. In subsequent inversion calculations, the data processing system only performs temperature and emissivity inversion calculations for pixels with V(r,c)=1.

[0057] The second stage is prior modeling, namely: S102, calculating the mean vector and covariance matrix of the emissivity samples in the emissivity sample library, and constructing a Gaussian prior model of emissivity.

[0058] To provide reasonable regularization constraints for the ill-conditioned inversion problem, the data processing system performs statistical parameter calculations on the M emissivity spectra in the sample database:

[0059] 1. Calculation of the mean vector μ: Take the arithmetic mean of the M spectra in each band to obtain an N-dimensional vector: That is, the i-th component of the mean vector is: .

[0060] 2. Calculation of the covariance matrix Σ (for simplicity, the covariance matrix will be abbreviated as Σ in the following text; please note the difference between it and the summation symbol): Calculate the outer product of the deviation vectors of the sample spectra and the mean vector, and take the average to obtain an N×N matrix: The (i,j)th element of the covariance matrix is: Among them, the diagonal elements The variance characterizing the emissivity of the i-th band, with off-diagonal elements. (i≠j) represents the covariance between the emissivity of the i-th band and the j-th band.

[0061] Based on the above parameters, a Gaussian prior model with a multivariate normal distribution is constructed, and the prior probability density function of the emissivity spectral vector ε is: That is, ε ~ N(μ, Σ) (following a normal distribution); this probability density reaches its maximum value at the mean vector μ, and the Mahalanobis distance defined by the covariance matrix Σ is used in the direction away from the mean vector. attenuation: .

[0062] In some embodiments, the Gaussian prior model can be constructed in various ways. Optionally, the data processing system calculates a uniform mean vector and covariance matrix for all samples in the emissivity sample library to obtain a global Gaussian prior model N(μ,Σ). This method is suitable for scenarios where the land cover types in the target area are unknown or have a high degree of mixing. Optionally, the data processing system uses auxiliary land cover classification information to divide the emissivity samples into K subsets according to category labels. , ,…, } Calculate the mean vector for each subset. Covariance Matrix (k=1,2,…,K), resulting in K categories of Gaussian prior models N( , During inversion, the corresponding prior model is selected based on the land cover category of the pixel. This method is suitable for scenarios where the land cover type of the target area is known and can be effectively classified. It is understood that other probabilistic models, such as Gaussian mixture models, can also be used to model the statistical distribution of the emissivity spectrum; this is not limited here.

[0063] The third stage: model building, namely: S103, based on the thermal infrared radiation transfer equation, observation noise distribution and Gaussian prior model, construct a Bayesian inversion model of hyperspectral observation data.

[0064] The thermal infrared radiation transfer equation is a physical equation describing the process of surface thermal radiation being transported through the atmosphere to the sensor. For the i-th band, the forward simulation expression for the sensor's entrance pupil radiance is: , i=1,2,…,N, where Let T be the surface emissivity of the i-th band, and T be the surface temperature. Let be the atmospheric transmittance of the i-th band. This represents the downward atmospheric radiation in the i-th band. Let B(λi,T) be the atmospheric upward radiation of the i-th band, and let B(λi,T) be the radiance of the Planck function at the center wavelength λi and temperature T. Where c1 = 1.191 × 10 8 W·μm 4 ·m⁻²·sr⁻¹ is the first radiation constant, c² = 1.4388 × 10⁻¹ 4 μm·K is the second radiation constant.

[0065] All N bands of forward simulation are written in vector form: The observation noise distribution is assumed to be a multivariate normal distribution with a mean of zero, i.e., the observation noise vector n ~ N(0, ),in The observation noise covariance matrix is ​​an N×N diagonal matrix, and its i-th diagonal element is... = Let be the noise variance of the i-th band. The observation data model is: .

[0066] Therefore, given a temperature T and an emissivity ε, the likelihood function of the observed data is: According to Bayes' theorem, the posterior probability density function is proportional to the product of the likelihood function and the prior probability density function: .

[0067] Taking the negative logarithm of the posterior probability density function and ignoring constant terms independent of the parameters, we obtain the objective function. The observation fitting term is: This is used to measure how well the current parameter estimates fit the observed data; the prior constraint term is... , is used to measure the degree to which the current emissivity estimate deviates from the prior distribution center, i.e., the square of the Mahalanobis distance.

[0068] The maximum a posteriori probability estimation model that minimizes the objective function J(ε,T) is called the Bayesian inversion model. .

[0069] In some embodiments, the data processing system refines the Bayesian inversion model construction process in step S103. When When the matrix is ​​diagonal, the calculation of the observation fitting term simplifies to dividing the square of the residuals for each band by the sum of the corresponding band noise variances: When the data processing system calculates the prior constraint terms, it has already performed Cholesky decomposition Σ=L on the prior covariance matrix. In this case, the calculation can be performed by solving Lx = ε − μ. = x is completed, avoiding explicit inverse calculation.

[0070] In some embodiments, the fusion of the observation fit term and the prior constraint term can be achieved in various ways. Optionally, the data processing system directly adds the observation fit term and the prior constraint term as the objective function, which corresponds to the standard form of Bayesian maximum a posteriori probability estimation. Optionally, the data processing system introduces a weighting coefficient α before the prior constraint term, and the objective function becomes: The strength of the prior constraint can be controlled by adjusting the value of α. When α > 1, the prior constraint is strengthened, which is suitable for high-noise scenarios. When α < 1, the prior constraint is weakened, which is suitable for scenarios with a high observation signal-to-noise ratio. The determination of α can be achieved through the L-curve method or the generalized cross-validation method. It is understandable that the objective function can also be constructed by adding the log-likelihood function and the log-prior function, which is not limited here.

[0071] In some embodiments, the mean vector μ of the Gaussian prior model may deviate significantly from the true emissivity spectra of certain land features in the target area (e.g., the target area contains new land feature types not included in the prior sample library), causing the prior constraint terms to excessively pull the inversion results in an incorrect direction. To address this, the data processing system monitors the ratio of the prior constraint terms to the observation fitting terms during the iterative solution process. .when It continuously increases in multiple iterations and If the rate of decrease fails to continue, the data processing system halves α, making the objective function rely more on the fit of the observed data, thus avoiding misleading inversion results due to mismatched prior knowledge.

[0072] The fourth stage is iterative optimization, namely: S104, alternatingly iterating and optimizing the temperature and emissivity in the Bayesian inversion model until convergence, to obtain the optimal temperature and optimal emissivity.

[0073] Since the joint optimization of the objective function J(ε,T) with respect to temperature and emissivity is a nonlinear problem, the data processing system employs a coordinate descent strategy for alternating iterative solutions. Let the temperature estimate in the k-th iteration be... The emissivity estimate is The iterative process is as follows:

[0074] Sub-step 1: Fix the temperature and optimize the emissivity.

[0075] Fix the current temperature At that time, the i-th band component of the positively simulated radiance affects the emissivity. It is a linear function. Define an N×N diagonal matrix. Its i-th diagonal element: Define the offset vector b, whose i-th component is: The forward simulation vector can then be uniformly written in matrix form as follows: .

[0076] Substituting this linear forward model into the objective function, taking the partial derivative with respect to ε and setting it to zero, yields the normal equation: Define the coefficient matrix. , is an N×N symmetric positive definite matrix; the vector on the right-hand side The data processing system performs Cholesky decomposition on A to obtain a lower triangular matrix L such that A = L. Then, Lz=y is solved using the previous substitution method and the back substitution method. ε=z, thus obtaining the optimal emission rate vector for the current iteration step. .

[0077] Sub-step 2: Fix the emissivity and optimize the temperature.

[0078] Fixed current emission rate When, the objective function J( (T) is a one-dimensional nonlinear function with respect to temperature T. The data processing system operates within the physical temperature range [ , Within ] (e.g., 200K to 400K) or a dynamic search range [ −ΔT, Search within +ΔT] for the temperature that minimizes the objective function, such as: The search method can employ either the golden ratio or Newton's method. When using Newton's method, the temperature update formula is: Where J'(T) = ∂J / ∂T is the first derivative of the objective function with respect to temperature, and J''(T) = ∂²J / ∂T² is the second derivative. When |J'( The temperature search will terminate when | < δT, where δT is a preset accuracy threshold (e.g., 0.001K).

[0079] Convergence Criterion: The data processing system calculates the temperature convergence index after each iteration. Emissivity convergence index: When Δ is satisfied simultaneously <θΔT and Δ The iteration terminates when θΔT < θΔε, where θΔT is the temperature convergence threshold (e.g., 0.01K) and θΔε is the emissivity convergence threshold (e.g., 10⁻⁻⁶). 4 The current temperature and emissivity are obtained as the optimal temperature T and optimal emissivity ε.

[0080] In some embodiments, alternating iterative optimization can be performed in various ways. Optionally, the data processing system applies physical range constraints to the solved emissivity vector in the emissivity update sub-step, such as: The components less than 0 are truncated to 0, and the components greater than 1 are truncated to 1. The constrained emissivity vector is then used in the temperature update sub-step. It is understood that other optimization algorithms, such as gradient descent or quasi-Newton methods, can also be used to perform alternating iterative optimization; this is not limited here.

[0081] In some embodiments, the iterative process may experience slow convergence or oscillation-induced non-convergence due to the initial values ​​deviating too far from the true solution. To address this, the data processing system sets a maximum limit on the number of iterations. (e.g., 50 times), when the number of iterations reaches If the convergence condition is not yet met, the data processing system records the sequence of objective function values ​​for the current iteration. , ,…, } Determine whether the objective function value shows a monotonically decreasing trend; if the objective function value oscillates (i.e., there exists a k such that...), > The data processing system narrows the temperature search range to ±5K of the current temperature estimate and reduces the emissivity update step size to 0.5 times the original step size, continuing to iterate in a damped manner until convergence or the upper limit of the second iteration is reached.

[0082] Phase 5: Uncertainty assessment, namely: S105, calculate the approximate covariance of the Bayesian inversion model under optimal temperature and optimal emissivity conditions to obtain the uncertainty measure.

[0083] After obtaining the optimal temperature T and optimal emissivity ε through inversion, the reliability of the results needs to be evaluated. The data processing system calculates the Hessian matrix H of the objective function at the optimal solution based on the Laplace approximation. Since the alternating iterative optimization uses a coordinate descent strategy to solve for temperature and emissivity independently, a block diagonal approximation is used here to ignore the cross-partial derivatives between temperature and emissivity, calculating only the emissivity-emissivity block and the temperature-temperature scalar block.

[0084] For the emissivity part, the emissivity-emissivity block of the Hessian matrix is: ,in Let be the Jacobian matrix of the forward model with respect to emissivity at the optimal solution.

[0085] Data processing system Inversely, the posterior covariance matrix of the emissivity is obtained as follows: The diagonal elements of this matrix represent the posterior variance of the emissivity for each band: The standard uncertainty of the emissivity for each band is the square root of the posterior variance, which is: = .

[0086] For the temperature component, the data processing system calculates the second-order partial derivative of the objective function with respect to temperature as follows: The posterior variance of temperature is the reciprocal of the second-order partial derivative, i.e.: The standard uncertainty of temperature is the square root of the posterior variance, i.e.: .

[0087] The data processing system calculates the standard uncertainty of the emissivity for each band. The standard uncertainty of (i=1,2,…,N) and temperature As an uncertainty measure, this measure provides a quantitative evaluation of the confidence level of the inversion results: a small uncertainty indicates that the inversion results are statistically reliable; a large uncertainty suggests that the inversion results may be significantly affected by observation noise or prior model bias, requiring comprehensive judgment in conjunction with other auxiliary information.

[0088] In some embodiments, the approximate covariance can be calculated in several ways. Optionally, the data processing system can directly calculate the posterior covariance matrix of emissivity by matrix inversion. =( ⁻¹, when the number of bands N is small (e.g., N≤64), the computational complexity of this method is acceptable. Optionally, when the number of bands N is large, the data processing system uses the Woodbury matrix identity to transform the inversion of the N×N matrix into an equivalent low-dimensional computation to reduce computational complexity. Understandably, Monte Carlo sampling methods can also be used to draw samples from the posterior distribution to estimate the posterior statistic; this is not limited here.

[0089] In some embodiments, a Hessian matrix may exist. The condition number is too large (e.g., more than 10). 6 This leads to a decrease in the numerical precision of matrix inversion. To address this, the data processing system... Before finding the inverse, calculate its condition number κ. )= / ,in and They are respectively The maximum and minimum eigenvalues; when κ( >10 6 At that time, the data processing system Perform truncated singular value decomposition, and select singular values ​​smaller than the maximum singular value by 10⁻. 6 The singular values ​​are truncated to zero, and the pseudo-inverse of the truncated matrix is ​​used to approximate the direct matrix inversion, so as to ensure the numerical stability of the posterior covariance matrix estimation.

[0090] Phase 6: Output of results, namely: S106, Based on the optimal temperature, optimal emissivity and uncertainty measure, generate temperature map, emissivity map and uncertainty map of the target area.

[0091] After performing the inversion and uncertainty assessment procedures S101 to S105 on each valid pixel (V(r,c)=1) in the target area, the data processing system organizes the inversion results into a three-dimensional data cube. Specifically:

[0092] 1. Temperature Map Generation: The data processing system arranges the optimal temperature T*(r,c) of each effective pixel into a two-dimensional matrix according to the spatial row and column positions (r,c) of the original pixels, generating a surface temperature distribution map with the same spatial resolution as the original hyperspectral image. Invalid pixel positions are filled with preset invalid value markers (e.g., −9999). The unit of each pixel value in the temperature map is Kelvin (K).

[0093] 2. Generation of Emissivity Map: The data processing system arranges the optimal emissivity ε*(r,c) of each effective pixel according to its spatial row and column positions, generating a three-dimensional emissivity data cube. The spatial dimension is the number of rows × the number of columns in the image, and the spectral dimension is the number of bands N. The value of each pixel in the emissivity map for each band is a dimensionless emissivity value (ranging from 0 to 1).

[0094] 3. Generation of uncertainty map: The data processing system will generate the standard uncertainty of temperature for each effective pixel. Arrange (r,c) into a two-dimensional matrix to generate a temperature uncertainty distribution map; calculate the standard uncertainty of emissivity for each effective pixel in each band. Arrange (r,c) into a three-dimensional data cube to generate an emissivity uncertainty distribution map. The uncertainty map has the same spatial resolution and number of bands as the temperature map and the emissivity map.

[0095] In some embodiments, the inversion results can be output in multiple ways. Optionally, the data processing system stores the temperature map, emissivity map, and uncertainty map as raster files in GeoTIFF format, retaining the geographic coordinate system, projection parameters, and spatial resolution information of the original imagery to facilitate subsequent spatial analysis and visualization in a geographic information system. Optionally, the data processing system stores the above results as scientific data files in HDF5 or NetCDF format, and attaches metadata to record inversion parameters (such as prior model type, number of iterations, convergence threshold, etc.) to support data traceability and repeatability verification. It is understood that the result data can also be stored in the ENVI standard format or a custom binary format, which is not limited here.

[0096] In some embodiments, the target area may have a large image size (e.g., containing millions of pixels), resulting in excessively long total computation time for pixel-by-pixel inversion. To address this, the data processing system employs a spatial block processing strategy in the result organization stage (S106). The system divides the entire image into multiple sub-blocks along the row direction. Each sub-block contains a fixed number of pixel data rows (e.g., 256 rows). Each sub-block independently executes the inversion process from S101 to S105 and then writes the results to its corresponding output data area. This block processing strategy reduces the memory requirements for a single computation, while the computations between each sub-block are independent, allowing for parallel acceleration using multi-core processors or computing clusters.

[0097] The above describes the basic workflow of the thermal infrared hyperspectral temperature and emissivity inversion method based on a Gaussian prior model. In this basic workflow, the method for setting the initial iteration values ​​and the spatial consistency processing of the inversion results are not yet specified in detail. In actual thermal infrared hyperspectral remote sensing data processing, the rationality of the initial iteration values ​​directly affects the convergence speed and results of the alternating iterative optimization. If the initial iteration values ​​deviate too far from the true solution, the iterative process will converge slowly or converge to a local extremum. Furthermore, since pixel-by-pixel independent inversion does not utilize the spatial correlation between adjacent pixels, spatial discontinuities and jumps may occur in the inversion results of adjacent pixels when observation noise is high.

[0098] The method provided in this embodiment will now be described in more detail. Please refer to [link / reference]. Figure 2 This is another flowchart illustrating the thermal infrared hyperspectral temperature and emissivity inversion method based on a Gaussian prior model in this application.

[0099] S201. Acquire hyperspectral observation data and a pre-constructed emissivity sample library.

[0100] Refer to step S101, which will not be repeated here.

[0101] S202. Calculate the mean vector and covariance matrix of the emissivity samples in the emissivity sample library, and construct a Gaussian prior model of emissivity.

[0102] Refer to step S102, which will not be repeated here.

[0103] In step S202, the data processing system uniformly calculates the mean vector and covariance matrix for all samples in the emissivity sample library to construct a Gaussian prior model. When the land cover types included in the emissivity sample library are diverse (e.g., simultaneously containing the spectra of quartz minerals, vegetation, and water bodies), a single Gaussian distribution may not accurately characterize the emissivity spectral distribution with multi-peak characteristics. To address this, the following process describes the specific implementation steps for constructing a refined Gaussian prior model through clustering and classification modeling.

[0104] In some embodiments, the data processing system refines the Gaussian prior model construction process in step S202. Specifically, the data processing system divides the emissivity samples in the emissivity sample library according to a clustering algorithm to obtain multiple land cover category subsets; calculates the mean vector and initial covariance matrix of the emissivity samples within the land cover category subsets; regularizes the initial covariance matrix by combining preset positive numbers and identity matrices to obtain a regularized covariance matrix; and combines the mean vector and the regularized covariance matrix to determine the Gaussian prior model of emissivity.

[0105] Clustering algorithms refer to unsupervised grouping of M emissivity spectra in an emissivity sample database based on spectral similarity. Commonly used clustering algorithms include K-means clustering and hierarchical clustering. The clustering result divides the M samples into K non-overlapping subsets. A land cover category subset is a set of emissivity samples belonging to the same category after clustering, where samples within each subset have similar spectral shapes and characteristics. The initial covariance matrix is ​​an N×N matrix calculated using the standard sample covariance formula within each land cover category subset. The preset positive number δ is a small positive real number used for regularization, typically ranging from 0.0001 to 0.01. The identity matrix I is an N×N identity matrix. The regularized covariance matrix is ​​the matrix obtained by adding δI to the initial covariance matrix, i.e. This operation ensures that the minimum eigenvalue of the covariance matrix is ​​not less than δ, satisfying the positive definiteness and numerical invertibility conditions.

[0106] Specifically, the data processing system first selects the clustering algorithm and the number of clusters K. Taking the K-means clustering algorithm as an example, the data processing system uses M N-dimensional emissivity spectra as input samples for the K-means algorithm, randomly initializes K cluster centers, and iteratively performs sample allocation (assigning each spectrum to the cluster center with the closest Euclidean distance) and center update (recalculating the mean spectrum of each cluster) until the cluster centers no longer change, thus obtaining K subsets of land cover categories. The number of categories K can be determined with the help of clustering evaluation indicators such as the silhouette coefficient or the elbow rule, and is usually in the range of 3 to 10. The data processing system then processes the k-th subset of land cover categories (containing...) Calculate the mean vector from the sample. (Taking the arithmetic mean of each band for all samples within the subset) and the initial covariance matrix (Calculated according to the sample covariance formula). When the number of samples in the subset When the number of bands is less than N, the initial covariance matrix Since it is a rank-deficient matrix, it cannot be directly inverted. Data processing systems... Adding δI yields the regularized covariance matrix. = +δI, ensuring its positive definiteness and invertibility. The data processing system will divide the K categories { , The combination is determined to be a Gaussian prior model of emissivity, and the corresponding prior parameters are selected according to the class assignment of the pixels in the subsequent inversion.

[0107] In some embodiments, clustering and regularization can be performed in various ways. Optionally, the data processing system uses the K-means clustering algorithm to divide the emissivity samples, setting K=5, performing 100 different random initializations, and selecting the clustering result with the smallest objective function value as the final division. After calculating the mean vector and initial covariance matrix for each subset, regularization is performed with δ=0.001. Optionally, the data processing system uses a hierarchical clustering algorithm (agglomerated, using Ward distance as the merging criterion), truncating the dendrogram at an appropriate distance threshold to obtain K subsets, and performing the same mean and covariance calculation and regularization on each subset. This method does not require pre-setting the K value and is suitable for scenarios where there is no prior knowledge of the number of categories. It is understood that the expectation-maximization algorithm of Gaussian mixture models can also be used to simultaneously complete clustering and parameter estimation, which is not limited here.

[0108] In some embodiments, there may be instances where the number of samples in certain subsets of land cover categories is too small (e.g., only 2 to 3 samples), leading to severely inaccurate covariance matrix estimation. To address this, the data processing system performs a sample size test after calculating the initial covariance matrix for each subset. When the variance is less than 2N, the data processing system replaces the covariance matrix of the subset with the result of regularization of the global covariance matrix (the covariance matrix calculated based on all M samples), and only retains the mean vector of the subset itself, so as to avoid the adverse effect of covariance estimation bias on the inversion result under small sample size.

[0109] S203. Based on the thermal infrared radiation transfer equation, observation noise distribution, and Gaussian prior model, a Bayesian inversion model for hyperspectral observation data is constructed.

[0110] Refer to step S103, which will not be repeated here.

[0111] In step S203, the data processing system constructs the overall framework of the Bayesian inversion model. The following process further describes the specific composition of the objective function in this Bayesian inversion model, including the specific calculation methods of the observation fitting term and the prior constraint term, as well as the fusion method between the two.

[0112] In some embodiments, the data processing system refines the Bayesian inversion model construction process in step S203. Specifically, the data processing system calculates the observation fitting term based on the hyperspectral observation data, the forward simulation data corresponding to the thermal infrared radiation transfer equation, and the covariance matrix of the observation noise distribution; calculates the prior constraint term based on the emissivity, the mean vector in the Gaussian prior model, and the covariance matrix in the Gaussian prior model; fuses the observation fitting term and the prior constraint term to obtain the objective function; and determines the probability estimation model containing the objective function as the Bayesian inversion model for the hyperspectral observation data.

[0113] Here, the forward simulation data refers to the vector M(ε, T) composed of simulated entrance pupil radiance values ​​of N band sensors, calculated using the thermal infrared radiative transfer equation under given temperature T and emissivity ε. The observation fitting term refers to the hyperspectral observation data vector. With forward simulation data vector The difference in the observation noise covariance matrix The quadratic form of the inverse matrix weighted by is expressed as follows: This term measures how well the current parameter estimates fit the observed data. The prior constraint term is the quadratic form of the difference between the emissivity vector ε and the mean vector μ of the Gaussian prior model, weighted by the inverse of the prior covariance matrix Σ. Its mathematical expression is: This term measures the degree to which the current emissivity estimate deviates from the prior distribution center, i.e., the square of the Mahalanobis distance. The objective function J(ε, T) is the sum of the observation fit term and the prior constraint term. A probabilistic estimation model refers to a model that minimizes the objective function. The maximum a posteriori probability estimation model is used as the criterion, and the optimal solution of this model corresponds to the maximum point of the posterior probability density function.

[0114] Specifically, when calculating the observation fitting term, the data processing system first uses the thermal infrared radiative transfer equation to calculate the forward simulated radiance band by band based on the current temperature estimate T and emissivity estimate ε. The simulated values ​​from N bands are combined to form a positive simulation vector M(ε, T). The data processing system calculates the residual vector. Then calculate the observation fitting term. .when When the matrix is ​​diagonal, the calculation simplifies to dividing the squared residuals of each band by the sum of the corresponding band noise variances. When calculating the prior constraint terms, the data processing system calculates the bias vector d = ε - μ, and then calculates the prior constraint terms. When Cholesky decomposition has been performed on the prior covariance matrix. In this case, the calculation can be performed by solving Lx=d. Complete, avoiding explicit inversion. The data processing system adds the observed fitting term and the prior constraint term to obtain the objective function. The probability estimation model that uses minimizing J(ε,T) as the solution criterion is determined to be the Bayesian inversion model.

[0115] In some embodiments, the fusion of the observation fitting term and the prior constraint term can be achieved in various ways. Optionally, the data processing system directly adds the observation fitting term and the prior constraint term as the objective function. This method corresponds to the standard form of Bayesian maximum a posteriori probability estimation and is suitable for the observation noise covariance matrix. This is suitable for scenarios where both the estimation of the prior covariance matrix Σ and the objective function are relatively accurate. Optionally, the data processing system introduces a weighting coefficient α before the prior constraint term, and the objective function becomes... The strength of the prior constraint can be controlled by adjusting the value of α. When α > 1, the prior constraint is strengthened, which is suitable for high-noise scenarios. When α < 1, the prior constraint is weakened, which is suitable for scenarios with a high observation signal-to-noise ratio. The determination of α can be achieved through the L-curve method or the generalized cross-validation method. It is understandable that the objective function can also be constructed by adding the log-likelihood function and the log-prior function, which is not limited here.

[0116] In some embodiments, the mean vector μ of the Gaussian prior model may deviate significantly from the true emissivity spectra of certain land features in the target area (e.g., the target area contains new land feature types not included in the prior sample library). This can cause the prior constraint term to excessively pull the inversion results in an incorrect direction. To address this, the data processing system monitors the ratio of the prior constraint term to the observation fitting term during the iterative solution process. When this ratio continuously increases in multiple iterations and the residual norm of the observation fitting term fails to continuously decrease, the data processing system reduces the weight of the prior constraint term (halving α), making the objective function rely more on the observation data fitting and avoiding misleading inversion results due to mismatched prior knowledge.

[0117] After constructing the Bayesian inversion model, the data processing system needs to determine reasonable initial values ​​for subsequent alternating iterative optimization. The reasonableness of the initial values ​​directly affects the convergence path and convergence speed of the coordinate descent method. Steps S204 to S207 describe the specific process of determining the initial values ​​based on physical constraints.

[0118] S204. Extract the reference band data that is least affected by atmospheric absorption from hyperspectral observation data.

[0119] The reference band refers to one or more bands within the thermal infrared atmospheric window that have the highest atmospheric transmittance. The observed radiance of these bands is minimally affected by atmospheric absorption and emission, and can directly reflect surface thermal radiation information. Reference band data refers to the observed radiance value corresponding to the reference band and its associated atmospheric parameters (atmospheric transmittance, upward atmospheric radiation, and downward atmospheric radiation). Minimal atmospheric absorption means that the atmospheric transmittance of this band is the highest among all bands, typically greater than 0.9.

[0120] Specifically, after constructing the Bayesian inversion model, the data processing system needs to determine a set of physically reasonable initial values ​​for iteration to initiate subsequent alternating iterations for optimization. The system reads atmospheric transmittance values ​​for each band from the atmospheric parameter data corresponding to the hyperspectral observation data, sorts the bands in descending order of transmittance, and selects the band with the highest transmittance as the reference band. When multiple bands have transmittance values ​​close to the highest value (difference less than 0.01), the system selects the band whose center wavelength is located at the center of the atmospheric window as the reference band to reduce the influence of residual atmospheric absorption at the window's edge. The system extracts the observed radiance value corresponding to the reference band, along with the band's atmospheric transmittance, upward atmospheric radiation, and downward atmospheric radiation, to form the reference band data.

[0121] In some embodiments, reference band data can be extracted in multiple ways. Optionally, the data processing system pre-sets the band number of the reference band in the sensor configuration file, and directly extracts the corresponding observed radiance and atmospheric parameters according to the fixed band number when processing each pixel. This method is suitable for situations where the atmospheric transmittance of each band of the sensor is relatively stable and does not change significantly with atmospheric conditions. Optionally, the data processing system independently calculates the atmospheric transmittance of each band for each pixel, and dynamically selects the band with the highest transmittance as the reference band based on the atmospheric state of the current pixel. This method can adapt to situations where the transmittance ranking of each band changes under different atmospheric conditions. The data processing system sorts the N band transmittances of the current pixel in descending order, takes the first band in the ranking as the reference band, and extracts its corresponding data. It is understood that a multi-band weighted average method can also be used to construct equivalent reference band data, which is not limited here.

[0122] In some embodiments, there may be abnormally high atmospheric water vapor content above the target area, resulting in atmospheric transmittance below 0.8 for all bands, making it impossible for any single band to provide high-quality reference data. To address this, the data processing system selects the three bands with the highest atmospheric transmittance, and uses the transmittance values ​​of each band as weights to perform a weighted average of the observed radiance and atmospheric parameters for each of the three bands. The weighted average data is then used as the equivalent reference band data and input into the subsequent initial temperature calculation step to reduce the impact of single-band noise on the initial temperature estimation.

[0123] S205. Based on the reference band data and the preset initial emissivity constant, the initial temperature is calculated by reversing the thermal infrared radiation transfer equation.

[0124] The initial emissivity constant refers to the emissivity value preset for the reference band during the initial temperature estimation stage. This value is typically a fixed value within the range of 0.95 to 0.98, reflecting the physical characteristic that the emissivity of most natural features is close to 1 within the thermal infrared atmospheric window. The inverted thermal infrared radiative transfer equation refers to the process of solving for T by treating the surface temperature T in the radiative transfer equation as the only unknown quantity, based on known observed radiance and atmospheric parameters. The initial temperature refers to the estimated surface temperature obtained through the above inverted calculation, which serves as the initial temperature value for subsequent iterative optimization.

[0125] Specifically, after acquiring reference band data, the data processing system inverts the thermal infrared radiation transfer equation for the reference band to solve for surface temperature. The data processing system first uses the observed radiance of the reference band... Subtract the atmospheric up-radiation in this band The radiance after atmospheric upward radiation correction is obtained; this result is then divided by the atmospheric transmittance τ of the reference band to obtain the surface emitted radiance; finally, the downward atmospheric radiation is subtracted from the surface emitted radiance. With (1- The product of (where) (Using a preset initial emissivity constant), the surface thermal radiation contribution term is obtained; this surface thermal radiation contribution term is divided by the initial emissivity constant. The Planck radiance value B(T) corresponding to the surface temperature was obtained. The data processing system then processed the Planck function, i.e., Perform the inverse operation, that is ,in and Let λ be the first and second radiation constants of the Planck function, λ be the center wavelength of the reference band, and B be the Planck radiance value calculated in the previous step. Solving for the initial temperature yields the solution. .

[0126] In some embodiments, the initial temperature can be calculated in multiple ways. Optionally, the data processing system performs the above-described inversion calculation using a single reference band, directly using the temperature value obtained from the single band as the initial temperature. This method requires the least computation and is suitable for situations where there are high-quality reference bands with atmospheric transmittance close to 1. Optionally, the data processing system selects multiple bands with high atmospheric transmittance, performs inversion calculations on each band to obtain multiple temperature estimates, and then takes the arithmetic mean of these temperature estimates as the initial temperature. This method reduces the impact of single-band observation noise on the initial temperature estimate through multi-band averaging. The data processing system performs the inversion of the radiative transfer equation and the inverse operation of the Planck function on each band to obtain the temperature estimate vector, and then calculates the arithmetic mean of the vector. It is understandable that the initial temperature can also be estimated based on the maximum brightness temperature, which is not limited here.

[0127] In some embodiments, there may be a significant deviation between the preset initial emissivity constant and the actual emissivity of the reference band (e.g., if the ground feature is a quartz mineral, the emissivity may be as low as 0.7 in some bands), resulting in an initial temperature estimation deviation exceeding 5K. To address this, the data processing system performs a physical plausibility check after calculating the initial temperature, comparing it with the climatic background temperature range of the target area. When the initial temperature exceeds the preset physical temperature range (e.g., 200K to 400K), the data processing system corrects the initial temperature to the boundary value of that range, and gradually corrects this initial deviation through alternating optimization in subsequent iterations.

[0128] In step S205, the data processing system calculates the initial temperature based on the reference band data and the initial emissivity constant inverted thermal infrared radiation transfer equation. The following process further describes how the surface thermal radiation term is constructed in this inverted calculation and the specific execution process of the Planck function inverse operation.

[0129] In some embodiments, the data processing system refines the initial temperature calculation process in step S205. Specifically, the data processing system acquires the radiance, atmospheric up-going radiation, and atmospheric down-going radiation corresponding to the reference band data; combines the radiance, atmospheric up-going radiation, atmospheric down-going radiation, and initial emissivity constant to construct the surface thermal radiation term; and solves the surface thermal radiation term based on the inverse function of the Planck function to obtain the initial temperature.

[0130] Among them, radiance This refers to the total radiance value received at the entrance pupil of the reference band sensor, measured in W / ( (·sr·μm). Atmospheric upward radiation This refers to the radiance component emitted by the atmosphere along the observation path in the reference band towards the sensor. Downward atmospheric radiation. This refers to the hemispherical integrated radiance emitted by the atmosphere towards the Earth's surface in the reference band, the portion of which, after being reflected by the surface, is transmitted through the atmosphere to the sensor. (Earth surface thermal radiation term) It refers to the equivalent Planck radiance value obtained by normalizing the observed radiance to emissivity after subtracting the contributions of atmospheric path radiation and surface reflection from the downward atmospheric radiation. Its physical meaning is the blackbody radiance corresponding to the temperature of the Earth's surface in the reference band.

[0131] Specifically, the data processing system extracts radiance from the reference band data. Atmospheric transmittance Atmospheric upward radiation and atmospheric downward radiation The data processing system constructs the surface thermal radiation term according to the following steps: First, calculate the surface emitted radiance. The first step is to subtract the upward atmospheric radiation from the observed radiance and divide by the atmospheric transmittance to obtain the radiance in the direction of surface emission; the second step is to calculate the contribution of downward atmospheric radiation reflected from the surface. ,in The first step is to set a preset initial emissivity constant; the third step is to calculate the contribution of the Earth's own thermal radiation. The fourth step is to divide the Earth's own surface thermal radiation contribution by the initial emissivity constant to obtain the surface thermal radiation term. The data processing system for Applying Planck's function (i.e. The inverse function of ) is given by: ,in The first radiation constant, The second radiation constant, The initial temperature was obtained using the center wavelength (in μm) of the reference band. .

[0132] In some embodiments, the surface thermal radiation term and the inverse Planck function can be constructed and performed in various ways. Optionally, the data processing system performs subtraction, division, and inverse function operations sequentially according to the above four steps, with intermediate variables stored as double-precision floating-point numbers to ensure numerical accuracy. This method is a direct calculation method and is suitable when all reference band data values ​​are within the normal range. Optionally, the data processing system verifies the consistency of atmospheric parameters before constructing the surface thermal radiation term, checking whether the upward atmospheric radiation, atmospheric transmittance, and atmospheric temperature profiles satisfy physical constraints. Only after the verification is passed can the surface thermal radiation term construction and the inverse Planck function operation be performed. This method adds a quality assurance step for the input data. It is understood that a lookup table method (pre-establishing a radiance-temperature reference table) can also be used to replace the analytical inverse function operation to improve calculation speed; this is not limited here.

[0133] In some embodiments, surface thermal radiation may be present. Negative or near-zero results may be caused by errors in atmospheric parameter estimation or noise in the reference band observations. To address this, the data processing system should... Then perform numerical verification, when When the radiance is less than the preset minimum radiance threshold (e.g., the radiance value corresponding to 150K blackbody radiation), the data processing system will... This threshold is set to ensure that the input value for the inverse operation of the Planck function is within the valid domain, thus avoiding calculation abnormalities caused by non-positive numbers in the logarithmic function.

[0134] S206. Substitute the initial temperature into the thermal infrared radiation transfer equation and perform forward deduction using hyperspectral observation data to calculate the initial emissivity vector.

[0135] Forward deduction refers to the process of calculating the emissivity of each band from the observed radiance by utilizing the relationships between the physical quantities in the thermal infrared radiative transfer equation under known temperature conditions. The initial emissivity vector refers to the N-dimensional emissivity estimate vector obtained through forward deduction, which serves as the initial emissivity value for subsequent iterative optimization.

[0136] Specifically, the data processing system obtains the initial temperature in step S205. After that, The initial emissivity vector is calculated by substituting the values ​​into the thermal infrared radiation transfer equations for each band. For the i-th band (i=1,2,...,N), the data processing system first calculates the temperature using the Planck function. Radiance B at the center wavelength of the i-th band , Then rewrite the radiative transfer equation in terms of emissivity. Explicit expression, ,in To observe the radiance in the i-th band, The atmospheric transmittance of the i-th band is... This represents the upward atmospheric radiation in the i-th band. This represents the downward atmospheric radiation in band i. The data processing system performs the above calculations sequentially on all N bands, combining the resulting N emissivity values ​​into an initial emissivity vector. .

[0137] In some embodiments, the initial emissivity vector can be processed in various ways. Optionally, the data processing system directly imposes physical range constraints on the calculated initial emissivity vector, setting components less than 0.5 to 0.5 and components greater than 1.0 to 1.0, to ensure that the initial emissivity vector is within a reasonable range of natural land cover emissivity. This method is suitable for scenarios where the target area is dominated by natural land cover types. Optionally, the data processing system calculates the Euclidean distance between the initial emissivity vector and the mean vector of the Gaussian prior model. When this distance exceeds a preset threshold, the initial emissivity vector is replaced with the mean vector μ of the Gaussian prior model to ensure that the iteration starting point is within the high probability density region of the prior distribution. This method is suitable for scenarios where the initial temperature estimation deviation is large, causing the initial emissivity vector to deviate significantly from a reasonable range. It is understood that regularization constraints or spectral smoothing can also be used to preprocess the initial emissivity vector, which is not limited here.

[0138] In some embodiments, certain bands may have extremely low atmospheric transmittance (e.g., transmittance less than 0.3), causing the denominator in the inversion calculation of the radiative transfer equation to approach zero, resulting in abnormally large or negative initial emissivity values ​​for these bands. To address this, the data processing system performs a threshold judgment on the atmospheric transmittance of each band. For bands with transmittance below 0.5, the system does not perform the inversion calculation of the radiative transfer equation, but instead sets the initial emissivity value of that band to the component value of the corresponding band in the mean vector of the Gaussian prior model, thus preventing the numerical instability of low-transmittance bands from being propagated to subsequent iterations.

[0139] S207. Use the initial temperature and initial emissivity vector as the initial values ​​for alternating iterative optimization to constrain the optimization space of the Bayesian inversion model.

[0140] Here, the initial values ​​for iteration refer to the temperature and emissivity vector used in the 0th iteration of the alternating iterative optimization process. These initial values ​​determine the starting position of the optimization algorithm in the parameter space. The optimization space refers to the search range of temperature and emissivity during the alternating iterative optimization process. A reasonable initial value for iteration ensures that the search range is constrained to a local region near the optimal solution.

[0141] Specifically, the data processing system obtains the initial temperature in steps S205 and S206 respectively. and initial emissivity vector After that, Assign the initial value to the temperature variable obtained through alternating iterative optimization. The initial value is assigned to the emissivity variable used in the alternating iterative optimization. The data processing system uses... Set the temperature search range around the center. -ΔT, +ΔT], where ΔT is a preset temperature search radius (e.g., 10K), and the optimal temperature is searched only within this range in each subsequent temperature update sub-step. Meanwhile, the data processing system uses... Emissivity optimization is initiated from the initial point. Due to the presence of the Gaussian prior constraint in the Bayesian inversion model, the search range of emissivity is implicitly constrained within the high probability density region of the prior distribution, i.e., the deviation from the mean vector μ does not exceed a certain number of standard deviations. By setting physically reasonable initial values ​​for iteration and combining them with prior constraints, the optimization space of the Bayesian inversion model is effectively constrained, preventing the optimization algorithm from diverging or converging to non-physical solutions in the high-dimensional parameter space.

[0142] Among them, the Gaussian prior model of emissivity is composed of a parameter set of K land cover categories { , When the data processing system constructs a set of parameters (k=1, 2, ..., K), it needs to determine the prior parameters to be used for the current pixel before initiating the alternating iterative optimization. The data processing system calculates the initial emissivity vector for the current pixel. with the mean vector of each category Mahalanobis distance between Select the category with the smallest Mahalanobis distance. As the prior class assignment of the current pixel, the prior mean vector and prior covariance matrix used in subsequent alternating iterative optimization are respectively... and During the iteration process, the data processing system recalculates the current emissivity estimate after each round of alternating iterations (i.e., one execution each for temperature update and emissivity update). The Mahalanobis distance to the mean vectors of each class, if the class corresponding to the minimum Mahalanobis distance changes (i.e., the new nearest class). ≠Currently Used Category The data processing system will switch the prior parameters to the mean vector of the new class. Covariance Matrix The iteration continues from the current temperature and emissivity estimates as a new starting point. The category switch is performed a maximum of a preset number of times (e.g., 3 times). After exceeding this number, the current category is fixed and no longer switched to avoid repeated oscillations at the category boundaries. This mechanism ensures that the prior constraints always best match the ground feature characteristics of the current pixel during the inversion process, thus providing the most effective regularization constraints.

[0143] In some embodiments, the optimization space can be constrained by iterative initial values ​​in various ways. Optionally, the data processing system dynamically adjusts the temperature search interval in the temperature update sub-step of each iteration. When the temperature change between two adjacent iterations is less than 1K, the temperature search radius is reduced to 0.5 times its original value to accelerate the fine-grained search in the convergence process. Optionally, the data processing system uses a wider temperature search interval in the first iteration (e.g., ...). (±20K) In subsequent iterations, the search interval is gradually narrowed according to the rate of decrease of the objective function value, so that the initial stage can escape possible local extrema, and the subsequent stages perform an accurate search near the global optimum. It is understood that a multi-starting-point strategy or simulated annealing method can also be used to expand the exploration range of the optimization space, which is not limited here.

[0144] In some embodiments, different land cover categories may correspond to different Gaussian prior models, and the initial temperature and initial emissivity vectors may not yet be associated with a specific land cover category. In response, the data processing system calculates the initial emissivity vector. The Mahalanobis distance between the pixel and the mean vector of the Gaussian prior model for each category is used to select the category with the smallest Mahalanobis distance as the initial category assignment for the current pixel. The Gaussian prior model parameters (mean vector and covariance matrix) corresponding to this category are then used to enter the subsequent alternating iteration optimization to ensure that the prior constraints match the land cover type of the current pixel.

[0145] S208. Iterate alternately to optimize the temperature and emissivity in the Bayesian inversion model until convergence, and obtain the optimal temperature and optimal emissivity.

[0146] Refer to step S104, which will not be repeated here.

[0147] S209. Calculate the approximate covariance of the Bayesian inversion model under optimal temperature and optimal emissivity conditions to obtain the uncertainty measure.

[0148] Refer to step S105, which will not be repeated here.

[0149] S210. Based on the uncertainty metric, the optimal temperature and optimal emissivity are screened to obtain the target temperature and target emissivity.

[0150] Refer to step S106, which will not be repeated here.

[0151] After completing the screening based on uncertainty metrics, the data processing system has obtained the target temperature and target emissivity for each valid pixel. Because the inversion process in steps S201 to S210 is performed pixel-by-pixel independently, it does not utilize the spatial correlation between adjacent pixels in the distribution of surface temperature and emissivity. In actual thermal infrared hyperspectral remote sensing images, surface temperature and emissivity typically exhibit localized smoothing characteristics in space, meaning that the temperature and emissivity differences between adjacent pixels are usually small under normal circumstances.

[0152] When observation noise causes abnormal deviations in the inversion results of individual pixels, spatial neighborhood information can be used to constrain and correct such anomalies. Steps S211 to S213 describe the process of jointly correcting the target temperature and target emissivity based on spatial neighborhood information.

[0153] S211. Obtain the set of spatial neighborhood pixels corresponding to the target temperature and target emissivity.

[0154] The spatial neighborhood pixel set refers to the set of surrounding pixels selected in the image space coordinate system according to a preset window size, centered on the currently processed pixel. The preset window size is typically (2W+1)×(2W+1), where W is the neighborhood radius, ranging from 1 to 5 pixels. Each pixel in the spatial neighborhood pixel set has completed the processing steps S201 to S210 and has the corresponding target temperature and target emissivity.

[0155] Specifically, after completing the inversion and filtering of all valid pixels in step S210, the data processing system extracts spatial neighborhood information for each valid pixel. Centered on the row and column coordinates (r, c) of the current pixel, the data processing system expands by W pixels in both the row and column directions to obtain the target temperature and target emissivity of all pixels within a (2W+1)×(2W+1) window centered at (r, c). The data processing system performs validity checks on the pixels within the window, retaining only those pixels marked as quality qualified in step S210, excluding invalid and low-confidence pixels, and forming a spatial neighborhood pixel set from the retained valid neighborhood pixels. When the current pixel is located in the image edge region, the data processing system does not process neighborhood positions outside the image range, only taking valid pixels within the image range to form the spatial neighborhood pixel set.

[0156] In some embodiments, the construction parameters of the spatial neighborhood cell set can be determined in various ways. Optionally, the data processing system uses a rectangular window of a fixed size (such as 5×5 or 7×7) as the neighborhood range, applying the same window size to all cells. This method is simple to implement, computationally efficient, and suitable for scenarios with uniform spatial resolution and sparse ground feature boundaries. Optionally, the data processing system dynamically adjusts the window size based on the uncertainty metric of the current cell. Cells with higher uncertainty use a larger window to obtain more neighborhood information for correction, while cells with lower uncertainty use a smaller window to retain more local details. The data processing system maps the posterior standard deviation of temperature to the window radius. This makes the correction strength proportional to the uncertainty of the result. Understandably, spatial neighborhood cell sets can also be constructed using methods such as circular windows or irregular neighborhoods based on ground feature segmentation, which are not limited here.

[0157] In some embodiments, the current pixel may be located at a boundary of abrupt change in land cover type (such as the boundary between a building and bare soil), resulting in the spatial neighborhood pixel set containing pixels of different land cover types. This can cause simple spatial averaging to blur the land cover boundaries. To address this, the data processing system introduces spectral similarity constraints when constructing the spatial neighborhood pixel set. It calculates the spectral angle between the target emissivity vector of each pixel in the neighborhood and the target emissivity vector of the current pixel, and only retains pixels with a spectral angle less than a preset threshold (such as 0.1 radians) in the spatial neighborhood pixel set. This ensures that the neighboring pixels participating in subsequent spatial correction belong to the same or similar land cover types as the current pixel.

[0158] S212. Calculate the spatial similarity weight between the target temperature and the temperatures of each neighboring pixel in the spatial neighborhood pixel set.

[0159] Spatial similarity weight refers to a normalized weight coefficient calculated based on the temperature difference between the current pixel and its neighboring pixels. Neighboring pixels with smaller temperature differences are assigned a higher weight. The calculation method for spatial similarity weight follows the mathematical form of the Gaussian kernel function. ,in The target temperature for the current pixel. For the target temperature of the j-th neighboring pixel, This is the bandwidth parameter of the range kernel function.

[0160] Specifically, after obtaining the spatial neighborhood pixel set in step S211, the data processing system calculates the spatial similarity weight of each neighborhood pixel relative to the current pixel. The data processing system first calculates the target temperature of the current pixel. Target temperature of each neighboring pixel The difference between Then calculate the range weights of each neighboring pixel. ,in The preset bandwidth parameter is used (e.g., 1K to 3K). Simultaneously, the data processing system calculates the spatial distance between the current pixel and its neighboring pixels. Calculate the spatial weight of each neighboring pixel using Euclidean distance (in pixels). ,in The preset spatial bandwidth parameter (e.g., 1 to 3 pixels) is used. The data processing system multiplies the range weight by the spatial weight to obtain the comprehensive weight for each neighboring pixel. Then, the combined weights of all neighboring pixels are normalized so that the sum of the weights equals 1. The normalized weights are the spatial similarity weights.

[0161] The normalization process specifically includes: the data processing system can calculate the comprehensive weight of all valid neighboring pixels in the spatial neighborhood pixel set. Summing yields the total weights. Then divide the combined weight of each neighboring pixel by . Obtain normalized weights In actual calculations, the following may occur: In cases of extremely small degradation (where the temperature difference between all neighboring pixels and the current pixel is significant, causing all Gaussian kernel function values ​​to approach zero), normalized division operations may result in numerical overflow or unstable results. The data processing system employs the following numerical protection strategy to address this: During calculations... Then it is determined whether it is less than the preset minimum weight and threshold. (like ),like This indicates that the temperature difference between the current pixel and all neighboring pixels exceeds the effective response range of the kernel function. In this case, the data processing system does not perform neighborhood-weighted temperature correction, but instead directly outputs the target temperature of the current pixel as the final temperature. The corresponding emissivity is not updated with inverse compensation; the target emissivity is directly output as the final emissivity. The physical meaning of this processing method is: when surrounding pixels do not possess similar temperature characteristics to the current pixel, spatial constraint information is unreliable, and retaining the independent inversion result of the current pixel is more reasonable than introducing irrelevant neighborhood information for correction. Furthermore, when the current pixel itself is also included in the weighted average calculation (i.e., its own weight is...),... =1, because the temperature difference is zero, so the Gaussian kernel function value is 1), the data processing system will use the current pixel's value as the base value. Included together In the calculation, even if the weights of all neighboring pixels approach zero, the weight of the current pixel itself approaches 1 after normalization, and the output is still dominated by the temperature of the current pixel itself.

[0162] In some embodiments, spatial similarity weights can be calculated in multiple ways. Optionally, the data processing system uses only the range weight (i.e., temperature difference weight) as the spatial similarity weight, without introducing a spatial distance weight. This approach is suitable when all neighboring pixels are equidistant (e.g., 8 neighbors in a 3×3 window) or when the influence of spatial distance is negligible. The data processing system calculates and normalizes the Gaussian kernel function value for the temperature difference of each neighboring pixel. Optionally, the data processing system introduces an additional uncertainty weight in addition to the range weight and spatial distance weight, using the reciprocal of the posterior standard deviation of the temperature of each neighboring pixel as the uncertainty weight factor. This allows neighboring pixels with lower uncertainty to contribute more to the weighted average, with the overall weight being... The normalized value is used as the final spatial similarity weight. It is understood that other kernel functions (such as uniform kernel or triangular kernel) can also be used to calculate spatial similarity weights, which are not limited here.

[0163] In some embodiments, there may be a situation where most of the neighboring pixels around the current pixel are marked as invalid due to cloud cover or sensor malfunction, resulting in an insufficient number of effective neighboring pixels (e.g., less than 3). To address this, when the number of effective neighboring pixels is below a preset minimum threshold, the data processing system skips the spatial joint correction step for that pixel and directly outputs the target temperature and target emissivity from step S210 as the final result, thus avoiding unreliable spatial correction when neighboring information is insufficient.

[0164] S213. Based on spatial similarity weights, perform spatial joint correction on the target temperature and target emissivity to obtain the final temperature and final emissivity.

[0165] Spatial joint correction refers to a joint processing procedure that uses spatial similarity weights to perform a weighted average correction on the target temperature of the current pixel, and then updates the emissivity accordingly based on the corrected temperature, ensuring physical consistency between temperature correction and emissivity updates. The final temperature and final emissivity refer to the surface temperature value and emissivity spectral vector output after spatial joint correction, serving as the final products of the entire inversion process.

[0166] Specifically, after obtaining the spatial similarity weights of each neighboring pixel in step S212, the data processing system performs a spatially weighted average of temperature and an inverse compensation update of emissivity. The data processing system calculates a weighted average of the target temperatures of the current pixel and all neighboring pixels in the spatial neighboring pixel set according to normalized spatial similarity weights. ,in This is the normalized weight of the current pixel (the temperature difference between the current pixel and itself is zero, corresponding to a Gaussian kernel function value of 1). The normalized weights for the j-th neighboring pixel are used to obtain the final temperature. The mathematical form of this weighted averaging operation is equivalent to the output of an adaptive bilateral filter, i.e., a locally weighted average constrained by both spatial distance and temperature range differences. This is used to obtain the final temperature. Afterwards, the data processing system will Substitute into the thermal infrared radiative transfer equation to calculate the simulated radiance of each band. And compare the simulated radiance with the observed radiance. Subtraction yields the radiance residual. Because temperature corrections cause changes in Planck radiance, these changes need to be compensated for by adjusting emissivity to maintain the accuracy of the observed data. The data processing system performs inverse compensation updates on the target emissivity based on the radiance residual ΔL(i). The final emissivity vector is obtained. .

[0167] In some embodiments, joint spatial correction can be performed in several ways. Optionally, the data processing system performs a single weighted average correction on temperature and a single inverse compensation update on emissivity to complete the correction process. This method is computationally efficient and suitable for cases where the temperature correction magnitude is small (e.g., less than 1K). Optionally, the data processing system performs multiple iterations (e.g., 2 to 3 times) on temperature correction and emissivity compensation updates. In each iteration, the temperature is first weighted and averaged using the current weights, then the radiance residual is calculated based on the new temperature and the emissivity is updated, and then the weights are recalculated based on the new emissivity until the temperature correction magnitude is less than 0.01K, at which point the iteration terminates. This method can handle cases where the accuracy of a single correction is insufficient when the temperature correction magnitude is large. It is understood that joint correction can also be performed using spatial constraint methods such as total variation regularization or Markov random fields, which are not limited here.

[0168] In some embodiments, after spatial joint correction, the final emissivity values ​​of some bands may exceed the physical range [0, 1]. To address this, the data processing system performs range truncation on the final emissivity vector after completing the inverse compensation update. Components less than 0 are set to the larger of the corresponding value in the Gaussian prior model mean vector and 0, while components greater than 1 are set to 1. Simultaneously, the identification information of the truncated bands is recorded and output to the quality flag bit for subsequent applications to determine the reliability of the pixel inversion result.

[0169] In step S213, the data processing system performs joint spatial correction on the target temperature and target emissivity based on spatial similarity weights. The following process further describes the construction method of the adaptive bilateral filter, the calculation method of the radiance residual, and the specific execution process of the inverse emissivity compensation update in this joint correction.

[0170] In some embodiments, the data processing system refines the spatial joint correction process in step S213. Specifically, the data processing system constructs an adaptive bilateral filter based on spatial similarity weights, performs local weighted averaging on the target temperature to obtain the final temperature, substitutes the final temperature back into the thermal infrared radiation transfer equation, and performs residual calculation based on the hyperspectral observation data of the corresponding pixel to obtain the radiance residual. Based on the radiance residual, the target emissivity is updated by inverse compensation to obtain the final emissivity.

[0171] The adaptive bilateral filter is a nonlinear spatial filter that uses spatial distance and temperature range differences as dual kernel functions and spatial similarity weights calculated in step S212 as filtering coefficients. This filter maintains the clarity of ground feature boundaries while locally smoothing the temperature. Local weighted averaging refers to the operation of weighted summation of the target temperature of the current pixel and its spatial neighbors using normalized spatial similarity weights. Radiance residual refers to the band-by-band difference vector between the simulated radiance obtained by substituting the final temperature and target emissivity into the thermal infrared radiative transfer equation and the actual observed radiance. Reverse compensation update refers to the parameter adjustment operation that adjusts the emissivity value in reverse based on the radiance residual, so that the updated emissivity and the final temperature, when substituted into the radiative transfer equation, can better fit the observed data.

[0172] Specifically, the data processing system first uses the normalized spatial similarity weights obtained in step S212 { Construct an adaptive bilateral filter. The filter outputs the temperature of the current pixel (r, c) as follows: The summation iterates through all valid neighboring pixels in the set of spatial neighboring pixels, excluding the current pixel. Let be the normalized spatial similarity weight of the j-th neighboring pixel. Let be the target temperature of the j-th neighboring pixel. The adaptability of this filter is reflected in the weights. The filter varies with temperature differences and spatial distance, resulting in a strong smoothing effect in uniform regions with flat temperature distribution, while at the boundaries of features with temperature gradients, the filter weakens the smoothing effect to maintain boundary sharpness.

[0173] The data processing system obtains the final temperature. Then, the radiance residuals are calculated for each of the N bands. For the i-th band, the data processing system uses the Planck function to calculate B( , Then, the radiative transfer equation is used to calculate the target emissivity. and final temperature Simulated radiance as parameters The residual radiance is The residual reflects the effect of temperature change from Revised to The resulting bias in the fitting of the observed data.

[0174] The data processing system performs inverse emissivity compensation updates based on the radiance residual. For the i-th band, the emissivity update amount is... This expression is obtained by inverting the partial derivative of the radiative transfer equation with respect to emissivity. Final emissivity After performing this compensation update on all N bands, the data processing system obtains the final emissivity vector. .

[0175] In some embodiments, adaptive bilateral filtering and inverse compensation updates can be performed in various ways. Optionally, the data processing system performs bilateral filtering and inverse compensation sequentially on all valid pixels of the full-frame image, with each pixel using a fixed spatial bandwidth parameter. Sum of bandwidth parameters After full-image processing, the final temperature map and final emissivity map are output. This method is suitable for scenes where the distribution of ground cover types in the image is relatively uniform. Optionally, the data processing system dynamically adjusts the value domain bandwidth parameter based on the uncertainty metric of each pixel. Pixels with greater uncertainty use larger [sizes / sizes]. (That is, more tolerant of temperature differences and stronger smoothing) to make better use of neighborhood information to correct unreliable inversion results, and use smaller pixels for pixels with less uncertainty. (That is, more sensitive to temperature differences and more protective) in order to retain its own more reliable inversion results. ,in As the reference value range bandwidth, The posterior standard deviation of the current pixel temperature. This represents the mean of the posterior standard deviation of temperature in the full-frame image. It is understood that other spatial filtering methods, such as guided filtering or nonlocal mean filtering, can be used to replace the bilateral filter; this is not limited here.

[0176] In some embodiments, there may be a denominator term in the reverse compensation update calculation. The near-zero value occurs when atmospheric transmittance is extremely low in certain wavelength bands or when the surface temperature is nearly equal to the equivalent temperature of downward atmospheric radiation. To address this, the data processing system performs a numerical check on the denominator before calculating the reverse compensation update. If the absolute value of the denominator is less than a preset threshold... At that time, the data processing system skips the reverse compensation update for that band and directly sets the final emissivity of that band to the target emissivity value. This is to avoid numerical overflow during division operations.

[0177] In this embodiment, by employing a multivariate Gaussian prior model constructed from the mean vector and covariance matrix of the emissivity sample library, integrating this Gaussian prior model as a regularization constraint into the Bayesian inversion framework and jointly constructing the objective function with the thermal infrared radiation transfer equation, and using alternating iterative optimization to solve for the optimal temperature and optimal emissivity and utilizing the Hessian matrix inverse approximation to achieve endogenous quantification of uncertainty, the statistical distribution characteristics of the emissivity spectrum (including the mean of each band and the covariance relationship between bands) are effectively integrated into the inversion optimization process to provide probabilistic constraints on the ill-conditioned equation set. This effectively solves the problem in the prior constraints based on simple empirical regression relationships in the prior art that cannot fully utilize high-spectral dimensional information and have large deviations in inversion results under complex noise environments. It achieves the technical goal of synchronously outputting high-precision temperature and emissivity inversion results and their pixel-level uncertainty quantification indicators under a unified Bayesian probabilistic framework.

[0178] The data processing system in this application embodiment is described below from a hardware processing perspective. Please refer to [link / reference]. Figure 3 This is a schematic diagram of the physical device structure of a data processing system in an embodiment of this application.

[0179] It should be noted that, Figure 3 The structure of the data processing system shown is merely an example and should not impose any limitations on the functionality and scope of use of the embodiments of this application.

[0180] like Figure 3 As shown, the data processing system includes a CPU 301, which can perform various appropriate actions and processes according to a program stored in ROM 302 or a program loaded from storage section 308 into RAM 303, such as executing the methods described in the above embodiments. RAM 303 also stores various programs and data required for system operation. The CPU 301, ROM 302, and RAM 303 are interconnected via bus 304. I / O interface 305 is also connected to bus 304.

[0181] The following components are connected to I / O interface 305: input section 306 including audio input devices, push-button switches, etc.; output section 307 including liquid crystal display (LCD) and audio output devices, indicator lights, etc.; storage section 308 including hard disks, etc.; and communication section 309 including network interface cards such as LAN (Local Area Network) cards, modems, etc. Communication section 309 performs communication processing via a network such as the Internet. Drive 310 is also connected to I / O interface 305 as needed. Removable media 311, such as disks, optical disks, magneto-optical disks, semiconductor memories, etc., are installed on drive 310 as needed so that computer programs read from them can be installed into storage section 308 as needed.

[0182] Specifically, according to embodiments of this application, the processes described above with reference to the flowcharts can be implemented as computer software programs. For example, embodiments of this application include a computer program product comprising a computer program carried on a computer-readable medium, the computer program including a computer program for performing the methods shown in the flowcharts. In such embodiments, the computer program can be downloaded and installed from a network via communication section 309, and / or installed from removable medium 311. When the computer program is executed by CPU 301, it performs the various functions defined in this application.

[0183] The flowcharts and block diagrams in the accompanying drawings illustrate the architecture, functionality, and operation of possible implementations of systems, methods, and computer program products according to various embodiments of this application. Each block in a flowchart or block diagram may represent a module, segment, or portion of code, which contains one or more executable instructions for implementing a specified logical function. It should also be noted that in some alternative implementations, the functions indicated in the blocks may occur in a different order than those shown in the drawings.

[0184] Specifically, the data processing system of this embodiment includes a processor and a memory. The memory stores a computer program. When the computer program is executed by the processor, it implements the thermal infrared hyperspectral temperature and emissivity inversion method based on the Gaussian prior model provided in the above embodiment.

[0185] In another aspect, this application also provides a computer-readable storage medium, which may be included in the data processing system described in the above embodiments; or it may exist independently and not assembled into the data processing system. The storage medium carries one or more computer programs that, when executed by a processor of the data processing system, cause the data processing system to implement the thermal infrared hyperspectral temperature and emissivity inversion method based on a Gaussian prior model provided in the above embodiments.

[0186] The above-described embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit it. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of this application.

[0187] As used in the above embodiments, depending on the context, the term "when..." can be interpreted as meaning "if...", "after...", "in response to determining...", or "in response to detecting...". Similarly, depending on the context, the phrase "when determining..." or "if (the stated condition or event) is interpreted as meaning "if determining...", "in response to determining...", "when (the stated condition or event) is detected", or "in response to detecting (the stated condition or event)".

Claims

1. A method for inverting temperature and emissivity in thermal infrared hyperspectral imaging based on a Gaussian prior model, characterized in that, Applied to a data processing system, the method includes: Acquire hyperspectral observation data and a pre-built emissivity sample library; Calculate the mean vector and covariance matrix of the emissivity samples in the emissivity sample library, and construct a Gaussian prior model of emissivity; Based on the thermal infrared radiative transfer equation, the observation noise distribution, and the Gaussian prior model, a Bayesian inversion model for the hyperspectral observation data is constructed. The temperature and emissivity in the Bayesian inversion model are iteratively optimized until convergence, yielding the optimal temperature and optimal emissivity. The approximate covariance of the Bayesian inversion model under the optimal temperature and optimal emissivity conditions is calculated to obtain an uncertainty measure; The optimal temperature and optimal emissivity are screened based on the uncertainty metric to obtain the target temperature and target emissivity.

2. The method according to claim 1, characterized in that, The steps of calculating the mean vector and covariance matrix of the emissivity samples in the emissivity sample library and constructing a Gaussian prior model of emissivity specifically include: The emissivity samples in the emissivity sample library are divided according to the clustering algorithm to obtain multiple subsets of land cover categories; Calculate the mean vector and initial covariance matrix of the emissivity samples within the subset of the land cover categories; The initial covariance matrix is ​​regularized by combining a preset positive number and an identity matrix to obtain a regularized covariance matrix; The mean vector and the regularized covariance matrix are combined to determine the Gaussian prior model of the emissivity.

3. The method according to claim 1, characterized in that, The steps for constructing a Bayesian inversion model of the hyperspectral observation data based on the thermal infrared radiative transfer equation, the observation noise distribution, and the Gaussian prior model specifically include: The observation fitting term is calculated based on the hyperspectral observation data, the forward simulation data corresponding to the thermal infrared radiative transfer equation, and the covariance matrix of the observation noise distribution. The prior constraint terms are calculated based on the emissivity, the mean vector in the Gaussian prior model, and the covariance matrix in the Gaussian prior model. The objective function is obtained by fusing the observed fitting term and the prior constraint term. The probability estimation model containing the objective function is determined as the Bayesian inversion model of the hyperspectral observation data.

4. The method according to claim 1, characterized in that, Before the step of iteratively optimizing the temperature and emissivity in the Bayesian inversion model until convergence to obtain the optimal temperature and optimal emissivity, the method further includes: Extract reference band data that is least affected by atmospheric absorption from the hyperspectral observation data; Based on the reference band data and the preset initial emissivity constant, the initial temperature is calculated by reversing the thermal infrared radiation transfer equation. Substituting the initial temperature into the thermal infrared radiation transfer equation and combining it with the hyperspectral observation data, a forward deduction is performed to calculate the initial emissivity vector; The initial temperature and the initial emissivity vector are used as the initial values ​​for the alternating iterative optimization to constrain the optimization space of the Bayesian inversion model.

5. The method according to claim 4, characterized in that, The step of calculating the initial temperature by reversing the thermal infrared radiation transfer equation based on the reference band data and a preset initial emissivity constant specifically includes: Obtain the radiance, atmospheric up-going radiation, and atmospheric down-going radiation corresponding to the reference band data; By combining the radiance, the upward atmospheric radiation, the downward atmospheric radiation, and the initial emissivity constant, a surface thermal radiation term is constructed; The initial temperature is obtained by solving the surface thermal radiation term using the inverse of the Planck function.

6. The method according to claim 1, characterized in that, After the step of filtering the optimal temperature and optimal emissivity based on the uncertainty metric to obtain the target temperature and target emissivity, the method further includes: Obtain the spatial neighborhood pixel set corresponding to the target temperature and the target emissivity; Calculate the spatial similarity weights between the target temperature and the temperatures of each neighboring pixel in the spatial neighborhood pixel set; Based on the spatial similarity weights, the target temperature and the target emissivity are jointly corrected spatially to obtain the final temperature and the final emissivity.

7. The method according to claim 6, characterized in that, The step of performing spatial joint correction on the target temperature and the target emissivity based on the spatial similarity weight to obtain the final temperature and the final emissivity specifically includes: An adaptive bilateral filter is constructed based on the spatial similarity weights, and the target temperature is locally weighted and averaged to obtain the final temperature. Substitute the final temperature back into the thermal infrared radiation transfer equation, and combine it with the hyperspectral observation data of the corresponding pixel to calculate the residual, and obtain the radiance residual. The target emissivity is updated by inverse compensation based on the radiance residual to obtain the final emissivity.

8. A data processing system, characterized in that, The data processing system includes: one or more processors and a memory; the memory is coupled to the one or more processors, the memory is used to store computer program code, the computer program code including computer instructions, and the one or more processors call the computer instructions to cause the data processing system to perform the method as described in any one of claims 1-7.

9. A computer-readable storage medium comprising instructions, characterized in that, When the instructions are executed on the data processing system, the data processing system performs the method as described in any one of claims 1-7.

10. A computer program product, characterized in that, When the computer program product is run on a data processing system, the data processing system performs the method as described in any one of claims 1-7.

Images