基于Landsat 8数据单通道算法反演地表温度

夏安全+齐建国+姜振飞+马津

摘要：当前地表温度反演的遥感数据源多为Landsat TM/ETM+、MODIS数据，Landsat 8热红外数据的使用还不是很多，许多针对Landsat 8热红外数据的地表温度反演算法虽然被提出，但是否能满足不同试验区的精度要求还有待考究，同时，经验模型的使用可能会对求解地表比辐射率和大气透过率等地表温度反演参数造成不同程度的影响。因此，本研究以山东省济南市为研究区，Landsat 8数据为数据源，分别利用大气校正法、JM_SC10算法、TIRS10_SC算法，结合分类回归树（classification and regression tree，简称CART）算法与中等光谱分辨率大气透过率算法计算机模型（moderate spectral resolution atmospheric transmittance algorithm and computer model，简称MODTRAN），构建适合济南地区的温度反演参数，实现地表温度的反演，并以济南市16个气象站的温度数据为基准进行精度验证。结果表明，3种算法反演温度平均误差为1.78 ℃，TIRS10_SC算法反演精度最高，其次是大气校正法、JM_SC10算法。在一定误差要求下，3种算法均可应用于济南地区的地表温度反演。

关键词：Landsat 8；地表温度反演；单通道算法；大气透过率；地表比辐射率

中图分类号： S127 文献标志码： A 文章编号：1002-1302（2017）20-0254-05

陆地表面温度（land surface temperature，简称LST）是地表能量平衡和资源环境变化的重要参数。地面监测站获取温度信息准确，但不能获取大范围温度值及其时空分布，热红外遥感的出现解决了这个难题。热红外遥感技术能获取热红外波段的辐射能量，根据地表物体的发射率特性反演得到其温度，实现大范围的温度信息获取，因而热红外遥感在地表温度反演方面有着重要的作用[1]。

目前地表温度的反演算法主要有单通道算法、劈窗算法、多通道多角度算法。单通道算法是指只利用1个热红外通道反演地表温度的方法，适用于几乎所有的热红外波段；劈窗算法是指利用相邻的2个热红外通道来进行地表温度反演的方法，是目前为止发展最为成熟的地表温度反演算法；多通道算法还在发展之中，目前还没有一个简便可行的多通道算法可以用来进行地表温度反演。Landsat 8有2个传感器：陆地成像仪（operational land imager，简称OLI）和热红外传感器（thermal infrared sensor，简称TRIS），其中TRIS有2个热红外波段，这虽然为劈窗算法的使用提供了条件，但一直以来美国地质调查局（United States Geological Survey，簡称USGS）对TIRS11波段的定标准确性把握不准，所以不鼓励使用劈窗算法来反演Landsat 8数据的LST[2]。因而本研究基于Landsat 8第10波段使用3种单通道算法，包括大气校正法[3]、JM_SC10 算法[4]、TIRS10_SC算法[5]进行LST反演。

本研究以山东省济南市为研究区，Landsat 8数据为数据源，根据济南市实时的大气水汽含量，选用合适的大气透过率模型进行大气透过率的估算；计算归一化植被指数（normalized difference vegetation index，简称NDVI）、改进型归一化水体指数（modified normalized difference water index，简称MNDWI）、压缩数据维建筑用地指数（index-based build-up index，简称IBI），采用基于分类回归树（classification and regression tree，简称CART）算法的数据挖掘技术，获取用于区分影像中不同类别地物的规则阈值，实现影像地物的精确分类；获取每一类别纯净植被像元与纯净裸土像元的NDVI，摒弃植被覆盖度计算的经验模型，实现研究区植被覆盖度及地表比辐射率的精确计算。在此基础上，对3种算法反演的地表温度进行比较分析，确定精度最高的单通道地表温度反演算法，为后期应用Landsat 8数据进行地表温度的反演提供参考。

1 材料与方法

1.1 研究区概况

济南市为山东省省会，位于山东省中西部，介于36°01′～37°32′N、116°11′～117°44′E。地处中纬度地带，季风明显，四季分明，年平均气温为13.8 ℃，夏季平均气温为26.7 ℃，冬季平均气温在1.0 ℃左右。济南三面环山，南依泰山，北跨黄河，地势南高北低，落差达500 m，平均海拔高度为118 m，这种地势构造令水汽和热空气回流聚集且不易扩散，一旦发生热污染将很难消除，因此实时监测大范围的济南市地表温度是很有必要的。

1.2 数据处理

遥感数据选用2015年4月25号覆盖研究区的Landsat 8 OLI、TIRS数据以及MODIS L1B 1KM数据。其中OLI数据用来计算研究区的地表比辐射率；TIRS数据用来反演LST；MODIS L1B 1KM数据用来计算研究区大气水汽含量；地面实测数据由济南市气象站提供。

数据处理包括MODIS L1B 1KM数据、Landsat 8多光谱数据辐射定标与大气校正、Landsat 8热红外数据辐射定标等，总体技术路线如图1所示。

2 地表温度反演算法

2.1 大气校正法

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

式中：λ为波长，μm；ελ为地表比辐射率；Ts为地表真实温度，K；B（Ts）为黑体热辐射亮度，W/（m2·sr·μm）；τλ为大气在波长λ处的透过率。变换（1）式得温度为Ts的黑体在热红外波段的辐射亮度B（Ts）为：

Ts可以用普朗克公式的函数获取：

对于TIRS 10，K1=774.89 W/（m2·sr·μm），K2=1 321.08 K，将式（3）中B（Ts）替换为Lλ即为亮度温度Tλ的计算表达式。

2.2 Jiménez-Muoz改进单通道算法

JM_SC10算法是由Jiménez-Muoz等于2014年在其原有单通道算法（single channel method，简称SC）的基础上提出来的，增加了针对Landsat 8的大气参数，JM_SC10算法如下：

式中：γ、δ是基于Planck函数的2个参数；bγ为常数，在TIRS10波段为1 324，Lλ表示热辐射亮度值，W/（m2·sr·μm）；Tλ表示亮度温度，K。ψ1、ψ2、ψ3与大气水汽含量（g/cm2）ω有关，当ω>3 g/cm2时，

2.3 TIRS10_SC算法

胡德勇等于2015年提出了TIRS10_SC算法[5]，该算法是在总结辐射传输方程和覃志豪单窗算法的基础上，专门针对Landsat 8 TIRS传感器开发的。

T10为TIRS10的亮温，K；K2=1 321.08 K；Ta为大气平均作用温度，K；根据影像的获取时间，采用中纬度夏季的大气平均作用温度估算方程进行Ta的计算，估算方程如下：

式中：T0为近地表温度（K），可根据当地气象资料获取。

3 大气透过率τ与地表比辐射率ε计算

3.1 大气透过率τ的计算

热红外波段大气透过率主要取决于大气水汽含量，可通过中等光谱分辨率大气透过率算法计算机模型（moderate spectral resolution atmospheric transmittance algorithm and computer model，简称MODTRAN）模拟计算出大气透过率与大气水汽含量的关系，根据获取的水汽含量信息来计算各像元大气透过率。Landsat 8数据很难进行大气水汽含量反演，但是MODIS 数据却可以。MODIS数据包含36个波段，其中17、18、19为大气吸收波段，第2、第5波段为大气窗口波段。毛克彪等研究发现，利用MODIS数据的第2、19波段可反演大气水汽含量[6]，大气水汽含量ω为：

式中：ρ19和ρ2 分别为MODIS数据的第19、第2波段的地表反射率；α、β为常数，分别为0.020 0、0.632 1。

不同研究区、不同水汽含量下，大气透过率的估算模型是不同的。当前针对Landsat 8数据的通用估算模型要求水汽含量在0.5～3.0 g/cm2，而研究区平均的水汽含量为 4.086 g/cm2，因此本研究区不能使用Landsat 8数据通用模型。鉴于Landsat 8数据的TIRS10与MODIS L1B数据的31波段中心波长和波宽相似（图2），那么应用于31波段的大气透过率估算方程，基本也适用于Landsat 8数据的TIRS10。根据中纬度大气剖面数据进行模拟[7]，得到Band31在各种水汽含量下的通用大气透过率估算方程（相关系数为 0.997 48）

通过式（12）实现TIRS10大气透过率的计算，并得到该波段大气透过率影像（图3）。

3.2 地表比辐射率ε计算

地表比辐射率的计算一直是地表温度反演的重点和难点。地表比辐射率别称地表发射率，与地表的物质结构有关。要实现对地表比辐射率的精确计算，前提就是对影像地物進行精确分类。基于CART算法的数据挖掘技术可以获取不同地物的规则阈值，大大提高决策树分类效率和精度。本研究将影像地物分为4类，分别是水体、自然表面（耕地、林地、草地）、城镇表面（建筑物、道路）、裸土，利用CART算法获取的各类规则阈值，确定决策树如图4所示，分类结果影像如图5所示。

根据ASTER光谱库和Nichol的研究成果[8]，获取各地类纯净像元在TIRS10波段的地表比辐射率，水体为0.996 83、建筑为0.964 885、植被为0.986 72、裸土为0.967 67。影像中虽然存在各地类的纯净像元，但也存在大量的混合像元，其中自然表面和裸土区域可以看作是植被与裸土的混合，城镇表面可以看作建筑物与植被的混合，考虑到混合像元这种情况，覃志豪等提出以下模型[9]计算地表比辐射率：

式中：ε自、ε城、ε裸分别为自然地表比辐射率、城镇地表比辐射率以及裸土地表比辐射率；RV、RB、RS分别为植被、建筑、裸土的温度比率；εV、εS、εB分别为纯净植被、裸土、建筑的地表辐射率；鉴于研究区南北高差相距500 m，可依据植被的构成比例简单估计dε，经验模型[10]为：

当0≤PV≤0.5时，dε=0.003 796PV；

当0.5< PV≤1时，dε=0.003 796（1- PV）

需要注意的是，如果应用公式计算出的ε自大于εV，取ε自=εV；ε城>εB，取ε城=εB；ε裸>εS，取ε裸=εS。PV为植被覆盖度，由改进像元二分模型计算：

NDVIS为纯净裸土或者建筑像元的NDVI值；NDVIV为纯净植被像元的NDVI值。温度比率Ri与植被覆盖度PV相关，覃志豪等根据各地表类型的温度差异进行模拟，确定出植被、建筑、裸土表面的温度比率[9]：

计算地表比辐射率需要参数多、运算复杂，总体波段运算表达式如下：

式中：b1、b6、b7、b9分别为植被、裸地、建筑、水体的掩膜影像；b2为植被覆盖度影像；b3、b4、b8分别为RV、RS、RB影像；b5为dε影像。地表比辐射率影像如图6所示。endprint

4 地表温度反演

4.1 地表温度反演结果

逐一计算各个参数之后，进行3种算法的地表温度反演，反演结果（取置信区间99.9%以去除异常值）如图7～9所示。

4.2 精度验证

根据济南市气象网统计，2015年4月25日济南市温度大致在16～31 ℃之间，部分地区温度达到36～40 ℃，与影像反演结果大致相同。结合反演影像得出高温天气主要发生在章丘市以及济南市区部分地区。究其原因是济南市重工业多集中在该地区，如三一重工、济南重工、枭龙重工等，重工业产生的热污染使章丘市与济南城区的部分地区温度居高不下，严重影响到人们的生活。

对3种算法反演的温度数据进行分类统计发现，耕地植被区域温度相对适中，温度在22～37 ℃之间，平均温度为 29 ℃；水体区域温度最低，温度在16～33 ℃之间，其中黄河水域的温度主要集中在17～21 ℃之间，各大水库如鹊山水库、卧虎山水库、玉清湖水库的温度在16～19 ℃之间，温度较高的水体区域多为零星分布的池塘以及海岸地带；裸地和建筑区域温度最高，在28～40 ℃之间，平均温度为33 ℃。整体来看，3种算法对济南市地表温度的反演比较合理。

为精确验证反演结果的精度，本研究选取济南市气象站16个站点温度数据作为实测数据，对3种算法反演的地表温度进行精度验证。各站点实测温度与应用3种算法反演的温度数据对比结果如表1所示。

由表1可看出，3种单通道算法中TIRS10_SC算法反演地表温度误差最小，平均误差为1.10 ℃，绝对平均误差为 1.24 ℃；其次是大气校正法，平均误差为1.40 ℃，绝对平均误差为1.45 ℃；最大误差来自JM_SC10算法反演的地表温度，平均误差为2.85 ℃，绝对平均误差为2.85 ℃。算法反演温度与温度站实测温度不尽相同，其中大部分反演温度略高于实测站温度，分析原因主要有：（1）温度站实测温度的获取时间为10：30，而影像的获取时间为10：48，在晴朗天气下，温度是逐渐上升的。（2）对于算法中应用到的大气向上和向下的辐射亮度L↑、L↓，采用美国航空航天局（National Aeronautics and Space Administration，简称NASA）提供的计算模型。要提高模型计算的精度，输入参数必须要有实时的气压、风速、海拔高度、相对湿度等参数，本研究没有这些数据，使用了模型的缺省值。（3）对JM_SC10算法，大气水汽含量是影响该算法精度的一个关键参数，当大气水汽含量降低到2 g/cm2时，地表温度的反演误差会降低1.5～3.0 ℃[4]。试验中研究区平均大气水汽含量为4.086 g/cm2，因此可以推断，一旦大气水汽含量降低到2 g/cm2，应用JM_SC10算法反演温度的误差有可能至少降低1.5 ℃， 使该算法反演的地表温度误差大大降低。

尽管上述原因造成了一定的误差，但3种算法反演的温度误差平均值为1.78 ℃，绝对误差的平均值为1.85 ℃，反演结果相对理想。

5 讨论

本研究选用3种地表温度反演算法进行地表温度的反演，根据济南地区的实际情况，构建适合济南地区的大气透过率和地表比辐射率参数。本研究创新点有：（1）对大气透过率的计算，考虑到各研究区不同大气状况、不同水汽含量的影响，选用与Landsat 8数据第10波段性质相近的MODIS第31波段的大气透过率估算模型，消除了含水量的影响，摒弃先前整景影像的单一大气透过率计算，实现影像上每像元的大气透过率求解。（2）基于CART算法的数据挖掘技术，结合多源影像，如NDVI影像、MNDWI影像、IBI影像，获取用于区分影像中不同类别地物的规则阈值，实现影像地物的精确分类。在此基础上，获取每一类别纯净植被像元与纯净裸土像元的NDVI值，摒弃植被覆盖度计算的经验模型，求解温度比率，实现研究区植被覆盖度以及地表比辐射率的精确计算。

研究展望：（1）考虑如果能根据济南地区的实际情况利用MODTRAN模型模拟出Landsat 8数据第10波段大气水汽含量与大气透过率的关系模型，对反演精度的提高会大有帮助；（2）获取济南市实时的气压、风速、海拔高度、相对湿度等参数，实现L↑、L↓精确计算；（3）JM_SC10算法提出的前提是ω>3 g/cm2，考虑对JM_SC10算法的使用进行改进，在原有算法的基础上 ±1.5 ℃，反演2幅温度影像，根据部分实测温度数据，确定出真实的温度反演影像，提高算法的精度。

6 结束语

不管是通用的大气校正法还是最近提出的地表温度反演算法，在实现济南市大气透过率和地表比辐射率的精确计算后，各算法均对济南市地表温度的反演有很好的适用性，尤其是TIRS10_SC算法精度相对较高，而当大气水汽含量低于 2 g/cm2 时，JM_SC10算法精度也可大大提高。

参考文献：

[1]宋 挺，段 峥，刘军志，等. Landsat 8数据地表温度反演算法对比[J]. 遥感学报，2015，19（3）：451-464.

[2]徐涵秋. 新型Landsat 8卫星影像的反射率和地表温度反演[J]. 地球物理学报，2015，58（3）：741-747.

[3]Qin Z H，Karnieli A，Berliner P. A mono-window algorithm for retrieving land surface temperature from Landsat TM data and its application to the Israel-Egypt border region[J]. lnternational Journal of Remote Sensing，2001，22（18）：3719-3746.

[4]Jiménez-Muoz J C，Sobrino J A，Skokovic D，et al. Land surface temperature retrieval methods from Landsat-8 thermal infrared sensor data[J]. IEEE Geoscience and Remote Sensing Letters，2014，11（10）：1840-1843.

[5]胡德勇，喬 琨，王兴玲，等. 单窗算法结合Landsat 8热红外数据反演地表温度[J]. 遥感学报，2015，19（6）：964-976.

[6]毛克彪，覃志豪，王建明，等. 针对MODIS数据的大气水汽含量反演及31和32波段透过率计算[J]. 国土资源遥感，2005（1）：26-29.

[7]Mao K B，Qin Z H，Shi J. A practical split-window algorithm for retrieving land surface temperature from MODIS data[J]. International Journal of Remote Sensing，2005，26（15）：3181-3204.

[8]Nichol J. An emissivity modulation method for spatial enhancement of thermal satellite images in urban heat island analysis[J]. Photogrammetric Engineering and Remote Sensing，2009，75（5）：547-556.

[9]覃志豪，李文娟，徐 斌，等. 陆地卫星TM 6波段范围内地表比辐射率的估计[J]. 国土资源遥感， 2004，16（3）：27-32.

[10]Sobrino J A，Jiménez-Munoz J C，Paolini L. Land surface temperature retrieval from Landsat TM 5[J]. Remote Sensing of Environment，2004，90（4）：434-440.endprint