机器学习在臭氧预测中的应用与优化

1. 臭氧预测问题概述

在德克萨斯州休斯顿地区，地面臭氧浓度超标是空气污染的主要表现形式之一。气象学家们发现，通过分析风速、温度等气象测量数据，可以预测次日地面臭氧浓度是否会达到触发公共空气污染警告的水平。这个问题在机器学习领域被称为"臭氧预测问题"，是一个典型的时间序列分类任务。

这个数据集包含了休斯顿地区1998-2004年共7年的每日气象观测数据，记录了是否出现臭氧浓度超标的情况。数据集有以下几个显著特点：

1. 包含72个气象变量，其中只有10个被环境科学家确认与臭氧形成直接相关
2. 数据存在大量缺失值，标记为"?"
3. 臭氧日(正样本)与非臭氧日(负样本)的比例严重不平衡(约6% vs 94%)
4. 臭氧日的出现具有明显的季节性特征，主要集中在夏季月份

在实际应用中，预测臭氧日的概率比简单二元分类更有价值，因为气象部门可以根据概率值决定是否发布预警以及预警级别。

2. 数据加载与初步分析

2.1 数据获取与加载

数据集可以从UCI机器学习仓库获取，我们使用8小时平均版本的数据。使用Pandas加载数据时需要注意以下几点：

from pandas import read_csv # 加载数据集 data = read_csv('eighthr.data', header=None, index_col=0, parse_dates=True, squeeze=True)

数据预处理的关键步骤包括：

1. 将日期列设为索引
2. 处理缺失值(初始用"?"标记)
3. 检查数据规模和类别分布

2.2 数据探索与可视化

通过简单的统计分析，我们可以看到数据的基本情况：

print(data.shape) # 输出数据维度 # (2534, 73) # 计算类别分布 counts = data.groupby(73).size() for i in range(len(counts)): percent = counts[i] / data.shape[0] * 100 print(f'Class={i}, total={counts[i]}, percentage={percent:.3f}%')

输出结果显示：

• 总样本数：2534天
• 非臭氧日(Class=0)：2374天(93.686%)
• 臭氧日(Class=1)：160天(6.314%)

绘制臭氧日的时间序列图可以明显看到臭氧日集中在每年夏季：

import matplotlib.pyplot as plt plt.plot(data.index, data.values[:,-1]) plt.show()

3. 数据预处理流程

3.1 缺失值处理

数据中存在大量缺失值，标记为"?"。我们初步采用以下策略：

1. 将"?"替换为0.0
2. 后续可以考虑更复杂的填补方法，如均值填补或基于时间序列的插值

values = data.values values[values=='?'] = 0.0

3.2 构建监督学习问题

我们将问题定义为：使用当天的气象观测数据预测次日是否为臭氧日。这需要将时间序列数据转换为监督学习格式：

from numpy import hstack, array supervised = [] for i in range(len(values) - 1): X, y = values[i, :-1], values[i + 1, -1] row = hstack((X, y)) supervised.append(row) supervised = array(supervised)

3.3 训练集与测试集划分

按照时间顺序划分数据集，使用前5年(约1803天)作为训练集，后2年(730天)作为测试集：

split = 365 * 2 train, test = supervised[:-split,:], supervised[-split:,:] train, test = train.astype('float32'), test.astype('float32')

4. 基准模型建立

4.1 朴素预测模型

作为基准，我们建立一个朴素模型：始终预测训练集中臭氧日的比例作为概率。

naive_prob = sum(train[:,-1]) / train.shape[0] print(f"朴素模型预测概率: {naive_prob:.4f}") # 输出: 0.0727

4.2 评估指标选择

我们使用Brier Score和Brier Skill Score作为评估指标：

1. Brier Score (BS): 概率预测的均方误差 $$BS = \frac{1}{N}\sum_{i=1}^N (f_i - o_i)^2$$

2. Brier Skill Score (BSS): 相对于基准模型的改进 $$BSS = 1 - \frac{BS_{model}}{BS_{ref}}$$

计算朴素模型的BS和BSS：

from sklearn.metrics import brier_score_loss # 朴素模型预测 yhat_naive = [naive_prob for _ in range(len(test))] bs_ref = brier_score_loss(test[:,-1], yhat_naive) # Brier Skill Score计算 bss_naive = (bs_ref - bs_ref) / (0 - bs_ref) print(f"Brier Score: {bs_ref:.6f}, Brier Skill Score: {bss_naive:.6f}")

5. 集成树模型比较

5.1 模型选择与配置

我们比较四种集成树模型，均设置100个基学习器：

1. BaggingClassifier
2. ExtraTreesClassifier
3. GradientBoostingClassifier
4. RandomForestClassifier

from sklearn.ensemble import (BaggingClassifier, ExtraTreesClassifier, GradientBoostingClassifier, RandomForestClassifier) models = { 'Bagging': BaggingClassifier(n_estimators=100), 'ExtraTrees': ExtraTreesClassifier(n_estimators=100), 'GBM': GradientBoostingClassifier(n_estimators=100), 'RandomForest': RandomForestClassifier(n_estimators=100) }

5.2 评估框架设计

由于集成方法具有随机性，我们进行10次重复实验取平均：

from sklearn.base import clone from numpy import mean def evaluate_model(model, trainX, trainy, testX, testy, n_repeats=10): scores = [] for _ in range(n_repeats): m = clone(model) m.fit(trainX, trainy) probs = m.predict_proba(testX)[:, 1] bs = brier_score_loss(testy, probs) bss = (bs_ref - bs) / bs_ref # 正数表示改进 scores.append(bss) return mean(scores), scores

5.3 模型性能比较

执行评估并可视化结果：

results = {} for name, model in models.items(): mean_bss, all_bss = evaluate_model(model, trainX, trainy, testX, testy) results[name] = all_bss print(f"{name:12s} | Mean BSS: {mean_bss:.4f}") # 绘制箱线图 plt.figure(figsize=(10,6)) plt.boxplot(results.values(), labels=results.keys()) plt.title('Brier Skill Score Comparison') plt.ylabel('BSS') plt.show()

典型输出结果：

Bagging | Mean BSS: 0.0698 ExtraTrees | Mean BSS: 0.1033 GBM | Mean BSS: 0.1198 RandomForest | Mean BSS: 0.1027

6. 梯度提升树调优

6.1 关键参数分析

基于初步结果，我们对梯度提升树(GBM)进行调优。关键参数包括：

1. 学习率(learning_rate)：控制每棵树的贡献
2. 最大深度(max_depth)：控制单棵树的复杂度
3. 子采样比例(subsample)：引入随机性防止过拟合
4. 特征采样比例(max_features)

6.2 网格搜索实现

使用网格搜索寻找最优参数组合：

from sklearn.model_selection import GridSearchCV param_grid = { 'learning_rate': [0.01, 0.05, 0.1], 'max_depth': [3, 5, 7], 'subsample': [0.6, 0.8, 1.0], 'max_features': ['sqrt', 'log2', None] } gsearch = GridSearchCV( estimator=GradientBoostingClassifier(n_estimators=100), param_grid=param_grid, scoring='neg_brier_score', cv=5, n_jobs=-1 ) gsearch.fit(trainX, trainy)

6.3 最优模型评估

获取最优参数并重新评估：

best_params = gsearch.best_params_ print("Best Parameters:", best_params) best_gbm = GradientBoostingClassifier(n_estimators=100, **best_params) mean_bss, _ = evaluate_model(best_gbm, trainX, trainy, testX, testy) print(f"Tuned GBM BSS: {mean_bss:.4f}")

7. 模型部署与监控

7.1 模型持久化

将训练好的模型保存为文件：

import joblib joblib.dump(best_gbm, 'ozone_prediction_gbm.pkl')

7.2 实时预测流程

构建预测流水线：

class OzonePredictor: def __init__(self, model_path): self.model = joblib.load(model_path) self.scaler = None # 如果有特征缩放 def preprocess(self, raw_data): # 实现与训练时相同的数据预处理 processed = raw_data.copy() processed[processed=='?'] = 0.0 return processed.astype('float32') def predict_proba(self, input_data): processed = self.preprocess(input_data) return self.model.predict_proba(processed)[:, 1]

7.3 性能监控指标

建议监控以下指标：

1. 每日预测概率与实际结果的Brier Score
2. 预测概率的校准曲线
3. 特征重要性变化
4. 预测结果的ROC曲线和PR曲线

8. 实际应用中的挑战与解决方案

8.1 类别不平衡问题

臭氧日占比仅约6%，我们采用以下策略：

1. 调整类别权重
2. 使用过采样技术(SMOTE)
3. 采用适合不平衡数据的评估指标(如F1-score)

from imblearn.over_sampling import SMOTE smote = SMOTE() trainX_res, trainy_res = smote.fit_resample(trainX, trainy)

8.2 特征工程优化

尝试以下特征工程方法：

1. 添加时间特征(月份、季节)
2. 创建滞后特征(前几天的气象数据)
3. 添加移动平均特征

# 示例：添加月份特征 trainX_with_month = np.hstack([trainX, data.index.month.values[:-split][:, None]])

8.3 模型集成策略

结合多个模型的优势：

1. 堆叠(Stacking)不同算法
2. 使用投票集成
3. 时间序列交叉验证

from sklearn.ensemble import StackingClassifier base_models = [ ('rf', RandomForestClassifier(n_estimators=100)), ('gbm', GradientBoostingClassifier(n_estimators=100)) ] stack = StackingClassifier( estimators=base_models, final_estimator=LogisticRegression() )

9. 扩展应用与未来方向

9.1 多步预测

将模型扩展为预测未来多天的臭氧概率：

def prepare_multi_step_data(data, n_steps): X, y = [], [] for i in range(len(data) - n_steps): X.append(data[i:i+n_steps, :-1].flatten()) y.append(data[i+n_steps, -1]) return np.array(X), np.array(y)

9.2 深度学习应用

尝试LSTM等深度学习模型：

from tensorflow.keras.models import Sequential from tensorflow.keras.layers import LSTM, Dense model = Sequential([ LSTM(64, input_shape=(n_steps, n_features)), Dense(1, activation='sigmoid') ]) model.compile(loss='binary_crossentropy', optimizer='adam')

9.3 可解释性分析

使用SHAP值分析模型决策：

import shap explainer = shap.TreeExplainer(best_gbm) shap_values = explainer.shap_values(trainX) shap.summary_plot(shap_values, trainX)

在实际部署这类空气质量预测系统时，需要特别注意模型的可解释性。气象部门和公共卫生决策者往往需要理解模型做出预测的依据，而不仅仅是接受预测结果。通过结合特征重要性分析和领域知识，我们可以构建既准确又可信赖的空气质量预警系统。