一种基于耦合作物模型同化光谱反射率的水稻估产方法与流程

本发明涉及一种基于耦合作物模型同化光谱反射率的水稻估产方法，属于农业遥感技术领域。

背景技术：

在作物估产应用过程中，基于光合、呼吸、蒸腾等机理过程的作物模型依靠其内在的物理过程和动力学机制，可以准确模拟作物在时间和空间上的连续变化，可以准确的模拟单点作物的生长状况和产量。但是运用到区域尺度时，由于地表环境的非均匀性，导致模型中的一些参数数据很难以获得。而卫星遥感方法具有空间连续和时间动态变化的优势，能够有效解决这个问题。但是遥感对地观测由于受卫星时空分辨率等因素的制约，还不能真正揭示环境气候条件与作物生长发育内在机理过程之间的关系，但这正是作物模型的优势所在。因此，通过同化的方法将二者结合起来，互相补充，发挥各自的优势，能够更加精确地反映作物的产量。

遥感信息与作物生长模型集成的方法分为驱动法和同化法，其中，同化法受到更多的关注也被更为广泛的采用。优化算法为同化法中的一个重要方法，目前比较通用的方法有以下两种：1、通过对比遥感反演和模型模拟获得的叶面积指数(leafareaindex，lai)；2、对比遥感反演和模型模拟获得的作物光谱信息或者植被指数。如acrm(atwo-layercanopyreflectancemodel)模型，需要与作物生长模型相耦合。

目前，中分辨率成像光谱仪(moderate-resolutionimagingspectroradiometer，modis)作为主要的遥感数据源被广泛应用于区域农作物长势监测和产量估算研究。不过在农田地块较小、分布较为零散、破碎化程度较高的区域，很难保证modis像元是纯像元，此时需要结合更高空间分辨率的遥感影像进行数据同化研究。将时序遥感观测信息加入到耦合模型中，通过调整模型的运行轨迹，使其更接近遥感观测的“客观”状态，模型模拟的lai能够反应作物生长的“真实”情况，使得模拟的作物产量更接近“真实值”。

技术实现要素：

本发明所要解决的技术问题是：提供一种基于耦合作物模型同化光谱反射率的水稻估产方法，将时序遥感观测信息加入到耦合作物模型中，有效改进了农田地块破碎化区域的农作物产量估算结果，使模拟值更接近真实值。

本发明为解决上述技术问题采用以下技术方案：

一种基于耦合作物模型同化光谱反射率的水稻估产方法，包括如下步骤：

步骤1，利用待估算区域的气象、土壤、植被数据对epic作物模型的参数进行适应性调整和修正，并采用待估算区域基准年描述作物品种的遗传参数对epic作物模型进行校验和标定，确定epic作物模型中描述作物品种的遗传参数的取值；

步骤2，获取待估算区域的遥感数据并与待估算区域的地面参数进行空间匹配，构建同化所需的时间序列遥感观测数据集，从数据集中获取地表反射率；所述遥感数据包括modis数据、lai数据、tm数据；

步骤3，利用叶面积指数耦合epic作物模型和acrm辐射传输模型，得到耦合模型，并对耦合模型进行参数敏感性分析；

步骤4，在待估算区域运行经参数敏感性分析后的耦合模型，得到待估算区域的模拟反射率；

步骤5，以模拟反射率和地表反射率之间的观测误差为权重，按以下五个待优化参数建立代价函数，所述五个待优化参数包括种植日期、种植密度、光周期敏感参数、叶片红光、近红外波段反射率，采用模拟退火算法对代价函数进行最小化，当代价函数满足收敛条件时，输出待优化参数的最优值；

步骤6，将待优化参数的最优值代入epic作物模型，计算得到水稻的产量结果。

作为本发明的一种优选方案，步骤1所述描述作物品种的遗传参数包括幼苗期生长特性参数、灌浆期特性参数、光周期敏感参数、出叶间隔特性参数。

作为本发明的一种优选方案，所述步骤2的具体过程如下：

将待估算区域的modis数据与待估算区域的土地利用图进行叠加，判断modis像元中旱地作物所占比例，利用modis数据分别在像元和亚像元尺度提取水稻作物种植面积和空间分布，结合tm数据，提取modis数据分别在像元和亚像元尺度的水稻冠层发射率信息，从而构建同化所需的时间序列遥感观测数据集。

作为本发明的一种优选方案，步骤3所述参数敏感性分析所采用的方法为直方图比较法。

作为本发明的一种优选方案，步骤5所述收敛条件为以下条件中的其中一个：

(1)连续5次循环后待优化参数的取值收缩到指定的范围；

(2)代价函数的目标函数值连续5次循环后提高的百分比小于0.0001％；

(3)计算代价函数的次数超过10000次。

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

本发明将时序遥感观测信息加入到耦合模型中，通过调整模型的运行轨迹，使其更接近遥感观测的“客观”状态，模型模拟的lai能够反应水稻生长的“真实”情况，使得模拟的水稻产量更接近“真实值”。

附图说明

图1是本发明一种基于耦合作物模型同化光谱反射率的水稻估产方法的流程图。

图2是采用本发明方法估算的水稻总产及单产示意图，其中，(a)为水稻单产，(b)为水稻总产量。

具体实施方式

下面详细描述本发明的实施方式，所述实施方式的示例在附图中示出。下面通过参考附图描述的实施方式是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。

以估算江苏省东南部的东台市研究区的水稻产量为例进一步阐述本发明的技术方案。本实施例的耦合作物模型同化光谱反射率的水稻估产方法的流程如图1所示，包括：

步骤s1，利用当地气象、土壤、植被等数据对epic(environmentalpolicy-integratedclimate)模型参数进行适应性调整、修正及模拟结果校验，以及对模型进行标定，确认模型中描述作物品种的遗传参数取值。

在进行敏感性分析前需要采集数据，选择江苏省东南部的东台市为研究区域，地理坐标位于东经120°07′～120°53′，北纬32°33′～32°57′，属于北亚热带南部季风气候区，土壤类型主要有黄白土、黄泥土、淤泥土、青泥土。获得以下数据：2012年30m分辨率的土地利用数据；2010、2011、2012三年包括日最高气温，日最低气温，降水量和日照时数在内的所需气象要素；从modis数据下载网站获取到2010、2011和2012年的陆地表面反射率产品mod09a1，空间分辨率为500m。将含有地理位置信息的数据统一坐标，完成空间匹配；以及2010、2011和2012三年水稻的地面实测数据。

将epic模型应用到区域尺度，需要对其进行区域校准，即对特定地区确定代表性品种的过程。考虑到一定地理范围内，受温度和日照条件影响，同种作物的种植习惯及生长发育过程均固定在一定时间段内。比如东台市水稻种植日期一般在5月中上旬左右、7月下旬左右开始吐丝、8月中上旬左右进行灌浆、10月上旬后收割，这就是一般年份的水稻生育时期。

在epic模型中，品种遗传参数：幼苗期生长特性参数(p1)、灌浆期特性参数(p5)、光周期敏感参数(p2)以及出叶间隔特性参数(phint)这四个参数决定了水稻的生育期，潜在灌浆速率参数(g3)和单株最大籽粒数(g2)影响作物的产量。以2010年为基准年，将种植日期、幼苗期生长特性参数(p1)、光周期敏感参数(p2)、灌浆期特性参数(p5)和出叶间隔特性参数(phint)分别在各自取值区间内均匀取值采样10000次，然后将不同的参数取值组合代入epic模型进行模拟，提取在水稻合理生育时期的参数组合，经过反复试验，确定出东台市水稻生育时期的遗传参数取值范围。将主推水稻品种最大籽粒数的平均值作为epic模型中单株最大籽粒数(g2)的取值，结合2010年实测产量数据，在幼苗期生长特性参数(p1)、光周期敏感参数(p2)、灌浆期特性参数(p5)和phint取区间中值、单株最大籽粒数(g2)固定的情况下，模拟水稻产量，然后和地面实测产量相比较，当模拟和实测产量基本吻合时，潜在灌浆速率参数(g3)的取值作为该参数的校准值。利用校准的epic模型分别预测2010、2011和2012年水稻lai和产量，比较模拟结果与水稻地面观测数据，评价模型参数取值的可靠度。

步骤s2，获取遥感数据并与地面参数进行空间匹配，构成时间序列遥感观测数据集。将土地利用图与modis数据进行叠加，判断modis像元中旱地作物所占比例。利用modis数据分别在像元和亚像元尺度提取水稻作物种植面积和空间分布，结合可用的tm遥感观测，提取modis像元和亚像元的水稻冠层反射率信息，构成时间序列遥感观测数据集。

从modis数据下载网站获取2010、2011和2012年陆地表面反射率产品mod09a1空间分辨率为500m生长季内获取了每年10次遥感观测数据，涉及到的年积日分别为168、178、186、198、202、208、220、232、243、256，这些日期均为mod09a1标识的观测日期。此外，下载了mcd43a4产品数据(校正到天顶观测的反射率数据)，空间分辨率为500m。采用mod09a1数据构建了时间序列遥感观测数据集，所用反射率为mod09a1红光和近红外波段反射率，用来进行同化估产；nbar数据主要用来判断modis像元是否被水稻作物覆盖，已校正到天顶观测，不存在因观测角度差异引起的地表反射率的变化；2010年doy202、243和256处tm数据主要连同mod09a1数据一起构成同化所需的时间序列遥感观测数据集。

步骤s3，通过lai耦合epic作物模型和acrm辐射传输模型，并且进行参数敏感性分析。

由于epic模型和acrm模型的结合点是lai，敏感性分析的主要目的是试图找出影响水稻lai随时间变化的关键参数，以及影响作物冠层反射率的相关参数，将这些参数作为同化过程中模型的重新初始化变量。本发明采用直方图比较法对耦合模型的作物参数进行全局敏感性分析。通过以下式子来计算敏感指数：

式中n表示montecarlo估算中的总样本数量。

计算得出输出结果如式(3)所示，并且通过式(4)、(5)、(6)计算敏感性指数：

ynor＝[y⁽¹⁾-min,y⁽²⁾-min,…,y⁽ⁿ⁾-min]^t(max-min)(4)

式中，ynor、代表了标准化的输出结果，min、max表示y的最小、最大值，min′、max′表示的最小、最大值，s代表敏感性指数，代表ynor频率直方图面积，代表ynor和直方图重叠部分面积。

步骤s4，在实验地区内运行s3中的耦合模型，得到当地的模拟反射率，以及预处理modis数据得到的高质量地表反射率。

步骤s5，按五个待优化参数建立代价函数，并以s4中得到的模拟反射率和高质量地表反射率之间的观测误差为权重求和得到总代价函数，对总代价函数进行最小化，使其不断重复初始化s3获得的需要同化的耦合模型的敏感性参数，在初始参数空间中搜索全局最优解，使得代价函数快速收敛，最终在收敛条件被满足时，同化过程结束，输出待优化参数最优值。

代价函数的表达式如下：

式中，x代表自由变量，i代表时间序列，l代表modis反射率总观测时间长度，yi代表来自mod09a1产品的地表反射率，hi代表耦合作物模型，α是除了x以外的其他输入参数，ri为观测误差协方差矩阵。

式中，x′为acrm模型中的自由变量，m为辐射传输模型acrm，α′为acrm模型中除了x′以外的其他输入参数。

待优化参数包括种植日期、种植密度、光周期敏感参数、叶片红光和近红外波段反射率。

选用vfsa(模拟退火)优化算法求解代价函数，优化成功后与最小代价函数值对应的种植日期、种植密度、光周期敏感参数(p2)、叶片红光和近红外波段反射率等5个参数取值称为“最优值”。将种植日期、种植密度和光周期敏感参数(p2)的“最优值”，连同其他输入参数一起输入到epic模型，模拟得到数据同化后的产量结果。

vfsa优化算法的收敛条件为以下条件之一：

(1)5次循环后待优化参数值已经收缩到指定的范围；

(2)目标函数值在五次循环后无法提高0.0001％；

(3)计算代价函数次数超过10000次。

本发明提出了基于多时相遥感观测地表反射率数据的区域水稻产量同化估算方案，分别在modis和tm尺度估算了水稻产量，有效改进了农田地块破碎化区域的水稻产量估算。利用江苏省东台市2010、2011和2012年的数据对该方案进行了验证，初步分析了遥感数据在年内产量及年际间产量差估测中的作用。如图2的(a)和(b)所示，东台市2010、2011和2012年遥感数据同化估算的水稻总产量误差分别为8.78％、13.59％和9.05％，平均单产误差分别为13.52％、7.53％和2.45％；与仅利用epic模型模拟得到的产量误差相比，2010、2011和2012年遥感估算总产量的误差分别减小了6.54％、2.21％和5.87％，平均单产误差分别减小了6.27％、2.35％和5.76％，表明经遥感数据同化后的产量优于仅运行epic模型模拟的结果。

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。