1. 赛题核心与破题思路从“全球变暖”到“量化评估”2022年亚太杯数学建模竞赛的C题题目是“全球是否变暖”。这个题目一出来很多同学的第一反应可能是这不是一个常识性问题吗全球变暖已经是科学界的共识这还用建模但恰恰是这种“常识性”的认知最容易让人掉入陷阱。数学建模竞赛尤其是亚太杯这种级别的赛事其核心从来不是让你去证明一个众所周知的结论而是考察你如何运用数学工具对一个复杂、模糊、充满争议的现实问题进行量化分析、数据挖掘和逻辑论证。这道题的精髓在于“是否”二字背后的不确定性。题目提供的是一系列全球范围内的温度数据但数据本身存在噪声、缺失、时空分布不均等问题。你的任务不是简单地喊一句“是”或“否”而是要通过严谨的数学模型去评估变暖趋势的显著性、量化其速率、分析其空间异质性并探讨其与人类活动如碳排放的关联性。这就像法官判案不能凭感觉必须拿出确凿的证据链。你的模型、算法和代码就是构建这条证据链的工具。因此破题的关键在于转换思维从回答一个“是非题”转变为完成一次“数据驱动的科学评估”。你需要系统地思考几个层面1趋势检测用什么统计方法如线性回归、Mann-Kendall检验从嘈杂的数据中提取出长期趋势2变化速率变暖的速度是多少是匀速、加速还是存在阶段性变化3空间格局变暖是全球均匀的吗哪些区域变暖更剧烈是否存在“变冷”的异常区域4归因分析如何将观测到的变暖趋势与可能的驱动因子如CO2浓度、太阳活动、火山喷发建立关联这每一步都需要具体的模型和代码来实现。2. 数据预处理与特征工程清洗“脏”数据提取有效信号拿到温度数据集通常是NC或CSV格式的栅格数据或站点数据第一步绝不是直接跑模型。原始数据往往“脏”得超乎想象直接使用会导致结果严重失真。数据预处理的质量直接决定了后续所有分析的可靠性。2.1 数据清洗处理缺失值与异常值温度数据常见的“脏”问题包括传感器故障导致的异常高/低值如-9999或9999这样的填充值、数据传输丢失形成的缺失值、以及由于站点迁移或仪器更换造成的非气候性跳变。对于缺失值简单的线性插值在时间序列中往往效果不佳因为温度有显著的季节性和日变化。更稳健的方法是时间序列插值对于短时间缺失如几天可以使用前后时刻的加权平均或者采用更高级的模型如STL分解季节性-趋势分解后对趋势项和季节项分别插值再合成。空间插值对于站点数据如果某些站点在特定时间段缺失可以利用邻近站点的数据通过克里金插值或反距离权重法进行空间插值。这需要你计算站点间的空间相关性。对于异常值不能武断删除需要结合物理意义判断。一个实用的方法是使用滑动窗口统计法计算每个数据点在某个时间窗口如30年内的均值和标准差如果该点偏离均值超过3倍标准差则标记为疑似异常值。然后需要查阅元数据如果提供或结合同期周边站点数据判断是真实极端天气事件如热浪还是仪器错误。如果是后者则按缺失值处理。import pandas as pd import numpy as np def detect_and_clean_anomalies(series, window365*30, n_sigmas3): 使用滚动统计量检测并清理时间序列中的异常值。 # 计算滚动均值和标准差 rolling_mean series.rolling(windowwindow, centerTrue, min_periods1).mean() rolling_std series.rolling(windowwindow, centerTrue, min_periods1).std() # 计算上下界 upper_bound rolling_mean n_sigmas * rolling_std lower_bound rolling_mean - n_sigmas * rolling_std # 标记异常值 anomalies (series upper_bound) | (series lower_bound) # 将异常值替换为NaN视为缺失 cleaned_series series.where(~anomalies, np.nan) # 可选对缺失值进行插值这里用简单线性插值实际可用更复杂方法 cleaned_series_interpolated cleaned_series.interpolate(methodtime) return cleaned_series_interpolated, anomalies # 示例假设temp_data是一个Pandas Series索引为时间 # cleaned_temp, anomaly_flags detect_and_clean_anomalies(temp_data)2.2 均一化处理消除非气候因素影响这是很多新手忽略但至关重要的一步。一个气象站可能在上世纪80年代从市中心迁到郊区或者更换了测温仪器这都会在数据中引入一个“阶跃”式的突变但这并非真实的气候变化。我们需要进行均一化检验与调整。常用的方法有双累积曲线法将目标站点的数据与一个参考站点或参考站点的平均值进行累积和对比如果曲线出现明显转折则说明可能存在不均一性。Pettitt检验一种非参数检验专门用于检测时间序列中的突变点。使用RHtests等专业软件包这是世界气象组织推荐的工具能自动检测并校正多个突变点。在比赛中如果时间有限一个简化的策略是使用区域平均序列作为参考。即计算目标站点所在气候区如华东地区所有站点的平均温度序列将单个站点的序列与之对比并用比值或差值法进行校正以削弱局部非均一性的影响。2.3 空间与时间聚合从点到面从高频到低频原始数据可能是每日甚至每小时的数据空间上可能是数万个网格点。我们需要将其聚合到有分析意义的尺度上。时间聚合通常聚合到年平均值以消除季节循环的干扰突出长期趋势。计算年均值时要确保每年有足够的数据如至少300天的有效数据否则年均值代表性不足。空间聚合如果是栅格数据可以计算全球平均温度序列这是衡量全球变暖最直接的指标。但要注意直接算术平均会高纬度地区的权重因为网格面积随纬度减小更科学的方法是进行面积加权平均每个网格点的值乘以该网格点所代表的地球表面积权重与cos(纬度)成正比。import xarray as xr import numpy as np # 假设已用xarray打开一个NetCDF文件ds # ds[tas] 是近地表空气温度变量维度为 (time, lat, lon) # 1. 计算每个网格点的年平均值 temp_annual ds[tas].resample(timeY).mean(dimtime) # Y代表日历年 # 2. 计算面积权重 lat ds[lat] lon ds[lon] # 将纬度从角度转换为弧度 lat_rad np.deg2rad(lat) # 每个纬度带的权重与其纬度的余弦成正比近似每个网格的面积权重 weight np.cos(lat_rad) # 将权重扩展到与数据相同的维度 (lat, lon) weight_2d np.tile(weight[:, np.newaxis], (1, len(lon))) # 或使用xarray的广播 # 3. 进行全球面积加权平均 # 确保权重和数据在非平均维度上对齐 temp_global_annual (temp_annual * weight_2d).sum(dim(lat, lon)) / weight_2d.sum() # 现在 temp_global_annual 就是一个时间维度的全球平均温度年序列注意在计算全球平均前务必检查数据是否有缺失值如海洋上的缺测并对权重进行相应处理避免缺失值参与计算影响结果。3. 趋势检测与显著性分析用数学说话得到清洗后的温度序列特别是全球平均序列后核心工作就是量化其变化趋势并检验该趋势是否显著。3.1 线性趋势拟合最直观的方法最简单的方法是使用最小二乘法线性回归将温度作为因变量时间如年份作为自变量拟合一条直线。斜率即为变暖速率单位°C/年或°C/10年。import numpy as np import statsmodels.api as sm from scipy import stats # 假设 years 是年份数组如 np.arange(1850, 2023) # global_temp 是对应的全球平均温度序列 X sm.add_constant(years) # 添加常数项截距 model sm.OLS(global_temp, X).fit() slope model.params[1] # 趋势斜率单位°C/年 slope_per_decade slope * 10 # 转换为°C/10年 p_value model.pvalues[1] # 趋势的p值 print(f线性趋势斜率: {slope:.4f} °C/年) print(f相当于 {slope_per_decade:.3f} °C/10年) print(f趋势显著性p值: {p_value:.4f}) if p_value 0.05: # 常用显著性水平0.05 print(在95%置信水平下变暖趋势是统计显著的。) else: print(变暖趋势在统计上不显著。)然而线性回归假设残差独立同分布但气候时间序列常有自相关性今年的温度与去年相关这会低估趋势的不确定性导致p值偏小可能错误地得出“显著”结论。因此必须考虑自相关的影响。3.2 考虑自相关的趋势检验Mann-Kendall与Sen‘s Slope对于存在自相关或非正态分布的数据Mann-Kendall趋势检验结合Sen‘s斜率估计是更稳健的非参数方法。Mann-Kendall检验不假设数据分布通过比较数据点相对顺序来检验趋势是否存在。它给出一个统计量Z和p值。Sen‘s斜率估计计算所有数据点对之间斜率的中位数作为趋势速率的稳健估计。from pymannkendall import original_test # 安装pip install pymannkendall result original_test(global_temp) print(fMann-Kendall趋势检验结果:) print(f 趋势: {result.trend}) # increasing, decreasing, no trend print(f p值: {result.p:.4f}) print(f Sen‘s斜率: {result.slope:.4f} °C/年) print(f 斜率置信区间(95%): ({result.slope_lower:.4f}, {result.slope_upper:.4f}))3.3 分段趋势与突变检测变暖在加速吗线性趋势假设变暖速率恒定。但很多研究指出变暖可能在加速。我们可以使用分段线性回归或滑动窗口趋势分析来探究。滑动窗口趋势定义一个窗口如30年在时间序列上滑动计算每个窗口内的线性趋势。这样可以可视化趋势速率随时间的变化。Piecewise Regression使用算法如pwlf库自动检测数据中的断点拟合不同阶段的不同趋势线。import matplotlib.pyplot as plt window_size 30 # 30年窗口 years np.array(years) global_temp np.array(global_temp) slopes [] mid_years [] for i in range(len(years) - window_size 1): y_sub years[i:iwindow_size] t_sub global_temp[i:iwindow_size] # 简单线性回归 slope, intercept np.polyfit(y_sub, t_sub, 1) slopes.append(slope) mid_years.append(y_sub[window_size//2]) # 窗口中间年份 slopes np.array(slopes) * 10 # 转换为°C/10年 plt.figure(figsize(10, 5)) plt.plot(mid_years, slopes, markero, linestyle-) plt.axhline(y0, colork, linestyle--, alpha0.3) plt.xlabel(年份窗口中心) plt.ylabel(变暖趋势 (°C/10年)) plt.title(f{window_size}年滑动窗口计算的全球变暖速率) plt.grid(True, alpha0.3) plt.show()如果滑动趋势图显示近几十年的斜率明显大于早期则为“加速变暖”提供了证据。4. 空间格局可视化与热点识别变暖并非均匀全球变暖不是“全球均匀变暖”。利用空间数据栅格数据我们可以制作温度趋势的空间分布图直观揭示变暖的“热点”和“冷点”。4.1 计算每个网格点的长期趋势对每个经纬度网格点的时间序列如1850-2021年的年数据进行线性回归或Sen‘s斜率估计得到该点的趋势斜率°C/年形成一个趋势斜率场。import xarray as xr import numpy as np # 假设 ds_annual 是年平均温度数据 xarray.Dataset, 维度 (time, lat, lon) def calc_trend(da): 计算一个数据数组时间序列的线性趋势斜率 time np.arange(len(da.time)) # 使用np.polyfit忽略NaN valid_mask ~np.isnan(da.values) if np.sum(valid_mask) 2: # 至少需要两个有效点 return np.nan slope, _ np.polyfit(time[valid_mask], da.values[valid_mask], 1) return slope # 对每个网格点应用函数 trend_slope xr.apply_ufunc( calc_trend, ds_annual[tas], input_core_dims[[time]], # 对time维度操作 output_core_dims[[]], # 输出是标量 vectorizeTrue ) trend_slope_per_decade trend_slope * 10 # 转换为°C/10年4.2 趋势显著性检验与制图同样对每个网格点进行趋势显著性检验如M-K检验得到p值场。通常将p0.05的区域视为趋势显著在制图时用打点或阴影突出显示。import cartopy.crs as ccrs import cartopy.feature as cfeature import matplotlib.pyplot as plt # 创建地图 fig plt.figure(figsize(14, 6)) ax plt.axes(projectionccrs.Robinson(central_longitude180)) ax.set_global() ax.add_feature(cfeature.LAND, facecolorlightgray) ax.add_feature(cfeature.COASTLINE, linewidth0.5) # 绘制趋势斜率 im ax.pcolormesh(trend_slope_per_decade.lon, trend_slope_per_decade.lat, trend_slope_per_decade, transformccrs.PlateCarree(), cmapRdBu_r, vmin-1, vmax1) # 设置对称色标 plt.colorbar(im, axax, orientationhorizontal, pad0.05, label温度趋势 (°C/10年)) # 叠加显著性打点假设有p_value场 # 假设 p_value 0.05 的区域显著 lon_grid, lat_grid np.meshgrid(trend_slope_per_decade.lon, trend_slope_per_decade.lat) sig_lon lon_grid[p_value 0.05] sig_lat lat_grid[p_value 0.05] ax.scatter(sig_lon, sig_lat, s1, colork, marker., transformccrs.PlateCarree(), alpha0.5) plt.title(全球地表温度变化趋势 (1850-2021) 及显著性 (p0.05), fontsize14) plt.show()通过这张图你可以清晰地论述变暖在高纬度地区特别是北极最为剧烈极地放大效应而某些海洋区域或个别地区趋势不明显甚至略有变冷。这比单纯说“全球变暖”要深刻得多。5. 归因分析与模型构建关联温室气体与温度变化为了增强论证的说服力可以尝试进行简单的归因分析探讨温度变化与人类活动以CO2浓度为代理变量的统计关系。5.1 数据对齐与处理收集全球年平均CO2浓度数据可从NOAA或CMIP6网站获取。将温度序列和CO2浓度序列在时间上对齐如都处理成1850-2021年的年数据。由于两者都有长期上升趋势直接相关可能会受“伪回归”问题干扰即两个趋势性序列即使无关相关性也很高。因此通常对数据进行去趋势或使用一阶差分后再计算相关性。5.2 格兰杰因果检验与回归模型格兰杰因果检验用于检验一个时间序列CO2是否有助于预测另一个时间序列温度。这比简单相关性更能暗示潜在的因果关系方向。可以使用statsmodels中的grangercausalitytests函数。多元线性回归构建一个以温度为因变量CO2浓度、太阳辐射指数、火山活动指数等为自变量的回归模型。通过分析标准化回归系数可以粗略比较各因子的相对贡献。import pandas as pd import statsmodels.api as sm from statsmodels.tsa.stattools import grangercausalitytests # 假设 df 是一个DataFrame包含year, temp_anomaly, co2_ppm列 df df.dropna() # 格兰杰因果检验检验CO2是否格兰杰引起温度变化 # 最大滞后阶数可以根据AIC准则选择这里示例用2 gc_res grangercausalitytests(df[[temp_anomaly, co2_ppm]], maxlag2, verboseFalse) # 查看滞后2阶的结果p值 p_value_granger list(gc_res.values())[0][0][ssr_chi2test][1] print(f格兰杰因果检验 (CO2 - Temp) p值 (lag2): {p_value_granger:.4f}) # 多元线性回归示例仅用CO2 X sm.add_constant(df[co2_ppm]) # 添加截距 y df[temp_anomaly] model sm.OLS(y, X).fit() print(model.summary())重要提示气候系统的归因是极其复杂的科学问题比赛中用到的简单统计模型只能提供初步的、启发性的证据不能作为确凿的因果结论。在论文中必须明确指出这是“统计关联分析”并讨论其局限性如遗漏变量、共同趋势等。6. 建模心得与代码避坑指南结合多次参赛和指导经验分享几个关键的实操心得这些在官方指南里往往不会写。第一数据源的选择与解释至关重要。亚太杯C题通常会提供或暗示数据源如NASA GISTEMP、Berkeley Earth。务必在论文中明确说明你使用的是哪个数据集并简要描述其特点如是否经过均一化处理、空间分辨率等。不同数据集的结果可能有细微差别这是正常的关键在于你能否自圆其说。如果题目没指定推荐使用Berkeley Earth的公开数据它提供了已融合和经过均一化处理的陆地与海洋温度数据且附带不确定性估计非常友好。第二代码的鲁棒性与可复现性是隐形加分项。评委可能不会运行你的代码但代码的结构清晰、注释完整、路径处理得当使用相对路径或明确的数据加载说明能极大体现你的专业素养。务必在代码开头用注释块说明运行环境Python 3.8所需库及版本requirements.txt。处理文件路径时建议使用pathlib或os.path.join避免硬编码绝对路径。# 良好的代码开头示例 2022 APMCM C题 - 全球变暖分析 作者Your_Team_Name 日期2023-XX-XX 环境要求 - Python 3.8 - 依赖库xarray, numpy, pandas, matplotlib, cartopy, pymannkendall, statsmodels 可通过 pip install -r requirements.txt 安装 数据 - 假设温度数据文件为 ./data/global_temperature.nc - CO2数据文件为 ./data/co2_annual.csv from pathlib import Path import xarray as xr import pandas as pd import numpy as np # ... 其他导入 DATA_DIR Path(./data) temp_file DATA_DIR / global_temperature.nc co2_file DATA_DIR / co2_annual.csv # 检查文件是否存在 if not temp_file.exists(): raise FileNotFoundError(f温度数据文件未找到: {temp_file}) # ... 加载数据第三可视化图表要“专业”且“信息丰富”。避免使用默认的艳丽配色。气候科学领域常用如viridis、plasma、RdBu_r用于正负异常等配色。每张图必须有清晰的标题、坐标轴标签带单位、图例。对于空间图务必添加比例尺、指北针或直接使用带经纬度的地图。使用Cartopy库绘制地图是专业性的体现。在图中用文字或箭头标注关键发现如“北极放大效应区域”、“变暖速率最快的海域”。第四不确定性分析是区分优秀与普通论文的关键。不要只给出一个趋势值0.08°C/十年。要给出其置信区间如95% CI: 0.07-0.09°C/十年。在计算全球平均时可以尝试不同的数据源或不同的空间插值方法进行简单的敏感性分析说明你的核心结论如变暖显著在不同处理下是否依然稳健。例如你可以写道“分别采用算术平均和面积加权平均计算全球温度序列两者得出的变暖趋势分别为0.082°C/十年和0.079°C/十年差异在不确定性范围内结论不变。”第五论文论述要围绕模型结果展开避免主观臆断。每一句结论性的表述最好都能对应到前文某张图或某个表格的输出。例如“如图5所示Mann-Kendall检验在99%的置信水平上拒绝了无趋势的原假设p0.01Sen‘s斜率估计显示变暖速率为0.083°C/十年”这样的论述就比“我们认为全球明显变暖”要有力得多。最后时间管理上数据预处理和基础分析趋势、空间图要稳扎稳打确保正确性这构成了论文的基石。归因分析等进阶内容可以作为亮点但如果时间紧张应优先保证基础部分的完整和准确。一个扎实的、没有错误的基础分析远比一个复杂但有漏洞的“高级”模型得分更高。