Design and Implementation of Evaluation Method for Spraying Coverage Region of Plant Protection UAV

Abstract

Plant protection UAVs are becoming the preferred plant protection method for agricultural pest control. At present, the evaluation of droplet distribution in aerial spraying is collected and evaluated after the completion of prevention and control operations, and there is a lack of real-time evaluation methods. Based on the flight parameter during the UAV plant protection process, real-time estimation of droplet distribution is the key to solving this problem and further improving the effectiveness of aerial spraying. This study proposes a merging algorithm for arbitrary polygonal regions, meshing the boundaries of the region, divide the mesh segments based on the overlapping meshes between the two regions, and connect the valid mesh connection segments of the two regions according to certain rules to obtain the intersection, union, and residual operation results between the regions. Afterwards, software based on this algorithm was developed and applied to generate spraying coverage regions, leakage spray regions, and repeated spray regions. The experimental results on theoretical and irregular routes show that the algorithm can accurately generate droplet distribution regions. The error of the calculation results with a mesh scale of 0.05 m is within 7‰, and the operating speed is above 30 Hz, meeting the real-time requirements. The smaller the mesh scale is, the higher the accuracy of the calculation results is, but the slower the calculation speed. Therefore, in practical applications, it is necessary to choose an appropriate mesh scale based on hardware computing power and accuracy level requirements. This study solves the problem of cumulative calculation of droplet distribution during the operation of plant protection UAVs, providing a basis for objectively evaluating the operation quality of plant protection UAVs and optimizing the setting of operation parameters.

1. Introduction

Diseases and pests are the major threat to agricultural production. According to the Food and Agriculture Organization of the United Nations, the loss of world food production caused by diseases and pests is approximately 24% [1,2]. Diseases and pests spread rapidly and have a wide range of spread, which can lead to crop yield reduction or even crop failure, and quality degradation [3]. However, manual control efficiency is low, and it is easy to threaten the health of operators. Ground mechanical control cannot be applied to special terrains [4]. UAVs are becoming the preferred method of plant protection due to their high operational efficiency, good prevention and control effects, good maneuverability, wide applicability, water and drug conservation, and low risk to operators [5,6,7,8]. However, the harm of leakage spraying is enormous, and the pests at the leakage spraying region may spread to other regions again, leading to a resurgence of pests. The phenomenon of repeated spraying can reduce the efficiency of pesticide spraying and pollute the soil [9,10].

In recent years, scholars have conducted extensive research on the spatial distribution of droplets in plant protection UAVs, aiming to clarify the drift patterns of droplets in different situations, and then take effective measures to reduce droplet drift and reduce leakage, and repeated spraying. The droplet distribution of aerial spraying is closely related to the UAV model for plant protection [11], flight load [12], spray parameters [13,14], droplet size [15,16], UAV rotor wind field [17,18], external wind field [19,20], flight operational parameters [21,22,23], and other factors. Based on these studies [24,25,26], it is found that different UAV models, mission loads, and spray parameters will cause a large change in the droplet distribution, while the speed and direction of the external wind field is the most important factor that causes droplet drift, and the flight operational parameters of plant protection UAV are the secondary factor that affects drift. Through research [27,28,29], a preliminary exploration was conducted to predict droplet drift during aerial spraying based on wind speed and direction parameters and to dynamically adjust the motion trajectory of plant protection UAV, which can effectively reduce leakage spraying and repeated spraying, and improve target deposition performance. However, these studies did not provide real-time evaluation of droplet deposition, and the specific droplet deposition effect still needs to be determined through subsequent measurements.

Predictably, applying the drift law of droplets to generate and analyze the spatial distribution of droplets throughout aerial spraying can effectively evaluate and improve the deposition performance of the target. Decomposing the region covered by droplets during the aerial spraying process into multiple time interval subregions, the spatial distribution map of droplets can essentially be simplified as a series of region merging operations between these subregions and with the region to be worked on; that is, performing a series of Boolean operations on these polygonal regions, including intersection, union, and difference. In recent years, the polygon Boolean operation methods published by researchers include algorithms that only support convex polygon operations [30,31] and algorithms that support both convex and concave polygon operations [32,33,34,35], to adapt to operations with different data structures, vector data structure algorithms [36,37,38], and raster data structure algorithms [39,40,41] are currently the main methods for solving polygon Boolean operations. These methods can solve polygon Boolean operation problems in various situations, but the dynamic construction of the region covered by the droplets includes the Boolean operations of polygons with holes, convex polygon, concave polygon, and two polygons with overlapping polygons, which lacks a general solution. In addition, the situation of polygon Boolean operations merging changes dynamically with the aerial spraying process of the plant protection UAV. The existing algorithms cannot effectively solve the problem of polygon region dynamic merging construction in all the above situations.

Therefore, this study proposes a new algorithm for polygon Boolean operation based on boundary meshing, which converts the operation between points and line segments into the operation of meshes and mesh connected segments, and solves the Boolean operation results of arbitrary polygons. Based on on-site measured flight data, real-time calculation of droplet deposition distribution is carried out to construct spraying coverage regions, leakage spraying regions, and repeated spraying regions, providing a reference basis for precise spraying, variable spraying, and improving operational effectiveness.

2. Materials and Methods

2.1. Instruments and Equipment

In the field, with coordinates (23°09′76″ N, 113°20′37″ E), first, manually handheld differential GPS positioning equipment is used to collect the boundary coordinates of the site as the boundary of the region to be sprayed. As shown in Figure 1, the self-developed data acquisition module is installed on an eight-rotor plant protection UAV. The flight controller visually controls the UAV to simulate pesticide spraying operations in the field, including forward, backward, left, right, turn, and short hover. Finally, the longitude and latitude coordinates of an irregular route are collected by self-developed PC software, which includes various flight conditions during the actual operation of the plant protection UAV.

The longitude and latitude coordinates are spherical coordinates and need to be converted into geodetic plane coordinates x and y through Gauss Kruger projection transformation.

$$x=X+\frac{{L_c}^2}{2}\cdot N\cdot \mathit{sin}B\cdot \mathit{cos}B+\frac{{L_c}^4}{24}\cdot N\cdot \mathit{sin}B\cdot {\mathit{cos}}^{3}B\cdot (5-t^2+9{\delta }^2+4{\delta }^4)+\frac{{L_c}^6}{720}\cdot N\cdot \mathit{sin}B\cdot {\mathit{cos}}^{5}B\cdot (61-58t^2+t^4)$$

(1)

$$y=N\cdot \mathit{cos}B+\frac{{L_c}^3}{6}\cdot N\cdot {\mathit{cos}}^{3}B\cdot (1-t^2+{\delta }^2)+\frac{{L_c}^5}{120}\cdot N{\cdot \mathit{cos}}^{5}B\cdot (5-18t^2+t^4+14{\delta }^4-58{\delta }^2\cdot t^2)+500000$$

(2)

L is the current longitude, L₀ is the longitude of the central meridian of the projection zone, and L_c = L − L₀. X is the meridian arc length from the equator to the parallel circle at latitude B, whose value is $X=c{\int }_0^B(1+e\cdot c^2{\mathit{cos}}^2{B}^{-3/2})d$. E is the first eccentricity of the ellipsoid, whose value is $e=2\partial -{\partial }^2$, and $\partial$ is the oblateness, whose value is $\partial =\sqrt{a^2/b^2-1}$. C is the radius of polar curvature, m; δ and t are auxiliary variables with values of $t=\mathit{tan}B$ and $\delta =e\cdot c\cdot \mathit{cos}B$. a and b are the long and short radii of the reference ellipsoid, respectively.

Based on the position coordinates, heading angle, and spray amplitude of the UAV, the operation coverage subregion is constructed. The construction rules are shown in Figure 2. The position coordinates of the drone are P_n(x_n, y_n), and the rotation angle when the coordinate axis X rotates in the same direction as the heading angle is α, $P_n^{\prime }{P}_n^{″}$ is perpendicular to the heading of UAV, the coordinates of $P_n^{\prime }$ and $P_n^{″}$ can be obtained through Equations (3) and (4).

When 0 ≤ α < π/2 or 3π/2 < α ≤ 2π:

$$\left[\begin{array}{cc}{x}_n^{\prime }& y_n^{\prime }\\ x_n^{″}& y_n^{″}\end{array}\right]=\left[\begin{array}{cc}{x}_n-\frac{d}{2}\cdot \mathit{sin}\alpha & y_n+\frac{d}{2}\cdot \mathit{cos}\alpha \\ x_n+\frac{d}{2}\cdot \mathit{sin}\alpha & y_n-\frac{d}{2}\cdot \mathit{cos}\alpha \end{array}\right]$$

(3)

When π/2 ≤ α < 3π/2:

$$\left[\begin{array}{cc}{x}_n^{\prime }& y_n^{\prime }\\ x_n^{″}& y_n^{″}\end{array}\right]=\left[\begin{array}{cc}{x}_n+\frac{d}{2}\cdot \mathit{sin}\alpha & y_n-\frac{d}{2}\cdot \mathit{cos}\alpha \\ x_n-\frac{d}{2}\cdot \mathit{sin}\alpha & y_n+\frac{d}{2}\cdot \mathit{cos}\alpha \end{array}\right]$$

(4)

2.2. Basic Principles of the Algorithm

The basic principle of the region merging algorithm based on boundary meshing is shown in Figure 3, using the union of region A and region B as an example for analysis. Step 1: as shown in Figure 3A, the two-dimensional drawing space is decomposed into countless meshes based on the original coordinate axis. The meshes penetrated by the region boundary are solved, and these meshes are connected in the order of the original region vertices to form a grid sequence. The mesh sequence is used to replace the original region boundary, which is the meshing of the region boundary.

Step 2: Classify the mesh sequence representing the boundary of the region, as shown in Figure 3B. The meshes are preliminarily divided into coincident meshes, overlapping meshes, valid meshes and invalid meshes. The valid mesh refers to the mesh representing the boundary of the new region in the region merging result. Therefore, the mesh of region A outside region B is the valid mesh, and the mesh of region B outside region A is also the valid mesh. Invalid mesh refers to mesh that do not belong to the boundaries of the new region, that is, the mesh of region A inside region B and the mesh of region B inside region A. Coincident mesh refers to the mesh where two regions coincide with each other. Using coincident meshes, the mesh sequence of each region is divided into several valid and invalid mesh segments. Then, extract the connection points of the valid mesh segments of two regions in the coincident meshes, which are called overlapping meshes.

Step 3: Connect the valid mesh segments of the two regions at the overlapping meshes according to certain rules to form a mesh sequence representing the new region, as shown in Figure 3C. Finally, extract the center point coordinates of the meshes as the region vertex coordinates to construct the new region boundary, which is the result of region merging.

2.3. Meshing of Region Boundary

The meshing of region boundaries is to rephrase them in the form of mesh coordinates. In order to reduce the number of operations, all vertices of the original region are first meshed, and then one of the meshes is used as the starting mesh. According to the order of vertex connections in the original region, grow the mesh to the next vertex until it returns to the starting mesh.

2.3.1. Meshing of Regional Vertices

The mesh coordinates (x, y) represent the internal region contained by the boundaries composed of vertices (x, y), (x, y + d_g), (x + d_g, y + d_g), and (x + d_g, y). Assuming that A_n−1(x_n−1, y_n−1) and A_n(x_n, y_n) are the vertices of the original region, and d_g is the scale value of the mesh division, the mesh coordinates ($x_{n-1}^{\prime },y_{n-1}^{\prime }$) and ($x_n^{\prime },y_n^{\prime }$) of the two vertices can be preliminarily calculated by using Equation (5).

$$\left[\begin{array}{ll}{x}_{n-1}^{\prime }& y_{n-1}^{\prime }\\ x_n^{\prime }& y_n^{\prime }\end{array}\right]=\left[\begin{array}{ll}⌊\frac{x_{n-1}}{d_g}⌋& ⌊\frac{y_{n-1}}{d_g}⌋\\ ⌊\frac{x_{n-1}}{d_g}⌋& ⌊\frac{y_{n-1}}{d_g}⌋\end{array}\right]\cdot d_g$$

(5)

Then the number of meshes n_x and n_y penetrated by line segment A_n−1A_n in the X and Y directions can be determined.

$$\left[n_x\hspace{1em}{n}_y\right]=\left[⌈\frac{\left|x_n-x_{n-1}^{\prime }\right|}{d_g}⌉⌈\frac{\left|y_n-y_{n-1}^{\prime }\right|}{d_g}⌉\right]$$

(6)

when using ($x_{n-1}^{\prime },y_{n-1}^{\prime }$) as the starting mesh, the main growth direction of the mesh is the coordinate axis of the side with higher n_x and n_y values, and the auxiliary growth direction of the grid is the coordinate axis of the other side. The mesh grows along the boundary line A_n−1A_n towards mesh ($x_n^{\prime },y_n^{\prime }$). All meshes penetrated by the region boundary A_n−1A_n along the way are sequentially added to the mesh sequence. Growth direction coefficient α determines the direction of mesh growth along the coordinate axis. When α = 1, mesh growth is along the positive direction of the coordinate axis; when α = −1, the mesh grows along the negative direction of the coordinate axis. α_x and α_y represents the growth direction coefficients in the X and Y directions, respectively.

$${\alpha }_x=\left\{\begin{array}{ll}\frac{x_n-x_{n-1}}{\left|x_n-x_{n-1}\right|}& \left(x_n-x_{n-1}\ne 0\right)\\ 1& \left(x_n-x_{n-1}=0\right)\end{array}\right.$$

(7)

$${\alpha }_y=\left\{\begin{array}{ll}\frac{y_n-y_{n-1}}{\left|y_n-y_{n-1}\right|}& \left(y_n-y_{n-1}\ne 0\right)\\ 1& \left(y_n-y_{n-1}=0\right)\end{array}\right.$$

(8)

2.3.2. Rules of Mesh Growth

The process of mesh growth is essentially searching for the mesh connected to the current mesh, and the connectivity rules of the mesh are its growth rules. This algorithm is based on a four mesh space composed of a starting mesh and three adjacent meshes. Firstly, the mesh connectivity rules within the four mesh space are determined, and then the last connected mesh is used as the starting mesh to search for the next four mesh space. Finally, all meshes penetrated by the boundary line are determined step by step. Assuming the starting mesh of the mesh growth is ($x_1^{\prime }$, $y_1^{\prime }$), the other three meshes of the four mesh space are ($x_1^{\prime }$, $y_1^{\prime }+{\alpha }_y\cdot d_g$), ($x_1^{\prime }+{\alpha }_x\cdot d_g$, $y_1^{\prime }+{\alpha }_y\cdot d_g$), and ($x_1^{\prime }+{\alpha }_x\cdot d_g$, $y_1^{\prime }$). As shown in Figure 4, the number of meshes that the boundary line penetrates within the four mesh space is 2 or 3, i.e., the meshes are two connected or three connected.

By calculating and analyzing the coordinate values in Figure 5, the rules of mesh growth can be obtained. When n_x ≥ n_y, a linear search should be performed along the direction of α_x. The first search should use a straight line of $x=x_1^{\prime }+\frac{1+{\alpha }_x}{2}\cdot d_g$. The intersection point between the straight line and the mesh boundary line is ($x_1^{\prime }+\frac{1+{\alpha }_x}{2}\cdot d_g$, $y_2$). When $y_2=y_1^{\prime }+\frac{1+{\alpha }_y}{2}\cdot d_g$, the meshes penetrated by A_n−1A_n are two connected, and the connected mesh is ($x_2^{\prime }$, $y_2^{\prime }$). When $y_2\ne y_1^{\prime }+\frac{1+{\alpha }_y}{2}\cdot d_g$, the number of meshes penetrated along the Y-axis direction should be calculate with $n_{y1}^{\prime }=⌈\frac{\left|y_2-y_1^{\prime }-\frac{1-{\alpha }_y}{2}\cdot d_g\right|}{d_g}⌉-1$; when $n_{y1}^{\prime }\le 0$, continue searching along the direction of α_x to find the intersection point between A_n-1A_n and the straight line $x=x_1^{\prime }+\frac{1+3{\alpha }_x}{2}\cdot d_g$, which is ($x_1^{\prime }+\frac{1+3{\alpha }_x}{2}\cdot d_g$, $y_3$). Then, the number of penetrating meshes along the Y-axis direction can be calculated with $n_{y2}^{\prime }=⌈\frac{\left|y_3-y_1^{\prime }-\frac{1-{\alpha }_y}{2}\cdot d_g\right|}{d_g}⌉-1$. When $n_{y2}^{\prime }\le 0$, the meshes penetrated by A_n−1A_n are two connected, and the connected mesh is ($x_2^{\prime }$, $y_1^{\prime }$); When $n_{y2}^{\prime }\ge 1$, the meshes penetrated by A_n−1A_n are three connected, and the connected meshes are ($x_2^{\prime }$, $y_1^{\prime }$) and ($x_2^{\prime }$, $y_2^{\prime }$). When $n_{y1}^{\prime }>0$, the meshes penetrated by A_n−1A_n are three connected, and the connected meshes are ($x_1^{\prime }$, $y_2^{\prime }$) and ($x_2^{\prime }$, $y_2^{\prime }$).

When n_x < n_y, a linear search should be performed along the direction of α_y. The first search should use a straight line of $y=y_1^{\prime }+\frac{1+{\alpha }_y}{2}\cdot d_g$. The intersection point between the straight line and the mesh boundary line is ($x_2$, $y_1^{\prime }+\frac{1+{\alpha }_y}{2}\cdot d_g$). When $x_2=x_1^{\prime }+\frac{1+{\alpha }_x}{2}\cdot d_g$, the meshes penetrated by A_n−1A_n are two connected, and the connected mesh is ($x_2^{\prime }$, $y_2^{\prime }$). When $x_2\ne x_1^{\prime }+\frac{1+{\alpha }_x}{2}\cdot d_g$, the number of meshes penetrated along the x-axis direction should be calculate with $n_{x1}^{\prime }=⌈\frac{\left|x_2-x_1^{\prime }-\frac{1-{\alpha }_x}{2}\cdot d_g\right|}{d_g}⌉-1$; when $n_{x1}^{\prime }\le 0$, continue searching along the direction of α_y to find the intersection point between A_n−1A_n and the straight line $y=y_1^{\prime }+\frac{1+3{\alpha }_y}{2}\cdot d_g$, which is ($x_3$, $y=y_1^{\prime }+\frac{1+3{\alpha }_y}{2}\cdot d_g$). Then, the number of penetrating meshes along the x-axis direction can be calculated with $n_{x2}^{\prime }=⌈\frac{\left|x_3-x_1^{\prime }-\frac{1-{\alpha }_x}{2}\cdot d_g\right|}{d_g}⌉-1$. When $n_{x2}^{\prime }\le 0$, the meshes penetrated by A_n−1A_n are two connected, and the connected mesh is ($x_1^{\prime }$, $y_2^{\prime }$); When $n_{x2}^{\prime }\ge 1$, the meshes penetrated by A_n−1A_n are three connected, and the connected meshes are ($x_1^{\prime }$, $y_2^{\prime }$) and ($x_2^{\prime }$, $y_2^{\prime }$). When $n_{x1}^{\prime }>0$, the meshes penetrated by A_n−1A_n are three connected, and the connected meshes are ($x_2^{\prime }$, $y_1^{\prime }$) and ($x_2^{\prime }$, $y_2^{\prime }$).

2.3.3. Processing of Original Mesh Sequence

Based on the above rules of mesh growth, starting from the mesh where the region’s vertex is located, it gradually grows within the four mesh space in the order of the original region’s vertex connection. When the mesh returns to the starting position again, all meshes that are penetrated by the region’s boundary line form a mesh sequence. There may be some coincident meshes in the original mesh sequence, which can cause errors in subsequent calculation results. Further approximation processing is required for the original grid sequence, removing coincident meshes, and adding new meshes based on mesh connectivity rules to ensure that all adjacent meshes meet mesh connectivity rules (Figure 6).

2.4. Classification Processing of Meshes

Before merging two regions, it is necessary to classify the meshes in the two region mesh sequences into coincident meshes, overlapping meshes, valid meshes, and invalid meshes, and then the valid meshes connection segments are segmented.

2.4.1. Extracting Coincident Meshes

Firstly, it is necessary to extract meshes with the same coordinate in both regions and extract mesh index numbers of the mesh sequence to construct an array. The program flowchart is shown in Figure 7. By comparing each mesh coordinates of region A and region B in a loop, the mesh index numbers of two regions that coincide are found. The coincident mesh index numbers of region A are stored in array N[m], and the coincident mesh index numbers of region B are also stored in array N’[m]. Therefore, there is a one-to-one correspondence between the coincident mesh numbers of region A and region B.

After extracting the coincident meshes, the mesh sequence shown on the left in Figure 8 can be obtained, which includes four typical cases. The first mesh A₀ of the mesh sequence in case 1 is coincident, and the last mesh is noncoincident, so it can be directly used for subsequent overlapping mesh extraction. In other cases, it is necessary to rearrange the mesh sequence so that the first mesh is the continuous coincident mesh with the minimum index number in the original sequence. All meshes before this mesh maintain the original sorting and are shifted to the last mesh of the original grid sequence to form a new mesh sequence, as shown on the right side of Figure 8.

2.4.2. Extracting Overlapping Meshes

After determining the coincident meshes of the two regions, it is necessary to extract overlapping meshes and binarize the noncoincident meshes. The flow chart of the binarization program for grid in region A is shown in Figure 9. The program needs to input the coincident mesh index number array N’[m] of region A, and mesh coordinates array of both region A and region B. By detecting the difference between adjacent coincident mesh index numbers, mesh index numbers with a difference of more than 1 are obtained, which are the boundary points between coincident and noncoincident meshes. These meshes are called overlapping meshes, and they all exist in pairs.

When $N^{\prime }\left[i+1\right]-N^{\prime }\left[i\right]\ne 1$, store i and i + 1 in the overlapping mesh serial number array $N_1\left[j\right]$, store $N^{\prime }\left[i+1\right]$ and $N^{\prime }\left[i\right]$ into the overlapping mesh index array $N_1^{\prime }\left[j\right]$. Finally, add 0 and m-1 as the first and last elements of $N_1\left[j\right]$, and $N^{\prime }\left[0\right]$ and $N^{\prime }\left[m-1\right]$ as the first and last elements of $N_1^{\prime }\left[j\right]$.

The noncoincident meshes between two adjacent overlapping meshes are called a connecting segment, and region A is divided into k connecting segments. When Equation $N_1^{\prime }\left[0\right]=0$ and $N_1^{\prime }\left[j\right]=n-1$ are true, that is, the first and last overlapping meshes are the first and last meshes of region A, $k=\frac{j+1}{2}$ can be obtained; otherwise, $k=\frac{j+1}{2}-1$. The connecting segment will be assigned a flag of 1 when all the meshes of this connecting segment are inside the other region; otherwise, the connecting segment will be assigned a flag of 0. The mesh segment with a flag of 1 is called a valid mesh segment.

Due to the fact that all mesh attributes within the same grid connecting segment are the same, this algorithm uniformly extracts the meshes after the odd overlapping mesh for binarizing, and stores the results of all connecting segments in the binary flag array $B_1\left[k\right]$ of region A. Similarly, the connecting segments of region B can be binarized.

2.4.3. Deletion Rules for Invalid Overlapping Meshes

Starting from overlapping meshes, connect all valid mesh segments of two regions one by one to form a new region. Therefore, the more overlapping grids, the more times the grid connection segments are connected. In order to reduce unnecessary calculations, overlapping meshes that cannot be directly connected to form new regions are removed from the overlapping mesh sequence in advance. These deleted overlapping meshes are called invalid overlapping meshes.

As shown in Figure 10, A_i–A_k+3 are the coincident meshes of region A, and B_m–B_m+13 are the coincident meshes of a region corresponding to region A in region B one by one. The overlapping meshes A_i and A_k+3 correspond to B_m and B_m+13, which are also overlapping meshes of region B, and such overlapping meshes are called valid overlapping meshes. The meshes of region B corresponding to A_i+3, A_j, A_j+5, and A_k are not overlapping meshes, and such overlapping meshes are called additional overlapping grids.

Therefore, there are three kinds of overlapping grid pairs. Two overlapping meshes that are both valid overlapping meshes are called valid overlapping mesh pairs, such as B_m and B_m+13; if only one of the two overlapping meshes is a valid overlapping mesh, it is called an additional valid overlapping mesh pair, such as A_i and A_i+3; two overlapping grids that are both additional overlapping meshes are called additional overlapping grid pairs, such as A_j and A _j+5.

The flag of the mesh segment before the valid overlapping mesh A_i is a₁. The flag of the mesh segment between the additional overlapping meshes A_i+3 and A_j is a₂. The flag of the mesh segment between the additional overlapping meshes A_k and A_k+3 is a₃. The flag of the mesh segment after the valid overlapping mesh A_k+3 is a₄; the flag of the mesh segment before the valid overlapping mesh B_m in region B is b₁. The flag of the mesh segment after valid overlapping mesh B_m+13 is b₂.

Taking A_j and A_j+5 as examples, analyze the deletion rules for additional overlapping mesh pairs. When a₃ = 1 and a₁ = a₄ = 0, retain additional overlapping grid pairs; otherwise, delete the additional overlapping mesh pairs and merge them into the mesh segment after A_j+5. Taking B_m and B_m+13 as examples, analyze the deletion rules for valid overlapping mesh pairs. When b₁ ≠ b₂ or b₁ = b₂ = 1 and a₁ = a₄ = 1, retain the valid overlapping mesh pairs; otherwise, delete the valid overlapping mesh pairs and merge the meshes between B_m and B_m+13 into the mesh segment before B_m and delete the corresponding overlapping meshes in region A. Taking A_i and A_i+3 as examples, analyze the deletion rules for additional valid overlapping mesh pairs. When a₁ = a₄ = 0 and a₂ = 1, retain additional valid overlapping mesh pairs; otherwise, delete additional valid overlapping meshes A_i and A_i+3.

2.4.4. Rearrange Mesh Segments

As shown in Figure 11, for the convenience of mesh connection, the mesh sequence after deleting invalid overlapping meshes is reordered, so that the first mesh segment is valid.

2.5. Connection Rules for Mesh Segments

Use the region with a large number of valid mesh segments as the active region for region merging. When both ends of the valid mesh segment are valid overlapping mesh pairs, and the other sides of the two valid overlapping mesh pairs are invalid mesh segments, it is called a general valid mesh segment. When there is at least one valid mesh segment at the other side of two overlapping mesh pairs, it is called a dual valid mesh segment. When one end of the valid mesh segment is an additional valid overlapping mesh pair, it is called an additional valid mesh segment. Analyze the mesh connection rules during the region merging based on the types and numbers of valid mesh segments.

2.5.1. Connection Rules for General Valid Mesh Segments

As shown in Figure 12, the overlapping mesh index number n₁, n₂, n₃, n₄ of the active region. First, find the overlapping mesh index numbers $n_1^{\prime }$ and $n_2^{\prime }$ of the passive region corresponding to n₁ and n₂, and solve for the minimum value $n_{min}^{\prime }$ in $n_1^{\prime }$ and $n_2^{\prime }$, then the validity flag B_min of the mesh segment after $n_{min}^{\prime }$ is used to search for another overlapping mesh pair $n_3^{\prime }$ and $n_4^{\prime }$. B = 0 indicates that the mesh segment of the passive region needs to be reversed during mesh connection, otherwise the original order will be maintained. B_C = 1 indicates that the merging result is a closed region, otherwise it is an unclosed region. The starting point of the active region’s valid mesh segment is s, and the length is l. The starting point of the passive region’s valid mesh segment is s′, and the length is l′. The merging situations of general valid mesh segments are shown in the Figure 13.

2.5.2. Connection Rules for Dual Valid Mesh Segments

The region formed by the dual valid mesh segment and another region must be closed, as shown in the program flowchart in Figure 14. It is necessary to extract the index numbers $m_1^{\prime }$, $m_2^{\prime }$, $m_3^{\prime }$, and $m_4^{\prime }$ of each overlapping mesh in the overlapping mesh array of the passive region. Equation $\left|m_3^{\prime }-m_2^{\prime }\right|=1$ represents that $n_3^{\prime }$ and $n_2^{\prime }$ are adjacent in the overlapping mesh array.

2.5.3. Connection Rules for Additional Valid Overlapping Mesh Segments

When the mesh segment of the active region is an additional valid, $n_2^{\prime }$ and $n_3^{\prime }$ are not overlapping meshes of the passive region. The valid mesh segment of the active region is connected to $n_2$ and $n_3$ to form a closed region, so $s=n_2 \mathrm{a}\mathrm{n}\mathrm{d} l=n_3-n_2+1$.

2.5.4. Connection Rules for Special Situations

When the active region only contains one pair of overlapping grids and one segment of invalid or valid meshes, the program flowchart for mesh connection is shown in Figure 15. For active region A, n is the total number of meshes, n₁ and n₂ are overlapping meshes and B_a is the validity flag of the mesh segment. For passive region B, m is the total number of meshes, $n_1^{\prime }$ and $n_2^{\prime }$ are overlapping meshes and B_b is the validity flag of the mesh segment. The total number of coincident meshes in two regions is n’. The number of new regions is B_n.

When all the elements of the mesh sequence are coincident grids, the active region and passive region are completely coincident, and the connection result is a closed region, which is an active region.

When there is no overlapping mesh between two regions, the subroutine is shown in Figure 16, the mesh connection result is two independent regions, or one of them, or an empty set.

3. Results

Design the software based on algorithm principles, and the software interface is shown in Figure 17. The software includes a data analysis mode selection area, which can be used for real-time online analysis or offline extraction of data documents for analysis. This interface selects the offline data document analysis mode. The mesh scale setting button can be used to set d_g.

3.1. Experiment of Theoretical Route

The polygon composed of points A₁(7.5,5), A₂(32.5,5), A₃(42.5,30), A₄(5,37.5), and A₅(2,20) is taken as the region to be worked on. The theoretically planned routes [42] are imported into the designed plant protection UAV operation map reconstruction software. The mesh scale is set to 0.05 m, and the software automatically generates the operation coverage region. Compare it with the theoretical calculation value to verify the software’s operation and calculation error under different starting operation directions. The experimental results are shown in

Table 1. Experimental results of theoretical route.

Starting Spraying Direction	Spraying Coverage Region (m2)		Repeated Spraying Region (m2)		Repeated Spraying Rate (%)
	Theoretical	Experimental	Theoretical	Experimental	Theoretical	Experimental
$\stackrel{⃑}{A_3{A}_2}$	1162.76	1162.84	70.66	71.07	6.08	6.11
$\stackrel{⃑}{A_3{A}_4}$	1297.74	1298.43	130.42	131.31	10.05	10.11
$\stackrel{⃑}{A_1{A}_2}$	1115.54	1115.67	58.37	58.15	5.23	5.21

. The experimental results show that the software calculation results are consistent with the theoretical calculation results. The maximum relative error of the spraying coverage region is 0.53 ‰, the maximum relative error of the repeated spraying coverage region is 6.8 ‰, the maximum relative error of the repeated spraying rate is 6 ‰, and the errors of the leakage spraying coverage region are all 0.

3.2. Experiment of Irregular Route

Import the irregular route into the software. The starting spraying direction angle is 180°, the spraying span of the UAV is 4 m, and the mesh scale is 0.05 m. The algorithm automatically calculates the spraying coverage region and leakage region. The operation results are shown in Figure 17. Export the data to obtain Figure 18. The solving process of the algorithm includes union, intersection, and remainder operations between arbitrary polygons.

From Figure 18a, it can be seen that the software has calculated four spraying coverage regions, where the black curve represents the outer boundary of the spraying coverage region, meaning that the interior of the region is the spraying coverage region; the other three are the inner boundaries of the spraying coverage region, and the interior of this region is the leakage spraying regions. From Figure 18b, it can be seen that the software has calculated six mutually independent leakage spraying regions, and the regions inside each boundary line are the leakage spraying regions. The results of region merging are all closed regions, and the connections between each vertex comply with mesh connectivity rules, resulting in accurate and reliable calculation results.

Comparing the running time of the software constructed in different regions under different data points can obtain

Table 2. Running time of constructing different regions of the software.

Number of Data Points	Constructing Spraying Coverage Region			Constructing Leakage Spraying Region
	Running Time (ms)	Number of Merging Regions	Maximum Number of Meshes in the Region	Running Time (ms)	Number of Merging Regions	Maximum Number of Meshes in the Region
100	783	1	888	1087	2	4882
200	1560	1	2562	5940	2	6542
300	2386	2	4340	16,982	3	9892
400	3274	3	5603	32,556	4	12,418
500	4267	4	6788	51,843	5	14,788
600	6004	4	9041	76,816	6	19,088
645	6663	4	9945	90,297	6	20,896

. From the table, it can be seen that with the increase in the number of regions and meshes, the running time of constructing the coverage region does not change significantly, while the running time of constructing the leakage spraying region significantly increases. As the number of regions and meshes increases, the difference in computational time between constructing spraying coverage region and constructing leakage spraying region becomes larger and larger. This is mainly due to the need to use the boundary of the pesticide spraying region to be worked on when calculating the leakage spraying region, which is relatively large and generates many meshes during the meshing process. Finding overlapping meshes and processing grid sequences require a lot of time.

By comparing the running time of the software constructing spraying coverage region with different mesh scales,

Table 3. Running time of constructing spraying coverage region with different mesh scales.

Mesh Scale (m)	Running Time (ms)
	100 Data	200 Data	300 Data	400 Data	500 Data	600 Data	645 Data
0.05	783	1560	2386	3274	4267	6004	6663
0.02	975	2274	4188	7198	12,166	18,445	21,021

can be obtained. From the data analysis in the table, it can be seen that when the mesh scale is 0.05 m, the average time for merging the regions of each data point is 10 milliseconds, which means that the algorithm can adapt to coordinate data update frequencies within 100 Hz; when the mesh scale is 0.02 m, the average time for merging regions of each data point is 33 milliseconds, which means that the algorithm can adapt to coordinate data update frequencies within 30 Hz.

Compared to the mesh scale of 0.02 m, the relative error of the spraying coverage region in the calculation result at the mesh scale of 0.05 m is 0.3‰, the error of the repeated spraying region is 5.9‰, and the relative error of the leakage spraying region is 6.8‰. Therefore, the error of the calculation result when the mesh scale is 0.05 m is within 7‰, but it increases the computational speed by three times.

Therefore, the performance of this algorithm is determined by the mesh scale. The larger the mesh scale, the simpler the boundary representation of the region, the smaller the storage space is occupied, and the shorter the calculation time of region merging, but the lower the accuracy of the calculation result. The smaller the mesh scale, the more accurate the boundary representation of the region, the larger the storage space occupied, and the longer the calculation time of region merging, but the higher the accuracy of the calculation results.

4. Discussion

In the process of solving the spraying coverage region, before forming a hollow region, the union calculation between regions is mainly carried out; after forming a hollow region, it is also necessary to complete the intersection operation between the boundary inside the hollow region and the new region. By solving the residual set of the spraying coverage region in the pesticide spraying region, the leakage coverage region can be obtained.

The algorithm in this study solves the dynamic calculation problem of union, intersection, and residual sets between arbitrary polygonal regions, and combines real-time flight parameter information of plant protection UAV to accumulate and calculate the distribution of fog droplets during the aerial spraying process in real time. Compared to the research results of [27,28,29], this study achieved the goal of real-time evaluation of droplet distribution, providing a reliable basis for precise planning, real-time replanning, and operation quality evaluation of plant protection UAV routes.

The performance of this algorithm is determined by the mesh scale. The larger the mesh scale, the simpler the boundary representation of the region, the smaller the storage space is occupied, and the shorter the calculation time of region merging, but the lower the accuracy of the calculation result. The smaller the mesh scale, the more accurate the boundary representation of the region, the larger the storage space occupied, and the longer the calculation time of region merging, but the higher the accuracy of the calculation results.

The spraying quality of plant protection UAV depends on the spraying coverage region and the uniformity of droplet distribution. The spraying coverage region in this manuscript can characterize the deposition area of the droplet on a two-dimensional plane, but it is not yet possible to characterize the uniformity of the deposition of the droplet in the spray coverage plane area and the deposition of the crop canopy at different heights. In further research, based on the experiments to obtain the uniformity data of droplet deposition under different parameters, the droplet deposition density model is established, and then the droplet deposition density is generated everywhere in the region while constructing the spraying coverage region, generating an accurate droplet distribution map.

5. Conclusions

This research proposes an algorithm for merging arbitrary polygonal regions based on boundary meshing, and develops software based on this algorithm. Combined with information collection hardware, a system for evaluating the spraying coverage region of plant protection UAV is constructed. The coverage region of pesticide droplets is constructed, displayed, and analyzed. The main conclusions were:

(1)

This algorithm can be used to dynamically construct the spraying coverage region of plant protection UAV, and it can provide a new idea for solving the intersection, union and difference of arbitrary polygons;

(2)

The system can calculate the spraying coverage region and leakage spraying region, with a calculation error of less than 7‰, generate spraying effect parameters, and mark the regions that need to be inspected;

(3)

Import the operational parameters of the plant protection UAV planning route, and use the algorithm in this paper to generate spraying coverage region, which can be used to estimate the operational effect of the planned route, thereby selecting the optimal operational parameters and guiding the precise planning of the operational route. In the process of plant protection UAV spraying, the quality of the operation can also be evaluated based on the real-time generated spraying coverage region generated by this algorithm. At the same time, the operation route can be replanned in real-time based on the generated leakage spraying region, which can assist in precise spraying and is expected to achieve the best control effect with the lowest pesticide consumption;

(4)

Further research should be carried out to clarify the distribution of droplet deposition density under different operating parameters, add it to the spraying coverage region to generate a complete droplet distribution map, and combine the droplet distribution map at different heights to generate droplet deposition at different heights of crop canopy in three-dimensional space.