, 谢美慧 , 田琨 , 谢晓晓 清华大学 土木工程系, 地球空间信息研究所, 北京 100084
收稿日期：2016-03-31
作者简介：刘钊(1967—), 男(汉), 副教授。E-mail:liuz@tsinghua.edu.cn
摘要：PM2.5空气污染问题目前是社会关注热点以及学术研究重点。该文对PM2.5污染的自然城市提取进行了研究，结合PM2.5的站点监测数据和气溶胶遥感数据并采用协同Kriging插值实现了PM2.5数据空间化，然后采用首尾分割分类方法实现了PM2.5污染分布的分类和污染自然城市的提取。对中国大陆PM2.5自然城市的提取结果进行了分析和讨论。结果表明：采用适当的分割阈值，首尾分割分类方法可以有效进行PM2.5污染自然城市提取工作，有助于决策者合理划分PM2.5联合治理的区域范围。
关键词：PM2.5空气污染协同Kriging插值首尾分割自然城市
Classification of PM2.5 for natural cities based on co-Kriging and head/tail break algorithms
LIU Zhao
, XIE Meihui, TIAN Kun, XIE XiaoxiaoInstitute of Geomatics, Department of Civil Engineering, Tsinghua University, Beijing 100084, China
Abstract: PM2.5 air pollution is now a hot topic in both social and academic circles. This study investigated the classification of natural cities based on PM2.5 concentrations in Mainland China. Firstly, the PM2.5 data obtained at monitoring stations and aerosol optical depths (AOD) obtained by remote sensing were fused to yield more accurate PM2.5 spatial distributions using a co-Kriging algorithm. Then, the PM2.5 concentrations were classified using the head/tail break clustering algorithm to identify natural cities with high PM2.5 pollution levels. Distribution of natural cities was also analyzed. The results show that the head/tail break algorithm with an appropriate segmentation threshold can efficiently identify natural cities with high PM2.5 concentrations. These classification results can guide policy makers to divide the country into several areas for pollution control.
Key words: PM2.5 air pollutionco-Kriging interpolationhead/tail breaknatural city
随着社会经济的发展与进步，人们的生活方式发生了巨大的变化，但在发展的同时也引发了诸多问题，其中环境污染问题已经成为社会关注热点，并越来越深刻地影响社会发展和人们的生活质量[1]。在中国大陆，PM2.5空气污染问题因其严重影响人们身体健康而引发社会广泛关注和高度重视。为实时掌握全国(本文特指中国大陆，不含香港和澳门特别行政区及台湾地区，下同)PM2.5污染情况，中国环境保护部先后在全国各主要城市建立了1 000多个空气监测站点，并在公开化平台(http://www.zhb.gov.cn/)实时发布PM2.5监测数据，众多****也通过该平台获取PM2.5监测数据用于环境污染相关研究。
目前，国内外****针对PM2.5空气污染时空分布问题进行了广泛探讨，但关于PM2.5空气污染空间分布特性的研究多采用传统统计方法，缺乏地理空间统计的分析[2]。本文采用协同Kriging空间插值方法对全国2015年PM2.5监测数据进行空间化分析，使用首尾分割方法提取全国PM2.5污染自然城市分布，最后针对污染自然城市提取结果，从地理空间统计层面分析和讨论全国PM2.5时空污染分布状况。
1 研究方法比较分析为精确研究PM2.5空气污染水平，研究者需要具有高空间分辨率的PM2.5分布数据，但PM2.5监测数据为稀疏的点数据，无法直接提供高空间分辨率的连续面数据。为了解决这个问题，多种点数据空间化的方法被提出，如大气扩散模型、土地使用回归模型和插值方法等[3]。其中，插值方法又可分为确定性内插方法(如反距离插值、全局多项式、径向基函数等)和地质统计学内插方法(如普通Kriging、协同Kriging、简单Kriging等)[4-5]。Kriging插值方法也被广泛应用于PM2.5分布数据空间化研究中，如Sampson等使用普通Kriging模型实现了美国地区PM2.5分布数据空间化[6]。
与PM2.5监测站数据相比较，MODIS气溶胶光学厚度(aerosol optical depth, AOD)数据具有广泛、连续的空间分布和长期监测等特性，且MODIS AOD数据与PM2.5数据具有一定的相关性，逐渐成为PM2.5空间化分析过程中的一个重要变量[7]。XIE等基于线性混合效应模型建立MODIS AOD与PM2.5的混合线性相关关系，实现了对北京在2013年3月至2014年4月间PM2.5数据的空间化，模型模拟精度达到0.8以上[8]。LONG等也基于MODIS AOD和PM2.5监测数据建立协同Kriging插值模型，得到了全国PM2.5污染空间化结果，空间化效果良好[9]。
在城市研究中，城市区域的提取至关重要，常用的城市提取方法有监督分类和非监督分类[10]。与传统城市边界相对应，JIANG在2010年提出了自然城市(natural cities)的概念[11]。在特定研究对象中，通过对均值进行切割，大于均值部分组成的集合即称为自然城市，其边界范围和边界形状与人为设定的边界(行政边界)具有很大程度上的不同，自然城市可以客观反映环境污染在空间上的分布情况。
自然界中表现出重尾分布特性的数据服从幂律法则。针对这种数据类型，JIANG提出了首尾分割的概念[12]，并使用ht-index来描述重尾分布特性，其大小等于1加上首尾分割次数[13-14]。LIU在2014年基于首尾分割方法对夜间灯光数据进行了自然城市提取，发现多数城市的自然城市边界与行政边界具有一定的重合度[10]。
本文基于自然城市的概念，提出使用污染自然城市来刻画PM2.5空气污染空间分布状况，收集并处理全国2015年PM2.5监测数据，获得2015年全国PM2.5空间分布，使用首尾分割分类法研究全国PM2.5污染自然城市的提取。
2 数据源及数据处理2.1 PM2.5数据中国环境保护部为公众实时提供了全国各监测站点的污染物监测状况，PM25.in网站(http://pm2.5.in/)则将环境保护部提供的数据进行统一整理并提供查询接口。本文通过PM25.in收集了全国2015年1月1日至12月31日369个城市的PM2.5监测站点每h浓度监测数据，共有1 497个监测站点，具体分布如图 1a所示。对所有PM2.5监测数据进行平均化处理得到每日PM2.5监测数据均值，将每日PM2.5监测数据均值统一转为ESRI Shapefile(shp)格式。
 |
| 图 1 全国PM2.5监测站分布及MYD04_3K分布示例 |
| 图选项 |
2.2 3 km MODIS AOD数据MODIS为搭载在Terra和Aqua太阳同步极轨卫星平台上的传感器，分别于1999年和2002年发射成功。MODIS标准数据共分为5级44种，包括陆地标准、大气标准和海洋标准数据[15]。标准MODIS AOD数据空间分辨率为10 km，但为了满足更高空间分辨率的需求，近期公布了空间分辨率为3 km的MODIS数据集006(MOD04_3K和MYD04_3K)[8, 16]。本文应用了全国2015年1月1日至12月31日的MYD04_3K数据。首先，使用ENVI二次开发工具并引入MODIS转换工具包(MODIS conversion toolkit, MCTK)进行几何校正，接着再经过拼接、提取等步骤获得全国范围内的气溶胶光学厚度分布。由于云、高反射率(例如，冰雪覆盖地区)和反演误差的原因，MODIS AOD数据有大量丢失数据的情况。如图 1b所示，MODIS AOD数据范围为-0.05~5。
2.3 MODIS AOD数据验证目前，有两种手段可以检测气溶胶光学厚度：地面监测、遥感卫星监测[17]。为验证MODIS AOD数据的有效性，多采用地面AOD监测数据进行验证。AERONET(aerosol robotic network)是由美国航空航天局NASA和法国的PHOTONS(photométrie pour le traitement opérationnel de normalisation satellitaire)联合在全球设立的气溶胶光学特性监测基地网络，在2002年全球就有155个永久性监测站点[15]。本次选择了level1.5的AERONET数据，且在全国范围内2015年1月1日至12月31日间的监测站点有9个。
AERONET监测站点提供多个波段的气溶胶光学厚度数据，但不包括0.55 μm波段处的气溶胶光学厚度数据。为与MODIS AOD数据保持一致，需根据Angstr?m公式[17]估计AERONET在0.55 μm处的气溶胶光学厚度值，
| ${\tau _\lambda } = \beta {\lambda ^\alpha }.$ |
最终得到304个有效的“AERONET-MODIS AOD”数据组。相关性分析表明，MODIS AOD数据与AERONET数据在0.55 μm波段处的数值具有较强的相关性，可认为MODIS AOD数据具有较高的准确性，如图 2所示，R2=0.83。
 |
| 图 2 AERONET与MODIS AOD气溶胶光学厚度数据的相关性分析 |
| 图选项 |
2.4 PM2.5估计MODIS AOD数据与PM2.5监测数据在2015年日均值的变化趋势如图 3所示，PM2.5值在春、冬季节较高，夏、秋季节较低，而MODIS AOD值则普遍偏高，但细化到月份中，则可以发现PM2.5与MODIS AOD数值变化趋势相近。
 |
| 图 3 2015年全国PM2.5与AOD均值变化 |
| 图选项 |
结合MODIS AOD数据和PM2.5监测站数据，采用协同Kriging插值方法进行PM2.5估计。以PM2.5监测数据为主变量，MODIS AOD数据为辅助变量，对应的协同Kriging方法的估计公式为
| $z_2^*\left( {{x_0}} \right) = \sum\limits_{i = 1}^{{N_1}} {{\lambda _{1i}}} {Z_1}\left( {{x_{1i}}} \right) + \sum\limits_{j = 1}^{{N_2}} {{\lambda _{2j}}} {Z_2}\left( {{x_{2j}}} \right).$ |
| $\sum\limits_{i = 1}^{{N_1}} {{\lambda _{1i}}} = 0, $ |
| $\sum\limits_{j = 1}^{{N_2}} {{\lambda _{2j}}} = 1.$ |
2.5 污染自然城市提取通过统计2015年全国PM2.5均值的空间分布情况，计算PM2.5均值在空间分布上所有像元平均值以及超过像元平均值的所有像元个数占总共像元个数的百分比。在本文中，以50%作为首尾分割阈值进行分割计算，若平均值以上的像元个数比例不超过50%，将继续对超过均值的所有像元进行循环计算，当均值以上的像元个数比例超过50%，循环终止。第n轮循环的计算公式如下：
| ${\text{mea}}{{\text{n}}_n} = \frac{{\sum\limits_{i = 1}^{{m_n}} {{C_{in}}} }}{{{m_n}}}, $ |
| $\begin{gathered} hea{d_n} = \frac{{\sum\limits_{i = 1}^{{m_n}} {({C_{in}} > {\text{mea}}{{\text{n}}_n})} }}{{{m_n}}}. \hfill \\ \end{gathered} $ |
假设计算至第N轮循环结束，则将meanN-1作为自然城市提取阈值，所有PM2.5数值大于meanN-1的像元集合即为将要提取的自然城市。
3 结果及讨论3.1 PM2.5污染分布描述分析图 4 2015年全国PM2.5均值分布以PM2.5监测数据为主变量，以MODIS AOD数据为协变量，采用协同Kriging估计方法得到空间分辨率为3 km的2015年全国PM2.5日污染分布，并对所有PM2.5日污染分布图作平均值计算得到2015年PM2.5分布图。图 4所示为2015年全国PM2.5均值分布图。其中，PM2.5均值为45.92 μg/m3，全国大约有40%的区域PM2.5均值超过该值。从图 4中可以看出，环渤海区、西北地区和东北地区均为PM2.5污染严重区域，其中西北地区PM2.5污染现象严重多为沙尘暴天气引起，而环渤海地区和东北地区的PM2.5污染现象则是由于重工业、汽车尾气等因素所造成[18]。
 |
| 图 4 2015年全国PM2.5均值分布 |
| 图选项 |
3.2 污染自然城市提取结果分析通过分析2015年全国PM2.5均值分布图(图 4)可以发现，低PM2.5均值范围远大于高PM2.5均值范围，其统计特征符合重尾分布，适合使用首尾分割分类方法对PM2.5质量浓度分布进行分类。表 1为2015年全国PM2.5均值分布图的首尾分割统计结果，计算中以50%作为首尾分割阈值，得到2个满足首尾分割条件的均值，分别为45.91 μg/m3和61.31 μg/m3。基于首尾分割分类结果可将2015年全国PM2.5均值分布划分为3层，其首尾分割指数ht-index=3。
表 1 2015年全国PM2.5均值分布图的首尾分割统计结果
| 序号 | 数值总和 | 数量 | 均值 | 首部像元数量 | 首部数量占比 | 尾部像元数量 | 尾部数量占比 |
| 1 | 77 088 104 | 1 678 858 | 45.92 | 669 703 | 0.40 | 1 009 155 | 0.60 |
| 2 | 41 061 720 | 669 703 | 61.31 | 236 540 | 0.35 | 433 163 | 0.65 |
| 3 | 18 253 270 | 236 540 | 77.17 | 107 857 | 0.46 | 128 683 | 0.54 |
| 4 | 9 550 114 | 107 857 | 88.54 | 61 478 | 0.57 | 46 379 | 0.43 |
表选项
为使得提取的PM2.5自然城市既能体现更丰富的城市信息，又能体现污染严重区域，需要选取合适的阈值进行自然城市提取。本文分别使用45.91 μg/m3和61.31 μg/m3进行PM2.5自然城市提取，并对二者提取结果进行对比分析，最终选择使用61.31 μg/m3作为PM2.5污染自然城市提取阈值。
PM2.5环境污染是一种区域化现象，呈现出明显的“区块”效应。图 5为在61.31 μg/m3阈值下所提取的两个PM2.5自然城市示例。其中：图 5a所示的PM2.5自然城市范围覆盖或部分覆盖新疆乌鲁木齐和吐鲁番两个行政城市，图 5b所示的PM2.5自然城市则覆盖或部分覆盖北京、河北、河南、山东和湖北行政区域。
 |
| 图 5 PM2.5自然城市提取结果示例 |
| 图选项 |
4 总结本文基于PM2.5地面监测数据和MODIS AOD数据，利用协同Kriging方法实现了PM2.5污染分布数据的空间化，使用首尾分割方法实现了全国PM2.5自然城市的提取工作。2015年全国PM2.5污染均值为45.91 μg/m3，从空间分布上看，大气污染在全国分布不均匀，在全国大部分地区污染现象较为轻微，而污染严重区主要集中在华北、东北和西北地区。
通过对PM2.5污染分布进行统计分析，发现2015年全国PM2.5污染分布符合“低PM2.5数值范围远大于高PM2.5数值范围”的分布规律。使用首尾分割方法进行PM2.5自然城市提取，提取结果表明，PM2.5环境污染是一种区域化现象。由于PM2.5污染为空气污染问题，其污染范围受到多种大气环境因素作用，因此PM2.5自然城市边界范围与行政边界范围具有很大的不同，且自然城市边界会横跨多个行政城市边界。大气污染是一种跨多个城市的区域性问题，其治理也不仅仅是某个具体城市的任务，而是应该由多个城市共同治理。
参考文献
| [1] | Journal of Central South University(Science and Technology), 41(2):649-654.-->郑思齐, 张晓楠, 宋志达, 等. 空气污染对城市居民户外活动的影响机制:利用点评网外出就餐数据的实证研究[J]. 清华大学学报(自然科学版), 2016, 56(1): 89–96.ZHENG Siqi, ZHANG Xiaonan, SONG Zhida, et al. Influence of air pollution on urban residents' outdoor activity: Empirical study based on dining-out data from the Dianping Website[J]. J Tsinghua Univ (Sci and Tech), 2016, 56(1): 89–96. (in Chinese) |
| [2] | Journal of Central South University(Science and Technology), 41(2):649-654.-->施益强, 王坚, 张枝萍. 厦门市空气污染的空间分布及其与影响因素空间相关性分析[J]. 环境工程学报, 2014(12): 5406–5412.SHI Yiqiang, WANG Jian, ZHANG Zhiping. Analysis on spatial distribution of air pollution and its spatial correlation with influencing factors in Xiamen City[J]. Chinese Journal of Environmental Engineering, 2014(12): 5406–5412. (in Chinese) |
| [3] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Hoek G, Beelen R, De Hoogh K, et al. A review of land-use regression models to assess spatial variation of outdoor air pollution[J]. Atmospheric Environment, 2008, 42(33): 7561–7578. DOI:10.1016/j.atmosenv.2008.05.057 |
| [4] | Journal of Central South University(Science and Technology), 41(2):649-654.-->易湘生, 李国胜, 尹衍雨, 等. 土壤厚度的空间插值方法比较:以青海三江源地区为例[J]. 地理研究, 2012, 31(10): 1793–1805.YI Xiangsheng, LI Guosheng, YIN Yanyu, et al. Comparison on soil depth prediction among different spatial interpolation methods: A case study in the Three-River Headwaters Regions of Qinghai Province[J]. Geographic Research, 2012, 31(10): 1793–1805. (in Chinese) |
| [5] | Journal of Central South University(Science and Technology), 41(2):649-654.-->张成才, 秦昆, 卢艳, 等. GIS空间分析理论与方法[M]. 武汉: 武汉大学出版社, 2004.ZHANG Chengcai, QIN Kun, LU Yan, et al. Theories and Methods of Spatial Analysis in GIS[M]. Wuhan: Wuhan University Press, 2004. (in Chinese) |
| [6] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Sampson P D, Richards M, Szpiro A A, et al. A regionalized national universal kriging model using partial least squares regression for estimating annual PM 2.5 concentrations in epidemiology[J]. Atmospheric Environment, 2013, 75: 383–392. DOI:10.1016/j.atmosenv.2013.04.015 |
| [7] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Van D A, Martin R V, Park R J. Estimating ground-level PM2.5 using aerosol optical depth determined from satellite remote sensing[J]. Journal of Geophysical Research: Atmospheres, 2006, 111: D21201. DOI:10.1029/2005JD006996 |
| [8] | Journal of Central South University(Science and Technology), 41(2):649-654.-->XIE Yuanyu, WANG Yuxuan, ZHANG Kai, et al. Daily estimation of ground-level PM2.5 concentrations over Beijing using 3 km resolution MODIS AOD[J]. Environmental Science & Technology, 2015, 49(20): 12280–12288. |
| [9] | Journal of Central South University(Science and Technology), 41(2):649-654.--> LONG Ying, WANG Jianghao, WU Kang, et al. Population Exposure to Ambient PM2.5 at the Subdistrict Level in China [Z/OL]. (2014-08-25) [2016-02-15]. https://ssrn.com/ab-stract=2486602 |
| [10] | Journal of Central South University(Science and Technology), 41(2):649-654.--> LIU Qingling. A Case Study on the Extraction of the Natural Cities from Nightlight Image of the United States of America [D]. G?vle, Sweden: University of G?vle, 2013. |
| [11] | Journal of Central South University(Science and Technology), 41(2):649-654.--> JIA Tao, JIANG Bing. Measuring Urban Sprawl Based on Massive Street Nodes and the Novel Concept of Natural Cities [Z/OL]. (2010-12-08) [2016-02-15]. https://arxiv.org/abs/1010.0541 |
| [12] | Journal of Central South University(Science and Technology), 41(2):649-654.-->JIANG Bing. Head/tail breaks: A new classification scheme for data with a heavy-tailed distribution[J]. The Professional Geographer, 2013, 65(3): 482–494. DOI:10.1080/00330124.2012.700499 |
| [13] | Journal of Central South University(Science and Technology), 41(2):649-654.-->JIANG Bing, YIN Junjun. Ht-index for quantifying the fractal or scaling structure of geographic features[J]. Annals of the Association of American Geographers, 2014, 104(3): 530–540. DOI:10.1080/00045608.2013.834239 |
| [14] | Journal of Central South University(Science and Technology), 41(2):649-654.-->GAO Peichao, LIU Zhao, XIE Meihui, et al. CRG index: A more sensitive ht-index for enabling dynamic views of geographic features[J]. The Professional Geographer, 2016, 68(4): 533–545. DOI:10.1080/00330124.2015.1099448 |
| [15] | Journal of Central South University(Science and Technology), 41(2):649-654.-->张小娟. 基于MODIS遥感DT和DB数据集的中国AOD分布与变化[D]. 南京: 南京信息工程大学, 2014. ZHANG Xiaojuan. The Distribution and Variation of AOD over China Based on MODIS Remote Sensing DT and DB Data Set [D]. Nanjing: Nanjing University of Information Science and Technology, 2014. (in Chinese) |
| [16] | Journal of Central South University(Science and Technology), 41(2):649-654.-->孙晓雷, 甘伟, 林燕, 等. MODIS 3 km气溶胶光学厚度产品检验及其环境空气质量指示[J]. 环境科学学报, 2015, 35(6): 1657–1666.SUN Xiaolei, GAN Wei, LIN Yan, et al. Validation of MODIS 3 km aerosol optical depth product and its air quality indication[J]. Acta Science Circumstantiae, 2015, 35(6): 1657–1666. (in Chinese) |
| [17] | Journal of Central South University(Science and Technology), 41(2):649-654.-->ZHOU Chunyan, LIU Qinhuo, TANG Yong, et al. Comparison between MODIS aerosol product C004 and C005 and evaluation of their applicability in the north of China[J]. Journal of Remote Sensing, 2009(5): 854–872. |
| [18] | Journal of Central South University(Science and Technology), 41(2):649-654.-->YOU Wei, ZANG Zengliang, ZHANG Lifeng, et al. National-scale estimates of ground-level PM2.5 concentration in China using geographically weighted regression based on 3 km resolution MODIS AOD[J]. Remote Sensing, 2016, 8(3): 184. DOI:10.3390/rs8030184 |
