1. 项目缘起与核心价值从一道赛题到一套可复用的分析框架去年带队参加数学建模竞赛队伍抽到了那道关于黄河水沙监测数据的E题。说实话当时看到题目要求对长时间序列的水沙数据进行建模分析和未来预测队伍里几个同学的第一反应是有点懵的。时间序列分析听起来高大上SARIMA模型更是带着一堆令人望而生畏的参数p, d, q, P, D, Q, m。网上能找到的教程要么是纯理论推导看得人云里雾里要么是拿一个完美处理过的AirPassengers数据集跑个Demo代码漂亮但一换自己的数据就各种报错对解决实际问题帮助有限。这道赛题的价值恰恰在于此。它不是一个虚构的完美数据集而是真实、复杂、充满挑战的黄河水沙监测数据。你需要处理的是可能存在的缺失值、明显的季节性波动、非平稳性以及“水”和“沙”这两个变量之间复杂的耦合关系。把这道题吃透你掌握的将不仅仅是如何在Python里调用statsmodels库拟合一个SARIMA模型而是一套应对真实世界时间序列数据的完整分析框架从数据审视、预处理、模型识别、参数调优到诊断预测的全流程。这远比学会一个孤立的模型更有用因为这套方法论可以迁移到销量预测、负荷预测、气象分析等无数个场景中。今天我就把当时我们解题的完整思路、踩过的坑、以及最终稳定可用的Python源码进行一次彻底的复盘和详解。我不会只给你一个“黑箱”代码而是会一步步拆解每个环节背后的“为什么”让你不仅能跑通代码更能理解每个决策背后的逻辑真正掌握时间序列建模的实战能力。2. 赛题数据深度剖析与预处理实战拿到数据的第一步绝对不是急着导入模型。很多新手会直接df pd.read_csv()然后就开始model.fit()这是最大的误区。真实数据尤其是水文监测数据往往“暗藏玄机”。2.1 数据加载与初步审视发现隐藏的问题我们假设数据文件为yellow_river_sediment.csv包含date日期、water_flow流量单位m³/s、sediment_concentration含沙量单位kg/m³等字段。import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns from statsmodels.tsa.stattools import adfuller from statsmodels.graphics.tsaplots import plot_acf, plot_pacf from statsmodels.tsa.seasonal import seasonal_decompose import warnings warnings.filterwarnings(ignore) plt.rcParams[font.sans-serif] [SimHei] # 用来正常显示中文标签 plt.rcParams[axes.unicode_minus] False # 用来正常显示负号 # 加载数据 df pd.read_csv(yellow_river_sediment.csv) df[date] pd.to_datetime(df[date]) df.set_index(date, inplaceTrue) # 初步查看 print(df.head()) print(df.info()) print(df.describe())关键审视点1时间索引与频率。df.info()会告诉你索引是否是DatetimeIndex以及频率freq。水文数据通常是日尺度或月尺度。如果显示freqNone说明数据可能存在缺失的日期。黄河冬季可能结冰某些站点数据会有规律性或非规律性的缺失。关键审视点2异常值与缺失值。df.describe()查看均值、标准差、最小最大值。如果最小流量为0或负数显然不合理或者最大含沙量高得离谱就需要警惕异常值。同时用df.isnull().sum()统计缺失值数量。踩坑记录1天真的填充方式。我们最初用df.fillna(methodffill)前向填充处理缺失值。但对于水文数据特别是流量相邻日期的值可能因降雨、调度发生剧变前向填充会严重扭曲数据的自相关结构导致模型误判。更合理的做法是对于短时间缺失如3-5天可以考虑用线性插值df.interpolate(methodlinear)对于长时间段缺失或明显的异常值如传感器故障可能需要结合历史同期均值或更复杂的方法进行估算甚至在建模前将其剔除并在报告中说明。2.2 序列平稳化检验与处理模型的前提SARIMA模型要求序列是“平稳”的即其统计特性均值、方差不随时间变化。水文序列通常有趋势和季节周期是非平稳的。第一步可视化判断。绘制原始序列图是第一步也是最直观的一步。fig, axes plt.subplots(2, 1, figsize(14, 10)) axes[0].plot(df[water_flow], label水流流量) axes[0].set_title(黄河某站水流流量原始序列) axes[0].legend() axes[0].grid(True) axes[1].plot(df[sediment_concentration], label含沙量, colororange) axes[1].set_title(黄河某站含沙量原始序列) axes[1].legend() axes[1].grid(True) plt.tight_layout() plt.show()从图上你很可能看到1) 长期趋势可能由于气候变化或水利工程2) 年度季节性周期雨季沙多水大旱季相反3) 波动幅度可能随时间变化异方差性。第二步ADF检验定量判断。Augmented Dickey-Fuller检验是一种严格的统计检验。def adf_test(timeseries): print(ADF检验结果:) dftest adfuller(timeseries, autolagAIC) # AIC准则自动选择滞后阶数 dfoutput pd.Series(dftest[0:4], index[ADF统计量, p-value, 滞后阶数, 观测数]) for key, value in dftest[4].items(): dfoutput[f临界值({key})] value print(dfoutput) if dftest[1] 0.05: print(- p值为{:.4f}小于0.05拒绝原假设序列平稳。.format(dftest[1])) else: print(- p值为{:.4f}大于0.05无法拒绝原假设序列非平稳。.format(dftest[1])) print(水流流量ADF检验:) adf_test(df[water_flow].dropna()) print(\n含沙量ADF检验:) adf_test(df[sediment_concentration].dropna())如果p-value 0.05则序列非平稳需要进行差分。第三步差分操作。SARIMA模型中的d和D参数就是用来做差分的。d是非季节性差分阶数D是季节性差分阶数。非季节性差分消除趋势。df[flow_diff1] df[water_flow].diff(1)。通常1阶或2阶差分足以消除趋势。季节性差分消除季节性。若数据是月度数据季节周期m12则季节性差分为df[flow_seasonal_diff] df[water_flow].diff(12)。实操心得不要盲目差分。先做1阶非季节性差分然后再次进行ADF检验和可视化观察序列是否变得“围绕0值上下波动”。如果还有趋势考虑2阶差分。过度差分会导致序列方差增大并引入不必要的相关性损害模型性能。对于同时有趋势和季节性的序列通常的做法是先进行季节性差分如果还不平稳再对季节性差分后的序列进行非季节性差分。在实践中可以借助statsmodels的plot_acf图辅助判断如果原始序列的ACF图衰减很慢说明非平稳差分后ACF图快速衰减到0则说明平稳性改善。3. SARIMA模型核心原理与定阶技巧当你的序列经过预处理变得平稳后就可以进入模型定阶环节。这是SARIMA建模最核心也是最考验经验的部分。3.1 SARIMA模型到底在干什么SARIMA(p,d,q)(P,D,Q)[m]模型可以拆解理解AR(p) 自回归当前值用过去p个时刻的历史值来回归。p越大考虑的历史“记忆”越长。I(d) 差分就是上一节我们做的让序列平稳。MA(q) 移动平均当前值用过去q个时刻的预测误差白噪声来回归。q越大模型对随机冲击的“反应”越灵敏。季节性部分 (P,D,Q)[m]在季节性周期m月度数据m12季度数据m4上重复上述AR、I、MA过程。P是季节性自回归阶数D是季节性差分阶数Q是季节性移动平均阶数。简单说SARIMA就是在同时捕捉序列的短期依赖关系通过p,d,q和季节性周期规律通过P,D,Q,m。3.2 利用ACF/PACF图进行初步定阶这是经典的时间序列定阶方法虽然对于复杂季节模型有时不够精确但能提供非常重要的初始参考。# 假设我们对水流流量进行了1阶非季节性差分和1阶季节性差分m12得到平稳序列 flow_stationary flow_stationary df[water_flow].diff(1).diff(12).dropna() fig, axes plt.subplots(1, 2, figsize(16, 4)) plot_acf(flow_stationary, lags40, axaxes[0]) # 绘制40阶自相关图 plot_pacf(flow_stationary, lags40, axaxes[1], methodywm) # 绘制40阶偏自相关图ywm方法更稳健 plt.show()解读规则经验法则确定非季节性阶数 (p, q):如果ACF图拖尾缓慢衰减PACF图在p阶后截尾p阶后突然降至不显著则建议AR(p)模型。如果PACF图拖尾ACF图在q阶后截尾则建议MA(q)模型。如果两者都拖尾则建议ARMA(p,q)或ARIMA(p,d,q)模型p和q可以尝试从PACF和ACF的显著阶数中选取。确定季节性阶数 (P, Q):观察ACF/PACF在季节周期倍数处如12, 24, 36...是否出现显著峰值。如果ACF在季节滞后处有显著峰值且拖尾PACF在季节滞后处截尾可能P0。如果PACF在季节滞后处有显著峰值且拖尾ACF在季节滞后处截尾可能Q0。重要提醒对于水文这种复杂序列ACF/PACF图可能非常混乱难以清晰判断。这时这个方法主要用来提供p, q, P, Q的大致范围例如p可能取0-3q取0-2为下一步的网格搜索缩小范围。3.3 网格搜索与AIC/BIC准则让数据说话当经验判断困难时最可靠的方法是让计算机通过网格搜索Grid Search来寻找最优参数组合。我们使用**AIC赤池信息准则或BIC贝叶斯信息准则**作为评价标准它们在衡量模型拟合优度的同时惩罚了模型复杂度防止过拟合。AIC倾向于选择更复杂的模型BIC的惩罚更重倾向于更简洁的模型。在样本量足够大时通常更信赖BIC。import itertools from statsmodels.tsa.statespace.sarimax import SARIMAX import warnings warnings.filterwarnings(ignore) # 定义参数搜索范围根据ACF/PACF和经验缩小范围 p d q range(0, 3) # 非季节性参数尝试0,1,2 P D Q range(0, 2) # 季节性参数尝试0,1 m 12 # 月度数据季节周期为12 # 生成所有参数组合 pdq list(itertools.product(p, d, q)) seasonal_pdq list(itertools.product(P, D, Q, [m])) best_aic np.inf best_order None best_seasonal_order None results_list [] print(开始网格搜索...) for param in pdq: for seasonal_param in seasonal_pdq: try: model SARIMAX(df[water_flow], orderparam, seasonal_orderseasonal_param, enforce_stationarityFalse, # 我们已经做了差分这里设为False enforce_invertibilityFalse) result model.fit(dispFalse, maxiter200) # dispFalse不显示迭代日志 current_aic result.aic results_list.append([param, seasonal_param, current_aic]) if current_aic best_aic: best_aic current_aic best_order param best_seasonal_order seasonal_param print(f发现更优模型: ARIMA{param}x{seasonal_param} - AIC:{current_aic:.2f}) except Exception as e: # 某些参数组合可能导致模型无法估计跳过 continue print(f\n搜索完成。) print(f最优模型参数: SARIMA{best_order}x{best_seasonal_order}) print(f最优AIC值: {best_aic:.2f})踩坑记录2网格搜索的代价。参数组合数量是len(p)*len(d)*len(q)*len(P)*len(D)*len(Q)。上面的例子就有33322*2216种组合。每种组合都要拟合一次模型非常耗时。务必在本地或算力充足的机器上运行并尽量根据先验知识缩小搜索范围。可以先固定季节性部分例如先尝试(1,1,1,12)搜索非季节性部分或者先使用auto_arimapmdarima库进行快速初步筛选再在其推荐值附近进行精细网格搜索。4. 模型诊断、预测与结果可视化找到最优参数后工作只完成了一半。我们必须检验这个模型是否真的“好”即它提取了数据的主要模式留下的残差应该是随机的白噪声。4.1 模型拟合与残差诊断# 使用最优参数拟合最终模型 final_model SARIMAX(df[water_flow], orderbest_order, seasonal_orderbest_seasonal_order, enforce_stationarityFalse, enforce_invertibilityFalse) final_result final_model.fit(dispFalse) print(final_result.summary())重点关注summary中的AIC/BIC、Log Likelihood以及各个系数coef的P|z|值。通常P值小于0.05认为该系数显著。残差诊断是必须的步骤# 获取残差 residuals final_result.resid fig, axes plt.subplots(2, 2, figsize(14, 10)) # 1. 残差时序图 axes[0, 0].plot(residuals) axes[0, 0].axhline(y0, colorr, linestyle--) axes[0, 0].set_title(残差序列图) axes[0, 0].set_xlabel(时间) axes[0, 0].set_ylabel(残差) # 2. 残差分布直方图 KDE axes[0, 1].hist(residuals, bins30, edgecolorblack, alpha0.7, densityTrue) residuals.plot(kindkde, axaxes[0, 1], secondary_yTrue, colorred) axes[0, 1].set_title(残差分布) axes[0, 1].set_xlabel(残差值) # 3. 残差Q-Q图检验正态性 from scipy import stats stats.probplot(residuals.dropna(), distnorm, plotaxes[1, 0]) axes[1, 0].set_title(Q-Q图) # 4. 残差自相关图 plot_acf(residuals.dropna(), lags40, axaxes[1, 1]) axes[1, 1].set_title(残差ACF图) plt.tight_layout() plt.show() # 林-博克斯检验Ljung-Box Test检验残差是否为白噪声 from statsmodels.stats.diagnostic import acorr_ljungbox lb_test acorr_ljungbox(residuals.dropna(), lags[10, 20, 30], return_dfTrue) # 检验滞后10,20,30阶 print(\nLjung-Box检验结果 (p-value):) print(lb_test)诊断标准残差序列图应围绕0随机波动无任何明显趋势或周期性。残差分布应近似正态分布钟形曲线。Q-Q图点应大致分布在45度参考线附近。残差ACF图各阶滞后自相关系数应基本落在置信区间内蓝色阴影区域无显著相关。Ljung-Box检验p-value应大于0.05说明无法拒绝“残差是白噪声”的原假设模型通过检验。如果诊断未通过例如残差ACF仍有显著峰值说明模型未能完全捕捉数据中的某些模式可能需要增加p、q、P或Q的阶数或者考虑更复杂的模型如加入外部变量。4.2 样本内拟合与样本外预测# 样本内拟合值动态预测 fitted_values final_result.get_prediction(start0, dynamicFalse).predicted_mean # 未来12步例如未来12个月的预测 forecast_steps 12 forecast_obj final_result.get_forecast(stepsforecast_steps) forecast_mean forecast_obj.predicted_mean forecast_ci forecast_obj.conf_int() # 置信区间 # 可视化 plt.figure(figsize(14, 7)) plt.plot(df[water_flow].index, df[water_flow], label观测值, colorblue, alpha0.7) plt.plot(df[water_flow].index, fitted_values, label样本内拟合, colorred, linestyle--, alpha0.9) plt.plot(pd.date_range(df.index[-1], periodsforecast_steps1, freqM)[1:], forecast_mean, label未来预测, colorgreen, markero) plt.fill_between(pd.date_range(df.index[-1], periodsforecast_steps1, freqM)[1:], forecast_ci.iloc[:, 0], forecast_ci.iloc[:, 1], colorgreen, alpha0.2, label95%置信区间) plt.title(黄河水流流量SARIMA模型拟合与预测) plt.xlabel(日期) plt.ylabel(流量 (m³/s)) plt.legend() plt.grid(True) plt.show() # 打印预测值 print(未来12个月的流量预测值:) print(forecast_mean.round(2))关于动态预测与静态预测dynamicFalse是静态一步向前预测即用直到t-1时刻的所有真实值来预测t时刻。dynamicTrue是动态预测即用模型自身的预测值来预测后续值。在评估样本内拟合时通常用静态预测在做长期样本外预测时只能用动态预测。4.3 针对含沙量序列的建模与思考对于含沙量序列重复上述2-4步流程即可。但需要特别注意两点水沙关系流量和含沙量高度相关。可以考虑建立向量自回归VAR或带外生变量的SARIMAX模型将流量作为预测含沙量的一个外生变量。这通常能显著提升含沙量预测的精度。SARIMAX模型在SARIMAX函数中通过exog参数传入外生变量。序列特性含沙量的波动可能更剧烈异方差性更强。如果残差诊断发现方差随时间变化可能需要考虑对序列进行对数变换np.log1p以稳定方差或者在模型层面使用GARCH族模型来刻画波动聚集效应。这在数学建模中是一个很好的加分点体现了对数据特性的深入思考。5. 竞赛应用延伸与源码整合在数学建模竞赛中仅仅跑出一个模型是远远不够的。你需要将分析过程、模型结果转化为有说服力的论文内容。论文写作要点问题重述与分析清晰定义你要解决的具体预测问题如未来一年月度水沙量预测。数据预处理图文并茂地展示缺失值、异常值处理过程并说明理由。模型建立阐述SARIMA模型的原理重点说明定阶过程展示ACF/PACF图说明网格搜索和AIC准则的选择。模型检验必须包含残差诊断图时序图、ACF图、Q-Q图和Ljung-Box检验结果证明模型的有效性。预测与评估展示拟合与预测图。如果数据允许可以将最后一部分数据留作“测试集”计算均方根误差RMSE、平均绝对百分比误差MAPE等指标来量化预测精度。模型对比与优化可以尝试不同模型如简单指数平滑、Holt-Winters、甚至LSTM进行对比说明SARIMA模型的优势。或者展示加入外生变量SARIMAX后的效果提升。结论与建议基于预测结果给出对黄河水沙调控、生态环境管理的简要建议。完整的、可运行的Python源码框架如下# -*- coding: utf-8 -*- 黄河水沙监测数据分析 - SARIMA模型实战 作者你的名字 日期2023年X月X日 import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns from statsmodels.tsa.stattools import adfuller from statsmodels.graphics.tsaplots import plot_acf, plot_pacf from statsmodels.tsa.seasonal import seasonal_decompose from statsmodels.tsa.statespace.sarimax import SARIMAX from statsmodels.stats.diagnostic import acorr_ljungbox import itertools import warnings warnings.filterwarnings(ignore) # 1. 数据加载与探索 def load_and_explore_data(filepath): df pd.read_csv(filepath) df[date] pd.to_datetime(df[date]) df.set_index(date, inplaceTrue) print(数据概览:) print(df.head()) print(df.info()) print(\n描述性统计:) print(df.describe()) print(\n缺失值统计:) print(df.isnull().sum()) return df # 2. 数据预处理函数 def preprocess_series(series, methodlinear): 处理缺失值默认线性插值 # 处理异常值例如负值置为NaN series[series 0] np.nan # 插值填充 series_filled series.interpolate(methodmethod, limit_directionboth) return series_filled # 3. 平稳性检验与差分函数 def test_stationarity_and_diff(series, seasonal_period12): 检验平稳性并返回差分后的序列及阶数建议 print(原始序列ADF检验:) adf_test(series.dropna()) # 可视化 fig, axes plt.subplots(2, 2, figsize(14, 8)) axes[0, 0].plot(series) axes[0, 0].set_title(原始序列) # 尝试1阶差分 diff1 series.diff(1).dropna() axes[0, 1].plot(diff1) axes[0, 1].set_title(1阶差分后序列) print(\n1阶差分后序列ADF检验:) adf_test(diff1) # 尝试季节性差分 seasonal_diff series.diff(seasonal_period).dropna() axes[1, 0].plot(seasonal_diff) axes[1, 0].set_title(f{seasonal_period}阶季节性差分后序列) print(f\n{seasonal_period}阶季节性差分后序列ADF检验:) adf_test(seasonal_diff) # 尝试1阶季节性差分 diff1_seasonal series.diff(1).diff(seasonal_period).dropna() axes[1, 1].plot(diff1_seasonal) axes[1, 1].set_title(1阶季节性差分后序列) print(f\n1阶{seasonal_period}阶季节性差分后序列ADF检验:) adf_test(diff1_seasonal) plt.tight_layout() plt.show() # 这里可以根据ADF检验p值自动判断但为了教学清晰我们手动选择。 # 通常选择p值显著小于0.05且差分阶数最小的序列。 return diff1_seasonal # 示例返回最平稳的序列 def adf_test(timeseries): 执行ADF检验并打印结果 dftest adfuller(timeseries, autolagAIC) p_value dftest[1] print(fADF统计量: {dftest[0]:.4f}) print(fp-value: {p_value:.4f}) print(临界值:) for key, value in dftest[4].items(): print(f\t{key}: {value:.4f}) if p_value 0.05: print(结论: 序列平稳 (p 0.05)) else: print(结论: 序列非平稳 (p 0.05)) return p_value # 4. 网格搜索最优SARIMA参数 def sarima_grid_search(series, p_range, d_range, q_range, P_range, D_range, Q_range, m, maxiter50): 网格搜索寻找最小AIC的参数组合 best_aic np.inf best_order None best_seasonal_order None all_results [] pdq list(itertools.product(p_range, d_range, q_range)) seasonal_pdq list(itertools.product(P_range, D_range, Q_range, [m])) print(f开始搜索 {len(pdq)*len(seasonal_pdq)} 种参数组合...) for order in pdq: for seasonal_order in seasonal_pdq: try: model SARIMAX(series, orderorder, seasonal_orderseasonal_order, enforce_stationarityFalse, enforce_invertibilityFalse) result model.fit(dispFalse, maxitermaxiter) current_aic result.aic all_results.append([order, seasonal_order, current_aic, result.bic]) if current_aic best_aic: best_aic current_aic best_order order best_seasonal_order seasonal_order print(f 当前最优: SARIMA{order}x{seasonal_order} - AIC:{current_aic:.2f}) except Exception as e: continue result_df pd.DataFrame(all_results, columns[order, seasonal_order, AIC, BIC]).sort_values(AIC) print(\n搜索完成。) print(f最优参数: SARIMA{best_order}x{best_seasonal_order}) print(f最优AIC: {best_aic:.2f}) return best_order, best_seasonal_order, result_df.head(10) # 返回前10个最佳结果 # 5. 模型诊断函数 def model_diagnostics(result): 绘制模型诊断图 residuals result.resid fig, axes plt.subplots(2, 2, figsize(14, 10)) axes[0, 0].plot(residuals) axes[0, 0].axhline(y0, colorr, linestyle--) axes[0, 0].set_title(残差序列) axes[0, 0].set_xlabel(时间) axes[0, 0].set_ylabel(残差) axes[0, 1].hist(residuals, bins30, edgecolorblack, alpha0.7, densityTrue) pd.Series(residuals).plot(kindkde, axaxes[0, 1], colorred) axes[0, 1].set_title(残差分布) axes[0, 1].set_xlabel(残差值) from scipy import stats stats.probplot(residuals.dropna(), distnorm, plotaxes[1, 0]) axes[1, 0].set_title(Q-Q图) plot_acf(residuals.dropna(), lags40, axaxes[1, 1]) axes[1, 1].set_title(残差ACF图) plt.tight_layout() plt.show() # Ljung-Box检验 lb_test acorr_ljungbox(residuals.dropna(), lags[10, 20, 30], return_dfTrue) print(Ljung-Box检验 (p-value):) print(lb_test) if (lb_test[lb_pvalue] 0.05).all(): print(- 残差序列在滞后10,20,30阶均无法拒绝白噪声原假设模型通过检验。) else: print(- 警告残差序列可能存在自相关模型可能未完全捕捉数据模式。) # 6. 预测与绘图函数 def plot_forecast(result, series, forecast_steps12): 绘制历史拟合与未来预测图 # 样本内动态拟合从起始点开始 fitted_values result.get_prediction(start0, dynamicFalse).predicted_mean # 未来预测 forecast_obj result.get_forecast(stepsforecast_steps) forecast_mean forecast_obj.predicted_mean forecast_ci forecast_obj.conf_int(alpha0.05) # 95%置信区间 # 生成未来日期索引 last_date series.index[-1] if isinstance(last_date, pd.Timestamp): freq pd.infer_freq(series.index) future_dates pd.date_range(startlast_date, periodsforecast_steps1, freqfreq)[1:] else: future_dates range(len(series), len(series)forecast_steps) plt.figure(figsize(14, 7)) plt.plot(series.index, series, label观测值, colorblue, alpha0.7, linewidth2) plt.plot(series.index, fitted_values, label样本内拟合, colorred, linestyle--, alpha0.9) plt.plot(future_dates, forecast_mean, labelf未来{forecast_steps}期预测, colorgreen, markero, linewidth2) plt.fill_between(future_dates, forecast_ci.iloc[:, 0], forecast_ci.iloc[:, 1], colorgreen, alpha0.2, label95% 置信区间) plt.title(SARIMA模型拟合与预测结果, fontsize16) plt.xlabel(时间, fontsize12) plt.ylabel(数值, fontsize12) plt.legend(locbest, fontsize10) plt.grid(True, alpha0.3) plt.tight_layout() plt.show() return forecast_mean, forecast_ci # 主程序执行 if __name__ __main__: # 步骤1: 加载数据 data_path yellow_river_sediment.csv # 请替换为你的文件路径 df load_and_explore_data(data_path) # 步骤2: 预处理 - 以水流流量为例 target_series_name water_flow original_series df[target_series_name] processed_series preprocess_series(original_series.copy()) # 步骤3: 平稳性检验与差分 print(f\n 对 {target_series_name} 进行平稳性检验 ) stationary_series test_stationarity_and_diff(processed_series, seasonal_period12) # 假设我们确定使用 1阶非季节性差分 1阶季节性差分 # stationary_series processed_series.diff(1).diff(12).dropna() # 步骤4: 网格搜索 (范围可根据ACF/PACF图调整) print(f\n 开始SARIMA模型网格搜索 ) p_range range(0, 3) # [0,1,2] d_range range(0, 2) # 差分阶数已确定这里d取0或1通常取1 q_range range(0, 3) # [0,1,2] P_range range(0, 2) # [0,1] D_range range(0, 2) # 季节性差分阶数已确定这里D取0或1通常取1 Q_range range(0, 2) # [0,1] m 12 # 月度数据 best_order, best_seasonal_order, top_models sarima_grid_search( processed_series, p_range, d_range, q_range, P_range, D_range, Q_range, m, maxiter100 ) print(f\n最优参数组合为: order{best_order}, seasonal_order{best_seasonal_order}) print(\nAIC最低的前10个模型:) print(top_models) # 步骤5: 用最优参数拟合最终模型 print(f\n 使用最优参数拟合最终模型 ) final_model SARIMAX(processed_series, orderbest_order, seasonal_orderbest_seasonal_order, enforce_stationarityFalse, enforce_invertibilityFalse) final_result final_model.fit(dispTrue, maxiter200) # dispTrue显示迭代信息 print(final_result.summary()) # 步骤6: 模型诊断 print(f\n 模型残差诊断 ) model_diagnostics(final_result) # 步骤7: 预测与绘图 print(f\n 生成预测结果 ) forecast_values, forecast_interval plot_forecast(final_result, processed_series, forecast_steps12) print(f\n未来12期的预测值:) for i, (date, val) in enumerate(zip(forecast_values.index, forecast_values.values), 1): print(f 第{i}期 ({date.date()}): {val:.2f}) # 步骤8: (可选)评估预测精度 - 如果需要测试集 # train_size int(len(processed_series) * 0.8) # train, test processed_series[:train_size], processed_series[train_size:] # ... 在训练集上重新拟合模型并预测测试集进行比较计算RMSE, MAPE等指标。这份源码提供了从数据加载到预测输出的完整管道并包含了详细的注释。你可以直接替换数据文件路径并根据自己数据的特性如季节周期m调整参数搜索范围。对于含沙量序列只需将target_series_name改为sediment_concentration重新运行主程序即可。通过这个项目你收获的不仅仅是一个比赛的解决方案更是一套应对现实世界时间序列问题的组合拳数据敏感度、模型原理理解、调参经验、诊断意识和完整的代码实现能力。这才是数学建模竞赛乃至后续科研或工作中真正宝贵的财富。 从数学建模到实战：基于SARIMA模型的黄河水沙时间序列分析与预测
1. 项目缘起与核心价值从一道赛题到一套可复用的分析框架去年带队参加数学建模竞赛队伍抽到了那道关于黄河水沙监测数据的E题。说实话当时看到题目要求对长时间序列的水沙数据进行建模分析和未来预测队伍里几个同学的第一反应是有点懵的。时间序列分析听起来高大上SARIMA模型更是带着一堆令人望而生畏的参数p, d, q, P, D, Q, m。网上能找到的教程要么是纯理论推导看得人云里雾里要么是拿一个完美处理过的AirPassengers数据集跑个Demo代码漂亮但一换自己的数据就各种报错对解决实际问题帮助有限。这道赛题的价值恰恰在于此。它不是一个虚构的完美数据集而是真实、复杂、充满挑战的黄河水沙监测数据。你需要处理的是可能存在的缺失值、明显的季节性波动、非平稳性以及“水”和“沙”这两个变量之间复杂的耦合关系。把这道题吃透你掌握的将不仅仅是如何在Python里调用statsmodels库拟合一个SARIMA模型而是一套应对真实世界时间序列数据的完整分析框架从数据审视、预处理、模型识别、参数调优到诊断预测的全流程。这远比学会一个孤立的模型更有用因为这套方法论可以迁移到销量预测、负荷预测、气象分析等无数个场景中。今天我就把当时我们解题的完整思路、踩过的坑、以及最终稳定可用的Python源码进行一次彻底的复盘和详解。我不会只给你一个“黑箱”代码而是会一步步拆解每个环节背后的“为什么”让你不仅能跑通代码更能理解每个决策背后的逻辑真正掌握时间序列建模的实战能力。2. 赛题数据深度剖析与预处理实战拿到数据的第一步绝对不是急着导入模型。很多新手会直接df pd.read_csv()然后就开始model.fit()这是最大的误区。真实数据尤其是水文监测数据往往“暗藏玄机”。2.1 数据加载与初步审视发现隐藏的问题我们假设数据文件为yellow_river_sediment.csv包含date日期、water_flow流量单位m³/s、sediment_concentration含沙量单位kg/m³等字段。import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns from statsmodels.tsa.stattools import adfuller from statsmodels.graphics.tsaplots import plot_acf, plot_pacf from statsmodels.tsa.seasonal import seasonal_decompose import warnings warnings.filterwarnings(ignore) plt.rcParams[font.sans-serif] [SimHei] # 用来正常显示中文标签 plt.rcParams[axes.unicode_minus] False # 用来正常显示负号 # 加载数据 df pd.read_csv(yellow_river_sediment.csv) df[date] pd.to_datetime(df[date]) df.set_index(date, inplaceTrue) # 初步查看 print(df.head()) print(df.info()) print(df.describe())关键审视点1时间索引与频率。df.info()会告诉你索引是否是DatetimeIndex以及频率freq。水文数据通常是日尺度或月尺度。如果显示freqNone说明数据可能存在缺失的日期。黄河冬季可能结冰某些站点数据会有规律性或非规律性的缺失。关键审视点2异常值与缺失值。df.describe()查看均值、标准差、最小最大值。如果最小流量为0或负数显然不合理或者最大含沙量高得离谱就需要警惕异常值。同时用df.isnull().sum()统计缺失值数量。踩坑记录1天真的填充方式。我们最初用df.fillna(methodffill)前向填充处理缺失值。但对于水文数据特别是流量相邻日期的值可能因降雨、调度发生剧变前向填充会严重扭曲数据的自相关结构导致模型误判。更合理的做法是对于短时间缺失如3-5天可以考虑用线性插值df.interpolate(methodlinear)对于长时间段缺失或明显的异常值如传感器故障可能需要结合历史同期均值或更复杂的方法进行估算甚至在建模前将其剔除并在报告中说明。2.2 序列平稳化检验与处理模型的前提SARIMA模型要求序列是“平稳”的即其统计特性均值、方差不随时间变化。水文序列通常有趋势和季节周期是非平稳的。第一步可视化判断。绘制原始序列图是第一步也是最直观的一步。fig, axes plt.subplots(2, 1, figsize(14, 10)) axes[0].plot(df[water_flow], label水流流量) axes[0].set_title(黄河某站水流流量原始序列) axes[0].legend() axes[0].grid(True) axes[1].plot(df[sediment_concentration], label含沙量, colororange) axes[1].set_title(黄河某站含沙量原始序列) axes[1].legend() axes[1].grid(True) plt.tight_layout() plt.show()从图上你很可能看到1) 长期趋势可能由于气候变化或水利工程2) 年度季节性周期雨季沙多水大旱季相反3) 波动幅度可能随时间变化异方差性。第二步ADF检验定量判断。Augmented Dickey-Fuller检验是一种严格的统计检验。def adf_test(timeseries): print(ADF检验结果:) dftest adfuller(timeseries, autolagAIC) # AIC准则自动选择滞后阶数 dfoutput pd.Series(dftest[0:4], index[ADF统计量, p-value, 滞后阶数, 观测数]) for key, value in dftest[4].items(): dfoutput[f临界值({key})] value print(dfoutput) if dftest[1] 0.05: print(- p值为{:.4f}小于0.05拒绝原假设序列平稳。.format(dftest[1])) else: print(- p值为{:.4f}大于0.05无法拒绝原假设序列非平稳。.format(dftest[1])) print(水流流量ADF检验:) adf_test(df[water_flow].dropna()) print(\n含沙量ADF检验:) adf_test(df[sediment_concentration].dropna())如果p-value 0.05则序列非平稳需要进行差分。第三步差分操作。SARIMA模型中的d和D参数就是用来做差分的。d是非季节性差分阶数D是季节性差分阶数。非季节性差分消除趋势。df[flow_diff1] df[water_flow].diff(1)。通常1阶或2阶差分足以消除趋势。季节性差分消除季节性。若数据是月度数据季节周期m12则季节性差分为df[flow_seasonal_diff] df[water_flow].diff(12)。实操心得不要盲目差分。先做1阶非季节性差分然后再次进行ADF检验和可视化观察序列是否变得“围绕0值上下波动”。如果还有趋势考虑2阶差分。过度差分会导致序列方差增大并引入不必要的相关性损害模型性能。对于同时有趋势和季节性的序列通常的做法是先进行季节性差分如果还不平稳再对季节性差分后的序列进行非季节性差分。在实践中可以借助statsmodels的plot_acf图辅助判断如果原始序列的ACF图衰减很慢说明非平稳差分后ACF图快速衰减到0则说明平稳性改善。3. SARIMA模型核心原理与定阶技巧当你的序列经过预处理变得平稳后就可以进入模型定阶环节。这是SARIMA建模最核心也是最考验经验的部分。3.1 SARIMA模型到底在干什么SARIMA(p,d,q)(P,D,Q)[m]模型可以拆解理解AR(p) 自回归当前值用过去p个时刻的历史值来回归。p越大考虑的历史“记忆”越长。I(d) 差分就是上一节我们做的让序列平稳。MA(q) 移动平均当前值用过去q个时刻的预测误差白噪声来回归。q越大模型对随机冲击的“反应”越灵敏。季节性部分 (P,D,Q)[m]在季节性周期m月度数据m12季度数据m4上重复上述AR、I、MA过程。P是季节性自回归阶数D是季节性差分阶数Q是季节性移动平均阶数。简单说SARIMA就是在同时捕捉序列的短期依赖关系通过p,d,q和季节性周期规律通过P,D,Q,m。3.2 利用ACF/PACF图进行初步定阶这是经典的时间序列定阶方法虽然对于复杂季节模型有时不够精确但能提供非常重要的初始参考。# 假设我们对水流流量进行了1阶非季节性差分和1阶季节性差分m12得到平稳序列 flow_stationary flow_stationary df[water_flow].diff(1).diff(12).dropna() fig, axes plt.subplots(1, 2, figsize(16, 4)) plot_acf(flow_stationary, lags40, axaxes[0]) # 绘制40阶自相关图 plot_pacf(flow_stationary, lags40, axaxes[1], methodywm) # 绘制40阶偏自相关图ywm方法更稳健 plt.show()解读规则经验法则确定非季节性阶数 (p, q):如果ACF图拖尾缓慢衰减PACF图在p阶后截尾p阶后突然降至不显著则建议AR(p)模型。如果PACF图拖尾ACF图在q阶后截尾则建议MA(q)模型。如果两者都拖尾则建议ARMA(p,q)或ARIMA(p,d,q)模型p和q可以尝试从PACF和ACF的显著阶数中选取。确定季节性阶数 (P, Q):观察ACF/PACF在季节周期倍数处如12, 24, 36...是否出现显著峰值。如果ACF在季节滞后处有显著峰值且拖尾PACF在季节滞后处截尾可能P0。如果PACF在季节滞后处有显著峰值且拖尾ACF在季节滞后处截尾可能Q0。重要提醒对于水文这种复杂序列ACF/PACF图可能非常混乱难以清晰判断。这时这个方法主要用来提供p, q, P, Q的大致范围例如p可能取0-3q取0-2为下一步的网格搜索缩小范围。3.3 网格搜索与AIC/BIC准则让数据说话当经验判断困难时最可靠的方法是让计算机通过网格搜索Grid Search来寻找最优参数组合。我们使用**AIC赤池信息准则或BIC贝叶斯信息准则**作为评价标准它们在衡量模型拟合优度的同时惩罚了模型复杂度防止过拟合。AIC倾向于选择更复杂的模型BIC的惩罚更重倾向于更简洁的模型。在样本量足够大时通常更信赖BIC。import itertools from statsmodels.tsa.statespace.sarimax import SARIMAX import warnings warnings.filterwarnings(ignore) # 定义参数搜索范围根据ACF/PACF和经验缩小范围 p d q range(0, 3) # 非季节性参数尝试0,1,2 P D Q range(0, 2) # 季节性参数尝试0,1 m 12 # 月度数据季节周期为12 # 生成所有参数组合 pdq list(itertools.product(p, d, q)) seasonal_pdq list(itertools.product(P, D, Q, [m])) best_aic np.inf best_order None best_seasonal_order None results_list [] print(开始网格搜索...) for param in pdq: for seasonal_param in seasonal_pdq: try: model SARIMAX(df[water_flow], orderparam, seasonal_orderseasonal_param, enforce_stationarityFalse, # 我们已经做了差分这里设为False enforce_invertibilityFalse) result model.fit(dispFalse, maxiter200) # dispFalse不显示迭代日志 current_aic result.aic results_list.append([param, seasonal_param, current_aic]) if current_aic best_aic: best_aic current_aic best_order param best_seasonal_order seasonal_param print(f发现更优模型: ARIMA{param}x{seasonal_param} - AIC:{current_aic:.2f}) except Exception as e: # 某些参数组合可能导致模型无法估计跳过 continue print(f\n搜索完成。) print(f最优模型参数: SARIMA{best_order}x{best_seasonal_order}) print(f最优AIC值: {best_aic:.2f})踩坑记录2网格搜索的代价。参数组合数量是len(p)*len(d)*len(q)*len(P)*len(D)*len(Q)。上面的例子就有33322*2216种组合。每种组合都要拟合一次模型非常耗时。务必在本地或算力充足的机器上运行并尽量根据先验知识缩小搜索范围。可以先固定季节性部分例如先尝试(1,1,1,12)搜索非季节性部分或者先使用auto_arimapmdarima库进行快速初步筛选再在其推荐值附近进行精细网格搜索。4. 模型诊断、预测与结果可视化找到最优参数后工作只完成了一半。我们必须检验这个模型是否真的“好”即它提取了数据的主要模式留下的残差应该是随机的白噪声。4.1 模型拟合与残差诊断# 使用最优参数拟合最终模型 final_model SARIMAX(df[water_flow], orderbest_order, seasonal_orderbest_seasonal_order, enforce_stationarityFalse, enforce_invertibilityFalse) final_result final_model.fit(dispFalse) print(final_result.summary())重点关注summary中的AIC/BIC、Log Likelihood以及各个系数coef的P|z|值。通常P值小于0.05认为该系数显著。残差诊断是必须的步骤# 获取残差 residuals final_result.resid fig, axes plt.subplots(2, 2, figsize(14, 10)) # 1. 残差时序图 axes[0, 0].plot(residuals) axes[0, 0].axhline(y0, colorr, linestyle--) axes[0, 0].set_title(残差序列图) axes[0, 0].set_xlabel(时间) axes[0, 0].set_ylabel(残差) # 2. 残差分布直方图 KDE axes[0, 1].hist(residuals, bins30, edgecolorblack, alpha0.7, densityTrue) residuals.plot(kindkde, axaxes[0, 1], secondary_yTrue, colorred) axes[0, 1].set_title(残差分布) axes[0, 1].set_xlabel(残差值) # 3. 残差Q-Q图检验正态性 from scipy import stats stats.probplot(residuals.dropna(), distnorm, plotaxes[1, 0]) axes[1, 0].set_title(Q-Q图) # 4. 残差自相关图 plot_acf(residuals.dropna(), lags40, axaxes[1, 1]) axes[1, 1].set_title(残差ACF图) plt.tight_layout() plt.show() # 林-博克斯检验Ljung-Box Test检验残差是否为白噪声 from statsmodels.stats.diagnostic import acorr_ljungbox lb_test acorr_ljungbox(residuals.dropna(), lags[10, 20, 30], return_dfTrue) # 检验滞后10,20,30阶 print(\nLjung-Box检验结果 (p-value):) print(lb_test)诊断标准残差序列图应围绕0随机波动无任何明显趋势或周期性。残差分布应近似正态分布钟形曲线。Q-Q图点应大致分布在45度参考线附近。残差ACF图各阶滞后自相关系数应基本落在置信区间内蓝色阴影区域无显著相关。Ljung-Box检验p-value应大于0.05说明无法拒绝“残差是白噪声”的原假设模型通过检验。如果诊断未通过例如残差ACF仍有显著峰值说明模型未能完全捕捉数据中的某些模式可能需要增加p、q、P或Q的阶数或者考虑更复杂的模型如加入外部变量。4.2 样本内拟合与样本外预测# 样本内拟合值动态预测 fitted_values final_result.get_prediction(start0, dynamicFalse).predicted_mean # 未来12步例如未来12个月的预测 forecast_steps 12 forecast_obj final_result.get_forecast(stepsforecast_steps) forecast_mean forecast_obj.predicted_mean forecast_ci forecast_obj.conf_int() # 置信区间 # 可视化 plt.figure(figsize(14, 7)) plt.plot(df[water_flow].index, df[water_flow], label观测值, colorblue, alpha0.7) plt.plot(df[water_flow].index, fitted_values, label样本内拟合, colorred, linestyle--, alpha0.9) plt.plot(pd.date_range(df.index[-1], periodsforecast_steps1, freqM)[1:], forecast_mean, label未来预测, colorgreen, markero) plt.fill_between(pd.date_range(df.index[-1], periodsforecast_steps1, freqM)[1:], forecast_ci.iloc[:, 0], forecast_ci.iloc[:, 1], colorgreen, alpha0.2, label95%置信区间) plt.title(黄河水流流量SARIMA模型拟合与预测) plt.xlabel(日期) plt.ylabel(流量 (m³/s)) plt.legend() plt.grid(True) plt.show() # 打印预测值 print(未来12个月的流量预测值:) print(forecast_mean.round(2))关于动态预测与静态预测dynamicFalse是静态一步向前预测即用直到t-1时刻的所有真实值来预测t时刻。dynamicTrue是动态预测即用模型自身的预测值来预测后续值。在评估样本内拟合时通常用静态预测在做长期样本外预测时只能用动态预测。4.3 针对含沙量序列的建模与思考对于含沙量序列重复上述2-4步流程即可。但需要特别注意两点水沙关系流量和含沙量高度相关。可以考虑建立向量自回归VAR或带外生变量的SARIMAX模型将流量作为预测含沙量的一个外生变量。这通常能显著提升含沙量预测的精度。SARIMAX模型在SARIMAX函数中通过exog参数传入外生变量。序列特性含沙量的波动可能更剧烈异方差性更强。如果残差诊断发现方差随时间变化可能需要考虑对序列进行对数变换np.log1p以稳定方差或者在模型层面使用GARCH族模型来刻画波动聚集效应。这在数学建模中是一个很好的加分点体现了对数据特性的深入思考。5. 竞赛应用延伸与源码整合在数学建模竞赛中仅仅跑出一个模型是远远不够的。你需要将分析过程、模型结果转化为有说服力的论文内容。论文写作要点问题重述与分析清晰定义你要解决的具体预测问题如未来一年月度水沙量预测。数据预处理图文并茂地展示缺失值、异常值处理过程并说明理由。模型建立阐述SARIMA模型的原理重点说明定阶过程展示ACF/PACF图说明网格搜索和AIC准则的选择。模型检验必须包含残差诊断图时序图、ACF图、Q-Q图和Ljung-Box检验结果证明模型的有效性。预测与评估展示拟合与预测图。如果数据允许可以将最后一部分数据留作“测试集”计算均方根误差RMSE、平均绝对百分比误差MAPE等指标来量化预测精度。模型对比与优化可以尝试不同模型如简单指数平滑、Holt-Winters、甚至LSTM进行对比说明SARIMA模型的优势。或者展示加入外生变量SARIMAX后的效果提升。结论与建议基于预测结果给出对黄河水沙调控、生态环境管理的简要建议。完整的、可运行的Python源码框架如下# -*- coding: utf-8 -*- 黄河水沙监测数据分析 - SARIMA模型实战 作者你的名字 日期2023年X月X日 import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns from statsmodels.tsa.stattools import adfuller from statsmodels.graphics.tsaplots import plot_acf, plot_pacf from statsmodels.tsa.seasonal import seasonal_decompose from statsmodels.tsa.statespace.sarimax import SARIMAX from statsmodels.stats.diagnostic import acorr_ljungbox import itertools import warnings warnings.filterwarnings(ignore) # 1. 数据加载与探索 def load_and_explore_data(filepath): df pd.read_csv(filepath) df[date] pd.to_datetime(df[date]) df.set_index(date, inplaceTrue) print(数据概览:) print(df.head()) print(df.info()) print(\n描述性统计:) print(df.describe()) print(\n缺失值统计:) print(df.isnull().sum()) return df # 2. 数据预处理函数 def preprocess_series(series, methodlinear): 处理缺失值默认线性插值 # 处理异常值例如负值置为NaN series[series 0] np.nan # 插值填充 series_filled series.interpolate(methodmethod, limit_directionboth) return series_filled # 3. 平稳性检验与差分函数 def test_stationarity_and_diff(series, seasonal_period12): 检验平稳性并返回差分后的序列及阶数建议 print(原始序列ADF检验:) adf_test(series.dropna()) # 可视化 fig, axes plt.subplots(2, 2, figsize(14, 8)) axes[0, 0].plot(series) axes[0, 0].set_title(原始序列) # 尝试1阶差分 diff1 series.diff(1).dropna() axes[0, 1].plot(diff1) axes[0, 1].set_title(1阶差分后序列) print(\n1阶差分后序列ADF检验:) adf_test(diff1) # 尝试季节性差分 seasonal_diff series.diff(seasonal_period).dropna() axes[1, 0].plot(seasonal_diff) axes[1, 0].set_title(f{seasonal_period}阶季节性差分后序列) print(f\n{seasonal_period}阶季节性差分后序列ADF检验:) adf_test(seasonal_diff) # 尝试1阶季节性差分 diff1_seasonal series.diff(1).diff(seasonal_period).dropna() axes[1, 1].plot(diff1_seasonal) axes[1, 1].set_title(1阶季节性差分后序列) print(f\n1阶{seasonal_period}阶季节性差分后序列ADF检验:) adf_test(diff1_seasonal) plt.tight_layout() plt.show() # 这里可以根据ADF检验p值自动判断但为了教学清晰我们手动选择。 # 通常选择p值显著小于0.05且差分阶数最小的序列。 return diff1_seasonal # 示例返回最平稳的序列 def adf_test(timeseries): 执行ADF检验并打印结果 dftest adfuller(timeseries, autolagAIC) p_value dftest[1] print(fADF统计量: {dftest[0]:.4f}) print(fp-value: {p_value:.4f}) print(临界值:) for key, value in dftest[4].items(): print(f\t{key}: {value:.4f}) if p_value 0.05: print(结论: 序列平稳 (p 0.05)) else: print(结论: 序列非平稳 (p 0.05)) return p_value # 4. 网格搜索最优SARIMA参数 def sarima_grid_search(series, p_range, d_range, q_range, P_range, D_range, Q_range, m, maxiter50): 网格搜索寻找最小AIC的参数组合 best_aic np.inf best_order None best_seasonal_order None all_results [] pdq list(itertools.product(p_range, d_range, q_range)) seasonal_pdq list(itertools.product(P_range, D_range, Q_range, [m])) print(f开始搜索 {len(pdq)*len(seasonal_pdq)} 种参数组合...) for order in pdq: for seasonal_order in seasonal_pdq: try: model SARIMAX(series, orderorder, seasonal_orderseasonal_order, enforce_stationarityFalse, enforce_invertibilityFalse) result model.fit(dispFalse, maxitermaxiter) current_aic result.aic all_results.append([order, seasonal_order, current_aic, result.bic]) if current_aic best_aic: best_aic current_aic best_order order best_seasonal_order seasonal_order print(f 当前最优: SARIMA{order}x{seasonal_order} - AIC:{current_aic:.2f}) except Exception as e: continue result_df pd.DataFrame(all_results, columns[order, seasonal_order, AIC, BIC]).sort_values(AIC) print(\n搜索完成。) print(f最优参数: SARIMA{best_order}x{best_seasonal_order}) print(f最优AIC: {best_aic:.2f}) return best_order, best_seasonal_order, result_df.head(10) # 返回前10个最佳结果 # 5. 模型诊断函数 def model_diagnostics(result): 绘制模型诊断图 residuals result.resid fig, axes plt.subplots(2, 2, figsize(14, 10)) axes[0, 0].plot(residuals) axes[0, 0].axhline(y0, colorr, linestyle--) axes[0, 0].set_title(残差序列) axes[0, 0].set_xlabel(时间) axes[0, 0].set_ylabel(残差) axes[0, 1].hist(residuals, bins30, edgecolorblack, alpha0.7, densityTrue) pd.Series(residuals).plot(kindkde, axaxes[0, 1], colorred) axes[0, 1].set_title(残差分布) axes[0, 1].set_xlabel(残差值) from scipy import stats stats.probplot(residuals.dropna(), distnorm, plotaxes[1, 0]) axes[1, 0].set_title(Q-Q图) plot_acf(residuals.dropna(), lags40, axaxes[1, 1]) axes[1, 1].set_title(残差ACF图) plt.tight_layout() plt.show() # Ljung-Box检验 lb_test acorr_ljungbox(residuals.dropna(), lags[10, 20, 30], return_dfTrue) print(Ljung-Box检验 (p-value):) print(lb_test) if (lb_test[lb_pvalue] 0.05).all(): print(- 残差序列在滞后10,20,30阶均无法拒绝白噪声原假设模型通过检验。) else: print(- 警告残差序列可能存在自相关模型可能未完全捕捉数据模式。) # 6. 预测与绘图函数 def plot_forecast(result, series, forecast_steps12): 绘制历史拟合与未来预测图 # 样本内动态拟合从起始点开始 fitted_values result.get_prediction(start0, dynamicFalse).predicted_mean # 未来预测 forecast_obj result.get_forecast(stepsforecast_steps) forecast_mean forecast_obj.predicted_mean forecast_ci forecast_obj.conf_int(alpha0.05) # 95%置信区间 # 生成未来日期索引 last_date series.index[-1] if isinstance(last_date, pd.Timestamp): freq pd.infer_freq(series.index) future_dates pd.date_range(startlast_date, periodsforecast_steps1, freqfreq)[1:] else: future_dates range(len(series), len(series)forecast_steps) plt.figure(figsize(14, 7)) plt.plot(series.index, series, label观测值, colorblue, alpha0.7, linewidth2) plt.plot(series.index, fitted_values, label样本内拟合, colorred, linestyle--, alpha0.9) plt.plot(future_dates, forecast_mean, labelf未来{forecast_steps}期预测, colorgreen, markero, linewidth2) plt.fill_between(future_dates, forecast_ci.iloc[:, 0], forecast_ci.iloc[:, 1], colorgreen, alpha0.2, label95% 置信区间) plt.title(SARIMA模型拟合与预测结果, fontsize16) plt.xlabel(时间, fontsize12) plt.ylabel(数值, fontsize12) plt.legend(locbest, fontsize10) plt.grid(True, alpha0.3) plt.tight_layout() plt.show() return forecast_mean, forecast_ci # 主程序执行 if __name__ __main__: # 步骤1: 加载数据 data_path yellow_river_sediment.csv # 请替换为你的文件路径 df load_and_explore_data(data_path) # 步骤2: 预处理 - 以水流流量为例 target_series_name water_flow original_series df[target_series_name] processed_series preprocess_series(original_series.copy()) # 步骤3: 平稳性检验与差分 print(f\n 对 {target_series_name} 进行平稳性检验 ) stationary_series test_stationarity_and_diff(processed_series, seasonal_period12) # 假设我们确定使用 1阶非季节性差分 1阶季节性差分 # stationary_series processed_series.diff(1).diff(12).dropna() # 步骤4: 网格搜索 (范围可根据ACF/PACF图调整) print(f\n 开始SARIMA模型网格搜索 ) p_range range(0, 3) # [0,1,2] d_range range(0, 2) # 差分阶数已确定这里d取0或1通常取1 q_range range(0, 3) # [0,1,2] P_range range(0, 2) # [0,1] D_range range(0, 2) # 季节性差分阶数已确定这里D取0或1通常取1 Q_range range(0, 2) # [0,1] m 12 # 月度数据 best_order, best_seasonal_order, top_models sarima_grid_search( processed_series, p_range, d_range, q_range, P_range, D_range, Q_range, m, maxiter100 ) print(f\n最优参数组合为: order{best_order}, seasonal_order{best_seasonal_order}) print(\nAIC最低的前10个模型:) print(top_models) # 步骤5: 用最优参数拟合最终模型 print(f\n 使用最优参数拟合最终模型 ) final_model SARIMAX(processed_series, orderbest_order, seasonal_orderbest_seasonal_order, enforce_stationarityFalse, enforce_invertibilityFalse) final_result final_model.fit(dispTrue, maxiter200) # dispTrue显示迭代信息 print(final_result.summary()) # 步骤6: 模型诊断 print(f\n 模型残差诊断 ) model_diagnostics(final_result) # 步骤7: 预测与绘图 print(f\n 生成预测结果 ) forecast_values, forecast_interval plot_forecast(final_result, processed_series, forecast_steps12) print(f\n未来12期的预测值:) for i, (date, val) in enumerate(zip(forecast_values.index, forecast_values.values), 1): print(f 第{i}期 ({date.date()}): {val:.2f}) # 步骤8: (可选)评估预测精度 - 如果需要测试集 # train_size int(len(processed_series) * 0.8) # train, test processed_series[:train_size], processed_series[train_size:] # ... 在训练集上重新拟合模型并预测测试集进行比较计算RMSE, MAPE等指标。这份源码提供了从数据加载到预测输出的完整管道并包含了详细的注释。你可以直接替换数据文件路径并根据自己数据的特性如季节周期m调整参数搜索范围。对于含沙量序列只需将target_series_name改为sediment_concentration重新运行主程序即可。通过这个项目你收获的不仅仅是一个比赛的解决方案更是一套应对现实世界时间序列问题的组合拳数据敏感度、模型原理理解、调参经验、诊断意识和完整的代码实现能力。这才是数学建模竞赛乃至后续科研或工作中真正宝贵的财富。