山西省粮食产量预测研究

doi:10.12677/SA.2022.116148

期刊菜单

山西省粮食产量预测研究
Study on Grain Yield Forecasting in Shanxi Province

DOI: 10.12677/SA.2022.116148, PDF, HTML, XML, 下载: 156 浏览: 298
作者: 汪路, 段婧霞：南京信息工程大学数学与统计学院，江苏南京
关键词: 灰色模型；残差修正；灰色线性组合模型；粮食产量；Gray Model； Residual Correction； Gray Linear Combination Model； Grain Yield

摘要: 本文基于2000~2017年山西省粮食产量数据分别构建了灰色预测模型、残差修正的灰色预测模型以及灰色线性组合模型，并预留2018~2021年的粮食产量数据进行预测精度检验。研究表明：灰色线性组合模型优于单一的灰色预测模型。灰色预测模型、残差修正的灰色预测模型以及灰色线性组合模型的平均绝对百分比预测分别为12.01%、11.99%、4.49%，拟合精度分别为87.99%、88.01%、95.51%，灰色线性组合模型的预测精度较灰色预测模型其预测精度提高了62.61%。

Abstract: This paper constructs a gray forecasting model, a residual-corrected gray forecasting model and a gray linear combination model based on 2000~2017 grain yield data in Shanxi Province, respectively, and sets aside 2018~2021 grain yield data for testing the forecasting accuracy. The study shows that the gray linear combination model is better than the single gray forecasting model. The average absolute percentage forecasts of the gray forecasting model, the residual-corrected gray forecasting model, and the gray linear combination model were 12.01%, 11.99%, and 4.49%, respectively, and the fitting accuracy was 87.99%, 88.01%, and 95.51%, respectively, and the forecasting accuracy of the gray linear combination model improved by 62.61% compared with that of the gray forecasting model.

文章引用：汪路, 段婧霞. 山西省粮食产量预测研究[J]. 统计学与应用, 2022, 11(6): 1421-1430. https://doi.org/10.12677/SA.2022.116148

1. 引言

古往今来，农业一直是人类生存发展的基础，其中粮食更是与人民生活息息相关。伴随着粮食产业的多元化发展，粮食生产行业的研究意义逐渐凸显，充足的粮食供应对国民经济长期稳定发展具有重大意义。同时，我国土地粮食产量预测一直都是农业问题中的重点，其工作意义重大，对国家规划国民经济发展起决定性作用 [1]，精准预测粮食产量也有利于国民经济等领域的发展。只有国家粮食储备充足，民众才能过上更好、更丰富的生活，追求更高层次的精神领域富足 [2]。山西省作为我国农业大省之一，相关政府部门十分重视粮食生产情况，并一直加大对农业产业的投资以保证粮食产量的稳速增长，故对山西省粮食产量的预测研究具备实际意义。

目前，国内外对相关粮食产量的预测方法主要有气象统计预测、遥感技术、动力学模拟、灰色系统、线性回归及时间序列等 [3]。例如，国外学者Hansen等通过反复测量冠层反射率后构建偏最小二乘回归预测冬小麦和春大麦的籽粒产量和蛋白质含量，预测结果较为精确，但此方法耗时较久 [4]；Bail等利用土壤、作物分析仪器(SPAD)的测量不同冬小麦品种的籽粒质量和产量值并进行预测，但此法潜在限制因素较多，结果可能差异较大 [5]；Mkhabela等下载美国宇航局研制的中分辨率成像光谱仪(MODIS)里的归一化植被指数数据(NDVI)后对加拿大大草原作物产量进行预测，其建立的模型在预测加拿大大草原作物产量方面显示出一定的潜力 [6]。以上介绍的方法都各有优劣，虽在其本身研究课题上具有一定的优越性，但仍存在如精确性、耗时长等局限性，有待进一步的改进。而国内主要将目光集中在粮食产量数据及其相关影响因素上，并基于此进行统计分析建模。郭晓婷通过多元回归和时间序列分析了安徽省粮食产量，并对比预测精度，选取多元回归模型进行预测，预测效果良好 [7]；李铭通过多元线性回归模型分析了15个可能影响山东省粮食生产的因素，并构建MLP神经网络模型，结果表明MLP模型具备稳定性和准确性且预测精度较高 [8]。也有不少研究将灰色预测应用到粮食产量上，如雷蕾等基于四川省粮食产量构建灰色模型且预测精度较高，表明其未来粮食产量呈增长趋势 [9]；孟凡琳基于河南省粮食产量构建灰色模型并进行改进，预测效果良好 [10]，但有相关研究显示使用灰色线性组合模型对河南省粮食产量进行预测的预测精度更优 [11]；薛晋芳对山西省的粮食产量进行了分析和预测，利用灰色线性组合模型进行预测，结果符合预期 [12]；陈婷怡等将二次指数平滑模型、灰色模型和支持向量回归模型进行组合优化，利用线性组合预测模型对云南省2021~2023年的粮食产量进行预测 [13]。戴剑勇等将灰色预测模型与逻辑斯蒂预测模型组合，并对相关数据进行预测，虽然不是对粮食产量进行预测，但其组合模型的预测思想值得借鉴 [14]。综上所述，组合模型的预测在总体上来说较单一模型的精度更优，也考虑了更多的可能性，才使得预测的结果更加接近现实，更有说服力，同样推测组合模型对山西省粮食产量的预测也具备较好的能力。

本文基于2000~2021年的山西粮食产量构建了灰色预测模型、灰色修正预测模型以及灰色线性组合模型，从中选取精度最优的预测模型对山西省粮食产量进行预测，并结合预测结果以及实际情况提出合理建议，以期为调整和发展我国粮食产业提供理论参考，为政府部门以及相关生产者们提供依据与决策。

2. 资料与方法

本文基于灰色系统理论对山西省粮食产量进行建模研究，其中粮食产量数据来自国家统计局 (http://www.stats.gov.cn/)。由于粮食产量受环境等多种未知因素的影响，符合灰色理论中的灰白不确定性特征，结合2001~2021年山西省的粮食产量呈上升趋势(图1)，经初步分析可建立灰色预测模型。

Figure 1. Histogram of year-on-year change of grain production in Shanxi Province from 2000 to 2021

图1. 2000~2021年山西省粮食产量逐年变化柱状图

1982年，邓聚龙教授为解决数据缺乏及不确定性问题提出了灰色系统理论，其相关模型又被称为灰色模型或灰色动态模型(简称GM模型)。构建灰色模型即可进行灰色预测，其核心是基于不确定背景下，通过数据处理分析建立模型，对发展趋势进行有效合理的预测评估。实现灰色预测首先需要确定处于某一范围内并与时间有关的灰色过程中的数据，再根据关联分析鉴别系统因素之间发展趋势的相异程度，其次对原始数据进行生成处理去寻找该系统内部的变动的规律，以最终生成的具有较强的规律性的数据序列建立相应的微分方程模型，最后，通过该灰色模型对事物未来的发展趋势进行预测。这里主要介绍灰色建模中数据检验与处理、建立模型以及模型检验三个部分的理论知识，并基于此进行山西省粮食产量灰色预测模型的构建。

2.1. 数据的检验与处理

首先，为保证建模方法的可行性，需要对已知数据列做必要的检验处理。设参考数据为时间序列 $x^{(0)} = (x^{(0)} (1), x^{(0)} (2), x^{0} (3), \dots, x^{(0)} (n))$ ，计算序列的级比：

$λ (k) = \frac{x^{(0)} (k)}{x^{(0)} (k - 1)}, k = 2, 3, \dots, n$ (1)

若(n − 1)个级比 $λ (k)$ 都落在 $Θ = (e^{- \frac{2}{n + 1}}, e^{\frac{2}{n + 1}})$ 范围内，则序列 $x^{(0)}$ 可以作为灰色模型的数据进行灰色预测。否则，需要对序列 $x^{(0)}$ 进行变换，使其落入 $Θ$ 范围内。变换方法为取适当的常数l，令

$y^{(0)} (k) = x^{(0)} (k) + l, k = 1, 2, \dots, n$ (2)

使序列 $y^{(0)} = (y^{(0)} (1), y^{(0)} (2), y^{0} (3), \dots, y^{(0)} (n))$ 的级比：

$λ_{(y)} (k) = \frac{y^{(0)} (k)}{y^{(0)} (k - 1)} \in Θ, k = 2, 3, \dots, n$ (3)

2.2. 灰色模型

已知的时间序列 $x^{(0)} = (x^{(0)} (1), x^{(0)} (2), x^{0} (3), \dots, x^{(0)} (n))$ ，对其一次累加生成新的序列:

$\begin{matrix} x^{(1)} = (x^{(1)} (1), x^{(1)} (2), x^{(1)} (3), \dots, x^{(1)} (n)) \\ = (x^{(0)} (1), x^{(0)} (1) + x^{(0)} (2), \dots, x^{(0)} (1) + \dots + x^{(0)} (n)) \end{matrix}$ (4)

$x^{(1)}$ 的均值生成序列：

$m^{(1)} = (m^{(1)} (2), m^{(1)} (3), \dots, m^{(1)} (n))$ ， (5)

式中，

$m^{(1)} (k) = 0.5 m^{(1)} (k) + 0.5 m^{(1)} (k - 1), k = 2, 3, \dots, n$ 。

建立灰分方程：

$x^{(0)} (k) + a m^{(1)} (k) = b, k = 2, 3, \dots, n$ (6)

相应得白化微分方程为

$\frac{d x^{(1)}}{d t} + a x^{(1)} (t) = b$ (7)

式中，a为发展灰数，b为内生控制灰数。

记a为待估参数向量，

$u = {[a, b]}^{T}$ ， $Y = {[x^{(0)} (2), x^{(0)} (3), \dots, x^{(0)} (n)]}^{T}$ ，

$B = [\begin{matrix} - m^{(1)} (2) & 1 \\ - m^{(1)} (3) & 1 \\ ⋮ & ⋮ \\ - m^{(1)} (n) & 1 \end{matrix}]$ ， (8)

利用最小二乘法，可得u的估计值为：

$\hat{u} = {[\hat{a}, \hat{b}]}^{T} = {(B^{T} B)}^{- 1} B^{T} Y$ (9)

求解白化微分方程，得

${\hat{x}}^{(1)} (k + 1) = (x^{(0)} (1) - \frac{\hat{b}}{\hat{a}}) e^{- \hat{a} k} + \frac{\hat{b}}{\hat{a}}, k = 0, 1, \dots, n - 1, \dots$ (10)

2.3. 灰色残差修正模型

灰色残差修正模型是对灰色预测模型进行残差检验后，根据实际值与预测值的误差，为提高灰色预测模型的预测精度而构建的模型，其主要目的是为了丰富和完善灰色预测模型。过程如下：

首先，根据原始序列得其残差序列为：

$e^{(0)} = (e^{(0)} (k), e^{(0)} (k + 1), \dots, e^{(0)} ( n ))$

为便于计算，改写为：

$e^{(0)} = (e^{(0)} (1^{'}), e^{(0)} (2^{'}), \dots, e^{(0)} (n^{'}))$ (11)

$e^{(0)}$ 的累加生成序列为：

$e^{(1)} = (e^{(1)} (1^{'}), e^{(1)} (2^{'}), \dots, e^{(1)} (n^{'})), (n^{'} = n - k)$

$e^{(1)}$ 可建立相应的GM (1, 1)模型：

${\hat{e}}^{(1)} (k + 1) = [e^{(0)} (1) - \frac{μ_{e}}{a_{e}}] e^{- a_{e} k} + \frac{μ_{e}}{a_{e}}$

${\hat{e}}^{(1)} (k + 1)$ 的导数为：

${\hat{e}}^{'} (k + 1) = (- a_{e}) [e^{(0)} (1) - \frac{μ_{e}}{a_{e}}] e^{- a_{e} (k - 1)}$

加上修正项，得到修正模型：

${\hat{x}}^{(1)} (k + 1) = [x^{(0)} (1) - \frac{μ}{a}] e^{- a k} \frac{μ}{a} + δ (k - 1) (- a_{e}) [e^{(0)} (1) - \frac{μ_{e}}{a_{e}}] e^{- a_{e} (k - 1)}$ (12)

其中， $δ (k + 1) = {\begin{matrix} 1, k \geq 2 \\ 0, k < 2 \end{matrix}$ 为修正系数。

2.4. 灰色线性组合模型

灰色线性组合模型适用同时具备指数增长趋势和线性趋势序列，利用线性回归方程和指数方程的有用信息来拟合之前的累加生成序列，可以将生成的序列写成：

${\hat{x}}^{(1)} (k) = C_{1} e^{v k} + C_{2} k + C_{3}$ (13)

v和 $C_{1}$ 、 $C_{2}$ 、 $C_{3}$ 是待定参数。

确定组合模型参数。设有参数序列：

$Z (k) = {\hat{x}}^{(1)} (k + 1) - {\hat{x}}^{(1)} (k) = C_{1} e^{v k} (e^{v} - 1) + C_{2}, k = 1, 2, \dots, n$ (14)

又假设

$Y_{m} (k) = Z (k + m) - Z (k) = C_{1} e^{v k} (e^{v m} - 1) (e^{v} - 1)$

$Y_{m} (k + 1) = C_{1} e^{v (k + 1)} (e^{v m} - 1) (e^{v} - 1)$ (15)

可得： $v = \ln [Y_{m} (k + 1) / Y_{m} (k)]$

将累加生成的序列代入公式(14)中可以得到v的近似值 $\hat{v}$ ，不同的m得到不同的 $\tilde{v}$ ，然后对其取平均得到估值 $\hat{v}$ 。即

$\hat{v} = \frac{\sum_{m = 1}^{n - 3} \sum_{k = 1}^{n - 2 - m} {\tilde{v}}_{m} (k)}{(n - 2) (n - 3) / 2}$ (16)

代入公式(2.1)中，利用最小二乘法可以解出 $C_{1}$ 、 $C_{2}$ 、 $C_{3}$ 的估计值。

GM (1, 1)模型对序列处理结果记为：

$x^{(1)} = [\begin{matrix} x^{(1)} (1) \\ x^{(1)} (2) \\ ⋮ \\ x^{(1)} (n) \end{matrix}]$ ， $C = [\begin{matrix} C_{1} \\ C_{2} \\ C_{3} \end{matrix}]$ ， $A = [\begin{matrix} e^{\hat{v}} & 1 & 1 \\ e^{2 \hat{v}} & 2 & 1 \\ ⋮ & ⋮ & ⋮ \\ e^{n \hat{v}} & n & 1 \end{matrix}]$ (17)

可得： $C = {(A^{T} A)}^{- 1} A^{T} x^{( 1 )}$

然后计算出累加生成序列的预测值，再通过累减生成进行还原，就可以得到原序列的预测序列。

2.5. 模型检验

对灰色预测模型得检验主要包括残差检验、关联度检验以及后验差检验。

1) 残差检验

令残差为 $ε (k)$ ，计算

$ε (k) = \frac{x^{(0)} (k) - {\hat{x}}^{(0)} (k)}{x^{(0)} (k)}, k = 1, 2, \dots, n$ (18)

这里 ${\hat{x}}^{(0)} (1) = x^{(0)} (1)$ ，若 $ε (k) < 0.2$ ，则认为达到一般要求；若 $ε (k) < 0.1$ ，则认为达到较高的要求。

2) 关联度检验

计算出 ${\hat{x}}^{(0)} (k)$ 与原始序列 $x^{(0)} (k)$ 的关联系数，然后计算出关联度，根据经验，当 $ρ = 0 .5$ 时，关联度大于0.6即可。

3) 后验差检验。计算原始序列的标准差：

$S_{1} = \sqrt{\frac{\sum {[x^{(0)} (k) - {\hat{x}}^{(0)}]}^{2}}{n - 1}}$ (19)

计算绝对误差序列的标准差：

$S_{2} = \sqrt{\frac{\sum {[Δ^{(0)} (k) - {\bar{Δ}}^{(0)}]}^{2}}{n - 1}}$ (20)

计算方差比：

$C = \frac{S_{2}}{S_{1}}$ (21)

计算小误差概率：

$P = p {| Δ^{(0)} (k) - {\bar{Δ}}^{(0)} | < 0.6745 S_{1}}$ (22)

令 $e_{k} = | Δ^{(0)} (k) - {\bar{Δ}}^{(0)} |$ ， $S_{0} = 0.6745 S_{1}$ ，则 $P = p {e_{k} < S_{0}}$ ，关于P与C的取值对应的拟合效果如下：

若 $P > 0.95, C < 0.35$ ，则拟合效果好；

若 $P > 0.80, C < 0.50$ ，则拟合效果合格；

若 $P \leq 0.70, C \geq 0.65$ ，则效果不合格。

3. 模型构建

这里仅以2000~2017年山西省粮食产量数据进行三类模型构建，并预留2018~2021年的粮食产量数据进行预测精度检验，选取最优模型。

利用R软件进行灰色预测模型的构建，可得到其待估参数向量为：

$\hat{u} = [- 0.0328, 810.8596]$

故其灰色预测模型为：

$x^{(1)} (k + 1) = 25574.68 \times e^{0.0328 k} - 24721 .33$

在构建灰色残差修正模型前，需要根据实际值和由灰色预测模型得出的预测值计算残差，并利用残差序列构建残差序列的GM模型，得到其待估参数向量为：

${\hat{u}}_{e} = [0.0797, 0.1122]$

基于残差序列得到的灰色预测模型为：

$e^{(1)} (k + 1) = - 1.4078 \times e^{- 0.0797 k} + 1 .4078$

其导数为

${\hat{e}}^{'}^{(1)} (k + 1) = - 0.0797 (- 1.4078 \times e^{- 0.0797 k - 1})$

故可得到灰色修正模型为：

${\hat{x}}^{(1)} (k + 1) = 25574.68 \times e^{0.0 328 k} - 24721 .33 + δ (k - 1) (0.1122 e^{- 0.0797 k - 1})$

其中， $δ (k + 1) = {\begin{matrix} 1, k \geq 2 \\ 0, k < 2 \end{matrix}$

为提高模型的预测精度，尝试构建灰色线性组合模型，首先将生成的序列写成： ${\hat{x}}^{(1)} (k) = C_{1} e^{v k} + C_{2} k + C_{3}$ 形式并确定组合模型参数。

设有参数序列：

$Z (k) = {\hat{x}}^{(1)} (k + 1) - {\hat{x}}^{(1)} (k) = C_{1} e^{v k} (e^{v} - 1) + C_{2}, k = 1, 2, \dots, n$

再假设

$Y_{m} (k) = Z (k + m) - Z (k) = C_{1} e^{v k} (e^{v m} - 1) (e^{v} - 1)$

$Y_{m} (k + 1) = C_{1} e^{v (k + 1)} (e^{v m} - 1) (e^{v} - 1)$

可得： $v = \ln [Y_{m} (k + 1) / Y_{m} (k)] = - 0.1391$

最后进行灰色线性组合模型预测，得到

$C = [\begin{matrix} C_{1} = 6102.6 \\ C_{2} = 1387.1 \\ C_{3} = - 5714.2 \end{matrix}]$

灰色线性组合模型的预测值为：

${\hat{x}}^{(1)} (k) = 6102.6 e^{- 0.1391 k} + 1387.1 k - 5714.2$

4. 模型检验与分析

4.1. 模型检验

模型检验主要包括残差检验、关联度检验以及后验差检验三部分，分别计算三个模型的平均误差、关联度、后验差(表1)。经对比分析可知，三类模型的平均误差均小于0.1，通过残差检验；关联度均大于0.6，通过关联度检验；同时，这三类模型的后验差比值检验C值均小于0.35，小误差概率大于0.95，故三类模型均通过后验差检验，且模型的拟合效果好。

Table 1. Three types of model test results

表1. 三类模型检验结果

4.2. 模型预测

由第3节和第4节的第1小节可知三类灰色预测模型均可适用，接下来分别计算相对误差(表2)，根据2018~2021年的实际值与绝对值的平均绝对百分比误差(MAPE)确定最优模型。

Table 2. Grey forecast of grain production in Shanxi Province from 2000 to 2017

表2. 2000~2017年山西省粮食产量灰色预测

如表3所示，GM (1, 1)模型的MAPE为12.01%，相对于残差修正后的GM (1, 1)模型(11.99%)和灰色线性组合模型(4.49%)的预测精度更低。残差修正后的灰色预测模型相较于与灰色预测模型而言只提高了较小精度，考虑到粮食产量的增长是呈现不完全指数增长，还有存在一定的线性关系，于是使用灰色线性组合模型，其MAPE相对于灰色预测模型提高了62.61%，这说明灰色线性组合模型的预测效果更优。同时，在预测模型相对误差表中可以观察到在2001、2004、2009年的粮食产量的相对误差较高，但是总体预测精度良好。

Table 3. 2018~2021 gray prediction accuracy comparison and MAPE

表3. 2018~2021灰色预测精度对比及MAPE

综上所述，结合三种模型检验以及MAPE，从灰色模型、灰色预测模型以及灰色线性组合模型中选取灰色线性组合模型对山西省的粮食产量进行预测，得到2018、2019、2020、2021年的粮食产量分别为1322.4、1330.8、1338.1、1344万吨，实际值与预测值散点图如图2所示。

Figure 2. Comparison chart of the actual and predicted changes in grain production in Shanxi Province from 2000 to 2021

图2. 2000~2021年山西省粮食产量实际值与预测值变化对比图

5. 结论与展望

本文基于2000~2021年的山西省粮食产量数据分别构建了灰色预测模型、残差修正的灰色预测模型以及灰色线性组合模型，三类模型的平均绝对百分比预测分别为12.01%、11.99%以及4.49%，拟合精度分别为87.99%、88.01%、95.51%，灰色线性组合模型的预测精度最优，较灰色预测模型其预测精度提高了62.61%。虽然构建的灰色预测模型未考虑到山西省粮食产量的环境影响因素以及品种因素，但仍可应用灰色系统理论的特性构建灰色预测模型或进行改进对数据进行预测。结合本文分析，虽然未分析影响山西省粮食产量因素，但在实际情况下，政府部门仍需加强对农业产业的调控，尽可能地控制认为因素，比如农作物的化肥使用量、农药使用量等，在后续工作中需要对影响因子进行分析，解决部分数据误差大的问题。

参考文献

[1]	姜长云. 影响我国粮食安全的新趋势新问题[J]. 人民论坛·学术前沿, 2022(4): 94-100.
[2]	李锁强. 粮食产量创历史新高农业生产结构进一步优化[N]. 中国信息报, 2020-01-20(002).
[3]	许德刚, 邢奎杰, 李凡, 郭奕欣. 粮食产量影响因素分析及预测方法研究[J]. 粮食与油脂, 2022, 35(10): 46-50.
[4]	Hansen, P.M., Jørgensen, J.R. and Thomsen, A. (2002) Predicting Grain Yield and Protein Content in Winter Wheat and Spring Barley Using Repeated Canopy Reflectance Measurements and Partial Least Squares Regression. The Journal of Agricultural Science, 139, 307-318. https://doi.org/10.1017/S0021859602002320
[5]	Le Bail, M., Jeuffroy, M.H., Bouchard, C. and Barbottin, A. (2005) Is It Possible to Forecast the Grain Quality and Yield of Different Varieties of Winter Wheat from Minolta SPAD Meter Measurements? European Journal of Agronomy, 23, 379-391. https://doi.org/10.1016/j.eja.2005.02.003
[6]	Mkhabela, M.S., Bullock, P., Raj, S., Wang, S. and Yang, Y. (2011) Crop Yield Forecasting on the Canadian Prairies Using MODIS NDVI Data. Agricultural and Forest Meteorology, 151, 385-393. https://doi.org/10.1016/j.agrformet.2010.11.012
[7]	郭晓婷. 基于多元回归分析的安徽省粮食产量预测方法研究[D]: [硕士学位论文]. 合肥: 安徽农业大学, 2017.
[8]	李铭, 方华. 山东省粮食产量影响因素探究及产量预测[J]. 农场经济管理, 2022(5): 34-39.
[9]	雷蕾, 李婷, 李鹄. 基于灰色GM(1,1)的四川粮食产量预测研究[J]. 粮食问题研究, 2022(2): 21-24.
[10]	孟凡琳. 基于改进灰色模型的河南省粮食安全量化分析[D]: [硕士学位论文]. 郑州: 河南农业大学, 2017.
[11]	陈祺琪, 李君, 梁保松. 基于灰色多元线性回归的河南省粮食产量预测[J]. 河南农业大学学报, 2012, 46(4): 449-452.
[12]	薛晋芳. 山西省粮食产量的分析及预测研究[D]: [硕士学位论文]. 太原: 山西大学, 2016.
[13]	陈婷怡, 王皓晔. 基于组合预测法对云南省粮食产量的预测[J]. 科学技术创新, 2022(17): 36-39.
[14]	戴剑勇, 郭睿. 基于组合预测模型的煤炭能源消费总量预测[J]. 南华大学学报, 2019, 20(3): 21-23.

为你推荐

友情链接