1. 引言
随着城市化进程的加速和工业化的迅猛发展,城市空气质量已成为全球关注的焦点问题。空气质量不仅直接关系到居民的身体健康和生活质量,还对城市的可持续发展产生深远影响。在中国,由于地域广阔且各城市的产业结构、能源利用方式以及自然地理条件存在显著差异,城市空气质量呈现出多样化的特征。
根据城市的发达程度、地理位置和人口密度等因素不同,选取北京、福州和乌鲁木齐这三个城市作为代表,绘制了2019年~2023年空气质量情况的折线图,如图1所示。从2019~2023这五年间,三个代表城市在各污染指标上呈现出不同的变化态势。北京作为我国的首都,污染物的平均浓度相对中等发达地区(福州)和一般水平地区(乌鲁木齐)要略高一些,而一年之中达到或好于二级的天数明显低于另外两个城市,其中达到或好于二级天数最多的城市是福州。三个城市中区别较为明显的污染物指标为O3和PM10。这五年北京在O3这一指标上日最大8小时第90百分位浓度高于150,而福州和乌鲁木齐均低于150。而PM10上却是乌鲁木齐的年平均浓度全部高于50,高于北京和福州的数据。三个城市其余污染物的年平均浓度基本都在50以下。由此可以看出在不同因素作用下的地区在空气质量的不同指标上同时兼具差异性和相似性。
Figure 1. Air quality in 2019~2023
图1. 2019~2023年空气质量
国内外学者用不同思路和方法对空气质量情况进行了研究。国外学者Kunwar P. Singha [1]等人利用主成分分析和决策树理论对空气质量问题进行了研究;国内孙越嘉[2]利用我国31个主要城市的空气质量影响因子SO2、NO2、PM10、CO、O3、PM2.5进行了聚类分析和因子分析;杨晓艳[3]等人使用模糊评价方法对各地空气质量进行评价。基于以上研究,本文使用SPSS对我国31个主要城市的7个主要指标SO2、NO2、PM10、CO、O3、PM2.5和二级天数进行主成分分析、聚类分析和判别分析。
2. 数据来源
本文数据为2022年我国主要城市空气质量情况,来源于2023年《中国统计年鉴》[4],选取了7个描述空气质量的相关指标,样本包含我国31个主要城市。
Figure 2. Histogram of indicators for different cities
图2. 不同城市指标柱状图
由图2可知,我国主要城市在SO2、NO2、PM10和PM2.5这些指标上差异较大,在CO、O3和二级天数这些指标上具有一定的相似性。
3. 方法原理介绍
3.1. 主成分分析
主成分分析法借助于正交变换将多个相关性指标转化为少数独立的综合指标,起到对数据降维的作用,每个综合指标都是原始指标的线性组合[5] [6]。综合指标即为主成分,其信息量用方差表示,这些主成分按照方差从大到小的顺序排列,方差越大,说明该主成分包含的原始数据信息越多[7]。通常保留较少的主成分来代替原来较多的变量。
3.2. 聚类分析
本文使用聚类分析中的K-均值聚类法,其主要用于将数据集中的样本划分为K个不同的聚类,通过迭代的方式,将n个数据对象划分为K个类,使得每个聚类内的数据对象之间的相似度较高。通常使用欧式距离作为相似度度量,目标是最小化每个数据点到其所属聚类中心的距离平方和。
4. 主要结果
为了消除量纲和数量级不同所可能带来的分析影响,先对数据进行标准化,然后对数据进行分析。
4.1. 相关性分析
本文通过对我国31个主要城市的六种空气污染物和二级天数共7个变量做Pearson相关性分析,得到变量之间的相关系数热力图如图3。相关系数取值范围为[−1,1],相关系数的绝对值越接近于1,变量之间的相关性越强。
根据热力图中方块不同颜色,可以判断出变量之间的相关性强弱。从图中可以看出,SO2,NO2,PM10,CO,PM2.5之间存在不同程度的正相关,其中NO2和PM10,PM10和PM2.5存在较强的相关性。二级天数与6种空气污染物均呈负相关,表明这些污染物浓度越高,空气质量达到二级的天数更少。
Figure 3. Correlation coefficient heat map
图3. 相关系数热力图
4.2. 主成分分析
为了进一步了解不同变量对城市空气质量的影响,我们对数据进行主成分分析。通过绘制方差解释柱状图和累计方差折线图来直观得到提取的主成分个数。
Figure 4. Cumulative contribution of variance
图4. 方差累计贡献率
从图4可以看出,前两个主成分方差解释能力显著高于之后的主成分,二者的累计方差贡献率达到了88.024%,充分保留了数据的主要特征。从第三个主成分往后,累计方差解释折线的变化速率减缓,表明其解释能力下降,因此,选择前两个主成分。
根据提取的2个主成分,计算因子载荷矩阵,通过对因子载荷矩阵除以各主成分特征值的平方根,得到主成分系数。关于主成分系数的雷达图如图5所示。可以看出第一主成分在6个污染物浓度上的系数均为正,第二主成分在CO和SO2上系数为正,其余指标上系数为负数或接近0。
Figure 5. Radar chart of principal component coefficients
图5. 主成分系数雷达图
两个主成分关于标准化原始变量的系数表达式为:
从第一主成分各指标的系数可知,第一主成分在所有污染物浓度变量中系数全部为正数,在二级以上天数上的系数为负数,因此第一主成分代表污染物浓度的综合水平,直接体现污染程度的高低。
从第二主成分各指标的系数可知,在SO2和CO上的系数为中等大的正数,在臭氧浓度上的系数为负数,表明第二主成分表示气体污染(CO和SO2的浓度)与光污染(O3浓度)的差异。
基于主成分分析得出的两个综合评估指标是对原始数据的数据降维,因此,通过以上两个综合指标评估主要城市的空气质量更具有代表性。
进一步计算出主成分得分。两个主成分得分的柱状图如图6所示。
Figure 6. Histogram of principal component scores
图6. 主成分得分柱状图
从图6中第一、二主成分得分可以看出,各城市在不同主成分上表现差异明显。主成分得分的正负反映了城市在对应主成分所代表的综合指标上相对于平均水平的偏离情况。
石家庄、太原、济南、郑州、西安、兰州这些城市在第一主成分得分较高且为正数,说明这些城市整体上污染程度较高。
哈尔滨、兰州、西宁这些城市在第二主成分得分较高且为正数,说明在第二主成分所代表的综合指标上表现突出,可能在某些关键空气质量指标或相关因素上具有优势。这些城市的气体污染物CO和SO2的含量较高,光污染程度较低。
4.3. 聚类分析
根据每个城市的主成分得分进行K均值聚类,通过手肘法得出应聚为3类。通过散点图来直观展示不同城市的主成分得分情况。
Figure 7. Scatter plot of principal component scores
图7. 主成分得分散点图
从图7可以看出,第一主成分解释了65.27%的方差,第二主成分解释了22.754%的方差。聚成的3类特点各异。第一类在第一主成分得分几乎全是负数,第二主成分得分有正有负。第二类在第一、二主成分得分均较高。第三类在第一主成分得分较高,第二主成分得分较低。
将每个类中的城市和各指标的平均值计算得出表1的数据,进一步对各城市的空气质量情况进行分析。
Table 1. Mean and variance for each category
表1. 各类平均值和差异性
聚类 |
城市 |
第一主成分平均得分 |
第二主成分平均得分 |
样本到类中心的平均距离 |
1 |
北京、长春、上海、南京、杭州、福州、南昌、长沙、广州、南宁、海口、贵阳、昆明、拉萨 |
−1.707 |
−0.484 |
1.648 |
2 |
呼和浩特、沈阳、哈尔滨、合肥、武汉、重庆、兰州、西宁、银川、乌鲁木齐 |
0.619 |
1.297 |
1.254 |
3 |
天津、石家庄、太原、济南、郑州、成都、西安 |
2.53 |
−0.889 |
1.065 |
第一类城市在两个主成分上平均得分均小于0,表明这类城市污染物综合浓度最低,且气体污染与光污染差异较小。这类城市以北京、上海、广州等经济发达的一线城市及福州、海口等沿海城市为主,具有较高的财政投入能力和技术创新优势。政策上,北京通过实行“国六”排放标准淘汰老旧柴油车,也是全国首个实施“挥发性有机物(VOCs)排污收费”,覆盖石化、印刷等行业。上海依托长三角一体化战略,外迁高污染企业并推广新能源汽车,降低了中心城区空气污染物的排放。沿海城市则利用季风加速污染物扩散。同时,昆明、贵阳等西南城市依托高植被覆盖率增强生态净化能力。经济实力支撑的长期治理、沿海地理的扩散优势与生态屏障的协同作用,使此类城市空气质量显著优于其他地区。
第二类城市在两个主成分上平均得分均大于0,第二主成分平均得分是三类中最高,表明这类城市污染物综合浓度中等,气体污染和光污染差异性较大,其中气体污染显著高于光污染水平,反映出这类城市以工业排放为主。这类城市包括天津、石家庄、太原等重工业城市以及兰州、乌鲁木齐等地处内陆和盆地地区。政策上,石家庄通过“煤改气”减少散煤污染,太原淘汰落后焦化产能以降低SO2排放,兰州则依赖冬季限产和人工增雪缓解污染。然而,由于重工业城市的工业产值占比高以及工业结构转型缓慢的现状,内陆地区由于其地理位置限制导致的污染物积累,政策多为短期应急措施,难以根治污染问题,导致这类城市污染程度较高且气体污染严重的特点。
第三类城市在第一主成分上平均得分最高,第二主成分平均得分最低,表明这些城市的整体污染水平最高,光污染问题显著。这类城市主要包括我国的华北地区和中西部地区。政策上,京津冀通过区域联防联控统一应急响应,西安划定“禁煤区”推广清洁能源,郑州加强扬尘智能监控。其中,郑州、西安及华北平原城市是典型代表,其工业化速度快、能源依赖煤炭,同时受地形气候影响(如华北冬季逆温、中西部盆地扩散受限)。但重工业排放强度高、平原逆温频繁及夏季高温加剧O3生成,导致PM2.5与O3协同升高。产业与气候的双重压力使此类城市污染治理难度最大,反映了其高污染与复杂污染类型的特征。
通过样本到每个类别中心的平均距离可以看出,第三类城市间的空气质量差异较小,第一类城市间的空气质量差异较大。
4.4. 判别分析
接下来用判别分析检验分类结果的准确率,判断之前的分类是否相对准确、可靠。关于主成分得分的分类结果判断如表2。
Table 2. Results of discriminant categorization based on principal components
表2. 基于主成分的判别归类结果
|
|
个案聚类编号 |
预测组成员信息 |
总计 |
|
|
|
1 |
2 |
3 |
|
原始 |
计数 |
1 |
14 |
0 |
0 |
14 |
2 |
1 |
9 |
0 |
10 |
3 |
0 |
0 |
7 |
7 |
% |
1 |
100 |
0 |
0 |
100 |
2 |
10 |
90 |
0 |
100 |
3 |
0 |
0 |
100 |
100 |
交叉验证b |
计数 |
1 |
13 |
1 |
0 |
14 |
2 |
1 |
8 |
1 |
10 |
3 |
1 |
0 |
6 |
7 |
% |
1 |
92.9 |
7.1 |
0 |
100 |
2 |
10 |
80 |
10 |
100 |
3 |
14.3 |
0 |
85.7 |
100 |
A. 正确地对96.8%个原始已分组个案进行了分类。B. 仅针对分析中的个案进行交叉验证。在交叉验证中,每个个案都由那些从该个案以外的所有个案派生的函数进行分类。C. 正确地对87.1%个进行了交叉验证的已分组个案进行了分类。
如表2所示,用K均值聚类法对城市进行分类后,分类器表现较好,交叉验证得到的正确率为87.1%,说明分类结果准确率较高。
5. 结论
针对我国主要城市的空气质量问题,通过主成分分析、聚类分析和判别分析三种方法对2022年数据进行了分析。得到了以下结论:
主成分分析得到的评估我国主要城市空气质量的两个主成分,第一主成分在所有污染物浓度变量(如SO2、NO2、PM10、CO、O3、PM2.5)上的系数均为正数,而在二级以上天数上的系数为负数,表明第一主成分主要反映了污染物浓度的综合水平,直接体现了城市空气污染程度的高低。第二主成分在SO2和CO上的系数为中等大小的正数,而在臭氧浓度上的系数为负数,表明第二主成分主要反映了气体污染(以SO2和CO为代表)与光污染(以O3为代表)之间的差异。综合来看,第一主成分可用于评估城市空气污染的整体严重程度,而第二主成分则有助于区分不同类型的污染特征(气体污染与光污染)。
本研究通过主成分分析与聚类分析,将中国城市空气质量划分为三类典型模式:
(1) 经济–自然协同治理型(第一类):以北京、上海及沿海城市为代表,经济实力与地理优势支撑了长期低污染水平,其政策创新(如VOCs排污收费、新能源汽车推广)与季风扩散、生态屏障形成协同效应,验证了“财力–技术–自然赋能”模式的有效性。
(2) 工业–地理约束型(第二类):以天津、石家庄等重工业城市及内陆盆地城市为主,工业惯性叠加地理扩散限制,导致气体污染突出,短期应急政策难以突破“高排放–低效益”路径依赖。
(3) 产业–气候耦合型(第三类):以郑州、西安及华北城市为典型,重工业与气候敏感性的叠加导致污染综合水平最高,PM2.5与O3协同升高,凸显“产业排放–气象条件”的深度耦合风险。
针对三类城市空气情况现状,提出以下治理建议:
第一类城市:利用经济优势推广碳交易和新能源技术,优化沿海工业布局借力季风扩散,通过生态补偿机制带动周边地区协同减排。
第二类城市:加速重工业循环化改造(如氢能炼钢),建设山地通风廊道缓解地理限制,推动“清洁能源 + 智能气象干预”替代短期管控。
第三类城市:建立污染弹性响应机制(静稳天限产、高温天控臭氧),强化PM2.5与O3协同治理,深化跨省联防联控与污染赔偿制度。
综上所述,相较于已有研究,本文将主成分分析与空间聚类结合,识别中国城市空气污染的“经济–地理–政策”耦合模式,突破了传统单维度分类。现有文献多基于单一指标划分污染物类别,本文综合主成分得分定义污染结构差异,分为气体污染和光污染,弥补了传统指标体系的片面性。根据聚类特点为城市空气质量分类治理提供了科学依据。
基金项目
北京信息科技大学2023年校级教学改革项目:基于模块化和案例的《多元统计分析》教学探索,5112310826。北京信息科技大学2024年研究生课程思政项目:思政“微故事”融入高等数理统计的教学模式探索,2024PYSZ04。