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.
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, 2013,2014):
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的函数,可用下式获得:鉴于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:
式中T10和T11为TIR 10、11波段的亮温(K); ε 为TIR 10、11波段的平均比辐射率,即ε= 0.5(ε10 + ε11),Δε为TIR 10、11波段的比辐射率差值,即Δε=(ε10-ε11), 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(θ)为对应太阳天顶角的大气透射率.从以上公式可以看出,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,并进行比较(表 3、4).由于2幅影像皆为同日影像,且前后间隔不到7 min,因此大气条件相同,可忽略公式(4)和(23)中的τ.另外,由于这2幅影像受大气影响的 只有可见光和近红外波段,因此表 4仅列出它们的对比结果.图 1以受大气影响最严重的蓝光波段为例,列出基于ρλ_COST的Landsat 7和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的大气校正.从表 3、4中可以看出,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可知,在单通道算法中,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. |