CN109829234B

CN109829234B - 一种基于高分辨率遥感数据和作物模型的跨尺度高精度动态作物长势监测和估产方法

Info

Publication number: CN109829234B
Application number: CN201910093078.5A
Authority: CN
Inventors: 张朝; 曹娟; 陶福禄; 骆玉川; 李子悦
Original assignee: Beijing Normal University
Current assignee: Beijing Normal University
Priority date: 2019-01-30
Filing date: 2019-01-30
Publication date: 2019-09-27
Anticipated expiration: 2039-01-30
Also published as: CN109829234A

Abstract

本发明公开了一种基于高分辨率卫星遥感数据和作物模型的跨尺度高精度动态作物长势监测和估产方法，包括：实现模型的本地化；利用GEE平台对遥感数据与地面数据进行空间匹配；设计多种模拟情景；将作物的生育期以返青期为节点分为前后两个时间窗口，并计算全生育期各种气象要素；构建回归方程，建立所有生育期内每天的回归方程；逐像元提取每年前后两个时间窗口的卫星遥感观测植被指数最大值及其对应的日期,并将提取到的植被指数最大值通过经验公式转化为回归方程的自变量LAI；以及以提取前后两个时间窗口观测日期为参照，利用对应组合日期的回归方程进行逐像元计算，在完成所有像元上的运行后，获得作物模拟产量。

Description

一种基于高分辨率遥感数据和作物模型的跨尺度高精度动态作物长势监测和估产方法

技术领域

本发明涉及农业遥感技术领域，更具体涉及一种基于高分辨率遥感数据和作物模型的跨尺度高精度动态作物长势监测和估产方法。

背景技术

传统的农作物估产方法主要分为基于农学统计数据、基于农业气象统计数据、基于经济统计数据和基于作物模型的预报方法等。这些方法由于高昂的时间、人力和物力以及作物模型本身的局限性都难以实现区域作物的高精度估产，尤其是田块尺度的估产。遥感技术具有快速、宏观、动态三个重要特点，可以快速获取大范围的地表信息，在大范围作物生长监测和产量预报发挥着巨大的作用，在区域估产和预测上具有得天独厚的优势，将遥感数据应用到作物估产模型上能够实现区域估产的目的。

国内外遥感作物估产模型主要分为遥感植被指数统计估产模型、光能利用率模型、遥感和作物生长耦合模型。遥感数据统计估产模型主要以遥感数据为基础，利用统计学的相关分析方法，建立遥感波段或波段组合数据与作物产量之间的关系，进而实现区域农作物估产；光能利用率模型主要是利用遥感估算生物量或作物的净初级生产力，加上修正的收获指数，进而得到作物的单产估算结果。目前遥感和作物生长耦合模型主要是以同化技术为主，将利用遥感技术获取的关键生物物理化学参数叶面积指数(LAI)、土壤水含量、叶片氮磷含量和地面生物量(AGB)等参数与基于作物光合、呼吸、蒸腾、营养等机理过程的作物模型结合进而达到区域高精度作物估产的目的，遥感与作物生长模型的耦合是目前国内外农业遥感估产的研究热点，这类模型能够结合作物模型的点模拟和遥感观测在面上的模拟优势。

基于遥感数据的统计估产模型虽然简单易用，但是局限于特定的时间，地点和年份，不具有时空外推性，且缺乏对作物生长的机理性认识。光能利用率参数模型虽对作物生长机理有一定的涉及，但部分参数量化(如光能利用效率、收获指数等)需要进一步加强研究；此外，光能利用率模型很难考虑到人为因素的影响，这在一定程度上制约了模型在农作物估产的精度。利用遥感数据同化模型虽然能够实现遥感数据和点尺度作物生长模型优势结合，但是存在着输入数据繁杂且难以获取，计算成本高和算法不确定性的缺点，该类研究主要利用中低分辨率遥感数据，研究结果难以实现田块间的估产研究。

因此，需要新的技术以至少部分解决现有技术中存在的局限。

发明内容

为解决现有作物估产存在的地面实测数据需求高、难以获取，模型空间应用能力较弱、同化算法不确定性和计算成本高等问题，本发明提供了一种基于Google EarthEngine(GEE)和作物生长模型高精度的农业估产方法，基于校准的CERES系列作物模型，根据当地实际的田间管理措施和作物环境设置不同的模拟情景，结合高时空分辨率遥感数据，构建产量与植被指数和气象要素的模型，逐像元模拟完成由点到面的产量估算。

根据本发明的一方面，提供一种基于高分遥感数据和作物模型的跨尺度高精度动态作物长势监测和估产方法，包括如下步骤：

S1:选择一定的研究区，对作物模型进行参数优化和验证，实现模型的本地化；

S2：利用GEE平台对遥感数据进行预处理，包括遥感数据的大气矫正、辐射矫正和去云处理；

S3:获取研究区内多年田间管理措施记录，依据该记录设计多种模拟情景，所述模拟情景包括种植日期、种植密度、灌溉和施肥；

S4：将作物的生育期以返青期为节点分为前后两个时间窗口，并且计算全生育期各种气象要素，所述气象要素包括研究区逐栅格点全生育期内的平均温度，平均降雨量和累积太阳辐射量；

S5:基于设定的田间管理模拟情景，多次运行作物生长模型以获得模拟的多条LAI曲线和每条LAI曲线对应的作物产量；以每次模拟的产量为因变量，以所述前后两个时间窗口中各取任意一个LAI以及生育期气象要素为自变量，构建回归方程，其中包括建立生育期内前后两个时间窗口内所有组合的回归方程，并存储全部的回归方程和相关系数；

S6：逐像元提取每年前后两个时间窗口的卫星遥感观测植被指数最大值及其对应的日期,并将提取到的最大植被指数通过经验公式转化为S5步骤中回归方程的自变量LAI；

S7:以S6步骤中提取前后两个时间窗口观测日期为参照，利用S6步骤中对应组合日期的回归方程进行逐像元计算，在完成所有像元上的运行后，获得区域作物模拟产量。

根据本发明的一个实施方案，其中步骤S1中所述作物模型为CERES(cropestimation through resource and environmental synthesis)作物模型。

根据本发明的一个实施方案，其中所述对作物模型进行参数优化和验证包括：

采集研究区内站点尺度的土壤参数、气象参数、生育期参数和田间管理参数；

对品种参数(遗传系数)进行逐步校准:(i)根据品种所需积温来优化影响开花和成熟的遗传系数；(ii)校准决定产量的遗传系数；

采用均方根误差(RMSE)和相对均方根误差(RRMSE)对生育期和产量的模拟值与实测值的偏差进行估计，其中所述均方根误差的计算公式为式(1)，相对均方根误差的计算公式为式(2)：

其中，O_i和S_i分别为观测值和模拟值,O_avg为观测的平均值；n为样本量。

根据本发明的一个实施方案，其中步骤S2中，所述大气矫正利用GEE平台的FLAASH模型进行，所述辐射矫正采用陆地卫星生态系统干扰自适应处理模型(landsat ecosystemdisturbance adaptive processing system,LEDAPS)进行，将辐射值转化为反射率；所述去云处理利用GEE平台的云概率计算(Simple Cloud Score)算法得出云概率，其中，所述云概率计算的公式如下:

其中，b1、b2、b3、b4、b5和b7分别代表遥感影像的蓝波段、绿波段、红波段、近红外波段、短波红外1波段和短波红外2波段，thermal表示热红外波段,clamp函数是指返回该像元是云的概率的函数，其中,返回的概率在0-1之间，大于0.2的像元被认为是云，需要将该类像元去掉。

根据本发明的一个实施方案，其中步骤S5中，所述回归方程为：

Yield＝β₀+β₁×W+β₂LAI₁+β₃×LAI₂ (4)

其中Yield代表每次模型模拟的产量；W代表生育内天气向量属性，包括降雨，温度和太阳辐射；d代表生育期内的日期，LAI₁代表前时间窗口所取得叶面积指数值，LAI₂代表后时间窗口所取得叶面积指数，，β₁，β₂和β₃分别表示天气因子、前时间窗口LAI和后时间窗口LAI对产量的贡献系数，β₀代表误差项。

根据本发明的一个实施方案，其中步骤S6中，所述植被指数是指遥感观测的宽动态范围植被指数(WDRVI)。

根据本发明的一个实施方案，其中所述WRDVI的计算公式为：

其中ρ_NIR为传感器载荷所测得的近红外波段反射率，ρ_red为红光波段反射率。

根据本发明的一个实施方案，其中所述经验公式为：

WDRVI＝-0.681+1.437(1-e^-0.351LAI) (6)

本发明结合高分辨率遥感观测数据例如哨兵2号(sentinel-2)，气象数据和作物模型，创建了一种由点扩展到区域的作物估产新方法，该方法可以通过大致了解该区域的种植环境和田间管理措施，不需要大量的地面实测数据，模型易操作且空间应用能力强，不受地面实测数据的难以获取的限制，能够估算县级甚至田块间的收益，可动态评估产量在时间和空间上的变异，为开展田块尺度的产量差异原因的研究提供数据支持。

附图说明

附图中相同的附图标记标示了相同或类似的部件或部分。本发明的目标及特征考虑到如下结合附图的描述将更加明显，附图中：

图1是根据本发明一个实施方案的基于高分辨率遥感数据和作物模型的跨尺度高精度动态作物长势监测和估产方法流程示意图。

图2为根据本发明一个实施方案的利用GEE平台对遥感数据进行处理获得的作物生长曲线图；

图3为根据本发明一个实施方案的模型模拟值与实测值对比图；

图4为根据本发明一个实施方案的模型模拟的100个不同管理情景下的逐日叶面积指数图；

图5为根据本发明一个实施方案的GCVI转化LAI的经验公式图。

图6为根据本发明一个实施方案的安徽省定州市估算的小麦产量结果图。

具体实施方式

为清楚的说明本发明中的方案，下面给出优选的实施例并结合附图详细说明。以下的说明本质上仅仅是示例性的而并不是为了限制本公开的应用或用途。

应该理解的是，本发明所引用的作物模型和遥感影像模型本身是已知的，例如模型的各个子模块、各种参数、运行机制等等，因此本发明重点阐述作物模型和遥感数据耦合，提供了一种结合各类模型的优点，更具有估产效率和精度优势的作物产量遥感估算方法。克服了现有技术中存在的缺陷。

图1是根据本发明一个实施方案的基于高分辨率遥感数据和作物模型的跨尺度高精度动态作物长势监测和估产方法流程示意图。如图1所示主要的模块包括情景的模拟①、气象数据的处理②、模型的建立③、遥感数据的提取④和逐像元产量的估算⑤，本发明的方法具体实施步骤如下：

S1:对作物模型进行参数优化和验证，实现模型的本地化，使其具备可靠的模拟能力。

S1步骤中所述的模型优选CERES(crop estimation through resource andenvironmental synthesis)系列作物模型，涉及的作物例如可以是大豆、小麦、水稻、玉米、花生等等。当然也可以采用其他适当的作物模型。

其中，当采用CERES作物模型时，采集研究区站点尺度的土壤参数、气象参数、生育期参数和田间管理参数，对品种参数进行逐步校准，例如可以基于CERES-DSSAT模型自带的GLUE工具来进行，或者利用其它工具:(i)根据品种所需积温优化影响开花和成熟的遗传系数；(ii)校准决定产量的遗传系数。验证可以采用均方根误差(RMSE)和相对均方根误差(RRMSE)对生育期和产量的模拟值和实测值的偏差进行估计。

其中所述均方根误差的计算公式(1)：

相对均方根误差的计算公式(2)：

Oi和Si分别为观测值和模拟值,O_avg为观测的平均值；n为样本量，大于1的自然数。

其中，所述的土壤参数可以包括土壤类型、颜色、坡度、渗透性、反射率、土层厚度、土壤水分蒸发限制(mm)、径流曲线数和土壤排水率(fraction day^-1)、光合作用因子(0-1)、土壤水下限或调萎点含水量(cm³cm^-3)、田间持水量(cm³cm^-3)，饱和含水量(cm³cm^-3)、土壤容量(g cm^-3)、土壤有机碳、氮(wt.％)、土壤PH值、粘粒含量(wt.％，<0.002mm粒径)和粉粒含量(wt.％，0.002-0.05mm粒径)。

其中，所述的气象参数包括日太阳辐射(MJ m^-2)19、日最高温度(℃)、日最低气温(℃)和日降雨量(mm)。

其中，所述的田间管理措施包括站点尺度的灌溉、品种、施肥、种植密度、播种方式。

其中，生育期参数包括播种到开花日期的间隔天数，播种到成熟日期的间隔天数。

其中，遗传系数即为品种参数，不同作物有不同的品种参数；如小麦有7个参数,玉米有6个参数，大豆有15个参数。玉米品种参数包含的6个与生长有关参数中，有2个控制玉米生长阶段(即P1：营养生长期积温，P5：抽雄-生理成熟生长的积温),有2个数控制产量(例如G2：单株玉米的籽粒数，G3:灌浆速率)。

本领域技术人员应该理解，上述土壤参数、气象参数、田间管理参数、生育期参数以及品种参数等本身为本领域技术人员所熟知，因此不做进一步的说明。

S2：利用GEE平台对遥感数据与地面的数据进行空间匹配，包括遥感数据的大气矫正、辐射矫正和去云去阴影处理。

S2步骤中所说的大气矫正可以采用的GEE平台自带的FLAASH模型，辐射矫正可以采用陆地卫星生态系统干扰自适应处理模型(landsat ecosystem disturbance adaptiveprocessing system，LEDAPS)，将辐射值转化为反射率；去云处理利用是GEE平台的云概率计算(SimpleCloudScore)算法，该算法是利用亮度、温度和归一化差分雪指数的组合计算出的每个像元是云概率。

其中，所述云概率计算的计算公式:

cloub_score＝clamp(min((b1-0.1)/0.2，(b1+b2+b3-0.2)/0.6，

(b4+b5+b7-0.3)/0.5，(300-thermal)/10，

(0.8-(b2-b5)/(b2+b5))/0.2，0.0，1.0))

(3)

其中，b1、b2、b3、b4、b5和b7分别代表遥感影像的蓝波段、绿波段、红波段、近红外波段、短波红外1波段和短波红外2波段，thermal表示热红外波段，clamp函数是指返回该像元是云的概率函数，其中，返回的概率在0-1之间，大于0.2的像元被认为是云，需要将该类像元去掉

S3:依据当地农气站记录的多年田间管理措施，设计不同的种植日期，种植密度，灌溉和施肥等模拟情景；

S4：将作物的生育期以返青期为节点分为前后两个时间窗口，并且计算全生育期各种气象要素。所述气象要素可以包括研究区逐栅格点全生育期内的平均温度(℃)，平均降雨量(mm)和累积太阳辐射量(MJ m^-2)。例如，可以利用地面气象站的观测记录数据为基础来进行插值计算。

S5:以每次模拟的产量为因变量，前后窗口任意两个日期的LAI以及生育期气象要素为自变量构建回归方程；

S5所说的回归方程的计算公式可以为：

Yield＝β₀+β₁×W+β₂LAI₁+β₃×LAI₂ (4)

Yield代表每次模型模拟的产量；W代表生育内天气向量属性(降雨，温度，太阳辐射)，由S3步骤获得；LAI₁代表前时间窗口所取得叶面积指数值，LAI₂代表后时间窗口所取得叶面积指数，，β₁,β₂和β₃分别表示天气因子、前窗口LAI和后窗口LAI对产量的贡献系数，β₀代表误差项。因为每个像元获取的LAI最大值影像日期不同，因此要建立所有生育期内每天的回归方程，最后将这些回归方程和相关系数储存起来，建立回归方程组并由此确定上述各项系数。

S6：提取逐像元提取每年前后两个时间窗口的卫星遥感观测植被指数最大值及其对应的日期,并将提取到的到植被指数通过经验公式转化为S5步骤中自变量LAI。

其中，S6所说的植被指数指的是遥感观测的各种植被指数，如归一化植被指数(NDVI)、增强性植被指数(EVI)和宽动态范围植被指数(WDRVI)等，如WRDVI的计算公式：

其中，S6步骤中所述的经验公式，主要是利用地面实测LAI与对应空间位置植被指数建立相关方程，本实施例中利用WDRVI与LAI的转化公式：

WDRVI＝-0.681+1.437(i-e^-0.351LAI) (6)

S7:以S6中提取前后窗口观测日期为参照，利用S5步骤中的对应组合日期的回归方程进行逐像元计算，在完成所有的像元上的运行后，最后将模拟产量综合到县级尺度进行精度评价，指导作物的生产。

实施例

下面以小麦为例，示例性阐述本发明的方法的具体实施过程。

步骤S1，选择华北平原安徽定州地区冬小麦主产区为研究区，在东经114°48'－115°15',北纬38°14'－38°40'之间,总面积127,500ha。地势平坦，沙河、孟良河和唐河横贯全境，主要作物为冬小麦、夏玉米、棉花和一些其他作物。该区以黄褐土为主，属半干旱季风气候，年平均温度和降水量分别为12.4℃和503.2mm，无霜期170－190d，常年积温4200－4800℃。

步骤S2，利用GEE平台对2015-2018年哨兵-2影像数据预处理，包括大气矫正、辐射矫正和集合矫正、去云处理等，并提取每年多时相的NDVI数据，构建SG滤波方法分析研究区夏玉米和冬小麦的NDVI时序特征曲线。通过分析作物的当地作物生育期和冬小麦的生长曲线(图2)，依据以下原则提取冬小麦种植面积：冬小麦5月中下旬处于乳熟期，NDVI最大值超过0.5；6月中旬冬小麦已收割，夏玉米刚播种，NDVI接近0.1；7月中旬，夏玉米处于吐丝乳熟期，NDVI值达到最大值超过0.7；8月下旬到9月夏玉米已收割,此时NDVI小于0.1。

步骤S3，基于2008-2013年冬小麦的生育期和产量数据，对研究区的农气站进行调参，调参结果表明CERES-Wheat可以很好地捕捉冬小麦的开花(ADAT)和成熟日期(MDAT)，平均RMSE(均方根误差)小于10天，RRMSE(相对均方根误差)低于6.5％；。产量的绝对误差为459.1kg/ha，RMSE为8.2％，具体参见图3。

步骤S4，基于校准的CERES-Wheat模型，根据实际管理实践设置不同播种日期和种植密度的组合生成一系列模拟情景。播种日期为9月15日至10月15日，间隔3天，种植密度在245至965×1000株/ha间随机产生10个，收集每年100个管理情景下的产量和逐日LAI，由此可以动态监测。

步骤S5,将小麦生育期以返青期为节点分为前后各30天的两个时间窗口(参见图4)，并且利用地面气象站的观测记录来计算全生育期的平均温度和辐射以及总降雨量。图4为根据本发明一个实施方案的模型模拟的100个不同管理情景下的逐日叶面积指数图，图中横纵坐标分别为日期和LAI，阴影部分为返青期节点(峰值)前后各30天，作为两个时间窗口。然后，利用这些LAI以及生育期气象要素值为自变量，构建回归方程。

步骤S6，利用GEE平台，提取每年每个像元前后两个时间窗口的最大WDRVI值及其日期并将WDRVI转换为LAI，经验关系如图5所示。

步骤S7，选取对应组合日期的回归关系进行逐像元计算，最后将模拟产量综合到县级尺度进行精度评价,估计的平均产量为5678.5千克/公顷(2016-2018年模拟产量分别见图6)，分别占2016年和2017年产量的76％和90％(表1)。

表1.估算和观测的产量

¹河北经济年鉴统计数据.

本实施例基于GEE平台，在像元尺度上结合哨兵-2影像和CERES-Wheat模型进行10m分辨率的冬小麦产量制图。估产结果精度较高且基本上可以表征产量的空间差异性。证明该方法可以很好地反映作物对气候和田间管理的响应，为作物产量以及灾害损失评估的业务化运行提供了新思路。

相比于现有技术，本发明实现了如下创新：

1、本发明提出了由点扩展到面的高精度动态农业估产方法技术。

2、充分利用了高分辨率遥感信息和作物生长作物模型信息的融合，结合辅助的背景信息，如气象数据、农气站数据，形成从宏观到微观的作物观测能力，相对于数据同化，该方法更加简单、可行且有效。

3、该方法不局限于作物估产，还可以用于农作物种子格局监测、作物长势监测、农作物灾害评估等各方面的农业应用。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的装置及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims

1.一种基于高分遥感数据和作物模型的跨尺度高精度动态作物长势监测和估产方法，包括如下步骤：

S2：利用GEE平台对遥感数据与地面数据进行数据预处理，包括遥感数据的大气矫正、辐射矫正和去云处理；

S3：获取研究区内多年田间管理措施记录，依据该记录设计多种模拟情景，所述模拟情景包括种植日期、种植密度、灌溉和施肥；

S5:基于设定的田间管理模拟情景，多次运行作物生长模型以获得模拟的多条LAI曲线和每条LAI曲线对应的作物产量，其中LAI为叶面积指数；以每次模拟的产量为因变量，以所述前后两个时间窗口中各取任意一个LAI以及生育期气象要素为自变量，构建回归方程，其中包括建立生育期内前后两个时间窗口内所有组合的回归方程，并存储全部的回归方程和相关系数；

S6：逐像元提取每年前后两个时间窗口的卫星遥感观测植被指数最大值及其对应的日期，并将提取到的植被指数最大值通过经验公式转化为S5步骤中回归方程的自变量LAI；

S7:以S6步骤中提取前后两个时间窗口观测日期为参照，利用S5步骤中对应组合日期的回归方程进行逐像元计算，在完成所有像元上的运行后，获得区域作物模拟产量；

其中，步骤S5中，所述回归方程为：

Yield＝β₀+β₁×W+β₂LAI₁+β₃×LAI₂ (4)

其中Yield代表每次模型模拟的产量；W代表生育内天气向量属性，包括降雨，温度和太阳辐射；LAI₁代表前时间窗口所取得叶面积指数，LAI₂代表后时间窗口所取得叶面积指数，β₁，β₂和β₃分别表示天气因子、前时间窗口LAI和后时间窗口LAI对产量的贡献系数，β₀代表误差项。

2.根据权利要求1所述的方法，其特征在于，步骤S1中所述作物模型为CERES(cropestimation through resource and environmental synthesis)作物模型。

3.根据权利要求2所述的方法，其特征在于，所述对作物模型进行参数优化和验证包括：

对品种参数进行逐步校准:(i)根据品种所需积温来优化影响开花和成熟的遗传系数；(ii)校准决定产量的遗传系数；

其中，O_i和S_i分别为观测值和模拟值，O_avg为观测平均值；n为样本量。

4.根据权利要求1所述的方法，其特征在于，步骤S2中，所述大气矫正利用GEE平台的FLAASH模型进行，所述辐射矫正采用陆地卫星生态系统干扰自适应处理模型(landsatecosystem disturbance adaptive processing system，LEDAPS)进行，将辐射值转化为反射率；所述去云处理利用GEE平台的云概率计算(Simple Cloud Score)算法得出云概率，其中，所述云概率计算的公式如下：

其中，b1、b2、b3、b4、b5和b7分别代表遥感影像的蓝波段、绿波段、红波段、近红外波段、短波红外1波段和短波红外2波段，thermal表示热红外波段，clamp函数是指返回该像元是云的概率的函数，其中，返回的概率在0-1之间，大于0.2的像元被认为是云，将所述大于0.2的像元去掉。

5.根据权利要求1所述的方法，其特征在于，步骤S6中，所述植被指数是指遥感观测的植被指数，包括归一化植被指数(NDVI)、增强性植被指数(EVI)和宽动态范围植被指数(WDRVI)。

6.根据权利要求5所述的方法，其特征在于，所述宽动态范围植被指数的计算公式为：

其中ρ_NIR为传感器载荷所测得的近红外波段反射率，ρ_red为红光波段反射率，α＝0.2。

7.根据权利要求5所述的方法，所述宽动态范围植被指数转化为LAI的经验公式为：

WDRVI＝-0.681+1.437(1-e^-0.351LAI) (6)。