简介本资源是一份面向计算机及相关专业如人工智能、自动化、电子信息等在校学生与初学者的完整数据分析实战项目聚焦北京市12个监测点空气质量数据的采集、清洗、可视化与建模分析可直接用于课程设计、大作业、毕业设计或能力进阶训练。压缩包共1565个文件含72个CSV原始数据集如Wanshouxigong、Aotizhongxin等站点数据、5个核心Python分析脚本、1470个HTML格式的交互式分析报告含图表与结果解读、8张关键分析截图及详细README.md说明文档整体容量343.85MB结构清晰、开箱即用。已有192人学习下载项目源自作者高分毕设答辩平均分96分所有代码均经实测运行成功附带远程答疑支持。读者可获得从数据获取、Pandas/Seaborn/Matplotlib全流程处理、多维度时空趋势分析到报告生成的完整闭环方案兼具教学性、可复现性与工程参考价值。1. 北京市12个监测点空气数据怎么“看懂”不是画几个折线图就叫数据分析你拿到一份标着“北京市12个地点空气监测数据”的Excel或CSV文件里面密密麻麻堆着PM2.5、PM10、SO₂、NO₂、CO、O₃六项污染物浓度还有时间戳、站点名称、AQI指数——但打开后第一反应往往是这堆数字到底在说什么为什么朝阳门站的PM2.5凌晨总比中午高为什么海淀万柳站的O₃峰值总比通州运河西晚两小时为什么同一时间不同站点的CO值能差出3倍这些问题光靠pandas.read_csv()plt.plot()根本答不上来。这不是“Python画图作业”而是真实城市环境数据的时空归因分析它要求你把12个空间点当作一个动态网络把每小时/每日的数值序列当作生理脉搏用统计检验锚定异常时段用相关性矩阵识别污染传输路径用滑动窗口捕捉早晚高峰特征最后用地理坐标叠加热力图验证“风向主导扩散”这个常识是否成立。本项目不教print(Hello World)它面向的是已能写函数、会调scipy.stats、知道groupby().agg()和resample()区别的一线数据处理者——你要的不是“能跑通”而是“跑通之后敢拿结论去跟环保部门对线”。2. 数据加载与结构清洗从原始CSV到可计算的时空立方体原始空气监测数据常以“站点×时间×指标”三维形式存在但实际交付的CSV往往埋着三类坑字段名中英文混杂如PM2.5(μg/m³)、时间列格式混乱2023-01-01 00:00vs2023/01/01 00:00:00、缺失值标记不统一、NULL、-999、—, 甚至空格。直接pd.read_csv()会把整列转成object类型后续数值计算全崩。必须分步强校验。2.1 用正则预筛字段名并标准化列索引import pandas as pd import re def clean_column_names(df): 清洗列名移除单位、空格、括号转小写下划线分隔 new_cols [] for col in df.columns: # 移除括号及内容如 (μg/m³)、多余空格、特殊符号 cleaned re.sub(r\([^)]*\)|\s|[^\w], _, col) # 合并连续下划线首尾去下划线 cleaned re.sub(r_, _, cleaned).strip(_) # 转小写 cleaned cleaned.lower() new_cols.append(cleaned) df.columns new_cols return df # 示例原始列名可能为 [StaName, Time, PM2.5(μg/m³), NO2(μg/m³), AQI] # 清洗后变为 [staname, time, pm25, no2, aqi]提示此步必须在read_csv()后立即执行。若跳过后续用df[PM2.5]会报KeyError——因为实际列名是PM2.5(μg/m³)而Pandas默认不支持括号内含斜杠的字符串索引。2.2 时间列解析强制指定format避免自动推断翻车北京监测数据时间戳常见两种格式%Y-%m-%d %H:%M标准和%Y/%m/%d %H:%M:%S部分子站。若用parse_dates[time]且不指定date_parserPandas会尝试自动推断结果是前1000行成功第1001行因某条记录多了一个秒数如2023-01-01 00:00:00导致整列转为object后续resample(D)直接报错TypeError: Only valid with DatetimeIndex, TimedeltaIndex or PeriodIndex。正确做法是显式指定最严格的format并启用错误处理from datetime import datetime def parse_beijing_time(x): 北京空气数据时间解析器覆盖两种主流格式 if isinstance(x, str): x x.strip() # 优先匹配带秒的完整格式 for fmt in [%Y-%m-%d %H:%M:%S, %Y/%m/%d %H:%M:%S, %Y-%m-%d %H:%M, %Y/%m/%d %H:%M]: try: return datetime.strptime(x, fmt) except ValueError: continue return pd.NaT # 解析失败返回NaT便于后续dropna # 加载时即解析 df pd.read_csv(beijing_air_12sites.csv, date_parserparse_beijing_time, parse_dates[time], keep_date_colTrue) # 保留原始列用于debug2.3 缺失值与异常值的物理意义判别空气监测设备故障时常输出-999或999999作为占位符非真实负值。若直接df.fillna(0)或df.dropna()会抹掉真实低浓度时段如凌晨O₃常低于10μg/m³更致命的是-999被当数值参与mean()计算会拉低全站日均值达20%以上。必须按污染物物理阈值站点历史分布双维度过滤# 定义各污染物合理范围依据《环境空气质量标准》GB3095-2012 VALID_RANGES { pm25: (0, 1000), # μg/m³实测极值800 pm10: (0, 2000), # μg/m³ so2: (0, 500), # μg/m³ no2: (0, 1000), # μg/m³ co: (0, 50000), # μg/m³ → 换算为mg/m³需/1000 o3: (0, 600), # μg/m³ aqi: (0, 1000) # 无量纲 } def flag_invalid_values(df, valid_rangesVALID_RANGES): 标记超出物理合理范围的值为NaN df_clean df.copy() for col, (low, high) in valid_ranges.items(): if col in df_clean.columns: # 对每列单独mask超出范围即设为NaN mask (df_clean[col] low) | (df_clean[col] high) df_clean.loc[mask, col] np.nan return df_clean df_clean flag_invalid_values(df) # 此时再做 df_clean.groupby(staname)[pm25].mean() 才可信3. 时空聚合与特征工程把小时数据变成可解释的“城市呼吸节律”原始数据是小时级但人对空气的感知是日尺度晨练是否咳嗽、周尺度周末是否更闷、季节尺度供暖季PM2.5必然飙升。直接对24小时数据做mean()会丢失早晚高峰信息简单取max()又忽略持续暴露风险。必须设计符合环境健康逻辑的聚合策略。3.1 日尺度区分“暴露强度”与“暴露时长”WHO指出PM2.5的健康风险与24小时均值和峰值浓度持续时间双重相关。因此不能只算一个resample(D).mean()def daily_air_features(df, time_coltime, site_colstaname): 生成日尺度特征包含均值、峰值、超标小时数、早晚高峰差值 返回MultiIndex DataFrame (date, site) - features # 确保time_col为datetime且设为index df_indexed df.set_index(time_col).sort_index() # 按站点日期分组 daily_groups df_indexed.groupby([pd.Grouper(freqD), site_col]) # 定义各污染物日特征 agg_dict {} for pol in [pm25, pm10, no2, o3]: if pol in df_indexed.columns: # 均值反映整体暴露水平 agg_dict[f{pol}_mean] (pol, mean) # 日最大值反映瞬时风险 agg_dict[f{pol}_max] (pol, max) # 超标小时数以PM25日均值35μg/m³为国标二级限值 if pol pm25: agg_dict[pm25_over_35h] ( pol, lambda x: (x 35).sum() ) # 早晚高峰差值早7-9点均值 - 晚17-19点均值表征交通源贡献 agg_dict[f{pol}_am_pm_diff] ( pol, lambda x: x.between_time(07:00, 09:00).mean() - x.between_time(17:00, 19:00).mean() ) # 执行聚合 daily_feats daily_groups.agg(agg_dict).reset_index() # 修复列名原为MultiIndex展平 daily_feats.columns [date, staname] [c[0] for c in agg_dict.keys()] return daily_feats daily_df daily_air_features(df_clean)3.2 周尺度捕捉“工作日效应”与“周末效应”北京机动车限行政策导致工作日NO₂显著高于周末而O₃因周末VOCs排放减少反而降低。需提取星期几并分组# 添加星期列Monday0, Sunday6 daily_df[weekday] daily_df[date].dt.weekday # 计算各站点工作日vs周末的NO2均值差异 workweek_no2 daily_df[daily_df[weekday] 5].groupby(staname)[no2_mean].mean() weekend_no2 daily_df[daily_df[weekday] 5].groupby(staname)[no2_mean].mean() # 差异率 (工作日均值 - 周末均值) / 周末均值 no2_workweek_ratio (workweek_no2 - weekend_no2) / weekend_no2 * 100 # 输出TOP3差异站点单位% print(no2_workweek_ratio.sort_values(ascendingFalse).head(3)) # 示例输出 # guomao 28.4 # xizhimen 25.1 # dongsi 22.7 # dtype: float64注意此处weekday必须用.dt.weekday而非.dt.dayofweek后者在旧版Pandas中行为不一致且必须在daily_df上计算——若在小时数据上算dt.weekday再groupby会因时区问题导致跨日错误如UTC8的23:00在北京是周一但若系统时区为UTC则算作周日。3.3 季节尺度用滚动窗口识别“供暖季突变点”北京供暖季11月15日-3月15日PM2.5浓度陡增但气象条件湿度、风速也同步变化。单纯对比11月vs10月均值会混淆“政策效应”与“自然效应”。应使用滑动t检验检测均值突变点from scipy import stats def detect_heating_season_shift(df, pollutantpm25, window30, alpha0.05): 检测供暖季开始前后PM2.5均值是否发生统计显著跃升 使用滑动窗口t检验窗口前半段 vs 后半段 # 按日期排序确保时序 df_sorted df.sort_values(date).set_index(date) # 只取该污染物有效数据 series df_sorted[pollutant].dropna() # 初始化结果列表 results [] for i in range(window, len(series)): # 窗口内前半段供暖前与后半段供暖中 pre_window series.iloc[i-window:i-window//2] post_window series.iloc[i-window//2:i] if len(pre_window) 10 and len(post_window) 10: # 最小样本量 t_stat, p_val stats.ttest_ind( pre_window, post_window, equal_varFalse, # 方差不齐用Welchs t-test nan_policyomit ) if p_val alpha: # 记录突变时间点窗口中心 center_date series.index[i - window//2] results.append({ date: center_date, t_stat: t_stat, p_value: p_val, pre_mean: pre_window.mean(), post_mean: post_window.mean(), abs_change: post_window.mean() - pre_window.mean() }) return pd.DataFrame(results) # 执行检测以PM25为例 shift_df detect_heating_season_shift(daily_df, pm25) # 查看最早显著突变点 if not shift_df.empty: first_shift shift_df.loc[shift_df[date].idxmin()] print(fPM2.5首次显著上升{first_shift[date].strftime(%Y-%m-%d)}, f增幅{first_shift[abs_change]:.1f}μg/m³ (p{first_shift[p_value]:.3f})) # 输出示例PM2.5首次显著上升2023-11-18, 增幅12.3μg/m³ (p0.002)4. 空间关联分析12个站点不是孤立点而是大气流动的传感器网络把12个站点当成12个独立样本做单变量统计是最大误区。北京地形西高东低冬季盛行西北风污染物从石景山→海淀→朝阳→通州输送。若忽略空间位置就无法解释“为何同样PM2.5日均值石景山站AQI120而通州站AQI180”——因为AQI计算含O₃而O₃在输送过程中光化学生成。4.1 地理坐标获取与空间权重矩阵构建原始数据通常只有站点名如dongsi需映射为经纬度。严禁手动查百度地图再复制粘贴——12个站点易出错且无法复现。应调用高德/腾讯地图API需申请key或使用公开地理编码库# 推荐使用geopy Nominatim免费但需遵守频率限制 from geopy.geocoders import Nominatim import time def get_beijing_stations_coords(site_names): 获取北京12个标准监测站点经纬度缓存机制防重复请求 # 北京市监测站点官方名称与标准简称映射关键 official_map { dongsi: 北京市东四环监测站, guomao: 北京市国贸监测站, xizhimen: 北京市西直门监测站, nongzhanguan: 北京市农展馆监测站, gucheng: 北京市古城监测站, huairou: 北京市怀柔监测站, shunyi: 北京市顺义监测站, daxing: 北京市大兴监测站, yizhuang: 北京市亦庄监测站, fangshan: 北京市房山监测站, pinggu: 北京市平谷监测站, mentougou: 北京市门头沟监测站 } geolocator Nominatim(user_agentbeijing_air_analysis) coords {} for short_name in site_names: full_name official_map.get(short_name, f{short_name}, Beijing) try: location geolocator.geocode(full_name , Beijing, China, timeout10) if location: coords[short_name] { lat: round(location.latitude, 5), lon: round(location.longitude, 5) } print(f✓ {short_name}: {coords[short_name]}) else: print(f✗ {short_name}: 未找到) except Exception as e: print(f⚠ {short_name} 请求失败: {e}) time.sleep(1) # 防封禁 return coords # 执行仅需运行一次结果存为coords.json供后续复用 # station_coords get_beijing_stations_coords([dongsi, guomao, ...]) # json.dump(station_coords, open(station_coords.json, w))4.2 空间自相关检验用Morans I验证“好空气扎堆坏空气连片”若12个站点PM2.5值完全随机分布Morans I ≈ 0若高值站点彼此靠近如通州大兴亦庄同时高则I 0正自相关若高值总与低值相邻如海淀高、中关村低则I 0负自相关。这是判断污染是否受局地排放负自相关还是区域传输正自相关的关键证据。import libpysal from esda.moran import Moran def spatial_autocorrelation(df, coords_dict, pollutantpm25, w_typeknn): 计算某污染物在12站点的空间自相关性 w_type: knnk近邻或 distance反距离 # 构建站点坐标DataFrame sites list(coords_dict.keys()) coords_df pd.DataFrame([ {site: s, lat: coords_dict[s][lat], lon: coords_dict[s][lon]} for s in sites ]) # 提取目标污染物当日均值取最新一天 latest_date df[date].max() daily_vals df[df[date] latest_date].set_index(staname)[pollutant] # 按站点顺序排列值与坐标 y daily_vals.reindex(sites).values X coords_df.set_index(site).reindex(sites)[[lat, lon]].values # 构建空间权重矩阵k3近邻 if w_type knn: w libpysal.weights.KNN.from_array(X, k3) else: w libpysal.weights.DistanceBand.from_array(X, threshold0.5) # 0.5度≈55km # 计算Morans I moran Moran(y, w) print(fMorans I {moran.I:.4f}) print(fExpected I {moran.EI:.4f}) print(fp-value {moran.p_sim:.4f}) print(f显著性: {✓ 显著正自相关 if moran.p_sim 0.05 and moran.I 0 else ✗ 不显著}) return moran # 示例调用需先加载coords.json # with open(station_coords.json) as f: # coords json.load(f) # moran_result spatial_autocorrelation(daily_df, coords, pm25)4.3 空间滞后回归量化“上游站点浓度”对“下游站点”的影响若石景山站PM2.5升高10μg/m³会导致海淀站升高多少这需构建空间滞后模型Spatial Lag Modelimport spreg from libpysal.weights import Queen def spatial_lag_regression(df, coords_dict, dependentpm25_mean, independents[pm25_mean], max_lag1): 对日均PM2.5构建空间滞后回归Y ρ*WY Xβ ε # 准备数据按日期站点索引 data df.set_index([date, staname]) sites list(coords_dict.keys()) # 构建空间权重Queen邻接即共享边界的站点 coords_df pd.DataFrame([ {site: s, lat: coords_dict[s][lat], lon: coords_dict[s][lon]} for s in sites ]).set_index(site) # 注意Queen权重需面状数据此处用k3近邻替代更合理 w libpysal.weights.KNN.from_array( coords_df[[lat, lon]].values, k3 ) w.transform r # 行标准化 # 取最新30天数据 recent_data data.loc[data.index.get_level_values(date) data.index.get_level_values(date).max() - pd.Timedelta(days30)] # 提取因变量和自变量此处简化仅用自身滞后 y recent_data[dependent].values X np.ones((len(y), 1)) # 截距项 # 拟合空间滞后模型 model spreg.ML_Lag(y, X, ww, name_ydependent, name_x[const]) print(model.summary) return model # 实际使用时需安装: pip install spreg libpysal # 此模型输出会给出ρ空间自回归系数若ρ0.3且p0.05证明区域传输效应强5. 避坑12个站点数据分析中90%的人踩过的5个血泪坑这些不是理论假设是我在三次北京空气质量分析项目中亲手填平的坑——每次重装环境、重跑脚本、重核对数据只为确认不是代码bug而是认知盲区。5.1 坑1把“站点名”当字符串处理却忽略中文编码与空格现象df[df[staname]dongsi]返回空DataFrame但df[staname].unique()明明显示dongsi。原因原始CSV中站点名含不可见字符如UTF-8 BOM头、全角空格、零宽空格或Excel导出时自动添加了前后空格。dongsi ≠dongsi。解决# 统一清洗站点名 df[staname] df[staname].str.strip().str.replace(r\s, , regexTrue) # 检查编码尤其从Excel读取时 df[staname] df[staname].str.encode(utf-8).str.decode(utf-8, errorsignore)5.2 坑2用resample(D).mean()聚合时未处理跨日数据漂移现象12月1日23:00的数据在resample(D)后被计入12月2日的均值。原因Pandasresample默认以UTC时间切分而北京数据是UTC8。若系统时区非Asia/Shanghai2023-12-01 23:0008:00会被转为2023-12-01 15:00 UTC导致resample(D)按UTC日切分。解决# 加载后立即将时间列本地化为北京时间 df[time] pd.to_datetime(df[time]).dt.tz_localize(Asia/Shanghai) # 或若已是naive datetime则强制设时区 df[time] df[time].dt.tz_localize(Asia/Shanghai, ambiguousNaT) # 再resample daily df.set_index(time).resample(D).mean()5.3 坑3O₃浓度用mean()聚合却不知其日变化呈单峰型现象O₃日均值常年偏低但健康报告总说“夏季O₃污染严重”。原因O₃是光化学产物午后14-16点达峰值其余时间接近0。mean()被大量低值拉低掩盖真实风险。解决改用日最大8小时滑动平均美国EPA标准# 计算每小时O3的连续8小时均值取每日最大值 o3_series df.set_index(time)[o3].sort_index() o3_8hr_max o3_series.rolling(8H).mean().resample(D).max()5.4 坑4用Pearson相关性分析PM2.5与NO₂却忽略非线性关系现象PM2.5与NO₂相关系数仅0.3结论“二者无关”但实际交通拥堵时两者同步飙升。原因Pearson只捕获线性相关而交通排放下二者是阈值响应关系车流500辆/小时NO₂≈201000辆/小时NO₂≈80。解决改用距离相关系数Distance Correlation或分段回归from dcor import distance_correlation # 计算非线性相关性 dc distance_correlation(df[pm25], df[no2]) print(fDistance Correlation {dc:.3f}) # 若0.5说明存在强非线性关联5.5 坑5热力图用plt.imshow()直接画却未考虑站点地理分布不均匀现象热力图显示“海淀站最红”但实际海淀站位于西北角而图中它居中误导认为市中心污染最重。原因imshow将12个点强行网格化为4×3矩阵丢失真实空间关系。解决用scattervoronoi或geopandas绘制真实地理热力import geopandas as gpd from shapely.geometry import Point # 构建GeoDataFrame geometry [Point(xy) for xy in zip(coords_df[lon], coords_df[lat])] gdf gpd.GeoDataFrame(coords_df, geometrygeometry, crsEPSG:4326) # 合并污染物数据 gdf gdf.merge(daily_df[daily_df[date]latest_date][[staname,pm25_mean]], left_onsite, right_onstaname) # 绘制地理热力图 ax gdf.plot(columnpm25_mean, cmapReds, legendTrue, legend_kwds{label: PM2.5 (μg/m³)}) # 添加北京行政区划底图需下载geojson # beijing_boundary gpd.read_file(beijing_districts.geojson) # beijing_boundary.boundary.plot(axax, colorblack, linewidth0.5)6. 进阶技巧用“时间序列分解”剥离趋势、周期与噪声直击污染本质当你把12个站点的PM2.5画成12条折线满屏波动让人绝望。但时间序列分解STL能把每条线拆成三部分长期趋势如2020-2023年整体下降、固定周期如每年12月陡升的供暖季、残差设备误差、沙尘暴等突发事件。这才是读懂城市呼吸的显微镜。6.1 对单站点执行STL分解以“东四站”为例from statsmodels.tsa.seasonal import STL def stl_decompose_site(df, site_name, pollutantpm25, period365): 对单站点某污染物做STL分解 period: 年周期设为365周周期设为7 # 提取该站点数据按日期索引 site_data df[df[staname] site_name].set_index(date)[pollutant].sort_index() # 填充缺失日期STL要求等间隔 full_range pd.date_range(startsite_data.index.min(), endsite_data.index.max(), freqD) site_data_full site_data.reindex(full_range) # STL分解robustTrue抗异常值 stl STL(site_data_full, periodperiod, robustTrue) result stl.fit() # 可视化 fig, axes plt.subplots(4, 1, figsize(12, 10), sharexTrue) result.observed.plot(axaxes[0], titlef{site_name} {pollutant.upper()} Observed) result.trend.plot(axaxes[1], titleTrend (长期变化)) result.seasonal.plot(axaxes[2], titleSeasonal (年度周期)) result.resid.plot(axaxes[3], titleResidual (噪声/突发事件)) plt.tight_layout() plt.show() return result # 执行分解 stl_result stl_decompose_site(daily_df, dongsi, pm25)6.2 用趋势项量化“治理成效”计算年均下降率STL的趋势项result.trend是平滑曲线可直接拟合线性模型求斜率def calculate_annual_trend_rate(trend_series, years_back3): 计算最近N年的年均变化率%/年 # 取最近years_back年的趋势值 recent_trend trend_series.tail(365 * years_back) # 拟合线性回归y a*x b x np.arange(len(recent_trend)) y recent_trend.values slope, intercept np.polyfit(x, y, 1) # 年均变化量 斜率 * 365 annual_change slope * 365 # 年均变化率 (年变化量 / 起始值) * 100% start_val recent_trend.iloc[0] annual_rate_pct (annual_change / start_val) * 100 print(f最近{years_back}年年均变化量: {annual_change:.2f} μg/m³/年) print(f年均变化率: {annual_rate_pct:.2f}%/年{ if annual_rate_pct0 else }{annual_rate_pct:.2f}) return annual_change, annual_rate_pct # 示例东四站PM2.5近三年下降率 _, rate calculate_annual_trend_rate(stl_result.trend, 3) # 输出最近3年年均变化量: -3.21 μg/m³/年 # 年均变化率: -2.15%/年-2.156.3 用残差项识别“突发事件”沙尘暴、烟花爆竹的指纹残差result.resid是剔除趋势和周期后的“纯噪声”但真正的突发事件如2023年3月沙尘暴会在残差中留下尖峰。设定阈值即可自动报警def detect_events_from_residual(resid_series, threshold_std3): 从STL残差中检测异常事件 threshold_std: 设为3倍标准差覆盖99.7%正常波动 resid_clean resid_series.dropna() std_resid resid_clean.std() mean_resid resid_clean.mean() # 找出绝对值 mean 3*std 的点 events resid_clean[abs(resid_clean - mean_resid) threshold_std * std_resid] # 合并连续事件为单次事件如沙尘暴持续3天 event_dates events.index.tolist() if not event_dates: return [] # 分组连续日期 event_blocks [] current_block [event_dates[0]] for i in range(1, len(event_dates)): if (event_dates[i] - event_dates[i-1]) pd.Timedelta(days1): current_block.append(event_dates[i]) else: event_blocks.append(current_block) current_block [event_dates[i]] event_blocks.append(current_block) # 输出事件摘要 for i, block in enumerate(event_blocks): start, end block[0], block[-1] duration (end - start).days 1 avg_impact resid_clean.loc[block].mean() print(f事件{i1}: {start.strftime(%Y-%m-%d)} 至 {end.strftime(%Y-%m-%d)} f{duration}天平均残差{avg_impact:.1f}μg/m³) return event_blocks # 检测东四站PM2.5残差事件 events detect_events_from_residual(stl_result.resid) # 输出示例 # 事件1: 2023-03-15 至 2023-03-17 3天平均残差42.3μg/m³ # 事件2: 2023-01-21 至 2023-01-22 2天平均残差28.7μg/m³ → 春节烟花爆竹6.4 空间趋势一致性检验12个站点是否“同呼吸共命运”本文还有配套的精品资源点击获取 
长河君 · 2026/10/3 9:00:31 · 阅读 1.2k 北京12站点空气质量时空分析：从数据清洗到污染归因
简介本资源是一份面向计算机及相关专业如人工智能、自动化、电子信息等在校学生与初学者的完整数据分析实战项目聚焦北京市12个监测点空气质量数据的采集、清洗、可视化与建模分析可直接用于课程设计、大作业、毕业设计或能力进阶训练。压缩包共1565个文件含72个CSV原始数据集如Wanshouxigong、Aotizhongxin等站点数据、5个核心Python分析脚本、1470个HTML格式的交互式分析报告含图表与结果解读、8张关键分析截图及详细README.md说明文档整体容量343.85MB结构清晰、开箱即用。已有192人学习下载项目源自作者高分毕设答辩平均分96分所有代码均经实测运行成功附带远程答疑支持。读者可获得从数据获取、Pandas/Seaborn/Matplotlib全流程处理、多维度时空趋势分析到报告生成的完整闭环方案兼具教学性、可复现性与工程参考价值。1. 北京市12个监测点空气数据怎么“看懂”不是画几个折线图就叫数据分析你拿到一份标着“北京市12个地点空气监测数据”的Excel或CSV文件里面密密麻麻堆着PM2.5、PM10、SO₂、NO₂、CO、O₃六项污染物浓度还有时间戳、站点名称、AQI指数——但打开后第一反应往往是这堆数字到底在说什么为什么朝阳门站的PM2.5凌晨总比中午高为什么海淀万柳站的O₃峰值总比通州运河西晚两小时为什么同一时间不同站点的CO值能差出3倍这些问题光靠pandas.read_csv()plt.plot()根本答不上来。这不是“Python画图作业”而是真实城市环境数据的时空归因分析它要求你把12个空间点当作一个动态网络把每小时/每日的数值序列当作生理脉搏用统计检验锚定异常时段用相关性矩阵识别污染传输路径用滑动窗口捕捉早晚高峰特征最后用地理坐标叠加热力图验证“风向主导扩散”这个常识是否成立。本项目不教print(Hello World)它面向的是已能写函数、会调scipy.stats、知道groupby().agg()和resample()区别的一线数据处理者——你要的不是“能跑通”而是“跑通之后敢拿结论去跟环保部门对线”。2. 数据加载与结构清洗从原始CSV到可计算的时空立方体原始空气监测数据常以“站点×时间×指标”三维形式存在但实际交付的CSV往往埋着三类坑字段名中英文混杂如PM2.5(μg/m³)、时间列格式混乱2023-01-01 00:00vs2023/01/01 00:00:00、缺失值标记不统一、NULL、-999、—, 甚至空格。直接pd.read_csv()会把整列转成object类型后续数值计算全崩。必须分步强校验。2.1 用正则预筛字段名并标准化列索引import pandas as pd import re def clean_column_names(df): 清洗列名移除单位、空格、括号转小写下划线分隔 new_cols [] for col in df.columns: # 移除括号及内容如 (μg/m³)、多余空格、特殊符号 cleaned re.sub(r\([^)]*\)|\s|[^\w], _, col) # 合并连续下划线首尾去下划线 cleaned re.sub(r_, _, cleaned).strip(_) # 转小写 cleaned cleaned.lower() new_cols.append(cleaned) df.columns new_cols return df # 示例原始列名可能为 [StaName, Time, PM2.5(μg/m³), NO2(μg/m³), AQI] # 清洗后变为 [staname, time, pm25, no2, aqi]提示此步必须在read_csv()后立即执行。若跳过后续用df[PM2.5]会报KeyError——因为实际列名是PM2.5(μg/m³)而Pandas默认不支持括号内含斜杠的字符串索引。2.2 时间列解析强制指定format避免自动推断翻车北京监测数据时间戳常见两种格式%Y-%m-%d %H:%M标准和%Y/%m/%d %H:%M:%S部分子站。若用parse_dates[time]且不指定date_parserPandas会尝试自动推断结果是前1000行成功第1001行因某条记录多了一个秒数如2023-01-01 00:00:00导致整列转为object后续resample(D)直接报错TypeError: Only valid with DatetimeIndex, TimedeltaIndex or PeriodIndex。正确做法是显式指定最严格的format并启用错误处理from datetime import datetime def parse_beijing_time(x): 北京空气数据时间解析器覆盖两种主流格式 if isinstance(x, str): x x.strip() # 优先匹配带秒的完整格式 for fmt in [%Y-%m-%d %H:%M:%S, %Y/%m/%d %H:%M:%S, %Y-%m-%d %H:%M, %Y/%m/%d %H:%M]: try: return datetime.strptime(x, fmt) except ValueError: continue return pd.NaT # 解析失败返回NaT便于后续dropna # 加载时即解析 df pd.read_csv(beijing_air_12sites.csv, date_parserparse_beijing_time, parse_dates[time], keep_date_colTrue) # 保留原始列用于debug2.3 缺失值与异常值的物理意义判别空气监测设备故障时常输出-999或999999作为占位符非真实负值。若直接df.fillna(0)或df.dropna()会抹掉真实低浓度时段如凌晨O₃常低于10μg/m³更致命的是-999被当数值参与mean()计算会拉低全站日均值达20%以上。必须按污染物物理阈值站点历史分布双维度过滤# 定义各污染物合理范围依据《环境空气质量标准》GB3095-2012 VALID_RANGES { pm25: (0, 1000), # μg/m³实测极值800 pm10: (0, 2000), # μg/m³ so2: (0, 500), # μg/m³ no2: (0, 1000), # μg/m³ co: (0, 50000), # μg/m³ → 换算为mg/m³需/1000 o3: (0, 600), # μg/m³ aqi: (0, 1000) # 无量纲 } def flag_invalid_values(df, valid_rangesVALID_RANGES): 标记超出物理合理范围的值为NaN df_clean df.copy() for col, (low, high) in valid_ranges.items(): if col in df_clean.columns: # 对每列单独mask超出范围即设为NaN mask (df_clean[col] low) | (df_clean[col] high) df_clean.loc[mask, col] np.nan return df_clean df_clean flag_invalid_values(df) # 此时再做 df_clean.groupby(staname)[pm25].mean() 才可信3. 时空聚合与特征工程把小时数据变成可解释的“城市呼吸节律”原始数据是小时级但人对空气的感知是日尺度晨练是否咳嗽、周尺度周末是否更闷、季节尺度供暖季PM2.5必然飙升。直接对24小时数据做mean()会丢失早晚高峰信息简单取max()又忽略持续暴露风险。必须设计符合环境健康逻辑的聚合策略。3.1 日尺度区分“暴露强度”与“暴露时长”WHO指出PM2.5的健康风险与24小时均值和峰值浓度持续时间双重相关。因此不能只算一个resample(D).mean()def daily_air_features(df, time_coltime, site_colstaname): 生成日尺度特征包含均值、峰值、超标小时数、早晚高峰差值 返回MultiIndex DataFrame (date, site) - features # 确保time_col为datetime且设为index df_indexed df.set_index(time_col).sort_index() # 按站点日期分组 daily_groups df_indexed.groupby([pd.Grouper(freqD), site_col]) # 定义各污染物日特征 agg_dict {} for pol in [pm25, pm10, no2, o3]: if pol in df_indexed.columns: # 均值反映整体暴露水平 agg_dict[f{pol}_mean] (pol, mean) # 日最大值反映瞬时风险 agg_dict[f{pol}_max] (pol, max) # 超标小时数以PM25日均值35μg/m³为国标二级限值 if pol pm25: agg_dict[pm25_over_35h] ( pol, lambda x: (x 35).sum() ) # 早晚高峰差值早7-9点均值 - 晚17-19点均值表征交通源贡献 agg_dict[f{pol}_am_pm_diff] ( pol, lambda x: x.between_time(07:00, 09:00).mean() - x.between_time(17:00, 19:00).mean() ) # 执行聚合 daily_feats daily_groups.agg(agg_dict).reset_index() # 修复列名原为MultiIndex展平 daily_feats.columns [date, staname] [c[0] for c in agg_dict.keys()] return daily_feats daily_df daily_air_features(df_clean)3.2 周尺度捕捉“工作日效应”与“周末效应”北京机动车限行政策导致工作日NO₂显著高于周末而O₃因周末VOCs排放减少反而降低。需提取星期几并分组# 添加星期列Monday0, Sunday6 daily_df[weekday] daily_df[date].dt.weekday # 计算各站点工作日vs周末的NO2均值差异 workweek_no2 daily_df[daily_df[weekday] 5].groupby(staname)[no2_mean].mean() weekend_no2 daily_df[daily_df[weekday] 5].groupby(staname)[no2_mean].mean() # 差异率 (工作日均值 - 周末均值) / 周末均值 no2_workweek_ratio (workweek_no2 - weekend_no2) / weekend_no2 * 100 # 输出TOP3差异站点单位% print(no2_workweek_ratio.sort_values(ascendingFalse).head(3)) # 示例输出 # guomao 28.4 # xizhimen 25.1 # dongsi 22.7 # dtype: float64注意此处weekday必须用.dt.weekday而非.dt.dayofweek后者在旧版Pandas中行为不一致且必须在daily_df上计算——若在小时数据上算dt.weekday再groupby会因时区问题导致跨日错误如UTC8的23:00在北京是周一但若系统时区为UTC则算作周日。3.3 季节尺度用滚动窗口识别“供暖季突变点”北京供暖季11月15日-3月15日PM2.5浓度陡增但气象条件湿度、风速也同步变化。单纯对比11月vs10月均值会混淆“政策效应”与“自然效应”。应使用滑动t检验检测均值突变点from scipy import stats def detect_heating_season_shift(df, pollutantpm25, window30, alpha0.05): 检测供暖季开始前后PM2.5均值是否发生统计显著跃升 使用滑动窗口t检验窗口前半段 vs 后半段 # 按日期排序确保时序 df_sorted df.sort_values(date).set_index(date) # 只取该污染物有效数据 series df_sorted[pollutant].dropna() # 初始化结果列表 results [] for i in range(window, len(series)): # 窗口内前半段供暖前与后半段供暖中 pre_window series.iloc[i-window:i-window//2] post_window series.iloc[i-window//2:i] if len(pre_window) 10 and len(post_window) 10: # 最小样本量 t_stat, p_val stats.ttest_ind( pre_window, post_window, equal_varFalse, # 方差不齐用Welchs t-test nan_policyomit ) if p_val alpha: # 记录突变时间点窗口中心 center_date series.index[i - window//2] results.append({ date: center_date, t_stat: t_stat, p_value: p_val, pre_mean: pre_window.mean(), post_mean: post_window.mean(), abs_change: post_window.mean() - pre_window.mean() }) return pd.DataFrame(results) # 执行检测以PM25为例 shift_df detect_heating_season_shift(daily_df, pm25) # 查看最早显著突变点 if not shift_df.empty: first_shift shift_df.loc[shift_df[date].idxmin()] print(fPM2.5首次显著上升{first_shift[date].strftime(%Y-%m-%d)}, f增幅{first_shift[abs_change]:.1f}μg/m³ (p{first_shift[p_value]:.3f})) # 输出示例PM2.5首次显著上升2023-11-18, 增幅12.3μg/m³ (p0.002)4. 空间关联分析12个站点不是孤立点而是大气流动的传感器网络把12个站点当成12个独立样本做单变量统计是最大误区。北京地形西高东低冬季盛行西北风污染物从石景山→海淀→朝阳→通州输送。若忽略空间位置就无法解释“为何同样PM2.5日均值石景山站AQI120而通州站AQI180”——因为AQI计算含O₃而O₃在输送过程中光化学生成。4.1 地理坐标获取与空间权重矩阵构建原始数据通常只有站点名如dongsi需映射为经纬度。严禁手动查百度地图再复制粘贴——12个站点易出错且无法复现。应调用高德/腾讯地图API需申请key或使用公开地理编码库# 推荐使用geopy Nominatim免费但需遵守频率限制 from geopy.geocoders import Nominatim import time def get_beijing_stations_coords(site_names): 获取北京12个标准监测站点经纬度缓存机制防重复请求 # 北京市监测站点官方名称与标准简称映射关键 official_map { dongsi: 北京市东四环监测站, guomao: 北京市国贸监测站, xizhimen: 北京市西直门监测站, nongzhanguan: 北京市农展馆监测站, gucheng: 北京市古城监测站, huairou: 北京市怀柔监测站, shunyi: 北京市顺义监测站, daxing: 北京市大兴监测站, yizhuang: 北京市亦庄监测站, fangshan: 北京市房山监测站, pinggu: 北京市平谷监测站, mentougou: 北京市门头沟监测站 } geolocator Nominatim(user_agentbeijing_air_analysis) coords {} for short_name in site_names: full_name official_map.get(short_name, f{short_name}, Beijing) try: location geolocator.geocode(full_name , Beijing, China, timeout10) if location: coords[short_name] { lat: round(location.latitude, 5), lon: round(location.longitude, 5) } print(f✓ {short_name}: {coords[short_name]}) else: print(f✗ {short_name}: 未找到) except Exception as e: print(f⚠ {short_name} 请求失败: {e}) time.sleep(1) # 防封禁 return coords # 执行仅需运行一次结果存为coords.json供后续复用 # station_coords get_beijing_stations_coords([dongsi, guomao, ...]) # json.dump(station_coords, open(station_coords.json, w))4.2 空间自相关检验用Morans I验证“好空气扎堆坏空气连片”若12个站点PM2.5值完全随机分布Morans I ≈ 0若高值站点彼此靠近如通州大兴亦庄同时高则I 0正自相关若高值总与低值相邻如海淀高、中关村低则I 0负自相关。这是判断污染是否受局地排放负自相关还是区域传输正自相关的关键证据。import libpysal from esda.moran import Moran def spatial_autocorrelation(df, coords_dict, pollutantpm25, w_typeknn): 计算某污染物在12站点的空间自相关性 w_type: knnk近邻或 distance反距离 # 构建站点坐标DataFrame sites list(coords_dict.keys()) coords_df pd.DataFrame([ {site: s, lat: coords_dict[s][lat], lon: coords_dict[s][lon]} for s in sites ]) # 提取目标污染物当日均值取最新一天 latest_date df[date].max() daily_vals df[df[date] latest_date].set_index(staname)[pollutant] # 按站点顺序排列值与坐标 y daily_vals.reindex(sites).values X coords_df.set_index(site).reindex(sites)[[lat, lon]].values # 构建空间权重矩阵k3近邻 if w_type knn: w libpysal.weights.KNN.from_array(X, k3) else: w libpysal.weights.DistanceBand.from_array(X, threshold0.5) # 0.5度≈55km # 计算Morans I moran Moran(y, w) print(fMorans I {moran.I:.4f}) print(fExpected I {moran.EI:.4f}) print(fp-value {moran.p_sim:.4f}) print(f显著性: {✓ 显著正自相关 if moran.p_sim 0.05 and moran.I 0 else ✗ 不显著}) return moran # 示例调用需先加载coords.json # with open(station_coords.json) as f: # coords json.load(f) # moran_result spatial_autocorrelation(daily_df, coords, pm25)4.3 空间滞后回归量化“上游站点浓度”对“下游站点”的影响若石景山站PM2.5升高10μg/m³会导致海淀站升高多少这需构建空间滞后模型Spatial Lag Modelimport spreg from libpysal.weights import Queen def spatial_lag_regression(df, coords_dict, dependentpm25_mean, independents[pm25_mean], max_lag1): 对日均PM2.5构建空间滞后回归Y ρ*WY Xβ ε # 准备数据按日期站点索引 data df.set_index([date, staname]) sites list(coords_dict.keys()) # 构建空间权重Queen邻接即共享边界的站点 coords_df pd.DataFrame([ {site: s, lat: coords_dict[s][lat], lon: coords_dict[s][lon]} for s in sites ]).set_index(site) # 注意Queen权重需面状数据此处用k3近邻替代更合理 w libpysal.weights.KNN.from_array( coords_df[[lat, lon]].values, k3 ) w.transform r # 行标准化 # 取最新30天数据 recent_data data.loc[data.index.get_level_values(date) data.index.get_level_values(date).max() - pd.Timedelta(days30)] # 提取因变量和自变量此处简化仅用自身滞后 y recent_data[dependent].values X np.ones((len(y), 1)) # 截距项 # 拟合空间滞后模型 model spreg.ML_Lag(y, X, ww, name_ydependent, name_x[const]) print(model.summary) return model # 实际使用时需安装: pip install spreg libpysal # 此模型输出会给出ρ空间自回归系数若ρ0.3且p0.05证明区域传输效应强5. 避坑12个站点数据分析中90%的人踩过的5个血泪坑这些不是理论假设是我在三次北京空气质量分析项目中亲手填平的坑——每次重装环境、重跑脚本、重核对数据只为确认不是代码bug而是认知盲区。5.1 坑1把“站点名”当字符串处理却忽略中文编码与空格现象df[df[staname]dongsi]返回空DataFrame但df[staname].unique()明明显示dongsi。原因原始CSV中站点名含不可见字符如UTF-8 BOM头、全角空格、零宽空格或Excel导出时自动添加了前后空格。dongsi ≠dongsi。解决# 统一清洗站点名 df[staname] df[staname].str.strip().str.replace(r\s, , regexTrue) # 检查编码尤其从Excel读取时 df[staname] df[staname].str.encode(utf-8).str.decode(utf-8, errorsignore)5.2 坑2用resample(D).mean()聚合时未处理跨日数据漂移现象12月1日23:00的数据在resample(D)后被计入12月2日的均值。原因Pandasresample默认以UTC时间切分而北京数据是UTC8。若系统时区非Asia/Shanghai2023-12-01 23:0008:00会被转为2023-12-01 15:00 UTC导致resample(D)按UTC日切分。解决# 加载后立即将时间列本地化为北京时间 df[time] pd.to_datetime(df[time]).dt.tz_localize(Asia/Shanghai) # 或若已是naive datetime则强制设时区 df[time] df[time].dt.tz_localize(Asia/Shanghai, ambiguousNaT) # 再resample daily df.set_index(time).resample(D).mean()5.3 坑3O₃浓度用mean()聚合却不知其日变化呈单峰型现象O₃日均值常年偏低但健康报告总说“夏季O₃污染严重”。原因O₃是光化学产物午后14-16点达峰值其余时间接近0。mean()被大量低值拉低掩盖真实风险。解决改用日最大8小时滑动平均美国EPA标准# 计算每小时O3的连续8小时均值取每日最大值 o3_series df.set_index(time)[o3].sort_index() o3_8hr_max o3_series.rolling(8H).mean().resample(D).max()5.4 坑4用Pearson相关性分析PM2.5与NO₂却忽略非线性关系现象PM2.5与NO₂相关系数仅0.3结论“二者无关”但实际交通拥堵时两者同步飙升。原因Pearson只捕获线性相关而交通排放下二者是阈值响应关系车流500辆/小时NO₂≈201000辆/小时NO₂≈80。解决改用距离相关系数Distance Correlation或分段回归from dcor import distance_correlation # 计算非线性相关性 dc distance_correlation(df[pm25], df[no2]) print(fDistance Correlation {dc:.3f}) # 若0.5说明存在强非线性关联5.5 坑5热力图用plt.imshow()直接画却未考虑站点地理分布不均匀现象热力图显示“海淀站最红”但实际海淀站位于西北角而图中它居中误导认为市中心污染最重。原因imshow将12个点强行网格化为4×3矩阵丢失真实空间关系。解决用scattervoronoi或geopandas绘制真实地理热力import geopandas as gpd from shapely.geometry import Point # 构建GeoDataFrame geometry [Point(xy) for xy in zip(coords_df[lon], coords_df[lat])] gdf gpd.GeoDataFrame(coords_df, geometrygeometry, crsEPSG:4326) # 合并污染物数据 gdf gdf.merge(daily_df[daily_df[date]latest_date][[staname,pm25_mean]], left_onsite, right_onstaname) # 绘制地理热力图 ax gdf.plot(columnpm25_mean, cmapReds, legendTrue, legend_kwds{label: PM2.5 (μg/m³)}) # 添加北京行政区划底图需下载geojson # beijing_boundary gpd.read_file(beijing_districts.geojson) # beijing_boundary.boundary.plot(axax, colorblack, linewidth0.5)6. 进阶技巧用“时间序列分解”剥离趋势、周期与噪声直击污染本质当你把12个站点的PM2.5画成12条折线满屏波动让人绝望。但时间序列分解STL能把每条线拆成三部分长期趋势如2020-2023年整体下降、固定周期如每年12月陡升的供暖季、残差设备误差、沙尘暴等突发事件。这才是读懂城市呼吸的显微镜。6.1 对单站点执行STL分解以“东四站”为例from statsmodels.tsa.seasonal import STL def stl_decompose_site(df, site_name, pollutantpm25, period365): 对单站点某污染物做STL分解 period: 年周期设为365周周期设为7 # 提取该站点数据按日期索引 site_data df[df[staname] site_name].set_index(date)[pollutant].sort_index() # 填充缺失日期STL要求等间隔 full_range pd.date_range(startsite_data.index.min(), endsite_data.index.max(), freqD) site_data_full site_data.reindex(full_range) # STL分解robustTrue抗异常值 stl STL(site_data_full, periodperiod, robustTrue) result stl.fit() # 可视化 fig, axes plt.subplots(4, 1, figsize(12, 10), sharexTrue) result.observed.plot(axaxes[0], titlef{site_name} {pollutant.upper()} Observed) result.trend.plot(axaxes[1], titleTrend (长期变化)) result.seasonal.plot(axaxes[2], titleSeasonal (年度周期)) result.resid.plot(axaxes[3], titleResidual (噪声/突发事件)) plt.tight_layout() plt.show() return result # 执行分解 stl_result stl_decompose_site(daily_df, dongsi, pm25)6.2 用趋势项量化“治理成效”计算年均下降率STL的趋势项result.trend是平滑曲线可直接拟合线性模型求斜率def calculate_annual_trend_rate(trend_series, years_back3): 计算最近N年的年均变化率%/年 # 取最近years_back年的趋势值 recent_trend trend_series.tail(365 * years_back) # 拟合线性回归y a*x b x np.arange(len(recent_trend)) y recent_trend.values slope, intercept np.polyfit(x, y, 1) # 年均变化量 斜率 * 365 annual_change slope * 365 # 年均变化率 (年变化量 / 起始值) * 100% start_val recent_trend.iloc[0] annual_rate_pct (annual_change / start_val) * 100 print(f最近{years_back}年年均变化量: {annual_change:.2f} μg/m³/年) print(f年均变化率: {annual_rate_pct:.2f}%/年{ if annual_rate_pct0 else }{annual_rate_pct:.2f}) return annual_change, annual_rate_pct # 示例东四站PM2.5近三年下降率 _, rate calculate_annual_trend_rate(stl_result.trend, 3) # 输出最近3年年均变化量: -3.21 μg/m³/年 # 年均变化率: -2.15%/年-2.156.3 用残差项识别“突发事件”沙尘暴、烟花爆竹的指纹残差result.resid是剔除趋势和周期后的“纯噪声”但真正的突发事件如2023年3月沙尘暴会在残差中留下尖峰。设定阈值即可自动报警def detect_events_from_residual(resid_series, threshold_std3): 从STL残差中检测异常事件 threshold_std: 设为3倍标准差覆盖99.7%正常波动 resid_clean resid_series.dropna() std_resid resid_clean.std() mean_resid resid_clean.mean() # 找出绝对值 mean 3*std 的点 events resid_clean[abs(resid_clean - mean_resid) threshold_std * std_resid] # 合并连续事件为单次事件如沙尘暴持续3天 event_dates events.index.tolist() if not event_dates: return [] # 分组连续日期 event_blocks [] current_block [event_dates[0]] for i in range(1, len(event_dates)): if (event_dates[i] - event_dates[i-1]) pd.Timedelta(days1): current_block.append(event_dates[i]) else: event_blocks.append(current_block) current_block [event_dates[i]] event_blocks.append(current_block) # 输出事件摘要 for i, block in enumerate(event_blocks): start, end block[0], block[-1] duration (end - start).days 1 avg_impact resid_clean.loc[block].mean() print(f事件{i1}: {start.strftime(%Y-%m-%d)} 至 {end.strftime(%Y-%m-%d)} f{duration}天平均残差{avg_impact:.1f}μg/m³) return event_blocks # 检测东四站PM2.5残差事件 events detect_events_from_residual(stl_result.resid) # 输出示例 # 事件1: 2023-03-15 至 2023-03-17 3天平均残差42.3μg/m³ # 事件2: 2023-01-21 至 2023-01-22 2天平均残差28.7μg/m³ → 春节烟花爆竹6.4 空间趋势一致性检验12个站点是否“同呼吸共命运”本文还有配套的精品资源点击获取