三亚市热异常遥感监测数据集 中华文库

摘要：针对三亚市热异常问题，本文以2008–2017年Landsat数据为主要数据源，基于辐射方程传输法反演地表温度和改进的箱线图法提取高温异常区，进行高温异常区频率统计，将频率大于60%的高温异常区判定为城市热异常区，得到10年间城市热污异常区。本数据集可应用在其他城市，进行多区域城市热异常区提取，从而发挥更高的应用价值。

关键词：热异常；辐射方程传输法；改进的箱线图法；高温异常区频率

Abstract: Aiming at the problem of thermal anomaly in Sanya, this paper takes Landsat data from 2008 to 2017 as the main data source, extracts the high temperature anomaly area based on the radiation equation transfer method and the improved box-plot method, carries out frequency statistics of the high temperature anomaly area, determines the high temperature anomaly area whose frequency is more than 60% as the urban thermal anomaly area, and obtains 10-year urban thermal pollution anomaly area. This data set can be applied to other cities to extract multi-regional urban thermal anomaly areas, so as to play a higher application value.

Keywords: thermal anomaly  ; radiation tra; fer equation; modified boxp; t method; frequency of; e high-temperature anomaly area

热异常是指由自然因素和人类活动中的热排放导致环境温度异常升高的现象^[1]。城镇化进程的快速推进可带动经济发展，但同时也对水资源、局部气候、生存环境等产生了一系列负面效应。据统计，近100 a来，全球的年平均气温上升了0.7–1℃，而大城市的平均气温则上升了2–3℃，北京、上海、东京等世界性大都市的年均气温比郊区高出2℃，夏季高出6℃左右^[2]。与物质污染带给环境的影响不同的是，物质污染是由污染物引起的，可以通过具体方法和措施得到控制，而热异常则不同，它主要来源于热源，有能源消耗则会产生热异常^[3]。因此，热异常是不可避免的。针对此类现象，国内外学者已经做了大量探索和研究，例如Wu和Chen等通过监测不同地区的核电站热流排放情况评估热异常强度和对周边环境的影响^[4][5]。Liu和Xia等基于热异常产品，通过利用不同方式提取工业热源^[6][7]；Zhang等利用Landsat数据对首钢搬迁前后的热环境进行探索和分析^[8]。还有一部分学者利用城市热岛强度与气温间的关系评估城市热异常^[9][10]。可见，热异常已被人们广泛关注，并成为研究热点。但以往关于热异常的大量文献多针对工业城市，或有工业污染的地方，针对非工业城市的热异常研究还较少。而热异常不仅存在于工业城市，还存在于非工业城市。因此，针对非工业城市，如何有效快速提取热异常区，如何缓解改善热异常具有重要研究意义。

三亚市作为我国唯一一个热带滨海旅游城市，其生态环境保护尤显重要。恰逢海南省推进“国际旅游岛”建设，作为海南旅游发展龙头的三亚正以国际化标准致力于打造“一港两地”的国际性滨海旅游城市。但伴随着三亚市城镇化的不断发展，城市热异常对城市环境的影响日益突出。因此，针对三亚城市建设和区域、环境、植被等特点，开展具有针对性的城市热异常研究，对三亚城市建设、规划布局、环境保护均具有重要意义。

综上，针对如何监测三亚市热异常区，并提取热异常区等问题，本文基于长时间序列的Landsat数据，利用辐射方程传输法获取地表温度信息，结合改进的箱线图法提取相应高温异常区，通过统计分析高温异常区的年度出现频率，得到城市热异常区。该方法不仅能够准确大面积获取城市热异常区，且相较于其他方法，可以有效避免潜在热异常的影响。

1.1 区域范围

三亚市位于北纬18°09′34″–18°37′27″，东经108°56′30″–109°48′28″，地处海南省最南端，主要以旅游业和农业为主，是海南省第二大城市。作为海南省最重要的旅游城市之一，其城市环境一直是广大群众关注的热点。但随着城市的不断发展与扩张，城市热异常现象日益加剧，因此，选取三亚市市辖区作为研究区，研究示意图如图1所示。

图1 三亚市边界示意图

1.2 基础数据准备

本文以2008–2017年间所有云量小于30%的Landsat 5、Landsat 7和Landsat 8数据为主要数据源（获取网址：https://glovis.usgs.gov）。本文温度反演所用到的波段及其空间分辨率如表1所示，其中，因Band11运行出现问题文中Landsat 8温度反演用到的是波段10。表2为本文所用影像的具体时间信息。

表1 温度反演所用波段信息

表2 数据集所用影像日期

1.3 数据处理

首先对Landsat数据进行拼接、裁剪、重投影等数据预处理。其中Landsat 7 ETM+影像因有条带，需在数据预处理基础上再进行去条带处理。其次用辐射方程传输法反演地表温度和改进的箱线图法提取高温异常区，得到高温异常区空间分布图。最后，为避免漏判现象，借鉴Liu等^[7]提取工业热源的思路，以频率大于60%为热异常区的标准统计分析热异常出现的频率，得到高温异常区频率空间分布图。同时，统计城市热异常区的总面积，具体数据处理流程如图2所示。

图2 数据处理流程图

1.3.1 地表温度反演

基于Landsat卫星影像，根据Landsat数据热红外波段的增益值与偏置值计算相应表观辐射亮度值；然后，利用影像植被覆盖度估算地表比辐射率；最后，根据Planck反函数和Landsat预设的定标系数计算地表温度。具体实现如下：

（1）表观辐亮度计算

表观辐亮度计算是将影像像元灰度值转换成相应热辐射强度的过程。利用ENVI软件中的辐射定标（Radiometric Calibration）工具将影像的像元灰度值转换为辐射亮度值。

（2）地表比辐射率估算

地表比辐射率表征地面物体的电磁波辐射的一种能力。基于影像的近红外波段和红光波段，根据公式（1）计算归一化植被指数NDVI^[11] 。 \[NDVI=\left({TM}_{\mathrm{N}\mathrm{I}\mathrm{R}}TM}_{\mathrm{R}\mathrm{e}\mathrm{d}}\right)/\left({TM}_{\mathrm{N}\mathrm{I}\mathrm{R}}+{TM}_{\mathrm{R}\mathrm{e}\mathrm{d}}\right)\]（1） 式中：\({TM}_{\mathrm{N}\mathrm{I}\mathrm{R}}\)为影像的近红外波段反射率；\({TM}_{\mathrm{R}\mathrm{e}\mathrm{d}}\)为影像的红光波段反射率。

其次，利用影像植被覆盖度估算地表比辐射率。由于地表物质结构具有差异性，因此针对不同地物类型，其相应的地表比辐射率计算方式也不同。大致可分为3种类型：水体、城镇和自然表面。而水体像元的地表比辐射率很高，与黑体的比辐射率较为相似，因此，在估算水体比辐射率时常被赋值为0.995，自然表面和城镇像元的地表比辐射率估算分别根据公式（2）和公式（3）计算^[12][13]： \[{e}_{\mathrm{S}\mathrm{u}\mathrm{r}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}}=0.9625+0.0614{P}_{\mathrm{V}\mathrm{e}\mathrm{g}\mathrm{e}0.0461{{P}_{\mathrm{V}\mathrm{e}\mathrm{g}\mathrm{e}}}^{2}\]（2）\[{e}_{\mathrm{B}\mathrm{u}\mathrm{i}\mathrm{l}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{g}}=0.9589+0.086{P}_{\mathrm{V}\mathrm{e}\mathrm{g}\mathrm{e}}-0.0671{{P}_{\mathrm{V}\mathrm{e}\mathrm{g}\mathrm{e}}}^{2}\]（3） 式中：\({e}_{\mathrm{S}\mathrm{u}\mathrm{r}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}}\)代表自然表面像元的地表比辐射率；\({e}_{\mathrm{B}\mathrm{u}\mathrm{i}\mathrm{l}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{g}}\)代表城镇像元的地表比辐射率；\({P}_{\mathrm{V}\mathrm{e}\mathrm{g}\mathrm{e}}\)代表植被覆盖度，计算方法如公式（4）所示^[14][15]： \[{P}_{\mathrm{V}\mathrm{e}\mathrm{g}\mathrm{e}}={\left[\frac{NDVINDVI}_{\mathrm{S}\mathrm{o}\mathrm{i}\mathrm{l}}}{{NDVI}_{\mathrm{V}\mathrm{e}\mathrm{g}\mathrm{e}{NDVI}_{\mathrm{S}\mathrm{o}\mathrm{i}\mathrm{l}}}\right]}^{2}\]（4） 式中：\(NDVI\)是影像的归一化植被指数，\({NDVI}_{\mathrm{V}\mathrm{e}\mathrm{g}\mathrm{e}}\)是影像中纯植被像元的NDVI值，\({NDVI}_{\mathrm{S}\mathrm{o}\mathrm{i}\mathrm{l}}\)是影像中纯土壤像元的NDVI值。当某个像元的\(NDVI\)大于0.70时，\({P}_{\mathrm{V}\mathrm{e}\mathrm{g}\mathrm{e}}\)取值为1；当\(NDVI\)小于0.05时，\({P}_{\mathrm{V}\mathrm{e}\mathrm{g}\mathrm{e}}\)取值为0；当\(NDVI\)介于0.05–0.7之间时，则分别将\({NDVI}_{\mathrm{V}\mathrm{e}\mathrm{g}\mathrm{e}}\)和\({NDVI}_{\mathrm{S}\mathrm{o}\mathrm{i}\mathrm{l}}\)取值0.70和0.05，并利用上式估算影像的植被覆盖度。

（3）地表温度计算

基于热红外波段的表观辐亮度和地表比辐射率计算结果，根据辐射传输方程反推地表热辐射亮度，如公式（5）所示： \[B\left({T}_{\mathrm{S}\mathrm{e}\mathrm{n}\mathrm{s}}\right)=\frac{LL}^{↑\tau \left(1-\epsilon \right){L}_{↓}}{\tau \epsilon }\]（5） 式中：\(B\left({T}_{\mathrm{S}\mathrm{e}\mathrm{n}\mathrm{s}}\right)\)是地表热辐射亮度；L是行星亮度温度值；\({L}^{↑}\)和\({L}_{↓}\)分别是大气上行辐射和大气下行辐射，可根据Landsat卫星过境时间、兴趣区中心经纬度、大气模式及传感器类型等获得，具体参数可通过NASA相关网页直接获得（http://atmcorr.gsfc.nasa.gov/）；\(\tau \)为大气路径透过率，可由NASA相关网页查询获取；\(\epsilon \)为地表比辐射率。

利用Plank反函数结合Landsat预设的定标常数K1 和K2 计算地表真实温度^[16][17]，如公式（6）所示： \[{T}_{\mathrm{S}\mathrm{u}\mathrm{r}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}}=\frac{{K}_{2}}{\mathrm{ln}\left(\frac{{K}_{1}}{\mathrm{B}\left({\mathrm{T}}_{\mathrm{S}\mathrm{e}\mathrm{n}\mathrm{s}}\right)}+1\right)}-273.15\]（6） 式中：\({T}_{\mathrm{S}\mathrm{u}\mathrm{r}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}}\)是地表真实温度；\({K}_{1}\)和\({K}_{2}\)为Landsat预设的定标常数，可从影像头文件中获取。

为避免多时相数据带来的影响，将多时相地表温度数据逐像元进行归一化处理，具体实现如公式（7）所示： \[Image=\frac{XX}_{\mathrm{m}\mathrm{i}\mathrm{n}}}{{X}_{\mathrm{m}\mathrm{a}\mathrm{x}{X}_{\mathrm{m}\mathrm{i}\mathrm{n}}}\]（7） 式中：\(Image\)代表归一化处理后的影像像元值，像元值处于0到1之间；\(X\)代表影像地表温度值；\({X}_{\mathrm{m}\mathrm{a}\mathrm{x}}\)和\({X}_{\mathrm{m}\mathrm{i}\mathrm{n}}\)分别代表影像地表温度的最大值与最小值。

1.3.2 基于改进的箱线图法提取热异常

改进的箱线图法共有5个统计量，分别为上下四分位数、上下非异常值截距线及中位数。其中，上下非异常值截距线用于区分正常值与异常值，处于非异常值截距线之内的数据为正常值，否则为异常值。该方法与统计学中最常见的箱线图法不同之处主要有两点：一是引入鲍利系数，将鲍利系数以“1”的形式作为样本数据的偏度调控乘数。当数据集是正态分布时，鲍利系数集中在0附近，当数据是右偏态分布时，鲍利系数取值为1；当数据是左偏态时，鲍利系数取值为−1。二是利用半四分位距计算非异常值截距线，可以更好地适用于具有偏态的数据^[18]。该方法的具体实现如下：

（1）获取上下四分位数以及半四分位距

通过IDL程序逐像元获取地表温度数据，然后利用SPSS软件获取相应数据的上下四分位数和中位数，即将数据从小到大排序后，分别处于25%、50%和75%位置的数值。

之后，基于根据四分位数相应的数值，利用EXCEL软件计算上下半四分位距，如公式（8）和公式（9）所示^[19]： \[{SIQR}_{\mathrm{d}\mathrm{o}\mathrm{w}\mathrm{n}}={Q}_{2}Q}_{1}\]（8）\[{SIQR}_{\mathrm{u}\mathrm{p}}={Q}_{3{Q}_{2}\]（9） 式中：\({SIQR}_{\mathrm{d}\mathrm{o}\mathrm{w}\mathrm{n}}\)和\({SIQR}_{\mathrm{u}\mathrm{p}}\)分别为箱体长度的上下界限，即改进的箱线图的上下半四分位距；\({Q}_{1}\)、\({Q}_{2}\)和\({Q}_{3}\)分别为样本数据的上四分位数、中位数和下四分位数。

（2）计算非异常值截距线

非异常值截距线包括非异常最大值截距线和非异常最小值截距线，用来判定样本数据是否为异常值。用非异常最大值截距线作为是否是热异常值的判断依据，若样本数据超出此值，则判定其为热异常值，即可能存在热异常现象。样本数据偏离极端异常值截距线外侧越远，说明其热异常现象就越严重。其实现过程如公式（10）、公式（11）和公式（12）所示： \[{B}_{\mathrm{c}}=\frac{{SIQR}_{\mathrm{u}\mathrm{p}}-{SIQR}_{\mathrm{d}\mathrm{o}\mathrm{w}\mathrm{n}}}{{SIQR}_{\mathrm{u}\mathrm{p}}+{SIQR}_{\mathrm{d}\mathrm{o}\mathrm{w}\mathrm{m}}}\]（10） 式中：Bc 是鲍利系数，介于−1到1之间；\({SIQR}_{\mathrm{u}\mathrm{p}}\)和\({SIQR}_{\mathrm{d}\mathrm{o}\mathrm{w}\mathrm{m}}\)分别为箱体的上下半四分位距。 \[{f}_{\mathrm{u}\mathrm{p}}={Q}_{3}+1.5×IQR×\left(\frac{1+{B}_{\mathrm{c}}}{1B}_{\mathrm{c}}}\right)\]（11）\[{f}_{\mathrm{d}\mathrm{o}\mathrm{w}\mathrm{n}}={Q}_{11.5×IQR×\left(\frac{1-{B}_{\mathrm{c}}}{1+{B}_{\mathrm{c}}}\right)\]（12） 式中：\({f}_{\mathrm{u}\mathrm{p}}\)和\({f}_{\mathrm{d}\mathrm{o}\mathrm{w}\mathrm{n}}\)分别为箱体的非异常最大值截距线和非异常最小值截距线，当样本数据集中某一数据大于\({f}_{\mathrm{u}\mathrm{p}}\)时，则归为热异常值；当样本数据集中某一数据小于\({f}_{\mathrm{d}\mathrm{o}\mathrm{w}\mathrm{n}}\)时，则归为冷异常值。

最后，利用ArcGIS软件逐像元筛选大于非异常最大值截距线的值。

1.3.3 基于热异常频率提取热异常区

热异常区的地表温度往往高于周边地表温度，形成局部热岛现象，且其出现的频率一般大于由其他因素引起的热异常出现的频率。因此，基于热异常提取结果，统计分析年度热异常区出现的频率。为避免漏判现象，当某一热异常区出现的频率大于60%时，则判定为城市热异常区，具体实现如公式（13）所示： \[{f}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}=\frac{{i}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}}{{n}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}}\]（13） 式中：\({f}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}\)代表热异常出现的频率；\({i}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}\)代表局部某一热异常区一年出现的次数；\({n}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}\)代表一年中整体热异常区出现的总次数。

为进一步探究城市热异常区的分布范围，统计分析2008–2017年热异常出现的频率，将提取频率大于60%的区域作为这10 a的城市热异常区，具体实现如公式（14）所示： \[{f}_{\mathrm{g}\mathrm{l}\mathrm{o}\mathrm{b}\mathrm{a}\mathrm{l}}=\frac{{i}_{\mathrm{g}\mathrm{l}\mathrm{o}\mathrm{b}\mathrm{a}\mathrm{l}}}{{n}_{\mathrm{g}\mathrm{l}\mathrm{o}\mathrm{b}\mathrm{a}\mathrm{l}}}\]（14） 式中：\({f}_{\mathrm{g}\mathrm{l}\mathrm{o}\mathrm{b}\mathrm{a}\mathrm{l}}\)代表10 a中热异常出现的频率；\({i}_{\mathrm{g}\mathrm{l}\mathrm{o}\mathrm{b}\mathrm{a}\mathrm{l}}\)代表局部某一热异常区10年中出现的总次数；\({n}_{\mathrm{g}\mathrm{l}\mathrm{o}\mathrm{b}\mathrm{a}\mathrm{l}}\)代表10年中整体热异常区出现的总次数，即取值为10。

使用上述步骤完成2008–2017年三亚热异常空间分布产品数据集，图3是三亚市10 a热异常区空间分布产品示例图。

图3 2008–2017年三亚市热异常区空间分布结果图

本数据集通过以下几方面进行质量控制：

为验证方法有效性，利用Google地球可以查看往年影像的优势进行统计2008–2017年逐年工厂区总个数和改进的箱线图法识别的工厂区个数及工厂区个数，以此间接验证方法的可靠性。具体方法为：首先，在Google地球上逐年圈出工厂区并统计个数；其次，把改进的箱线图法识别的热异常转为kmz格式数据，加载到Google地球里并逐年统计其在工厂区的个数；最后，一个工厂区包含一个或多个工厂，将改进的箱线图法识别出的工厂统计到所属工厂区，统计工厂区个数，该工厂区个数与Google地球里圈出的工厂区个数之比为精度。由于有些工厂停产或排放的热量不足以达到热异常的范围，或某些人为原因致使小领域温度达到热异常范围，以及人工判读工厂的误差，都会影响最终精度，其具体逐年精度如表3。

表3 2008–2017年三亚市热异常区提取精度表

本数据集以三亚市为研究区，基于长时间序列的Landsat数据，结合改进的箱线图法进行热异常提取的成果可为城市规划与建设和生态质量评价等提供数据产品服务，为三亚市的居住环境舒适度提供参考标准。

此外，通过利用遥感技术对城市热异常进行监测和评价，可有助于减轻城市发展对环境的压力，促进环境保护和治理，为三亚国际岛建设提供技术参考。该数据集的长时间序列城市热异常产品和技术也可在其他城市中进行推广和应用。

感谢城市陆表环境团队提供技术支撑服务，感谢USGS和中国空气质量在线监测分析平台提供数据支撑服务。

孟庆岩, 谷艳春, 郝丽春, 等. 三亚市热异常遥感监测数据集[DB/OL]. Science Data Bank, 2018. (2019-03-23). DOI: 10.11922/sciencedb.700.