近年来全球多地遭遇罕见极端天气, 给人们的生产生活造成了极大影响[1, 2].全球极端天气频发的主要原因之一是气候变化, 自工业革命以来人类社会排放的大量温室气体正在重新塑造地球气候, 过量的温室气体使得地球气候系统的“能量”过剩, 而这些过剩“能量”最终以极端天气的形式进行释放[3].中国是受全球气候变化影响最为显著的国家之一, 干旱和洪涝等灾害频发, 但同时中国也是世界上最大的能源消费国, 二氧化碳排放总量居世界第一位, 单位国内生产总值能耗为世界平均水平的2倍, 人均碳排放量则超过世界平均水平的40%[4, 5].面对巨大的气候环境压力和减排压力, 近年来中国制定了一系列的降碳减排计划, 习近平总书记也向世界做出了力争在2030年前碳排放达峰, 2060年前实现碳中和的“中国承诺”.基于这一目标, 如何科学准确地评估中国碳排放情况成为了实施“碳达峰”“碳中和”相关政策的前提和迫切需求.

近年来, 学者们就中国的碳排放情况进行了大量科学研究.主要有以能源消费统计数据为基础估算国家尺度的碳排放并计算碳排放强度[6~9]和人均碳排放强度等[10]; 基于社会经济数据、产业结构数据和能源统计数据等, 研究影响中国能源碳排放的驱动因子[11~15]; 此外, 随着遥感技术的发展, 大量的遥感数据产品也被用于对碳排放估算的研究, 其中夜间灯光作为对人类社会活动最直接的体现, 已被证明可用于经济社会指标估算、城镇化进程研究和重大事件影响评估等, 因此不少学者在碳排放的研究中引入了夜间灯光数据, 诸如基于夜间灯光、土地利用数据估算中国栅格尺度碳排放[16], 如苏泳娴等[17]基于1992~2010年DMSP/OLS夜间灯光结合省市级的统计数据估算了中国大陆碳排放量, 并分析了碳排放密度以及影响机制, 发现经济增长是二氧化碳排放增长的主要因素, 能源结构、能源利用效率和产业结构则是碳排放强度变化的主要因素; 王少剑等[18]基于1992~2013年DMSP/OLS夜间灯光数据反演的中国地级能源碳排放, 揭示了经济增长和第二产业占比升高对碳排放增长呈现正向效应, 人口集聚、技术水平、对外开放度和公路运输强度则表现为负向效应; 除此之外还有从流域视角研究碳排放格局和驱动因素[19]以及基于碳排放格局、变化趋势和主要驱动因素进行的碳排放交易[20]、碳减排政策[21~24]和减排路径[25~28]的研究.

上述关于中国能源碳排放的相关研究, 为摸清碳排放家底和实施具体的碳减排措施提供了思路和方法, 但仍存在3点不足：①单一的基于能源消费统计数据的研究, 由于缺乏市级和更小尺度的数据以及不同尺度统计数据的统计口径、计算方法和统计误差等因素影响, 使得现有研究大多聚焦于国家和省级尺度, 对于市级以及更小尺度碳排放情况研究则较为稀缺; ②基于夜间灯光数据的碳排放研究由于DMSP/OLS数据与NPP/VIIRS数据时间范围不同且两种数据存在较大差异, 使得长时间连续的碳排放监测研究面临困难; ③单一的基于夜间灯光数据的研究, 由于夜间灯光与碳排放本身并无直接因果关系, 在碳排放模拟结果尤其是在较小尺度上的模拟结果受到了部分学者的质疑.

因此, 本文在前人研究的基础上, 通过一系列科学方法力争解决上述问题.首先, 通过线性回归的方法重构了长时间序列可相互比较的夜间灯光数据集; 其次, 通过人口栅格数据的引入以及基于人均碳排放和单位DN值所承载的碳排放进行数据分组, 最终估算了2000、2005、2010、2015和2018年中国碳排放, 并对其时空特征、演变特征和城市碳排放与经济增长的关联性进行了分析, 以期为实施碳减排措施和实现“碳达峰”和“碳中和”目标提供科学支撑.

本文中用到的数据主要有5类, 各类数据信息及来源如表 1所示.

夜光数据包含DMSP/OLS和NPP/VIIRS两套数据的处理, 其中DMSP/OLS稳定夜间灯光数据去除了月光、火光和油气燃烧的噪声影响, 但存在像元饱和, 本文采用增强型植被指数对其进行去饱和度处理[39], 并进行年内合成校正与年际序列校正.对于NPP/VIIRS数据本文采用Zhao等[40]提出的方法对其进行异常值去除、年内合成校正与年际序列校正.从原理上讲, 不同的传感器对同一时间同一地点采集的夜光数据应大致相同, 因此本文以DMSP/OLS数据为基准, 把NPP/VIIRS数据作为校正对象, 对2012年和2013年中国省级尺度的DMSP/OLS与NPP/VIIRS数据总亮度值(total digital number, TDN)进行回归分析.由图 1可知两组数据存在较高的相关性(R2=0.958 9), 可以利用该方法进行数据拟合构建长时间序列连续的夜间灯光数据集.LandScan人口数据主要采用了图像分析以及多变量不对称的建模方法来分解区域内的人口普查数据.

本文基于省级的能源消费统计数据计算碳排放量.具体计算方法是利用美国能源部橡树岭国家实验室(ORNL)提出的化石能源燃烧产生二氧化碳的计算公式[41]以及政府间气候变化专门委员会(IPCC)提供的各类能源碳排放系数, 选取了9种主要能源的消费统计数据计算CO2排放量, 计算公式如下：

式中, i表示能源类型, 共计9种; Ei表示能源i的消费量, 按标准煤计(万t); Ki表示能源i的碳排放系数(万t碳/万t标准煤).各类能源的换算标准煤系数和碳排放系数见表 2.

(1) 统计碳排放数据、夜光数据和人口数据分组由于各省份的社会经济发展程度不一, 导致各地人均碳排放以及单位DN值所承载的碳排放均有较大差异[42, 43], 利用全国所有的数据进行一次碳排放的回归分析会造成一定程度的失真, 因此本文对统计碳排放数据、夜光数据和人口数据进行二次分析, 以实现对原始数据符合实际的分组来减小统一回归分析的失真程度.具体方法是统计各个年份省级的DN值总值及总人口, 分别计算各年份各省份的人均碳排放和单位DN值所承载的碳排放, 计算公式如下：

式中, i表示不同年份, j表示不同省份; CO2ij表示i年j省的统计碳排放量, POPij表示i年j省的总人口数, DNij表示i年j省的DN值总值; Aij表示i年j省的人均碳排放量, Bij表示i年j省的单位DN值所承载的碳排放量.

利用自然断裂法分别对Aij与Bij进行两组断裂分组, 综合Aij与Bij可以将5 a 30个省份共计150组数据分为3类(命名为ab型, a表示单位DN值所承载的碳排放, b表示人均碳排放, 1表示低值碳排放, 2表示高值碳排放), 分组情况见表 3.

(2) 能源消费碳排放模拟利用分组后的3组原始数据导入SPSS进行不含截距的线性回归分析, 结果显示统计碳排放数据与夜光数据和人口数据具有良好的线性相关关系, R2均大于0.9.碳排放回归方程参数见表 4.

因此可以对中国的能源消费碳排放进行估算, 计算公式如下：

式中, y表示碳排放的估算值, x1表示夜间灯光DN值, x2表示人口数据, a和b分别为各个自变量的回归参数.

(1) 全局空间自相关使用Moran's I刻画属性值在整个区域内的空间相关性, Moran's I值介于-1~1, 大于0表示正相关, 即该属性值在区域内存在强的相似性; 小于0则为负相关.计算公式如下[44]：

式中, n表示地级行政单元数量; xi和xj分别表示地级行政单元i和j的碳排放量; x表示各地级单元碳排放量的平均值; wij表示邻近地级单元i和j的空间权重矩阵.

(2) 局部空间自相关一般利用LISA刻画在整体空间中可能存在的局部空间变异性, 计算公式如下：

式中, 当Ii>0时, 表示邻近区域为正的空间自相关, 为“高高”或“低低”型; 当Ii < 0时, 表示邻近区域为负的空间自相关, 为“高低”或“低高”型.

为反映2000年以来中国大陆地区的能源碳排放变化特征, 本文引入数字编码法, 对5期碳排放估算结果进行两两差值处理, 再对4期差值数据进行数学编码运算, 最终获得一幅体现2000~2018年中国大陆地区能源碳排放变化频次的栅格地图.具体方法是利用栅格计算器分别对2000、2005、2010、2015和2018年的5期碳排放估算结果进行做差运算, 将等于0的赋值为0, 不等于0的赋值为1, 最后将4期差值图按照公式(6)叠加[45].计算公式如下：

式中, code2000~2005、code2005~2010、code2010~2015和code2015~2018分别表示4期差值的编码(0或1), code2000~2018表示2000~2018年中国大陆地区能源碳排放变化频次代码, 变化频次代码共计5类16种(表 5).

本研究结果发现, 基于统计数据的中国能源碳排放存在明显的阶段性差异, 2000~2005年间碳排放增长率为77.27%, 2005~2010年间碳排放增长率为49.03%, 2010~2015年间碳排放增长率为17.49%, 2015~2018年间碳排放的增长率为9.60%, 其中2000~2005年与2005~2010年为碳排放快速增长的两个阶段, 因此选择2005年和2010年两个关键时间节点, 选择深圳、珠海、厦门、忻州和晋城等29个地级市的碳排放统计数据进行模型精度验证.将估算结果与统计结果进行回归发现两者存在明显的线性关系(图 2), 在2005年R2为0.92, 均方根误差RMSE为1 145.56, 在2010年R2为0.90, RMSE为1 272.77.由此可见, 综合夜间灯光与人口数据的碳排放估算模型有着良好的精度, 误差在可控范围内, 可以利用该模型进行中国能源碳排放估算以及时空演变分析的相关研究.

如表 6所示, 从2000~2018年我国的CO2排放总量不断增长, 但增长速率出现下降趋势.碳排放总量从2000年的48.2亿t增长至2018年的153.5亿t, 碳排放年均增长速率则由2000~2010年的9.9%下降至2010~2018年的7.4%.该项数据也反映出由于经济社会发展对能源的需求, 我国的碳排放总量在2000~2018年间快速增长, 但由于后期国家层面开始调整能源结构, 实施降碳减排政策, 并同时研究相关的技术政策, 使得碳排放增速开始下降.

为更好地分析我国碳排放的空间分布格局, 本文将估算结果分为3类：无碳排放(< 0万t)、低碳排放(0~1.16万t)、高碳排放(>1.16万t).结果表明, 无碳排放地区主要分布在西北及东北部分地区, 低碳排放主要分布在广大的中小城镇地区, 高碳排放则集中分布于华北、华中、东部沿海以及西部的省会城市及城市群附近(图 3).西北地区地广人稀, 存在大量的无人区, 如荒漠戈壁、雪山冰川等, 因此碳排放少, 主要的碳排放地区分布在重点城市、交通线附近的节点城市以及人口聚居区, 东北地区的无碳排放主要分布在大兴安岭、小兴安岭和长白山地区; 低碳排放主要分布在广大的中小城镇地区, 因为这里分布了我国94.1%的人口, 产生了不少的低碳排放; 高碳排放分布的华北、华中地区是我国重要的能源基地、工业基地和粮食主产区, 除此之外还分布着京津冀城市群、太原城市群、中原城市群、山东半岛城市群以及长江中游城市群[46], 东部沿海地区是全国经济最发达的地区, 人类活动强烈, 能源消耗量较大且密集分布着大量的制造业, 同时还有长江三角洲城市群、海峡西岸城市群、珠江三角洲城市群, 西部的高碳排放则集中在省会地区附近, 其中尤以成渝城市群和关中城市群最为显著.

通过Moran's I刻画中国能源碳排放空间自相关特征, 发现该指数值均大于0, 表明中国大陆的碳排放呈现较为显著的正相关, 即存在高值聚集或低值聚集现象.具体变化为2000年莫兰指数值为0.22, 2005年降为0.19; 在此之后总体为上升趋势, 2010年为0.23, 2015年值为0.23, 2018年则略有下降为0.22, 从指数变化看, 整体趋于稳定, 表明中国大陆的碳排放高值聚集或低值聚集现象趋于稳定.利用LISA指数对中国能源碳排放进行局部变异性分析发现(图 4), 中国碳排放的低低集聚区正在破碎化, 规模在减小而同时高高集聚区正在加强巩固, 规模开始扩大; 除此之外, 分析集聚区的地域分布还发现地区的社会经济发展水平与碳排放存在密切联系[47].具体表现为2000年在西部形成了新疆-青海-内蒙古西部-甘肃-四川西部连片的低低集聚区, 海南岛也是一个低低集聚区, 此外还有湖南西部、江西和安徽接壤处也形成了低低集聚区; 而在京津冀城市群-太原城市群、潍坊-淄博-泰安-济宁、南通-苏州-嘉兴和中山-东莞-惠州-清远形成了4个高高集聚区; 到2010年, 西部地区的低低集聚区规模大幅减小, 仅剩甘肃西北-青海-四川西部-云南西北部; 高高集聚区中太原城市群加强为忻州-吕梁-临汾-长治-晋中, 长三角地区加强为泰州-南通-无锡-苏州-嘉兴-上海; 2015年, 西部地区的低低集聚区规模进一步减小, 湖南西部的低低集聚区消失; 高高集聚区中长三角地区进一步加强为泰州-南通-无锡-苏州-湖州-嘉兴-绍兴-上海; 珠三角地区加强为清远-惠州-深圳-东莞-中山; 到2018年, 在怀化、天门-荆门-潜江出现新的低低集聚区; 高高集聚区中辽宁鞍山以及长三角地区规模进一步扩大; 其余地区的碳排放空间依赖性则无明显变化.

为进一步探究2000~2018年中国能源碳排放的变化情况, 本文在能源碳排放估算结果的基础上, 利用数字编码法对碳排放变化情况做了量化分析(图 5).本研究结果表明我国大部分地区的碳排放比较稳定, 没有大幅增加或减少, 结合碳排放的空间分布结果可以发现, 这部分地区主要是无碳排放地区和低碳排放地区, 值得注意的是在高碳排放地区也存在未发生变化的地方, 这些地方主要是省会城市和重点城市的中心城区, 例如北京、天津、西安、成都和上海等城市的中心城区, 出现这种现象主要是由于这些省会城市和重点城市的中心城区已经过多年的发展, 空间结构、产业结构相对稳定, 而且中心城区基本没有高碳排放企业, 所以中心城区的碳排放一直相对稳定.碳排放发生1~2次变化和3~4次变化的地区主要分布在省会城市和重点城市的外围城区以及郊区, 其中以发生1~2次变化为主, 例如西安、成都和郑州等地, 这些城市中心城区的外围主要是1~2次的碳排放变化, 出现这种情况主要与城市扩张、人口、产业和企业向市郊县区流动有关; 发生3~4次变化的地区主要由山西、河北、武汉和长沙这些地区的城市外围, 山西和河北的变化主要与能源的直接消耗有关, 例如山西在2000年之后开始大力关停私营煤矿和黑煤窑, 淘汰和合并产能落后的煤矿, 合并整顿国营煤炭企业, 研发降碳技术等; 河北在重度霾天气后开始大力淘汰钢铁产能, 由于煤炭和钢铁这些高碳排放产业的变动, 导致了山西和河北出现了碳排放的多次变化.武汉和长沙这些地区的城市外围碳排放的变化一方面来自于自身城市的扩张、中心城区人口企业的外移, 另一方面则来自于东部沿海地区制造业的内迁, 武汉和长沙等地作为靠近东部沿海地区, 经济基础又相对较好的内陆城市是这些内迁制造业的主要承接城市, 综合这些原因使得武汉、长沙等城市外围出现了碳排放的多次变化.

由于中国的高碳排放地区主要集中在重点城市和城市群附近, 而城市的快速发展对能源的需求是产生碳排放的根本原因, 因此本文在获取了城市碳排放估算结果的基础上, 根据城市地区生产总值和城市碳排放量两个不同维度划分, 运用四分位图刻画了2000~2018年中国城市在发展中能源消耗与经济增长的关联性(图 6).该图横坐标为地区生产总值, 代表经济发展, 纵坐标为碳排放量, 代表能源消耗.其中第三象限的点代表“低能耗-低收入”城市, 第二象限的点代表“高能耗-低收入”城市, 第一象限的点代表“高能耗-高收入”城市, 第四象限的点代表“低能耗-高收入”城市.根据结果可发现, 在2000~2018年中国城市发展进程中, 遵循了由第三象限到第二象限再到第一象限、第四象限的一种顺时针发展规律, 说明中国各城市的发展经历了由“低排放-低收入”转向“高排放-低收入”转向“高排放-高收入”最终转向“低排放-高收入”的发展过程.由此可见, 城市的发展与碳排放有着密切的联系, 同时相关研究也显示城市在无干预与规划干预两种情况下, 城市的发展扩张以及对能源的需求有着较大的差异[48].在当前“碳达峰”和“碳中和”大背景下, 如何科学地干预城市发展并使“低排放-低收入”、“高排放-低收入”和“高排放-高收入”城市快速步入“低排放-高收入”城市行列, 实现城市低碳发展是需要深入探讨的一个命题.

(1) 从2000~2018年, 我国的CO2排放总量持续增长, 但增长速率呈现放缓的趋势.碳排放从空间分布来看, 无碳排放地区主要分布在西北及东北地区, 低碳排放主要分布在除西北以外的广大地区, 高碳排放则集中分布于华北、华中、东部沿海以及西部的省会城市及城市群附近.无论从总量还是空间分布来看, 未来几年东部、中部地区都将面临更大的减排压力, 在碳交易市场建立后东部、中部地区也将向西部和东部地区支付更多资金获取排放额, 此举也将倒逼高排放产业进行减排技术革新.

(2) 碳排放在地级市水平上呈现较为显著的全局空间正相关, 存在高值聚集或低值聚集现象, 且该聚集现象整体趋于稳定.

(3) 中国大部分地区的碳排放量相对稳定, 碳排放发生变化的地区主要分布在省会城市和重点城市的外围地区, 即存在中心城区无变化, 外围区域碳排放变化的圈层结构.城市作为减排工作的基本单元, 中心城区并不是工作中心, 城市外围地区才是降低城市碳排放的重点地区.

(4) 在2000~2018年中国城市碳排放与经济关联性发展的历程中, 整体上遵循一个由“低排放-低收入”转向“高排放-低收入”再转向“高排放-高收入”最终转向“低排放-高收入”的发展规律.