# PM2.5浓度预测实战：从数据清洗到随机森林模型调优全解析 空气质量预测一直是环境科学和公共健康领域的重要课题。PM2.5作为大气污染的主要指标之一，其浓度变化直接影响着城市居民的生活质量和健康水平。传统统计方法在PM2.5预测中存在诸多局限，而机器学习技术，特别是集成学习方法，为解决这一问题提供了新的思路。 ## 1. 环境准备与数据清洗 在开始构建预测模型前，我们需要搭建合适的工作环境并准备数据。Python生态提供了丰富的工具链，从数据清洗到模型训练都能找到成熟的解决方案。 首先安装必要的依赖库： ```bash pip install pandas numpy scikit-learn matplotlib seaborn ``` 空气质量数据通常包含多种污染物指标和气象参数。以北京市空气质量数据为例，原始数据集可能包含以下字段： | 字段名 | 描述 | 单位 | 处理方式 | |--------|------|------|----------| | PM2.5 | 细颗粒物浓度 | μg/m³ | 目标变量 | | SO2 | 二氧化硫浓度 | μg/m³ | 异常值处理 | | NO2 | 二氧化氮浓度 | μg/m³ | 缺失值填充 | | CO | 一氧化碳浓度 | mg/m³ | 单位统一 | | O3 | 臭氧浓度 | μg/m³ | 时间序列对齐 | | TEMP | 温度 | ℃ | 标准化处理 | | RH | 相对湿度 | % | 范围检查 | 数据清洗是建模过程中最耗时的环节之一，但也是确保模型效果的基础。常见的数据问题及处理方法包括： - **异常值处理**：PM2.5浓度理论上不应为负值，也不应超过合理上限 - **缺失值填充**：可采用前后时间点插值或同类站点均值填充 - **单位统一**：确保所有数值在同一量纲下比较 - **时间对齐**：多源数据的时间戳需要统一 ```python import pandas as pd import numpy as np def clean_air_quality_data(df): # 处理异常值 df = df[(df['PM2.5'] > 0) & (df['PM2.5'] <= 500)] # 填充缺失值 for col in ['SO2', 'NO2', 'CO', 'O3']: df[col] = df[col].interpolate() # 单位转换：CO从mg/m³转为μg/m³ df['CO'] = df['CO'] * 1000 # 去除仍存在缺失值的行 df = df.dropna() return df ``` ## 2. 特征工程与数据探索 高质量的特征工程往往比模型选择更能提升预测性能。对于PM2.5预测，我们需要考虑以下几类特征： 1. **基础污染物特征**：SO2、NO2、CO、O3等污染物的当前浓度 2. **气象特征**：温度、湿度、风速、气压等气象条件 3. **时间特征**：小时、星期、月份等时间维度信息 4. **空间特征**：监测站点的地理位置信息 5. **衍生特征**：污染物比值、滑动平均值等统计特征 ```python from sklearn.preprocessing import StandardScaler from sklearn.model_selection import train_test_split def feature_engineering(df): # 添加时间特征 df['hour'] = df['time'].dt.hour df['weekday'] = df['time'].dt.weekday df['month'] = df['time'].dt.month # 添加污染物比值特征 df['SO2/NO2'] = df['SO2'] / df['NO2'] df['CO/O3'] = df['CO'] / df['O3'] # 添加滑动窗口特征 window_size = 6 # 6小时滑动窗口 for col in ['PM2.5', 'SO2', 'NO2', 'CO', 'O3']: df[f'{col}_rolling_mean'] = df[col].rolling(window_size).mean() df[f'{col}_rolling_std'] = df[col].rolling(window_size).std() # 去除因滑动窗口产生的缺失值 df = df.dropna() # 特征与目标变量分离 X = df.drop(['PM2.5', 'time'], axis=1) y = df['PM2.5'] # 数据标准化 scaler = StandardScaler() X_scaled = scaler.fit_transform(X) # 划分训练集和测试集 X_train, X_test, y_train, y_test = train_test_split( X_scaled, y, test_size=0.2, random_state=42 ) return X_train, X_test, y_train, y_test, scaler, X.columns ``` 数据探索是理解数据分布和特征间关系的重要步骤。通过可视化可以直观发现数据规律： ```python import seaborn as sns import matplotlib.pyplot as plt def explore_data(df): # PM2.5浓度分布 plt.figure(figsize=(10, 6)) sns.histplot(df['PM2.5'], bins=50, kde=True) plt.title('PM2.5 Concentration Distribution') plt.xlabel('PM2.5 (μg/m³)') plt.ylabel('Frequency') # 污染物相关性热力图 plt.figure(figsize=(12, 8)) corr_matrix = df[['PM2.5', 'SO2', 'NO2', 'CO', 'O3', 'TEMP', 'RH']].corr() sns.heatmap(corr_matrix, annot=True, cmap='coolwarm') plt.title('Pollutants Correlation Heatmap') # PM2.5时间序列趋势 plt.figure(figsize=(14, 6)) df.set_index('time')['PM2.5'].plot() plt.title('PM2.5 Time Series') plt.ylabel('PM2.5 (μg/m³)') ``` ## 3. 随机森林原理与实现 随机森林是一种基于Bagging思想的集成学习算法，通过构建多棵决策树并综合它们的预测结果来提高模型性能。与传统决策树相比，随机森林在两个方面引入了随机性： 1. **数据随机性**：通过自助采样法(Bootstrap)为每棵树生成不同的训练子集 2. **特征随机性**：在节点分裂时，从全部特征中随机选取部分特征进行最优分裂 这种双重随机性机制有效降低了模型的方差，提高了泛化能力。随机森林的主要优势包括： - 对高维数据适应性强 - 能够处理非线性关系 - 对异常值和噪声不敏感 - 提供特征重要性评估 在scikit-learn中实现基础随机森林模型： ```python from sklearn.ensemble import RandomForestRegressor from sklearn.metrics import mean_squared_error, r2_score def train_baseline_rf(X_train, X_test, y_train, y_test): # 初始化随机森林回归器 rf = RandomForestRegressor( n_estimators=100, # 树的数量 max_depth=None, # 树的最大深度 min_samples_split=2, # 分裂所需最小样本数 random_state=42, n_jobs=-1 # 使用所有CPU核心 ) # 训练模型 rf.fit(X_train, y_train) # 预测与评估 y_pred = rf.predict(X_test) mse = mean_squared_error(y_test, y_pred) r2 = r2_score(y_test, y_pred) print(f'Baseline Random Forest - MSE: {mse:.2f}, R2: {r2:.2f}') return rf ``` 随机森林提供了特征重要性评估，这有助于我们理解哪些因素对PM2.5浓度影响最大： ```python def plot_feature_importance(rf, feature_names): importances = rf.feature_importances_ indices = np.argsort(importances)[::-1] plt.figure(figsize=(12, 8)) plt.title('Feature Importance') plt.bar(range(len(importances)), importances[indices], align='center') plt.xticks(range(len(importances)), feature_names[indices], rotation=90) plt.tight_layout() ``` ## 4. 模型调优与验证 虽然随机森林在默认参数下通常表现不错，但通过调优可以进一步提升模型性能。主要调优参数包括： - **n_estimators**：树的数量，增加可提升性能但会增加计算成本 - **max_depth**：树的最大深度，控制模型复杂度 - **min_samples_split**：节点分裂所需最小样本数 - **max_features**：寻找最佳分裂时考虑的特征数量比例 使用网格搜索进行参数调优： ```python from sklearn.model_selection import GridSearchCV def tune_rf_parameters(X_train, y_train): param_grid = { 'n_estimators': [50, 100, 200], 'max_depth': [None, 10, 20, 30], 'min_samples_split': [2, 5, 10], 'max_features': ['auto', 'sqrt', 0.5] } rf = RandomForestRegressor(random_state=42, n_jobs=-1) grid_search = GridSearchCV( estimator=rf, param_grid=param_grid, cv=5, scoring='neg_mean_squared_error', verbose=2 ) grid_search.fit(X_train, y_train) print(f'Best parameters: {grid_search.best_params_}') print(f'Best MSE: {-grid_search.best_score_:.2f}') return grid_search.best_estimator_ ``` 随机森林特有的OOB(Out-of-Bag)误差可以作为模型性能的快速评估： ```python def evaluate_with_oob(X_train, y_train): rf_oob = RandomForestRegressor( n_estimators=200, max_depth=20, min_samples_split=5, max_features=0.5, oob_score=True, random_state=42, n_jobs=-1 ) rf_oob.fit(X_train, y_train) print(f'OOB Score: {rf_oob.oob_score_:.2f}') return rf_oob ``` 模型验证阶段，除了常规的MSE和R2指标外，还应该关注预测值与真实值的散点图，检查模型在不同浓度区间的表现： ```python def plot_predictions(y_test, y_pred): plt.figure(figsize=(10, 10)) plt.scatter(y_test, y_pred, alpha=0.5) plt.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'k--', lw=2) plt.xlabel('Actual PM2.5') plt.ylabel('Predicted PM2.5') plt.title('Actual vs Predicted PM2.5') # 添加误差统计信息 errors = y_pred - y_test plt.text(0.05, 0.9, f'MAE: {np.mean(np.abs(errors)):.2f}', transform=plt.gca().transAxes) plt.text(0.05, 0.85, f'MSE: {np.mean(errors**2):.2f}', transform=plt.gca().transAxes) plt.text(0.05, 0.8, f'R2: {r2_score(y_test, y_pred):.2f}', transform=plt.gca().transAxes) ``` ## 5. 模型部署与实时预测 训练好的模型可以部署到生产环境中进行实时PM2.5预测。以下是一个简单的预测流水线示例： ```python import joblib from datetime import datetime class PM25Predictor: def __init__(self, model_path, scaler_path): self.model = joblib.load(model_path) self.scaler = joblib.load(scaler_path) self.feature_names = None # 应在初始化时设置 def preprocess(self, input_data): """处理实时输入数据""" # 转换输入数据为DataFrame df = pd.DataFrame([input_data]) # 添加时间特征 current_time = datetime.now() df['hour'] = current_time.hour df['weekday'] = current_time.weekday() df['month'] = current_time.month # 确保特征顺序与训练时一致 df = df[self.feature_names] # 数据标准化 scaled_data = self.scaler.transform(df) return scaled_data def predict(self, input_data): """预测PM2.5浓度""" processed_data = self.preprocess(input_data) prediction = self.model.predict(processed_data) return prediction[0] ``` 模型保存与加载： ```python def save_model_artifacts(model, scaler, feature_names, model_path, scaler_path, features_path): joblib.dump(model, model_path) joblib.dump(scaler, scaler_path) joblib.dump(feature_names, features_path) def load_model_artifacts(model_path, scaler_path, features_path): model = joblib.load(model_path) scaler = joblib.load(scaler_path) feature_names = joblib.load(features_path) return model, scaler, feature_names ``` 在实际应用中，可以考虑将模型部署为REST API服务： ```python from flask import Flask, request, jsonify app = Flask(__name__) predictor = None # 应在启动时初始化 @app.route('/predict', methods=['POST']) def predict(): data = request.get_json() try: prediction = predictor.predict(data) return jsonify({ 'status': 'success', 'prediction': float(prediction) }) except Exception as e: return jsonify({ 'status': 'error', 'message': str(e) }), 400 if __name__ == '__main__': # 初始化预测器 predictor = PM25Predictor('models/rf_model.pkl', 'models/scaler.pkl') predictor.feature_names = joblib.load('models/feature_names.pkl') app.run(host='0.0.0.0', port=5000) ``` ## 6. 模型优化进阶 基础随机森林模型虽然强大，但仍有优化空间。以下是几种进阶优化策略： ### 6.1 特征选择优化 通过递归特征消除(RFE)选择最有价值的特征子集： ```python from sklearn.feature_selection import RFECV def feature_selection_rfe(X_train, y_train): rf = RandomForestRegressor( n_estimators=100, random_state=42, n_jobs=-1 ) rfecv = RFECV( estimator=rf, step=1, cv=5, scoring='neg_mean_squared_error', min_features_to_select=10, n_jobs=-1 ) rfecv.fit(X_train, y_train) print(f'Optimal number of features: {rfecv.n_features_}') plt.figure(figsize=(10, 6)) plt.plot(range(1, len(rfecv.cv_results_['mean_test_score']) + 1), -rfecv.cv_results_['mean_test_score']) plt.xlabel('Number of features selected') plt.ylabel('Cross validation score (MSE)') plt.title('Recursive Feature Elimination') return rfecv ``` ### 6.2 集成其他算法 将随机森林与其他算法集成，如梯度提升树(XGBoost, LightGBM)： ```python from xgboost import XGBRegressor from lightgbm import LGBMRegressor from sklearn.ensemble import VotingRegressor def ensemble_models(X_train, X_test, y_train, y_test): rf = RandomForestRegressor( n_estimators=200, max_depth=20, min_samples_split=5, max_features=0.5, random_state=42, n_jobs=-1 ) xgb = XGBRegressor( n_estimators=200, max_depth=6, learning_rate=0.1, random_state=42, n_jobs=-1 ) lgbm = LGBMRegressor( n_estimators=200, max_depth=6, learning_rate=0.1, random_state=42, n_jobs=-1 ) ensemble = VotingRegressor( estimators=[('rf', rf), ('xgb', xgb), ('lgbm', lgbm)], weights=[0.4, 0.3, 0.3] ) ensemble.fit(X_train, y_train) y_pred = ensemble.predict(X_test) mse = mean_squared_error(y_test, y_pred) r2 = r2_score(y_test, y_pred) print(f'Ensemble Model - MSE: {mse:.2f}, R2: {r2:.2f}') return ensemble ``` ### 6.3 时间序列特征优化 对于PM2.5这类时间序列数据，可以引入更复杂的时间特征： ```python def enhanced_time_features(df, time_col='time'): # 原始时间特征 df['hour'] = df[time_col].dt.hour df['weekday'] = df[time_col].dt.weekday df['month'] = df[time_col].dt.month # 周期性编码 df['hour_sin'] = np.sin(2 * np.pi * df['hour']/24) df['hour_cos'] = np.cos(2 * np.pi * df['hour']/24) df['weekday_sin'] = np.sin(2 * np.pi * df['weekday']/7) df['weekday_cos'] = np.cos(2 * np.pi * df['weekday']/7) df['month_sin'] = np.sin(2 * np.pi * df['month']/12) df['month_cos'] = np.cos(2 * np.pi * df['month']/12) # 滞后特征 for lag in [1, 2, 3, 6, 12, 24]: # 1小时到24小时的滞后 df[f'PM2.5_lag_{lag}'] = df['PM2.5'].shift(lag) # 滑动统计量 for window in [3, 6, 12, 24]: # 3小时到24小时的窗口 df[f'PM2.5_rolling_mean_{window}'] = df['PM2.5'].rolling(window).mean() df[f'PM2.5_rolling_std_{window}'] = df['PM2.5'].rolling(window).std() # 差分特征 df['PM2.5_diff_1'] = df['PM2.5'].diff(1) df['PM2.5_diff_3'] = df['PM2.5'].diff(3) # 去除因滞后和滑动窗口产生的缺失值 df = df.dropna() return df ```