A Three-Dimensional Full-Coverage Operation Path Planning Method for Plant Protection Unmanned Aerial Vehicles Based on Energy Consumption Modeling

: With the rapid development of precision agriculture technology worldwide, plant protection unmanned aerial vehicles (UAVs) have been widely applied in the prevention and control of agricultural pests and diseases due to their advantages such as unrestricted terrain, strong ma-neuverability, hover capability, and high operational efﬁciency. In view of the fact that most of the current path planning methods for plant protection UAVs are mainly applicable to large ﬁelds with relatively ﬂat terrain, and there is a lack of three-dimensional path planning methods applicable to hilly orchards, as well as a lack of intuitive evaluation of the advantages and disadvantages of the paths by utilizing the amount of energy consumption, this paper proposes a three-dimensional (3D) path planning method for plant protection UAVs in hilly orchards. The method utilizes the three-dimensional (3D) coordinates of the operational area obtained from a realistic 3D model of the area, and determines the 3D operational path through the use of Boustrophedon Cellular Decomposition and the height variations of the operational area. Based on the induced power used by the plant protection UAV for ﬂight and the power needed to overcome wind resistance, combined with the change of the load of the plant protection UAV and the supply situation, the relative energy consumption model of the plant protection UAV is constructed. Then, by optimizing the operational heading angle (1–180 ◦ ), the 3D path with the minimum relative energy consumption is obtained. Through a simulation test, it is demonstrated that the path with the minimum relative energy consumption (operational heading angle of 133 ◦ ) reduces the average relative energy consumption by 62.44% and 47.21% compared to the path with the maximum relative energy consumption (operational heading angle of 51 ◦ ) and the average relative energy consumption, respectively. This proves the feasibility of the proposed path planning method in hilly orchard environments and provides a reference for the application of plant protection UAVs in such orchards.


Introduction
Pest and disease prevention is one of the most crucial elements of orchard management.Due to terrain and planting patterns, ground plant protection machinery has difficulty maneuvering between rows in hilly or mountainous orchard terrains [1,2].As such, pest and disease prevention still mainly relies on manual pesticide application, which is labor-intensive, inefficient, leads to significant pesticide wastage, and frequently results in incidents of pesticide poisoning among workers [3].Compared to manual application, ground-based plant protection machinery, and fixed-wing plant protection UAVs [4], plant protection UAVs, with their capabilities for autonomous take-off and landing, high mobility and agility, intelligence, low carbon footprint, and minimal environmental impact, have garnered widespread attention from researchers.They have become an effective means of pest and disease prevention in hilly regions, and their development continues to show rapid growth [5].Hilly orchards have significant terrain undulations, coupled with issues such as low standardization of orchards and inconsistent planting patterns [6][7][8].When planning the operation path for plant protection drones, it is necessary to adjust the flight height over time according to the requirements of plant protection operations and the terrain undulations and height changes of fruit trees.The radar fitted to existing plant protection UAVs is unable to provide timely feedback to ensure the operational safety of plant protection UAVs in the face of sudden changes in terrain.At present, the operation of plant protection UAVs in hilly orchards largely relies on manual remote control, which cannot guarantee the precision, efficiency, and effect of the operations.
As the development of plant protection UAVs is relatively recent, path planning for these UAVs is still in the stage of algorithm research; that is, coverage path planning based on the terrain changes and obstacle distribution in the operation area, combined with the battery capacity, liquid capacity, and spraying width of the plant protection UAV.Researchers have conducted a large amount of research on the comprehensive coverage path planning of plant protection UAVs.In terms of algorithm optimization, researchers have explored the coverage of the operational area based on the Boustrophedon Cellular Decomposition or the Internal Spiral method [9,10].Then, they have optimized the algorithm on the planned path using the gravitational search algorithm [11], ant colony algorithm [12], genetic algorithm [13], and annealing algorithm [14].However, the above studies mainly focus on two-dimensional path planning for large fields with flat terrain, and there are fewer studies on three-dimensional paths [15].In addition, most of the objectives of path planning are to obtain the path with the shortest operation length and the shortest operation time, and there is no energy consumption modeling for the plant protection UAV itself and its real-time change of load during operation.However, the shortest length and the shortest time cannot directly evaluate the energy consumption of the plant protection UAV.In the operation environment of hilly orchards, compared with large fields, the terrain changes more frequently, and the altitude difference within the operation area is often large; if the two-dimensional operation path is still used, the plant protection UAV can only operate on the higher horizontal plane to ensure the safety of the plant protection UAV.However, this mode of operation is susceptible to the influence of natural wind, increasing the degree of droplet drift, and cannot guarantee the effect of plant protection operations in the operation area at a lower altitude.Therefore, it is necessary to perform three-dimensional (3D) path planning [16] based on two-dimensional (2D) path planning to ensure the effectiveness of plant protection operations [17] and the safety of plant protection UAVs.In order to realize the affine 3D operation of plant protection UAVs, with the development of UAV technology and image recognition technology, the plant protection UAV operation method based on the heart of trees has attracted the attention of scholars because it can reduce the operating distance of plant protection UAVs [18].However, due to the precision of the 3D model or image recognition technology, the effect of plant protection operations cannot be guaranteed, and problems such as missed spraying and uneven spraying are likely to occur.The 3D reality model of the operation area can provide complete 3D information of the operation area, and coverage 3D path planning based on this can ensure the spraying effect of the plant protection UAV and avoid leakage of spraying caused by inaccurate image recognition.
In this paper, a 3D path planning method for plant protection UAVs in hilly orchards is proposed.The method consists of three parts: (1) 3D coordinate acquisition of the operation area.We used an aerial survey UAV (Phantom 4-RTK) to collect orchard images and establish a 3D model of a hilly orchard.Then, based on the 3D model, we divided the operation area to obtain the latitude, longitude, and elevation coordinate information of the area.(2) 3D path planning.In order to obtain the full-coverage 3D path of the operation area, we first imported the obtained 3D coordinate information into Matlab, extracted the boundary coordinates of the operation area, and used the Boustrophedon Cellular Decomposition to carry out 2D path planning.Then, on the basis of the 2D path, we combined the terrain change of the operation area and the flight altitude of the plant protection UAV to obtain the 3D path points, and expanded the 2D path to the 3D path.Finally, in order to obtain the plant planting UAV operation path with the minimum energy consumption in the operation area, we changed the operation heading angle of the plant planting UAV (1-180 • ), and obtained 180 operation paths in the operation area in 1 • steps.
(3) The energy consumption model of the plant protection UAV.The DJI T30 six-rotor plant protection UAV was taken as the object, which is widely used at present, to decompose the motion of the plant protection UAV in the horizontal and vertical directions, and combine the real-time change of the load of the plant protection UAV in the operation to construct the energy consumption model.Using the above model, the energy consumption of 3D paths under heading angles of 1-180 • was calculated separately, so as to obtain the 3D path with the smallest relative energy consumption.
The contributions of this paper mainly include three aspects: (1) A method of path planning based on a realistic 3D model is proposed.The method can obtain more accurate 3D coordinate information, which can guarantee the accuracy and safety of plant protection tasks performed by plant protection UAVs in subsequent research.(2) Expanding the 2D full-coverage planting path to a 3D path, which will be more suitable for planting UAVs to carry out planting operations in hilly orchards.(3) Evaluating the advantages and disadvantages of the operation paths by using the energy consumption model instead of path length, operation time, and other indexes usually used in the existing path planning, in order to obtain a better plant protection path.

Materials and Methods
The path planning of a plant protection UAV is a form of global path planning [10], i.e., with the known coordinate information of the operational area, combined with constraints such as operational width, battery capacity, flight height, flight energy consumption model, resupply distance, etc.The goal is to achieve the minimum energy consumption and shortest time under the premise of fully covering the operational area based on actual needs.
In this study, the relatively accurate 3D coordinate information of the operation area was obtained by constructing a realistic 3D model of the operation area, and the operation heading angle of the plant protection UAV optimized by establishing the relevant energy consumption model with the plant protection UAV (DJI T30) as the object, so as to obtain the operation path with the lowest relative energy consumption.The research process is shown in Figure 1.
boundary coordinates of the operation area, and used the Boustrophedon Cellular Decomposition to carry out 2D path planning.Then, on the basis of the 2D path, we combined the terrain change of the operation area and the flight altitude of the plant protection UAV to obtain the 3D path points, and expanded the 2D path to the 3D path.Finally, in order to obtain the plant planting UAV operation path with the minimum energy consumption in the operation area, we changed the operation heading angle of the plant planting UAV (1-180°), and obtained 180 operation paths in the operation area in 1° steps.
(3) The energy consumption model of the plant protection UAV.The DJI T30 six-rotor plant protection UAV was taken as the object, which is widely used at present, to decompose the motion of the plant protection UAV in the horizontal and vertical directions, and combine the real-time change of the load of the plant protection UAV in the operation to construct the energy consumption model.Using the above model, the energy consumption of 3D paths under heading angles of 1-180° was calculated separately, so as to obtain the 3D path with the smallest relative energy consumption.
The contributions of this paper mainly include three aspects: (1) A method of path planning based on a realistic 3D model is proposed.The method can obtain more accurate 3D coordinate information, which can guarantee the accuracy and safety of plant protection tasks performed by plant protection UAVs in subsequent research.(2) Expanding the 2D full-coverage planting path to a 3D path, which will be more suitable for planting UAVs to carry out planting operations in hilly orchards.(3) Evaluating the advantages and disadvantages of the operation paths by using the energy consumption model instead of path length, operation time, and other indexes usually used in the existing path planning, in order to obtain a better plant protection path.

Materials and Methods
The path planning of a plant protection UAV is a form of global path planning [10], i.e., with the known coordinate information of the operational area, combined with constraints such as operational width, battery capacity, flight height, flight energy consumption model, resupply distance, etc.The goal is to achieve the minimum energy consumption and shortest time under the premise of fully covering the operational area based on actual needs.
In this study, the relatively accurate 3D coordinate information of the operation area was obtained by constructing a realistic 3D model of the operation area, and the operation heading angle of the plant protection UAV optimized by establishing the relevant energy consumption model with the plant protection UAV (DJI T30) as the object, so as to obtain the operation path with the lowest relative energy consumption.The research process is shown in Figure 1.

Creation of the 3D Model of the Work Area
To plan the path of plant protection UAVs in hilly orchards, it is necessary to map the work area to build a 3D model of the work area based on the actual 3D terrain data.This serves as the data foundation for the path planning of plant protection UAVs [19].
The test area is located in Ding Yangling, Yuan Dong Township, Lanxi City, Jinhua City, Zhejiang Province (119.84282479• E-119.84375015• E, 29.27923194 • N-29.28023677 • N).The aerial photography of the orchard was taken using a DJI PHANTOM 4 RTK UAV (DJI Innovation Technology Co., Ltd., Shenzhen, China) to obtain images.Based on DJI Terra for 3D terrain modeling, the final high-resolution 3D real-world model was obtained (Figure 2).The operational boundary points were determined in the model, and the json file covering the longitude, latitude, and elevation information of the operational area (accuracy 0.5 m) was obtained through Smart 3D processing, which provided 3D coordinate information of the operational area for route planning.

Creation of the 3D Model of the Work Area
To plan the path of plant protection UAVs in hilly orchards, it is necessary to map the work area to build a 3D model of the work area based on the actual 3D terrain data.This serves as the data foundation for the path planning of plant protection UAVs [19].
The test area is located in Ding Yangling, Yuan Dong Township, Lanxi City, Jinhua City, Zhejiang Province (119.84282479°E-119.84375015°E, 29.27923194° N-29.28023677°N).The aerial photography of the orchard was taken using a DJI PHANTOM 4 RTK UAV (DJI Innovation Technology Co., Ltd., Shenzhen, China) to obtain images.Based on DJI Terra for 3D terrain modeling, the final high-resolution 3D real-world model was obtained (Figure 2).The operational boundary points were determined in the model, and the json file covering the longitude, latitude, and elevation information of the operational area (accuracy 0.5 m) was obtained through Smart 3D processing, which provided 3D coordinate information of the operational area for route planning.

Coordinate System Construction
The operation of a plant protection UAV is a reciprocal covering movement, i.e., after determining each operational boundary point, the plant protection UAV performs spraying operations in a polygonal area consisting of each boundary point [20,21].The work area in this study is located in the eastern longitude and northern latitude region (119.84282479°E-119.84375015°E, 29.27923194° N-29.28023677°N), and the longitude, latitude, and altitude values of each vertex of the work area were obtained separately for path planning to facilitate the calculation.By applying Miller Projection, the longitude and latitude values of the work area were converted into Cartesian coordinates, which were established with the point O as the origin, the east direction as the positive direction of the horizontal coordinate X axis, and the north direction as the positive direction of the vertical coordinate axis Y.The polygon work area is guaranteed to be located in the first quadrant of the OXY coordinate system, as shown in Figure 3.

Construction and Conversion of Coordinate System for the Work Area 2.2.1. Coordinate System Construction
The operation of a plant protection UAV is a reciprocal covering movement, i.e., after determining each operational boundary point, the plant protection UAV performs spraying operations in a polygonal area consisting of each boundary point [20,21].The work area in this study is located in the eastern longitude and northern latitude region (119.84282479• E-119.84375015• E, 29.27923194 • N-29.28023677 • N), and the longitude, latitude, and altitude values of each vertex of the work area were obtained separately for path planning to facilitate the calculation.By applying Miller Projection, the longitude and latitude values of the work area were converted into Cartesian coordinates, which were established with the point O as the origin, the east direction as the positive direction of the horizontal coordinate X axis, and the north direction as the positive direction of the vertical coordinate axis Y.The polygon work area is guaranteed to be located in the first quadrant of the OXY coordinate system, as shown in Figure 3.

Coordinate System Conversion
To facilitate the calculation of the energy consumption of the plant protection UAV on the operational path under different headings, it is necessary to perform a coordinate transformation so that the positive direction of the X' axis in the transformed coordinate system is parallel to the operational heading [13].The angle  between the transverse axis of the transformed coordinate system and the original coordinate axis is the heading angle, and the operational area is still located in the first quadrant of the new coordinate system.The method of coordinate transformation is as follows (Figure 4).In the new coordinate system X′OY′, the coordinates (xn, yn) of the vertex An are transformed into ( n x ′ , n y ′ ), i.e., cos sin sin cos where  is the angle (rad) between the transformed horizontal coordinate axis and the original horizontal coordinate axis, which is calculated as below:

Coordinate System Conversion
To facilitate the calculation of the energy consumption of the plant protection UAV on the operational path under different headings, it is necessary to perform a coordinate transformation so that the positive direction of the X' axis in the transformed coordinate system is parallel to the operational heading [13].The angle θ between the transverse axis of the transformed coordinate system and the original coordinate axis is the heading angle, and the operational area is still located in the first quadrant of the new coordinate system.The method of coordinate transformation is as follows (Figure 4).

Coordinate System Conversion
To facilitate the calculation of the energy consumption of the plant protection UAV on the operational path under different headings, it is necessary to perform a coordinate transformation so that the positive direction of the X' axis in the transformed coordinate system is parallel to the operational heading [13].The angle  between the transverse axis of the transformed coordinate system and the original coordinate axis is the heading angle, and the operational area is still located in the first quadrant of the new coordinate system.The method of coordinate transformation is as follows (Figure 4).In the new coordinate system X′OY′, the coordinates (xn, yn) of the vertex An are transformed into ( n x ′ , n y ′ ), i.e., cos sin sin cos where  is the angle (rad) between the transformed horizontal coordinate axis and the original horizontal coordinate axis, which is calculated as below: In the new coordinate system X OY , the coordinates (x n , y n ) of the vertex An are transformed into (x n , y n ), i.e., where θ is the angle (rad) between the transformed horizontal coordinate axis and the original horizontal coordinate axis, which is calculated as below: After the path planning under this heading angle is completed, the coordinates of the path points are transformed back to the original coordinate system.The rotation angle of the coordinate system becomes −θ, which is used to calculate the 2D coordinates of the internal waypoints within the path.The conversion equation is as follows: x n y n (5) After determining the operating heading angle, in the operation area, the boustrophedon is applied to carry out two-dimensional path planning, i.e., in accordance with the spraying width of the plant protection drone as an interval, forming a complete coverage of the operation area, and the center line of the spraying range is the two-dimensional flight path of the plant protection drone.When the coverage is completed, the adjacent center lines are connected to each other at the head and tail, and a complete two-dimensional path is obtained.
Let the angle between the operation path of the plant protection drone and the positive direction of the X axis be the heading angle θ (1-180 • ), the operation width be w, and the distance between the fruit trees in the work area be d.The steps of generating a 2D operation path are as follows: After the path planning under this heading angle is completed, the coordinates of the path points are transformed back to the original coordinate system.The rotation angle of the coordinate system becomes −, which is used to calculate the 2D coordinates of the internal waypoints within the path.The conversion equation is as follows: y y cos( ) sin( ) sin( ) cos( )

2D Path Planning
After determining the operating heading angle, in the operation area, the boustrophedon is applied to carry out two-dimensional path planning, i.e., in accordance with the spraying width of the plant protection drone as an interval, forming a complete coverage of the operation area, and the center line of the spraying range is the two-dimensional flight path of the plant protection drone.When the coverage is completed, the adjacent center lines are connected to each other at the head and tail, and a complete twodimensional path is obtained.
Let the angle between the operation path of the plant protection drone and the positive direction of the X axis be the heading angle  (1-180°), the operation width be w, and the distance between the fruit trees in the work area be d.The steps of generating a 2D operation path are as follows:

3D Path Planning
As shown in Figure 6, for the same work area, the change of the operation heading angle does not affect the total path length if the 3D terrain is not considered.However, when the 3D terrain is considered, the length of the operation path under different operation heading angles is different and the operation energy consumption changes [11].In addition, in a hilly orchard environment, due to the varying and frequent changes in terrain and fruit tree heights, plant protection UAVs are prone to uneven fog droplet settlement under the 2D operation path, which affects the plant protection operation effect and even fails to ensure the operational safety of the plant protection UAVs.Therefore, in a hilly orchard environment, it is necessary to carry out 3D path planning of plant protection UAVs.

3D Path Planning
As shown in Figure 6, for the same work area, the change of the operation heading angle does not affect the total path length if the 3D terrain is not considered.However, when the 3D terrain is considered, the length of the operation path under different operation heading angles is different and the operation energy consumption changes [11].In addition, in a hilly orchard environment, due to the varying and frequent changes in terrain and fruit tree heights, plant protection UAVs are prone to uneven fog droplet settlement under the 2D operation path, which affects the plant protection operation effect and even fails to ensure the operational safety of the plant protection UAVs.Therefore, in a hilly orchard environment, it is necessary to carry out 3D path planning of plant protection UAVs.As the continuity between each coordinate of the 3D data in the above work area cannot be guaranteed due to the influence of accuracy, it is necessary to import the 3D data of the work area into Matlab and carry out 3D terrain surface interpolation to obtain the 3D coordinates of the complete area.We decomposed the 2D path (Section 2.3.1)into multiple path points at intervals of fruit tree spacing in the orchard, and extracted the height value of each 2D path point in the 3D model.Based on the height of the plant protection operation, we obtained the 3D operation path under this heading angle.
To ensure effective screening of three-dimensional paths in terms of energy consumption, taking into account the algorithm computation time as well as the results of previous research, we chose to change the operating heading angle from the minimum value (1°) to the maximum value (180°) in steps of 1° to obtain 180 operating paths.Combined with the energy consumption model, the relative energy consumption required for the 180 paths was calculated, resulting in an operational heading angles for the path with the lowest energy consumption.The flow chart of the path planning algorithm is shown in Figure 7. Figure 8 shows the 3D operation path for the heading angle of 60° planned using this method.
Through the above method, it is possible to obtain the 3D operation path under a certain heading angle  .To study the change in energy consumption under different heading angles, the algorithm generates with the heading angle  as the variable.The range of  is 1° ≤  ≤ 180° (the 3D path under 180-360° is the same as the 3D path under 1-180°).Because the range of the variable changes is small, all the 3D operation paths in the range of 1° ≤  ≤ 180° are obtained with a step length of 1°.As the continuity between each coordinate of the 3D data in the above work area cannot be guaranteed due to the influence of accuracy, it is necessary to import the 3D da-ta of the work area into Matlab and carry out 3D terrain surface interpolation to obtain the 3D coordinates of the complete area.We decomposed the 2D path (Section 2.3.1)into multiple path points at intervals of fruit tree spacing in the orchard, and extracted the height value of each 2D path point in the 3D model.Based on the height of the plant protection operation, we obtained the 3D operation path under this heading angle.
To ensure effective screening of three-dimensional paths in terms of energy consumption, taking into account the algorithm computation time as well as the results of previous research, we chose to change the operating heading angle from the minimum value (1 • ) to the maximum value (180 • ) in steps of 1 • to obtain 180 operating paths.Combined with the energy consumption model, the relative energy consumption required for the 180 paths was calculated, resulting in an operational heading angles for the path with the lowest energy consumption.The flow chart of the path planning algorithm is shown in Figure 7. Figure 8 shows the 3D operation path for the heading angle of 60 • planned using this method.
Through the above method, it is possible to obtain the 3D operation path under a certain heading angle θ.To study the change in energy consumption under different heading angles, the algorithm generates with the heading angle θ as the variable.The range of θ is 1 • ≤ θ ≤ 180 • (the 3D path under 180-360 • is the same as the 3D path under 1-180 • ).Because the range of the variable changes is small, all the 3D operation paths in the range of 1 • ≤ θ ≤ 180 • are obtained with a step length of 1 • .

Calculation of Energy Consumption of Plant Protection UAVs
The energy consumption of a plant protection UAV includes flight power consumption, which ensures that the UAV remains at high altitude and supports its motion and communication energy consumption, such as signal processing, circuitry, and power amplification.The most important part of the energy consumption required by a plant protection UAV to accomplish the plant protection task is to carry the liquid for flight, so we chose to ignore the energy consumption used for signal processing, communication, and control components when constructing the energy consumption model.

Calculation of Energy Consumption of Plant Protection UAVs
The energy consumption of a plant protection UAV includes flight power consumption, which ensures that the UAV remains at high altitude and supports its motion and communication energy consumption, such as signal processing, circuitry, and power amplification.The most important part of the energy consumption required by a plant protection UAV to accomplish the plant protection task is to carry the liquid for flight, so we chose to ignore the energy consumption used for signal processing, communication, and control components when constructing the energy consumption model.
To compare the relative operating energy consumption of the plant protection UAV at different heading angles, its energy consumption was simplified according to the horizontal and altitude variations of adjacent operating path points [22].An energy consumption model was constructed, and the plant protection UAV was assumed to maintain a uniform speed during operation to facilitate the calculation.

Horizontal and Vertical Flight Energy Consumption Model
The flight energy consumption of a rotor UAV consists of the power component for flight induction and the power component for overcoming air resistance [23,24].
The power of flight  at a horizontal speed  is: where  is the gravity of the plant protection UAV, N; S is the total rotor area of the plant protection UAV, m 2 ;  is the speed scalar of the plant protection UAV in horizontal flight;  is the power parameter for the plant protection UAV when hovering; and ρ is the air density, kg/m 3 .
The plant protection UAV needs to spray pesticides while working, and its own mass changes in real time, so  can be calculated as below: To compare the relative operating energy consumption of the plant protection UAV at different heading angles, its energy consumption was simplified according to the horizontal and altitude variations of adjacent operating path points [22].An energy consumption model was constructed, and the plant protection UAV was assumed to maintain a uniform speed during operation to facilitate the calculation.

Horizontal and Vertical Flight Energy Consumption Model
The flight energy consumption of a rotor UAV consists of the power component for flight induction and the power component for overcoming air resistance [23,24].
The power of flight P xy at a horizontal speed V xy is: where G RW is the gravity of the plant protection UAV, N; S RW is the total rotor area of the plant protection UAV, m 2 ; V xy is the speed scalar of the plant protection UAV in horizontal flight; P h is the power parameter for the plant protection UAV when hovering; and ρ is the air density, kg/m 3 .The plant protection UAV needs to spray pesticides while working, and its own mass changes in real time, so G RW can be calculated as below: where m s is the unloaded weight of the plant protection UAV, kg; m l is the maximum liquid load, kg; v l is the liquid spraying speed, kg/s; t is the flight time, s; and g is the acceleration of gravity, m/s 2 .
The power of flying P Z at the vertical speed V Z is: Electronics 2023, 12, 4051 10 of 16 Wind resistance power P RW is: where C D 0 is the wind resistance coefficient; V is flight speed, m/s.When the plant protection UAV is flying horizontally, its flight power P TD is: When the plant protection UAV is flying horizontally, its flight power P v is:

.2. Energy Consumption Model
In the operating environment of hilly orchards, plant protection UAVs usually fly in 3D scenarios, and their flight states are not limited to horizontal or vertical flight, such as simultaneous rising and forward flight states.Therefore, we further studied the power consumption issues of drones in 3D flight and to establish a corresponding general flight power consumption model.Speed and mass are key factors determining UAV power consumption [22] and the total power required for multi-rotor drones can be studied by analyzing the corresponding vertical and horizontal power consumptions [25].
When a multi-rotor UAV is flying in a 3D scenario, its flight speed can be decomposed into three different directional speeds, namely V x , V y , and V z , as shown in Figure 9a.In this study, the agricultural drone needs to traverse 3D path points to perform protection operation.Therefore, according to the position changes between two adjacent work points, the decomposed speed of the plant protection UAV is divided into horizontal speed and vertical speed, as shown in Figure 9b.
acceleration of gravity, m/s 2 .
The power of flying  at the vertical speed  is: Wind resistance power  is: where C is the wind resistance coefficient; is flight speed, m/s.When the plant protection UAV is flying horizontally, its flight power  is: When the plant protection UAV is flying horizontally, its flight power  is:

.2. Energy Consumption Model
In the operating environment of hilly orchards, plant protection UAVs usually fly in 3D scenarios, and their flight states are not limited to horizontal or vertical flight, such as simultaneous rising and forward flight states.Therefore, we further studied the power consumption issues of drones in 3D flight and to establish a corresponding general flight power consumption model.Speed and mass are key factors determining UAV power consumption [22] and the total power required for multi-rotor drones can be studied by analyzing the corresponding vertical and horizontal power consumptions [25].
When a multi-rotor UAV is flying in a 3D scenario, its flight speed can be decomposed into three different directional speeds, namely  ,  , and  , as shown in Figure 9a.In this study, the agricultural drone needs to traverse 3D path points to perform protection operation.Therefore, according to the position changes between two adjacent work points, the decomposed speed of the plant protection UAV is divided into horizontal speed and vertical speed, as shown in Figure 9b.The calculation methods for horizontal speed components V xy and vertical speed components V z are: Combined with the changes in the mass of the plant protection UAV, flight time, etc., the energy consumption is calculated.The calculation formula for its work energy consumption W w is as follows: where t i is the time when the UAV moves to the i point, s; and t i+1 is the time when the UAV moves to the i + 1 point.After the pesticide is exhausted in a single operation, the plant protection UAV needs to return to the replenishment point.The replenishment energy consumption W R consists of the energy consumption used by the plant protection drone to return to the replenishment point from the breakpoint when it is empty, and the energy consumption used by the plant protection UAV to return to the breakpoint from the replenishment point when it is fully loaded.Suppose the breakpoint coordinates are (x 1 , y 1 , z 1 ), and the replenishment point coordinates are (x 2 , y 2 , z 2 ).Then, the time used for both full-load and no-load flights is as follows: where v is the flight speed, m/s.According to the energy consumption model described in Section 2.4.1, the energy consumption for a single recharge is calculated for a full-load flight (W f ) and an empty-load flight (W e ), respectively: Then, the energy required to resupply the plant protection UAV, W r , is: The total power consumption of the plant protection UAV when completing the operation with a heading angle of θ is shown in Equation (20):

Results
The DJI T30 plant protection UAV was selected as the theoretical model with a maximum flight time of 20 min.Since the maximum drug load is 30 kg and the pesticide spraying speed is 0.1 kg/s, the time for one application is 300 s, which is much less than the maximum flight time of 20 min.Therefore, it is assumed that when replenishing the liquid, the UAV battery is replaced to maintain sufficient power.The relevant energy consumption calculation parameters are shown in Table 1.10.

Results
The DJI T30 plant protection UAV was selected as the theoretical model with a max imum flight time of 20 min.Since the maximum drug load is 30 kg and the pesticide spray ing speed is 0.1 kg/s, the time for one application is 300 s, which is much less than th maximum flight time of 20 min.Therefore, it is assumed that when replenishing the liquid the UAV battery is replaced to maintain sufficient power.The relevant energy consump tion calculation parameters are shown in Table 1.Based on the above algorithm, the path planning results are shown in Table 2, with the objective of searching for the heading angle that has the minimum energy consumption.As seen from Table 2, the minimum energy consumption required for operation in the above-mentioned work area is at a heading angle of 133 • , which is 98.09 KJ.There is a difference of 163.09KJ and 87.73 KJ compared to the maximum energy consumption path (heading angle 51 • ) and average energy consumption (heading angle 1-180 • ), respectively.
The work time under the minimum energy consumption path (heading angle 133 • ) is 53.84 s and 37.46 s less than the maximum energy consumption (heading angle 51 • ) path and average energy consumption (heading angle 1-180 • ) work time, reducing the required time by 7.62% and 5.43%, respectively.Among them, the working time is reduced by 87 s (20.07%) and 43.37 s (11.12%), respectively; the replenishment time is increased by 28.16 s (10.32%) and 8.65 s (2.96%), respectively.
In paths with different heading angles, work energy consumption and work time always constitute the main components of total energy consumption and total time.Taking the operation path with a heading angle of 133 • as an example, the work energy consumption and work time account for 80.47% and 53.11% of the total energy consumption and total time, respectively.
From the above results, in the simulation test of this study, the working time under different operating heading angles has some changes, as mentioned in Section 2.3.2,unlike the 2D path planning, the path length under different operating heading angles will be changed due to the change of terrain in the 3D path planning.In addition, both the path with maximum energy consumption and the path with minimum energy consumption have a small difference between the working time and the replenishment time, but the difference between the relative energy consumption required for working and the energy consumption required for replenishment is huge, and the working energy consumption accounts for 80% of the total energy consumption required, so the change of the load during the working time is the main factor that affects the total energy consumption of the operation.The above conclusion verifies that the main factors of energy consumption of plant protection UAVs are the change of altitude during flight and the change of load carried.

Discussion
In this study, for the 3D path planning of plant protection UAVs in irregularly shaped operational areas, after optimizing the operational heading angles from 1 to 180 • , the relative energy consumption and total time required for plant protection operations were significantly reduced.The path with the lowest relative energy consumption (133 • ) reduced the required energy consumption by 62.44% compared to the path with the highest relative energy consumption (51 • ).Compared to the average energy consumption (heading angle 1-180 • ), the path with the lowest relative energy consumption reduced the required energy consumption by 47.21%.The gap between the path with maximum energy consumption and the path with minimum energy consumption is mainly due to the fact that the plant protection UAV needs to make frequent climbs on the path with maximum energy consumption (51 • ) to realize the change of flight altitude, which leads to such a large gap in the operating energy consumption (63.80% reduction), which further illustrates the necessity of performing three-dimensional path planning, and optimization of the operational heading angle.
In the operation area selected for this experiment, the plant protection UAV only needs one replenishment to complete the plant protection operation, and the energy consumption of replenishment accounts for about 10% of the total energy consumption.However, this does not mean that the energy consumption of replenishment can be ignored.In the face of a large area to be operated within or the continuous operation of plant protection tasks in multiple operation areas, the number of replenishments required by the plant protection UAV to perform the operation will be increased, and the energy consumption required for replenishment will also rise.
We believe that further research can be divided into two directions: on the one hand, the optimization of the obtained 3D path, and on the other hand, the further reduction of replenishment energy consumption.In the optimization of the 3D path, such as the processing of path smoothness, combined with the related contents of plant protection UAV dynamics, such as turning, variable speed flight, etc., the optimization of the path can be used to reduce the energy consumption of the flight during the whole operation process.When facing a large area to be operated upon, there will be multiple interruption points, and the plant protection UAV needs to fly back from the interruption point to the supply point for the replenishment of consumables such as liquid and batteries, so the location of the interruption point affects the amount of energy consumption of the replenishment.For the optimization of the location of the interruption point, ant colony algorithm, genetic algorithm, annealing algorithm, etc., can be added to make the replenishment in advance when the remaining amount of liquid and battery is low, and the interruption point can be closer to the replenishment point, so as to reduce the energy consumption of replenishment.
In actual applications, when plant protection UAVs operate according to the path of minimum relative energy consumption, deviations from the predetermined path are inevitable.This is related to the positioning system used by the UAV itself and the operational environment.UAVs require a more precise navigation positioning system and stronger anti-interference ability to achieve the expected operational effect.In future research, the study of the navigation positioning system and anti-interference ability of plant protection UAVs will be of great significance for the application of UAVs' path planning.

Conclusions
(1) We used aerial drones for aerial triangular oblique photography to construct a realistic 3D model, which was used to determine the operational area, and obtain 3D coordinate information of the operational area with an accuracy of 0.5 m to ensure the reliability of the coordinate information.We imported the 3D coordinate information of the operational area into Matlab, and performed fitting and reconstruction of the operational area, which provides an information basis for path planning.
(2) Based on the reconstructed 3D model, we extracted the boundary coordinates of the operational area.After setting relevant parameters, such as the operation width of the plant protection UAV and the fruit tree spacing, we applied the Boustrophedon method for 2D path planning, and combined it with the height change of the operational area for 3D path planning.We constructed an energy consumption model for the six-rotor plant protection UAVs, which are widely used currently, to calculate the energy needed for the plant protection drone to complete the plant protection task along the planned 3D path.
(3) We carried out case tests, determined the operational area, and planned the path.With the minimum energy required to complete the plant protection task as the goal and the operation heading angle of the plant protection UAVs as the variable, we optimized over the range of the 1-180 • operation heading angle.The results show that in the path with the minimum energy consumption (operation course angle 133 • ), compared to the path with the maximum energy consumption (operation course angle 51 • ), the average energy consumption decreased by 163.09KJ (62.44%) and 87.73 KJ (47.21%), respectively, which validates the rationality and feasibility of this method.

Figure 1 .
Figure 1.The research process of plant protection UAV path planning based on minimum energy consumption.

Figure 1 .
Figure 1.The research process of plant protection UAV path planning based on minimum energy consumption.

( 1 )
Import the 3D data of the work area into Matlab (ver.2020b, The MathWorks, Natick, MA, USA); (2) Obtain the coordinates of the boundary vertices of the work area; (3) Rotate the coordinate system with the heading angle θ to obtain the boundary vertices of the work area after rotation; (4) Cover the work area with parallel lines that are spaced w apart and parallel to the rotated X axis; (5) Obtain path points in the work area, located on the above parallel lines, spaced d apart; (6) Connect the head and tail of the path points located on two adjacent parallel lines in sequence to obtain the 2D operation path under this heading angle.An example of a 2D operation path is shown in Figure 5 (taking a heading angle of 60 • as an example).

( 1 )
Import the 3D data of the work area into Matlab (ver.2020b, The MathWorks, Natick, MA, USA); (2) Obtain the coordinates of the boundary vertices of the work area; (3) Rotate the coordinate system with the heading angle  to obtain the boundary vertices of the work area after rotation; (4) Cover the work area with parallel lines that are spaced w apart and parallel to the rotated X axis; (5) Obtain path points in the work area, located on the above parallel lines, spaced d apart; (6) Connect the head and tail of the path points located on two adjacent parallel lines in sequence to obtain the 2D operation path under this heading angle.An example of a 2D operation path is shown in Figure 5 (taking a heading angle of 60° as an example).

Figure 7 .
Figure 7. Flow chart of path planning algorithm.

Figure 7 .
Figure 7. Flow chart of path planning algorithm.

Figure 9 .
Figure 9. Schematic diagram of speed decomposition.Schematic diagram of speed decomposition.(a) for decomposition of speed in three directions; (b) for decomposition of speed in two directions.
The selected work area (119.84282479°E-119.84375015°E, 29.27923194° N 29.28023677°N) has a total area of 6027.3295m 2 with irregular boundary shape.The loca tion of the work area and the recharge point (119.84367639°E, 29.28023268° N) are show in Figure10.

Figure 10 .
Figure 10.The location of the work area and the recharge point.

Table 1 .
Parameters related to energy consumption calculation.The selected work area (119.84282479• E-119.84375015• E, 29.27923194 • N-29.28023677 • N) has a total area of 6027.3295m 2 with irregular boundary shape.The location of the work area and the recharge point (119.84367639• E, 29.28023268 • N) are shown in Figure

Table 1 .
Parameters related to energy consumption calculation.

Table 2 .
The results of path planning.