发明内容
为了解决现有技术中存在的至少一个技术问题,本申请提供了一种基于时序遥感数据的作物单产及长势评估方法。
第一方面,本申请公开了一种基于时序遥感数据的作物单产评估方法,包括如下步骤:
步骤1、构建作物单产预测模型
步骤101、获取研究区域内多个年份范围内的时序遥感数据(例如MODIS遥感数据);
步骤102、对时序遥感数据进行归一化植被指数提取,以构建NDVI时序数据集;
步骤103、提取植被指数特征值及其对应的时间节点:结合历年时序遥感数据提取作物生长过程特征曲线,对播种期、成熟期最低谷值和生长期最高峰值时点进行位置识别,对应时段分别记为
并提取相应NDVI
min_s、NDVI
min_m、NDVI
max_g值;并基于已经提取的最低谷值和最高峰值,分别用算术平均法得到生长前期中值NDVI
mid_gf,其对应时段记为
和生长后期中值NDVI
mid_gb,其对应时段记为
步骤104、计算不同时间点的累积植被指数变化率
先计算各个时间节点的植被指数变化率;然后依次计算
到
至
之间任一时间节点的累积植被指数变化率;
其中,植被指数变化率和累积植被指数变化率通过如下公式获得:
式中,CRVIi为i时段的植被指数变化率,NDVIi和NDVIi-1分别为i时段和i-1时段的植被指数,i表示时序产品对应的时间节点,i-1为i紧邻的前一个时间节点;ACRVIi为截止i时间节点的累积植被指数变化率,t0为播种期最低值出现时间,T为特定时间节点;
步骤105、获取所述多个年份所对应的作物单产统计数据,以形成作物单产统计数据集;
步骤106、将所述累积植被指数变化率时序产品集和对应的作物单产统计数据集进行回归分析,得出不同特定时间节点上针对研究区域作物单产预测模型;
步骤2、单产评估
获取待评估年份的遥感数据,根据步骤103的方法提取植被指数特征值及其对应的估产时间节点,根据步骤104的方法计算待评估年份不同时间节点的植被指数变化率及累积植被指数变化率,将其代入步骤106中得出的对应时间节点的作物单产预测模型中,经计算输出目标作物单产评估结果。
根据本申请的至少一个实施方式,如果估产时没有完整的作物生长季相关遥感数据,采用历史多年对应时间节点的有效均值特征影像作为未结束年份的准同期数据,随着作物生长期发展,进行数据的迭代更新。
需要说明的是,根据本发明的思路,采用的遥感数据源除MODIS外,还可以是其他卫星遥感数据;另外,除归一化植被指数(NDVI)外,可以采用其它反映作物生长状况的特征参量数据,例如叶面积指数(LAI)、增强植被指数(EVI)。
根据本申请的至少一个实施方式,在所述步骤102中,对MODIS遥感数据进行归一化植被指数提取之前,还包括:
对MODIS遥感数据进行预处理,所述预处理包括拼接处理、剪裁处理以及质量控制处理。
根据本申请的至少一个实施方式,在所述步骤102中,对预处理后的MODIS遥感数据进行归一化植被指数提取之前,还包括:
利用中高分辨率遥感数据进行监测区域目标作物识别,得到目标作物分布数据;
利用目标作物分布数据对预处理后的MODIS遥感数据进行掩膜处理。
根据本申请的至少一个实施方式,在所述步骤102中,对掩膜处理的MODIS遥感数据进行归一化植被指数提取之后,还包括:
进行滤波处理,本发明的实施方案中优选采用Savitzky-Golay滤波处理。
根据本申请的至少一个实施方式,数据时间步长为8d,所述多个年份指8年以上。
根据本申请的至少一个实施方式,在所述步骤2中进行目标作物估产时,目标作物的估产时间节点选自
其中,m取1-9之间自然数。
第二方面,本申请还公开了一种基于时序遥感数据的作物长势评估方法,包括如下步骤:
步骤301、采用上述第一方面中任一项所述的基于时序遥感数据的作物单产评估方法,得到目标作物的单产评估结果;
步骤302、根据得到的单产评估结果,结合单产等级划分标准进行单产等级评估,其中,单产等级包括高产、中等水平、低产;
步骤303、根据单产等级划分标准和产量预测模型测算各特定时间节点长势植被指数特征参量的等级标准后,评估待评估区域目标作物的长势。其中,单产等级的高产、中等水平、低产分别对应长势等级的好、中等、差。
根据本申请的至少一个实施方式,在所述步骤302中,单产等级以多年单产均值μ和标准差σ确定,将(μ-0.5σ)≤评估产量≤(μ+0.5σ)定级为中等水平,评估产量大于(μ+0.5σ)定级为高产,评估产量小于(μ-0.5σ)定级为低产。
本申请的创新及有益技术效果:
本发明提出了一个新的作物估产与长势监测评估的特征参量:累积植被指数变化率。应用试验显示,该参量较累积NDVI、累积植被变化速率在作物估产和长势评估方面更具优势;另外,确立了基于产量的相对客观且应用便捷的长势遥感监测评估模式。
具体创新点如下:
1)本申请的基于时序遥感数据的作物单产及长势评估方法中,以累积植被指数变化率作为作物估产与长势监测评估的特征参量,一方面利于反映作物在特定生长时段内生物量累积变化状况,另一方面可以降低作物时空异质性影响,增强不同时空范围的可比性;
2)本申请提出了累积植被指数变化率特征参量的提取方法,提出了合理的累积时段范围并明确了适宜的产量评估时间区间(比如本发明在示例性实施例中发现,合理的累积时段范围是自
到按照数据时间步长依次划分的包含
至
在内的任一时间节点;适宜的产量评估时间区间是生长前期中值对应时间节点+8d至生长前期中值对应时间节点+72d的区间),从而提高了累积植被指数变化率特征参量用于作物估产和长势监测评估的有效性;
3)本申请中提出了基于预测产量的长势量化评估方法,基于该模式的长势监测评估更加客观合理,且易于实现。
具体实施方式
为使本申请实施的目的、技术方案和优点更加清楚,下面将结合本申请实施例中的附图,对本申请实施例中的技术方案进行更加详细的描述。
1.研究方法
1.1遥感数据处理
1、基于GIS软件或相关平台,对MODIS遥感数据进行预处理,包括拼接、剪裁、质量控制等;
2、利用中高分辨率遥感数据,结合物候信息、解译标志数据等,进行监测区域作物识别,得到作物分布数据,本发明中以大豆作物为例;
3、利用作物分布数据对MODIS时序数据进行掩膜处理,得到目标作物区域时序数据集;
4、生产归一化植被指数(NDVI)产品,进行Savitzky-Golay滤波处理,构建NDVI 时序数据集;
NDVI是由红光和近红外两个波段得到的,由于该两个波段是植物光谱中反映光合作用和呼吸作用的最重要波段,该指数能够反映植被生长差异特征,且指数计算简单,易于实现,被广泛用于作物长势监测中。其形式如下:
式中,Rnir为近红外波段反射率,Rred为红光波段反射率。
1.2累积植被指数变化率(Accumulated Changing Rate of VI,ACRVI)
以变化指数表征的植被生长过程通常包括正增长阶段和负增长阶段(即下降阶段),利用前后植被指数差值可以表征其增长趋势,为突出表征某一时点植被变化程度和趋势,本研究提出植被指数变化率,形式见公式(2);为表征植被生长过程变化程度,提出累积植被指数变化率,形式见公式(3):
式中,CRVIi为i时段的植被指数变化率,NDVIi和NDVIi-1分别为i时段和i-1时段的植被指数,i表示时序产品对应的时间节点,i-1为i紧邻的前一个时间节点;
ACRVIi为截止i时间节点的累积植被指数变化率,t0为初始时间,本研究以播种期最低值出现时间为t0,T为特定时间节点。
1.3累积植被指数变化率特征值提取
结合作物生长过程曲线,提取5个植被指数特征值,包括播种期最低谷值(记为NDVImin_s)、生长期最高峰值(记为NDVImax_g)、生长前期最低-最高的中值(简称生长前期中值,记为NDVImid_gf)、生长后期最高-最低的中值(简称生长后期中值,记为NDVImid_gb) 和成熟期最低谷值(记为NDVImin_m)。
首先,根据研究区大豆种植习惯,确定监测起止时段(以月记);
其次,采用二次差分法对播种期、成熟期最低谷值和生长期最高峰值时点进行位置识别,对应时段分别记为
再提取相应NDVI值,基于已经提取的最低谷值和最高峰值,分别用算术平均法得到生长前期中值(其对应时段记为
)和生长后期中值(其对应时段记为
);
最后,以
为累计植被指数变化率统计起点,以
为统计终点,依次计算
到上述统计区间内按数据时间步长划分的不同时间节点的累积植被指数变化率ACRVI,从而得到不同年份截止不同时间点的累积植被指数变化率序列,数据时间步长为8d。
需要指出的是,通常作物产量预测和长势评估是在作物生长季中进行的,即可能没有完整的作物生长季相关遥感数据,因此,本研究参考Hoolst R V等[1]的方法得到多年不同时点的有效均值特征影像并以此作为未结束年份的准同期数据,随着作物生长期发展,进行数据的迭代更新。
1.4其他特征参量提取
为客观评估基于累积植被指数变化率的估产和长势评估效果,本研究同步计算得到较为常用的两个特征参量,即累积NDVI和累积植被指数变化速率。计算公式如下:
式中,ANDVIi为截止i时间节点的累积植被指数,t0为初始时间,本研究以播种期最低值出现时间为t0;SRVIi为i时段时间节点的植被指数变化速率,t为i时段时间节点和i-1时段时间节点的时间间隔;ASRVIi为截止i时间节点的累积植被指数变化速率。
1.5精度分析方法
为了评价模型的精度,采用常用的模型检验指标对其进行分析评价。通过计算模型模拟值和实测值之间的相对误差、决定系数对它们之间的符合度进行分析。计算公式如下:
式中,Re为相对误差,yi和y分别为预测值和实际值;R2为决定系数,n为样本数。R2取值范围为(0,1),其值越接近1,表示模拟值和实测值的符合度越高,相反,越接近 0表明符合度越低。
示例性应用:
1)研究区域
本研究区为位于巴西巴拉那州西部的托莱多市,其经度范围是-54°3′~-53°32′,纬度范围是-24°27′~-24°56′,海拔在254~690m之间,平均海拔为503m,全县面积为1196.2km2(图3)。托莱多气温通常在9.4℃~30℃之间。肥沃的土壤和相对平坦的地势使该市成为农业发达区。它是该州重要粮食生产基地之一,其农业GDP在巴拉那州和南部地区排名首位,农业增加值在全国位列第三[2]。主要农作物类型为大豆、玉米、烟叶、咖啡等。其中,大豆种植生长季为9月~次年3月。
2)数据来源与处理
2.1)统计调查数据
托莱多市大豆年度单产统计数据和苗情进度与长势调查数据来源于巴西巴拉那州农村经济部(DERAL)。其中,使用的单产统计数据年限为2010/11~2020/21,苗情进度与长势调查数据年限为2019/20~2020/21。
2.2)遥感及基础地理数据
本研究使用中分辨率成像光谱仪(Moderate Resolution ImagingSpectroradiometer, MODIS)的MOD09A1(地表反射率)8d合成产品,空间分辨率为500m,数据时间为2010 年9月~2021年3月。对其进行拼接、剪裁、质量控制等预处理。由于使用农作物掩膜之后的预报结果准确率最高[3],因此,利用Landsat8数据基于机器学习方法对研究区进行大豆种植区域提取,并以此为掩膜数据,完成掩膜处理。再进行归一化植被指数(NDVI) 提取,经过Savitzky-Golay[4]滤波处理后,生成大豆产区NDVI时序产品。
本研究中NDVI时序产品的时间表达以8d合成的截止日期为准。
2.3)提取植被指数特征值及其对应的特定时间节点
结合2010年9月~2021年3月时序遥感数据提取作物生长过程特征曲线,对各年的播种期、成熟期最低谷值和生长期最高峰值时点进行位置识别,对应时段分别记为
并提取相应NDVI
min_s、NDVI
min_m、NDVI
max_gNDVI值;并基于已经提取的最低谷值和最高峰值,分别用算术平均法得到生长前期中值NDVI
mid_gf,其对应时段记为
和生长后期中值NDVI
mid_gb,其对应时段记为
2.4)计算不同时间节点的累积植被指数变化率
累积植被指数变化率通过如下公式获得:
式中,CRVIi为i时段的植被指数变化率,NDVIi和NDVIi-1分别为i时段和i-1时段的植被指数,i表示时序产品对应的时间节点,i-1为i紧邻的前一个时间节点;ACRVIi为截止i时间节点的累积植被指数变化率,t0为播种期最低值出现时间,T为特定时间节点。
步骤106、将所述累积植被指数变化率时序产品集和对应的作物单产统计数据集进行回归分析,得出不同特定时间点上针对研究区域作物单产预测模型。
3)研究结果与分析
3.1)累积植被指数特征参量与产量关系分析
对2010/11-2018/19年的
的累积植被指数变化率与大豆单产进行回归分析,如图4(a)所示,累积植被指数变化率与大豆单产呈现正向线性关系。
利用Pearson相关分析计算相关系数为0.88,决定系数为0.77,并通过显著性水平为0.01的极显著性检验;利用线性回归模型,对2019/20和2020/21两个年份单产进行预测,得到相对误差分别为-2.51%和2.14%,平均相对误差为-0.18%。
累积NDVI与大豆单产的相关系数为0.79,决定系数为0.62,通过显著性水平为0.05 的显著性检验。利用线性回归模型,得到2019/20和2020/21两个年份预测单产结果,其相对误差分别为13.97%和27.20%,平均相对误差为20.59%。
累积植被指数变化速率与大豆单产的相关系数为0.89,决定系数为0.8,通过显著性水平为0.01的极显著性检验。利用线性回归模型,得到2019/20和2020/21两个年份预测单产结果,其相对误差分别为10.74%和0.78%,平均相对误差为-5.76%。
经比较认为,在该时段,基于累积植被指数变化率进行产量预测效果最佳,累积植被指数变化速率预测效果次之,累积NDVI预测效果不佳。
对2010/11-2018/19年的
的累积植被指数变化率与大豆单产进行回归分析,如图4(b)所示,累积植被指数变化率与大豆单产呈现负向线性关系。
利用Pearson相关分析计算相关系数为0.54,决定系数为0.288,未通过显著性水平为0.05的检验。利用线性回归模型,对2019/20和2020/21两个年份单产进行预测,得到相对误差分别为-21.48%和-4.48%,平均相对误差为-12.98%。
累积NDVI与大豆单产的相关系数为0.79,决定系数为0.63,未通过显著性水平为0.05的显著性检验。利用线性回归模型,得到2019/20和2020/21两个年份预测单产结果,其相对误差分别为-8.82%和4.65%,平均相对误差为-2.09%。
累积植被指数变化速率与大豆单产的相关系数为0.77,决定系数为0.59,通过显著性水平为0.05的显著性检验。利用线性回归模型,得到2019/20和2020/21两个年份预测单产结果,其相对误差分别为-18.62%和-0.91%,平均相对误差为-9.77%。
该阶段,仅累积植被指数变化速率通过显著性检验,但平均相对误差偏高,虽然累积NDVI平均相对误差较低,好于另外两个特征参量拟合效果,但未通过显著性检验,因此认为,在该阶段三种特征参量均不适宜进行产量预测。
3.12)多时段累积植被变化特征参量与单产关系分析
以生长前期中值对应时间节点
为起始时间,分别以其顺次加8d,直至加 80d为计算累积植被变化特征参量的终止时间,共计11个监测时间模式,对 2010/11-2018/19共计9年的大豆单产与11个监测时间模式的三个累积植被指数变化特征参量进行相关性分析,结果如图5-图7所示(其中,**为极显著相关,*为显著相关)。
累积植被指数变化率和累积植被指数变化速率与单产相关系数时序尺度上分布较为相似,呈抛物线形状,在前期和后期均较低,相关关系表现为不显著,随着生长阶段发展,显著性水平表现为显著或极显著,相关关系在生长前期中值对应时间节点+40d时达到峰值,分别为0.96和0.97。累积NDVI除生长前期中值对应时间节点+24d时出现低值外,总体呈现线性递增趋势,即随生长阶段发展,相关系数逐渐上升,在生长后期达到高值,为0.89,且表现为极显著性相关。
生长前期中值主要为大豆苗期-分枝期,刚进入该阶段时,累积植被指数变化率、累积NDVI及累积植被指数变化速率三个特征参量与单产相关系数均较低,且相关关系未达到显著水平。随着大豆生长,累积植被指数变化率和累积植被指数变化速率与单产呈现显著或极显著相关关系,前者相关系数总体高于后者,可认为累积植被指数变化率在生长前期中值对应时间节点+8d至生长前期中值对应时间节点+72d期间,即TNDVImid_gf+8d~TNDVImid_gf+72d期间,对大豆单产均有比较好的反映,可以用于大豆估产,是三个特征参量中的优选指标。另外,累积NDVI与大豆单产相关系数在生长前期中值对应时间节点 +72d及以后表现好于另外两个特征参量,可以作为全生育期产量估测参量。
进一步,对2010/11-2018/19年的
的累积植被指数特征参量与大豆单产进行回归分析,建立不同时期的线性拟合模型,并基于模型对2019/20和2020/21两个年份单产进行预测,得到基于三个累积植被指数特征参量的单产预测相对误差结果,如图8所示:
基于累积NDVI的预测结果总体偏高,且相对误差较高,其中,2019/20预测相对误差在-13.53%~36.63%之间,2020/21预测相对误差在-0.07%~43.71%之间。
基于累积植被指数变化速率的预测结果总体偏低,在大豆生长中后期相对误差两年份差异较大,其中,2019/20预测相对误差在-17.76%~4.33%之间,2020/21预测相对误差在 -7.35%~5.3%之间。
基于累积植被指数变化率的预测结果在前期偏高、后期偏低,其中,2019/20预测相对误差在-13.52%~14.65%之间,2020/21预测相对误差在-7.43%~6.98%之间。
对2019/20和2020/21两年相对误差进行平均,得到三个特征参量的相对误差均值时序曲线,其中,基于累积NDVI的相对误差均值在-6.8%~40.2%之间,基于累积植被指数变化率的相对误差均值在-6.8%~6.6%之间,基于累积植被指数变化速率的相对误差均值在-11.0%~-0.6%之间。
进一步,对2019/20和2020/21两个年份基于三个累积植被指数特征参量的单产预测相对误差进行分时段平均,一个时段为
一个时段为
结果如表1所示。
表1多时段累积植被指数特征参量不同时段平均相对误差对比
可以看出:基于累积NDVI的预测相对误差总体最高,2019/20年度多时段平均相对误差为23.28%,2020/21年度平均相对误差为32.74%。
平均相对误差最小的为累积植被指数变化率,其中,2019/20年度多时段平均相对误差为0.48%,2020/21年度平均相对误差为0.71%。基于累积植被指数变化速率的预测效果次之,2019/20的为-6.96%,2020/21的为-1.6%。
由于
多为大豆的开花结荚-鼓粒期,是大豆产量形成关键阶段,因此对该阶段预测平均相对误差进行比较分析。
基于累积NDVI的预测相对误差总体最高,2019/20年度的平均相对误差为33.22%, 2020/21年度平均相对误差为40.65%。平均相对误差总体最小的为累积植被指数变化率,其中,2019/20年度的平均相对误差为0.11%,2020/21年度的平均相对误差为3.69%。基于累积植被指数变化速率的预测效果次之,2019/20的为-6.42%,2020/21的为1.17%。
3.2)基于累积植被指数变化特征参量的典型年长势分析
基于对累积植被指数变化特征参量与单产关系分析结果,认为在特定预测精度范围内,由于该累积植被指数变化特征参量可以反映作物未来产量水平和趋势,因此可以利用相关参量进行不同时段的长势监测评估。
选取2010/11~2020/21间典型年,包括2016/17(最高产年)、2019/20(上一年)、2020/21(当年)三个年份的最低谷值-最高峰值时段进行应用分析。
各特定时间节点长势等级需基于单产等级和单产预测模型进行推算,其中,单产等级以多年单产均值μ和标准差σ确定,本研究中,产量在[μ-0.5σ,μ+0.5σ]区间定级为中等水平(其对应的长势等级为中等),大于(μ+0.5σ)定级为好(其对应的长势等级为好),小于(μ-0.5σ)定级为差(其对应的长势等级为差)。
基于该标准得到基于累积植被指数变化特征参量的三个年份最高峰值期的长势评价结果,如图9所示。
利用DERAL2019/20和2020/21两年同期调查统计数据进行对比分析和相对误差分析,如下表2所示。
表2长势等级与精度统计
根据统计数据,2019/20年为3846kg/ha,2020/21年大豆单产为3350kg/ha。按年份,基于三个累积植被指数变化特征参量的长势监测图(图9)总体上都能体现现出2020/21年长势不及前两年份,根据表2看出,该年度苗情为好的比例低于2019/20。这是由于该年度大豆种植和出苗期出现了较重干旱,导致部分地区播种偏晚、苗情偏弱,影响了生物量积累,长势表现偏差。
为评估基于三个累积植被指数变化特征参量的长势评价效果,结合统计调查数据,进行了相对误差测算。
由表2可知,基于累积植被指数变化率的2019/20年度的相对误差分别为0.5%、1.3%和-4.4%;2020/21年度的相对误差分别为-1.9%、-0.2%和23.1%。两年份中,好和中等等级的平均相对误差分别为-0.7%和0.5%,差等级偏低,为9.4%,总体上,相对误差低于10%,尤其是好和中等的相对误差非常低,在-1.9%~1.3%之间。
基于累积NDVI的2019/20年度的相对误差分别为7.6%、4.3%和-53.4%;2020/21年度的相对误差分别为20.9%、-47.8%和-52.3%。两年份中,不同等级的相对误差差异较大,尤其是中等和差的等级,这是由于基于累积NDVI的估产在常年估计时有偏高、高产年估计时有偏低的趋势,因此导致在长势评估中也出现了较大的偏差。
基于累积植被指数变化速率的2019/20年度的相对误差分别为-3.4%、29.1%和3.0%; 2020/21年度的相对误差分别为-3.6%、9.8%和2.6%。两年份中,好和差等级的平均相对误差分别为-3.5%和2.8%,中等等级偏低,为19.5%。其评估效果介于前两特征参量之间。
综上,总体看,累积植被指数变化率是三个特征参量中最具可行性和有效性的长势评估参量。
以上所述,仅为本申请的具体实施方式,但本申请的保护范围并不局限于此,任何熟悉本技术领域的技术人员在本申请揭露的技术范围内,可轻易想到的变化或替换,都应涵盖在本申请的保护范围之内。因此,本申请的保护范围应以所述权利要求的保护范围为准。
参考文献:
[1]Hoolst R V,Eerens H,Haesen D,et al.FAO's AVHRR-based AgriculturalStress Index System(ASIS)for global drought monitoring[J].InternationalJournal of Remote Sensing,2016,37(1-2):418-439.
[2]Manfrin,J.,Ceraso,A.F.F.,&
Jr,A.C.Evaluation of OdoursPerception by Residents of the Municipality of Toledo/ParanáState/Brazil[J].Journal of Scientific Research and Reports,2018,21(2),1-9.https://doi.org/10.9734/JSRR/2018/44890.
[3]李新星.基于时间序列植被指数的水稻发育期提取和估产[D].浙江大学,2017.
[4]边金虎,李爱农,宋孟强,et al.MODIS植被指数时间序列Savitzky-Golay滤波算法重构[J]. 遥感学报,2010,14(04):725-741.