地球物理学报  2015, Vol. 58 Issue (3): 741-747   PDF    
新型Landsat 8卫星影像的反射率和地表温度反演
徐涵秋    
福州大学环境与资源学院、福州大学遥感信息工程研究所, 福建省水土流失遥感监测评价重点实验室, 福州 350108
摘要:Landsat 8卫星自2013年2月发射以来,其影像的定标参数经过了不断调整和完善,针对Landsat 8开发的各种算法也相继问世.本文采用最新的参数、算法和引入COST算法建立的大气校正模型,对Landsat 8多光谱和热红外波段进行了处理,反演出它们的反射率和地表温度,并与同日的Landsat 7数据和实测地表温度数据进行了对比.结果表明,现有Landsat 8多光谱数据的定标参数和大气顶部反射率反演算法已有很高的精度,本文引入COST算法建立的Landsat 8大气校正模型也与Landsat 7的COST模型所获得的结果几乎相同,相关系数可高达0.99.但是现有针对Landsat 8提出的地表温度反演算法仍不理想,已提出的劈窗算法误差都较大.鉴于TIRS 11热红外波段的定标参数仍不理想,因此在现阶段建议采用单通道算法单独反演TIRS 10波段来求算地表温度,但要注意根据大气水汽含量的情况选用正确的大气参数计算公式.
关键词Landsat 8     影像处理     地表温度     反射率     遥感    
Retrieval of the reflectance and land surface temperature of the newly-launched Landsat 8 satellite
XU Han-Qiu    
College of Environment and Resources; Institute of Remote Sensing Information Engineering; Fujian Provincial Key Laboratory of Remote Sensing of Soil Erosion; Fuzhou University, Fuzhou 350108, China
Abstract: The quality of Landsat 8 sensor data has been improved after a series of corrections and data reprocessing since its successful launch on February 11, 2013. These include the correction of all calibration parameters, the improvement of the radiance conversion coefficients for the all Operational Land Imager (OLI) sensor bands, the refinement of the OLI detector linearization, the radiometric offset correction for the two Thermal Infrared Sensor (TIRS) bands, the slight improvement to the geolocation of the TIRS data, and the reprocessing of all Landsat 8 data held in the USGS archives. In addition, several algorithms specially developed for Landsat 8 data have also been proposed over the pass year. This paper aims to assess the accuracy of the retrieved the reflectance of the OLI sensor and the land surface temperature (LST) of the TIRS sensor of the new satellite with those of the well-calibrated Landsat 7 ETM+ sensor and the ground-measured LST. This study retrieved the top of the atmosphere (TOA) reflectance of the OLI multispectral bands and the LST of the TIRS thermal bands with the most recent calibration parameters, algorithms and USGS-processed Landsat 8 image data. In addition, the Chavez's COST model has been introduced to correct the atmospheric effects on the OLI multispectral bands. To examine the performance of the calibration parameters and the developed algorithms, the retrieved TOA reflectance of each multispectral band and the computed LST of the thermal infrared bands have been compared with that of the corresponding band of the synchronized, well-calibrated Landsat 7 data and the in situ LST, respectively.
The results show that the current Landsat 8 calibration parameters of multispectral bands can achieve high accuracy for the retrieval of TOA reflectance. The proposed COST-based atmospheric correction algorithm can also have a nearly identical performance when compared with the Landsat 7 COST model's result. Nevertheless, two recently-proposed split window algorithms for computing the LST from Landsat 8 thermal infrared bands did not perform well, as they offered a large difference between the algorithm-modeled LST and the ground-measured LST. Given the scaling parameters of the TIRS thermal infrared band 11 is still unstable, as announced by the Landsat 8 project team, it is recommended that at this stage users might use the single channel (SC) algorithm of Jiménez-Muñoz and Sobrino to retrieve the LST from Landsat 8 thermal band 10 (like working on Landsat TM/ETM+ band 6) rather than attempt a split-window algorithm using both TIRS bands 10 and 11. However, care should be taken in the selection of correct atmosphere parameters for the SC-based LST computing, especially when a very high atmospheric water vapor condition occurs.
Key words: Landsat 8     Image processing     Land surface temperature     Reflectance     Remote sensing    

1 引言

在诸多的卫星对地观测系统中,Landsat系列卫星以其40多年的观测历史而成为应用最为广泛的系列卫星之一.其最新一代的Landsat 8卫星自2013年2月11日发射以来,迄今已在轨运行了一年多,它以每天550幅的速率向地面传回了海量的对地观测数据,并与Landsat-7 ETM+一起,构成了8天间隔的Landsat重复观测周期,为地球表面的生态环境变化监测提供了极为宝贵的数据.与之前的Landsat系列卫星对比,Landsat 8在许多方面进行了明显的改进(USGS,2013).发射一年多来,Landsat 8卫星影像的各种处理文件和参数又经过了多次修改,同时还有一些针对Landsat 8开发的算法也相继问世.这些参数和算法究竟是否有效?它们的应用潜力究竟如何?当前国内外尚无相应的研究.因此,本文将对这些参数变化、相关算法及其应用处理进行梳理.在此基础上,根据Landsat 8最新的定标参数、算法和引入COST算法建立的大气校正模型,对其多光谱波段反射率和热红外波段地表温度的反演进行探讨,查明它们的应用情况,以使用户可以有效、正确地对Landsat 8的多光谱和热红外数据进行处理. 2 Landsat 8发射以来与影像处理有关的重要事件

Landsat 8于2013年2月11日发射后,经过了一系列的测试,于当年3月18日发回了第一幅影像.运行一年多来,Landsat 8的影像处理经历了数次重要的变化(USGS, 20132014):

2013年4月10日,Landsat 8卫星到达了预期的705 km运行高度(称为WRS-2),而在此之前接收的数据由于未达到运行高度,其空间分辨率不是设计的30 m,所以被统称为Pre-WRS-2数据,这点在下载数据和数据处理时必须注意;

5月17日,Landsat 8官方网站公布了Landsat 8影像的大气顶部反射率转换公式和初始定标参数;

9月16日,Landsat项目团队发现TIRS热红外传感器的10波段需要进行偏差校正;

10 月25 日,Landsat 8 网站发布公告称,其之前发布的热红外数据由于定标参数有误差,因此必须对TIRS 10 波段求出的光谱辐射值减去0.32 W/(m2·sr·μm),而TIRS 11 波段由于误差更大而未给出修正方案,USGS也因此暂不鼓励用劈窗算法来反演地表温度,而是建议仍采用TM/ETM+的单波段方式来计算地表温度;

11月14日的公告又将TIRS 10波段的调整参 数由上述的0.32改为0.29,同时给出了TIRS 11 波段的调整参数,要求将所求出的辐射值减去0.51 W/(m2·sr·μm);

2014年1月29日,Landsat网站发布公告称:由于一年多来定标参数的不断改进,USGS 将于2月3日撤下所有已经上网的Landsat 8数据,在用新的定标参数进行重新处理后,再提供给用户下载,这一数据重新处理预计将耗时50天;

2月3日,USGS开始对存档数据进行重新处理.这一处理对OLI传感器多光谱数据的影响很小(辐射值<2 %,反射率<0.8 %),但对TIRS传感器热红外数据的影响很大,根据上述2013年11月14日发布的定标调整参数,重新处理后的TIRS 10波段的温度会降低2.1 K,11波段会降低4.4 K.但可能是对TIRS 11波段定标的准确性仍无把握,USGS还是不鼓励用劈窗算法来反演TIRS的地表温度.

以上一系列的数据和参数更动,必须引起用户在影像处理时的足够重视,如果用户手上的Landsat 8影像是下载于2014年2月3日之前,则最好重新下载经最新定标参数处理过的数据. 3 Landsat 8影像处理 3.1 多光谱影像处理

Landsat 8 具有2个传感器,分别为多光谱OLI传感器和热红外TIRS传感器.对于多光谱的OLI数据,用户可根据Landsat 8网站提供的公式和参数(USGS,2013),将影像的亮度值转换为大气顶部的光谱辐射值,即

式中 Lλ为波段λ的大气顶部光谱辐射值(TOA spectral radiance),ML为波段λ的调整因子,可从头文件(MTL)中获得,为语句Radiance_Mult_B and _x后的数值,x代表波段号,AL为波段λ的调整参数,在MTL文件中为语句Radiance_Add_ B and _x后的数值,Qcal 为影像以16位量化的亮度值(DN).

Landsat 8数据也可以通过下式转换为大气顶部的反射率(TOA Reflectance),即

式中ρ′ λ为波段λ未经太阳角度纠正的大气顶部反射率,Mρ为波段λ的反射率调整因子,在MTL文件中为语句Reflectance_Mult_B and _x后的数值,Aρ为波段λ的反射率调整参数,在MTL文件中为语句Reflectance_Add_B and _x后的数值.

ρ′ λ可经过太阳天顶角的进一步纠正转换为大气顶部反射率ρλ

式中θZ为影像中心的太阳天顶角,也可采用太阳高度角θE,但要采用正弦函数sinθE.

由于所获得的ρλ并未经过大气校正,因此本文引入在Landsat TM/ETM+中广为应用的Chavez(1996)的COST大气校正算法来建立Landsat 8的大气校正模型(记为ρλ_COST).COST模型所需大气参数少,主要以影像自身的数据进行校正(Xu,2006),即

式中Qh是大气影响的修正值,可以通过最暗像元 法获得(Chavez,1996);τ = cos [(90-θZ/180], 为基于θZ估算的大气透射率.但由于这种计算经常造成τ的高估,特别是在晴空无云、太阳天顶角很大 或北方高纬度地区这3种情况(RS/GIS Laboratory,2014),因此在实际应用中τ常被忽略(Ramsey et al., 2004; Lowry et al., 2005). 3.2 热红外影像处理 3.2.1 亮温反演

Landsat 8的TIRS热红外传感器具有2个波段,分别编号为TIRS 10和TIRS 11,它们可以通过下式转换为亮温:

式中T为传感器处的亮温(K),L为公式(1)求出的热红外波段的辐射值,K1和K2为热红外波段的定标常数,对于TIRS 10波段,K1=774.89 W/(m2·sr·μm),K2=1321.08 K;对于TIRS 11波 段,K1=480.89 W/(m2·sr·μm),K2=1201.14 K. 3.2.2 地表温度反演

当前,针对Landsat 8热红外波段地表温度(LST)反演所提出的算法有Jiménez-Muñoz等(2014)的单通道算法(以下简称SC).尽管USGS在现阶段不提倡用劈窗算法来反演LST,但他们仍同时提出了针对Landsat 8热红外波段的劈窗算法; 提出劈窗算法的还有Rozenstein等(2014).

(1)单通道算法(Jiménez-Muñoz et al., 2014): Jiménez-Muñoz等(2014)在其原有的SC算法的基础上,增加了针对Landsat 8的大气参数,其算法为

式中ε为地表比辐射率,可通过ASTER光谱库获得主要地物在TIRS 10、11波段的比辐射率,表 1为本文综合ASTER光谱库和Nichol(2009)的研究成果计算的比辐射率;参数bγ分别为:TIRS 10=1324 K,TIRS 11=1199 K;L和T分别由公式(1)和(5)获得,ψ1、ψ2、ψ3为大气水汽含量w的函数,可用下式获得:
表 1 Landsat 8 TIR 10、11波段的比辐射率 Table 1 Land surface emissivity for Landsat 8 TIRS B and s 10 and 11

式中的pij(i=1,2,3,j=1,2,3)是与大气水汽含量w相关的大气参数.可能是由于TIRS 11波段的定标仍不稳定,Jiménez-Muñoz等只给出了针对TIRS 10波段的特定参数:

鉴于SC算法在w>3 g·cm-2时精度会受到明显影响,因此Jiménez-Muñoz等(2014)建议此时不用公式(8)—(9)计算ψ1ψ2ψ3,而是用它们的原始推导公式来计算,即

式中τ为大气透过率,L和L为大气上行和下行辐射强度.

(2)Jiménez-Muñoz等(2014)的劈窗算法:该算法计算简单,在求出TIR 10、11波段的亮温后,可由下式计算出LST:

式中T10T11为TIR 10、11波段的亮温(K); ε 为TIR 10、11波段的平均比辐射率,即ε= 0.5(ε10 + ε11),Δε为TIR 10、11波段的比辐射率差值,即Δε=(ε1011), w为总大气水汽含量(g·cm-2); 参数 c0 ~ c6分别为:-0.268、1.378、0.183、54.3、-2.238、 -129.2、16.4.

(3)Rozenstein等(2014)的劈窗算法:该算法基于Qin等(2001)的算法提出,其基本公式为

式中A0、A1、A2是由大气透射率和地表比辐射率所确定的TIR 10、11波段的参数,它们可以由以下公式计算:

式中ai、bi分别为TIRS 10、11波段根据不同温度范围确定的回归系数(表 2),Ci、Di是由大气透射率和比辐射率所确定的参数,由下式计算:

式中εi为TIRS 10、11波段对应的比辐射率;τi(θ)为对应太阳天顶角的大气透射率.
表 2 不同温度范围内的a、b取值 Table 2 Coefficient a and b for different temperature ranges
4 对比和讨论 4.1 多光谱波段

从以上公式可以看出,Landsat 8影像与之前的Landsat 7影像无论是在多光谱或热红外波段的处理上都有一定的区别.其中多光谱波段最主要的区别表现在大气顶部反射率ρλ的计算上.Landsat 7的ρλ及其基于COST算法的大气校正模型为(Chander et al., 2009; Ramsey et al., 2004; Lowry et al., 2005):

将以上公式分别与公式(3)、(4)对比可以看出,Landsat 8的ρλ的计算已大为简化,除了太阳天顶角外,已经不再需要大气顶部平均太阳辐照度(ESUN)、日-地天文单位距离(d)等参数.但可能是出于惯性思维,仍有很多用户在许多相关网站上讨论是否需要ESUN等参数,有的还模拟计算了针对 Landsat 8的ESUN参数(http://www.researchgate.net).

由于Landsat 7的ETM+数据是经过了严格定标的准确数据,因此,为了验证Landsat 8究竟是否需要ESUN等参数,以及本文引入COST模型建立的大气校正公式(4)是否有效,本次研究采用Landsat 8网站提供的2对位于美国西部、于2013年3月29日同日过空的Landsat 7和Landsat 8的影像对,分别计算出它们的大气顶部反射率ρλ和经过COST模型大气校正的反射率ρλ_COST,并进行比较(表 34).由于2幅影像皆为同日影像,且前后间隔不到7 min,因此大气条件相同,可忽略公式(4)和(23)中的τ.另外,由于这2幅影像受大气影响的 只有可见光和近红外波段,因此表 4仅列出它们的对比结果.图 1以受大气影响最严重的蓝光波段为例,列出基于ρλ_COST的Landsat 7和Landsat 8的对比影像.

图 1 经过ρλ_COST模型大气校正的Landsat 7和Landsat 8的蓝光波段大气顶部反射率对比影像(a—d) 以及3种地表温度算法反演的福州2013-8-4的地表温度影像(e—g) Fig.1 Comparison of the ρλ_COST-corrected TOA reflectance of blue b and s between Landsat 7 and Landsat 8(a—d) and the LST images of Fuzhou(8-4-2013)retrieved using SC-10,SW-JM, and SW-R algorithms(e—g)

表 3 基于ρλ的Landsat 7和Landsat 8比较 Table 3 Comparison of ρλ between the corresponding b and s of Landsat 7 and Landsat 8

表 4 基于ρλ-COST的Landsat 7和Landsat 8比较 Table 4 Comparison of ρλ-COST between the corresponding b and s of Landsat 7 and Landsat 8

表 3、4可以看出,无论是大气顶部反射率(ρ)或经大气校正的反射率(ρλ_COST),Landsat 8各波段的均值与Landsat 7各对应波段的均值都非常接近,其各波段均值的平均差最大只有0.003,误差最大也只有2.7%,相关系数r最小的也都接近0.95.显然,Landsat 8已经不再需要ESUN等参数,其本质差别就是Landsat 8用公式(3)的ρ′ λ取代了 Landsat 7公式(22)中的Lλ,而计算ρ′ λMρ和Aρ 参数(见公式(2))本身已经是基于反射率定标的调整系数.

表 4表明,引入COST算法的Landsat 8大气校正模型计算的结果与Landsat 7 COST模型的结果非常接近,特别是影像对2的最小平均误差几乎为0,相关系数高达0.99,说明该模型可以和Landsat 7一样,用于Landsat 8的大气校正.从表 34中可以看出,Landsat 8和Landsat 7各对应波段的差值在近红外、中红外波段通常要大于可见光波段,这是由于Landsat 8在近红外和中红外波段的波长设置范围与Landsat 7的对应波段有较明显的差别所致(徐涵秋、唐菲,2013). 4.2 热红外波段

对于热红外波段地表温度反演算法的验证,我

们用上述2种劈窗算法和SC算法对2幅福州Landsat 8 影像的地表温度进行反演(日期为2013-7-3、2013-8-4),然后与当日同时的福州国家基准气候站实测的地表温度进行对比.根据实验区所处地理方位与季相,将中纬度夏季大气廓线作为MODTRAN大气剖面数据输入,模拟得到大气参数.与计算有关的大气水汽含量w、大气透射率τ、大气上行辐射强度L、大气下行辐射强度L分别为:7月3日:4.35,0.3742,4.2094,6.0665;8月4日:4.18,0.4494,4.1208,6.1377.表 5是3种算法反演的地表温度与实测地表温度对比的结果,图 1e—1g是以2013-8-4影像为例的3种算法反演的地表温度影像.

表 5 各种地表温度反演算法的精度比较 Table 5 Comparison between the algorithm-retrieved LST and in situ LST data

表 5可知,在单通道算法中,SC-10所计算的地表温度LST的误差为-1.53~-2.77 ℃,平均为-2.15 ℃,而SC-9的误差很大.分析其原因可能与2幅实验影像当天的大气水汽含量w过高(都超过了4 g·cm-2)有关.Jiménez-Muñoz等指出,当w>3 g·cm-2时,就不宜使用SC-9提供的大气参数(公式(9))来反演LST,而必须用SC-10的算法.从本次实验的结果来看,当w过高时,SC-10算法的精度确实大大高于SC-9的.

就劈窗算法而言,2种算法的精度都很低.劈窗算法原本可借助Landsat 8提供的TIRS 10、11这两个波段来消除大气的影响(Irons et al., 2012),但本次实验获得的精度却很不理想,这可能与TIRS 11波段的定标参数仍不稳定有关.也正因为此,USGS的多次公告一直不鼓励用劈窗算法来反演TIRS的LST,而是建议仍采用Landsat TM/ETM+的单波段方法,通过单独反演TIRS 10波段来求取LST.

总的来看,本次实验的结果与USGS的建议是吻合的,使用劈窗算法估算的Landsat 8的LST误差较大,而用SC-10单独反演TIRS 10波段所获得的LST的误差最小.虽然本次实验的SC-10算法仍平均低估了-2.15 ℃,但前已述及,这可能是大气水汽含量过高引起的.根据Jiménez-Muñoz等(2009)的实验结果,当w下降到2 g·cm-2时,LST的估算误差会降低1.5~3 ℃.因此可以预测,当大气水汽含量较低时,SC算法的精度可能会明显提高,从而使我们看到了SC在反演Landsat 8地表温度上的潜力. 5 结论

Landsat 8发射一年多来,其影像的定标参数经过了不断调整,所有已上网的影像也经过了重新处理.本次研究表明,现有Landsat 8多光谱数据的定标参数和大气顶部反射率反演算法已有较高的精度,与经过严格定标的Landsat 7同步数据对比没有明显的差别.本文引入COST算法建立的大气校正模型也与Landsat 7 COST模型所获得的结果几乎一致.与Landsat 7对比,Landsat 8 影像大气顶部反射率的计算已大为简化,也不再需要ESUN等参数.

但是现有针对Landsat 8提出的地表温度劈窗算法的反演精度并不理想,由于TIRS 11热红外波段的定标参数仍不理想,因此已提出的2个劈窗算法的误差都较大.在现阶段,建议采用SC算法单独反演TIRS 10波段来求算地表温度,待TIRS 11热红外波段的定标参数稳定后,再采用劈窗算法.但采用SC算法时,要根据大气水汽含量的情况选用正确的大气参数计算公式.

参考文献
[1] Chander G, Markham B L, Helder D L. 2009. Summary of current radiometric calibration coefficients for Landsat MSS, TM, ETM+, and EO-1 ALI sensors. Remote Sens. Environ., 113(5): 893-903.
[2] Chavez P S Jr. 1996. Image-based atmospheric corrections-Revisited and revised. Photogram. Eng. Remote Sens., 62: 1025-1036.
[3] Irons J R, Dwyer J L, Barsi J A. 2012. The next Landsat satellite: the Landsat data continuity mission. Remote Sens. Environ., 122: 11-21.
[4] Jiménez-Muñoz J C, Sobrino C J, Soria J A, et al. 2009. Revision of the single-channel algorithm for land surface temperature retrieval from Landsat thermal-infrared data. IEEE Trans. Geosci. Remote Sens., 47(1): 339-349.
[5] Jiménez-Muñoz J C, Sobrino J A, Skokovic D, et al. 2014. Land surface temperature retrieval methods from landsat-8 thermal infrared sensor data. IEEE Geosci. Remote Sens. Lett., 11(10): 1840-1843, doi: 10.1109/LGRS.2014.2312032.
[6] Lowry J, Kirby J, Langs L. 2005. SWReGAP Land Cover Mapping Methods Documentation. http://earth.gis.usu.edu/swgap/data/landcover/map_methods/AZ/AZ4_Methods.pdf.[2014-6-3].
[7] Nichol J. 2009. An emissivity modulation method for spatial enhancement of thermal satellite images in urban heat island analysis. American Society for Photogrammetry and Remote Sensing, 75(5): 547-556.
[8] Qin Z H, Dall'Olmo G, Karnieli A, et al. Derivation of split window algorithm and its sensitivity analysis for retrieving land surface temperature from NOAA-advanced very high resolution radiometer data. J. Geophys. Res.: Atmos., 2001, 106(D19): 22655-22670.
[9] Ramsey R D, Wright D L Jr, McGinty C. 2004. Evaluating the use of Landsat 30m Enhanced Thematic Mapper to monitor vegetation cover in shrub-steppe environments. Geoc. Int., 19(2): 39-47.
[10] Rozenstein O, Qin Z H, Derimian Y, et al. 2014. Derivation of land surface temperature for Landsat-8 TIRS using a split window algorithm. Sensors, 14(4): 5768-5780.
[11] RS/GIS Laboratory at Utah State University. 2005. Southwest Regional Gap Analysis Project Land Cover Map. http://www.gis.usu.edu/current_proj/swregap_landcover.html.[2014-6-3].
[12] USGS. 2013. Landsat 8 Data. http://landsat.usgs.gov/landsat8.php..
[13] USGS. 2014. Landsat Headlines. http://landsat.usgs.gov/mission_headlines2014.php.[2014-6-3].
[14] Xu H Q, Tang F. 2013. Analysis of new characteristics of the first Landsat 8 image and their eco-environmental significance. Acta Ecol. Sinica (in Chinese), 33(11): 3249-3257.
[15] Xu H Q. 2006. Evaluation of two absolute radiometric normalization algorithms for pre-processing of Landsat imagery. J. China Univ. Geosci., 17(2): 146-150, 157.
[16] 徐涵秋, 唐菲. 2013. 新一代Landsat 系列卫星: Landsat 8遥感影像新增特征及其生态环境意义. 生态学报, 33(11): 3249-3257.