CN109029243B - 一种改进的农机作业面积测量终端和方法 - Google Patents

Abstract

本发明涉及一种能较准确地测量农机作业面积，可进行工况数据采集与保存，并将相关信息发送至监控中心的面积测量终端。本发明在硬件设计上由主控芯片STM32F103ZET6、ATK‑S1216F8‑BD GPS/北斗模块、电源模块、人机交互模块(LCD触摸显示屏和键盘)、GPRS模块、存储模块、NRF24L01模块和CAN总线接口等构成，其价格便宜且功能多用，解决了以往测亩仪功能单一和成本较高的问题。

Description

一种改进的农机作业面积测量终端和方法

技术领域

本发明旨在做出一种能较准确地测量农机作业面积，可进行工况数据采集与保存，并将相关信息发送至监控中心的车载终端。

背景技术

在“精准农业”的潮流下，由于田块的大小直接决定了生产资料的投入量，因此，准确地计算农机作业的面积对指导农机高效作业、控制生产资料的投入，以及计算农作物的单位面积产量等具有很大意义。农机机主和农户都希望给农用机械安装一个“计亩器”，就像城市里的出租车有个“记价器”，耕作了多少面积的土地能实时计算显示出来。国内学者对测亩仪表的研究简要情况如下。

在农业机械化的大背景下，测亩仪的发展方向是与机械同步使用，进行实时测量；且如今我们的生活与网络也息息相关，如若测亩仪具备了联网与采集农机运行状态信息的功能，这也有利于实现监控中心对农机的远程监控；另一方面，较高的测量精度和较合适的价格对于农机机主来说也是非常重要的考虑因素。因此，本发明主要是基于STM32F103ZET6和北斗/GPS双模接收机等模块，结合相应算法，设计出一种低成本测亩仪，在此基础上，进一步扩展其数据采集、存储与传输的功能，提高其性价比，更好地适应智能化农业的发展。

发明内容

本发明旨在做出一种能较准确地测量农机作业面积，可进行工况数据采集与保存，并将相关信息发送至监控中心的车载终端。

为了实现上述目的，本发明一种改进的农机作业面积测量终端，其特征在于，所述终端由手持式部分和基座两部分构成，所述手持式部分包括：主控芯片、GPS/北斗模块、电源模块、人机交互模块和NRF24L01模块，所述基座部分包括电源转换模块、GPRS 模块、存储模块和CAN总线接口。

更进一步的，所述基座底部为吸盘结构，所述手持部分与基座通过插片和插孔连接。

更进一步的，本终端使用位置差分技术，即使用两个北斗/GPS接收机，其中一个作为固定基站，通过NRF24L01射频模块向另一个测量终端发送实时经纬度误差。

更进一步的，本终端使用自适应交互多模型卡尔曼滤波算法，即根据农机工作状态使用不同卡尔曼滤波模型，对北斗/GPS接收机接收回来的经纬度进行修正。

一种改进的农机作业面积测量方法，其特征在于，包括以下步骤：

步骤1：坐标转换和面积计算；本发明通过北斗GPS双模接收机模块获取当前农机的经纬度信息，经过高斯克吕格投影方式进行坐标转换，再利用完全离散格林公式对不规则农田面积进行计算。

步骤2：提高低成本接收机的定位精度；本发明采用差分的方法，消除电离层误差等公共误差；结合卡尔曼滤波的算法，对经纬度数据进行预测与修正，进一步减弱误差；在卡尔曼滤波上增加自适应因子，避免因为量测设备故障或者信息传输间隔中断而导致的不能正常提供状态测量信息的问题，提高卡尔曼滤波的鲁棒性；构造多模型卡尔曼滤波算法，针对目标不同的运动状态使用不同的模型匹配，与各个模型对应的滤波器并行工作，实现不同模型的切换，避免因为使用单一模型而造成的偏差；

步骤3：运用两种面积测量方法；一是航迹法，此时农机利用回耕法进行作业，即农机先沿农田边界作业，再逐次往内作业，在农机沿外边界绕行一周结束后，利用测亩仪计算出这一周的周长和所围面积；二是定幅法，当农田形状比较规则，农机按照梭形耕法进行作业时，默认农机的作业幅宽为固定值，用幅宽乘以农机走过的距离得到面积；

步骤4：数据的存储与传输；本发明通过CAN总线获取农机工况信息，一方面通过本地数据存储设备(SD卡)将数据保存在车上，另一方面通过GPRS网络发送至指定IP 的监控中心平台上，实现监控中心对车辆工况的实时监测和对突发状况的处理。

更进一步的，步骤1所述经过高斯克吕格投影方式进行坐标转换，再利用完全离散格林公式对不规则农田面积进行计算，具体为：

高斯投影是一种横切椭圆柱投影，假想有一个椭圆柱与地球椭球体上某一经线相切，并且将椭球筒的中心轴放置在赤道平面上，按等角条件将地球表面投影到椭球筒上，然后将椭球筒展开，此投影为高斯投影；在高斯投影变换中，中央经线投影后的直线为X 轴，赤道投影后的直线为Y轴，它们的交点为原点，为使转换之后的Y轴坐标不为负值，通常在Y轴上加上500km，具体变换公式如下:

式中，B—纬度；l＝L-L₀,L为经度，L₀为轴子午线的经度；N—卯酉圈曲率半径；t＝tanB；η＝e'cosB；e'—椭球第二偏心率；X₀—当l＝0时，从赤道起算的子午线弧长；参数计算如下： e²＝(a²-b²)/a²e'²＝(a²-b²)/b²

X₀＝a(1-e²)(A₀B+A₂sin2B+A₄sin4B+A₆sin6B+A₈sin8B)

其中：

为减少投影变形，保证投影变换的精度，通常采用分带投影的方法，即按经度把椭球分为很多带，各带分别投影；在精细农业实践中，通常采用3度分带投影，从东经1°30’开始，每3°为一带，将全球划分为120个投影带，利用高斯投影的公式，编写代码求得转换后的平面坐标X，Y；

面积测量利用完全离散格林公式，其具体内容如下：

对于一个封闭的区域，如果获得其边界上n个坐标点，如(X₁,Y₁),(X₂,Y₂),…(X_n,Y_n)则其面积可以由下面的格林公式计算获得：

其中i＝1,2,…n，当i＝n时，X_i+1＝X₁,Y_i+1＝Y₁

周长可由以下公式得到：

更进一步的，步骤2所述差分的方法和卡尔曼滤波的算法具体为：

因为农机的工作面积相对较小，两台接收机接收到的误差基本相同，因此本发明借鉴位置差分的原理，设计测量装置时采用两台接收机，一台作为固定参考站，一台作为移动站；测量面积前，参考站放在固定位置多次接收定位信号，取平均值获得相对参考坐标(B_r,L_r)作为标准值，测量面积时，参考站实时接收定位信号(B_x,L_x)，并与之前获得的相对参考坐标相减计算误差，得出实时误差值(ΔB_x,ΔL_x)，参考站将实时误差值发送给移动站，移动站利用实时误差值对经纬度数据进行修正，从而减小外界的干扰，消除参考站和移动站的公共误差，例如卫星轨道误差、电离层误差等，实现定位精度的提高；

卡尔曼滤波的实质是由量测值重构系统的状态向量，它以“预测—实测—修正”的顺序递推，根据系统的量测值来消除随机干扰，再现系统的状态，或根据系统的量测值从被污染的系统中恢复系统的本来面目；

卡尔曼滤波的算法流程如下：

(1)设定状态与协方差初值X₀和P₀；

(2)时间更新；根据上一时刻的状态变量估算值和误差协方差阵估算值，结合线性化的状态方程，更新当前时刻的状态变量和相应误差协方差阵；公式如下：

P_k-＝FP_k-1F^T+Q

式中：F为状态传递矩阵，

为k-1时刻的状态量，K为过程噪声增益矩阵，v为过程噪声扰动量，其相应的协方差为Q；

(3)测量状态更新；

首先由以下公式计算卡尔曼增益

K_k＝P_k ^-H^T(HP_k ^-H^T+R)^-1

式中，H为量测矩阵，R为量测噪声协方差矩阵；

接着由观测变量z_k更新估计，

再更新误差协方差，

P_k＝(1-K_kH)P_k ^-

在原有卡尔曼滤波的基础上，我们增加了一个由比例因子组成的比例矩阵S_k+1，来调整量测协方差矩阵，让理论值接近真实值，关系式如下：

其中，δ是移动窗口的宽度，上式的左边实际上代表实际的新息协方差，右边代表理论的新息协方差，经调整可以得到：

在某一个传感器发生量测故障，而其他传感器正常工作的情况下，由比例因子组成的比例矩阵对应的项将变得相对更大，将增大传感器在R矩阵的量测噪声的协方差，最终故障测量将不会被滤波器理会；另一方面，比例因子只有在测量故障的情况下，才会影响系统的估算过程，在正常的运作情况下，这个比例矩阵将是一个单位矩阵 S_k+1＝I_z×z，其中z为更新向量的维数；

矩阵S_k+1可以用上式进行计算，但有时对角元素可能为负或者远小于1，因为δ是一个有限的数，且计算过程中会有估算和圆滑误差；矩阵S_k+1是一个对角矩阵，因为只有对角元素才会对调整有影响，此外矩阵S_k+1不能有元素小于1，因为R矩阵不能针对特殊问题快速递减，因此，按照如下公式计算比例矩阵：

根据观测量对当前时刻的状态变量预估值进行修正，得到当前状态的估计值及相应的误差协方差矩阵；

说明交互多模型的算法思想；在实际工作过程中，农机并非一直定常速行驶，在转弯或者其他状况下时，速度会出现一定的变化，因此卡尔曼滤波中的定常速模型会出现一定的偏差，交互多模型算法的基本思想是针对目标不同的运动状态使用不同的模型匹配，与各个模型对应的滤波器并行工作，通过多个模型的概率加权构建更加符合实际的系统模型，滤波器最终的估计结果为各个单独模型估计的加权和；

本发明在卡尔曼滤波中引入了交互多模型，构造了一种多模型卡尔曼滤波算法，即 IMM-EKF算法，具体流程如下：

(1)计算IMM-EKF算法的滤波初值；初始状态向量

初始协方差矩阵P_i(0|0)和初始模型概率μ_i(0)；

(2)状态估计的交互式作用

IMM-EKF算法状态估计的交互式作用后得到

和P_oi(k|k)，具体计算公式如下

其中，

(3)EKF算法

将

和P_oi(k|k)作为k+l时刻第i个模型的输入，进行EKF，得到滤波输出

和P_i(k+1|k+1)；具体步骤如下：

计算状态和协方差一步预测：

P_i(k+1|k)＝F_i(k)P_oi(k|k)F_i(k)^T+G_i(k)Q_i(k)G_i(k)^T

计算量测预测值：

计算新息协方差矩阵：

计算滤波增益矩阵：

计算状态和协方差更新值：

其中，

是非线性方程的雅可比矩阵，如下所示：

(4)模型概率更新

若模型i的滤波新息为v_i(k+1)，相应的新息协方差为S_i(k+1)，并假定服从高斯分布，那么模型i的可能性为

其中，

则模型i的概率更新为：

其中，

(5)滤波交互输出

本发明的有益效果为：本发明一种改进的农机作业面积测量终端和方法可进行工况数据采集与保存，并将相关信息发送至监控中心的面积测量终端。本发明在硬件设计上由主控芯片STM32F103ZET6、ATK-S1216F8-BD GPS/北斗模块、电源模块、人机交互模块、GPRS模块、存储模块、NRF24L01模块和CAN总线接口等构成，其价格便宜且功能多用，解决了以往测亩仪功能单一和成本较高的问题。

附图说明

图1是本发明的总体系统组成；

图2是本发明的思维导图；

图3是本发明的上层控制算法程序流程；

图4是本发明高斯克吕格投影原理图；

图5是本发明基于格林公式的面积计算示意图；

图6是本发明位置差分的基本原理图；

图7是本发明卡尔曼滤波的基本流程图；

图8是本发明鲁棒卡尔曼滤波的流程图；

图9是本发明利用多模型进行滤波的流程图；

图10是本发明两种面积测量方法；

图11是本发明车载式与手持式结构简要外观图；

图12是本发明基于CAN总线和GPRS的农机远程监控系统的示意图。

具体实施方式

下面参考附图，详细说明本发明的具体实施方式。

1、参考图1，说明本发明的总体系统组成。本发明主要分为硬件系统和软件系统两大子系统，硬件系统主要由主控芯片、北斗/GPS双模接收机、人机交互模块(LCD触摸显示屏和键盘)、电源模块、存储模块、NRF24L01无线传输模块、GPRS模块等构成，软件系统主要由底层硬件驱动程序和上层控制算法构成。

2、参考图2，说明本发明的设计思维导图。本发明主要实现面积测量与数据传输两项功能，面积测量方面，本发明利用高斯克吕格投影方法和格林公式实现不规则土地面积的计算，利用卡尔曼滤波、自适应因子和交互多模型提高面积测量的精度；数据传输方面，本发明利用CAN总线获取农机工况信息，利用SD卡存储相关信息，并利用GPRS 模块实现农机与监控中心的数据传输。

3、参考图3，说明软件系统中上层控制算法的程序流程。此程序主要完成三个部分的功能，分别为面积测量，数据存储与数据传输。用户使用时，在人机交互界面下达指令，控制程序执行实现相应的功能。

4、参考图4，说明基于高斯克吕格投影方法的坐标转换。高斯投影是一种横切椭圆柱投影，假想有一个椭圆柱与地球椭球体上某一经线相切，并且将椭球筒的中心轴放置在赤道平面上，按等角条件将地球表面投影到椭球筒上，然后将椭球筒展开，此投影为高斯投影。在高斯投影变换中，中央经线投影后的直线为X轴，赤道投影后的直线为Y 轴，它们的交点为原点，为使转换之后的Y轴坐标不为负值，通常在Y轴上加上500km，具体变换公式如下:

式中，B—纬度；l＝L-L₀,L为经度，L₀为轴子午线的经度；N—卯酉圈曲率半径；t＝tanB；η＝e'cosB；e'—椭球第二偏心率；X₀—当l＝0时，从赤道起算的子午线弧长；参数计算如下： e²＝(a²-b²)/a²e'²＝(a²-b²)/b²

X₀＝a(1-e²)(A₀B+A₂sin2B+A₄sin4B+A₆sin6B+A₈sin8B)

其中：

为减少投影变形，保证投影变换的精度，通常采用分带投影的方法，即按经度把椭球分为很多带，各带分别投影。在精细农业实践中，通常采用3度分带投影，从东经1°30’开始，每3°为一带，将全球划分为120个投影带，利用高斯投影的公式，编写代码求得转换后的平面坐标X，Y。

5、参考图5，说明利用格林公式的面积测量方法。面积测量利用完全离散格林公式，其具体内容如下：

对于一个封闭的区域，如果获得其边界上n个坐标点，如(X₁,Y₁),(X₂,Y₂),…(X_n,Y_n)则其面积可以由下面的格林公式计算获得：

其中i＝1,2,…n，当i＝n时，X_i+1＝X₁,Y_i+1＝Y₁

周长可由以下公式得到：

6、参考图6，说明差分的原理与实现。差分实际上是在一个测站对两个目标的观测量、两个测站对一个目标的观测量或一个测站对一个目标的两次观测量之间进行求差，其目的在于消除公共项，包括公共误差和公共参数。因为农机的工作面积相对较小，两台接收机接收到的误差基本相同，因此本发明借鉴位置差分的原理，设计测量装置时采用两台接收机，一台作为固定参考站，一台作为移动站。测量面积前，参考站放在固定位置多次接收定位信号，取平均值获得相对参考坐标(B_r,L_r)作为标准值，测量面积时，参考站实时接收定位信号(B_x,L_x)，并与之前获得的相对参考坐标相减计算误差，得出实时误差值(ΔB_x,ΔL_x)，参考站将实时误差值发送给移动站，移动站利用实时误差值对经纬度数据进行修正，从而减小外界的干扰，消除参考站和移动站的公共误差，例如卫星轨道误差、电离层误差等，实现定位精度的提高。

7、参考图7，说明实现卡尔曼滤波的基本流程。卡尔曼滤波的实质是由量测值重构系统的状态向量。它以“预测—实测—修正”的顺序递推，根据系统的量测值来消除随机干扰，再现系统的状态，或根据系统的量测值从被污染的系统中恢复系统的本来面目。

卡尔曼滤波的算法流程如下：

(1)设定状态与协方差初值X₀和P₀。

(2)时间更新。根据上一时刻的状态变量估算值和误差协方差阵估算值，结合线性化的状态方程，更新当前时刻的状态变量和相应误差协方差阵。公式如下：

P_k ^-＝FP_k-1F^T+Q

式中：F为状态传递矩阵，

为k-1时刻的状态量，K为过程噪声增益矩阵，v为过程噪声扰动量，其相应的协方差为Q。

(3)测量状态更新。

首先由以下公式计算卡尔曼增益

K_k＝P_k ^-H^T(HP_k-H^T+R)^-1

式中，H为量测矩阵，R为量测噪声协方差矩阵。

接着由观测变量z_k更新估计，

再更新误差协方差，

P_k＝(1-K_kH)P_k ^-

8、参考图8，说明鲁棒卡尔曼滤波的基本流程。在故障情况下，卡尔曼滤波器没有较强的鲁棒性，在实际系统中，难免会出现接收机设备故障或者信息传输中断等问题，导致不能正常提供测量信息，此时传统的卡尔曼滤波器准确性就会大大地降低，甚至会发散。因此用于状态估计的卡尔曼滤波算法需使用一种自适应的方法改善其鲁棒性。抑制测量误差的自适应处理过程其实就是比较新息方差矩阵，当系统中出现测量故障时，真实误差就会超过理论误差，在这种情况下，可以通过调整量测噪声协方差对角矩阵R，确保滤波器在传感器发生故障的情况下的鲁棒性。修正过程中最基本的目的就是为测量噪声协方差矩阵R找到合适的修正系数矩阵，以确保新息协方差真实值和理论值相匹配。

在原有卡尔曼滤波的基础上，我们增加了一个由比例因子组成的比例矩阵S_k+1，来调整量测协方差矩阵，让理论值接近真实值，关系式如下：

其中，δ是移动窗口的宽度，上式的左边实际上代表实际的新息协方差，右边代表理论的新息协方差，经调整可以得到：

在某一个传感器发生量测故障，而其他传感器正常工作的情况下，由比例因子组成的比例矩阵对应的项将变得相对更大，将增大传感器在R矩阵的量测噪声的协方差，最终故障测量将不会被滤波器理会。另一方面，比例因子只有在测量故障的情况下，才会影响系统的估算过程，在正常的运作情况下，这个比例矩阵将是一个单位矩阵 S_k+1＝I_z×z，其中z为更新向量的维数。

矩阵S_k+1可以用上式进行计算，但有时对角元素可能为负或者远小于1，因为δ是一个有限的数，且计算过程中会有估算和圆滑误差。矩阵S_k+1是一个对角矩阵，因为只有对角元素才会对调整有影响，此外矩阵S_k+1不能有元素小于1，因为R矩阵不能针对特殊问题快速递减，因此，按照如下公式计算比例矩阵：

根据观测量对当前时刻的状态变量预估值进行修正，得到当前状态的估计值及相应的误差协方差矩阵。

9、参考图9，说明交互多模型的算法思想。在实际工作过程中，农机并非一直定常速行驶，在转弯或者其他状况下时，速度会出现一定的变化，因此卡尔曼滤波中的定常速模型会出现一定的偏差，交互多模型算法的基本思想是针对目标不同的运动状态使用不同的模型匹配，与各个模型对应的滤波器并行工作，通过多个模型的概率加权构建更加符合实际的系统模型，滤波器最终的估计结果为各个单独模型估计的加权和。

本发明在卡尔曼滤波中引入了交互多模型，构造了一种多模型卡尔曼滤波算法，即 IMM-EKF算法，具体流程如下：

(1)计算IMM-EKF算法的滤波初值。初始状态向量

初始协方差矩阵P_i(0|0)和初始模型概率μ_i(0)；

(2)状态估计的交互式作用

IMM-EKF算法状态估计的交互式作用后得到

和P_op(k|k)，具体计算公式如下

其中，

(3)EKF算法

将

和P_oi(k|k)作为k+l时刻第i个模型的输入，进行EKF，得到滤波输出

和P_i(k+1|k+1)。具体步骤如下：

计算状态和协方差一步预测：

P_i(k+1|k)＝F_i(k)P_oi(k|k)F_i(k)^T+G_i(k)Q_i(k)G_i(k)^T

计算量测预测值：

计算新息协方差矩阵：

计算滤波增益矩阵：

计算状态和协方差更新值：

其中，

是非线性方程的雅可比矩阵，如下所示：

(4)模型概率更新

若模型i的滤波新息为v_i(k+1)，相应的新息协方差为S_i(k+1)，并假定服从高斯分布，那么模型i的可能性为

其中，

则模型i的概率更新为：

其中，

(5)滤波交互输出

10、参考图10，说明两种面积测量方法。一是航迹法，此时农机利用回耕法进行作业，即农机先沿农田边界进行作业，再逐次往内作业。在农机沿外边界绕行一周结束后，由北斗/GPS接收机获取农田边界的经纬度数据，结合相应的算法计算出这一周的周长和所围面积；二是当农田形状比较规则，农机按照梭形耕法进行作业时，通过转弯模型，结合对农机刀具的状态分析，利用农机内部控制刀具升降的电路，引出相关标志，判别出农机运行的边界点，将这些边界点进行滤波与估计，得出边界线，再进行面积计算。

11、参考图11，说明本发明车载式与手持式的两用型结构，本发明由可手持式部分和固定基座两部分构成。可手持式部分上有主控制器、GPS模块、NRF24L01无线传输模块、显示屏及SD卡插口，可实现面积测量的基本功能；固定基座上有CAN总线接口、电源接口以及GPRS模块，可通过底部的四个吸盘吸附在农机驾驶舱门上。正常工作时，手持式部分固定在基座上，由基座进行充电，进行面积测量，并由基座采集CAN 总线数据，利用GPRS模块实现与监控中心的数据传输；若在不需要作业的情况下想单独测量地块的面积，则可以直接取下手持式部分，将其当作普通的手持式测亩仪进行使用。这种两用型结构一方面扩展了单纯手持式结构所不具备的采集CAN总线数据信息的功能，另一方面又考虑了必要情况下所需的农机耕作时测量农田面积或农机耕作后测量农田面积，更便于用户使用。

12、说明利用SD卡进行数据存储的原理。在实际过程中，需要将作业位置、作业面积等信息保存下来，便于一些突发情况时的分析，由于数据比较多，需要较大的存储空间，因此本发明采用SD卡作为存储介质。为了对SD卡进行管理，使用了开源的FATFS 文件系统模块，使用FATFS可以方便地对SD卡进行读写操作。同时，在保存数据时是以文件的形式进行保存，文件名是以时间命名的，这有利于对数据进行区分。

13、参考图12，说明利用GPRS实现通信的原理。本发明要把移动终端的定位数据以及面积测量数据实时传输到远程监控中心，采用无线网络传输系统。本系统选用GPRS 网络来完成移动终端的信息传输。SIM800A的无线数据传输是整个系统的通信核心。 SIM800A中内嵌有TCP/IP协议，开发人员可以通过AT指令操作SIM800A，大大简化了开发的过程。在使用GPRS网络进行数据传输过程中，我们的系统采用TCP通信传输协议。在建立移动终端与PC的GPRS连接时，PC机必须有一个固定的公网IP地址，移动终端的GPRS模块通过AT指令完成模块的初始化，建立TCP/IP网络连接，与远程的监控中心的PC机连接起来，最后通过AT指令集发送数据。在移动终端上网络的连接和数据的传输接收都是通过串口向SIM800A发送数据完成的，整个过程中运行正常的话，SIM800A将回传数据“OK"。

14、参考图12，说明基于CAN总线和GPRS技术的农机远程监控系统。本发明拟在测亩仪的基础功能上，附加与监控中心通信的功能。通过CAN总线获取农机运行状态、油量、温度等数据，通过GPRS网络发送至指定IP的监控中心上。监控中心通过获取农机发送来的数据，在监控平台端PC机上，完成监控中心对车辆的监测功能。在实际应用中，遇到紧急情况时，农机驾驶员可以通过专用报警按钮发出警报，该信号通过CAN总线发送给车载终端，车载终端通过GPRS网络向监控中心发出求救信号。同时，监控中心也可以通过网络发送控制指令到车载终端，车载终端再通过CAN总线和农机其他部件进行通信，实现对农机一定权限内的控制，如切断总电源等，以此来保障农机的安全。

Claims (1)

1.一种改进的农机作业面积测量方法，其特征在于，包括以下步骤：

步骤1：坐标转换和面积计算：通过北斗GPS双模接收机模块获取当前农机的经纬度信息，经过高斯克吕格投影方式进行坐标转换，再利用完全离散格林公式对不规则农田面积进行计算，具体为：

在高斯投影变换中，中央经线投影后的直线为X轴，赤道投影后的直线为Y轴，它们的交点为原点，为使转换之后的Y轴坐标不为负值，在Y轴上加上500km，具体变换公式如下：

式中，B-纬度；l＝L-L₀，L为经度，L₀为轴子午线的经度；N-卯酉圈曲率半径；t＝tanB；η＝e′cosB；e′-椭球第二偏心率；X₀-当l＝0时，从赤道起算的子午线弧长；参数计算如下：

e²＝(a²-b²)/a²e′²＝(a²-b²)/b²

X₀＝a(1-e²)(A₀B+A₂sin2B+A₄sin4B+A₆sin6B+A₈sin8B)

其中：

在精细农业实践中，采用3度分带投影，从东经1°30’开始，每3°为一带，将全球划分为120个投影带，利用高斯投影的公式，求得转换后的平面坐标X，Y；

面积测量利用完全离散格林公式，其具体内容如下：

对于一个封闭的区域，获得其边界上n个坐标点，(X₁，Y₁)，(X₂，Y₂)，…(X_n，Y_n)，则其面积由格林公式计算获得：

其中i＝1，2，…n，当i＝n时，X_i+1＝X₁，Y_i+1＝Y₁

周长由以下公式得到：

步骤2：采用差分的方法，消除电离层误差；结合卡尔曼滤波的算法，对经纬度数据进行预测与修正，具体为：

测量面积前，参考站放在固定位置多次接收定位信号，取平均值获得相对参考坐标(B_r，L_r)作为标准值，测量面积时，参考站实时接收定位信号(B_x，L_x)，并与之前获得的相对参考坐标相减计算误差，得出实时误差值(ΔB_x，ΔL_x)，参考站将实时误差值发送给移动站，移动站利用实时误差值对经纬度数据进行修正；

卡尔曼滤波的算法流程如下：

(1)设定状态与协方差初值X₀和P₀；

(2)时间更新；根据上一时刻的状态变量估算值和误差协方差阵估算值，结合线性化的状态方程，更新当前时刻的状态变量和相应误差协方差阵；公式如下：

P_k ^-＝FP_k-1F^T+Q

式中：F为状态传递矩阵，

为k-1时刻的状态量，K为过程噪声增益矩阵，v为过程噪声扰动量，其相应的协方差为Q；

(3)测量状态更新；

首先由以下公式计算卡尔曼增益

K_k＝P_k ^-H^T(HP_k ^-H^T+R)^-1

式中，H为量测矩阵，R为量测噪声协方差矩阵；

接着由观测变量z_k更新估计，

再更新误差协方差，

P_k＝(1-K_kH)P_k ^-

在原有卡尔曼滤波的基础上，增加了一个由比例因子组成的比例矩阵S_k+1，来调整量测协方差矩阵，让理论值接近真实值，关系式如下：

其中，δ是移动窗口的宽度，上式的左边实际上代表实际的新息协方差，右边代表理论的新息协方差，经调整可以得到：

在某一个传感器发生量测故障的情况下，由比例因子组成的比例矩阵对应的项将变得相对更大，将增大传感器在R矩阵的量测噪声的协方差，最终故障测量将不会被滤波器理会；另一方面，比例因子只有在测量故障的情况下，才会影响系统的估算过程，在正常的运作情况下，这个比例矩阵将是一个单位矩阵S_k+1＝I_z×z，其中z为更新向量的维数；

矩阵S_k+1用上式进行计算，δ是一个有限的数，且计算过程中会有估算和圆滑误差；矩阵S_k+1是一个对角矩阵，按照如下公式计算比例矩阵：

根据观测量对当前时刻的状态变量预估值进行修正，得到当前状态的估计值及相应的误差协方差矩阵；

在实际工作过程中，农机并非一直定常速行驶，在转弯时，速度会出现变化，因此卡尔曼滤波中的定常速模型会出现偏差，交互多模型算法是针对目标不同的运动状态使用不同的模型匹配，与各个模型对应的滤波器并行工作，通过多个模型的概率加权构建更加符合实际的系统模型，滤波器最终的估计结果为各个单独模型估计的加权和；

在卡尔曼滤波中引入了交互多模型，构造了一种多模型卡尔曼滤波算法，即IMM-EKF算法，具体流程如下：

(1)计算IMM-EKF算法的滤波初值；初始状态向量

初始协方差矩阵P_i(0|0)和初始模型概率μ_i(0)；

(2)状态估计的交互式作用

IMM-EKF算法状态估计的交互式作用后得到

和P_oi(k|k)，具体计算公式如下

其中，

(3)EKF算法

将

和P_oi(k|k)作为k+l时刻第i个模型的输入，进行EKF，得到滤波输出

和P_i(k+1|k+1)；具体步骤如下：

计算状态和协方差一步预测：

P_i(k+1|k)＝F_i(k)P_oi(k|k)F_i(k)^T+G_i(k)Q_i(k)G_i(k)^T

计算量测预测值：

计算新息协方差矩阵：

计算滤波增益矩阵：

计算状态和协方差更新值：

其中，

是非线性方程的雅可比矩阵，如下所示：

(4)模型概率更新

若模型i的滤波新息为v_i(k+1)，相应的新息协方差为S_i(k+1)，并假定服从高斯分布，那么模型i的可能性为

其中，

则模型i的概率更新为：

其中，

(5)滤波交互输出

步骤3：运用两种面积测量方法；一是航迹法，此时农机利用回耕法进行作业，即农机先沿农田边界作业，再逐次往内作业，在农机沿外边界绕行一周结束后，利用测亩仪计算出这一周的周长和所围面积；二是定幅法，当农田形状比较规则，农机按照梭形耕法进行作业时，默认农机的作业幅宽为固定值，用幅宽乘以农机走过的距离得到面积；

步骤4：数据的存储与传输；通过CAN总线获取农机工况信息，一方面通过本地数据存储设备将数据保存在车上，另一方面通过GPRS网络发送至指定IP的监控中心平台上，实现监控中心对车辆工况的实时监测和对突发状况的处理。

Images