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

Abstract

A Landsat 8 satellite data land surface temperature inversion method. The method is completely based on Landsat 8 data and does not need any external data source, thereby overcomes the limitation that conventional Landsat data land surface temperature inversion must rely on an external data source, and the method provided by the invention is of great practical significance to realization of using Landsat 8 data to produce land surface temperature products in an operational manner.

Description

A kind of Landsat 8 satellite data Surface Temperature Retrieval method

Technical field

The present invention relates to a kind of method from Landsat 8 satellite data inverting surface temperature, the industry departments such as forestry, agricultural, meteorology, ecologic environment can be applied in.

Background technology

Surface temperature is the key parameter of survey region energy exchange and Water Cycle, is also the important input parameter of of the process model such as ecological, the hydrology and weather.Obtain the important content that region surface temperature is region resource environmental dynamic monitor.Thermal infrared satellite remote sensing technology obtains one of region surface temperature very important channel.Landsat 8 data are a kind of novel satellite data sources, and compared with traditional Landsat series of satellites (Landsat5,7), Landsat 8 improves on the radiometric resolution of the quantity of wave band, the spectral range of wave band and image.Landsat8 carries two sensors: 1) Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS).OLI sensor is in the data of visible ray, near infrared and short-wave infrared areas accept nine spectral bands; TIRS sensor by original Landsat5,7 Thermal infrared bands be divided into two, be arranged to two Detection Using Thermal Infrared Channel (Band 10:10.6-11.19 μm; Band11:11.5-12.51 μm).For Landsat5, 7 Surface Temperature Retrieval, usually single channel Surface Temperature Retrieval algorithm is utilized, this algorithm at least needs two input parameters: atmospheric water vapor content and emissivity, emissivity can utilize NDVI threshold method to obtain from Landsat data itself, atmospheric water vapor content then must rely on external data source, usually indirectly obtained by weather data or MODIS data, but no matter utilizing weather data or MODIS data all have obvious limitation: weather data is a kind of point data, and remotely-sensed data is a kind of face data, the mode of weather data Points replacing surfaces can cause larger error, and for remote districts or historical archive satellite data, obtain corresponding weather data just very difficult, there is larger difference in MODIS data and Landsat data, the geometrical registration between two kinds of data and projection transform also can bring error on imaging time and spatial resolution.The more important thing is, for most of China area, the geographical overlapping region between Landsat data and MODIS data is very little (being less than 1/3rd) often, even can not find the MODIS data corresponding with Landsat data.These defects cause very large difficulty to traditional Landsat Surface Temperature Retrieval above.Fortunately, the band setting of Landsat8 is given based on Landsat8 data itself, is carried out inverting surface temperature bring possibility without using external data source.For Landsat8 data, emissivity can utilize NDVI threshold method to obtain from Landsat8 data itself equally, and atmospheric water vapor content can utilize two of Landsat8 Detection Using Thermal Infrared Channels to carry out inverting based on splitting window covariance-variance ratio algorithm, so just can realize carrying out inverting surface temperature based on Landsat8 data itself without any need for outsource data completely.This invention has important practical significance for realizing with utilizing Landsat8 data service production surface temperature product.

Summary of the invention

The object of the present invention is to provide a kind of Landsat 8 satellite data Surface Temperature Retrieval method, the method completely based on Landsat8 data itself without any need for external data, practicality is very strong.

For achieving the above object, the method that the present invention proposes comprises the following steps:

The first step, calculate Landsat 8 the 10th wave band and the 11st wave band star on brightness temperature on radiance and star

L _sen＝M _LQ _cal+A _L

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

Wherein, L _senradiance on star, T _senbrightness temperature on star, M _lfor the gain of wave band, A _lfor being biased of wave band, Q _calfor image DN value, K1 and K2 is constant, M _l, A _land K1 and K2 obtains from Landsat 8 header file;

Second step, utilize NDVI (Normalized Difference Vegetation Index) threshold method to obtain emissivity ε:

$\mathrm{NDVI}=\frac{{\mathrm{DN}}_{\mathrm{band}5}-D{N}_{\mathrm{band}4}}{{\mathrm{DN}}_{\mathrm{band}5}+{\mathrm{DN}}_{\mathrm{band}4}}$

Wherein DN _band5and DN _band4represent the DN value of Landsat8 the 5th wave band and the 4th wave band image respectively;

As NDVI < NDVI _stime, ε=ε _s, wherein NDVI _sthe NDVI in pure exposed soil region, ε _sit is the emissivity of soil;

As NDVI > NDVI _vtime, ε=ε _v, wherein NDVI _vthe NDVI of pure vegetation area, ε _vit is the emissivity of vegetation;

Work as NDVI _s≤ NDVI≤NDVI _vtime, ε=ε _s(1-FVC)+ε _vfVC

FVC is vegetation coverage:

$\mathrm{FVC}={\left[\frac{\mathrm{NDVI}-\mathrm{NDV}{I}_s}{{\mathrm{NDVI}}_v-{\mathrm{NDVI}}_s}\right]}^2$

NDVI _sand NDVI _vthe exposed soil region of homogeneous can be chosen and vegetation area obtains from image; ε _sand ε _vcalculated by MODIS UCSB emissivity storehouse and Landsat 8 TIRS spectral response function;

3rd step: calculate 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}$

Wherein, τ _ifor the atmospheric transmittance of i wave band, τ _jfor the atmospheric transmittance of j wave band, ε _ifor the emissivity of i wave band, ε _jfor the emissivity of j wave band, k represents a kth pixel, T _{i, k}for brightness temperature on the star of a kth pixel i wave band, T _{j, k}for brightness temperature on the star of a kth pixel j wave band, for brightness temperature on the average star of N number of pixel i wave band, for brightness temperature on the average star of N number of pixel j wave band, for Landsat8 data, i, j are respectively 10,11, N and represent window size, get 20 pixel * 20 pixels;

Coefficient a and b utilizes MODTRAN4.0 atmospheric radiation transmission and TIGR database to carry out simulated atmosphere vapour content w and Landsat8 Thermal infrared bands atmospheric transmittance ratio τ ₁₁/ τ ₁₀between relation obtain:

w＝-18.973(τ ₁₁/τ ₁₀)+19.13 R ²＝0.9663，τ ₁₁/τ ₁₀＞0.9

w＝-13.412(τ ₁₁/τ ₁₀)+14.158 R ²＝0.9366，τ ₁₁/τ ₁₀＜0.9

4th step: calculate surface temperature

$T_s=\mathrm{\&gamma;}[\frac{1}{\mathrm{\&epsiv;}}({\mathrm{\&psi;}}_1{L}_{\mathrm{sen}}+{\mathrm{\&psi;}}_2)+{\mathrm{\&psi;}}_3]+\mathrm{\&delta;}$

Wherein T _sbe surface temperature, ε is emissivity, L _senbe Landsat8 the 10th wave band star on radiance, (γ, δ) can be expressed as:

$\mathrm{\&gamma;}\&ap;\frac{T_{\mathrm{sen}}^2}{b_{\mathrm{\&gamma;}}{L}_{\mathrm{sen}}},\mathrm{\&delta;}\&ap;T_{\mathrm{sen}}-\frac{T_{\mathrm{sen}}^2}{b_{\mathrm{\&gamma;}}}$

Wherein T _senbe Landsat8 the 10th wave band star on brightness temperature, b _γequal 1324K, ψ ₁, ψ ₂, and ψ ₃be air function, following formula can be utilized to carry out approximate obtaining from atmospheric water vapor content (w):

ψ ₁＝0.04019w ²+0.02916w+1.01523

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

ψ ₃＝0.00918w ²+1.36072w-0.27514

Accompanying drawing explanation

The relation of Fig. 1 Landsat 8 Thermal infrared bands atmospheric transmittance ratio and atmospheric water vapor content

Embodiment

The present invention utilizes single channel Surface Temperature Retrieval method to carry out inverting surface temperature from Landsat 8 the 10th wave band, and single channel method obtains based on the simplification of heat wave section radiation transfer equation, can be expressed as:

$T_s=\mathrm{\&gamma;}[\frac{1}{\mathrm{\&epsiv;}}({\mathrm{\&psi;}}_1{L}_{\mathrm{sen}}+{\mathrm{\&psi;}}_2)+{\mathrm{\&psi;}}_3]+\mathrm{\&delta;}$

Wherein T _sbe surface temperature, ε is emissivity, L _senbe radiance on star, (γ, δ) can be expressed as:

$\mathrm{\&gamma;}\&ap;\frac{T_{\mathrm{sen}}^2}{b_{\mathrm{\&gamma;}}{L}_{\mathrm{sen}}},\mathrm{\&delta;}\&ap;T_{\mathrm{sen}}-\frac{T_{\mathrm{sen}}^2}{b_{\mathrm{\&gamma;}}}$

Wherein T _senbrightness temperature on star, b _γequal 1324K, ψ ₁, ψ ₂, and ψ ₃be air function, following formula can be utilized to be similar to from atmospheric water vapor content (w) and obtain (Jim é nez- j.C., Sobrino, J.A., 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 _LQ _cal+A _L

M _lfor the gain of wave band, A _lfor being biased of wave band, M _land A _lobtain from Landsat 8 header file, Q _calfor image DN value.

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

K1 and K2 is constant, obtains from Landsat 8 header file.

Emissivity utilizes NDVI (Normalized Difference Vegetation Index) threshold method to obtain:

$\mathrm{NDVI}=\frac{{\mathrm{DN}}_{\mathrm{band}5}-D{N}_{\mathrm{band}4}}{{\mathrm{DN}}_{\mathrm{band}5}+{\mathrm{DN}}_{\mathrm{band}4}}$

Wherein DN _band5and DN _band4represent the DN value of Landsat8 the 5th wave band and the 4th wave band image respectively.

As NDVI < NDVI _stime, ε=ε _s, wherein NDVI _sthe NDVI in pure exposed soil region, ε _sit is the emissivity of soil;

As NDVI > NDVI _vtime, ε=ε _v, wherein NDVI _vthe NDVI of pure vegetation area, ε _vit is the emissivity of vegetation;

Work as NDVI _s≤ NDVI≤NDVI _vtime, ε=ε _s(1-FVC)+ε _vfVC

FVC is vegetation coverage:

$\mathrm{FVC}={\left[\frac{\mathrm{NDVI}-\mathrm{NDV}{I}_s}{{\mathrm{NDVI}}_v-{\mathrm{NDVI}}_s}\right]}^2$

NDVI _sand NDVI _vthe exposed soil region of homogeneous can be chosen and vegetation area obtains from image.ε _sand ε _vcalculated by MODIS UCSB emissivity storehouse and Landsat 8 TIRS spectral response function.

Atmospheric water vapor content (w) is based on splitting 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.) inverting is carried out, this algorithm hypothesis is under cloudless condition, (for Landsat 8 in N number of adjacent picture elements region, N can value be 20, namely window size is 20 pixel * 20 pixels), atmospheric conditions and emissivity do not change, only surface temperature changes, w is calculated as follows:

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)$

Wherein, τ _ifor the atmospheric transmittance of i wave band, τ _jfor the atmospheric transmittance of j wave band, ε _ifor the emissivity of i wave band, ε _jfor the emissivity of j wave band, k represents a kth pixel, T _{i, k}for brightness temperature on the star of a kth pixel i wave band, T _{j, k}for brightness temperature on the star of a kth pixel j wave band, for brightness temperature on the average star of N number of pixel i wave band, for brightness temperature on the average star of N number of pixel j wave band.For Landsat8 data, i, j are respectively 10, and 11.

For Landsat8 TIRS data, employing formula (1) and (2) Retrieval of Atmospheric Water Vapor content, need to determine coefficient a and b, coefficient a and b can be solved by the relation of atmospheric radiation transmission simulated atmosphere vapour content and Thermal infrared bands atmospheric transmittance ratio and obtain.

MODTRAN4.0 atmospheric radiation transmission and TIGR (Thermodynamic Initial Guess Retrieval, TIGR) database is utilized to carry out simulated atmosphere vapour content w and Landsat8 Thermal infrared bands atmospheric transmittance ratio τ ₁₁/ τ ₁₀between relation.TIGR database is a meteorogical phenomena database be made up of 2311 atmospheric profiles; Wherein every bar cross-sectional data all contains from earth's surface to the air pressure at every layer, atmospheric envelope top, temperature, moisture content and ozone content.872 tropical atmosphere sections are included, 742 mid latitude atmosphere sections and 697 high latitude atmospheric profiles in TIGR database.Contain in TIGR database one widely atmospheric water vapor content range (from 0.066 to 7.833g/cm ²).TIGR database is carried out simulated atmosphere vapour content w and Thermal infrared bands atmospheric transmittance ratio τ as the input of MODTRAN4.0 model ₁₁/ τ ₁₀between relation.Fig. 1 represents the relation of Landsat 8 Thermal infrared bands atmospheric transmittance ratio and the atmospheric water vapor content obtained based on 2311 TIGR atmospheric profiles and MODTRAN4.0 atmospheric radiation transmission.

The relation of Fig. 1 Landsat 8 Thermal infrared bands atmospheric transmittance ratio and atmospheric water vapor content

As shown in Figure 1, Landsat8 data 11 wave band and 10 wave band atmospheric transmittance ratios and atmospheric water vapor content have good correlativity.As can be seen from Figure 1, be that 0.9 place exists a flex point at transmitance ratio, in order to the relational expression better between matching atmospheric transmittance ratio and atmospheric water vapor content, carry out matching with 0.9 for atmospheric transmittance ratio is divided into two sections by separation, obtain the relational expression between atmospheric transmittance ratio and atmospheric water vapor content:

w＝-18.973(τ ₁₁/τ ₁₀)+19.13 R ²＝0.9663，τ ₁₁/τ ₁₀＞0.9 (3)

w＝-13.412(τ ₁₁/τ ₁₀)+14.158 R ²＝0.9366，τ ₁₁/τ ₁₀＜0.9 (4)

Coefficient a and b can be obtained by formula 3 and formula 4.

Claims (1)

1. a Landsat 8 satellite data Surface Temperature Retrieval method, the steps include:

The first step, calculate Landsat 8 the 10th wave band and the 11st wave band star on brightness temperature on radiance and star

L _sen＝M _LQ _cal+A _L

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

Wherein, L _senradiance on star, T _senbe brightness temperature on star, ML is the gain of wave band, A _lfor being biased of wave band, Q _calfor image DN value, K1 and K2 is constant, ML, A _land K1 and K2 obtains from Landsat 8 header file;

Second step, utilize NDVI (Normalized Difference Vegetation Index) threshold method to obtain emissivity ε:

Wherein DN _band5and DN _band4represent the DN value of Landsat8 the 5th wave band and the 4th wave band image respectively;

As NDVI < NDVI _stime, ε=ε _s, wherein NDVI _sthe NDVI in pure exposed soil region, ε _sit is the emissivity of soil;

As NDVI > NDVI _vtime, ε=ε _v, wherein NDVI _vthe NDVI of pure vegetation area, ε _vit is the emissivity of vegetation;

Work as NDVI _s≤ NDVI≤NDVI _vtime, ε=ε _s(1-FVC)+ε _vfVC

FVC is vegetation coverage:

NDVI _sand NDVI _vthe exposed soil region of homogeneous can be chosen and vegetation area obtains from image; ε _sand ε _vcalculated by MODIS UCSB emissivity storehouse and Landsat 8TIRS spectral response function;

3rd step: calculate atmospheric water vapor content w

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

and

Wherein, τ _ifor the atmospheric transmittance of i wave band, τ _jfor the atmospheric transmittance of j wave band, ε i is the emissivity of i wave band, ε _jfor the emissivity of j wave band, k represents a kth pixel, T _{i, k}for brightness temperature on the star of a kth pixel i wave band, T _{j, k}for brightness temperature on the star of a kth pixel j wave band, for brightness temperature on the average star of N number of pixel i wave band, for brightness temperature on the average star of N number of pixel j wave band, for Landsat8 data, i, j are respectively 10,11, N and represent window size, get 20 pixel * 20 pixels;

Coefficient a and b utilizes MODTRAN4.0 atmospheric radiation transmission and TIGR database to carry out simulated atmosphere vapour content w and Landsat8 Thermal infrared bands atmospheric transmittance ratio τ ₁₁/ τ ₁₀between relation obtain:

w＝-18.973(τ ₁₁/τ ₁₀)+19.13 R ²＝0.9663，τ ₁₁/τ ₁₀＞0.9

w＝-13.412(τ ₁₁/τ ₁₀)+14.158 R ²＝0.9366，τ ₁₁/τ ₁₀＜0.9

4th step: calculate surface temperature

Wherein T _sbe surface temperature, ε is emissivity, L _senbe Landsat8 the 10th wave band star on radiance, (γ, 6) can be expressed as:

Wherein T _senbe Landsat8 the 10th wave band star on brightness temperature, b _γequal 1324K, ψ ₁, ψ ₂, and ψ ₃be air function, following formula can be utilized to carry out approximate obtaining from atmospheric water vapor content (w):

ψ ₁＝0.04019w ²+0.02916w+1.01523

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

ψ ₃＝0.00918w ²+1.36072w-0.27514 。