夏安全+齊建國+姜振飛+馬津

摘要：當(dāng)前地表溫度反演的遙感數(shù)據(jù)源多為Landsat TM/ETM+、MODIS數(shù)據(jù)，Landsat 8熱紅外數(shù)據(jù)的使用還不是很多，許多針對Landsat 8熱紅外數(shù)據(jù)的地表溫度反演算法雖然被提出，但是否能滿足不同試驗(yàn)區(qū)的精度要求還有待考究，同時(shí)，經(jīng)驗(yàn)?zāi)Ｐ偷氖褂每赡軙η蠼獾乇肀容椛渎屎痛髿馔高^率等地表溫度反演參數(shù)造成不同程度的影響。因此，本研究以山東省濟(jì)南市為研究區(qū)，Landsat 8數(shù)據(jù)為數(shù)據(jù)源，分別利用大氣校正法、JM_SC10算法、TIRS10_SC算法，結(jié)合分類回歸樹（classification and regression tree，簡稱CART）算法與中等光譜分辨率大氣透過率算法計(jì)算機(jī)模型（moderate spectral resolution atmospheric transmittance algorithm and computer model，簡稱MODTRAN），構(gòu)建適合濟(jì)南地區(qū)的溫度反演參數(shù)，實(shí)現(xiàn)地表溫度的反演，并以濟(jì)南市16個(gè)氣象站的溫度數(shù)據(jù)為基準(zhǔn)進(jìn)行精度驗(yàn)證。結(jié)果表明，3種算法反演溫度平均誤差為1.78 ℃，TIRS10_SC算法反演精度最高，其次是大氣校正法、JM_SC10算法。在一定誤差要求下，3種算法均可應(yīng)用于濟(jì)南地區(qū)的地表溫度反演。

關(guān)鍵詞：Landsat 8；地表溫度反演；單通道算法；大氣透過率；地表比輻射率

中圖分類號： S127 文獻(xiàn)標(biāo)志碼： A 文章編號：1002-1302（2017）20-0254-05

陸地表面溫度（land surface temperature，簡稱LST）是地表能量平衡和資源環(huán)境變化的重要參數(shù)。地面監(jiān)測站獲取溫度信息準(zhǔn)確，但不能獲取大范圍溫度值及其時(shí)空分布，熱紅外遙感的出現(xiàn)解決了這個(gè)難題。熱紅外遙感技術(shù)能獲取熱紅外波段的輻射能量，根據(jù)地表物體的發(fā)射率特性反演得到其溫度，實(shí)現(xiàn)大范圍的溫度信息獲取，因而熱紅外遙感在地表溫度反演方面有著重要的作用[1]。

目前地表溫度的反演算法主要有單通道算法、劈窗算法、多通道多角度算法。單通道算法是指只利用1個(gè)熱紅外通道反演地表溫度的方法，適用于幾乎所有的熱紅外波段；劈窗算法是指利用相鄰的2個(gè)熱紅外通道來進(jìn)行地表溫度反演的方法，是目前為止發(fā)展最為成熟的地表溫度反演算法；多通道算法還在發(fā)展之中，目前還沒有一個(gè)簡便可行的多通道算法可以用來進(jìn)行地表溫度反演。Landsat 8有2個(gè)傳感器：陸地成像儀（operational land imager，簡稱OLI）和熱紅外傳感器（thermal infrared sensor，簡稱TRIS），其中TRIS有2個(gè)熱紅外波段，這雖然為劈窗算法的使用提供了條件，但一直以來美國地質(zhì)調(diào)查局（United States Geological Survey，簡稱USGS）對TIRS11波段的定標(biāo)準(zhǔn)確性把握不準(zhǔn)，所以不鼓勵使用劈窗算法來反演Landsat 8數(shù)據(jù)的LST[2]。因而本研究基于Landsat 8第10波段使用3種單通道算法，包括大氣校正法[3]、JM_SC10 算法[4]、TIRS10_SC算法[5]進(jìn)行LST反演。

本研究以山東省濟(jì)南市為研究區(qū)，Landsat 8數(shù)據(jù)為數(shù)據(jù)源，根據(jù)濟(jì)南市實(shí)時(shí)的大氣水汽含量，選用合適的大氣透過率模型進(jìn)行大氣透過率的估算；計(jì)算歸一化植被指數(shù)（normalized difference vegetation index，簡稱NDVI）、改進(jìn)型歸一化水體指數(shù)（modified normalized difference water index，簡稱MNDWI）、壓縮數(shù)據(jù)維建筑用地指數(shù)（index-based build-up index，簡稱IBI），采用基于分類回歸樹（classification and regression tree，簡稱CART）算法的數(shù)據(jù)挖掘技術(shù)，獲取用于區(qū)分影像中不同類別地物的規(guī)則閾值，實(shí)現(xiàn)影像地物的精確分類；獲取每一類別純凈植被像元與純凈裸土像元的NDVI，摒棄植被覆蓋度計(jì)算的經(jīng)驗(yàn)?zāi)Ｐ?，?shí)現(xiàn)研究區(qū)植被覆蓋度及地表比輻射率的精確計(jì)算。在此基礎(chǔ)上，對3種算法反演的地表溫度進(jìn)行比較分析，確定精度最高的單通道地表溫度反演算法，為后期應(yīng)用Landsat 8數(shù)據(jù)進(jìn)行地表溫度的反演提供參考。

1 材料與方法

1.1 研究區(qū)概況

濟(jì)南市為山東省省會，位于山東省中西部，介于36°01′～37°32′N、116°11′～117°44′E。地處中緯度地帶，季風(fēng)明顯，四季分明，年平均氣溫為13.8 ℃，夏季平均氣溫為26.7 ℃，冬季平均氣溫在1.0 ℃左右。濟(jì)南三面環(huán)山，南依泰山，北跨黃河，地勢南高北低，落差達(dá)500 m，平均海拔高度為118 m，這種地勢構(gòu)造令水汽和熱空氣回流聚集且不易擴(kuò)散，一旦發(fā)生熱污染將很難消除，因此實(shí)時(shí)監(jiān)測大范圍的濟(jì)南市地表溫度是很有必要的。

1.2 數(shù)據(jù)處理

遙感數(shù)據(jù)選用2015年4月25號覆蓋研究區(qū)的Landsat 8 OLI、TIRS數(shù)據(jù)以及MODIS L1B 1KM數(shù)據(jù)。其中OLI數(shù)據(jù)用來計(jì)算研究區(qū)的地表比輻射率；TIRS數(shù)據(jù)用來反演LST；MODIS L1B 1KM數(shù)據(jù)用來計(jì)算研究區(qū)大氣水汽含量；地面實(shí)測數(shù)據(jù)由濟(jì)南市氣象站提供。

數(shù)據(jù)處理包括MODIS L1B 1KM數(shù)據(jù)、Landsat 8多光譜數(shù)據(jù)輻射定標(biāo)與大氣校正、Landsat 8熱紅外數(shù)據(jù)輻射定標(biāo)等，總體技術(shù)路線如圖1所示。

2 地表溫度反演算法

2.1 大氣校正法

衛(wèi)星傳感器接收到的熱紅外輻射亮度值Lλ由3部分組成：大氣向上輻射亮度L↑、地面的真實(shí)輻射亮度經(jīng)過大氣層之后到達(dá)衛(wèi)星傳感器的能量、大氣向下輻射能量L↓到達(dá)地面后反射的能量。衛(wèi)星傳感器接收到的熱紅外輻射亮度值Lλ的表達(dá)式（輻射傳輸方程）可寫為：endprint

式中：λ為波長，μm；ελ為地表比輻射率；Ts為地表真實(shí)溫度，K；B（Ts）為黑體熱輻射亮度，W/（m2·sr·μm）；τλ為大氣在波長λ處的透過率。變換（1）式得溫度為Ts的黑體在熱紅外波段的輻射亮度B（Ts）為：

Ts可以用普朗克公式的函數(shù)獲取：

對于TIRS 10，K1=774.89 W/（m2·sr·μm），K2=1 321.08 K，將式（3）中B（Ts）替換為Lλ即為亮度溫度Tλ的計(jì)算表達(dá)式。

2.2 Jiménez-Muoz改進(jìn)單通道算法

JM_SC10算法是由Jiménez-Muoz等于2014年在其原有單通道算法（single channel method，簡稱SC）的基礎(chǔ)上提出來的，增加了針對Landsat 8的大氣參數(shù)，JM_SC10算法如下：

式中：γ、δ是基于Planck函數(shù)的2個(gè)參數(shù)；bγ為常數(shù)，在TIRS10波段為1 324，Lλ表示熱輻射亮度值，W/（m2·sr·μm）；Tλ表示亮度溫度，K。ψ1、ψ2、ψ3與大氣水汽含量（g/cm2）ω有關(guān)，當(dāng)ω>3 g/cm2時(shí)，

2.3 TIRS10_SC算法

胡德勇等于2015年提出了TIRS10_SC算法[5]，該算法是在總結(jié)輻射傳輸方程和覃志豪單窗算法的基礎(chǔ)上，專門針對Landsat 8 TIRS傳感器開發(fā)的。

T10為TIRS10的亮溫，K；K2=1 321.08 K；Ta為大氣平均作用溫度，K；根據(jù)影像的獲取時(shí)間，采用中緯度夏季的大氣平均作用溫度估算方程進(jìn)行Ta的計(jì)算，估算方程如下：

式中：T0為近地表溫度（K），可根據(jù)當(dāng)?shù)貧庀筚Y料獲取。

3 大氣透過率τ與地表比輻射率ε計(jì)算

3.1 大氣透過率τ的計(jì)算

熱紅外波段大氣透過率主要取決于大氣水汽含量，可通過中等光譜分辨率大氣透過率算法計(jì)算機(jī)模型（moderate spectral resolution atmospheric transmittance algorithm and computer model，簡稱MODTRAN）模擬計(jì)算出大氣透過率與大氣水汽含量的關(guān)系，根據(jù)獲取的水汽含量信息來計(jì)算各像元大氣透過率。Landsat 8數(shù)據(jù)很難進(jìn)行大氣水汽含量反演，但是MODIS 數(shù)據(jù)卻可以。MODIS數(shù)據(jù)包含36個(gè)波段，其中17、18、19為大氣吸收波段，第2、第5波段為大氣窗口波段。毛克彪等研究發(fā)現(xiàn)，利用MODIS數(shù)據(jù)的第2、19波段可反演大氣水汽含量[6]，大氣水汽含量ω為：

式中：ρ19和ρ2 分別為MODIS數(shù)據(jù)的第19、第2波段的地表反射率；α、β為常數(shù)，分別為0.020 0、0.632 1。

不同研究區(qū)、不同水汽含量下，大氣透過率的估算模型是不同的。當(dāng)前針對Landsat 8數(shù)據(jù)的通用估算模型要求水汽含量在0.5～3.0 g/cm2，而研究區(qū)平均的水汽含量為 4.086 g/cm2，因此本研究區(qū)不能使用Landsat 8數(shù)據(jù)通用模型。鑒于Landsat 8數(shù)據(jù)的TIRS10與MODIS L1B數(shù)據(jù)的31波段中心波長和波寬相似（圖2），那么應(yīng)用于31波段的大氣透過率估算方程，基本也適用于Landsat 8數(shù)據(jù)的TIRS10。根據(jù)中緯度大氣剖面數(shù)據(jù)進(jìn)行模擬[7]，得到Band31在各種水汽含量下的通用大氣透過率估算方程（相關(guān)系數(shù)為 0.997 48）

通過式（12）實(shí)現(xiàn)TIRS10大氣透過率的計(jì)算，并得到該波段大氣透過率影像（圖3）。

3.2 地表比輻射率ε計(jì)算

地表比輻射率的計(jì)算一直是地表溫度反演的重點(diǎn)和難點(diǎn)。地表比輻射率別稱地表發(fā)射率，與地表的物質(zhì)結(jié)構(gòu)有關(guān)。要實(shí)現(xiàn)對地表比輻射率的精確計(jì)算，前提就是對影像地物進(jìn)行精確分類?；贑ART算法的數(shù)據(jù)挖掘技術(shù)可以獲取不同地物的規(guī)則閾值，大大提高決策樹分類效率和精度。本研究將影像地物分為4類，分別是水體、自然表面（耕地、林地、草地）、城鎮(zhèn)表面（建筑物、道路）、裸土，利用CART算法獲取的各類規(guī)則閾值，確定決策樹如圖4所示，分類結(jié)果影像如圖5所示。

根據(jù)ASTER光譜庫和Nichol的研究成果[8]，獲取各地類純凈像元在TIRS10波段的地表比輻射率，水體為0.996 83、建筑為0.964 885、植被為0.986 72、裸土為0.967 67。影像中雖然存在各地類的純凈像元，但也存在大量的混合像元，其中自然表面和裸土區(qū)域可以看作是植被與裸土的混合，城鎮(zhèn)表面可以看作建筑物與植被的混合，考慮到混合像元這種情況，覃志豪等提出以下模型[9]計(jì)算地表比輻射率：

式中：ε自、ε城、ε裸分別為自然地表比輻射率、城鎮(zhèn)地表比輻射率以及裸土地表比輻射率；RV、RB、RS分別為植被、建筑、裸土的溫度比率；εV、εS、εB分別為純凈植被、裸土、建筑的地表輻射率；鑒于研究區(qū)南北高差相距500 m，可依據(jù)植被的構(gòu)成比例簡單估計(jì)dε，經(jīng)驗(yàn)?zāi)Ｐ蚚10]為：

當(dāng)0≤PV≤0.5時(shí)，dε=0.003 796PV；

當(dāng)0.5< PV≤1時(shí)，dε=0.003 796（1- PV）

需要注意的是，如果應(yīng)用公式計(jì)算出的ε自大于εV，取ε自=εV；ε城>εB，取ε城=εB；ε裸>εS，取ε裸=εS。PV為植被覆蓋度，由改進(jìn)像元二分模型計(jì)算：

NDVIS為純凈裸土或者建筑像元的NDVI值；NDVIV為純凈植被像元的NDVI值。溫度比率Ri與植被覆蓋度PV相關(guān)，覃志豪等根據(jù)各地表類型的溫度差異進(jìn)行模擬，確定出植被、建筑、裸土表面的溫度比率[9]：

計(jì)算地表比輻射率需要參數(shù)多、運(yùn)算復(fù)雜，總體波段運(yùn)算表達(dá)式如下：

式中：b1、b6、b7、b9分別為植被、裸地、建筑、水體的掩膜影像；b2為植被覆蓋度影像；b3、b4、b8分別為RV、RS、RB影像；b5為dε影像。地表比輻射率影像如圖6所示。endprint

4 地表溫度反演

4.1 地表溫度反演結(jié)果

逐一計(jì)算各個(gè)參數(shù)之后，進(jìn)行3種算法的地表溫度反演，反演結(jié)果（取置信區(qū)間99.9%以去除異常值）如圖7～9所示。

4.2 精度驗(yàn)證

根據(jù)濟(jì)南市氣象網(wǎng)統(tǒng)計(jì)，2015年4月25日濟(jì)南市溫度大致在16～31 ℃之間，部分地區(qū)溫度達(dá)到36～40 ℃，與影像反演結(jié)果大致相同。結(jié)合反演影像得出高溫天氣主要發(fā)生在章丘市以及濟(jì)南市區(qū)部分地區(qū)。究其原因是濟(jì)南市重工業(yè)多集中在該地區(qū)，如三一重工、濟(jì)南重工、梟龍重工等，重工業(yè)產(chǎn)生的熱污染使章丘市與濟(jì)南城區(qū)的部分地區(qū)溫度居高不下，嚴(yán)重影響到人們的生活。

對3種算法反演的溫度數(shù)據(jù)進(jìn)行分類統(tǒng)計(jì)發(fā)現(xiàn)，耕地植被區(qū)域溫度相對適中，溫度在22～37 ℃之間，平均溫度為 29 ℃；水體區(qū)域溫度最低，溫度在16～33 ℃之間，其中黃河水域的溫度主要集中在17～21 ℃之間，各大水庫如鵲山水庫、臥虎山水庫、玉清湖水庫的溫度在16～19 ℃之間，溫度較高的水體區(qū)域多為零星分布的池塘以及海岸地帶；裸地和建筑區(qū)域溫度最高，在28～40 ℃之間，平均溫度為33 ℃。整體來看，3種算法對濟(jì)南市地表溫度的反演比較合理。

為精確驗(yàn)證反演結(jié)果的精度，本研究選取濟(jì)南市氣象站16個(gè)站點(diǎn)溫度數(shù)據(jù)作為實(shí)測數(shù)據(jù)，對3種算法反演的地表溫度進(jìn)行精度驗(yàn)證。各站點(diǎn)實(shí)測溫度與應(yīng)用3種算法反演的溫度數(shù)據(jù)對比結(jié)果如表1所示。

由表1可看出，3種單通道算法中TIRS10_SC算法反演地表溫度誤差最小，平均誤差為1.10 ℃，絕對平均誤差為 1.24 ℃；其次是大氣校正法，平均誤差為1.40 ℃，絕對平均誤差為1.45 ℃；最大誤差來自JM_SC10算法反演的地表溫度，平均誤差為2.85 ℃，絕對平均誤差為2.85 ℃。算法反演溫度與溫度站實(shí)測溫度不盡相同，其中大部分反演溫度略高于實(shí)測站溫度，分析原因主要有：（1）溫度站實(shí)測溫度的獲取時(shí)間為10：30，而影像的獲取時(shí)間為10：48，在晴朗天氣下，溫度是逐漸上升的。（2）對于算法中應(yīng)用到的大氣向上和向下的輻射亮度L↑、L↓，采用美國航空航天局（National Aeronautics and Space Administration，簡稱NASA）提供的計(jì)算模型。要提高模型計(jì)算的精度，輸入?yún)?shù)必須要有實(shí)時(shí)的氣壓、風(fēng)速、海拔高度、相對濕度等參數(shù)，本研究沒有這些數(shù)據(jù)，使用了模型的缺省值。（3）對JM_SC10算法，大氣水汽含量是影響該算法精度的一個(gè)關(guān)鍵參數(shù)，當(dāng)大氣水汽含量降低到2 g/cm2時(shí)，地表溫度的反演誤差會降低1.5～3.0 ℃[4]。試驗(yàn)中研究區(qū)平均大氣水汽含量為4.086 g/cm2，因此可以推斷，一旦大氣水汽含量降低到2 g/cm2，應(yīng)用JM_SC10算法反演溫度的誤差有可能至少降低1.5 ℃， 使該算法反演的地表溫度誤差大大降低。

盡管上述原因造成了一定的誤差，但3種算法反演的溫度誤差平均值為1.78 ℃，絕對誤差的平均值為1.85 ℃，反演結(jié)果相對理想。

5 討論

本研究選用3種地表溫度反演算法進(jìn)行地表溫度的反演，根據(jù)濟(jì)南地區(qū)的實(shí)際情況，構(gòu)建適合濟(jì)南地區(qū)的大氣透過率和地表比輻射率參數(shù)。本研究創(chuàng)新點(diǎn)有：（1）對大氣透過率的計(jì)算，考慮到各研究區(qū)不同大氣狀況、不同水汽含量的影響，選用與Landsat 8數(shù)據(jù)第10波段性質(zhì)相近的MODIS第31波段的大氣透過率估算模型，消除了含水量的影響，摒棄先前整景影像的單一大氣透過率計(jì)算，實(shí)現(xiàn)影像上每像元的大氣透過率求解。（2）基于CART算法的數(shù)據(jù)挖掘技術(shù)，結(jié)合多源影像，如NDVI影像、MNDWI影像、IBI影像，獲取用于區(qū)分影像中不同類別地物的規(guī)則閾值，實(shí)現(xiàn)影像地物的精確分類。在此基礎(chǔ)上，獲取每一類別純凈植被像元與純凈裸土像元的NDVI值，摒棄植被覆蓋度計(jì)算的經(jīng)驗(yàn)?zāi)Ｐ?，求解溫度比率，?shí)現(xiàn)研究區(qū)植被覆蓋度以及地表比輻射率的精確計(jì)算。

研究展望：（1）考慮如果能根據(jù)濟(jì)南地區(qū)的實(shí)際情況利用MODTRAN模型模擬出Landsat 8數(shù)據(jù)第10波段大氣水汽含量與大氣透過率的關(guān)系模型，對反演精度的提高會大有幫助；（2）獲取濟(jì)南市實(shí)時(shí)的氣壓、風(fēng)速、海拔高度、相對濕度等參數(shù)，實(shí)現(xiàn)L↑、L↓精確計(jì)算；（3）JM_SC10算法提出的前提是ω>3 g/cm2，考慮對JM_SC10算法的使用進(jìn)行改進(jìn)，在原有算法的基礎(chǔ)上 ±1.5 ℃，反演2幅溫度影像，根據(jù)部分實(shí)測溫度數(shù)據(jù)，確定出真實(shí)的溫度反演影像，提高算法的精度。

6 結(jié)束語

不管是通用的大氣校正法還是最近提出的地表溫度反演算法，在實(shí)現(xiàn)濟(jì)南市大氣透過率和地表比輻射率的精確計(jì)算后，各算法均對濟(jì)南市地表溫度的反演有很好的適用性，尤其是TIRS10_SC算法精度相對較高，而當(dāng)大氣水汽含量低于 2 g/cm2 時(shí)，JM_SC10算法精度也可大大提高。

參考文獻(xiàn)：

[1]宋 挺，段 崢，劉軍志，等. Landsat 8數(shù)據(jù)地表溫度反演算法對比[J]. 遙感學(xué)報(bào)，2015，19（3）：451-464.

[2]徐涵秋. 新型Landsat 8衛(wèi)星影像的反射率和地表溫度反演[J]. 地球物理學(xué)報(bào)，2015，58（3）：741-747.

[3]Qin Z H，Karnieli A，Berliner P. A mono-window algorithm for retrieving land surface temperature from Landsat TM data and its application to the Israel-Egypt border region[J]. lnternational Journal of Remote Sensing，2001，22（18）：3719-3746.

[4]Jiménez-Muoz J C，Sobrino J A，Skokovic D，et al. Land surface temperature retrieval methods from Landsat-8 thermal infrared sensor data[J]. IEEE Geoscience and Remote Sensing Letters，2014，11（10）：1840-1843.

[5]胡德勇，喬 琨，王興玲，等. 單窗算法結(jié)合Landsat 8熱紅外數(shù)據(jù)反演地表溫度[J]. 遙感學(xué)報(bào)，2015，19（6）：964-976.

[6]毛克彪，覃志豪，王建明，等. 針對MODIS數(shù)據(jù)的大氣水汽含量反演及31和32波段透過率計(jì)算[J]. 國土資源遙感，2005（1）：26-29.

[7]Mao K B，Qin Z H，Shi J. A practical split-window algorithm for retrieving land surface temperature from MODIS data[J]. International Journal of Remote Sensing，2005，26（15）：3181-3204.

[8]Nichol J. An emissivity modulation method for spatial enhancement of thermal satellite images in urban heat island analysis[J]. Photogrammetric Engineering and Remote Sensing，2009，75（5）：547-556.

[9]覃志豪，李文娟，徐 斌，等. 陸地衛(wèi)星TM 6波段范圍內(nèi)地表比輻射率的估計(jì)[J]. 國土資源遙感， 2004，16（3）：27-32.

[10]Sobrino J A，Jiménez-Munoz J C，Paolini L. Land surface temperature retrieval from Landsat TM 5[J]. Remote Sensing of Environment，2004，90（4）：434-440.endprint