1. 项目概述从一道赛题到一套可复现的解决方案去年带队打完亚太杯数学建模竞赛APMCM后一直有学弟学妹来问C题的代码和思路。与其零散回复不如系统整理出来。这份代码分享不仅仅是几个脚本的堆砌它完整记录了我们团队从题目解析、模型构建、编程实现到结果分析的全过程。对于正在备赛亚太杯、美赛MCM/ICM甚至国赛的同学来说这类“全球变暖与极端天气”主题的题目具有很强的代表性和前瞻性。你拿到手的将是一个可以直接运行、模块清晰、注释详细的代码包更重要的是你能看到我们当时每一个技术决策背后的“为什么”——为什么用LSTM而不是ARIMA为什么特征工程要那样处理模型调参的坑是怎么踩过来的这篇文章我就以2022年APMCM C题为蓝本拆解我们的完整解题链条并附上经过优化和注释的MATLAB与Python双版本核心代码。无论你是建模新手还是有一定基础想提升代码实现能力的同学都能从中找到可以直接“抄作业”的模块和值得深思的建模思想。2. 赛题核心与解题框架设计2.1 题目回顾与问题本质剖析2022年APMCM C题的标题通常围绕“全球气候变化背景下极端天气事件的预测与评估”展开。题目会提供长时间序列的全球或区域气候数据如温度、降水、海平面压力等要求参赛者分析历史变化趋势建立模型预测未来极端天气事件如热浪、暴雨、干旱的频率与强度并评估其社会经济影响。这本质上是一个时间序列预测与风险评估相结合的复合型问题。面对这类问题新手容易陷入两个误区一是急于找代码、套模型忽视了对数据本身和问题定义的深刻理解二是试图用一个“超级模型”解决所有子问题导致结构臃肿解释性差。我们的策略是“分而治之阶段推进”。将整个问题分解为三个核心阶段数据洞察与特征工程阶段理解数据故事构造对预测有效的特征。核心预测模型构建阶段针对不同预测目标如趋势、极端事件选择合适的模型。综合评估与可视化阶段将预测结果转化为风险评估指标并用直观图表呈现。这个框架的优点是逻辑清晰每个阶段的目标明确且便于团队分工。例如数据处理能力强的同学主攻第一阶段擅长算法的同学聚焦第二阶段而写作和可视化能力突出的同学负责第三阶段。2.2 工具选型为什么是MATLAB与Python双线并行在工具选择上我们采用了MATLAB与Python并行的策略这并非炫技而是基于实际需求的最优解。MATLAB其强项在于强大的内置数学函数、优雅的矩阵运算以及专业的工具箱如Curve Fitting, Statistics and Machine Learning Toolbox。在前期数据探索、快速原型验证比如尝试不同的回归模型、以及生成出版级图表时MATLAB的效率非常高。对于题目中可能涉及的信号处理如滤波去噪或地理信息展示MATLAB也有现成的工具。Python其生态在深度学习、复杂机器学习流程以及大规模数据处理方面更具优势。当我们需要实现LSTM、GRU等循环神经网络或者进行复杂的特征交叉、使用XGBoost等集成学习模型时Python的TensorFlow/PyTorch和scikit-learn库是更自然的选择。此外网络爬虫如果需要补充数据和自动化报告生成也是Python的拿手好戏。我们的实操心得是在72小时的竞赛中“用什么工具最快实现想法”是第一原则。对于趋势拟合、统计检验等任务我们用MATLAB快速搞定对于需要调参、迭代的深度学习模型则用Python构建。代码仓库里两个文件夹/matlab_code和/python_code分别存放相应代码并在README.md中明确说明每个脚本的输入输出和用途。3. 数据预处理与特征工程全流程解析3.1 数据清洗与异常值处理实战竞赛提供的数据通常来自公开气候数据集如NASA GISS、CRU但并非“开箱即用”。原始数据可能包含缺失值记为-999.9或NaN、明显的录入错误以及由于仪器更替造成的非气候性跳变。我们的处理流程如下并附上关键代码片段缺失值处理对于时间序列简单删除可能导致序列断裂。我们采用线性插值对短期缺失和季节性插值对具有年周期性的数据如温度相结合的方法。% MATLAB: 使用fillmissing函数进行线性插值 temperature_data raw_data.Temperature; temperature_data_filled fillmissing(temperature_data, linear); % 对于具有年周期的数据可以考虑基于同期历史均值插值简化示例 monthly_means grpstats(temperature_data_filled, month(raw_dates), mean); % ... 根据月份索引用月均值填充仍存在的缺失值如需# Python (Pandas): 插值功能更灵活 import pandas as pd df[Temperature] df[Temperature].interpolate(methodlinear) # 线性插值 # 或者使用时间序列友好的方法如‘time’ df[Temperature] df[Temperature].interpolate(methodtime)异常值检测与修正我们使用滑动窗口统计法结合物理常识进行判断。例如计算每个数据点在一个滑动窗口如15年内的均值和标准差将超出均值 ± 3倍标准差范围的点视为统计异常。然后结合历史气候记录如已知的极端热浪年份判断是真正的极端事件还是错误数据。对于确认为错误的数据用插值结果替代。# Python 异常值检测示例 def detect_outliers_rolling(series, window180, n_sigmas3): rolling_mean series.rolling(windowwindow, centerTrue).mean() rolling_std series.rolling(windowwindow, centerTrue).std() outliers (series - rolling_mean).abs() (n_sigmas * rolling_std) return outliers注意事项处理气候数据异常值时必须谨慎直接删除或修正一个“异常高温”点可能会抹掉一次真实的极端热浪事件这恰恰可能是题目分析的重点。我们的原则是除非有充分证据表明是仪器错误否则优先保留并在后续分析中将其作为极端事件样本进行专门标记。3.2 构造对预测极端事件有效的特征原始的温度、降水序列是基础但直接将其丢进模型效果往往不佳。我们需要构造一些更能反映问题本质的特征Feature Engineering这是提升模型性能的关键。滞后特征这是时间序列预测的标配。不仅包含t-1,t-2时刻的值我们还构造了季节性滞后特征如t-12对于月度数据即去年同月、t-24等这对捕捉年周期规律至关重要。# Python 创建滞后特征 for lag in [1, 2, 3, 12, 24]: df[f‘Temperature_lag_{lag}’] df[‘Temperature’].shift(lag)滚动统计特征计算滑动窗口内的统计量可以平滑噪声并提取趋势信息。滚动均值反映短期背景状态。滚动标准差反映短期波动性或变率。波动性增大本身可能就是气候不稳定的信号。滚动最大值/最小值直接与极端事件阈值相关。df[‘Temp_rolling_mean_10y’] df[‘Temperature’].rolling(window120).mean() df[‘Temp_rolling_std_10y’] df[‘Temperature’].rolling(window120).std()交互特征与衍生指标温雨复合指数例如高温 × 干旱可能加剧热浪影响。可以构造温度异常 × (1 - 降水百分位数)这样的简单指数。极端事件阈值特征定义一个极端高温阈值如历史第95百分位数生成一个布尔序列Is_Extreme表示每一天是否超过该阈值。这个序列本身可以作为预测目标也可以作为特征。气候态特征计算每个日历日或月的长期如30年气候平均值然后用每日数据减去该平均值得到异常值序列。这个序列过滤掉了强烈的年周期更能突出长期趋势和异常信号。实操心得特征不是越多越好。过多的特征会导致维度灾难特别是对于数据量有限的气候序列。我们采用的方法是先基于气候学知识构造一批特征然后使用递归特征消除RFE或观察特征与目标变量的互信息筛选出最重要的10-15个特征用于最终模型训练。4. 核心预测模型从传统时序到LSTM的递进4.1 基准模型季节性自回归积分滑动平均模型在尝试复杂模型前建立一个稳健的基准模型是必要的。SARIMA模型是处理具有明显季节性的时间序列的标准方法。我们用它来预测未来几十年的平均气候态趋势如全球平均温度。这一步的目的有三个提供一个可解释的、稳健的趋势预测基线。其残差分析可以帮助我们理解数据中未被线性模型捕捉到的模式这些模式可能是非线性、非平稳的需要用更复杂的模型来处理。与后续复杂模型的结果进行对比量化机器学习模型带来的提升。在MATLAB中可以使用Econometrics Toolbox的arima和estimate函数。在Python中statsmodels库的SARIMAX是首选。关键步骤是确定(p,d,q)和(P,D,Q,s)参数我们通过观察自相关图、偏自相关图并结合AIC/BIC信息准则网格搜索来确定。# Python SARIMA示例 (需安装statsmodels) import statsmodels.api as sm # 假设已处理好数据series model sm.tsa.SARIMAX(series, order(1,1,1), seasonal_order(1,1,1,12)) result model.fit(dispFalse) forecast result.get_forecast(steps120) # 预测未来10年月度数据踩过的坑SARIMA对数据的平稳性要求高。即使做了差分和季节性差分有时残差仍不白噪声。这时不要强行调参到过拟合而应意识到数据的非线性果断转向机器学习模型。4.2 主力模型长短期记忆网络模型详解对于预测极端事件的发生频率和强度LSTM是我们的主力模型。因为极端事件往往由复杂的、非线性的气候系统内部相互作用和长期记忆效应所驱动LSTM非常适合捕捉这种长期依赖关系。4.2.1 数据准备与序列构造这是LSTM建模中最容易出错的一步。我们需要将时间序列转化为监督学习问题的格式(samples, timesteps, features)。samples训练样本数。timesteps回溯的时间步长如用过去10年的数据预测下一年。features每个时间步上的特征数量即我们之前构造的所有特征。import numpy as np from sklearn.preprocessing import StandardScaler def create_dataset(data, look_back120, forecast_horizon12): X, Y [], [] for i in range(len(data) - look_back - forecast_horizon): X.append(data[i:(i look_back), :]) # 取look_back步长的特征 Y.append(data[i look_back forecast_horizon - 1, target_idx]) # 预测未来forecast_horizon步的某个值如温度 return np.array(X), np.array(Y) # 假设df_features是包含所有特征的DataFrame scaler StandardScaler() scaled_data scaler.fit_transform(df_features) X, y create_dataset(scaled_data, look_back120, forecast_horizon12)4.2.2 网络结构设计与调参我们使用了相对经典的堆叠LSTM结构from tensorflow.keras.models import Sequential from tensorflow.keras.layers import LSTM, Dense, Dropout, Input model Sequential() model.add(Input(shape(X.shape[1], X.shape[2]))) # (timesteps, features) model.add(LSTM(units100, activation‘relu’, return_sequencesTrue)) model.add(Dropout(0.2)) # 防止过拟合 model.add(LSTM(units50, activation‘relu’)) model.add(Dropout(0.2)) model.add(Dense(25, activation‘relu’)) model.add(Dense(1)) # 回归问题输出一个值如温度 model.compile(optimizer‘adam’, loss‘mse’, metrics[‘mae’])关键调参经验units神经元数从50开始尝试根据数据量和复杂度增加。太大易过拟合太小欠拟合。我们通过验证集损失曲线来判定。Dropout在竞赛数据量有限的情况下Dropout是防止过拟合的利器。率设置在0.2到0.5之间。优化器Adam是默认且有效的选择。学习率lr可以使用回调函数ReduceLROnPlateau动态调整。早停必须使用EarlyStopping回调函数监控验证集损失耐心值patience设为15-20个epoch避免无效训练。4.2.3 损失函数与评估指标的选择对于极端事件预测我们不仅关心整体误差更关心对“尾部”极端值的预测能力。损失函数除了均方误差MSE我们还尝试了Huber Loss它对异常值的敏感度低于MSE在预测极端值时可能更稳健。评估指标整体RMSE均方根误差,MAE平均绝对误差。极端事件专项我们定义了Extreme Event Capture Rate在测试集中真实值超过极端阈值的样本里有多少被模型预测也超过了阈值或一个较低的预警阈值。这比单纯的误差指标更有业务意义。5. 模型集成、评估与结果可视化5.1 多模型结果对比与集成策略我们不会只依赖单一模型。除了SARIMA和LSTM我们还尝试了梯度提升树如XGBoost来处理特征间复杂的非线性关系。最终提交的预测结果往往来自多个模型的加权平均或分位数集成。集成方法简单加权平均根据各个模型在验证集上的RMSE倒数作为权重。final_forecast (w1 * forecast_lstm w2 * forecast_xgb) / (w1 w2)分位数回归森林/集成为了提供预测区间而不仅仅是点预测我们训练多个模型预测不同分位数如5% 50% 95%。这能给出未来气候值的可能范围对于风险评估至关重要。scikit-learn的QuantileRegressor或GradientBoostingRegressor的loss‘quantile’参数可以实现。模型评估我们严格在时间序列交叉验证的框架下进行评估。即按时间顺序划分训练集和验证集避免未来信息泄露。例如用1950-2000年的数据训练预测2001-2010年评估然后用1950-2010年训练预测2011-2020年如此滚动。5.2 可视化让结果自己说话在数学建模论文中一图胜千言。我们精心设计了以下几类图表历史序列与预测序列对比图将历史观测数据、多个模型的预测曲线以及预测区间阴影表示画在同一张图上。使用MATLAB的plot和fill函数或Python的matplotlib可以做出非常专业的图表。% MATLAB 绘制预测区间示例 figure; plot(historical_time, historical_data, ‘k-’, ‘LineWidth‘, 1.5); hold on; plot(forecast_time, forecast_mean, ‘b-’, ‘LineWidth‘, 2); fill([forecast_time; flipud(forecast_time)], [forecast_lower; flipud(forecast_upper)], ‘b’, ‘FaceAlpha‘, 0.2, ‘EdgeColor‘, ‘none’); legend(‘观测数据‘, ‘预测均值‘, ‘95%预测区间‘); xlabel(‘年份‘); ylabel(‘全球平均温度异常 (℃)’);极端事件频率变化图计算历史期和预测期内每年超过极端阈值的天数/次数绘制其变化趋势。可以使用滑动窗口计算频率并用条形图或线图展示清晰显示极端事件是否变得更多。空间分布图如果数据支持如果题目提供了格点数据可以用填色图展示未来极端温度或降水变化的空间分布。MATLAB的m_map工具箱或Python的cartopy库是绘制地理地图的利器。注意事项所有图表必须清晰标注坐标轴、单位、图例。颜色选择要直观如用红色表示增暖蓝色表示变冷并考虑黑白印刷时的可区分度。6. 代码使用指南与常见问题排查6.1 代码仓库结构与运行说明我们提供的代码包结构如下APMCM2022_C_Code/ ├── README.md # 项目总说明环境配置快速开始 ├── data/ # 数据存放目录原始数据需自行按说明放置 │ ├── raw/ # 原始数据 │ └── processed/ # 预处理后的数据由脚本生成 ├── matlab_code/ # MATLAB代码 │ ├── 01_data_preprocessing.m │ ├── 02_feature_engineering.m │ ├── 03_sarima_modeling.m │ ├── 04_visualization.m │ └── utils/ # 自定义函数 ├── python_code/ # Python代码 │ ├── requirements.txt # Python依赖包列表 │ ├── 01_data_preprocess.py │ ├── 02_feature_engineer.py │ ├── 03_lstm_model.py │ ├── 04_xgboost_model.py │ ├── 05_ensemble_evaluation.py │ └── utils.py └── reports/ # 生成的分析图表运行顺序根据README.md配置Python环境pip install -r requirements.txt和MATLAB工具箱。将竞赛数据放入data/raw/。按编号顺序运行脚本。Python和MATLAB代码相对独立但处理后的数据格式我们设计为可以互通如保存为.csv或.mat文件。6.2 常见报错与解决方案实录在复现过程中你可能会遇到以下问题这里是我们踩坑后的解决方案问题运行LSTM代码时出现内存不足错误。原因构造的序列数据(samples, timesteps, features)维度太大一次性加载进内存导致溢出。解决使用TensorFlow的tf.data.DatasetAPI或生成器函数进行批量化数据加载而不是一次性将全部数据转为NumPy数组。def data_generator(X, y, batch_size32): num_samples len(X) while True: indices np.random.permutation(num_samples) for i in range(0, num_samples, batch_size): batch_idx indices[i:ibatch_size] yield X[batch_idx], y[batch_idx]问题LSTM模型训练损失不下降预测结果是一条直线。原因a数据未进行标准化。LSTM对输入数据的尺度敏感。解决确保对每个特征列进行了标准化如StandardScaler并且必须用训练集的scaler去变换验证集和测试集避免数据泄露。原因b学习率可能太高或网络结构不合适。解决降低学习率如从1e-3调到1e-4尝试更简单的网络如先只用一层LSTM确保激活函数使用正确ReLU或tanh。问题SARIMA模型拟合报错提示“矩阵奇异”或“非平稳”。原因参数(p,d,q)或(P,D,Q,s)选择不当或数据确实不满足平稳性要求。解决首先通过差分d和季节性差分D确保序列平稳。可以使用statsmodels的adfuller函数检验平稳性。然后从低阶参数开始尝试如(0,1,0)或(1,1,1)逐步增加复杂度。问题特征工程后特征数量太多模型训练慢且效果差。解决进行特征筛选。除了之前提到的RFE还可以计算特征与目标变量的皮尔逊相关系数或最大信息系数剔除相关性极低的特征。也可以使用主成分分析进行降维但要注意这会损失特征的可解释性。问题预测结果与历史数据相比波动性方差明显偏小。原因这是机器学习模型特别是使用MSE损失的回归模型的常见问题模型倾向于预测“安全”的平均值低估了极端值。解决在损失函数中增加对极端样本的权重。使用分位数回归直接预测高分位数和低分位数。在特征中引入更多的波动性指标如滚动标准差。这份代码和文档是我们团队心血的结晶它不仅仅是为了解决一道赛题更希望能为你提供一套处理时间序列预测、特别是气候与极端事件预测问题的完整方法论和工具箱。数学建模竞赛中清晰的思路、合理的假设、稳健的模型和令人信服的可视化比追求极致的模型精度更重要。希望你在接下来的比赛中能灵活运用这些工具和方法取得理想的成绩。如果在使用代码中有任何问题欢迎在项目仓库中提出Issue我们会持续维护和更新。