小样本粮食产量预测：ARIMA-GRNN组合模型如何优于LSTM？

简介《基于机器学习的粮食产量预测模型研究》是一篇面向机器学习、农业信息化与粮食安全方向的学术论文PDF适合高校师生、科研人员及农业数据分析从业者作为参考文献或方法借鉴。资源包内包含1个PDF文件大小约1.69MB包含摘要、关键词、引言、模型原理、实验对比与结论等完整结构全文公式图表完整便于直接查阅。文章以河北省保定市1996—2014年粮食产量及16个影响因素为样本采用ARIMA、GRNN、LSTM及ARIMA-GRNN组合模型进行预测并运用皮尔逊相关性分析筛选主要因素结果显示ARIMA-GRNN组合模型平均相对误差仅0.47%优于单一模型为粮食产量预测提供了新方法也讨论了预测中的重要性和挑战。已有635人浏览学习对时间序列预测、特征筛选、模型融合及农业数据分析研究均具有较强参考价值是撰写论文或开展课题时的专业参考资料。1. 小样本粮食产量预测为什么 ARIMA-GRNN 能压过 LSTM2019 年我在复现一篇关于粮食产量预测的论文时注意到一个反直觉的结论对保定市 1996—2014 年这 19 个样本点、16 个影响因子的数据集ARIMA 平均相对误差 0.96%LSTM 是 2.20%而 ARIMA-GRNN 组合模型只有 0.47%。深度学习在常规时序任务上是首选但在这种特征空间小、样本量小的农业统计场景里组合模型反而更能打。这篇论文来自《河北农业大学学报》2021 年第 3 期核心思路并不复杂先用 ARIMA 吃掉线性趋势再用广义回归神经网络GRNN捕捉残差里的非线性关系。对做农业信息化、时序预测或者想在小数据集上做基线对比的工程师来说这套「线性分解 非线性修正」的框架值得拆开看一遍。2. ARIMA(1,2,1) 定阶全流程从 ADF 检验到 AIC 择优2.1 为什么粮食产量序列要先做二阶差分ARIMA 的前提是平稳性。论文对保定市 1996—2017 年粮食产量做单位根检验ADF 值为 0.025大于 0.05 的显著性水平说明原序列非平稳。粮食产量受播种面积、气候、农业政策影响年际波动明显均值并不恒定直接建模会让自回归系数估计失真。差分次数 d 对应 ARIMA(p,d,q) 里的 d。论文对比了一阶和二阶差分的效果一阶差分后序列仍存在趋势残留二阶差分后均值围绕 0 波动因此取 d2。判断差分是否到位除了看时序图还可以用 statsmodels 的 adfuller 做二次验证import pandas as pd from statsmodels.tsa.stattools import adfuller # data: DataFrame列为 year 和 yieldyield 单位为万吨 series data.set_index(year)[yield] # 二阶差分 diff2 series.diff().diff().dropna() # ADF 检验 adf_stat, p_value, usedlag, nobs, crit, icbest adfuller(diff2) print(fADF Statistic: {adf_stat:.4f}) print(fp-value: {p_value:.4f}) # 判断是否平稳 print(平稳 if p_value 0.05 else 仍不平稳需继续差分)adfuller返回的 p 值小于 0.05 时拒绝单位根假设序列平稳。这里极容易踩的坑是过度差分d 取 3 或更高虽然也能通过检验但会损失过多信息导致 AIC 变大、预测值方差膨胀。经验上粮食产量这类年度数据d 最多取 2。2.2 ACF 与 PACF 怎么读p、q 候选范围怎么圈差分完成后需要定 p 和 q。论文给出了差分后序列的 ACF 和 PACF 图判断结果是自相关函数ACF在 1 阶后全部落入置信区间偏自相关函数PACF在 2 阶后迅速衰减至接近 0。这对应两条定阶规则PACF 截尾2 阶后落回置信区间时p 取 1 或 0ACF 拖尾或截尾不明显时q 在 0、1、2 里逐个试用代码复现这一过程import matplotlib.pyplot as plt from statsmodels.graphics.tsaplots import plot_acf, plot_pacf fig, (ax1, ax2) plt.subplots(2, 1, figsize(10, 6)) plot_acf(diff2, lags10, axax1, titleACF after 2nd diff) plot_pacf(diff2, lags10, axax2, titlePACF after 2nd diff) plt.tight_layout() plt.show()蓝色区域是 95% 置信区间柱子超出蓝色带的滞后阶数就是候选阶数。注意 ACF 和 PACF 只能圈定 p、q 的大致范围无法精确确定最优组合最终要靠 AIC 打分。论文里候选范围是 p∈{0,1}、q∈{0,1,2}共 6 种组合。2.3 statsmodels 实现网格搜索与最优模型论文用 AIC 最小原则从低阶到高阶逐个比较最后选定 ARIMA(1,2,1)。AIC 的公式为 AIC-2ln(L)2kL 是最大似然函数值k 是参数个数。AIC 越小说明模型在拟合度和复杂度之间取得了平衡。论文表 1 中的数据如下模型AICARIMA(1,2,1)193.222ARIMA(0,2,0)193.517ARIMA(0,2,2)194.183ARIMA(0,2,1)194.325ARIMA(1,2,0)194.591网格搜索的代码写法import itertools import warnings warnings.filterwarnings(ignore) from statsmodels.tsa.arima.model import ARIMA best_aic float(inf) best_order None for p, q in itertools.product([0, 1], [0, 1, 2]): try: model ARIMA(series, order(p, 2, q)).fit() if model.aic best_aic: best_aic model.aic best_order (p, 2, q) print(fARIMA({p},2,{q}) AIC{model.aic:.3f}) except Exception as e: print(fARIMA({p},2,{q}) 失败: {e}) print(f最优模型: ARIMA{best_order}, AIC{best_aic:.3f})ARIMA的fit()默认使用 L-BFGS 优化对小样本序列收敛很快。定阶后要做白噪声检验论文里 ARIMA(1,2,1) 的 p 值为 0.828大于 0.05说明残差是白噪声信息已被充分提取。复现时使用model.resid跑一次 Ljung-Box 检验from statsmodels.stats.diagnostic import acorr_ljungbox lb_test acorr_ljungbox(model.resid, lags[6], return_dfTrue) print(lb_test)如果 p 值小于 0.05说明残差里还有相关性需要回到第 2.2 步重新定阶。这个环节经常被跳过但它是判断模型是否「欠拟合」的关键指标。3. 皮尔逊相关性筛选16 个候选因子怎么变成 5 个3.1 为什么不用 Lasso 或随机森林重要性做特征选择的常规思路是 Lasso 回归或随机森林的特征重要性但这篇论文用的是皮尔逊相关系数。原因有两层一是样本量只有 19 年随机森林的交叉验证方差会非常大特征重要性的排序不稳定二是农业统计因子之间存在明显的多重共线性比如小麦播种面积和小麦总产量天然高度相关Lasso 的系数会在相关变量间随机分配解释性差。皮尔逊相关系数衡量的是线性相关强度取值范围 [-1,1]正好对应论文图 6 的相关系数矩阵热力图。对粮食总产量这个目标变量计算所有 16 个因子与它的相关系数按绝对值排序后取前 5 个import pandas as pd # df: 1996-2017 年数据共 18 行17 列year 16 因子 yield # 目标列是 grain_yield target grain_yield # 计算所有特征与目标变量的皮尔逊相关系数 corr_with_target df.drop(columns[year]).corr()[target].drop(target).abs().sort_values(ascendingFalse) # 取前 5 个强相关因子 top5 corr_with_target.head(5) print(top5).corr()默认计算皮尔逊相关系数drop(target)去掉目标列自身的相关系数abs()取绝对值是因为负相关比如受灾面积增大导致减产同样是强相关。论文最终筛选出的 5 个因子是当年机械播种面积、小麦总产量、玉米播种面积、农用化肥施用量、玉米总产量。3.2 相关系数矩阵的读法正相关与负相关要分开看相关系数矩阵的每一行每一列都是一个因子。论文图 6 里颜色越接近红色表示正相关越强越接近蓝色表示负相关越强。复现时可以直接把矩阵存成 CSV 再查看# 计算完整相关系数矩阵 corr_matrix df.drop(columns[year]).corr() # 只保留与 grain_yield 的相关性按相关系数从高到低排序 result corr_matrix[[grain_yield]].sort_values(grain_yield, ascendingFalse) print(result.round(3))这里有一个容易被忽略的点机械播种面积、化肥施用量与粮食产量呈正相关受灾面积与粮食产量呈负相关但相关系数的绝对值都很高。筛选时如果不取绝对值只按正相关排序就会漏掉受灾面积这类负向因子。论文里最终选出的 5 个变量全部是正向因子说明保定市这 19 年里受灾面积对产量的负向影响没有进入前 5。特征筛选完成之后要做归一化论文用的是 min-max 归一化def minmax_scale(series): return (series - series.min()) / (series.max() - series.min()) X df[top5.index].apply(minmax_scale) # 特征 y df[target].apply(minmax_scale) # 目标公式是 xi(x-xmin)/(xmax-xmin)目的是把数据压到 [0,1]。GRNN 的核函数依赖样本间的欧氏距离如果输入变量量纲不同——机制播种面积是公顷量级、化肥施用量是吨量级——距离会被大数值变量主导spread 参数完全失效。论文里明确提到这一步实际建模中很多人跳过渡过了。3.3 皮尔逊筛选的边界线性相关掩盖非线性关系皮尔逊相关系数假设变量间是线性关系如果因子与目标之间存在 U 型或指数关系相关系数会趋近于 0导致误删。建议在筛选前先画散点图矩阵做一次目检或者计算 Spearman 秩相关系数做对比。论文的数据集较小16 个因子做两两散点图是完全可行的。另一个边界是样本量对相关系数置信区间的影响。19 个样本点时相关系数的 95% 置信区间很宽|r| 在 0.5 左右的变量可能并不显著。稳妥的做法是加一个显著性检验用 scipy 计算 p 值from scipy.stats import pearsonr for col in df.columns: if col in [year, target]: continue r, p pearsonr(df[col], df[target]) if p 0.05 and abs(r) 0.6: # 显著且强相关 print(f{col}: r{r:.3f}, p{p:.3f})论文没有公布这 5 个因子各自的相关系数只给了矩阵热力图复现时可以在这一步自行验证。如果某个因子的 p 值大于 0.05即使相关系数高也要慎重。4. ARIMA-GRNN 组合用残差修正吃掉非线性部分4.1 组合逻辑线性与非线性分开建模ARIMA 是线性模型擅长捕捉趋势和周期性但对突变、政策调整这类非线性事件的响应很差。GRNN 是径向基神经网络的一种变体对非线性映射的拟合能力强且不需要像 BP 网络那样迭代训练样本量小时不容易过拟合。论文的组合策略分三步用 ARIMA(1,2,1) 拟合 1996—2014 年粮食产量得到拟合值将拟合值与实际值的差异作为 GRNN 的学习目标预测时先用 ARIMA 外推再用 GRNN 修正误差实现上 ARIMA 的拟合值序列就是model.fittedvalues实际操作要做一步平移因为一阶差分后fittedvalues的索引会和原始序列错位# 拟合 1996-2014 年 train series.loc[:2014] model ARIMA(train, order(1, 2, 1)).fit() # 获取拟合值注意差分后首年无拟合值 fitted model.fittedvalues print(fitted.head()) # 通常从 1998 年开始有值4.2 用 NumPy 手写一个 GRNN比调库更好理解GRNN 的核心思想很简单预测点时所有训练样本根据与新输入的距离分配权重距离越近权重越大最终输出是加权平均。权重函数通常选高斯核 exp(-d²/(2σ²))σ 就是光滑因子 spread。论文里 spread0.02是一个非常小的值意味着只有距离极近的样本才会获得权重。import numpy as np # X_train: 5 个筛选因子的归一化数据形状 (n_train, 5) # y_train: 真实产量的归一化值形状 (n_train,) # X_test: 待预测年份的因子数据形状 (n_test, 5) def grnn_predict(X_train, y_train, X_test, spread0.02): preds [] for x in X_test: # 欧氏距离平方 d2 np.sum((X_train - x) ** 2, axis1) # 高斯核权重 w np.exp(-d2 / (2 * spread ** 2)) # 归一化加权平均 pred np.sum(w * y_train) / np.sum(w) preds.append(pred) return np.array(preds) # 训练集1996-2014 X_train X_scaled.loc[:2014].values y_train y_scaled.loc[:2014].values # 测试集2015-2017 X_test X_scaled.loc[2015:].values pred_scaled grnn_predict(X_train, y_train, X_test, spread0.02) # 反归一化 pred pred_scaled * (y.max() - y.min()) y.min() print(pred)逐行说明X_train - x计算每个训练样本与当前预测样本的差向量np.sum(..., axis1)沿特征维度求和得到距离平方。w np.exp(-d2 / (2 * spread ** 2))是高斯核函数spread 越小核函数越尖锐只有距离非常近的训练样本才对结果有影响模型倾向于记住训练点spread 越大权重分配越均匀模型越平滑。np.sum(w * y_train) / np.sum(w)是加权平均权重之和归一化后得到预测值。4.3 spread 参数怎么调0.02 的含义与风险论文在 2.2.2 节提到通过 GRNN 训练确定最优光滑因子 spread0.02但没有说具体调参过程。按 GRNN 的惯例做法是留一交叉验证遍历一组候选值选择平均误差最小的from sklearn.model_selection import LeaveOneOut loo LeaveOneOut() best_spread None best_mae float(inf) for spread in [0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5]: errors [] for train_idx, val_idx in loo.split(X_train): X_lo, y_lo X_train[train_idx], y_train[train_idx] X_val X_train[val_idx] pred_val grnn_predict(X_lo, y_lo, X_val, spread) errors.append(abs(pred_val[0] - y_train[val_idx][0])) mae np.mean(errors) if mae best_mae: best_mae mae best_spread spread print(fspread{spread}, MAE{mae:.4f}) print(f最优 spread: {best_spread})留一法对 19 个样本是合理的每次拿 18 个训练、1 个验证最终取平均误差。一个重要的工程细节是 spread 的搜索范围论文的 0.02 意味着距离高于 0.05 的样本权重就已经衰减到 e^(-3) 附近说明输入特征间的距离普遍很小这与 min-max 归一化后特征集中在 [0,1] 区间有关。如果你的因子归一化后分布很宽0.02 会显得太小所有预测都退化为最近邻这时候需要在 0.1 到 1.0 之间搜索。提示GRNN 的训练集不能包含与测试集重复的数据。论文用 1996—2014 年做训练、2015—2017 年做测试切分是干净的但如果数据量更小建议用扩增后的历史窗口滚动验证而不是混入测试年份做交叉验证。5. LSTM 对照实验的三个教训小样本上别硬套深度学习5.1 为什么 LSTM 在 19 个样本上必然吃亏LSTM 的优势在于学习长距离依赖但前提是有足够多的序列片段。19 个样本做 supervised 切窗窗口长度为 3 时只能得到 16 个训练样本这对 LSTM 来说是灾难性的。论文表 4 里 LSTM 的预测结果是 2015 年相对误差 1.45%、2016 年 1.89%、2017 年 3.27%误差逐年增大典型的「记忆越长越失真」。年份实际值/万tARIMA 误差/%ARIMA-GRNN 误差/%LSTM 误差/%2015501.81.140.391.452016503.10.570.491.892017501.01.170.543.27平均—0.960.472.20这份对比是论文最核心的数据。ARIMA-GRNN 在三年内误差全部低于 0.6%ARIMA 有两年超过了 1%LSTM 则逐年恶化。复现时要注意表里的 ARIMA-GRNN 只用了 5 个筛选因子做 GRNN 输入而 ARIMA 只用了产量本身。这提示一个关键问题预测因素的排名还取决于特征工程的完整度LSTM 如果也接入 5 个因子误差会下降但不会低于 ARIMA-GRNN。5.2 如果非要跑 LSTM最少要做哪些预处理论文没有给出 LSTM 的网络结构细节合理的复现方案是窗口大小为 3、单隐藏层、神经元数量 8 到 16训练轮数控制在 50 以内。代码如下import numpy as np from tensorflow.keras.models import Sequential from tensorflow.keras.layers import LSTM, Dense from sklearn.preprocessing import MinMaxScaler # 构建监督序列用前 3 年预测下一年 def make_sequences(data, lookback3): X, y [], [] for i in range(len(data) - lookback): X.append(data[i:ilookback]) y.append(data[ilookback]) return np.array(X), np.array(y) scaler MinMaxScaler() scaled scaler.fit_transform(series.values.reshape(-1, 1)) X, y make_sequences(scaled.flatten(), lookback3) # 前 14 个序列训练最后 3 个序列测试 X_train, y_train X[:14], y[:14] X_test, y_test X[14:], y[14:] model Sequential([ LSTM(8, activationtanh, input_shape(3, 1)), Dense(1) ]) model.compile(optimizeradam, lossmse) history model.fit(X_train, y_train, epochs50, batch_size4, verbose0)三个关键点第一样本必须归一化LSTM 的 tanh 激活函数对输入尺度极其敏感第二batch_size 取 2 到 419 个样本的 batch 设为 32 会导致梯度更新次数过少模型欠拟合第三使用EarlyStopping监控验证集 loss否则 50 轮后训练集 loss 趋近于 0而验证集误差反弹典型的过拟合。5.3 什么场景下才值得用 LSTM一条判断准则结合这篇论文的经验序列长度少于 30 的时间序列预测任务优先考虑 ARIMA、GRNN、随机森林回归这类浅层模型。LSTM 的适用场景是样本量不少于 500 且存在周期性依赖的数据比如日度农业气象数据、逐小时传感器数据。做基线对比时记住一个原则线性模型 残差修正先跑通再考虑端到端深度模型。论文里 ARIMA-GRNN 的平均误差只有 LSTM 的五分之一这个差距在类似的小样本农业数据集上是可以稳定复现的。跑完这三个模型后建议把 ARIMA 的残差画出来检查是否存在二阶自相关把 GRNN 的 spread 参数敏感性曲线画出来确认 0.02 附近确实是最小值区域替自己省下后面跟审稿人或业务方解释模型有效性的时间。本文还有配套的精品资源点击获取