A Multi-UAV Cooperative Path Planning Method for Mountain Monitoring

Abstract

The present invention provides a multi-UAV collaborative path planning method for mountain monitoring. This method uses a set of UAV clusters to collect images of the mountain to be measured, extracts feature points of the mountain to be measured by using the Harris corner detection algorithm, marks the peak feature points, valley feature points and foot-of-mountain feature points, calculates the coordinates of each feature point in the world coordinate system, and then obtains the three-dimensional model of the mountain to be measured. Then, the greedy algorithm is used to assign the feature points as task detection points to the UAV clusters, and based on the three-dimensional model of the mountain to be measured, the A-Star algorithm is used to perform three-dimensional path planning for each UAV assigned with detection tasks. The present invention uses automatic formation of UAVs to collaboratively collect images of the mountain to be measured, reduces manual input, can adapt to different terrains and environmental conditions at the same time, has a wider data collection coverage range, improves work efficiency, and provides strong support for the monitoring of mountain natural disasters.

Description

Technical Field

[0001] The present invention relates to the technical field of unmanned aerial vehicle path planning, and particularly relates to a multi-unmanned aerial vehicle collaborative path planning method for mountain monitoring. Background Art

[0002] A geological disaster refers to a geological process or geological phenomenon formed under the action of natural or human factors, which causes losses to human life and property and damages to the environment. In the prior art, artificial simple detection methods or instrument detection methods are usually used to detect and warn geological disasters. For example, personnel are dispatched to conduct regular screening and exploration at geological disaster hidden danger points, or multiple types of sensors are set at geological disaster hidden danger points to collect and process data according to a set period. Landslide is a typical type of geological disaster, which is characterized by strong suddenness, great destructiveness and high prediction difficulty, and it is difficult to make timely and accurate predictions through conventional artificial means.

[0003] In order to reduce the potential threat of landslides, professionals need to conduct regular inspections and construction protection on mountains. At present, the detection of landslides mainly relies on manual screening and ground exploration technologies. However, manual screening not only has a large labor intensity, low efficiency, and limited mountainous areas covered by each detection, but also the ground exploration technology is costly and is easily restricted by external factors such as weather and environment, making it difficult to achieve long-term and continuous monitoring of mountains. Summary of the Invention

[0004] Object of the Invention: The object of the present invention is to provide a multi-unmanned aerial vehicle collaborative path planning method for mountain monitoring that can automatically obtain a mountain model.

[0005] Technical Solution: A multi-unmanned aerial vehicle collaborative path planning method for mountain monitoring includes the following steps:

[0006] S1. Select several unmanned aerial vehicles with the same structure to form an unmanned aerial vehicle cluster. Each unmanned aerial vehicle is equipped with a camera and a radar, and control signals for each unmanned aerial vehicle in the unmanned aerial vehicle cluster are generated based on a set unmanned aerial vehicle kinematic model;

[0007] S2. Construct an unmanned aerial vehicle cluster formation based on a consensus control protocol. Set one unmanned aerial vehicle in the unmanned aerial vehicle cluster as a leader unmanned aerial vehicle, and the remaining unmanned aerial vehicles as follower unmanned aerial vehicles. During flight, the follower unmanned aerial vehicles are arranged on the same side of the leader unmanned aerial vehicle in a fixed order in turn;

[0008] S3. Control the UAV cluster formation set in step S2 to start from the starting point, fly at the set speed and approach the mountain to be measured. Use the cameras carried by each UAV to collect images of the mountain to be measured at a preset time interval, and mark the corresponding position information on each image of the mountain to be measured. Use the radars carried by each UAV to measure the distance between each UAV and the mountain to be measured in real time, and obtain the inclination angle and inclination direction of the mountain to be measured. By setting a distance threshold between the UAV cluster and the mountain to be measured, when the distance between any UAV in the UAV cluster and the mountain to be measured is less than or equal to the distance threshold, the UAV cluster stops collecting images of the mountain to be measured and returns to the starting point.

[0009] S4. Convert the images of the mountain to be measured obtained in step S3 into grayscale images, use the Harris corner detection algorithm to extract feature points of the mountain to be measured, and mark the extracted feature points as mountaintop feature points, mountain depression feature points and mountain foot feature points respectively according to the gradient change direction. Calculate the depth from each feature point to the corresponding UAV according to the camera imaging principle, and obtain the average coordinates of each feature point in the world coordinate system using the transformation matrix.

[0010] S5. Perform multivariate polynomial fitting based on the coordinates of each feature point in the world coordinate system obtained in step S4 to obtain the ridge curve function of the mountain to be measured. Set the ridge curve generation range, extract the ridge curve segments and connect the endpoints of each ridge curve segment to form a closed cross-section diagram of the mountain to be measured. Then, according to the inclination angle and inclination direction of the mountain to be measured obtained by the radar, perform coordinate transformation on the feature points on the cross-section diagram of the mountain to be measured, and sweep the transformed feature points along the inclination direction to obtain the three-dimensional model of the mountain to be measured.

[0011] S6. Use the mountaintop feature points and mountain depression feature points as task detection points, calculate the Manhattan distance between each UAV and each task detection point, and use the greedy algorithm to perform detection task allocation with the goal of minimizing the total detection distance, so that each task detection point is assigned to a UAV, and the UAVs not assigned to the detection task stay at the starting point.

[0012] S7. Use the three-dimensional model of the mountain to be measured obtained in step S5, select the Manhattan distance as the evaluation function, and use the A-Star algorithm to perform three-dimensional path planning for each UAV assigned to the detection task in step S6.

[0013] Specifically, step S1 includes:

[0014] Obtain the second-order consensus control input of UAV i using the second-order motion equation of UAV i :

[0015] ,

[0016] where: 、 , are the position coordinates of the UAV i in the x, y, and z directions in the world coordinate system , respectively. , , are the velocity vectors of the UAV i in the x, y, and z directions, respectively. , , are the components of the second-order consensus control input of the UAV i in the x, y, and z directions.

[0017] Specifically, step S2 includes the following sub-steps:

[0018] S21. Calculate the relative distance between each follower UAV and the leader UAV;

[0019] S22. The leader UAV and each follower UAV take off from the starting point in sequence and fly to the preset altitude, and then fly at the preset altitude;

[0020] S23. Set the speeds of each follower UAV and the leader UAV to be the same, and control the leader UAV to fly on the leftmost side of the formation. Each follower UAV is arranged on the right side of the leader UAV in sequence at a fixed distance according to the three-dimensional space-consensus control protocol with spacing. The control input of the follower UAV is:

[0021] ,

[0022] where: is the control input of the follower UAV i at time t; is the fixed gain; is the damping gain; is the connection strength between the follower UAV i and the follower UAV j; i and j are the UAV numbers; is the relative distance between the follower UAV i and the leader UAV at time t; is the relative distance between the follower UAV j and the leader UAV at time t; is the expected relative distance between any two UAVs; is the position distance between the follower UAV i and the leader UAV on the z-axis at time t in the world coordinate system ; is the position distance between the follower UAV j and the leader UAV on the z-axis at time t in the world coordinate system ; , are the velocity amplitudes of the follower UAV i and the follower UAV j at time t, respectively; is the total number of UAVs in the UAV cluster;

[0023] S24. Under the action of the control input of the following drone, the following drone realizes the tracking of the desired position by tracking the desired relative distance until the speed and position states of the drone cluster satisfy the consistency control protocol with spacing in three-dimensional space, and the formation of the drone cluster is completed.

[0024] Specifically, step S4 includes the following sub-steps:

[0025] S41. Convert the image of the mountain to be measured into a grayscale image by using the maximum value method;

[0026] S42. Apply the Harris corner detection algorithm to segment the grayscale image converted in step S41, detect the feature points of the segmented grayscale image, obtain all the feature points in the grayscale image, and then mark the identified feature points as mountaintop feature points, mountain depression feature points and mountain foot feature points respectively according to the gradient change direction;

[0027] S43. Optionally select two grayscale images collected by two different drone cameras at the same moment, use the stereo matching algorithm based on Census transform to perform stereo matching on the same feature points in the two grayscale images, and obtain the disparity of the same feature points in the two grayscale images;

[0028] S44. Use the disparity of the same feature points in the two grayscale images obtained in step S43, and according to the camera imaging principle and the law of similar triangles, calculate the depth of the feature point from the camera carried by one of the two drones in step S43, and then obtain the three-dimensional coordinates of the feature point in the camera coordinate system, and then use the camera parameters to convert the three-dimensional coordinates of the feature point in the camera coordinate system into the coordinates in the world coordinate system;

[0029] S45. Repeat steps S43 and S44 to obtain the coordinate set of the world coordinate system of the same feature point obtained from the grayscale images collected by any two drones at each acquisition moment, and calculate the average coordinate of the feature point in the world coordinate system;

[0030] S46. Repeat steps S43 to S45 to obtain the average coordinates of all feature points in the world coordinate system.

[0031] Specifically, step S43 includes: using the stereo matching algorithm based on Census transform to match the grayscale images collected by two cameras, selecting the pixel corresponding to a certain feature point in one grayscale image , calculate the Hamming distance between the Census value of each feature point pixel in the other grayscale image and the Census value of the pixel , and obtain the feature point pixel with the smallest Hamming distance from the pixel , perform stereo matching on the same feature points in the grayscale images collected by two cameras, and then obtain the pixel points of the feature points on the two grayscale images 、 The disparity in the horizontal direction is as follows:

[0032] ,

[0033] In the formula: is the Census value of the pixel point corresponding to the feature point in the k-th acquisition by the camera carried by the exploration UAV i; is the bit-by-bit concatenation operation of bits; is the neighborhood window of the pixel point corresponding to the feature point ; q is the pixel in the neighborhood window; 、 are the grayscale values of the pixel corresponding to the feature point and the neighborhood pixel q respectively; is the comparison operation;

[0034] ,

[0035] In the formula: is the Hamming distance between the pixels and ; is the Census value of the pixel point corresponding to the feature point in the k-th acquisition by the camera carried by the exploration UAV j; E is the number of bits of the Census transform; is the exclusive OR operation; It is stipulated that in the two grayscale images to be matched, the exploration UAV i should be located on the left side of the exploration UAV j; is the comparison operation;

[0036] ,

[0037] In the formula: is the disparity between the matching pixel points of the feature point and in the images acquired by the cameras carried by the exploration UAVs i and j in the k-th acquisition; 、 are the horizontal pixel coordinates of the matching points in the image pixel coordinate system 、 respectively.

[0038] Specifically, step S44 includes:

[0039] Calculate the depth of the feature point from the camera carried by the exploration UAV i and its three-dimensional coordinates in the camera coordinate system according to the camera imaging principle and the law of similar triangles. The formula is:

[0040] ,

[0041] In the formula: , , are the coordinates of the feature point in the x, y, and z directions in the camera coordinate system obtained from the grayscale images captured by cameras i and j in the k-th acquisition; i and j are the numbers of the two UAVs respectively; f is the camera focal length; is the expected relative distance between the UAVs; , are the x and y direction coordinates of the feature point in the image pixel coordinate system of camera i; is the depth of the feature point from the camera carried by exploration UAV i deduced according to the camera imaging principle and the law of similar triangles in the images captured by the cameras carried by exploration UAVs i and j in the k-th acquisition; , are the principal point coordinates of camera i;

[0042] Construct a transformation matrix using the rotation matrix and translation matrix of the camera to convert the three-dimensional coordinates of the feature point in the camera coordinate system into the coordinates in the world coordinate system:

[0043] ,

[0044] In the formula: is the coordinate of the feature point in the world coordinate system obtained from the grayscale images captured by cameras i and j in the k-th acquisition; is the coordinate of the feature point in the camera coordinate system obtained from the images captured by cameras i and j in the k-th acquisition; is the rotation matrix of camera i, is the translation matrix of camera i, obtained from the camera parameters.

[0045] Specifically, step S5 includes the following sub-steps:

[0046] S51. Perform multivariate polynomial fitting using the average coordinates of each feature point in the world coordinate system, obtain the optimal parameters of the multivariate polynomial through the least error method, and construct a ridge curve function.

[0047] S52. Set the generation range of the ridge curve according to the average coordinates of the mountaintop feature point and the foot-of-mountain feature point in the world coordinate system:

[0048] ,

[0049] where: , are the x-axis and y-axis coordinates of the mountain to be measured in the world coordinate system respectively; is the average coordinate in the x-axis direction of the foot-of-mountain feature point 1, is the average coordinate in the x-axis direction of the foot-of-mountain feature point 2, is the average coordinate in the y-axis direction of the mountaintop feature point i; S is the number of mountaintop feature points;

[0050] S53. Extract the ridge curve segments according to the ridge curve function and the ridge curve generation range, and connect the endpoints of the ridge curve segments to obtain a closed cross-sectional view of the mountain to be measured.

[0051] S54. According to the inclination angle and inclination direction of the mountain to be measured obtained by the radar, establish a reduction matrix and a rotation matrix about the y-axis to transform the feature points on the cross-sectional view of the mountain to be measured. The rotation matrix and the reduction matrix are as follows:

[0052] ,

[0053] ,

[0054] where: is the rotation matrix about the y-axis; is the inclination angle of the mountain to be measured obtained by the radar; is the reduction matrix; , and are the reduction ratios in the x, y, and z-axis directions respectively;

[0055] S55. Sweep the feature points on the transformed cross-sectional view of the mountain to be measured along the inclination direction to obtain a three-dimensional model of the mountain to be measured.

[0056] Preferably, step S5 further includes:

[0057] S56. Use Laplacian smoothing to reduce the unnatural transition at the edge of the three-dimensional model of the mountain to be measured. The formula is:

[0058] ,

[0059] In the formula: is the position of the point that needs to be smoothed; is the position of the point after smoothing; is the point adjacent point set; is the set size of.

[0060] Specifically, step S7 includes the following sub-steps:

[0061] S71. Grid the task three-dimensional space of a certain drone assigned with the detection task, and determine the initial position coordinates and task detection point coordinates of the drone;

[0062] S72. Use the Manhattan distance as the evaluation function of the A-Star algorithm. Starting from the initial position coordinates of the drone, select the node with the minimum value of the evaluation function among its adjacent nodes as the next central node, and repeat this selection method until the task detection point coordinates are reached;

[0063] S73. Perform backward tracing from the task detection point coordinates to construct a flight path from the initial position coordinates of the drone to the task detection point coordinates;

[0064] S74. Repeat steps S71 to S73 until all drones assigned with the detection task complete the corresponding complete path planning.

[0065] Preferably, step S73 further includes:

[0066] Smooth the flight path using a third-order Bezier curve, use the nodes on the flight path obtained by the A-Star algorithm as the control points of the third-order Bezier curve, and use 4 consecutive nodes as a group of control points. The third-order Bezier curve formula is as follows:

[0067] ,

[0068] In the formula: is the output point of the third-order Bezier curve; , , and are 4 consecutive Bezier control points; is the variation parameter.

[0069] Beneficial effects: Compared with the prior art, the remarkable effect of the present invention is:

[0070] 1. Use drones to automatically form a formation and cooperate to collect images of the mountain to be measured. Compared with manual measurement or single-plane measurement, it reduces the investment of manpower. The drone formation can adapt to different terrains and environmental conditions, can perform collection operations in areas that are difficult for humans to reach, and has a wider data collection coverage, effectively shortening the data collection time and improving work efficiency.

[0071] 2. Use the cameras carried by the drone formation to collect images of the mountain to be measured multiple times, providing more image sets for visual ranging for target point ranging. Moreover, the images obtained by different drones from different angles are conducive to reducing the parallax error caused by a single perspective and improving the accuracy of distance calculation to ensure the reliability of the three-dimensional modeling of the mountain to be measured.

[0072] 3. By identifying, modeling, and using the mountaintop, mountain depression, and mountain foot as the characteristic points of the mountain and using them as the end points for drone detection, and matching them one by one with the drones for monitoring, the efficiency of drone image collection is achieved.

[0073] 4. This method can be applied to the monitoring of various mountain natural disasters. For example, it can achieve long-term and continuous landslide monitoring of the mountain, and combine machine learning models to predict and warn of mountain landslides, improving the automation level of natural disaster monitoring. BRIEF DESCRIPTION OF THE DRAWINGS

[0074] Figure 1 is the flowchart of the method of the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0075] The following further illustrates a preferred embodiment of the present invention with reference to the accompanying drawings.

[0076] Please refer to Figure 1 as shown, the present invention provides a multi-drone cooperative path planning method for mountain monitoring, including the following steps:

[0077] S1. Select several drones with the same structure to form a drone cluster. Each drone is equipped with the same camera and radar, satisfies the same kinematic model, and establish a set of N drones:

[0078] ,

[0079] Generate control signals for each drone in the drone cluster based on the drone kinematic model. Specifically, use the second-order motion equation of drone i to obtain the second-order consensus control input of drone i:

[0080] :

[0081] ,

[0082] Where: , , are the position coordinates of the UAV i in the x, y, and z directions in the world coordinate system respectively, , , are the velocity vectors of the UAV i in the x, y, and z directions respectively, , , are the components of the second-order consensus control input of the UAV i in the x, y, and z directions respectively.

[0083] According to the kinematic geometric relationship of the UAV, it can be expressed as:

[0084] ,

[0085] where: is the velocity amplitude of the UAV i; is the velocity of the UAV i; is the heading angle of the UAV i; is the pitch angle of the UAV i.

[0086] S2. Construct a UAV swarm formation based on the consensus control protocol, adopting the "leader-follower" method framework. During the entire exploration mission, the exploration UAV formation needs to maintain a specific flight altitude, with the nose facing the mountain to be explored, and the formation arranged in a linear "one" shape. Set one of the UAVs in the UAV swarm as the leader UAV, and the rest as follower UAVs. During flight, the follower UAVs are arranged in a fixed order on the same side of the leader UAV. In this embodiment, the leader UAV is located on the leftmost side of the entire formation, and the rest of the follower UAVs are arranged on the right side of the leader exploration UAV in order of number, at a fixed distance. This formation can not only effectively relieve the communication and computing pressure but also ensure the reliability of subsequent visual ranging of each feature point of the mountain.

[0087] This step specifically includes the following sub-steps:

[0088] S21. Calculate the relative distance between each follower UAV and the leader UAV. Let the leader UAV be , and its position coordinate in the world coordinate system is , and the velocity is . Let the follower UAV i be , and the relative distance between the follower UAV i and the leader UAV is:

[0089] ,

[0090] In the formula: is the relative distance between the following drone i and the leading drone 0 (0 is the number) at time t; and are respectively the position coordinates of the leading drone 0 and the following drone i in the world coordinate system at time t.

[0091] S22. The leading drone and each following drone take off from the starting point in sequence and fly to the preset flight altitude, and then maintain that altitude for flight.

[0092] S23. It is set that the speeds of each following drone and the leading drone are kept consistent, and the leading drone is controlled to fly on the leftmost side of the formation. Each following drone is arranged on the right side of the leading drone in sequence at a fixed distance according to the number order. The coordinate position is described as that each following drone needs to have the same z-axis coordinate as the leading drone, and the x-axis coordinates are equally spaced:

[0093] ,

[0094] In the formula: i and j are drone numbers; and are respectively the speed amplitudes of the following drones i and j at time t; and are the z-axis positions of the following drones i and j in the world coordinate system at time t; and are respectively the distances of the following drones i and j relative to the leading drone at time t; is the expected relative distance between each drone.

[0095] According to the consistency control protocol with spacing in three-dimensional space, the control input of the following drone is:

[0096] ,

[0097] In the formula: is the control input of the following drone i at time t; is a fixed gain; is a damping gain; is the connection strength between the following drone i and the following drone j; i and j are drone numbers; is the relative distance between the following drone i and the leading drone at time t; is the relative distance between the following drone j and the leading drone at time t; is the expected relative distance between any two drones; is the following drone i and the leading drone at time t in the world coordinate system The position distance on the lower z-axis; is the position distance on the lower z-axis at time t of the follower UAV j and the leader UAV in the world coordinate system The position distance on the lower z-axis; 、 are the velocity amplitudes of the follower UAV i and the follower UAV j at time t, respectively; is the total number of UAVs in the UAV cluster.

[0098] S24. Under the action of the control input of the follower UAV the follower UAV realizes the tracking of the desired position by tracking the desired relative distance until the velocity and position states of the UAV cluster satisfy the consistency control protocol with spacing in three-dimensional space, and the UAV cluster formation is completed.

[0099] S3. The UAV cluster formation set in control step S2 starts from the starting point, flies at the set speed and approaches the mountain to be measured. The cameras carried by each UAV collect images of the mountain to be measured at a preset time interval and mark the corresponding position information on each image of the mountain to be measured, so as to realize the accurate association of spatial position and image data. The radar carried by each UAV measures the distance between each UAV and the mountain to be measured in real time, and obtains the inclination angle and inclination direction of the mountain to be measured. To prevent the exploration UAV from colliding with the mountain during the forward process, it is necessary to set a distance threshold between the UAV cluster and the mountain to be measured. When the distance between any UAV in the UAV cluster and the mountain to be measured is less than or equal to the distance threshold, the UAV cluster stops collecting images of the mountain to be measured and returns to the starting point, and executes step S4.

[0100] S4. Convert the image of the mountain to be measured obtained in step S3 into a grayscale image, extract the feature points of the mountain to be measured by using the Harris corner detection algorithm, and mark the extracted feature points as mountaintop feature points, mountain depression feature points and mountain foot feature points respectively according to the gradient change direction. Calculate the depth of each feature point from the corresponding UAV according to the camera imaging principle, and obtain the average coordinates of each feature point in the world coordinate system by using the transformation matrix.

[0101] This step specifically includes the following sub-steps:

[0102] S41. Based on the internal parameter matrix and distortion coefficients obtained in camera calibration, perform distortion correction on all the collected images to eliminate the image distortion problem caused by the camera lens. The color images after distortion correction need to be further pre-processed. First, use the median filtering technique to filter the color images, which can better protect the edge characteristics of the images while eliminating noise. Second, use the maximum value method to grayscale the color images, which can effectively convert the color images into grayscale images while reducing noise and improving image quality, and reducing the complexity of image processing.

[0103] S42. Define the key feature positions of the mountain to be measured, including "the top of the mountain", "the mountain depression" and "the foot of the mountain", as the feature points of the mountain to be measured. In the two-dimensional image, these feature points are manifested as the intersection points of two mountain ranges and belong to the boundaries of two mountain ranges in different directions in the local area. In the grayscale image, the gradient values and the rate of change of the gradient direction of these points are extremely significant. Apply the Harris corner detection algorithm to perform image segmentation on the grayscale image transformed in step S41, and perform feature point detection on the segmented grayscale image to obtain all the feature points in the grayscale image. Then, according to the gradient change direction, mark the identified feature points as the top-of-the-mountain feature points , the mountain-depression feature points , the foot-of-the-mountain feature points and .

[0104] Among them, the top-of-the-mountain feature points are the highest peak points of the local gradient change. The direction of the gradient radiates out from the top-of-the-mountain feature points and points to the downward direction of the mountain slope; the mountain-depression feature points are the lowest valley points of the local gradient change. The direction of the gradient converges inward and points to the mountain-depression feature points along the downward direction of the mountain slope; the foot-of-the-mountain feature points are the critical value points of the local gradient change. The direction of the gradient is arranged along the mountain slope direction. Among them, when the gradient vector points to the lower left , and when the gradient vector points to the lower right . The number of each feature point is determined by the detection of the grayscale image by the Harris corner detection algorithm.

[0105] S43. Arbitrarily select two grayscale images collected by two different UAV cameras at the same moment, and use the stereo matching algorithm based on the Census transform to perform stereo matching on the same feature points in the two grayscale images, and obtain the disparity of the same feature points in the two grayscale images.

[0106] Specifically, in the kth collection, arbitrarily select two images, use the stereo matching algorithm based on the Census transform to complete the stereo matching of each feature point in the images, and obtain the disparity of each point. Take the top-of-the-mountain feature points For example, arbitrarily select two cameras i and j from N cameras (it is stipulated that for the two selected cameras, the exploration UAV carrying camera i should be located to the left of the exploration UAV carrying camera j, that is, i < j). Use the stereo matching algorithm based on Census transform to match the grayscale images collected by the two cameras at the kth time, and select the mountaintop feature points in the image obtained by the left camera i The corresponding pixel , by calculating the Hamming distance between the Census value of each feature point pixel p in the image obtained by the right camera j and the Census value of the pixel , the feature point pixel with the smallest Hamming distance from the pixel is obtained , realizing the stereo matching of the mountaintop feature points in the pictures collected by cameras i and j at the kth time , and then obtaining the pixel points of the mountaintop feature points on the two images , The parallax in the horizontal direction , the transformation formula function of Census is:

[0107] ,

[0108] In the formula: is the Census value of the pixel point corresponding to the feature point in the image collected by the camera carried by exploration UAV i at the kth time; is the bit-by-bit concatenation operation of bits; is the neighborhood window of the feature point corresponding to the pixel point ; q is the pixel in the neighborhood window; , are the grayscale values of the pixel corresponding to the feature point and the neighborhood pixel q respectively; is the comparison operation;

[0109] The Hamming distance formula is:

[0110] ,

[0111] In the formula: is the Hamming distance between the pixels and ; is the Census value of the pixel point corresponding to the feature point in the image collected by the camera carried by exploration UAV j at the kth time; E is the number of bits of the Census transform; is the exclusive OR operation; Among the two grayscale images for matching as specified, the exploration UAV i should be located to the left of the exploration UAV j; is a comparison operation;

[0112] The parallax formula is:

[0113] ,

[0114] In the formula: are the matching pixel points of the feature point in the images acquired by the cameras carried by the exploration UAVs i and j in the k-th acquisition, and is the parallax; , are respectively the horizontal pixel coordinates of the matching points in the image pixel coordinate system , .

[0115] S44. Using the parallax of the mountaintop feature points obtained in the two grayscale images in step S43, according to the camera imaging principle and the law of similar triangles, calculate the depth of the mountaintop feature point from the camera carried by the UAV i in step S43, and then obtain the three-dimensional coordinates of the mountaintop feature point in this camera coordinate system. Then, use the camera parameters to convert the three-dimensional coordinates of the mountaintop feature point in the camera coordinate system to the coordinates in the world coordinate system.

[0116] The three-dimensional coordinates in the camera coordinate system are:

[0117] ,

[0118] In the formula: , , are respectively the x, y, and z direction coordinates of the mountaintop feature point obtained from the grayscale images acquired by the cameras i and j in the k-th acquisition in the camera coordinate system ; i and j are the numbers of the two UAVs respectively; f is the camera focal length; is the expected relative distance between the UAVs; , are respectively the x and y direction coordinates of the mountaintop feature point in the image pixel coordinate system of the camera i ; is the mountaintop feature point deduced according to the camera imaging principle and the law of similar triangles for the cameras carried by the exploration UAVs i and j in the images acquired in the k-th acquisition The depth from the camera carried by the exploration drone i; and are the principal point coordinates of camera i.

[0119] Using the rotation matrix and translation matrix of the camera to construct a transformation matrix, convert the three-dimensional coordinates of the mountaintop feature point in the camera coordinate system to the coordinates in the world coordinate system:

[0120] ,

[0121] where: is the coordinate of the mountaintop feature point in the k-th acquisition, obtained from the grayscale images collected by cameras i and j, in the world coordinate system; is the coordinate of the mountaintop feature point in the camera coordinate system obtained from the images collected by cameras i and j in the k-th acquisition; is the rotation matrix of camera i, is the translation matrix of camera i, obtained from the parameters of the camera.

[0122] S45. Repeat steps S43 and S44 to obtain the coordinate set of the mountaintop feature point in the world coordinate system obtained from the grayscale images collected by any two drones in the k-th acquisition:

[0123] .

[0124] To improve the accuracy of the measurement model, use neighborhood filtering to screen and eliminate the feature point coordinates with large errors, sum and average the coordinate set of the mountaintop feature points after elimination, and obtain the average coordinate of the mountaintop feature point in the world coordinate system .

[0125] S46. Repeat steps S43 to S45 to obtain the mountaintop feature point , the mountain depression feature point , the mountain foot feature point and the average coordinates in the world coordinate system.

[0126] S5. Perform multivariate polynomial fitting based on the coordinates of each feature point in the world coordinate system obtained in step S4 to obtain the ridge curve function of the mountain to be measured. Set the range for generating the ridge curve, extract the ridge curve segments, and connect the endpoints of each ridge curve segment to form a closed cross-sectional view of the mountain to be measured. Then, according to the inclination angle and inclination direction of the mountain to be measured obtained by the radar, perform coordinate transformation on the feature points on the cross-sectional view of the mountain to be measured, and sweep the transformed feature points along the inclination direction to obtain the three-dimensional model of the mountain to be measured.

[0127] This step specifically includes the following sub-steps:

[0128] S51. Use the average coordinates of each feature point in the world coordinate system for multivariate polynomial fitting. Since the number of feature points such as "mountain top", "mountain depression", and "mountain foot" of the mountain to be measured is limited, in order to avoid overfitting, a quadratic polynomial function is used as the fitting model for the ridge line of the mountain to be measured. This model can fully describe the distribution trend of each feature point in the world coordinate system. Obtain the optimal parameters of the quadratic polynomial through the least squares error method to construct the ridge curve function.

[0129] S52. Set the range for generating the ridge curve according to the average coordinates of the mountain top feature point and the mountain foot feature point in the world coordinate system:

[0130] ,

[0131] In the formula: , are the x-axis and y-axis coordinates of the mountain to be measured in the world coordinate system respectively; is the average coordinate in the x-axis direction of the mountain foot feature point 1 in the world coordinate system, is the average coordinate in the x-axis direction of the mountain foot feature point 2 in the world coordinate system, is the average coordinate in the y-axis direction of the mountain top feature point i in the world coordinate system; S is the number of mountain top feature points.

[0132] S53. According to the ridge curve function and the range for generating the ridge curve, extract the ridge curve segments and connect the endpoints of the ridge curve segments to obtain a closed cross-sectional view of the mountain to be measured.

[0133] S54. According to the inclination angle and inclination direction of the mountain to be measured obtained by the radar, establish a reduction matrix and a rotation matrix about the y-axis to transform the feature points on the cross-sectional view of the mountain to be measured. The rotation matrix and the reduction matrix are as follows:

[0134] ,

[0135] ,

[0136] In the formula: is the rotation matrix about the y-axis; is the inclination angle of the mountain to be measured obtained by the radar; is the reduction matrix; , and are the reduction ratios in the x, y, and z-axis directions respectively;

[0137] S55. Sweep the feature points on the cross-sectional view of the mountain to be measured after transformation along the inclination direction to obtain the three-dimensional model of the mountain to be measured.

[0138] S56. There are some points with unnatural transitions at the edge of the three-dimensional model of the mountain to be measured obtained by sweeping. Fix the positions of the feature points at the mountain top , the feature points in the mountain depression , and the feature points at the mountain foot and in the model. By detecting the edge of the three-dimensional model of the mountain to be measured, use Laplacian smoothing to reduce the unnatural transition at the edge of the three-dimensional model of the mountain to be measured. The formula is:

[0139] ,

[0140] In the formula: is the position of the point to be smoothed; is the position of the point after smoothing; is the set of adjacent points of point ; is the size of the set .

[0141] Output the three-dimensional model of the mountain to be measured in STL format after smoothing.

[0142] S6. Use the feature points at the mountain top and the feature points in the mountain depression as task detection points, calculate the Manhattan distance between each drone and each task detection point, and use the greedy algorithm to perform detection task allocation with the goal of minimizing the total detection distance, so that each task detection point is assigned to a drone, and the drones not assigned to the detection task stay at the starting point.

[0143] This step specifically includes the following sub-steps:

[0144] S61. Establish a set of N drones :

[0145] ,

[0146] Set of task detection points:

[0147] ,

[0148] In the formula: is the set of task detection points; is the mean coordinate of the top feature points of the mountain to be measured, is the mean coordinate of the characteristic points of the depressions in the mountain to be measured.

[0149] S62. Calculate the Manhattan distance between each drone and each detection target point:

[0150] ,

[0151] Where: For drones Task Checkpoint The Manhattan distance between and For drones and task checkpoints In the world coordinate system The location coordinates of the following.

[0152] S63. Taking the minimum total detection distance as the performance indicator, according to the greedy strategy, select one from the unassigned task detection points, and preferentially assign the selected task detection point to the drone with the shortest Manhattan distance from the point, and add the drone to the set of assigned drones. Repeat the assignment steps until all tasks are assigned.

[0153] S64. UAVs that are not assigned to inspection tasks remain at the starting point and wait.

[0154] S7. Using the three-dimensional model of the mountain to be tested obtained in step S5, selecting Manhattan distance as the evaluation function, and using the A-Star algorithm (A* algorithm) to perform three-dimensional path planning for each UAV assigned to the detection task in step S6.

[0155] This step specifically includes the following sub-steps:

[0156] S71, rasterizing the three-dimensional space of a certain UAV assigned to the detection task, and determining the initial position coordinates and task detection point coordinates of the UAV;

[0157] S72, using Manhattan distance as the evaluation function of the A-Star algorithm, starting from the initial position coordinates of the UAV, selecting the node with the minimum evaluation function among its adjacent nodes as the next central node, and repeating this selection method until the task detection point coordinates are reached;

[0158] S73. Perform reverse tracing from the task detection point coordinates to construct the flight path between the initial position coordinates of the UAV and the task detection point coordinates; use a third-order Bézier curve to smooth the flight path. Take the nodes on the flight path obtained by the A-Star algorithm as the control points of the third-order Bézier curve, and take 4 consecutive nodes as a group of control points. By adjusting the control point coordinates, change the shape of the flight path to make the output flight path smoother, ensure the stability of the detection UAV during flight, reduce energy consumption, and avoid problems such as position deviation caused by too large turning angles. The formula for the third-order Bézier curve is as follows:

[0159] ,

[0160] In the formula: is the output point of the third-order Bézier curve; , , and are 4 consecutive Bézier control points; is the variation parameter.

[0161] S74. Repeat steps S71 to S73 until all UAVs assigned with detection tasks have completed the corresponding complete path planning.

[0162] This method can be applied to the monitoring of various mountain natural disasters. For example, it can be applied to the monitoring of landslides, and combined with a machine learning model to predict and give early warnings for landslides. The following is a brief description taking the landslide monitoring as an example.

[0163] S8. The UAV executes the detection of the mountain to be measured along the optimized flight path in step S7, and uses the equipped camera to collect image data of the mountain. During the flight process, the UAV real-time monitors its own position and speed, and compensates for the position, attitude and speed deviations caused by environmental factors through control and adjustment. For the needs of later image stitching, the overlap rate between adjacent images on the same flight line needs to reach 25%, and the camera shooting speed needs to match the flight speed of the detection UAV to ensure full coverage and repetition rate of the images. During the task process, the collected image data is transmitted in real time to the central processing unit in the cabin of the UAV cluster through the image transmission module. When the UAV flies to the end point of the flight path, the camera stops collecting data of the mountain to be measured, and the UAV automatically returns to the cabin along the flight path.

[0164] The central processing unit in the cabin of the UAV cluster performs preprocessing such as cropping and image enhancement on the received images, and uses the pre-trained landslide crack recognition model to perform crack feature recognition and detection on the preprocessed images. If the collected images contain preset crack features, it is determined that there are precursors of landslides in the mountain, and enter step S10. If there are no preset crack features, continue to execute step S9.

[0165] S10. Stitch the remaining collected images on the route where the crack is located, mark the area where the crack is located, and promptly issue a disaster warning to the outside.

Claims

1. A multi-UAV collaborative path planning method for mountain monitoring, characterized in that: The following steps are involved: S1. Select several drones with the same structure to form a drone cluster. Each drone is equipped with a camera and a radar. The control signal of each drone in the drone cluster is generated based on the set drone kinematic model. S2. Construct a drone cluster formation based on the consistency control protocol, set a drone in the drone cluster as the pilot drone, and set the remaining drones as follower drones. When flying, the follower drones are arranged in a fixed order on the same side of the pilot drone; S3, the drone cluster formation set in control step S2 starts from the starting point, flies at a set speed and approaches the mountain to be measured, uses the camera carried by each drone to collect images of the mountain to be measured at preset time intervals, and marks the corresponding position information on each image of the mountain to be measured, uses the radar carried by each drone to measure the distance between each drone and the mountain to be measured in real time, and obtains the inclination angle and inclination direction of the mountain to be measured; by setting a distance threshold between the drone cluster and the mountain to be measured, when the distance between any drone in the drone cluster and the mountain to be measured is less than or equal to the distance threshold, the drone cluster stops collecting images of the mountain to be measured and returns to the starting point; S4, converting the image of the mountain to be measured obtained in step S3 into a grayscale image, extracting feature points of the mountain to be measured using the Harris corner detection algorithm, and marking the extracted feature points as mountain top feature points, valley feature points and mountain foot feature points according to the direction of gradient change, calculating the depth of each feature point to the corresponding drone according to the camera imaging principle, and obtaining the average coordinates of each feature point in the world coordinate system using the conversion matrix; S5, performing multivariate polynomial fitting according to the coordinates of each feature point obtained in step S4 in the world coordinate system to obtain a ridge curve function of the mountain to be measured, setting a ridge curve generation range, extracting ridge curve segments and connecting the endpoints of each ridge curve segment to form a closed profile of the mountain to be measured, and then performing coordinate transformation on the feature points on the profile of the mountain to be measured according to the inclination angle and inclination direction of the mountain to be measured obtained by the radar, and scanning the transformed feature points along the inclination direction to obtain a three-dimensional model of the mountain to be measured; S6. Take the mountain top feature points and valley feature points as task detection points, calculate the Manhattan distance between each UAV and each task detection point, take the minimum total detection distance as the goal, use the greedy algorithm to allocate detection tasks, so that each task detection point is assigned to a UAV, and the UAVs that are not assigned to the detection task stay at the starting point; S7. Using the three-dimensional model of the mountain to be tested obtained in step S5, selecting Manhattan distance as the evaluation function, and using the A-Star algorithm to perform three-dimensional path planning for each UAV assigned to the detection task in step S6.

2. The multi-UAV collaborative path planning method according to claim 1, characterized in that: The step S1 comprises: The second-order consistent control input of UAV i is obtained by using the second-order motion equation of UAV i. : ， Where: , , They are respectively UAV i in the world coordinate system The position coordinates in the x, y, and z directions, , , are the velocity vectors of drone i in the x, y, and z directions, respectively. , , are the second-order consistency control inputs of UAV i Components in the x, y, and z directions.

3. The multi-UAV collaborative path planning method according to claim 1, characterized in that: The step S2 comprises the following sub-steps: S21, calculating the relative distance between each follower drone and the pilot drone; S22, the pilot drone and each follower drone take off from the starting point in sequence and fly to a preset altitude, and then maintain the preset altitude; S23, setting the speed of each follower drone to be consistent with that of the pilot drone, and controlling the pilot drone to fly on the far left of the formation, and the follower drones are arranged in sequence on the right side of the pilot drone at a fixed distance in the order of numbers. According to the consistency control protocol with spacing in three-dimensional space, the control input of the follower drone is: ， Where: is the control input of the following drone i at time t; is a fixed gain; is the damping gain; is the connection strength between following drone i and following drone j; i and j are drone numbers; is the relative distance between the following UAV i and the pilot UAV at time t; is the relative distance between the following UAV j and the pilot UAV at time t; is the expected relative distance between any two UAVs; The world coordinate system of the following drone i and the pilot drone t at time The position distance on the lower z axis; The world coordinate system of the following drone j and the pilot drone t at all times The position distance on the lower z axis; , are the velocity amplitudes of following UAV i and following UAV j at time t respectively; is the total number of drones in the drone swarm; S24. Under the control input of the following drone, the following drone tracks the desired position by tracking the desired relative distance until the speed and position state of the drone cluster meets the consistency control protocol with spacing in three-dimensional space, thus completing the drone cluster formation.

4. The multi-UAV collaborative path planning method according to claim 1, characterized in that: The step S4 comprises the following sub-steps: S41, converting the image of the mountain to be measured into a grayscale image using a maximum value method; S42, applying the Harris corner detection algorithm to perform image segmentation on the grayscale image converted in step S41, and performing feature point detection on the segmented grayscale image to obtain all feature points in the grayscale image, and then marking the identified feature points as mountain top feature points, valley feature points and foot feature points according to the gradient change direction; S43, randomly selecting two grayscale images captured by two different drone cameras at the same time, using a stereo matching algorithm based on Census transformation to perform stereo matching on the same feature point in the two grayscale images, and obtaining the disparity of the same feature point in the two grayscale images; S44, using the disparity of the same feature point in the two grayscale images obtained in step S43, according to the camera imaging principle and the triangle similarity law, calculate the depth of the feature point from the camera carried by one of the two drones in step S43, and then obtain the three-dimensional coordinates of the feature point in the camera coordinate system, and then use the camera parameters to convert the three-dimensional coordinates of the feature point in the camera coordinate system into coordinates in the world coordinate system; S45, repeating steps S43 and S44, obtaining a coordinate set of the same feature point in the world coordinate system obtained by grayscale images collected by any two drones at each collection time, and calculating the average coordinate of the feature point in the world coordinate system; S46. Repeat steps S43 to S45 to obtain the average coordinates of all feature points in the world coordinate system.

5. The multi-UAV collaborative path planning method according to claim 4, characterized in that: The step S43 includes: using a stereo matching algorithm based on Census transformation to match the grayscale images collected by the two cameras, selecting a pixel corresponding to a feature point in a grayscale image , calculate the Census value and pixel value of each feature point in another grayscale image The Hamming distance of the Census value is obtained with the pixel The feature point pixel with the smallest Hamming distance , realize stereo matching of the same feature point in the grayscale images collected by two cameras, and then obtain the pixel point of the feature point on the two grayscale images , The parallax in the horizontal direction is given by the following formula: ， Where: To explore the feature points of the camera carried by drone i in the kth acquisition The corresponding pixel Census value; It is a bit-by-bit connection operation; For feature points Corresponding pixel Neighborhood window; q is the pixel in the neighborhood window; , The feature points Corresponding pixels and the gray value of the neighborhood pixel q; For comparison operations; ， Where: Pixel and The Hamming distance of To explore the feature points of the camera carried by UAV j in the kth acquisition The corresponding pixel Census value; E is the number of digits of Census transformation; is an XOR operation; It is stipulated that in the two grayscale images to be matched, the exploration drone i should be located on the left side of the exploration drone j; For comparison operations; ， Where: To explore the feature points in the image acquired by the camera carried by UAV i and j for the kth time The matching pixels and parallax; , The image pixel coordinate system matching point under , The horizontal pixel coordinate of .

6. The multi-UAV collaborative path planning method according to claim 5, characterized in that: The step S44 comprises: According to the camera imaging principle and the triangle similarity law, the depth of the feature point from the camera carried by the exploration drone i and its three-dimensional coordinates in the camera coordinate system are calculated. The formula is: ， Where: , , The feature points The grayscale image acquired by cameras i and j in the kth acquisition is obtained in the camera coordinate system The x, y, and z coordinates of the drone. i and j are the numbers of the two drones respectively. f is the camera focal length; is the expected relative distance between each UAV; , The feature points At camera i in the image pixel coordinate system The x and y coordinates of the following: To explore the feature points in the image acquired by the camera carried by UAV i and j for the kth time according to the camera imaging principle and the triangle similarity law The depth of the camera carried by the exploration drone i; , is the principal point coordinate of camera i; The camera's rotation matrix and translation matrix are used to construct a transformation matrix to transform the three-dimensional coordinates of the feature points in the camera coordinate system into coordinates in the world coordinate system: ， Where: For feature points In the kth acquisition, the grayscale image acquired by cameras i and j is obtained in the world coordinate system The coordinates below; For feature points The image acquired by cameras i and j in the kth acquisition is obtained in the camera coordinate system The coordinates below; is the rotation matrix of camera i, is the translation matrix of camera i, obtained from the camera parameters.

7. The multi-UAV collaborative path planning method according to claim 1, characterized in that: The step S5 comprises the following sub-steps: S51, performing multivariate polynomial fitting using the average coordinates of each feature point in the world coordinate system, obtaining optimal parameters of the multivariate polynomial by minimizing the error, and constructing a ridge curve function; S52, setting the range of generating the ridge curve according to the average coordinates of the feature points at the top of the mountain and the feature points at the foot of the mountain in the world coordinate system: ， Where: , They are the x-axis and y-axis coordinates of the mountain to be measured in the world coordinate system; is the average coordinate of the feature point 1 at the foot of the mountain in the x-axis direction in the world coordinate system, is the average coordinate of the feature point 2 at the foot of the mountain in the x-axis direction in the world coordinate system, is the average coordinate of the feature point i on the top of the mountain in the y-axis direction in the world coordinate system; S is the number of feature points on the top of the mountain; S53, extracting ridge curve segments according to the ridge curve function and the ridge curve generation range, and connecting the end points of the ridge curve segments to obtain a closed profile of the mountain to be measured; S54, according to the inclination angle and inclination direction of the mountain to be measured obtained by the radar, a reduction matrix and a rotation matrix about the y-axis are established to transform the feature points on the profile of the mountain to be measured. The rotation matrix and the reduction matrix are as follows: ， ， Where: is the rotation matrix about the y-axis; The inclination angle of the mountain to be measured is obtained by the radar; To reduce the matrix; , and They are the reduction ratios in the x, y, and z axis directions respectively; S55, scanning the characteristic points on the transformed profile of the mountain to be measured along the inclined direction to obtain a three-dimensional model of the mountain to be measured.

8. The multi-UAV collaborative path planning method according to claim 7, characterized in that: The step S5 further comprises: S56. Use Laplace smoothing to reduce the unnatural transition of the edge of the three-dimensional model of the mountain to be measured. The formula is: ， Where: is the point position that needs to be smoothed; is the point position after smoothing; For point The set of adjacent points of ; For collection size.

9. The multi-UAV collaborative path planning method according to claim 1, characterized in that: The step S7 comprises the following sub-steps: S71, rasterizing the three-dimensional space of a certain UAV assigned to the detection task, and determining the initial position coordinates and task detection point coordinates of the UAV; S72, using Manhattan distance as the evaluation function of the A-Star algorithm, starting from the initial position coordinates of the UAV, selecting the node with the minimum evaluation function among its adjacent nodes as the next central node, and repeating this selection method until the task detection point coordinates are reached; S73, tracing back from the coordinates of the mission detection point to construct a track path from the initial position coordinates of the UAV to the coordinates of the mission detection point; S74. Repeat steps S71 to S73 until all drones assigned to the detection task have completed the corresponding complete path planning.

10. The multi-UAV collaborative path planning method according to claim 9, characterized in that: The step S73 further includes: The track path is smoothed using a third-order Bezier curve. Each node on the track path obtained by the A-Star algorithm is used as the control point of the third-order Bezier curve. Four consecutive nodes are used as a group of control points. The formula of the third-order Bezier curve is as follows: ， Where: is the output point of the third-order Bezier curve; , , and are 4 consecutive Bezier control points; is the changing parameter.

Images