CN104897289A - Landsat 8 satellite data land surface temperature inversion method - Google Patents

Abstract

A Landsat 8 data surface temperature inversion method, which is completely based on the Landsat 8 data itself does not require any external data sources, overcomes the limitations caused by the traditional Landsat data surface temperature inversion must rely on external data sources, the invention is useful for realizing the use of Landsat 8 It is of great practical significance to commercially produce land surface temperature products from data.

Description

A method for inversion of land surface temperature from Landsat 8 satellite data

technical field

The invention relates to a method for retrieving surface temperature from Landsat 8 satellite data, which can be applied in forestry, agriculture, meteorology, ecological environment and other industry sectors.

Background technique

Surface temperature is a key parameter for studying regional energy exchange and water cycle, and it is also an important input parameter for process models such as ecology, hydrology and climate. Obtaining regional surface temperature is an important part of dynamic monitoring of regional resources and environment. Thermal infrared satellite remote sensing technology is an important way to obtain regional surface temperature. Landsat 8 data is a new type of satellite data source. Compared with the traditional Landsat series satellites (Landsat5, 7), Landsat 8 has improved the number of bands, the spectral range of the bands and the radiometric resolution of the images. Landsat8 carries two sensors: 1) Operational LandImager (OLI) and Thermal Infrared Sensor (TIRS). The OLI sensor receives data from nine spectral bands in the visible light, near-infrared and short-wave infrared regions; the TIRS sensor divides the original Landsat5 and 7 thermal infrared bands into two, and sets them into two thermal infrared channels (Band 10: 10.6-11.19μm ; Band 11: 11.5-12.51 μm). For Landsat5 and 7 surface temperature retrieval, a single-channel surface temperature retrieval algorithm is usually used. This algorithm requires at least two input parameters: atmospheric water vapor content and specific emissivity. The specific emissivity can be obtained from the Landsat data itself using the NDVI threshold method. However, atmospheric water vapor content must rely on external data sources, usually obtained indirectly through meteorological data or MODIS data. However, both meteorological data and MODIS data have obvious limitations: meteorological data is a kind of point data, and Remote sensing data is a kind of surface data, and the method of replacing meteorological data with points will lead to large errors, and it is very difficult to obtain corresponding meteorological data for remote areas or historical archived satellite data; There is a large difference in spatial resolution and geometric registration and projection transformation between the two data will also bring errors. What's more, for most of China, the geographic overlap area between Landsat data and MODIS data is often very small (less than one-third), and even MODIS data corresponding to Landsat data cannot be found. These defects have caused great difficulties to the traditional Landsat surface temperature retrieval. Fortunately, the band setting of Landsat8 makes it possible to invert the land surface temperature based on the Landsat8 data itself without using external data sources. For Landsat8 data, the specific emissivity can also be obtained from the Landsat8 data itself using the NDVI threshold method, and the atmospheric water vapor content can be inverted using the two thermal infrared channels of Landsat8 based on the split window covariance-variance ratio algorithm, so that The realization is completely based on Landsat8 data itself does not need any external data to invert the surface temperature. This invention has important practical significance for realizing the commercial production of land surface temperature products using Landsat8 data.

Contents of the invention

The purpose of the present invention is to provide a method for inversion of land surface temperature from Landsat 8 satellite data, which is completely based on Landsat 8 data itself and does not require any external data, and has very strong practicability.

To achieve the above object, the method proposed by the present invention comprises the following steps:

The first step is to calculate the on-board radiance and on-board brightness temperature of the 10th and 11th bands of Landsat 8

L _sen =M _L Q _cal +A _L

T _sen =K ₂ /ln(1+K ₁ /L _sen )

Among them, L _sen is the radiance of the star, T _sen is the brightness temperature of the star, M _L is the gain of the band, _AL is the offset of the band, Q _cal is the DN value of the image, K1 and K2 are constants, M _L , A _L and K1 and K2 are obtained from the Landsat 8 header file;

The second step is to use the NDVI (Normalized Difference Vegetation Index) threshold method to obtain the specific emissivity ε:

$\mathrm{<span\; class="notranslate">\; NDVI</span>}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; DN</span>}}_{\mathrm{<span\; class="notranslate">\; band</span>}\mathrm{<span\; class="notranslate">\; 5</span>}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; D.</span>}{\mathrm{<span\; class="notranslate">\; N</span>}}_{\mathrm{<span\; class="notranslate">\; band</span>}\mathrm{<span\; class="notranslate">\; 4</span>}}}{{\mathrm{<span\; class="notranslate">\; DN</span>}}_{\mathrm{<span\; class="notranslate">\; band</span>}\mathrm{<span\; class="notranslate">\; 5</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; DN</span>}}_{\mathrm{<span\; class="notranslate">\; band</span>}\mathrm{<span\; class="notranslate">\; 4</span>}}}$

Among them, DN _band5 and DN _band4 represent the DN values of Landsat8's 5th band and 4th band images respectively;

When NDVI<NDVI _s , ε=ε _s , where NDVI _s is the NDVI of the pure bare soil area, and ε _s is the specific emissivity of the soil;

When NDVI > NDVI _v , ε = ε _v , where NDVI _v is the NDVI of the pure vegetation area, and ε _v is the specific emissivity of vegetation;

When NDVI _s ≤ NDVI ≤ NDVI _v , ε = ε _s (1-FVC) + ε _v FVC

FVC is vegetation coverage:

$\mathrm{<span\; class="notranslate">\; FVC</span>}<span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; [</span>\frac{\mathrm{<span\; class="notranslate">\; NDVI</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; NDV</span>}{\mathrm{<span\; class="notranslate">\; I</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}}{{\mathrm{<span\; class="notranslate">\; NDVI</span>}}_{\mathrm{<span\; class="notranslate">\; v</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; NDVI</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}}<span\; class="notranslate">\; ]</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}$

NDVI _s and NDVI _v can be obtained by selecting homogeneous bare soil areas and vegetation areas from the image; ε _s and ε _v are calculated by MODIS UCSB specific emissivity library and Landsat 8 TIRS spectral response function;

Step 3: Calculate the atmospheric water vapor content w

w=a(τ _j /τ _i )+b

$\frac{{\mathrm{\&tau;}}_j}{{\mathrm{\&tau;}}_i}=\frac{{\mathrm{\&epsiv;}}_i}{{\mathrm{\&epsiv;}}_j}{R}_{j,i}$ and $R_{j,i}=\frac{\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}(T_{i,k}-{\stackrel{\&OverBar;}{T}}_i)(T_{j,k}-{\stackrel{\&OverBar;}{T}}_j)}{\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{(T_{j,k}-{\stackrel{\&OverBar;}{T}}_i)}^2}$

Among them, τ _i is the atmospheric transmittance of the i-band, τ _j is the atmospheric transmittance of the j-band, ε _i is the specific emissivity of the i-band, ε _j is the specific emissivity of the j-band, k represents the kth image element, T _{i, k} is the on-board brightness temperature of the k-th pixel in band i, T _{j, k} is the on-board brightness temperature of the k-th pixel in j-band, is the average on-star brightness temperature of the i-band of N pixels, is the average on-star brightness temperature of N pixel j-band, for Landsat8 data, i, j are 10, 11 respectively, N represents the window size, take 20 pixels*20 pixels;

The coefficients a and b use the MODTRAN4.0 atmospheric radiative transfer model and the TIGR database to simulate the relationship between the atmospheric water vapor content w and the Landsat8 thermal infrared band atmospheric transmittance ratio τ ₁₁ /τ ₁₀ to get:

w=-18.973(τ ₁₁ /τ ₁₀ )+19.13 R ² =0.9663, τ ₁₁ /τ ₁₀ >0.9

w=-13.412(τ ₁₁ /τ ₁₀ )+14.158 R ² =0.9366, τ ₁₁ /τ ₁₀ <0.9

Step 4: Calculate the surface temperature

${\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; \&gamma;</span>}<span\; class="notranslate">\; [</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&psi;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; L</span>}}_{\mathrm{<span\; class="notranslate">\; sen</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&psi;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&psi;</span>}}_{\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; ]</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; \&delta;</span>}$

Where T _s is the surface temperature, ε is the specific emissivity, L _sen is the radiance of the star in the 10th band of Landsat8, (γ, δ) can be expressed as:

$\mathrm{<span\; class="notranslate">\; \&gamma;</span>}<span\; class="notranslate">\; \&ap;</span>\frac{{\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; sen</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{{\mathrm{<span\; class="notranslate">\; b</span>}}_{\mathrm{<span\; class="notranslate">\; \&gamma;</span>}}{\mathrm{<span\; class="notranslate">\; L</span>}}_{\mathrm{<span\; class="notranslate">\; sen</span>}}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; \&delta;</span>}<span\; class="notranslate">\; \&ap;</span>{\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; sen</span>}}<span\; class="notranslate">\; -</span>\frac{{\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; sen</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{{\mathrm{<span\; class="notranslate">\; b</span>}}_{\mathrm{<span\; class="notranslate">\; \&gamma;</span>}}}$

Where T _sen is the brightness temperature of the star in the 10th band of Landsat8, b _γ is equal to 1324K, ψ ₁ , ψ ₂ , and ψ ₃ are atmospheric functions, which can be approximated from the atmospheric water vapor content (w) by the following formula:

ψ ₁ =0.04019w ² +0.02916w+1.01523

ψ ₂ =-0.38333w ² -1.50294w+0.20324

ψ ₃ =0.00918w ² +1.36072w-0.27514

Description of drawings

Figure 1 The relationship between the ratio of atmospheric transmittance in the thermal infrared band of Landsat 8 and the content of atmospheric water vapor

Detailed ways

The present invention uses the single-channel surface temperature inversion method to invert the surface temperature from the 10th band of Landsat 8. The single-channel method is obtained based on the simplification of the thermal band radiation transfer equation, which can be expressed as:

${\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; \&gamma;</span>}<span\; class="notranslate">\; [</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&psi;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; L</span>}}_{\mathrm{<span\; class="notranslate">\; sen</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&psi;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&psi;</span>}}_{\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; ]</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; \&delta;</span>}$

Where T _s is the surface temperature, ε is the specific emissivity, L _sen is the radiance of the star, (γ, δ) can be expressed as:

$\mathrm{<span\; class="notranslate">\; \&gamma;</span>}<span\; class="notranslate">\; \&ap;</span>\frac{{\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; sen</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{{\mathrm{<span\; class="notranslate">\; b</span>}}_{\mathrm{<span\; class="notranslate">\; \&gamma;</span>}}{\mathrm{<span\; class="notranslate">\; L</span>}}_{\mathrm{<span\; class="notranslate">\; sen</span>}}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; \&delta;</span>}<span\; class="notranslate">\; \&ap;</span>{\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; sen</span>}}<span\; class="notranslate">\; -</span>\frac{{\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; sen</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{{\mathrm{<span\; class="notranslate">\; b</span>}}_{\mathrm{<span\; class="notranslate">\; \&gamma;</span>}}}$

where T _sen is the brightness temperature on the star, b _γ is equal to 1324K, ψ ₁ , ψ ₂ , and ψ ₃ are atmospheric functions, which can be approximated from the atmospheric water vapor content (w) using the following formula (Jiménez- JC, Sobrino, JA, D, Mattar C, and J. (2014). Land Surface Temperature Retrieval Methods From Landsat-8 Thermal Infrared Sensor Data. IEEE Geoscience and Remote Sensing Letters, 11(10), 1840-1843.):

ψ ₁ =0.04019w ² +0.02916w+1.01523

ψ ₂ =-0.38333w ² -1.50294w+0.20324

ψ ₃ =0.00918w ² +1.36072w-0.27514

L _sen =M _L Q _cal +A _L

_ML is the gain of the band, _{AL is the offset of the band, ML and AL} are obtained from the Landsat 8 header file, and Q _cal is the DN value of the image.

T _sen =K ₂ /ln(1+K ₁ /L _sen )

K1 and K2 are constants obtained from Landsat 8 header files.

The specific emissivity is obtained using the NDVI (Normalized Difference Vegetation Index) threshold method:

$\mathrm{<span\; class="notranslate">\; NDVI</span>}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; DN</span>}}_{\mathrm{<span\; class="notranslate">\; band</span>}\mathrm{<span\; class="notranslate">\; 5</span>}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; D.</span>}{\mathrm{<span\; class="notranslate">\; N</span>}}_{\mathrm{<span\; class="notranslate">\; band</span>}\mathrm{<span\; class="notranslate">\; 4</span>}}}{{\mathrm{<span\; class="notranslate">\; DN</span>}}_{\mathrm{<span\; class="notranslate">\; band</span>}\mathrm{<span\; class="notranslate">\; 5</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; DN</span>}}_{\mathrm{<span\; class="notranslate">\; band</span>}\mathrm{<span\; class="notranslate">\; 4</span>}}}$

Among them, DN _band5 and DN _band4 represent the DN values of Landsat8's 5th band and 4th band images, respectively.

When NDVI<NDVI _s , ε=ε _s , where NDVI _s is the NDVI of the pure bare soil area, and ε _s is the specific emissivity of the soil;

When NDVI > NDVI _v , ε = ε _v , where NDVI _v is the NDVI of the pure vegetation area, and ε _v is the specific emissivity of vegetation;

When NDVI _s ≤ NDVI ≤ NDVI _v , ε = ε _s (1-FVC) + ε _v FVC

FVC is vegetation coverage:

$\mathrm{<span\; class="notranslate">\; FVC</span>}<span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; [</span>\frac{\mathrm{<span\; class="notranslate">\; NDVI</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; NDV</span>}{\mathrm{<span\; class="notranslate">\; I</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}}{{\mathrm{<span\; class="notranslate">\; NDVI</span>}}_{\mathrm{<span\; class="notranslate">\; v</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; NDVI</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}}<span\; class="notranslate">\; ]</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}$

NDVI _s and NDVI _v can be obtained by selecting homogeneous bare soil areas and vegetation areas from the image. ε _s and ε _v are calculated by MODIS UCSB specific emissivity library and Landsat 8 TIRS spectral response function.

Atmospheric water vapor content (w) based on split window covariance-variance ratio algorithm (SOBRINO J A, Li Z L, Stoll MP, et al. Improvements in the split-window technique for land surface temperature determination[J]. Geoscience and Remote Sensing, IEEE Transactions on, 1994, 32(2): 243-253.) to invert, the algorithm assumes that under cloudless conditions, within N adjacent pixel areas (for Landsat 8, N can take a value of 20 , that is, the window size is 20 pixels*20 pixels), the atmospheric conditions and specific emissivity do not change, only the surface temperature changes, and w is calculated by the following formula:

w=a(τ _j /τ _i )+b (1)

$\frac{{\mathrm{\&tau;}}_j}{{\mathrm{\&tau;}}_i}=\frac{{\mathrm{\&epsiv;}}_i}{{\mathrm{\&epsiv;}}_j}{R}_{j,i}$ and $R_{j,i}=\frac{\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}(T_{i,k}-{\stackrel{\&OverBar;}{T}}_i)(T_{j,k}-{\stackrel{\&OverBar;}{T}}_j)}{\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{(T_{j,k}-{\stackrel{\&OverBar;}{T}}_i)}^2}---\left(2\right)$

Among them, τ _i is the atmospheric transmittance of the i-band, τ _j is the atmospheric transmittance of the j-band, ε _i is the specific emissivity of the i-band, ε _j is the specific emissivity of the j-band, k represents the kth image element, T _{i, k} is the on-board brightness temperature of the k-th pixel in band i, T _{j, k} is the on-board brightness temperature of the k-th pixel in j-band, is the average on-star brightness temperature of the i-band of N pixels, is the average on-star brightness temperature of the j-band of N pixels. For Landsat8 data, i, j are 10, 11 respectively.

For the Landsat8 TIRS data, formulas (1) and (2) are used to invert the atmospheric water vapor content, and the coefficients a and b need to be determined. The coefficients a and b can simulate the atmospheric water vapor content and the thermal infrared band atmospheric transmittance through the atmospheric radiation transfer model The relationship of the ratio is solved.

The MODTRAN4.0 atmospheric radiative transfer model and the TIGR (Thermodynamic Initial Guess Retrieval, TIGR) database were used to simulate the relationship between the atmospheric water vapor content w and the Landsat8 thermal infrared band atmospheric transmittance ratio τ ₁₁ /τ ₁₀ . The TIGR database is a meteorological database consisting of 2311 atmospheric profiles; each profile data contains the pressure, air temperature, water vapor content and ozone content of each layer from the surface to the top of the atmosphere. The TIGR database includes 872 tropical atmospheric profiles, 742 mid-latitude atmospheric profiles and 697 high-latitude atmospheric profiles. A wide range of atmospheric water vapor contents (from 0.066 to 7.833 g/cm ² ) is included in the TIGR database. The TIGR database is used as the input of the MODTRAN4.0 model to simulate the relationship between the atmospheric water vapor content w and the thermal infrared band atmospheric transmittance ratio τ ₁₁ /τ ₁₀ . Figure 1 shows the relationship between the Landsat 8 thermal infrared band atmospheric transmittance ratio and atmospheric water vapor content based on 2311 TIGR atmospheric profiles and the MODTRAN4.0 atmospheric radiative transfer model.

Figure 1 The relationship between the ratio of atmospheric transmittance in the thermal infrared band of Landsat 8 and the content of atmospheric water vapor

As shown in Figure 1, there is a good correlation between the atmospheric transmittance ratio of the 11th band and the 10th band of Landsat8 data and the atmospheric water vapor content. It can be seen from Figure 1 that there is an inflection point at the transmittance ratio of 0.9. In order to better fit the relationship between the atmospheric transmittance ratio and atmospheric water vapor content, the The rate ratio is divided into two sections for fitting, and the relationship between the atmospheric transmittance ratio and the atmospheric water vapor content is obtained:

w=-18.973(τ ₁₁ /τ ₁₀ )+19.13 R ² =0.9663, τ ₁₁ /τ ₁₀ >0.9 (3)

w=-13.412(τ ₁₁ /τ ₁₀ )+14.158 R ² =0.9366, τ ₁₁ /τ ₁₀ <0.9 (4)

The coefficients a and b can be obtained from formula 3 and formula 4.

Claims (1)

1. A Landsat 8 satellite data surface temperature retrieval method, the steps are: The first step is to calculate the on-board radiance and on-board brightness temperature of the 10th and 11th bands of Landsat 8 L _sen =M _L Q _cal +A _L T _sen =K ₂ /ln(1+K ₁ /L _sen ) Among them, L _sen is the radiation brightness on the star, T _sen is the brightness temperature on the star, _ML is the gain of the band, AL is the bias of the band, Q _cal is the image DN value, K1 and K2 are constants, _ML , AL and K1 and K2 are obtained from the Landsat 8 header; The second step is to use the NDVI (Normalized Difference Vegetation Index) threshold method to obtain the specific emissivity ε: Among them, DN _band5 and DN _band4 represent the DN values of Landsat8's 5th band and 4th band images respectively; When NDVI<NDVI _s , ε=ε _s , where NDVI _s is the NDVI of the pure bare soil area, and ε _s is the specific emissivity of the soil; When NDVI > NDVI _v , ε = ε _v , where NDVI _v is the NDVI of the pure vegetation area, and ε _v is the specific emissivity of vegetation; When NDVI _s ≤ NDVI ≤ NDVI _v , ε = ε _s (1-FVC) + ε _v FVC FVC is vegetation coverage: NDVI _s and NDVI _v can be obtained by selecting homogeneous bare soil areas and vegetation areas from the image; ε _s and ε _v are calculated by MODIS UCSB specific emissivity library and Landsat 8TIRS spectral response function; Step 3: Calculate the atmospheric water vapor content w w=a(τ _j /τ _i )+b and Among them, τ _i is the atmospheric transmittance of the i-band, τ _j is the atmospheric transmittance of the j-band, εi is the specific emissivity of the i-band, ε _j is the specific emissivity of the j-band, and k represents the kth pixel , T _{i, k} is the on-board brightness temperature of the k-th pixel in band i, T _{j, k} is the on-board brightness temperature of the k-th pixel in j-band, is the average on-star brightness temperature of the i-band of N pixels, is the average on-star brightness temperature of N pixel j-band, for Landsat8 data, i, j are 10, 11 respectively, N represents the window size, take 20 pixels*20 pixels; The coefficients a and b use the MODTRAN4.0 atmospheric radiative transfer model and the TIGR database to simulate the relationship between the atmospheric water vapor content w and the Landsat8 thermal infrared band atmospheric transmittance ratio τ ₁₁ /τ ₁₀ to get: w=-18.973(τ ₁₁ /τ ₁₀ )+19.13 R ² =0.9663, τ ₁₁ /τ ₁₀ >0.9 w=-13.412(τ ₁₁ /τ ₁₀ )+14.158 R ² =0.9366, τ ₁₁ /τ ₁₀ <0.9 Step 4: Calculate the surface temperature where T _s is the surface temperature, ε is the specific emissivity, L _sen is the onboard radiance of the 10th band of Landsat8, (γ, 6) can be expressed as: Where T _sen is the brightness temperature of the star in the 10th band of Landsat8, b _γ is equal to 1324K, ψ ₁ , ψ ₂ , and ψ ₃ are atmospheric functions, which can be approximated from the atmospheric water vapor content (w) by the following formula: ψ ₁ =0.04019w ² +0.02916w+1.01523 ψ ₂ =-0.38333w ² -1.50294w+0.20324 ψ ₃ =0.00918w ² +1.36072w-0.27514.