两种多时相遥感地表温度正规化处理方法对比分析

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

(1.昆明理工大学 国土资源工程学院， 云南 昆明 650093；2.中国科学院 深圳先进技术研究院， 广东 深圳 518052)

城市热岛(Urban Heat Ialand，UHI)是指由城市化引起的城市地表及大气温度高于周边非城市环境的一种现象[1]。随着全球城市化进程的不断加快，近几十年城市热岛现象不断加剧，引起了社会的广泛关注[2-3]。当前城市热岛研究通常采用气象观测、遥感图像温度反演等方法，对热岛的时空分布及演变规律进行研究。其中，气象观测的方法是基于有限的地面气象站数据，覆盖范围有限且分布不均匀，很难全面地反映城市热岛的空间分布情况。而遥感图像具有覆盖范围广、时效性强、数据综合性和可比性强等优势，使其在城市热岛的研究中得到了越来越广泛的应用[4]。

为了开展城市热岛效应演变研究，需要将不同时相遥感图像的热红外波段加以比较[3]。然而，不同时相遥感图像在获取时，大气条件、传感器辐射定标、太阳辐射强度等的不同，导致图像中同一目标地物具有较大的辐射差异[5]。即便在数据预处理过程中已经普遍采用FLAASH模型或ATCOR模型等方法对获取的遥感图像进行辐射校正[6]，时相差异却仍然对地表温度反演结果产生着不同程度的影响。如同一地点地表覆盖不变，反演得到的温度却有较大差异，即地表温度的“伪变化”[7]，从而对热岛现象演变分析产生较大干扰。因此，要想进行基于遥感的城市热岛的动态变化分析，就必须考虑对多时相遥感地表温度图像的正规化处理。将不同时相的地表温度图像进行正规化处理后，解决客观条件对地表温度反演的干扰问题，从而进一步研究时间序列的城市热岛演变。

目前，城市热岛演变研究中不同时相地表温度图像正规化处理方法效果较好，被广泛应用的主要有两种：典型归一化方法[3,8]和伪不变特征(Pseudo Invariant Feature，PIF)法[9]。前者由徐涵秋于2003年首次应用于多时相地表温度图像的归一化处理，后被学者广泛应用于城市热岛演变的相关研究。后者由Schottetal于1988年首次提出，目的是使多时相图像具有相同的辐射尺度(或辐射分辨率)。多用于对遥感影像的相对辐射校正，目前有少数学者将其应用于多时相地表温度处理[10]，并取得了良好效果。

深圳市近20年来经济、文化等飞速发展，使城市建成区扩展的速度迅速，规模较大，同时也使得其城市热岛现象明显加剧。因此，本文以深圳市陆域范围为例，选用Landsat系列数据，在不同年份、不同季节的时间序列上城市热岛效应演变研究中，讨论上述两种方法对地表温度图像正规化处理效果的差异。

1 研究区与数据

深圳市位于珠江三角洲东岸，陆域范围22°26′59″—22°51′49″N，113°45′44″—114°37′21″，东临大亚湾与惠州相邻，西至珠江口，南望香港特别行政区，北接莞、惠两地，整体呈狭长型。需要说明的是，本文所提到的深圳市均表示深圳市陆域范围，不包括海域。深圳市作为快速城市化的典型代表，30余年间从最初的边陲小镇发展为下辖6个行政区和4个新区内的60个街道办事处的国际化都市。城市空间结构以包括福田、罗湖、南山3个行政区在内的中心城区为核心，形成“南北贯通，西联东拓”的区域空间协调发展框架。

深圳市处于121/43、122/44共2景Landsat图像上，121/43包含深圳东半部分区域，122/44包含西半部分区域。本文选用了Landsat 5 TM、Landsat 8 OLI/TIRS数据，在2000、2005、2010、2015年4个时间段各获取的2景遥感图像，云量均为0，共获取8景，数据信息如表1所示。由于中国南部沿海常年多云多雨，深圳为典型区域，Landsat数据缺失、质量差的现象尤其明显，同一地区2年内可能只有一景图像可用。因此，本文选取的遥感图像日期并非全为当地同一季节。基于此，就更需要对地表温度图像用正规化处理方法进行分析探讨，解决客观条件对地表温度反演的干扰问题。

表1 Landsat数据信息

另外，本文数据处理过程中还选用了深圳市2000、2005、2010、2015年4期的土地利用分类产品，由中国科学院深圳先进技术研究院空间信息研究中心提供。该数据产品的生产基于Landsat系列遥感卫星，其中1980年基于80 m分辨率Landsat MSS数据，1990、2000、2010年基于30 m分辨率Landsat TM数据，2015年基于30 m分辨率Landsat 8 OLI数据。数据来源于美国地质勘探局(United States Geological Survey，USGS)网站(https://landsatlook.usgs.gov)。数据产品的生产过程：(1)基于ENVI对所有图像进行辐射定标，消除传感器误差；(2)基于FLAASH模块进行大气校正；(3)基于易康软件建立分类规则集，进行监督分类；(4)以人机交互的方式对分类结果进行修正。该数据产品的分类系统包括建设用地、草地、耕地、林地、湿地、未利用地和近海海域，分类精度达到95%以上。

2 研究方法

2.1 地表温度反演

本文选用覃志豪的单窗算法[10]，基于ENVI 5.1软件平台，进行Landsat 5、Landsat 8图像地表温度反演，得到地表温度图像。普朗克函数用于定量描述黑体辐射的能量与波长(频率)、温度的客观规律，已知物体的能量谱和辐射强度，即可推算出亮度温度。

根据普朗克函数将黑体辐射亮度转化为亮度温度[10]，公式为

(1)

式中Tb是亮度温度(K)；Lb为同温度下的黑体辐射亮度，由对热红外波段进行辐射定标后的绝对辐射亮度L和地表比辐射率ε计算得到；k1、k2为热红外波段的参数，Landsat 8第10波段的参数为k1=774.89、k2=1 321.08，Landsat 5第6波段的参数为k1=607.76、k2=1 260.56。

得到亮度温度Tb后，用覃志豪的单窗算法[10]计算地表温度(Land Surface Temperature，LST)，公式为

(2)

式中LST为地表温度(℃)；a、b为常数，地表温度在0～70 ℃内，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为气象局发布的地表1.5～2.0 m左右的地面温度，即气温(K)；C、D为中间变量，C=ετ，D=(1-τ)×[1+(1-ε)τ]，ε为比辐射率，τ为大气透射率。

将本研究选取的8景Landsat图像分别通过上述方法，计算得到8幅地表温度图像。然后将表1中同一年份的2景地表温度图像进行地表温度的融合拼接、裁剪[11]等一系列处理，最后得到深圳市陆域范围2001年3月1日、2005年11月23日、2009年11月2日、2015年10月18日的地表温度图像。下文提到的2000、2005、2010、2015年分别指上述日期的LST图像。

此外，本文在反演地表温度的过程中，同时获取了深圳市归一化植被指数(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 地表温度图像正规化方法

2.2.1 典型归一化方法

归一化是一种无量纲处理手段，在城市热岛的研究中，使用此方法主要针对研究热岛现象相对强弱的空间分布。当侧重地表温度相对高低的空间分布特征时，时间不同改变的只是地表温度值的大小，并不改变地表温度的空间分布[3]。因此，徐涵秋等采用典型归一化方法对反演得到的LST图像进行处理，令不同时相的LST图像分布范围统一到0～1之间。计算公式为

(4)

式中Ni为第i个像元温度归一化后的值，TSi为第i个像元地表温度值，TSmin为LST图像像元灰度值的最小值，TSmax为LST图像像元灰度值的最大值。将4幅LST图像分别经过公式(4)处理后，像元灰度值范围均被归一化至0～1之间，此时的图像即为典型归一化方法正规化处理后的图像。

2.2.2 伪不变特征法

伪不变特征法(Pseudo Invariant Feature，PIF)在遥感图像处理中也被称为相对辐射校正(或相对辐射归一化)方法，是基于同一区域，多时相图像相同波段的像元灰度值存在线性关系的假设提出来的。首先要确定参考图像和实验图像，然后选取伪不变特征点(或称为辐射地面控制点)，确定伪不变特征点后，提取两期图像伪不变特征点的像元灰度值，建立参考图像与实验图像之间的回归模型[8]。将多时相图像的像元灰度值逐渐归一化到参考图像，使多时相图像具有相同的辐射尺度(或辐射分辨率)。该方法的具体步骤如下。

2.2.2.1 选取参考图像

参考图像应该目视质量较好，变动范围较大，以便于信息量尽可能多地保留下来[14]。本文共有2000、2005、2010、2005年共4个时间的地表温度图像，前3幅由Landsat 5数据反演得到，后1幅由Landsat 8数据反演得到。依据文章2.1小节，4幅LST图像的具体日期分别为2001年3月1日(春)、2005年11月23日(秋)、2009年11月2日(秋)和2015年10月18日(夏)，深圳市一年四季植物长势良好，春秋季与夏季图像的目视质量不相上下。

图1为上述4幅LST图像的像元灰度值曲线图。从图中可以看出，2015年的像元灰度值变动范围最大，2010和2000年与之相近，参考图像可从这3幅图中选取；另外，同一季节下图像的地表温度“伪变化”会相对较小。因此，本文选取2010年的图像作为参考图像，其余3幅图为实验图像。

图1 LST像元灰度值曲线图

2.2.2.2 选取伪不变特征点

伪不变特征点是指不同时相的遥感图像中，地表不随时间发生变化或变化很小可以忽略不计的地物，深水体(如水库)、常年裸岩裸土、大屋顶等典型地物[7]。对于样本点的选取，样本容量应当适中，样本容量过大或过小都会直接影响回归模型的结果。通常伪不变特征点的选取主要为人工目视的方法，得到的点数有限，且由于个体差异导致特征点的质量水平参差不齐。本文选取伪不变特征点的过程主要分为两个步骤：第1步，基于深圳市4期的土地利用分类产品，分别提取2000、2005、2015年3幅实验图像与2010年的参考图像之间土地利用未发生变化的区域，在土地利用不变的区域上，再分别提取3000个随机点，记为S1、S2、S3；第2步，依次提取S1在2000年和2010年遥感图像的NDVI值S1a和S1b、S2在2005年和2010年遥感图像的NDVI值S2a和S2b、S3在2015年和2010年遥感图像的NDVI值S3a和S3b。选取伪不变特征点Pi的规则如下：

Pi∶|Sia-Sib|<0.06，

(5)

式中i=1、2、3，满足上述条件的随机点，则为选取的伪不变特征点。这样得到的伪不变特征点点数适中，分布均匀，且具有客观、定量化的选取标准。

2.2.2.3 建立回归模型

回归模型是以参考图像伪不变特征点的像元灰度值为纵坐标y，以实验图像伪不变特征点的像元灰度值为横坐标x进行点绘，然后对这些点做回归分析，求解模型参数。首先从已有的LST图像中提取伪不变特征点的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，参数α和β采用最小二乘法求解。计算得到3组模型参数，分别代入2000、2005、2015年的LST实验图像进行计算，得到的3幅图像即为伪不变特征法正规化处理后的图像，2010年LST图像为参考图像，不需要再做处理。

以2010年LST图像为参考图像，利用上述方法提取出的3组伪不变特征点P1、P2、P3，建立3组伪不变特征回归模型，如表2所示。

表2 伪不变特征法回归模型

表2中3组回归模型均显示出了较强的相关性，确定系数0.5

3 城市热岛演变分析

3.1 热岛强度等级划分

城市热岛的研究中，有关地表热场分布、热岛结构、热岛演变等的研究，都涉及城市热岛的界定问题，即地表温度的等级划分[14]。常用的灰度图像等级划分方法有等间距法、分位数划分法、自然间断点分级法、几何间隔法以及均值-标准差法等。本文综合各分级方法对2010年LST图像典型归一化方法正规化和伪不变特征法正规化的分级效果，采用目视判别，发现均值-标准差法划分等级，可视化处理的效果最佳。因此本文对正规化后的LST图像分别以μ(均值)-σ(标准差)设置分割点，将热岛强度等级划分为6个级别：低温区、次低温区、中温区、次高温区、高温区和特高温区。等级划分标准如表3所示，表中Ts为需要划分等级的图像。

表3 均值-标准差法划分热岛强度等级

将分级结果按照时间顺序排列，如图2所示。对图2进行目视分析，观察两种正规化方法处理过后的热岛强度等级划分结果，可以看出：

图2 两种正规化处理方法下的热岛强度等级划分结果

(1)2000—2015年间，两组结果均体现了深圳市热岛现象随时间逐渐增强的现象，即“特高温区”逐年增多、“低温区”逐年减少。且两者各个热岛级别的空间分布移动特征相似：高温区、特高温区均表现出了随时间逐渐往中西部、西北和西南部转移的特征，其中西部与西北部热岛聚集现象最为严重。

(2)两组结果在不同年份的热岛强度增强程度存在差异。经过典型归一化方法处理的热岛强度等级划分结果中，2000、2005、2010年的热岛转移趋势明显，但增强程度相对不明显；2010—2015年的热岛增强程度相对较明显。经过伪不变特征法处理的热岛强度等级划分结果中，2000、2005、2010、2015年各阶段热岛转移趋势和增强程度都较为明显。以2010年为分界，典型归一化法2000、2005年的热岛强度比伪不变特征法的强；前者2015年的高温区与特高温区分布与后者相差不大，低温区分布明显大于后者。

两组结果在不同年份的热岛强度增强程度存在差异，表明典型归一化方法和伪不变特征法对多时相遥感地表温度图像的处理效果不同，也就是说，对消除地表温度的“伪变化”的作用有差异。在本研究中，两种方法孰优孰劣，还需进一步对比分析。

3.2 城市热岛演变与土地利用变化关系

相关研究表明，城市环境效应产生的内在驱动力是土地利用与土地覆被的变化，城市热环境与土地利用类型密切相关。因此，本文列出4个时间段深圳市不同土地利用类型的分布情况，如图3所示，以进一步对两种方法的处理效果进行分析评价。

图3 深圳市2000—2015年土地利用类型分布

分别统计两种正规化方法处理过后的LST图像上各热岛强度等级的面积及所占比例，统计结果如表4、表5所示。统计深圳市2000—2015年不同土地利用类型的面积，如表6所示。

表4 深圳市2000—2015年热岛强度等级面积——典型归一化方法

表5 深圳市2000—2015年热岛强度等级面积——伪不变特征法

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

依据图3、表6土地利用类型分布的变化，可见深圳市人工表面与林地占地面积最多。2000、2005、2010年间人工表面占地面积大量增加，2010—2015年增加不明显；2000—2015年林地、耕地和湿地面积不断减少，草地占地面积少且缓慢增加。结合图2、图3进行目视对比，整体来看，热岛强度级别为高温区及以上的区域变化与人工表面的分布变化接近，低温区的变化与林地的分布变化接近。已有研究[15]表明，人工表面的占比对热岛现象的强弱有很大影响，而林地对城市热岛效应具有减缓作用。因此，本文依据表4—表6，以人工表面与林地为主，将高温区与特高温区作为热岛区域、次低温区与低温区作为热岛缓解区域，绘制不同区域面积变化与土地利用变化对比折线图(图4)。

(a) 热岛区域 (b) 热岛缓解区域

对比图4(a)热岛区域和人工表面的面积变化趋势，随着人工表面的显著增加。经过伪不变特征法处理的热岛区域面积变化趋势与人工表面面积变化趋势相同；而经过典型归一化方法处理的热岛区域面积却并没有明显变化，2015年热岛区域面积甚至少于2010年，这是不合理的。

对比图4(b)热岛缓解区域和人工表面的面积变化趋势，林地面积在2000—2015年间不断减少。经过典型归一化方法处理的热岛缓解区域面积在2010年之前的变化趋势与林地变化相近，但相反地，在2015年热岛缓解区域面积出现回升，整体表现为2000年与2015年几乎没有变化。经过伪不变特征法处理的热岛缓解区域面积在2000—2015年间表现为先减少、再增加、又减少，但整体是减少的，在2000年与2015年两个时间节点与林地面积减少量相对应；对于2005年热岛缓解区域面积的大幅减少，结合表5，发现该现象正好对应了2005年中温区面积的大幅增加。

从以上两方面来看，在本研究中，经过伪不变特征法处理过后的效果更好。

4 结语

实际上，两种正规化方法在城市热岛研究中都已有大量的论证和应用，各自都具有其运用价值。通过本研究，说明了典型归一化处理方法对年际变化的城市热岛效应演变研究适用性有限，特别是针对深圳市这样发展迅速、城市环境年际变化剧烈的地区，采用伪不变特征法处理后所得的结果相对较好。宏观上，从目视分析的结果来看，经过典型归一化方法处理所得的热岛等级划分结果在热岛分布、热岛的空间转移方面，与经过伪不变特征法处理所得的结果差别不大。但微观上，在城市热岛演变与土地利用变化关系的量化分析中，两组结果差异明显，伪不变特征法处理得到的结果更合理。

本研究也存在不足之处：(1)较典型归一化方法而言，伪不变特征法处理过程相对繁琐。若是针对日期相近、或一定时间范围内研究区土地利用类型变化不大的地表温度图像，用典型的归一化方法处理结果同样准确，且更方便。(2)伪不变特征法的精度较难以控制，特别是伪不变特征点的选取。目前常采用人工目视选点的方法，人工判断失误不可控且不适合较大区域。本文虽利用已有的土地利用分类产品设立了客观、定量的选取方法，但回归模型精度有所下降，方法有待改进。