A method and system for reconstructing a vegetation-shielded slope DEM by fusing multi-source remote sensing data

Abstract

The application discloses a kind of fusion multi-source remote sensing data's vegetation sheltering side slope DEM reconstruction method and system, belong to remote sensing surveying and mapping and digital terrain modeling technical field.Getting unmanned aerial vehicle laser radar point cloud, orthographic image and existing digital elevation model, after data pre-processing and initial ground filtering, to initial non-ground point is carried out multi-scale terrain continuity analysis, comprehensive point to local reference surface distance, normal angle, curvature and ground connection relationship, it is judged as ground point to form revised ground point set and extract terrain structure line in the non-ground point that meets continuity condition under multiple scales;According to the direction distribution of missing area peripheral ground point and the consistency with existing DEM terrain trend, it can be determined that the region can be reconstructed;Terrain sub-region is divided using terrain structure line, and the objective function including boundary elevation constraint, DEM trend constraint and local smoothing constraint is established in each sub-region to reconstruct ground surface;The revised ground point is fused with reconstructed point cloud to generate high-precision digital elevation model.

Description

Technical Field

[0001] This invention relates to the fields of remote sensing mapping, 3D point cloud information processing, and digital terrain modeling. Specifically, it relates to a method and system for generating a vegetation-covered slope DEM by using UAV lidar point clouds as the main observation data, fusing orthophotos (DOM) with existing digital elevation models (DEMs), and through candidate ground structure point back-identification, reconstructable area determination, and structurally constrained surface reconstruction. It is applicable to digital terrain modeling and engineering surveying applications in complex terrain areas with vegetation cover, such as road cut slopes, high fill slopes, and natural slopes. Background Technology

[0002] Mountainous highways are characterized by numerous road cut slopes, natural slopes, and engineering protective structures. Due to factors such as terrain cutting, slope gradient, vegetation growth, and scanning perspective, traditional manual surveying struggles to continuously acquire the spatial morphology of slope surfaces. UAV-based lidar can rapidly acquire three-dimensional point clouds of slopes and has some penetration capability into vegetation canopy; orthophotos can provide information on surface cover such as vegetation, bare rock, road platforms, and artificial structures; and existing digital elevation models (DEMs) can provide regional topographic trends. Therefore, generating digital elevation models of slopes using multiple types of observation data has become an important technical approach in highway slope surveying and monitoring.

[0003] Existing methods for generating digital elevation models (DEMs) for slopes typically employ the following technical approach: First, methods such as CSF (Content Simulation Filtering), PMF (Progressive Morphological Filtering), or PTD (Progressive Density Irregular Triangular Network Filtering) are used to separate ground points from non-ground points in the original point cloud. Then, methods such as moving least squares, Kriging interpolation, spline interpolation, radial basis functions, or irregular triangular networks are used to perform surface fitting or cavity repair in sparse and missing areas of the point cloud. Existing methods often uniformly treat areas that meet the point density or cavity scale conditions as the interpolation region, without distinguishing between the reliability of surface observations and the reconstructability of missing areas. Furthermore, some methods introduce local terrain feature preservation constraints or existing DEM trend constraints during the interpolation process to improve the modeling effect for steep slope areas.

[0004] However, the existing technologies still have the following shortcomings: First, existing point cloud filtering methods mainly classify points based on geometric features such as elevation difference, local slope, or distance from the point to the fitted surface. In steep slope areas, these methods can easily misclassify real surface structures such as bare rock surfaces, steep slopes, ridges, and slope breaks as non-ground points. Furthermore, when using continuous surface interpolation, smoothing may be performed across engineering structure boundaries such as graded platforms, drainage ditches, and retaining walls, weakening key natural terrain and artificial engineering structures. Second, sparse point clouds or missing surfaces can be caused not only by vegetation shading but also by various other reasons such as excessive scanning incident angle, insufficient coverage at the flight path edge, shading by artificial structures, or abrupt changes in the surface itself. The reliability and reconstructability of surrounding observation data vary depending on the cause of the missing area. Existing methods typically treat areas that meet the point density or cavity scale conditions uniformly as interpolation areas, without fully considering the directional distribution of reliable surrounding ground points and their consistency with the existing DEM topographic trend. This results in some insufficiently constrained areas being interpolated and filled. Third, while orthophotos and existing DEMs can provide auxiliary information, their reliability varies across different locations: orthophotos primarily reflect land cover types and cannot directly provide the true ground elevation under vegetation cover; existing DEMs may have smoothing or timeliness errors in areas such as steep slopes, valleys, and ridges. How to distinguish between point cloud missing areas of different origins and determine their reconstructability, reducing misclassification of true terrain structure points and terrain distortion caused by unfounded interpolation, remains a pressing technical problem to be solved in the generation of DEMs for vegetation-covered slopes. Summary of the Invention

[0005] The purpose of this invention is to provide a method and system for reconstructing vegetation-shaded slopes by fusing multi-source remote sensing data, addressing the aforementioned problems. This method uses UAV laser point clouds as the main observation data, orthophotos and existing DEMs as auxiliary data, and can select ground laser scanning point clouds or measured control points for local supplementation, registration or accuracy verification.

[0006] The technical solution of the present invention is as follows: A method for reconstructing a vegetation-shaded slope DEM by fusing multi-source remote sensing data includes the following steps: Acquire multi-source remote sensing data, including UAV lidar point clouds, orthophotos, and existing digital elevation models; By unifying all data sources to the same geographic coordinate system, initial ground filtering is performed on the UAV point cloud to obtain the initial ground point set. and the initial non-ground point set ; A multi-scale three-dimensional neighborhood is established based on the median nearest neighbor distance of the initial ground point set. At each scale, a local terrain reference surface is established using the initial ground points within the neighborhood. Multi-scale terrain continuity analysis is performed on non-ground points. Non-ground points that meet the continuity condition at at least two scales are reclassified as ground points and incorporated into the initial ground point set to form a modified ground point set. And extract terrain structure lines based on the corrected ground point set. ; The corrected ground point set is divided into a regular grid, and candidate missing regions are extracted by combining orthophoto vegetation cover information. Reconstructable regions are determined based on the directional distribution of ground points around the missing regions and their consistency with the existing DEM topographic trend. ; Utilizing terrain structure lines The reconstructable area is divided into several terrain sub-regions. An objective function is established within each sub-region, incorporating boundary elevation constraints, existing DEM trend constraints, and local smoothing constraints. Solving this function yields the reconstructed surface point cloud. ; Correct the ground point set With reconstructing surface point clouds Integration, based on terrain structure lines As a forced edge, an irregular triangular network is constructed to generate a digital elevation model, and areas that do not meet the reconstruction conditions are retained as null values.

[0007] This step establishes the overall technical framework. By performing multi-scale terrain continuity analysis on the non-ground points after initial filtering and re-identifying the real ground points, the shortcomings of initial filtering in misidentifying real terrain structures such as bare rock surfaces, slope ridges, and slope toes in steep slope areas are compensated. By determining reconstructable areas based on directional distribution and DEM trend consistency, the unfounded terrain distortion caused by uniform interpolation of all sparse areas in traditional methods is avoided. By dividing terrain sub-regions by structural lines and blocking smoothing across structural lines, the reconstructed surface is prevented from crossing structural boundaries such as platforms, drainage ditches, and retaining walls, significantly improving the overall quality and terrain fidelity of the DEM for vegetation-covered slopes.

[0008] Furthermore, the multiple progressively scaled three-dimensional neighborhoods are specifically as follows: Calculate the initial ground point set Median of nearest neighbor distances for each point ,by As the distance between feature points, establish a radius of... , and The three three-dimensional spherical neighborhoods, where the scaling coefficients satisfy ; The local terrain reference surface is established as follows: when the number of initial ground points in the neighborhood meets the threshold and the orientation distribution is uniform, a robust weighted least squares or RANSAC plane fitting method is used to establish the local terrain reference surface. And calculate the fitting residual of the local terrain reference surface. ;in The unit normal vector of the reference surface. It is a plane constant. The coordinate vector of the spatial point is used. When the fitting residual is less than the threshold, it is considered a single continuous slope. When the fitting residual is large and the normal of the points in the neighborhood shows a significant bimodal distribution, two local reference surfaces are fitted respectively for the judgment of the slope break structure.

[0009] By establishing three progressively larger three-dimensional neighborhoods based on the median nearest neighbor distance of the initial ground point, a multi-scale comprehensive assessment of small-scale local details, medium-scale slope extension, and large-scale overall trends is achieved. Compared to single-scale fixed-scale analysis, this method can take into account terrain features of different scales, avoiding potential omissions or misjudgments that may occur at a single scale. Furthermore, by using fitted residuals to determine single continuous slopes or slope break structures, subsequent structural line extraction has a clear geological basis, improving the reliability of ground point re-identification and the accuracy of structural line positioning.

[0010] Furthermore, the multi-scale terrain continuity analysis includes: Calculate the non-ground points to be judged at each effective scale Distance to local topographic reference surface Local normal angle and local curvature Specifically: non-ground points Distance to local topographic reference surface : , If there are two local reference surfaces in the neighborhood, then calculate the non-ground points separately. The distances to the two reference planes are taken, and their minimum values ​​are calculated: , in, Let be the unit normal vector of the first local reference surface. Let be the plane constant of the first local reference surface. The unit normal vector of the second local reference surface. The plane constant of the second local reference surface; At non-ground points Construct the covariance matrix in the neighborhood And calculate non-ground points Local normal vector Local normal vector unit normal vector to the reference plane The spatial angle is the local normal angle. The local normal angle is expressed as: ; Let non-ground points The eigenvalues ​​of the neighborhood covariance matrix are The local curvature is then expressed as: ; Adaptive tolerance threshold for distance ,in and For adjustment coefficients, The median of the nearest neighbor distances of the initial set of ground points. The fitting residuals for the local terrain reference surface; When a non-ground point Simultaneously satisfying at least two effective scales The included normal angle is no greater than a set threshold, the curvature is no greater than a set threshold, and it can be connected to the initial ground point set through a continuous point set. At that time, the point will be reclassified as a ground point.

[0011] By comprehensively evaluating four dimensions—distance from the point to the reference surface, normal angle, curvature, and ground connectivity—a multi-angle, all-round evaluation system for non-ground points is formed. The adaptive distance threshold is dynamically adjusted based on local point spacing and fitting residuals, ensuring stable judgment capabilities across regions with varying point cloud density and surface roughness. These methods ensure that only points with continuous relationships to adjacent terrain in spatial location, surface orientation, local morphology, and topological connectivity are reclassified, effectively excluding canopy points and isolated noise points, thus improving the accuracy of ground point reclassification.

[0012] Furthermore, the extraction of terrain structure lines The specific method is as follows: Based on the modified ground point set In this study, candidate structural points are identified through abrupt changes in normal direction and local maxima of curvature. Abrupt changes in normal direction are determined by the rate of change of the angle between the normal vectors of neighboring points, while local maxima of curvature are determined by the local curvature of each grid node. Non-maximum suppression extraction is performed; connected component analysis is conducted on candidate structural points to remove discrete line segments with lengths less than a set threshold, resulting in ridge lines, toe lines, upper and lower edges of steep slopes, and platform edge lines.

[0013] By identifying abrupt changes in normal and local maxima of curvature in a modified ground point set, candidate points for terrain structures were identified. Connectivity analysis and length thresholding were then used to remove discrete noise segments, enabling accurate extraction of ridge lines, toe lines, upper and lower edges of steep slopes, and platform edges. These structural lines serve as the boundary criteria for subsequent terrain sub-region division, directly determining the effectiveness of cross-structure smoothing and providing a data foundation for preserving the morphology of key terrain and engineering structures.

[0014] Furthermore, the utilization of terrain structure lines The specific method for dividing the reconstructable region into several terrain sub-regions is as follows: Terrain structure lines Projected onto the reconstructable area The regular grid uses four-neighbor or eight-neighbor domains to establish node adjacency relationships. When the line connecting two adjacent nodes coincides with a terrain structure line... When two regions intersect, the adjacency relationship is deleted, thereby dividing the reconstructable region into several topologically independent terrain sub-regions.

[0015] By projecting structural lines onto a regular grid and deleting node adjacencies across structural lines, topological partitioning of the reconstructable region is achieved. This method ensures that subsequent smoothing constraints are established only within the same terrain sub-region, guaranteeing independent reconstruction of different structural regions such as the sides of the hierarchical platform, inside and outside the drainage ditch, and above and below the retaining wall. This fundamentally avoids the drawback of traditional continuous surface interpolation smoothing across structural boundaries, which weakens key terrain features.

[0016] Furthermore, the objective function Represented as: , The discrete calculation formulas for each constraint term are as follows: Boundary connection constraints : Used to constrain the connection between the reconstructed ground surface and the surrounding measured ground points; the calculation formula is: , in, This represents the set of meshes at the boundary of the reconstructable region. The boundary grid elevation is to be determined. The known measured elevations at the boundary; Topographic trend constraints Used to limit the deviation of the reconstructed surface. The overall trend is provided, and the calculation formula is as follows: , In the formula, For reconstructable areas, For the desired grid elevation, For the corrected The elevation value at this grid point; Local smoothing constraints Used to maintain reasonable continuity within the same structural region, the calculation formula is: , in, To avoid crossing terrain structure lines The set of adjacent grid pairs; and To balance the weighting coefficients, To and Elevation of adjacent grid nodes to be determined.

[0017] By constructing an objective function comprising three constraints—boundary connectivity, DEM trend, and local smoothing—the reconstructed surface is simultaneously constrained by fixed measured elevations at the boundary, the macroscopic trend of the existing DEM, and continuous smoothing within sub-regions. The contributions of each constraint under different data conditions are adjusted through a balancing weight coefficient. This ensures seamless integration between the reconstructed results and surrounding measured data, limits deviations of the reconstructed surface from the macroscopic topographic trend, and maintains reasonable continuity within sub-regions, thus achieving global optimization of the reconstructed surface.

[0018] Furthermore, the determined reconstructable region The specific method is as follows: Divide the candidate missing region into eight directional sectors, and statistically analyze the directional distribution of the surrounding corrected ground points within the set search distance. When the surrounding corrected ground points are distributed in no less than three directional sectors, or distributed in two opposite directional sectors and the difference between the surrounding slope aspect and the existing DEM slope aspect does not exceed the set threshold, the region is identified as a reconstructable region; otherwise, it is marked as a non-reconstructable region and null values ​​are retained.

[0019] By dividing the missing area into eight directional sectors and statistically analyzing the directional distribution of surrounding ground points, the system quantitatively determines whether the area meets multi-directional constraints. Simultaneously, it introduces verification based on slope aspect consistency with the existing DEM to ensure that the surrounding topographic trend does not conflict with the macro-topography. Only areas satisfying multi-directional constraints or relative direction plus DEM consistency constraints are considered reconstructable areas. Areas with insufficient constraints or conflicting trends are left null, fundamentally eliminating topographic errors caused by unfounded interpolation in areas with insufficient constraints, a problem inherent in traditional methods.

[0020] Furthermore, the initial ground filtering adopts cloth simulation filtering, progressively encrypted irregular triangular mesh filtering, or progressive morphological filtering methods; the establishment of local terrain reference surfaces includes: when the number of initial ground points in a certain scale neighborhood is lower than the fitted set threshold, or when the initial ground points are only concentrated in a single direction, the scale is marked as insufficiently constrained and does not participate in the ground point back judgment.

[0021] By explicitly stating that various existing filtering methods can be used for initial filtering, and by setting thresholds for the number of neighborhood points and directional distribution criteria in multi-scale analysis, scales with insufficient constraints are automatically marked and excluded. This approach makes the entire method independent of any specific filtering algorithm, adaptable to application scenarios with different data quality and terrain features, while ensuring that only scales with sufficient constraint information participate in ground point re-judgment, thus improving the method's universality and the reliability of the judgment.

[0022] Furthermore, the specific method for generating the digital elevation model is as follows: modifying the ground point set... With reconstructing surface point clouds Integration, based on terrain structure lines As a forced edge to construct a constrained irregular triangular mesh, the triangular mesh is sampled to a regular grid through bilinear interpolation; the grid marked as non-reconstruction area is kept null and the existing DEM is not used for high-resolution filling.

[0023] By constructing a constrained triangular mesh using structural lines as forced edges, the terrain structural lines are accurately represented in the final DEM. By keeping null values ​​in non-reconstruction areas without high-resolution filling, the propagation of terrain errors caused by forced interpolation in areas that do not meet the reconstruction conditions, as in traditional methods, is avoided, thus ensuring the authenticity and reliability of the DEM results in each region.

[0024] This application also includes a vegetation-shaded slope DEM reconstruction system that integrates multi-source remote sensing data, and implements a method for vegetation-shaded slope DEM reconstruction that integrates multi-source remote sensing data, including: The data acquisition module is used to acquire UAV lidar point clouds, orthophotos, and existing digital elevation models. The data preprocessing module is used to unify multi-source data into the same geographic coordinate system and perform initial ground filtering to obtain an initial ground point set and a non-ground point set. The ground point reclassification module is used to perform multi-scale terrain continuity analysis on non-ground points, reclassify non-ground points that meet the conditions as ground points, form a corrected ground point set and extract terrain structure lines. The reconstructable area determination module is used to determine the reconstructable area based on the directional distribution of ground points around the missing area and its consistency with the existing DEM topographic trend. The structural constraint reconstruction module is used to divide the terrain into sub-regions using terrain structure lines, establish and solve an objective function that includes boundary elevation constraints, DEM trend constraints and local smoothing constraints, and obtain the reconstructed surface point cloud. The DEM generation module is used to fuse and correct ground point sets with reconstructed surface point clouds to generate a digital elevation model.

[0025] Compared with existing technologies, the advantages of this invention are: 1. This application identifies and back-judges real ground points with multi-scale continuous relationships with adjacent terrain from the initial non-ground point set, which makes up for the defect of the initial filtering in high and steep slope areas that easily misclassifies real surface structure points such as bare rock surfaces, steep slopes, ridges and slope breaks as non-ground points, and provides a more complete ground point data foundation for subsequent structure line extraction and DEM generation.

[0026] 2. This application determines the reconstructable area based on the directional distribution and distance of reliable ground points around the missing area and their consistency with the existing DEM topographic trend. For areas with insufficient constraints or obvious conflicts in topographic trends, null values ​​are retained instead of being forcibly filled, thus avoiding the unfounded topographic distortion caused by the traditional method of uniformly interpolating all sparse areas.

[0027] 3. This application uses structural lines such as ridges, toes, steep slopes and platform edges to divide the terrain into sub-regions and eliminates the grid connection relationship across structures, so that the surface reconstruction is carried out only within the same terrain sub-region. This avoids the reconstruction surface from crossing different terrain structures such as graded platforms, drainage ditches and retaining walls to make it smooth, which is conducive to maintaining the morphological characteristics of key terrain structures and artificial engineering structures on the slope. Attached Figure Description

[0028] Figure 1 This is a flowchart of a method for reconstructing a vegetation-shaded slope DEM by integrating multi-source remote sensing data.

[0029] Figure 2 Figure 1 shows the orthophoto image of the target side wave, the lidar point cloud, and the schematic diagram of the control point layout. Figure 2 shows the orthophoto image of the slope, Figure 3 shows the overall distribution of the lidar point cloud, and Figure 4 shows the schematic diagram of the control point layout.

[0030] Figure 3 This is a comparative diagram of the initial ground point set and the initial non-ground point set in this application; wherein, Figure (a) shows the initial ground point set. Schematic diagram; Figure (b) shows the initial non-ground point set. Schematic diagram.

[0031] Figure 4 This is a schematic diagram illustrating the process of candidate ground point re-judgment, reconstructable region determination, and structural constraint DEM generation in this application. Detailed Implementation

[0032] It should be noted that relational terms such as "first" and "second" are used merely to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element.

[0033] The features and performance of the present invention will be further described in detail below with reference to embodiments.

[0034] Please see Figures 1-4 A method for reconstructing DEMs of vegetation-shaded slopes by fusing multi-source remote sensing data, such as... Figure 1 As shown, it includes the following steps: Multi-source data acquisition: Acquire UAV lidar point clouds and corresponding orthophotos covering the target slope; acquire existing digital elevation models covering the target slope to provide macroscopic topographic trends of the slope; for areas with insufficient UAV point cloud coverage or requiring accuracy verification, acquire supplementary data from ground 3D laser scanning, total station or other ground measurement methods.

[0035] Multi-source data preprocessing and spatial unification: Convert all data sources to a unified geographic coordinate system and elevation datum, and establish a unified spatial index; use UAV integrated navigation data, ground control points or known coordinates to complete the initial registration of multi-source data; use ground filtering methods to perform initial classification of UAV point clouds to obtain initial ground point sets and initial non-ground point sets.

[0036] Candidate ground structure point identification and reclassification: Based on the median nearest neighbor distance of the initial ground point set, a three-dimensional neighborhood with three progressive scales (small, medium, and large) is established; at each effective scale, a local terrain reference surface is established using the initial ground points in the neighborhood, and the multi-scale geometric features of non-ground points are calculated; for non-ground points that meet the terrain continuity condition at at least two effective scales, they are reclassified as ground points and incorporated into the initial ground point set to form a modified ground point set; terrain structure lines are extracted based on the normal abrupt change and curvature abrupt change of the modified ground points.

[0037] Determining the area to be reconstructed: Divide the corrected ground point set into regular grids and count the ground point density and spatial coverage of each grid; extract candidate sparse and missing areas by combining orthophoto vegetation cover information; determine the reconstructable areas based on the directional distribution, distance and consistency with the existing DEM topographic trend of the missing areas, and retain the areas that do not meet the conditions as null values.

[0038] Structurally constrained surface reconstruction: The reconstructable area is divided into several topologically independent sub-regions using topographic structure lines; an objective function is established in each sub-region, including boundary elevation constraints, existing DEM trend constraints, and local continuity constraints; the reconstructed surface point cloud is obtained by solving the objective function.

[0039] Vegetation occlusion area identification: Based on the local point density, neighborhood elevation continuity, and orthophoto vegetation cover information of the fused ground point cloud, identify areas with missing surface points caused by vegetation occlusion, remove non-vegetation occlusion areas such as man-made structures, and determine the areas to be restored. .

[0040] DEM generation: The corrected ground point set is merged with the reconstructed surface point cloud, and an irregular triangular network or regular grid is constructed using terrain structure lines as forced edges to generate a digital elevation model of vegetation-shaded slopes; areas that do not meet the reconstruction conditions are retained as null values.

[0041] The present application will now be described clearly and completely in conjunction with the technical solutions and accompanying drawings in the embodiments of the present application.

[0042] Example 1: Scenario of a steep rock-cut slope: This embodiment uses a vegetation-covered road cut slope on a highway in Yunnan Province as an example to illustrate the specific implementation of the present invention. The selected slope is approximately 150m long along the route, with a maximum slope height of approximately 50m and a slope angle of approximately 45° to 65°. The slope surface includes graded platforms, drainage ditches, local bare rock steep slopes, and areas covered by trees and shrubs. The vegetation coverage is approximately 65%, which is a typical high and steep rock-soil mixed slope under vegetation cover conditions.

[0043] Multi-source data acquisition: UAV LiDAR and Optical Image Data Acquisition: A UAV equipped with a LiDAR and optical camera is used to acquire data on the target slope, obtaining UAV LiDAR point clouds. and digital orthophotos The drone platform used was the DJI M300 RTK, equipped with a Zenmuse L1 LiDAR (wavelength 905nm, maximum range 450m @ reflectivity 80%, range accuracy ±3cm, point frequency 240kHz, scanning angle range 70°×30°) and a P1 full-frame camera (45 million effective pixels, pixel size 4.4μm, equipped with a 35mm fixed-focus lens). Flight parameters were set as follows: relative flight altitude 120m (based on the top of the slope), flight speed 8m / s, lateral overlap 80%, lateral overlap 60%, and flight paths parallel to the slope direction, with a total of 4 parallel flight paths. Quality checks showed that the airborne point cloud density was approximately 80 points / m², with a point cloud coverage area of ​​approximately 0.08km²; the orthophoto ground resolution was approximately 0.05m.

[0044] Ground-based supplementary data acquisition: Obtaining an existing digital elevation model with a 12.5m grid covering the target slope. The DEM is derived from the 1:10000 basic survey results of the region and serves as prior data for the macro-topographic trend of the slope.

[0045] Twenty Global Navigation Satellite System (GNSS) control points were deployed in accessible areas such as the foot of the slope and the graded platform, of which 12 were used for spatial registration and 8 served as independent checkpoints.

[0046] In areas with severe vegetation obstruction at the foot of the slope, a Faro Focus S150 terrestrial 3D laser scanner was used to acquire supplementary ground point clouds. ,like Figure 2 As shown. The main parameters of the scanner are: ranging distance 150m @ reflectivity 90%, ranging accuracy ±1mm, angle accuracy ±19″, scanning speed 2 million points / second, field of view 300° horizontally and 360° vertically. The scanning parameters are set as follows: scanning resolution 1 / 2 (i.e., angular interval of 0.036° per point), quality coefficient 3×, and scanning time per station is about 8 minutes. A total of 7 ground laser scanner stations were set up at the bottom of the slope (3 stations), at each level of platform (2 stations), and in areas with dense vegetation cover (2 stations), with a station spacing of about 50-80m. The ground point cloud density is about 5000 points / m², and the point cloud coverage area within the effective scanning radius is about 0.02km².

[0047] Multi-source data preprocessing and spatial unification: All data sources (drone laser point cloud) Ground-based supplemental point clouds Digital orthophotos Existing digital elevation models The coordinates were uniformly converted to the 2000 National Geodetic Coordinate System (CGCS2000, ellipsoid parameters: a=6378137m, f=1 / 298.257222101) and the Gauss-Kruger projection (central meridian 108°E, 3° zone, zone number 36) system, and the elevation datum adopted was the 1985 National Elevation Datum.

[0048] Multi-source data spatial registration: Initial spatial registration of multi-source data is completed using deployed GNSS control points. The UAV integrated navigation system (GNSS / IMU) employs three satellites with seven frequencies (GPS L1 / L2, GLONASS L1 / L2, BDS B1 / B2), with an IMU sampling frequency of 200Hz, attitude accuracy of 0.015°, and heading accuracy of 0.08°. The ground-based 3D laser scanner, equipped with dual-frequency GNSS (horizontal accuracy ±5mm +0.5ppm, vertical accuracy ±10mm +0.5ppm), performs resection with the total station (angle measurement accuracy ±2″, ranging accuracy ±2mm +2ppm) to complete the initial spatial registration of each station. The registration residual for each station is controlled within ±5mm.

[0049] Corresponding points were selected in stable, non-vegetated areas such as road platforms and bare rock surfaces. Fine-grained registration was performed using the Iterative Closest Point (ICP) method from point to plane, with a maximum corresponding distance of 1.0m, a maximum number of iterations of 100, and a pose parameter change of less than 10. -6 The iteration was stopped at that point. Verification showed that the initial spatial registration error between the airborne point cloud and the ground point cloud was approximately 0.15m.

[0050] Unified Spatial Index Construction: Construct an octree unified spatial index structure, setting the minimum voxel side length (leaf node size) to 0.5m and the maximum recursion depth to 10 levels. This index structure supports... Achieve fast neighborhood retrieval and spatial association of multi-source data within time complexity. This represents the total number of point clouds.

[0051] Initial ground filtering classification: Noise points are removed using a statistical outlier method. Specifically, each point searches for 20 neighboring points, and points whose distance exceeds twice the standard deviation of the neighborhood average distance are considered outliers and removed.

[0052] The cloth simulation filtering (CSF) algorithm was used for initial ground classification. The cloth grid resolution was set to 0.5m (consistent with the node size of the octagonal leaf), the classification distance threshold (rigidity threshold) was set to 0.30m, the cloth stiffness parameter was set to 2 (medium rigidity, suitable for vegetated slopes), the maximum number of iterations was set to 500, the time step was 0.65 (to ensure stable convergence of the cloth simulation), and the gravitational acceleration was 9.8m / s². The filtered initial airborne ground point set was then output. and the initial non-ground point set ,like Figure 3 As shown. After filtering, about 45% (about 2.7 million points) of the airborne point cloud were classified as ground points, and about 55% were classified as non-ground points (mainly vegetation canopy points); about 85% (about 4.25 million points) of the ground supplementary point cloud were classified as ground points, and about 15% were a small number of low vegetation points or noise points.

[0053] Candidate ground structure point identification and back-judgment: This step solves the problem that the initial ground filtering may misclassify real surface points such as steep slopes, bare rock surfaces, slope ridges, slope toes and steep slopes as non-ground points.

[0054] Multi-scale neighborhood construction: Calculating the initial ground point set Median of nearest neighbor distances for each point According to statistics, in this embodiment... Non-ground points to be determined Centered on, establish radii of and respectively. , and The scale progression relationship of the three three-dimensional spherical neighborhoods is as follows: , , , Among them, the proportionality coefficient satisfies : Small-scale neighborhood ( ): Used to determine non-ground points Local spatial connections with neighboring ground points; Mesoscale neighborhood ( ): Used to determine the local geometric extension and continuity of a slope; Large-scale neighborhood ( ): Used to determine non-ground points Does it conform to the overall topographical trend of the slope?

[0055] In this embodiment, (Approximately 0.36m) (approximately 0.72m) and Three three-dimensional spherical neighborhoods (approximately 1.44m) are defined. These three scales are used to determine local spatial connectivity (small scale), local geometric extension of the slope (medium scale), and overall topographic trend consistency of the slope (large scale), respectively. If the number of initial ground points in a neighborhood at a certain scale is less than 12, or if the initial ground points are concentrated only in a single-directional sector (i.e., the number of points in the largest directional sector accounts for more than 80% of the total), then that scale is marked as insufficiently constrained and will not participate in subsequent judgments.

[0056] Local terrain reference surface establishment: For each effective neighborhood, initial ground points are extracted, and a local terrain reference surface is established using the Random Sample Consensus (RANSAC) plane fitting method. The maximum number of samplings is set to 200, and the interior point distance threshold is set to... .

[0057] In the Each scale ( The implicit equation for the local terrain reference surface under () is expressed as: ; in The unit normal vector of the reference surface. It is a plane constant. is the coordinate vector of a point in space.

[0058] Calculate the fitting residual of the local topographic reference surface : ; in This represents the initial number of ground points used in the fitting process. For the first A spatial point coordinate vector.

[0059] When the fitting residual of the local terrain reference surface Less than the threshold (In this embodiment) When the fitted residual is large and the normal of the ground points in the neighborhood shows a significant bimodal distribution, the ground points in the neighborhood are divided into two point sets to fit two local reference surfaces. and It is used to determine the structure of subsequent slope bends, ridges, toes, or steep slopes.

[0060] Multi-scale geometric feature calculation for non-ground points: For non-ground points to be judged In each effective scale Calculate the following geometric features: (1) Distance from non-ground points to local topographic reference surfaces If the neighborhood is a single continuous slope, then: ; If there are two local reference surfaces in the neighborhood, then calculate the non-ground points separately. Distances to two reference planes and taking the minimum value: ; Distance adaptive tolerance threshold Based on the median of the nearest neighbor distances of the initial ground points Fitting residuals with local topographic reference surface Dynamically determined: In this embodiment , .

[0061] in, and This is an adjustment factor. Flat and dense point cloud regions correspond to a smaller threshold, while coarse or sparse regions automatically have a wider threshold.

[0062] (2) Local normal angle At non-ground points Construct the covariance matrix in the neighborhood Then perform eigenvalue decomposition to calculate non-ground points. Local normal vector Then calculate the local normal vector. unit normal vector to the reference plane The spatial angle is the local normal angle. : ; In this embodiment, the normal angle threshold is set to 20°. The smaller the value, the more likely it is to be a non-ground point. The more consistent the orientation of the local micro-surface with the adjacent ground, the better.

[0063] (3) Local curvature Let the eigenvalues ​​of the neighborhood covariance matrix be... ,but In this embodiment, the curvature threshold is set to 0.10. Points with small curvature and gentle changes in normal are selected as candidate points for continuous slopes; points with large curvature or abrupt changes in normal are retained in the logic for determining slope break points and are not directly identified as noise.

[0064] (4) Ground connectivity. Establish local adjacency relationships based on the distance between points and the included normal angle, and determine the non-ground points. Is it possible to connect a continuous set of points to an initial set of ground points? Isolated point clusters that are only connected to other non-ground points and cannot be connected to known ground points are not directly reclassified as ground points.

[0065] (5) Land cover category. Non-ground points... Projected onto orthophoto This process retrieves the land cover category for the corresponding location. The image classification results are used to help exclude non-ground targets such as vegetation canopies, buildings, vehicles, and pole-shaped structures.

[0066] When the image location is within a vegetation-covered area, the judgment is made in conjunction with the vertical layering relationship of the point cloud: points located inside the vegetation canopy or branch clusters remain as non-ground points; low-level echo points located below vegetation clusters and with a continuous relationship to the adjacent ground can still be included in the ground point re-judgment.

[0067] (6) Spatial Relationship of Engineering Structures. The road platform, the bottom of the drainage ditch, and the surface of the protective works that are continuously connected to the slope are retained as part of the ground surface; the vertical facade of the retaining wall, guardrails, vehicles, and pole-like facilities are not incorporated into the ground point set. The top of the retaining wall, the foot of the wall, the edge of the platform, and the edge of the drainage ditch serve as the structural boundaries for subsequent ground surface reconstruction.

[0068] Ground point return judgment: when it is not a ground point A point is identified as a candidate point for a continuous slope when all of the following conditions are met simultaneously at at least two effective scales: (1) The distance from the point to the local terrain reference surface is not greater than the corresponding adaptive distance threshold; (2) The angle between the local normal of a point and the normal of the reference surface is no greater than 20°; (3) The direction of elevation change of the point is consistent with the direction of extension of the adjacent slope; (4) Points can be connected to the initial ground point set through a continuous set of points. ; (5) The point does not belong to the vegetation canopy, vehicles, guardrails, pole-shaped facilities or buildings.

[0069] Connectivity analysis is performed on candidate points to eliminate point clusters with fewer than 20 points, spatially discrete points, or points lacking neighborhood support, thus obtaining the set of ground points for back-judgment. Incorporate it into the initial ground point set. The corrected ground point set is obtained: ; according to By connecting candidate points of the structure based on the normal abrupt change, local maxima of curvature, and the overlap of the two reference planes, discrete line segments with a length less than 1.0m are eliminated to obtain the structure line. ,like Figure 4 As shown. The vertical facade of the retaining wall, railings, and pole-like structures are not included in the ground point set, and the top and base of the retaining wall can be used as structural boundaries.

[0070] Areas to be rebuilt identified: Extraction of sparse and missing regions: The ground point set is corrected using a grid resolution of 0.5m. The grid is divided into regular grids. The median density of grid points in the exposed area is used as the reference density (approximately 85 points / m² in this embodiment). Grids with no ground points or a point density lower than 25% of the reference density (approximately 21 points / m²) are marked as candidate sparse or missing grids.

[0071] Exclusion of non-surface areas: combining orthophotos Based on vegetation cover information and artificial structure identification results (based on deep learning semantic segmentation, with an accuracy of approximately 92%), buildings, retaining wall facades, road facilities, and other non-surface areas are excluded to obtain candidate areas for reconstruction.

[0072] Determining Reconstructable Regions: Divide the candidate missing region into 8 directional sectors (each sector 45°). Count effective surface points within a 5.0m search distance. If reliable surface points are distributed across at least 3 directional sectors, it is considered to meet multi-directional constraints; if they are distributed across two opposing directional sectors (difference 180°±45°) and the surrounding slope aspect is consistent with the existing digital elevation model... When the slope difference is ≤20°, it is determined that the relative direction constraint condition is met.

[0073] Areas that meet any of the above conditions and do not exhibit significant conflicts in local topographic trends are identified as reconstructable areas. ,like Figure 4 As shown; otherwise, mark it as a non-reconstructable area and leave it blank. In this embodiment, the reconstructable area accounts for approximately 22% of the total slope area.

[0074] Structurally constrained surface reconstruction: Topographic sub-region division: dividing structural lines Projected onto the reconstructable area Regular grid (grid resolution) ), and uses four-neighborhood (up, down, left, and right directions) to establish node adjacency relationships. When the line connecting two adjacent nodes intersects with the line connecting the two nodes, the node is considered an adjacency. When they intersect, delete the adjacency relationship, thus... Divided into several topologically independent terrain sub-regions, such as Figure 4 As shown.

[0075] Boundary constraints and trend correction: A width of [value missing] is applied to the perimeter of each sub-region. Extract corrected ground points within a buffer zone of 2.5m. Calculate the elevation of ground points and The median of the elevation residuals at the same location, and corrected accordingly. Trend elevation, resulting in the corrected trend surface An inverse distance weighting method (weight index) with a search radius of 2h (i.e., 1.0m) is adopted. According to the surrounding area Determine the boundary observation elevation.

[0076] Objective function establishment: Using grid node elevations as unknowns, establish an objective function that includes boundary elevation constraints, trend constraints, and local smoothing constraints. : ; The discrete calculation formulas for each constraint term are as follows: Boundary connection constraints : Used to constrain the connection between the reconstructed ground surface and the surrounding measured ground points; the calculation formula is: , in, This represents the set of meshes at the boundary of the reconstructable region. The boundary grid elevation is to be determined. The known measured elevations at the boundary are given.

[0077] Topographic trend constraints Used to limit the deviation of the reconstructed surface. The overall trend is provided, and the calculation formula is as follows: , in, For reconstructable areas, For the desired grid elevation, For the corrected The elevation value at this grid point.

[0078] Local smoothing constraints Used to maintain reasonable continuity within the same structural region, the calculation formula is: , in, To avoid crossing terrain structure lines The set of adjacent grid pairs; and These are the balancing weighting coefficients.

[0079] For grids separated by ridges, steep slopes, platform edges, or retaining wall boundaries, the connecting lines cross terrain structure lines. Therefore, it is not included in the adjacent grid set. This means that a smooth relationship is not established across structural lines to avoid unconditional interpolation of the reconstructed surface across different structural regions. Each item is set after normalization based on the number of nodes. , , Objective function solution: The preconditioned conjugate gradient (PCG) method is used to solve the above objective function. The preconditioning matrix uses incomplete Cholesky decomposition (IC(0)) to improve convergence speed. The convergence condition is set as the relative change in the objective function value between two adjacent iterations being less than 10. -6 (Right now The maximum number of iterations is 1000. It converges after approximately 120 iterations, outputting the surface point cloud of the reconstructable region. .

[0080] DEM generation: Constrained DEM generation: Correcting the ground point set With reconstructing surface point clouds Integration, with structural lines A constrained delaunay triangulation (DEM) was constructed using forced edges. The maximum triangle side length was limited to 3m (to prevent excessively large triangles in sparse data areas), and the maximum triangle angle was limited to 120° (to prevent overly sharp triangles from affecting interpolation quality). The TIN elevation was sampled to a 0.5m regular grid using bilinear interpolation to generate a digital elevation model (DEM) for vegetation-shaded slopes. For areas that did not meet the reconstruction criteria, null values ​​(NoData) were retained, and high-resolution filling was performed directly without using existing DEMs.

[0081] Accuracy verification: The mean absolute error (MAE) and root mean square error (RMSE) were calculated using 8 independent checkpoints and ground points reserved from the exposed area (a total of 30 verification points), and compared with those using only the initial ground point set. The results of ordinary triangulation interpolation were compared. Checkpoints and reserved ground points were used only for verification and did not participate in registration, ground point re-identification, or surface reconstruction. Verification showed that the final DEM RMSE at the verification points in this embodiment was approximately 0.16m, an improvement of approximately 62% compared to the method using only initial filtered ground points for ordinary TIN interpolation (RMSE approximately 0.42m).

[0082] Output: Outputs a high-resolution digital elevation model (DEM) of the slope with a spatial resolution of 0.5m. Elevation values ​​are stored in GeoTIFF floating-point format. Auxiliary outputs such as slope maps, aspect maps, and curvature maps are also generated.

[0083] The above implementation process can be achieved through point cloud spatial indexing, random sampling consistency plane fitting, connected component analysis, grid segment intersection judgment, and sparse linear equation solving program. The input and output data of each step are the same as those of the above steps. Figure 1 The overall process shown corresponds to this. The structural lines extracted in the candidate ground structure point identification and back-judgment step... The data is passed to the structurally constrained surface reconstruction step for terrain subregion division and smooth blocking across structural lines; the candidate ground structure point identification and back-judgment step generates a corrected ground point set. Passed to the DEM generation step, and The slope DEM is generated by fusion.

[0084] Example 2: High embankment slope scenario: The difference between this embodiment and Embodiment 1 lies in the application scenario and some parameter settings.

[0085] This embodiment uses a high embankment slope as the application scenario. The slope is about 300m long along the route, with a maximum embankment height of about 28m and a slope of about 30° to 40°. The embankment is mainly composed of gravelly soil, with a rough surface and a large number of gravel particles with a diameter of 0.05m to 0.2m. The vegetation coverage is about 40%, and the main vegetation types are low-growing herbs and sparse shrubs. The slope is equipped with a three-tiered platform and a longitudinal drainage ditch. The significant differences between this scenario and Embodiment 1 are: (1) the terrain is relatively gentle (the slope is about 60% of that of a rock slope); (2) the vegetation coverage is low and mainly low-growing plants (which have little impact on point cloud penetration); (3) the surface is loose deposits, and the geometric features are not as regular as those of a rock slope.

[0086] The drone's flight altitude was adjusted to 80m (the fill slope is relatively low), the flight speed was adjusted to 6m / s, the forward overlap rate was maintained at 80%, the lateral overlap rate was maintained at 60%, the airborne point cloud density was approximately 120 points / m², and the ground resolution of the orthophoto was set to 0.03m.

[0087] The number of ground-based 3D laser scanner stations was reduced to four (due to the relatively open terrain of the fill site with less obstruction), with a station spacing of approximately 100-150m. The ground point cloud density was approximately 3000 points / m². Twelve GNSS control points were deployed, eight for registration and four for independent inspection.

[0088] The cloth simulation filter (CSF) has its cloth mesh resolution adjusted to 0.3m (the surface undulation of the fill slope is small, and a finer cloth mesh can better fit the surface details), the classification distance threshold (rigid threshold) adjusted to 0.20m (the surface roughness of loose deposits is high, and a more flexible cloth can better track the surface undulations), and the maximum number of iterations is kept at 500.

[0089] Median of the nearest neighbor distances of the initial ground point The three neighborhood radii are respectively set to (Approximately 0.24m) (approximately 0.48m) and (Approximately 0.96m). The threshold for the number of initial ground points in the neighborhood at each scale was adjusted to no less than 8 (the ground points on the fill slope are relatively evenly distributed).

[0090] The interior point distance threshold for RANSAC plane fitting is set to .

[0091] The normal angle threshold is adjusted to 25° (due to the high surface roughness of the fill slope), and the curvature threshold is adjusted to 0.15.

[0092] The minimum segment length threshold for structural line extraction was adjusted to 0.5m (due to the smaller size of the fill slope structure).

[0093] The resolution of the regular grid is adjusted to 0.3m, and the reference density threshold is adjusted to 20% of the reference density (approximately 100 points / m²) (i.e., 20 points / m²).

[0094] The search distance for determining the reconstructable area is adjusted to 3.0m (for areas with relatively gentle slope changes), and the aspect difference threshold is kept at 20°.

[0095] Grid resolution The outer buffer width of each sub-region is 5h (i.e., 1.5m).

[0096] After normalizing each term of the objective function according to the number of nodes, set the parameters accordingly. , (The existing DEM trend reliability of the fill slope is moderate, so the trend constraint strength should be appropriately reduced.) The convergence threshold should be kept at 10. -6 The maximum number of iterations is 1000. The DEM grid resolution is set to 0.3m.

[0097] The remaining steps are the same as in Embodiment 1, and will not be repeated here. This embodiment shows that, under different types of slope scenarios, this application can effectively achieve high-precision reconstruction of slope DEM under vegetation shading conditions by adjusting the corresponding filtering parameters, neighborhood scale, and judgment threshold.

[0098] Example 3: Densely vegetated slope along the reservoir bank: The difference between this embodiment and Embodiment 1 lies in the application scenario and some parameter settings.

[0099] This embodiment uses a reservoir bank slope of a water conservancy project as the application scenario. The slope is approximately 1200m long, 60m high, and has a slope of approximately 30°–50°. The vegetation coverage is approximately 70%, mainly consisting of a mixed forest of trees (primarily pine and eucalyptus, with tree heights of 8m–15m and canopy closure of approximately 0.75) and shrubs. The vegetation cover density is significantly higher than in Embodiment 1. Typical characteristics of the reservoir bank slope include: the slope is affected by the periodic rise and fall of the reservoir water level (water level fluctuation of approximately 10-15m), the surface is covered with a thick layer of residual slope soil, the soil has a high water content (volume water content of approximately 25%–35%), the vegetation has a well-developed root system, and surface point cloud acquisition is difficult. This scenario is used to verify the applicability of the invention under extreme vegetation obstruction and complex surface conditions.

[0100] Due to the high vegetation coverage (approximately 70%) on the reservoir bank slopes, the UAV's flight altitude was adjusted to 100m (20m lower than in Example 1 to increase the probability of point clouds penetrating the vegetation canopy), and the flight speed was adjusted to 5m / s (to increase the number of scan lines per unit area). The forward overlap rate was adjusted to 85%, and the lateral overlap rate was adjusted to 65%, increasing the airborne point cloud density to approximately 60 points / m² to obtain more point cloud information below the vegetation canopy. The lidar pulse emission frequency remained unchanged at 240kHz, but by reducing the flight speed, the actual ground point cloud density increased from 50 points / m² to 60 points / m².

[0101] The number of ground-based 3D laser scanner stations has increased to 10 (due to severe obstruction by dense vegetation, more ground stations are needed to ensure coverage), with a station spacing of approximately 30m to 50m. The ground point cloud density is approximately 8000 points / m² (due to the increased station density, the scanning resolution per station can be appropriately reduced, using a scanning resolution of 1 / 2 and a quality coefficient of 2×, with a single station scanning time of approximately 5 minutes, to balance data acquisition efficiency).

[0102] The cloth simulation filtering (CSF) was adjusted with a cloth mesh resolution of 0.5m and a rigidity threshold of 0.6 (due to the low penetration rate of point clouds in densely vegetated areas and the small proportion of ground points, a more rigid cloth parameter (rigidity threshold of 0.6) can effectively cross the gaps in the vegetation canopy, preventing the cloth from getting stuck inside the vegetation canopy and reducing the misclassification of vegetation points as ground points), and the maximum number of iterations was kept at 500. After filtering, approximately 30% of the points in the airborne point cloud were classified as ground points (due to the low penetration rate caused by dense vegetation), and approximately 70% were vegetation canopy points.

[0103] The neighborhood radius of DBSCAN clustering is adjusted to (Trees with large canopies (canopy diameter of about 3m to 6m) require a larger neighborhood radius to correctly cluster the complete canopy point cloud), the minimum number of neighborhood points MinPts is adjusted to 15 (to suppress small noise clusters).

[0104] The resolution of the regular grid remains at 0.5m, and the reference density threshold is adjusted to 20% of the reference density (approximately 55 points / m²) (i.e., 11 points / m²).

[0105] The search distance for reconstructable areas has been adjusted to 8.0m (dense vegetation makes the distribution of known ground points sparser, requiring a larger search range).

[0106] Grid resolution The outer buffer width of each sub-region is 5h (i.e., 2.5m).

[0107] After normalizing each term of the objective function according to the number of nodes, set the parameters accordingly. , (In densely vegetated areas, DEM trend constraints need to be strengthened to compensate for the lack of point cloud information, and the smoothing constraint strength should be appropriately reduced to avoid over-smoothing).

[0108] The preprocessed conjugate gradient method (PCG) converged after about 200 iterations (the area to be restored accounts for about 40% of the total slope area, about 0.24 km², and the solution scale is about 3.2 times that of Example 1).

[0109] The remaining steps are basically the same as in Example 1. The purpose of this example is to demonstrate that even under extreme conditions with a vegetation coverage of up to 70%, this application can still stably complete surface reconstruction and DEM generation by adaptively adjusting the combination of coefficients, demonstrating broad scene adaptability and robustness.

[0110] The embodiments described above merely illustrate specific implementation methods of this application, and while the descriptions are detailed and specific, they should not be construed as limiting the scope of protection of this application. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the technical solution of this application, and these modifications and improvements all fall within the scope of protection of this application.

Claims

1. A method for reconstructing a vegetation-shaded slope DEM by fusing multi-source remote sensing data, characterized in that, Includes the following steps: Acquire multi-source remote sensing data, including UAV lidar point clouds, orthophotos, and existing digital elevation models; By unifying all data sources to the same geographic coordinate system, initial ground filtering is performed on the UAV point cloud to obtain the initial ground point set. and the initial non-ground point set ; A multi-scale three-dimensional neighborhood is established based on the median nearest neighbor distance of the initial ground point set. At each scale, a local terrain reference surface is established using the initial ground points within the neighborhood. Multi-scale terrain continuity analysis is performed on non-ground points. Non-ground points that meet the continuity condition at at least two scales are reclassified as ground points and incorporated into the initial ground point set to form a modified ground point set. And extract terrain structure lines based on the corrected ground point set. ; The corrected ground point set is divided into a regular grid, and candidate missing regions are extracted by combining orthophoto vegetation cover information. Reconstructable regions are determined based on the directional distribution of ground points around the missing regions and their consistency with the existing DEM topographic trend. ; Utilizing terrain structure lines The reconstructable area is divided into several terrain sub-regions. An objective function is established within each sub-region, incorporating boundary elevation constraints, existing DEM trend constraints, and local smoothing constraints. Solving this function yields the reconstructed surface point cloud. ; Correct the ground point set With reconstructing surface point clouds Integration, based on terrain structure lines As a forced edge, an irregular triangular network is constructed to generate a digital elevation model, and areas that do not meet the reconstruction conditions are retained as null values.

2. The method for reconstructing a vegetation-shaded slope DEM by fusing multi-source remote sensing data according to claim 1, characterized in that, The three-dimensional neighborhoods with multiple progressive scales are specifically: Calculate the initial ground point set Median nearest neighbor distance ,by As the distance between feature points, establish a radius of... , and The three three-dimensional spherical neighborhoods, where the scaling coefficients satisfy ; The local terrain reference surface is established as follows: when the number of initial ground points in the neighborhood meets the threshold and the orientation distribution is uniform, a robust weighted least squares or RANSAC plane fitting method is used to establish the local terrain reference surface. And calculate the fitting residual of the local terrain reference surface. ,in Let be the unit normal vector of the reference surface. It is a plane constant. The coordinate vector of the spatial point is used. When the fitting residual is less than the threshold, it is considered a single continuous slope. When the fitting residual is large and the normal of the points in the neighborhood shows a significant bimodal distribution, two local reference surfaces are fitted respectively for the judgment of the slope break structure.

3. The method for reconstructing a vegetation-shaded slope DEM by fusing multi-source remote sensing data according to claim 2, characterized in that, The multi-scale terrain continuity analysis includes: Calculate the non-ground points to be judged at each effective scale Distance to local topographic reference surface Local normal angle and local curvature Specifically: non-ground points Distance to local topographic reference surface : ， If there are two local reference surfaces in the neighborhood, then calculate the non-ground points separately. The distances to the two reference planes are taken, and their minimum values ​​are calculated: ， in, Let be the unit normal vector of the first local reference surface. Let be the plane constant of the first local reference surface. The unit normal vector of the second local reference surface. The plane constant of the second local reference surface; At non-ground points Construct the covariance matrix in the neighborhood And calculate non-ground points Local normal vector Local normal vector unit normal vector to the reference plane The spatial angle is the local normal angle. The local normal angle is expressed as: ； Let non-ground points The eigenvalues ​​of the neighborhood covariance matrix are The local curvature is then expressed as: ； Adaptive tolerance threshold for distance ,in and For adjustment coefficients, The median of the nearest neighbor distances of the initial set of ground points. The fitting residuals for the local terrain reference surface; When a non-ground point Simultaneously satisfying at least two effective scales The included normal angle is no greater than a set threshold, the curvature is no greater than a set threshold, and it can be connected to the initial ground point set through a continuous point set. At that time, the point will be reclassified as a ground point.

4. The method for reconstructing a vegetation-shaded slope DEM by fusing multi-source remote sensing data according to claim 1, characterized in that, The extraction of terrain structure lines The specific method is as follows: Based on the modified ground point set In this study, abrupt changes in normal direction and local maxima of curvature are used to identify candidate structural points. Abrupt changes in normal direction are determined by the rate of change of the angle between the normal vectors of points in the neighborhood, while local maxima of curvature are determined by the local curvature of each grid node. Non-maximum suppression extraction is performed; connected component analysis is conducted on candidate structural points to remove discrete line segments with lengths less than a set threshold, resulting in ridge lines, toe lines, upper and lower edges of steep slopes, and platform edge lines.

5. The method for reconstructing a vegetation-shaded slope DEM by fusing multi-source remote sensing data according to claim 1, characterized in that, The use of terrain structure lines The specific method for dividing the reconstructable region into several terrain sub-regions is as follows: Terrain structure lines Projected onto the reconstructable area The regular grid uses four-neighbor or eight-neighbor domains to establish node adjacency relationships. When the line connecting two adjacent nodes coincides with a terrain structure line... When two regions intersect, the adjacency relationship is deleted, thereby dividing the reconstructable region into several topologically independent terrain subregions.

6. The method for reconstructing a vegetation-shaded slope DEM by fusing multi-source remote sensing data according to claim 1, characterized in that, The objective function Represented as: ， The discrete calculation formulas for each constraint term are as follows: Boundary connection constraints : Used to constrain the connection between the reconstructed ground surface and the surrounding measured ground points; the calculation formula is: ， in, This represents the set of meshes at the boundary of the reconstructable region. The boundary grid elevation is to be determined. The known measured elevations at the boundary; Topographic trend constraints Used to limit the deviation of the reconstructed surface. The overall trend is provided, and the calculation formula is as follows: ， in, For reconstructable areas, For the desired grid elevation, For the corrected The elevation value at this grid point; Local smoothing constraints Used to maintain reasonable continuity within the same structural region, the calculation formula is: ， in, To avoid crossing terrain structure lines The set of adjacent grid pairs; and To balance the weighting coefficients, To and Elevation of adjacent grid nodes to be determined.

7. The method for reconstructing a vegetation-shaded slope DEM by fusing multi-source remote sensing data according to claim 1, characterized in that, The area that can be reconstructed is determined. The specific method is as follows: Divide the candidate missing region into eight directional sectors, and statistically analyze the directional distribution of the surrounding corrected ground points within the set search distance. When the surrounding corrected ground points are distributed in no less than three directional sectors, or distributed in two opposite directional sectors and the difference between the surrounding slope aspect and the existing DEM slope aspect does not exceed the set threshold, the region is identified as a reconstructable region; otherwise, it is marked as a non-reconstructable region and null values ​​are retained.

8. The method for reconstructing a vegetation-shaded slope DEM by fusing multi-source remote sensing data according to claim 1, characterized in that, The initial ground filtering adopts cloth simulation filtering, progressively encrypted irregular triangular mesh filtering, or progressive morphological filtering methods; the establishment of local terrain reference surfaces includes: when the number of initial ground points in a certain scale neighborhood is lower than the fitted set threshold, or when the initial ground points are only concentrated in a single direction, the scale is marked as insufficiently constrained and does not participate in the ground point back judgment.

9. The method for reconstructing a vegetation-shaded slope DEM by fusing multi-source remote sensing data according to claim 1, characterized in that, The specific method for generating the digital elevation model is as follows: modify the ground point set... With reconstructing surface point clouds Integration, based on terrain structure lines As a forced edge to construct a constrained irregular triangular mesh, the triangular mesh is sampled to a regular grid through bilinear interpolation; the grid marked as non-reconstruction area is kept null and the existing DEM is not used for high-resolution filling.

10. A vegetation-shaded slope DEM reconstruction system integrating multi-source remote sensing data, characterized in that... The method for reconstructing a vegetation-shaded slope DEM by fusing multi-source remote sensing data as described in any one of claims 1-9 includes: The data acquisition module is used to acquire UAV lidar point clouds, orthophotos, and existing digital elevation models. The data preprocessing module is used to unify multi-source data into the same geographic coordinate system and perform initial ground filtering to obtain an initial ground point set and a non-ground point set. The ground point reclassification module is used to perform multi-scale terrain continuity analysis on non-ground points, reclassify non-ground points that meet the conditions as ground points, form a corrected ground point set and extract terrain structure lines. The reconstructable area determination module is used to determine the reconstructable area based on the directional distribution of ground points around the missing area and its consistency with the existing DEM topographic trend. The structural constraint reconstruction module is used to divide the terrain into sub-regions using terrain structure lines, establish and solve an objective function that includes boundary elevation constraints, DEM trend constraints and local smoothing constraints, and obtain the reconstructed surface point cloud. The DEM generation module is used to fuse and correct ground point sets with reconstructed surface point clouds to generate a digital elevation model.

Images