最近在整理空气污染相关的数据拿到手的就是 CnOpenData 的站点空气质量表。这种表在环境数据分析里太常见了一行是一个监测站点在某个时间点的记录后面跟着AQI和六项常规污染物浓度。字段看着不复杂但真正用起来会发现从原始数据到能直接分析的标准表中间要过的关卡比想象中多得多。这篇文章我就以这份站点空气质量表为中心把完整处理流程摊开讲一遍从字段认知、数据清洗、空间分析到常见问题排查每个环节都给出可以直接参考的代码和踩坑提醒。如果你正准备拿这类数据写论文、做课程设计或者给业务方产出一份空气质量分析报告这篇文章应该能让你少走不少弯路。1. 站点空气质量表的数据结构认知拿到表先别急着跑模型1.1 一张典型的站点空气质量表长什么样站点空气质量表的核心逻辑一句话就能讲完每一个观测站点在每一个观测时刻都对应一行记录。我接触过日度版本也有小时版本粒度不同但结构大同小异。日度版本大概是全国上千个站点乘365天几十万行小时版本直接上百万甚至上千万行。后面的分析方式会有差别但底层数据结构是同一套。这类表的典型字段如下面这张表所示字段名类型说明station_code字符串站点编码全局唯一标识station_name字符串站点名称可能存在重复province / city字符串站点所在的省份和城市longitude / latitude浮点数站点经纬度date / hour时间观测时间部分表是单独两列AQI整数或浮点数空气质量指数PM2_5 / PM10 / SO2 / NO2 / O3 / CO浮点数六项污染物浓度level字符串或整数空气质量等级为什么要强调长表这个概念因为清洗和聚合阶段长表用 groupby 非常顺手。我之前见过有人拿到表之后先 pivot把每个站点变成一列结果后面新增站点时脚本全崩。正确习惯是保持长表需要宽表时再临时转换。1.2 字段单位和小细节弄错一个全盘皆输字段名都很直白真正需要小心的是单位。PM2.5、PM10、SO2、NO2、O3 的单位通常是微克每立方米而 CO 则是毫克每立方米两个量纲差了1000倍。做相关性分析之前最好先用 describe 确认每一列最大值。如果 CO 平均浓度到了几百基本就是单位写错了要么源数据有误要么读表时列错位。另一个容易忽略的是 AQI。AQI 不是六项污染物浓度的简单平均而是先算出每个污染物的 IAQI 分指数然后取最大值。空气质量等级对应 AQI 区间比如优、良、轻度污染等。如果表里的 AQI 值和用浓度算出来的不一致别急着怀疑自己源表可能混入了错误行这部分我在第五章会给出一个校验函数。时间字段也值得多看一眼。有的表用 2023-01-01 00:00:00有的拆成 date 和 hour 两列还有可能混入 2023/1/1 0:00 这种写法。读进来之后先看 info 里时间这一列是不是 datetime 类型如果不是先统一格式再继续。经纬度的话大部分公开数据用 WGS84 坐标系但也有少量来源会做偏移后面画地图时如果发现站点整体偏移通常就是坐标系问题。1.3 拿到数据后的第一步先给数据做个体检我的习惯是任何站点空气质量表到手先不做统计建模花一分钟做数据体检。用 pandas 读进来打一个 info 和 describe比什么都有用。import pandas as pd df pd.read_csv(station_air_quality.csv) print(df.shape) print(df.columns.tolist()) print(df.head(10)) print(df.info()) print(df.describe(includeall))shape 看行数和列数columns 看列名是否有重复head 快速看前几行info 能看出每列有没有缺失、类型是否正确describe 会暴露所有数值列的分布。这一步最关键的检查点有三个是否有负浓度、AQI 最大值是否离谱、date 列是不是被读成了 object。我见过有人在数据体检都没做的情况下直接 groupby 画图后来发现站点编号有重复画出来的曲线断断续续回去排查花了两个小时。还要顺手看一眼站点数、时间范围print(df[station_code].nunique()) print(df[date].min(), df[date].max())这一步做完你对这张表的健康状况就有了基本判断。如果发现 date 还是字符串别犹豫继续往下清洗。2. 数据清洗从原始表到标准表重点是别把数据洗坏2.1 统一时间格式并建立索引如果体检发现时间列是 object 类型第一步就是把它转成 DatetimeIndex。这里有个小技巧尽量指定 format转换速度快也能避免自动推断带来的意外结果。df[date] pd.to_datetime(df[date], format%Y-%m-%d %H:%M:%S, errorscoerce)如果时间来源是分开的日期和小时两列先把两列拼起来再转df[datetime] pd.to_datetime( df[date].astype(str) df[hour].astype(str), errorscoerce ) df df.dropna(subset[datetime]) df df.sort_values([station_code, datetime]).set_index(datetime)设置索引之后立刻检查索引是否唯一。很多分析操作比如 resample、merge_asof 都要求索引干净重复索引会报错或者结果错乱。print(df.index.is_unique)如果发现重复索引保留第一条即可df df[~df.index.duplicated(keepfirst)]这个操作看似简单但能避免后面大量报错。有一次我处理小时级数据直接 set_index 后做 resample结果报 “cannot reindex from a duplicate axis”排查半天才知道是数据源里同一个站点同一小时出现了两条记录。2.2 缺失值处理能补才补不能补就留着缺失值处理没有万能公式我的原则是能通过重算得到的不插值短连续缺失可以插值长段缺失宁可留着。先把各列缺失比例跑出来missing df.isna().mean().sort_values(ascendingFalse) print(missing)如果 PM2.5 缺失率超过20%这个站点基本不值得补如果只是零星几个小时的缺失线性插值就够了。注意插值要按站点分组进行不同站点的浓度水平差很多不能把全国所有站点放在一起插值否则会把高值站点的浓度拉低。cols [PM2_5, PM10, SO2, NO2, O3, CO] for col in cols: df[col] df.groupby(station_code, group_keysFalse)[col].apply( lambda s: s.interpolate(limit6) if s.isna().mean() 0.2 else s )有一点必须提醒插值后的数据带有平滑效果会压低真实波动。用来做图表展示没问题但如果后续要做趋势检验或者浓度超标统计建议保留原始缺失只在绘图时临时补全。这是很多刚接触数据分析的同学容易忽略的点。2.3 异常值和重复值比缺失值更隐蔽异常值常见两种浓度出现负值以及浓度高到完全不可能。负值直接过滤掉不用心疼。离谱高值要结合站点和季节判断比如 PM10 超过 5000 微克每立方米在普通城市站点几乎不可能出现多半是单位错误或者传感器故障。df df[df[cols].ge(0).all(axis1)]重复值方面除了整行重复还有“同一站点同一时间出现两条不同浓度”的更隐蔽情况。前者用 drop_duplicates 处理后者要先用 groupby 判断是否有重复dup_mask df.reset_index().duplicated(subset[station_code, datetime], keepFalse) print(dup_mask.sum())真实情况里我遇到过一种很坑的场景两个表拼接时CO 浓度一列是 mg/m³另一列却是 μg/m³拼起来之后均值和相关性全乱。排查了半个小时才通过 max 值发现量纲不一致。所以体检里的 describe每个字段的最大值一定要亲自看一眼。3. 站点与空间维度把点位信息变成可分析的资产3.1 从长表里抽出站点维度表站点空气质量表是典型的长表站点信息会在每一行重复出现。与其每次分析都去筛选站点字段不如先抽一张站点维度表后面城市匹配、画分布图、做地图联动都直接查表。station_info df[[station_code, station_name, province, city, longitude, latitude]] station_info station_info.drop_duplicates(station_code).reset_index(dropTrue)站点编号一定是唯一标识不要用站点名称中文名可能重复还可能包含特殊字符。抽完维度表后顺手检查经纬度缺失值和站点数量。如果同一站点编码对应多个经纬度说明源数据里存在重复且信息不一致需要按时间最新的一条保留。3.2 把监测点画到地图上站点分布可视化看起来只是画个图但在实际工作中能帮你快速判断数据覆盖范围。我常用 plotly一张图搞定交互import plotly.express as px fig px.scatter_mapbox( station_info, latlatitude, lonlongitude, colorcity, zoom3, mapbox_styleopen-street-map, ) fig.show()不需要花哨的话用 matplotlib 的静态散点图也行。全国范围直接用经纬度画散点能看出站点集中在哪些区域。要特别注意如果叠加在线地图上发现站点整体偏移几百米大概率是坐标系问题。公开站表多为 WGS84但有些间接来源会变成 GCJ-02需要先转换再叠加。3.3 站点数量对城市均值的代表性影响站点空气质量表最常用的统计口径是城市日均浓度。但如果一个城市只有一两个站点它的日均值波动会很大和周边城市做横向对比时容易出现“异常高”或“异常低”。处理前先统计每个城市的站点数量city_station_count station_info.groupby(city)[station_code].nunique().sort_values() print(city_station_count)如果某个城市站点数太少在报告里最好备注一下代表性有限。做政策评估或归因分析时站点密度会影响结论的可靠性不能只盯着数值。这个细节很容易被带过但实际影响很大。4. 实战从站点空气质量表提炼时空规律4.1 计算城市日均浓度分组聚合的细节站点维度的日均浓度直接按城市和日期分组代码非常短df[date_only] df.index.date city_daily df.groupby([city, date_only])[PM2_5].mean().reset_index() city_daily.columns [city, date, pm25_daily]这里有个问题简单平均默认把每个站点等权重处理。如果站点在城市内分布不均匀简单平均会偏向站点密集的片区。标准做法里可以考虑按区县人口加权但多数公开研究直接用算术平均。我的建议是先跑简单均值同时统计每日参与计算的站点数如果某天站点数骤减那天的均值不能作为可靠结果。4.2 污染物相关性矩阵一张热力图快速看同源性把污染物之间的相关性算出来是做来源分析前很常用的动作。做法是先选一个时间范围去掉缺失较多的站点再算相关系数矩阵pivot df.pivot_table( indexdf.index, columnsstation_code, valuesPM2_5 ) corr pivot.corr(methodpearson)PM2.5 和 PM10 通常高度相关因为它们来自燃烧、扬尘等共同源NO2 和 PM2.5 相关性高往往说明交通源影响明显。相关分析不能证明因果但能帮你发现异常如果两个应该高相关的污染物相关系数很低要么是站点本身数据质量差要么是你在处理异常值时把方向弄反了。做热力图的时候可以用 seaborn 一行出图import seaborn as sns import matplotlib.pyplot as plt sns.heatmap(corr, annotTrue, cmapRdBu_r, center0) plt.show()4.3 时间序列趋势滑动平均是最快的平滑方式站点级浓度噪声很大日度数据画出来全是毛刺直接看趋势很难。我的习惯是先按城市聚合再做7日滑动平均city_daily city_daily.set_index(date).sort_index() smooth city_daily[pm25_daily].rolling(7, centerTrue, min_periods3).mean()如果数据有完整一年以上可以用 statsmodels 的季节分解把趋势、季节、残差分开from statsmodels.tsa.seasonal import seasonal_decompose result seasonal_decompose( city_daily[pm25_daily].dropna(), modeladditive, period365 )注意 seasonal_decompose 要求时间索引等间隔插值后的数据更适合做这种分解。但分解结果要谨慎解读特别是对短序列趋势部分很容易被边界效应影响。5. 常见问题与排查技巧实录5.1 先看速查表这一节我直接把平时最常遇到的问题整理成了一张表适合遇到问题对照着查现象可能原因处理办法中文乱码文件编码不匹配用 encodinggbk 或 gb18030 重新读取时间解析后大量 NaT源表时间格式混用先看 unique 值统一格式后转 datetimeresample 报 duplicate axis索引有重复先用 duplicated 查重再保留唯一值站点画在地图上整体偏移坐标系不是 WGS84确认源数据坐标系做坐标转换AQI 与浓度对不上源表脏数据用 IAQI 反算标记异常行读大文件内存不足默认读取方式占用太高指定 dtype、category、chunksize 分段读入下面挑三个我反复踩过的坑展开讲。5.2 AQI不靠谱用IAQI反算来校验站点空气质量表里AQI 和浓度之间理应有严格的计算关系。如果 AQI 和污染物浓度对不上很有可能表里某些行是历史旧版本或被人为改过。校验思路就是根据浓度反算每个污染物的 IAQI取最大值再和表里的 AQI 比较。以 PM2.5 为例日平均浓度和 IAQI 存在分段线性关系可以写一个简化版函数def iaqi_pm25(v): breaks [ (0, 35, 0, 50), (35, 75, 50, 100), (75, 115, 100, 150), (115, 150, 150, 200), (150, 250, 200, 300), (250, 350, 300, 400), (350, 500, 400, 500), ] for lo, hi, ilo, ihi in breaks: if lo v hi: return (v - lo) / (hi - lo) * (ihi - ilo) ilo return None实际项目中完整校验需要六项污染物各自的浓度区间表写法都一样只是区间不同。反算之后把计算值和表里的 AQI 做差超过1就标出来。我曾经用这个办法在几百万行的表里找出好几百条浓度正确但 AQI 写错的脏数据修完再做分析结果干净很多。5.3 大文件读取与内存优化小时级的站点空气质量表很容易上千万行直接 read_csv 不是不行只是内存占用很高。先指定 dtype把站点编码、城市这类重复值很多的列设为 category时间列用 parse_dates 解析内存能省下一大截df pd.read_csv( station_air_quality_hour.csv, dtype{ station_code: category, city: category, }, parse_dates[datetime], )如果文件实在太大就分段读取每段只保留你要的聚合结果最后再合并chunks [] for chunk in pd.read_csv(station_air_quality_hour.csv, chunksize500000): chunks.append( chunk.groupby(city)[PM2_5].mean() ) city_mean pd.concat(chunks).groupby(level0).mean()这一步放在清洗之前也不会影响后续流程只是读表方式变了。我平时处理这类长表基本都会先评估数据规模和内存否则后面 pivot_table 或者 merge 很容易把进程卡死。6. 扩展玩法把站点空气质量表用得更深6.1 与气象数据打通分析污染和天气的关系站点表本身只有浓度数据很多研究还要结合气象条件。最简单的做法是用 merge_asof让每个站点时间点匹配到最近一小时内的气象观测tolerance 设成 1 小时即可。注意匹配前两边的数据都要按时间排序并按站点编码分组。weather pd.read_csv(weather.csv, parse_dates[datetime]) merged pd.merge_asof( df.sort_index(), weather.sort_index(), ondatetime, bystation_code, tolerancepd.Timedelta(1h), )这样匹配后就能分析风速低时 PM2.5 浓度是否更高、湿度对二次污染物生成的影响等。做这类分析时相关性不代表因果还要考虑站点周边排放源、地形等因素。6.2 我处理这类表的一些个人习惯最后分享几个我自己长期用下来的习惯没什么高深理论但能省不少时间。第一站点信息永远单独存一张维度表不混在明细数据里第二所有清洗步骤都写成函数输入输出明确方便同一份数据反复更新第三任何聚合结果都要保留“参与计算的记录数”否则无法判断数据质量波动对结论的影响。站点空气质量表看起来简单真正做好清洗、空间匹配、时间聚合每一步都需要对数据的理解。这篇文章里写的都是我在实际处理过程中反复用到的操作你可以把自己的数据套进去试一遍。如果遇到上面表格里没有的问题也建议顺着这个思路继续排查先查明字段含义再检查数据粒度最后再怀疑数据本身。