基于Landsat-8數(shù)據(jù)和劈窗算法的地表溫度反演及城市熱島效應研究

宋挺， 段崢， 劉軍志，嚴飛，黃君，吳蔚

(1. 無錫市環(huán)境監(jiān)測中心站，江蘇 無錫 214121;2. 代爾夫特理工大學,代爾夫特,2628 CN,荷蘭)

0 引言

城市熱島效應是指城市地區(qū)整體或局部溫度高于周圍地區(qū)，溫度較高的城市地區(qū)被溫度較低的郊區(qū)所包圍或部分包圍的現(xiàn)象[1]。熱量平衡是城市熱島形成的能量基礎，城市化改變了城市下墊面的性質(zhì)和結(jié)構(gòu)，增加了人為熱，從而影響了城市熱量平衡各分量的變化[2]。城市熱島效應通常被視為一種城市公害，特別是中低緯度地區(qū)的夏季，熱島增溫效應加劇城市的酷熱，致使城市氣候舒適度變差，加重了空氣污染和資源消耗，影響了城市生態(tài)環(huán)境質(zhì)量[1]。

衛(wèi)星遙感的出現(xiàn)為城市熱島效應的研究提供了新的數(shù)據(jù)源，其中衛(wèi)星影像的熱紅外波段在城市熱島效應研究中應用尤為廣泛。已有學者利用熱紅外波段數(shù)據(jù)對城市熱島進行了研究，取得了不少成績[3-6]。目前，在城市熱島遙感研究中，地表溫度(Land Surface Temperature, LST)數(shù)據(jù)主要以遙感圖像灰度值、亮度值和遙感反演等方式獲取。地物灰度值、亮度溫度與地表實際溫度相差較大，不能有效支撐城市熱島的定量研究。而基于LST的反演算法則被證明是一種精度較高的地表溫度獲取方法，適用于城市熱島的定量研究。

在LST反演的眾多算法中，劈窗算法是一種精度較高的代表性算法[7-8]。該算法利用10～13 μm大氣窗口內(nèi)兩個相鄰熱紅外通道(一般為10.5～11.5 μm和11.5～12.5 μm)對大氣吸收作用的不同,通過兩個通道測量值的各種組合來剔除大氣的影響,進行大氣和地表比輻射率的修正。劈窗算法已成為目前由熱紅外遙感數(shù)據(jù)獲取地表溫度的主要方法。

劈窗算法最初主要是針對NOAA/AVHRR(The National Oceanic and Atmospheric Administration/Advanced Very High Resolution Radiometer)與MODIS(Moderate Resolution Imaging Spectroradiometer)等遙感產(chǎn)品而開發(fā)，經(jīng)過20多年的發(fā)展，目前公開發(fā)表的劈窗算法已將近20個。覃志豪等[9]對12種劈窗算法進行了詳細的討論。劈窗算法的一種表現(xiàn)形式如下:

Ts=A0+A1Ti-A2Tj

(1)

式中:Ts——地表溫度;Ti和Tj——熱通道i和j的亮度溫度;A0，A1和A2——參數(shù),由大氣透過率和地表比輻射率等因子確定,不同算法中它們的值各不相同，也可以改寫成如下形式:

Ts=Ti+A(Ti-Tj)+B

(2)

式中:Ts，Ti和Tj含義同上;系數(shù)A和B也由大氣透過率和地表比輻射率等因子確定。

NOAA/AVHRR與MODIS等遙感產(chǎn)品的空間分辨率較低，往往不能滿足城市熱島定量研究對精細LST空間分布信息的需求。因此利用劈窗算法和較高分辨率的遙感產(chǎn)品進行LST反演，將有力促進城市熱島的定量研究。2013年2月11號，美國宇航局NASA 成功發(fā)射了 Landsat-8 衛(wèi)星，為走過了40年輝煌歲月的 Landsat系列衛(wèi)星計劃重新注入新鮮血液。Landsat-8上攜帶有2個主要載荷：陸地成像儀(Operational Land Imager, OLI)和熱紅外傳感器(Thermal Infrared Sensor，TIRS)。Landsat-8數(shù)據(jù)的2個熱紅外波段(10和11波段)與MODIS熱紅外波段的31和32波段的波寬以及中心波長基本一致，具體對比細節(jié)如表1所示。

表1 Landsat-8(第10和11波段)與MODIS熱波段(31和32波段)數(shù)據(jù)對比

根據(jù)表1中的數(shù)據(jù)參數(shù)，基于劈窗算法和Landsat-8 TIRS數(shù)據(jù)進行LST反演是可行的。因此，筆者利用Landsat-8的雙通道的熱紅外數(shù)據(jù)，采用劈窗算法對無錫地區(qū)的地表溫度進行反演，并根據(jù)地表溫度對無錫市熱島的空間分布進行分析和研究。

1 數(shù)據(jù)源與工作流程

1.1 數(shù)據(jù)源

1.1.1 遙感數(shù)據(jù)

Landsat-8屬于NASA的Landsat系列，屬于近極地太陽同步軌道衛(wèi)星。筆者利用TIRS的10和11波段估算星上亮度溫度；利用OLI的3、4、5和6波段生成研究區(qū)的歸一化植被指數(shù)(Normalized Difference Vegetation Index，NDVI)、改進型歸一化水體指數(shù)(Modified Normalized Difference Water Index ，MNDWI)以及歸一化建筑指數(shù)(Normalized Difference Barren Index，NDBI)用來參與比輻射率的估算。TIRS與OLI要分別進行輻射定標，筆者還對OLI進行了FLAASH大氣校正，各項定標參數(shù)可從Landsat-8影像的頭文件里獲取。

使用2013年12月20日與2014年3月16日包含無錫地區(qū)的Landsat-8數(shù)據(jù)，圖幅號分別為LC81190382013344和LC81190382014075。由于缺少夏季的Landsat-8數(shù)據(jù)參與熱島分析，研究還收集了2013年8月12日的Landsat-7的ETM+數(shù)據(jù)(圖幅號為LE71190382013224)，對陸地表面溫度反演所使用的方法為單通道算法(SC)[10]。另外，對應時段的MODIS1B數(shù)據(jù)和MOD11-L2地表溫度產(chǎn)品被用于計算水汽含量和對Landsat-8反演的結(jié)果進行精度驗證。

1.1.2 實測數(shù)據(jù)

選取無錫市環(huán)境監(jiān)測中心站在無錫太湖水域設立的16個浮標點位(圖1)的水溫數(shù)據(jù)作為地面實測數(shù)據(jù)對反演的LST進行比較和驗證。在全部的16個站點，監(jiān)測頻率均為30 min/次。選取了與Landsat-8過境時間最為接近的實測站點數(shù)據(jù)用于遙感反演結(jié)果的驗證。

圖1 選取的太湖水域內(nèi)16個實測水溫水質(zhì)自動監(jiān)測點位分布

1.2 工作流程

首先利用Landsat-8的TIRS數(shù)據(jù)、OLI數(shù)據(jù)和MODIS1B數(shù)據(jù)分別進行亮溫、地表比輻射率和大氣水汽含量的計算，然后再以上述中間結(jié)果作為輸入，利用劈窗算法進行LST反演。在反演的地表溫度基礎上，利用實測水溫數(shù)據(jù)和MODIS溫度產(chǎn)品對反演結(jié)果進行驗證與對比分析，最后用得到的地表溫度反演結(jié)果進行熱島效應的分級與評價。工作流程如圖2所示。

圖2 Landsat-8地表溫度反演與熱島效應分級評價流程

2 針對Landsat-8影像數(shù)據(jù)的劈窗算法

2.1 地表熱輻射傳輸方程的構(gòu)建

地表熱輻射傳輸方程是描述熱輻射傳播通過介質(zhì)時與介質(zhì)發(fā)生相互作用(吸收、散射、發(fā)射等)而使熱輻射能按一定規(guī)律傳輸?shù)姆匠?。遙感器所接受到的熱輻射主要有由地表熱輻射經(jīng)大氣衰減后被遙感器接受的熱輻射(即被測目標本身的熱輻射)、大氣向上的熱輻射(大氣直接熱輻射)和大氣向下熱輻射(大氣向地面的熱輻射)經(jīng)地表反射后又被大氣衰減最終被遙感器接受的熱輻射3部分[11-12]，故MODIS遙感器所接收到的熱輻射計算公式如下:

Bi(Ti)=τi(θ)·[εi·Bi(Ts)+(1-εi)·Li↓]+Li↑

(3)

式中:Ti——Landsat-8通道i的亮度溫度；Ts——地表溫度;εi——地表比輻射率;τi(θ)——通道i在遙感器視角為θ下從地面到遙感器的大氣透射率;Bi(Ti)——遙感器所接收到的輻射強度；Bi(Ts)——在地表溫度為Ts時的黑體輻射強度；Li↑和Li↓——大氣向上和向下的輻射強度。

此外，由于用大氣向上平均輻射強度代替大氣向下輻射強度不會產(chǎn)生較大的地表溫度誤差[9]。

Li↑≈[(1-τi(θ)]·Bi(Ta)

(4a)

Li↓≈[(1-τi(θ)]·Bi(Ta↓)

(4b)

用Ta替代Ta↓對方程的計算不產(chǎn)生實質(zhì)性的影響，將Li↑與Li↓代入(3)式，簡化得到：

Bi(Ti)=εi·τi(θ)·Bi(Ts)+[1-τi(θ)]·[1+(1-εi)·τi(θ)]Bi(Ta)

(5)

式中:Ti和Ts含義同上；Ta——大氣向上平均作用溫度;εi和τi(θ)含義同上;Bi(Ti)和Bi(Ts)含義同上；Bi(Ta)——通道i的大氣向上平均輻射強度。

針對Landsat-8數(shù)據(jù)的10和11波段，則有：

B10(T10)=ε10·τ10(θ)·B10(Ts)+[1-τ10(θ)]·[1+(1-ε10)·τ10(θ)]B10(Ta)

(6a)

B11(T11)=ε11·τ11(θ)·B11(Ts)+[1-τ11(θ)]·[1+(1-ε11)·τ11(θ)]B11(Ta)

(6b)

式中:B10(T10)、B10(Ts)和B10(Ta)——波段10對應于亮度溫度、地表溫度和大氣向上平均溫度的輻射強度;B11(T11)、B11(Ts)和B11(Ta)——波段11對應于亮度溫度、地表溫度和大氣向上平均溫度的輻射強度;τ10和τ11——波段10和11的地表比輻射率;τ10(θ)和τ11(θ)——波段10和11從地面到遙感器的大氣透過率。

2.2 Planck輻射函數(shù)的線性展開

Planck方程是復雜的非線性函數(shù),要想從方程組(6)中直接求出地表溫度很難。而用其Taylor展開式的前兩項來表示它的近似值,是推導劈窗算法的通用做法[11-13]。依據(jù)Planck方程中輻射強度與溫度之間的線性關系,通過對Planck方程的泰勒展開并取前兩項得到下式：

Bi(Tj)=Bi(T)+(Tj-T)ΔBi(T)/ΔT=(Li+Tj-T)ΔBi(T)/ΔT

(7)

針對Landsat-8波段10或11數(shù)據(jù),若上式中i=10或11，j=10或11,則Tj表示波段10或11的亮度溫度;若j=s,則Tj表示所要求解的地表溫度;若j=a,則Tj表示大氣向上平均作用溫度;T是任意確定的一個亮度溫度;Bi(Tj)和Bi(T)分別是溫度為Tj和T時的第i波段的輻射強度。參數(shù)Li是輻射強度與其隨溫度變化率的比值。

針對Landsat-8波段10和11,建立下式的計算Li的表達式：

Li=Bi(T)/[ΔBi(T)/ΔT]

(8)

L10=a10+b10T10

(9a)

L11=a11+b11T11

(9b)

在式(7)中命令T=T10、i=10和j=10，則有:

B10(T10)=(L10+T10-T10)ΔB10(T10)/ΔT=L10B10(T10)/ΔT

(10a)

同理,針對波段10和11其他各個溫度,有如下的表達式:

B10(Ts)=(L10+Ts-T10)ΔB10(T10)/ΔT

(10b)

B10(Ta)=(L10+Ta-T10)ΔB10(T10)/ΔT

(10c)

B10(T11)=(L10+T11-T10)ΔB11(T10)/ΔT

(10d)

B10(Ts)=(L10+Ts-T10)ΔB11(T10)/ΔT

(10e)

B10(Ta)=(L10+Ta-T10)ΔB11(T10)/ΔTi

(10f)

2.3 劈窗算法的推導

文中的劈窗算法,是建立在輻射傳輸方程組(6)的基礎上。為了簡化起見,將方程組(6)中的系數(shù)定義如下:

Ci=εiτi(θ)

(11)

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

(12)

據(jù)此，將方程組(6)改寫成下式：

B10(T10)=C10B10(Ts)+D10B10(Ta)

(13a)

B11(T11)=C11B11(Ts)+D11B11(Ta)

(13b)

將式(10)代入式(13)中，再將求解Li的式(9)代入，兩式相減并消解方程，最終可得到如下的反演地表溫度的公式，它就是所要推導的劈窗算法：

Ts=A0+A1T10-A2T11

(14)

式中:Ts、T10、T11含義同上,單位K。A0、A1和A2——系數(shù)。表2給出了不同的溫度范圍時，系數(shù)ai和bi的值[14]?？梢钥闯?，溫度范圍減小，精度較高。為了獲得準確的地表反演溫度，最好是根據(jù)溫度范圍選擇系數(shù)。

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

最終公式推演如下表示：

A0=a10E1-a11E2

(15a)

A1=1+A+b10E1

(15b)

A2=A+b11E2

(15c)

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

(15d)

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

(15e)

A=D10/E0

(15f)

E0=D11C10-D10C11

(15g)

式中：C10、C11、D10和D11可由式(11)和式(12)得出。由式(14)和(15)可知,求解地表溫度的關鍵是確定亮度溫度、地表比輻射率、大氣透過率和遙感器視角。亮度溫度利用Planck方程可直接計算。Landsat-8的TIR最大天頂角僅為7.5°，為窄視場傳感器，可以不用考慮天頂觀測角的影響。因此下文接著討論反演地表溫度所需的2個參數(shù)(大氣透射率τ和地表比輻射率ε)的計算。

3 亮度溫度、地表比輻射率、大氣透過率的計算

3.1 亮度溫度

亮度溫度是指輻射出與觀測物體相等輻射能量的黑體的溫度。將影像DN值定標為熱輻射強度(Ii)之后,可用Planck函數(shù)求解出星上亮度溫度,計算公式如下：

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

(16)

Landsat-8在第10波段的中心波長λ10=10.9 μm；在第11波段的中心波長λ11=12.0 μm，根據(jù)Planck公式得出Landsat- 8的Ki1和Ki2常量。

對于第i=10波段,分別為K10 1= 774.89 W/(m2·sr·μm),K10 2=1 321.08 K。

對于第i=11波段,分別為K11 1= 480.89 W/(m2·2sr·μm),K11 2=1 201.14 K。

3.2 地表比輻射率的估計方法

地表比輻射率是物體與黑體在同溫度、同波長下的輻射出射度的比值,它受物體的表面狀態(tài)、介電常數(shù)、含水量、溫度、物體輻射能的波長、觀測角度等多種因素的影響。目前求地表比輻射率的方法主要有差值法、獨立溫度光譜指數(shù)法(TISI)和NDVI門檻值(NDVITHM)等方法[15-16]。

對于Landsat系列影像圖像，地表比輻射率主要取決于地表的物質(zhì)結(jié)構(gòu)和遙感的波段區(qū)間。地球表面不同區(qū)域的地表結(jié)構(gòu)雖然很復雜，但從衛(wèi)星像元的尺度來看，可以大體視作由4種類型構(gòu)成：水面、建筑、裸土和植被[17]。水面結(jié)構(gòu)簡單；建筑包括城市和農(nóng)莊，主要由道路、各種建筑和房屋組成；植被包括林地和農(nóng)田等。

3.2.1 典型地物的比輻射率

水體在熱波段范圍內(nèi)的比輻射率很高，接近于黑體。植被的比輻射率也很高，一般認為在0.981～0.983。HUMES等[18]測量結(jié)果表明，全覆蓋灌木葉冠的熱波段輻射比率為0.986左右。典型土壤的平均比輻射率為0.972 15左右。對于城鎮(zhèn)的各種建筑表面而言，在11.5～12.5 μm波段范圍內(nèi)，其比輻射率一般在0.960～0.980之間變動，取其中值，比輻射率大約在0.970左右，略小于自然土壤的地表比輻射率。

根據(jù)如上判定，確定landsat-8的10與11通道比輻射率值如下：εi w、εi v、εi s和εm分別是水體、植被、裸土和建筑在第i(i同上)波段的地表比輻射率。分別取ε10w=0.996 83,ε11w=0.992 54,ε10v=0.986 72,ε11v=0.989 90,ε10s=0.967 67，ε11s=0.977 90,ε10m=0.964 885,ε11m=0.975 115。

3.2.1.1 自然表面比輻射率的估計方法

實際上，組成自然表面的象元可以簡單看作是由不同比例的植被葉冠和裸土組成，即混合象元。對于面積較大的100%植被或裸土表面，可直接用這兩種類型的地表比輻射率來表示象元的比輻射率，因此，當Pv=1時，ε=εv；當Pv=0時，ε=εs。但是，通常很難有100%的植被覆蓋或裸土表面，因此，覃志豪等[17]通過下式來估計混合象元的地表比輻射率。

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

(17)

式中，Pv——植被覆蓋度，即植被的構(gòu)成比例:Pv=(NDVI-NDVIs)/(NDVIv-NDVIs)，其中NDVIv與NDVIs分別為植被和裸土的NDVI值，文中取NDVI的5%置信區(qū)間來定義。Rv和Rs——植被和裸土的溫度比率。在地表相對平整的情況下，一般可取dε=0；在地表高低相差較大的情況下，dε可以根據(jù)植被的構(gòu)成比例簡單估計。由于熱輻射相互作用在植被與裸土各占一半時達到最大，所以，覃志豪等[17]提出如下經(jīng)驗公式來估計dε：

當Pv<0.5時，dε=0.003 8Pv

當Pv>0.5時，dε=0.003 8(1-Pv)

當Pv=0.5時，dε最大，dε=0.001 9

值得指出的是，用公式算出的ε如果大于εv，則取ε=εv。

3.2.1.2 城鎮(zhèn)地表比輻射率的估計方法

城鎮(zhèn)在地表溫度反演中也很重要，尤其是在村莊較密集的地區(qū)或城市郊區(qū)。一般來說，城鎮(zhèn)像元的地表比輻射率與確定自然表面比輻射率的方法類似。由于城鎮(zhèn)主要是由各種建筑表面和分布其中的綠化植被所組成，因此有：

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

(18)

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

3.2.1.3 典型地物溫度比率估計

比較準確的地表比輻射率估計需要對典型地物的溫度比率進行估計。為了確定典型地物的溫度比率，根據(jù)主要地表類型的溫度差異進行模擬，得出Ri隨植被覆蓋度變化而變化的狀況。覃志豪等[17]提出用如下公式估計植被、裸土和建筑表面的溫度比率：

Rv=0.933 2+0.058 5Pv

(19a)

Rs=0.990 2+0.106 8Pv

(19b)

Rm=0.988 6+0.128 7P

(19c)

這一估計雖然沒有考慮溫度變化的影響，但基本上已經(jīng)能夠滿足地表溫度反演對地表比輻射率估計的要求。

3.2.2 地表構(gòu)成比例的確定

文中使用分類回歸決策樹(Classification And Regression Tree,CART)對地表進行分類，決策樹模型分類法結(jié)構(gòu)明確，能自動選取特征并融入影像以外的各種知識，是人工智能遙感信息的主要提取方法之一，廣泛用于地表分類的研究中。

結(jié)合地物光譜特征和光譜影像信息，建立分類規(guī)則，將地表分為植被、建筑物、裸土和水體4個部分。文中使用NDVI、MNDWI與NDBI 3種遙感信息指數(shù)對地表信息進行決策樹分類。

根據(jù)NDVI指數(shù)、MNDWI指數(shù)和NDBI指數(shù)結(jié)合人工判定閾值構(gòu)建無錫地區(qū)地表分類的決策樹模型。以2014年3月16日為例，根據(jù)10%置信度，NDVI閾值取0.65，MNDWI取0.17，NDBI取-0.05,構(gòu)建模型流程如圖3所示。

圖3 決策樹模型分類流程

分類后，對水體、植被、裸土、建筑物分別計算其地表比輻射率。其中，10波段水體的比輻射率ε10w=0.996 83,11波段水體的比輻射率ε11w=0.992 54；植被和裸土等自然表面的比輻射率根據(jù)公式(14)得出；建筑物的比輻射率根據(jù)公式(15)得出。

3.3 大氣透過率的計算

在地表反演過程中，水汽是估計大氣透過率的主要考慮因素。通常的做法是通過MODTRAN、6S和LOWTRAN等大氣模型軟件模擬大氣透過率與水汽含量的關系[19]。因為實時的大氣剖面資料很難獲得，所以這種模擬結(jié)果精度有時難以保證。

覃志豪等[20]提出了求算透過率的方法，即通過MODTRAN來計算大氣水汽含量和透過率的關系，然后通過大氣水汽含量來求算透過率。對于Landsat-8數(shù)據(jù)而言，大氣水汽含量很難從影像反演得到，文中使用與Landsat-8影像同一天的MODIS數(shù)據(jù)，根據(jù)毛克彪[21]的方法，可以用MODIS第2和第19波段來反演大氣水分含量，推算公式如下：

ω={[α-ln(ρ19/ρ2)]/β}2

(20)

式中:ω——大氣水汽含量，g/cm2；α和β——常量，α= 0.02，β=0.632 1；ρ19和ρ2——MODIS第19和2波段的地面反射率。

大氣透過率與大氣水汽含量的關系主要通過大氣模擬來確定。根據(jù)中緯度大氣剖面數(shù)據(jù)進行模擬[22]，得出大氣透過率估算方程：

τ31=2.897 98-1.883 66e-(ω/-21.227 04),R2=0.997 48

(21a)

τ32=-3.592 89+4.604 14e-(ω/-32.706 39),R2=0.996 85

(21b)

由于Landsat-8第10和11通道的中心波長與MODIS第31和32通道的中心波長基本相對應，所以式(21)推得的31通道與32通道的大氣透過率，基本適用于Landsat-8的第10和11通道的大氣透過率。

4 無錫市地表溫度反演結(jié)果分析

4.1 溫度反演的比較驗證

可根據(jù)實測值或其他標準值對一個反演算法進行對比驗證，文中利用無錫市環(huán)境監(jiān)測站在太湖水域建立的水質(zhì)自動監(jiān)測站的實測水溫數(shù)據(jù)進行比對，并結(jié)合美國宇航局提供的MODIS地表溫度數(shù)據(jù)集進行相對驗證。

文中選取16個太湖水域的水質(zhì)自動監(jiān)測站的水溫數(shù)據(jù)，監(jiān)測數(shù)據(jù)的時間精度可以達到30 min，完全可以滿足對遙感數(shù)據(jù)進行比對的數(shù)據(jù)同步的要求。

同時MOD11-L2級MODIS溫度產(chǎn)品也被用于驗證基于Landsat-8反演的LST結(jié)果的精度。該MODIS產(chǎn)品根據(jù)WAN等[23]提出的多通道反演計算，能同時反演地表溫度和比輻射率。在晴朗天氣條件下，MODIS溫度產(chǎn)品誤差在1K以下[24]。該MODIS溫度產(chǎn)品過境時間為當?shù)貢r間11:35，Landsat-8衛(wèi)星過境對應的當?shù)貢r間為10:33，兩者差距接近1 h，MODIS過境時間的實際地表溫度應稍高于Landsat-8過境時間的實際地表溫度。得到的水溫比對的結(jié)果如表3所示(其中，Landsat-8統(tǒng)計數(shù)據(jù)為空間分辨率100 m的象元，MODIS地溫產(chǎn)品空間分辨率為1 000 m)。

表3 基于Landsat-8反演得出的LST、MODIS地溫產(chǎn)品與實測數(shù)據(jù)比對結(jié)果 K

從表3可以看出，2014年3月16日Landsat- 8反演的對照點的平均水溫是284.93 K，MODIS溫度產(chǎn)品的平均水溫是286.09 K，而實測的平均水溫為284.96 K。從平均值來看，Landsat-8反演的水溫比MODIS溫度產(chǎn)品較實測的誤差更小，僅為0.03 K，而MODIS溫度產(chǎn)品較實測誤差有1.13 K。從絕對誤差來看，Landsat-8反演的水溫的平均絕對誤差為0.52 K，較MODIS溫度產(chǎn)品的1.13 K也明顯要小。2013年12月20日的結(jié)果也類似， Landsat-8反演的水溫較實測水溫的平均絕對誤差為 0.62 K，較MODIS溫度產(chǎn)品的1.29 K也明顯要小。

由于Landsat-8熱紅外數(shù)據(jù)空間分辨率是100 m，而MODIS地表溫度產(chǎn)品空間分辨率是1 km，所以Landsat-8熱紅外數(shù)據(jù)需重采樣到1 km，以便與MODIS數(shù)據(jù)保持一致。重采樣為1 km分辨率的Landsat-8反演結(jié)果和MODIS溫度產(chǎn)品的影像對比見圖4。

圖4 2014年3月16日Landsat-8地溫反演結(jié)果與MODIS地溫產(chǎn)品對比

文中采用最小值、最大值和平均值這些統(tǒng)計指標對反演結(jié)果進行評價。分別對Landsat-8地表反演結(jié)果與MODIS地溫產(chǎn)品設置0.5%的置信區(qū)間，以去除異常值。統(tǒng)計結(jié)果如表4所示。

表4 2014年3月16日Landsat-8劈窗算法反演結(jié)果與MODIS地溫產(chǎn)品對比統(tǒng)計 K

由于不可避免地存在一定幾何校正方面的誤差，統(tǒng)計區(qū)域范圍會有略微的偏差。從無錫區(qū)域的統(tǒng)計來看，利用劈窗算法對Landsat-8溫度數(shù)據(jù)反演值相對于MODIS地溫產(chǎn)品最小值誤差為-1.2 K，最大值誤差為2.7 K，平均值誤差為-0.1 K。雖然平均值誤差很小，但是由于MODIS的溫度產(chǎn)品過境時間為11:35，Landsat-8衛(wèi)星過境時間為10:33，兩者差距接近1 h，在這段時差里，溫度變化大概為1K左右，且水域溫度上文已做過比較，MODIS溫度產(chǎn)品的水溫比Landsat-8反演的水溫要高1 K左右。所以兩者陸地區(qū)域的實際平均誤差應該在1K左右，利用劈窗算法對Landsat-8陸地溫度數(shù)據(jù)反演值要略高于MODIS地溫產(chǎn)品，而水域溫度數(shù)據(jù)較MODIS溫度產(chǎn)品更接近于實測值，表明該方法具有一定的可信精度。

還需要指出的是，Landsat-8的2013年12月10日的數(shù)據(jù)，無錫區(qū)域基本無云，數(shù)據(jù)質(zhì)量較高，無錫區(qū)域反演的平均溫度為284.4 K。但是MODIS比Landsat-8晚過境1 h，在MODIS過境的時候，無錫陸地區(qū)域上空已經(jīng)有了薄云層，造成MODIS地溫數(shù)據(jù)較實際要低，與Landsat- 8反演數(shù)據(jù)相比，平均值要低2.3 K，再考慮到MODIS比Landsat-8晚過境1 h，且水溫數(shù)據(jù)比Landsat-8反演水溫又要略高，陸地表面溫度的實際差別可能要達到3 K以上，綜合MODIS產(chǎn)品陸域上空薄云層的影響與統(tǒng)計結(jié)果差異較大的原因，這一天MODIS地溫產(chǎn)品陸地表面溫度數(shù)據(jù)不建議用作比對。

4.2 敏感性分析

敏感性分析通常是先假定某一參數(shù)有一微小誤差，其他參量不變或者在指定的范圍內(nèi)變化，分析不同情況下由這一誤差帶來的最終結(jié)果的變化。通過下式可以分析參數(shù)誤差引起反演地表溫度的誤差[25-26]:

式中:△Ts——地表溫度的反演誤差；△x——參數(shù)x的估計誤差；Ts(x+△x)和Ts (x)——用(x+△x)和x進行反演得到的地表溫度。

文中以2014年3月16日的Landsat-8數(shù)據(jù)為例，進行敏感性分析。

4.2.1 比輻射率敏感性分析

首先對比輻射率的影響進行分析，結(jié)果表明：給定溫度，大氣越干燥，比輻射率影響越大；一定的大氣狀況，溫度越高，影響反演結(jié)果越大；同樣的比輻射率誤差，通道的波長越長，所引起的溫度誤差越大。

通過對地表比輻射率不同改變值的設定，同時其他參量不變，進行地表比輻射率的敏感性分析。將2014年3月16日第10與11波段的比輻射率都增加與減去一定的量，值分別為0.01，0.02，0.04和0.06。最終得到的結(jié)果如表5所示。

表5 地表比輻射率變化時，地表溫度反演結(jié)果

從表5可以看出，隨著比輻射率的降低，反演得到的LST基本是等量升高的，但升高的幅度略微增加。地表比輻射率每改變0.01，反演得到的LST改變0.66 K，劈窗算法對地表比輻射的敏感性較低。

4.2.2 水汽含量敏感性分析

由于大氣透過率是基于大氣中的水汽含量，預計大氣透過率的偏差會在波段10與波段11都會發(fā)生。因此，敏感性分析是對水汽含量作為模型的自變量。

通過對水汽含量不同濃度的設定，同時其他參量不變，進行水汽含量的敏感性分析。

將水汽含量設為0.5，1，1.5，2，2.5，3，4，5和6 g/cm2，對2014年3月16日的數(shù)據(jù)分別代入進行計算，最終得到的結(jié)果如表6和圖5所示。

表6 水汽含量為變化時，地表溫度反演結(jié)果

圖5 水汽含量變化時候反演LST趨勢

由表6和圖5可以看出：劈窗算法反演得到的地表溫度結(jié)果呈一元多次函數(shù)上升，在低水汽含量(ω<1 g/cm2)情況下，對水汽的敏感性較高，在中等水汽含量(1 g/cm2<ω<3 g/cm2)的情況下， 敏感性較低，在高水汽含量(ω>3 g/cm2)情況下，敏感性又略微升高，總體來說，敏感性不是很高，水汽含量每變化1 g/cm2，劈窗算法的地表溫度反演結(jié)果平均變化0.4 K左右。

由于地物的比輻射率誤差較小，一般在1%以內(nèi)，而水汽含量使用MODIS的實時數(shù)據(jù)進行反演，誤差也不會很大，從而表明該算法具有較高的精度，可滿足大多數(shù)研究和應用需求。

5 城市熱島效應與分析

在反演LST的基礎上，計算城市熱場變異指數(shù)以定量分析城市熱島效應[27]。某點的LST與研究區(qū)域平均LST的差值同研究區(qū)域平均LST之比，可用于描述該點的熱場變異情況：

HI(T)=(T-TMEAN)/TMEAN

(22)

式中：HI(T)——熱場變異指數(shù)；T——城市某點的遙感反演LST；TMEAN——城市研究區(qū)域的平均LST。

為了更直觀地描述城市熱場變化情況。進一步采用閾值法將熱場變異指數(shù)分為6級[25]。該指標的閾值劃分以及對應的生態(tài)意義見表7。

表7 熱島效應評價指標的閾值劃分

以上文中反演的無錫及其周邊區(qū)域的LST為基礎，計算該地區(qū)的城市熱場變異指數(shù)，并且按照生態(tài)評價指標的閾值進行分級，分別制作出2013年12月10日、2014年3月16日和2013年8月12日，分別代表冬季、春季和夏季的分級圖(圖6)。

圖6 冬季、春季與夏季無錫市熱島效應生態(tài)評價分級

表8為熱島效應等級面積比例與對應的LST平均值。

無錫市地表溫度在不同的季節(jié)變化較大，具體表現(xiàn)為：冬季(12月)地表最高溫度為290.4 K，最低地表溫度為278.9 K，平均溫度為284.7 K；春季(3月)地表最高溫度為304.2 K，最低地表溫度為283.4 K，平均溫度為294.7 K；夏季(8月)地表最高溫度為324.9 K，最低地表溫度為292.5 K，平均溫度為312.3 K，季節(jié)變化明顯。

表8 熱島效應等級面積比例與對應的LST平均值

從熱島效應等級面積比例來看，冬季相對高溫區(qū)所占面積比例較少，強和極強所占面積比例僅為2.52%；春季與夏季相對高溫區(qū)所占面積比例相近，強和極強所占比例分別為18.07%和17.31%。冬季相對低溫區(qū)面積比例最高，無與弱所占比例為66.52%；夏季次之，無與弱所占比例為62.56%；春季相對低溫區(qū)面積比例最高，無與弱所占比例為47.18%。

從熱島效應強度來看，冬季熱島效應等級為無時，LST的平均溫度為282.15 K，熱島效應等級為極強時，LST的平均溫度為291.74 K，差異 9.56 K；春季熱島效應熱島效應等級為無時，LST的平均溫度為289.73 K，熱島效應等級為極強時，LST的平均溫度為302.35 K，差異為12.62 K；冬季熱島效應等級為無時，LST的平均溫度為306.32 K，熱島效應等級為極強時，LST的平均溫度為2 323.52 K，差異為17.2 K。夏季的熱島效應強度最強，春季次之，冬季最弱。

利用之前完成的土地利用分類結(jié)果，結(jié)合地面溫度反演結(jié)果，得到不同地表覆蓋類型的LST平均值，如表9所示。

表9 各地表覆蓋類型對應的LST平均值 K

從表9可以看出，不同地表覆被類型的LST相差較大，城鎮(zhèn)和裸土表面溫度相對較高，植被和水體溫度相對較低。在春夏季，不同地表覆蓋類型之間的地表溫度差異尤其明顯。這是由于地表性質(zhì)的差異對地表能量交換有重要影響。在同樣的太陽輻射條件下，水體因反射率高且比熱較大，升溫慢，植被因蒸騰作用吸收熱量降溫，所以這兩類區(qū)域平均LST較低，對熱島起到緩解作用。而城鎮(zhèn)區(qū)和裸土反射率低且比熱較小，因此LST明顯高于其他地表覆蓋類型?；谏鲜龇治?，可以通過合理規(guī)劃城市建筑密集度以及增加城市綠化帶和水域面積，來改善城市熱島效應。

6 結(jié)語

文中基于Landsat-8的光學和熱紅外數(shù)據(jù)，采用劈窗算法反演了地表溫度LST，并采用同步的MODIS溫度產(chǎn)品和地面實測的16個站點水溫數(shù)據(jù)進行了比較和驗證。結(jié)果表明，基于Landsat-8數(shù)據(jù)的劈窗算法誤差<1 K，比MODIS 溫度產(chǎn)品更為接近地面實測的水溫數(shù)據(jù)。文中進一步對劈窗算法涉及的關鍵輸入?yún)?shù)(大氣水汽含量和地表比輻射率)進行了敏感性分析。結(jié)果表明，在正常的大氣水汽含量和地表比輻射率變化的情況下，其誤差較低(1 K左右)。因此，利用Landsat-8進行劈窗算法反演地表溫度，具備較高的精度，可滿足大多數(shù)研究應用需求。

在LST反演結(jié)果的基礎上，計算了可直觀反映城市熱島效應的HI指數(shù)，發(fā)現(xiàn)無錫市春夏季熱島效應較冬季更為明顯。不同的地表覆蓋類型對熱島效應的貢獻不同，植被和水體因具有明顯的降溫作用，對熱島起到了緩解作用。

筆者在基于劈窗算法和Landsat-8數(shù)據(jù)的LST反演方面取得了初步結(jié)論，但僅使用了一種劈窗算法進行反演計算。為進一步提高地表溫度反演的精度，在今后的工作中，將基于Landsat-8數(shù)據(jù)采用多種反演方法進行比對分析，以更好地輔助城市熱島效應的定量研究。

[1] 葉柯,覃志豪.基于MODIS數(shù)據(jù)的南京市夏季城市熱島分析[J].遙感技術與應用,2006(5):426-431.

[2] 楊士弘.城市生態(tài)環(huán)境學[M].北京:科學出版社，2002.

[3] 劉文淵,謝亞楠,萬智龍,等.不同地表參數(shù)變化的上海市熱島效應時空分析[J].遙感技術與應用,2012(5):797-803.

[4] 劉玉安,唐志勇,程濤,等.基于HJ-1B數(shù)據(jù)的武漢市LST反演及熱環(huán)境分析[J].長江流域資源與環(huán)境,2014(4):526-533.

[5] 吳傳慶,王橋,王文杰,等.利用TM影像監(jiān)測和評價大亞灣溫排水熱污染[J].中國環(huán)境監(jiān)測,2006(6):80-84.

[6] 白楊,孟浩,王敏,等.上海市城市熱島景觀格局演變特征研究[J].環(huán)境科學與技術,2013(3):196-201.

[7] 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.

[8] LI Z L, TANG B H,WU H,et al.Satellite-derived land surface temperature: Current status and perspectives[J]. Remote Sensing of Environment,2013,131:14-37.

[9] 覃志豪,ZHANG M H, ARNON K. 用NOAA-AVHRR熱通道數(shù)據(jù)演算地表溫度的劈窗算法[J].國土資源遙感,2001(2):33-42.

[10] JUAN C J,SOBRINO J A,NINYEROLA M,et al.Revision of the single-channel algorithm for land surface temperature retrieval from Landsat thermal-infrared data[J].IEEE Trans.Geosci.Remote Sens,2009,47(1):339-349.

[11] COLL C,CASSELLES V,SOBRINO A. On the Atmospheric Dependence of the Spit-window Equation for Land Surface Temperature[J].International Journal of Remote Sensing,1994(15):1105-1112.

[12] 柳欽火,徐希孺.遙測地表溫度與比輻射率的迭代反演方法[J].遙感學報,1998,2(1):1-9.

[13] FRANCA G B,GRACKNELL A P.Retrieval of Land and Sea Surface Temperature Using NOAA-11 AVHRR Data in Northeastern Brazil[J].International Journal of Remote Sensing,1994(15):1695-1712.

[14] 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.

[15] 劉玉潔,楊忠東.MODIS遙感信息信息處理原理與方法[M].北京：科學出版社，2001.

[16] 俞宏,石漢青.利用分裂窗算法反演陸地表面溫度的研究進展[J].氣象科學,2002,22(4):494-500.

[17] 覃志豪,李文娟,徐斌,等.陸地衛(wèi)星TM6波段范圍內(nèi)地表比輻射率的估計[J].國土資源遙感,2004(3):28-32.

[18] HUMES K S,KUSTAS W P,MORAN M S.Variability of emissivity and surface temperature over a sparsely vegetated surface[J].Water Resources Research,1994,30:1299-1310.

[19] 毛克彪,覃志豪, 施建成,等.針對MODIS影像的劈窗算法研究[J].武漢大學學報：信息科學版,2005(8):703-707.

[20] 覃志豪,LI W J,ZHANG M H.單窗算法的大氣參數(shù)估計方法[J].國土資源遙感,2004(5):37-43.

[21] 毛克彪.用于MODIS數(shù)據(jù)的地表溫度反演方法研究[D].南京：南京大學,2004.

[22] MAO K,QIN Z,SHI J.A practical split-window alogorithm. for retrirving L and surface temperature from MODIS data[J].Intemational Joumal of Remote Sensing,2005,15(26):3181-3204.

[23] WAN Z M,LI Z L.A Physics-Based Algorithm for Retrieving Land-Surface Emissivity and Temperature From EOS/MODIS Data[J].IEEE Transaction on Geoscience and Remote Sensing,1997,35(2):980-996.

[24] WAN Z M,ZHANG Y L,ZHANG Q C.Quality Assessment and Validation of the MODIS Global Land Surface Temperature[J].International Journal of Remote Sensing,2004,25(1):261-274.

[25] 高懋芳,覃志豪,劉三超.MODIS數(shù)據(jù)反演地表溫度的參數(shù)敏感性分析[J].遙感信息,2005(6):3-6.

[26] QIN Z H,DALL O G.Derivation of spit window algorithm and its sensitivity analysis for retrieving land surface temperature from NOAA-advanced very high resolution radiometer data[J].Geophys Research, 2001,106(D19):22655-22670.

[27] 張勇,余濤,顧行發(fā). CBERS-02 IMRSS熱紅外數(shù)據(jù)地表溫度反演及其在城市熱島效應定量分析中的應用[J].遙感學報,2006,10(5):789-797.