摘要: 城市气候环境和下垫面性质对人类健康具有重要的影响. 为了研究福州主城区不透水面与城市热岛效应之间的时空变化关系，本研究以 1994 年、2000 年、2009 年、2018 年 4 期 Landsat 遥感影像数据作为数据源，基于线性光谱分解原理提取福州主城区不同时期不透水面信息和利用单窗算法反演出地表温度，分别从空间维度和时间维度进行对比分析，结果表明: 不透水面与城市热岛的分布具有显著的正相关关系，R2 达到 0. 70 以上. 1994 年 - 2018 年，福州主城区不透水面增加了 89. 58 km2 ，年平均增速为 3. 73 km2 / 年. 就福州主城区而言，不透水面丰度每增加 0. 1，可导致局部地表温度上升 1 ℃ 以上，并且不透水面的空间分布与热岛的分布特征具有较强的一致性.
关键词: 不透水面; 热岛效应; 遥感; 福州
城市热岛效应是指城市因大量的人工产热、建筑物、厂房和道路等高蓄热体的增多及绿地和植被的减少等因素，造成城市 “高温化”，城市内部的气温明显高于外围郊区的现象[1]. 城市热岛是城市热环境最常见的气候分布特征，其中城市下垫面性质和城市热岛的形成密切相关，类型和结构的不同通过影响城市局部热环境的变化，进而影响到城市的气候环境和人类的健康生活[2]. 下垫面可根据透水性简单分为不透水面 ( impervious surface area，ISA) 和透水面两类. 文中不透水面是指人为制造的能够阻止水体渗透到地表以下的人工地貌类型[3]，而从生态学的角度来理解，则通常指那些相对于植被和土壤渗透率更小的建筑区域 ( 例如屋顶、公路、人行道、停车场等) [4]. 随着城市化的不断发展，土地城市化也日渐显现，特别是耕地或水体向建筑用地的转变，使原本不易蓄热的地面表层转化为温度更易升高的不透水层，从而导致了城市热岛效应的不断加剧[5].
目前，遥感技术在很多方面得到广泛应用，其中研究利用卫星对地观测影像如何快速准确地获取不透水面受到国内外众多学者的关注. 如 RIDD[6]提出了研究地表覆盖的植被-不透水面-土壤模型 ( vegetable-impervious-soil model，简称 V-I-S 模型) ，认为城市的地表覆盖类型是由植被、不透水面、土壤组成或是 3 种地物的线性组合. CARLSON 等[7]利用植被覆盖度与不透水面的关系，提出适合城市建成区不透水面信息提取的方法. YANG[8]基于回归思想，利用不同季节 Landsat 影像各波段和穗帽变化结果的多种组合分别作为分类回归树的输入，分别在不同的空间尺度上提取了不透水面盖度，结果表明平均误差在 8% ～ 11. 4% 之间. WU[9]为减少同种地物端元光谱亮度差异而产生不透水面盖度提取误差，发展了亮度归一化的线性光谱混合分解 ( LSMA) 模型，研究结果表明不透水面盖度估算的标准误差为 10. 1% . 徐涵秋[10 - 11]提出利用归一化差值不透水面指数 ( normalized difference impervious surface index，NDISI) 提取不透水面的方法，并将其用于快速、自动提取大区域范围内的不透水面信息，结果具有很高的精度. 林云杉等[12]利用植被覆盖度与不透水面的负相关关系，提取并研究了泉州地区不透水面与城市热岛效应之间的关系. 岳文泽等[13]基于混合光谱分解对上海城市不透水面分布进行估算，结果表明上海市不透水面分布比率高，不透水面分布的空间差异进一步揭示了城市土地覆被空间结构以及城市空间扩展的差异性.
本文针对 Landsat 影像数据特点，选择线性光谱分解模型提取不透水面和单窗算法反演地表温度，研究福州市主城区不透水面与热岛效应的空间分布格局及时间变化，并使用数学方法分析二者之间的关系.
1 研究区概况
福州 市 ( 25° 15' ～ 26° 39' N， 118°08' ～ 120° 31' E) 是福建省的省会 ( 见图 1) ，位于中国东南沿海、福建省中东部的闽江口，与台湾省隔海相望. 地貌类型属典型的河口盆地，盆地四周被群山峻岭所环抱，境内地势自西向东倾斜，其海拔多在 600 ～ 1 000 m 之间，该地区属典型的亚热带季风气候，年平均降水量为 900 ～ 2 100 mm，年平均气温为 20 ～ 25 ℃，全市陆地总面积 11 968 km2 ，其中市区面积 1 786 km2 ，2016 年末全市常住人口 757 万人.
2 基本原理和方法
2. 1 数据来源及预处理
本文在地理空间数据云网站上获取 4 期 Landsat 遥感影像，影像数据信息见表 1. 由于原始影像为 DN 值，需要对原始影像进行一系列的预处理，具体包括几何精校正、辐射定标和 FLAASH 大气校正和，由于研究区内水体的反射率较低，会对不透水面的提取带来干扰，因此在提取不透水面的过程中要对水体进行掩膜处理. 水体掩膜的方法本文采用修正后的归一化水体指数 ( MNDWI) .
2. 2 建成区不透水面信息提取
2. 2. 1 线性光谱分析原理
影像中一个像元的覆盖范围包含多种地物，像元的光谱特征其实是像元内几种纯净地物光谱值的组合，这样的像元称为混合像元. 混合像元内的这些纯净地类被称为端元 ( endmembers) . 目前，比较流行的光谱混合模型有线性、概率、几何光学、随机几何及模糊模型等，线性模型因其简单和原理清晰受到广泛欢迎，首先假设每一光谱波段中单一像元的反射率，为各端元组分特征反射率与它们各自比率的线性组合，其数学模型如式 ( 1) 所示[14]: Ρ( λi ) = ∑ Fj Ρj ( λi ) + ε( λi ) ， ( 1) 式中，i = 1，…，m( 光谱波段数) ; j = 1，…，n( 端元数) ; λi 为第 i 波段的波长; Ρ( λi ) 为第 i 波段像元反射率，其包含一种或多种端元组分; Fj 为一个像元内第 j 个端元的反射率所占比率; Ρj ( λi ) 表示一个像元内端元 j 在第 i 波段上的反射率; ε( λi ) 是第 i 波段的误差. 在求解端元所占比率Fj 的过程中，经常采用有约束条件的最小二乘方法，要求Fj 同时满足∑ n i = 1 Fj = 1 且Fj ≥ 0 这两个条件，模式的正确性还必须通过检验影像中每个波段的残差的均方根来确定，如式( 2) 所示: RMS = ( ∑ m i = 1 ε ( λ ) i 2 ) 槡 m ， ( 2) 其中 RMS 代表残差 ε 的均方根; m 是影像中的波段数. 2. 2. 1 最小噪声分离变换变换和端元选择最小噪声分离变换 ( Minimum NoiseFraction，MNF) 是 一种空间图位变换方法，是能对多光谱和高光谱遥感影像进行降维、消除波段间的相关性和压缩数据量 的 有 效 方 法[15]，它的 应用可以提高光谱混合分解精度[16]. 本研究以 2000 年 ETM + 影像为例进行 MNF 变换，生成 6 个变换分量 ( 见图 2) ，前 3 个变换分量空间纹理清晰，特征值占总 信 息 量 的 87. 71% ( 见 表 2) ，而后 3 个变换分量空间纹理模糊，噪声较大，占总信息量的比重较小. 因此选取前 3 个变换分量进行端元的选，研究选取 4 类端元，包括高反射率、低反射率、土壤和植被.
2. 2. 2 线性光谱分解和精度验证
在对研究区 ETM + 反射率影像进行线性分解后，得到一个包含 5 层数据的影像，前 4 层分别为高反照率、植被、土壤和低反照率类型在每个像元中所占面积比例，第 5 层是分解结果的均方差统计. 分解结果的均方差统计影像能很好地反映分解结果的准确性，能直接用来对分解结果进行评价，分解结果的总体均方根 ( RMS) 平均值必须小于 0. 02[17]. 根据解混中得到的 RMS 频率分布与空间图像，本研究区域 RMS 最大值为 0. 043 65，最小值为 0，平均值为 0. 007 35，绝大部分像元的 RMS 值都小于 0. 015，影像分解的 RMS 远小于 0. 02，达到了精度要求，从定量的角度说明了本次研究选取的终端地类数目合适，地类光谱值准确，分解精度高，分解结果可靠. 反演流程见图 3. 2. 3 地表温度反演地表温度反演算法有很多种，结合本文研究的数据类型 Landsat TM /ETM 数 据 以 及 Landsat8 OLI /TIRS 数 据，采 用 覃 志豪[18]的单窗算法反演地表温度流程如图 4.覃志豪单窗算法是根据地表热辐射方程推导出来的，计算公式如下: T = [a( 1 - C - D) + ( b( 1 - C - D) + C + D) TSensor - D Ta]/C ， C = τε ， ( 3)
式中，ε 是地表比辐射率，τ 是大气透过率， TSensor 是星上辐射亮度对应的亮度温度，Ta 是大气平均作用温度，a、b 为系数.辐射量温是假定地球为绝对黑体的条件下求得的，利用定标系数将影像 DN 值转换成星上辐射亮度( Ii ) 后，可用普朗克公式求解出星上亮温，计算公式如下: Ii = Ki，2 /ln( 1 + Ki，1 / Ii ) ， ( 4) 式中，Ki，1 和Ki，2 是常量，Ki，1 = 2hc2 /λ5 i， Ki，2 = hc /( k λi ) ，h 为 普 克 常 数，约 为 6. 626 068 96 × 10 -34 J·s，c 为光速约为 2. 997 92 × 108 m·s -1，k 为玻尔兹曼常数，约为1. 380 650 5 × 10 -23 J·K-1，λi 为第 i 通道内的中心波长. Landsat8 第 10 波段的中心波长λ10 = 10. 9 μm，根据 Planck 公式得出K10，1 = 774. 89 W·m-2·sr -1·μm-1，K10，2 = 1 321. 08 K. Landsat TM 第 6 波段的中心波长 λ6 = 11. 435 μm，根据 Planck 公式得出K6，1 = 607. 76 W·m-2·sr -1·μm-1，K6，2 = 1 260. 56 K. Landsat ETM 第 6 波段的中心波长 λ6 = 11. 435 μm，根据 Planck 公式得出K6，1 = 666. 09 W·m-2·sr -1·μm-1，K6， 2 = 1 282. 71 K. 因为本文研究区处于中纬度地区，遥感影像云量都在 1% 以下，由于缺乏过境时大气实测数据，1994 年选择大气水汽含量经验值 ω = 2. 0 g·cm-2 进行计算. 2001 年、2009 年、2018 年根据 NASA 网站 ( http: / /atmcorr. gsfc. nasa. gov /) 查询得到大气透过率分别为 0. 82、0. 90、0. 85.
根据覃志豪[19]等人提出的 NDVI 阈值改进算法求取地表比辐射率时除了考虑自然表面之外，还考虑了水面和城镇这 2 种地表覆盖类型. 通过监督分类提取水体，其它 2 类分为城镇像元和自然表面像元，其中组成自然表面的像元可以看做是比例不同的植被和裸土组成的混合像元，而城镇像元则是建筑物和植被组成的混合像元. 采用植被覆盖度公式计算获得，其估算公式如下: 自然表面像元 ε = PV RV εv + ( 1 - Pv) Rsεs + dε， ( 5) 城镇像元 ε = PV RV εv + ( 1 - Pv) Rm εm + dε， ( 6) 式中，PV 是植被覆盖度，RV、RS 和Rm 分别是纯植被、纯裸土和纯建筑表面像元的温度比率，εv、εs 和 εm 分别是纯植被、纯裸土和纯建筑表面像元的比辐射率，分别取 εv = 0. 987，εs = 0. 973，εm = 0. 971，RV = 0. 933 4 + 0. 058 2 PV，RS = 0. 991 2 + 0. 107 8 PV，Rm = 0. 988 7 + 0. 128 6 PV . dε 为地形形状，在地表相对平整情况下，一般可取 0，覃志豪根据植被的构成比例提出了计算地表形状的经验公式，当PV ≤ 0. 5 时，dε = 0. 003 8Pv; 当PV ≥ 0. 5 时，dε = 0. 003 8( 1 - PV ) .
3 结果和分析
3. 1 城市不透水面分布特征
利用线性光谱分解法，得福州市建成区在 1994 年、2000 年、2009 年、2018 年不透水面分布图 ( 见图 5) . 图中的颜色的差异表现为不透水面信息的强弱，红亮色代表较高数值，此时不透水面的信息较强. 反之，蓝绿色代表较低数值，此时不透水面信息较弱. 主城区不透水面大体上分布在鼓楼区、台江区、仓山区，并从 1994 - 2018 年增加明显，以 2018 年为例，数值在 0. 4 以下区域为植被区 ( 包括森林和耕地等) ; 数值在 0. 6 以上区域为城镇用地或农村住宅建筑用地，特别是新开发用地不透水面数值可达到0. 9 以上. 数值位于0. 4 ～ 0. 6 之间的区域，多为植被和不透水面的混合区. 由图5 可以看出增加区域主要为仓山区及福州主城区周边地区. 表 3 统计出研究区不同年份不透水面的面积及年际之间的差值.
3. 2 地表温度分布特征
根据单窗算法地表温度反演结果 ( 见图 6) ，由于本文的 4 景影像获取的时相不同，因此不能使用反演的绝对温度进行探究不同时相的温度变化. 为使 2 期影像温度数值尽量相同，首先将地表温度标准化，随后将其拉伸到 0 ～ 40 ℃之间. 表达式如下: T = 40( TS - Tmin Tmax - Tmin ) ， ( 7) 式中，Tmin 为地表温度的最小值; Tmax 为地表温度的最大值.地表温度作为衡量城市热环境的重要指标，可以用来模拟城市热环境的分布格局，并可初步了解城市热岛效应的分布特征. Landsat 卫星在福州地区过境时间为上午 10 点半前后，查找历史天气资料可知反演出的地表温度与真实温度相差不大，可用于分析与不透水面的关系. 其中水体表现为低温，植被要略高于水体，高温区主要分布在城市主城区.
本研究利用 Jenks[20]提出的自然断点法 ( Natural Breaks) 分类. 分类的原则就是差不多的放在一起，分成若干类. 统计上可以用方差来衡量，通过计算每类的方差，再计算这些方差之和，用方差和的大小来比较分类的好坏. 因而需要计算各种分类的方差和，其值最小的就是最优的分类结果 ( 但并不唯一) . 利用 Arcgis10. 2 软件中的 Natural Breaks 分类工具将地表温度分为 7 个等级，分别为低温、次低温、次中温、中温、次高温、高温、特高温 ( 见图 7) .
在研究区内，茂密植被覆盖区表现为低温和次低温，城市居住区表现为中温及以上. 总体上，低温和高温分布都很集中，其中高温大部分集中在建成区，由城区向郊区逐渐减弱过度到低温，福州城区位于研究区的中间，也是高温最为集中的区域. 其中仓山区 ( 图中岛状位置) 从 1994 - 2018 年地表温度变化强度和面积最大，地表温度整体表现出东进南拓的趋势，与不透水面的变化趋势大致相同. 表 4 统计各个温度等级所占面积及百分比。
3. 3 城市热岛分布特征
地表温度作为衡量城市热环境的重要指标，可以用来模拟城市热岛分布格局，并对热岛效应的分布特征进行分析. 热岛强度是反应城市温度差异的一个指数，可用地表温度计算得出，并对其进行阈值分类，得到不同强度的城市热岛的分布[21]，如下: THI = ( T - Tmean ) /Tmean， ( 8) 式中，T 为地表温度，Tmean 为研究区域地表温度的平均值，THI 为热岛强度. 由于本文的4 景影像获取的时相不同，不能直接将式( 8) 计算结果进行分级，本研究将计算结果归一化处理，并分为 6 个等级分别为强绿岛、中绿岛、弱绿岛、弱热岛、中热岛和强热岛( 见图 8) .
在研究区内，茂密植被覆盖区表现为高等级的绿岛，城市居住区表现为中级及以上热岛. 总体上，热岛和绿岛分布都很集中，其中热岛大部分集中在建成区，由城区向郊区逐渐减弱过度为绿岛，福州城区位于研究区的中间位置，也是热岛最为集中的区域. 其中仓山区在 1994 - 2018 年热岛强度变化最显著、面积增加最多. 表 5 为统计各个热岛等级所占面积及百分比.
3. 4 不透水面与城市热岛关系研究
城市不透水面的存在导致下垫面环境的变化，进而导致城市局部热环境的改变，主要体现在温度的明显升高，通过对比图 5 和图 6，可以发现城市热岛空间分布和变化与不透水面的空间分布和变化有明显的一致性. 为了研究其相互关系，利用求得的 1994 - 2018 年不透水面与温度进行相关性分析，从定量的角度分析两者之间的关系. 结果表明，不透水面分布与地比温度具有显著的正相关关系，R2 均在 0. 70 以上. 表明不由水面率较大的区域温度也较高，从二者的回归方程( 图 9) 可知，不透水面率每 加 0. 1，地表温度分别升高 1. 56℃ ( 1994 年) 、1. 52℃ ( 2000 年) 、1. 58℃ ( 2009 年) 、1. 49℃ ( 2018 年) . 由此可见两者关系密切，大面积的城市不透水面的增加会引起城市热环境的变化，进而产生了热岛效应.
4 结论
研究结果表明: ( 1) 福州市主城区不透水面面积从 1994 - 2018 年 24a 间增加了 89. 58 km2 ，平均每年增加 3. 73 km2 . 增加区域主要是福州市南部的仓山区和靠近主城区的周边地区. ( 2) 不透水面的增加会引起城市局部地表温度的升高，且具有显著的正相关关系，R2 达 0. 70 以上，不透水面的分布特征和发展趋势与城市热岛的分布规律具有较强的一致性. 就福州市而言，不透水面每增加 0. 1，就会导致局部地表温度上升 1℃以上. ( 3) 福州主城区主要位于闽江沿岸，城市北面有山脉阻隔，地势相对较陡，城市扩张困难. 城市南部仓山区附近地势相对平缓，适宜城市的发展，东部靠近沿海发展相对较快，这也符合福州东进南拓的发展方针，造成了城市南部不透水面的增加进而导致热岛效应更加明显. ( 4) 本文针对 Landsat 系列影像数据，采用覃志豪单窗算法反演地表温度分析热岛效应是可行的，这对福州以后的城市发展规划具有一定的参考意义.
参考文献:
[1] 李延明，徐佳，张济和，等. 城市绿化对北京城市热岛效应的缓解作用 [J]. 北京园林，2002，18 ( 4) : 12 - 16.
[2] 薛志成. 城市热岛效应威胁人类健康 [J]. 安全与健康，2002 ( 13) : 15 - 16.
[3] 徐永明，刘勇洪. 基于 TM 影像的北京市热环境及其与不透水面的关系研究 [J]. 生态环境学报，2013，22 ( 4) : 639 - 643.
[4] 徐涵秋. 城市不透水面与相关城市生态要素关系的定量分析 [J]. 生态学报，2009，29 ( 5) : 2456 - 2462.
福州主城区热岛效应与不透水面的关系及时空变化分析相关论文期刊你还可以了解：《闽都文化福州市文化期刊》
