Split-window remote sensing inversion method and system for land surface temperature based on virtual channel

Abstract

The present invention belongs to the field of surface temperature inversion technology, and relates to a surface temperature split window remote sensing inversion method and system based on a virtual channel. The method includes constructing a simulation data set, constructing a virtual channel, determining a split window coefficient, and determining a target surface temperature. The present invention eliminates the influence of the surface temperature angle effect caused by the different projection areas of vegetation coverage of ground mixed pixels under different observation zenith angles by normalizing the angle of the brightness temperature on the satellite, so that the observation data under different observation zenith angles are comparable; through the PSO-BP neural network model, the time correction and spectral conversion of the brightness temperature on the satellite are realized to realize the construction of the virtual channel of the 11th channel of Landsat-8 TIRS; the present invention reselects the split window coefficient of the appropriate sub-interval through the surface temperature approximation, the atmospheric water vapor content and the average surface emissivity, and can estimate a more accurate surface temperature.

Description

Technical Field

[0001] The present invention belongs to the technical field of surface temperature inversion, and in particular, relates to a surface temperature split-window remote sensing inversion method and system based on a virtual channel. Background Art

[0002] Land surface temperature is a key driver of the interaction between the land and atmosphere system, and is also one of the most important physical parameters in the water and energy balance from regional to global scales. Since the 1970s, scholars have developed and used various land surface temperature inversion algorithms to produce land surface temperature products. These methods mainly include: single channel algorithm, split window algorithm and temperature emissivity separation algorithm. The thermal infrared sensor carried by Landsat-8 (Landsat 8 Satellite) satellite has a total of 2 thermal infrared channels with a spatial resolution of 100m. As an important part of the Landsat series of satellites, it provides important support for time-continuous 100-meter-level land surface temperature data. Due to the influence of stray light outside the field of view, the calibration accuracy of the Landsat-8 thermal infrared sensor data did not meet the design target, and the uncertainty of the 11th channel was much greater than that of the 10th channel, which made the Landsat-8 thermal infrared sensor have only one available thermal infrared channel in the atmospheric window, and it was impossible to use the existing split window algorithm to invert the land surface temperature, which caused trouble in obtaining time-continuous 100-meter-level land surface temperature data.

[0003] The prior art discloses a method for generating a split window coefficient corresponding to a satellite according to a spectral response function corresponding to the satellite's thermal infrared dual-channel data and a set simulation scheme, which is used for the inversion of surface temperature, and specifically discloses the following technical features: updating the tp5 file according to the spectral response function of the satellite's thermal infrared dual-channel data and the simulation settings of different atmospheric transmission models; generating multiple groups of dual-channel brightness temperature data according to the updated tp5 file and the atmospheric radiation transmission model; performing group regression according to the multiple groups of dual-channel brightness temperature data, and establishing a multi-dimensional regression coefficient lookup table; generating the split window coefficient corresponding to the satellite according to the multi-dimensional regression coefficient lookup table. This method requires that the satellite sensor is provided with at least two thermal infrared earth observation channels in the atmospheric window, and its applicability is limited when the data of a certain earth observation channel cannot be used. Summary of the invention

[0004] In order to solve the above technical problems, the present invention provides a surface temperature split window remote sensing inversion method and system based on virtual channels.

[0005] In a first aspect, the present invention provides a split-window remote sensing inversion method for land surface temperature based on a virtual channel, including constructing a simulation data set, constructing a virtual channel, determining a split-window coefficient, and determining a target land surface temperature;

[0006] Construct a simulated data set, including:

[0007] The atmospheric profile data under clear sky conditions were selected from the TIGR database, and the atmospheric profile data were input into the radiation transfer model to simulate the atmospheric parameters;

[0008] Convolve the ground object spectrum data with the sensor spectral response function to simulate the channel effective surface emissivity;

[0009] The atmospheric parameters and surface emissivity were input into the thermal infrared radiation transfer model to obtain the Terra MODIS on-board brightness temperature simulation data and the Landsat-8 TIRS on-board brightness temperature simulation data as simulation data sets;

[0010] Build a virtual channel, including:

[0011] Normalize the brightness temperature angle on the satellite based on MODIS products, LAI products, vegetation aggregation index and MODIS L1B data;

[0012] A PSO-BP neural network model was constructed for time correction and spectral conversion. The brightness temperature simulation data of the 31st and 32nd channels of the Terra MODIS on-board brightness temperature simulation data at several observation zenith angles and the brightness temperature simulation data of the 10th channel of Landsat-8 TIRS at an observation zenith angle of 0° were used as input, and the brightness temperature simulation data of the 11th channel of Landsat-8 TIRS at the corresponding observation zenith angle was used as output. The PSO-BP neural network model was trained, and the virtual channel of the 11th channel of Landsat-8 TIRS was obtained based on the trained PSO-BP neural network model.

[0013] Determine the split window coefficient, including:

[0014] Based on the brightness temperature simulation data on Landsat-8 TIRS, the split window coefficient is determined according to the split window algorithm and the grouping rate of average surface emissivity, surface temperature and atmospheric water vapor content;

[0015] An approximate value of the surface temperature is determined by a split-window algorithm;

[0016] Determine target land surface temperature, including:

[0017] The splitting window coefficient is reselected according to the average surface emissivity, the approximate value of the surface temperature and the atmospheric water vapor content to determine the target surface temperature.

[0018] In a second aspect, the present invention provides a surface temperature split window remote sensing inversion system based on a virtual channel, comprising a simulation data set construction unit, a virtual channel construction unit, a split window coefficient determination unit and a target surface temperature determination unit;

[0019] Simulated dataset building blocks for:

[0020] The atmospheric profile data under clear sky conditions were selected from the TIGR database, and the atmospheric profile data were input into the radiation transfer model to simulate the atmospheric parameters;

[0021] Convolve the ground object spectrum data with the sensor spectral response function to simulate the channel effective surface emissivity;

[0022] The atmospheric parameters and surface emissivity were input into the thermal infrared radiation transfer model to obtain the Terra MODIS on-board brightness temperature simulation data and the Landsat-8 TIRS on-board brightness temperature simulation data as simulation data sets;

[0023] Virtual channel building blocks for:

[0024] Normalize the brightness temperature angle on the satellite based on MODIS products, LAI products, vegetation aggregation index and MODIS L1B data;

[0025] A PSO-BP neural network model was constructed for time correction and spectral conversion. The brightness temperature simulation data of the 31st and 32nd channels of the Terra MODIS on-board brightness temperature simulation data at several observation zenith angles and the brightness temperature simulation data of the 10th channel of Landsat-8 TIRS at an observation zenith angle of 0° were used as input, and the brightness temperature simulation data of the 11th channel of Landsat-8 TIRS at the corresponding observation zenith angle was used as output. The PSO-BP neural network model was trained, and the virtual channel of the 11th channel of Landsat-8 TIRS was obtained based on the trained PSO-BP neural network model.

[0026] The split window coefficient determination unit is used to:

[0027] Based on the brightness temperature simulation data on Landsat-8 TIRS, the split window coefficient is determined according to the split window algorithm and the grouping rate of average surface emissivity, surface temperature and atmospheric water vapor content;

[0028] An approximate value of the surface temperature is determined by a split-window algorithm;

[0029] The target surface temperature determination unit is used to:

[0030] The splitting window coefficient is reselected according to the average surface emissivity, the approximate value of the surface temperature and the atmospheric water vapor content to determine the target surface temperature.

[0031] Based on the above technical solution, the present invention can also be improved as follows.

[0032] Furthermore, the atmospheric parameters include channel transmittance, atmospheric upwelling radiation and atmospheric downwelling radiation.

[0033] Furthermore, when training the PSO-BP neural network model, a set proportion of the simulated data set is used as the training set, and the rest is used as the test set; the input training set parameters are normalized to between 0 and 1; the parameters of the PSO-BP neural network model are initialized, and the PSO-BP neural network model trains the training set; the number of hidden layers of the PSO-BP neural network model is set, each hidden layer contains a number of hidden layer nodes, the learning rate, the number of training times and the target minimum error are set, and the quantized conjugate gradient method is used for optimization; the PSO algorithm is used to iteratively search for the best, the number of particle populations is 5, and the number of iterations, cognitive learning factor, social learning factor, inertia factor, maximum speed, minimum speed, maximum boundary and minimum boundary are set.

[0034] Furthermore, for mixed pixels containing both soil and vegetation components, the directional effect of surface temperature is corrected according to vegetation coverage, including:

[0035] Assume that the vegetation coverage is , the brightness temperature of the sensor's i-th channel is , Indicates the brightness temperature contributed by the vegetation component of the sensor's i-th channel, Represents the surface emissivity of the soil part in the mixed pixel, represents the brightness temperature contributed by the soil component of the sensor's channel i, represents the atmospheric transmittance, represents the surface emissivity of the vegetation component in the mixed pixel, represents the part of the atmospheric downward radiation in the i-th channel reflected by the soil components, represents the part of the atmospheric downward radiation in the i-th channel reflected by the vegetation components, It represents the contribution of vegetation when the atmosphere of the i-th channel is divided into vegetation and soil according to vegetation coverage. It represents the soil contribution when the atmosphere of the i-th channel is divided into vegetation and soil according to vegetation coverage, then:

[0036]

[0037]

[0038]

[0039] ;

[0040] set up is the base of natural logarithms, is the vegetation aggregation index, is the leaf area index, is the leaf area orientation function; the leaf area orientation function is the average projection of the unit leaf area in the viewing direction, which is used to describe the effect of leaf inclination on radiation attenuation, then:

[0041] ;

[0042] Under the condition that the vegetation aggregation index and leaf area index are known, the value of vegetation coverage depends on the observation zenith angle. Assuming that the pixels in the 3×3 spatial window adjacent to the target pixel have the same component temperature, including vegetation temperature and soil temperature, and the atmospheric environment does not change, for the nine pixels in the spatial window, the equation group is constructed:

[0043] ;

[0044] According to the vegetation coverage of each pixel at the current observation zenith angle, the brightness temperature contributed by each component is determined, and then the brightness temperature is normalized to the brightness temperature at the observation zenith angle of 0° by observing the vegetation coverage at the observation zenith angle of 0°.

[0045] Furthermore, according to the leaf area index, vegetation aggregation index data and the observation zenith angle, the brightness temperatures of the 31st and 32nd channels of Terra MODIS were angle-normalized to obtain the brightness temperature when the observation zenith angle was 0°.

[0046] Furthermore, the surface emissivity of the 10th and 11th channels of Landsat-8 TIRS was calculated using the FR97 algorithm, and the atmospheric water vapor content was estimated using the atmospheric relative humidity provided by the ERA5 reanalysis data.

[0047] Furthermore, the brightness temperature simulation data on Landsat-8 are divided into two groups according to the average emissivity. The simulation data below the first set threshold are the first emissivity group, and the simulation data above the second set threshold are the second emissivity group. The first set threshold is greater than the second set threshold, and the second emissivity group and the second emissivity group containing overlapping areas are obtained. The surface temperature and atmospheric water vapor content are grouped.

[0048] The beneficial effects of the present invention are:

[0049] (1) The present invention eliminates the influence of the surface temperature angle effect caused by the different projection areas of vegetation coverage at different observation zenith angles on ground mixed pixels by normalizing the angle of brightness temperature on the satellite, making the observation data at different observation zenith angles comparable;

[0050] (2) The present invention uses the PSO-BP neural network model to achieve the time correction and spectral conversion of the on-board brightness temperature to realize the construction of the virtual channel of the 11th channel of Landsat-8 TIRS;

[0051] (3) The present invention can estimate the surface temperature more accurately by reselecting the split window coefficient of the appropriate sub-interval through the approximate value of the surface temperature, the atmospheric water vapor content and the average surface emissivity. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] Figure 1 A schematic diagram of a surface temperature split window remote sensing inversion method based on a virtual channel provided in Example 1 of the present invention;

[0053] Figure 2 It is a flow chart of the surface temperature split window remote sensing inversion method based on virtual channel of the present invention;

[0054] Figure 3 It is the relationship curve between the true value and the predicted value; Figure 3 (a) is the relationship curve between the true value and the predicted value of the training set. Figure 3 (b) is a curve chart showing the relationship between the true value and the predicted value of the test set;

[0055] Figure 4 The simulation diagram shows the root mean square error between the inverted surface temperature and the actual surface temperature in different atmospheric water vapor content sub-intervals. Figure 4 (a) is a simulation diagram of the root mean square error corresponding to the Terra MODIS observation zenith angle of 0°. Figure 4 (b) is the simulation diagram of the root mean square error corresponding to the Terra MODIS observation zenith angle of 0°;

[0056] Figure 5 A simulation diagram showing the comparison between the inverted surface temperature and the ground station temperature; Figure 5 (a) is the simulation diagram of the surface temperature at ground station 1. Figure 5 (b) is a simulation diagram of the surface temperature at ground station 1;

[0057] Figure 6 This is a schematic diagram of a surface temperature split-window remote sensing inversion system based on a virtual channel provided in Example 2 of the present invention. DETAILED DESCRIPTION

[0058] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Generally, the components of the embodiments of the present invention described and shown in the drawings here can be arranged and designed in various different configurations.

[0059] Example 1

[0060] As an example, Figure 1 As shown, in order to solve the above technical problems, this embodiment provides a split window remote sensing inversion method for land surface temperature based on a virtual channel, including constructing a simulation data set, constructing a virtual channel, determining a split window coefficient and determining a target land surface temperature;

[0061] Construct a simulated data set, including:

[0062] The atmospheric profile data under clear sky conditions were selected from the TIGR database, and the atmospheric profile data were input into the radiation transfer model to simulate the atmospheric parameters;

[0063] Convolve the ground object spectrum data with the sensor spectral response function to simulate the channel effective surface emissivity;

[0064] The atmospheric parameters and surface emissivity were input into the thermal infrared radiation transfer model to obtain the Terra MODIS (Terra Moderate - Resolution Imaging Spectroradiometer) on-board brightness temperature simulation data and the Landsat-8 TIRS on-board brightness temperature simulation data as simulation data sets;

[0065] Build a virtual channel, including:

[0066] Normalize the brightness temperature angle on board according to MODIS products, LAI (Leaf Area Index) products, vegetation aggregation index and MODISL1B data (Moderate - Resolution Imaging Spectroradiometer Level 1B data);

[0067] A PSO-BP (Particle Swarm Optimization - Back Propagation) neural network model was constructed for time correction and spectral conversion. The on-board brightness temperature simulation data of the 31st and 32nd channels of the Terra MODIS on-board brightness temperature simulation data at several observation zenith angles and the on-board brightness temperature simulation data of the 10th channel of the Landsat-8 TIRS (Thermal Infrared Sensor) at an observation zenith angle of 0° were used as input, and the on-board brightness temperature simulation data of the 11th channel of the Landsat-8 TIRS at the corresponding observation zenith angle was used as output. The PSO-BP neural network model was trained, and the virtual channel of the 11th channel of the Landsat-8 TIRS was obtained based on the trained PSO-BP neural network model.

[0068] Determine the split window coefficient, including:

[0069] Based on the brightness temperature simulation data on Landsat-8 TIRS, the split window coefficient is determined according to the split window algorithm and the grouping rate of average surface emissivity, surface temperature and atmospheric water vapor content;

[0070] An approximate value of the surface temperature is determined by a split-window algorithm;

[0071] Determine target land surface temperature, including:

[0072] The splitting window coefficient is reselected according to the average surface emissivity, the approximate value of the surface temperature and the atmospheric water vapor content to determine the target surface temperature.

[0073] The flowchart of this method is shown in the attached Figure 2 The present invention proposes a method for constructing a simulated data set, namely, a method for constructing a virtual channel of the 11th channel of Landsat-8TIRS and a method for constructing a split window algorithm.

[0074] The simulated data set is used for the construction of virtual channels and split window algorithms. Through the construction of virtual channels, the split window remote sensing inversion of the surface temperature of Landsat-8 data is realized.

[0075] In actual application, the simulation data set contains the on-board brightness temperature simulation data of Terra MODIS channel 31 and channel 32 and Landsat-8 TIRS channel 10 and channel 11.

[0076] The on-board brightness temperature simulation data is simulated based on the thermal infrared radiation transfer equation. The surface emissivity is obtained by convolving the ECOSTRESS (Ecosystem Spaceborne Thermal Radiometer Experiment on Space Station) ground object spectrum library with the sensor spectral response function. The atmospheric parameters are obtained by the radiation transfer model MODTRAN5.2. For example, by inputting 978 atmospheric profiles under clear sky conditions screened from the TIGR database into the model, the atmospheric transmittance, atmospheric uplink radiation and atmospheric downlink radiation at each wavelength are obtained.

[0077] After obtaining the atmospheric parameters, the channel transmittance, atmospheric upwelling radiation, and atmospheric downwelling radiation of the 31st and 32nd channels of Terra MODIS and the 10th and 11th channels of Landsat-8 TIRS are obtained by convolution with the sensor spectral response function. In the actual application process, considering the influence of the observation angle of Terra MODIS, the atmospheric parameters are simulated at 6 different observation zenith angles (0°, 33.56°, 44.42°, 51.32°, 56.25°, and 60°). In order to make the simulated data more consistent with the actual situation, the possible changes in the surface temperature depend on the boundary layer temperature of the input atmospheric profile. When the boundary layer temperature T0 of the atmospheric profile is lower than 290 K, the surface temperature is set to 3 values: T0–5 K, T0, and T0+5 K. When the boundary layer temperature T0 of the atmospheric profile is higher than 290 K, the surface temperature is set to 5 values: T0–15 K, T0–10 K, T0–5 K, T0, and T0+5 K. Taking into account the different transit times of Terra and Landsat-8, when simulating the MODIS on-board brightness temperature to set the surface temperature, the surface temperature set according to the Landsat-8 temperature is perturbed. Each Landsat-8 surface temperature is set with 7 MODIS temperatures, from T-2K to T+1K, with a step size of 0.5K.

[0078] Optionally, the atmospheric parameters include channel transmittance, atmospheric up-radiation and atmospheric down-radiation.

[0079] Optionally, when training the PSO-BP neural network model, a set proportion of the simulated data set is used as the training set, and the rest is used as the test set; the input training set parameters are normalized to between 0 and 1; the parameters of the PSO-BP neural network model are initialized, and the PSO-BP neural network model trains the training set; the number of hidden layers of the PSO-BP neural network model is set, each hidden layer contains a number of hidden layer nodes, the learning rate, the number of training times and the target minimum error are set, and the quantized conjugate gradient method is used for optimization; the PSO algorithm is used for iterative optimization, the number of particle populations is 5, and the number of iterations, cognitive learning factor, social learning factor, inertia factor, maximum speed, minimum speed, maximum boundary and minimum boundary are set.

[0080] In the actual application process, PSO-BP neural network is used to realize time correction and spectral conversion, and a virtual channel of the 11th channel of Landsat-8 TIRS is constructed. The training data set is the simulated brightness temperature data of MODIS and Landsat-8. The input parameters are: the brightness temperature of Terra MODIS channels 31 and 32 at 6 observation zenith angles, the secant value of MODIS observation zenith angle, and the brightness temperature of the 10th channel of Landsat-8 TIRS when the observation zenith angle is 0°; the output parameter is the brightness temperature of the 11th channel of Landsat-8 TIRS at the corresponding observation zenith angle. 2 / 3 of the simulated data set is used as the training set, and the remaining 1 / 3 is used as the test set. In order to facilitate the comparison of parameters of different dimensions, the input training set parameters are normalized to between 0 and 1. Then the parameters of the BP neural network are initialized, and the PSO-BP neural network model is trained using the training set. The hidden layer is set to 1 layer, with a total of 7 hidden layer nodes, the learning rate is set to 0.01, the number of training times is 1000, and the target minimum error is 1×10 -6 , the quantized conjugate gradient method was used for optimization; then the PSO algorithm was used for iterative optimization, the particle population size was 5, the number of iterations was 30, the cognitive learning factor and the social learning factor were both set to 2, the inertia factor was set to 1, the maximum speed was set to 1, the minimum speed was set to -1, the maximum boundary was set to 1, and the minimum boundary was set to -1.

[0081] The construction of the virtual channel includes the angle normalization, time correction and spectral conversion of the brightness temperature on the satellite. For mixed pixels containing soil components and vegetation components, the proportion of the projected area of ​​each component in the pixel will change with the change of the observation zenith angle, resulting in the total outgoing radiation brightness of the pixel changing with the change of the observation angle. The brightness temperature at the sensor is decomposed into the soil contribution part and the vegetation contribution part. The directionality of the brightness temperature on the satellite is mainly determined by the projection ratio, that is, the vegetation coverage.

[0082] Optionally, for mixed pixels containing both soil and vegetation components, the directional effect of surface temperature is corrected according to vegetation cover, including:

[0083] Assume that the vegetation coverage is , the brightness temperature of the sensor's i-th channel is , Indicates the brightness temperature contributed by the vegetation component of the sensor's i-th channel, Represents the surface emissivity of the soil part in the mixed pixel, represents the brightness temperature contributed by the soil component of the sensor's channel i, represents the atmospheric transmittance, represents the surface emissivity of the vegetation component in the mixed pixel, represents the part of the atmospheric downward radiation in the i-th channel reflected by the soil components, represents the part of the atmospheric downward radiation in the i-th channel reflected by the vegetation components, It represents the contribution of vegetation when the atmosphere of the i-th channel is divided into vegetation and soil according to vegetation coverage. It represents the soil contribution when the atmosphere of the i-th channel is divided into vegetation and soil according to vegetation coverage, then:

[0084]

[0085]

[0086]

[0087] ;

[0088] set up is the base of natural logarithms, is the vegetation aggregation index, is the leaf area index, is the leaf area orientation function; the leaf area orientation function is the average projection of the unit leaf area in the viewing direction, which is used to describe the effect of leaf inclination on radiation attenuation, then:

[0089] ;

[0090] Under the condition that the vegetation aggregation index and leaf area index are known, the value of vegetation coverage depends on the observation zenith angle. Assuming that the pixels in the 3×3 spatial window adjacent to the target pixel have the same component temperature, including vegetation temperature and soil temperature, and the atmospheric environment does not change, for the nine pixels in the spatial window, the equation group is constructed:

[0091] ;

[0092] According to the vegetation coverage of each pixel at the current observation zenith angle, the brightness temperature contributed by each component is determined, and then the brightness temperature is normalized to the brightness temperature at the observation zenith angle of 0° by observing the vegetation coverage at the observation zenith angle of 0°.

[0093] Optionally, the brightness temperatures of the 31st and 32nd channels of Terra MODIS are angle-normalized based on the leaf area index, vegetation aggregation index data and the observation zenith angle to obtain the brightness temperature when the observation zenith angle is 0°.

[0094] Optionally, the surface emissivity of the Landsat-8 TIRS Channels 10 and 11 is calculated using the FR97 algorithm; the atmospheric water vapor content is estimated using the atmospheric relative humidity provided by the ERA5 reanalysis data.

[0095] Optionally, the brightness temperature simulation data on Landsat-8 is divided into two groups according to the average emissivity, the simulation data below a first set threshold is the first emissivity group, and the simulation data above a second set threshold is the second emissivity group; the first set threshold is greater than the second set threshold, and the second emissivity group and the second emissivity group containing overlapping areas are obtained; the surface temperature and the atmospheric water vapor content are grouped.

[0096] In actual application, the brightness temperature simulation data on Landsat-8 is divided into two groups according to the average emissivity: the simulated data below 0.96 is the low emissivity group, and the simulated data above 0.94 is the high emissivity group. An overlapping area of ​​0.02 is set between the high emissivity group and the low emissivity group, which can reduce the grouping error caused by the error in emissivity estimation and thus affect the accuracy of the algorithm. Similarly, the surface temperature and atmospheric water vapor content are grouped. The surface temperature is divided into 5 groups and an overlapping interval of 5 K is set: less than 280 K, between 275 K and 295 K, between 290 K and 310 K, between 305 K and 325 K, and greater than 320 K. Assume that the atmospheric water vapor content is WVC (Water Vapor Content, referred to as atmospheric water vapor content), the atmospheric water vapor content is divided into 6 groups and an overlapping interval of 0.5 is set. Overlap space: WVC less than 1.5 、WVC is located at 1 to 2.5 Between, WVC is located at 2 To 3.5 Between, WVC is located at 3 to 4.5 Between, WVC is located at 4 To 5.5 Between and WVC greater than Considering that the constructed Landsat-8 Channel-11 virtual channel has the same radiation path fitting coefficient as MODIS, let 𝑇 𝑖 The simulated brightness temperature of Channel-10 when the zenith angle is 0°, To simulate the brightness temperature of Channel-11 at the observed zenith angle 𝜃, is the average of the surface emissivity of channels i and j, is the difference between the surface emissivity of channel i and the surface emissivity of channel j. When inverting the surface temperature, The actual sensor brightness temperature of Channel-10 is obtained by interpolating the algorithm coefficients at 6 observation zenith angles to determine the algorithm coefficients of the target pixel at different observation zenith angles. The auxiliary parameter atmospheric water vapor content is estimated by the atmospheric relative humidity provided by the ERA5 reanalysis data, and the emissivity of the mixed pixel is estimated by the FR97 algorithm.

[0097] The process of inverting the surface temperature from the Landsat-8 TIRS thermal infrared data using the split window algorithm can be achieved by the following steps:

[0098] Firstly, the brightness temperature of the 31st and 32nd channels of Terra MODIS was normalized according to the leaf area index, vegetation aggregation index data and the observation zenith angle to obtain the brightness temperature when the observation zenith angle was 0°. The brightness temperature of the 31st and 32nd channels of Terra MODIS, the secant value of the observation zenith angle and the brightness temperature of the 10th channel of Landsat-8 TIRS when the observation zenith angle was 0° were input into the trained PSO-BP neural network to construct a virtual channel of the 11th channel of Landsat-8 TIRS. Then, the surface emissivity of the 10th and 11th channels of Landsat-8 TIRS was estimated by the FR97 algorithm, and the atmospheric water vapor content was estimated by the ERA reanalysis data.

[0099] Under the condition of known atmospheric water vapor content and surface emissivity, the split window coefficients at the 6 observation zenith angles in the sub-interval are determined only according to the grouping of average emissivity and atmospheric water vapor content. The split window coefficients at the observation zenith angle of the Terra MODIS pixel corresponding to the target Landsat-8 TIRS pixel are determined by interpolation, and an approximate value of the surface temperature is obtained through the split window algorithm formula. Finally, according to the approximate value of the surface temperature obtained in the previous step, the atmospheric water vapor content and the average surface emissivity, the split window coefficients of the appropriate sub-interval are reselected to estimate a more accurate surface temperature.

[0100] The present invention eliminates the influence of the surface temperature angle effect caused by different vegetation coverage projection areas at different observation zenith angles on ground mixed pixels by normalizing the angle of brightness temperature on the satellite, so that the observation data at different observation zenith angles are comparable.

[0101] The present invention realizes the time correction and spectrum conversion of the brightness temperature on the satellite through the PSO-BP neural network model to realize the construction of the virtual channel of the 11th channel of Landsat-8 TIRS.

[0102] Attached Figure 3 is the relationship curve between the true value and the predicted value. Figure 3 (a) is the relationship curve between the true value and the predicted value of the training set. Figure 3(b) is the relationship curve between the true value and the predicted value of the test set, with the horizontal axis being the predicted value in K and the vertical axis being the true value in K. The root mean square error of the training set is 0.2962 K, and the root mean square error of the test set is 0.2968 K. By constructing a virtual channel, the Landsat-8 satellite has two thermal infrared earth observation channels in the atmospheric window, and the split window algorithm can be used to invert the surface temperature.

[0103] The present invention adopts a group fitting method in the process of fitting the split window algorithm coefficients, which reduces the coefficient fitting error. In addition, considering the influence of different virtual channel observation paths, the fitting is performed according to the different observation zenith angles of the virtual channels to reduce the error.

[0104] As attached Figure 4 The simulation diagram of the root mean square error between the surface temperature inverted from different atmospheric water vapor content sub-intervals and the actual surface temperature is shown. Figure 4 (a) is a simulation diagram of the root mean square error corresponding to the Terra MODIS observation zenith angle of 0°. Figure 4 (b) is a simulation diagram of the root mean square error corresponding to the Terra MODIS observation zenith angle of 0°. The horizontal axis is the atmospheric water vapor content sub-interval, and the vertical axis is the root mean square error, unit: K. The root mean square of the coefficients of the fitted split window algorithm is less than 1K when the virtual channel observation zenith angle is 0°, and less than 1.2K when the observation zenith angle is 60°. Figure 5 The simulation diagram of the comparison between the inverted surface temperature and the ground station temperature is shown. Figure 5 (a) is a simulation diagram of the surface temperature of ground station 1. The horizontal axis is the inverted surface temperature, and the vertical axis is the surface temperature of station 1. Figure 5 (b) is a simulation diagram of the surface temperature at ground station 1. The horizontal axis is the inverted surface temperature, and the vertical axis is the surface temperature at station 2. The root mean square errors of the inverted surface temperatures are 3.23 K and 3.09 K, respectively.

[0105] Example 2

[0106] Based on the same principle as the method shown in Example 1 of the present invention, as shown in the attached Figure 6 As shown, an embodiment of the present invention also provides a ground surface temperature split window remote sensing inversion system based on a virtual channel, including a simulation data set construction unit, a virtual channel construction unit, a split window coefficient determination unit and a target ground surface temperature determination unit;

[0107] Simulated dataset building blocks for:

[0108] The atmospheric profile data under clear sky conditions were selected from the TIGR database, and the atmospheric profile data were input into the radiation transfer model to simulate the atmospheric parameters;

[0109] Convolve the ground object spectrum data with the sensor spectral response function to simulate the channel effective surface emissivity;

[0110] The atmospheric parameters and surface emissivity were input into the thermal infrared radiation transfer model to obtain the Terra MODIS on-board brightness temperature simulation data and the Landsat-8 TIRS on-board brightness temperature simulation data as simulation data sets;

[0111] Virtual channel building blocks for:

[0112] Normalize the brightness temperature angle on the satellite based on MODIS products, LAI products, vegetation aggregation index and MODIS L1B data;

[0113] A PSO-BP neural network model was constructed for time correction and spectral conversion. The brightness temperature simulation data of the 31st and 32nd channels of the Terra MODIS on-board brightness temperature simulation data at several observation zenith angles and the brightness temperature simulation data of the 10th channel of Landsat-8 TIRS at an observation zenith angle of 0° were used as input, and the brightness temperature simulation data of the 11th channel of Landsat-8 TIRS at the corresponding observation zenith angle was used as output. The PSO-BP neural network model was trained, and the virtual channel of the 11th channel of Landsat-8 TIRS was obtained based on the trained PSO-BP neural network model.

[0114] The split window coefficient determination unit is used to:

[0115] Based on the brightness temperature simulation data on Landsat-8 TIRS, the split window coefficient is determined according to the split window algorithm and the grouping rate of average surface emissivity, surface temperature and atmospheric water vapor content;

[0116] An approximate value of the surface temperature is determined by a split-window algorithm;

[0117] The target surface temperature determination unit is used to:

[0118] The splitting window coefficient is reselected according to the average surface emissivity, the approximate value of the surface temperature and the atmospheric water vapor content to determine the target surface temperature.

[0119] Optionally, the atmospheric parameters include channel transmittance, atmospheric up-radiation and atmospheric down-radiation.

[0120] Optionally, when training the PSO-BP neural network model, a set proportion of the simulated data set is used as the training set, and the rest is used as the test set; the input training set parameters are normalized to between 0 and 1; the parameters of the PSO-BP neural network model are initialized, and the PSO-BP neural network model trains the training set; the number of hidden layers of the PSO-BP neural network model is set, each hidden layer contains a number of hidden layer nodes, the learning rate, the number of training times and the target minimum error are set, and the quantized conjugate gradient method is used for optimization; the PSO algorithm is used for iterative optimization, the number of particle populations is 5, and the number of iterations, cognitive learning factor, social learning factor, inertia factor, maximum speed, minimum speed, maximum boundary and minimum boundary are set.

[0121] Optionally, for mixed pixels containing both soil and vegetation components, the directional effect of surface temperature is corrected according to vegetation cover, including:

[0122] Assume that the vegetation coverage is , the brightness temperature of the sensor's i-th channel is , Indicates the brightness temperature contributed by the vegetation component of the sensor's i-th channel, Represents the surface emissivity of the soil part in the mixed pixel, represents the brightness temperature contributed by the soil component of the sensor's channel i, represents the atmospheric transmittance, represents the surface emissivity of the vegetation component in the mixed pixel, represents the part of the atmospheric downward radiation in the i-th channel reflected by the soil components, represents the part of the atmospheric downward radiation in the i-th channel reflected by the vegetation components, It represents the contribution of vegetation when the atmosphere of the i-th channel is divided into vegetation and soil according to vegetation coverage. It represents the soil contribution when the atmosphere of the i-th channel is divided into vegetation and soil according to vegetation coverage, then:

[0123]

[0124]

[0125]

[0126] ;

[0127] set up is the base of natural logarithms, is the vegetation aggregation index, is the leaf area index, is the leaf area orientation function; the leaf area orientation function is the average projection of the unit leaf area in the viewing direction, which is used to describe the effect of leaf inclination on radiation attenuation, then:

[0128] ;

[0129] Under the condition that the vegetation aggregation index and leaf area index are known, the value of vegetation coverage depends on the observation zenith angle. Assuming that the pixels in the 3×3 spatial window adjacent to the target pixel have the same component temperature, including vegetation temperature and soil temperature, and the atmospheric environment does not change, for the nine pixels in the spatial window, the equation group is constructed:

[0130] ;

[0131] According to the vegetation coverage of each pixel at the current observation zenith angle, the brightness temperature contributed by each component is determined, and then the brightness temperature is normalized to the brightness temperature at the observation zenith angle of 0° by observing the vegetation coverage at the observation zenith angle of 0°.

[0132] Optionally, the brightness temperatures of the 31st and 32nd channels of Terra MODIS are angle-normalized based on the leaf area index, vegetation aggregation index data and the observation zenith angle to obtain the brightness temperature when the observation zenith angle is 0°.

[0133] Optionally, the surface emissivity of the Landsat-8 TIRS Channels 10 and 11 is calculated using the FR97 algorithm; the atmospheric water vapor content is estimated using the atmospheric relative humidity provided by the ERA5 reanalysis data.

[0134] Optionally, the brightness temperature simulation data on Landsat-8 is divided into two groups according to the average emissivity, the simulation data below a first set threshold is the first emissivity group, and the simulation data above a second set threshold is the second emissivity group; the first set threshold is greater than the second set threshold, and the second emissivity group and the second emissivity group containing overlapping areas are obtained; the surface temperature and the atmospheric water vapor content are grouped.

[0135] The above are only preferred embodiments of the present invention and are not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A split-window remote sensing inversion method for land surface temperature based on virtual channels, characterized in that: It includes constructing simulation data set, constructing virtual channel, determining split window coefficient and determining target surface temperature; Construct a simulated data set, including: The atmospheric profile data under clear sky conditions were selected from the TIGR database, and the atmospheric profile data were input into the radiation transfer model to simulate the atmospheric parameters; Convolve the ground object spectrum data with the sensor spectral response function to simulate the channel effective surface emissivity; The atmospheric parameters and surface emissivity were input into the thermal infrared radiation transfer model to obtain the Terra MODIS on-board brightness temperature simulation data and the Landsat-8 TIRS on-board brightness temperature simulation data as simulation data sets; Build a virtual channel, including: The brightness temperature angle is normalized on the satellite based on MODIS products, LAI products, vegetation aggregation index and MODIS L1B data; for mixed pixels containing soil components and vegetation components, the directional effect of the surface temperature is corrected according to the vegetation coverage, including: Assume that the vegetation coverage is , the brightness temperature of the sensor's i-th channel is , Indicates the brightness temperature contributed by the vegetation component of the sensor's i-th channel, Represents the surface emissivity of the soil part in the mixed pixel, represents the brightness temperature contributed by the soil component of the sensor's channel i, represents the atmospheric transmittance, represents the surface emissivity of the vegetation component in the mixed pixel, represents the part of the atmospheric downward radiation in the i-th channel reflected by the soil components, represents the part of the atmospheric downward radiation in the i-th channel reflected by the vegetation components, It represents the contribution of vegetation when the atmosphere of the i-th channel is divided into vegetation and soil according to vegetation coverage. It represents the soil contribution when the atmosphere of the i-th channel is divided into vegetation and soil according to vegetation coverage, then: ； set up is the base of natural logarithms, is the vegetation aggregation index, is the leaf area index, is the leaf area orientation function; the leaf area orientation function is the average projection of the unit leaf area in the viewing direction, which is used to describe the effect of leaf inclination on radiation attenuation, then: ； Under the condition that the vegetation aggregation index and leaf area index are known, the value of vegetation coverage depends on the observation zenith angle. Assuming that the pixels in the 3×3 spatial window adjacent to the target pixel have the same component temperature, including vegetation temperature and soil temperature, and the atmospheric environment does not change, for the nine pixels in the spatial window, the equation group is constructed: ； According to the vegetation coverage of each pixel at the current observation zenith angle, the brightness temperature contributed by each component is determined, and then the brightness temperature is normalized to the brightness temperature at the observation zenith angle of 0° by observing the vegetation coverage at the observation zenith angle of 0°; A PSO-BP neural network model was constructed for time correction and spectral conversion. The on-board brightness temperature simulation data of the 31st and 32nd channels of the TerraMODIS on-board brightness temperature simulation data at several observation zenith angles and the on-board brightness temperature simulation data of the 10th channel of Landsat-8 TIRS at an observation zenith angle of 0° were used as input, and the on-board brightness temperature simulation data of the 11th channel of Landsat-8 TIRS at the corresponding observation zenith angle was used as output. The PSO-BP neural network model was trained, and the virtual channel of the 11th channel of Landsat-8 TIRS was obtained based on the trained PSO-BP neural network model. Determine the split window coefficient, including: Based on the brightness temperature simulation data on Landsat-8 TIRS, the split window coefficient is determined according to the split window algorithm and the grouping rate of average surface emissivity, surface temperature and atmospheric water vapor content; An approximate value of the surface temperature is determined by a split-window algorithm; Determine target surface temperature, including: The splitting window coefficient is reselected according to the average surface emissivity, the approximate value of the surface temperature and the atmospheric water vapor content to determine the target surface temperature.

2. The method for remote sensing inversion of land surface temperature by split window based on virtual channel according to claim 1 is characterized in that: Atmospheric parameters include channel transmittance, atmospheric upwelling radiation and atmospheric downwelling radiation.

3. The method for remote sensing inversion of land surface temperature using split windows based on virtual channels according to claim 1 is characterized in that: When training the PSO-BP neural network model, a set proportion of the simulated data set is used as the training set, and the rest is used as the test set; the input training set parameters are normalized to between 0 and 1; the parameters of the PSO-BP neural network model are initialized, and the PSO-BP neural network model trains the training set; the number of hidden layers of the PSO-BP neural network model is set, each hidden layer contains a number of hidden layer nodes, the learning rate, the number of training times and the target minimum error are set, and the quantized conjugate gradient method is used for optimization; the PSO algorithm is used to iteratively search for the best, the number of particle populations is 5, and the number of iterations, cognitive learning factor, social learning factor, inertia factor, maximum speed, minimum speed, maximum boundary and minimum boundary are set.

4. The method for remote sensing inversion of land surface temperature using split windows based on virtual channels according to claim 1 is characterized in that: According to the leaf area index, vegetation aggregation index data and the observation zenith angle, the brightness temperature of the 31st and 32nd channels of Terra MODIS was angle-normalized to obtain the brightness temperature when the observation zenith angle was 0°.

5. The method for remote sensing inversion of land surface temperature using split windows based on virtual channels according to claim 1, characterized in that: The surface emissivity of the 10th and 11th channels of Landsat-8 TIRS was calculated using the FR97 algorithm. The atmospheric water vapor content was estimated using the atmospheric relative humidity provided by the ERA5 reanalysis data.

6. The method for remote sensing inversion of ground surface temperature using split windows based on virtual channels according to claim 1 is characterized in that: The brightness temperature simulation data on Landsat-8 are divided into two groups according to the average emissivity. The simulation data below the first set threshold are the first emissivity group, and the simulation data above the second set threshold are the second emissivity group. When the first set threshold is greater than the second set threshold, the second emissivity group and the second emissivity group containing overlapping areas are obtained. The surface temperature and atmospheric water vapor content are grouped.

7. The surface temperature split window remote sensing inversion system based on virtual channel is characterized by: It includes a simulation data set construction unit, a virtual channel construction unit, a split window coefficient determination unit and a target surface temperature determination unit; Simulated dataset building blocks for: The atmospheric profile data under clear sky conditions were selected from the TIGR database, and the atmospheric profile data were input into the radiation transfer model to simulate the atmospheric parameters; Convolve the ground object spectrum data with the sensor spectral response function to simulate the channel effective surface emissivity; The atmospheric parameters and surface emissivity were input into the thermal infrared radiation transfer model to obtain the Terra MODIS on-board brightness temperature simulation data and the Landsat-8 TIRS on-board brightness temperature simulation data as simulation data sets; Virtual channel building blocks for: The brightness temperature angle is normalized on the satellite based on MODIS products, LAI products, vegetation aggregation index and MODIS L1B data; for mixed pixels containing soil components and vegetation components, the directional effect of the surface temperature is corrected according to the vegetation coverage, including: Assume that the vegetation coverage is , the brightness temperature of the sensor's i-th channel is , Indicates the brightness temperature contributed by the vegetation component of the sensor's i-th channel, Represents the surface emissivity of the soil part in the mixed pixel, represents the brightness temperature contributed by the soil component of the sensor's channel i, represents the atmospheric transmittance, represents the surface emissivity of the vegetation component in the mixed pixel, represents the part of the atmospheric downward radiation in the i-th channel reflected by the soil components, represents the part of the atmospheric downward radiation in the i-th channel reflected by the vegetation components, It represents the contribution of vegetation when the atmosphere of the i-th channel is divided into vegetation and soil according to vegetation coverage. It represents the soil contribution when the atmosphere of the i-th channel is divided into vegetation and soil according to vegetation coverage, then: ； set up is the base of natural logarithms, is the vegetation aggregation index, is the leaf area index, is the leaf area orientation function; the leaf area orientation function is the average projection of the unit leaf area in the viewing direction, which is used to describe the effect of leaf inclination on radiation attenuation, then: ； Under the condition that the vegetation aggregation index and leaf area index are known, the value of vegetation coverage depends on the observation zenith angle. Assuming that the pixels in the 3×3 spatial window adjacent to the target pixel have the same component temperature, including vegetation temperature and soil temperature, and the atmospheric environment does not change, for the nine pixels in the spatial window, the equation group is constructed: ； According to the vegetation coverage of each pixel at the current observation zenith angle, the brightness temperature contributed by each component is determined, and then the brightness temperature is normalized to the brightness temperature at the observation zenith angle of 0° by observing the vegetation coverage at the observation zenith angle of 0°; A PSO-BP neural network model was constructed for time correction and spectral conversion. The on-board brightness temperature simulation data of the 31st and 32nd channels of the TerraMODIS on-board brightness temperature simulation data at several observation zenith angles and the on-board brightness temperature simulation data of the 10th channel of Landsat-8 TIRS at an observation zenith angle of 0° were used as input, and the on-board brightness temperature simulation data of the 11th channel of Landsat-8 TIRS at the corresponding observation zenith angle was used as output. The PSO-BP neural network model was trained, and the virtual channel of the 11th channel of Landsat-8 TIRS was obtained based on the trained PSO-BP neural network model. The split window coefficient determination unit is used to: Based on the brightness temperature simulation data on Landsat-8 TIRS, the split window coefficient is determined according to the split window algorithm and the grouping rate of average surface emissivity, surface temperature and atmospheric water vapor content; An approximate value of the surface temperature is determined by a split-window algorithm; The target surface temperature determination unit is used to: The splitting window coefficient is reselected according to the average surface emissivity, the approximate value of the surface temperature and the atmospheric water vapor content to determine the target surface temperature.

Images