Landsat 8地表溫度產(chǎn)品驗(yàn)證與分析

李勝林,劉波,*,趙犇,李大成,王宇恒,郭鳳云

1.中國電子科技集團(tuán)公司第二十二研究所,新鄉(xiāng) 453002 2.中國農(nóng)業(yè)科學(xué)院 農(nóng)田灌溉研究所,新鄉(xiāng) 453002 3.太原理工大學(xué) 礦業(yè)工程學(xué)院,太原 030024

1 引言

地表溫度(land surface temperature,LST)是影響不同尺度地表物理過程的一個(gè)關(guān)鍵參數(shù),影響著大氣與陸地之間相互作用[1-2],對(duì)熱島研究[3-6]、水資源管理[7]、蒸散發(fā)[8]等有重要的意義。熱紅外遙感是在區(qū)域甚至全球范圍內(nèi)檢索 LST 的有效方法[9-10]。Landsat系列衛(wèi)星自 1982 年以來一直在捕獲熱紅外圖像,空間分辨率為 60～120m,時(shí)間分辨率為16～18d[11],再加上其校準(zhǔn)良好的傳感器和數(shù)據(jù)存檔,使得 Landsat地表溫度在長(zhǎng)時(shí)間序列以及中等空間分辨率的研究中起到了重要的作用[12-13]。

Landsat 8是Landsat系列衛(wèi)星的最新任務(wù),與之前的Landsat系列不同,它有兩個(gè)熱紅外波段(波段10和波段11),空間分辨率為100m。兩個(gè)或更多熱紅外波段使得使用多通道算法獲得LST成為可能。然而,由于雜散光影響,波11段誤差過大,使得傳統(tǒng)的劈窗算法難以從Landsat 8熱紅外傳感器(thermal infrared sensor,TIRS)的兩個(gè)熱紅外(thermal infrared,TIR)波段中反演LST[14]。即使雜散光校正算法[15]于2017年在美國地質(zhì)調(diào)查局 (United States Geological Survey,USGS)地面系統(tǒng)中實(shí)施,Landsat 團(tuán)隊(duì)仍然不建議使用波段11,直到完成更多的工作來評(píng)估這種校正是否足夠[16]。

因此USGS建議用戶開發(fā)單通道算法進(jìn)行LST檢索。雖然文獻(xiàn)[17]表明基于劈窗算法對(duì)雜散光校正后的Landsat 8進(jìn)行溫度反演精度可靠,但是仍然需要更多的探索。目前常用的單通道算法主要為輻射傳輸方程法和經(jīng)驗(yàn)單通道算法[18],輻射傳輸方程法利用輻射傳輸模型與大氣廓線進(jìn)行大氣校正,最終反演得到地表溫度。文獻(xiàn)[19]開發(fā)了一種基于網(wǎng)絡(luò)的工具,使用美國環(huán)境預(yù)測(cè)中心 (National Centers for Environmental Prediction,NCEP) 大氣剖面和大氣輻射傳輸軟件MODTRAN來推導(dǎo)大氣參數(shù),但是該方法只能提供單點(diǎn)大氣參數(shù),對(duì)于整景數(shù)據(jù)的LST反演將會(huì)帶來誤差。經(jīng)驗(yàn)單通道算法主要包括單窗算法[20]和普適性單通道算法[21]。單窗算法通過構(gòu)建平均作用溫度(一般由近地表氣溫通過經(jīng)驗(yàn)關(guān)系得到)與大氣參數(shù)的經(jīng)驗(yàn)關(guān)系反演地表溫度,然而平均作用溫度與近地表氣溫的經(jīng)驗(yàn)關(guān)系并不適合所有區(qū)域[22],在不同區(qū)域應(yīng)用時(shí)需要考慮。普適性單通道算法通過構(gòu)建大氣參數(shù)與水汽的經(jīng)驗(yàn)關(guān)系反演地表溫度,但是在水汽濃度較大的情況下,得到反演結(jié)果精度較低[23-24]。另一方面,目前常用的發(fā)射率獲取方法主要包括NDVI閾值法、植被覆蓋度法、分類賦值法等,這些方法適合本地研究,在全球應(yīng)用時(shí)可能會(huì)出現(xiàn)問題,特別是在發(fā)射率隨光譜和空間上變化較大的干旱與半干旱地區(qū)[25]。因此,上述方法無法用來檢索全球范圍內(nèi)所有土地覆蓋類型的LST。

針對(duì)上述問題,文獻(xiàn)[26]提出了一種為所有Landsat傳感器生成 Landsat LST 的算法,并且應(yīng)用于全球地表溫度的生成,得到了全球標(biāo)準(zhǔn)的LST產(chǎn)品。該算法由USGS實(shí)施,LST產(chǎn)品在Landsat集合2,2級(jí)(Landsat collection 2 level 2)數(shù)據(jù)中獲得。該算法包括3個(gè)步驟:①大氣校正:使用輻射傳輸模型和再分析數(shù)據(jù)對(duì)傳感器觀測(cè)到的熱輻射進(jìn)行大氣校正;②發(fā)射率獲取:基于ASTER GEDv3發(fā)射率進(jìn)行光譜調(diào)整,然后結(jié)合Landsat可見光、短波紅外數(shù)據(jù)進(jìn)行修改,以解決植被物候變化和積雪變化帶來的影響;③溫度檢索:使用查找表方法對(duì)大氣校正和發(fā)射率校正后的Landsat輻射數(shù)據(jù)進(jìn)行LST檢索。Landsat collection2 level2數(shù)據(jù)在2020年中開始提供,包含了全球標(biāo)準(zhǔn)的Landsat LST 產(chǎn)品,為L(zhǎng)andsat熱波段在全球的廣泛應(yīng)用提供了基礎(chǔ)。但是,目前使用的Landsat LST多數(shù)是基于傳統(tǒng)的單通道算法進(jìn)行反演得到的,直接使用該溫度產(chǎn)品的應(yīng)用較少,也缺乏足夠的驗(yàn)證。成熟的地表溫度產(chǎn)品往往需要長(zhǎng)時(shí)間序列的驗(yàn)證,例如MODS LST產(chǎn)品、VIIRS LST產(chǎn)品等均經(jīng)過了大量的長(zhǎng)時(shí)間序列驗(yàn)證[27-28]。針對(duì)目前對(duì)Landsat LST產(chǎn)品進(jìn)行精度驗(yàn)證的工作較少,特別是針對(duì)長(zhǎng)時(shí)間序列數(shù)據(jù)的驗(yàn)證,本研究基于黑河流域生態(tài)-水文過程綜合遙感觀測(cè)聯(lián)合試驗(yàn)(HiWATER)[29-30]中的濕地站(SD)、戈壁站(GB)、大滿超級(jí)站(CJZ)的地面測(cè)量數(shù)據(jù)對(duì)Landsat 8 LST產(chǎn)品進(jìn)行了長(zhǎng)時(shí)間序列驗(yàn)證,并將該產(chǎn)品與廣泛流行的普適性單通道算法反演結(jié)果進(jìn)行了對(duì)比。

2 研究區(qū)以及數(shù)據(jù)介紹

研究區(qū)位于甘肅省張掖市,屬黑河流域中游,地理范圍為38°41'15.33"-39°02'37"(N),100°10'26"-100°43'11"(E)。該地區(qū)氣候干旱,年平均氣溫6℃, 區(qū)域內(nèi)地勢(shì)平坦,主要地貌有戈壁、農(nóng)田、沙地、建筑等。

遙感數(shù)據(jù)為研究區(qū)內(nèi)2013-04-09—2016-11-02清晰無云的31景Landsat 8數(shù)據(jù)以及1km分辨率的MODIS水汽產(chǎn)品(MOD05),MOD05用來提供普適性單通道算法中需要的水汽含量。Landsat 8為collection2 level2數(shù)據(jù),該數(shù)據(jù)包括反射率、地表溫度、發(fā)射率等19個(gè)波段,所有波段以30m分辨率存儲(chǔ),其中地表溫度由TIR 10波段反演得到,實(shí)際空間分辨率為100m。在研究區(qū)內(nèi)Landsat 8過境時(shí)間為11:55左右,選擇與Landsat 8過境時(shí)間相近的MOD05,過境時(shí)間為11:20左右。

地面驗(yàn)證數(shù)據(jù)使用黑河流域中游SD、GB和CJZ三個(gè)自動(dòng)氣象站的地面測(cè)量數(shù)據(jù)作為驗(yàn)證數(shù)據(jù),自動(dòng)氣象站每30s觀測(cè)一次,輸出為10min觀測(cè)的平均值,站點(diǎn)詳細(xì)信息如表1所示。

表1 驗(yàn)證站點(diǎn)信息

根據(jù)衛(wèi)星過境時(shí)間對(duì)10min內(nèi)的觀測(cè)值進(jìn)行線性插值,得到衛(wèi)星過境時(shí)刻的實(shí)測(cè)值。真實(shí)地表溫度Ts通過架設(shè)在自動(dòng)氣象站上的SI-111紅外輻射計(jì)和四分量輻射儀測(cè)量值計(jì)算:

式中:Tr為SI-111觀測(cè)溫度(單位K);L↓為通過四分量輻射儀測(cè)量到的下行長(zhǎng)波輻射(單位W·m-2);ε為地表發(fā)射率,由ASTER GEDv3計(jì)算得到;σ為斯蒂芬-玻爾茲曼常數(shù)(5.67×10-8W·m-2·K-4)。

3 LST反演算法

3.1 普適性單通道算法

普適性單通道算法又稱JMS算法[21],并且該算法在后續(xù)的研究中進(jìn)行了改進(jìn)[23],本文使用改進(jìn)后的算法。從 Landsat 8 TIR 10波段檢索地表溫度Ts,該算法表示為:

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

式中:L為輻射值,根據(jù)TIR 10波段觀測(cè)值進(jìn)行輻射定標(biāo)得到;ε為發(fā)射率;γ,δ的計(jì)算公式如下:

(1)

式中:Tb為TIR 10的亮度溫度;bγ=C2/λ,其中C2=1.43877×104μm·K,λ為TIR 10的有效波長(zhǎng),對(duì)于TIR 10,bγ=1324K;ψ1、ψ2和ψ3為大氣功能參數(shù),是大氣水汽含量ω的函數(shù)[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進(jìn)行地表溫度的檢索中需要兩個(gè)重要的參數(shù),地表發(fā)射率ε和大氣水汽含量ω,在本文中使用的地表發(fā)射率來自Landsat8 collection2 Level2數(shù)據(jù)中的發(fā)射率波段,將發(fā)射率波段乘縮放系數(shù)0.0001即得到實(shí)際發(fā)射率值。大氣水汽含量來自MOD05數(shù)據(jù),在溫度反演前將MOD05數(shù)據(jù)投影至WGS84_UTM_47N,并重采樣至30m分辨率,重采樣方法使用最鄰近法。

3.2 Landsat 8 LST 產(chǎn)品反演算法

首先使用發(fā)射率、大氣參數(shù)和Landsat 8 TIR 10波段的輻射數(shù)據(jù)計(jì)算地面輻射Ls:

(2)

式中:L↑為上行輻射;τ為大氣透過率;B為普朗克公式。

理論上來講Ts可由如下公式求得:

式中:λ為TIR10波段的有效波長(zhǎng);C1=1.19104×108W·μm4·m-2·sr-1;C2=14387.7μm·K。然而,對(duì)于偏離狄拉克δ函數(shù)的傳感器光譜響應(yīng),該公式將越來越不準(zhǔn)確。因此,使用查找表(LUT)來計(jì)算Landsat 8 LST。對(duì)地面所有可能的溫度(通常為150～380K)進(jìn)行模擬,計(jì)算Landsat 8傳感器在0.01K溫度范圍內(nèi)光譜響應(yīng)的預(yù)期輻射,生成一張Landsat 8傳感器接收輻射值與相應(yīng)溫度值的對(duì)應(yīng)表。Landsat 8 LST可以通過使用式(2)求得的地面輻射Ls和生成的表格來檢索地表溫度。對(duì)于大氣參數(shù)、發(fā)射率的計(jì)算以及算法詳細(xì)信息,請(qǐng)參考文獻(xiàn)[26]。

對(duì)于collection2 Level2中的Landsat 8 LST產(chǎn)品,觀測(cè)值(observed)與實(shí)際溫度之間K的關(guān)系可由如下公式計(jì)算:

LST(K)=0.00341802×observed+149.0

4 研究結(jié)果與分析

本研究基于站點(diǎn)實(shí)測(cè)數(shù)據(jù)對(duì)2013—2016年研究區(qū)內(nèi)清晰無云的31景Landsat 8 LST產(chǎn)品進(jìn)行了驗(yàn)證,并將Landsat 8 LST產(chǎn)品與廣泛流行的JMS算法的反演結(jié)果進(jìn)行了對(duì)比。評(píng)價(jià)指標(biāo)采用相關(guān)系數(shù)R2,平均偏差Bias和均方根誤差RMSE。

圖1給出了2015-08-28的地表溫度反演結(jié)果以及差值圖。從圖1可以看出,兩類LST空間分布一致,都反映出了不同地類的溫度差異,可以明顯看出植被區(qū)域溫度較低,周邊的戈壁和沙地溫度較高,并且兩類LST偏差均在-1～1K,具有較小的誤差分布。

圖1 2015-08-28地表溫度反演結(jié)果及差值Fig.1 Land surface temperature inversion results and difference(2015-08-28)

為了更加詳細(xì)地分析兩類LST的差異,將研究區(qū)分為植被、城鎮(zhèn)、戈壁、沙地4個(gè)地類,分析不同地類的LST平均偏差(LST產(chǎn)品與JMS LST差值的平均值),每個(gè)地類選擇一個(gè)80×80像元的區(qū)塊進(jìn)行研究,被選擇的區(qū)塊位置以及地表覆蓋狀況如圖2所示。1為植被區(qū)域,2為戈壁區(qū)域,3為城鎮(zhèn)區(qū)域,4為沙地。兩類LST不同地類平均偏差的時(shí)間序列如圖3所示。各地類大部分日期平均偏差大于0,反演溫度LST 產(chǎn)品整體上高于JMS LST。對(duì)于各個(gè)地類間平均偏差對(duì)比,大部分日期城鎮(zhèn)地區(qū)平均偏差最小,沙地平均偏差最大,但整體來看,同一日期影像不同地類間的平均偏差大小接近,隨著地類的不同平均偏差變化較小。各地類的平均偏差隨著日期的不同均有較大的變化,這可能與水汽含量相關(guān)。從圖3可以看出平均偏差的變化趨勢(shì)與水汽含量變化趨勢(shì)(研究區(qū)內(nèi)水汽含量的平均值)比較一致,在平均偏差較大的幾個(gè)峰值,研究區(qū)內(nèi)水汽含量大小也在峰值附近。圖4顯示了研究區(qū)內(nèi)兩類LST的平均偏差與水汽含量變化的散點(diǎn)圖,當(dāng)水汽含量較低時(shí),特別是水汽含量小于1g·cm-2時(shí),平均偏差幾乎在0.5K以內(nèi),當(dāng)水汽含量增大時(shí),兩類LST的平均偏差整體也出現(xiàn)增大的趨勢(shì)。

上述對(duì)LST產(chǎn)品和JMS LST進(jìn)行了對(duì)比分析,為了進(jìn)一步對(duì)兩類LST的絕對(duì)精度進(jìn)行評(píng)估,下面基于SD、GB、CJZ三個(gè)站點(diǎn)實(shí)測(cè)數(shù)據(jù)對(duì)兩類LST進(jìn)行驗(yàn)證。驗(yàn)證前需要對(duì)地面站點(diǎn)的異質(zhì)性進(jìn)行分析。一般選擇站點(diǎn)周圍3×3像素進(jìn)行站點(diǎn)異質(zhì)性分析,由于Landsat 8 LST產(chǎn)品被采樣為30m存儲(chǔ) ,而實(shí)際分辨率為100m,因此本文采用站點(diǎn)周圍9×9像素進(jìn)行異質(zhì)性分析。統(tǒng)計(jì)站點(diǎn)周圍9×9像素溫度,并計(jì)算31景LST產(chǎn)品平均標(biāo)準(zhǔn)差作為異質(zhì)性分析依據(jù)。得到統(tǒng)計(jì)結(jié)果顯示: SD站為0.85K,GB站為0.412K,CJZ站為0.491K,三個(gè)站點(diǎn)的平均標(biāo)準(zhǔn)差均小于1K,滿足理想地表溫度驗(yàn)證點(diǎn)的不確定度應(yīng)在1K以內(nèi)的要求[31]。

圖2 不同地類的位置以及地表覆蓋Fig.2 The location and surface coverage of different land types

圖3 不同地類平均偏差對(duì)比以及與水汽含量關(guān)系Fig.3 Comparison of bias of different land types and their relationship with water vapor content

圖4 研究區(qū)內(nèi)兩類LST的平均偏差與水汽含量散點(diǎn)圖Fig.4 Scatter plot of the bias and water vapor content of the two types of LST in the study area

圖5為站點(diǎn)處兩類LST與站點(diǎn)實(shí)測(cè)數(shù)據(jù)的散點(diǎn)圖,從圖中可以看出兩類LST與地面測(cè)量數(shù)據(jù)具有很高的一致性,在每個(gè)站點(diǎn)處R2均高于0.949。圖6顯示了兩類LST在不同站點(diǎn)處與地面測(cè)量LST的偏差。與地面測(cè)量LST比較,兩類LST在SD站與GB站大部分日期偏差為負(fù),低估了地表溫度。對(duì)于濕地站,站點(diǎn)異質(zhì)性強(qiáng),地面測(cè)量數(shù)據(jù)作為像元尺度的驗(yàn)證值時(shí),可能會(huì)出現(xiàn)站點(diǎn)實(shí)際觀測(cè)值比像元尺度實(shí)際值高的情況,造成兩類LST低估。對(duì)于GB站,站點(diǎn)處溫度隨時(shí)間變化敏感,而站點(diǎn)觀測(cè)溫度是根據(jù)衛(wèi)星過境時(shí)間對(duì)10min內(nèi)的觀測(cè)值進(jìn)行線性插值得到,因此無法與衛(wèi)星過境時(shí)間時(shí)的溫度準(zhǔn)確對(duì)應(yīng),造成低估的情況。在GB站兩類LST RMSE在3.500K以上,明顯高于SD站與CJZ站,這可能與GB站驗(yàn)證數(shù)據(jù)較少有關(guān)。在CJZ站,兩類LST相比站點(diǎn)實(shí)測(cè)數(shù)據(jù)在大部分日期出現(xiàn)了高估的情況,因?yàn)樵撜军c(diǎn)下墊面為玉米,發(fā)射率變化隨季節(jié)變化敏感,高估的情況可能是溫度反演中發(fā)射率計(jì)算不準(zhǔn)確造成。表2為兩類LST在各個(gè)站點(diǎn)處的精度指標(biāo)。對(duì)比LST 產(chǎn)品與JMS LST,在SD站與GB站,LST 產(chǎn)品精度指標(biāo)Bias與RMSE均小于JMS LST, LST產(chǎn)品精度優(yōu)于JMS LST。在CJZ站,JMS LST(Bias=-0.051K,RMSE=2.028K)精度優(yōu)于LST產(chǎn)品(Bias=0.578K,RMSE=2.230K)。在所有站點(diǎn)處,LST產(chǎn)品 Bias為-0.714K,RMSE為2.865K, JMS LST Bias為-1.252K,RMSE為2.958K。整體來講, LST產(chǎn)品與JMS LST均有較高的精度,RMSE均在3K以內(nèi),LST產(chǎn)品的精度稍高于JMS LST。

從圖4可看出,兩類LST的偏差隨著水汽含量的不同而有明顯的變化,在水汽含量增大時(shí)偏差有增大趨勢(shì),基于所有站點(diǎn)對(duì)研究區(qū)內(nèi)水汽大于1g·cm-2和小于1g·cm-2兩種情況進(jìn)行兩類LST驗(yàn)證,如表3所示。水汽含量小于1g·cm-2時(shí),LST產(chǎn)品和JMS LST分別低估了1.617K和1.909K,同時(shí)LST產(chǎn)品RMSE小于JMS LST,LST產(chǎn)品精度稍高。在水汽大于1g·cm-2時(shí),LST產(chǎn)品高估了0.730K, JMS LST低估了0.201K,JMS LST具有更好的平均偏差,但是RMSE稍高于LST產(chǎn)品。整體來看,在本研究區(qū)的水汽范圍內(nèi),基于站點(diǎn)的驗(yàn)證結(jié)果兩類LST精度相似。由于研究區(qū)內(nèi)水汽含量整體較低,所有日期水汽含量均小于3g·cm-2,因此無法對(duì)高水汽含量情況下兩類LST的精度進(jìn)行評(píng)估。

圖5 各站點(diǎn)反演溫度與驗(yàn)證溫度散點(diǎn)圖Fig.5 Scatter plots of retrieved temperature and verified temperature at each site

圖6 各站點(diǎn)反演溫度與驗(yàn)證溫度偏差Fig.6 The bias of the retrieved temperature and the verified temperature at each site

表2 基于各站點(diǎn)的精度評(píng)估結(jié)果

表3 不同水汽含量情況下所有站點(diǎn)精度評(píng)估結(jié)果

5 結(jié)論

Landsat 8 LST產(chǎn)品具有重要的應(yīng)用價(jià)值,其精度評(píng)估亦需要展開。本文以黑河流域中游為研究區(qū),基于JMS 反演結(jié)果和地面站測(cè)量數(shù)據(jù)對(duì)2013-2016研究區(qū)內(nèi)31景清晰無云的Landsat 8 LST產(chǎn)品進(jìn)行了評(píng)估與分析,可以得出以下結(jié)論:

1)Landsat 8 LST產(chǎn)品與JMS反演結(jié)果相似,大部分日期偏差在1K以內(nèi),整體上反演溫度Landsat 8 LST產(chǎn)品稍高于JMS反演結(jié)果。

2)通過對(duì)比兩類LST之間的平均偏差,不同的地類對(duì)偏差大小影響較小,水汽含量的大小是影響兩類LST偏差的主要因素。 當(dāng)水汽含量小于1g·cm-2時(shí),大部分偏差在0.5K以內(nèi),水汽含量大于1g·cm-2,偏差有增大的趨勢(shì),最大偏差達(dá)到了2.5K。

3)基于站點(diǎn)數(shù)據(jù)驗(yàn)證,兩類LST整體上低估了地表溫度,RMSE均在3K以內(nèi),Landsat 8 LST產(chǎn)品的評(píng)價(jià)指標(biāo)(Bias=-0.714K,RMSE=2.865K)優(yōu)于JMS LST(Bias=-1.252K,RMSE=2.958K),LST 產(chǎn)品的精度稍高。

本研究區(qū)內(nèi)水汽含量較低,均小于3g·cm-2,因此缺乏高水汽含量情況下的驗(yàn)證與評(píng)估,本文的驗(yàn)證精度不夠全面,今后需要在高水汽含量地區(qū)進(jìn)行實(shí)驗(yàn)。