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

李外賓， 湯 軍， 高賢君, 李志遠

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

中國的城市化進程隨著經濟的發展正在穩步推進。城市建設用地擴張、人口規模增長與空間格局的調整導致大量自然地表被人為改變，城市熱島效應成為當今全球面臨最嚴峻的生態環境問題之一[1]。熱島效應由英國氣候學家Hawke[2]于19世紀初在《倫敦的氣候》中首次提到，1985年Manley首次提出了城市熱島這一概念詞[3]。城市熱島是指城市中心溫度明顯高出周邊郊區溫度的現象[4]，這一現象的出現伴隨著城市病的產生，例如空氣污染、城市高溫等，不僅影響人類的正常工作學習，還會危害身體健康。

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

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

1 研究區概況

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

圖1 研究區概況圖Fig.1 Survey map of the study area

2 數據及處理方法

2.1 數據源

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

表1 遙感數據參數Table 1 Remote sensing data parameters

2.2 地表溫度的反演

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

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

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

(1)

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

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

(2)

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

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

(3)

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

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

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

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

(4)

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

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

ε=0.004Pv+0.986

(5)

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

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

(6)

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

2.3 熱島強度的計算

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

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

(7)

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

3 結果與分析

3.1 地表溫度變化分析

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

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

表2 2004—2020年熱島強度面積統計Table 2 Statistics of heat island intensity area from 2004 to 2020

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

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

3.2 土地利用格局演變

基于支持向量機對研究區進行分類，土地利用變化明顯，呈現主城區向外蔓延擴張[22]，結果如圖3所示。結合轉移矩陣(表3～表5)，可以看出2004—2020年石家莊市主城區的土地利用的變化。建筑用地面積增加了552.77 km2，其中大部分由耕地和裸地及其他轉變而來，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年的時間，耕地減少了549.75 km2，裸地及其他減少了134.67 km2；林地在2004—2010年多被砍伐或用來耕作等造成面積減少，到2016年增加到108.40 km2，2020年達到了244.70 km2，這與2015年啟動的太行山生態綠化工程和2018年的重點造林綠化工程有很大關系；另外石家莊市是一個缺水型城市，因此水體面積變化微小，并且總面積小。石家莊市城市化進程，大量耕地、裸地及其他等被人為建筑取代，建筑用地大量增添，城市熱島嚴重加劇。

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

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

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

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

3.3 熱島擴張的方向特征

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

表6 熱島變化標準差橢圓信息表Table 6 Ellipse information table of heat island variation standard deviation

表7 建筑用地變化標準差橢圓信息表Table 7 Ellipse information table of standard deviation of building land change

圖4 熱島變化標準差橢圓圖Fig.4 Standard deviation ellipse of heat island variation

圖5 建筑用地變化標準差橢圓圖Fig.5 Standard deviation ellipse of heat island variation

向西北，主城區與鹿泉區之間交通路網發達，與其他主城區通行便捷。鹿泉在持續發展的進程中，建筑設施也逐漸完善發展，然而，鹿泉區位于石家莊的“上風口”，也是重要的生態區，主城區向西持續發展的空間會在未來受限。

向東南，隨著高新區的開發，以及高新區與主城區之間緊密的城區發展格局，將主城區向此擴張，并將太行大街以東、南二環東延線以南的區域作為“處女地”不斷開發，東南方向將成為未來主要繼續延伸發展的方向。

3.4 城市熱島效應與下墊面關系

3.4.1 溫度與NDVI關系

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

圖6 地表溫度與NDVI的相關性Fig.6 Correlation between surface temperature and NDVI

3.4.2 溫度與NDBI關系

NDBI建筑信息提取公式為

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

(8)

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

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

圖7 地表溫度與NDBI的相關性Fig.7 Correlation between surface temperature and NDBI

3.4.3 溫度與TVDI關系

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

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

(9)

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

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

Tmin=a1+b1NDVI

(10)

Tmin=a2+b2NDVI

(11)

式中：a1、b1、a2、b2為方程系數。

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

圖8 地表溫度與TVDI的相關性Fig.8 Correlation between surface temperature and TVDI

4 結論

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

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

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

(3)建筑用地是城市地表升溫的重要因素，植被覆蓋度和水體能夠對地表溫度調控起到有效作用。可以通過改變城市格局及建筑材料來降低城市熱島效應，比如使用太陽光反射率好或能降溫節能材料等[25]，做到“適地適樹”的原則，選擇適合植被種植，增加城市綠地面積，城市定時灑水、建造公園濕地等方法解決石家莊的熱島效應問題[30]。