亚太赛C题实战复盘：数据驱动下的全球变暖建模与预测

1. 项目概述：一次完整的亚太赛C题实战复盘

去年带队打完亚太赛，C题“全球是否在变暖？”给我留下了挺深的印象。这题看起来是个老生常谈的话题，但组委会出得相当巧妙，它不是一个简单的“是或否”判断题，而是一个典型的数据驱动型政策分析题。你需要从海量的、多源异构的气候数据里，挖出证据，构建模型，最后还要给出有说服力的评估和建议。很多新手队伍一看到题目里又是温度又是二氧化碳的，容易一头扎进复杂的气候模型里，结果时间耗尽，模型也没跑通。其实，这道题的核心在于数据处理的工程能力和统计建模的逻辑链条，而不是去搞一个能预测百年气候的超级模型。

这篇文章，我就以2022年亚太赛C题为例，完整复盘我们当时的解题思路、数据处理的关键步骤、模型构建的具体实现，以及最后论文写作的要点。我会附上核心的Python代码（基于Pandas、Scikit-learn、Statsmodels等库），并重点分享我们在实战中踩过的坑和总结出的技巧。无论你是正在备赛的数模新手，还是想提升数据分析实战能力的朋友，相信这篇从“战场”上带回来的经验，会比单纯的赛题解析更有参考价值。

2. 赛题核心剖析与解题策略设计

2.1 题目要求拆解：我们到底要回答什么？

拿到题目，第一步不是急着找数据、写代码，而是必须把题目要求一字一句地拆解清楚。2022年C题的要求大致可以归纳为以下几个核心任务：

1. 趋势分析与归因：分析全球及关键区域（如大陆、海洋）的温度变化趋势（1880-2022年），并定量评估主要因素（如CO2浓度、太阳辐射、火山活动等）对变暖的贡献。
2. 预测与情景分析：预测未来几十年（例如到2050或2100年）的温度变化，并分析在不同温室气体排放情景下的可能结果。
3. 影响评估与政策建议：评估全球变暖对自然环境（如海平面、冰川）和人类社会（如农业、健康）的潜在影响，并提出缓解或适应措施。

这里的关键洞察在于，题目虽然叫“全球是否在变暖”，但它的期望答案远不止一个结论。它考察的是你用数据证明趋势的能力、用模型量化关系的能力、用逻辑推演未来的能力以及将复杂科学问题转化为可管理子问题的能力。因此，我们的策略必须覆盖“描述现状、解释原因、预测未来、提出建议”这一完整链条。

2.2 整体解题思路与模型选型

基于以上拆解，我们制定了“数据奠基-统计验证-模型预测-综合评估”的四步走策略。

第一步：数据奠基与预处理。这是整个项目最耗时但也最决定性的环节。我们需要收集权威的长时间序列数据，包括全球平均温度异常、CO2浓度、太阳辐照度、气溶胶光学深度（代表火山活动）等。数据来源如NASA GISS、NOAA、Mauna Loa观测站、CMIP6模型输出等。预处理包括对齐时间分辨率（年均值）、处理缺失值、进行标准化等。

第二步：统计验证趋势与相关性。在构建复杂模型前，先用稳健的统计方法揭示初步规律。我们计划使用：

• Mann-Kendall趋势检验：一种非参数检验，对数据分布没有要求，非常适合检验温度等气候数据是否存在单调上升或下降趋势。
• Sen‘s Slope估计：与M-K检验配套，用于估计趋势的斜率，即变暖的速率。
• 皮尔逊/斯皮尔曼相关系数：分析温度与各潜在驱动因子之间的线性或单调关系强度。

第三步：构建核心预测模型。这是答题的“技术硬核”部分。我们放弃了试图模拟完整地球系统物理过程的复杂模型（时间不允许），选择了多元线性回归和时间序列分析的结合方案。

• 多元线性回归模型：用于归因分析。以温度为因变量，CO2、太阳辐射等为自变量，建立统计关系。通过标准化回归系数比较各因素的贡献度。注意：这里必须考虑变量的共线性（如CO2和其他温室气体可能高度相关），我们计划使用方差膨胀因子(VIF)进行诊断，或采用岭回归(Ridge Regression)来缓解。
• 时间序列模型（ARIMA/Prophet）：用于温度自身的预测。考虑到温度序列的非平稳性和季节性，我们计划使用季节性ARIMA模型。同时，也可以尝试Facebook的Prophet模型，它对于缺失值和趋势变化点处理得比较好，且更易于使用。

第四步：情景分析与综合评估。利用建立好的回归模型，代入IPCC等机构提供的未来不同排放情景（如SSP1-2.6, SSP2-4.5, SSP5-8.5）下的CO2浓度预测数据，从而得到未来温度的可能范围。影响评估部分则主要基于文献调研，将预测的温度上升幅度与已知的研究结论（如每升温1℃海平面上升约X毫米）相结合，进行量化或半量化推断。

策略选择心得：在数模比赛中，模型“足够好”比“理论上最优”更重要。我们的方案没有用深度学习，因为数据量有限且解释性要求高。线性回归和ARIMA虽然传统，但结果稳健、解释性强，而且能在有限时间内实现、调试并写出清晰的论文。这是比赛策略的关键。

3. 数据获取、处理与探索性分析实战

3.1 关键数据源与获取方式

可靠的数据是分析的基石。以下是我们当时使用的主要数据集及其来源，这些源都是公开且权威的：

1. 全球温度数据：来自NASA戈达德空间研究所的GISTEMP数据集。它提供了自1880年以来的全球月度、年度平均地表温度异常（相对于1951-1980年基线）。可以直接从其官网下载CSV或NetCDF格式文件。NOAA的类似数据集也可作为交叉验证。
2. 大气CO2浓度：夏威夷莫纳罗亚观测站的月度平均数据。这是衡量全球背景CO2浓度的黄金标准。数据可以从ESRL网站获取，包含从1958年至今的连续记录。对于更早的年份（1880-1958），需要使用冰芯重建数据（如Law Dome冰芯数据），这些数据通常包含在综合数据集中。
3. 太阳辐照度：使用太阳物理实验室提供的重建的总太阳辐照度序列。这是一个相对缓慢变化的因子。
4. 火山活动指数：使用气溶胶光学深度数据或火山爆发指数。我们采用了来自卫星和地面观测重建的全球平均平流层气溶胶光学深度数据集，它能较好地反映大型火山爆发对全球温度的冷却效应。

实操技巧：建议将所有数据统一处理为年度平均值，并以共同的时间轴（如1880-2022年）进行对齐。对于早期部分缺失的数据（如1958年前的精确CO2），可以采用插值或使用已发表的重建序列。务必在论文中说明数据来源、处理方法和任何假设。

3.2 数据预处理与清洗代码示例

这里用Python的Pandas库演示核心的数据加载与合并步骤。

import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns # 1. 加载温度数据 (假设已下载为CSV，包含‘Year’和‘Temp_Anomaly’列) df_temp = pd.read_csv('gistemp_global_annual.csv') # 可能需要进行一些清洗，比如重命名列、选择时间段 df_temp = df_temp[['Year', 'Temp_Anomaly']].copy() df_temp.set_index('Year', inplace=True) # 2. 加载CO2数据 (MLO数据，需要处理) # MLO数据通常是月度，需要计算年度平均 df_co2_monthly = pd.read_csv('co2_mlo_monthly.csv', comment='#', delim_whitespace=True, names=['year', 'month', 'decimal_date', 'co2_ppm']) # 计算年度平均 df_co2 = df_co2_monthly.groupby('year')['co2_ppm'].mean().reset_index() df_co2.rename(columns={'year':'Year', 'co2_ppm':'CO2_ppm'}, inplace=True) df_co2.set_index('Year', inplace=True) # 3. 加载太阳辐照度数据 (示例) df_solar = pd.read_csv('solar_irradiance_annual.csv', index_col='Year') # 4. 数据合并 # 使用外部连接，确保所有年份都保留，缺失值后续处理 df_merged = pd.concat([df_temp, df_co2, df_solar], axis=1, join='outer') # 重命名列以便清晰 df_merged.columns = ['Temp_Anomaly', 'CO2_ppm', 'Solar_Irradiance'] # 5. 处理缺失值 # 对于CO2早期数据缺失，可以使用前向填充或插值，但更科学的方法是使用冰芯数据填充。 # 这里为演示，对缺失值进行线性插值（需谨慎，实际情况可能更复杂） df_merged['CO2_ppm'] = df_merged['CO2_ppm'].interpolate(method='linear') # 检查是否有剩余缺失值 print(df_merged.isnull().sum()) # 6. 可视化初步查看 fig, axes = plt.subplots(3, 1, figsize=(12, 10)) df_merged['Temp_Anomaly'].plot(ax=axes[0], title='Global Temperature Anomaly (℃)') df_merged['CO2_ppm'].plot(ax=axes[1], title='Atmospheric CO2 Concentration (ppm)', color='orange') df_merged['Solar_Irradiance'].plot(ax=axes[2], title='Total Solar Irradiance (W/m²)', color='red') plt.tight_layout() plt.show()

3.3 探索性分析与统计检验

在建模前，通过可视化和统计检验对数据有一个直观认识至关重要。

from scipy import stats import pymannkendall as mk # 需要安装 pymannkendall 库 # 1. 计算温度与CO2的滚动相关性（例如30年滚动窗口） window_size = 30 df_merged['Rolling_Corr'] = df_merged['Temp_Anomaly'].rolling(window=window_size).corr(df_merged['CO2_ppm']) # 2. Mann-Kendall趋势检验和Sen's Slope估计 # 对温度序列进行检验 result_temp = mk.original_test(df_merged['Temp_Anomaly'].dropna()) print(f"Temperature Trend Test:") print(f" Trend: {result_temp.trend}") print(f" p-value: {result_temp.p:.6f}") print(f" Sen's Slope: {result_temp.slope:.6f} per year") # 对CO2序列进行检验 result_co2 = mk.original_test(df_merged['CO2_ppm'].dropna()) print(f"\nCO2 Trend Test:") print(f" Trend: {result_co2.trend}") print(f" p-value: {result_co2.p:.6f}") # 3. 绘制温度与CO2的散点图与拟合线 plt.figure(figsize=(8,6)) sns.scatterplot(data=df_merged, x='CO2_ppm', y='Temp_Anomaly', alpha=0.6) # 添加线性拟合线 z = np.polyfit(df_merged['CO2_ppm'].dropna(), df_merged['Temp_Anomaly'].dropna(), 1) p = np.poly1d(z) plt.plot(df_merged['CO2_ppm'], p(df_merged['CO2_ppm']), "r--", linewidth=2) plt.xlabel('CO2 Concentration (ppm)') plt.ylabel('Temperature Anomaly (℃)') plt.title('Temperature vs CO2 with Linear Fit') plt.grid(True) plt.show() # 计算相关系数 corr_pearson, p_val_pearson = stats.pearsonr(df_merged['CO2_ppm'].dropna(), df_merged['Temp_Anomaly'].dropna()) print(f"\nPearson Correlation between Temp and CO2: {corr_pearson:.4f} (p={p_val_pearson:.4e})")

数据处理避坑指南：

1. 时间对齐是魔鬼：确保所有数据的时间戳精确对齐到同一年份。有时数据集使用的“年份”可能是观测年份，也可能是发布年份，务必查清。
2. 缺失值处理要谨慎：对于气候数据，尤其是早期数据，简单的前向填充或均值填充可能引入巨大偏差。尽可能使用科学界公认的重建数据来填补空白，或在论文中明确说明处理方式及其潜在局限性。
3. 单位统一：温度异常是相对于基线的变化值，CO2单位是ppm，太阳辐照度是W/m²。确保理解每个数据的物理意义，避免错误比较。
4. 保存中间结果：将清洗合并后的最终数据集保存为新的CSV文件（如climate_data_processed_1880_2022.csv），这样后续所有模型都基于同一份数据，保证可复现性。

4. 核心模型构建：归因分析与预测

4.1 多元线性回归与归因分析

我们使用多元线性回归来量化各因素对温度变化的贡献。假设温度异常T是CO2浓度C、太阳辐照度S和火山气溶胶指数V的线性函数（这里V需要另外获取数据）。

import statsmodels.api as sm from statsmodels.stats.outliers_influence import variance_inflation_factor from sklearn.preprocessing import StandardScaler # 假设 df_merged 已经包含了 'Volcanic_AOD' 列 # 为了评估贡献度，我们通常对自变量进行标准化，使回归系数具有可比性 features = ['CO2_ppm', 'Solar_Irradiance', 'Volcanic_AOD'] X = df_merged[features].dropna() # 确保没有缺失值 y = df_merged.loc[X.index, 'Temp_Anomaly'] # 标准化特征 scaler = StandardScaler() X_scaled = scaler.fit_transform(X) X_scaled = pd.DataFrame(X_scaled, columns=features, index=X.index) # 添加常数项（截距） X_scaled_with_const = sm.add_constant(X_scaled) # 构建OLS模型 model_ols = sm.OLS(y, X_scaled_with_const).fit() # 打印详细的回归结果 print(model_ols.summary()) # 检查多重共线性 - 计算VIF vif_data = pd.DataFrame() vif_data['feature'] = X_scaled_with_const.columns vif_data['VIF'] = [variance_inflation_factor(X_scaled_with_const.values, i) for i in range(X_scaled_with_const.shape[1])] print("\nVariance Inflation Factor (VIF):") print(vif_data) # 解读：标准化后的系数大小反映了该变量对温度变化的相对贡献度。 # 例如，如果CO2的系数为0.65，太阳辐射为0.05，则表明在当前模型框架下， # CO2的变化对温度变化的解释力远强于太阳辐射。

结果解读与注意事项：

• 统计显著性：查看每个系数的P值（P>|t|）。通常P<0.05认为该变量对模型有显著贡献。
• 系数大小：在特征标准化的前提下，系数的绝对值大小直接反映了该因素对温度变化的相对影响强度。这是我们进行归因分析的主要依据。
• 多重共线性警告：如果VIF值大于10（有些严格标准是大于5），说明变量间存在较强的共线性，这会使系数估计不稳定。我们的数据中，CO2和其他温室气体代理变量可能共线性高。解决方案包括：1) 剔除高度相关的变量之一；2) 使用主成分回归(PCR)或岭回归。
• 模型诊断：务必检查残差是否符合正态分布、是否独立（无自相关）。对于时间序列数据，残差自相关是常见问题，会破坏OLS的假设。可以使用Durbin-Watson检验（summary里有）查看，如果DW统计量远偏离2，则存在自相关。

4.2 时间序列预测模型：季节性ARIMA

对于温度自身的时间序列预测，我们采用SARIMA模型。首先需要确定模型的阶数(p,d,q)和季节性阶数(P,D,Q,s)。

from statsmodels.tsa.stattools import adfuller from statsmodels.graphics.tsaplots import plot_acf, plot_pacf from statsmodels.tsa.statespace.sarimax import SARIMAX import warnings warnings.filterwarnings('ignore') # 使用温度异常序列 ts_temp = df_merged['Temp_Anomaly'].dropna() # 1. 平稳性检验 - Augmented Dickey-Fuller Test result_adf = adfuller(ts_temp) print(f'ADF Statistic: {result_adf[0]:.4f}') print(f'p-value: {result_adf[1]:.4f}') if result_adf[1] > 0.05: print("Series is non-stationary. Differencing may be needed.") else: print("Series is stationary.") # 2. 观察ACF和PACF图，初步判断阶数 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12,8)) plot_acf(ts_temp, lags=40, ax=ax1) plot_pacf(ts_temp, lags=40, ax=ax2, method='ywm') plt.tight_layout() plt.show() # 根据ACF衰减慢，判断需要差分；根据PACF在滞后1、2处显著，初步判断p=1或2。 # 3. 尝试差分并再次检验 ts_temp_diff1 = ts_temp.diff().dropna() result_adf_diff1 = adfuller(ts_temp_diff1) print(f'\nAfter 1st differencing, p-value: {result_adf_diff1[1]:.4f}') # 4. 手动尝试或使用网格搜索确定最佳参数（这里演示手动设定一个模型） # 假设我们通过观察和简单尝试，确定阶数为 (1,1,1) 季节性 (1,1,1,10)？ 气候数据的季节性周期是年，但年数据没有季节性，这里用非季节性ARIMA。 # 对于年度数据，通常不考虑季节性，或季节性周期很长。我们使用非季节性ARIMA。 order = (1, 1, 1) # (p,d,q) # 为了演示，我们分割训练集和测试集（例如用2020年之前的数据训练，预测之后） train = ts_temp[ts_temp.index < 2020] test = ts_temp[ts_temp.index >= 2020] model_arima = SARIMAX(train, order=order, trend='c') results_arima = model_arima.fit(disp=False) print(results_arima.summary()) # 5. 进行预测 forecast_steps = len(test) + 20 # 预测测试期+未来20年 forecast_obj = results_arima.get_forecast(steps=forecast_steps) forecast_mean = forecast_obj.predicted_mean forecast_ci = forecast_obj.conf_int() # 6. 绘制结果 plt.figure(figsize=(12,6)) plt.plot(train.index, train, label='Training Data') plt.plot(test.index, test, label='Actual Test Data', color='gray') plt.plot(forecast_mean.index, forecast_mean, label='ARIMA Forecast', color='red') plt.fill_between(forecast_ci.index, forecast_ci.iloc[:, 0], forecast_ci.iloc[:, 1], color='red', alpha=0.2, label='95% Confidence Interval') plt.xlabel('Year') plt.ylabel('Temperature Anomaly (℃)') plt.title('ARIMA Model Forecast for Global Temperature') plt.legend() plt.grid(True) plt.show()

模型调参心得：

1. ARIMA参数选择：p（自回归阶数）和q（移动平均阶数）不宜过大，对于年度数据，(0,1,1)、(1,1,0)、(1,1,1)是常见的起点。可以使用pmdarima库的auto_arima函数进行自动搜索，但在比赛中要写明你的选择依据。
2. 平稳性处理：气候温度序列通常是非平稳的（有上升趋势）。一阶差分（d=1）通常足以使其平稳。务必用ADF检验验证。
3. 预测不确定性：一定要画出置信区间！这能直观展示预测的不确定性，随着预测时间变长，区间会越来越宽，这符合常识，也是论文中的关键信息。
4. 模型诊断图：运行results_arima.plot_diagnostics()来检查残差是否近似白噪声（无自相关、正态分布）。如果残差有问题，说明模型还有改进空间。

5. 情景分析、影响评估与论文写作要点

5.1 基于回归模型的情景预测

我们利用前面建立的多元线性回归模型进行情景预测。关键在于获得未来不同排放情景下自变量的预测值。

# 假设我们已经有一个训练好的回归模型 `model_ols` 和对应的标准化器 `scaler` # 以及未来情景数据框 df_future_scenarios，包含年份和不同情景下的CO2等变量预测值。 # 1. 加载未来情景数据（例如SSP2-4.5情景下的CO2预测） df_future = pd.read_csv('ssp245_co2_projections.csv') # 示例文件，包含'Year', 'CO2_ppm'等列 # 这里需要同样处理太阳辐射和火山活动的未来假设。通常假设太阳辐射遵循已知周期，火山活动随机或取历史平均。 # 2. 使用训练时相同的scaler来标准化未来特征 # 注意：这里必须使用训练集的均值和标准差来转换未来数据，不能重新拟合！ future_features_scaled = scaler.transform(df_future[features]) # 使用之前的scaler future_features_scaled = pd.DataFrame(future_features_scaled, columns=features, index=df_future['Year']) # 3. 添加常数项并预测 future_features_scaled_with_const = sm.add_constant(future_features_scaled, has_constant='add') future_predictions = model_ols.get_prediction(future_features_scaled_with_const) # 获取预测均值及置信区间 future_mean = future_predictions.predicted_mean future_ci = future_predictions.conf_int(alpha=0.05) # 95%置信区间 # 4. 将预测结果整合 df_future['Predicted_Temp_Anomaly'] = future_mean.values df_future['Predicted_Temp_Lower'] = future_ci.iloc[:, 0].values df_future['Predicted_Temp_Upper'] = future_ci.iloc[:, 1].values # 5. 可视化比较不同情景 plt.figure(figsize=(12,6)) # 绘制历史数据 plt.plot(df_merged.index, df_merged['Temp_Anomaly'], label='Historical (Observed)', color='black', linewidth=2) # 绘制预测情景 plt.plot(df_future['Year'], df_future['Predicted_Temp_Anomaly'], label='SSP2-4.5 Projection', color='blue') plt.fill_between(df_future['Year'], df_future['Predicted_Temp_Lower'], df_future['Predicted_Temp_Upper'], color='blue', alpha=0.2, label='95% CI') plt.xlabel('Year') plt.ylabel('Temperature Anomaly (℃)') plt.title('Global Temperature Projection under SSP2-4.5 Scenario') plt.legend() plt.grid(True) plt.show()

5.2 影响评估与政策建议框架

这部分更多是定性或半定量的，基于文献和预测结果进行逻辑推导。

1. 海平面上升：可以引用IPCC报告中的关系，例如全球平均温度每升高1℃，海平面可能上升X米（考虑热膨胀和冰川融化）。将我们的温度预测值代入，估算未来海平面上升范围。
2. 极端天气事件：引用研究指出，温度升高会增加热浪、强降水等极端事件的频率和强度。可以建立简单的统计关系（如温度与某地区极端降水指数的历史关系），并外推。
3. 生态系统与农业：分析温度升高对主要农作物生长季、病虫害的影响。可以使用“积温”等概念进行粗略估算。
4. 政策建议：必须紧扣模型结果。例如，模型显示CO2是主导因子，那么建议就应聚焦于减排。可以量化说明，如果将排放路径从SSP5-8.5（高排放）切换到SSP1-2.6（低排放），到本世纪末可能避免多少度的升温，从而避免多少海平面上升或经济损失。

5.3 数模论文写作核心要点

论文是最终交付物，其清晰度和逻辑性直接决定成绩。

1. 摘要：用一段话概括全部工作。模板：“针对全球变暖问题，本文建立了基于多元线性回归的归因模型和基于时间序列的预测模型。首先，我们收集并处理了1880-2022年的全球温度、CO2浓度等数据，通过Mann-Kendall检验证实了显著的变暖趋势。其次，构建多元线性回归模型，发现CO2浓度变化可解释约XX%的温度变异。进而，建立ARIMA(1,1,1)模型预测未来温度，并在SSP2-4.5情景下，预计2100年全球温度将比工业化前升高X℃。最后，评估了相关影响并提出了减排建议。”
2. 问题重述与分析：用自己的话精炼题目要求，并阐述总体解决思路。
3. 模型假设与符号说明：清晰列出所有重要假设（如“忽略臭氧变化的影响”）和文中使用的符号。
4. 数据分析与预处理：详细展示数据来源、处理步骤、可视化图表。这是体现工作量的重要部分。
5. 模型建立与求解：分小节阐述每个模型（趋势检验、回归、时间序列）。一定要解释为什么选这个模型，参数怎么定的。配上核心代码片段或流程图。
6. 结果分析与讨论：展示所有关键结果图表（趋势图、回归系数表、预测图、置信区间）。对结果进行深入讨论：例如，为什么回归模型中火山活动的系数是负的？预测的不确定性主要来自哪里？
7. 模型评价与推广：客观评价自己模型的优缺点（例如：回归模型简单易解释，但可能忽略了非线性关系和滞后效应；ARIMA模型只依赖历史数据，未考虑未来强制力变化）。提出可能的改进方向。
8. 参考文献：规范引用数据源和关键方法文献。
9. 附录：可以放上完整的、注释良好的核心代码。

最后一点个人体会：数学建模比赛，三分靠建模，七分靠“表达”。这里的表达不仅指论文写作，更指整个解题过程的逻辑叙事能力。你的论文需要像一个引人入胜的故事，从提出问题、分析数据、建立模型、解读结果到给出建议，环环相扣，让评委即使不看细节，也能清晰地跟上你的思路。代码要整洁，图表要美观专业，这些“软实力”往往在队伍水平接近时成为制胜关键。