閆李月， 左小清*， 王占峰， 李洪忠

(1.昆明理工大學(xué) 國土資源工程學(xué)院， 云南 昆明 650093；2.中國科學(xué)院 深圳先進(jìn)技術(shù)研究院， 廣東 深圳 518052)

城市熱島(Urban Heat Ialand，UHI)是指由城市化引起的城市地表及大氣溫度高于周邊非城市環(huán)境的一種現(xiàn)象[1]。隨著全球城市化進(jìn)程的不斷加快，近幾十年城市熱島現(xiàn)象不斷加劇，引起了社會的廣泛關(guān)注[2-3]。當(dāng)前城市熱島研究通常采用氣象觀測、遙感圖像溫度反演等方法，對熱島的時空分布及演變規(guī)律進(jìn)行研究。其中，氣象觀測的方法是基于有限的地面氣象站數(shù)據(jù)，覆蓋范圍有限且分布不均勻，很難全面地反映城市熱島的空間分布情況。而遙感圖像具有覆蓋范圍廣、時效性強(qiáng)、數(shù)據(jù)綜合性和可比性強(qiáng)等優(yōu)勢，使其在城市熱島的研究中得到了越來越廣泛的應(yīng)用[4]。

為了開展城市熱島效應(yīng)演變研究，需要將不同時相遙感圖像的熱紅外波段加以比較[3]。然而，不同時相遙感圖像在獲取時，大氣條件、傳感器輻射定標(biāo)、太陽輻射強(qiáng)度等的不同，導(dǎo)致圖像中同一目標(biāo)地物具有較大的輻射差異[5]。即便在數(shù)據(jù)預(yù)處理過程中已經(jīng)普遍采用FLAASH模型或ATCOR模型等方法對獲取的遙感圖像進(jìn)行輻射校正[6]，時相差異卻仍然對地表溫度反演結(jié)果產(chǎn)生著不同程度的影響。如同一地點(diǎn)地表覆蓋不變，反演得到的溫度卻有較大差異，即地表溫度的“偽變化”[7]，從而對熱島現(xiàn)象演變分析產(chǎn)生較大干擾。因此，要想進(jìn)行基于遙感的城市熱島的動態(tài)變化分析，就必須考慮對多時相遙感地表溫度圖像的正規(guī)化處理。將不同時相的地表溫度圖像進(jìn)行正規(guī)化處理后，解決客觀條件對地表溫度反演的干擾問題，從而進(jìn)一步研究時間序列的城市熱島演變。

目前，城市熱島演變研究中不同時相地表溫度圖像正規(guī)化處理方法效果較好，被廣泛應(yīng)用的主要有兩種：典型歸一化方法[3,8]和偽不變特征(Pseudo Invariant Feature，PIF)法[9]。前者由徐涵秋于2003年首次應(yīng)用于多時相地表溫度圖像的歸一化處理，后被學(xué)者廣泛應(yīng)用于城市熱島演變的相關(guān)研究。后者由Schottetal于1988年首次提出，目的是使多時相圖像具有相同的輻射尺度(或輻射分辨率)。多用于對遙感影像的相對輻射校正，目前有少數(shù)學(xué)者將其應(yīng)用于多時相地表溫度處理[10]，并取得了良好效果。

深圳市近20年來經(jīng)濟(jì)、文化等飛速發(fā)展，使城市建成區(qū)擴(kuò)展的速度迅速，規(guī)模較大，同時也使得其城市熱島現(xiàn)象明顯加劇。因此，本文以深圳市陸域范圍為例，選用Landsat系列數(shù)據(jù)，在不同年份、不同季節(jié)的時間序列上城市熱島效應(yīng)演變研究中，討論上述兩種方法對地表溫度圖像正規(guī)化處理效果的差異。

1 研究區(qū)與數(shù)據(jù)

深圳市位于珠江三角洲東岸，陸域范圍22°26′59″—22°51′49″N，113°45′44″—114°37′21″，東臨大亞灣與惠州相鄰，西至珠江口，南望香港特別行政區(qū)，北接莞、惠兩地，整體呈狹長型。需要說明的是，本文所提到的深圳市均表示深圳市陸域范圍，不包括海域。深圳市作為快速城市化的典型代表，30余年間從最初的邊陲小鎮(zhèn)發(fā)展為下轄6個行政區(qū)和4個新區(qū)內(nèi)的60個街道辦事處的國際化都市。城市空間結(jié)構(gòu)以包括福田、羅湖、南山3個行政區(qū)在內(nèi)的中心城區(qū)為核心，形成“南北貫通，西聯(lián)東拓”的區(qū)域空間協(xié)調(diào)發(fā)展框架。

深圳市處于121/43、122/44共2景Landsat圖像上，121/43包含深圳東半部分區(qū)域，122/44包含西半部分區(qū)域。本文選用了Landsat 5 TM、Landsat 8 OLI/TIRS數(shù)據(jù)，在2000、2005、2010、2015年4個時間段各獲取的2景遙感圖像，云量均為0，共獲取8景，數(shù)據(jù)信息如表1所示。由于中國南部沿海常年多云多雨，深圳為典型區(qū)域，Landsat數(shù)據(jù)缺失、質(zhì)量差的現(xiàn)象尤其明顯，同一地區(qū)2年內(nèi)可能只有一景圖像可用。因此，本文選取的遙感圖像日期并非全為當(dāng)?shù)赝患竟?jié)?；诖?，就更需要對地表溫度圖像用正規(guī)化處理方法進(jìn)行分析探討，解決客觀條件對地表溫度反演的干擾問題。

表1 Landsat數(shù)據(jù)信息

另外，本文數(shù)據(jù)處理過程中還選用了深圳市2000、2005、2010、2015年4期的土地利用分類產(chǎn)品，由中國科學(xué)院深圳先進(jìn)技術(shù)研究院空間信息研究中心提供。該數(shù)據(jù)產(chǎn)品的生產(chǎn)基于Landsat系列遙感衛(wèi)星，其中1980年基于80 m分辨率Landsat MSS數(shù)據(jù)，1990、2000、2010年基于30 m分辨率Landsat TM數(shù)據(jù)，2015年基于30 m分辨率Landsat 8 OLI數(shù)據(jù)。數(shù)據(jù)來源于美國地質(zhì)勘探局(United States Geological Survey，USGS)網(wǎng)站(https://landsatlook.usgs.gov)。數(shù)據(jù)產(chǎn)品的生產(chǎn)過程：(1)基于ENVI對所有圖像進(jìn)行輻射定標(biāo)，消除傳感器誤差；(2)基于FLAASH模塊進(jìn)行大氣校正；(3)基于易康軟件建立分類規(guī)則集，進(jìn)行監(jiān)督分類；(4)以人機(jī)交互的方式對分類結(jié)果進(jìn)行修正。該數(shù)據(jù)產(chǎn)品的分類系統(tǒng)包括建設(shè)用地、草地、耕地、林地、濕地、未利用地和近海海域，分類精度達(dá)到95%以上。

2 研究方法

2.1 地表溫度反演

本文選用覃志豪的單窗算法[10]，基于ENVI 5.1軟件平臺，進(jìn)行Landsat 5、Landsat 8圖像地表溫度反演，得到地表溫度圖像。普朗克函數(shù)用于定量描述黑體輻射的能量與波長(頻率)、溫度的客觀規(guī)律，已知物體的能量譜和輻射強(qiáng)度，即可推算出亮度溫度。

根據(jù)普朗克函數(shù)將黑體輻射亮度轉(zhuǎn)化為亮度溫度[10]，公式為

(1)

式中Tb是亮度溫度(K)；Lb為同溫度下的黑體輻射亮度，由對熱紅外波段進(jìn)行輻射定標(biāo)后的絕對輻射亮度L和地表比輻射率ε計算得到；k1、k2為熱紅外波段的參數(shù)，Landsat 8第10波段的參數(shù)為k1=774.89、k2=1 321.08，Landsat 5第6波段的參數(shù)為k1=607.76、k2=1 260.56。

得到亮度溫度Tb后，用覃志豪的單窗算法[10]計算地表溫度(Land Surface Temperature，LST)，公式為

(2)

式中LST為地表溫度(℃)；a、b為常數(shù)，地表溫度在0～70 ℃內(nèi)，a=-67.355 351，b=0.458 606；TC為大氣平均作用溫度，估算方程在中緯度夏季TC=16.011 0+0.926 21×T0，在熱帶大氣下TC=17.976 9+0.917 15×T0，T0為氣象局發(fā)布的地表1.5～2.0 m左右的地面溫度，即氣溫(K)；C、D為中間變量，C=ετ，D=(1-τ)×[1+(1-ε)τ]，ε為比輻射率，τ為大氣透射率。

將本研究選取的8景Landsat圖像分別通過上述方法，計算得到8幅地表溫度圖像。然后將表1中同一年份的2景地表溫度圖像進(jìn)行地表溫度的融合拼接、裁剪[11]等一系列處理，最后得到深圳市陸域范圍2001年3月1日、2005年11月23日、2009年11月2日、2015年10月18日的地表溫度圖像。下文提到的2000、2005、2010、2015年分別指上述日期的LST圖像。

此外，本文在反演地表溫度的過程中，同時獲取了深圳市歸一化植被指數(shù)(Normalized Difference Vegetation Index，NDVI)圖像。計算方法為近紅外波段的反射值與紅光波段的反射值之差比兩者之和，公式為

NDVI=(NIR-R)/(NIR+R)，

(3)

式中NIR為近紅外波段，R為紅光波段。Landsat 5紅光波段(0.63～0.69 μm)為3，近紅外波段(0.76～0.90 μm)為4；Landsat 8紅光波段(0.630～0.680 μm)為4，近紅外波段(0.845～0.885 μm)為5。

2.2 地表溫度圖像正規(guī)化方法

2.2.1 典型歸一化方法

歸一化是一種無量綱處理手段，在城市熱島的研究中，使用此方法主要針對研究熱島現(xiàn)象相對強(qiáng)弱的空間分布。當(dāng)側(cè)重地表溫度相對高低的空間分布特征時，時間不同改變的只是地表溫度值的大小，并不改變地表溫度的空間分布[3]。因此，徐涵秋等采用典型歸一化方法對反演得到的LST圖像進(jìn)行處理，令不同時相的LST圖像分布范圍統(tǒng)一到0～1之間。計算公式為

(4)

式中Ni為第i個像元溫度歸一化后的值，TSi為第i個像元地表溫度值，TSmin為LST圖像像元灰度值的最小值，TSmax為LST圖像像元灰度值的最大值。將4幅LST圖像分別經(jīng)過公式(4)處理后，像元灰度值范圍均被歸一化至0～1之間，此時的圖像即為典型歸一化方法正規(guī)化處理后的圖像。

2.2.2 偽不變特征法

偽不變特征法(Pseudo Invariant Feature，PIF)在遙感圖像處理中也被稱為相對輻射校正(或相對輻射歸一化)方法，是基于同一區(qū)域，多時相圖像相同波段的像元灰度值存在線性關(guān)系的假設(shè)提出來的。首先要確定參考圖像和實(shí)驗(yàn)圖像，然后選取偽不變特征點(diǎn)(或稱為輻射地面控制點(diǎn))，確定偽不變特征點(diǎn)后，提取兩期圖像偽不變特征點(diǎn)的像元灰度值，建立參考圖像與實(shí)驗(yàn)圖像之間的回歸模型[8]。將多時相圖像的像元灰度值逐漸歸一化到參考圖像，使多時相圖像具有相同的輻射尺度(或輻射分辨率)。該方法的具體步驟如下。

2.2.2.1 選取參考圖像

參考圖像應(yīng)該目視質(zhì)量較好，變動范圍較大，以便于信息量盡可能多地保留下來[14]。本文共有2000、2005、2010、2005年共4個時間的地表溫度圖像，前3幅由Landsat 5數(shù)據(jù)反演得到，后1幅由Landsat 8數(shù)據(jù)反演得到。依據(jù)文章2.1小節(jié)，4幅LST圖像的具體日期分別為2001年3月1日(春)、2005年11月23日(秋)、2009年11月2日(秋)和2015年10月18日(夏)，深圳市一年四季植物長勢良好，春秋季與夏季圖像的目視質(zhì)量不相上下。

圖1為上述4幅LST圖像的像元灰度值曲線圖。從圖中可以看出，2015年的像元灰度值變動范圍最大，2010和2000年與之相近，參考圖像可從這3幅圖中選取；另外，同一季節(jié)下圖像的地表溫度“偽變化”會相對較小。因此，本文選取2010年的圖像作為參考圖像，其余3幅圖為實(shí)驗(yàn)圖像。

圖1 LST像元灰度值曲線圖

2.2.2.2 選取偽不變特征點(diǎn)

偽不變特征點(diǎn)是指不同時相的遙感圖像中，地表不隨時間發(fā)生變化或變化很小可以忽略不計的地物，深水體(如水庫)、常年裸巖裸土、大屋頂?shù)鹊湫偷匚颷7]。對于樣本點(diǎn)的選取，樣本容量應(yīng)當(dāng)適中，樣本容量過大或過小都會直接影響回歸模型的結(jié)果。通常偽不變特征點(diǎn)的選取主要為人工目視的方法，得到的點(diǎn)數(shù)有限，且由于個體差異導(dǎo)致特征點(diǎn)的質(zhì)量水平參差不齊。本文選取偽不變特征點(diǎn)的過程主要分為兩個步驟：第1步，基于深圳市4期的土地利用分類產(chǎn)品，分別提取2000、2005、2015年3幅實(shí)驗(yàn)圖像與2010年的參考圖像之間土地利用未發(fā)生變化的區(qū)域，在土地利用不變的區(qū)域上，再分別提取3000個隨機(jī)點(diǎn)，記為S1、S2、S3；第2步，依次提取S1在2000年和2010年遙感圖像的NDVI值S1a和S1b、S2在2005年和2010年遙感圖像的NDVI值S2a和S2b、S3在2015年和2010年遙感圖像的NDVI值S3a和S3b。選取偽不變特征點(diǎn)Pi的規(guī)則如下：

Pi∶|Sia-Sib|<0.06，

(5)

式中i=1、2、3，滿足上述條件的隨機(jī)點(diǎn)，則為選取的偽不變特征點(diǎn)。這樣得到的偽不變特征點(diǎn)點(diǎn)數(shù)適中，分布均勻，且具有客觀、定量化的選取標(biāo)準(zhǔn)。

2.2.2.3 建立回歸模型

回歸模型是以參考圖像偽不變特征點(diǎn)的像元灰度值為縱坐標(biāo)y，以實(shí)驗(yàn)圖像偽不變特征點(diǎn)的像元灰度值為橫坐標(biāo)x進(jìn)行點(diǎn)繪，然后對這些點(diǎn)做回歸分析，求解模型參數(shù)。首先從已有的LST圖像中提取偽不變特征點(diǎn)的LST值：依次提取P1在2000年和2010年的LST值T1a和T1b，P2在2005年和2010年的LST值T2a和T2b，P3在2015年和2010年的LST值T3a和T3b。則回歸模型表示為

Tib=αi+βi×Tia，

(6)

式中i=1、2、3，參數(shù)α和β采用最小二乘法求解。計算得到3組模型參數(shù)，分別代入2000、2005、2015年的LST實(shí)驗(yàn)圖像進(jìn)行計算，得到的3幅圖像即為偽不變特征法正規(guī)化處理后的圖像，2010年LST圖像為參考圖像，不需要再做處理。

以2010年LST圖像為參考圖像，利用上述方法提取出的3組偽不變特征點(diǎn)P1、P2、P3，建立3組偽不變特征回歸模型，如表2所示。

表2 偽不變特征法回歸模型

表2中3組回歸模型均顯示出了較強(qiáng)的相關(guān)性，確定系數(shù)0.5

3 城市熱島演變分析

3.1 熱島強(qiáng)度等級劃分

城市熱島的研究中，有關(guān)地表熱場分布、熱島結(jié)構(gòu)、熱島演變等的研究，都涉及城市熱島的界定問題，即地表溫度的等級劃分[14]。常用的灰度圖像等級劃分方法有等間距法、分位數(shù)劃分法、自然間斷點(diǎn)分級法、幾何間隔法以及均值-標(biāo)準(zhǔn)差法等。本文綜合各分級方法對2010年LST圖像典型歸一化方法正規(guī)化和偽不變特征法正規(guī)化的分級效果，采用目視判別，發(fā)現(xiàn)均值-標(biāo)準(zhǔn)差法劃分等級，可視化處理的效果最佳。因此本文對正規(guī)化后的LST圖像分別以μ(均值)-σ(標(biāo)準(zhǔn)差)設(shè)置分割點(diǎn)，將熱島強(qiáng)度等級劃分為6個級別：低溫區(qū)、次低溫區(qū)、中溫區(qū)、次高溫區(qū)、高溫區(qū)和特高溫區(qū)。等級劃分標(biāo)準(zhǔn)如表3所示，表中Ts為需要劃分等級的圖像。

表3 均值-標(biāo)準(zhǔn)差法劃分熱島強(qiáng)度等級

將分級結(jié)果按照時間順序排列，如圖2所示。對圖2進(jìn)行目視分析，觀察兩種正規(guī)化方法處理過后的熱島強(qiáng)度等級劃分結(jié)果，可以看出：

圖2 兩種正規(guī)化處理方法下的熱島強(qiáng)度等級劃分結(jié)果

(1)2000—2015年間，兩組結(jié)果均體現(xiàn)了深圳市熱島現(xiàn)象隨時間逐漸增強(qiáng)的現(xiàn)象，即“特高溫區(qū)”逐年增多、“低溫區(qū)”逐年減少。且兩者各個熱島級別的空間分布移動特征相似：高溫區(qū)、特高溫區(qū)均表現(xiàn)出了隨時間逐漸往中西部、西北和西南部轉(zhuǎn)移的特征，其中西部與西北部熱島聚集現(xiàn)象最為嚴(yán)重。

(2)兩組結(jié)果在不同年份的熱島強(qiáng)度增強(qiáng)程度存在差異。經(jīng)過典型歸一化方法處理的熱島強(qiáng)度等級劃分結(jié)果中，2000、2005、2010年的熱島轉(zhuǎn)移趨勢明顯，但增強(qiáng)程度相對不明顯；2010—2015年的熱島增強(qiáng)程度相對較明顯。經(jīng)過偽不變特征法處理的熱島強(qiáng)度等級劃分結(jié)果中，2000、2005、2010、2015年各階段熱島轉(zhuǎn)移趨勢和增強(qiáng)程度都較為明顯。以2010年為分界，典型歸一化法2000、2005年的熱島強(qiáng)度比偽不變特征法的強(qiáng)；前者2015年的高溫區(qū)與特高溫區(qū)分布與后者相差不大，低溫區(qū)分布明顯大于后者。

兩組結(jié)果在不同年份的熱島強(qiáng)度增強(qiáng)程度存在差異，表明典型歸一化方法和偽不變特征法對多時相遙感地表溫度圖像的處理效果不同，也就是說，對消除地表溫度的“偽變化”的作用有差異。在本研究中，兩種方法孰優(yōu)孰劣，還需進(jìn)一步對比分析。

3.2 城市熱島演變與土地利用變化關(guān)系

相關(guān)研究表明，城市環(huán)境效應(yīng)產(chǎn)生的內(nèi)在驅(qū)動力是土地利用與土地覆被的變化，城市熱環(huán)境與土地利用類型密切相關(guān)。因此，本文列出4個時間段深圳市不同土地利用類型的分布情況，如圖3所示，以進(jìn)一步對兩種方法的處理效果進(jìn)行分析評價。

圖3 深圳市2000—2015年土地利用類型分布

分別統(tǒng)計兩種正規(guī)化方法處理過后的LST圖像上各熱島強(qiáng)度等級的面積及所占比例，統(tǒng)計結(jié)果如表4、表5所示。統(tǒng)計深圳市2000—2015年不同土地利用類型的面積，如表6所示。

表4 深圳市2000—2015年熱島強(qiáng)度等級面積——典型歸一化方法

表5 深圳市2000—2015年熱島強(qiáng)度等級面積——偽不變特征法

表6 深圳市2000—2015年不同土地利用類型面積 km2

依據(jù)圖3、表6土地利用類型分布的變化，可見深圳市人工表面與林地占地面積最多。2000、2005、2010年間人工表面占地面積大量增加，2010—2015年增加不明顯；2000—2015年林地、耕地和濕地面積不斷減少，草地占地面積少且緩慢增加。結(jié)合圖2、圖3進(jìn)行目視對比，整體來看，熱島強(qiáng)度級別為高溫區(qū)及以上的區(qū)域變化與人工表面的分布變化接近，低溫區(qū)的變化與林地的分布變化接近。已有研究[15]表明，人工表面的占比對熱島現(xiàn)象的強(qiáng)弱有很大影響，而林地對城市熱島效應(yīng)具有減緩作用。因此，本文依據(jù)表4—表6，以人工表面與林地為主，將高溫區(qū)與特高溫區(qū)作為熱島區(qū)域、次低溫區(qū)與低溫區(qū)作為熱島緩解區(qū)域，繪制不同區(qū)域面積變化與土地利用變化對比折線圖(圖4)。

(a) 熱島區(qū)域 (b) 熱島緩解區(qū)域

對比圖4(a)熱島區(qū)域和人工表面的面積變化趨勢，隨著人工表面的顯著增加。經(jīng)過偽不變特征法處理的熱島區(qū)域面積變化趨勢與人工表面面積變化趨勢相同；而經(jīng)過典型歸一化方法處理的熱島區(qū)域面積卻并沒有明顯變化，2015年熱島區(qū)域面積甚至少于2010年，這是不合理的。

對比圖4(b)熱島緩解區(qū)域和人工表面的面積變化趨勢，林地面積在2000—2015年間不斷減少。經(jīng)過典型歸一化方法處理的熱島緩解區(qū)域面積在2010年之前的變化趨勢與林地變化相近，但相反地，在2015年熱島緩解區(qū)域面積出現(xiàn)回升，整體表現(xiàn)為2000年與2015年幾乎沒有變化。經(jīng)過偽不變特征法處理的熱島緩解區(qū)域面積在2000—2015年間表現(xiàn)為先減少、再增加、又減少，但整體是減少的，在2000年與2015年兩個時間節(jié)點(diǎn)與林地面積減少量相對應(yīng)；對于2005年熱島緩解區(qū)域面積的大幅減少，結(jié)合表5，發(fā)現(xiàn)該現(xiàn)象正好對應(yīng)了2005年中溫區(qū)面積的大幅增加。

從以上兩方面來看，在本研究中，經(jīng)過偽不變特征法處理過后的效果更好。

4 結(jié)語

實(shí)際上，兩種正規(guī)化方法在城市熱島研究中都已有大量的論證和應(yīng)用，各自都具有其運(yùn)用價值。通過本研究，說明了典型歸一化處理方法對年際變化的城市熱島效應(yīng)演變研究適用性有限，特別是針對深圳市這樣發(fā)展迅速、城市環(huán)境年際變化劇烈的地區(qū)，采用偽不變特征法處理后所得的結(jié)果相對較好。宏觀上，從目視分析的結(jié)果來看，經(jīng)過典型歸一化方法處理所得的熱島等級劃分結(jié)果在熱島分布、熱島的空間轉(zhuǎn)移方面，與經(jīng)過偽不變特征法處理所得的結(jié)果差別不大。但微觀上，在城市熱島演變與土地利用變化關(guān)系的量化分析中，兩組結(jié)果差異明顯，偽不變特征法處理得到的結(jié)果更合理。

本研究也存在不足之處：(1)較典型歸一化方法而言，偽不變特征法處理過程相對繁瑣。若是針對日期相近、或一定時間范圍內(nèi)研究區(qū)土地利用類型變化不大的地表溫度圖像，用典型的歸一化方法處理結(jié)果同樣準(zhǔn)確，且更方便。(2)偽不變特征法的精度較難以控制，特別是偽不變特征點(diǎn)的選取。目前常采用人工目視選點(diǎn)的方法，人工判斷失誤不可控且不適合較大區(qū)域。本文雖利用已有的土地利用分類產(chǎn)品設(shè)立了客觀、定量的選取方法，但回歸模型精度有所下降，方法有待改進(jìn)。