基于Landsat8数据的地表温度反演及分析研究——以武汉市为例

吴志刚，江 滔，樊艳磊，陈联君

(1.武汉市测绘研究院，湖北 武汉 430022；2.中国地质大学 信息工程学院，湖北 武汉 430074)



基于Landsat8数据的地表温度反演及分析研究
——以武汉市为例

吴志刚1，江 滔1，樊艳磊1，陈联君2

(1.武汉市测绘研究院，湖北 武汉 430022；2.中国地质大学 信息工程学院，湖北 武汉 430074)

地表温度的反演方法经过多年发展，多种算法相继被提出。为了研究新型遥感Landsat8影像条件下，各地表温度反演算法的优劣，利用2013年9月Landsat8 数据，以武汉市武昌区为研究对象，使用辐射传导方程法(大气校正法)、单窗算法对研究区进行了地表温度反演，并与其亮度温度进行了对比研究，结果表明：这两种算法反演的地表温度和亮度温度的空间分布情况大体相同，但也存在着差异，单窗法和辐射传导法均值温度相差1.3K，与亮度温度均值相比，单窗方法高1.6K，辐射传导方程法高2.9K。单窗算法精度优于辐射传导方程法。

辐射传导方程法；单窗算法；地表温度

1 引 言

地表温度在地—气相互作用过程中有着非常重要的作用，是地球系统水热平衡的关键因子，在灾害监测预防、气候变化、天气预报和军事目标识别等众多方面有着重大的意义[1,2]。获得地表温度的传统方法是利用地面气象站提供的地表温度观测资料，但由于观测站密度较低等问题，不能全面反映地面的热状况[3]。

伴随着遥感技术的发展进步，其在较多领域有着广泛的应用[4-6]。特别是热红外遥感技术的进步，使其能够大范围快速地监测地表温度。20世纪60年代起，遥感数据开始用于地表温度的反演，国内外学者相继提出了很多温度反演算法，主要有：单波段算法、多波段算法、单波段多角度法和多波段多角度法等[7-11]。

单波段的方法，主要包括辐射传导方程法(大气校正法)、单通道算法和单窗算法。辐射传导方程法需要较为精确的大气剖面数据；单窗算法和单通道算法相对精确。这些算法都需要进行地表辐射率的计算[12]。

由于辐射传导方程法(大气校正法)和单窗方法都只需要一个热红外波段就可以进行相关研究，简单易行；多波段算法需较多气象资料反演相关参数，较为复杂。因此本文选择了单波段进行反演。

本文以武汉市武昌区为研究区域，利用最新的Landsat8数据，采用辐射传导方程法和单窗算法对研究区地表温度进行反演，并对结果进行分析比较。

2 研究区概况

武汉市是湖北省会，位于江汉平原东部，地处113°41′E～115°05′E，29°58′N～31°22′N，属北亚热带季风性(湿润)气候，常年雨量充沛，热量丰富、四季分明。夏季最长约135天；春秋二季各约60天。全年的平均气温为17.6℃，1月平均气温最低，为0.4℃，7、8月份平均气温最高，可达28.7℃。因武汉地处北纬30°，夏天太阳高度大，地处内陆离海洋较远，地形如盆地容易集热，河湖多故夜间水汽较多，加上城市热岛效应和副热带高气压带控制，天气闷热，极端气温可达44.5℃。武昌区是武汉市一个市辖区，位于武汉城区东南，是旧武昌市的主体部分。研究区的位置详见图1。

图1 研究区位置及Landsat8影像Fig.1 Location of the study area and Landsat8 data

3 数据源及预处理

本文选用的是2013年9月17日的Landsat8遥感数据，PATH/ROW是123/039。Landsat8数据来源于美国地质调查局相关网站[13]，该影像无云，质量较好。该卫星于2013年2月11日发射，具有两个传感器Operational Land Imager(OLI)和Thermal Infrared Sensor(TIRS)[14]，OLI除了具有Landsat7所有光谱波段之外，还增加了一个深蓝波段(Deep blue)、海岸/气溶胶波段(Coastal/Aerosol)和用于检测卷云的卷云波段(Crirus)，同时收窄了全色波段和近红外波段的光谱范围，Landsat8各个波段的辐射分辨率从8bit提高至12bit，增大了影像的灰度量化级，信噪比提高[15]。TIRS有两个热红外波段，比ETM+增加了一个热波段，分辨率为100 m，这使其大气校正更加容易，可用劈窗算法进行热红外校正[16]或者ENVI的Thermal Atm Correction工具。Landsat7与Landsat8卫星参数对比见表1。

表1 Landsat7和Landsat8数据参数对比

研究区温度、相对湿度、气压等数据，主要通过历史气象分享网站(http://weather.bsyan.com/)获得。

文中使用的研究区行政边界矢量图是从全国县级行政区矢量图提取得到。Landsat8数据用的UTM/WGS84 投影/坐标系，数据产品Level 1T已进行了基于地形的几何校正，可直接使用。

OLI波段数据的大气校正主要使用ENVI的FLAASH模块进行。因TIRS有两个热红外波段，可使用ENVI Thermal Atm Correction工具进行热红外波段的大气校正。在大气校正之前需要进行辐射定标，将像元灰度值转换为热辐射强度值(Radiance)。辐射定标的公式如下：

Lλ=MLQcal+AL

(1)

ML为增益参数，AL为偏移参数，两参数都可以从影像元数据文件中获得：ML对应元数据中RADIANCE_MULT_BAND_x，AL对应RADIANCE_ADD_BAND_x，x为相应波段数，Qcal为相应波段的灰度值[17]。

波段10参数RADIANCE_MULT_BAND_10 = 3.3420E-04，RADIANCE_ADD_BAND_10 = 0.10000；波段11的参数相同，因此两个热红外波段辐射亮度公式为：

Lλ=3.3420E-04×Qcal+0.1

(2)

TIRS有两个热红外波段B10和B11，文中所使用的是B10波段。

4 地表温度反演算法

4.1 辐射传导方程法

在地—气的辐射传输中，卫星接收到的热红外辐射能量Lλ包含3部分的内容：地面真实辐射经大气衰减之后被卫星传感器接收到的热辐射能量、大气的上行辐射亮度L↑、大气下行辐射亮度L↓(大气向地面热辐射)。卫星接收到的热红外辐射亮度值表达即辐射传导方程：

Lλ=[ε·B(TS)+(1-ε)L↓]·τ+L↑

(3)

式(3)中ε地表辐射率；TS是地面真实温度单位为k；τ为大气在热红外波段的透过率。大气下行辐射亮度L↓和大气上行辐射亮度L↑单位均为W·m-2·sr-1·μm-1。

辐射传导方程法(Radioactive Transfer Equation)，又被称作大气校正法，其主要的原理是根据实时的大气探测数据、大气廓线数据(或大气模型，如6S、MODTRAN等)来评估大气对地表热辐射的影响，并从卫星观测的热辐射总量中去除此部分大气影响，得到真实的地面热辐射强度并转换成相应的地表温度[18]。假设地表、大气对热辐射具有朗伯体性质，则根据辐射传导方程可得出与地表真实温度相同的黑体的辐射亮度B(TS)，公式为

(4)

式(4)中的透过率τ、大气上行辐射亮度L↑(W·m-2·sr-1·μm-1)、大气下行辐射亮度L↓(W·m-2·sr-1·μm-1)三个参数可以通过NASA官网(http://atmcorr.gsfc.nasa.gov/)输入影像的成像时间、中心经纬度、相关地区气压等相关信息生成，黑体的辐射亮度B(TS)单位为W·m-2·sr-1·μm-1。本文参与计算的影像成像时间为2013年9月17日02时58分、影像中心经纬度113.8192E和30.3011N，气压为1011百帕，相对湿度为46%，最后得到的大气参数透过率τ为0.66，大气上行辐射亮度L↑为3.10 W·m-2·sr-1·μm-1，大气下行辐射亮度L↓为4.86 W·m-2·sr-1·μm-1。由于NASA官网暂时只能获取B10波段的参数，而不能获取B11波段的参数，为了和下文中的单窗算法进行严谨的比较，因此统一使用B10波段。

估算出地表真实温度相同的黑体的辐射亮度B(TS)后，根据普朗克定律反函数，得出地面真实温度。公式如下：

(5)

K1,K2为常数，Landsat8 TIRS波段10的K1值为K1_CONSTANT_BAND_10，波段10的K2值为K2_CONSTANT_BAND_10[19]，波段11的参数与其类似：

K1_CONSTANT_BAND_10=774.89，K2_CONSTANT_BAND_10 = 1 321.08

K1_CONSTANT_BAND_11= 480.89，K2_CONSTANT_BAND_11 = 1 201.14

4.2 单窗算法

辐射传导方程法由于对大气剖面数据、探空数据的依赖，覃志豪等提出了一种基于TM数据的地表温度反演算法，即单窗算法(MonoK-window Algorithm，简称为MW算法)[20]。公式为：

TS={a(1-C-D)+[b(1-C-D)+C+D]Tb-DTa}/C

(6)

式(6)中TS是地表反演温度(K)，Tb为亮度温度(K)，Ta为大气平均作用温度(K)，a、b数值取-67.355 35和0.458 61；C、D为中间量，通过公式(7)和(8)计算得到。Ta估算根据中纬度夏季平均大气廓线，由公式(9)得到[18]。

C=ε×τ

(7)

D=(1-τ)[1+(1-ε)τ]

(8)

Ta=16.011 0+0.926 21T0

(9)

公式中ε为地表比辐射率，τ为大气透射率，T0为近地面气温，单位为K。

根据历史气象资料，武汉市当时近地面气温为31℃，得出Ta为297.72 K。

通常的大气透射率τ结合表2和大气水分含量w进行估算[21]。

表2 大气透射率估算方程

大气水分含量的计算，需要获取研究区域的气温、相对湿度数据，使用公式(10)估算出e(绝对水汽压)。可表示[22]为：

(10)

式(10)中RH为相对湿度，T0为气温(K)，e的单位为千帕。

利用绝对水汽压e和杨景梅[23]等研究确定的大气水分含量和地面水汽压的关系式，估算出大气水汽含量。公式如下：

W=a0+a1e

(11)

式(11)中a0、a1为经验系数(分别取0.178 8和0.197 8)，e为绝对水汽压，本式中e单位为百帕，根据公式参数得出W=3.16 g/cm2。

大气透射率τ利用杨槐[24]Landsat8数据的水汽和透过率关系进行估算：

τ10=-0.106 7w+1.040 2 R2=0.994 8

(12)

亮温的估算，利用公式(2)先将像元灰度值转化为热辐射强度值Lλ，然后根据普朗克定理将其转化为亮度温度，公式表示为：

(13)

式(13)中K1、K2值与公式(5)中相同。

4.3 地表比辐射率计算

地表比辐射率是计算地表温度的重要参数之一，是不同地表温度反演方法的共同参数[25]。其对地表温度精度影响较大。估算方法主要有经验公式法和混合像元法。

经验公式方法为Van de Griend等人于1993年发现。经研究他认为地表的比辐射率和归一化植被指数有较好的正相关关系，相关系数为0.94，因此地表比辐射率可用NDVI进行计算[26]，公式如下：

ε=1.009 4+0.047ln(NDVI)

(14)

混合像元方法主要有Sobrino[27]的基于地表覆盖类型的加权混合模型和覃志豪[28]等的混合模型。

Sobrino混合模型认为地表由植被和裸地构成[25]。用NDVI进行地表分类：

1)当NDVI<0.2时，则认为全部由裸地覆盖，地表比辐射率取裸地典型发射率值0.973;

2)当0.2≤NDVI≤0.5时，认为像元是由植被和裸地构成的混合像元，则地表比辐射率由简化公式(15)计算：

ε=0.004PV+0.986

(15)

3)当NDVI>0.5时，认为完全为植被覆盖地面，地表比辐射率取植被典型发射率0.986。

覃志豪等的混合模型认为，计算地表比辐射率除了考虑自然表面，还应考虑城镇、水体这两种地表类型。城镇可看作城镇和植被的混合，自然表面看作植被和裸土的混合。水体、植被、土壤和城镇的比辐射率[29]为0.995、0.986、0.973和0.970。

城镇和自然表面比辐射率分别用化简的公式(16)和(17)估算：

(16)

(17)

式(15)、(16)、(17)中：PV为植被覆盖度，其估算用公式(18)进行：

(18)

式(18)中NDVIV、NDVIS分别为植被和裸地的NDVI值。

本文估算地表辐射率，是对Sobrino混合模型进行了改进，认为地表由植被、裸地和水体构成，根据像元统计，得到植被茂密区NDVI值约为0.72，裸地均值为0.12，Landsat8 B10波段的水体、裸地，植被比辐射率取0.992、0.973 1、0.984 4，利用NDVI进行地表分类：

1)当NDVI≤0，认为像元由水体覆盖，比辐射率取0.992；

2)当0

3)当0.12

ε=0.004PV+0.984 4

(19)

4)当NDVI≥0.72，认为像元完全由植被构成，比辐射率为0.984 4，由公式(18)计算。

5 温度反演专题制图及分析

利用研究区域Landsat8数据(热波段使用B10波段)结合辐射传导方程法、单窗算法分别得到武汉市武昌区地表温度的反演结果，对其进行密度分割，并在ARCGIS中进行专题制图。图2为研究区域亮度温度影像(简称BT影像)；图3为辐射传导方程法地表温度影像(简称LST-RTE影像)；图4为单窗算法地表温度影像(简称LST-MW影像)。

由图2～图4可知，亮度温度得到的地表温度和其他两种算法反演的地表温度空间分布情况基本相同，只是变化的幅度稍有差异。其主要的地表温度序列为：城市居民区温度>裸地温度>植被温度>水体温度。城市由于其下垫面比自然地表太阳吸收率高以及大气污染、人工热源等因素，温度表现较高；裸地和绿地相比，因为含水量、比热容的不同导致温度的差异。水体温度在白天日照时，温度较低；在夜晚时由于水比热容大，温度较高。

图2 武昌区亮度温度影像(℃)Fig.2 Lightness temperature image of Wuchang districts

图4 武昌区单窗算法地表温度影像(℃)Fig.4 Land surface temperature image based on the mono-window in Wuchang districts

图3 武昌区辐射传导方程法地表温度影像(℃)Fig.3 Land surface temperature image based on the RTE in Wuchang districts

根据反演温度的影像统计得出：图2为亮度温度转化的地表温度影像(BT影像)，最高温度为43℃(316.2K)，最低温度为16℃(289.1K)，均值温度为28.8℃(301.9K)；图3为LST-RTE影像，最高温度为48℃(321.2K)，最低温度为15℃(288.2K)，均值温度为31.7℃(304.8K)；图4为LST-MW影像，最高温度为45℃(318.5 K)，最低温度为17℃(290.2K)，均值温度为30.4℃(303.5K)(表3)。LST-MW影像均值温度和亮温均值相差1.2℃，LST-RTE影像均值温度和亮温均值相差2.4℃，LST-RTE影像均值温度和LST-MW影像均值温度相差3.6℃(表4)。

上述分析可以看出，辐射传导方程法(大气校正法)和单窗方法所得均温都高于亮度温度。两者与亮度温度比较可以看出，单窗算法反演精度优于辐射传输算法，主要是由于辐射传输算法对大气剖面数据、探空数据较为依赖，且准确获取数据较为困难，本文所使用大气剖面数据是NASA网站获取数据，精度不高；单窗算法的参数主要依赖于大气水分，精度较高。

表3 研究区各算法影像温度范围及均值

表4 研究区各算法温度差值范围和均值差值

为研究各算法对不同地表覆盖类型温度反演的差异，提取城镇、自然表面(植被裸地混合)和水体3种地表覆盖类型对各影像进行掩膜处理和分析的结果。

经过统计可以得到：BT影像的水体均值温度为25.6℃，植被裸地混合均值温度为29.5℃，城镇均值温度为30.2℃；LST-MW影像水体均值温度为26.2℃，植被裸地混合均值温度为31.7℃，城镇均值温度为32.8℃；LST-RTE影像水体均值温度为26.4℃，植被裸地混合均值温度为32.8℃，城镇均值温度为33.9℃(表5)。

表5 研究区水体、自然表面、城镇各算法均值温度

6 结 论

本文基于Landsat8数据分别使用了辐射传导方程法(大气校正法)、单通道算法对研究区进行了地表温度的反演研究，并对其结果和研究区亮度温度进行了差值比较，还对不同的地表覆盖类型地表温度进行了分类统计，以便对不同算法结果和亮度温度、不同地物类型温度差异和分布规律进行研究。以下几点结论：

1)两种算法反演的地表温度和亮度温度的空间分布情况大体相同，只是变化的幅度有差异；

2)反演的地表均值温度与亮度温度相比，单窗方法高1.6K，辐射传导方程法高2.9K，单窗法和辐射传导法均值温度相差1.3K；

3)不同类型地表在不同反演方法下均值相差不大，辐射传导法反演值高，单窗法反演值较为接近亮度温度值。

本文使用的温度和绝对水汽压等参数值是武汉市大范围的历史记录值，和研究区值存在一定偏差。Landsat8 TIRS具有两个热波段数据，本文只使用了B10数据，如利用劈窗算法两个波段数据都参与反演，反演精度会进一步提升。

[1]徐希孺，柳钦火，陈家宜.遥感陆面温度[J].北京大学学报(自然科学版)，1998，34(2-3)：248-253.

[2]孟宪红，吕世华，张宇，等.使用LANDSAT-5TM数据反演金塔地表温度量衡[J].高原气象，2005，24(5):721-726.

[3]谭桂容，蔡哲，徐永明.基于Landsat影像的南京地区热岛效应[J].安徽农业科学，2009,37(13):6 050-6 052.

[4]刘建华.遥感技术在土地变更调查中的应用研究[J].工程地球物理学报，2014，11(6)：901-904.

[5]郑杨琳，赵楠，王华伟，等.QuickBird影像在低中放废物处置场选取中的应用研究[J].工程地球物理学报，2013,10(3):415-419.

[6]黄伦春，肖玉环，张红英.基于PCA变换的多元遥感影像融合方法研究[J].工程地球物理学报，2014，11(2)：266-270.

[7]Qin Z, 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]. International Journal of Remote sensing，2001，22:3 719-3 746.

[8]Jimenez-Munoz J C, Sobrino J A. A Generalized Single channel Method for Retrieving Land Surface Temperature from Remote Sensing Data[J]. Journal of Geophysics Research，2003，108:4 688-4 697.

[9]俞红,石汉青.利用分裂窗算法反演陆地表面温度的研究进展[J].气象科学，2002，22(4):494-500.

[10]覃志豪, Zhang M H, Kamieli A,等.用陆地卫星TM6数据演算地表温度的单窗算法[J].地理学报，2001，56 (4) :456-466.

[11]So brino J A, Li ZL, Stoll MP, et al. Multi-channel and multi-angle algorithms for estimating sea and land surface temperature with ATSR data [J]. Internati nal Journal of Remote Sensing，1996，17(11): 2 089-2 114.

[12]樊辉.基于Landsat TM 热红外波段反演地表温度的算法对比分析[J].遥感应用，2009(1)：36-40.

[13]United States Geological Survey/[EB/OL]. http://glovis.usgs.gov，2014-5-4.

[14]USGS. Landsat Data Continuity Mission(LDCM)-Landsat8[EB/OL]. http://landsat.usgs.gov/LDCM_Landsat8.php，2014-5-4.

[15]徐涵秋,唐菲.新一代Landsat系列卫星：Landsat8遥感影像新增特征及其生态环境意义[J].生态学报，2013，33(11)：3 249-3 257.

[16]Irons J R, Dwyer J L, Barsi J A. The next Landsat satellite:The Landsat Data Continuty Misson[J].Remote Sensing of Environment，2012，122:11-21.

[17]Using the USGS Landsat 8 Product[EB/OL]. http://landsat.usgs.gov/Landsat8_Using_Product.php，2014-5-4.

[18]游绚,晏路明.基于ETM+影像辐射传导方程算法的地表温度反演[J].科技情报开发与经济，2009，19(27)：134-136.

[19]陈云.基于Landsat 8的城市热岛效应研究初探——以厦门市为例[J].测绘与空间地理信息，2014,37(2)：123-128.

[20]Qin Z H, Karnieli A, Berliner P. A Mono-window Algorithm of Retrieving Land Surface Temperature from Landsat and TM Data and its Application to the Israel-Egypt Border Region[J]. International Journal of Remote Sensing，2001,18(2):3 719-37 46.

[21]覃志豪,Li Wenjuan,Zhang Minghua,等.单窗算法的大气参数估计方法[J].国土资源遥感，2003(2)：37-43.

[22]乔建民.不同土地覆盖类型对城市地表温度的影响——以山东省龙口市为例[J].宁波农业科技，2012(3)：21-25.

[23]杨景梅,邱金桓.用地面湿度参量计算我国整层大气可降水量及有效水汽含量方法的研究[J].大气科学，2002,26(1)：9-22.

[24]杨槐.从Landsat8影像反演地表温度的劈窗算法研究[J].测绘地理信息，2014，39(4)：73-77.

[25]丁凤,徐涵秋.TM 热波段图像的地表温度反演算法与实验分析[J].地球信息科学，2006， 8(3)：125-130.

[26]Van de Griend A A，Owe M. On the relationship between thermal emissivity and the normalized different vegetation index for natural surfaces.[J].International Journal of Remote Sensing，1993,14(6)：1 119-1 131.

[27]Sobrino J A, Raissouni N, Li Zhao Liang. A comparative study of land surface emissivity retrieval from NOAA data[J].Remote Sensing Environment,2001,75(2):256-266.

[28]覃志豪,李文娟,徐斌,等.陆地卫星TM6波段范围内地表比辐射率的估计[J].国土资源遥感，2004,61(3)：28-32.

[29]郑国强,鲁敏,张涛,等.地表比辐射率求算对济南市地表温度反演结果的影响[J].山东建筑大学学报，2010,25(5)：519-523.

Land Surface Temperature Retrieval and Result Analysis Based on Landsat8 Data in Wuhan City

Wu Zhigang1，Jiang Tao1，Fan Yanlei1，Chen Lianjun2

(1.WuhanGeomaticInstitute,WuhanHubei430022,China; 2.CollegeofInformationEngineering,ChinaUniversityofGeosciences,WuhanHubei430074,China)

After years of development of land surface temperature inversion method, many algorithms have been proposed. In order to study the advantages and disadvantages of surface temperature inversion algorithm in the new remote sensing Landsat8 images, the Landsat8 image data are used in the paper acquired on September 17, 2013. Taking Wuchang district, in Wuhan as an example, a comparison between the LSTs retrieved from the radioactive transfer equation algorithm(RTE) and mono-window algorithm(MW) is made. The research results show that the LSTs from two algorithms have similar overall temperature distributions, but the difference of the average LST also exists. When compared to the land surface temperature, the difference of mean value between the RTE and MW is 1.3 K, and when compared to the mean brightness temperature, the LST retrieved from the MW is about 1.6 K higher than the brightness temperature, from the RTE is 2.9 K higher. mono-window algorithm(MW) is superior to the radioactive transfer equation algorithm(RTE).

radioactive transfer equation algorithm; mono-window algorithm; land surface temperature

1672—7940(2016)01—0135—08

10.3969/j.issn.1672-7940.2016.01.023

吴志刚(1986-)，男，助理工程师，主要从事资源环境与遥感研究。E-mail: 457903327@qq.com

P237

A

2015-06-13