Landsat-8影像地表溫度反演算法對比分析

朱江山， 楊麗萍*， 白宇興， 蘇志強， 孔金玲

(1.長安大學(xué)地質(zhì)工程與測繪學(xué)院， 西安 710054； 2.長安大學(xué)地球科學(xué)與資源學(xué)院， 西安 710054)

地表溫度(land surface temperature，LST)是區(qū)域和全球尺度地表物理過程的關(guān)鍵因子之一[1]，作為地表和大氣能量交換的直接驅(qū)動因素，在環(huán)境變化、水文循環(huán)和災(zāi)害監(jiān)測等多個研究領(lǐng)域得到廣泛應(yīng)用[2-4]。

Landsat系列衛(wèi)星數(shù)據(jù)是反演地表溫度的重要數(shù)據(jù)源之一。目前，針對該系列衛(wèi)星數(shù)據(jù)，前人已經(jīng)提出多種LST反演算法。其中，普適性單通道(generalized single-channel, GSC)算法是應(yīng)用最為廣泛的一種算法。2003年，Jiménez-Muoz等[5]首次提出GSC算法，該算法分析了大氣水汽含量(ω)、大氣上行輻射、大氣下行輻射和大氣透過率(τ)之間的關(guān)系，將大氣水汽含量與其余3個參量建立相關(guān)關(guān)系，在已知地表比輻射率的條件下，將GSC算法表達(dá)成以大氣水汽含量為變量的函數(shù)，由于該算法所需大氣校正參數(shù)較少，所以得到了廣泛的應(yīng)用。2009年，Jiménez-Muoz等[6]對原算法做了修正，使之不僅能適用于Landsat-5 TM數(shù)據(jù)，也可適用于Landsat-7 ETM+數(shù)據(jù)，從而擴大了該算法的適用范圍。2013年，Landsat-8衛(wèi)星成功發(fā)射，該衛(wèi)星攜帶有運行陸地成像儀(operational land imager，OLI)和熱紅外傳感器(thermal infrared sensor，TIRS)。與Landsat系列的其他衛(wèi)星相比，Landsat-8衛(wèi)星采用了量子學(xué)紅外光子檢測技術(shù)，增強了TIRS對熱紅外波段的敏感性。同時，其所攜帶的TIRS可以獲取2個長波熱紅外波段的數(shù)據(jù)，使得地表溫度與大氣溫度的分離易于實現(xiàn)。因此，基于Landsat-8數(shù)據(jù)反演地表溫度表現(xiàn)出較大優(yōu)勢。2014年，Jiménez-Muoz等[7]進一步提出了適用于Landsat-8熱紅外波段的GSC算法。與此同時，中外學(xué)者也相繼提出了多種針對Landsat-8數(shù)據(jù)的劈窗(split-window, SW)算法。例如，2014年，Rozenstein等[8]提出了一種以大氣透過率 (atmospheric transmittance) 與地表發(fā)射率 (land surface emissivity, LSE) 為輸入?yún)?shù)的針對Landsat-8數(shù)據(jù)的SW算法，分析了輸入?yún)?shù)的多種選擇方案及其敏感性，得出該SW算法的均方根誤差(root mean square error，RMSE)為0.93 ℃。同年，Jiménez-Muoz等[7]也提出了針對Landsat-8數(shù)據(jù)的SW算法，采用大氣廓線數(shù)據(jù)庫與地表比輻射率光譜庫的模擬數(shù)據(jù)進行驗證，結(jié)果表明該算法的平均絕對誤差(mean absolute error，MAE)小于1.5 K，當(dāng)大氣水汽含量增加時有較好的精度。隨后，Yu等[9]、Jin等[10]和Du等[11]分別提出了不同的SW算法，并發(fā)現(xiàn)各種SW算法均具有相對較高的精度。綜上可見，針對Landsat-8數(shù)據(jù)的地表溫度反演，前期的GSC算法主要利用了其TIRS的第10波段，而SW算法則綜合利用了其Band 10 和 Band 11 兩個熱紅外波段。但是由于Landsat-8 TIRS Band 11的定標(biāo)不確定性問題，SW算法在Landsat-8 LST反演中的適用性及反演精度一直以來備受關(guān)注，已開展了很多關(guān)于GSC與SW算法反演精度的對比研究。宋挺等[12]比較了SW算法和GSC算法在無錫周邊地區(qū)的反演精度，發(fā)現(xiàn)SW算法精度較高，在0.7 K左右，GSC算法的精度在 1.4 K 左右。孫靜等[13]在南京市的研究表明， GSC算法的反演誤差在1 K以內(nèi)，結(jié)果較SW算法更為可靠。岳輝等[14]以神東礦區(qū)為例，基于MODIS地溫產(chǎn)品，比較了輻射傳輸方程法(radiative transfer equation，RTE)、基于影像(image-based method，IB)算法、單窗(Mono-Window，MW)算法與GSC算法的反演精度，結(jié)果表明GSC算法的反演精度為0.58 K，MW算法精度最高。段金饋[15]以北京市為研究區(qū)，對比了RTE法與GSC算法，得出GSC算法的精度要高于RTE法。同年，侯宇初等[16]在大興安嶺地區(qū)漠河市，基于野外實測數(shù)據(jù)對比了RTE法、SW算法與GSC算法的LST反演精度，得出這3種算法精度較高，均在1 K以內(nèi)。

以上工作為Landsat-8數(shù)據(jù)地表溫度反演研究走向?qū)嵱没於嘶A(chǔ)，然而由于驗證站點數(shù)量有限，加之研究時間和地表覆蓋的差異，已有結(jié)果尚存爭議，難以實現(xiàn)相互之間的綜合評價與對比分析，因而亟需基于長期、準(zhǔn)確、連續(xù)的實測數(shù)據(jù)開展進一步的評價研究，以提供更為系統(tǒng)和全面的分析結(jié)果。與此同時，針對TIRS Band 11的不確定性，USGS對TIRS波段的定標(biāo)參數(shù)進行了不斷調(diào)整[17]。2017年，USGS提出了Landsat-8數(shù)據(jù)的雜散光校正算法 (stray light correction algorithm, SLCA)，以期改善數(shù)據(jù)質(zhì)量，因此，其熱紅外波段的精度也有待進一步評價。針對上述問題，基于NOAA地表輻射通量網(wǎng)(surface radiation budget network, SURFRAD) 的Fort peck(FPK)站點、Sioux falls(SXF)站點與Goodwin creek(GWN)站點，選取37景質(zhì)量較好的Landsat-8熱紅外影像，分別采用GSC算法與SW算法進行地表溫度反演，以站點觀測的氣象參數(shù)為基礎(chǔ)，以實測地溫數(shù)據(jù)作為驗證數(shù)據(jù)，對這2種算法的反演精度進行評價與對比分析。

1 研究數(shù)據(jù)

采用37景Landsat-8影像，獲取時間為2018年1月1日—2019年6月30日，云量均低于10%，影像來源于https://earthexplorer.usgs.gov/。驗證站點采用SURFRAD網(wǎng)的FPK站點、SXF站點與GWN站點，采用與成像時間同步的站點實測數(shù)據(jù)進行驗證，數(shù)據(jù)源自ftp://aftp.cmdl.noaa.gov/data/radiation/surfrad/。1995年，SURFRAD站點在NOAA的支持下開始運行。它提供上行輻射通量、下行輻射通量、相對濕度和近地表氣溫等氣象參數(shù)[18-19]。測站的設(shè)計、選址、數(shù)據(jù)處理方法、儀器校準(zhǔn)原理和高質(zhì)量的操作程序確保SURFRAD站點可為衛(wèi)星數(shù)據(jù)驗證分析、氣候變化等研究提供長期、準(zhǔn)確、連續(xù)的地面輻射實測數(shù)據(jù)[20]。站點信息如表1所示。

表1 驗證站點信息

2 研究方法

2.1 GSC算法

Ts=γ[(φ1Lsen+φ2)/ε+φ3]+δ

(1)

(2)

(3)

式中：ε為地表比輻射率；γ為由式(2)確定的參數(shù)；δ為由式(3)確定的參數(shù)；bγ為常量，對于TIRS Band 10、Band 11，bγ分別為1 324 K、1 199 K；Lsen為影像像元對應(yīng)的光譜輻射值(W·m-2·sr-1·μm-1)；Tsen表示由熱紅外波段轉(zhuǎn)換后的亮溫值；φ1、φ2、φ3分別為大氣功能參數(shù)，其原始推導(dǎo)公式為

φ1=1/τ

(4)

φ2=-L↓-L↑

(5)

φ3=L↓

(6)

式中：τ為大氣透過率；L↓為下行輻射通量，W/m2，L↑為上行輻射通量，W/m2。

(7)

式(7)中：ω為大氣水汽含量，g/cm2。

2.2 SW算法

Ts=Ti+c1(Ti-Tj)+c2(Ti-Tj)2+c0+

(c3+c4ω)(1-ε)+(c5+c6ω)Δε

(8)

式(8)中：Ts為反演得到的LST(K)；Ti和Tj分別為波段i與j的亮度溫度；ε為波段i與j的平均地表比輻射率；Δε為波段i與j的地表比輻射率之差；由模擬數(shù)據(jù)確定的系數(shù)c0～c6見表2。

表2 由模擬數(shù)據(jù)確定的系數(shù)c0～c6對應(yīng)數(shù)值

2.3 重要參數(shù)計算

2.3.1 LSE的計算

LSE是物體與黑體在同溫度、同波長下的輻射出射度的比值，是確定地表熱輻射的重要參數(shù)。LSE受水分含量、化學(xué)組成、地表結(jié)構(gòu)、粗糙度和觀測條件(例如波長、分辨率和觀測角度)等多種因素影響。

考慮到各個站點的地表覆蓋類型情況，研究采用歸一化差值植被指數(shù)(normalized difference vegetation index，NDVI)閾值法計算LSE。即當(dāng)NDVI≤0時，認(rèn)為像元全部由水體覆蓋；當(dāng)00.7時，認(rèn)為像元完全由植被覆蓋[22]。對于混合像元時，LSE的計算公式[23]為

εimix=εivpv+εin(1-pv)+(1-εin)(1-pv)Fεiv

(9)

(10)

式中：εimix為Band 10和Band 11的混合像元的LSE；εiv為相應(yīng)波段植被的LSE；εin為相應(yīng)波段裸土的LSE；Pv為植被覆蓋度；F取決于地表的幾何分布，值為0.55；NDVImin和NDVImax分別為裸土和植被全覆蓋時的歸一化植被指數(shù)。

采用宋挺等[12]根據(jù)ASTER提供的常用地物光譜庫和Landsat-8的熱紅外波段波譜響應(yīng)函數(shù)確定的Band 10與Band 11分別對應(yīng)的水體、植被和裸土的LSE計算各個站點的LSE(表3)。

表3 Landsat-8 Band 10、Band 11對應(yīng)的水體、植被和 裸土的LSE

2.3.2 大氣透過率的計算

大氣透過率對地表熱輻射在大氣中的傳導(dǎo)具有重要影響，大氣水汽含量的變化直接影響大氣透過率。根據(jù)SURFRAD站點實測氣象參數(shù)計算大氣水汽含量，公式如下[24]：

e*=(1.007+3.46×10-6P)×

(11)

e=e*RH

(12)

式中：e*為飽和水汽壓，hPa；P為站點壓強，hPa；T為近地表溫度，℃；RH為相對濕度，%；e為實際水汽壓，hPa。

運用Cole提出[25]的大氣水汽含量與實際水汽壓之間的關(guān)系式計算大氣水汽含量

(13)

參考Rozenstein等[8]基于MODTRAN 4.0大氣模擬程序得到的大氣水汽含量與大氣透過率之間的關(guān)系式計算大氣透過率(表4)。

表4 大氣水汽含量與大氣透過率的關(guān)系式

2.3.3 實測溫度的計算

基于SURFRAD的站點實測氣象參數(shù)，采用式(14)計算得到實測溫度[26]：

(14)

式(14)中：R↑與R↓分別為上行、下行熱紅外輻射通量(W·m-2)；σ為斯蒂芬波爾茲曼常數(shù)，取值為5.670 51×10-8(W/m2)/K4；ε為地表寬波段輻射率，計算公式為[27]

ε=0.212 2ε29+0.385 9ε31+0.402 9ε32

(15)

式(15)中：ε29、ε31、ε32分別為MODIS數(shù)據(jù)波段29、波段31和波段32的地表窄波段輻射率，該數(shù)據(jù)從MODIS地溫產(chǎn)品(MOD11C1、MOD11A1)直接獲得。

現(xiàn)運用中心最接近站點的像元值代表算法反演得到的LST，采用NOAA 地球系統(tǒng)研究實驗室(earth system research laboratories，ESRL)提供的與SURFRAD站點影像獲取時間相同(相同分鐘)的氣象參數(shù)計算實測溫度。

3 結(jié)果與分析

3.1 反演結(jié)果

基于FPK站點、GWN站點和SXF站點，運用GSC算法與SW算法反演得到地表溫度，并結(jié)合站點實測數(shù)據(jù)進行算法精度驗證，圖1為LST反演結(jié)果。圖1中，LST_G為實測溫度，LST_GSC為GSC算法反演結(jié)果，LST_SW為SW算法反演結(jié)果，圖1(a)、圖1(b)與圖1(c)分別為FPK站點、GWN站點與SXF站點的反演結(jié)果。

圖1 LST反演結(jié)果Fig.1 The retrieved LST

3.2 對比分析

3.2.1 站點對比分析

為直觀對比算法精度，分別計算了算法的MAE與RMSE。由于衛(wèi)星數(shù)據(jù)反演得到的地表溫度精度具有季節(jié)性變化的特點[28]，加之每一站點地表結(jié)構(gòu)與類型的多樣性，首先對每個站點、每個季節(jié)進行算法精度的對比分析，結(jié)果如圖2所示。

由圖2(a)與圖2(d)可見，F(xiàn)PK站點整體反演精度在1～4 K，RMSE幾乎均小于4 K。夏季的反演精度最低，MAE在2～4 K，RMSE較大；冬季的反演精度最高，在1 K左右。與此同時，春季與夏季，GSC算法的反演精度要高于SW算法的精度。春季，GSC算法比SW算法誤差小0.35 K左右，夏季，GSC算法比SW算法誤差小1.2 K左右。數(shù)據(jù)量的大小影響精度評定，由于秋季的Landsat-8影像數(shù)據(jù)量較少，最后的評定結(jié)果可能略欠說服力。

由圖2(b)與圖2(e)可見，GWN站點整體反演精度在1～2 K之間，RMSE幾乎均小于2 K，精度較高。其中，春季與夏季，GSC算法在該站點的反演精度比在FPK站點高0.5 K左右；SW算法在該站點的精度較FPK站點的更高，這可能是由于海拔的影響。FPK站點海拔634 m，GWN站點海拔僅 98 m?？梢?，海拔越高，地表溫度反演誤差越大。該站點秋季與冬季Landsat-8影像分別為1景與3景，數(shù)據(jù)量較少，因此，采用春季與夏季的統(tǒng)計結(jié)果進行精度分析更具說服力。

圖2 站點精度對比Fig.2 Accuracy comparison within the validation sites

由圖2(c)與圖2(f)可見，在SXF站點，冬季2種算法的反演精度均小于1 K，并且RMSE小于1.6 K，因此，反演精度以及可信度都較高。在春季、夏季與秋季，GSC算法與SW算法的MAE值均小于 2 K。

綜上，由于夏季大氣水汽含量較高，大氣透過率下降，使得3個站點GSC算法與SW算法的反演精度均較低，但在FPK站點的精度最低，可能的原因是FPK站點海拔最高，海拔高度影響地表溫度算法反演精度。同時，由于該站點為高密度雜草覆蓋，夏季草地長勢較為茂盛，由此可能引起較大的地表比輻射率誤差，進而導(dǎo)致LST的估算誤差。

3.2.2 算法對比分析

為了整體評價2種典型地表溫度反演算法的精度，對所有站點地表溫度反演結(jié)果進行統(tǒng)計分析，結(jié)果如表5所示。

由表5可見，對于GSC算法，在春季與冬季，該算法的反演精度在1 K左右，夏季與秋季則為1～2 K。春季與冬季的反演精度要高于夏季與秋季的反演精度。出現(xiàn)該結(jié)果的原因可能是大氣水汽含量的影響。春季與冬季，水汽含量較??；夏季與秋季，水汽含量較大。由于GSC算法中，3個大氣函數(shù)均與大氣水汽含量有關(guān)，當(dāng)水汽含量較小時，對反演結(jié)果的影響較小；反之，當(dāng)水汽含量較大時，對反演結(jié)果的影響較大。因此，水汽含量較小的春季與冬季的反演精度要高于水汽含量較大的夏季與秋季的反演精度，即GSC算法在春季與冬季的反演精度要好于夏季與秋季的反演精度。

表5 2種算法精度統(tǒng)計

對于SW算法，整體而言，該算法的反演精度小于1.5 K，RMSE為1～2.5 K。夏季的反演精度最低，冬季的反演精度較高。究其原因可能是大氣水汽含量與地表比輻射率的影響。首先，在夏季，水汽含量較大，冬季則較小。大氣水汽含量的增加會導(dǎo)致大氣透過率的下降，這是SW算法地表溫度反演中的1個重要誤差來源。其次，在夏季，當(dāng)站點地表為高密植被覆蓋時，衛(wèi)星觀測到的為植被的莖與冠，較難觀測到地表，而實測得到的為地表溫度，因此夏季地表溫度反演結(jié)果誤差稍大。在FPK站點與SXF站點地表覆蓋類型分別為高密草地與農(nóng)田，夏季的反演精度稍低較為合理。

由統(tǒng)計分析可知，2種算法反演精度對比而言，在春季，GSC算法的MAE比SW算法的MAE低約0.53 K，GSC算法的RMSE比SW算法低約0.37 K，GSC算法精度較高；在夏季，GSC算法的MAE比SW算法高約0.34 K，SW算法精度較高；在秋季，GSC算法的MAE比SW算法高約0.24 K，GSC算法的RMSE比SW算法高約0.24 K，SW算法精度較高；在冬季，GSC算法的MAE比SW算法低約0.14 K，GSC算法的RMSE比SW算法低約0.13 K，GSC算法精度較高。

由此可見，春季與冬季，GSC算法的反演精度明顯高于SW算法?？赡艿脑蚴荊SC算法運用1個熱紅外波段進行地表溫度反演，該算法需要輸入的參數(shù)較少，不確定性因素較小。而SW算法運用2個熱紅外波段進行地表溫度反演，輸入的參數(shù)較多，不確定性因素增加，因此，GSC算法反演精度較高。反之，在夏季與秋季，SW算法的反演精度要高于GSC算法，可能的原因為大氣水汽含量較高時，SLCA算法的效果得到了較好的體現(xiàn)，其對2個熱紅外波段重新標(biāo)定的貢獻(xiàn)使得結(jié)合2個熱紅外波段反演地表溫度的效果要好于僅采用單一的熱紅外波段?？傊?，MAE是從準(zhǔn)確性方面評價算法的反演精度，2種算法的MAE在2 K以內(nèi)，準(zhǔn)確性較高。RMSE則從精確性方面來評價反演精度，2種算法的RMSE都在2.5 K以內(nèi)，精確性較高。

3.2.3 相關(guān)性分析

將2種算法的地表溫度反演結(jié)果與通過各個站點氣象參數(shù)計算得到的實測地表溫度進行相關(guān)性分析，結(jié)果如圖3和圖4所示。由圖可知，2種算法的反演結(jié)果與實測溫度具有較高的一致性，比較其相關(guān)系數(shù)，R2分別為0.990 7與0.992 2，均與實測溫度具有強烈的相關(guān)性，SW算法與實測地表溫度的相關(guān)性更強，說明由于SLCA的引入，較好地改善了Landsat-8熱紅外波段的數(shù)據(jù)質(zhì)量，使得SW算法在Landsat-8數(shù)據(jù)LST反演中同樣具有較強的適用性。綜合2種算法的精度對比得出，春季與冬季，GSC算法的反演精度較高，而夏季與秋季，SW算法的反演精度較高。

圖3 GSC算法反演值與實測溫度的相關(guān)性Fig.3 The correlation between GSC algorithm retrieved LST and in-situ LST

圖4 SW算法反演值與實測溫度的相關(guān)性Fig.4 The correlation between SW algorithm retrieved LST and in-situ LST

4 結(jié)論

基于37景不同季節(jié)、不同地表覆蓋類型及不同海拔地區(qū)的Landsat-8數(shù)據(jù)，采用SURFRAD的3個觀測站點，對比了Landsat-8數(shù)據(jù)地表溫度反演中的2種典型算法，即GSC算法和SW算法的反演精度，主要結(jié)論如下。

(1)GSC算法與SW算法的反演精度呈現(xiàn)出一定的季節(jié)依賴性。夏季的反演精度要低于春季與冬季的反演精度。由于夏季水汽含量較高，大氣透過率下降，導(dǎo)致所使用的2種算法夏季誤差稍大，冬季反演誤差較小。同時，地表覆蓋類型對反演結(jié)果也有一定的影響。

(2)春季與冬季，GSC算法的精度高于SW算法精度；夏季與秋季，SW算法的精度要高于GSC算法的精度。其中，春季，GSC算法的精度比SW算法精度高約0.53 K，且GSC算法的RMSE比SW算法的RMSE低約0.37 K；冬季，GSC算法的精度比SW算法精度高約0.14 K， 且GSC算法的RMSE比SW算法的RMSE低約0.13 K；夏季，GSC算法的精度比SW算法精度低約0.35 K；秋季，GSC算法的精度比SW算法精度低約0.25 K， 且GSC算法的RMSE比SW算法的RMSE高約0.24 K。在大氣水汽含量較低時，輸入?yún)?shù)的多少可能影響著地表溫度反演精度。GSC算法輸入?yún)?shù)較少，不確定因素較少，春季與冬季該算法精度較高。在大氣水汽含量較高時，可能是由于SLCA算法對2個熱紅外波段重新標(biāo)定的貢獻(xiàn)使得結(jié)合2個熱紅外波段反演地表溫度的效果要好于僅采用單一的熱紅外波段。

(3)GSC算法和SW算法與實測溫度的相關(guān)系數(shù)R2分別為0.990 7與0.992 2，與實測溫度具有強烈相關(guān)性，SW算法與實測地表溫度的相關(guān)性更強，表明由于SLCA的引入，較好地改善了Landsat-8熱紅外波段的數(shù)據(jù)質(zhì)量，使得SW算法在Landsat-8數(shù)據(jù)LST反演中同樣具有較強的適用性。總體而言，春季與冬季，GSC算法的反演精度較高，而夏季與秋季，SW算法的反演精度較高。

綜上可見，GSC算法與SW算法的反演精度雖然具有一定的季節(jié)依賴性，但2種算法的反演結(jié)果均具有較高的可信度。但是如何改善輸入?yún)?shù)如大氣水汽含量和地表比輻射率的不確定性是進一步提高地表溫度反演精度的關(guān)鍵問題，今后尚需深入探討。