A method for rapid extraction of vegetation canopy height based on UAV RGB camera

Abstract

This invention specifically refers to a method for rapidly extracting vegetation canopy height based on an RGB camera from a UAV. It utilizes orthophotos collected by the UAV flying along a set flight path, employs a motion structure algorithm to generate a 3D dense point cloud, and combines an aggressive two-step filtering method with a binary image obtained from Otsu segmentation (ExG-ExR) to achieve ground point filtering. The ground points are then processed into a regularized grid. Based on the fundamental idea of ​​least-squares quadratic fitting and prediction of profile lines, rapid interpolation is performed on the gridded ground points along both rows and columns to obtain two interpolated ground models. The average of the two ground models yields the final ground model. The original dense point cloud is then gridded at the same resolution. The vegetation canopy height model is obtained by subtracting the final ground model from the original gridded point cloud. This invention simplifies the parameter selection difficulty of the two-step filtering by utilizing the RGB information of the point cloud in ground point filtering. The proposed interpolation algorithm is less time-consuming than traditional algorithms while achieving similar accuracy.

Description

Technical Field

[0001] This invention relates to the field of ecological environment vegetation canopy extraction technology, specifically to a method for rapid extraction of vegetation canopy height based on a drone RGB camera. Background Technology

[0002] Vegetation, as a crucial component of ecosystems and a primary producer, is most closely linked to other elements (such as soil water and groundwater), playing a vital role in mitigating desertification and protecting biodiversity. Vegetation canopy height is one of the most commonly used indicators for assessing regional vegetation biomass and health.

[0003] Drone-based remote sensing overcomes the inefficiency of traditional manual measurement. Currently, the main method for obtaining vegetation height using drones is LiDAR (Light Detection and Ranging) technology, leading to a focus on LiDAR data for ground filtering and interpolation algorithms. In recent years, simultaneous improvements in computer hardware performance and image matching software have enabled faster and more accurate dense matching of stereo images using moving structures. These latest technological advancements have made photogrammetry a competitive alternative to LiDAR.

[0004] Compared to LiDAR, the 3D point cloud format generated by SFM is similar to LiDAR data, but due to its production principle, there are still differences between the two. These differences are mainly reflected in: the inability to penetrate dense vegetation canopies; fewer ground points in vegetated areas with complex canopy structures; and a lower ability to distinguish low-lying vegetation compared to LiDAR; the number of point clouds is tens or even hundreds of times greater than that of LiDAR point clouds, making it a truly "dense point cloud"; and the point clouds contain RGB information. Current ground filtering algorithms mostly rely on precise structural threshold parameter settings, without considering the RGB information of the point clouds; when non-ground points are not completely removed, the remaining non-ground points often cause traditional interpolation algorithms based on certain neighborhood ground points to overestimate the surface of vegetated areas, resulting in an underestimation of canopy height; efficient interpolation algorithms have not yet been introduced, and computational complexity is a common problem in interpolation algorithms, making traditional algorithms too inefficient, especially when it comes to quickly obtaining regional canopy heights. Summary of the Invention

[0005] This invention addresses the shortcomings of existing technologies by providing a method for rapid extraction of vegetation canopy height based on a drone RGB camera.

[0006] A method for rapid extraction of vegetation canopy height based on a drone RGB camera includes the following steps:

[0007] Step 1: Based on the image data acquired by the UAV, the motion structure method is used to process and generate a 3D dense point cloud with RGB band information;

[0008] Step 2: Perform ground point filtering on the dense point cloud based on the aggressive two-step filtering, and calculate the point cloud (ExG-ExR) index at the same time. Obtain ground points based on the Otsu segmentation (ExG-ExR) index. Take the intersection of the ground points obtained by the aggressive two-step filtering and Otsu to obtain the ground point cloud.

[0009] Step 3: Based on the ground point cloud obtained in Step 2, perform regular meshing;

[0010] Step 4: Based on the regular gridded ground point cloud in Step 3, extract profile lines according to the row and column directions, perform least squares quadratic fitting and prediction on each profile line to fill in all ground holes, and obtain two surface models filled in according to the row and column directions. Perform a weighted average on the two surface models to obtain a weighted average ground model.

[0011] Step 5: Perform the same regular meshing on the 3D dense point cloud obtained in Step 1 as in Step 3 to obtain the original surface model. Subtract the weighted average ground model obtained in Step 4 from the original surface model to obtain the canopy structure model.

[0012] Furthermore, in step 1, the method for the UAV to acquire image data is as follows: based on the shape of the vegetation canopy height zone, two sets of mutually perpendicular flight strips are preset to cover the entire area, with the overlap in both the flight direction and the lateral direction reaching 80%.

[0013] In step 1, the method for generating a three-dimensional dense point cloud is as follows: feature point extraction and matching are performed based on the SIFT operator and the sample consistency algorithm RANSC, the bundle adjustment is used to optimize the reconstruction results, and finally the dense point cloud is obtained based on the three-dimensional multi-view stereo vision algorithm PMVS based on the patch.

[0014] Further, (i) the original surface model obtained in step 5 is obtained by: determining the initial filter window size based on the maximum plant canopy size in the vegetation canopy height area, applying it to the three-dimensional dense point cloud generated in step 1 to identify the lowest point in each window, classifying it into several initial ground points, and then performing linear interpolation on several initial ground points to obtain the initial terrain model.

[0015] (ii) Evaluate remaining unclassified points: Points that meet the following two conditions can be added to the ground category.

[0016] (a) The distance from the remaining unclassified points to the initial terrain surface does not exceed the maximum distance threshold;

[0017] (b) The angle difference between the nearest ground point and the initial terrain surface, and the angle difference between the nearest initial ground point and the point being evaluated, are less than the maximum angle threshold;

[0018] Among them, the filtering window, maximum distance threshold and maximum angle threshold are gradually reduced using an aggressive strategy to achieve step-by-step removal of vegetation from large to small. For complex scenes with mixed trees, shrubs and grasses, three aggressive filters are used to obtain ground point clouds.

[0019] Secondly, the (ExG-ExR) index of the three-dimensional dense point cloud generated in step 1 is calculated using equation (1). Based on the Otsu automated segmentation (ExG-ExR) index, a binary map of the ground classification result is obtained. The intersection of this map with the ground point cloud obtained by progressive filtering is used to obtain the ground point cloud described in step 2.

[0020] ExG-ExR=3*G-2.4*RB (1).

[0021] Furthermore, the specific method for obtaining the weighted average ground model in step 4 includes:

[0022] (i) Extract all profile lines of the regular meshed point cloud from step 3 by row and column respectively;

[0023] (ii) Detect the number and length of void areas for each profile line, and sort them from largest to smallest void length;

[0024] (III) Calculate the length L of the region with the largest cavity length. If it is even, add one. Then take 3(L-1) / 2 points from the left boundary to the left and 3(L-1) / 2 points from the right boundary to the right. Divide the 3(L-1) points into three segments: left, middle and right. Take the first 3(L-1) / 2 points from the left of the left segment and perform a least squares quadratic fitting. Predict (L-1) / 2 points to the right. Take the first 3(L-1) / 2 points from the right of the right segment and perform a least squares quadratic fitting. Predict (L-1) / 2 points to the left. Take (L-1) points from the middle segment and perform a least squares quadratic fitting. Predict the cavity length L points. The three segments overlap by L points.

[0025] (iv) Using y (i) (d1), y (i+1) (d2), d1, d2 = 1, 2, ..., L represent the fitting polynomials of the i-th segment and the (i+1)-th segment, respectively. The weighted average of the left-middle and middle-right overlapping regions is calculated using equation (2). The weight decreases linearly with the distance between the point and the center of the line segment, and effectively eliminates any jumps or discontinuities near the boundary of adjacent parts.

[0026]

[0027] (v) Replace the hollow areas with the weighted values ​​of the hollow areas obtained in step (iv), and repeat steps (ii) to (iv) until all hollow areas are eliminated, to obtain two terrain models extracted by row and column;

[0028] (vi) The final terrain model is obtained by weighted averaging of the two terrain models. Compared with the prior art, the present invention has the following advantages:

[0029] (1) This invention introduces the aggressive approach into two-step filtering, gradually reducing the threshold to achieve step-by-step removal of trees, shrubs and grasses. This avoids the misclassification caused by the heterogeneity of slope and vegetation type in complex scenarios when filtering with a fixed threshold. It can filter vegetation point clouds more effectively. The Otsu method is used to automatically classify the (ExG-ExR) index calculated by the RGB band of the point cloud and take the intersection with the two-step filtering results to obtain the final ground point cloud. This can effectively assist in the selection of two-step filtering parameters. Because the Otsu method is good at distinguishing low vegetation, this can avoid the loss of ground points caused by the overly strict threshold selection of the aggressive two-step filtering. The reduction of ground points will affect the accuracy of the interpolation algorithm.

[0030] (2) The least squares quadratic fitting and prediction algorithm based on profile lines proposed in this invention adaptively extracts ground points within a certain length based on the size of the missing window in one dimension and performs fitting and prediction in three segments. By adopting a reasonable weighting method, the local terrain trend is preserved. At the same time, the introduction of a trend term effectively suppresses the unavoidable misclassification points in the ground point filtering process. In addition, because the proposed algorithm has low computational complexity and is easy to implement in parallel, the running time is greatly shortened compared with traditional interpolation algorithms such as inverse distance weighted interpolation and Kriging interpolation. Attached Figure Description

[0031] Figure 1 This is a flowchart of a method for rapidly extracting vegetation canopy height based on an RGB camera from a drone, provided by the present invention. Detailed Implementation

[0032] The present invention will now be described in further detail with reference to the accompanying drawings.

[0033] The technical solution adopted in this invention is a method for rapid extraction of vegetation canopy height based on UAV RGB camera. This method mainly consists of two parts: first, an aggressive two-step filtering method based on point cloud RGB band information to extract ground points; second, addressing data gaps in the filtered ground point cloud, proposing a fast interpolation algorithm based on profile line least-squares quadratic fitting and prediction. This algorithm has low complexity and high parallelizability, while effectively suppressing the influence of misclassified points on the interpolation results. The method is implemented according to the following steps:

[0034] Combination Figure 1 The present invention will be further described as follows:

[0035] A method for rapid extraction of vegetation canopy height based on a drone RGB camera, the method comprising:

[0036] Step 1: Based on the image data acquired by the UAV, the motion structure method is used to process and generate a 3D dense point cloud with RGB band information;

[0037] Step 2: Perform ground point filtering on the dense point cloud based on the aggressive two-step filtering, and calculate the point cloud (ExG-ExR) index at the same time. Obtain ground points based on the Otsu segmentation (ExG-ExR) index. Take the intersection of the ground points obtained by the aggressive two-step filtering and Otsu to obtain the ground point cloud.

[0038] Step 3: Based on the ground point cloud obtained in Step 2, perform regular meshing;

[0039] Step 4: Based on the regular gridded ground point cloud in Step 3, extract profile lines according to the row and column directions, perform least squares quadratic fitting and prediction on each profile line to fill in all ground holes, and obtain two surface models filled in according to the row and column directions. Perform a weighted average on the two surface models to obtain a weighted average ground model.

[0040] Step 5: Perform the same regular meshing on the 3D dense point cloud obtained in Step 1 as in Step 3 to obtain the original surface model. Subtract the weighted average ground model obtained in Step 4 from the original surface model to obtain the canopy structure model.

[0041] In step 1, the flight route should be designed as two mutually perpendicular flight strips covering the entire study area, with a heading and lateral overlap of 80% to expand the camera's field of view, thereby increasing the number of ground points and improving the accuracy of aerial triangulation. If conditions permit, oblique photography should be introduced. The selection of flight altitude needs to be calculated based on camera parameters and the minimum ground resolution required for the operation. The time for data acquisition should be selected during periods of low light variation and low wind speed, such as 2 pm. Multiple flights should be scheduled at the same time. In addition, to ensure the reliability of the results, a certain number of ground control points need to be set up according to the accuracy requirements.

[0042] The UAV image sequence was used to extract and match feature points based on the SIFT operator and the sample consistency algorithm RANSC. The bundle adjustment was used to optimize the reconstruction results. Finally, a dense point cloud with RGB band information was obtained based on the patch-based three-dimensional multi-view stereo vision algorithm PMVS.

[0043] In step 2, the initial filter window size is determined based on the maximum canopy size of the plant in the study area. The window must be larger than the maximum canopy area to ensure that there is at least one ground point in each window. This is then applied to the point cloud generated in step 1 to identify the lowest point in each window and classify it as an initial ground point. Then, linear interpolation is performed on these initial ground points to obtain an approximate initial terrain model. Secondly, all remaining unclassified points are evaluated, and they are added to the ground category if they meet both of the following conditions: (i) their distance to the initial terrain surface does not exceed the maximum distance threshold, and (ii) the angle difference between the nearest ground point and the initial terrain surface and the angle difference between the nearest initial ground point and the evaluated point are less than the maximum angle threshold. Among them, the filtering window, maximum distance threshold and maximum angle threshold are gradually reduced using an aggressive strategy to achieve stepwise removal of vegetation from large to small. For complex scenarios with mixed trees, shrubs and grasses, it is recommended to use three aggressive filters to obtain the ground point cloud. Taking the mixed tree, shrub and grass area in the semi-arid region as an example, the parameters of the three aggressive filters are: window size: 10m, 2m, 0.5m, maximum distance threshold: 0.6m, 0.3m, 0.1m, maximum angle threshold: 7°, 3°, 1°.

[0044] The (ExG-ExR) index of the point cloud generated in step 1 is calculated using equation (1). The ground classification result binary map is obtained based on Otsu automatic segmentation (ExG-ExR). The intersection of this map with the ground point cloud obtained by progressive two-step filtering is used to obtain the final ground point cloud.

[0045] ExG-ExR=3*G-2.4*RB (1).

[0046] In step 3, the selection of grid size in the process of regularizing ground point clouds follows the principle of effectively identifying the research object. For example, if the main research is on the height of shrub canopy, a 5cm grid can be selected. It is recommended that the grid be slightly larger than the ground resolution of the UAV image. The point cloud located in the grid takes its average or maximum elevation value as the grid value. If there are no points in the grid, it is set to a null value.

[0047] The specific steps for step 4 are as follows:

[0048] (i) Extract all profile lines of the regular meshed point cloud from step 3 by row and column respectively;

[0049] (ii) Detect the number and length of void areas for each profile line, and sort them from largest to smallest void length;

[0050] (III) Calculate the length L of the region with the largest cavity length. If it is even, add one. Then take 3(L-1) / 2 points from the left boundary to the left and 3(L-1) / 2 points from the right boundary to the right. Divide the 3(L-1) points into three segments: left, middle and right. Take the first 3(L-1) / 2 points from the left of the left segment and perform a least squares quadratic fitting. Predict (L-1) / 2 points to the right. Take the first 3(L-1) / 2 points from the right of the right segment and perform a least squares quadratic fitting. Predict (L-1) / 2 points to the left. Take (L-1) points from the middle segment and perform a least squares quadratic fitting. Predict the cavity length L points. The three segments overlap by L points.

[0051] (iv) Using y (i) (d1), y (i+1) (d2), d1, d2 = 1, 2, ..., L represent the fitting polynomials of the i-th segment and the (i+1)-th segment, respectively. The weighted average of the left-middle and middle-right overlapping regions is calculated using equation (2). The weight decreases linearly with the distance between the point and the center of the line segment, and effectively eliminates any jumps or discontinuities near the boundary of adjacent parts.

[0052]

[0053] (v) Replace the hollow areas with the weighted values ​​of the hollow areas obtained in step (iv), and repeat steps (ii) to (iv) until all hollow areas are eliminated, to obtain two terrain models extracted by row and column;

[0054] (vi) The final terrain model is obtained by weighted averaging of the two terrain models.

[0055] Finally, it should be noted that in step 4, taking the profile lines of the gridded point cloud as the basic research object according to the rows and columns can be understood as a four-directional interpolation of ground points. If a more accurate and smooth vegetation canopy model is desired, technicians in the relevant field can introduce the two diagonal directions of the gridded point cloud to achieve eight-directional interpolation, or even introduce sixteen or thirty-two directions. Finally, after removing error points from the interpolation results of multiple directions in the missing area, the average is taken as the final interpolation result. This can ensure that the results are more reliable and restore the hydrological continuity to a certain extent.

[0056] The present invention and its embodiments have been described above. This description is not restrictive, and the accompanying drawings are only one embodiment of the present invention; the actual structure is not limited thereto. In conclusion, if those skilled in the art are inspired by this description and design methods and embodiments similar to this technical solution without departing from the spirit of the present invention, such methods and embodiments should fall within the protection scope of the present invention.

Claims

1. A method for rapid extraction of vegetation canopy height based on a drone RGB camera, characterized in that, Includes the following steps: Step 1: Based on the image data acquired by the UAV, the motion structure method is used to process and generate a 3D dense point cloud with RGB band information; Step 2: Perform ground point filtering on the dense point cloud based on the aggressive two-step filtering, and calculate the point cloud (ExG-ExR) index at the same time. Obtain ground points based on the Otsu segmentation (ExG-ExR) index. Take the intersection of the ground points obtained by the aggressive two-step filtering and Otsu to obtain the ground point cloud. Step 3: Based on the ground point cloud obtained in Step 2, perform regular meshing; Step 4: Based on the regular gridded ground point cloud in Step 3, extract profile lines according to the row and column directions, perform least squares quadratic fitting and prediction on each profile line to fill in all ground holes, and obtain two surface models filled in according to the row and column directions. Perform a weighted average on the two surface models to obtain a weighted average ground model. Step 5: Perform the same regular meshing on the 3D dense point cloud obtained in Step 1 as in Step 3 to obtain the original surface model. Subtract the weighted average ground model obtained in Step 4 from the original surface model to obtain the canopy structure model.

2. The method for rapid extraction of vegetation canopy height based on a UAV RGB camera according to claim 1, characterized in that, In step 1, the method for the UAV to acquire image data is as follows: based on the shape of the vegetation canopy height zone, two sets of mutually perpendicular flight strips are preset to cover the entire area, with the overlap in both the flight direction and the lateral direction reaching 80%. In step 1, the method for generating a three-dimensional dense point cloud is as follows: feature point extraction and matching are performed based on the SIFT operator and the sample consistency algorithm RANSC, the bundle adjustment is used to optimize the reconstruction results, and finally the dense point cloud is obtained based on the three-dimensional multi-view stereo vision algorithm PMVS based on the patch.

3. The method for rapid extraction of vegetation canopy height based on a UAV RGB camera according to claim 1, characterized in that: (i) Obtaining the original surface model in step 5: Based on the size of the largest plant canopy in the vegetation canopy height area, determine the size of the initial filtering window, apply it to the three-dimensional dense point cloud generated in step 1 to identify the lowest point in each window, and classify it into several initial ground points. Then, perform linear interpolation on several initial ground points to obtain the initial terrain model. (ii) Evaluate remaining unclassified points: Points that meet the following two conditions can be added to the ground category. (a) The distance from the remaining unclassified points to the initial terrain surface does not exceed the maximum distance threshold; (b) The angle difference between the nearest ground point and the initial terrain surface, and the angle difference between the nearest initial ground point and the point being evaluated, are less than the maximum angle threshold; Among them, the filtering window, maximum distance threshold and maximum angle threshold are gradually reduced using an aggressive strategy to achieve step-by-step removal of vegetation from large to small. For complex scenes with mixed trees, shrubs and grasses, three aggressive filters are used to obtain ground point clouds. Secondly, the (ExG-ExR) index of the three-dimensional dense point cloud generated in step 1 is calculated using equation (1). Based on the Otsu automated segmentation (ExG-ExR) index, a binary map of the ground classification result is obtained. The intersection of this map with the ground point cloud obtained by progressive filtering is used to obtain the ground point cloud described in step 2. ExG-ExR=3*G-2.4*RB (1).

4. The method for rapid extraction of vegetation canopy height based on a UAV RGB camera according to claim 1, characterized in that, The specific method for obtaining the weighted average ground model in step 4 includes: (i) Extract all profile lines of the regular meshed point cloud from step 3 by row and column respectively; (ii) Detect the number and length of void areas for each profile line, and sort them from largest to smallest void length; (III) Calculate the length L of the region with the largest cavity length. If it is even, add one. Then take 3(L-1) / 2 points from the left boundary to the left and 3(L-1) / 2 points from the right boundary to the right. Divide the 3(L-1) points into three segments: left, middle and right. Take the first 3(L-1) / 2 points from the left of the left segment and perform a least squares quadratic fitting. Predict (L-1) / 2 points to the right. Take the first 3(L-1) / 2 points from the right of the right segment and perform a least squares quadratic fitting. Predict (L-1) / 2 points to the left. Take (L-1) points from the middle segment and perform a least squares quadratic fitting. Predict the cavity length L points. The three segments overlap by L points. (iv) Using y (i) (d1), y (i+1) (d2), d1, d2 = 1, 2, ..., L represent the fitting polynomials of the i-th segment and the (i+1)-th segment, respectively. The weighted average of the left-middle and middle-right overlapping regions is calculated using equation (2). The weight decreases linearly with the distance between the point and the center of the line segment, and effectively eliminates any jumps or discontinuities near the boundary of adjacent parts. (v) Replace the hollow areas with the weighted values ​​of the hollow areas obtained in step (iv), and repeat steps (ii) to (iv) until all hollow areas are eliminated, to obtain two terrain models extracted by row and column; (vi) The final terrain model is obtained by weighted averaging of the two terrain models.

Images