Landsat 8地表温度产品验证与分析

李胜林,刘波,*,赵犇,李大成,王宇恒,郭凤云

1.中国电子科技集团公司第二十二研究所,新乡 453002 2.中国农业科学院 农田灌溉研究所,新乡 453002 3.太原理工大学 矿业工程学院,太原 030024

1 引言

地表温度(land surface temperature,LST)是影响不同尺度地表物理过程的一个关键参数,影响着大气与陆地之间相互作用[1-2],对热岛研究[3-6]、水资源管理[7]、蒸散发[8]等有重要的意义。热红外遥感是在区域甚至全球范围内检索 LST 的有效方法[9-10]。Landsat系列卫星自 1982 年以来一直在捕获热红外图像,空间分辨率为 60～120m,时间分辨率为16～18d[11],再加上其校准良好的传感器和数据存档,使得 Landsat地表温度在长时间序列以及中等空间分辨率的研究中起到了重要的作用[12-13]。

Landsat 8是Landsat系列卫星的最新任务,与之前的Landsat系列不同,它有两个热红外波段(波段10和波段11),空间分辨率为100m。两个或更多热红外波段使得使用多通道算法获得LST成为可能。然而,由于杂散光影响,波11段误差过大,使得传统的劈窗算法难以从Landsat 8热红外传感器(thermal infrared sensor,TIRS)的两个热红外(thermal infrared,TIR)波段中反演LST[14]。即使杂散光校正算法[15]于2017年在美国地质调查局 (United States Geological Survey,USGS)地面系统中实施,Landsat 团队仍然不建议使用波段11,直到完成更多的工作来评估这种校正是否足够[16]。

因此USGS建议用户开发单通道算法进行LST检索。虽然文献[17]表明基于劈窗算法对杂散光校正后的Landsat 8进行温度反演精度可靠,但是仍然需要更多的探索。目前常用的单通道算法主要为辐射传输方程法和经验单通道算法[18],辐射传输方程法利用辐射传输模型与大气廓线进行大气校正,最终反演得到地表温度。文献[19]开发了一种基于网络的工具,使用美国环境预测中心 (National Centers for Environmental Prediction,NCEP) 大气剖面和大气辐射传输软件MODTRAN来推导大气参数,但是该方法只能提供单点大气参数,对于整景数据的LST反演将会带来误差。经验单通道算法主要包括单窗算法[20]和普适性单通道算法[21]。单窗算法通过构建平均作用温度(一般由近地表气温通过经验关系得到)与大气参数的经验关系反演地表温度,然而平均作用温度与近地表气温的经验关系并不适合所有区域[22],在不同区域应用时需要考虑。普适性单通道算法通过构建大气参数与水汽的经验关系反演地表温度,但是在水汽浓度较大的情况下,得到反演结果精度较低[23-24]。另一方面,目前常用的发射率获取方法主要包括NDVI阈值法、植被覆盖度法、分类赋值法等,这些方法适合本地研究,在全球应用时可能会出现问题,特别是在发射率随光谱和空间上变化较大的干旱与半干旱地区[25]。因此,上述方法无法用来检索全球范围内所有土地覆盖类型的LST。

针对上述问题,文献[26]提出了一种为所有Landsat传感器生成 Landsat LST 的算法,并且应用于全球地表温度的生成,得到了全球标准的LST产品。该算法由USGS实施,LST产品在Landsat集合2,2级(Landsat collection 2 level 2)数据中获得。该算法包括3个步骤:①大气校正:使用辐射传输模型和再分析数据对传感器观测到的热辐射进行大气校正;②发射率获取:基于ASTER GEDv3发射率进行光谱调整,然后结合Landsat可见光、短波红外数据进行修改,以解决植被物候变化和积雪变化带来的影响;③温度检索:使用查找表方法对大气校正和发射率校正后的Landsat辐射数据进行LST检索。Landsat collection2 level2数据在2020年中开始提供,包含了全球标准的Landsat LST 产品,为Landsat热波段在全球的广泛应用提供了基础。但是,目前使用的Landsat LST多数是基于传统的单通道算法进行反演得到的,直接使用该温度产品的应用较少,也缺乏足够的验证。成熟的地表温度产品往往需要长时间序列的验证,例如MODS LST产品、VIIRS LST产品等均经过了大量的长时间序列验证[27-28]。针对目前对Landsat LST产品进行精度验证的工作较少,特别是针对长时间序列数据的验证,本研究基于黑河流域生态-水文过程综合遥感观测联合试验(HiWATER)[29-30]中的湿地站(SD)、戈壁站(GB)、大满超级站(CJZ)的地面测量数据对Landsat 8 LST产品进行了长时间序列验证,并将该产品与广泛流行的普适性单通道算法反演结果进行了对比。

2 研究区以及数据介绍

研究区位于甘肃省张掖市,属黑河流域中游,地理范围为38°41'15.33"-39°02'37"(N),100°10'26"-100°43'11"(E)。该地区气候干旱,年平均气温6℃, 区域内地势平坦,主要地貌有戈壁、农田、沙地、建筑等。

遥感数据为研究区内2013-04-09—2016-11-02清晰无云的31景Landsat 8数据以及1km分辨率的MODIS水汽产品(MOD05),MOD05用来提供普适性单通道算法中需要的水汽含量。Landsat 8为collection2 level2数据,该数据包括反射率、地表温度、发射率等19个波段,所有波段以30m分辨率存储,其中地表温度由TIR 10波段反演得到,实际空间分辨率为100m。在研究区内Landsat 8过境时间为11:55左右,选择与Landsat 8过境时间相近的MOD05,过境时间为11:20左右。

地面验证数据使用黑河流域中游SD、GB和CJZ三个自动气象站的地面测量数据作为验证数据,自动气象站每30s观测一次,输出为10min观测的平均值,站点详细信息如表1所示。

表1 验证站点信息

根据卫星过境时间对10min内的观测值进行线性插值,得到卫星过境时刻的实测值。真实地表温度Ts通过架设在自动气象站上的SI-111红外辐射计和四分量辐射仪测量值计算:

式中:Tr为SI-111观测温度(单位K);L↓为通过四分量辐射仪测量到的下行长波辐射(单位W·m-2);ε为地表发射率,由ASTER GEDv3计算得到;σ为斯蒂芬-玻尔兹曼常数(5.67×10-8W·m-2·K-4)。

3 LST反演算法

3.1 普适性单通道算法

普适性单通道算法又称JMS算法[21],并且该算法在后续的研究中进行了改进[23],本文使用改进后的算法。从 Landsat 8 TIR 10波段检索地表温度Ts,该算法表示为:

Ts=γ[ε-1(ψ1L+ψ2)+ψ3]+δ

式中:L为辐射值,根据TIR 10波段观测值进行辐射定标得到;ε为发射率;γ,δ的计算公式如下:

(1)

式中:Tb为TIR 10的亮度温度;bγ=C2/λ,其中C2=1.43877×104μm·K,λ为TIR 10的有效波长,对于TIR 10,bγ=1324K;ψ1、ψ2和ψ3为大气功能参数,是大气水汽含量ω的函数[24]:

ψ1=0.04019ω2+0.02916ω+1.01523

ψ2=-0.38333ω2-1.50294ω+0.20324

ψ3=0.00918ω2+1.36072ω-0.27514

上述基于Landsat 8 TIR 10进行地表温度的检索中需要两个重要的参数,地表发射率ε和大气水汽含量ω,在本文中使用的地表发射率来自Landsat8 collection2 Level2数据中的发射率波段,将发射率波段乘缩放系数0.0001即得到实际发射率值。大气水汽含量来自MOD05数据,在温度反演前将MOD05数据投影至WGS84_UTM_47N,并重采样至30m分辨率,重采样方法使用最邻近法。

3.2 Landsat 8 LST 产品反演算法

首先使用发射率、大气参数和Landsat 8 TIR 10波段的辐射数据计算地面辐射Ls:

(2)

式中:L↑为上行辐射;τ为大气透过率;B为普朗克公式。

理论上来讲Ts可由如下公式求得:

式中:λ为TIR10波段的有效波长;C1=1.19104×108W·μm4·m-2·sr-1;C2=14387.7μm·K。然而,对于偏离狄拉克δ函数的传感器光谱响应,该公式将越来越不准确。因此,使用查找表(LUT)来计算Landsat 8 LST。对地面所有可能的温度(通常为150～380K)进行模拟,计算Landsat 8传感器在0.01K温度范围内光谱响应的预期辐射,生成一张Landsat 8传感器接收辐射值与相应温度值的对应表。Landsat 8 LST可以通过使用式(2)求得的地面辐射Ls和生成的表格来检索地表温度。对于大气参数、发射率的计算以及算法详细信息,请参考文献[26]。

对于collection2 Level2中的Landsat 8 LST产品,观测值(observed)与实际温度之间K的关系可由如下公式计算:

LST(K)=0.00341802×observed+149.0

4 研究结果与分析

本研究基于站点实测数据对2013—2016年研究区内清晰无云的31景Landsat 8 LST产品进行了验证,并将Landsat 8 LST产品与广泛流行的JMS算法的反演结果进行了对比。评价指标采用相关系数R2,平均偏差Bias和均方根误差RMSE。

图1给出了2015-08-28的地表温度反演结果以及差值图。从图1可以看出,两类LST空间分布一致,都反映出了不同地类的温度差异,可以明显看出植被区域温度较低,周边的戈壁和沙地温度较高,并且两类LST偏差均在-1～1K,具有较小的误差分布。

图1 2015-08-28地表温度反演结果及差值Fig.1 Land surface temperature inversion results and difference(2015-08-28)

为了更加详细地分析两类LST的差异,将研究区分为植被、城镇、戈壁、沙地4个地类,分析不同地类的LST平均偏差(LST产品与JMS LST差值的平均值),每个地类选择一个80×80像元的区块进行研究,被选择的区块位置以及地表覆盖状况如图2所示。1为植被区域,2为戈壁区域,3为城镇区域,4为沙地。两类LST不同地类平均偏差的时间序列如图3所示。各地类大部分日期平均偏差大于0,反演温度LST 产品整体上高于JMS LST。对于各个地类间平均偏差对比,大部分日期城镇地区平均偏差最小,沙地平均偏差最大,但整体来看,同一日期影像不同地类间的平均偏差大小接近,随着地类的不同平均偏差变化较小。各地类的平均偏差随着日期的不同均有较大的变化,这可能与水汽含量相关。从图3可以看出平均偏差的变化趋势与水汽含量变化趋势(研究区内水汽含量的平均值)比较一致,在平均偏差较大的几个峰值,研究区内水汽含量大小也在峰值附近。图4显示了研究区内两类LST的平均偏差与水汽含量变化的散点图,当水汽含量较低时,特别是水汽含量小于1g·cm-2时,平均偏差几乎在0.5K以内,当水汽含量增大时,两类LST的平均偏差整体也出现增大的趋势。

上述对LST产品和JMS LST进行了对比分析,为了进一步对两类LST的绝对精度进行评估,下面基于SD、GB、CJZ三个站点实测数据对两类LST进行验证。验证前需要对地面站点的异质性进行分析。一般选择站点周围3×3像素进行站点异质性分析,由于Landsat 8 LST产品被采样为30m存储 ,而实际分辨率为100m,因此本文采用站点周围9×9像素进行异质性分析。统计站点周围9×9像素温度,并计算31景LST产品平均标准差作为异质性分析依据。得到统计结果显示: SD站为0.85K,GB站为0.412K,CJZ站为0.491K,三个站点的平均标准差均小于1K,满足理想地表温度验证点的不确定度应在1K以内的要求[31]。

图2 不同地类的位置以及地表覆盖Fig.2 The location and surface coverage of different land types

图3 不同地类平均偏差对比以及与水汽含量关系Fig.3 Comparison of bias of different land types and their relationship with water vapor content

图4 研究区内两类LST的平均偏差与水汽含量散点图Fig.4 Scatter plot of the bias and water vapor content of the two types of LST in the study area

图5为站点处两类LST与站点实测数据的散点图,从图中可以看出两类LST与地面测量数据具有很高的一致性,在每个站点处R2均高于0.949。图6显示了两类LST在不同站点处与地面测量LST的偏差。与地面测量LST比较,两类LST在SD站与GB站大部分日期偏差为负,低估了地表温度。对于湿地站,站点异质性强,地面测量数据作为像元尺度的验证值时,可能会出现站点实际观测值比像元尺度实际值高的情况,造成两类LST低估。对于GB站,站点处温度随时间变化敏感,而站点观测温度是根据卫星过境时间对10min内的观测值进行线性插值得到,因此无法与卫星过境时间时的温度准确对应,造成低估的情况。在GB站两类LST RMSE在3.500K以上,明显高于SD站与CJZ站,这可能与GB站验证数据较少有关。在CJZ站,两类LST相比站点实测数据在大部分日期出现了高估的情况,因为该站点下垫面为玉米,发射率变化随季节变化敏感,高估的情况可能是温度反演中发射率计算不准确造成。表2为两类LST在各个站点处的精度指标。对比LST 产品与JMS LST,在SD站与GB站,LST 产品精度指标Bias与RMSE均小于JMS LST, LST产品精度优于JMS LST。在CJZ站,JMS LST(Bias=-0.051K,RMSE=2.028K)精度优于LST产品(Bias=0.578K,RMSE=2.230K)。在所有站点处,LST产品 Bias为-0.714K,RMSE为2.865K, JMS LST Bias为-1.252K,RMSE为2.958K。整体来讲, LST产品与JMS LST均有较高的精度,RMSE均在3K以内,LST产品的精度稍高于JMS LST。

从图4可看出,两类LST的偏差随着水汽含量的不同而有明显的变化,在水汽含量增大时偏差有增大趋势,基于所有站点对研究区内水汽大于1g·cm-2和小于1g·cm-2两种情况进行两类LST验证,如表3所示。水汽含量小于1g·cm-2时,LST产品和JMS LST分别低估了1.617K和1.909K,同时LST产品RMSE小于JMS LST,LST产品精度稍高。在水汽大于1g·cm-2时,LST产品高估了0.730K, JMS LST低估了0.201K,JMS LST具有更好的平均偏差,但是RMSE稍高于LST产品。整体来看,在本研究区的水汽范围内,基于站点的验证结果两类LST精度相似。由于研究区内水汽含量整体较低,所有日期水汽含量均小于3g·cm-2,因此无法对高水汽含量情况下两类LST的精度进行评估。

图5 各站点反演温度与验证温度散点图Fig.5 Scatter plots of retrieved temperature and verified temperature at each site

图6 各站点反演温度与验证温度偏差Fig.6 The bias of the retrieved temperature and the verified temperature at each site

表2 基于各站点的精度评估结果

表3 不同水汽含量情况下所有站点精度评估结果

5 结论

Landsat 8 LST产品具有重要的应用价值,其精度评估亦需要展开。本文以黑河流域中游为研究区,基于JMS 反演结果和地面站测量数据对2013-2016研究区内31景清晰无云的Landsat 8 LST产品进行了评估与分析,可以得出以下结论:

1)Landsat 8 LST产品与JMS反演结果相似,大部分日期偏差在1K以内,整体上反演温度Landsat 8 LST产品稍高于JMS反演结果。

2)通过对比两类LST之间的平均偏差,不同的地类对偏差大小影响较小,水汽含量的大小是影响两类LST偏差的主要因素。 当水汽含量小于1g·cm-2时,大部分偏差在0.5K以内,水汽含量大于1g·cm-2,偏差有增大的趋势,最大偏差达到了2.5K。

3)基于站点数据验证,两类LST整体上低估了地表温度,RMSE均在3K以内,Landsat 8 LST产品的评价指标(Bias=-0.714K,RMSE=2.865K)优于JMS LST(Bias=-1.252K,RMSE=2.958K),LST 产品的精度稍高。

本研究区内水汽含量较低,均小于3g·cm-2,因此缺乏高水汽含量情况下的验证与评估,本文的验证精度不够全面,今后需要在高水汽含量地区进行实验。