2000年—2017年中国PM2.5污染的健康影响及经济损失时空分布模式

周梦杰，胡文庆，杨思远，朱 翔*，赵琛琛，白紫月

(1.湖南师范大学资源与环境科学学院，长沙 410006；2.地理空间大数据挖掘与应用湖南省重点实验室，长沙 410006)

随着工业化和城市化的发展，空气污染日趋严重，而PM2.5作为大气颗粒物的首要组成成分，近年来广泛受到公众的关注.国内外已有的流行病学研究表明，大气中PM2.5及其成分可引起人体炎症、心率变异、凝血障碍等急性效应[1-3]，与循环、呼吸系统等疾病发病率、死亡率的增加成正相关.美国ASC研究显示PM2.5浓度升高10 μg·m-3，会使全死因死亡、心肺疾病死亡、肺癌死亡风险分别增加4%、6%、8%[4]；Raaschou-Nielsen等[5]一项涉及9个欧洲国家、31万人，平均随访12.8年的研究显示，PM2.5每增加5 μg·m-3，肺癌发病风险增加18%；Bowe等[6]发现PM2.5浓度高于2.4 μg·m-3时，糖尿病风险随PM2.5浓度增加显著上升，PM2.5达到10 μg·m-3后，糖尿病风险随PM2.5浓度增加平稳上升；李娜等[7]发现PM2.5、PM10、NO2三种污染物浓度每上升1 μg·m-3，研究人群中高血压患病风险会分别增加2.9%、1.5%、6.9%.PM2.5引起人体健康水平下降，对人群健康及社会经济造成的损失不容忽视.Cohen等[8]研究发现2015年全球范围内PM2.5的暴露导致约420万例死亡，占全球总死亡人数的7.6%；李慧娟等[9]估算2015年我国62个环保重点检测城市PM2.5污染约造成12.51万人早逝及1 009.59万人次患病、门诊和住院，约占这些城市市区总人口的3.53%；黄德生等[10]估计2009年京津冀地区PM2.5浓度若达到2012年的《空气质量标准》，将产生1 729亿元的健康效益，相当于2009年GDP的4.68%.测算大气PM2.5对健康的影响并进行货币化分析，可以为制定大气污染防控政策和相关措施的效益—成本分析提供数据依据，因而具有重要意义.

目前已有的大气污染健康经济效应研究中，对包括PM2.5在内的细颗粒物浓度数据的选取多采用监测站数据.如谢鹏等[11]基于珠三角地区16个站点观测的大气污染浓度，评价了2006年珠三角地区大气污染的人群健康影响；谢元博等[12]利用北京市23个大气环境监测站点的PM2.5浓度数据，评估了高浓度PM2.5污染的健康影响.由于我国监测站点分布不均衡，大部分位于东部和中部地区，城市测站也比农村更密集，因此以少量站点数据估计总体PM2.5浓度分布存在不合理性.此外，现有的研究尺度多为行政区划单元.如李慧娟等[9]对62个环保重点监测城市开展了PM2.5污染的人群健康影响和经济损失研究；陈仁杰等[13]对我国113个主要城市进行了大气颗粒物污染的健康经济损失评价.仅以城市为研究尺度，忽略了人口经济空间分布的异质性和复杂性.此外，现有研究多是对单一年份的静态估算，如曾贤刚[14]等估算了2017年我国PM2.5污染的健康效应和经济损失；王桂芝等[15]基于CGE模型估算了2013年北京市PM2.5污染的健康影响.仅对单一或少量年份的健康影响和经济损失进行评估，会使结果缺乏动态变化分析.近年来，遥感卫星的蓬勃发展为PM2.5大范围、长时间的连续观测提供了可能，大气PM2.5遥感定量反演技术顾及了卫星观测气溶胶光学厚度(AOD)与PM2.5关系的时空异质性[16].本研究在5 km×5 km空间格网尺度下，采用遥感反演PM2.5暴露浓度数据，通过相关健康效应模型和经济损失模型，对2000年—2017年全国PM2.5污染的人群健康影响及经济损失进行了估算，获得了PM2.5污染的健康影响及经济损失的时空动态分析结果.

1 数据处理与分析方法

1.1 PM2.5暴露浓度和暴露人口估计

在以往的研究中，PM2.5暴露浓度数据的选择多采用地面监测站数据，然而监测站建设成本高、实施难度较大、数量较少，难以实现大范围连续监测，遥感卫星能提供大范围观测数据.鉴于此本文采用大气PM2.5遥感反演浓度数据[17-18].该遥感反演数据集是基于大气化学传输模型(GEOS-Chem)和地理加权回归模型(GWR)，估算的中国地区2000年—2017年0.01°×0.01°年平均PM2.5浓度分布.以2015年为例，PM2.5浓度空间分布如图1.本文在ArcMap中对数据进行转换后，得到5 km×5 km空间尺度的PM2.5浓度分布，作为计算所需要的基准浓度面.

在PM2.5污染引起的健康影响的估算中，需要拟定作为参考的浓度，即控制浓度，认为污染物浓度超过控制浓度时，引起人群的健康效应.一般来说，控制浓度可以选择零值、流行病学研究观察到的最低作用浓度值或最高无作用浓度值、政府机构制定的卫生标准等[19].本文选用我国环境空气质量标准(GB 3095-2012)中规定的PM2.5二级浓度限值35 μg·m-3作为控制浓度.

图1 2015年全国PM2.5浓度空间分布Fig.1 Spatial distribution of national PM2.5 concentration in 2015

本研究所采用的暴露人口数据来源于中国科学院资源环境科学数据中心[20]，该数据是在全国分县人口统计数据的基础上，综合考虑了与人口密切相关的土地利用类型、夜间灯光亮度、居民点密度等多因素，利用多因子权重分配法，将以行政区为基本统计单元的人口数据展布到空间格网上，从而实现人口的空间化.该数据分辨率为1 km.该数据集从1995年起始，每五年统计一次，其他年份数据则由已有的人口数据和人口增长率换算得到.此外，该数据集不包含香港、澳门、台湾地区.以2015年为例，人口分布如图2.为了获得与污染物相同空间尺度的暴露人口分布，本文将人口公里格网分布转换为了人口5 km×5 km格网分布.

图2 2015年全国人口1 km×1 km空间分布Fig.2 Spatial distribution of national population in 2015

1.2 健康效应终点及其基线发生率

基于国际疾病分类标准，考虑到各疾病发病率和死亡率资料的有限性，本文选用的健康终点分为死亡、住院两类.其中，死亡健康终点包含心脑血管疾病死亡、呼吸系统疾病死亡和肺癌死亡；住院健康终点包含心脑血管疾病住院、呼吸系统疾病住院.2000年基线发生率源于《2001年卫生与计划生育统计摘要》，2003至2017年基线发生率源于当年《卫生与计划生育统计年鉴》，且由于未发表《2002年卫生与计划生育统计摘要》，故2001年的基线发生率以2002年的统计数据代替.此外，由于统计标准更改，2000年相关健康终点的疾病别少于2003及之后的年份.以2015年为例，基线发生率数据如表1.

1.3 选择暴露—反应系数

大气暴露—反应系数用于定量描述大气污染物浓度变化和人体健康效应的关系，是大气污染健康影响研究的关键之一.选取流行病学暴露—反应关系研究文献应有2方面的考虑：一是随时间推移我国大气污染物浓度水平和成分会发生变化；二是国内外人群中暴露—反应关系存在差异.鉴于此本文尽量选取研究年份接近本文研究年份，且由国内研究报道的暴露—反应关系分析结果.最终选取的各健康终点的暴露—反应系数平均值及其95%置信区间如表1所示[10，21-22].此外，本文取用暴露—反应系数的平均值参与健康影响模型的计算.

表1 2015年各健康终点的基线发生率和暴露—反应关系系数Tab.1 Reference incidence rates and dose-response relationship coefficients of different health outcomes

1.4 人群健康影响估算

基于暴露—反应关系和泊松回归的比例模型，可以使用以下公式计算各个健康终点归因于PM2.5污染的人群健康影响[23]：

ΔYi=Yi(1-e-βiΔX)×POP,

(1)

其中，ΔYi为健康终点i的人群归因健康影响(例数)；Yi为健康终点i的基线发生率；βi为健康终点i的污染物浓度变化和健康效应的关系系数(即暴露—反应系数)；ΔX为PM2.5基准浓度与控制浓度之差；POP为暴露人口数.计算出暴露人群各个健康终点的健康影响后，各健康终点的健康影响相加即为PM2.5污染的人群总健康影响.

1.5 经济损失估算

健康终点的人群健康影响乘以该健康终点的单位经济损失，即该健康终点的经济损失，各健康终点的经济损失相加即PM2.5污染的总经济损失.对于早逝的单位经济损失估计，本文通过统计生命价值法(value of statistical life，VSL)计算.VSL是支付意愿法(WTP)的一种，VSL并非估计一个具体的人的价值，而是衡量人群为规避一定的死亡风险的而愿意付出的价值[24].由于VSL是人们对风险和财富之间的权衡，因此会受到个体收入水平的影响，在不同省份和同一省份的不同年份，VSL值有所不同.本文以谢旭轩[25]对2010年北京市统计生命价值的研究成果168万元为基础，结合省份之间的人均可支配收入比值，通过下式折算得到其他年份各省的统计生命价值：

(2)

式中，VSLit为省份i在t年的单位统计生命价值；Iit为省份i在t年的人均可支配收入；e为收入弹性系数，本文取1.

对于住院的单位经济损失估计采用疾病成本法(cost of illness，COI).单位疾病成本等于单位医疗费用加上误工损失，其基本的计算公式如下：

Ci=(CPi+GDPi·TLi)·Ii，

(3)

式中，Ci为PM2.5污染造成的健康终点i的单位疾病成本；CPi为健康终点i的单位病例的医疗费用；GDPi为地区人均生产日均值；TLi为因健康终点i导致的误工时间(d)；Ii为PM2.5污染造成的健康终点i的健康效应变化量.

本文参考北京市2013年相关疾病住院的单位疾病成本[26]，通过其他年份各省份人均可支配收入与2013年北京市人均可支配收入的比值系数，得到其他年份各省份的健康终点的单位经济损失.

表2 北京市各健康终点单位经济价值Tab.2 Unit economic cost of different health outcomes in Beijing

1.6 EMD距离

EMD(Earth Mover’s distance)距离是度量图像之间的相似性的一种方法.一般来说，一张图像可以转换成一个直方图分布，在EMD方法中，直方图中的每一根柱称为特征量，柱高称为特征的权重.EMD就是计算将一个直方图分布转换为另一个直方图分布所需要的最小代价.本文利用EMD距离比较不同年份的健康影响和经济损失空间分布的相似度.EMD的计算可用如下公式表示：

(4)

式中，P和Q个直方图分布，包括特征量pi、qi特征权重wpi、wqi；dij是一个矩阵，每一项表示pi、qi差值的绝对值；fij也是一个矩阵，每一项表示从pi到qi的移动量.以最小工作量为目的，求出每个两特征量之间的最小移动量fij，再进行归一化，就可以得到两个分布之间的EMD距离.在图像对比分析中，EMD值越小，表明两个直方图分布差异越小，图像越相似.本研究中，EMD值越接近零，表示不同年份的健康影响和经济损失空间分布模式越近似.

2 结果及分析

通过健康效应模型及经济损失模型，可以计算出各格网内不同健康终点的健康影响和经济损失，对不同健康终点的健康影响及经济损失分别求和，得到该格网的总健康影响和经济损失，再经统计得到2000年—2017年各年份的全国健康影响和经济损失，如图3所示.将全国健康影响和经济损失划分为9个区间，计算各区间内格网数量占比，如图4所示.各年份的全国健康影响空间分布图，如图5所示；各年份的全国经济损失空间分布图，如图6所示.此外，本文利用EMD距离定量分析不同年份之间健康影响和经济损失的空间分布差异，计算得到2001年—2017各年与2000年的健康影响和经济损失的空间分布相似度，以及2001年—2017各年与前一年之间的健康影响和经济损失的空间分布相似度，结果如图7.

2.1 人群健康影响历年变化

由图3分析可知，2000年—2017年，PM2.5污染的健康影响整体可分为两个阶段：2000年—2013年呈波动上升趋势；2013年—2017年基本为逐年下降趋势(2017年上升).健康影响人数最小为2000年的9.07万例，最大为2013年的71.35万例，2000年—2017年的平均值为38.51万例.如图4所示，在2000年—2017年间，健康影响人数为0例的格网达50%以上，侧面反映出我国一半以上地区年均PM2.5浓度未超过二级标准的限值35 μg·m-3；此外，绝大部分格网的健康影响人数位于0～30例/25 km2的区间内，大于30例的格网占比极少.结合统计数据来看：1) 健康影响人数为0例的格网占比从2000年的73.01%降至2017年的54.18%，平均每年占比62.39%.2) 健康影响人数大于5例的格网占比从2000年的1.15%增至2017年的9.15%；健康影响人数大于30例的格网占比从2000年的0.05%增至2017年的0.95%；3) 健康影响人数在局部年份现了超高值，2010至2017年均出现了超过300例的格网，2013和2015年出现了超过900例的格网.

图3 2000年—2017年全国PM2.5污染的人群健康影响和经济损失变化趋势Fig.3 PM2.5-related health impact and economic loss in China from 2000 to 2017

图4 2000年—2017年全国分区间PM2.5污染的健康影响和经济损失变化趋势Fig.4 PM2.5-related health impact and economic loss by range in China from 2000 to 2017

得到每一年的健康影响空间分布(图5)后，对不同年份的分布差异进行定量衡量.图7为两种EMD距离的变化趋势，其中，从基于2000年的EMD值来看，2001年—2007年的健康影响空间分布变化大致可分为两阶段：1) 2001年—2009年与2000年的健康影响的分布差异基本逐年增大(2002年稍下降)；2) 2010年—2017年与2000年的健康影响分布差异呈波动变化趋势.期间，2010年—2012年和2015年—2016年两段时间与2000年的健康影响分布差异逐渐减小；2013年—2014年和2017年两段时间与2000年的健康影响分布差异逐渐增大.从与前一年的EMD值来看，大部分年份健康影响分布的变化较为平缓，在短时间内健康影响分布变化大的年份有2001年、2013年、2016年和2017年.结合图5可知，2001年、2013年和2017年相较于前一年的健康影响均大幅增大；2016年的健康影响相较于2015年则大幅减小.

通过图5进一步分析局部的健康影响分布的变化：1) 2000年—2009年多地区健康影响呈增长趋势.2009年东北地区南部、山东、山西、安徽、湖南、江西、江苏、浙江、上海、重庆、贵州、两广地区的健康影响的分布范围较2000年明显扩大；2009年京津冀南部、河南省东部、山东省西部地区的健康影响对比2000年有较大增幅，高密度区域也逐渐扩大.2) 2011和2012年多地区的健康影响相较于2009年逐渐下降：东北地区南部、湖北省西南部、广东省中部地区的健康影响范围明显缩小；京津冀南部、河南省、山东省的健康影响的高密度区域减小.3) 与2012年相比，2013年北方地区的健康影响增长显著：东北地区南部、山西、陕西、宁夏、甘肃、内蒙古西南部地区的健康影响范围扩大；京津冀南部、山东省西部、河南省东北部地区的健康影响的高密度区域明显扩大.南方地区则有长三角地区的高密度区域扩大.但2013年贵州省东部和广西省北部地区的健康影响范围较2012年缩小.

图5 2000年—2017年全国PM2.5污染的健康影响空间分布(5 km×5 km)Fig.5 Spatial distribution of national health impact of PM2.5 pollution in from 2000 to 2017(5 km×5 km)

图6 2000年—2017年全国PM2.5污染的经济损失空间分布(5 km×5 km)Fig.6 Spatial distribution of national economic loss of PM2.5 pollution in from 2000 to 2017(5 km×5 km)

Fig.7 2000年—2017年全国PM2.5污染的健康影响分布和经济损失分布EMD变化趋势Fig.7 The EMD values of national PM2.5-related and economic loss from 2000 to 2017

此外，2013年人群健康影响是2000年—2017年的首位，占18年来健康影响总数的10.29%.4)2014年—2016年健康影响总数逐渐下降，2015年和2016年较2000年的健康影响分布差异连续缩小.2016年陕甘宁三省、湖南和湖北省西部、两广地区的健康影响范围相比于2013年明显缩小，华北平原和长三角地区健康影响的高密度区域缩少.2013年也是“大气十条”政策的推行之年，可认为该政策对我国大气污染的防治起到了一定程度积极效用.5) 值得注意的是，2017年多地区健康影响范围有所扩大：西北地区、东北地区南部和两广地区的健康影响范围扩大较明显；长三角和安徽省西北部地区的高密度区域也有所扩大；但京津冀地区的高密度区域有所缩小.PM2.5污染问题仍然值得关注.

2.2 经济损失历年变化

由图3可知，2000年—2017年，全国PM2.5污染的经济损失总值呈现出波动上升趋势.经济损失最小值为2000年的61.93亿元，占当年GDP总值的0.06%；最大值为2017年的1 961.1亿元，占当年GDP总值的0.24%.若按经济损失占当年GDP的比重排序，则最小值为2000年的0.06%，最大值为2014年的0.26%.由2000年—2017年经济损失的各区间格网占比(图4)可看出：由于经济损失模型和健康效应模型相关联，每年经济损失为0的格网也超过了50%；18年间绝大多数格网的经济损失分布在0～500万元/25 km2的区间内；大于500万元的格网占比极少.结合具体统计数据可知：1) 经济损失为0的格网占比与健康损失为0的格网占比相同；2) 经济损失高于10万元的格网占比从2000年的4.14%增至2017年的21.5%；经济损失高于100万元的格网占比从2000年0.1%增至2017年的9.89%；3) 经济损失也在局部地区出现了超高值，2013至2017年均出现了经济损失超过4亿元的格网.

经可视化后，5 km分辨率的经济损失空间分布如图6所示，通过两种EMD距离定量分析不同年份的经济损失分布差异如图7所示.从基于2000年的EMD可以看到：首先，仍然可以2009年为界，将2001年—2017年与2000年的经济损失分布差异划分为两个变化阶段.然后，我们也注意到个别年份相较于2000年的经济损失分布差异与健康影响分布差异不一致的情况，如2003年和2013年差异更大，2008年差异较小.这是因为不同地区间的经济水平差距拉大了这些地区的经济损失的差距，造成同一年份不同地区的经济损失差异与健康影响的差异的不同.从各年份与前一年的经济损失分布EMD值可以看出，与前一年分布差异大的年份仍然包括2001年、2012年、2016年和2017年.此外，2003年与2002年、2004年与2003年的分布差异也比较大.较于2002年，2003年华北平原和长三角地区的经济损失有较大增幅；较于2003年，2004年长江以南多地区经济损失范围扩大.造成以上年份与前一年经济损失差异与健康损失差异有所不同的原因是，随时间推移，经济水平的增长使同一地区在不同年份的经济损失差距增大，以至于经济损失差异大于健康影响差异.

2000年—2017年，经济损失局部的空间分布的变化可从图6可知：一方面，由于经济损失模型和健康效应模型之间存在关联，经济损失的空间分布一定程度上与健康影响的空间分布是一致的，因此历年各地区经济损失范围的阶段性变化可以参考健康影响范围的变化.另一方面，随时间推移，各地区经济水平的增长使单位经济损失逐渐增大，造成各地区经济损失高值区域向低值区域的扩张.据统计，18年内，经济损失大于50万元/25 km2的区域基本呈逐年增大趋势(2004年、2012年、2016年有所下降)，这些高值区域主要位于华北平原、长三角、湖南和湖北省东部、四川盆地地区.

2.3 健康影响及经济损失的空间分布格局

由图5、图6可看出，中国PM2.5污染的健康影响和经济损失空间分布具有相似性，两者空间分布的区域性显著，总体上表现为由东南及西北地区向西南、东北地区递减的格局.

健康影响的低密度区主要位于福建、云南、海南、四川省西部、青海、西藏自治区、天山以北地区、内蒙古自治区东部及黑龙江省北部；健康影响的高密度区则表现出明显集聚特征，三大平原及其周围地区、四川盆地、塔里木盆地和两广地区是主要集聚区域.经济损失的低值区域与健康影响的低密度区域基本一致，经济损失的高值区域则更集中分布在京津冀南部、华东、华中、四川盆地及珠三角地区，这也与我国经济空间格局相适应.基于健康影响和经济损失分布的空间集聚性，按地理区划统计各区域的健康影响及经济损失，由图8可知：1) 2000年—2017年，华北、华中、华东地区的健康损失始终位居前三，三区域健康损失平均每年占比74.21%.随时间推移，各区域的人群健康影响差异不断扩大，18年内华东地区增幅最大，其次为华中地区.2) 对于经济损失，2000年—2017年仍然是华中、华北、华东地区位居前三，三区域经济损失和平均每年占比80.78%.18年内华东地区增幅最大，其次为华北地区.

Fig.8 2000年—2017年各区域健康影响和经济损失变化趋势Fig.8 PM2.5-related health impact and economic loss in different regions from 2000 to 2017

3 讨论

在本研究中，健康影响及经济损失估算结果与暴露人口数、PM2.5暴露浓度、健康终点的发生率和暴露—反应关系系数的选取均有关系，因此也存在如下不确定性.1) 本文所采用的PM2.5浓度数据是基于室外环境，但现实中人们不会长时间暴露于室外的PM2.5污染下，因此，对于PM2.5室内污染的健康经济效应的评估也有必要.此外，选取不同的控制浓度会导致估算结果的不同.2) 在健康终点的选取中，出于数据获取难度，本文未予包括门诊、急性和慢性支气管炎、哮喘等健康终点，同样未考虑PM2.5污染引起认知功能下降、产生不良生殖等.3) 暴露—反应系数的选择是评估健康影响的关键要素之一，但出于资料限制，本文使用部分区域暴露—反应系数的研究结果来代替全国，可能使估算结果产生一定偏差，但参数本身的不确定性对各年份的影响是一致的，因此对各年份PM2.5污染损失结果的相对大小无影响.4) 对于社会而言，老年人的过早死亡给社会带来的预期损失小于青年人或小孩.穆泉等[27]对中国13年内PM2.5重污染损失的研究显示，未考虑过早死亡的年龄结构对健康损失的影响，有可能高估PM2.5污染所致的健康损失额2倍以上，因此年龄结构对健康损失评估的影响应纳入考虑.