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

热红外遥感（Infrared Remote Sensing）是指传感器工作波段限于红外波段范围之内的遥感。即利用星载或机载传感器收集、记录地物的热红外信息，并利用这种热红外信息来识别地物和反演地表参数如温度、湿度和热惯量等。

目前有很多的卫星携带了热红外传感器，包括ASTER、AVHRR、MODIS、TM/ETM+/ TIRS等。

为了后面阅读方法，表15.5列举了常见的几个名词解释。

表15.5热红外遥感中常见名词解释

名词	说明
辐射出射度	单位时间内从单位面积上辐射的辐射能量称为辐射出射度，单位一般是W/m2 (瓦特/平方米)
辐射亮度（Radiance）	辐射源在某一方向上单位投影表面、单位立体角内的辐射通量，称为辐射亮度，单位一般是W/(m2*µm*sr)（瓦特/平方米.微米.球面度）
比辐射率( Emissivity)	也称发射率，物体的辐射出射度与同温度黑体辐射出射度的比值。如果物体指的是地表，称为地表比辐射率。
大气透射率	通过大气(或某气层)后的辐射强度与入射前辐射强度之比。
亮度温度 （Brightness Temperature）	当一个物体的辐射亮度与某一黑体的辐射亮度相等时，该黑体的物理温度就被称之为该物体的“亮度温度”，所以亮度温度具有温度的量纲，但是不具有温度的物理含义，它是一个物体辐射亮度的代表名词。

15.4.1 地表温度反演模型概述

目前，地表温度反演算法主要有以下三种：大气校正法（也称为辐射传输方程：Radiative Transfer Equation——RTE）、单通道算法和分裂窗算法。

1. 大气校正法

基本原理：首先估计大气对地表热辐射的影响, 然后把这部分大气影响从卫星传感器所观测到的热辐射总量中减去, 从而得到地表热辐射强度, 再把这一热辐射强度转化为相应的地表温度。

具体实现为：卫星传感器接收到的热红外辐射亮度值L_λ由三部分组成：大气向上辐射亮度L↑，地面的真实辐射亮度经过大气层之后到达卫星传感器的能量；大气向下辐射到达地面后反射的能量。卫星传感器接收到的热红外辐射亮度值L_λ的表达式可写为（辐射传输方程）：

          L_λ = [εB(T_S) + (1-ε)L↓]τ + L↑                           （15.1）

式中，ε为地表比辐射率，T_S为地表真实温度(K)，B(T_S)为黑体热辐射亮度，τ为大气在热红外波段的透过率。则温度为T的黑体在热红外波段的辐射亮度B(T_S)为：

         B(T_S) = [L_λ - L↑- τ(1-ε)L↓]/τε                            （15.2）

T_s可以用普朗克公式的函数获取。

            T_S = K₂/ln(K₁/ B(T_S)+ 1)                              （15.3）

对于TM，K₁ =607.76 W/(m²*µm*sr)，K₂ =1260.56K。

对于ETM+，K₁=666.09 W/(m²*µm*sr)，K₂ =1282.71K。

对于TIRS Band10，K₁= 774.89 W/(m²*µm*sr)，K₂ = 1321.08K。

从上可知此类算法需要2个参数：大气剖面参数和地表比辐射率。大气剖面参数 在NASA提供的网站(http://atmcorr.gsfc.nasa.gov/)中，输入成影时间以及中心经纬度可以获取大气剖面参数。适用于只有一个热红外波段的数据，如Landsat TM /ETM+/TIRS数据。

2. 单窗算法

    单窗算法(Mono- window Algorithm) 是覃志豪（2004）等根据地表热辐射传导方程, 推导出的一个利用Landsat TM /ETM+第六波段数据反演地表温度的算法，该算法计算公式如下：

              T_S =[a(1- C- D)+(b(1- C- D)+C+D) T₆- DT_a]/C             （15.4）

式中，T_S是为地表真实温度(K) ; a 和 b 是常量，分别为-67.355351 和 0.458606； C 和 D 是中间变量， C =ετ，D=(1-τ) ([1+(1-ε) τ], 其中，ε是地表比辐射率, τ是大气透射率；T₆ 是卫星高度上传感器所探测到的像元亮度温度(K) ；大气平均作用温度(T_a)与地面附近(一般为2 m处)气温(T0)存在如下线性关系( Ta与T₀的单位为K) ：

• 热带平均大气(北纬 15°，年平均)

Ta= 17.9769+0.91715 T₀

• 中纬度夏季平均大气(北纬 45°，7 月)

Ta= 16.0110+0.92621 T₀

• 中纬度冬季平均大气(北纬 45°，1 月)

Ta= 19.2704+0.91118 T₀

从上可知此类算法需要3个参数：大气平均作用温度、大气透过率和地表比辐射率。

3. 分裂窗算法

分裂窗算法（也叫劈窗算法）最初是为反演海面温度开发的，具体地说是针对NOAA/AVHRR的4和5通道设计的，后来也被用来反演地表温度，这种算法较成熟，精度也高。分裂窗算法以地表辐射传导方程为基础，利用10~13μm 大气窗口内，两个相邻热红外通道(一般为10.5～11.5μm、11.5～12.5μm)对大气吸收作用的不同，通过两个通道测量值的各种组合来剔除大气的影响，进行大气和地表比辐射率的修正。表达式为：

                T _S= T ₄+ A (T ₄- T ₅) + B                              （15.4）

式中，T _S为地表真实温度，T ₄和 T ₅分别为AVHRR的4和5通道，A和B为常量。

AVHRR 的通道4(10.15~11.13μm) 和通道5 (11.15~12.15μm) 恰与MODIS 的第31 波段 (10.178~ 11.128μm) 和32 波段 (11.177~ 12.127μm ) 的中心波长相对应, 可将MODIS 的31、32 波段用于分裂窗算法进行地表温度计算。

很多学者对这个算法进行了推演，得到很多新的算法，如覃志豪、毛克彪等人。如下为覃志豪（2005）等人改进的分裂窗法模型用于MODIS数据，算法为：

T _s = A ₀ + A ₁T ₃₁ - A ₂T ₃₂                                 （15.5）

其中：T_s是地表温度( K ) , T ₃₁和T₃₂分别是MODIS第31和32波段的亮度温度；A₀、A₁和 A ₂是分裂窗算法的参数，分别定义如下：

A₀=[D₃₂(1-C₃₁-D₃₁)/(D₃₂C₃₁-D₃₁C₃₂) ] a₃₁-[D₃₁(1-C₃₂-D₃₂)/(D₃₂C₃₁-D₃₁C₃₂)]a₃₂  （15.6）

A ₁=1+D₃₁/(D ₃₂C₃₁-D ₃₁C₃₂)+[D ₃₂(1-C₃₁-D₃₁)/(D₃₂C₃₁-D₃₁C₃₂)]b₃₁            （15.7）

A₂=D₃₁/(D₃₂C₃₁ -D₃₁C₃₂)+[D₃₁(1-C₃₂-D ₃₂)/(D₃₂C₃₁-D₃₁C₃₂)]b₃₂               （15.8）

式中，a ₃₁，b₃₁，a ₃₂ 和 b ₃₂ 是常量, 根据 MODIS的波段特征确定, 在地表温度 0 ~ 5 e 范围内，这些常量分别可取 a ₃₁ = -64.60363，b₃₁ = 0.440817，a ₃₂ = -68.72575，b₃₂ = 0.473453。

上述公式的中间参数分别计算如下：

C_i = Ɛ _iτ_i (ɵ)                                           （15.9）

D _i= [1-τ_i (ɵ)][1+(1-Ɛ _i)τ_i (ɵ)]                              （15.10）

式中，i是指 MODIS的第31和32波段，分别为 i= 31或 32；τ_i (ɵ)是视角为ɵ的大气透过率；Ɛ _i 是波段i的地表比辐射率。

T ₃₁和T₃₂的亮度温度可用普朗克函数计算：

T_i= K_i2 /ln(1+K _i1/I_i)                                    （15.11）

式中，K _i1和K _i2是常量，对于第i =31波段，分别为 K _31,1 = 729.541636W/(m²*µm* sr)，K_{31, 2} =1304.413871K；对于第 i = 32 波段，为 K _{32 ,1} = 474.6 84780 W/(m²*µm*sr)，K _{31 , 2} =1196 . 978785 K；I_i是31或者32波段的辐射亮度值。

上述三种方法，都需要两个参数：大气透过率和地表比辐射率。

4. 大气透过率

一般可以通过大气水汽含量来计算大气透过率，对于Landsat可以在NASA公布的网站查询（http://atmcorr.gsfc.nasa.gov）。

5. 地表比辐射率

比较常用的是其中一种方法是先对遥感影像进行分类, 将地表分为不同的覆盖类型, 再根据实测或者经验值的地物比辐射率给各个地表覆盖类型赋予不同的值, 从而生成地表比辐射率影像。目前一些比辐射率数据库，如MODIS UCSB 比辐射率库（http://www.icess.ucsb.edu/modis/EMIS/html/em.html）等。

另外还有利用归一化植被指数（NDVI）计算地表比辐射率，如NDVI的对数与地表比辐射率存在线性相关性；利用NDVI的阈值对地表进行分类，然后给各个地表覆盖类型赋予不同的值。

15.4.2 反演流程

本实例是基于大气校正法，利用Landsat8 TIRS反演地表温度。

注：可直接使用ENVI扩展工具：Landsat 8 LST 进行地表温度反演，本文为了让大家更好的理解整个反演过程。

主要内容就是使用BandMath工具计算公式（15.2）和公式（15.3），处理流程如图15.8所示。

图15.8 基于大气校正法的TIRS反演流程图

15.4.3 图像辐射定标和大气校正

     在主界面中，选择File→Open，在文件选择对话框中选择“txt”文件，ENVI自动按照波长分为五个数据集：多光谱数据（1-7波段），全色波段数据（8波段），卷云波段数据（9波段），热红外数据（10，11波段）和质量波段数据（12波段）。

（1）在Toolbox工具箱中，选择Radiometric Correction/Radiometric Calibration。在File Selection对话框中，选择数据LC81230322013132LGN02_MTL_Thermal，单击Spectral Subset选择Thermal Infrared1（9），打开Radiometric Calibration面板。

（2）在Radiometric Calibration面板中，设置以下参数：

（3）定标类型（Calibration Type）：辐射亮度值（radiance）。

（4）其他选择默认参数。

（5）选择输出路径和文件名，单击OK按钮执行定标处理。得到Band10辐射亮度图像。

 

15.4.4 地表比辐射率计算

TIRS的Band10热红外波段与TM/ETM+6热红外波段具有近似的波谱范围，本例采用TM/ETM+6相同的地表比辐射率计算方法。使用Sobrino提出的NDVI阈值法计算地表比辐射率。

                    ε=0.004Pv+0.986                                   （15.12）

其中， Pv是植被覆盖度，用以下公式计算：

        Pv = [(NDVI- NDVI_Soil)／(NDVI_Veg - NDVI_Soil)]                      （15.13）

其中，NDVI为归一化植被指数，NDVI_Soil为完全是裸土或无植被覆盖区域的NDVI值，NDVI_Veg则代表完全被植被所覆盖的像元的NDVI值，即纯植被像元的NDVI值。取经验值NDVI_Veg = 0.70和NDVI_Soil = 0.05，即当某个像元的NDVI大于0.70时，Pv取值为1；当NDVI小于0.05，Pv取值为0。

（1）在Toolbox工具箱中，双击Spectral/Vegetation/NDVI工具，在文件输入对话框中，选择Landsat8 OLI大气校正结果。

    提示：覃志豪（2004）提出使用原始的DN值图像计算NDVI对反演结果影响不大。

（2）在NDVI Calculaton parameters对话框中，选择NDVI计算波段：Red：4，Near IR：5。

（3）选择输出文件名和路径。

（4）在Toobox中，选择Band Ratio/Band Math，输入表达式：

(b1 gt 0.7)*1+(b1 lt 0.05)*0+(b1 ge 0.05 and b1 le 0.7)*((b1-0.05)/(0.7-0.05))

      其中，b1：NDVI。

计算得到植被覆盖度图像。

（5）在Toobox中，选择Band Ratio/Band Math，输入表达式：

0.004*b1+0.986

      其中，b1：植被覆盖度图像。

     计算得到地表比辐射率图像。

提示：为了得到更精确的地表比辐射率数据，可使用覃志豪（2004）等提出的先将地表分成水体、自然表面和城镇区，分别针对三种地表类型计算地表比辐射率：

• 水体像元比辐射率：995
• 自然表面像元比辐射率：ε_surface = 0.9625 + 0.0614Pv - 0.0461Pv²
• 城镇区像元比辐射率：ε_building = 0.9589 + 0.086Pv - 0.0671Pv²

15.4.5 黑体辐射亮度与地表温度计算

在NASA公布的网站查询（http://atmcorr.gsfc.nasa.gov），输入成影时间：2013-10-03  02:55和中心经纬度（Lat：40.32899857，Lon：116.70610046），以及其他相应的参数，得到大气剖面信息为：

• 大气在热红外波段的透过率τ：90
• 大气向上辐射亮度L↑：75 W/(m²·sr·μm)
• 大气向下辐射亮辐射亮度L↓：29W/(m²·sr·μm)

    提示：由于缺少地表相关参数（气压、温度、相对湿度等信息），得到的结果是基于模型计算的结果。

（1）依据公式（2），在Toolbox工具箱中，双击Band Ratio/Band Math工具，输入表达式：

(b2-0.75-0.9*(1-b1)*1.29)/(0.9*b1)

       其中，b1：地表比辐射率图像

             b2：Band10辐射亮度图像

   计算得到同温度下的黑体辐射亮度图像。

（2）依据公式（3），在Toobox中，双击Band Ratio/Band Math工具，输入表达式：

            (1321.08)/alog(774.89/b1+1)-273

      其中，b1：同温度下的黑体辐射亮度图像

得到单位为摄氏度的地表温度图像。

提示：公式（15.3中），TIRS Band10的K₁和K₂ 是从*_MTL.txt元数据文件中获取。

（3）在图层管理器（Layer Manager）中的地表温度图像图层，右键选择 Raster Color Slices。将温度划分为四个区间：

• 25℃以上
• 22℃至25℃
• 20℃至22℃
• 15℃至20℃
• 低于15℃

（4）分别浏览几个温度区间的空间分布范围。

（5）统计反演结果得出81%区域的温度集中在15~22°区间。

提示：缺少同步温度测量数据用于验证反演结果，查询2013年10月3号北京市最低气温为10°，最高气温为22°。本示例反演结果大部分在这个区间内，反演结果有一定的参考价值。