石家莊市主城區(qū)熱島特征及變化分析

李外賓， 湯 軍， 高賢君, 李志遠(yuǎn)

(長江大學(xué)地球科學(xué)學(xué)院，武漢 430100)

中國的城市化進(jìn)程隨著經(jīng)濟(jì)的發(fā)展正在穩(wěn)步推進(jìn)。城市建設(shè)用地?cái)U(kuò)張、人口規(guī)模增長與空間格局的調(diào)整導(dǎo)致大量自然地表被人為改變，城市熱島效應(yīng)成為當(dāng)今全球面臨最嚴(yán)峻的生態(tài)環(huán)境問題之一[1]。熱島效應(yīng)由英國氣候?qū)W家Hawke[2]于19世紀(jì)初在《倫敦的氣候》中首次提到，1985年Manley首次提出了城市熱島這一概念詞[3]。城市熱島是指城市中心溫度明顯高出周邊郊區(qū)溫度的現(xiàn)象[4]，這一現(xiàn)象的出現(xiàn)伴隨著城市病的產(chǎn)生，例如空氣污染、城市高溫等，不僅影響人類的正常工作學(xué)習(xí)，還會危害身體健康。

城市熱島的研究大多是依靠氣象觀測獲取或者遙感技術(shù)的處理，分析相關(guān)因素對熱島影響、熱島空間格局演變和時(shí)序演變規(guī)律。在相關(guān)因素研究上，王寶強(qiáng)等[5]采用相關(guān)分析方法證明轉(zhuǎn)變的產(chǎn)業(yè)結(jié)構(gòu)化、能源消耗和機(jī)動車數(shù)目增多是改城區(qū)增溫重要表現(xiàn)因素；Niu等[6]基于行政邊界方法，揭示了GDP和人口等人為因素是城市熱島增強(qiáng)的重要因素；黃初冬等[7]基于Pearson相關(guān)性分析、灰色綜合關(guān)聯(lián)分析和回歸分析等方法，證明海綿城市建設(shè)可以緩解城市熱島效應(yīng)。在針對特定城市熱島的空間格局演變上，Liu等[8]基于決策樹算法，探討了深圳地區(qū)城市熱島環(huán)境與土地利用格局演變及相互關(guān)系；錢敏蕾等[9]揭示了城市熱島分布與城市擴(kuò)張的規(guī)律；李恒凱等[10]基于混合像元線性波譜分離算法，證明城區(qū)熱島形態(tài)與平均不透水面指數(shù)和平均植被指數(shù)關(guān)系密切。時(shí)序演變規(guī)律的研究上，蘇玥等[11]基于MODIS遙感衛(wèi)星數(shù)據(jù)，發(fā)現(xiàn)內(nèi)蒙古大部分區(qū)域夜間城市熱島效應(yīng)強(qiáng)于白天；潘瑩等[12]采用定性和定量得出重慶市夏季和冬季的熱島效應(yīng)較為顯著，并且在2001—2011年時(shí)間上有愈發(fā)嚴(yán)重趨勢。

石家莊建設(shè)用地?cái)U(kuò)展、人口增加和人為熱源的增加，導(dǎo)致熱島效應(yīng)日益嚴(yán)重，夏季尤為明顯。為了有效減緩熱島效應(yīng)影響，研究石家莊城市熱島空間格局變化、擴(kuò)張情況和下墊面相關(guān)分析。單窗算法和劈窗算法對地表比輻射率誤差不敏感，劈窗算法反演Landsat8的地表溫度又有較大誤差[13-14]，本文選取2004、2010、2016、2020年的Landsat數(shù)據(jù)用輻射傳導(dǎo)方程法對石家莊市的7個(gè)主城區(qū)(橋西區(qū)、新華區(qū)、長安區(qū)、裕華區(qū)、鹿泉區(qū)、欒城區(qū)和藁城區(qū))進(jìn)行地表溫度的反演。根據(jù)反演結(jié)果，計(jì)算地表溫度的變化情況，結(jié)合土地利用和標(biāo)準(zhǔn)差橢圓分析城市熱島的變遷方向和擴(kuò)張趨勢，研究地表溫度與歸一化植被指數(shù)(normalized difference vegetation Index，NDVI)、建筑指數(shù)(normalized difference built-up index，NDBI)和溫度植被干旱指數(shù)(temperature vegetation dryness index，TVDI)之間的相關(guān)性，以此為石家莊城市國土空間優(yōu)化布局提供理論依據(jù)依據(jù)。

1 研究區(qū)概況

石家莊市地處河北省中南部，位于北緯37°27′～38°47′，東經(jīng)113°30′～115°20′，海拔多處于30～100 m，常住人口數(shù)量有千萬多。石家莊地處太行山東麓，東進(jìn)的氣流抵達(dá)太行山后下沉隨后溫度增加，呈“西風(fēng)效應(yīng)”，使得石家莊市溫度高于周邊，年平均氣溫偏高，夏季氣溫經(jīng)?？蛇_(dá)到35 ℃，太行山阻擋，靜風(fēng)天氣常有，污染物和熱量難以散去，加劇熱島效應(yīng)。石家莊市轄8個(gè)區(qū)包括新華區(qū)、橋西區(qū)、長安區(qū)、裕華區(qū)、藁城區(qū)、鹿泉區(qū)、欒城區(qū)、井陘礦區(qū)，選取除井陘礦區(qū)外的7個(gè)區(qū)為研究對象，其包含了絕大部分的城市建設(shè)和人為活動，研究區(qū)概況如圖1所示。

圖1 研究區(qū)概況圖Fig.1 Survey map of the study area

2 數(shù)據(jù)及處理方法

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

以日期相差不大的夏季Landsat影像為主要數(shù)據(jù)源，選用了專題制圖儀(thematic mapper，TM)在2004年8月和2010年8月的成像數(shù)據(jù)，熱紅外傳感器(thermal infrared sensor,TIRS)在2016年8月和2020年5月成像的四期數(shù)據(jù)(表1)，每期均由行列號124/33和124/34拼接而成，拼接前每幅影像均進(jìn)行過輻射定標(biāo)和大氣校正預(yù)處理。數(shù)據(jù)可從地理空間數(shù)據(jù)云(http://www.gscloud.cn/)中獲取研究區(qū)內(nèi)無云的數(shù)據(jù)，同時(shí)，為了獲取大氣水汽含量信息，可以在NASA(https://atmcorr.gsfc.nasa.gov/)網(wǎng)站獲取[3]。

表1 遙感數(shù)據(jù)參數(shù)Table 1 Remote sensing data parameters

2.2 地表溫度的反演

選取石家莊市主城區(qū)2004、2010、2016、2020年的四期影像作為數(shù)據(jù)源，用輻射傳導(dǎo)方程法反演地表溫度。輻射傳導(dǎo)方程法就是先獲取到大氣水汽含量等信息后預(yù)估大氣對地表熱輻射的影響，再把這些大氣影響從衛(wèi)星傳感器觀測到的熱輻射總量中減去,最終得到地表熱輻射強(qiáng)度，最后把這一熱輻射強(qiáng)度轉(zhuǎn)化為相應(yīng)的地表溫度[15]。

衛(wèi)星收到的熱紅外輻射亮度值Lλ由大氣向上輻射亮度、地面的真實(shí)輻射亮度通過大氣層之后到衛(wèi)星傳感器的能量、大氣向下輻射到地面后反射的能量共同組成[16]。表達(dá)式為

Lλ=[εB(Ts)+(1-ε)Ldown]τ+Lup

(1)

式(1)中：ε為地表比輻射率；Ts為地表真實(shí)溫度，K；B(Ts)為黑體熱輻射亮度；τ為大氣在熱紅外波段的透過率；Lup為大氣向上輻射亮度；Ldown為大氣向下輻射亮度[16-17]。B(Ts)表達(dá)式為

B(Ts)=[Lλ-Lup-τ(1-ε)Ldown]/(τε)

(2)

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

Ts=K2/ln[K1/B(Ts)+1]-273.15

(3)

式(3)中：K1、K2為常數(shù)，Landsat5的K1=607.76，K2=1 260.56；Landsat8的K1=774.89，K2=1 321.08。

式(1)中所需的大氣剖面信息可以通過輸入影像獲取的年月日、中心經(jīng)緯度、對應(yīng)數(shù)據(jù)的衛(wèi)星系列和郵箱在NASA網(wǎng)站(https://atmcorr.gsfc.nasa.gov/)獲取。

計(jì)算地表比輻射率還需先將NDVI算出，計(jì)算公式為

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

(4)

式(4)中：NIR為近紅外波段的反射率值；R為紅外波段反射率值。

地表比輻射率用到Sobrino等[18]提出的NDVI閾值法計(jì)算，公式為

ε=0.004Pv+0.986

(5)

式(5)中：Pv是植被覆蓋度[19]，計(jì)算公式為

Pv=(NDVI-NDVIsoil)/(NDVIveg-NDVIsoil)

(6)

式(6)中：NDVIsoil表示完全裸土或無植被覆蓋區(qū)域的NDVI值，通過經(jīng)驗(yàn)值選取NDVI為2%的值，NDVIveg表示完全被植被覆蓋的區(qū)域的NDVI值，選取NDVI為98%的值[3]。理論上低植被覆蓋區(qū)域的NDVI值應(yīng)該趨近于0，全被植被覆蓋的區(qū)域應(yīng)該趨近于1，但是由于土壤類型、地表濕度等因素導(dǎo)致誤差的存在，因此選擇這種經(jīng)驗(yàn)值的辦法盡可能保證數(shù)值趨近于理論值。

2.3 熱島強(qiáng)度的計(jì)算

將地表溫度反演出來以后，為更清晰直觀地看出地溫的區(qū)域分布情況以及時(shí)間跨度上的變化。首先用地表溫度正歸化的方法，使得地溫的值域區(qū)間在0～1，然后用ArcGIS中的自然間隔法、標(biāo)準(zhǔn)差法、等間隔法等，發(fā)現(xiàn)分出來的都會有部分誤差較大。自然間隔方法在低溫區(qū)的劃分上有一定誤差，根據(jù)付瑜等[20]的研究使用0.1劃分后，無論在2004年、2010年、2016年還是2020年劃分的效果都不錯(cuò)，最后選擇自然間隔方法，把低溫區(qū)與次低溫區(qū)以0.1劃分。溫度劃分為五個(gè)級別，分別為高溫區(qū)、次高溫區(qū)、中溫區(qū)、次低溫區(qū)、低溫區(qū)，中溫區(qū)可視為平均溫度，高溫區(qū)、次高溫區(qū)可視為城市的熱島區(qū)域。地表溫度正規(guī)化[21]公式為

Ni=(Ti-Tmin)/(Tmax-Tmin)

(7)

式(7)中：Ni第i個(gè)像元正規(guī)化后得到的溫度；Ti是第i個(gè)像元原始地表溫度；Tmin、Tmax分別是原始地表溫度統(tǒng)計(jì)得到的最小值、最大值。

3 結(jié)果與分析

3.1 地表溫度變化分析

將Landsat數(shù)據(jù)處理后得到地表溫度，計(jì)算出2004年的最低溫是21.6 ℃、最高溫50.44 ℃和平均溫度31.1 ℃；2010年的最低溫與最高溫溫差達(dá)到28.9 ℃、平均溫度為33.25 ℃；2016年的最低溫25.73 ℃、最高溫51.5 ℃、平均溫度33.94 ℃；2020年的最低溫23.64 ℃、最高溫57.06 ℃、平均溫度35.11 ℃。2004—2020年平均溫度從31.1 ℃升高到了35.11 ℃，上升較為明顯；最低溫基本是水體的溫度，雖然水體一直處于低溫區(qū)，但是2004—2020年低溫也在波動性升高。

除了溫度明顯升高，熱島強(qiáng)度占比也發(fā)生了較大變化，2004—2020年低溫區(qū)基本持平，變化微?。淮蔚蜏貐^(qū)總體占比縮小最多，2004—2010年從66.74%上升到67.53%，2010—2016年下降到52.91%，面積縮減了303.02 km2，到2020年繼續(xù)下降到41.82%，面積縮減到916.43 km2；2004—2020年中溫區(qū)上升，先從17.61%下降到15.59%再上升到22.99%，2020年占比達(dá)到24.30%，總體呈上升趨勢；次高溫區(qū)逐年上升，2004年次高溫區(qū)占10.87%，2010年占12.99%，2016年上升到17.55%，2020年達(dá)到23.03%；高溫區(qū)增長較多，2004年高溫區(qū)占比為4.26%，2010年下降成3.32%，下降原因是因?yàn)?004年反演的地溫中，一部分高溫區(qū)來自滹沱河，滹沱河沉積物的泥沙含量較高，在嚴(yán)重缺水下地表溫度較高，2016年又上升到5.96%，面積從2004年的93.39 km2增長到2016年的130.65 km2，到2020年高達(dá)227.82 km2，總體呈現(xiàn)上升。具體面積的詳細(xì)信息變化如表2所示。

表2 2004—2020年熱島強(qiáng)度面積統(tǒng)計(jì)Table 2 Statistics of heat island intensity area from 2004 to 2020

從圖2可以直接看出石家莊市城區(qū)基本都處于熱島狀態(tài)，溫差界限在逐漸外擴(kuò)，碎斑化的區(qū)域也逐漸合并，熱島效應(yīng)存在且加強(qiáng)。低溫區(qū)主要在西北部的黃壁莊水庫，2004—2020年一直呈現(xiàn)著最低溫；次低溫區(qū)覆蓋面積最廣，因?yàn)檫x擇的影像在五月底和八月中下旬，正處于農(nóng)作物生長茂盛的時(shí)期，很明顯可以看出有農(nóng)作物的地方基本都處在次低溫區(qū)溫度區(qū)間；2004—2016年新華區(qū)、橋西區(qū)、長安區(qū)、裕華區(qū)始終是熱島主要分布位置并保持不斷向外擴(kuò)散；另外還存在四個(gè)小面積熱島區(qū)域，分別是鹿泉區(qū)、藁城區(qū)和欒城區(qū)的城鎮(zhèn)點(diǎn)，以及銅冶鎮(zhèn)，到2020年熱島區(qū)域的持續(xù)外擴(kuò)，離散的鹿泉區(qū)、銅冶鎮(zhèn)和欒城區(qū)與城區(qū)熱島連接呈片，出現(xiàn)大面積熱島區(qū)。從石家莊市主城區(qū)地表溫度分布圖上可以清晰看出分布情況(圖2)。

圖2 2004—2020年石家莊市主城區(qū)地表溫度分布情況Fig.2 Distribution of surface temperature in main urban area of Shijiazhuang City from 2004 to 2020

3.2 土地利用格局演變

基于支持向量機(jī)對研究區(qū)進(jìn)行分類，土地利用變化明顯，呈現(xiàn)主城區(qū)向外蔓延擴(kuò)張[22]，結(jié)果如圖3所示。結(jié)合轉(zhuǎn)移矩陣(表3～表5)，可以看出2004—2020年石家莊市主城區(qū)的土地利用的變化。建筑用地面積增加了552.77 km2，其中大部分由耕地和裸地及其他轉(zhuǎn)變而來，2004—2010年，55.61 km2的耕地和50.16 km2的裸地及其他變更成建筑用地；2010—2016年，97.40 km2的耕地和116.96 km2的裸地及其他變更成建筑用地；2016—2020年，207.74 km2的耕地和108.04 km2的裸地及其他變更成建筑用地。其次變化比較明顯的是耕地和裸地及其他，16年的時(shí)間，耕地減少了549.75 km2，裸地及其他減少了134.67 km2；林地在2004—2010年多被砍伐或用來耕作等造成面積減少，到2016年增加到108.40 km2，2020年達(dá)到了244.70 km2，這與2015年啟動的太行山生態(tài)綠化工程和2018年的重點(diǎn)造林綠化工程有很大關(guān)系；另外石家莊市是一個(gè)缺水型城市，因此水體面積變化微小，并且總面積小。石家莊市城市化進(jìn)程，大量耕地、裸地及其他等被人為建筑取代，建筑用地大量增添，城市熱島嚴(yán)重加劇。

表3 2004—2010年土地利用轉(zhuǎn)移矩陣Table 3 Land use transfer matrix from 2004 to 2010

表4 2010—2016年土地利用轉(zhuǎn)移矩陣Table 4 Land use transfer matrix from 2010 to 2016

表5 2016—2020年土地利用轉(zhuǎn)移矩陣Table 5 Land use transfer matrix from 2016 to 2020

圖3 2004—2020年石家莊市主城區(qū)土地利用Fig.3 Land use in Shijiazhuang City from 2004 to 2020

3.3 熱島擴(kuò)張的方向特征

為了進(jìn)一步可以看出熱島擴(kuò)張的方向性，采用標(biāo)準(zhǔn)差橢圓的方法進(jìn)行熱島的空間分析(圖4、圖5)。結(jié)果顯示，熱島的整體向著西北和東南方向擴(kuò)張，建筑用地整體同樣向著西北和東南方向擴(kuò)張，結(jié)合表6、表7可以看到方向角度近似，中心點(diǎn)基本不變，長軸趨于增長趨勢，這也是說明熱島與建筑用地具有良好的一致性，建筑用地在空間擴(kuò)張的情況下會帶動熱島的擴(kuò)張[23-24]。

表6 熱島變化標(biāo)準(zhǔn)差橢圓信息表Table 6 Ellipse information table of heat island variation standard deviation

表7 建筑用地變化標(biāo)準(zhǔn)差橢圓信息表Table 7 Ellipse information table of standard deviation of building land change

圖4 熱島變化標(biāo)準(zhǔn)差橢圓圖Fig.4 Standard deviation ellipse of heat island variation

圖5 建筑用地變化標(biāo)準(zhǔn)差橢圓圖Fig.5 Standard deviation ellipse of heat island variation

向西北，主城區(qū)與鹿泉區(qū)之間交通路網(wǎng)發(fā)達(dá)，與其他主城區(qū)通行便捷。鹿泉在持續(xù)發(fā)展的進(jìn)程中，建筑設(shè)施也逐漸完善發(fā)展，然而，鹿泉區(qū)位于石家莊的“上風(fēng)口”，也是重要的生態(tài)區(qū)，主城區(qū)向西持續(xù)發(fā)展的空間會在未來受限。

向東南，隨著高新區(qū)的開發(fā)，以及高新區(qū)與主城區(qū)之間緊密的城區(qū)發(fā)展格局，將主城區(qū)向此擴(kuò)張，并將太行大街以東、南二環(huán)東延線以南的區(qū)域作為“處女地”不斷開發(fā)，東南方向?qū)⒊蔀槲磥碇饕^續(xù)延伸發(fā)展的方向。

3.4 城市熱島效應(yīng)與下墊面關(guān)系

3.4.1 溫度與NDVI關(guān)系

由圖6可以看出地表溫度與NDVI相關(guān)性較高并呈負(fù)相關(guān)，NDVI越高地表溫度越低，NDVI越低地表溫度越高，這與崔林林等[25]研究結(jié)果統(tǒng)一。樣本點(diǎn)在NDVI越高，與趨勢線的聚集程度越高，較高的NDVI代表該像元大概率為植被，植被有明顯調(diào)控作用[26]。

圖6 地表溫度與NDVI的相關(guān)性Fig.6 Correlation between surface temperature and NDVI

3.4.2 溫度與NDBI關(guān)系

NDBI建筑信息提取公式為

NDBI=(SWIR-NIR)/(SWIR+NIR)

(8)

式(8)中：SWIR為短波紅外的反射率值；NIR為近紅外的反射率值。

由圖7可以看出地表溫度與NDBI高度相關(guān)并呈正相關(guān)，NDBI越高地表溫度越高，NDBI越低地表溫度越低[27]。樣本點(diǎn)在NDBI大于0時(shí)，與趨勢線的離散程度較大，這與城市土地利用決策規(guī)劃和管理有很大關(guān)系[28]。

圖7 地表溫度與NDBI的相關(guān)性Fig.7 Correlation between surface temperature and NDBI

3.4.3 溫度與TVDI關(guān)系

TVDI是土壤水分反演的一種方法[29]。TVDI的值域?yàn)閇0，1]，值越大土壤越干旱缺水，值越小，土壤水分含量越高，公式為

TVDI=(Ts-Tmin)/(Tmax+Tmin)

(9)

式(9)中：Ts是任意影像像元地表溫度值；Tmin、Tmax分別是統(tǒng)計(jì)得到NDVI的最低和最高地表溫度值。

在簡化特征空間中，對干邊(Tmax)、濕邊(Tmin)進(jìn)行線性擬合，公式為

Tmin=a1+b1NDVI

(10)

Tmin=a2+b2NDVI

(11)

式中：a1、b1、a2、b2為方程系數(shù)。

由圖8可以看出地表溫度與TVDI正相關(guān)。整體的樣本點(diǎn)都與趨勢線較高的趨近，這說明地表溫度與地表的濕潤程度密切關(guān)聯(lián)，可以用來重點(diǎn)考慮成調(diào)控城市熱島問題方法。

圖8 地表溫度與TVDI的相關(guān)性Fig.8 Correlation between surface temperature and TVDI

4 結(jié)論

通過輻射傳導(dǎo)方程法計(jì)算出石家莊市主城區(qū)的地表溫度后，用了自然分割法結(jié)合0.1值劃分低溫區(qū)的方法分析熱島的空間格局、強(qiáng)度變化，以及土地利用格局變化對熱島的影響，再研究地表溫度與NDVI、NDBI和TVDI的相關(guān)性，得出結(jié)論并給出建議。

(1)2004—2020年反演得到的地表溫度的平均溫度從31.1 ℃升高到了35.11 ℃，最低溫度從21.6 ℃升高到了23.64 ℃，最高溫度均達(dá)到50 ℃以上，增溫趨勢明顯。

(2)2004—2020年熱島(次高溫區(qū)和高溫區(qū))面積增加，新華區(qū)、長安區(qū)、橋西區(qū)和裕華區(qū)是熱島區(qū)域聚集程度最高的地方，另外存在鹿泉區(qū)、藁城區(qū)和欒城區(qū)和銅冶鎮(zhèn)四個(gè)孤立熱島，并且不斷向外擴(kuò)張，擴(kuò)張情況與建筑用地趨于一致，整體上向著西北和東南方向擴(kuò)張。

(3)建筑用地是城市地表升溫的重要因素，植被覆蓋度和水體能夠?qū)Φ乇頊囟日{(diào)控起到有效作用?？梢酝ㄟ^改變城市格局及建筑材料來降低城市熱島效應(yīng)，比如使用太陽光反射率好或能降溫節(jié)能材料等[25]，做到“適地適樹”的原則，選擇適合植被種植，增加城市綠地面積，城市定時(shí)灑水、建造公園濕地等方法解決石家莊的熱島效應(yīng)問題[30]。