基于机器学习的重庆市粮食产量预测及影响因素分析*

邬粒，邹黎敏，周科

(1.重庆工商大学数学与统计学院,重庆市,400067; 2.重庆市统计局,重庆市,401147)

0 引言

“民以食为天”,粮食是国家发展的战略物资,是人民生活的必需品,是一切消费品生产的基础。近几年国家多次强调国家粮食安全问题,并出台众多关于粮食和农业发展的政策。解决粮食生产问题、提高粮食产量是保障粮食安全首要关注的事情。重庆市政府积极响应国家号召,为解决国家粮食安全问题作出贡献。虽然重庆境内山高谷深,沟壑纵横,山地面积占76%,使得农用可开发地不多,但是重庆市地处中国西南部、长江上游地区,且重庆属亚热带季风性湿润气候,年降雨量充沛,四季分明,无霜期短,日照时间短,少霜雪,这样的气候资源为重庆的农业提供有利条件。

近年来,已有许多学者通过不同研究方法和研究尺度对国家[1-5]或各省[6-11]的粮食产量影响因素进行了分析,也有学者针对粮食产量预测进行研究[12-16]。王国敏[1]、王双进[2]等采用因子分析法对我国粮食产量及其影响因素进行实证分析,黄坤等[3]从定性和定量分别分析了我国粮食产量的宏观影响因素,江松颖[4]、刘守义[5]等采用柯布-道格拉斯生产函数对我国粮食产量影响因素进行研究,以上学者的研究显示成灾面积、化肥施用量、农业机械总动力和播种面积等因素是影响我国粮食产量的关键因素。江松颖等[4]指出在指标选取时指出劳动力投入、机械总动力等影响因素的数据是对整个农业或农林牧渔业而言的,直接代替为粮食的相关数据不妥,要乘以相应的权重。

在省市层面,各学者采用偏最小二乘模型[6]、主成分分析法[7]、C-D函数[8-9]、滤波技术与分位数回归法[10]、GM模型[11]等方法对不同省市粮食产量影响因素进行分析,其中效灌溉面积、农用机械总动力、粮食播种面积、成灾面积、化肥施用量、劳动力投入等因素是影响各省市粮食产量的重要影响因素。马云倩等[12]采用了LASSO-GM(1,N)、GM(1,1)、GM(1,N)以及LASSO这四种模型对我国粮食产量进行预测,分析发现LASSO-GM(1,N)预测效果更好。陈全润等[13]根据河南省粮食产量的变动特点建立了不同的预测模型,再使用模型平均的方法将不同模型的结果集成,结果表明该方法的预测精度较高。贾梦琦[14]、孟国庆[15]等基于神经网络、灰色预测等方法的组合模型对河北省粮食产量进行预测,都得出组合模型预测效果更好的结论。Nosratabadi等[16]以伊朗的相关数据为例,运用了自适应神经模糊系统(ANFIS)和多层感知器(MLP)的机器学习方法对粮食产量进行预测,在伊朗的数据案例中ANFIS模型比MLP模型预测效果更好。

重庆市粮食产量的提高有助于促进国家粮食安全的发展,也有学者对重庆市粮食产量及其影响因素进行了研究。亓永静等[17]选取了粮食播种面积、化肥使用量、粮食单产、劳动力及成灾面积五个变量为解释变量,依据重庆市1990—2007的数据,建立的线性生产函数,对粮食产量影响因素进行了分析,得到化肥施用量、粮食播种面积、粮食单产是影响粮食产量的重要因素。李岩岩等[18]选取农用机械总动力、农用化肥施用量、农作物播种面积、粮食单产,基于切片逆回归降维方法对重庆市粮食产量进行预测。粮食单产是由粮食产量和粮食播种面积确定后计算的结果,并不能作为粮食总产量的一个影响因素进行分析。同时文献[18]中缺少了成灾面积这一影响粮食产量的重要因素[19-21],重庆市粮食产量易受旱灾和洪涝的影响,例如1998年大洪灾和2006年的旱灾,对重庆市粮食产量有较大的负面影响。

在进行粮食产量影响因素分析时,大多数文献都选择了粮食播种面积、农用机械总动力、成灾面积、有效灌溉面积等指标;在进行粮食产量预测时,近期文献都采用了机器学习及其组合模型的方法,为研究粮食产量影响因素的指标选取和粮食产量预测方法提供了思路。鉴于此,同时考虑到数据的可获取性、指标量化性、主导因素重点凸显性等指标选取原则,本文选取8个指标,以重庆市粮食总产量为目标函数,利用不同核函数的高斯过程回归(GPR)的组合模型对重庆市粮食产量进行预测,并借助于多元函数泰勒公式及偏最小二乘回归(PLSR)对重庆市粮食产量影响因素进行分析。

1 数据采集与数据预处理

1.1 数据采集

本文选取8个指标,分别为粮食播种面积X1(104hm2)、有效灌溉面积X2(104hm2)、农用机械总动力X3(104hm2)、劳动力投入X4(万人)、农村用电量X5(104kW·h)、农用化肥施用量(折纯)X6(104t)、农药使用量X7(104t)、成灾面积X8(104hm2),并选用重庆市1997—2021年数据进行研究,如图1所示,其中数据X1～X8来源于《重庆市统计年鉴》和《中国农业统计年鉴》。

(a) X1

由于统计年鉴中无粮食生产的直接统计数据,以农业的相关指标代替会存在一定误差,因此有效灌溉面积、农用机械总动力、农村用电量、折纯化肥施用量、农药使用量、成灾面积以当年粮食播种面积与农作物播种面积之比为权数[4]。由于统计年鉴中没有粮食及农业劳动力投入数据,但有第一产业的劳动力投入数据,第一产业是指农、林、牧、渔业(不含农、林、牧、渔专业及辅助性活动),因此粮食劳动力投入按如下方式计算。

第一产业就业人员

1.2 数据预处理

根据图2的箱线图,观察各变量是否有离群值。变量X1～X8的值均无离群值,取值均在箱线图的下边缘与上边缘的范围内,而粮食总产量Y存在一个离群值,该离群值为808.4,为重庆市2006年的粮食总产量值。为了研究结果不受离群值的影响,本文删除2006年的相关数据。

图2 各变量的箱线图

从理论上分析,这些因素与粮食产量间不是简单的线性关系,例如农药、化肥使用量与粮食总产量之间的关系,随着农药、化肥使用量的增加,粮食总产量也会增加,但在某个临界点之后,农药、化肥使用量的增加并不会提高粮食产量[22]。

从数据上分析,根据图1中的X1～X8与重庆市粮食总产量Y的关系图可知,各变量与重庆市粮食产量之间并不是简单的线性关系,例如农用机械总动力逐年较稳定地增加,粮食产量呈先下降再上升,后又逐渐下降的趋势,表明农用机械总动力与重庆市粮食产量间不是简单的线性关系。因此,使用非线性模型拟合粮食产量与其影响因素之间的关系更合理。

2 重庆市粮食产量预测分析

2.1 高斯过程回归(GPR)

高斯过程是指所有有限维分布都是多元正态分布的随机过程。高斯过程回归是将回归模型所对应的函数空间视为高斯过程:y=f(X,ω),f～GP(μ,κ)来学习样本,然后估计出回归模型参数的过程,其中μ为均值函数,κ为协方差函数。

现有n组训练集数据X={x1,x2,…,xn},y={y1,y2,…,yn},其中n为训练样本数量。X与y回归的一般模型如式(1)所示。

y=f(X)+ε,ε～N(0,σn2I)

(1)

式中:ε——噪声,且服从独立同分布的0均值正态分布;

σn2——残差方差;

I——单位矩阵。

假设f是一个高斯过程,f(X)～GP[μ,κ(X,X′)],则协方差函数κ(X,X′)公式如式(2)所示。

κ(X,X′)=E[f(X)f(X′)]

(2)

由式(1)和式(2)得到真实输出y与预测值f*的先验分布:y～N(0,K+σn2I),f*～[0,κ(X*,X*)]。

由此得二者的联合高斯先验分布为

(3)

式中:K——协方差矩阵。

高斯过程由其均值函数与协方差函数完全决定,通常选择零均值函数,协方差函数是一个核函数,一般选用径向基核函数

(4)

r=‖xp-xq‖

(5)

式中:xp、xq——数据集中任意两个样本;

δ2——超参数。

对于训练样本,根据式(4)得到最优超参数和协方差矩阵,然后可得与测试样本X*对应的最大概率的预测值,X*的预测值f*的最大概率分布为

P(f*|y,X,X*)=N(f*,μ*,σ*2)

(6)

μ*=K(X*,X)(K+σn2I)-1y

(7)

σ*2=κ(X*,X*)-K(X*,X)(K+σn2I)-1

K(X,X*)

(8)

(9)

(10)

(11)

式中:μ*——预测值f*的均值;

σ*2——预测值f*的方差。

本文使用以下三种核函数的高斯过程回归。

1) Matern核函数

(12)

式中:l、v——核函数的超参数;

Kv——修正贝塞尔函数。

2) 指数核函数

(13)

3) 径向基核函数

(14)

尽管GPR有白噪声分布假设及计算量大的缺点,但GPR有非参数推断灵活、超参数自适应获取等优点,能很好地适应高维数、非线性、小样本等复杂问题,泛化能力强,且在核函数和指数集满足特定条件时,GPR可以近似任意函数[23]。

2.2 GPR模型的预测效果分析

使用留一交叉验证法,分别训练Matern GPR、指数GPR、径向基GPR模型,并得到每一年的粮食产量真实值与预测值之间的绝对百分比误差,如图3所示。

从图3可知,这三种GPR模型的绝对百分比误差大多数都在5%以下,其最大值分别为6.578 7%、7.112 8%、6.717 5%,在10%以下,且通过计算,三种模型的平均绝对百分比误差分别为1.909 9%,1.744 3%,2.013 9%,都在5%以下,根据易丹辉[24]的研究,这三种模型的预测效果都很好。

2.3 组合预测

组合模型可以降低单一模型预测存在的一些风险。第一,单一模型样本假设空间往往很大,会使得结果被误选,降低模型的泛化能力,组合模型能有效减少这一风险。第二,局部极小解不一定是全局极小解,如果在求解算法的最优解时,落入某个局部极小,可能会使得模型的泛化能力很差,这时组合模型能有效地降低这种风险。第三,当前学习算法考虑的空间可能不包含某些学习任务的真实假设,这时使用单一模型无效,而多个模型的组合使得假设空间有所扩大,可能会学习到更好的模型[25]。组合预测就是利用加权平均的形式将不同的单一预测模型组合,然后进行组合模型的预测[25-26]。

2.3.1 方差分析

若要组合预测效果比单一模型更好,则需单一模型的预测能力相当,因此使用方差分析,检验这三种模型的能力是否相当。记Matern GPR、指数GPR、径向基GPR模型的绝对百分比误差的正态总体均值分别为μ1,μ2,μ3,正态总体方差分别表示为σ12,σ22,σ32。

先对以上三种模型的绝对百分比误差两两进行正态总体方差假设检验,原假设为σi2=σj2,备择假设为σi2≠σj2,i,j=1,2,3。在显著性水平为0.05下计算检验统计量

(15)

式中:si2——第i种模型的绝对百分比误差的样本方差;

sj2——第j种模型的绝对百分比误差的样本方差。

F值计算结果如表1所示,查F分布表,总体方差假设检验的拒绝域临界值为F0.025(23,23)=2.311 6。由于表1中的F值均小于2.311 6,故不能拒绝原假设,即任意两个模型的绝对百分比误差的方差无显著性差异。

表1 F值计算结果矩阵

由此可以对三种模型的绝对百分比误差进行方差分析,其原假设和备择假设分别为

H0:μ1=μ2=μ3↔H1:μ1,μ2,μ3不全相等

选择显著水平为α=0.05,得到方差分析的结果如表2所示。可以看出,在0.05的显著性水平下,F值为0.11,小于F0.05(2,66)=3.135 9且P值为0.896 2,大于0.05,故接受原假设,即认为这三种模型的绝对百分比误差的均值无显著性差异。

表2 三种模型的绝对百分比误差方差分析表

2.3.2 组合预测

在个体学习器性能相近时宜使用等权平均法[25],因此本文选用等权平均的方法进行模型的组合。记Matern GPR、指数GPR、径向基GPR模型下的预测值分别为y1、y2、y3,组合预测模型的预测值为ycomb,则

(16)

对这三种模型进行泛化能力的测试,以1997—2019年的数据作为训练集,2020—2021年数据为测试集,得到其泛化结果如表3所示。

表3 2020和2021年预测结果

由表3可知,组合预测模型对2020年和2021年数据预测的绝对百分比误差分别为0.074 4%和0.632 4%,都小于5%,由此说明组合模型的泛化能力非常好,可以作为预测模型,且预测效果比较好。另一方面,虽然组合预测模型的绝对百分比误差略高于指数GPR模型,但是组合预测模型具有更高的鲁棒性[25-26]。

3 重庆市粮食产量影响因素分析

3.1 偏最小二乘回归(PLSR)

使用第2节的非线性模型进行影响因素分析是有困难的。为了使问题简化,通常将目标函数在某点附近展开为泰勒多项式来逼近原函数。多元函数在点X0处的展开式

f(X)=f(X0)+(X-X0)T∇f(X0)+o(ρ)

(17)

其中ρ=‖X-X0‖,因此非线性模型可以近似地表示为线性模型,故可以考虑用线性模型进行影响因素分析。

用条件数方法诊断影响因素之间的多重共线性,计算得条件数为1 996,远大于100,因此变量间存在严重的多重共线性[24],而偏最小二乘回归可以有效地解决这一问题。

3.2 PLSR模型建立

对样本数据进行PLSR建模,计算预测值的绝对百分比误差,并绘制其随主成分个数变化的箱线图,如图4所示。除去个别离群点,其余各个主成分的绝对误差率都集中在2%左右,其中主成分个数为6～8时的分散程度较高,出现了过拟合现象。主成分个数为2时自变量的方差贡献率达到95%以上,绝对百分比误差都集中在4%以下。为防止过拟合,选择2个主成分为最佳主成分个数。

图4 绝对百分比误差的箱线图

使用2个主成分进行偏最小二乘回归,得到如下的回归方程

Y=921.811 1+0.446 5X1-0.031 1X2+

0.169 3X3+0.270 2X4-0.000 2X5+

0.013 4X6-0.003 4X7-0.798 0X8

(18)

使用留一交叉验证方法计算该模型粮食产量真实值与预测值的绝对百分比误差,如图5所示,绝对百分比误差的最大值为2001年的10.101 7%,其余年份的都在10%以下,且除了1997年、2001年及2005年的绝对百分比误差是大于5%的,其余都是小于5%的,平均绝对百分比误差为2.951 7%,根据易丹辉[24]的研究,此线性回归模型也有较好的拟合效果。

将第2章节中的三种GPR模型与PLSR模型的绝对百分比误差进行方差分析,在显著水平α=0.05下,计算得F值为1.36,小于临界值F0.05(3,88)=2.708 2,且P值为0.259 4,大于0.05,故接受原假设,即认为这四种模型的绝对百分比误差的均值无显著性差异。

综合上述分析可知,PLSR模型与三种GPR模型的预测能力大致相当,因此可以使用PLSR模型进行影响因素分析。

3.3 影响因素分析

回归结果表明,在上述的8个影响重庆市粮食产量的因素中,粮食播种面积、农用机械总动力、劳动力投入及受灾面积等因素对重庆市粮食产量具有较大的影响。其中,粮食播种面积每增加一个单位,粮食产量平均增加0.446 5个单位;农用机械总动力每增加一个单位,粮食产量平均增加0.169 3个单位;劳动力投入每增加一个单位,粮食产量平均增加0.270 2个单位;成灾面积每增加一个单位,粮食产量平均减少0.798 0个单位。

3.3.1 粮食播种面积

由图6可以看出,总体上,重庆市粮食总产量和粮食播种面积之间呈现出较强的正相关关系,都呈现下降的趋势。

图6 重庆市粮食产量和粮食播种面积变化趋势

由PLSR方程可知粮食播种面积的回归系数为0.446 5,是对该市粮食产量影响较大的一个因素。重庆市粮食播种面积从1997年以来呈现出逐年递减趋势,由1997年的2.882×106hm2下降至2021年的2.013×106hm2,下降了30.147 3个百分点,且在1997到2007年间的减少幅度较大。这与退耕还林的政策、重庆市城市化水平的提高及农业结构的调整等因素有关。耕地面积逐渐缩小,进而导致粮食播种面积和粮食产量下降,粮食产量也由1997年的1.185×107t下降至2021年的1.092×107t,下降了7.748 4个百分点。重庆市要确保粮食安全,就要有一定粮食播种面积保证粮食产量。

3.3.2 农用机械总动力

由PLSR方程可知农用机械总动力的回归系数为0.169 3,即粮食产量是随农用机械总动力的增加而增加的。但从图7可以看出,从1997年至2021年,重庆市农用机械总动力以一个较稳定的状态增加,1997年至2010年以平均每年5.888 6%的速度增加,2011年至2021年间以平均每年2.135 8%的速度增加,而粮食产量总体上呈下降趋势,即重庆市粮食产量与农用机械总动力间呈负相关关系,出现这一现象的主要原因是粮食播种面积的大幅减少,导致粮食减产。农用机械总动力对粮食产量的正向影响使得粮食产量没有随播种面积的大幅减少而也大幅减产。

图7 重庆市粮食产量和农用机械总动力变化趋势

3.3.3 劳动力投入

劳动力投入也是对重庆市粮食产量比较重要的影响因素,且影响是正向的,劳动力投入每增加1单位,粮食产量就会增加0.270 2个单位。从图8可以看出,劳动力投入与粮食产量总体上都是下降的趋势,在1997—2000年和2010—2014年,重庆市粮食产量与劳动力投入有较强的正相关关系,在2001—2009年和2015—2021年,粮食产量与劳动力投入呈负相关关系。随着重庆市城镇化提高,产业结构调整,且农村里从事农业的收入也比较微薄,越来越多农民放弃从事农业,外出务工,使得粮食以及农业的劳动力逐年减少。

图8 重庆市粮食产量和劳动力投入变化趋势

3.3.4 成灾面积

重庆市粮食产量与成灾面积之间具有一定的负相关关系,受灾面积每增加一个单位,粮食产量平均减少0.798 0个单位,是影响力最大的一个因素。从图9中可以看出,重庆市成灾面积在1997—2021年期间呈现出波动下降的趋势,其中1998年以及2001年成灾面积大,对重庆市粮食产量产生了较为严重的负面影响。1998年,重庆市成灾面积为8.097×105hm2,同比增长75.568 5%;2001年,重庆市成灾面积为6.871×105hm2,同比增长104.971 0%。

图9 重庆市粮食产量和成灾面积变化趋势

1998年的特大洪灾对重庆市造成一定的影响,致使粮食产量同比减少2.470 8%;2001年重庆市遭受春旱、夏旱、特大伏旱和秋旱的袭击,且地下水资源也比较少,造成全市40个区县(自治县、市)普遍受灾,粮食减产9.586×105t,同比减少8.474 1%。可以看出,重庆市粮食产量受自然灾害的影响较大,需要加强自然灾害的防范,尽量降低成灾面积带来的负面影响。

4 结论与建议

4.1 结论

本文选取粮食播种面积、有效灌溉面积、农用机械总动力等指标,基于重庆市1997—2021年度数据,使用三种不同核函数的GPR模型和PLSR模型,对重庆市粮食产量进行预测及影响因素分析。

1) 不同核函数的GPR模型均具有良好的拟合效果和泛化能力,把预测值与实际值进行比较,绝对百分比误差均在10%以下,大多都在5%以下,可作为预测模型进行预测,且其预测能力相当,采用等权平均对三种核函数的GPR模型进行组合预测,能有效提高预测精度。

2) PLSR模型平均绝对百分比误差较低,为2.951 7%,与非线性模型GPR模型的预测能力大致相当,因此可代替非线性模型进行重庆市粮食产量因素分析。通过PLSR模型,发现粮食播种面积、农用机械总动力、劳动力投入以及成灾面积对重庆市粮食产量的影响较大,其中影响最大的因素是成灾面积。粮食播种面积、农用机械总动力和劳动力投入对粮食产量的影响是正向的,成灾面积的影响是负向的。

3) 重庆市粮食产量目前处于低产量区,在整体上呈负增长,粮食播种面积减少、劳动力投入降低、洪涝及干旱灾害对重庆市粮食产量具有严重的负面影响,但农用机械总动力的增加降低了粮食播种面积减少等带来的负面影响。

4.2 建议

1) 保护耕地面积,提高耕地利用率。2021年重庆市人民政府办公厅发布的《关于防止耕地“非粮化”稳定粮食生产的实施意见》明确要求“十四五”时期,重庆粮食播种面积不低于2.003×106hm2,粮食产量不低于1.081×107t。保护耕地面积是提高粮食产量的关键,也是保障粮食安全的关键。严格实行耕地保护政策,才能进一步有效保障粮食播种面积。

2) 农业科技创新,提高粮食产量。农用机械总动力在重庆市粮食产量上的投入还未达到饱和,增加农用机械总动力的投入,能有效降低粮食播种面积下降以及部分自然灾害带来的负面影响。

3) 稳定粮食价格,推广科学种植。保持粮食价格的稳定,出台农业扶持政策,从而保持农民的收入,吸引更多人员回乡就业创业。鼓励乡村人员回乡就业创业的同时也要教授农民科学种植的方法,提高粮食产量。

4) 加强气候监测,及时灾害预警。重庆市气候多变,自然灾害频发,加快完善自然灾害预警预报系统,提高预警精度,能在一定程度上减少自然灾害对农业的影响,减少农业经济损失。维护水利设施,抵御洪涝或旱灾带来的负面影响。