CN113409240B - A method of agricultural machinery behavior analysis and operation area statistics based on Beidou positioning data - Google Patents

Abstract

A method of agricultural machinery behavior analysis and operation area statistics based on Beidou positioning data, including automatic identification of agricultural machinery operation area, calculation area, analysis of overlapping area and missing area; the automatic identification method of agricultural machinery operation area is based on spatial clustering. An automatic identification algorithm for agricultural machinery operation area; the method for calculating the area includes a grid-based area calculation method and a contour-based area calculation method. In order to improve the accuracy and efficiency of the area statistics of the agricultural machinery operation area, reduce the input of manpower, material resources and time, and meet the needs of modern agricultural development, the invention provides a method for automatically analyzing the agricultural machinery behavior and counting the area of the operation area through the Beidou positioning data, It can automatically identify each sub-area of agricultural machinery operations, and can be applied to area statistics when overlapping operations and missing operations exist at the same time.

Description

A method of agricultural machinery behavior analysis and operation area statistics based on Beidou positioning data

technical field

The invention relates to an agricultural machinery behavior analysis and operation area statistics method based on Beidou positioning data, and belongs to the technical field of agricultural machinery operation data collection and analysis.

Background technique

Most of my country's traditional agriculture is still done through manual intervention, and the degree of automation is not high. For example, in the statistics and calculation of agricultural machinery operation area, most of them are still realized by manual division and tape measurement, which consumes a lot of manpower and time. Moreover, there are certain errors in manual measurement, and they are all static measurement methods, which cannot calculate the working area of agricultural machinery in real time.

Precision agriculture is the top priority of agricultural science and technology development in the 21st century. Its high technology content and strong integration and comprehensiveness greatly improve agricultural production efficiency and become an important goal of modern agricultural production management. As one of the core supporting technologies, Beidou positioning and navigation system has played an irreplaceable role in agricultural fine sowing and harvesting, real-time positioning and monitoring of agricultural machinery, automatic operation and navigation of agricultural machinery, and remote command and dispatch. With the popularization of large-scale, intensive, and regional agricultural production processes, it is of great significance to grasp the operation area of agricultural machinery and the completion of tasks in time for overall efficiency evaluation and subsequent operation scheduling. Most of the traditional agricultural machinery operation area statistical methods are implemented manually by agricultural machinery operators or by entrusting a third-party on-site surveying and mapping, which involves many human factors, large errors, and consumes a lot of manpower, material resources and time. In response to the above problems, the technical problem raised by those skilled in the art is: how to realize automatic and accurate measurement of agricultural machinery operation area by means of dynamic measurement?

At present, in the dynamic measurement of agricultural machinery operation area, there are mainly wheel revolution measurement method, ultrasonic measurement method, laser measurement method and measurement method based on positioning and navigation systems (such as GPS and Beidou). The wheel revolution measurement method mainly uses the speed sensor to count the number of wheel revolutions of the agricultural machinery, indirectly calculates the working distance of the agricultural machinery according to the wheel radius, and then obtains the working area by multiplying the working width of the agricultural machinery by the working distance. This method is simple in principle and low in cost, but agricultural machinery Wheel slippage often occurs in the field, resulting in large errors in area measurement, and it is impossible to accurately measure the operating area when the agricultural machinery has overlapping operations. Ultrasonic and laser measurement methods are similar in principle, but the accuracy of laser measurement method is higher than that of ultrasonic measurement method. It mainly measures the actual working width in real time through ultrasonic or laser sensors, and then combines the speed sensor to calculate the working area. This method can effectively measure The area when the agricultural machinery overlaps, but the actual working width can be obtained only when the crops are used as a reference, which is mostly used for the calculation of the working area during harvesting. The measurement methods based on positioning and navigation mainly include boundary-based measurement methods and trajectory-based measurement methods. Among them, most of the boundary-based measurement methods are to circle around the agricultural machinery operation area, and then use the triangle segmentation algorithm or the Simpson algorithm to calculate the area.

The boundary-based measurement method can be used for any shape of operation area. Generally, the larger the area, the higher the calculation accuracy, but the area of the missing operation area cannot be counted.

The measurement methods based on the trajectory mainly include the calculation methods such as the width method and the buffer vector method. The breadth method is to obtain the area of the operation area according to the length of the operation track of the agricultural machinery multiplied by the operation width. The operation length is regarded as the cumulative sum of the distances between two adjacent data points. This method can dynamically and real-time calculate the operation area of the agricultural machinery. However, when there are overlapping operation areas, the operation area cannot be accurately calculated. This method is mostly used for the area calculation when the agricultural machinery with accurate autonomous navigation performs full-width operation. The buffer vector method is mainly to construct the buffer area of the trajectory line entity, which is essentially based on the running trajectory of the agricultural machinery, and translates the distance of half the working width along the vertical direction on both sides of the trajectory (assuming that the positioning terminal is on the central axis of the agricultural machinery) to get Two parallel lines are fitted with smooth curves at both ends or one end of the parallel lines, and the final closed area is the buffer area of the line entity.

The buffer vector method can be used for any shape of operation area, and it can count the area of effective operation and missing operation area, but the calculation complexity of the design is relatively high, involving multiple and multiple intersection operations.

In view of the above technical problems, Chinese patent documents disclose the following technical contents:

Chinese patent document CN107036572B discloses a method and device for obtaining agricultural machinery operation area, including: receiving agricultural machinery operation trajectory data sent by an agricultural machinery positioning device; improving the neighborhood radius determination method of the dbscan clustering algorithm based on the operation speed; using the improved dbscan clustering algorithm Filter the road driving points and field transition points in the agricultural machinery operation trajectory data; determine the number of agricultural machinery operation fields according to the filtered agricultural machinery operation trajectory data; use the distance method to calculate the area of each agricultural machinery operation field. However, the method described in this patent document is not suitable for the acquisition requirements of the existing agricultural machinery operation trajectory data: the agricultural machinery operation trajectory data needs to include speed and heading information, so the distance method is not suitable for area statistical analysis during overlapping operations.

Chinese patent document CN107462208A discloses an agricultural machinery and an agricultural machinery working area measuring device and a measuring method. The latitude and longitude data is collected in real time through the working area measuring device on the agricultural machinery, and the running track of the agricultural machinery is recorded; the offset latitude and longitude information points are removed to determine the working track of the agricultural machinery. ; Determine the contour line of the agricultural machinery operation track graph as the boundary line of the operation area; calculate the area of the operation area based on a plurality of measurement reference points and the boundary line of the operation area. The invention accurately determines the boundary line of the operation area by removing the offset points and supplementing the blank points, and further calculates the area of the operation area through multiple measurement reference points, and has the advantages of fast measurement speed and high precision. The technology in this document essentially adopts the area calculation method based on the boundary, so it is impossible to analyze and count the area of the missing work area.

Chinese patent document CN107843228B discloses a method for obtaining track area in multi-layer scanning time series, which includes: performing Gauss-Kruger projection on the running track points on the agricultural machinery operation track; obtaining the first circumscribed rectangle of the coordinates of the running track points; The coordinates of the two running track points are generated, and line buffers are generated respectively; each line buffer is scanned and rasterized, and the sum of the grid areas covered by each line buffer is calculated to obtain the first working area; each line buffer is scanned again. There are line buffers, re-rasterize the incompletely covered grids in each line buffer, calculate the sum of the areas of the grids covered by each line buffer, and obtain the second working area; the second working area is the same as When the absolute value of the difference of the first working area is smaller than the set error threshold, the second working area is taken as the actual working area of the agricultural machinery. This patent document calculates the area by a combination of the buffer vector method and the grid method, and involves two measurements. The calculation complexity is high, and the operation area cannot be automatically identified, and the area of the missing operation area cannot be analyzed and counted.

To sum up, the present invention utilizes the Beidou positioning terminal to realize the real-time collection of the driving trajectories of agricultural machinery, and further proposes a method for analyzing the behavior of agricultural machinery, which can automatically identify the operating area of agricultural machinery only by relying on the latitude and longitude information, time and operation width of the operation of the agricultural machinery. The area of effective operation, missing operation, and overlapping operation of agricultural machinery is counted.

SUMMARY OF THE INVENTION

Aiming at the deficiencies of the prior art, the present invention discloses a method for analyzing the behavior of agricultural machinery and statistical operation area based on Beidou positioning data.

The purpose of the present invention is: in order to improve the accuracy and efficiency of the area statistics of agricultural machinery operation areas, reduce the input of manpower, material resources and time, and meet the needs of modern agricultural development, it provides an automatic analysis of agricultural machinery behavior and statistical operation areas through Beidou positioning data. The area method can automatically identify each sub-area of agricultural machinery operations, and can be applied to area statistics when overlapping operations and missing operations exist at the same time.

The detailed technical scheme of the present invention is as follows:

A method of agricultural machinery behavior analysis and operation area statistics based on Beidou positioning data, characterized in that it includes automatic identification of agricultural machinery operation areas, calculation of area, analysis of overlapping area and omission area;

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

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

The method for analyzing the overlapping area is to calculate the area obtained by multiplying the track length by the width minus the area obtained by the grid-based area calculation method:

where d(Q _i , Q _i+1 ) represents the distance between adjacent track points of agricultural machinery operations;

The method for analyzing the missing area is:

S _miss =S _contour -S _grid ;

Wherein, the Scontour is the area based on the contour; the Sgrid is the area based on the calculation grid.

Preferably according to the present invention, the automatic identification algorithm of agricultural machinery operation area based on spatial clustering includes:

1-1) Acquisition of agricultural machinery operation data

The track data set P={P ₁ (t ₁ , lat ₁ , long ₁ ), P ₂ (t ₂ , lat ₂ , long ₂ ), ..., P _n is obtained through the vehicle-mounted Beidou positioning terminal and the GPRS mobile communication device. (t _n , lat _n , long _n ), where t represents time, lat represents latitude, long represents longitude, and n represents the total number of trajectory points;

1-2) Data preprocessing

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

1-3) Projection

A set of data points in the UTM coordinate system Q={Q ₁ (t ₁ , x ₁ , y ₁ ), Q ₂ (t ₂ , x ₂ , y ₂ ), ..., Q _n (t _n , x _n , y ) _n );

1-4) Spatial clustering

Using spatial clustering to identify the work area, the specific steps of identification are as follows:

1-4-1) Take each data point after preprocessing as the center of the circle, draw a circle with a certain radius r, and how many adjacent data points in the circle become the density value of the point;

1-4-2) If the density value of the point is less than the set threshold min_pts, then mark the point as a low density point, otherwise it is a high density point;

1-4-3) If a high-density point is inside the circle of another high-density point, connect the two points; if a low-density point is inside the circle of another high-density point, also connect it to On the nearest high-density point, it becomes the boundary point;

1-4-4) Repeat steps 1-4-2) and 1-4-3), remove the low-density points that are not in the circle of any high-density points, and retain the high-density point set as the trajectory points of the agricultural machinery operation area;

1-5) Calculate the area Scontour based on the contour;

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

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

Calculated as the area obtained by multiplying the track length by the width minus the area obtained by the grid-based area calculation method:

where d(Q _i , Q _i+1 ) represents the distance between adjacent track points of agricultural machinery operations;

The method for analyzing the missing area is:

S _miss = S _contour - S _grid .

Preferably according to the present invention, the method of step 1-2) data preprocessing includes, but does not limit the order of excluding data:

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

lat∈[-90°, 90°]

long∈[-180°, 180°]

Remove data points that do not satisfy the above formula as outliers;

1-2-2) Eliminate drift points: for two adjacent track points P _i (t _i , lat _i , long _i ), P _i+1 (t _i+1 , lat _i+1 , long _i+1 ) to calculate the operating speed of the agricultural machinery:

where d(P _i , P _i+1 ) represents the distance between adjacent trajectory points P _i and P _i+1 , and the trajectory points with v(P _i P _i+1 )>v _max are eliminated, where v _max is The maximum operating speed of the agricultural machinery;

1-2-3) Eliminate the stopping point: Calculate the average speed for k consecutive agricultural machinery running track points:

Eliminate trajectory points whose average speed is less than a certain threshold δ;

1-2-4) Eliminate random noise points: for two adjacent trajectory points P _i (t _i , lat _i , long _i ), P _i+1 (t _i+1 , lat _i+1 , long _{i+ 1} ) Calculate its orientation:

Convert it to the representation method of unit vector: θ _{i, i+1} → (cos(θ _{i, i+1} ), sin(θ _{i, i+1} )), then for k consecutive agricultural machinery running trajectory points to calculate Its direction mean is:

Calculate its standard deviation as:

Eliminate points with standard deviation greater than a certain threshold.

Preferably according to the present invention, described 15) calculate the method for the area Scontour based on the outline, the area calculation method based on the outline comprises the following steps:

1-5-1) Concave Packet Calculation

According to the data points obtained by clustering, the concave hull calculation is performed for each type of data point, and the specific steps are as follows:

(1) Find the point with the smallest y value. If the y value is the same, take the point with the largest x value as the starting point O;

(2) Starting from the starting point O, with (1, 0) as the reference vector, first find an edge smaller than the radius R as the initial edge, and this point is A at this time;

(3) Loop to find the next edge. Assuming that the previous edge is AB, then the next edge must start from point B and connect to a point C in the R neighborhood of B. Use the following rules to find point C: The points in the R neighborhood are sorted in the polar coordinate direction with B as the center and the BA vector as the reference, and then a circle with BC _i as the chord is established for the points C ₀ ~ C _n in the R neighborhood of B in turn, and then checks whether it contains If other neighbor points do not exist, the chord is a new edge and jumps out of the loop;

(4) Find all the edges in turn until no new edge is found or a point that has been used as an edge before is encountered;

1-5-2) Calculate the area

Calculate the concave hull of each category to obtain multiple polygons, and then use the triangle segmentation algorithm or the Simpson algorithm to calculate the polygon area to obtain S _{contour_in} , that is

where w represents the working width, and d(Q _i , Qi ₊₁ ) _represents the distance between adjacent boundary points Qi and Qi ₊₁ .

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

1-6-1) Regional expansion, regional expansion is to rasterize the actual operation area according to the agricultural machinery track data and width. The specific steps are as follows:

(a-1) Find the minimum value x _min , y _min and the maximum value x _max , y _max of x, y;

(a-2) Determine the size of the grid matrix to be developed according to the ratio μ of the pixel to the actual size:

Where ε represents an additional boundary to ensure that the data points are all located in the matrix; open up a two-dimensional array representing the grid matrix and initialize it to 0;

(a-3) Calculate the area to be expanded according to the trajectory and width of the agricultural machinery:

Knowing the adjacent agricultural machinery running track points Q _i (t _i , x _i , y _i ), Q _i+1 (t _i+1 , x _i+1 , y _i+1 ) and the working width w, carry out Area expansion is actually to find the coordinates of four points Q′ _i , Q″ _i , Q′ _i+1 , Q″ _i+1 :

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 the rectangle generated by the four points Q′ _i , Q″ _i , Q′ _i+1 , Q″ _i+1 and the rectangle generated by the two points Q′ _i+1 , Q″ _i+1

Superimpose the semicircle with the radius and the grid matrix: get the rasterized agricultural machinery running track map;

1-6-2) Calculate the area

According to the above grid matrix, count the number of grid cells with a value of 1, and then calculate the area according to the ratio of the cell to the actual size:

S _grid =N*w ²

where N is the number of grid cells with a value of 1.

The technical advantages of the present invention are:

1. The present invention can automatically identify the agricultural machinery operation area, reducing the input of manpower, material resources and time, and the identification results are shown in the attached drawings.

2. The present invention can make real-time statistical analysis of the effective farming area, the missing area and the overlapping farming area of the agricultural machinery operation area, which provides a basis for precision agricultural analysis.

3. The present invention is applicable to low-cost positioning terminals, and is insensitive to noise and drift in positioning data.

4. The present invention only needs the all-day positioning data and the working width of the agricultural machinery to effectively analyze the working area of the agricultural machinery, and does not need other hardware to inform the algorithm model whether the agricultural machinery is in the working state.

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

Description of drawings

Fig. 1 is the flow chart of the scheme of the present invention;

Figure 2 is a schematic diagram of the agricultural machinery supervision platform;

Fig. 3 is the schematic diagram of the driving track of agricultural machinery;

Figure 4 is the outer contour of the driving track of the agricultural machinery;

Fig. 5 is the schematic diagram of area expansion;

Fig. 6 is the schematic diagram of the running track of grid agricultural machinery;

Fig. 7 is the principle explanatory diagram of finding new edge during concave hull calculation;

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

9a-9f are images of the agricultural machinery operation area obtained after eliminating the interference of the overlapping area and the missing area in the embodiment of the present invention.

Detailed ways

The present invention will be described in detail below with reference to the embodiments and the accompanying drawings, but is not limited thereto.

example,

As shown in Figure 1, a method of agricultural machinery behavior analysis and operation area statistics based on Beidou positioning data, including automatic identification of agricultural machinery operation area, calculation area, analysis of overlapping area and omission area;

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

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

The method for analyzing the overlapping area is to calculate the area obtained by multiplying the track length by the width minus the area obtained by the grid-based area calculation method:

where d(Q _i , Q _i+1 ) represents the distance between adjacent track points of agricultural machinery operations;

The method for analyzing the missing area is:

S _miss =S _contour -S _grid ;

Wherein, the Scontour is the area based on the contour; the Sgrid is the area based on the calculation grid.

The automatic identification algorithm of agricultural machinery operation area based on spatial clustering includes:

1-1) Acquisition of agricultural machinery operation data

As shown in Figure 2, the track data set P={P ₁ (t ₁ , lat ₁ , long ₁ ), P ₂ (t ₂ , lat ₂ , long ₂ ) is obtained through the vehicle-mounted Beidou positioning terminal and the GPRS mobile communication device. ), ..., P _n (t _n , lat _n , long _n ), where 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 eliminating 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, it is necessary to convert the latitude and longitude information (WGS84 coordinate system) into data points in the plane Cartesian coordinate system (UTM coordinate system), and obtain the set of data points in the UTM coordinate system Q={ Q ₁ (t ₁ , x ₁ , y ₁ ), Q ₂ (t ₂ , x ₂ , y ₂ ), ..., Q _n (t _n , x _n , y _n );

1-4) Spatial clustering

According to the division of agricultural machinery driving data points, there are differences in the spatial distribution of road driving points and field transfer driving points and actual operating points: as shown in Figures 3 and 4, the distribution of trajectory points in the operating area is relatively dense, while road driving and The distribution of track points in the field transfer area is relatively sparse, and spatial clustering is used to identify the operation area. The specific steps of identification are as follows:

1-4-1) Take each data point after preprocessing as the center of the circle, draw a circle with a certain radius r, and how many adjacent data points in the circle become the density value of the point;

1-4-2) If the density value of the point is less than the set threshold min_pts, then mark the point as a low density point, otherwise it is a high density point;

1-4-3) If a high-density point is inside the circle of another high-density point, connect the two points; if a low-density point is inside the circle of another high-density point, also connect it to On the nearest high-density point, it becomes the boundary point;

1-4-4) Repeat steps 1-4-2) and 1-4-3), remove the low-density points that are not in the circle of any high-density points, and keep the high-density point set as the trajectory points of the agricultural machinery operation area. The points that can be connected together form a class, and the low-density points that are not in the circle of any high-density points are abnormal points (ie road driving points or field transfer points), which can be eliminated to obtain high-density points (ie agricultural machinery). track points in the work area);

1-5) Calculate the area Scontour based on the contour;

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

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

Calculated as the area obtained by multiplying the track length by the width minus the area obtained by the grid-based area calculation method:

where d(Q _i , Q _i+1 ) represents the distance between adjacent track points of agricultural machinery operations;

The method for analyzing the missing area is:

Since the contour-based area calculation method may take the missed job area into account, the area of the missed area is:

S _miss = S _contour - S _grid .

Described step 1-2) the method for data preprocessing includes, but does not limit the order of eliminating data:

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

lat∈[-90°, 90°]

long∈[-180°, 180°]

Remove data points that do not satisfy the above formula as outliers;

1-2-2) Eliminate drift points: for two adjacent track points P _i (t _i , lat _i , long _i ), P _i+1 (t _i+1 , lat _i+1 , long _i+1 ) to calculate the operating speed of the agricultural machinery:

where d(P _i , P _i+1 ) represents the distance between adjacent trajectory points P _i and P _i+1 , and the trajectory points with v(P _i P _i+1 )>v _max are eliminated, where v _max is The maximum operating speed of the agricultural machinery;

1-2-3) Eliminate the stop point: Since the Beidou positioning terminal still uploads data continuously when the agricultural machine is in the stop state, it is necessary to remove the stop point. When the agricultural machinery turns, the speed will also be close to 0, so the average speed is used when eliminating the stop point, and the average speed is calculated for k consecutive agricultural machinery running trajectory points:

Eliminate the trajectory points whose average speed is less than a certain threshold δ (a very small value, which can be preset by the user according to the actual working conditions of agricultural machinery);

1-2-4) Eliminate random noise points: Because the Beidou positioning and navigation system is inevitably subject to various interferences, certain random noises appear, especially when the agricultural machinery is in a stopped state, the agricultural machinery data points are not always fixed in one position, Instead, it walks randomly around a certain center point, so it is necessary to remove these noises to reduce the impact of noise on subsequent clustering algorithms. The characteristic of random noise is that the direction of data points changes randomly within a certain period of time, that is, the variance is large; For two adjacent trajectory points P _i (t _i , lat _i , long _i ), P _i+1 (t _i+1 , lat _i+1 , long _i+1 ) calculate their directions:

Convert it to the representation method of unit vector: θ _{i, i+1} → (cos(θ _{i, i+1} ), sin(θ _{i, i+1} )), then for k consecutive agricultural machinery running trajectory points to calculate Its direction mean is:

Calculate its standard deviation as:

The points whose standard deviation is greater than a certain threshold (the user can preset according to the actual working conditions of agricultural machinery) are eliminated.

According to the preferred method of the present invention, the method 1-5) calculates the area Scontour based on the contour, the actual operation area of the agricultural machinery can be obtained by spatial clustering, because there may be multiple blocks in the operation area (that is, the clustering results are multiple categories) , so the area is calculated for each work area according to the category. The contour-based area calculation method includes the following steps:

1-5-1) Concave Packet Calculation

According to the data points obtained by clustering, the concave hull calculation is performed for each type of data point, as shown in Figure 4. The specific steps are as follows:

(1) Find the point with the smallest y value. If the y value is the same, take the point with the largest x value as the starting point O, where the x and y points are the coordinate values of the data points, and this point must be on the concave envelope;

(2) Starting from the starting point O, with (1, 0) as the reference vector, first find an edge smaller than the radius R as the initial edge, and this point is A at this time;

(3) Loop to find the next edge. Assuming that the previous edge is AB, then the next edge must start from point B and connect to a point C in the R neighborhood of B. Use the following rules to find point C: The points in the R neighborhood are sorted in the polar coordinate direction with B as the center and the BA vector as the reference, and then a circle with BC _i as the chord is established for the points C ₀ ~ C _n in the R neighborhood of B in turn, and then checks whether it contains If other neighborhood points do not exist, then the chord is a new edge and jumps out of the loop. Figure 7 shows how to find the C point, that is, to establish a circle with BC as the chord, and then determine whether it contains other neighborhood points. exists, that is, the point is found;

(4) Find all the edges in turn until no new edge is found or a point that has been used as an edge before is encountered;

1-5-2) Calculate the area

Calculate the concave hull of each category to obtain multiple polygons similar to those shown in Figure 8, and then use the triangle segmentation algorithm or the Simpson algorithm to calculate the polygon area to obtain S _{contour_in} , assuming that the Beidou positioning terminal is installed on the central axis of the agricultural machinery, then the actual The area also includes the area of the outer contour length multiplied by the width/2, that is

where w represents the working width, and d(Q _i , Qi ₊₁ ) _represents the distance between adjacent boundary points Qi and Qi ₊₁ .

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

The grid data structure is an array data composed of equal-sized, uniformly distributed, and closely connected pixels (grid cells), which can be used to represent the distribution of spatial objects or phenomena. calculation, and it is easier to deal with overlapping work areas. The specific steps are as follows:

1-6-1) Regional expansion, regional expansion is to rasterize the actual operation area according to the agricultural machinery track data and width. The specific steps are as follows:

(a-1) Find the minimum value x _min , y _min and the maximum value x _max , y _max of x, y;

(a-2) Determine the size of the grid matrix to be developed according to the ratio μ of the pixel to the actual size:

Where ε represents an additional boundary to ensure that the data points are located in the matrix; open up a two-dimensional array representing the grid matrix and initialize it to 0;

(a-3) Calculate the area to be expanded according to the trajectory and width of the agricultural machinery:

As shown in Fig. 5, it is known that the adjacent agricultural machinery running track points Q _i (t _i , x _i , y _i ), Q _i+1 (t _i+1 , x _i+1 , y _i+1 ) and the working width The width w, the area expansion of it is actually to find the coordinates of the four points Q′ _i , Q″ _i , Q′ _i+1 , Q″ _i+1 :

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 )

为半径的半圆(为了作业区域的平滑)与栅格矩阵进行叠加处理：把此矩形以及半圆范围内的栅格单元的值修改为1，如果之前已经是1(重叠作业)，无需修改，得到栅格化的农机运行轨迹图，如图6所示；According to the rectangle generated by the four points Q′ _i , Q″ _i , Q′ _i+1 , Q″ _i+1 and the rectangle generated by the two points Q′ _i+1 , Q″ _i+1

Superimpose the semicircle with the radius (for the smoothing of the work area) and the grid matrix: modify the value of this rectangle and the grid cells within the semicircle to 1. If it was already 1 before (overlapping work), there is no need to modify it, and you will get The rasterized agricultural machinery running trajectory is shown in Figure 6;

1-6-2) Calculate the area

According to the above grid matrix, count the number of grid cells with a value of 1, and then calculate the area according to the ratio of the cell to the actual size:

S _grid =N*w ²

where N is the number of grid cells with a value of 1.

application comparison,

A method of agricultural machinery behavior analysis and operation area statistics based on Beidou positioning data as described in the embodiment, the method is applied to 6 different plots to detect the running trajectories of agricultural machinery when operating in the above-mentioned different plots. 9a to 9f.

In order to verify the accuracy of the analysis and operation area statistics method described in the present invention, the actual operation area of the above 6 plots was manually measured to obtain the "actual operation area".

At the same time, the "boundary method" and "width method" described in the prior art are introduced to process the agricultural machinery operation data respectively, and obtain corresponding area data respectively.

The above-mentioned agricultural machinery operation area obtained by the method of the present invention, using the "boundary method" and "width method" is compared, as shown in Table 1:

Table 1:

Among them, "195001018747_0, 195001018747_1, 195001018747_2, 195001018747_3" refer to the three identification areas in the plot 195001018747;

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

"206003470996_0, 206003470996_1" refers to the two identified areas in the parcel 206003470996.

It can be seen from Table 1 that the "boundary method" has a large error in the area statistics when there are missing operations; while the "width method" takes the overlapping operations into account, and there is a large error; the present invention is based on contour-based and grid-based methods. The area calculation method can count the overlapping area and the missing area separately, and then eliminate the interference of the two, and the actual effective area error obtained is smaller than the error of the above comparison algorithm.

Claims (2)

1. a kind of agricultural machinery behavior analysis and operation area statistical method based on Beidou positioning data, it is characterized in that, comprise the automatic identification of agricultural machinery operation area, calculate area, analyze overlapping area and omission area; The method for automatic identification of the agricultural machinery operation area is an automatic identification algorithm of the agricultural machinery operation area based on spatial clustering; The method for calculating the area includes a grid-based area calculation method and a contour-based area calculation method; The method for analyzing the overlapping area is to calculate the area obtained by multiplying the track length by the width minus the area obtained by the grid-based area calculation method:

where d(Q _i , Q _i+1 ) represents the distance between adjacent track points of agricultural machinery operations; The method for analyzing the missing area is: S _miss =S _contour -S _grid ; Wherein, the Scontour is the area based on the contour; the Sgrid is the area based on the calculation grid; The automatic identification algorithm of agricultural machinery operation area based on spatial clustering includes: 1-1) Acquisition of agricultural machinery operation data The track data set P={P ₁ (t ₁ , lat ₁ , long ₁ ), P ₂ (t ₂ , lat ₂ , long ₂ ), ..., P _n is obtained through the vehicle-mounted Beidou positioning terminal and the GPRS mobile communication device. (t _n , lat _n , long _n ), where 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 eliminating data abnormal points, drift points, stop points and random noise points; 1-3) Projection Obtain a set of data points in the UTM coordinate system Q={Q ₁ (t ₁ , x ₁ , y ₁ ), Q ₂ (t ₂ , x ₂ , y ₂ ), ..., Q _n (t _n , x _n , y ) _n ); 1-4) Spatial clustering Using spatial clustering to identify the work area, the specific steps of identification are as follows: 1-4-1) Take each data point after preprocessing as the center of the circle, draw a circle with a certain radius r, and how many adjacent data points in the circle become the density value of the point; 1-4-2) If the density value of the point is less than the set threshold min_pts, then mark the point as a low density point, otherwise it is a high density point; 1-4-3) If a high-density point is inside the circle of another high-density point, connect the two points; if a low-density point is inside the circle of another high-density point, also connect it to On the nearest high-density point, it becomes the boundary point; 1-4-4) Repeat steps 1-4-2) and 1-4-3), remove the low-density points that are not in the circle of any high-density points, and retain the high-density point set as the trajectory points of the agricultural machinery operation area; 1-5) Calculate the area Scontour based on the contour; 1-6) Calculate the grid-based area Sgrid; 1-7) The method for analyzing the overlapping area is: Calculated as the area obtained by multiplying the track length by the width minus the area obtained by the grid-based area calculation method:

where d(Q _i , Q _i+1 ) represents the distance between adjacent track points of agricultural machinery operations; The method for analyzing the missing area is: S _miss =S _contour -S _grid ; Described 1-5) the method based on the area Scontour of the outline, the area calculation method based on the outline may further comprise the steps: 1-5-1) Concave Packet Calculation According to the data points obtained by clustering, the concave hull calculation is performed for each type of data point, and the specific steps are as follows: (1) Find the point with the smallest y value. If the y value is the same, take the point with the largest x value as the starting point O; (2) Starting from the starting point O, with (1, 0) as the reference vector, first find an edge smaller than the radius R as the initial edge, and this point is A at this time; (3) Loop to find the next edge. Assuming that the previous edge is AB, then the next edge must start from point B and connect to a point C in the R neighborhood of B. Use the following rules to find point C: The points in the R neighborhood are sorted in the polar coordinate direction with B as the center and the BA vector as the reference, and then a circle with BC _i as the chord is established for the points C ₀ ~ C _n in the R neighborhood of B in turn, and then checks whether it contains If other neighbor points do not exist, the chord is a new edge and jumps out of the loop; (4) Find all the edges in turn until no new edge is found or a point that has been used as an edge before is encountered; 1-5-2) Calculate the area Calculate the concave hull of each category to obtain multiple polygons, and then use the triangle segmentation algorithm or the Simpson algorithm to calculate the polygon area to obtain S _{contour_in} , that is

where w represents the working width, and d(Q _i , Qi ₊₁ ) _represents the distance between adjacent boundary points Qi and Qi ₊₁ ; Described step 1-6) the method for calculating grid-based area Sgrid is as follows: 1-6-1) Regional expansion, regional expansion is to rasterize the actual operation area according to the agricultural machinery track data and width. The specific steps are as follows: (a-1) Find the minimum value x _min , y _min and the maximum value x _max , y _max of x, y; (a-2) Determine the size of the grid matrix to be developed according to the ratio μ of the pixel to the actual size:

Where ε represents an additional boundary to ensure that the data points are all located in the matrix; open up a two-dimensional array representing the grid matrix and initialize it to 0; (a-3) Calculate the area to be expanded according to the trajectory and width of the agricultural machinery: Knowing the adjacent agricultural machinery running track points Q _i (t _i , x _i , y _i ), Q _i+1 (t _i+1 , x _i+i , y _i+1 ) and the working width w, carry out Area expansion is actually to find the coordinates of four points Q′ _i , Q″ _i , Q′ _i+1 , Q″ _i+1 :

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 the rectangle generated by the four points Q′ _i , Q″ _i , Q′ _i+1 , Q″ _i+1 and the rectangle generated by the two points Q′ _i+1 , Q″ _i+1

Superimpose the semicircle with the radius and the grid matrix: get the rasterized agricultural machinery running track map; 1-6-2) Calculate the area According to the above grid matrix, count the number of grid cells with a value of 1, and then calculate the area according to the ratio of the cell to the actual size: S _grid =N*μ ² where N is the number of grid cells with a value of 1. 2. a kind of agricultural machinery behavior analysis and operation area statistics method based on Beidou positioning data according to claim 1, is characterized in that, described step 1-2) method of data preprocessing comprises the order of excluding data: 1-2-1) Eliminate abnormal points: any agricultural machinery trajectory point should satisfy P(t, lat, long): lat∈[-90°, 90°] long∈[-180°, 180°] Remove data points that do not satisfy the above formula as outliers; 1-2-2) Eliminate drift points: for two adjacent track points P _i (t _i , lat _i , long _i ), P _i+1 (t _i+1 , lat _i+1 , long _i+1 ) to calculate the operating speed of the agricultural machinery:

where d(P _i , P _i+1 ) represents the distance between adjacent trajectory points P _i and P _i+1 , and the trajectory points with v(P _i P _i+1 )>v _max are eliminated, where v _max is The maximum operating speed of the agricultural machinery; 1-2-3) Eliminate the stopping point: Calculate the average speed for k consecutive agricultural machinery running track points:

Eliminate trajectory points whose average speed is less than a certain threshold δ; 1-2-4) Eliminate random noise points: for two adjacent trajectory points P _i (t _i , lat _i , long _i ), P _i+1 (t _i+1 , lat _i+1 , long _{i+ 1} ) Calculate its orientation:

The representation method of converting it into a unit vector: θ _{i, i+1} → (cos(θ _{i, i+1} ), sin(θ _{i, i+1} )), then for k consecutive agricultural machinery running trajectory points to calculate Its direction mean is:

Calculate its standard deviation as:

Eliminate points with standard deviation greater than a certain threshold.

Images