于夢(mèng)馨 劉 波，2，3

（1.湖南農(nóng)業(yè)大學(xué)信息科學(xué)技術(shù)學(xué)院 長(zhǎng)沙 410128）（2.邵陽學(xué)院湘西南農(nóng)村信息化服務(wù)湖南省重點(diǎn)實(shí)驗(yàn)室 邵陽 422000）（3.湖南省農(nóng)村農(nóng)業(yè)信息化工程技術(shù)研究中心 長(zhǎng)沙 410128）

1 引言

地表溫度反演研究是熱紅外遙感領(lǐng)域的重要研究方向，隨著衛(wèi)星遙感技術(shù)的發(fā)展、熱紅外波段應(yīng)用的拓寬和新的反演算法的出現(xiàn)，熱紅外遙感地表溫度反演將為全球能量變化研究、生態(tài)環(huán)境研究和數(shù)值氣象預(yù)報(bào)研究提供更準(zhǔn)確的地面溫度信息。在諸多的衛(wèi)星對(duì)地觀測(cè)系統(tǒng)中，Landsat系列衛(wèi)星數(shù)據(jù)具有長(zhǎng)期連續(xù)，較高空間分辨率，以及免費(fèi)使用等眾多優(yōu)勢(shì)，成為全球應(yīng)用最為廣泛的遙感數(shù)據(jù)之一。2013年2月11日，其最新一代Landsat8衛(wèi)星發(fā)射成功［1］，其攜帶的熱紅外傳感器TIRS（Thermal Infrared Sensor）具有兩個(gè)熱紅外波段，分別 是 第 10（10.6～11.9μm） 波 段 和 第 11（11.5～12.51μm）波段，這為地表溫度反演提供了新的數(shù)據(jù)。而基于Landsat8的地表溫度反演的實(shí)際應(yīng)用研究還相對(duì)較少，就這一問題，本文以長(zhǎng)沙市為研究區(qū)域，利用Landsat8數(shù)據(jù)，采用兩種劈窗算法［2～9］和覃志豪單窗算法［10～11］分別對(duì)研究區(qū)進(jìn)行地表溫度的反演計(jì)算，并結(jié)合MODIS MYD11A1地表溫度產(chǎn)品，對(duì)三種典型算法的反演結(jié)果進(jìn)行分析對(duì)比，旨在為L(zhǎng)andsat8數(shù)據(jù)的地表溫度反演應(yīng)用提供參考。

2 三種典型算法介紹

本文選擇以下三種典型算法：分別是Jua-n C.Jimênez-Mu?oz劈窗算法、Offer Rozenstein劈窗算法和覃志豪單窗算法。

1）Juan C.Jimênez-Mu?oz劈窗算法

Juan C.Jimenez-Munoz等（2014）提出的適用于Landsat8數(shù)據(jù)的劈窗算法（Split-Window Alorithm，簡(jiǎn)稱為SWJ），該算法計(jì)算較為簡(jiǎn)單，其計(jì)算模型如下

式中T10、T11為第10、11波段的亮度溫度（K）；ω為大氣水汽含量（g·cm-2）；ε為第10、11波段的平均地表比輻射率，即 ε=0.5(ε10+ε11)，Δε為第10、11波段的比輻射率差值，即 Δε=(ε10-ε11)；相關(guān)系數(shù)C0=-0.268，C1=1.378，C2=-0.183，C3=-54.30，C4=-2.238，C5=-129.20，C6=-16.4。

2）Offer Rozenstein劈窗算法

Rozenstein等（2014）提出的適用于Landsat8數(shù)據(jù)的劈窗算法（Split-Window Alorithm，簡(jiǎn)稱為SWR），其基本計(jì)算公式如下

式中，T10，T11為第10、11波段的亮度溫度（K）；A0、A1、A2是該算法的反演系數(shù)，它們由以下公式計(jì)算得到：

ai和bi分別是TIRS10、11波段根據(jù)不同溫度范圍確定的回歸系數(shù)，如表1所示。為了反演結(jié)果的精確性，需要根據(jù)研究區(qū)所在的地理位置以及影像數(shù)據(jù)選取的時(shí)間，估算出大致的溫度范圍，然后選取合適的系數(shù)。式中C、D是由地表比輻射率ε和大氣透射率τ所確定的參數(shù)，計(jì)算公式如下

表1 不同溫度范圍內(nèi)TIRS的反演回歸系數(shù)

3）覃志豪單窗算法

覃志豪單窗算法（Mono-window Algorithm，簡(jiǎn)稱為MW），是一種基于TM熱紅外波段的地表溫度反演算法，公式為

式中，T10是亮度溫度（K），ε是地表比輻射率，τ是大氣透過率，Ta是大氣平均作用溫度（K）。a、b為系數(shù)，可由表1查找。綜合上述三種算法可以看出，每種算法都需要三個(gè)參數(shù)，亮度溫度、大氣透過率和地表比輻射率，而單窗算法還需要計(jì)算大氣平均作用溫度。

3 長(zhǎng)沙市地表溫度反演實(shí)例

3.1 研究方法

本文以Landsat8影像數(shù)據(jù)為基礎(chǔ)，在ENVI和GIS的支持下，對(duì)影像數(shù)據(jù)進(jìn)行處理、分析、統(tǒng)計(jì)等操作，實(shí)現(xiàn)長(zhǎng)沙市地表溫度分布情況，并對(duì)三種算法進(jìn)行分析對(duì)比，以期為L(zhǎng)andsat8數(shù)據(jù)的地表溫度反演應(yīng)用提供參考。技術(shù)路線如圖1所示。

3.2 數(shù)據(jù)源及預(yù)處理

1）數(shù)據(jù)源

本文選用的數(shù)據(jù)來源于美國(guó)地質(zhì)調(diào)查局相關(guān)網(wǎng)站，是2016年7月23日的長(zhǎng)沙市Landsat8遙感數(shù)據(jù)，圖幅號(hào)為L(zhǎng)C81230402016205。該影像無云，質(zhì)量較好。Landsat8衛(wèi)星有兩個(gè)傳感器TIRS和OLI。陸地成像儀OLI具有8個(gè)多光譜波段和1個(gè)全色波段，分辨率分別為30m和15m。熱紅外傳感器TIRS具有兩個(gè)熱紅外波段，分辨率為100m，是有史以來最先進(jìn)，性能最好的TIRS，可以監(jiān)測(cè)陸地和水資源消耗情況。

圖1 地表溫度反演流程圖

2）數(shù)據(jù)預(yù)處理

本文對(duì)熱紅外影像和多光譜影像進(jìn)行一些預(yù)處理。

（1）輻射定標(biāo)：對(duì)于Landsat8 TIRS的兩個(gè)熱紅外波段，利用輻射定標(biāo)系數(shù)將其像元灰度值（DN）轉(zhuǎn)化為輻射亮度值，輻射定標(biāo)的公式如下

ML為增益參數(shù)，AL為偏移參數(shù)，兩個(gè)參數(shù)都可以從影像元數(shù)據(jù)中文件中獲取，Qcal是像元灰度值，也已知。ML對(duì)應(yīng)元數(shù)據(jù)中的RANDIANCE_MULT_BAND_X，AL對(duì)應(yīng)元數(shù)據(jù)中的RANDIANCE_ADD_BAND_X。

（2）反射率定標(biāo)：對(duì)于Landsat8 OLI多光譜波段，利用反射率定標(biāo)系數(shù)將其像元灰度值（DN）轉(zhuǎn)換為對(duì)應(yīng)的反射率值，公式如下

其中增益值Mρ和偏移值A(chǔ)ρ都可以從影像元數(shù)據(jù)文件中獲取，Qcal是像元灰度值，也已知。Mρ對(duì)應(yīng)元數(shù)據(jù)中的REFLECTANCE_MULT_BAND_X，Aρ對(duì)應(yīng)元數(shù)據(jù)REFLECTANCE_ADD_BAND_X。接著按太陽高度角進(jìn)行驗(yàn)證，公式如下

其中θSE為太陽高度角，可直接從影像元數(shù)據(jù)文件中獲取，語句為SUN_ELEVATION后的數(shù)值。

3.3 亮度溫度反演

根據(jù)普朗克定律將Landsat8第10和11波段的輻射亮度值轉(zhuǎn)化為星上亮度溫度，其計(jì)算公式如下

式中，T為星上亮度溫度，單位為K；Lλ為輻射亮度值；K1和 K2參數(shù)可從影像頭文件中獲取。通常T可以轉(zhuǎn)換為攝氏度溫度T′來表示：T′=T-273.15。下圖所分別為第10波段和11波段的星上亮度溫度。星上亮度溫度不能夠準(zhǔn)確地反應(yīng)地表溫度值，所以還需要進(jìn)一步的反演。

圖2 第10波段亮度溫度

圖3 第11波段亮度溫度

3.4 地表比輻射率的估算

地表比輻射率主要取決于地表的物質(zhì)結(jié)構(gòu)，是反演地表溫度的重要參數(shù)之一。為了得到更精確的地表比輻射率，本文選用覃志豪提出的地表比輻射率估算方法。其認(rèn)為地球表面不同區(qū)域的地表結(jié)構(gòu)雖然復(fù)雜，但可以大體認(rèn)為由三種類型構(gòu)成：水面、建筑物、和自然表面。其中自然表面可以看成是由植被和裸土組成的混合像元。為了確定研究區(qū)各類地物的比輻射率，首先利用Landsat8多光譜波段對(duì)地表進(jìn)行監(jiān)督分類，將研究區(qū)土地覆蓋類型分為水面、建筑物、和自然表面。然后分別計(jì)算比輻射率。自然表面的比輻射率由式（13）得出，建筑物的比輻射率由式（14）得出，水體的比輻射率較高，所以需要單獨(dú)提取計(jì)算。

注：如果用公式算出的ε＞εv，則取ε=εv

式中：i為L(zhǎng)andsat8熱紅外的第10和11波段。εiv、εis、εim分別為i波段處植被、裸土和建筑物的比輻射率。根據(jù)ASTER提供的常用地物比輻射率光譜庫(kù)（來源：http：//speclib.jpl.nasa.gov）和 Landsat8的第10和11波段區(qū)間，對(duì)第10波段的水體、植被、裸土和建筑物的比輻射率分別取0.9908、0.9867、0.9676、0.9648；對(duì)第11波段的水體、植被、裸土建筑物的比輻射率分別取0.9902、0.9899、0.9779、0.9751。其中PV為植被覆蓋度，用以下公式計(jì)算：

式中，NDVI為歸一化植被指數(shù)，NDVIs為完全裸土或無植被覆蓋區(qū)域的NDVI值，NDVIv則代表完全被植被所覆蓋的像元的NDVI值，即純植被像元NDVI值。本文取 NDVIv=0.7，NDVIs=0.05，當(dāng)NDVI≥NDVIv時(shí)，PV=1，認(rèn)為該地區(qū)完全被植被覆蓋；當(dāng)NDVI≤NDVIs時(shí)，PV=0，認(rèn)為該地區(qū)完全裸露。

RV，RS，Rm分別為植被、裸土和建筑表面的溫度比率，不僅與溫度變化有關(guān)，還與植被覆蓋度有關(guān)計(jì)算公式如下：

dε為熱輻射校正項(xiàng)，可以根據(jù)經(jīng)驗(yàn)公式進(jìn)行計(jì)算：

3.5 大氣透過率的估算

在地表溫度反演過程中，水汽含量是估算大氣透過率需要考慮的主要因素。本文的大氣水汽含量是通過楊景梅等［12］提出的經(jīng)驗(yàn)?zāi)Ｐ陀傻孛嫠麎汗浪愕玫?。本文使?016年7月23日Landsat8數(shù)據(jù)，上下景影像過境時(shí)間分別為10點(diǎn)56分和10點(diǎn)57分。從中國(guó)氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)查詢到衛(wèi)星過境當(dāng)天11：00的水汽壓為28.4h/Pa，估算水汽含量為3.4625g/cm2。前后相差3min，可認(rèn)為這兩個(gè)時(shí)刻的大氣水汽含量近似相等。大氣透過率與大氣水汽含量的關(guān)系主要是通過大氣模擬來確定。本文使用覃志豪等人根據(jù)大氣輻射傳輸軟件LOWTRAN模擬得到的大氣透過率估算方程（表2）進(jìn)行估算，采用中緯度夏季的計(jì)算方程來估算大氣透過率。

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

3.6 大氣平均作用溫度的估算

大氣平均作用溫度可由近地面氣溫T0估算，本文選擇中緯度夏季大氣的估算方程進(jìn)行計(jì)算，Ta為大氣平均作用溫度，計(jì)算公式如下：

4 反演結(jié)果對(duì)比分析

利用研究區(qū)的Landsat8數(shù)據(jù)，結(jié)合上述兩種劈窗算法和覃志豪單窗算法反演得到長(zhǎng)沙市地表溫度的反演結(jié)果圖，然后進(jìn)行密度分割，如下圖所示。從長(zhǎng)沙區(qū)域的統(tǒng)計(jì)來看，SWJ，SWR，MW得到的地表溫度反演結(jié)果的空間分布情況基本相同，只是溫度值稍有差異。將反演結(jié)果與之前進(jìn)行決策樹分類得到的土地利用類型圖進(jìn)行對(duì)比，可以得出，研究區(qū)的低溫區(qū)主要集中在湘江流域及研究區(qū)分布的其他水體，由于水的比熱容比較大，接受太陽輻射后升溫較慢，所以在研究區(qū)形成低溫區(qū)；中溫度區(qū)域，主要是2011年由望城縣劃為望城區(qū)的絕大部分區(qū)域，該區(qū)域植被覆蓋度較高；高溫度區(qū)域，主要是長(zhǎng)沙市主體城區(qū)部分，由于城市的下墊面比自然地表太陽吸收率高以及人工熱源等因素，溫度較高。

表3 各種地表溫度反演算法的比較

如表3所示：兩種劈窗算法反演得到的地表溫度高于單窗算法反演得到的地表溫度。兩種劈窗算法的最低溫度分別為15.87℃和17.36℃，最高溫度分別為51.53℃和50.49℃，分別相差1.49k，1.04k。單窗算法的最低溫度為15.24℃，最高溫度為47.97℃，與兩種劈窗算法相比較，最低溫度分別相差0.63K，2.12K，最高溫度分別相差3.56K，1.52K，可以看出兩種劈窗算法的反演結(jié)果更為相近。此外，本文還采用當(dāng)天的MYD11A1地表溫度產(chǎn)品進(jìn)行相對(duì)驗(yàn)證，由于研究區(qū)市中心影像的缺失，最終得到的溫度結(jié)果為16.57℃～35.89℃，由溫度的空間分布情況可知，市中心溫度最高，缺失影像的溫度沒有辦法進(jìn)行驗(yàn)證，但三種反演算法得到的最低溫度與MYD11A1溫度產(chǎn)品的最低溫度值分別相差0.7K，-0.79K，1.33K，與兩種劈窗算法的溫度差別較小，與單窗算法溫度差別較大。

圖4 長(zhǎng)沙市SWJ 地表溫度反演結(jié)果

圖5 長(zhǎng)沙市SWR地表溫度反演結(jié)果

圖6 長(zhǎng)沙市MW地表溫度反演結(jié)果

通過對(duì)反演結(jié)果的驗(yàn)證分析，發(fā)現(xiàn)兩種劈窗算法的反演結(jié)果更為接近，且精度均較高，單窗算法的反演精度較低，對(duì)于Landsat8影像數(shù)據(jù)，兩種劈窗算法的實(shí)際精度優(yōu)于單窗算法。

5 結(jié)語

本文以長(zhǎng)沙市為研究區(qū)域，選用2016年7月23日Landsat8衛(wèi)星影像數(shù)據(jù)，采用3種典型算法對(duì)研究區(qū)進(jìn)行了地表溫度反演，并對(duì)這3種算法的反演結(jié)果進(jìn)行了對(duì)比分析及精度驗(yàn)證。得出以下幾點(diǎn)結(jié)論：

1）三種算法的反演結(jié)果的總體趨勢(shì)和空間分布情況大體相同，只是溫度值稍有差異。就本文所使用的數(shù)據(jù)可以看出，兩種劈窗算法反演得到的地表溫度高于單窗算法反演得到的地表溫度。

2）兩種劈窗算法的溫度反演結(jié)果更為相近，最低溫度相差1.49k，最高溫度相差1.04k。單窗算法的溫度反演結(jié)果與兩種劈窗算法相比較，最低溫度分別相差0.63K，2.12K，最高溫度分別相差3.56K，1.52K，誤差較大。與MYD11A1溫度產(chǎn)品相比，最低溫度值分別相差0.7K，-0.79K，1.33K，兩種劈窗算法的反演精度優(yōu)于單窗算法。

3）本文使用的MYD11A1溫度產(chǎn)品，市中心有部分影像的缺失，不能更好地進(jìn)行溫度反演結(jié)果的驗(yàn)證，這是本文的不足之處。