基于作物模型、历史和气象预报数据的作物产量预报方法与流程

本发明属于农业遥感领域，具体涉及一种基于作物模型、历史和气象预报数据的作物产量预报方法。

背景技术：

随着作物生长模型应用与研究的不断深入，目前作物模型已从最早的单点尺度模拟扩展到了区域尺度的模拟，而且在气候变化情景模拟和遥感数据同化方面应用尤为广泛。而农业气象领域对于作物估产，尤其是大区域作物产量的准确预报一直都有应用需求，因此实现基于作物模型同化的产量预报是一个亟待解决的问题。

在作物模型应用方面，由于作物模型对输入气象要素的时间分辨率要求较高，而且需要输入全生育期的完整气象要素才能模拟得到较好的产量结果。基于作物模型的产量预报方法主要还是受制于预报期内的气象资料的获取。在现有的研究中，部分研究直接利用区域气候模型生成的数据作为月度或季度的气象预报数据，由于生成的气象数据存在一定的误差，还需要对气候模型的输出进行纠偏和订正，同时低分辨率的气候模型模拟的逐日天气数据也会带来时间和空间尺度不匹配的问题。

也有部分研究将气象资料多年平均值作为预报期的驱动数据。由于我国大部地区为典型的季风气候，年际变率大，且季节分布不均，因此，这种方法并不能很好地体现气象要素的波动性，而且灾害性的异常天气会被平滑，这对模拟潜在的产量损失十分不利；而如果将历史同期所有年份的气象要素都作为集合成员进行驱动模拟，则可以重现本地历史天气情景下的产量，通过集合预报的技术手段得出产量预报的不确定性，将以往单一数值的预报提升为概率预报并提升其预报能力。

技术实现要素：

本发明提供一种基于作物模型、历史和气象预报数据的作物产量预报方法。

本发明提供一种基于作物模型、历史和气象预报数据的作物产量预报方法，具体步骤如下：

s1、在研究区以县界为基准分区标定敏感参数，实现wofost(worldfoodstudies：世界粮食研究)模型和prosail(是耦合prospect叶片模型和sailh冠层结构模型得到的一种辐射传输模型)模型的本地化；

s2、筛选生育期内modis(eos系列卫星上的中分辨率成像光谱仪)反射率产品中的高质量反射率数据，用于输入耦合模型进行数据同化；

s3、耦合prosail和wofost模型，并基于全局敏感性分析确定耦合模型的待优化参数，构建代价函数搭建数据同化框架，以实现区域参数优化；

s4、将ecmwf(欧洲中期天气预报中心)中提供的50组气象预报数据作为wofost模型未来15d气象输入；

s5、利用s4中各组数据的日最高温构成的向量，使用广义向量夹角法计算得到该像元历史年份同时段气象数据中相似度最高的年份，将该年份的气象数据作为wofost模型未来时段的气象输入；

s6、用s5中得到的50组不同的气象输入数据，驱动耦合模型后得到50个不同的预报产量，实现概率预报；

s7、逐格网运行，得到区域作物产量概率预报统计图。

步骤s1所述分区标定敏感参数，方法如下：

将每个县设为一个分区，根据多年历史气象数据，再将产量水平与气候条件均较为相似的县合并为一个分区；在每个分区的基础上，选择敏感的或在空间上具有变异性的参数，分别对wofost模型和prosail模型进行校准。

所述气候条件指气候带、降水量、温度等。

所述参数为tdwi(作物初始干物质量)，span(35℃下叶片存活天数)和slatb(比叶面积)。

所述校准，按常规方法进行，例如按照efast方法校准。

步骤s2所述筛选生育期内modis反射率产品中的高质量反射率数据，具体如下：筛选生育期范围内的modis反射率产品mod09ga，并剔除其中state_1km波段超过阈值的像元，其余像元作为modis反射率观测数据。所述阈值设定为72(此阈值设定为经验值，设定标准是使经过筛选后生育期内仍有5个左右反射率数据可用)。

步骤s3所述耦合prosail和wofost模型的方法如下：

使用lai(叶面积指数)、比叶面积作为耦合变量，wofost模拟的lai和比叶面积(通过lai与叶片干物质重计算得到)作为prosail的输入参数，驱动prosail模型来模拟红光、近红外(nir)和短波红外(swir)波长处的光谱反射率。

步骤s3所述分析确定耦合模型的待优化参数，构建代价函数搭建数据同化框架，具体方法如下：

采用4dvar(四维变分算法)数据同化程序将modis反射率数据和步骤s3耦合得到的wofost-prosail模型的模拟反射结合起来，通过最小化4dvar代价函数来重新初始化四个参数，包括tsum1(出苗期至开花期积温)，tdwi(作物初始干物质量)和span(35℃下叶片存活天数)，lidfa(叶倾角函数系数)(dvs＝0.8，拔节期)；代价函数j(x)的构造如公式(1)所示：

其中，x表示重新初始化参数的向量tsum1、tdwi、span，lidfa，x表示重新初始化参数的tsum1、tdwi、span，lidfa的矢量集合；x^b表示关于这四个参数的先验知识；b是四个参数的误差协方差矩阵；t代表观察时间的集合，t为观测日期，t为t的最大取值；yt代表在观测日期t对特定红光、近红外和短波红外的遥感反射向量；ht表示观测算子，即耦合的wofost-prosail模型；c是平衡观测在同化过程中的影响的一个常数值；q代表不同波长的观测误差协方差矩阵；基于多样本实测地面光谱数据确定：红光、近红外和短波红外波长的观测协方差分别为0.05、0.03和0.04。

步骤s4所述气象预报数据作为wofost模型未来15d气象输入的处理方法如下：

ecmwf提供的气象预报数据中包含七个要素：6h内最高温、6h内最低温、总降水量、露点温度、陆地表面太阳净辐射量、u/v方向(u和v分别代表经度方向和纬度方向，是两个要素)的风速；

其中，wofost气象输入要求的水汽压可用露点温度换算得到，换算公式见公式(2)：

式中，ea为水汽压，t为露点温度；e为自然常数；

其中，wofost所需风速由u/v方向风速矢量合成得到。

其中，s5所述使用广义向量夹角法计算得到该像元历史年份同时段气象数据中相似度最高的年份，具体方法如下：

气象预报数据中日最高温构成的15维向量x＝(x1,x2,…xn)，与历史年份同时段日最高温构成的15维向量y＝(y1,y2,…yn)，n为正整数，其它们的广义夹角θ定义见公式(3)(4)：

即：

求得夹角θ最小的年份为相似度最高的年份。

步骤s6所述用s5中得到的50组不同的气象输入数据，驱动耦合模型后得到50个不同的预报产量，具体如下：

认为50组气象数据在未来15d出现的概率相等，即模型输出的50个不同的预报产量出现的概率相等，以此可得到研究区高产或地产的概率(即高于或者低于某个具体产量的概率)(例如：求某地产量高于6000kg/hm²的概率，且预报产量中高于6000kg/hm²的结果共有18个，则概率为18/50×100％＝36％)。

所述作物为适用于wofost模型和prosail模型的作物，例如小麦、水稻、玉米等大面积种植的作物，优选冬小麦。

本发明还提供所述一种基于作物模型、历史和气象预报数据的作物产量预报方法在指导作物生产中的应用。

本发明与现有技术相比，有益效果为：基于向量相似性的思想建立了历史气象数据与未来气象数据的联系，并结合ecmwf提供的多种气象预报数据实现了概率预报，考虑了未来时段气象不确定性对产量的影响，提高了夏玉米产量预报的时效性和准确性。

附图说明

图1为本发明实施例1的流程示意图。

图2为实施例1中研究区5月16日的区域集合中位数产量预报图与产量高于8000kg/hm²的概率预报统计图。

具体实施方式

下面结合实施例，对本范明的具体实施方式做进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

实施例1

以下面以在2009年河北省衡水市的冬小麦产量预报为例进一步阐述本发明的技术方案，具体流程见图1。

s1、在研究区以县界为基准分区标定敏感参数，实现wofost模型和prosail模型的本地化；

所述区标定敏感参数，方法如下：

将每个县设为一个分区，根据多年历史气象数据，再将产量水平与气候条件均较为相似的县合并为一个分区；在每个分区的基础上，选择敏感的或在空间上具有变异性的参数，分别对wofost模型和prosail模型进行校准；

s2、筛选生育期内modis反射率产品中的高质量反射率数据，用于输入耦合模型进行数据同化；

所述筛选生育期内modis反射率产品中的高质量反射率数据，具体如下：筛选生育期范围内的modis反射率产品mod09ga，并剔除其中state_1km波段超过阈值的像元，其余像元作为modis反射率观测数据；所述阈值设定为72。

s3、耦合prosail和wofost模型，并基于全局敏感性分析确定耦合模型的待优化参数，构建代价函数搭建数据同化框架，以实现区域参数优化；

所述耦合prosail和wofost模型的方法如下：

使用lai、比叶面积作为耦合变量，wofost模拟的lai和比叶面积作为prosail的输入参数，驱动prosail模型来模拟红光、近红外和短波红外波长处的光谱反射率；

所述分析确定耦合模型的待优化参数，构建代价函数搭建数据同化框架，具体方法如下：

采用4dvar数据同化程序将遥感反射率和步骤s3耦合得到的wofost-prosail模型的模拟反射结合起来，通过最小化4dvar代价函数来重新初始化四个参数，包括tsum1，tdwi和span，lidfa；代价函数j(x)的构造如公式(1)所示：

其中，x表示重新初始化参数的向量tsum1、tdwi、span，lidfa，x表示重新初始化参数的tsum1、tdwi、span，lidfa的矢量集合；x^b表示关于这四个参数的先验知识；b是四个参数的误差协方差矩阵；t代表观察时间的集合，t为观测日期，t为t的最大取值；yt代表在观测日期t对特定红光、近红外和短波红外的遥感反射向量；ht表示观测算子，即耦合的wofost-prosail模型；c是平衡观测在同化过程中的影响的一个常数值；q代表不同波长的观测误差协方差矩阵；基于多样本实测地面光谱数据确定：红光、近红外和短波红外波长的观测协方差分别为0.05、0.03和0.04；

s4、将ecmwf中提供的50组气象预报数据作为wofost模型未来15d气象输入，方法如下：

ecmwf提供的气象预报数据中包含七个要素：6h内最高温、6h内最低温、总降水量、露点温度、陆地表面太阳净辐射量、u/v方向的风速；

其中，wofost气象输入要求的水汽压可用露点温度换算得到，换算公式见公式(2)：

式中，ea为水汽压，t为露点温度；e为自然常数；

其中，wofost所需风速由u/v方向风速矢量合成得到；

s5、利用s4中各组数据的日最高温构成的向量，使用广义向量夹角法计算得到该像元历史年份同时段气象数据中相似度最高的年份，将该年份的气象数据作为wofost模型未来时段的气象输入；

所述使用广义向量夹角法计算得到该像元历史年份同时段气象数据中相似度最高的年份，具体方法如下：

气象预报数据中日最高温构成的15维向量x＝(x1,x2,…xn)，与历史年份同时段日最高温构成的15维向量y＝(y1,y2,…yn)，n为正整数，其它们的广义夹角θ定义见公式(3)(4)：

即：

求得夹角θ最小的年份为相似度最高的年份；

s6、用s5中得到的50组不同的气象输入数据，驱动耦合模型后得到50个不同的预报产量，实现概率预报；

所述用s5中得到的50组不同的气象输入数据，驱动耦合模型后得到50个不同的预报产量，具体如下：

认为50组气象数据在未来15d出现的概率相等，即模型输出的50个不同的预报产量出现的概率相等，以此可得到研究区高产或地产的概率；

s7、逐格网运行，得到区域作物产量概率预报结果，经实验得到研究区5月16日的区域集合中位数产量预报图与产量高于8000kg/hm²的概率预报统计图，如图2所示。

虽然，上文中已经用一般性说明及具体实施方案对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。