一种基于BDS/GPS定位数据时空特征的农机作业面积计算方法与流程

本发明涉及农业应用技术领域，具体涉及一种基于bds/gps定位数据时空特征的农机作业面积计算方法。

背景技术：

随着农业机械化、现代化、智能化的发展，智能农机技术广受关注。在农机跨区进行耕种管收作业时，作业面积作为关键数据至关重要。在精细农业中，作业面积也是一个很重要的指标，不仅决定了种子、农药、化肥等资源的投入量，还用于确定最终作业工作量和最终收费。但是农机作业地形具有极大的复杂性，很多地块面积无法准确测量统计，导致农民的相关利益受损，政府部门也无法准确统计地块面积进而影响到决策和管理。

为此，在农机上安装地理信息系统、空间定位系统等数字技术导航装置，结合后台服务系统的模式，可以满足农机的深松、深翻、旋耕、联合收割、秸秆还田、农机植保等多种类型作业的远程监控、管理需求，实现农业生产耕种收的精准作业以及精度高、可靠性高、实时便捷的农机作业面积计算结果。

现有的农业平台计算农机作业面积的方式概述为：通过农机上安装的传感器的bds/gps定位模块实现对农机的实时定位，上传轨迹点至平台；通过农机上安装的传感器实时获取作业的状态，并上传平台；作业结束后，平台根据各点的作业状态，对合格的轨迹点根据其数据特征进行画线及分地块，再计算作业面积。

目前，农业平台在作业面积计算上一直存在三个主要问题：1、非耕地作业区域的路上无效面积。有些农机在路上的轨迹线也被判定成作业状态，会产生无效面积，而且轨迹线越长，产生的无效面积越大。2、轨迹点漂移产生的无效面积。有些农机作业的轨迹点会出现漂移，导致地块形状失真、非地块区域有作业面积，会不同程度的增加无效面积，同时也会导致部分正常作业面积的丢失。3、地块划分不准确，如一个地块被划分成多个，或者多个地块被划分成一个，导致判断重耕和漏耕现象时的结果不准确。

针对上述三个技术问题，现有技术中：对于路上行驶或地块间行驶且标记为正常作业状态的轨迹数据，无法做到有效识别或识别效率低下，在现有的一种农机作业面积获取方法及装置中，虽然对dbscan算法的邻域半径做了优化，但仅考虑速度，未考虑农具幅宽对算法的影响，未对密度阈值做优化，也未对算法的效率问题提出有效改进措施；对于漂移数据，现有方法仅依靠速度是否异常来判别，如速度过大、速度为0或1km/h识别零点漂移，无法识别因终端原因导致漂移数据中速度与正常作业时差异较小的情况，去除不彻底；因漂移导致作业数据会有部分丢失，从而影响面积计算，现有方法只补偿转弯点，未进行此类数据的补偿；地块划分不准确或划分准确但执行效率很低。

技术实现要素：

本发明要解决的技术问题是提供一种基于bds/gps定位数据时空特征的农机作业面积计算方法，可以有效识别路上轨迹、地块间轨迹，并进行去除，还可以识别不同原因导致的多种类型的漂移数据，并做到有效去除，提升原地块划分算法的执行效率和效果，准确划分地块并补偿缺失数据以得到精确的作业面积。

为解决上述技术问题，本发明采用的技术方案是：

一种基于bds/gps定位数据时空特征的农机作业面积计算方法，包括以下步骤：

步骤a、获取农机作业的轨迹数据集d；

步骤b、数据的预处理：去除轨迹数据集d中速度不正常的漂移点、停车点以及位置完全相同的冗余轨迹点数据，并将轨迹数据集d中的轨迹点按定位时间排序；

步骤c、地块划分以及道路轨迹点的去除：

步骤c.1、将轨迹数据集d切分为k个子集，每个子集中具有m个数据，同时记录子集di的分割点，i＝1,...,k，pi为di的最后一个点；

步骤c.2、对每个数据子集di分别建立kd树并保存数据在树中的结构；

步骤c.3、对每个数据子集di运用dbscan+kd树搜索算法对子集内的数据进行聚类、保留核心点并去除噪声点；对每个类进行地块编号并统计地块数量；

步骤c.4、更新分割点集：聚类完成后，若某一个pi已不存在，则按排序的时间选取原di中与pi最近的点作为pi；判断pi与di+1第一个点的时间间隔，若小于轨迹数据上报时间间隔t的a倍则合并pi与di+1第一个点所属的类，然后更新每个类的标号及地块数量；

步骤d、漂移数据二次过滤：遍历每个地块，运用农机航向差分数理统计判定算法识别漂移地块并进行去除，更新各地块信息；

步骤e、地块缺失数据还原：遍历每个地块，采用农机轨迹线性插值算法将地块中丢失的轨迹数据进行还原；

步骤f、根据过滤后的农机作业轨迹数据确定农机作业田块数量并分别计算各个农机作业田块的面积。

本发明的有益效果是：1、提升原地块划分算法的执行效率和效果，准确划分地块；2、运用算法有效识别路上轨迹、地块间轨迹，并进行去除；3、识别并去除各种类型的漂移数据，如停车静态漂移、作业过程中的动态漂移等，并对无法直接去除的漂移数据进行二次过滤。

下面结合附图对本发明进行详细说明。

附图说明

图1是本发明农机作业面积的计算方法的流程图；

图2是正常作业、路上有轨迹、停车有漂移现象都存在的实例图；

图3是还原图2的原始轨迹数据仿真图；

图4是运行本发明的改进聚类算法划分地块及去除路上轨迹的效果图；

图5是农机在原地停车产生漂移、正常作业的情况的角度差分情况图；

图6是500条漂移数据的分布图；

图7是正常作业500条数据的分布图；

图8是漂移情况下航向差分数据的频数直方图；

图9是正常作业的航向差分数据的频数直方图；

图10是运行本发明的农机航向差分数理统计算法去除漂移的效果图；

图11是农机全天轨迹均是路上行驶时采用本发明的方法处理前后的对比图；

图12是农机正常作业时有路上轨迹采用本发明的方法处理前后的对比图；

图13是农机正常作业过程中有作业轨迹但出现大量漂移点的情况，采用本发明的方法处理前后的对比图；

图14是农机全天轨迹均是原地停车出现漂移的情况，采用本发明的方法处理前后的对比图。

具体实施方式

参见附图1，本发明提供了一种基于bds/gps定位数据时空特征的农机作业面积计算方法，包括以下步骤。

步骤a、获取农机作业的轨迹数据集d。轨迹数据集d一般是农机一天所采集的数据。

目前农机终端以设定的时间间隔t(如每2s一次)向平台发送农机的轨迹数据，以形成农机每日的轨迹日报表。农机在全天不关机的情况下，一天最多可上报43200条数据。针对本方法所要解决的几个问题，所选取的数据字段主要有每辆农机的标识号、每次定位的经纬度、实时速度、实时航向、农具幅宽等数据。

需要说明的是：1、实时速度与实时航向为bds/gps定位模块直接提供，终端未做处理，其数据在卫星信号较差的情况下会失真；2、受bds/gps定位精度影响，经纬度数据会存在一定偏差，农机在同一位置停留，每次的经纬度数据也可能会不同，因而会产生小范围的漂移，而此时轨迹点的速度会在正常范围内(0～15km/h)，航向也没有明显变化，会被判定为正常作业状态，从而产生无效面积；3、不同作业类型的农具幅宽差异明显，而同一作业类型，不同型号的农具幅宽也不相同。由于农具幅宽和作业轨迹点的稀疏程度直接相关，这个字段会是本方法中很重要的一个参数。

步骤b、数据的预处理：去除轨迹数据集d中速度不正常(小于0或大于100km/h)的漂移点、停车点(速度为0)以及位置完全相同的冗余轨迹点数据，并将轨迹数据集d中的轨迹点按定位时间排序。

被过滤的数据不仅对算法没有用处，还会对执行效率和处理效果造成很大影响，过滤后不仅保证了数据质量，还可大幅提高后续算法处理的效率。这里以某一农机某一天的作业情况为例，如图2，将数据量压缩后，可从之前的43200条骤减至约7520条。

步骤c、地块划分以及道路轨迹点的去除。

目前的dbscan算法聚类效果较好，原理简单较易实现；能去路上面积，能划分地块，还能识别少部分的漂移点，用途很广。但缺点也很明显：1、执行效率非常低，每个点要与所有其他点进行计算比较，如果有n个点，那么复杂度就是n²，经测试，对一辆车一天产生的43200条数据运算的时间约20分钟，在现有资源条件下无法实现工程化；2、算法中的两个参数ε和minpts太固定，设置大了会存在将路上轨迹识别成正常作业轨迹的问题，还会把相邻但不同的地块划分到一起，设置小了正常作业轨迹会存在被识别成路上轨迹而去除的风险，自适应性很差。

基于kd树的搜索方法能够明显提高效率，使原来的n²变成nlogn，本方法在此基础上，借鉴大数据“分而治之”的并行处理思想，做了进一步改进：对数据进行分段处理，即在一开始就对数据进行等间隔的切分，将n个数据分成k份，每份有m个数据，即n＝km，再在每份数据中做基于kd树的dbscan聚类，总的执行情况为：能够最大程度的减少算法中最耗时的查询与比对的工作。由于切分数据会使原本属于同一类的数据被分成了两类，那么当聚类完成后，只需判断每个分割点与它相邻的上一个点的时间间隔即可将原本属于同一类的两类点再聚成一类。

详细的，本步骤包括：

步骤c.1、将轨迹数据集d切分为k个子集，每个子集中具有m个数据，本例中，将预处理后的7520条数据可切分为10个子集，每个子集由752个数据组成，同时记录子集di的分割点pi，i＝1,...,k，pi为di的最后一个点。本步骤中是直接先将数据硬性分割成10类，在此基础上运用聚类算法对每个分割的数据内部再进行准确聚类完成地块划分，之后再通过步骤c.4把这第一步不该分割的地块再重新连接起来，这样能最大程度提高计算效率。

步骤c.2、对每个数据子集di分别建立kd树并保存数据在树中的结构。kd树方法分为建立kd树和遍历kd树两步，本步骤将数据以树形结构存储，之后在每次聚类过程中，只需按照此结构查找待查找点的相邻点即可。

步骤c.3、对每个数据子集di运用dbscan+kd树搜索算法对子集内的数据进行聚类、保留核心点并去除噪声点；对每个类进行地块编号并统计地块数量。

搜索算法中的半径参数ε和密度阈值minpts参数设定为：

ε＝max(k1ε1,k2ε2)，

其中，ε1为前后两个轨迹点的间隔，ε1＝vt，v为轨迹点的速度，t为轨迹数据上报时间间隔，本例中为2s，ε2为农机的农具幅宽，k1、k2为放大系数，取值范围都在区间[1,2]，可根据实际情况做微调，η和μ为基础参数，一般根据实际应用场景确定，η取值范围在区间[3,10]，μ取值范围在区间[2,4]，在本实施例中η＝5,μ＝2，也可根据实际情况做相应调整。两个参数的确定方式就是按照地块中的点的特征去设计的，这样保证了地块内的轨迹点会满足核心点条件，即当k1ε1≥k2ε2时，地块内的每个点周围至少有5个点，而路上的每个轨迹点在自己的邻域内只包含不多于3个点，从而被识别成噪声而去除；当k1ε1＜k2ε2时，地块内的核心点周围会按照两者的比值等比例放大密度阈值，而路上的每个轨迹点邻域半径内虽然也会增加相同的点数，但依然不会达到核心点的要求。本例中切分数据后聚类结果为13类。

步骤c.4、更新分割点集：聚类完成后，若某一个pi已不存在，则按排序的时间选取原di中与上述pi最近的点作为pi；判断pi与di+1第一个点的时间间隔，若小于轨迹数据上报时间间隔t的a倍(a的取值范围为3～5，可以设为5，认为相隔10s以上就不是一个地块)则合并pi与di+1第一个点所属的类，然后更新每个类的标号及地块数量。数据被切分会导致同一地块被划分成多块，这一步用来判断哪些地块本应属于一块，并重新合并成同一块。重新更新地块编号与地块数量。本例中10个分割点有7个与其相邻点时间差小于5t，故重新合并点集后聚类结果为6类，所剩点的总数为5024。

步骤d、漂移数据二次过滤：遍历每个地块，运用农机航向差分数理统计判定算法，识别漂移地块，并进行去除，更新各地块信息。

通过步骤b和步骤c已将路上轨迹有效去除，形成形状规则的地块，如图4，并且去除一部分速度数据异常的漂移点，但是一些停在房屋附近受遮挡导致停车漂移的现象无法用上述方法滤除，因该部分数据中速度和航向体现的特征有时与正常状态下无异。通过观察发现，这种漂移数据通常与正常作业地块相分离，会被识别为单独的地块，如图2和图3所示。正常作业和原地漂移的本质区别在于：正常作业轨迹一般是沿着直线或近似直线行驶，偶尔转弯或掉头，其前后两点的航向角度差通常在15°以内，而漂移则可能每时每刻航向都在无规律的变化，航向角度差会在[0,180]°之间近似于均匀分布。如图5，前500条数据发生在农机在原地停车产生漂移，后面是农机正常作业的情况的角度差分情况，可以看出，前后差异很明显。图6放大显示了500条漂移数据的情况，图7选取了正常地块中500条数据，与图6做更鲜明的对比；图8和图9分别以频数直方图的方式展示了漂移和正常作业的航向差分数据的区别，据此，根据数据的数理统计规律，可用方差来作为区分漂移和正常作业的指标。

本方法在上述聚类算法划分地块后，以地块为对象，考虑到发生漂移现象时，实施航向数据有可能失真，因此用墨卡托坐标转换方法将经纬度转换为直角坐标后，用每两点之间求取的航向代替原来的实时航向，再对航向进行差分，运用农机航向差分数理统计判定算法识别漂移地块并进行去除，具体步骤如下：

步骤d.1、遍历地块，对任一地块，将位置数据的经纬度坐标用墨卡托法转换成直角坐标：[经度，纬度]-->[x,y]；

步骤d.2、统计该地块的轨迹点数l，求每相邻两点间斜率：ki＝(x(i+1)-x(i))/(y(i+1)-y(i))，i＝1,...,l-1；

步骤d.3、将斜率转换成航向角度：heading(i)＝atan(ki),i＝1,...,l-1,并且需要根据x(i+1)-x(i)、y(i+1)-y(i)的正负转换成[0,360)°范围内；

步骤d.4、对相邻两点做角度差分：

heading_diff(i)＝min(abs(heading(i+1)-head(i)),360-abs(heading(i+1)-heading(i)))，i＝1,...,l-1；其中，abs()表示对数据取绝对值；

步骤d.5、将[0,180]°划分成q个区间，每个区间宽度为统计角度差分在每个区间的个数并生成一组q个数的数组。

具体地，对[0,180]°以每15°做区间划分，共12个区间，统计heading_diff在每个区间的个数，生成一组12个数的数组num；

步骤d.6、对数组归一化处理：h＝num/(l-1)；

步骤d.7、对归一化的数组求方差s＝var(h)，s<δ(本实施例中δ取0.01)则判定为漂移地块，删除该地块信息，对剩余地块更新编号。本例中6个地块有一个被判定为漂移地块，并进行了过滤，剩余5个地块。

步骤e、地块缺失数据还原：遍历每个地块，采用农机轨迹线性插值算法将地块中丢失的轨迹数据进行还原。

步骤d已将漂移地块进行了有效去除，如图10，此时所有无效数据均已去除，只需对各地块进行数据补偿即可计算面积。具体包括以下步骤：

步骤e.1、对任一地块，轨迹点数为l，将轨迹点集的上报时间做差分t_diff(i)，i＝1,...,l-1；

步骤e.2、若某t_diff(i)>3t，t＝2s为轨迹数据上报时间间隔，按照坐标位置插入个坐标点。

步骤f、根据过滤后的农机作业轨迹数据确定农机作业田块数量并分别计算各个农机作业田块的面积。

以下各实施例为除本例外，本发明对各种不同情况所做处理后的效果：

实施例一：农机全天轨迹均是路上行驶的情况，参见附图11。采用本发明的方法后已完全去除路上轨迹。

实施例二：农机正常作业时，既有路上行驶轨迹，也有作业轨迹，参见附图12。采用本发明的方法后路上无效轨迹也已被去除。

实施例三：农机正常作业过程中，有作业轨迹，但出现大量漂移点，参见附图13。采用本发明的方法后已完全去除漂移，保留了正常作业轨迹。

实施例四：农机全天轨迹均是原地停车出现漂移的情况，参见附图14。采用本发明的方法后已完全去除整个漂移地块。

最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制；尽管参照较佳实施例对本发明进行了详细的说明，所属领域的普通技术人员应当理解：依然可以对本发明的具体实施方式进行修改或者对部分技术特征进行等同替换；而不脱离本发明技术方案的精神，其均应涵盖在本发明请求保护的技术方案范围当中。