Landsat 计算LST（地表温度）——没有大气剖面参数计算器怎么办

最新推荐文章于 2024-03-22 09:00:00 发布

承伤成殇

最新推荐文章于 2024-03-22 09:00:00 发布

阅读量1.5w

点赞数 107

分类专栏：遥感算法文章标签：经验分享算法

本文链接：https://blog.csdn.net/weixin_43953542/article/details/135849420

版权

遥感算法专栏收录该内容

1 篇文章 11 订阅

订阅专栏

本文探讨了在AtmosphericCorrectionParameterCalculator网站停止服务后，如何使用USGS的LandsatCollection-2SurfaceTemperature数据来反演地表温度，包括使用Level-2SurfaceTemperatureBands数据、数据转换、精度验证以及替代计算方法，如利用ASTERGED数据和VFC计算地表比辐射率。

摘要由CSDN通过智能技术生成

Landsat 计算LST（地表温度）——没有大气剖面参数计算器怎么办

背景

地表温度反演的主要方法有大气校正法（辐射传输方程）、单窗算法和劈窗算法（分裂窗）三种，而这三种方法都至少需要知道大气剖面参数（主要是大气透过率）和地表比辐射率。最常用的大气剖面参数获取来源于Atmospheric Correction Parameter Calculator (nasa.gov)，通过几个参数即可得到大气剖面参数（大气上行辐射、大气下行辐射和波段大气平均透过率），但2024年1月开始上述网站停止服务，不能再获取到大气剖面参数，这对地表温度的反演有很大影响。

解决方案

根据网站描述，USGS提供新的USGS Landsat Collection-2 Surface Temperature 产品，按条件选择后获取数据集。根据网站提供的描述其提供的文件信息，选择Level-2 Surface Temperature Bands数据集。

在这里插入图片描述

根据Landsat Collection 2 表面温度 |美国地质调查局 (usgs.gov)网站提供的https://www.usgs.gov/media/files/landsat-8-collection-2-level-2-science-product-guide数据使用手册。根据命名可以找到ProductID_ST_B10的单波段.tif，根据数据使用手册分析，该文件为USGS已经反演好的地表温度产品。

根据数据规范，该文件单位为开尔文（开式温度），一般我们常用的是摄氏度（℃），所以需要转换。
在这里插入图片描述

由于USGS存储数据都是整形，所以与真实结果往往存在一定的缩放系数和偏移量，通过上表可知Multiplicative Scale Factor = 0.00341802；Additive Offset=149；在ENVI中使用Band Math进行计算b1*0.00341802+149，即可得到真实地表温度，但要注意此时的单位仍然是开尔文，所以需要在减去273.15才能得到摄氏度。

在这里插入图片描述

该数据为2022年7月某地区的结果

测试分析

恰巧在大气剖面参数网站停止服务前，曾经处理过某地区的LST，所以用这种方法和历史的LST产品进行对比，查看一下反演水平。

在这里插入图片描述

从图中可以看出，两者趋势基本一致，均值也比较接近，说明该方案可行。

同时根据Modis数据反演的某区域LST和Landsat8的结果做对比，发现该方案的结果通常比Modis的结果高2-5℃，而使用大气剖面参数的Landsat8反演结果和Modis结果非常接近；
在这里插入图片描述

在这里插入图片描述

从结果来看偏差还是存在的，但整体还是在一个比较合理的区间。下面就来探究一下USGS怎么生成ST_B10文件的。

原理探究

根据USGS的描述，他们使用了如下数据集来生产

Landsat 4-9 Surface Temperature products are generated from Landsat Collection 2 Level-1 thermal infrared bands, Top of Atmosphere (TOA) reflectance, TOA Brightness Temperature, Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Emissivity Database (GED) data, ASTER Normalized Difference Vegetation Index (NDVI) data, and atmospheric profiles of geopotential height, specific humidity, and air temperature extracted from the following inputs:

January 1, 2024 to present: Goddard Earth Observing System (GEOS) for Instrument Teams (IT)
January 1, 2000 to December 31, 2023: Forward Processing for Instrument Teams (FP-IT)
August 22, 1982 to December 31, 1999: Modern Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2)

反演LST的论文为《Atmospheric Compensation for a Landsat Land Surface Temperature Product》根据该论文描述，其采用大气校正法来计算地表温度。

在这里插入图片描述

上图很好的说明了大气校正法的原理，对于卫星Sensor接收到的热红外辐射量，总共有三个来源：

A（upwelled radiance）：来源于目标（Target）物体上方的的大气，属于大气自发向上辐射的能量，通常被称为上行辐射或上涌辐射；
B（downwelled radiance）：来源于目标上方的大气先对目标（Target）辐射，然后经过目标的反射到达传感器，通常被称为下行辐射；
C（radiance due to temperature）：是目标本身的辐射，也是我们最后需要确定的结果。
D（background）：属于背景辐射，背景辐射到目标然后反射到传感器，但根据参考论文的描述，当背景物体较少或只遮挡一小部分天空时，这种能量路径产生的光子通常可以忽略不计；

基于上述情况，可以认为传感器接收到的能量减去上行和下行，就是目标的能量，基于这个观点，提出了大气校正法。
$L_{obs} = (L_T\epsilon+(1-\epsilon)L_d)\tau + L_u$

其中： $L_{obs}$ 表示卫星传感器接收到的辐射亮度， $L_T$ 表示黑体亮度， $\epsilon$ 表示比辐射率， $L_d$ 表示大气下行辐射， $\tau$ 表示大气透过率， $L_u$ 表示大气上行辐射。

显然，根据公式（1），可以推出 $L_T$ 同温度下的黑体辐射亮度为：
$L_T = \frac{L_{obs}-L_u-\tau(1-\epsilon)L_d}{\tau\epsilon}$

由普朗克公式，可以计算得到目标的温度 $T_s$
$T_s = \frac{K_2}{\ln(\frac{K_1}{L_T}+1)}$

对于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。

这个信息可以在数据的元数据*_MTL.txt中读取的。
在这里插入图片描述

此时的得到的结果单位为开尔文，需要在减去273.15才能得到摄氏温度

原理实操

如果理论正确的话，那么根据理论和完备的数据可以推导出结果。根据USGS提供的C2数据集中的过程文件，尝试进行手动处理。

先说明一下每个文件的作用，详细可见数据使用手册的附录A:

产品名称	产品描述	含义
ProductID_ST_QA	Surface Temperature Uncertainty	QA质量波段主要反映的是产品的质量，说明数据的精确度
ProductID_ST_TRAD	Thermal Radiance	表示传感器接收到的热红外辐射亮度， $L_{obs}$
ProductID_ST_URAD	Upwelled Radiance	表示大气上行辐射， $L_u$
ProductID_ST_DRAD	Downwelled Radiance	表示大气下行辐射， $L_d$
ProductID_ST_ATRAN	Atmospheric Transmittance	表示大气透过率， $\tau$
ProductID_ST_EMIS	Emissivity estimated from ASTER GED	表示表面发射率，也意味着比辐射率*， $\epsilon$
ProductID_ST_EMSD	Emissivity standard deviation	表示表面发射率的预期偏差，
ProductID_ST_CDIST	Pixel distance to cloud	表示每个像元到最近的云层的距离，单位是 $km$
ProductID_ST_B10	Band 10 ST	表示利用Landsat8的第10波段反演的表面温度，单位开尔文

***注：**发射率即在已知温度下材料的辐射与理想黑体的发射之间的比率，也称为比辐射率。

根据缩放比例计算真实结果

由于数据存在比例缩放系数和偏移量，所以首先需要分别处理以获得各波段的真实结果，而且由于各波段的缩放系数不尽相同，所以要先进行处理，不能等到最后在处理。

按照数据使用手册中的波段文件说明，对所用到的波段进行波段运算。**注意：**表面发射率结果EMIS还有个EMSD表示预计的偏差，在这里我选择将EMIS和EMSD结果先加起来在进行0.0001的比例缩放。

计算黑体辐射

由上述公式（2）可以计算出同温度下的黑体辐射亮度，可以选择把所用的五个波段文件进行Layer Stack，然后在ENVI中进行波段运算(b1-b2-b3*(1-b4)*b5)/(b3*b4)其中b1~b5分别代表传感器接受辐射亮度（TRAD）、大气上行辐射（UPAD)、大气透射率（ATRAN）、比辐射率（表面发射率 EMIS+EMSD）、大气下行辐射（DRAD）的真实结果。

计算表面温度

由上述公式（3）代入Landat8 TIRS所使用的K1= 774.89 W/(m2*µm*sr)，K2 = 1321.08K两个参数，使用ENVI的band math进行波段运算，参考公式为(1321.08)/alog(774.89/b1+1)；b1代表黑体辐射亮度结果（单位为开尔文）。对该结果减去273.15转化为℃；

检验结果

在这里插入图片描述

可以看到这两个方法得到的结果非常接近，平均相差0.1℃，误差可能来源于浮点数计算、累计误差、QA波段影响，实际应用过程中这个偏差不会有太大的问题，此举主要说明C2数据集中ST_B10结果的来源和计算方法，验证大气校正法的实际过程。

后记

曾经最常使用的大气剖面参数计算器的上行辐射、下行辐射和平均波段大气透过率均为数据，该网站并没有标示数据来源和处理方法；一般来说有两种可能，一是使用整个影像范围的均值来代表，二是使用中心经纬度的像元值来代表，回想一下，曾经的大气剖面参数计算器是不是要输入影像的中心经纬度呢？当然这只是一种猜测，根据新的C2数据集中提供的上下行辐射和透过率数据，简单测试了一下网站计算的三个值并不等于C2数据集中对应数据的均值或者中心点值。所以这个猜测可能是无效的，我也没有任何依据。可能原来纯纯是理论模型算出来的吧。

**但是！！**神奇就神奇在了，根据C2数据集中的上下行辐射、大气透射率三个数据，我将其均值作为结果，带入到LST的计算过程（应用的是ENVI的插件），结果竟然相差不大。
在这里插入图片描述