一种农作物季节早期多源遥感产量数据处理方法及预测方法

摘要

本发明公开了一种农作物季节早期多源遥感产量数据处理方法及预测方法，其中数据处理方法包括以下步骤：对多光谱图像、激光雷达数据进行预处理，得到对应目标区域的多个多光谱特征图、激光雷达特征图；将所有的多光谱特征图、激光雷达特征图合并成一个原数据矩阵，对原数据矩阵进行双变量参数的非负矩阵分解，得到新特征权重矩阵；从新特征权重矩阵提取构成每个新特征所占比重最大的原始特征序号，从原数据矩阵中筛选得到新数据矩阵，由此完成数据处理。

说明书

技术领域

本发明涉及农业的遥感图像处理技术领域，具体为一种农作物季节早期多源遥感产量数据处理方法及预测方法。

背景技术

快速、精确的作物产量预测，对国家制定粮食政策，进行价格宏观调控，实现乡村经济的振兴以及粮食进出口贸易等都具有重要意义。传统的作物产量预测方法主要有农学估产法、统计估产法以及气象统计法等，而这些方法往往只适用于小面积估产。自20世纪70年代以来，遥感技术因其宏观、动态和快速的特点在农业生产上得到了广泛的应用，为农作物长势的监测以及籽粒产量的预测提供了有效的途径。

近年来，一方面获取光谱及其他结构信息的传感器逐渐小型化并不断提升可集成性，另一方面搭载这些传感器的无人机平台的载荷与稳定性亦有了长足进步，这些变化使得多种传感器联合、同步数据获取越来越成为主流的遥感探测方案之一。将多模式传感数据结合应用于产量预测工作，可以结合不同传感器系统所包含的丰富的光谱、空间和结构信息的优点，获得更全面、更详细的当前作物信息，有助于改善产量预测的预测精度。

现有的作物遥感估产技术大体可分为3个步骤：通过分析遥感影像数据估算作物种植面积；通过分析遥感影像数据提取作物相应的植被指数，以进行作物长势状况的监测；构建植被指数和作物产量以及其它气象、农学参数等资料的单产估算模型，再经过进一步的计算得到总产。而在利用多源遥感数据特征和作物产量构建单产估算模型时，为避免冗余特征或无效特征降低估产模型的预测性能，往往会采用PCA算法对数据进行降维处理，然后针对降维后的遥感数据与产量数据采用多元线性回归的方法构建关系函数，得到简单的产量预测模型，但是所获得的模型可解释性差且预测精度不足。

发明内容

本发明为了解决以上现有技术存在不足与缺陷的问题，提供了一种农作物季节早期多源遥感产量数据处理方法及预测方法。

为实现上述本发明目的，采用的技术方案如下：

一种农作物季节早期多源遥感产量数据处理方法，包括以下步骤：

对多光谱图像、激光雷达数据进行预处理，得到对应目标区域的多个多光谱特征图、激光雷达特征图；

将所有的多光谱特征图、激光雷达特征图合并成一个原数据矩阵，对原数据矩阵进行双变量参数的非负矩阵分解，得到新特征权重矩阵；

从新特征权重矩阵提取构成每个新特征所占比重最大的原始特征序号，从原数据矩阵中筛选得到新数据矩阵，由此完成数据处理。

优选地，对多光谱图像进行预处理，得到多光谱特征图，具体步骤如下：

对多光谱图像的多个通道图，利用能反映农作物生长态势的N1个不同的植被指数进行处理，得到N1个植被指数图；

以最终产量收成时生物量采样区的大小为标准，对N1个植被指数图同步进行网格划分；

然后在每幅植被指数图的每个网格单位内提取N2个不同的衍生物形成N1×N2个个多光谱特征图，其中N1、N2均为正整数。

进一步地，所述的植被指数共N1＝10项，分别为归一化植被指数、归一化差异红色边缘指数、绿光归一化差值植被指数、增强型植被指数、改性花青素含量指数、优化土壤调节植被指数、简化冠层叶绿素含量指数、转化叶绿素吸收和反射率指数、三角绿色指数、可见大气阻力指数。

进一步地，所述的衍生物共N2＝7项，分别为包括最大值、最小值、平均值、标准值、低四分位值、中位值、高四分位值。

优选地，对激光雷达图像进行预处理得到激光雷达特征图，具体步骤如下：

以最终产量收成时生物量采样区的大小为标准，在各个区域单位内提取所需植被数据，生成多个激光雷达特征图。

进一步地，所述的植被数据包括最大高度、平均高度、平均平方高度、高度标准差、高度偏态、高度峰度、第5至95位高度百分位、第5至95位从地面到高度百分位数的分数、4个高度点数的占比。

优选地，将多光谱特征图和激光雷达特征图按同样的位置顺序各自设置成行向量X

然后将多光谱特征图和激光雷达特征图共计m个特征图转化的行向量进行合并，得到集合多光谱特征图和激光雷达特征图的原数据矩阵X

进一步地，对原数据矩阵进行双变量参数的非负矩阵分解，得到新特征权重矩阵，具体如下：

根据原数据矩阵X

为了使原数据矩阵X

采用交替最小二乘方法，先随机初始化参数p，利用p和式(1)计算得到参数q，再利用参数q和式(1)求解参数p，不断重复这个过程直至收敛于局部最小值；

通过不断交替优化P

再进一步地，从原数据矩阵中筛选得到新数据矩阵，具体如下：

初始化一个空的新数据矩阵Y

对所述的新特征权重矩阵Q

在原数据矩阵X

在经过k次遍历后，得到经过关键特征提取后的新数据矩阵Y

一种农作物季节早期多源遥感产量数据预测方法，包括以下步骤：

通过所述的农作物季节早期多源遥感产量数据处理方法获得新数据矩阵Y

通过采集对应的最终产量收成时生物量采样区的产量数据Z

利用最终的产量预测模型对农作物产量进行预测。

本发明的有益效果如下：

本发明针对能提供丰富光谱信息和作物三维结构信息的多光谱及激光雷达两种不同模式的遥感数据，进行预处理并在同一网格划分单位内进行衍生物提取，挖掘更多对产量预测有潜在效益的预测因子。然后通过对数据矩阵进行双变量参数的非负矩阵分解来筛选出比较关键的预测因子，在消除原特征之间的冗余性的同时使降维得到的新特征具备实际意义。最后采用神经网络对预测因子和产量结果进行非线性回归拟合，得到的产量预测模型，可以在线性回归模型上进一步提高产量的预测精度。

附图说明

图1是本发明所述的数据处理方法的流程图。

图2是特征图转化为一个数据矩阵的示意图。

图3是新特征权重矩阵得到新数据矩阵的示意图。

图4是本发明的所述预测方法的示意图。

图5是利用本发明得到的甘蔗的产量预测模型效果图。

图6是传统的甘蔗的产量预测模型效果图。

具体实施方式

下面结合附图和具体实施方式对本发明做详细描述。

实施例1

如图1所示，一种农作物季节早期多源遥感产量数据处理方法，包括以下步骤：

对多光谱图像、激光雷达数据进行预处理，得到对应目标区域的多个多光谱特征图、激光雷达特征图；

将所有的多光谱特征图、激光雷达特征图合并成一个原数据矩阵，对原数据矩阵进行双变量参数的非负矩阵分解，得到新特征权重矩阵；

从新特征权重矩阵提取构成每个新特征所占比重最大的原始特征序号，从原数据矩阵中筛选得到新数据矩阵，由此完成数据处理。

在一个具体的实施例中，对多光谱图像进行预处理，得到多光谱特征图，具体步骤如下：

对多光谱图像的多个通道图，利用能反映农作物生长态势的N1个不同的植被指数进行处理，得到N1个植被指数图(VIs)；

以甘蔗遥感数据集为例，所述的植被指数共N1＝10项，分别为归一化植被指数(NDVI)、归一化差异红色边缘指数(NDRE)、绿光归一化差值植被指数(GNDVI)、增强型植被指数(EVI)、改性花青素含量指数(MACI)、优化土壤调节植被指数(OSAVI)、简化冠层叶绿素含量指数(SCCCI)、转化叶绿素吸收和反射率指数(TCARI)、三角绿色指数(TGI)、可见大气阻力指数(VARI)。

所述的植被指数的计算公式具体如下：

NDVI＝(NIR-R)/(NIR+R)

NDRE＝(NIR-RE)/(NIR+RE)

GNDVI＝(NIR-G)/(NIR+G)

EVI＝2.5(NIR-R)/(NIR+6R-7.5B+1)

MACI＝NIR/G

OSAVI＝(1+0.16)(NIR-R)/(NIR+R+0.16)

SCCCI＝NDRE/NDVI

TCARI＝3[RE-R-0.2(RE/G)(RE/R)]/OSAVI

TGI＝-0.5[(668-475)(R-G)-(668-560)(R-B)]

VARI＝(G*R)/(G+R-B)

式中，RE表示传感器的红边通道值、NIR表示传感器的近红外通道值、R表示传感器的标准红色通道值、G表示传感器的标准绿色通道值、B表示传感器的标准蓝色通道值；

以最终产量收成时生物量采样区的大小为标准，对10个植被指数图同步进行网格划分；比如多光谱图像分辨率为0.5，指定区域为2m*2m，则每个网格单位包含16个像素点。

然后在每幅植被指数图的每个网格单位内提取N2＝7个不同的衍生物形成N1×N2个个多光谱特征图，其中N1、N2均为正整数。比如取每个单位中NDVI的最大值形成max NDVI特征图。所述的衍生物共7项，分别为最大值max，最小值min，平均值avg，标准值std，低四分位值p25，中位值p50，高四分位值p75。以衍生物为多光谱数据特征使得每幅VIs图各自生成7个多光谱特征图，共计得到70个多光谱特征图。

在一个具体的实施例中，对激光雷达图像进行预处理得到激光雷达特征图，具体步骤如下：

以最终产量收成时生物量采样区的大小为标准，使用LasTools等软件在各个区域单位内提取所需植被数据，生成多个激光雷达特征图。

所述的植被数据包括最大高度max，平均高度avg，平均平方高度qav，高度标准差std，高度偏态ske，高度峰度kur，第5至95位高度百分位p05～p95，第5至95位从地面到高度百分位数的分数b05～b95，四个高度点数的占比d00～d03，共计48个预测变量生成48个激光雷达特征图。

在一个具体的实施例中，将多光谱特征图和激光雷达特征图按同样的位置顺序各自设置成行向量X

然后将多光谱特征图和激光雷达特征图共计m个特征图转化的行向量进行合并，得到集合多光谱特征图和激光雷达特征图的原数据矩阵X

如图2所示，以甘蔗预测实验为例，上述提取有70个多光谱特征图和48个激光雷达特征图，也就是一共有118个特征图，而随机的生物量采样点(2m*2m为单位且有记录最终产量数据的重点区块)有56个，所以将这56个区块的118个特征放在一起并转化成行向量会得到118个格式为[1，56]的行向量，最后合并成一个原数据矩阵X[118，56]，行数118表示特征个数，列数56表示样本数。

在一个具体的实施例中，对原数据矩阵进行双变量参数的非负矩阵分解，得到新特征权重矩阵，具体如下：

根据原数据矩阵X

为了使原数据矩阵X

采用交替最小二乘方法，先随机初始化参数p，利用p和式(1)计算得到参数q，再利用参数q和式(1)求解参数p，不断重复这个过程直至收敛于局部最小值。

通过不断交替优化P

在一个具体的实施例中，从原数据矩阵中筛选得到新数据矩阵，具体如下：

初始化一个空的新数据矩阵Y

对所述的新特征权重矩阵Q

在原数据矩阵X

在经过k＝15次遍历后，得到经过关键特征提取后的新数据矩阵Y

本实施例首先对所有的甘蔗种植区域进行2m*2m的网格划分，然后在甘蔗开始生长早期利用无人机搭载传感器采集可以精确定位网格位置的多光谱图像和激光雷达数据，并且在最后收成时记录下各个产量采样区域的产量值(就是记录下黑色区域的甘蔗有多少kg)

所述的多光谱图像和激光雷达数据以像素点的形式保存(其分辨率要小于2m*2m)。

实施例2

如图4所示，一种农作物季节早期多源遥感产量预测方法，包括以下步骤：

对多光谱图像、激光雷达数据进行预处理，得到对应目标区域的多个多光谱特征图、激光雷达特征图；

将所有的多光谱特征图、激光雷达特征图合并成一个原数据矩阵，对原数据矩阵进行双变量参数的非负矩阵分解，得到新特征权重矩阵；

从新特征权重矩阵提取构成每个新特征所占比重最大的原始特征序号，从原数据矩阵中筛选得到新数据矩阵Y

通过采集对应的最终产量收成时生物量采样区的产量数据Z

利用最终的产量预测模型对农作物产量进行预测。

将新数据矩阵Y[15，56]和步骤一中同一网格位置的产量数据Z[1,56]送入到深度神经网络中进行非线性回归拟合训练，构建遥感数据特征与最终产量的潜在关系函数，可以获得预测精度较高的产量预测模型，同时保留模型输入特征参数的实际意义，提高产量预测方法的可解释性。

在本实施例中，要把数据预处理成网格单位的特征图，而不是原来的像素点为单位的特征图，是为了与产量数据保持在一个尺度上，以网格为单位的特征值可以一一对应同一网格位置上的产量值，以此来构建输入与输出的关系。

在本实施例中，所用数据是农田遥感数据，如图5所示，其是通过本发明对于甘蔗得到的产量预测模型性能。如图6所示，其传统方法得到的产量预测模型性能的拟合效果比较示意图。由此，根据对比，可以看出本发明所得到的产量预测模型的效果更优。

实施例3

一种计算机设备，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，所述的处理器执行所述的计算机程序时，实现如下步骤：

对多光谱图像、激光雷达数据进行预处理，得到对应目标区域的多个多光谱特征图、激光雷达特征图；

将所有的多光谱特征图、激光雷达特征图合并成一个原数据矩阵，对原数据矩阵进行双变量参数的非负矩阵分解，得到新特征权重矩阵；

从新特征权重矩阵提取构成每个新特征所占比重最大的原始特征序号，从原数据矩阵中筛选得到新数据矩阵，由此完成数据处理。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。