在GEE中实现Landsat C2L2 大气校正法地表温度反演方法对比

一、大气校正法

1.地表温度反演公式

(1)得到B(Ts)​同温度下的黑体辐射亮度：

Lλ：卫星传感器接收到的辐射亮度（[W·m−2·sr−1·μm−1]）

B(Ts)​：温度Ts对应的黑体辐射亮度（[W·m−2·sr−1·μm−1]）。B是普朗克函数，TS是地表温度LST。

ε：地表比辐射率（0-1，无量纲）

τ：大气透射率（0-1，无量纲），可以用大气水分含量来估计。

Ld：大气下行辐射（[W·m−2·sr−1·μm−1]）

Lu：大气上行辐射（[W·m−2·sr−1·μm−1]）

(2)再由普朗克公式，可求得真实地表温度：

对于TM，K1 =607.76 W/(m2*µm*sr)，K2 =1260.56K。
对于ETM+，K1=666.09 W/(m2*µm*sr)，K2 =1282.71K。
对于TIRS Band10，K1= 774.89 W/(m2*µm*sr)，K2 = 1321.08K。

因此，只要能知道地表比辐射率，就可求得地表温度。

2.地表比辐射率计算

（1）Sobrino等人提出的NDVI阈值法：（简化公式）

Pv为植被覆盖度。

也可以不用简化公式，用NDVI0.2、0.5分段求取，结果会更精确一些。具体见文章《Land surface temperature retrieval from LANDSAT TM 5》，

（二）覃志豪等人提出的针对不同地表的经验公式

覃志豪于2004年发表在《国土资源遥感》的《陆地卫星TM6波段范围内地表比辐射率的估计》的方法：

根据文章提出的经验公式，根据NDVI阈值计算水体、自然表面、城镇三种不同地表的比辐射率：

这也是简化算法。

二、算法实现

（1）数据及波段选取：

使用landsat8 collection2-level2数据集反演地表温度，用到的波段如图：

SR_B4：红色波段 (Red)

SR_B5：近红外波段 (NIR)

ST_TRAD：热辐射亮度

ST_URAD：大气上行辐射

ST_DRAD：大气下行辐射

ST_ATRAN：大气透过率

代码思路可参考：

使用GEE实现温度反演的三种方法_gee反演地表温度-CSDN博客

Landsat 计算LST（地表温度）——没有大气剖面参数计算器怎么办_大气校正参数计算器用不了-CSDN博客

（2）算法设计

根据第一节地表温度计算原理，分为4步进行：

植被覆盖度Pv——地表比辐射率ε——亮温BT——地表温度Ts：

第1步：利用red、NIR波段计算NDVI，用NDVI反演出植被覆盖度Pv:

var Pv = ndvi.expression('(b1 <= 0.05) ? 0 : ((b1 >= 0.7) ? 1 : (b1 - 0.05) / 0.65)', {b1: ndvi});

第2步：利用覃志豪的经验公式，反演地表比辐射率ε：

var LSE = ndvi.expression( 'b1 < 0 ? 0.995 : (b1 > 0.7 ? (0.9625 + 0.0614*Pv - 0.0461*Pv*Pv) : (0.9589 + 0.0860*Pv - 0.0671*Pv*Pv ))', {'Pv': Pv, 'b1': ndvi} );

第3步：利用反演出的比辐射率ε，以及4个波段：ST_TRAD（热辐射亮度）、ST_URAD（大气上行辐射）、ST_DRAD（大气下行辐射）、ST_ATRAN（大气透过率），计算B(Ts)​同温度下的黑体辐射亮度，采用前文的公式1：

式中，Lλ=ST_TRAD ，Lu=ST_URAD，Ld=ST_DRAD ，t = ST_ATRAN。将4波段和ε代入。

第4步：由前文公式2，B(Ts)利用普朗克公式反演真实地表温度，其中，K1、 K2根据传感器选取。

三、实验结果

使用landsat8 collection2-level2级别数据，相同时间段实验结果

1.ST B10：

仅需将开尔文换算为摄氏度

2.大气校正法：

3.MOD11A2

顺便对比了下MODIS的结果：

1.ST B10波段由于官方使用Aster 比辐射率数据进行反演，因比辐射率数据本身的缺失，导致地表温度局部有缺失。

2.landsat大气校正法反演结果总体与ST B10官方反演结果相差不多，但不会出现缺失情况（因去云操作导致的缺失值除外）。

3.MODIS地表温度数据比landsat大气校正结果低1°C左右，对比了其他年份和其他研究区，也有低于2°C的，但基本都是MODIS地表温度数据较低。研究区南部也存在明显差异

4 如果不想自行反演，也可直接使用官方collection2-level2 ST波段地表温度数据，用MODIS数据填补缺失值。结果如下：

但是，如果研究区ST波段有大范围缺失，用MODIS填补，因二者数据值本身的差异，填补区域会与周围区域有明显边界。