石家庄市主城区热岛特征及变化分析

李外宾， 汤 军， 高贤君, 李志远

(长江大学地球科学学院，武汉 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]。