基于改进主成分分析的LSTM水稻产量预测模型研究

Research on LSTM Rice Yield Prediction Model Based on Improved Principal Component Analysis

摘要: 作为农业大国，农业一直是我国经济发展的重要基石，而水稻作为我国主要的粮食作物，随着人口不断增长，其需求量也不断攀升，因此水稻的产量预测对我国农业的发展建设以及保障粮食安全具有重要意义。长短期记忆(LSTM)循环神经网络因其不仅能够较好地处理各因素间的非线性关系，且适合处理时间序列数据的预测问题，在作物产量预测领域应用前景良好。本文提出一种基于改进主成分分析(IPCA)的LSTM循环神经网络，对神经网络的输入进行数据降维，旨在提高神经网络训练的收敛速度，并消除输入数据间由于相关性导致的信息冗余，从而提高预测精度。

Abstract: As a major agricultural country, agriculture has always been an important cornerstone of our country’s economic development. As the main food crop in our country, with the continuous growth of the population, the demand for rice is also rising. Therefore, the prediction of rice production is very important for the development and construction of agriculture and guarantee of food security. Long short-term memory (LSTM) recurrent neural network has a good application prospect in the field of crop yield prediction because it can not only deal with the nonlinear relationship between various factors, but also is suitable for dealing with the prediction problem of time series data. In this paper, a LSTM recurrent neural network based on improved principal component analysis is proposed, which reduces the data dimension of the input of the neural network, aiming to improve the convergence speed of neural network and eliminate the information redundancy caused by the correlation between the input data, thereby improving prediction accuracy.

1. 引言

水稻是一种主要的谷类作物，随着全球人口的持续增长，对其产量的需求也随之大幅增加。我国水稻的产量占农业总产量的五分之一 [1]，因此预测其产量对农业发展与建设具有非常重要的意义 [2]。同时全球变暖已经成为了一个不争的客观事实 [3]，其带来的影响不可避免地对农作物生产产生了冲击 [4]，结合各种气象因素来对水稻产量进行预测，不仅有助于了解气象因素对水稻产量产生的影响，还有利于提高产量预测的精度 [5]。

由于受多种因素影响，使用传统方法进行水稻产量的准确预测较为困难。常规的产量预测主要包括基于统计回归原理预测、灰色理论生长模拟 [6]、卫星遥感预测 [7] 等方法，这些方法虽然简便易行，可以用于预测产量，但实际生产中水稻产量往往与多种因素之间存在较为复杂的非线性关系 [8]，使用这些传统的预测方法往往无法实现较高的预测精度。随着人工智能相关技术的不断发展与研究 [9] [10]，神经网络因其对分析因素之间复杂非线性关系的优势，在作物产量预测方面得到了广泛的应用 [11]，而长短期记忆(LSTM)循环神经网络结合了处理时间序列数据的能力，能够考虑到作物产量随时间波动的特点。同时影响作物产量的因素多且复杂，且各因素之间往往存在着较强的相关性，彼此之间既相互促进又相互制约，共同决定了作物的产量 [12]，这为作物产量预测的建模带来了困难 [13]，主成分分析是解决该问题的方法之一 [14]。主成分分析是一种数据降维的方法 [15]，可以将原始变量转换为一组基于原始变量的线性组合，其目标是在低维子空间中表示高维数据，使得在最小误差平方和的意义下，低维能够较好地描述原始高维数据 [16]。

本文采用LSTM对水稻的产量进行预测，同时提出了一种改进的主成分分析方法对神经网络模型的输入进行数据降维，去除原始数据中的噪声和冗余，同时尽可能多地保留原始输入所包含的信息，提高预测精度。

2. 相关工作

2.1. LSTM循环神经网络

LSTM是一种特殊的循环神经网络(RNN)，通过在标准RNN中引入门控单元的概念来解决标准RNN中可能存在的梯度消失问题 [17]。LSTM循环神经网络的基本结构如图1所示。其中 $f_t$ 为遗忘门， $i_t$ 为输入门， $o_t$ 为输出门， ${\tilde{c}}_t$ 是标准RNN单元中的结构。 $x_t$ 表示一个输入节点，用来对应一个特征参数， $h_t$ 表示t时刻该单元的信息输出， $C_t$ 表示t时刻记忆单元的状态。 $\sigma$ 一般选择Sigmoid作为激励函数，主要起门控作用，由于Sigmoid函数的输出为0-1，当输出接近0或1时符合物理意义上的关和开。tanh函数作为生成候选记忆单元c的选项，因其输出在−1到1之间，符合多数场景下中心为0的特征分布，且梯度(求导)在接近0处的收敛速度快于Sigmoid函数。

Figure 1. Basic structure of LSTM recurrent neural network

图1. LSTM循环神经网络基本结构

由图1可见LSTM的隐藏层通过增加神经元的复杂性，在前向计算中增加了遗忘门、输入门、输出门和记忆单元，以此保持和控制信息流动，一定程度上结局了梯度消失以及短时记忆的问题。

观察遗忘门的输出，可以得到公式(1)：

$f_t=\sigma \left(W_{xf}{x}_t+W_{hf}{h}_{t-1}+b_f\right)$ (1)

其中 $W_{xf}$ 表示权重，下标中的f表示这是输入到遗忘门的权重，x表示是 $x_t$ 的权重，故这是 $x_t$ 输入到遗忘门中的权重，其余权重同理， $b_f$ 表示遗忘门偏置值，因此可以得到输入门和输出门以及 ${\tilde{c}}_t$ 的输出，分别见公式(2)、公式(3)和公式(4)：

$i_t=\sigma \left(W_{xi}{x}_t+W_{hi}{h}_{t-1}+b_i\right)$ (2)

$o_t=\sigma \left(W_{xo}{x}_t+W_{ho}{h}_{t-1}+b_o\right)$ (3)

${\tilde{c}}_t=\tanh \left(W_{x\tilde{c}}{x}_t+W_{h\tilde{c}}{h}_{t-1}+b_{\tilde{c}}\right)$ (4)

由公式和图1不难理解，遗忘门反映的是输入x和上一层隐藏层输出h被遗忘的程度的大小；输入门反映的是输入x和当前计算的状态更新到记忆单元的程度大小；输出门反映了输入x和当前输出取决于当前记忆单元的程度大小。将这些输出通过运算可得到 $C_t$ 和 $h_t$，作为下一时刻的输入节点和隐藏层输入信息。计算方法见公式(5)和公式(6)：

$C_t=C_{t-1}\ast f_i+{\tilde{c}}_t\ast i_t$ (5)

$h_t=o_t\ast \tanh \left(C_t\right)$ (6)

对于一个训练好的LSTM模型，由遗忘门、输入门、输出门和内部记忆单元共同控制LSTM模型输出h的设计可以使整个网络更好地把握序列信息之间的关系。

2.2. 主成分分析

主成分分析是一种对数据进行降维的统计方法 [18]，它借助于一个正交变换，将分量相关的原随机变量转化为分量不相关的新随机变量，以方差来衡量所包含原始变量信息的多少。在线性代数上表现为将原随机变量的协方差矩阵变换成对角矩阵，在几何上的意义则是将原坐标系通过移动和旋转变换成新的坐标系，使得所有样本点距离新的坐标轴最近，且新的坐标轴指向样本点散布最开的方向，即方差最大的方向。假设F1为方差最大的一个线性组合，那么称F1为第一主成分，如果F1不足以表征原变量的全部信息，则考虑选取剩下的线性组合中方差最大的线性组合F2作为第二主成分，依此类推可以根据需要构造出第三、第四以及更多主成分。

记原始变量 $y={\left(Y_1,Y_2,Y_3,\cdots ,Y_P\right)}^{\prime }$，其协方差矩阵为 $\Sigma$。主成分分析的目的即是定义一组互不相关的变量，成为 $Y_1,Y_2,Y_3,\cdots ,Y_P$ 的主成分(Principal component, PC)，记为 $Z_1,Z_2,Z_3,\cdots ,Z_P$，每一个主成分都是 $Y_1,Y_2,Y_3,\cdots ,Y_P$ 的线性组合：

$\begin{array}{l}{Z}_1={a^{\prime }}_{1}y=a_{11}{Y}_1+a_{12}{Y}_2+a_{13}{Y}_3+\cdots +a_{1p}{Y}_p\\ Z_2={a^{\prime }}_{2}y=a_{21}{Y}_1+a_{22}{Y}_2+a_{13}{Y}_3+\cdots +a_{2p}{Y}_p\\ ⋮\\ Z_p={a^{\prime }}_{p}y=a_{p1}{Y}_1+a_{p2}{Y}_2+a_{p3}{Y}_3+\cdots +a_{pp}{Y}_p\end{array}$

由于使用线性组合的方差作为指标来衡量主成分所包含的原变量信息的多少，因此问题转化为最大化每一个主成分的方差，即 $\mathrm{var}\left(Z_j\right)={a^{\prime }}_j\Sigma a_j$ 最大，同时主成分之间应满足不相关性，否则达不到充分降维的目的，故需要同时保证 $\mathrm{cov}\left(Z_j,Z_k\right)={a^{\prime }}_j\Sigma a_k=0$。

依次对每一个主成分按照方差贡献度的指标进行排序并分析可知，第一主成分 $Z_1={a^{\prime }}_{1}y$ 的方差是最大的，即 ${a^{\prime }}_1\Sigma a_1$ 最大，观察其方差的表达式不难看出如果 $a_1$ 的长度越大，则会导致 ${a^{\prime }}_1\Sigma a_1$ 也越大，但讨论方差贡献度时并不希望方差的大小是依赖 $a_1$ 的长度的，因此加上限定条件 ${a^{\prime }}_1{a}_1=1$，则 $a_1$ 成为一个单位向量。同理可得第二主成分 $Z_2={a^{\prime }}_{2}y$ 需要在满足 ${a^{\prime }}_2{a}_2=1$ 的同时最大化方差 ${a^{\prime }}_2\Sigma a_2$，但由于主成分之间需要满足不相关性，因此 $\mathrm{cov}\left(Z_1,Z_2\right)={a^{\prime }}_1\Sigma a_2=0$。以此类推，第j主成分 $Z_j={a^{\prime }}_{j}y$ 则需要满足 ${a^{\prime }}_j{a}_j=1$，且 ${a^{\prime }}_j\Sigma a_k=0$，其中 $k=1,2,\cdots ,j-1$ 的同时，最大化方差 ${a^{\prime }}_j\Sigma a_j$。

根据线性代数定理，对于任意向量 $a$ 和 $p\times p$ 的对称矩阵 $\Sigma$，记 $\left({\lambda }_j,e_j\right)$ 为 $\Sigma$ 的特征对，其中 $j=1,2,3,\cdots ,p$ 且 ${\lambda }_1\ge {\lambda }_2\ge {\lambda }_3\ge \cdots \ge {\lambda }_p\ge 0$，则有 $\left(a^{\prime }\Sigma a\right)/\left(a^{\prime }a\right)$ 在 $a=e_1$ 时取得最大值，且最大值为 ${\lambda }_1$。当 $a$ 与 $e_1,e_2,\cdots ,e_{j-1}$ 都正交时， $\left(a^{\prime }\Sigma a\right)/\left(a^{\prime }a\right)$ 在 $a=e_j$ 时取得最大值，且最大值为 ${\lambda }_j$，其中 $j=2,3,\cdots ,p-1$。

由之前的条件可知，以 ${a^{\prime }}_2$ 为例，当满足 ${a^{\prime }}_2\Sigma a_1=0$ 时，由于 $a_1=e_1$，故可推出 ${a^{\prime }}_2\lambda e_1=0$，即 ${a^{\prime }}_2$ 与 $e_1$ 是垂直的，可以看出该定理可以用于解决本文的最大化方差问题。则变量 $Y_1,Y_2,Y_3,\cdots ,Y_P$ 的第j个主成分公式如式(7)：

$Z_j={e^{\prime }}_{j}y=e_{j1}{Y}_1+e_{j2}{Y}_2+e_{j3}{Y}_3+\cdots +e_{jp}{Y}_p,j=1,2,\cdots ,p$ (7)

其中 $\left({\lambda }_j,e_j\right)$ 为协方差矩阵 的特征对， ${\lambda }_1\ge {\lambda }_2\ge {\lambda }_3\ge \cdots \ge {\lambda }_p\ge 0$ 且特征向量均为标准化特征向量，而且有如下性质：

$\mathrm{var}\left(Z_j\right)={e^{\prime }}_j\Sigma e_j={\lambda }_j$ (8)

$\mathrm{cov}\left(Z_j,Z_k\right)={e^{\prime }}_j\Sigma e_k=0$ (9)

进一步地可以得到式(10)：

${\displaystyle {\sum }_{j=1}^p\mathrm{var}\left(Z_j\right)}={\displaystyle {\sum }_{j=1}^p\mathrm{var}\left(Y_j\right)}$ (10)

由(10)式可以看出，得到的主成分可以完全描述原始变量的方差。

实际生产中，变量的方差可能会因为量纲的不同而导致部分变量方差过大或过小的问题，为了消除量纲带来的影响，可以将原始变量除以标准差，从而构造出一个新的标准化矩阵来进行主成分分析，即：

$W_j=\left(Y_j-u_i\right)/\sqrt{{\sigma }_{jj}}$ (11)

其中分子的作用是进行中心化。而新变量的协方差矩阵即为原始变量的相关系数矩阵，因此可以通过求原始变量的相关系数矩阵的特征值和特征向量来消除量纲的影响，达到标准化的目的。但是需要注意的是，这两种方法分析得到的两组主成分的数值和方向是没有关系的，二者之间并不能通过简单的线性变换进行转化。

2.3. 主成分分析存在的问题及改进

特征与特征之间的影响程度可以由特征间的相关系数来刻画，而特征的变异程度则由变异系数(Coefficient of Variation, CV)来刻画，变异系数计算方式如公式(12)所示，其中分子为标准差，分母为平均值：

$\text{CV}=\frac{\sigma }{\mu }$ (12)

如公式(11)所示的标准化过程虽然可以消除量纲对方差带来的影响，但经过标准化后的变量计算是基于原始变量的相关系数矩阵来进行的，观察相关系数矩阵可以发现其对角线元素都是1，即标准化使得变量的方差都变成了1，也就忽略了特征本身的变异程度，从而不能完整刻画原始变量的全部信息。

为了解决上述问题，采用如公式(13)所示的标准化方法：

$y_{ij}=\frac{x_{ij}}{\overline{x_j}}$ (13)

其中 $x_{ij}$ 为原始变量n × p矩阵中的元素， $y_{ij}$ 为标准化后的矩阵中的元素， $\overline{x_j}$ 为原始变量矩阵第j列的均值，即第j个特征的均值，计算方法如公式(14)所示：

$x_j=\frac{1}{n}{\displaystyle {\sum }_{k=1}^n{x}_{kj}}$ (14)

设标准化后的矩阵的协方差矩阵为 $C={\left(c_{ij}\right)}_{n\times p}$，其计算方法如公式(15)所示：

$c_{ij}=\frac{1}{n}{\displaystyle {\sum }_{k=1}^n\left(y_{ki}-\overline{y_i}\right)\left(y_{kj}-\overline{y_j}\right)}$ (15)

由公式(13)和公式(14)易推得标准化后的矩阵中每个特征的均值都为1，即 $\overline{y_i}=1$，可得：

$c_{ij}=\frac{1}{n}{\displaystyle {\sum }_{k=1}^n\left(\frac{x_{ki}}{\overline{x_i}}-1\right)\left(\frac{x_{kj}}{\overline{x_j}}-1\right)}$ (16)

化简式(16)可得：

$c_{ij}=\frac{\frac{1}{n}{\displaystyle {\sum }_{k=1}^n\left(x_{ki}-\overline{x_i}\right)\left(x_{kj}-\overline{x_j}\right)}}{\overline{x_i}\cdot \overline{x_j}}$ (17)

不难看出式(17)中分子为原始变量的协方差矩阵中的元素，特别地，矩阵C的对角线元素为原始变量的变异系数的平方，因此该标准化没有丢失变量变异程度的信息。同时观察标准化后矩阵的相关系数，设原始变量的相关系数为 ${\rho }_{ij}$，经过标准化后为 $r_{ij}$， $s_{ij}$ 为原始变量协方差矩阵中的元素：

$r_{ij}=\frac{c_{ij}}{\sqrt{c_{ii}}\cdot \sqrt{c_{jj}}}$ (18)

由式(17)可知 $c_{ij}=\frac{s_{ij}}{\overline{x_i}\cdot \overline{x_j}}$，代入式(18)并化简，得到公式(19)如下：

$r_{ij}=\frac{s_{ij}}{\sqrt{s_{ii}}\cdot \sqrt{s_{jj}}}$ (19)

而 $\frac{s_{ij}}{\sqrt{s_{ii}}\cdot \sqrt{s_{jj}}}={\rho }_{ij}$，因此该标准化也不会丢失原始变量的相关系数信息，综上所述，该标准化方法可以更完整地刻画原始变量的信息。

3. 模型构建

对原始数据的输入进行主成分分析，将输入数据降维后作为LSTM的输入节点，模型结构如图2所示。其中 $X_1,X_2,\cdots ,X_7$ 为原始输入变量，经过主成分分析后降维为两个变量，即PC1和PC2，再以这两个主成分作为神经网络的输入节点进行预测并输出产量(Yield)。

Figure 2. Structure of LSTM based on PCA

图2. 基于主成分分析的LSTM结构

3.1. 主成分分析相关计算

首先计算数据集的KMO (Kaiser-Meyer-Olkin)检验值和Bartlett球形检验值，以判断该数据集是否适合进行主成分分析。KMO值大于0.5且越接近1时就越适合进行主成分分析。Bartlett球形检验值小于显著水平0.05或者0.01时适合进行主成分分析。

在本文提出的改进主成分分析中，先将数据进行标准化后使用其协方差矩阵的特征值作为指标来选取主成分，选取主成分一般遵循以下原则：特征值大于1；累计方差贡献率大于85%。

3.2. 神经网络模型构建

根据Kolmogorov定理，含有输入层、隐含层和输出层的三层神经网络可以精确表征任意一个连续函数 [19]，因此本文同样采用输入层、隐含层和输出层的三层结构来构建神经网络，隐含层激活函数采用Softmax函数，输出层激活函数采用恒等式。训练的参数通过比较不同取值对预测结果的影响来确定，最终结果如表1所示。

Table 1. Parameters

表1. 参数设置

由最终确定的节点个数可以看出，BP神经网络的最优节点数显著小于LSTM的最优节点数，原因在于过于复杂的网络拓扑结构虽然可以提高拟合的精度，但也可能导致这两种神经网络的过拟合问题，而过拟合问题对于BP神经网络而言较为突出，LSTM因其本身独特的门控结构使得其的权重更新不会受到所有数据的影响，而是可以选择记住或者遗忘前时刻的特征，因此一定程度上降低了产生过拟合的可能性 [20]。

4. 实验结果及分析

4.1. 数据集

为了验证模型的有效性，本文使用的某地水稻产量数据集中原始变量包含降雨量(Precipitation)、最低温度(Minimum Temperature)、平均温度(Mean Temperature)、最高温度(Maximum Temperature)、相关作物需水量(Reference Crop Evapotranspiration)、种植面积(Area)和总产量(Production)。

4.2. 相关参数计算及数据预处理

由表2数据可知，该数据集KMO > 0且Sig < 0.01，因此变量之间相关性较高，存在信息重叠，故可以进行主成分分析。

Table 2. KMO and Bartlett’s test of sphericity

表2. KMO和巴特利特球形检验

将数据集按改进的标准化方法处理后计算其协方差矩阵及其特征值来进行主成分的选取，计算结果如表3所示。

Table 3. Variance contribution rate

表3. 方差贡献率

根据选取主成分的原则可知，该数据集应选取的主成分个数为2，设第一主成分为F1，第二主成分为F2，其成分矩阵如表4所示。

Table 4. Component matrix

表4. 成分矩阵

由表4中主成分系数可得主成分表达式如(20)式和(21)式所示。

$\text{F}1=-0.781X_1+0.796X_2+0.949X_3+0.991X_4+0.925X_5+0.058X_6+0.076X_7$ (20)

$\text{F}2=0.490X_1+0.325X_2+{0.1573}_1+0.023X_4-0.192X_5+0.973X_6+0.963X_7$ (21)

根据(12)式和(13)式计算主成分F1和F2，则将原变量由7维降到了2维。

4.3. 模型对比

通过在相同数据集上对比标准BP神经网络的预测效果和经过改进主成分分析后的LSTM神经网络(IPCA-LSTM)预测效果，以验证本文模型的有效性。评价指标采用均方根误差(RMSE)和平均绝对百分比误差(MAPE)，计算结果见公式(22)和公式(23)，其中 $h\left(x_i\right)$ 为预测值， $y_i$ 为观测值。图3和图4分别为BP神经网络和PCA-LSTM在训练集上随着迭代次数增加的损失(Loss)变化情况，横坐标为迭代次数，纵坐标为损失。实验结果如表5所示。

$\text{RMSE}\left(X,h\right)=\sqrt{\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{\left(h\left(x_i\right)-y_i\right)}^2}}$ (22)

$\text{MAPE}\left(X,h\right)=100\%\times \frac{1}{n}{\displaystyle {\sum }_{i=1}^n\left|h\left(x_i\right)-y_i\right|}$ (23)

Figure 3. Loss of BP neural network on training set

图3. BP神经网络训练模型损失

Figure 4. Loss of IPCA-LSTM on training set

图4. IPCA-LSTM训练模型损失

Table 5. System results of different models on the dataset

表5. 不同模型在数据集上的结果

从实验结果不难看出，本文提出的基于改进主成分分析的LSTM神经网络模型相较于仅使用BP神经网络的模型具有更高的预测精度，误差更小且模型训练收敛更快。

5. 结论

本文针对作物产量预测研究中数据维度高和包含时间序列信息的问题，提出了基于改进主成分分析的LSTM神经网络预测模型，在能够处理复杂输入变量之间的非线性关系的同时，可以有效降低原始输入变量的维度，消除其冗余性，提高模型训练时的收敛速度，提高预测精度，应用前景良好。