1. An agricultural machinery behavior analysis and operation area statistical method based on Beidou positioning data is characterized by comprising the steps of automatically identifying an agricultural machinery operation area, calculating the area, analyzing the overlapped area and the omitted area;

the method for automatically identifying the agricultural machinery operation area is an agricultural machinery operation area automatic identification algorithm based on spatial clustering;

the method for calculating the area comprises a grid-based area calculation method and a contour-based area calculation method;

the method for analyzing the overlapping area is calculated by subtracting the area obtained by the grid-based area calculation method from the area obtained by multiplying the track length by the width:

wherein d (Q)_i，Q_i+1) Representing the distance between adjacent track points of the agricultural machinery operation;

the method for analyzing the missing area comprises the following steps:

S_miss＝S_contour-S_grid；

wherein the Scontour is a contour-based area; the Sgrid is to calculate the grid-based area.

2. The agricultural machinery behavior analysis and working area statistical method based on Beidou positioning data as set forth in claim 1, wherein the agricultural machinery working area automatic identification algorithm based on spatial clustering comprises:

1-1) agricultural machinery operation data acquisition

Track data set P ═ P of agricultural machinery operation is obtained through vehicle-mounted Beidou positioning terminal and GPRS mobile communication equipment₁(t₁，lat₁，long₁)，P₂(t₂，lat₂，long₂)，…，P_n(t_n，lat_n，long_n) Wherein t represents time, lat represents latitude, long represents longitude, and n represents the total number of track points;

1-2) data preprocessing

The preprocessing refers to removing data abnormal points, drift points, stop points and random noise points;

1-3) projection

Obtaining a data point set Q ═ Q under a UTM coordinate system₁(t₁，x₁，y₁)，Q₂(t₂，x₂，y₂)，…，Q_n(t_n，x_n，y_n)；

1-4) spatial clustering

Identifying the operation area by utilizing spatial clustering, wherein the identification comprises the following specific steps:

1-4-1) drawing a circle by a certain radius r with each preprocessed data point as the center of the circle, wherein the density value of the point is formed by how many adjacent data points in the circle;

1-4-2) if the density value of the point is less than a set threshold value min _ pts, marking the point as a low density point, otherwise, marking the point as a high density point;

1-4-3) connecting two points if a high density point is within the circle of another high density point; if a certain low-density point is in the circle of another high-density point, connecting the low-density point to the high-density point closest to the low-density point to form a boundary point;

1-4-4) repeating the steps 1-4-2) and 1-4-3), eliminating low-density points which are not in the circle of any high-density point, and reserving a high-density point set as track points of an agricultural machinery operation area

1-5) calculating the contour-based area Scontour;

1-6) calculating the grid-based area Sgrid;

1-7) the method of analyzing the analysis overlap area is:

the area obtained by subtracting the area obtained by the grid-based area calculation method from the area obtained by multiplying the track length by the width is calculated as follows:

wherein d (Q)_i，Q_i+1) Representing the distance between adjacent track points of the agricultural machinery operation;

the method for analyzing the missing area comprises the following steps:

S_miss＝S_contour-S_grid。

3. the agricultural machinery behavior analysis and working area statistical method based on Beidou positioning data as set forth in claim 2, wherein the data preprocessing method of step 1-2) comprises, but is not limited to, the sequence of data elimination:

1-2-1) rejecting abnormal points: any agricultural machinery track point should satisfy P (t, lat, long):

lat∈[-90°，90°]

long∈[-180°，180°]

data points which do not meet the formula are taken as abnormal points to be removed;

1-2-2) removing drift points: for 2 adjacent track points P_i(t_i，lat_i，long_i)，P_i+1(t_i+1，lat_i+1，long_i+1) Calculating the running speed of the agricultural machine:

wherein d (P)_i，P_i+1) Representing adjacent track points P_i、P_i+1A distance between, v (P)_iP_i+1)＞v_maxTrace point elimination of (v), wherein_maxThe maximum operating speed of the agricultural machine;

1-2-3) elimination stop point: calculating the average speed of k continuous agricultural machine operation track points:

removing track points with the average speed smaller than a certain threshold value delta;

1-2-4) eliminating random noise points: for 2 adjacent track points P_i(t_i，lat_i，long_i)，P_i+1(t_i+1，lat_i+1，long_i+1) Calculate its direction:

the expression method for converting the vector into the unit vector comprises the following steps: theta_i，i+1→(cos(θ_i，i+1)，sin(θ_i，i+1) Then, the mean direction value of k continuous agricultural machine operation track points is calculated as:

the standard deviation was calculated as:

and eliminating points with standard deviation larger than a certain threshold value.

4. The agricultural machinery behavior analysis and working area statistical method based on Beidou positioning data as set forth in claim 2, wherein 1-5) the method for calculating the contour-based area Scontour comprises the following steps:

1-5-1) foveal bag calculation

And performing concave packet calculation on each type of data points according to the data points obtained by clustering, wherein the specific steps are as follows:

(1) finding out the point with the minimum y value, and taking the point with the maximum x value as a starting point O if the y values are the same;

(2) starting from an initial point O, taking (1, 0) as a reference vector, firstly finding an edge with a radius of R smaller than that of the initial edge, and taking the point as A;

(3) looping to find the next edge, assuming the previous edge is AB, then the next edge must start at point B and connect to a point C in the R neighborhood of B, using the following rule: firstly, the points in the R neighborhood of B are sorted in the polar coordinate direction by taking B as the center and BA vector as the reference, and then the points C in the R neighborhood of B are sequentially sorted₀～C_nSet up with BC_iThe circle is a chord, whether other neighborhood points are included is checked, if the other neighborhood points do not exist, the chord is a new edge, and a cycle is jumped out;

(4) finding all edges in sequence until no new edge can be found or a point which is used as an edge before is encountered;

1-5-2) calculated area

Calculate the value of each classCalculating the area of the polygon by adopting a triangle segmentation algorithm or Simpson algorithm to obtain S_{contour_in}I.e. by

Wherein w represents the working width, d (Q)_i，Q_i+1) Indicating adjacent boundary points Q_i、Q_i+1The distance between them.

5. The agricultural machinery behavior analysis and working area statistical method based on Beidou positioning data as set forth in claim 2, wherein the method of calculating grid-based area Sgrid in step 1-6) is as follows:

1-6-1), wherein the area expansion is to perform rasterization processing on an actual operation area according to agricultural machinery track data and breadth, and the specific steps are as follows:

(a-1) finding the minimum x of x, y_min，y_minAnd a maximum value x_max，y_max；

(a-2) determining the size of the grid matrix to be opened up according to the ratio mu of the pixel elements to the actual size:

where ε represents an additional boundary that ensures that the data points are all located within the matrix; opening up a two-dimensional array representing a grid matrix and initializing to 0;

(a-3) calculating the area to be expanded according to the agricultural machine track and the width:

known adjacent agricultural machinery operation track point Q_i(t_i，x_i，y_i)，Q_i+1(t_i+1，x_i+1，y_i+1) And an operating width w, the region expansion of which is substantially Q'_i，Q″_i，Q′_i+1，Q″_i+1Coordinates of four points:

Q′_i(x′_i，y′_i)＝(x_i+Δ_x，y_i-Δ_y)

Q″_i(x″_i，y″_i)＝(x_i-Δ_x，y_i+Δ_y)

Q′_i+1(x′_i+1，y′_i+1)＝(x_i+1+Δ_x，y_i+1-Δ_y)

Q″_i+1(x″_i+1，y″_i+1)＝(x_i+1-Δ_x，y_i+1+Δ_y)

according to Q'_i，Q″_i，Q′_i+1，Q″_i+1Rectangles and Q 'generated from these four points'_i+1，Q″_i+1Generated by the two pointsSuperposing the semi-circle with the radius and the grid matrix: obtaining a rasterized agricultural machinery operation track diagram

1-6-2) calculated area

According to the grid matrix, counting the number of grid unit values of 1, and then calculating the area according to the proportion of the pixel elements to the actual size:

S_grid＝N*w²

where N is the number of grid cell values 1.

Background

Most of the traditional agriculture in China is still completed in a manual intervention mode, the automation degree is not high, for example, in the aspects of statistics and calculation of the operating area of the agricultural machine, most of the traditional agriculture is still realized in a manual division and tape measure mode, a large amount of manpower and time are consumed, certain errors exist in manual measurement, the manual measurement is a static measurement mode, and the operating area of the agricultural machine cannot be calculated in real time.

The precision agriculture is the most important part in the development of the science and technology of the agriculture in the twenty-first century, has high science and technology content and strong integration comprehensiveness, greatly improves the agricultural production efficiency, and becomes an important target of modern agricultural production management. The Beidou positioning and navigation system is one of core supporting technologies, and plays an irreplaceable role in the aspects of fine agricultural seeding and harvesting, real-time agricultural machine positioning and monitoring, automatic agricultural machine operation navigation, remote command and scheduling and the like. With the popularization of large-scale, intensive and regional agricultural production processes, timely mastering of agricultural machine operation areas and task completion conditions is of great significance to overall efficiency assessment and subsequent operation scheduling. Most of traditional agricultural machinery operation area statistical methods are realized by manually reporting by agricultural machinery operating personnel or entrusting a third party to carry out on-site surveying and mapping, and relate to the problems of more human factors, larger errors and consumption of a large amount of manpower, material resources and time. In view of the above problems, the technical problems proposed by those skilled in the art are: how to use dynamic metering to achieve automatic and accurate measurement of the working area of agricultural machinery?

At present, in the aspect of dynamic measurement of the working area of agricultural machinery, wheel revolution measurement methods, ultrasonic measurement methods, laser measurement methods and measurement methods based on a positioning navigation system (such as GPS and Beidou) are mainly used. The wheel revolution measuring method mainly utilizes a speed sensor to count the revolution of wheels of the agricultural machine, the operation distance of the agricultural machine is indirectly calculated according to the radius of the wheels, and then the operation area is obtained by multiplying the operation width of the agricultural machine by the operation distance. The method can effectively measure the area of the agricultural machinery during overlapping operation, but needs crops as reference to obtain the actual operation width, and is mostly used for calculating the operation area during harvesting. The measurement method based on positioning and navigation mainly comprises a boundary-based measurement method and a track-based measurement method. Most boundary-based measurement methods are that the agricultural machinery operation area is wound for a circle, and then a triangle segmentation algorithm or a Simpson algorithm is adopted to realize area calculation.

The boundary-based measurement method can be used for a working area with an arbitrary shape, and the larger the area is, the higher the calculation accuracy is.

The track-based measuring method mainly comprises a width method, a buffer area vector method and other calculation methods. The breadth method is to obtain the area of an operation area according to the multiplication of the length of an operation track of the agricultural machine and the operation breadth, the operation length is regarded as the sum of the distances between two adjacent data points, the method can dynamically calculate the operation area of the agricultural machine in real time, but the operation area cannot be accurately calculated when the overlapped operation area exists, and the method is mainly used for calculating the area when the agricultural machine with accurate autonomous navigation performs full-breadth operation. The vector method of the buffer area is mainly to construct the buffer area of the track line entity, which essentially translates the distance of half the operation width along the vertical direction to two sides of the track (assuming that the positioning terminal is on the central axis of the agricultural machinery) to obtain two parallel lines according to the running track of the agricultural machinery, and fits the two ends or one end of the parallel lines by adopting a smooth curve, and the finally obtained closed area is the buffer area of the line entity.

The buffer area vector method can be used for operation areas in any shapes, effective operation and missing operation area can be counted, but the design has high calculation complexity and relates to multiple intersection calculation.

In order to solve the above technical problems, the chinese patent literature discloses the following technical contents:

chinese patent document CN107036572B discloses a method and a device for obtaining the working area of an agricultural machine, which includes: receiving agricultural machinery operation track data sent by an agricultural machinery positioning device; improving a neighborhood radius determination method of a dbscan clustering algorithm based on the operation speed; filtering road driving points and field transition points in the agricultural machinery operation track data by adopting an improved dbscan clustering algorithm; determining the number of agricultural machinery operation field pieces according to the filtered agricultural machinery operation track data; and respectively calculating the area of each agricultural machine operation field by using a distance method. However, the method described in the patent document is not suitable for the acquisition requirement of the existing agricultural machinery operation track data: the agricultural machinery operation track data is required to contain speed and course information, so that the distance method is not suitable for area statistical analysis during overlapping operation.

Chinese patent document CN107462208A discloses an agricultural machine and an agricultural machine working area measuring device and measuring method, which collects longitude and latitude data in real time by the working area measuring device on the agricultural machine, and records the running track of the agricultural machine; removing the longitude and latitude information points with the deviation, and determining the operation track of the agricultural machine; determining the contour line of the agricultural machinery operation track graph as an operation area boundary line; the work area is calculated based on the plurality of measurement reference points and the work area boundary line. The method accurately determines the boundary line of the operation area by removing the offset point and supplementing the blank point, further calculates the area of the operation area by the multiple measurement reference points, and has the advantages of high measurement speed and high precision. The technique of this document essentially uses a boundary-based area calculation method, and thus cannot analyze and count the area of the missing work area.

Chinese patent document CN107843228B discloses a method for acquiring a spatial trajectory area of a multi-layer scanning timing sequence, which includes: performing Gauss-Krueger projection on a running track point on an agricultural machinery operation track; acquiring a first external rectangle of the coordinates of the operation track points; respectively generating line buffer areas aiming at the coordinates of every two adjacent running track points; scanning and rasterizing each line buffer area, and calculating the sum of the grid areas covered by each line buffer area to obtain a first operation area; scanning each line buffer area again, regrooving the grids which are not completely covered in each line buffer area, and calculating the sum of the areas of the grids covered by each line buffer area to obtain a second operation area; and when the absolute value of the difference value between the second working area and the first working area is smaller than a set error threshold value, taking the second working area as the actual working area of the agricultural machine. The method for calculating the area in the patent document adopts a mode of combining a buffer vector method and a grid method, and relates to twice measurement, the calculation complexity is high, the operation area cannot be automatically identified, and the area of the omitted operation area cannot be analyzed and counted.

In conclusion, the invention utilizes the Beidou positioning terminal to realize real-time acquisition of the travel track of the agricultural machine, and further provides a method for analyzing the behavior of the agricultural machine, which can automatically identify the operation area of the agricultural machine only by depending on the longitude and latitude information, the time and the operation width of the operation of the agricultural machine, and analyze and count the areas of the agricultural machine during effective operation, missed operation and overlapped operation.

Disclosure of Invention

Aiming at the defects of the prior art, the invention discloses an agricultural machinery behavior analysis and operation area statistical method based on Beidou positioning data.

The invention aims to: in order to improve the accuracy and efficiency of the area statistics of the agricultural machinery operation area, reduce the investment of manpower, material resources and time and adapt to the requirements of modern agricultural development, the method for automatically analyzing the agricultural machinery behaviors and counting the area of the operation area through Beidou positioning data is provided, each subarea of the agricultural machinery operation can be automatically identified, and the method is suitable for area statistics under the conditions of overlapping operation and missing operation.

The detailed technical scheme of the invention is as follows:

an agricultural machinery behavior analysis and operation area statistical method based on Beidou positioning data is characterized by comprising the steps of automatically identifying an agricultural machinery operation area, calculating the area, analyzing the overlapped area and the omitted area;

the method for automatically identifying the agricultural machinery operation area is an agricultural machinery operation area automatic identification algorithm based on spatial clustering;

the method for calculating the area comprises a grid-based area calculation method and a contour-based area calculation method;

the method for analyzing the overlapping area is calculated by subtracting the area obtained by the grid-based area calculation method from the area obtained by multiplying the track length by the width:

wherein d (Q)_i，Q_i+1) Representing the distance between adjacent track points of the agricultural machinery operation;

the method for analyzing the missing area comprises the following steps:

S_miss＝S_contour-S_grid；

wherein the Scontour is a contour-based area; the Sgrid is to calculate the grid-based area.

According to the invention, the agricultural machinery operating area automatic identification algorithm based on spatial clustering preferably comprises the following steps:

1-1) agricultural machinery operation data acquisition

Move through on-vehicle big dipper positioning terminal and GPRSThe mobile communication equipment obtains a track data set P ═ { P ═ P of agricultural machinery operation₁(t₁，lat₁，long₁)，P₂(t₂，lat₂，long₂)，…，P_n(t_n，lat_n，long_n) Wherein t represents time, lat represents latitude, long represents longitude, and n represents the total number of track points;

1-2) data preprocessing

The preprocessing refers to removing data abnormal points, drift points, stop points and random noise points;

1-3) projection

Obtaining a data point set Q ═ Q under a UTM coordinate system₁(t₁，x₁，y₁)，Q₂(t₂，x₂，y₂)，…，Q_n(t_n，x_n，y_n)；

1-4) spatial clustering

Identifying the operation area by utilizing spatial clustering, wherein the identification comprises the following specific steps:

1-4-1) drawing a circle by a certain radius r with each preprocessed data point as the center of the circle, wherein the density value of the point is formed by how many adjacent data points in the circle;

1-4-2) if the density value of the point is less than a set threshold value min _ pts, marking the point as a low density point, otherwise, marking the point as a high density point;

1-4-3) connecting two points if a high density point is within the circle of another high density point; if a certain low-density point is in the circle of another high-density point, connecting the low-density point to the high-density point closest to the low-density point to form a boundary point;

1-4-4) repeating the steps 1-4-2) and 1-4-3), removing low-density points which are not in the circle of any high-density point, and reserving a high-density point set as a track point of an agricultural machinery operation area;

1-5) calculating the contour-based area Scontour;

1-6) calculating the grid-based area Sgrid;

1-7) the method of analyzing the analysis overlap area is:

the area obtained by subtracting the area obtained by the grid-based area calculation method from the area obtained by multiplying the track length by the width is calculated as follows:

wherein d (Q)_i，Q_i+1) Representing the distance between adjacent track points of the agricultural machinery operation;

the method for analyzing the missing area comprises the following steps:

S_miss＝S_contour-S_grid。

according to the invention, the method for preprocessing the data in the step 1-2) comprises, but is not limited to, removing the data in the following order:

1-2-1) rejecting abnormal points: any agricultural machinery track point should satisfy P (t, lat, long):

lat∈[-90°，90°]

long∈[-180°，180°]

data points which do not meet the formula are taken as abnormal points to be removed;

1-2-2) removing drift points: for 2 adjacent track points P_i(t_i，lat_i，long_i)，P_i+1(t_i+1，lat_i+1，long_i+1) Calculating the running speed of the agricultural machine:

wherein d (P)_i，P_i+1) Representing adjacent track points P_i、P_i+1A distance between, v (P)_iP_i+1)＞v_maxTrace point elimination of (v), wherein_maxThe maximum operating speed of the agricultural machine;

1-2-3) elimination stop point: calculating the average speed of k continuous agricultural machine operation track points:

removing track points with the average speed smaller than a certain threshold value delta;

1-2-4) eliminating random noise points: for 2 adjacent track points P_i(t_i，lat_i，long_i)，P_i+1(t_i+1，lat_i+1，long_i+1) Calculate its direction:

the expression method for converting the vector into the unit vector comprises the following steps: theta_i，i+1→(cos(θ_i，i+1)，sin(θ_i，i+1) Then, the mean direction value of k continuous agricultural machine operation track points is calculated as:

the standard deviation was calculated as:

and eliminating points with standard deviation larger than a certain threshold value.

According to a preferred embodiment of the present invention, the method 15) for calculating the contour-based area Scontour includes the following steps:

1-5-1) foveal bag calculation

And performing concave packet calculation on each type of data points according to the data points obtained by clustering, wherein the specific steps are as follows:

(1) finding out the point with the minimum y value, and taking the point with the maximum x value as a starting point O if the y values are the same;

(2) starting from an initial point O, taking (1, 0) as a reference vector, firstly finding an edge with a radius of R smaller than that of the initial edge, and taking the point as A;

(3) looping to find the next edge, assuming the previous edge is AB, then the next edge must start at point B and connect to a point C in the R neighborhood of B, using the following rule: firstly, the points in the R neighborhood of B are sorted in the polar coordinate direction by taking B as the center and BA vector as the reference, and then the points C in the R neighborhood of B are sequentially sorted₀～C_nSet up with BC_iThe circle is a chord, whether other neighborhood points are included is checked, if the other neighborhood points do not exist, the chord is a new edge, and a cycle is jumped out;

(4) finding all edges in sequence until no new edge can be found or a point which is used as an edge before is encountered;

1-5-2) calculated area

Calculating the concave packet of each category to obtain a plurality of polygons, and then calculating the polygon area by adopting a triangle segmentation algorithm or Simpson algorithm to obtain S_{contour_in}I.e. by

Wherein w represents the working width, d (Q)_i，Q_i+1) Indicating adjacent boundary points Q_i、Q_i+1The distance between them.

Preferably, according to the present invention, the method for calculating the grid-based area Sgrid in step 16) is as follows:

1-6-1), wherein the area expansion is to perform rasterization processing on an actual operation area according to agricultural machinery track data and breadth, and the specific steps are as follows:

(a-1) finding the minimum x of x, y_min，y_minAnd a maximum value x_max，y_max；

(a-2) determining the size of the grid matrix to be opened up according to the ratio mu of the pixel elements to the actual size:

where ε represents an additional boundary that ensures that the data points are all located within the matrix; opening up a two-dimensional array representing a grid matrix and initializing to 0;

(a-3) calculating the area to be expanded according to the agricultural machine track and the width:

known adjacent agricultural machinery operation track point Q_i(t_i，x_i，y_i)，Q_i+1(t_i+1，x_i+1，y_i+1) And an operating width w, the region expansion of which is substantially Q'_i，Q″_i，Q′_i+1，Q″_i+1Coordinates of four points:

Q′_i(x′_i，y′_i)＝(x_i+Δ_x，y_i-Δ_y)

Q″_i(x″_i，y″_i)＝(x_i-Δ_x，y_i+Δ_y)

Q′_i+1(x′_i+1，y′_i+1)＝(x_i+1+Δ_x，y_i+1-Δ_y)

Q″_i+1(x″_i+1，y″_i+1)＝(x_i+1-Δ_x，y_i+1+Δ_y)

according to Q'_i，Q″_i，Q′_i+1，Q″_i+1Rectangles and Q 'generated from these four points'_i+1，Q″_i+1Generated by the two pointsSuperposing the semi-circle with the radius and the grid matrix: obtaining a rasterized agricultural machinery operation track diagram;

1-6-2) calculated area

According to the grid matrix, counting the number of grid unit values of 1, and then calculating the area according to the proportion of the pixel elements to the actual size:

S_grid＝N*w²

where N is the number of grid cell values 1.

The technical advantages of the invention are as follows:

1. the agricultural machinery operation area can be automatically identified, the investment of manpower, material resources and time is reduced, and the identification result is shown in the attached drawing.

2. The invention can carry out real-time statistical analysis on the effective farming area, the area of the lost and lost area and the overlapping farming area of the agricultural machinery operation area, and provides a foundation for accurate agricultural analysis.

3. The invention is suitable for the positioning terminal with low cost and is insensitive to the noise and drift of the positioning data.

4. The agricultural machinery operation area can be effectively analyzed only by positioning data and operation width of the agricultural machinery all day long, and other hardware is not needed to inform the algorithm model whether the agricultural machinery is in an operation state or not.

5. The method can not only count the total effective operation area of the agricultural machinery all day, but also effectively count the area of each sub-operation area.

Drawings

FIG. 1 is a flow chart of the scheme of the present invention;

FIG. 2 is a schematic view of an agricultural machinery supervision platform;

FIG. 3 is a schematic view of a travel track of an agricultural machine;

FIG. 4 is a peripheral profile of an agricultural machine travel track;

FIG. 5 is a schematic illustration of zone expansion;

FIG. 6 is a schematic diagram of a rasterized agricultural machine trajectory;

FIG. 7 is a schematic illustration of finding a new edge during the computation of a notch;

FIG. 8 is a schematic coordinate diagram of a boundary-based work area calculation method;

fig. 9a to 9f are images of agricultural machinery working areas obtained after eliminating interference of overlapping areas and missing areas in the embodiment of the invention.

Detailed Description

The invention is described in detail below with reference to the following examples and the accompanying drawings of the specification, but is not limited thereto.

Examples

As shown in fig. 1, an agricultural machinery behavior analysis and operation area statistical method based on Beidou positioning data comprises automatic identification of an agricultural machinery operation area, area calculation, overlapped area analysis and missing area analysis;

the method for automatically identifying the agricultural machinery operation area is an agricultural machinery operation area automatic identification algorithm based on spatial clustering;

the method for calculating the area comprises a grid-based area calculation method and a contour-based area calculation method;

the method for analyzing the overlapping area is calculated by subtracting the area obtained by the grid-based area calculation method from the area obtained by multiplying the track length by the width:

wherein d (Q)_i，Q_i+1) Representing the distance between adjacent track points of the agricultural machinery operation;

the method for analyzing the missing area comprises the following steps:

S_miss＝S_contour-S_grid；

wherein the Scontour is a contour-based area; the Sgrid is to calculate the grid-based area.

The agricultural machinery operation area automatic identification algorithm based on spatial clustering comprises the following steps:

1-1) agricultural machinery operation data acquisition

As shown in fig. 2, through the vehicle-mounted big dipperThe positioning terminal and the GPRS mobile communication equipment obtain a track data set P ═ { P ] of agricultural machinery operation₁(t₁，lat₁，long₁)，P₂(t₂，lat₂，long₂)，…，P_n(t_n，lat_n，long_n) Wherein t represents time, lat represents latitude, long represents longitude, and n represents the total number of track points;

1-2) data preprocessing

The preprocessing refers to removing data abnormal points, drift points, stop points and random noise points;

1-3) projection

In order to facilitate the calculation of the subsequent area and the distance between two points, the longitude and latitude information (WGS84 coordinate system) needs to be converted into data points in a plane rectangular coordinate system (UTM coordinate system), and a data point set Q ═ Q in the UTM coordinate system is obtained₁(t₁，x₁，y₁)，Q₂(t₂，x₂，y₂)，…，Q_n(t_n，x_n，y_n)；

1-4) spatial clustering

According to the division of the agricultural machinery driving data points, the road driving points and the field transfer driving points are different from the actual operation points in the spatial distribution: as shown in fig. 3 and 4, the track points in the working area are distributed more densely, while the track points in the road driving and field transfer areas are distributed sparsely, and the working area is identified by using spatial clustering, wherein the identification comprises the following specific steps:

1-4-1) drawing a circle by a certain radius r with each preprocessed data point as the center of the circle, wherein the density value of the point is formed by how many adjacent data points in the circle;

1-4-2) if the density value of the point is less than a set threshold value min _ pts, marking the point as a low density point, otherwise, marking the point as a high density point;

1-4-3) connecting two points if a high density point is within the circle of another high density point; if a certain low-density point is in the circle of another high-density point, connecting the low-density point to the high-density point closest to the low-density point to form a boundary point;

1-4-4) repeating the steps 1-4-2) and 1-4-3), eliminating low-density points which are not in the circle of any high-density point, reserving a high-density point set as track points of an agricultural machine operation area, wherein all the points which can be connected together form a class, and the low-density points which are not in the circle of any high-density point are abnormal points (namely road driving points or field transfer points) and can be eliminated, so that the high-density points (namely the track points of the agricultural machine operation area) are obtained;

1-5) calculating the contour-based area Scontour;

1-6) calculating the grid-based area Sgrid;

1-7) the method of analyzing the analysis overlap area is:

the area obtained by subtracting the area obtained by the grid-based area calculation method from the area obtained by multiplying the track length by the width is calculated as follows:

wherein d (Q)_i，Q_i+1) Representing the distance between adjacent track points of the agricultural machinery operation;

the method for analyzing the missing area comprises the following steps:

since the contour-based area calculation method may take into account the missing operation region, the area of the missing region is:

S_miss＝S_contour-S_grid。

the data preprocessing method of the step 1-2) comprises, but is not limited to, the data removing sequence:

1-2-1) rejecting abnormal points: any agricultural machinery track point should satisfy P (t, lat, long):

lat∈[-90°，90°]

long∈[-180°，180°]

data points which do not meet the formula are taken as abnormal points to be removed;

1-2-2) removing drift points: for 2 adjacent track points P_i(t_i，lat_i，long_i)，P_i+1(t_i+1，lat_i+1，long_i+1) Calculating the running speed of the agricultural machine:

wherein d (P)_i，P_i+1) Representing adjacent track points P_i、P_i+1A distance between, v (P)_iP_i+1)＞v_maxTrace point elimination of (v), wherein_maxThe maximum operating speed of the agricultural machine;

1-2-3) elimination stop point: because the agricultural machinery is when the stall condition, the big dipper positioning terminal is the upload data that still does not stop, so need eliminate the stall point, when the agricultural machinery is stopped, the velocity of motion lasts for 0 in certain time, but when the agricultural machinery turns, speed also can be close to 0, so adopt average speed when eliminating the stall point, calculate its average speed to k continuous agricultural machinery operation track points:

removing track points with the average speed smaller than a certain threshold value delta (a small value can be preset by a user according to the actual working condition of the agricultural machine);

1-2-4) eliminating random noise points: because the Beidou positioning and navigation system inevitably suffers from various interferences and generates certain random noise, particularly when the agricultural machinery is in a stopped state, the data points of the agricultural machinery are not always fixed at one position but randomly walk around a certain central point, the noise needs to be eliminated, the influence of the noise on a subsequent clustering algorithm is reduced, and the random noise is characterized in that the direction of the data points changes randomly within a certain time, namely the variance is large; for 2 adjacent track points P_i(t_i，lat_i，long_i)，P_i+1(t_i+1，lat_i+1，long_i+1) Calculate its direction:

the expression method for converting the vector into the unit vector comprises the following steps: theta_i，i+1→(cos(θ_i，i+1)，sin(θ_i，i+1) Then, the mean direction value of k continuous agricultural machine operation track points is calculated as:

the standard deviation was calculated as:

and eliminating points with standard deviation larger than a certain threshold (preset by a user according to the actual working condition of the agricultural machine).

According to the present invention, preferably, the method 1-5) of calculating the contour-based area Scontour can obtain the actual working area of the agricultural machinery through spatial clustering, because there may be a plurality of working areas (i.e. the clustering result is a plurality of categories), the area calculation is performed on each working area according to the category, and the contour-based area calculation method includes the following steps:

1-5-1) foveal bag calculation

Performing the dip-bag calculation on each type of data points according to the data points obtained by clustering, as shown in fig. 4, the specific steps are as follows:

(1) finding out the point with the minimum y value, and taking the point with the maximum x value as the starting point O if the y values are the same, wherein the x and y points are coordinate values of a data point and the point is fixed on the concave bag;

(2) starting from an initial point O, taking (1, 0) as a reference vector, firstly finding an edge with a radius of R smaller than that of the initial edge, and taking the point as A;

(3) the next edge is found in a loop, and if the previous edge is AB, the next edge must start from B and connect to B in the R neighborhoodFor point C, the following rule is used to find point C: firstly, the points in the R neighborhood of B are sorted in the polar coordinate direction by taking B as the center and BA vector as the reference, and then the points C in the R neighborhood of B are sequentially sorted₀～C_nSet up with BC_iThe circle is a chord circle, whether other neighborhood points are included in the circle is checked, if the circle does not exist, the chord is a new edge, a loop is formed, fig. 7 shows how to find a point C, namely, a circle with BC as the chord is established, and then whether other neighborhood points are included is judged, and the point is found if the other neighborhood points do not exist;

(4) finding all edges in sequence until no new edge can be found or a point which is used as an edge before is encountered;

1-5-2) calculated area

Calculating the concave packet of each category to obtain a plurality of polygons similar to those shown in FIG. 8, and then calculating the polygon area by using a triangle segmentation algorithm or Simpson algorithm to obtain S_{contour_in}If the Beidou positioning terminal is installed on the central axis of the agricultural machinery, the actual area also comprises the area of multiplying the length of the peripheral outline by the width/2, namely

Wherein w represents the working width, d (Q)_i，Q_i+1) Indicating adjacent boundary points Q_i、Q_i+1The distance between them.

Preferably, according to the present invention, the method for calculating the grid-based area Sgrid in step 16) is as follows:

the grid data structure is array data composed of pixels (grid units) which are equal in size, uniform in distribution and closely connected, can be used for representing the distribution of spatial ground objects or phenomena, can be used for calculating the operation area of an agricultural machine, and is easier to process the overlapped operation area, and the specific steps are as follows:

1-6-1), wherein the area expansion is to perform rasterization processing on an actual operation area according to agricultural machinery track data and breadth, and the specific steps are as follows:

(a-1) finding the minimum x of x, y_min，y_minAnd a maximum value x_max，y_max；

(a-2) determining the size of the grid matrix to be opened up according to the ratio mu of the pixel elements to the actual size:

where ε represents an additional boundary that ensures that the data points are all located within the matrix; opening up a two-dimensional array representing a grid matrix and initializing to 0;

(a-3) calculating the area to be expanded according to the agricultural machine track and the width:

as shown in figure 5, the running track points Q of the known adjacent agricultural machinery_i(t_i，x_i，y_i)，Q_i+1(t_i+1，x_i+1，y_i+1) And an operating width w, the region expansion of which is substantially Q'_i，Q″_i，Q′_i+1，Q″_i+1Coordinates of four points:

Q′_i(x′_i，y′_i)＝(x_i+Δ_x，y_i-Δ_y)

Q″_i(x″_i，y″_i)＝(x_i-Δ_x，y_i+Δ_y)

Q′_i+1(x′_i+1，y′_i+1)＝(x_i+1+Δ_x，y_i+1-Δ_y)

Q″_i+1(x″_i+1，y″_i+1)＝(x_i+1-Δ_x，y_i+1+Δ_y)

according to Q'_i，Q″_i，Q′_i+1，Q″_i+1Rectangles and Q 'generated from these four points'_i+1，Q″_i+1Generated by the two pointsThe half circle of radius (for the smoothing of the working area) is superimposed with the grid matrix: modifying the value of the grid unit in the rectangular and semicircular range to 1, if the value is 1 (overlapping operation), obtaining a rasterized agricultural machinery operation track diagram without modification, as shown in FIG. 6;

1-6-2) calculated area

According to the grid matrix, counting the number of grid unit values of 1, and then calculating the area according to the proportion of the pixel elements to the actual size:

S_grid＝N*w²

where N is the number of grid cell values 1.

Application of comparative example,

According to the agricultural machinery behavior analysis and operation area statistical method based on Beidou positioning data, the method is applied to 6 different plots and used for detecting the operation tracks of agricultural machinery during operation of the different plots, and the method is formed as shown in fig. 9a to 9 f.

In order to verify the accuracy of the analysis and operation area statistical method, the actual operation area of the 6 plots is artificially measured to obtain the actual operation area.

Meanwhile, a boundary method and a breadth method in the prior art are introduced to process the agricultural machinery operation data respectively to obtain corresponding area data respectively.

The working areas of the agricultural machinery obtained by the method of the invention and the boundary method and the breadth method are compared, and are shown in table 1:

table 1:

wherein "195001018747 _0, 195001018747_1, 195001018747_2, 195001018747_ 3" refers to three identification areas in a parcel 195001018747;

"195001018965 _0, 195001018965_1, 195001018965_2, 195001018965_3, 195001018965_4, 195001018965_ 5" refers to six identification areas in the parcel 195001018965;

"206003470996 _0, 206003470996_ 1" refers to two identification regions in parcel 206003470996.

As can be seen from table 1, the "boundary method" has a large error in area statistics when there is a missing operation; the 'breadth method' takes overlapping operation into consideration, and has a large error; according to the area calculation method based on the contour and the grid, the overlapping area and the missing area can be respectively counted, then the interference of the overlapping area and the missing area is eliminated, and the obtained actual effective area error is smaller than the error of the comparison algorithm.