1. 项目背景与分析目标拆解1.1 为什么偏偏是 2014-2024 这十年做空气质量数据分析最怕的就是样本区间太短。只看两三年数据你分不清某个变化到底是长期趋势还是单纯的气候波动。所以当我把目光锁定在 2014 到 2024 这个窗口时其实是在刻意选择一段有特殊意义的完整周期。这个时间窗口的起点很有意思。在那之前公众对 PM2.5 的感知已经很强烈但系统的、大范围的公开监测数据还不够完整。从 2014 年前后开始国家层面的环境监测网络逐步加密公开数据服务接口陆续上线数据质量也慢慢稳定下来。也就是说这十年恰好覆盖了一个从大家都在谈空气污染到手里真正有完整数据可用的转折期。另一方面十年的时间跨度足够长可以平滑掉单年气象条件带来的随机波动让结构性变化浮出水面——比如某类污染物是在真正下降还是仅仅因为那年风大。就 AQI空气质量指数本身而言它并不是一个可以直接拿来就用的数字。AQI 是把六项污染物——PM2.5、PM10、SO2、NO2、CO、O3——按照各自的分段 IAQI 换算后取最大值得到的综合指数。每项污染物都有一套线性插值公式比如 PM2.5 的 24 小时平均浓度在 35 到 75 微克每立方米之间时对应的 IAQI 段是 50 到 100。取所有污染物的最大 IAQI 作为当日 AQI同时把贡献最大那个污染物标记为首要污染物。这个设计逻辑决定了AQI 下降不等于所有污染物都在下降它只代表当天最突出的那个问题缓解了。理解这一点是后面所有分析的地基。如果你没搞清楚 AQI 的取最大值机制就很容易犯一个错误看到 AQI 年平均值下降了就得出空气质量全面好转的结论。实际上臭氧单独拉高的夏天AQI 可能一点不低。1.2 这个项目到底想回答什么问题我搭这个分析项目不是纯粹为了画几张漂亮的趋势图。整理数据之前先问自己三个问题后面所有处理步骤都是围绕这些问题展开的。第一个问题这十年里整体空气质量到底有没有结构性改善注意结构性这三个字。单看某一年降低说明不了什么要看趋势线是否持续向下、下降斜率在哪个阶段变缓、不同污染物之间的下降幅度是否同步。第二个问题季节规律有没有发生变化传统认知里秋冬季是重污染高发期因为采暖排放叠加不利扩散条件。但十年里采暖结构、工业布局都在调整这个规律是否被打破夏季臭氧问题是不是在逐年加重这些都需要从十年尺度的数据里找答案。第三个问题极端污染事件是变多了还是变少了这一点对普通人最实际。年均 AQI 从 100 降到 80听起来是好事但如果你关心的是家里老人冬天能不能出门锻炼那你更该关注的是重度污染天数这个指标而不是年均值。年均值会被大量好天气稀释掉而极端天数才是真正影响健康防护决策的变量。这三个问题分别对应趋势分析、季节分解、极值统计三套方法。整个项目就是围绕这三条线展开的后面每一节都会回到这三条线上来。注意这个分析适合谁如果你手里有一份多年环境数据想做回顾分析或者你想验证某个城市空气质量是否真的变好了这种日常争论又或者你想练手时间序列分析与数据可视化这套流程都能直接抄作业。需要的基础只有 Python 和 pandas 的基本操作不需要机器学习背景。2. 数据获取与预处理栈2.1 数据源选型公开平台批量下载不做爬虫数据是整个项目的地基。这一步选错了后面全白搭。我采用的是某国家级环境监测数据发布平台提供的公开历史数据以及部分第三方整理好的开放数据仓库。两套来源都测过最终以官方平台为主、第三方数据为辅交叉验证后再入库。为什么不建议自己写爬虫两个原因。一是大部分这类平台的页面反爬机制虽然不复杂但逐日、逐站点的抓取量太大容易被限制没必要为这点数据冒风险二是历史数据本身就支持打包下载批量获取的效率和稳定性远高于实时爬取。直接下载 CSV 压缩包解压后是结构化表格省掉了解析 HTML 的时间。数据粒度上我选择了日尺度。原始数据里有逐小时的但考虑到 AQI 本身就以日为单位发布而且十年 × 数百站点 × 逐小时的数据量对个人项目来说太大日粒度已经足够回答前面那三个问题。如果你要研究臭氧的日内变化规律再回头取小时数据就好没必要一开始就把摊子铺大。字段信息大致如下站点编码、站点名称、所属城市日期AQI、首要污染物PM2.5、PM10、SO2、NO2、CO、O3 各自浓度需要注意不同年份、不同批次的文件字段名和格式会有细微差异。有的叫quality有的叫aqi有的时间字段是字符串有的已经是时间戳。统一口径是预处理阶段最琐碎、也最不能省的工作。2.2 清洗流程缺失值、异常值和口径统一拿到手的数据绝对不是干净的实测下来主要有三类脏数据。第一类是站点变更问题。十年里有些监测站点被撤销、有些是新增的、还有些站点编码在某一年的文件中突然换了编号。如果不处理直接按站点编码聚合会得到一条断裂的时间序列看起来像空气质量突变其实是站点换了马甲。我的处理策略是先用站点名称和经纬度做去重和匹配把相同物理位置的站点合并为一个实体再按城市维度聚合分析不去纠结单个站点的连续性。城市层面数据更稳定因为单个站点的变动不会对城市均值造成跳跃式影响。第二类是缺失值。日尺度数据里某些站点的某些污染物字段会整体缺测特别是 O3 在早几年的数据里缺得比较多。对于缺失率低于 20% 的字段我用前后均值 季节性修正的方式补齐也就是取前后 7 天的均值再乘以该月份的历史均值比例做一个缩放。缺失率超过 20% 的字段直接放弃因为硬补出来的数据反而是噪音。第三类是异常值。AQI 体系里有一个特殊设计当浓度超过一定阈值后AQI 最高只能报到 500。这就导致真实情况可能是 AQI 800但记录里只有一个 500。这类封顶值不能当普通数值处理否则年均值会被低估。另外数据里偶尔会出现负浓度、AQI 大于 500 但浓度字段异常等明显录入错误直接剔除并标记为缺失走缺失值处理流程。还有一个特别容易被忽略的坑统计口径。PM2.5 用的是 24 小时平均浓度O3 用的是日最大 8 小时滑动平均CO 用的是 24 小时平均但不同批次文件里单位可能不一样——有的给微克/立方米有的给毫克/立方米。合并之前必须统一单位否则你会看到 CO 的数值比其他污染物大几百倍然后误判成CO 污染最严重。提示清洗完成后请务必做一次眼睛校验。随机抽出几个城市、几个日期去官方平台上人工核对一下数值。这一步不耗时但能提前拦截掉大批你根本想不到的脏数据问题。3. 核心分析方法与可视化实现3.1 时间序列趋势分解把十年数据拆成三层看拿到干净的日度数据后第一件事不是直接画年均值曲线而是做时间序列分解。我用的方法是 STLSeasonal-Trend decomposition using LoessPython 里 statsmodels 库直接支持。选择 STL 而不是简单移动平均是因为它能把时间序列拆成三部分长期趋势、季节性波动、残差项。趋势项回答空气质量在变好还是在变坏季节项回答一年内什么时候污染重残差项回答有没有什么突发事件打破了常规规律比如某年某月突然出现大面积异常高值。代码非常简单import pandas as pd from statsmodels.tsa.seasonal import STL # df 为按城市聚合后的日度数据包含 date 和 aqi 两列 df[date] pd.to_datetime(df[date]) df df.set_index(date).sort_index() # 周期设为 365因为我们要提取的是年内的季节规律 stl STL(df[aqi], period365, robustTrue) result stl.fit() # result.trend 是趋势项result.seasonal 是季节项result.resid 是残差 df[trend] result.trend df[seasonal] result.seasonal df[resid] result.resid三个参数值得展开说。period365表示以一年为周期做季节分解这个选择符合空气质量以年为周期的物理事实。robustTrue是关键它让分解过程对极端值不那么敏感——否则某几天爆表的 AQI 会把整条趋势线往上拉掩盖真实走势。STL默认的 Loess 平滑窗口不需要手动调但如果你发现趋势线太曲折可以适当增大seasonal的平滑参数。分解完之后画一张三行图上面是原始序列中间是趋势下面是季节性。这张图的信息量远超一张堆积折线图。我最直观的感受是趋势项不是一条均匀下降的直线而是前段下降快、中段趋缓、后段又有新变化的阶梯状曲线季节项则非常规律地呈现冬季高、夏季低的形态但不同城市、不同污染物的季节性相位差得很远。3.2 空间分布与城市画像不要只画一张全国均值图全国均值曲线有一个大问题它会抹平巨大的地域差异。华北某重点城市的空气质量变化轨迹和南方沿海城市的轨迹几乎不重叠硬要平均成一条线意义不大。所以我在全国视角之外单独挑了三类典型城市做对比华北某重点城市、南方某沿海城市、中西部某工业城市。用折线图对比这三座城市的年度 AQI 趋势能看到完全不同的故事线。华北某重点城市十年里 AQI 年均值下降幅度明显尤其是冬季改善最显著南方某沿海城市本来基数就好下降空间有限但这两年夏季 AQI 反而在抬升中西部某工业城市的下降拐点来得比前两者晚后程发力但波动大。这就是城市画像的雏形。我按主要污染物画像 季节性特征把所有城市粗略分成了三类冬季燃烧型污染物以 PM2.5 为主季节性强烈冬季峰值极高典型如华北、东北部分城市。夏季光化学型首要污染物在夏季频繁切换为 O3AQI 的季节曲线呈现双峰一个在夏季、一个在冬季。全年扬尘型PM10 贡献长期偏高季节差异不如前两类明显多见于西北、内陆地市。分类不需要搞复杂模型用首要污染物出现频次统计 AQI 季节曲线形态就能做。这步分类的价值在于后面所有对比分析都要在同类城市内部进行否则就是拿苹果比橘子。3.3 极端污染事件统计衡量健康风险的关键指标前面说过年均值会骗人。所以第三个分析模块专门聚焦极端事件把 AQI 大于 200 的天数定义为重度及以上污染日统计每个城市每年的这类天数同时记录每次极端过程的持续时长和发生时段。具体实现上用 pandas 的分组统计就能完成# 构造一个标记列 df[severe] df[aqi] 200 # 按城市、年份统计重度污染天数 df[year] df[date].dt.year city_year df.groupby([city, year])[severe].sum().unstack() # 查找连续极端过程用分组技巧标记连续段 df[group] (df[severe] ! df[severe].shift()).cumsum() episodes df[df[severe]].groupby([city, group]).agg( start(date, min), end(date, max), duration(date, count) )统计结果非常有意思。大多数北方城市的重度污染天数从早期的每年几十天压缩到近年的个位数这是趋势分析里看不到的细节。同时还发现极端事件的发生时段越来越集中从曾经的整个采暖季频繁爆发收敛为偶发于 12 月到 1 月的少数几次过程。这个规律如果只盯着年均值是完全没有办法发现的。极端事件的统计还有一个实用价值它能直接回答空气质量是不是变好了这种日常提问。普通人感受到的不是年均值的波动而是冬天还有几天要戴口罩。用极端天数这个指标来衡量往往比官方年均值更有说服力。4. 实操链路与关键代码解读4.1 数据读取与日期规整这一步的目标是把散落在多个 CSV 文件里的数据统一加载到一个 DataFrame 里。因为文件多我用pathlib遍历目录按年份循环读取。读取时注意两个细节一是用dtype参数把城市和站点编码指定为字符串防止前导零丢失二是用parse_dates直接解析日期列避免之后二次转换。import pandas as pd from pathlib import Path data_dir Path(./aqi_data) frames [] for f in sorted(data_dir.glob(*.csv)): df_tmp pd.read_csv( f, encodingutf-8, dtype{city: string, station_code: string}, parse_dates[date] ) frames.append(df_tmp) df pd.concat(frames, ignore_indexTrue) df df.drop_duplicates(subset[station_code, date], keeplast)drop_duplicates这行是实测中加上去的。不同年份文件的边界可能存在重复记录同样的站点、同一天可能出现两次直接取最后一次避免后续聚合时把同一天的数值算两次。4.2 城市级聚合与数据透视站点级别的数据太碎而且站点变动会影响连续性所以我统一按城市聚合。日均值用简单平均但更稳健的做法是先按站点平均、再按城市平均避免站点数量不均衡导致某些城市被站点多的区域主导。# 先按站点、日期聚合此时已是日度数据主要做字段检查 station_day df.groupby([city, station_code, date]).agg( aqi(aqi, mean), pm25(pm25, mean), pm10(pm10, mean), o3(o3, mean), ).reset_index() # 再按城市、日期聚合 city_day station_day.groupby([city, date]).agg( aqi(aqi, mean), pm25(pm25, mean), pm10(pm10, mean), o3(o3, mean), ).reset_index() # 生成城市 × 年份的透视表 city_year city_day.pivot_table( indexcity, columnscity_day[date].dt.year, valuesaqi, aggfuncmean )透视表的好处是快速查看哪些城市哪一年明显偏高但如果要看完整时间序列透视表反而碍事。所以这个表我只用来做快速扫描正式绘图和分析还是用长表city_day。4.3 可视化从折线图到热力图可视化是整个项目里最出效果、也最坑的一步。先说坑matplotlib 默认不显示中文坐标轴上的城市名全变成方框。解决办法很简单指定一个支持中文的字体即可import matplotlib.pyplot as plt plt.rcParams[font.sans-serif] [SimHei] # 或用你系统里的中文字体 plt.rcParams[axes.unicode_minus] False核心可视化我做了四类图。第一类是趋势图用全国均值加置信带的方式展示整体走势置信带用seaborn.lineplot的ci参数即可。第二类是城市对比图用双 y 轴同时画 AQI 和 PM2.5观察两者是否同步变化。第三类是季节箱型图把所有年份的月度 AQI 分布画成箱型能直接看出冬季中位数和离散度是否逐年收窄——这是判断极端污染少了没最直观的图。第四类是热力图横轴是月份、纵轴是年份、颜色代表 AQI一张图就能看尽十年的季节规律演变。import seaborn as sns import matplotlib.pyplot as plt # 月份-年份热力图 city_data city_day[city_day[city] 华北某重点城市].copy() city_data[month] city_data[date].dt.month city_data[year] city_data[date].dt.year pivot_monthly city_data.pivot_table( indexyear, columnsmonth, valuesaqi, aggfuncmean ) fig, ax plt.subplots(figsize(12, 8)) sns.heatmap(pivot_monthly, cmapYlOrRd, axax) plt.title(华北某重点城市 AQI 月度均值热力图) plt.xlabel(月份) plt.ylabel(年份) plt.show()这张热力图是整份分析里信息密度最高的一张图。我能从中快速读到十年前的颜色从 11 月一直深到次年 3 月近年深色区域明显收缩但夏季颜色并没有变浅太多——氧气化型污染在悄悄抬头。单看年均值曲线你永远发现不了这个规律。4.4 季节性分解的进阶用法前面 STL 分解用的是 AQI 序列但更细的分析需要把每个污染物单独做一遍分解。比如 SO2 的趋势下降反映的是能源结构调整O3 的季节项峰值在夏季且逐年抬升反映的是光化学反应增强。这比只看 AQI 综合指标准确得多。把分解函数封装一下可以批量跑def decompose_trend(series, period365): stl STL(series, periodperiod, robustTrue) result stl.fit() return result.trend # 以城市和污染物为分组批量计算趋势项 trend_result {} for city in city_day[city].unique(): tmp city_day[city_day[city] city].set_index(date).sort_index() trend_result[city] { col: decompose_trend(tmp[col]) for col in [aqi, pm25, pm10, o3] }这一步输出的是趋势数据之后可以计算十年末相比十年初的降幅百分比用作城市排名的核心指标。我实测下来用趋势项的降幅而非原始年均值的降幅做排名更能反映真实改善水平——因为趋势项滤掉了气象波动带来的逐年噪声。5. 常见问题与排查技巧实录5.1 数据断裂站点编码变动导致的假突降做聚合时最隐蔽的坑就是站点编码变更。我遇到过一个城市的数据2014 到 2019 年用的是老编码2020 年开始换成新编码。如果直接按站点编码分组会看到两个完全不同的站点各自身下的数据都不连续。更麻烦的是这种断裂不会报错你会把站点换了误判成空气质量突然变化来写结论。排查技巧做站点级数据完整性检查逐年统计每个站点全年的非空天数。正常站点一年应该有 350 天以上数据。如果一个站点在某一年突然只剩半年数据、或者编码消失、次年又出现新编码就去查站点名称和经纬度确认是不是同一个物理站点换了个 ID。确认后做站点合并或者直接放弃站点级分析、只保留城市级聚合——后者更省心。5.2 封顶值对统计量的污染AQI 封顶 500 的问题在做极端事件统计时不明显但做年均值、趋势分解时会被放大。原因很简单假设某城市冬季某天真实 AQI 是 700但因为封顶只记录了 500。这个被压低的数值参与年均值计算导致该城市那一年年均值偏低几十甚至上百。如果封顶只发生在一年里的几天影响有限如果某城市连续多年冬季频繁爆表累加起来的低估效应足以改变趋势判断。我的处理方式统计阶段增加一个是否封顶标记列单独计算封顶天数指标把它和年均值并列呈现。这样做既不需要修改原始数据又不会让封顶值污染均值分析。5.3 热力图颜色映射的假可视化热力图看似简单但默认的vmin和vmax会自动按数据范围映射。当数据里有极端高值时比如某年冬季 AQI 达到 500颜色映射会自动拉伸导致大多数年份的中低值全部挤在浅色区域颜色对比失真。解决办法是手动指定色标范围sns.heatmap(pivot_monthly, cmapYlOrRd, axax, vmin0, vmax250) # 手动设定范围防止极端值拉伸这个细节值得多说一句可视化不是为了把图做得好看而是为了不误导读者。自动色标给出的图很可能让你得出这些年变化不大的错误结论。5.4 一个经典分析误区只看年均值整个项目做下来我认为最值得强调的误区是年均值是分析起点不是终点。年均 AQI 从 100 降到 80听起来改善很大但如果极端污染天数的下降曲线在近几年趋于平缓说明常规日子变好了但极端情况改善不大。反过来有的城市年均值降幅不大但极端事件显著减少对居民健康防护的意义反而更大。所以我建议在看任何城市时同时看三个数年均值、冬季均值、重污染天数。三个数一起看才能勾勒出完整画面。单拎一个数出来说事都容易带偏结论。提示如果你做的是城市间对比记得先用人口密度、地形条件、气候背景做粗分类再在同类城市里比较。否则把沿海城市和内陆工业城市放在同一张排名表里除了引发争论得不到任何有价值的信息。6. 事后复盘与扩展方向这次分析项目做完之后有几个数据之外的体会想分享。第一环境数据是典型的杂而全数据。它能支撑的结论很多但每一层结论都要付出相应的清洗代价。城市级聚合分析大概能覆盖 80% 的常见问题但如果你把尺度下探到站点级别所有脏数据问题都会放大十倍。对于个人项目来说宁可做粗一点、做对一点也不要一开始就追求精细化。第二时间序列分析的方法论是通用的。STL 分解、极端事件统计、季节箱型图这些方法不只适用于空气质量数据用在气温、降水、用电负荷、车流量等任何有时间戳的数据上都能得出有价值的结论。学会这套流程等于多了一个通用的数据复盘工具。第三关于数据源的另一个思考。这次用的是官方公开数据数据质量相对可控。如果换一个数据源质量参差不齐建议建立一套抽样人工校验的流程从数据里随机抽 20 个点逐条去原始页面核对。我这次就通过抽样发现了一个批次文件的臭氧浓度字段整体偏移约 10% 的问题这种问题靠代码很难自动发现。最后说一个后续可以做的方向把气象数据叠加进来做关联分析。AQI 的波动很大程度上受气象条件驱动温度、风速、湿度、边界层高度都影响污染物扩散。单纯看 AQI 时间序列你只能描述发生了什么叠加气象数据后才能初步讨论为什么发生。再往后走还可以用机器学习模型做污染事件的归因分析区分排放减少的贡献和气象变好的贡献。但那是另一个项目的事了。如果这个十年回顾分析能给你带来一条可复用的数据处理思路或者让你对空气质量到底变好没有这个问题有了新的看法那这次折腾数据就没有白费。