Digital elevation model hole repairing method of airborne laser radar

Abstract

The invention relates to a digital elevation model hole repairing method for an airborne laser radar, and the method comprises the following steps: S1, extracting ground point data through the airborne laser radar, and creating a TIN (Triangulated Irregular Network) mesh model; s2, setting a side length threshold value and an angle threshold value, rejecting unqualified triangular patches in the TIN mesh model by using the side length threshold value and the angle threshold value, and performing monitoring to obtain a polygonal hole; s3, filling the polygonal holes in the feature plane to obtain the feature plane with the filled holes; s4, according to the TIN irregular triangulation network mesh model, constructing an implicit curved surface which is close to the original terrain and is smoothly connected with the original terrain; and S5, fitting the implicit curved surface with the hole-filled feature plane by using a gradient descent method to complete hole repair. According to the method, the situation that the quality of a digital elevation model (DEM) is poor due to local data missing of the airborne laser radar is relieved, and smooth transition between hole repairing data and original data is achieved.

Description

Technical Field

[0001] The invention belongs to the technical field of hole repair, and in particular relates to a method for repairing holes in a digital elevation model of an airborne laser radar. Background Art

[0002] In recent years, airborne LiDAR technology, due to its unique advantage of multi-echo penetration of vegetation, can quickly, efficiently and safely obtain high-density, high-precision ground point cloud data, thereby generating a high-precision digital elevation model (DEM). This technology is widely used in topographic mapping, geological disaster monitoring and other fields, and its performance is particularly outstanding in complex terrain environments. However, due to factors such as ground object occlusion, laser scattering and terrain complexity, it is difficult to avoid the generation of data holes during LiDAR scanning. These holes will cause the local quality of the generated DEM to deteriorate, causing surface distortion and blurred local landform details in the subsequent three-dimensional model construction or topographic map drawing, which not only reduces the recognition of geological disasters, but also affects the visual effect of the three-dimensional model and the aesthetics of the map. Therefore, repairing data holes to restore the original terrain as much as possible is of great practical significance for improving the quality of DEM.

[0003] In the prior art, hole repair methods are mainly divided into two types: direct repair of scattered point clouds and data hole repair based on grid models. However, these two methods have obvious defects in practical applications. The method of directly repairing scattered point clouds usually ignores the influence of hole size on the repair results, resulting in large amount of calculation and limited effect. Although the data hole repair method based on the grid model can improve the repair efficiency to a certain extent, it is not effective when dealing with large-area data holes under complex terrain, especially in complex environments such as mountainous surfaces, and its applicability is limited. These defects seriously restrict the application of DEM in complex terrain, and it is urgent to propose a more efficient and adaptable data hole repair method. Summary of the invention

[0004] In view of the shortcomings of the prior art, the present invention provides a method for repairing holes in a digital elevation model of an airborne laser radar, which alleviates the poor quality of the digital elevation model DEM caused by the missing of local data of the airborne laser radar, realizes a smooth transition between the hole repair data and the original data, and has a good repair effect.

[0005] To achieve the above object, the technical solution of the present invention is as follows:

[0006] A method for repairing holes in a digital elevation model of an airborne laser radar comprises the following steps:

[0007] S1, using airborne laser radar to extract ground point data and create a TIN irregular triangulated mesh model;

[0008] S2, set the side length threshold and the angle threshold, use the side length threshold and the angle threshold to eliminate the unqualified triangular patches in the TIN irregular triangular mesh model, and perform monitoring to obtain polygon holes;

[0009] S3, fill the polygon holes on the feature plane to obtain a feature plane with the holes filled;

[0010] S4, according to the TIN irregular triangular mesh model, construct an implicit surface that is similar to the original terrain and smoothly connected;

[0011] S5, use the gradient descent method to fit the implicit surface with the feature plane with the holes filled to complete the hole repair.

[0012] The present invention repairs holes on the feature plane. The repair algorithm is fast and simple. At the same time, combined with the implicit surface of the radial basis function, it makes full use of the hole boundary points and their neighborhood information, and can achieve a smooth transition between the hole repair data and the original model. After repair, the terrain features are enhanced and become more obvious, and the repair effect is good.

[0013] Preferably, the step S2 includes the following steps:

[0014] S201, set the side length threshold and the angle threshold, and eliminate the unqualified triangular patches in the TIN irregular triangular mesh model;

[0015] S202, select an edge of the TIN irregular triangular mesh model as the starting calculation edge;

[0016] S203, according to the starting calculation edge, calculate the topological attributes of each edge;

[0017] S204, according to the topological attributes of each edge, judge that one side of the edge is a triangle and the other side is a polygon. If so, the edge is a hole boundary edge and enter step S205. Otherwise, the edge is not a hole boundary edge, and return to step S202;

[0018] S205: According to the hole boundary edges, perform closed connection to obtain polygon holes.

[0019] The present invention eliminates the unqualified triangular patches in the TIN irregular triangular mesh model, ensures the formation of closed polygon holes, and improves the hole repair effect.

[0020] Preferably, the side length threshold is twice the average side length of the TIN irregular triangular mesh model, the angle threshold is 15°, and the angle of the interior angle of the triangular patch cannot be less than 15°.

[0021] According to the set side length threshold and angle threshold, the present invention removes triangular patches with unsatisfactory rendering effects, ensuring the subsequent formation of closed polygon holes and improving the hole repair effect.

[0022] Preferably, step S3 includes the following steps:

[0023] S301, deleting the isolated triangular patches in the TIN irregular triangular network grid model;

[0024] S302, using the space projection formula to project the hole onto the feature plane and calculating the interior angle size of the polygon hole;

[0025] S303, starting from the smallest interior angle of the polygon hole and constructing triangular patches according to different value ranges of the interior angles;

[0026] S304, judging whether the number of remaining boundary points on the feature plane for constructing the triangular patches is equal to 3 and the side length is not greater than the side length threshold. If so, the feature plane with the hole filled is obtained; otherwise, return to step S303.

[0027] The present invention performs hole repair on the feature plane, and can obtain a more perfect repair effect.

[0028] Preferably, in step S303:

[0029] Starting from the smallest interior angle of the polygon hole, when the interior angle α satisfies α ≤ 90°, connect the adjacent points of the interior angle. If the side formed by connecting the adjacent points of the interior angle is greater than twice the average side length in the TIN irregular triangular network grid model, add a point at the midpoint of the side formed by the adjacent points, and connect the interior angle with the added point to form two triangular patches;

[0030] Starting from the smallest interior angle of the polygon hole, when the interior angle α satisfies 90° < α ≤ 120°, add a point at the same distance as the average side length in the TIN irregular triangular network grid model on the angle bisector of the interior angle, and connect the added point with the two adjacent points of the interior angle respectively to form two triangular patches;

[0031] Starting from the smallest interior angle of the polygon hole, when the interior angle α satisfies 120° < α ≤ 180°, add two points at the same distance as the average side length in the TIN irregular triangular network grid model on the trisector of the interior angle, and connect the adjacent points of the interior angle and the added two points in sequence to form three triangular patches.

[0032] The present invention adopts different triangular patch filling methods according to different interior angles, which is convenient for quickly constructing triangular patches and restoring the original shape of the model as much as possible.

[0033] Preferably, in step S4, according to the TIN irregular triangular network grid model, the grid vertices 3 to 5 grids extending outward from the hole are fitted by using radial basis functions to obtain an implicit surface that is similar to the original terrain and smoothly connected.

[0034] The present invention obtains a patched grid that is similar to the original terrain and smoothly connected, which is convenient for subsequent operations.

[0035] Preferably, the expression of the implicit surface is as follows:

[0036]

[0037] P(r) = p0 + p1x + p2y + p3z

[0038] where r is the grid vertex passed through when establishing the implicit surface, n is the number of grid vertices in space, c j is the jth sampling point required for establishing the implicit surface, w j is the jth weight corresponding to each sampling point, is the radial basis function, F(r) is the implicit surface, P(r) is the polynomial, p0, p1, p2, and p3 are all polynomial coefficients, and x, y, and z are all coordinate axes.

[0039] The present invention uses radial basis functions to construct an implicit surface that is similar to the original terrain and smoothly connected, enabling it to perfectly fit the feature plane that has been patched.

[0040] Preferably, step S5 includes the following steps:

[0041] S501: Using the gradient descent method, the iteration points in the feature plane with the hole filled are iteratively approximated to the implicit surface along the surface gradient direction;

[0042] S502: Determine whether the iteration is completed. If the iteration is completed, the fitting of the implicit surface and the feature plane with the hole filled ends, and the hole repair is completed. Otherwise, return to S501.

[0043] The present invention obtains a patched TIN irregular triangular network grid model, which is convenient for subsequent analysis and research of the terrain.

[0044] Preferably, the expression of the iterative approximation is as follows:

[0045]

[0046] where C1 is the newly calculated iteration point, C0 is the iteration starting point, γ is the learning efficiency, is the gradient of the surface, is the representation symbol of the partial derivative, is the partial derivative of the function F with respect to x, is the partial derivative of function F with respect to y, is the partial derivative of function F with respect to z, where x, y, and z are all coordinate axes.

[0047] The present invention uses iterative approximation to adjust the vertices of the filled triangular patches to the implicit surface, ensuring smooth connection between the feature plane with the hole filled and the implicit surface, and improving the hole repair effect.

[0048] Preferably, the expression for the completion of the iteration is as follows:

[0049] ||C1 - C0|| ≤ ε

[0050] where ε is a preset limiting error, || || is the norm operator, C1 is the newly calculated iteration point, and C0 is the iteration starting point.

[0051] The present invention uses the limiting error to limit the iterative approximation, ensuring the fitting completion between the feature plane with the hole filled and the implicit surface, and improving the hole repair effect.

[0052] Compared with the prior art, the advantages of the present invention are as follows:

[0053] The present invention provides a method for repairing holes in the digital elevation model of an airborne lidar, which has the advantages of fast and simple repair algorithms, can significantly reduce the computational complexity, improve the processing efficiency, and is applicable to the repair of large-scale point cloud data. By combining the implicit surface technology of radial basis functions, it fully utilizes the hole boundary points and their neighborhood information to ensure a smooth transition between the repaired data and the original model, making the repaired terrain model seamlessly connect with the original terrain, enhancing the detailed features of the terrain, including surface boundary features, deformation features, and local landform details, and improving the clarity and integrity of the digital elevation model; the present invention has strong adaptability and can achieve excellent results in hole repair in complex terrain environments and different areas, especially effectively solving the deficiencies of the prior art in mountainous areas and large holes. By adopting a step-by-step optimized repair strategy, constructing triangular patches through the feature plane and combining the gradient descent method, the vertices of the filled triangular patches are adjusted to the implicit surface to achieve high-precision repair, significantly improving the quality and three-dimensional visual effect of the terrain model; the present invention has important application value in the fields of topographic surveying and mapping, geological disaster monitoring, three-dimensional modeling, and topographic map drawing, can enhance the recognition of geological disasters and the mapping quality, and provides an excellent solution for the application of the digital elevation model in complex terrains. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 is the flowchart of the method for repairing holes in the digital elevation model of the airborne lidar of the present invention;

[0055] Figure 2It is a diagram of a TIN irregular triangular network grid model of the present invention;

[0056] Figure 3 It is a schematic diagram of the newly added triangular patches in the structure of the present invention;

[0057] Figure 4 It is a schematic diagram of adjusting the triangular grid vertices to the implicit surface of the present invention;

[0058] Figure 5 It is a diagram of the original point filtering processing result and hole repair result of the present invention;

[0059] Figure 6 It is a repair example of the present invention;

[0060] Figure 7 It is the point cloud holes and test point arrangement in the experimental area of the present invention;

[0061] Figure 8 It is the repair accuracy analysis area with different terrain complexities of the present invention;

[0062] Figure 9 It is the mountain shadow diagram of the comparison of the repair of various hole areas with different terrain complexities of the present invention. Detailed implementation manners

[0063] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art belong to the present invention.

[0064] In addition, the described features, structures or characteristics can be combined in any suitable manner in one or more embodiments. In the following description, many specific details are provided to give a full understanding of the embodiments of the present application. However, those skilled in the art will realize that the technical solutions of the present application can be practiced without one or more of the specific details, or other methods, components, devices, steps, etc. can be adopted. In other cases, well-known methods, devices, implementations or operations are not shown or described in detail to avoid obscuring various aspects of the present application.

[0065] As Figure 1 shown, this embodiment discloses a method for repairing holes in a digital elevation model of an airborne lidar, including the following steps:

[0066] S1, extracting ground point data by using an airborne lidar and creating a TIN irregular triangular network grid model;

[0067] S2. Set the side length threshold and the angle threshold, use the side length threshold and the angle threshold to eliminate the unqualified triangular patches in the TIN irregular triangular mesh model, and conduct monitoring to obtain polygon holes;

[0068] S3. Fill the polygon holes on the feature plane to obtain a feature plane with the holes filled;

[0069] S4. According to the TIN irregular triangular mesh model, construct an implicit surface that is similar to the original terrain and is smoothly connected;

[0070] S5. Use the gradient descent method to fit the implicit surface with the feature plane with the holes filled to complete the hole repair.

[0071] As Figure 2 shown, it is a TIN irregular triangular mesh model, and the area with darker color is the hole that needs to be repaired. As Figure 2 shown, there is a lack of point cloud data in the hole area, the triangular network in the area is sparse, and features such as longer triangle side lengths or too small angles appear. Such triangular patches will cause problems such as streaking and fuzzing, and misrepresentation of terrain and landform information. They are unqualified triangular patches. Therefore, certain side length and angle thresholds are set to eliminate the unqualified triangular patches at the hole positions. The side length threshold is determined by twice the average side length of the triangular mesh model. The specific steps are S2.

[0072] The specific steps of S2 are as follows:

[0073] S201: Set the side length threshold and the angle threshold, and eliminate the unqualified triangular patches in the TIN irregular triangular mesh model;

[0074] S202: Select an edge of the TIN irregular triangular mesh model as the starting calculation edge;

[0075] S203: Calculate the topological attributes of each edge according to the starting calculation edge;

[0076] S204: According to the topological attributes of each edge, judge that one side of the edge is a triangle and the other side is a polygon. If so, the edge is a hole boundary edge and go to step S205. Otherwise, the edge is not a hole boundary edge, and return to step S202;

[0077] S205: Make a closed connection according to the hole boundary edge to obtain a polygon hole.

[0078] Among them, the side length threshold is twice the average side length of the TIN irregular triangular network grid model, the angle threshold is 15°, and the interior angle of the triangular patch cannot be less than 15°. The unqualified triangular patches in the TIN irregular triangular network grid model are removed according to the threshold to ensure the formation of a closed polygon hole subsequently and improve the hole repair effect.

[0079] As Figure 3 shown, it is a schematic diagram for constructing new triangular patches. Figure 3 (a) is the repair method when the interior angle α satisfies α ≤ 90°. Figure 3 (b) is the repair method when the interior angle α satisfies 90° < α ≤ 120°. Figure 3 (c) is the repair method when the interior angle α satisfies 120° < α ≤ 180°. Let the interior angle of the polygon hole be α, and the two adjacent sides of the interior angle be P i-1 P i and P i P i+1 . Analyze that there are two types of polygon holes on the hole feature plane: convex polygons and concave polygons. The interior angles α of convex polygons all satisfy 0° < α ≤ 180°, and the smallest interior angle α of concave polygons must satisfy 0° < α ≤ 180°. Moreover, with the addition of new triangular patches, the interior angles greater than 180° will also be gradually divided and finally satisfy 0° < α ≤ 180°. Therefore, the smallest angle α in the triangle interior angles must be less than 180°.

[0080] Furthermore, step S3 includes the following steps:

[0081] S301, delete the isolated triangular patches in the TIN irregular triangular network grid model;

[0082] S302, use the space projection formula to project the hole onto the feature plane and calculate the size of the interior angle of the polygon hole;

[0083] S303, take the smallest interior angle of the polygon hole as the starting interior angle and construct triangular patches according to different value ranges of the interior angles;

[0084] S304, determine whether the number of remaining boundary points on the feature plane where the triangular patches are constructed is equal to 3 and the side length is not greater than the side length threshold. If so, obtain the feature plane with the hole filled. Otherwise, return to step S303.

[0085] In this embodiment, the specific method of step S303 is:

[0086] Taking the interior angle of the smallest polygon hole as the starting interior angle, when the interior angle α satisfies α ≤ 90°, connect the adjacent points of the interior angle. If the side formed by connecting the adjacent points of the interior angle is greater than twice the average side length in the TIN irregular triangular network grid model, then add a point at the midpoint of the side formed by the adjacent points, connect the interior angle with the added point, and form two triangular patches;

[0087] Taking the interior angle of the smallest polygon hole as the starting interior angle, when the interior angle α satisfies 90° < α ≤ 120°, then add a point at a position on the angle bisector of the interior angle where the distance is the same as the average side length in the TIN irregular triangular network grid model, connect the added point with the two adjacent points of the interior angle respectively, and form two triangular patches;

[0088] Taking the interior angle of the smallest polygon hole as the starting interior angle, when the interior angle α satisfies 120° < α ≤ 180°, then add two points respectively at positions on the trisector of the interior angle where the distance is the same as the average side length in the TIN irregular triangular network grid model, and connect the adjacent points of the interior angle and the added two points in sequence to form three triangular patches.

[0089] Furthermore, in step S4, according to the TIN irregular triangular network grid model, use the radial basis function to fit the grid vertices that extend 3 to 5 grids outward from the hole, and obtain an implicit surface that is similar to the original terrain and smoothly connected.

[0090] Among them, the expression of the implicit surface is as follows:

[0091]

[0092] P(r) = p0 + p1x + p2y + p3z

[0093] Among them, r is the grid vertex passed through when establishing the implicit surface, n is the number of grid vertices in space, c j is the jth sampling point required for establishing the implicit surface, w j is the jth weight value corresponding to each sampling point, is the radial basis function, F(r) is the implicit surface, P(r) is the polynomial, p0, p1, p2, and p3 are all polynomial coefficients, and x, y, and z are all coordinate axes.

[0094] In order to obtain a patched grid that is similar to the original terrain and smoothly connected, use the grid vertices that extend 3 to 5 rings outward from the hole boundary to fit the implicit surface of the radial basis function, and use the gradient descent method to gradually approximate all the newly added triangular patch vertices along the surface gradient direction to the implicit surface until the newly generated triangular patch vertices are adjusted to the fitted hole surface.

[0095] The said step S5 includes the following steps:

[0096] S501: Using the gradient descent method, iteratively approximate the iterative points in the feature plane with holes filled along the surface gradient direction towards the implicit surface;

[0097] S502: Determine whether the iteration is completed. If the iteration is completed, the fitting of the implicit surface and the feature plane with holes filled ends, and the hole repair is completed. Otherwise, return to S501.

[0098] As Figure 4 shown, it is a schematic diagram of the adjustment of triangular mesh vertices towards the implicit surface, that is, the process of fitting the feature plane with holes filled and the implicit surface using the gradient descent method.

[0099] Among them, the expression for iterative approximation is as follows:

[0100]

[0101] Among them, C1 is the newly calculated iterative point, C0 is the iterative starting point, γ is the learning efficiency, is the gradient of the surface, is the representation symbol of the partial derivative, is the partial derivative of the function F with respect to x, is the partial derivative of the function F with respect to y, is the partial derivative of the function F with respect to z, and x, y, and z are all coordinate axes.

[0102] Furthermore, the expression for the completion of the iteration is as follows:

[0103] ||C1 - C0|| ≤ ε

[0104] Among them, ε is the preset limit error, || || is the norm operator, C1 is the newly calculated iterative point, and C0 is the iterative starting point.

[0105] Intercept a point cloud data at Songlinping, Luding County. The original data format is las format, and a total of 6,810,492 point clouds are read. The point cloud density is 37.72 pts / m 2 . As Figure 5 shown, it is the result diagram of the original point filtering process and the hole repair result. The original points are as Figure 5 (a) shown. The main vegetation features are some shrubs and tall trees. Due to vegetation occlusion or complex terrain, after filtering the original points, 775,689 ground points are extracted, as Figure 5 (b) shown. There are some holes in the ground points. Using the method of the present invention, the holes are repaired. The effect before hole repair is as Figure 5 (c) shown, and the effect after hole repair is as Figure 5 (d) shown.

[0106] As Figure 6As shown, it is a patching example. To more intuitively see the patching effect, the data before and after patching are used to generate DEM (Digital Elevation Model) with a unified resolution, and hillshade visualization is performed. The vegetation coverage in this area is dense as Figure 6 (a). The vegetation coverage in this area is dense, and there are holes in the ground point cloud, resulting in streaks in the generated hillshade map and missing local geomorphic features, as Figure 6 (b) shows; after patching the point cloud holes, a hillshade map as Figure 6 (c) is obtained. Compared with the unpatched one, it is more beautiful and complete, and the landslide details shown are richer; Figure 6 (e) and Figure 6 (g) are both hillshade maps of the holes, Figure 6 (f) and Figure 6 (h) are both hillshade maps of the repaired area. By comparison, it can be seen that after patching, the secondary sliding of the landslide is restored, the gully on the right side of the landslide is obvious after repair, and the boundary features of the landslide are enhanced. After patching using the method of this article, the landslide contour is shown more clearly, and the deformation characteristics of the landslide form are more obvious.

[0107] In order to detect the patching accuracy of the method of the present invention under different terrain conditions, according to the terrain complexity index of different terrains, slope, terrain roughness, and terrain undulation are selected as the constituent factors of the terrain complexity index. The larger the terrain complexity index value, the more complex the terrain. After determining the constituent factors of the terrain complexity index, the principal component analysis method is used to allocate the weights of each factor. After obtaining the weights, the terrain complexity index is obtained by weighted overlay using the raster calculator in the GIS software.

[0108] As Figure 7 shown, it is the layout of the point cloud holes and test points in the experimental area; an area is selected for the experiment. The terrain complexity index range of this area is 0 - 0.667, and the average is 0.402, as Figure 7 (a) shows; in this area, point cloud holes with three different terrains, namely flat area, slope area, and ridge area, are constructed, and the hole area is 2000 m 2 ; first, the data before and after patching in the hole area are extracted. The DEM elevation constructed from the original ground point cloud data is used as the true value, and then the DEM elevation generated from the patched triangular mesh model is used as the predicted value. The constructed DEM resolution is 0.2 m. Then, the known value and the predicted value are compared and analyzed, and a detection point is arranged every 5 m on the DEM products in the three experimental areas; as Figure 6 (b) shows, it is the DEM product of the ridge area, as Figure 6 (c) shows, it is the DEM product of the slope area, as Figure 6 (d) shows, it is the DEM product of the flat area.

[0109] Compare the difference between the elevation value measured at the detection point on the original DEM and the elevation value of the DEM after patching one by one. To more intuitively reflect the accuracy of point cloud hole patching, the root mean square error of elevation is selected as the evaluation index. The smaller this value is, the smaller the error. The calculation formula of the root mean square error of elevation is as follows:

[0110]

[0111] where RMSE is the root mean square error of elevation, and Z k is the true elevation value of the detection point k; z k is the measured elevation value of the detection point k; n is the number of detection points, and k is the kth detection point.

[0112] The test results in different terrain experimental areas are as follows. In the flat area, the maximum elevation difference is 0.5 m, and the average elevation difference is 0.193 m; in the slope area, the maximum elevation difference is 0.95 m, and the average elevation difference is 0.378 m; in the ridge area, the maximum elevation difference is 1.98 m, and the average difference is 0.547 m. Through the statistical analysis of the elevation differences of the detection points, in the flat area, the terrain complexity is 0.06, the terrain change is small, the elevation differences of 100% of the test points are less than 0.5 m, the RMSE is 0.23, the terrain after patching is similar to the original, and the patching accuracy is high; in the slope area, 64.6% of the test points have elevation differences in the range of 0 - 0.5 m, and 35.4% are in the range of 0.5 - 1 m. Compared with the flat area, the terrain complexity increases by 0.207, and the patching accuracy decreases by 50%; in the ridge experimental area, 49.2% of the test points have elevation differences in the range of 0 - 0.5 m, 26.1% are in the range of 0.5 - 1 m, and 24.7% of the test points have differences greater than 1 m. As the terrain complexity increases, the RMSE value becomes larger, and the root mean square error of elevation reaches 1.21, which cannot meet the accuracy requirements of 1:2000 DEM for mountain production. To sum up, the hole patching accuracy under different terrains decreases with the increase of terrain complexity. At the same time, all three experimental areas show that the closer the hole patching accuracy is to the hole boundary, the better the patching effect and the smaller the DEM elevation difference.

[0113] To monitor the influence of different micro - landforms on the patching accuracy, as Figure 8 shown, for the analysis area of patching accuracy with different terrain complexities, based on the terrain complexity index, four gradient holes of 200 m 2 , 500 m 2 , 1000 m 2 and 2000 m 2 are respectively constructed at the rear scarp, middle mound, and boundary gully of the landslide, focusing on analyzing the influence of different terrain complexities on the patching accuracy of point cloud holes.

[0114] As Figure 9As shown in the figure, it is a mountain shadow map for comparing the repair of various hole areas with different terrain complexities. From the comparison before and after repairing different hole areas, it can be seen that as the hole area increases, the repair effect gradually decreases. When the hole area is 200m 2 , the mountain shadows after repairing at three different landforms have a small difference from the original shadow map in terms of image display and are almost the same as the original landform. When the hole area reaches 1000m 2 , the repair marks of the holes gradually appear. The shadow map at the hole area is finer than the original map, and the difference from the surrounding area of the hole gradually increases. However, the mountain shadow with the hole has obvious streaks, covering up the landform information. After repairing the hole, the original landform can be restored to a certain extent. When the hole area reaches 2000m 2 , the repair effect significantly decreases, and the repaired mountain shadow map is quite different from the original map. Although the hole can be repaired, the repair deviation of the hole is large, and the repair effect is not ideal.

[0115] In summary, the terrain complexity and the size of the enlarged area are the key factors affecting the terrain repair accuracy. Using the method of the present invention to repair the terrain, for the terrain in flat areas and slope areas, the terrain and landform can be restored to a certain extent, the quality of the DEM digital elevation model can be improved, and it has a high repair accuracy. For the mountainous terrain, when the hole area ≤ 200m 2 , it meets the accuracy requirements of the 1:500 DEM digital elevation model. When the hole area ≤ 500m 2 , it meets the accuracy requirements of the 1:1000 DEM digital elevation model. When the hole area ≤ 1000m 2 , it meets the accuracy requirements of the 1:2000 DEM digital elevation model.

[0116] In summary, the present invention discloses a method for repairing holes in a digital elevation model of an airborne lidar, which includes extracting ground point data using the airborne lidar and creating a TIN irregular triangular network grid model; removing unqualified triangular patches by setting side length thresholds and angle thresholds, and identifying and closing to form polygon holes; filling the hole area on the feature plane, constructing an implicit surface and achieving smooth connection with the original terrain; finally completing hole repair through the gradient descent method. This method adopts the implicit surface technology combined with radial basis functions and an efficient repair algorithm, making full use of the hole boundary points and their neighborhood information to ensure a smooth transition between the repaired data and the original model. The repaired terrain model is clearer and more complete, capable of enhancing terrain boundary features and geomorphic detail features. The present invention has significant advantages such as fast, efficient, and highly adaptable repair algorithms, and can achieve good results in repairing complex terrains and large-area holes. Especially in complex environments such as mountainous areas, it solves the technical problems that cannot be effectively handled by existing technologies. Through a high-precision repair strategy, it significantly improves the quality and accuracy of the digital elevation model, while enhancing the recognition of geological disasters and the visual effect of the 3D model, providing an efficient solution for the application of the digital elevation model in terrain mapping, geological disaster monitoring, 3D modeling, topographic map drawing and other fields. The present invention not only overcomes the limitations of existing technologies in dealing with holes in complex terrains, but also makes important contributions to improving the accuracy and application scope of the digital elevation model, has wide industrial applicability and practical value, and is of great significance for promoting the technological development of the digital elevation model in the field of terrain monitoring and modeling.

[0117] The above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit them; under the idea of the present invention, the technical features in the above embodiments or different embodiments can also be combined, and the steps can be implemented in any order, and there are many other variations in different aspects of the present invention as described above. For the sake of brevity, they are not provided in detail; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements on some of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for repairing holes in a digital elevation model of an airborne laser radar, characterized in that: The following steps are involved: S1, using airborne laser radar to extract ground point data and create a TIN irregular triangulated mesh model; S2, setting the edge length threshold and the angle threshold, using the edge length threshold and the angle threshold to remove unqualified triangular facets in the TIN irregular triangulated mesh model, and monitoring them to obtain polygonal holes; S3, filling the polygonal holes on the feature plane to obtain a feature plane with completed hole filling; S4, based on the TIN irregular triangulated network mesh model, construct an implicit surface that is close to the original terrain and smoothly connected; S5, using the gradient descent method to fit the implicit surface with the feature plane on which the hole filling has been completed, to complete the hole repair.

2. The method for repairing holes in a digital elevation model of an airborne laser radar according to claim 1, characterized in that: The step S2 comprises the following steps: S201, setting edge length threshold and angle threshold, and removing unqualified triangular facets in the TIN irregular triangulated mesh model; S202, selecting an edge of the TIN irregular triangulated network model as a starting edge for calculation; S203, calculating the topological attribute of each edge according to the starting edge calculation; S204, judging whether one side of the edge is a triangle and the other side is a polygon according to the topological attributes of each edge, if so, the edge is a hole boundary edge, and proceeding to step S205, otherwise, the edge is not a hole boundary edge, and returning to step S202; S205: Perform a closed connection according to the boundary edges of the hole to obtain a polygonal hole.

3. The method for repairing holes in a digital elevation model of an airborne laser radar according to claim 2, characterized in that: The edge length threshold is twice the average edge length of the TIN irregular triangulated mesh model, the angle threshold is 15°, and the internal angle of the triangular facet cannot be less than 15°.

4. The method for repairing holes in a digital elevation model of an airborne laser radar according to claim 3, characterized in that: The step S3 comprises the following steps: S301, deleting isolated triangular facets in the TIN irregular triangulated mesh model; S302, projecting the hole onto the feature plane using a spatial projection formula, and calculating the size of the inner angle of the polygonal hole; S303, taking the minimum polygon hole inner angle as the starting inner angle, and constructing a triangular face according to the value range of different inner angles; S304, determining whether the number of remaining boundary points on the feature plane for constructing the triangular face is equal to 3 and the side length is not greater than the side length threshold, if so, obtaining the feature plane for which the hole filling has been completed, otherwise, returning to step S303.

5. The method for repairing holes in a digital elevation model of an airborne laser radar according to claim 4, characterized in that: In step S303: Taking the minimum polygon hole inner angle as the starting inner angle, when the inner angle α satisfies α≤90°, connecting the adjacent points of the inner angle, if the edge formed by connecting the adjacent points of the inner angle is greater than twice the average edge length in the TIN irregular triangulated mesh model, then adding a point at the midpoint of the edge formed by the adjacent points, connecting the inner angle and the added point, to form two triangular facets; Taking the minimum polygon hole inner angle as the starting inner angle, when the inner angle α satisfies 90°<α≤120°, a point is added on the angle bisector of the inner angle at a distance equal to the average side length in the TIN irregular triangulated mesh model, and the added point is connected to two adjacent points of the inner angle to form two triangular facets; Taking the inner angle of the minimum polygonal hole as the starting inner angle, when the inner angle α satisfies 120°<α≤180°, two points are added on the trisection line of the inner angle at a distance equal to the average side length in the TIN irregular triangulated mesh model, and the adjacent points of the inner angle and the two added points are connected in sequence to form three triangular facets.

6. The method for repairing holes in a digital elevation model of an airborne laser radar according to claim 5, characterized in that: In step S4, according to the TIN irregular triangulated network mesh model, radial basis functions are used to fit the mesh vertices extending 3 to 5 grids outward from the hole to obtain an implicit surface that is close to the original terrain and smoothly connected.

7. The method for repairing holes in a digital elevation model of an airborne laser radar according to claim 6, characterized in that: The expression of the implicit surface is as follows: P(r)=p0+p1x+y2y+y3z Among them, r is the mesh vertex passed when establishing the implicit surface, n is the number of mesh vertices in space, and c is the number of mesh vertices in space. j To establish the jth sampling point required for the implicit surface, w j is the jth weight corresponding to each sampling point, is the radial basis function, F(r) is the implicit surface, P(r) is a polynomial, p0, p1, p2 and p3 are all polynomial coefficients, and x, y, z are all coordinate axes.

8. The method for repairing holes in a digital elevation model of an airborne laser radar according to claim 7, characterized in that: The step S5 comprises the following steps: S501: using the gradient descent method, iteratively approximate the iterative points in the feature plane where the hole filling has been completed to the implicit surface along the surface gradient direction; S502: Determine whether the iteration is completed. If the iteration is completed, the implicit surface is fitted with the feature plane with the hole filling completed, and the hole repair is completed. Otherwise, return to S501.

9. The method for repairing holes in a digital elevation model of an airborne laser radar according to claim 8, characterized in that: The expression of the iterative approximation is as follows: Among them, C1 is the calculated new iteration point, C0 is the iteration starting point, γ is the learning efficiency, is the gradient of the surface, is the symbol for partial derivatives, is the partial derivative of function F with respect to x, is the partial derivative of function F with respect to y, is the partial derivative of the function F with respect to z, where x, y, and z are coordinate axes.

10. The method for repairing holes in a digital elevation model of an airborne laser radar according to claim 9, characterized in that: The expression of the iteration is as follows: ||C1-C0||≤ε Among them, ε is the preset limit error, || || is the norm operator, C1 is the calculated new iteration point, and C0 is the iteration starting point.

Images