基于三次样条插值的农机作业面积计算

李亚硕，贾晓峰，伏伟浩，赵 博，王瑞雪，苏仁忠，徐名汉，庞在溪

(1.中国农业机械化科学研究院集团有限公司，北京 100083； 2.土壤植物机器系统技术国家重点实验室，北京100083；3.中机建工有限公司，北京100083； 4.湖北省农业机械化技术推广总站，湖北 武汉 430017)

0 引言

随着国内农机化水平的提高和农机规模的扩大，田间作业全过程已基本实现机械化。农机作业面积是补贴发放、费用收取、效率评估、工资发放等的依据，精准有效的作业面积计算是农机管理部门和农机合作社管理人员的迫切所需。

在农机上安装机载定位设备可以记录作业轨迹信息，每隔几秒向平台发送一次数据，利用北斗定位信息提取行驶轨迹点[1-2]。如何利用这些散点还原农机真实行驶轨迹，并计算出真实作业面积，是值得关注和需要解决的问题。鲁植雄等[3]利用图像处理方法计算作业面积，同时标注重耕、漏耕区域，但该方法受轨迹点密集度和分辨率影响较大。有学者使用栅格法通过计算轨迹覆盖区域，减去空白区域和重叠区域统计作业面积[4-6]。但栅格大小不好确定，占用内存空间较大，并且边界会有毛刺。传统距离测量法是简单易行的作业面积计算方法，将轨迹点连成折线，计算折线总长度，根据农机作业轨迹长度和幅宽计算。当农机沿直线行驶时，该方法计算面积误差不大。当农机出现转向或角度发生变化时，转向处用折线计算会产生误差。有研究采用矢量缓冲区算法，将拐角处用圆弧代替，减少了拐角误差[7-8]。缓冲区法以幅宽作为转角直径不能真实代表所有转向幅度，对行驶过程中偶然转向有效，当地块是弧形时，缓冲区法仍将两点间轨迹用直线代替。

本研究基于农机运行轨迹特点及作业地块不确定性，将农机轨迹用弧线表示，代替传统直线法计算，更能反映农机行驶的真实情况。本研究设计了基于三次样条插值的作业面积计算方法，对轨迹点进行插值，计算得到的轨迹曲线长度，结合作业幅宽得到作业面积，准确率较缓冲区法和距离测量法都有提高。

1 农机作业面积计算

农机作业轨迹分道路行驶、田间行驶和原点停留，有效的作业面积为田间行驶轨迹，如图1所示。农机在田间作业多采用往返作业形式，即到地块边界后掉头反向作业。因此，计算作业面积时，除了意外停车轨迹点，也应将掉头转向点去除。保留的轨迹点为连续单向的线段，将这些线段分别插值计算再累加。

图1 农机作业轨迹Fig.1 Operation track of agricultural machinery

1.1 轨迹点预处理

机载设备在开机情况下，会定时发送定位信息。有一些特殊点会影响作业面积计算，应提前去除。

(1)农机停止散点。农机在停止时，轨迹点不是集中在一个点，而是散落在农机附近小范围内，称为农机停止散点，速度为0或者很小值，如图2所示，将这些轨迹点直接去除。

图2 停止散点和轨迹漂移点Fig.2 Dwell scatter point and trajectory drift point

(2)断点。由于设备故障问题，两点之间时间间隔会大于正常时间间隔，称为数据点丢失。将前一断点作为前段轨迹的终点，后一断点作为后段轨迹的起点。

(3)漂移点。由于定位精度原因，个别轨迹点与前后两点距离远大于正常距离间隔，称该点为漂移点。去除该漂移点，将漂移点前一点作为前段轨迹的终点，后一点作为后段轨迹的起点。

(4)农机在田间作业一般是往返行驶，到地块边界后会掉头行驶，地块边界转向部分不属于作业面积。某轨迹点与前一轨迹点和后一轨迹点分别连成直线，两直线夹角范围为0°～180°。当夹角大于一定阈值时，认为该轨迹点处于地块边界调转，将该点去除。

所有田间异常点处理后，农机田间作业轨迹被划分为多条沿地垄的单向独立轨迹，如图3所示。图3a为原始轨迹点图，图3b为根据轨迹点速度和转向角度，将地块转向轨迹点与田间行驶轨迹点标记为不同类型。按照时间序列，计算每条作业行地头之间作业轨迹点的长度，虽然由于信号问题导致个别轨迹点丢失，不影响整体轨迹长度计算。将计算得到的各段轨迹累加，得到所有轨迹总长度。

图3 作业轨迹预处理效果Fig.3 Effect of job track preprocessing

1.2 三次样条插值算法

两个点确定一条直线，最简单的方法就是每两个点之间连一条线，最后得到一种折线图。但是，这种方法显然具有很明显的缺陷：曲线不够光滑，连接点处的斜率变化过大。

如图4所示，将折线连接处用弧形代替，虽然解决了拐角处误差，但对于其余折线P1P2、Pn-1Pn、PnPn+1等处，若地块形状如图所示为弧形，用直线来代替弧形计算轨迹长度，也会产生误差。

图4 折线连接处用弧线代替Fig.4 Broken line connection is replaced by an arc

如图5所示，用折线P1P2+P2P3+P3P4+P4P5代替弧P1P5，必然会出现计算误差。因此需要尽可能还原出农机作业的真实轨迹，才能更准确计算行驶距离和作业面积。

图5 两点之间直线与弧线对比Fig.5 Comparison of straight line and arc between two points

如图6所示，离散的坐标点为轨迹点，穿过每个轨迹点的平滑曲线为基于三次样条插值得到的弧线，未完全经过每个轨迹点的弧形线为最小二乘法拟合出的曲线。计算农机作业面积，首先要计算农机轨迹线长度，每个采集的轨迹点都是真实点，数据点拟合曲线只能统计轨迹点趋势和规律，不能保证所有点都在统计内，因此三次样条插值算法得到的弧线更代表农机真实轨迹[9-11]。通过插值细分线段，将两点之间的线段平滑，得到的弧线作为作业轨迹线。

图6 拟合和插值方法效果对比Fig.6 Comparison of fitting and interpolation methods

对于区间[a，b]，每个点坐标为[xi，yi]，i=0，1，…，n。将区间分成x0=a

S0(x)=y0+b0(x-x0)+c0(x-x0)2+d0(x-x0)3
S1(x)=y1+b1(x-x1)+c1(x-x1)2+d1(x-x1)3

(1)

……

Sn(x)=yn+bn(x-xn)+cn(x-xn)2+dn(x-xn)3

其中，S0(x)在区间[x0，x1]上，S1(x)在区间[x1，x2]上，……，Sn(x)在区间[xn-1，xn]上。函数S(x)若为区间[a，b]上的三次样条插值函数，需满足以下条件。

(1)在每个子区间[xi-1，xi]上S(x)是三次多项式。

(2)S(xi)=yi，i=0，1，…，n。

(3)在区间[a，b]上S(xi)的二阶导数S″(xi)连续。

Si(x)=aix3+bix2+cix+di，i=0，1，…，n

(2)

S(xi)=yi，S(xi-0)=S(xi+0)，i=0，1，…，n

(3)

S′(xi-0)=S′(xi+0)

(4)

S″(xi-0)=S″(xi+0)，i=0，1，…，n

(5)

式(2)中的ai、bi、ci、di待定，需满足式(3)、式(4)和式(5)。

式(3)、式(4)和式(5)共给出4n-6个条件，需要待定4(n-1)个系数，因此要唯一确定三次插值函数，还要附加2个边界条件。

S′(x0)=y1′，S′(xn)=yn′

(6)

S″(x0)=y1″，S″(xn)=yn″

(7)

2 试验及结果分析

试验采用弧形地块农机作业数据作为数据集，地块内作业面积为122 544 m2。在农机上安装采集设备，每5 s采集一次作业点数据，通过无线发送至后端服务器。每个轨迹点包含农机作业点经纬度、车辆行驶速度等信息。通过三次样条插值方法计算农机作业轨迹长度，乘以农机幅宽得到农机作业面积。对比试验为采用折线累加、矢量缓冲区算法和最小二乘法拟合算法3种方法计算作业轨迹长度。

试验步骤如下。

(1)轨迹点预处理。

(2)利用三次样条插值法获取农机作业轨迹线。

(3)计算单条轨迹线长度。

(4)累计轨迹线总长度，乘以幅宽，得到总作业面积。

准确率计算方式为

(8)

式中S——真实测量面积，m2

S0——计算面积，m2

采用不同算法计算的作业面积结果如表1所示。本文提出的三次样条插值法在该地块的作业面积计算中准确率高于其他几种方法。作业缓冲区算法在处理有转向弧度的作业轨迹时，准确率要优于折线累加法，而曲线拟合方法准确率低于折线累加法。

表1 不同算法计算作业面积结果Tab.1 Calculation results of working area by different algorithms

3 结论

针对农机在弧形地块作业时，作业面积依靠轨迹点连成折线方式计算误差较大的问题，提出了利用三次样条插值法对轨迹点进行插值，将轨迹线用弧线表示。该方法既能保证计算涵盖所有轨迹点，也根据已有轨迹点信息还原真实农机运行轨迹。试验结果表明，本文方法面积计算准确率为97.83%，高于其他面积计算方法。