基于Landsat 8的地表溫度反演研究

郭亞茹，張繼賢，杜漫飛

(1.遼寧工程技術(shù)大學(xué)，測繪與地理科學(xué)學(xué)院，遼寧 阜新 123000； 2.中國測繪科學(xué)研究院，北京 100830；3.國家測繪產(chǎn)品質(zhì)量檢驗(yàn)測試中心，北京 100830)

基于Landsat 8的地表溫度反演研究

郭亞茹1，2*，張繼賢3，杜漫飛1

(1.遼寧工程技術(shù)大學(xué)，測繪與地理科學(xué)學(xué)院，遼寧 阜新 123000； 2.中國測繪科學(xué)研究院，北京 100830；3.國家測繪產(chǎn)品質(zhì)量檢驗(yàn)測試中心，北京 100830)

為研究Landsat 8衛(wèi)星 TIRS傳感器的第10通道和第11通道對地表溫度反演的不同表現(xiàn)，以北京市為研究區(qū)，將單窗算法和普適性單通道算法分別應(yīng)用到TIRS傳感器的第10、11通道，并將這兩種算法的反演結(jié)果與綜合運(yùn)用兩個(gè)通道的劈窗算法反演結(jié)果進(jìn)行比較，結(jié)果表明：①三種算法反演的結(jié)果總體趨勢較為一致，中心高，邊緣低；②兩種單通道算法的反演精度普適性單通道算法優(yōu)于單窗算法；③第10通道的反演結(jié)果要優(yōu)于第11通道。

地表溫度反演；單窗算法；普適性單通道算法；劈窗算法；地表比輻射率

1 引 言

地表溫度與光合作用、近地氣溫、土壤溫度、風(fēng)形成、蒸散發(fā)、火災(zāi)危險(xiǎn)等都有較為緊密的關(guān)系，是地表能量平衡的重要參數(shù)，也是資源環(huán)境動(dòng)態(tài)變化的主要影響因素，地表溫度反演已經(jīng)成為遙感地學(xué)分析的一個(gè)重要研究領(lǐng)域。

隨著遙感技術(shù)的快速發(fā)展，利用熱紅外信息來反演地表溫度等參數(shù)已越來越受到重視且日漸發(fā)展成熟。2001年，覃志豪等根據(jù)地表熱輻射傳導(dǎo)方程，推導(dǎo)出簡單易行并且精度較高的單窗算法，均取得了較好的結(jié)果[1]；2003年，Jimenez-Munoz和Sobrino提出了適用于任何傳感器的普適性單通道算法，并利用TM6數(shù)據(jù)對地表溫度進(jìn)行了反演得到了較高精度的結(jié)果[2]。隨后出現(xiàn)了以各種單通道算法為對比的研究，孟憲紅等基于輻射傳輸方程算法、單窗算法和普適性單通道算法，利用TM6數(shù)據(jù)對甘肅金塔地區(qū)的地表溫度進(jìn)行反演[3]。近年來劈窗算法[4]以其綜合考慮雙通道的優(yōu)勢成為目前獲取地表溫度的主要方法。

本文以單通道算法分別對兩個(gè)熱紅外波段反演，定性地分析了三種方法反演結(jié)果，并利用北京氣象局提供的實(shí)測數(shù)據(jù)定量地驗(yàn)證了結(jié)果精度。

2 算法原理

2.1 劈窗算法[5]

劈窗算法在地表輻射傳導(dǎo)方程的基礎(chǔ)上，利用 10 μm～13 μm大氣窗口內(nèi)兩個(gè)相鄰熱紅外通道(一般為 10.5 μm～ 11.5 μm和 11.5 μm～12.5 μm)所攜帶信息來間接計(jì)算地表溫度。推導(dǎo)公式如下：

Ts=A0+A1T10-A2T11

(1)

式中：Ts即為地表溫度(K)；Ti表示熱通道i的亮度溫度(K)；A0，A1和A2表示系數(shù)，由大氣透過率和地表比輻射率等因子確定。其推導(dǎo)公式如下：

A0=a10E1-a11E2

(2)

A1=1+A+b10E1

(3)

A2=A+b11E2

(4)

E1=D11(1-C10-D10)/E0

(5)

E2=D10(1-C11-D11)/E0

(6)

A=D10/E0

(7)

E0=D11C10-D10C11

(8)

式中：

Ci=εi×τi

(9)

Di=[1-τi][1+(1-εi)×τi]

(10)

式中，εi為i通道的地表比輻射率；τi為i通道大氣透過率。

對于方程中的系數(shù)a10，a11，b10，b11根據(jù)如表1獲取。

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

2.2 單窗算法[1]

覃志豪等在地表熱輻射傳導(dǎo)方程以及一系列假設(shè)的基礎(chǔ)上，避免了對于實(shí)時(shí)探空數(shù)據(jù)的依賴，建立了適用于僅有1個(gè)熱紅外通道的反演方法——單窗算法。由于本次試驗(yàn)數(shù)據(jù)源為Landsat 8，該數(shù)據(jù)具有兩個(gè)熱通道，所以推導(dǎo)公式修改如下推導(dǎo)出單窗算法：

Tis={a(1-Ci-Di)+[b(1-Ci-Di)+Ci+Di]Ti-DiTa}/Ci

(11)

其中，Ts是地表溫度(K)；a，b為回歸系數(shù)；Ci，Di為中間變量，求解同上；Ta是大氣平均作用溫度(K)，在標(biāo)準(zhǔn)大氣狀態(tài)下(天氣晴朗、沒有旋渦作用)，其與近地面溫度(T0)存在如下線性關(guān)系：

美國1976平均大氣

Ta=25.9396+0.88045T0

(12)

熱帶平均大氣(北緯15°，年平均)

Ta=17.9769+0.91715T0

(13)

中緯度冬季平均大氣(北緯45°，冬季)

Ta=19.2704+0.91118T0

(14)

2.3 普適性單通道算法[2]

普適性單通道算法由JIMENEZ-MUNOZ等于2003年提出，該方法所需的參數(shù)較少，僅需通道有效波長、大氣水汽含量和傳感器輻射亮度3個(gè)參數(shù)，操作簡便且結(jié)果較為精確，該法主要優(yōu)點(diǎn)是可以使用相同方程和系數(shù)應(yīng)用于不同熱傳感器。推導(dǎo)公式如下：

(15)

(16)

δi=-γ×Li+Ti

(17)

λ為有效波長；εi為i通道的地表比輻射率；Li表示輻射亮度；Ti表示亮度溫度；γ和δ均為中間變量，根據(jù)Sobrino給出的大氣水分含量與傳感器響應(yīng)函數(shù)之間的關(guān)系，可以確定大氣參數(shù)ψ1，ψ2，ψ3：

ψ1=0.14714ω2-0.15583ω+1.1234

(18)

ψ2=-1.1836ω2-0.37607ω-0.52894

(19)

ψ3=-0.04554ω2+1.8719ω-0.39071

(20)

2.4 參數(shù)估計(jì)

2.4.1 亮度溫度

亮度溫度是在普朗克定律下傳感器獲取的波長與通過均衡測量的輻射率積分而獲得的定向溫度，計(jì)算公式如下：

Ti=Ki2/ln(1+Ki1/Li)

(21)

式中，Ki1和Ki2為普朗克常量，Li為輻射亮度。

2.4.2 地表比輻射率[6]

地表比輻射率是物體與黑體在同溫度、同波長下的輻射出射度的比值。研究表明：在地表溫度反演過程中，比輻射率有0.01的誤差，就會(huì)導(dǎo)致反演出的地表溫度產(chǎn)生1℃～2℃的誤差，并且發(fā)射率糾正引起的誤差將是所有大氣糾正所引起誤差之和[7]。由此可以看出，地表比輻射率的精度必須得到保證。目前應(yīng)用較為廣泛且考慮較為全面的是覃志豪等所提出的算法，本文將使用其法，利用可見光和近紅外波段來估計(jì)研究區(qū)的地表比輻射率：TIRS傳感器的分辨率為 100 m，從影像中我們可以將地表大體分為水體、城鎮(zhèn)和自然表面3種類型[8]，其中組成自然表面的像元可以看作是由不同比例的植被葉冠和裸土組成的混合像元，而城鎮(zhèn)像元也可看成是由建筑物和綠化植被組成的混合像元，計(jì)算過程如下：

(1)自然表面混合像元的地表比輻射率估算公式如下：

εi=PvRvεiv+(1-Pv)Rsεis+dε

(22)

式中，εiv，εis分別表示植被，裸土的比輻射率；Pv是植被占混合像元的比例：Pv=(NDVI-NDVImin)/(NDVImax-NDVImin)，其中NDVI為植被歸一化指數(shù)。Rv和Rs表示植被和裸土的溫度比率，公式如下：

(23)

dε表示混合像元之間的熱輻射相互作用量，估計(jì)式如下：

(24)

(2)城鎮(zhèn)混合像元的地表比輻射率計(jì)算公式：

εi=PvRvεiv+(1-Pv)Rmεim+dε

(25)

式中：Rm表示建筑表面的溫度比率；εim表示建筑表面的比輻射率。

3 實(shí)例驗(yàn)證與分析

3.1 數(shù)據(jù)概況

本文數(shù)據(jù)為2017年1月31日10：53的Landsat 8影像(由于影像已覆蓋研究區(qū)域的絕大部分，故僅選取了一幅)。研究區(qū)域?yàn)楸本┦械貐^(qū)(覆蓋率非100%)，該地平均海拔為 435 m，為典型的北溫帶半濕潤大陸性季風(fēng)氣候，冬季寒冷干燥，全年平均太陽輻射量為112～136×103k/cm2。兩個(gè)高值區(qū)分別分布在延慶盆地及密云縣西北部至懷柔東部一帶；低值區(qū)位于房山區(qū)的霞云嶺附近。

本文預(yù)處理由ENVI5.1軟件完成，主要包括以下幾步：

(1)輻射定標(biāo)：即將灰度值(digital number，DN)圖像轉(zhuǎn)換為輻射亮度圖像的處理過程。

(2)大氣校正：采用FLAASH大氣校正模型對Landsat 8 OLI傳感器的多光譜波段進(jìn)行大氣校正，以消除大氣和光照等因素對地物反射的影響，消除大氣分子和氣溶膠散射的影響，該過程也是反演地物真實(shí)反射率的過程。

(3)幾何校正：將幾何誤差控制在1個(gè)像元內(nèi)。

3.2 結(jié)果分析

利用ENVI5.1軟件中的band math工具，逐步實(shí)現(xiàn)上述理論方法，前述三種不同方法反演分別得到北京市地表溫度空間分布圖如圖1所示。

圖1北京市地表溫度反演結(jié)果

從三種算法的反演結(jié)果圖中可以看出，不同算法得到的北京市地區(qū)溫度空間分布圖變化趨勢基本一致，市中心范圍內(nèi)溫度高，向外逐漸降低，市中心及海淀區(qū)明顯高于郊區(qū)，這主要是由于市中心區(qū)人口密度較大，汽車尾氣排放量、生活中的能量消耗等因素。一些高耗型、高排放量的工廠主要分布在西部，這與圖像上呈現(xiàn)的溫度較高一致。城區(qū)水體附近溫度均低于零下且變化不大。

本文驗(yàn)證數(shù)據(jù)來自于北京氣象局提供的2017年1月31日10：55區(qū)局實(shí)況數(shù)據(jù)，在反演結(jié)果圖上以區(qū)局為單位隨機(jī)選點(diǎn)，對比情況如下。

北京市Landsat 8第10、11波段地表溫度反演結(jié)果驗(yàn)證表 表2

三種算法與實(shí)測溫度的比較 表3

表中，Ts表示實(shí)測溫度；T1，T2，T3，T4，T5分別表示單窗算法的第10、11波段、普適性單通道算法的10、11波段和劈窗算法的反演結(jié)果；△1，△2，△3，△4，△5分別表示單窗算法的第10、11波段、普適性單通道算法的10、11波段以及劈窗算法反演結(jié)果與實(shí)測溫度之差。

根據(jù)表格我們可以看到每一種算法中的第10、11波段反演結(jié)果具有良好的一致性，通過標(biāo)準(zhǔn)差的比較，我們可以看出第10波段的反演精度要高于第11波段。在表2中，通過與實(shí)測溫度的比較，在觀測區(qū)域內(nèi)，基于Landsat 8的3種地表溫度反演模型中，我們可以發(fā)現(xiàn)3種模型跟實(shí)測氣溫具有較好的一致性，不同觀測區(qū)域的實(shí)測溫度和劈窗算法反演結(jié)果都很接近，精度最高的是普適性單通道算法，其次是劈窗算法，單窗算法的精度較低。

4 結(jié) 論

本文針對Landsat 8衛(wèi)星的熱紅外數(shù)據(jù)特點(diǎn)，應(yīng)用3種模型反演了北京市地表溫度。得到如下結(jié)論：

(1)3種反演模型溫度空間分布趨勢大致相同，表現(xiàn)為市中心范圍內(nèi)溫度高，向外逐漸降低，水體的溫度變化較小。普適性單通道算法與地表溫度實(shí)際溫度最為接近，其次是劈窗算法，單窗算法的反演結(jié)果較差。本次試驗(yàn)結(jié)果顯示第10通道的反演結(jié)果要優(yōu)于第11通道的反演結(jié)果。

(2)本次試驗(yàn)總體來說，僅僅基于Landsat 8數(shù)據(jù)比較了劈窗算法、單窗算法和普適性單通道算法的精度，這為后續(xù)的與地表溫度相關(guān)的研究提供了較好的基礎(chǔ)。同時(shí)實(shí)測溫度與衛(wèi)星過境時(shí)實(shí)際地表溫度也是存在著差異的，故使得算法的絕對精度驗(yàn)證受到了一定的限制。

[1] 覃志豪，Zhang M，Karnieli A等. 用陸地衛(wèi)星TM6數(shù)據(jù)演算地表溫度的單窗算法[J]. 地理學(xué)報(bào)，2001，56(4)：456～466.

[2] Jiménez-Munoz JC，Sobrino J A. A generalized single-channel method for retrieving land surface temperature from remote sensing data[J]. Journal of Geophysical Research，2003，108(D22)：ACL(2-1)～(2-9).

[3] 孟憲紅，呂世華，張宇等. 使用LANDSAT 5 TM數(shù)據(jù)反演金塔地表溫度[J]. 高原氣象，2005，24(5)：721～726.

[4] GAO C X，TANG B H，WU H et al. A Generalized Split-Window Algorithm for Land Surface Temperature Estimation from MSG-2 /SEVIRI Data[J]. International Journal of Remote Sensing，2013，34(12) ：4182～4199.

[5] OFFER R，QIN Z H，YEVGENY D et al. Derivation of Land Surface Temperature for landsat 8 TIRS Using a Split Window Algorithm[J]. Sensors，2014( 14)：5768～5780.

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

[7] 柳菲，王新生，徐靜等. 基于NDVI閾值法反演地表溫度比輻射率的參數(shù)敏感性分析[J]. 遙感信息，2012，27(4)：3～12.

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

ResearchofRetrievalofLandSurfaceTemperaturefromLandsat8

Guo Yaru1，2，Zhang Jixian3，Du Manfei1
(1.School of Geomatics，Liaoning Technical University，F(xiàn)uxin 123000，China；2.Chinese Academy of Surveying and Mapping，Beijing 100830，China；3.National Quality Inspection and Testing Center for Surveying and Mapping Products，Beijing 100830，China)

In order to study the different performance of the 10th and 11th channels of the Landsat 8 TIRS sensor to the surface temperature inversion，this paper use the two single window algorithms and the split-window algorithm into the TIRS sensor. The inversion results of the single algorithms are compared with those of the split-window algorithm. The results show that：①The overall trend of the three algorithms is more consistent，showing the center is high and the edge is low；②The inversion accuracy of the generalized single-channel algorithm is superior to the accuracy of single - window algorithm；③The inversion results of the 10th channel are better than the 11th channel.

surface temperature inversion；single window algorithm；generalized single-channel method；split-window algorithm；surface emissivity

1672-8262(2017)06-74-05

TP722.5

A

2017—03—20

郭亞茹(1991—)，女，碩士研究生，研究方向：攝影測量與遙感。

國家基礎(chǔ)測繪科技計(jì)劃(2016KJ0103)；中國博士后科學(xué)基金資助項(xiàng)目(2016M591219)