一种基于R语言的小麦生育期特征参数估算方法

摘要

本发明提出一种基于R语言的小麦生育期特征参数估算方法，利用R语言对小麦生长模拟模型进行重新编写，通过拉丁超立方设置参数初始值，运算得到小麦品种特征参数的先验概率分布，根据初始值和先验概率分布选取候选参数，计算生育期的概率密度函数，判断是否接受新的参数，最终得到了品种各特征参数的后验概率密度分布。本发明的方法采用非线性最小二乘法提高了抽样的效率，避免了参数陷入局部最优；基于贝叶斯理论框架，有效地考虑了参数的先验分布，通过转化原始数据解决了站点之间误差相关性的问题，使得参数调试结果准确高效，在同类模型小麦品种特征参数估算当中具有普遍适用性。

说明书

技术领域

本发明属于农业信息技术领域，尤其是一种基于R语言的小麦生育期特征参数估算方法。

背景技术

自工业革命以来，全球气候变化已经成为一个不争的事实。全球范围的气候变化不仅影响人类赖以生存的环境，也给农业生产带来了严峻的挑战。随着对全球气候变化研究的不断深入，有关气候变化对作物生长发育及产量形成影响的研究也开始受到广泛的关注。在未来气候变化情景下，及时准确的预测作物生长情况，对指导及调控农业生产、提高管理水平、制定适宜的粮食生产政策等意义重大。

作物生产系统是受气候、土壤、品种以及栽培管理措施等因素影响的系统，从而使农业生产者难以综合考虑多因子互作、预测农业生产趋势并量化各因子对农作物生产的影响。作物生长模拟模型为定量描述作物生长发育过程及其与环境和技术的动态关系提供了新的方法和手段。但作物模拟模型在应用时存在着一定的不确定性。这些不确定性主要来源于三个方面，即输入变量导致的不确定性，参数值导致的不确定性以及模型公式导致的不确定性。

不确定性分析是基于模型的风险分析和决策的重要组成部分，它能够给风险分析人员和决策支持者提供模型输出结果的准确性程度。参数与作物模型密切相关，模型中的参数值一般难以直接获得，并且由于获得参数值的方法不同以及输入数据的不同，从而导致模拟结果存在较大的不确定性，增加了使用者基于模型决策的不确定性。然而，已有对作物模拟模型不确定性的研究大都将参数值、输入变量及模型公式这三种不确定性来源混合在一起进行分析，没有研究参数不确定性对模拟结果不确定性产生的影响。

发明内容

本发明所解决的技术问题在于提供一种基于R语言的小麦生育期特征参数估算方法，利用R语言对小麦生长模拟模型进行重新编写，结合非线性最小二乘法，判断开花期和成熟期的实测值与模拟值差方的和的最小值，计算得到差方的和的最小值对应的参数调试结果，再利用R语言计算得到小麦品种的参数范围。

实现本发明目的的技术解决方案为：

一种基于R语言的小麦生育期特征参数估算方法，包括以下步骤：

步骤1：采集N组小麦生长模拟模型WheatGrow的生育期子模块所需的生育期数据和气象数据；

步骤2：利用拉丁超立方体为小麦生育期相关的小麦品种特征参数调试设置10组不同的初始值；

步骤3：采用R语言对小麦生长模拟模型的生育期子模块进行重新编写，并将输出的儒历日期连续化；

步骤4：将田间观测的某一品种Nsites个种植站点的气象数据通过R语言读取到小麦生长模拟模型中，调用第i组小麦品种特征参数初始值，运行小麦生长模拟模型N次，得到模拟的N次小麦主要生育期结果，其中，i＝1,2，…，10，Nsites是数据的站点数；

步骤5：将田间实测的N次小麦主要生育期数据和模拟的N次小麦主要生育期结果相结合，计算生育期中开花期和成熟期实测值与模拟值差方的和

步骤6：若达到最小值，计算得到参数调试结果并保存，利用R语言中vcov()函数得到参数间的协方差，利用cov2cor(vcov())函数得到参数间的相关系数，若i＜10则令i＝i+1并转到步骤4，若i＝10则转到步骤7；若未达到最小值，则在第i组参数的基础上产生新的第i组参数，返回步骤4；

步骤7：计算i站点j年的模型误差和其中为播种-开花的天数误差，为开花-成熟的天数误差，再利用R语言中的Kendall秩相关函数计算和的相关性大小P；

步骤8：若P＜0.5，则和是相关独立的，步骤6中得到的参数调试结果为最终结果；若P≥0.5，则对各站点的原始数据进行转化，使转化后的模型误差相关性P＜0.5，再转到步骤5进行参数调试并得到最终调试结果；

步骤9：选取10套参数调试结果中最小的参数组作为最优值，结合参数间协方差，利用R语言中的正太分布函数mvrnorm(1000,mu,sigma)得到1000套参数值并计算其概率分布，取80％区间为置信区间，该区间内的参数值为该品种的参数范围，其中，mu为最优值，sigma为参数间协方差。

进一步的，本发明的基于R语言的小麦生育期特征参数估算方法，步骤2中的小麦生育期相关的小麦品种特征参数包括：生理春化天数、光周期敏感性、基本早熟性、温度敏感性和灌浆期因子。。

进一步的，本发明的基于R语言的小麦生育期特征参数估算方法，步骤3中采用在开花期和成熟期用一个带小数位的天数替代整数天数来实现小麦生长模拟模型输出的儒历日期连续化。

进一步的，本发明的基于R语言的小麦生育期特征参数估算方法，步骤5中开花期和成熟期实测值与模拟值差方的和的计算公式为：

其中，n_i是i站点的数据年份数，和分别是i站点j年的播种-开花天数和开花-成熟天数的实测值，和分别是i站点j年的播种-开花天数和开花-成熟天数的模拟值。

进一步的，本发明的基于R语言的小麦生育期特征参数估算方法，步骤6中采用非线性最小二乘法判断是否达到最小值。

进一步的，本发明的基于R语言的小麦生育期特征参数估算方法，步骤7中的和的计算公式为：

其中，μ^ant和μ^mat分别是导致播种-开花和开花-成熟模拟误差的固定值，和分别是由站点效应导致播种-开花和开花-成熟模拟误差的随机值，和分别是播种-开花和开花-成熟的残差。

进一步的，本发明的基于R语言的小麦生育期特征参数估算方法，步骤8中的对各站点的原始数据进行转化的步骤包括：

步骤8-1：将和的值利用R语言中的lmer函数计算出误差的方差-协方差矩阵∑；

步骤8-2：采用R语言中的Cholesky对∑进行分解：∑＝R^TR，得到R^T，利用R^T对原始数据Y进行转化：Y^*＝(R^T)^-1Y，得到转化后的数据Y^*。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

1、本发明的方法通过基于R语言的非线性最小二乘法和小麦生长模型相结合，有效反演出小麦品种特征参数；

2、本发明的方法在计算生育期时，不仅将各个生育期看作独立分布的，而且考虑到各生育期之间的相互作用，从而确保参数调试结果的准确性；

3、本发明的方法采用的非线性最小二乘法提高了抽样的效率，避免了参数陷入局部最优；

4、本发明的方法基于贝叶斯理论框架，以较大值为搜索起始点，不断从两端向中间寻找误差最小的值，有效地考虑了参数的先验分布，通过转化原始数据解决了站点之间误差相关性的问题，使得参数调试结果准确高效。

附图说明

图1是本发明的方法流程图；

图2是小麦生长模拟模型的框架结构图；

图3是京冬8号小麦的生育期特征参数概率分布图；

图4是京冬8号小麦的生育期模拟值与实测值的1：1关系图。

具体实施方式

下面详细描述本发明的实施方式，所述实施方式的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施方式是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。

利用本发明的基于R语言的小麦生育期特征参数估算方法对1980-2010年在12个站点种植63次的冬小麦品种京冬8号的生育期数据进行有效地估计，并结合图1的流程进行详细说明，如图1所示具体包括以下步骤：

步骤1：采集63组小麦生长模拟模型WheatGrow(如图2所示)所需的生育期数据和气象数据，包括12个站点。其中，采集气象数据：搜集整理1980-2010年研究站点的逐日气象数据，涉及日最高气温、日最低气温、日平均气温、日照时数和日降水量5个气象要素，来自于中国气象科学数据共享服务网(http://cdc.cma.gov.cn/)；采集小麦生育期数据：包括1980-2010京冬8号生育期资料(包括播种期、开花期、成熟期)。

数据均来自国家气象中心资料室。冬小麦京冬8号基本种植信息见表1所示：表1京冬8号实测生育期

步骤2：利用拉丁超立方体为小麦生育期相关的小麦品种特征参数调试设置10组不同的初始值，小麦品种特征参数具体包括：生理春化天数(PVT)、光周期敏感性(PS)、基本早熟性(IE)、温度敏感性(TS)和灌浆期因子(FDF)，具体公式为：

ParaTable[i,1]＜-runif(1,0.6+C1[i]/(10/0.9)-0.12,0.6+C1[i]/(10/0.9))

ParaTable[i,2]＜-runif(1,C2[i]/(10/3)-0.3,C2[i]/(10/3))

ParaTable[i,3]＜-runif(1,C3[i]/(10/0.01)-0.001,C3[i]/(10/0.01))

ParaTable[i,4]＜-runif(1,C4[i]/(10/60)-6,C4[i]/(10/60))

ParaTable[i,5]＜-runif(1,0.8+C5[i]/(10/0.7)-0.12,0.8+C5[i]/(10/0.7))

ParaTable[i,1]为IE，ParaTable[i,2]为TS，ParaTable[i,3]为PS，ParaTable[i,4]为PVT，ParaTable[i,5]为FDF，i＝1,2,…,10。

得到的初始值如表2所示：

表2 京冬8号10组初始参数值

步骤3：采用R语言对小麦生长模拟模型的生育期子模块进行重新编写，并采用在开花期和成熟期用一个带小数位的天数替代整数天数将小麦生长模拟模型输出的儒历日期连续化。在本模型中，存在着几个临界温度变化的转折点，临界温度变化导致模型公式变化，从而导致模拟值输出的不连续性，据此进行连续化的点见表3所示：

表3 公式改变时的PDT值

步骤4：将田间观测的某一品种12个种植站点的气象数据通过R语言读取到小麦生长模拟模型中，调用第i组小麦品种特征参数初始值，运行小麦生长模拟模型63次，得到模拟的63次小麦主要生育期结果，其中，i＝1,2，…，10。

步骤5：将田间实测的63次小麦主要生育期数据和模拟的63次小麦主要生育期结果相结合，如图4所示，计算生育期中开花期和成熟期实测值与模拟值差方的和的计算公式为：

其中，n_i是i站点的数据年份数，和分别是i站点j年的播种-开花天数和开花-成熟天数的实测值，和分别是i站点j年的播种-开花天数和开花-成熟天数的模拟值。

步骤6：采用非线性最小二乘法判断是否达到最小值，若达到最小值，计算得到参数调试结果并保存，如表4所示：

表4 京冬8号10组参数估计结果

然后利用R语言中vcov()函数得到参数间的协方差，如表5所示：

表5 京冬8号参数间协方差

利用cov2cor(vcov())函数得到参数间的相关系数，如表6所示：

表6 京冬8号参数间相关系数

若i＜10则令i＝i+1并转到步骤4，若i＝10则转到步骤7；若未达到最小值，则在第i组参数的基础上产生新的第i组参数，返回步骤4；

步骤7：计算i站点j年的模型误差和其中，为播种-开花的天数误差，为开花-成熟的天数误差，再利用R语言中的Kendall秩相关函数计算和的相关性大小P；和的计算公式为：

其中，μ^ant和μ^mat分别是导致播种-开花和开花-成熟模拟误差的固定值，和分别是由站点效应导致播种-开花和开花-成熟模拟误差的随机值，和分别是播种-开花和开花-成熟的残差。

R语言中Kendall的秩相关函数表达为：

和的计算结果如表7所示：

表7 京冬8号实测值与模拟值的差值

步骤8：若P＜0.5，则和是相关独立的，步骤6中得到的参数调试结果为最终结果，即表4所示；若P≥0.5，则对各站点的原始数据进行转化，使转化后的模型误差相关性P＜0.5，再转到步骤5进行参数调试并得到最终调试结果，具体的转化步骤包括：

步骤8-1：将和的值利用R语言中的lmer函数计算出误差的方差-协方差矩阵∑；

步骤8-2：采用R语言中的Cholesky对∑进行分解：∑＝R^TR，得到R^T，利用R^T对原始数据Y进行转化：Y^*＝(R^T)^-1Y，得到转化后的数据Y^*。

京冬8号品种的相关性大小P＝0.26<0.5，因此不需要对该品种原始数据进行转化，以上通过非线性最小二乘法调试出的结果即为最终的结果，即表4所示。从中选取为最小时的参数值作为最优值，如表8所示：

表8 京冬8号最优值

步骤9：选取10套参数调试结果中最小的参数组作为最优值，结合参数间协方差，利用R语言中的正太分布函数mvrnorm(1000,mu,sigma)得到1000套参数值并计算其概率分布，取80％区间为置信区间，该区间内的参数值为该品种的参数范围，如图3所示，其中，mu为最优值，即表8所示，sigma为参数间协方差，即表5所示。

以上所述仅是本发明的部分实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进，这些改进应视为本发明的保护范围。