基于MapReduce和神经网络的农作物产量遥感估算方法

技术领域

本发明涉及遥感处理领域，尤其涉及一种基于MapReduce和神经网 络的农作物产量遥感估算方法。

背景技术

传统的农作物产量估算存在调查范围小、人力物力耗费巨大的问题。 遥感技术的发展，为农作物的产量估算提供了有力的工具。

公开号为102162850A的专利文献公开了一种作物遥感估产的方法。 该方法，基于遥感信息获取的瞬时性与广域性，结合小麦产量形成过程及 其与气候环境的关系，建立了较为简化的小麦产量预测模型；通过组件化 的设计方法实现了遥感信息和估产模型的耦合，即，利用抽穗期遥感影像 反演的LAI和生物量及时替换小麦估产模型对应参数变量，进而实现对 单点小麦产量的估测，估产精度能达到90％以上；进一步，采用“点”(样 点产量)与“面”(遥感区域)尺度转换的方法，进行区域小麦产量遥感分 级监测预报，制作区域小麦产量遥感监测分级预报专题图，具有直观、具 体、时效性好的特点，对县级农业技术人员获取区域小麦布局信息或指导 生产管理具有较好的实用性。

在这类方法中，需要从遥感图像中提取NDVI值。

NDVI是Normalized Difference Vegetation Index的简写，中文名为归 一化植被指数，又被称为标准化植被指数，它对植被空间分布密度以及 植物生长状态有极佳的指示作用。NDVI与植被分布密度有线性相关性。 根据实验结果显示，归一化植被指数对土壤背景的变化比较敏感；其 次，归一化植被指数是单位像元内的植被覆盖形态、植被类型、植被生 长状况等的综合反映，它的数值大小由植被覆盖度和叶面积指数这两个 要素所决定；NDVI在植被覆盖度的检测领域应用较广，主要原因是它对 植被覆盖度的检测幅度较宽，有比较好的空间和时间适应性。归一化植 被指数NDVI在植被指数中有着非常重要的位置，它较其他植被指数主要 有以下几个方面的突出优点：1.植被覆盖度的检测范围较大；2.植被检测 灵敏度较高；3.能削弱太阳高度角以及大气所带来的噪音；4.能消除地 形和群落结构的阴影以及辐射干扰。

NDVI计算通常需要将红色可见光通道（波长范围为0.6-0.7μm）和 近红外光谱通道（波长范围为0.7-1.1μm）进行组合，用来设计NDVI， 具体的计算公式如下：

NDVI＝（Rn－Rr）／（Rn＋Rr）

上式中，Rn表示的是近红外波段的反射率，Rr表示的是红光波段的反 射率。

而现有技术利用遥感图像进行农作物产量估算时，存在着遥感数据量 庞大，处理速度慢，NDVI提取效率低等问题。

MapReduce是一个编程框架，它为程序员提供了一种快速开发海量数 据处理程序的编程环境，并能够让基于这种机制开发出的处理程序以稳 定、容错的方式并行运行于由大量商用硬件组成的集群上。同时， MapReduce又是一个运行框架，它需要为基于MapReduce机制开发出的 程序提供一个运行环境，并透明管理运行中的各个细节。每一个需要由 MapReduce运行框架运行的MapReduce程序也称为一个MapReduce作业 (mapreduce job)，它需要由客户端提交，由集群中的某专门节点负责接收 此作业，并根据集群配置及待处理的作业属性等为其提供合适的运行环 境。其运行过程分为两个阶段：map阶段和reduce阶段，每个阶段都根据 作业本身的属性、集群中的资源可用性及用户的配置等启动一定数量的任 务(也即进程)负责具体的数据处理操作。

如何利用MapReduce来提高遥感数据的处理效率，从而提高农作物 的遥感产量估算效率，是亟需解决的问题。

发明内容

为了解决遥感图像巨大的数据量问题、切图的效率问题、NDVI值提 取的效率问题、农作物估产的准确性和稳定性问题等实际问题，本发明结 合MapReduce程序，将神经网络进行优化，提出一种农作物产量遥感估 算的方法。

一种基于MapReduce和神经网络的农作物产量遥感估算方法，包括 如下步骤：

步骤1，对输入的遥感图像进行多线程并发切图，得到若干瓦片，各 个瓦片以其顶点的经纬度数据命名；

步骤2，根据瓦片名称所提取的经纬度数据以及遥感图像中各个地区 的边界经纬度数据，对所有瓦片进行MapReduce操作，得到遥感图像中 各个地区的NDVI值；

步骤3，对于各个地区，将其NDVI值输入至经过训练的神经网络中， 得到该地区的农作物产量估算值。

在步骤1中，各个瓦片以其顶点的经纬度数据命名，是指每个瓦片以 其顶点的经度和纬度组合命名，且顶点为瓦片四个顶点的其中一者，例如 以瓦片左下角顶点的经纬度数据命名，方便步骤2中的经纬度数据提取。 其中瓦片的尺寸由用户预设，例如为512像素*512像素。该方法并行切图 和MapReduce的运用，大大优化了光谱遥感图像的处理效率，使本发明 具有高效的特点。

在步骤1中，多线程并发切图的步骤如下：

步骤1-1，由一个Dispatch线程计算切割任务，并判断是否还有切割 任务：是，则将所得切割任务插入到Task队列；否则，发送消息通知Task 线程已无切割任务；

步骤1-2，由若干个Task线程从Task队列依次获取切割任务进行切 割，每个Task线程在完成当前切割任务后判断Task队列中是否还有切割 任务：是，则获取下一个切割任务；否则，判断是否接收到已无切割任务 的消息：是，结束切割；否则，等待Task队列中插入的切割任务。

单线程的图像切割算法，在每一层次的图像切割时只负责将图像切割 成每一片分辨率为512*512的小图（即瓦片），每一个瓦片的切割与前一 个瓦片或者后一个瓦片的切割是没有直接关系，也即，两个瓦片的切割不 存在逻辑上的关系，它们之间的关系仅仅是这两个瓦片在规整结果中可能 是相临的，两个瓦片的切割是完全是可以独立开。而且，现在的计算机硬 件配置中，CPU大多为多核，可以并发地执行多个任务。因而，传统的流 程切割图片，属于线性地切割每一个瓦片，每一个时刻只有一个核在工作， 在一定程度上浪费了一定的计算资源。既然任意两个瓦片的切割在逻辑上 没有直接关系，且现代计算机CPU大多为多核，可以通过多线程技术手 段，将瓦片的切割任务分配给多个线程进行，从而更进一步的提高算法性 能。

本发明的实现方案主要两种不同角色的线程组成，Dispatch线程与 Task线程。在系统中只有一个Dispatch线程，有若干个Task线程。Task 线程的具体数量根据当前系统CPU核数而定，默认为4。Dispatch线程负 责计算各个瓦片的左下角经纬度以及各个瓦片左下角像素点在遥感图像 中的坐标值。Task线程负责从Task队列中取若干个切割任务，切割任务 具体数量视系统实际情况而定，默认为4，进行图像切割，计算文件名， 并保存。当Dispatch线程将计算完成所有的任务描述信息并已插入到Task 队列，并且，Task队列中的所有任务描述都已经被Task线程计算完成， 此层次下的影像切割流程完成。

在步骤1-1中，Dispatch线程计算切割任务的方法为：计算出该切割 任务中的瓦片左下角顶点位置处的经纬度，以及左下角处像素点在遥感图 像中的坐标值。

这些信息描述了一个切割任务，所有计算任务组成了Task队列，并 由Dispatch线程维护。

在步骤1-2中，每个Task线程切割的方法为：根据步骤1-1所得的瓦 片左下角像素点在遥感图像中的坐标值，按照预设的瓦片尺寸从遥感图像 中进行相应切割。

瓦片左下角像素点在遥感图像中的坐标值包括横坐标值和纵坐标值， 将横坐标值作为瓦片的左边边界线的横坐标值，将纵坐标值作为瓦片的下 方边界线的纵坐标值，再加上预设的瓦片尺寸（例如512*512），从而得到 遥感图像中每个瓦片的切割区域。

步骤2的具体步骤如下：

步骤2-1，获取遥感图像中各个地区的边界经纬度数据，并按照地理 位置从西向东，从北向南的顺序排列边界经纬度数据；

步骤2-2，对各个瓦片进行Map操作，得到各个瓦片的NDVI值，并 输入对应的Reduce节点；

步骤2-3，对输入的NDVI值进行Reduce操作，得到遥感图像中各个 地区的NDVI值。

其中遥感图像中各个地区的边界经纬度数据为预先获取的数据，例如 可以由用户从历史数据中得到。其中将从西向东方向作为第一排序字段， 将从北向南方向作为第二排序字段。排序完成后的边界经纬度数据顺序 为，越靠西北方位的边界经纬数据次序越靠前，且同样纬度下，越靠西的 边界经纬数据次序越靠前。

在步骤2-2中，对各个瓦片进行Map操作的步骤如下：

步骤2-21，确定当前瓦片所在的地区，其中各个地区具有对应的地区 编号；

步骤2-22，得到当前瓦片的NDVI值；

步骤2-23，将所获得的NDVI值按对应的地区编号输入至相应的 Reduce节点。

其中Reduce节点与地区编号对应，具有相同地区编号的不同瓦片将 自身的NDVI值输入同一Reduce节点。

在步骤2-3中，在每个Reduce节点进行的Reduce操作的步骤如下：

步骤2-31，将该Reduce节点中输入的NDVI值相加，得到该节点的 NDVI值之和；

步骤2-32，得到该Reduce节点对应的地区编号与该节点的NDVI值 之和构成的数值对，从而得到各个地区的NDVI值。

由于地区编号与Reduce节点对应，因此通过将Reduce节点中的NDVI 值相加，所得到的NDVI值就是该地区编号对应的NDVI值。

在步骤2-2中，当前瓦片所在的地区确定方法为：

步骤a，根据瓦片四角的经纬度计算瓦片中心点C的经纬度（Lng，Lat）， 其中Lng表示经度，Lat表示纬度；

步骤b，从经过排序的边界经纬度数据中，自最西北角的地区开始， 将纬度值为Lat的边界点按经度值从大到小排序；

步骤c，用二分查找法找到最后一个经度值大于Lng的边界点，该边界 点所在的经度线就是当前瓦片所属的经度线；

步骤d，沿步骤c所确定的经度线向南查找，直到找到最后一个纬度 大于Lat的边界点，该边界点所在的地区就是当前瓦片所在的地区。

其中越靠西经度值越大，且越靠北纬度值越大。其中瓦片四个顶点的 经纬度通过瓦片左下角顶点的经纬度以及瓦片尺寸获得。

在步骤3中，神经网络为经过遗传算法优化的BP神经网络。

用神经网络进行产量估算可以省去对NDVI值进行平滑处理的步骤， 具有方便快捷的特点，而用遗传算法对传统的BP神经网络进行优化，则使 得神经网络更易得到全局最优解，且收敛速度更快。

在步骤3中，对每个地区的进行产量估算的方法，用包含NDVI-农作 物产量数值对的样本数据对神经网络进行训练，并将待估算的地区NDVI 数值输入经过训练的神经网络中，经过训练的神经网络输出即为该地区的 农作物产量估算值。

其中NDVI-农作物产量数值对是指遥感图像中所包含的各个地区 NDVI值及其对应农作物产量值的各个数值对。这些样本数据来自历史数 据，例如，过去一年所采集的包含各个地区的NDVI值-农作物产量数值对 的样本数据。

本发明的实现过程中，已经将光谱遥感图像巨大的数据量问题、切图 的效率问题、NDVI值提取的效率问题、农作物估产的准确性和稳定性问 题等实际问题考虑在内，为农作物产量的遥感估算提供了高效、可靠的解 决方案。

附图说明

图1为本发明一个实施例中并发切割图像的示意图；

图2为本发明当前实施例中并发切割图像的流程图；

图3为本发明方法当前实施例中神经网络训练相对误差随时间变化的 示意图；

图4为本发明实施例中利用GA-BP模型所得误差与利用简单的BP模型 所得误差的对比图；

图5为本发明方法的步骤流程图。

具体实施方式

现结合实施例和附图对本发明方法进行详细解释。

本发明采用小麦主产区河南省南阳市在2008年CentOS提供的 EOS/MODIS光谱遥感图像。如图5所示，本发明方法的步骤如下：

步骤1，对输入的遥感图像进行多线程并发切图，得到若干瓦片，各 个瓦片以其顶点的经纬度数据命名。

步骤1由1个Dispatch线程和多个Task线程共同完成。如图1所示，在本 发明实施例中，Task线程的个数为4。Task线程负责从Task队列中取4个具 体的Task进行图像切割，计算文件名，保存文件。

多线程并发切图的具体流程如图2所示：

由Dispatch线程分配切割任务，并判断是否还有切割任务：是，则将 所得切割任务插入到Task队列末端；否则，发送消息通知Task线程已无切 割任务。每个切割任务的计算方法为：计算出该切割任务对应的瓦片在左 下角位置处的经纬度，以及左下角处像素点在原遥感图像中的坐标值。

由Task线程进行切割，每个Task线程执行切割任务的步骤如下：从Task 队列的队头获取切割任务进行切割，完成当前切割任务后判断队列中是否 还有等待完成的切割任务：是，则获取下一个切割任务；否则，判断是否 接收到Dispatch线程所发送已无切割任务的消息：接收到消息则结束切割， 否则等待Dispatch线程插入切割任务。

完成全部切割任务后，得到尺寸为512像素*512像素的瓦片，每个瓦 片均含有自身左上、左下、右上以及右下四角的经纬度信息。将每个瓦片 以左下角顶点的经纬度数据命名，由于瓦片尺寸已知，因此在后续步骤中， 只需要从瓦片名称中提取提取左下角经纬度数据即可得到对应瓦片的四 角顶点经纬度数据。

步骤2，根据瓦片名称所提取的经纬度数据以及遥感图像中各个地区 的边界经纬度数据，对所有瓦片进行MapReduce操作，得到遥感图像中各 个地区的NDVI值。

步骤2-1，进行数据准备，获取南阳市各地区的边界经纬度数据，按 地理位置按从西向东，从北向南的地理位置排列数据的顺序。

对于各个瓦片，输入由步骤1所得到的瓦片以及各个瓦片四角顶点的 经纬度数据，采用NDVI计算软件ERDAS IMAGINE的调用脚本进行Map。

Map的输出为每个瓦片的NDVI值以及该瓦片中所包含的南阳市各个 地区的编号。

步骤2-2，对于每一个瓦片，Map的具体步骤如下：

步骤2-21，确定当前瓦片所在的地区。当前实施例中，所要确定的地 区为南阳市的下属区县。确定的方式如下：

步骤a，根据瓦片四角的经纬度计算瓦片中心点C的经纬度（Lng，Lat）, 其中Lng表示经度，Lat表示纬度。

步骤b，从准备数据中已排好序的边界经纬度数据中，自最西北角的 地区开始，将纬度＝Lat的边界点按经度值从小到大排序。

步骤c，用二分查找法找到最后一个经度值<Lng的边界点，该边界点 所在的经度线就是当前瓦片所属的经度线。

步骤d，沿步骤c所确定的经度线向南查找，直到找到最后一个纬 度>Lat的边界点，该边界点所在的地区就是当前瓦片所在的地区。

确定了每一个瓦片所在地区后，进入步骤2-22。

步骤2-22，计算各个瓦片的NDVI值。

调用NDVI计算软件ERDAS IMAGINE的“interpreter/Spectral  Enhancement/Indice”脚本，得到当前瓦片的NDVI值。

步骤2-23，将所得NDVI值进行归并和划分。

将获得的NDVI值按照区域编号输入相应的Reduce节点。这样就使得 同一个地区的NDVI值输入到了相同的节点进行处理。

Map得到各个瓦片的NDVI值之后，对输入的NDVI值进行Reduce操 作。

步骤2-3，Reduce的具体步骤为：

对于每个Reduce节点，将该节点输入的NDVI值相加；

将地区编号与该地区所对应的Reduce节点NDVI值之和所构成的数值 对输出，计算得到南阳市下属13个县以及县级市各自的NDVI值。

步骤3，对于各个地区，将其NDVI值输入至经过训练的神经网络中， 得到该地区的农作物产量估算值。

用遗传算法（GA）对传统的BP神经网络（ANN）进行优化，将优化 过的BP神经网络（GA-BP神经网络）进行训练，将训练后所得的神经网络 对南阳市各个地区的农作物产量进行估算。在当前实施例中，农作物为小 麦。

具体方法如下：

首先用GA（遗传算法）对传统的BP神经网络的拓扑结构、连接边的 权重、各个节点的阈值进行初始化。其中，为了保持种群的多样性，当前 实施例对传统的GA直接用子代替代父代生成下一代繁殖种群的方法进行 改造和优化，优化的方式用如下伪代码表示：

if(f_i>f_avg)

将该个体加入新繁殖种群中

上述伪代码表示对于第i个农作物的个体，如果其评价函数值f_i高于神 经网络的平均评价函数值f_avg，则将该个体加入到新繁殖种群中；否则进行 进一步判断：设定将丢弃个体的阈值=min{1，exp（-|f_i-f_avg|/T_k）}，即取1 和exp（-|f_i-f_avg|/T_k）中较小者，将随机生成的概率p与该阈值进行比较， 若随机生成的概率p小于该阈值，则将该个体加入新繁殖种群中，否则将 该个体丢弃，其中，k表示的是进化的代数，T_k表示的是一个随着k递减的 变量。

用各个地区历史的“NDVI－农作物产量”数值对作为训练集合对优 化过的神经网络进行训练，训练持续进行直到误差达到预设的误差阈值。

训练过程中的相对误差随时间的变化情况如图3所示，其中纵坐标为 误差，横坐标为训练的时间。用2008年的NDVI值作为ANN的输入，得到 的输出即为当年小麦产量的估算值。本发明得到的GA-BP模型与简单的BP 模型对2008南阳小麦产量估算所得到的相对误差的比较如图4所示。

图4中以方块表示的折线代表普通的BP模型预测的相对误差，以圆点 表示的折线代表GA-BP模型预测的相对误差。横坐标表示各个地区，纵坐 标表示相对误差。由图4可以明显地看到，无论是从产量估算的准确性来 看还是从估算结果的稳定性来看，基于GA-BP的模型都比简单的BP模型要 好。

本发明的实现过程中，已经将光谱遥感图像巨大的数据量问题、切图 的效率问题、NDVI值提取的效率问题、农作物估产的准确性和稳定性问 题等实际问题考虑在内，为农作物的遥感估产提供高效、可靠的解决方案。