GEE学习篇：Landsat地表温度反演的方法

最近在学习使用GEE平台反演地表温度，遇到了一些困惑也学到了很多方法，与大家共勉~

一、原理

地表温度反演：主要是基于卫星的热红外传感器观测到的地表热辐射，之后将其中的大气影响减去从而得到地表热辐射强度，最后经过热辐射强度转换得到地表温度。

卫星传感器所观测到的辐射亮度值主要包括三个部分：地表辐射经大气到达卫星传感器的能量；大气向上辐射亮度以及大气向下辐射并由地表反射后的能量。

二、汇总一些使用GEE进行地表温度反演的方法

1、GEE平台中的数据集LANDSAT/LC08/C02/T1_L2及LANDSAT/LT05/C02/T1_L2数据集中的ST_B10波段即为处理好的LST波段（This dataset contains atmospherically corrected surface reflectance and land surface temperature derived from the data produced by the Landsat 8 OLI/TIRS sensors.），可以使用公式（" B1-273.15 "）转变为地表温度后，直接进行调用。

https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C02_T1_L2

https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LT05_C02_T1_L2

GEE代码：

//现实矢量图层
var styling = {color:'red',fillColor:'00000000'}//设置显示样式
Map.centerObject(table,10)
Map.addLayer(table2.style(styling),{},'table')//可视化地理数据
var roi=table

//定义数据筛选时间
var star_date = '2014-01-1'//定义起始时间
var end_date  = '2014-12-31'//定义终止时间

//GEE平台对该数据集进行了缩放偏移，使用前需要进行预处理
function maskL8sr(image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
  return image.addBands(opticalBands, null, true)
              .addBands(thermalBands, null, true);
}

//影像筛选
var dataset = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
    .filterBounds(roi)
    .filterDate(star_date, end_date)//时间过滤
    .filter(ee.Filter.lt("CLOUD_COVER",20))    //云量小于10%
    .map(maskL8sr)
    .median()//计算均值    
var dataset1 = dataset.clip(roi)
//地表温度计算
var img = dataset1.select("ST_B10")
var lst = img.expression(
'B1-273.15',
{
B1:img.select('ST_B10'),
});

print("LST处理后直方图", ui.Chart.image.histogram(lst, table, 100, 258))
print(lst)

Map.addLayer(dataset1, {bands: ['SR_B4', 'SR_B3', 'SR_B2'], min: 0, max: 0.3});
Map.addLayer(lst.clip(roi), {}, "LST");

但是该方法有两个缺点：

一是数据可能会存在少量的异常值；二是不知道什么原因导致该数据集中的ST_B10波段有残缺（影像中的其它波段没有这个问题，只有B10波段存在这个问题），在使用前要先检查自己的研究区有没有这个问题（目前我知道的是广东省的部分地区存在这个问题），如果B10波段存在残缺就不能使用这个方法了，可以选择其它数据集。

Landsat/LT05/C01/T1_SR和Landsat/LC08/C01/T1_SR数据集虽然已经弃用（仅限于2021年之前）但是B10波段不会出现这个问题

2.大气校正法（也称为辐射传输方程：Radiative Transfer Equation——RTE）

3.单通道算法

4.单窗算法

可以参考这篇博客：基于Google Earth Engine的Landsat单窗算法地表温度（LST）反演_前面几篇博客介绍了基于landsat这一多光谱遥感图像数据的多种地表温度(lst)反演方-优快云博客

5.地表绝对温度算法（目前再GEE、PIE、Python等代码中使用最多的计算地表温度方法）

Malaret等人的地表绝对温度算法即把DN值转变为辐射温度（也可以称为亮度温度），再应用Artis等人提出的绝对表面温度表达式，公式如下：

Ts=Tb/(1+(λ*Tb/ρ)lnε)

式中:Tb为辐射温度；热红外波段的中心波长λ=11.5 μm；ρ=h×c/σ (1.439×10^-2mK) , 其中, 光速c=2.998×10^8 m/s;普朗克常数h=6.626×10^-34Js;波耳兹曼常数σ=1.38×10^-23J/K；ε为地表比辐射率。已知公式中只有ε为未知量

采用NDVI阈值法计算地表比辐射率，主要步骤如下：

1）计算NDVI

2）计算植被覆盖度PV

3）计算地表比辐射率：x=0.004PV+0.986

GEE代码：

//地表温度反演
var thermal = images.select('B10').multiply(0.1)
var Emissivity = function(image, nir, red) {
  var ndvi_e = "(nir-red)/(nir+red)"
  var ndvi =  ee.Image(0).expression(ndvi_e, {'nir': image.select(nir), 'red': image.select(red)}).rename('ndvi')
  var dic = ndvi.reduceRegion(ee.Reducer.percentile([5, 95]), roi,1000);
  var ndvi_min  = ee.Number(ee.Dictionary(dic).get("ndvi_p5"));
  var ndvi_max  = ee.Number(ee.Dictionary(dic).get("ndvi_p95"));
  var fvc = (ndvi.subtract(ndvi_min)).divide(ndvi_max.subtract(ndvi_min));
var e_building = ndvi.expression(
    '0.9589+0.0860*F-0.0671*(F**2)', {F: fvc});
var e_surface = ndvi.expression(
    '0.9625+0.0614*F-0.0461*(F**2)',{F: fvc});
var e_water = 0.995;
//针对不同的地物，计算地表比辐射率
var exp = e_building.multiply(ndvi.gt(ndvi_min).and(ndvi.lt(ndvi_max))).add(e_surface.multiply(ndvi.gte(ndvi_max))).add(ndvi.lte(ndvi_min).multiply(e_water));
  return exp.rename('EM')}
var EMM1 = Emissivity(images, 'B5', 'B4')
var EMM = EMM1.clip(roi);
var LST = thermal.expression(
    '(Tb/(1 + (0.0010904 * (Tb / 1.438))*log(Ep)))-273.15', {
      'Tb': thermal.select('B10'),
      'Ep': EMM.select('EM'),
});
print("LST处理后直方图", ui.Chart.image.histogram(LST, roi, 100, 258))

如有问题换用大家指出~

参考文章：

【硬核】Landsat8地表温度反演 - 知乎

利用Landsat8 TIRS反演地表温度实例

GEE-计算Landsat8地表温度LST - 知乎

PIE-Engine 基于影像的通州区地表温度反演代码解析_m0_56819632的博客-优快云博客

张建平,王崇倡.利用TM图像反演翁牛特旗地表温度[J].水资源与水工程学报,2009,20(04):62-66.