摘要：针对三亚城市热污染问题，本文以2007–2017年Landsat数据为主要数据源，基于辐射方程传输法反演地表温度和改进的箱线图法提取高温异常区，进行高温异常区频率统计，将频率大于60%的高温异常区判定为城市热污染区，得到11年城市热污染区。本数据集可应用在其他城市，进行多区域城市热污染区提取，从而发挥更高的应用价值。
关键词：城市热污染；辐射方程传输法；改进的箱线图法；高温异常区频率
海南资源环境遥感产品数据集专题
I 区论文（评审中）
•
版本
ZH1
Abstract: Targeting at urban thermal pollution in Sanya, we used radiation transfer equation method to perform surface temperature inversion. Thern based on the Landsat data 2007-2007, modified boxplot method was used together to extract the high temperature anomaly area and summarize the frequency statistics of high temperature anomaly. Areas with an anomaly frequency greater than 60% during the 11 years were identified as thermal pollution areas. This data set is of a high applicability in extraction of urban thermal pollution areas.
Keywords: urban thermal pollution; radiation transfer equation; modified boxplot method; frequency of the high-temperature anomaly area
| 数据库（集）名称 | 三亚市城市热污染遥感监测数据集 |
| 数据作者 | 郝丽春、孟庆岩、胡蝶、张颖、谷艳春、张琳琳 |
| 数据通信作者 | 孟庆岩（mengqy@radi.ac.cn） |
| 数据时间范围 | 2007–2017年 |
| 地理区域 | 三亚市（18°09′34″–18°37′27″N、108°56′30″–109°48′28″E） |
| 空间分辨率 | 100 m |
| 数据量 | 5.24 MB |
| 数据格式 | *.vsdx，*.jpg |
| 数据服务系统网址 | http://www.sciencedb.cn/dataSet/handle/683 |
| 基金项目 | 海南省财政科技计划资助——三亚市城市热污染遥感监测技术（2017CZ07）；基于大数据机器学习与冠层反射率模型结合的水稻叶面积指数提取技术（2018JZ0054）；海南省重大科技计划（ZDKJ2017009）。 |
| 数据库（集）组成 | 本数据集主要包括5个文件夹，每个文件夹的内含文件格式和命名方式如下： 1、封面图文件夹：内含命名为封面图*.jpg格式的数据； 2、三亚市Landsat 8 OLI遥感影像图（真彩色合成）文件夹：内含三亚市Landsat 8 OLI遥感影像图（真彩色合成）*.jpg格式的数据； 3、数据处理流程图文件夹：内含命名为数据处理流程图*.vsdx格式的数据； 4、2007–2017年三亚市年度城市热污染区空间产品示例图文件夹：内含6幅*.jpg格式的数据，分别是2007年三亚热污染区分布图.jpg、2009年三亚热污染区分布图.jpg、2011年三亚热污染区分布图.jpg、2013年三亚热污染区分布图.jpg、2015年三亚热污染区分布图.jpg、2017年三亚热污染区分布图.jpg。 |
| Title | Data set of remote sensing monitoring for urban thermal pollution in Sanya city |
| Data corresponding author | Meng Qingyan (mengqy@radi.ac.cn) |
| Data authors | Hao Lichun, Meng Qingyan, Hu Die, Zhang Ying, Gu Yanchun, Zhang Linlin |
| Time range | 2007–2017 |
| Geographical scope | Sanya (18°09′34″–18°37′27″N, 108°56′30″–109°48′28″E) |
| Spatial resolution | 100 m |
| Data volume | 5.24 MB |
| Data format | *.vsdx, *.jpg |
| Data service system | <http://www.sciencedb.cn/dataSet/handle/683> |
| Sources of funding | Science and Technology Fund of Hainan Provincial Department of Finance - Remote Sensing Monitoring Technology for Urban Thermal Pollution in Sanya (2017CZ07); Technology for Extraction of rice LAI Based on Large Data Machine Learning and Canopy Reflectance Model (2018JZ0054); Major Science and Technology Program of Hainan Province (ZDKJ2017009). |
| Dataset composition | This data set mainly includes five folders comprised of the following: 1. Cover image, which contains a file named “cover image.jpg”; 2. Landsat 8 OLI remote sensing image images (true color synthesis) of Sanya city in JPG format. 3. Data processing flow chart in vsdx format; 4. Spatial product photos of thermal pollution area in Sanya city for 2007-2017, which contain six *. jpg format. They are, respectively, 2007 Sanya thermal pollution layout.jpg, 2009 Sanya thermal pollution layout. jpg, 2011 Sanya thermal pollution layout. jpg, 2013 Sanya thermal pollution layout. jpg, 2015 Sanya thermal pollution layout. jpg, 2017 Sanya thermal pollution layout. jpg. |
热污染是指由自然因素和人类活动中的热排放导致环境温度异常升高的现象[1]。城镇化进程的快速推进可带动经济发展，但同时也对水资源、局部气候、生存环境等产生了一系列负面效应。据统计，近100年来，全球的年平均气温上升了0.7–1℃，而大城市的平均气温则上升了2–3℃，北京、上海、东京等世界性大都市的年均温比郊区高出2℃，夏季高出6℃左右[2]。与物质污染不同的是，物质污染是由污染物引起的，可以通过具体方法和措施得到控制，而热污染则不同，它主要来源于热污染源。只要有能源消费，即使一点污染物也不产生，也会产生热污染[3]。因此，热污染是不可避免的。针对此类现象，国内外学者已经做了大量探索和研究，例如Wu、Chen等通过监测不同地区的核电站热流排放情况评估热污染强度和对周边环境的影响[4-5]。Liu、Xia等基于热异常产品，通过利用不同方式提取工业热源[6-7]；Zhang等利用Landsat数据对首钢搬迁前后的热环境进行探索和分析[8]。还有一部分学者利用城市热岛强度与气温间的关系评估城市热污染[9-10]。可见，热污染已被人们广泛关注，并成为研究热点。但以往关于热污染的大量文献多针对工业城市，或有工业污染的地方，针对非工业城市的热污染研究还较少。而热污染不仅存在于工业城市，还存在于非工业城市。因此，针对非工业城市，如何有效快速提取热污染区，如何缓解改善热污染具有重要研究意义。
三亚市作为我国唯一一个热带滨海旅游城市，其生态环境保护尤显重要。恰逢海南省推进“国际旅游岛”建设，作为海南旅游发展龙头的三亚正以国际化标准致力于打造“一港两地”的国际性滨海旅游城市。但伴随着三亚市城镇化的不断发展，城市热污染对城市环境的影响日益突出。因此，针对三亚城市建设和区域、环境、植被等特点，开展具有针对性的城市热污染研究，对三亚城市建设、规划布局、环境保护均具有重要意义。
综上，针对如何监测三亚城市热污染区，并提取热污染区及验证等问题，本文基于长时间序列的Landsat数据，利用辐射方程传输法获取地表温度信息，结合改进的箱线图法提取相应高温异常区，通过统计分析高温异常区的年度出现频率，得到城市热污染区。该方法不仅能够准确大面积获取城市热污染区，且相较于其他方法，可以有效避免潜在热异常的影响。通过对比分析三亚的城市空气质量指数AQI与相应热污染区面积间的相关性评估产品的可靠性。
1.1
区域范围
1.2
基础数据准备
本文以2007–2017年所有云量小于30%的Landsat 4-5 TM、Landsat 7 ETM+和Landsat 8 OLI数据为主要数据源（获取网址：https://glovis.usgs.gov），2014–2017年空气质量指数AQI为辅助数据源（获取网址：https://www.aqistudy.cn/），具体影像参数如表1所示。
| 数据源 | 发射时间 | 覆盖周期（天） | 波段数 | 分辨率（m） | 星载传感器 | 数据获取途径 |
|---|---|---|---|---|---|---|
| Landsat数据 | Landsat 4-5 TM | 1984年 | 6 | 7 | 波段1-5和波段7的空间分辨率为30 m，波段6（热红外波段）的空间分辨率为120 m。 | USGS官网 |
| Landsat 7 ETM+ | 1999年 | 6 | 8 | 波段1-5和波段7的空间分辨率为30 m，波段8（全色波段）的空间分辨率为15 m，波段6（热红外波段）的空间分辨率为60 m。 | ||
| Landsat 8 OLI | 2013年 | 6 | 11 | 波段1-7和波段9的空间分辨率为30 m，波段8（全色波段）的空间分辨率为15 m，波段10和波段11（热红外波段）的空间分辨率为100 m。 | ||
| 空气质量指数（AQI） | AQI | - | - | - | - | 中国空气质量在线监测分析平台 |
1.3
数据处理
首先对Landsat数据进行拼接、裁剪、重投影等数据预处理。其中Landsat 7 ETM+影像因有条带，需在数据预处理基础上再进行去条带处理。其次用辐射方程传输法反演地表温度和改进的箱线图法提取高温异常区，得到高温异常区空间分布图。最后，以频率大于60%为城市热污染区的标准统计分析热异常出现的频率，得到高温异常区频率空间分布图。同时，统计城市热污染区的总面积，分析该面积与空气质量指数AQI间的相关性来分析数据集的可靠性。具体数据处理流程如图2所示。
1.3.1
地表温度反演
基于Landsat卫星影像，根据Landsat数据热红外波段的增益值与偏置值计算相应表观辐射亮度值；然后，利用影像植被覆盖度估算地表比辐射率；最后，根据Plank反函数和Landsat预设的定标系数计算地表温度。具体实现如下：
（1）表观辐亮度计算
表观辐亮度计算是将影像像元灰度值转换成相应热辐射强度的过程。利用ENVI软件中的辐射定标（Radiometric Calibration）工具将影像的像元灰度值转换为辐射亮度值。
（2）地表比辐射率估算
地表比辐射率表征地面物体的电磁波辐射的一种能力。基于影像的近红外波段和红光波段，根据公式（1）计算归一化植被指数NDVI[11]。
\[\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}=\left({TM}_{NIR}-{TM}_{Red}\right)/\left({TM}_{NIR}+{TM}_{Red}\right)\]
(1)
式中，\({TM}_{NIR}\)为影像的近红外波段，\({TM}_{Red}\)为影像的红光波段。
其次，利用影像植被覆盖度估算地表比辐射率。由于地表物质结构具有差异性，因此针对不同地物类型，其相应的地表比辐射率计算方式也不同。大致可分为3种类型：水体、城镇和自然表面。而水体像元的地表比辐射率很高，与黑体的比辐射率较为相似，因此，在估算水体比辐射率时常被赋值为0.995，自然表面和城镇像元的地表比辐射率估算分别根据公式（2）和公式（3）计算[12-14]：
式中，\({\mathrm{e}}_{Surface}\)代表自然表面像元的地表比辐射率，\({\mathrm{e}}_{Building}\)代表城镇像元的地表比辐射率，\({P}_{Vege}\)代表植被覆盖度，计算方法如公式（4）所示[15-16]：
\[{\mathrm{P}}_{Vege}={\left[\frac{NDVI-{NDVI}_{Soil}}{{NDVI}_{Vege}-{NDVI}_{Soil}}\right]}^{2}\]
（4）
式中，\(NDVI\)是影像的归一化植被指数，\({NDVI}_{Vege}\)是影像中纯植被像元的NDVI值，\({NDVI}_{Soil}\)是影像中纯土壤像元的NDVI值。当某个像元的\(NDVI\)大于0.70时，\({P}_{Vege}\)取值为1；当\(NDVI\)小于0.05时，\({P}_{Vege}\)取值为0；当\(NDVI\)介于0.05–0.7之间时，则分别将\({NDVI}_{Vege}\)、\({NDVI}_{Soil}\)取值0.70和0.05，并利用上式估算影像的植被覆盖度。
（3）地表温度计算
基于热红外波段的表观辐亮度和地表比辐射率计算结果，根据辐射传输方程反推地表热辐射亮度，如公式（5）所示。利用Plank反函数结合Landsat预设的定标常数K1、K2计算地表真实温度[17-18]，如公式（6）所示：
\[B\left({T}_{Sens}\right)=\frac{L-{L}^{↑}-\tau \left(1-\epsilon \right){L}_{↓}}{\tau \epsilon }\]
（5）
式中，\(B\left({T}_{Sens}\right)\)是地表热辐射亮度，L是行星亮度温度值，\({L}^{↑}\)、\({L}_{↓}\)分别是大气上行辐射和大气下行辐射，可根据Landsat卫星过境时间、兴趣区中心经纬度、大气模式及传感器类型等获得，具体参数可通过NASA相关网页直接获得（http://atmcorr.gsfc.nasa.gov/）。\(\tau \)为大气路径透过率，可由NASA相关网页查询获取；\(\epsilon \)为地表比辐射率。
\[{T}_{Surface}=\frac{{K}_{2}}{\mathrm{ln}\left(\frac{{K}_{1}}{\mathrm{B}\left({\mathrm{T}}_{Sens}\right)}+1\right)}-273.15\]
（6）
式中，\({T}_{Surface}\)是地表真实温度，\({\mathrm{K}}_{1}\)、\({\mathrm{K}}_{2}\)为Landsat预设的定标常数，可从影像头文件中获取。
1.3.2
基于改进的箱线图法提取热异常
为避免多时相数据带来的影响，将多时相地表温度数据逐像元进行归一化处理，具体实现如公式（7）所示：
\[\mathrm{I}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}=\frac{\mathrm{X}-{\mathrm{X}}_{min}}{{X}_{max}-{X}_{min}}\]
（7）
式中，\(\mathrm{I}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}\)代表归一化处理后的影像，像元值处于0到1之间。\(\mathrm{X}\)代表影像地表温度值，\({X}_{max}\)、\({X}_{min}\)分别代表影像地表温度的最大值与最小值。
改进的箱线图法共有5个统计量，分别为上下四分位数，上下非异常值截距线及中位数。其中，上下非异常值截距线用于区分正常值与异常值，处于非异常值截距线之内的数据为正常值，否则为异常值。该方法与统计学中最常见的箱线图法不同之处主要有两点：一是引入鲍利系数，将鲍利系数以“1”的形式作为样本数据的偏度调控乘数。当数据集是正态分布时，鲍利系数集中在0附近，当数据是右偏态分布时，鲍利系数取值为1；当数据是左偏态时，鲍利系数取值为−1。二是利用半四分位距计算非异常值截距线，可以更好地适用于具有偏态的数据[19]。该方法的具体实现如下：
（1）获取上下四分位数以及半四分位距
通过IDL程序逐像元获取地表温度数据，然后利用SPSS软件获取相应数据的上下四分位数和中位数，即将数据从小到大排序后，分别处于25%、50%和75%位置的数值。
之后，基于根据四分位数相应的数值，利用EXCEL软件计算上下半四分位距，如公式（8）和公式（9）所示[20]：
式中，\({SIQR}_{down}\)、\({SIQR}_{up}\)分别为箱体长度的上下界限，即改进的箱线图的上下半四分位距，\({Q}_{1}\)、\({Q}_{2}\)、和\({Q}_{3}\)分别为样本数据的上四分位数、中位数和下四分位数。
（2）计算非异常值截距线
非异常值截距线包括非异常最大值截距线和非异常最小值截距线，用来判定样本数据是否为异常值。用非异常最大值截距线作为是否是热异常值的判断依据，若样本数据超出此值，则判定其为热异常值，即可能存在城市热污染现象。样本数据偏离极端异常值截距线外侧越远，说明其热异常现象就越严重。其实现过程如公式（10）、公式（11）和公式（12）所示：
式中，Bc 是鲍利系数，介于−1到1之间，\({SIQR}_{up}\)、\({SIQR}_{dowm}\)分别为箱体的上下半四分位距。
\[{f}_{\mathrm{u}\mathrm{p}}={Q}_{3}+1.5×\mathrm{I}\mathrm{Q}\mathrm{R}×\left(\frac{1+{B}_{\mathrm{c}}}{1-{B}_{\mathrm{c}}}\right)\]
（11）
\[{f}_{\mathrm{d}\mathrm{o}\mathrm{w}\mathrm{n}}={Q}_{1}-1.5×\mathrm{I}\mathrm{Q}\mathrm{R}×\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}_{local}\)代表热异常出现的频率，\(i\)代表局部某一热异常区一年出现的次数，\(n\)代表一年中整体热异常区出现的总次数。
为进一步探究城市热污染区的分布范围，统计分析2007–2017年热异常出现的频率，将提取频率大于60%的区域作为这11年的城市热污染区，具体实现如公式（14）所示：
式中，\({f}_{global}\)代表11年中热异常出现的频率，\(i\)代表局部某一热异常区11年中出现的总次数，\(n\)代表11年中整体热异常区出现的总次数，即取值为11。
本数据成果有利于更合理地对三亚城市功能模块进行配置，减轻城市发展对环境的压力，为进一步提高三亚市人居环境的舒适度提供参考。同时为园林绿地规划和质量评价、市政管理、城市规划等提供数据产品服务。
此外，通过利用遥感技术对城市热污染进行监测和评价，可有助于合理规划城市建设、促进环境保护和治理，为三亚国际岛建设提供技术参考。长时间序列城市热污染产品可在城市建设和环境保护部门推广应用。
[4]
Wu C , Wang Q , Yang Z , et al. Monitoring heated water pollution of the DaYaWan nuclear power plant using TM images[J]. International Journal of Remote Sensing, 2007, 28(5):885-890.
[5]
ChuqunChen, PingShi, QingwenMao. Application of Remote Sensing Techniques for Monitoring the Thermal Pollution of Cooling-Water Discharge from Nuclear Power Plant[J]. Environmental Letters, 2003, 38(8):10.
[6]
Liu Y , Hu C , Zhan W , et al. Identifying industrial heat sources using time-series of the VIIRS Nightfire product with an object-oriented approach[J]. Remote Sensing of Environment, 2017:S0034425717304820.
[7]
Xia H, Chen Y, Quan J. A simple meathod based on the thermal anomaly index to detect industrial heat sources[J]. International Journal of Applied Earth Observation & Geoinformation 2018 (73): 627-637.
[8]
Zhang L , Meng Q , Sun Z , et al. Spatial and Temporal Analysis of the Mitigating Effects of Industrial Relocation on the Surface Urban Heat Island over China[J]. ISPRS International Journal of Geo-Information, 2017, 6(121):121.
[9]
Murata A , Sasaki H , Hanafusa M , et al. Estimation of urban heat island intensity using biases in surface air temperature simulated by a nonhydrostatic regional climate model[J]. Theoretical & Applied Climatology, 2013, 112(1-2):351-361.
[10]
Yang X , Leung L R , Zhao N , et al. Contribution of urbanization to the increase of extreme heat events in an urban agglomeration in east China[J]. Geophysical Research Letters, 2017, 44.
[15]
Sobrino J A , Jimenez-Muoz J C , Soria G , et al. Land Surface Emissivity Retrieval From Different VNIR and TIR Sensors[J]. IEEE Transactions on Geoscience and Remote Sensing, 2008, 46(2):316-327.
[16]
Carlson T N , Ripley D A . On the relation between NDVI, fractional vegetation cover, and leaf area index[J]. Remote Sensing of Environment, 1997, 62(3):241-252.
[17]
Jimenezmunoz J C , Cristobal J , Sobrino, José A, et al. Revision of the Single-Channel Algorithm for Land Surface Temperature Retrieval From Landsat Thermal-Infrared Data[J]. IEEE Transactions on Geoscience & Remote Sensing, 2008, 47(1):339-349.
[18]
Chander G , Markham B . Revised landsat-5 tm radiometric calibration procedures and postcalibration dynamic ranges[J]. IEEE Transactions on Geoscience and Remote Sensing, 2003, 41(11):2674-2677.
[19]
Walker M L , Dovoedo Y H , Chakraborti S , et al. An Improved Boxplot for Univariate Data[J]. American Statistician, 2018:1-13.
郝丽春, 孟庆岩, 胡蝶, 等. 三亚市城市热污染遥感监测数据集[DB/OL]. Science Data Bank, 2018. (2018-12-05). DOI: 10.11922/sciencedb.700.
论文引用格式
郝丽春, 孟庆岩, 胡蝶, 等. 三亚市城市热污染遥感监测数据集[J/OL]. 中国科学数据, 2018. (2018-12-07). DOI: 10.11922/csdata.2018.0077.zh.
出版历史
I区发布时间：2018年12月12日
（
版本ZH1）
