Topographic mapping system based on three-dimensional laser scanning

Abstract

The present application belongs to the technical field of topographic geometry measurement, and discloses a topographic surveying and mapping system based on three-dimensional laser scanning, which adopts a scanning and multi-view coaxial acquisition scheme driven by topographic geometry entropy value, constructs a topographic geometry information entropy evaluation function based on a preset survey area and an initial rough scanning point cloud, calculates regional topographic information entropy by fusing slope, curvature and elevation coefficient of variation, quantifies the disorder degree of topographic features and divides topographic grades into three levels; taking topographic elevation measurement uncertainty as a constraint, the optimal laser emission frequency and scanning step angle are matched for different grade regions, low resolution and large step angle sampling is adopted for gentle areas, point cloud density is improved and scanning lines are encrypted for areas with severe changes, characteristic lines such as cliff edge lines and gully axes are directionally sampled with high density, and high complexity areas are preferentially traversed in combination with planned paths; three groups of laser transceiving units with different fields of view are configured, a coaxial reference is established by sharing a rotating main shaft, and phase difference time triggering and coordinate conversion are performed.

Description

Technical Field

[0001] This invention belongs to the field of terrain geometry measurement technology, specifically a terrain mapping system based on three-dimensional laser scanning. Background Technology

[0002] Traditional topographic surveying primarily utilizes single-point positioning equipment such as total stations and GNSS-RTK to collect elevation data of topographic feature points point by point and then interpolate the data to generate contour lines and digital elevation models. This approach suffers from limitations such as low fieldwork efficiency, high labor costs, and significant data gaps in complex terrain areas like steep slopes and canyons. Three-dimensional laser scanning technology, based on the principle of laser pulse ranging, acquires massive amounts of three-dimensional point cloud data of the surveyed area, enabling a complete record of the surface's geometric shape and texture information. However, existing three-dimensional laser scanning topographic surveying systems still have the following technical problems:

[0003] Existing systems mostly use fixed angular resolution or fixed point spacing to perform uniform scanning across the entire area. They cannot adjust the sampling density according to the degree of geometric changes in local terrain. This results in a large amount of redundant point cloud data in flat terrain areas, increasing the computational burden on subsequent data storage and processing. In critical terrain areas with geometric features such as steep slopes, gullies, and abrupt slope changes, insufficient point cloud density leads to the loss of terrain geometric features.

[0004] Due to terrain obstruction and limited scanning field of view, areas with dramatic undulations such as mountains and canyons generally have scanning blind spots and data gaps, and single-view scanning cannot obtain complete surface geometric information. Multi-field acquisition systems mostly use independently installed laser units, lack a unified coaxial geometric measurement benchmark, and the time-sharing acquisition and scanning movement are not synchronized. Multi-station data stitching requires manual placement of reflective targets, resulting in a large workload in the field and easy introduction of human registration errors, which affects the geometric accuracy of the measurement results.

[0005] Existing filtering algorithms mostly distinguish between surface points and non-surface points based on echo intensity and elevation statistical characteristics, resulting in digital elevation models that generally have elevation deviations. Furthermore, they do not introduce differential geometric constraints to address data gaps caused by scanning blind spots, making it difficult to restore the natural geometric patterns of terrain undulations and easily leading to terrain morphology distortion. Summary of the Invention

[0006] The purpose of this invention is to provide a terrain mapping system based on three-dimensional laser scanning to solve one or more of the problems mentioned in the background art.

[0007] To achieve the above objectives, the present invention provides the following technical solution: a terrain mapping system based on three-dimensional laser scanning, comprising an intelligent scanning module, a multi-view acquisition module, a point cloud filtering module, a standard-free registration module, a surface separation module, a cavity filling module, and a terrain reconstruction module;

[0008] Furthermore, the intelligent scanning module constructs a terrain geometric entropy value evaluation function based on the preset survey area range and the initial coarse scan point cloud data. It calculates the terrain geometric information entropy of the corresponding area by comprehensively considering three indicators: local slope, average curvature, and elevation variation coefficient. The terrain geometric entropy value evaluation function quantifies the disorder of the terrain geometric shape and divides the survey area into three levels: gentle curved surface area, general undulation area, and drastic change area.

[0009] The intelligent scanning module uses a preset uncertainty threshold for terrain elevation measurement as a constraint to solve for the optimal laser emission frequency and scanning step angle for different levels of areas. In gently curved areas, it adopts low-resolution scanning with a large step angle to control the sampling density at the lowest level that meets the uncertainty requirements for terrain elevation measurement. In areas with drastic changes, it increases the point cloud sampling density and densifies the scanning lines. It also performs directional densification sampling on terrain geometric feature lines such as steep slope edges and gully axes.

[0010] The intelligent scanning module has a built-in scanning path planning algorithm that prioritizes traversing highly complex areas and sequentially performs coarse scanning assessment, terrain level classification, constraint matching, parameter solving, and fine scanning acquisition processes.

[0011] Furthermore, the multi-view acquisition module is configured with three independent laser transceiver units, corresponding to three fixed field of view directions: frontal view, forward downward view, and side downward view. The three units share the same set of rotating scanning spindle mechanism to establish a coaxial geometric measurement reference with the spindle rotation center as the origin. The multi-view acquisition module pre-calibrates the geometric parameters of the field of view units and obtains geometric measurement parameters such as the spatial angle between the optical axis and the spindle and the optical center offset of each unit, as well as time measurement parameters such as the pulse delay.

[0012] The front-viewing unit performs mid-to-long-range main scanning, the front-downward-viewing unit acquires point cloud data of the near-range steep slope and ridge area, and the side-downward-viewing unit acquires terrain occlusion blind area data on the side of the scanning field of view;

[0013] The multi-view acquisition module adopts phase difference time-division trigger control synchronized with the rotation angle of the rotating scanning spindle mechanism. Three sets of laser pulses are emitted sequentially according to a fixed phase difference. The echo data acquired by each unit is uniformly converted to the coaxial measurement coordinate system according to the pre-calibrated geometric parameters. The single-station scan outputs a multi-view fused point cloud.

[0014] Furthermore, the point cloud filtering module receives the original point cloud data and uses a filtering strategy based on local differential geometric consistency to filter and remove noise points, isolated points, and gross error points.

[0015] The point cloud filtering module fits the local surface of each point and calculates its first and second basic forms based on the spatial neighborhood distribution characteristics of the point cloud. It extracts differential geometric parameters such as normal vector, principal curvature, and Gaussian curvature, calculates the deviation between the differential geometric parameters of a single point and the differential geometric parameters of the fitted neighborhood surface, and identifies discrete noise and jump points.

[0016] The point cloud filtering module combines the echo intensity and waveform broadening attribute parameters of the point cloud to remove invalid echo data formed by specular reflection and water mist scattering. The point cloud filtering module retains terrain geometric feature points and rigid ground feature points during the filtering process.

[0017] Furthermore, the unstandardized registration module receives multi-station surface point cloud data and adopts a registration scheme based on topological geometric invariants. The unstandardized registration module extracts viewpoint-invariant topological geometric feature primitives from the station cloud. These topological geometric feature primitives include point features such as topological saddle points, steep slope turning points, and ridgeline curvature extrema, as well as linear features such as ridgelines, valley lines, and steep slope edges. The unstandardized registration module completes the initial coarse registration through the topological adjacency relationship of the topological geometric feature primitives to obtain the initial transformation matrix between stations.

[0018] The unstandardized registration module uses the curvature-constrained iterative nearest point algorithm for fine registration, selects rigid points on the ground surface for calculation, and removes point cloud data from easily deformable areas such as vegetation and loose soil.

[0019] The standard-free registration module performs overall adjustment optimization of multi-station data, using the topological consistency of terrain geometric feature primitives as a constraint, identifies and eliminates mismatched point pairs, and outputs multi-station cloud stitching results.

[0020] Furthermore, the surface separation module receives waveform laser echo data and extracts geometric feature information such as the multiple echo positions, wavefront rising slope, waveform half width at half maximum and peak spacing of each laser pulse through a waveform geometric parameter decomposition algorithm.

[0021] The surface separation module distinguishes between canopy echoes, branch echoes, and ground echoes based on the number of echoes and vertical distribution characteristics. It constructs a surface geometric continuity discrimination model by combining the differential geometric distribution characteristics of local point clouds. It comprehensively evaluates the consistency of the elevation neighborhood of the point, the stability of the normal vector, the variation law of surface curvature, and the matching degree of the echo geometric parameters, and calculates the probability value of the corresponding point location belonging to the real surface surface.

[0022] The surface separation module extracts a clean set of surface points through threshold segmentation, and separates the vegetation canopy and branch data.

[0023] Furthermore, the cavity filling module receives surface point cloud data and uses an interpolation reconstruction strategy based on terrain differential geometric constraints to fill data cavities formed by scanning blind spots and vegetation occlusion. The cavity filling module performs cavity boundary detection and geometric classification on the surface point cloud, classifying cavities into four types: smooth curved surface cavities, feature line crossing cavities, steep slope boundary cavities, and composite cavities.

[0024] For smooth curved surface cavities, a local terrain surface fitting method with boundary normal vector and curvature constraints is used for interpolation filling; for feature lines crossing cavities, the geometric trend of terrain geometric feature lines such as ridge lines and gully axes is first continued, and then regional surface interpolation filling is performed with terrain geometric feature lines as hard constraints.

[0025] For cavities at the boundary of steep slopes, a piecewise interpolation method constraining the curvature of the upper and lower edges of the steep slopes is used for filling; for complex cavities formed by large-area shading, redundant observation data from overlapping areas of multiple stations are extracted for complementary filling, and geometric consistency verification is performed on the supplemented data before filling.

[0026] Furthermore, the terrain reconstruction module receives a clean and complete surface point cloud, and uses a geometric feature-constrained irregular triangular mesh construction algorithm to generate a digital elevation model (DEM). The terrain reconstruction module adjusts the triangular mesh size according to the point cloud density and terrain curvature, using a larger triangular mesh in flat terrain areas and a denser triangular mesh in areas with abrupt terrain changes.

[0027] The terrain reconstruction module extracts terrain geometric feature lines such as ridge lines, valley lines, and steep slope edges from the surface point cloud, and embeds them as hard constraint edges into the triangular mesh, forcing the edges of the triangular mesh not to cross the terrain geometric feature lines. The terrain reconstruction module has a built-in contour line tracking algorithm to extract multi-level contour lines from the triangular mesh, and performs morphological optimization and line break repair on the contour lines in areas such as steep slopes and gullies based on the differential geometric curvature change law.

[0028] The terrain reconstruction module outputs standard-format DEM raster data, vector contour lines, terrain curvature distribution maps, cross-section data, and other terrain survey results.

[0029] The beneficial effects of this invention are as follows:

[0030] 1. This invention employs a terrain geometric entropy-driven scanning and multi-view coaxial acquisition scheme. Based on a preset survey area and initial coarse-scan point cloud, a terrain geometric information entropy evaluation function is constructed. The regional terrain information entropy is calculated by integrating slope, curvature, and elevation variation coefficients, quantifying the disorder of terrain features and classifying them into three terrain levels. With the uncertainty of terrain elevation measurement as a constraint, the optimal laser emission frequency and scanning step angle are matched for different level areas. Low-resolution, large-step-angle sampling is used in flat areas, while point cloud density is increased and scanning lines are densified in areas with drastic changes. Characteristic lines such as steep slope edges and gully axes are targeted and densified for sampling. High-complexity areas are traversed first in combination with the planned path. Three sets of laser transceiver units with different fields of view are configured, sharing a common rotating spindle to establish a coaxial reference. Through phase differential time-triggered and coordinate transformation, the main scan, close-range steep slopes, and lateral blind areas are covered, reducing redundant sampling, eliminating scanning blind areas, and improving the efficiency and coverage integrity of single-station terrain data acquisition.

[0031] 2. This invention constructs a processing link from raw point clouds to clean surface points. Based on a local differential geometric consistency strategy, it extracts parameters such as normal vectors, principal curvature, and Gaussian curvature through neighborhood surface fitting. It identifies discrete noise and abrupt points by comparing the deviation of single-point parameters with those of the neighborhood. Combined with echo intensity and waveform broadening attributes, it removes invalid echoes such as specular reflection and water mist scattering, while simultaneously preserving topographic feature points and rigid ground features. It extracts view-invariant feature primitives such as topographic saddle points, steep slope inflection points, and ridgelines, and completes coarse registration through topological adjacency matching. It completes fine registration using a curvature-constrained iterative nearest-point algorithm, and combines overall adjustment optimization to remove mismatched point pairs. It decomposes waveform geometric parameters, constructs a surface continuity discrimination model based on differential geometric features, calculates the surface attribution probability, and extracts the surface point set through threshold segmentation. This improves the quality of point cloud data and the accuracy of multi-station registration, and removes vegetation interference.

[0032] 3. This invention employs a categorized cavity filling and feature-constrained terrain reconstruction scheme. First, cavity boundaries are detected and geometrically classified on the surface point cloud, categorizing cavities into four types: gentle curved surfaces, feature line crossings, steep slope boundaries, and composite cavities. For gentle curved surface cavities, local surface fitting interpolation with boundary normal vectors and curvature constraints is used. For feature line crossing cavities, the trend of ridgelines and gully axes is followed by regional interpolation. For steep slope boundary cavities, segmented interpolation with upper and lower edge curvature constraints is used. Composite cavities are filled by fusing redundant point clouds from multiple stations and undergoing geometric consistency verification. A feature-constrained irregular triangular mesh algorithm is used, adjusting the mesh size based on point cloud density and curvature, embedding terrain geometric feature lines as hard constraints, and using contour line tracking and morphological optimization algorithms to generate multi-level contour lines. This repairs data cavities, restores terrain details, and outputs terrain results. Attached Figure Description

[0033] Figure 1 This is a flowchart of the topographic mapping point cloud processing of the present invention;

[0034] Figure 2 This is a flowchart of the point cloud cavity filling process of the present invention. Detailed Implementation

[0035] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0036] The specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. This embodiment takes a small-scale watershed soil and water conservation topographic mapping project in the hilly and mountainous areas of southern China as the application scenario. The survey area covers areas with gentle terraces, undulating slopes, and dense gullies and steep slopes, with vegetation mainly consisting of shrubs and young trees. All parameter values ​​are corresponding results obtained through calibration calculations under this scenario and are not fixed technical limitations of the present invention.

[0037] like Figures 1 to 2 As shown, the terrain mapping system based on three-dimensional laser scanning described in this embodiment consists of an intelligent scanning module, a multi-view acquisition module, a point cloud filtering module, a standard-free registration module, a surface separation module, a cavity filling module, and a terrain reconstruction module connected in sequence, which complete the processing from field scanning to output of indoor results according to the process.

[0038] In this embodiment, the intelligent scanning module first performs a rapid coarse scan of the entire area based on the preset 1.2km² survey area to obtain initial low-resolution point cloud data. Then, it constructs a terrain geometric information entropy evaluation function and calculates the terrain geometric information entropy of the corresponding area by comprehensively considering three indicators: local slope, average curvature, and elevation variation coefficient. This quantifies the degree of disorder in the terrain geometry and ultimately divides the survey area into three levels: gentle curved surface area, general undulation area, and drastic change area.

[0039] Topographic geometric information entropy evaluation function:

[0040] ;

[0041] H represents the terrain geometric entropy value corresponding to a single grid cell; the higher the entropy value, the more severe the terrain undulation and the higher the geometric complexity within the grid.

[0042] i and j represent the serial numbers of the terrain evaluation indicators, which include a total of 3 core terrain indicators;

[0043] x1 represents the normalized local slope value, which is obtained by linearly normalizing the original slope value calculated from the plane fitted by the neighborhood points within the grid.

[0044] x2 represents the normalized average curvature value, which is obtained by linearly normalizing the original average curvature value calculated from the quadratic surface fitting within the grid.

[0045] x3 represents the normalized elevation variation coefficient value, which is obtained by linearly normalizing the original elevation variation coefficient calculated from the grid point elevation sequence.

[0046] This represents the sum of the values ​​of the three normalized topographic indicators.

[0047] The grid cell size is determined based on the initial coarse sweep point cloud density and the minimum topographic feature scale of the survey area. The minimum topographic feature of the survey area is a small gully with a width of 2m. The grid side length is taken as 1m. The number of sampling points in a single grid is not less than 6, which meets the minimum number of points required for quadratic surface fitting in the grid. Statistical calculations are carried out on the point cloud data in each grid.

[0048] For each grid, the local slope value is calculated by fitting a local plane based on the neighborhood point set, the average curvature value is calculated by fitting a quadratic surface, and the elevation variation coefficient is calculated based on the elevation sequence of points within the grid. The three statistics are linearly normalized respectively, and the normalized three indicators are simultaneously input into the terrain geometric entropy value evaluation function to calculate and output the terrain geometric entropy value corresponding to each grid.

[0049] The natural breakpoint method was used to divide the global topographic entropy value sequence into three levels. The calculated thresholds for the two levels were 0.25 and 0.55, respectively. Grids with entropy values ​​less than 0.25 were classified as gently curved areas, those between 0.25 and 0.55 as general undulating areas, and those greater than 0.55 as areas of drastic change.

[0050] The relationship between the entropy value of terrain geometric information and the two-level threshold is compared grid by grid. Based on the comparison results, each grid is assigned to the corresponding terrain level area. The gentle curved surface area is mainly terraced fields, the generally undulating area is mainly gentle slope forest, and the drastic change area is mainly dense gullies and steep slopes.

[0051] The intelligent scanning module uses a terrain elevation measurement uncertainty threshold of ±5cm as a constraint to solve for the optimal laser emission frequency and scanning step angle for different levels of areas. In gently curved areas, it adopts low-resolution scanning with a large step angle to control the sampling point spacing to the minimum level that meets the elevation measurement accuracy requirements. In areas with drastic changes, it reduces the sampling point spacing and densifies the scanning line density according to the elevation measurement uncertainty constraint. At the same time, it performs directional densification sampling on terrain geometric feature lines such as steep slope edges and gully axes.

[0052] Vector trajectory data of steep slope edge lines and gully axis are extracted from the initial coarse point cloud using curvature change detection and edge tracking algorithms. The extracted feature lines are then discretized at nodes. The density step size is set according to 2 / 5 of the baseline sampling interval of the drastic change area along the feature line direction, and sampling control points are set sequentially.

[0053] Calculate the tangent direction of the feature line at each sampling control point, adjust the horizontal and vertical scanning step angles of the scanning device so that the direction of the laser scanning line intersects perpendicularly with the tangent direction of the feature line, simultaneously increase the laser emission frequency of the corresponding scanning line, reduce the spacing between point cloud sampling points along the feature line to within 2 / 5 of the reference sampling spacing in the area of ​​drastic change, and simultaneously record the spatial range of the encrypted sampling area and the corresponding scanning parameters.

[0054] A mapping table of the relationship between terrain grade, sampling point spacing and elevation measurement uncertainty was established in advance through calibration tests. Using the elevation measurement uncertainty threshold of ±5cm as the constraint boundary, the maximum allowable sampling point spacing under the corresponding terrain grade was obtained by looking up the mapping table.

[0055] Combining the scanning distance parameters between the current station and the target area, the maximum sampling point spacing is converted into laser emission frequency and scanning step angle parameters, which are then output to the scanning control unit. When the scanning distance is 100m, based on the elevation transfer error model and accuracy constraints, the scanning step angle corresponding to the smooth curved surface area is 0.12° and the laser emission frequency is 50kHz, the scanning step angle corresponding to the general undulating area is 0.06° and the laser emission frequency is 100kHz, and the scanning step angle corresponding to the drastic change area is 0.03° and the laser emission frequency is 200kHz.

[0056] Identify the boundary vectors of adjacent terrain level regions, and extend them to both sides along the boundary lines to generate transition zones. The width of the transition zone is two to three times the reference sampling interval of the adjacent level regions. Divide the transition zone into sub-intervals at equal intervals along the direction perpendicular to the boundary lines. The number of sub-intervals ensures that the sampling density change rate between adjacent intervals does not exceed 20%. Based on the reference scanning parameters of the two level regions, calculate the laser emission frequency and scanning step angle parameters corresponding to each sub-interval using linear interpolation.

[0057] During the scanning process, the scanning parameters of the corresponding sub-interval are called in real time according to the location of the sub-interval where the scanning field of view is located, so as to achieve a smooth transition of sampling density between adjacent level areas and synchronously record the coordinates of parameter switching nodes within the transition zone.

[0058] The intelligent scanning module has a built-in scanning path planning algorithm that prioritizes traversing highly complex terrain areas such as areas with drastic changes. Each station sequentially performs a complete workflow of coarse scanning assessment, terrain level classification, constraint matching, parameter solving, and fine scanning data acquisition.

[0059] Based on the topographic classification results of the survey area, the spatial distribution range and boundary coordinates of each level area are identified. With the nominal effective scanning radius of the equipment and the minimum overlap requirement of adjacent stations as constraints, candidate station locations are initially deployed in each level area according to the principle of uniform coverage.

[0060] With the optimization objectives of minimizing the total traversal path length and prioritizing coverage of drastically changing areas, a genetic path planning algorithm was used to solve for the optimal traversal order of all candidate stations, generating a scanning path scheme that includes station number, station coordinates, and scanning order. The path scheme was then output to the scanning control unit. A total of 12 stations were deployed, prioritizing the scanning of the upstream drastically changing areas with dense gullies.

[0061] In this embodiment, the multi-view acquisition module is configured with three independent laser transceiver units, which correspond to three fixed field of view directions: frontal view, front downward view, and side downward view. The three units share the same set of rotating scanning spindle mechanism to establish a coaxial geometric measurement reference with the spindle rotation center as the origin.

[0062] Before field operations, the multi-view acquisition module completes the geometric parameter calibration of the field of view units in an open area of ​​the survey area, and obtains geometric measurement parameters such as the spatial angle between the optical axis and the principal axis and the optical center offset of each unit, as well as time measurement parameters such as the pulse delay.

[0063] The standard three-dimensional target was placed sequentially at three preset calibration positions within the scanning field of view: 50m near distance, 150m medium distance, and 300m far distance. Five calibration points with different heights and horizontal deflection angles were set up at each calibration position. For each calibration point, 20 frames of measurement data of the target feature points were collected by three sets of laser transceiver units. The coordinates of the target feature points measured by the total station were taken as the true values.

[0064] A spatial coordinate transformation error equation including the optical axis angle and optical center offset is constructed. The least squares method is used to iteratively solve the three-dimensional offset parameters of the optical axis and principal axis spatial angle and the optical center relative to the principal axis rotation center of each laser transceiver unit. The calibration residual is controlled within 2mm. The timestamps of the laser emission trigger signal and the echo reception signal are synchronously recorded by a synchronous acquisition device. Statistical calculations are performed on multiple sets of timestamp data to obtain the system pulse delay parameters of each unit. All calibration parameters are stored in the system parameter configuration file for subsequent point cloud coordinate transformation.

[0065] The frontal view unit performs a mid-to-long-range main scan within a range of 100-300m, covering the main terrain of the survey area; the front downward view unit collects point cloud data of the near-distance steep slope and ridge area within a range of 50m, filling the near-distance bottom blind zone of the frontal view unit; the side downward view unit collects data of the mountain-occupied blind zone on the side of the scanning field of view, supplementing the lateral data gap of the single-sided scan.

[0066] The multi-view acquisition module adopts phase difference time-division trigger control synchronized with the rotation angle of the rotating scanning spindle mechanism. Three sets of laser pulses are emitted sequentially with a fixed phase difference of 120° spindle rotation angle, and the emission time interval corresponding to the spindle speed of 60 rpm is about 333 ms. The echo data collected by each unit is uniformly converted to the coaxial measurement coordinate system according to the pre-calibrated geometric parameters, and the single-station scan outputs a multi-view fused point cloud.

[0067] For the raw echo ranging data output by each laser transceiver unit, the real-time rotation angle data of the rotating spindle at the current moment is read, and the optical axis spatial angle and optical center offset calibration parameters of the corresponding unit in the system parameter configuration file are called. The three-dimensional coordinate values ​​of the corresponding sampling point in the coaxial measurement coordinate system with the spindle rotation center as the origin are calculated through spatial rigid body coordinate transformation. The coordinates of all sampling points output by the three sets of units are uniformly stored to form single-station multi-view fused point cloud data.

[0068] Identify the overlapping coverage area of ​​the point clouds of three sets of laser transceiver units, and construct a three-dimensional voxel grid in the overlapping area according to the spatial voxel grid size not less than the average point spacing of the overlapping area; traverse all sampling points in the overlapping area and assign each sampling point to the corresponding voxel unit;

[0069] For voxel cells containing multiple sampling points, priority is given to retaining labeled terrain geometric feature points and rigid feature points; if there are no labeled points, the sampling point closest to the voxel center is selected as the retention point, and other duplicate sampling points within the same voxel are removed; no downsampling processing is performed on point cloud data in non-overlapping areas, and the sampling points are directly retained.

[0070] The rotational scanning spindle mechanism's encoder's rotational pulse signal is acquired in real time to obtain the spindle's current real-time rotational angle data. Based on the preset 120° phase difference parameter, three independent trigger pulse signals are generated sequentially with the spindle rotational angle as the reference. The three trigger pulses are output to the corresponding laser transceiver unit's transmission control terminal to control the three sets of laser pulses to be emitted sequentially with a fixed phase difference.

[0071] In this embodiment, the point cloud filtering module receives the raw point cloud data after single-station fusion and adopts a filtering strategy based on local differential geometric consistency to screen and remove noise points, isolated points and gross error points formed by equipment noise and water mist scattering.

[0072] The point cloud filtering module fits the local quadratic surface of each point through the spatial neighborhood distribution characteristics of the point cloud, extracts differential geometric parameters such as normal vector, principal curvature, and Gaussian curvature, calculates the deviation between the differential geometric parameters of a single point and the corresponding parameters of the neighborhood fitted surface, and identifies discrete noise and jump points.

[0073] For each point to be processed, a spherical neighborhood is constructed using four times the average point spacing of the grid point cloud. The search radius is 2m for the smooth curved surface area, 1m for the general undulating area, and 0.4m for the drastic change area. The neighborhood point set of the point is obtained by filtering within the neighborhood range.

[0074] When the number of neighborhood points meets the minimum requirement of being no less than two to three times the number of fitting parameters for the quadratic surface, the coefficients of the local quadratic surface are solved by using the three-dimensional coordinates of the neighborhood point set as input and the least squares fitting method. Based on the fitted local quadratic surface, the normal vector, principal curvature, and Gaussian curvature differential geometric parameters are further extracted.

[0075] A pre-set differential geometric deviation threshold is used. Based on the statistical results of differential geometric parameters in a flat and homogeneous region of the survey area, three times the standard deviation of the parameter sequence is taken as the noise judgment threshold. In this implementation, the normal vector deviation threshold is 8°, and the principal curvature deviation threshold is 0.2m. -1 The Gaussian curvature deviation threshold is 0.05m. -2 For each point to be detected, calculate the angle between its normal vector and the normal vector of the corresponding point on the fitted surface in the neighborhood, as well as the numerical difference between the principal curvature and the Gaussian curvature; compare the difference of each parameter with the preset deviation judgment threshold one by one, and mark the point as a discrete noise point or a jump point when the difference of any parameter exceeds the judgment threshold.

[0076] The point cloud filtering module combines the echo intensity and waveform broadening properties of the point cloud to remove invalid echo data formed by water surface reflection and post-rain fog scattering; during the filtering process, it actively retains topographic geometric feature points such as steep slope edges and field ridges, as well as rigid ground features such as rocks and artificial markers.

[0077] The upper limit threshold, lower limit threshold, and waveform broadening judgment threshold of echo intensity are preset. The effective intensity range is determined according to the standard hard target echo intensity range calibrated by the equipment. The waveform broadening threshold is 1.5 to 2 times the half-width of the standard hard target echo.

[0078] The echo intensity value and waveform broadening parameter value attached to the original point cloud are read point by point; the echo intensity value of a single point is compared with the upper and lower thresholds of the intensity, and points whose echo intensity exceeds the threshold range are marked as suspected invalid echo points; the waveform broadening parameter of a single point is compared with the waveform broadening judgment threshold, and points whose waveform broadening is greater than the preset threshold are marked as suspected invalid echo points; the final set of invalid echo points is confirmed by combining the two judgment results.

[0079] For each of the remaining valid points, the principal curvature and Gaussian curvature parameters are calculated. Points that meet the extreme value conditions are marked as topographic geometric feature points according to the curvature extreme value judgment rule. At the same time, rigid ground features such as rocks and artificial markers are judged by the echo intensity variation coefficient and the compactness of the point cloud circumscribed cube: points of the same material in the neighborhood of a single point are judged to be intensity stable if the echo intensity variation coefficient is less than 0.1, and points with a morphological regularity if the ratio of the volume of the point cloud circumscribed cube to the number of points is less than a threshold. Points that meet both conditions are marked as rigid ground features. Attribute labels are added to the marked topographic geometric feature points and rigid ground features and stored separately for subsequent registration and reconstruction.

[0080] Spatial connectivity analysis is performed on valid points. The connectivity between points is constructed using twice the average point spacing of the point cloud as the neighborhood radius, and the number of points contained in each connectivity is counted. Connectivity with fewer than 1 / 3 of the minimum number of points for quadratic surface fitting is identified as isolated point clusters, and all points in the isolated point clusters are removed from the valid point set. After single-station filtering, a preliminary surface point set is quickly extracted based on the statistical characteristics of the lowest point and elevation of the single echo, and then output to the standardless registration module. After full-domain registration is completed, fine-grained surface separation is performed.

[0081] In this embodiment, the unmarked registration module receives the surface point cloud data after initial separation from a single station and adopts an unmarked registration scheme based on topological geometric invariants. The unmarked registration module extracts viewpoint-invariant topological geometric feature primitives from the station cloud. Point features include topographic saddle points, steep slope turning points, and ridgeline curvature extreme points, while linear features include ridgelines, valley lines, and steep slope edges. Initial coarse registration is completed through the topological adjacency relationship of the topological geometric feature primitives to obtain the initial transformation matrix between stations.

[0082] For each feature primitive, extract local neighborhood point clouds within a radius of 10 to 15 times the average point spacing of the point cloud. Perform scale normalization on the neighborhood point clouds to uniformly map them to a standard scale space. Calculate the local principal direction based on the principal component analysis results of the neighborhood point clouds, and perform orientation normalization on the descriptors of the feature primitives using the principal direction as a reference. Generate a 32-dimensional feature descriptor vector with viewpoint invariance, containing 16 dimensions each of the curvature distribution histogram and the normal vector angle distribution histogram of the neighborhood point clouds. Bind and store the descriptors with the corresponding feature primitives.

[0083] For each site, the topological adjacency descriptor for the terrain geometric feature primitives extracted from the cloud is constructed. The descriptor contains the type, relative spatial position, number of adjacent features, and adjacency direction information of the feature primitive. Descriptor similarity matching is performed between the feature primitive sets of two sites to generate an initial candidate matching pair set.

[0084] A random sampling consensus algorithm is used to screen candidate matching pairs. After removing erroneous matching pairs, valid matching pairs are retained. The initial transformation matrix between the two stations is calculated based on the valid matching pairs. After screening by the random sampling consensus algorithm, the number of valid matching pairs is no less than 30.

[0085] The unstandard registration module uses the curvature constraint iterative nearest point algorithm for fine registration. Before the calculation, it first selects rigid points with stable geometric shape based on local differential geometric features and echo waveform attributes, and eliminates non-rigid points with deformation features. Only rigid points are selected to participate in the calculation.

[0086] In the corresponding point search phase of each iteration of the iterative nearest point algorithm, a constraint condition of curvature difference between point pairs is added. For each set of candidate corresponding points, the curvature difference between the two points is calculated. Only point pairs whose curvature difference is less than twice the standard deviation of curvature of the same surface area are included in the effective corresponding point set of this iteration. Iterative calculation is continuously performed until the difference of transformation obtained by two adjacent iterations is less than 1 / 10 of the single-point measurement accuracy of the system. After reaching the convergence threshold, the iteration stops and the final fine registration transformation matrix is ​​output.

[0087] The standard-free registration module performs overall adjustment optimization of multi-station data across the entire survey area. Using the topological consistency of terrain geometric feature primitives as a constraint, it further identifies and eliminates mismatched point pairs, and finally outputs the overall stitched global point cloud data.

[0088] All spatial transformation parameters of the stations involved in the splicing are incorporated into a unified adjustment calculation system. Feature elements such as topographic saddle points and steep slope inflection points with the same name in the common overlapping area of ​​multiple stations are extracted. The optimization objective is to minimize the sum of squared coordinate deviations of the feature elements with the same name in the common area, and the topological consistency constraint is that the relative deviation of the spatial distance between adjacent feature elements is less than 1%. An indirect adjustment optimization equation with constraints is constructed.

[0089] The optimal spatial transformation parameters for each station were obtained through iterative solution. Based on the optimal parameters, a unified coordinate transformation was performed on the point cloud data of all stations. The results were verified using 20 GNSS-RTK checkpoints in the survey area. After adjustment, the mean square error of the planar registration of the point cloud over the entire area was 1.5 cm, and the mean square error of the elevation registration was 1.8 cm.

[0090] Multi-station global adjustment optimization model with topological constraints:

[0091] ;

[0092] ;

[0093] Let represent the objective function of adjustment optimization, characterize the sum of squared coordinate deviations after transformation of the same feature primitives at multiple stations, and the optimization objective is to minimize its value;

[0094] This represents the set of spatial transformation parameters for all stations, including the rotation matrix and translation vector for each station, which are the parameters to be solved in the adjustment.

[0095] , Let represent the three-dimensional rigid body transformation matrix corresponding to the cloud at station s and station t;

[0096] , This represents the original three-dimensional coordinate vector of the i-th terrain feature primitive with the same name in the s-th and t-th stations;

[0097] This indicates the total number of terrain feature primitives with the same name within the common overlapping area of ​​multiple stations;

[0098] , Let represent the original three-dimensional coordinate vector of the i-th and j-th adjacent terrain feature primitives in the s-th station;

[0099] The reference spatial distance between the i-th and j-th adjacent feature primitives is predetermined by the topological adjacency relationship of the original point cloud;

[0100] This represents the topological consistency constraint threshold, which corresponds to the constraint requirement that the relative deviation of the spatial distance between adjacent feature primitives does not exceed 1%.

[0101] In this embodiment, the surface separation module receives the full-domain point cloud and its associated full-waveform echo data after standard-free registration stitching. For the shrub and young forest covered areas of the survey area, the geometric feature information such as the multiple echo positions, wavefront rising slope, waveform half-width and peak spacing corresponding to each laser pulse is extracted by the waveform geometric parameter decomposition algorithm.

[0102] DC baseline correction processing is performed on the original acquired laser echo waveform sequence to remove the baseline drift component in the waveform signal; the waveform sampling rate is 1GS / s, the inherent rise time of the laser pulse is 3ns, and the sliding window length is 10 sampling points, covering the complete pulse rise time and not less than 3 times the sampling interval;

[0103] In the smoothed waveform sequence, all local peak points are identified, and each peak point is decomposed into an independent sub-waveform corresponding to each peak. For each sub-waveform, the echo position, wavefront rising slope, waveform half-width and peak spacing geometric feature parameters are extracted.

[0104] The surface separation module initially screens the ground echo candidate point set based on the echo number and vertical distribution characteristics. For shrubland with a height of 1-3m in the survey area, it constructs a surface geometric continuity discrimination model by combining the differential geometric distribution characteristics of the local point cloud. It comprehensively evaluates the consistency of the candidate point's elevation neighborhood, the stability of the normal vector, the curvature variation law of the surface and the matching degree of the echo geometric parameters, and calculates the probability value of the corresponding point location belonging to the real surface surface.

[0105] All echoes corresponding to the same laser pulse are sorted from high to low according to their elevation values, and then classified into levels according to the vertical distribution interval and waveform characteristics of the echoes; the echoes with the highest elevation and waveform broadening greater than 6ns are classified as canopy echoes.

[0106] Echoes with a vertical spacing greater than 0.5m between intermediate echo points and a waveform peak amplitude less than 50% of the ground echo peak are classified as branch echoes; echoes with the lowest elevation and a waveform half-width variation coefficient less than 0.2 are classified as ground echoes; and canopy point cloud subsets, branch point cloud subsets, and ground candidate point cloud subsets are generated according to the hierarchical classification results.

[0107] Four indicators—elevation neighborhood consistency, normal vector stability, surface curvature change rate, and echo waveform parameters—were selected as inputs for the surface geometric continuity discrimination model. Each of the four indicators was quantified and normalized. The entropy weight method was used to calculate the weight coefficients based on the information contribution of each indicator in the survey area. The four normalized indicators were then weighted and summed, and the weighted calculation results were mapped to probability values ​​in the range of 0 to 1, which were used as the surface attribution probability of the corresponding points.

[0108] The surface separation module extracts a clean set of surface points through threshold segmentation, and separates the vegetation canopy and branch data.

[0109] The probability values ​​of surface attribution for all points within the survey area are statistically analyzed to generate a probability distribution histogram. The Otsu method is used to calculate the probability distribution histogram to obtain a segmentation threshold. The initial segmentation of surface points and non-surface points is completed using the segmentation threshold as the boundary. Spatial connectivity analysis is performed on the surface point set obtained from the initial segmentation, and the coverage area and number of points in each connectivity are statistically analyzed. Non-surface connectivity areas with an area less than twice the area of ​​a single grid cell are removed. For scattered non-surface point holes within surface connectivity areas, interpolation using neighboring surface points is used for filling and correction.

[0110] In this embodiment, the cavity filling module receives the separated clean surface point cloud data and uses an interpolation reconstruction strategy based on terrain differential geometric constraints to fill the data cavities formed by steep slopes and tree occlusions. The cavity filling module first performs cavity boundary detection and geometric classification on the surface point cloud, classifying the cavities into four types: smooth curved surface cavities, feature line crossing cavities, steep slope boundary cavities, and composite cavities.

[0111] Construct a triangular mesh topology for the input surface point cloud data. The maximum side length is taken as 2 to 3 times the average point spacing of the surface point cloud. Traverse all the edges of the triangular mesh and count the number of times each edge is shared by a triangular face. Mark the edges that are shared by only a single triangular face as free boundary edges and perform head-to-tail connection tracing on all free boundary edges.

[0112] Spatial determination is made by combining the boundary of the survey area with the boundary of the ground features. Free boundary areas corresponding to ground features such as buildings and large rocks are eliminated. Areas enclosed by free boundary edges that are not caused by ground features and are located inside the survey area are determined as point cloud cavities. The boundary point set and boundary ring geometric parameters corresponding to each cavity are extracted, and the spatial location and coverage of the cavity are recorded.

[0113] For each detected cavity, the curvature distribution features of its boundary point set, the boundary closure morphology features, and the correlation features of the surrounding terrain geometric feature lines are extracted as classification criteria. According to the preset classification rules, if the cavity boundary and the terrain geometric feature line have two or more intersections and the connecting line passes through the cavity, it is determined that the feature line passes through the cavity.

[0114] If the length of the cavity boundary coinciding with the edge of the steep slope is more than 30% of the total boundary length, it is judged as a cavity at the edge of the steep slope; if it meets two or more of the above judgment conditions, it is judged as a composite cavity; the rest are judged as cavities on a gentle curved surface, mainly cavities on a gentle curved surface formed by cavities at the edge of the steep slope and cavities on a gentle curved surface formed by trees blocking the view.

[0115] For smooth curved surface cavities, a local terrain surface fitting method with boundary normal vector and curvature constraints is used for interpolation filling; for feature lines crossing cavities, the geometric trend of terrain geometric feature lines such as ridge lines and gully axes is first continued, and then regional surface interpolation filling is performed with terrain geometric feature lines as hard constraints.

[0116] Extract the coordinates of the endpoints of the existing topographic geometric feature lines on both sides of the cavity, and calculate the tangent direction and local curvature value of the feature lines at the two endpoints respectively. Using the coordinates, tangent direction and curvature value of the endpoints on both sides as constraints, use the cubic spline curve interpolation method to generate continuous feature line extensions inside the cavity, so that the extensions and the original feature lines on both sides achieve a second-order smooth connection of tangent and curvature at the endpoints. The generated feature line extensions are used as hard constraints for cavity filling.

[0117] For cavities at the boundary of steep slopes, a piecewise interpolation method constraining the curvature of the upper and lower edges of the steep slopes is used for filling; for complex cavities formed by large-area shading, redundant observation data from overlapping areas of multiple stations are extracted for complementary filling, and geometric consistency verification is performed on the supplemented data before filling.

[0118] Common reference points are uniformly selected in the surrounding boundary area of ​​the cavity to be filled. The number of reference points ensures that the local registration accuracy meets the requirements. The reference points are evenly distributed in different directions of the cavity boundary. Based on the common reference point set, local coordinate registration is performed on the supplementary point cloud from other stations so that the supplementary point cloud and the original point cloud are under the same coordinate reference.

[0119] Point-by-point comparison is performed between the added point cloud and the original landmark point cloud to determine the planar deviation, elevation difference, and curvature difference at corresponding locations. The planar deviation should be less than 2cm, the elevation difference less than 10cm, and the curvature difference less than 0.15m. -1 When filling in the gaps, retain the corresponding data and include it in the gap-filling point set.

[0120] The average point spacing of the original point cloud within a range of 3 times the average point spacing of the point cloud around the hole is statistically analyzed. Using this average point spacing as the benchmark, the filling point set inside the hole is uniformly resampled. The boundary point set between the filling region and the original point cloud is extracted, and a connection region of 3 to 5 times the average point spacing of the point cloud is taken on both sides of the boundary. Neighborhood interpolation smoothing is performed on the original points and filling points in the connection region to adjust the spatial distribution density of the points in the connection region so that the point cloud density on both sides of the boundary is consistent.

[0121] Identify the boundary point sets of the upper and lower edges of the steep slope surrounding the cavity, and perform smooth fitting on the upper and lower edge point sets to obtain the geometric shapes of the top and bottom surfaces of the slope. For the slope wall area, using the boundary points of the upper and lower edges as constraints, a linear interpolation method is used along the normal direction of the slope to generate the transition point cloud of the slope wall. The filling point clouds of the top, wall, and bottom areas of the slope are stitched together to form a complete result of filling the cavity at the boundary of the steep slope.

[0122] In this embodiment, the terrain reconstruction module receives the clean and complete surface point cloud after filling, and uses a geometric feature-constrained irregular triangular mesh construction algorithm to generate a digital elevation model; the terrain reconstruction module adjusts the triangular mesh size according to the point cloud density and terrain curvature, the maximum side length of the triangular mesh is inversely proportional to the terrain curvature, and matches the corresponding side length threshold according to the terrain curvature value.

[0123] Three threshold values ​​for the maximum side length of the triangulation network are pre-set, corresponding to three numerical ranges of terrain curvature from low to high, and corresponding to terrain levels of gentle curved areas, general undulation areas, and drastically changing areas. Based on the relationship model between DEM elevation error and triangulation network side length, with ±5cm elevation error as a constraint, the three threshold values ​​for side length are derived as 5m, 2m, and 0.8m. During the construction of the irregular triangulation network, the terrain curvature value of the current point is calculated point by point, and the corresponding threshold value for side length is matched according to the interval to which the curvature value belongs. This threshold is used to constrain the side length of adjacent nodes during the generation of the triangulation network.

[0124] The terrain reconstruction module extracts topographic geometric feature lines such as ridge lines, valley lines, and steep slope edges from the surface point cloud, and embeds them as hard constraint edges into the triangulation network, forcing the edges of the triangulation network not to cross the topographic geometric feature lines; the terrain reconstruction module has a built-in contour line tracking algorithm, which extracts multi-level contour lines with corresponding contour intervals according to the topographic relief level and mapping scale specifications of the survey area, and optimizes the shape and repairs broken lines of contour lines in areas such as steep slopes and gullies based on the differential geometric curvature change law;

[0125] Based on the 1m contour interval parameter, an elevation value sequence of contour lines at each level is generated. For each level of elevation value, all triangular mesh edges are traversed. The intersection points of the contour line elevation surface and the triangular mesh edges are detected, and the coordinates of the intersection points and the mesh information are recorded. All intersection points are sequentially traced and connected according to the mesh adjacency relationship to form continuous contour line segments. The contour line segments are then processed by closing the breakpoints and simplifying redundant nodes.

[0126] During the construction of irregular triangular meshes, all line segments of terrain geometric feature lines are inserted as hard constraint edges into the triangular mesh generation system; mesh edges that cross hard constraint edges in the initially generated triangular mesh are detected, and vertex splitting and mesh reconstruction are performed on the crossed mesh edges to adjust the vertex connection relationship of the triangular mesh so that the edges of the finally generated triangular mesh coincide with the terrain geometric feature lines or are distributed on both sides of the terrain geometric feature lines.

[0127] Based on the differential geometric parameters of the surface point cloud, the D8 flow direction algorithm is used to calculate the cumulative runoff at each point. Raster lines are extracted from grids with a cumulative runoff greater than 100, and ridge lines are extracted from grids with a cumulative runoff less than 10. Complete vector data is generated through region growing. A surface curvature abrupt change detection algorithm is used to identify steep slope distribution areas, with the curvature abrupt change threshold set at 0.3m. -1 We used edge tracking and line segment fitting methods to extract the vector data of steep slope edges; and performed breakpoint connection and shape smoothing on all extracted terrain geometric feature lines to generate a complete terrain geometric feature line dataset.

[0128] For contour line breaks in steep slopes and gullies, the tangent direction and local curvature change trend at the break point are extracted. The contour line segment is extended to the other side of the break point along the tangent direction and curvature change trend until it connects with the break point on the opposite side. After the break line is connected, the nodes of the contour line connection segment are smoothed by local node interpolation.

[0129] Before outputting the results, 50 GNSS-RTK checkpoints evenly distributed within the survey area were used to verify the elevation accuracy. The overall elevation error of the survey area was 4.2 cm, of which the elevation error in the flat area was 2.8 cm and the elevation error in the area with drastic changes was 4.9 cm. The terrain reconstruction module outputs topographic survey results such as DEM raster data, vector contour lines, terrain curvature distribution maps, and typical cross-section line data in the mapping reference format.

[0130] The grid resolution is set at 1.5 times the average point spacing of the point cloud, and the DEM grid resolution is 0.5m. A regular grid covering the survey area is constructed, and the planar coordinates of each grid point are determined. The generated irregular triangular network is used as the elevation benchmark, and the elevation value corresponding to each grid point is calculated by the triangular face linear interpolation method to generate DEM grid data in standard GeoTIFF format.

[0131] Based on DEM raster data, the average curvature and Gaussian curvature values ​​of the terrain are calculated grid by grid using the 3×3 neighborhood window quadratic surface fitting method to generate a terrain curvature distribution map. According to the preset starting and ending coordinates of the main channel cross-section line of the watershed, the cross-section line trajectory is generated, and the elevation sequence of the corresponding position is extracted from the DEM raster along the trajectory to generate terrain cross-section line data containing plane coordinates and elevation.

[0132] It should be noted that, in this document, relational terms such as "first" and "second" are used only 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 process, method, article, or apparatus.

[0133] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A terrain mapping system based on three-dimensional laser scanning, characterized in that, Includes the following modules: The intelligent scanning module is used to adjust the laser scanning parameters according to the terrain level of the survey area and plan the scanning path to complete the sampling control; The multi-view acquisition module is used to acquire multi-field echo data through multiple sets of coaxially arranged laser transceiver units and output a fused point cloud. The point cloud filtering module is used to remove noise from the original point cloud based on local differential geometric consistency and to remove invalid echoes by combining echo attributes. The unstandardized registration module is used to extract terrain geometric feature primitives to complete coarse registration, and combined with the curvature constraint iterative nearest point algorithm to complete multi-site cloud fine registration; The surface separation module is used to decompose the laser echo waveform and extract the surface point set by combining the surface geometric continuity judgment; The hole filling module is used to detect and classify point cloud holes, fill them using the corresponding methods, and verify the filled data. The terrain reconstruction module is used to generate digital elevation models by constructing irregular triangular networks using geometric feature constraints, and outputs terrain survey results.

2. The terrain mapping system based on three-dimensional laser scanning according to claim 1, characterized in that, The intelligent scanning module constructs a terrain geometric entropy evaluation function based on the preset survey area and the initial coarse scan point cloud data. It calculates the terrain geometric information entropy by comprehensively considering three indicators: local slope, average curvature, and elevation variation coefficient. This quantifies the disorder of the terrain geometric shape and divides the survey area into three levels: gentle curved surface area, general undulation area, and drastic change area. With a preset uncertainty threshold for terrain elevation measurement as a constraint, the optimal laser emission frequency and scanning step angle are solved for different levels of areas, and directional densification sampling is performed on terrain geometric feature lines, including steep slope edges and gully axes.

3. The terrain mapping system based on three-dimensional laser scanning according to claim 2, characterized in that, The intelligent scanning module has a built-in scanning path planning algorithm that prioritizes traversing the drastically changing areas and sequentially executes the complete process of coarse scanning assessment, terrain classification, elevation measurement constraint matching, scanning parameter solving, and fine scanning acquisition.

4. The terrain mapping system based on three-dimensional laser scanning according to claim 3, characterized in that, The multi-view acquisition module is configured with three sets of laser transceiver units with different field of view directions: frontal view, front-downward view, and side-downward view. They share a common rotating scanning spindle mechanism to establish a coaxial geometric measurement reference. Using a phase difference time-division triggering control method synchronized with the spindle rotation angle, three sets of laser pulses are emitted sequentially with a fixed phase difference. Based on the pre-calibrated geometric and time measurement parameters, the multi-field echo data is uniformly converted to the coaxial measurement coordinate system, and a single-station multi-view fused point cloud is output.

5. The terrain mapping system based on three-dimensional laser scanning according to claim 4, characterized in that, The point cloud filtering module fits a local surface by the spatial neighborhood distribution characteristics of the point cloud, extracts differential geometric parameters including normal vector, principal curvature, and Gaussian curvature, calculates the deviation between the differential geometric parameters of a single point and the corresponding parameters of the local surface obtained by fitting the neighborhood point set, and identifies discrete noise and jump points. Simultaneously, by combining echo intensity and waveform broadening properties, invalid echoes formed by specular reflection and water mist scattering are eliminated, while the filtering process retains terrain geometric feature points and rigid ground features.

6. The terrain mapping system based on three-dimensional laser scanning according to claim 5, characterized in that, The unstandardized registration module extracts the view-invariant terrain geometric feature primitives from the multi-site cloud, completes feature primitive matching based on the topological adjacency relationship of the terrain geometric feature primitives, completes the initial coarse registration, and obtains the initial transformation matrix between sites. The curvature-constrained iterative nearest point algorithm is used to select surface rigid points to participate in the calculation to complete the fine registration, and multi-station data overall adjustment optimization is performed.

7. The terrain mapping system based on three-dimensional laser scanning according to claim 6, characterized in that, The surface separation module extracts the multiple echo positions and waveform geometric feature parameters of a single laser pulse through a waveform geometric parameter decomposition algorithm. A surface geometric continuity discrimination model is constructed by combining the differential geometric features of local point clouds, the surface belonging probability of each point is calculated, and the pure surface point set is extracted by threshold segmentation, and the vegetation canopy and branch data are stripped away.

8. The terrain mapping system based on three-dimensional laser scanning according to claim 7, characterized in that, The cavity filling module performs point cloud cavity boundary detection and geometric classification on the surface point cloud, classifying the point cloud cavity into four types: smooth curved surface cavity, feature line crossing cavity, steep slope boundary cavity, and compound cavity. Different types of voids are filled using surface fitting, feature constraint interpolation, piecewise interpolation, and multi-station redundancy fusion methods. After geometric consistency verification, the filled data is retained as valid data.

9. The terrain mapping system based on three-dimensional laser scanning according to claim 8, characterized in that, The terrain reconstruction module uses the geometric feature-constrained irregular triangular mesh construction algorithm to adjust the triangular mesh size according to the point cloud density and the average curvature of the terrain. Larger meshes are used in flat terrain areas, and denser triangular meshes are used in areas with abrupt terrain changes. The terrain geometric feature lines, including ridge lines, valley lines, and steep slope lines, are extracted and embedded as hard constraint edges into the irregular triangular mesh, forcing the mesh edges not to cross the terrain geometric feature lines.

10. The terrain mapping system based on three-dimensional laser scanning according to claim 9, characterized in that, The terrain reconstruction module has a built-in contour line tracking algorithm to extract multi-level contour lines from the irregular triangular network, and optimizes the shape and repairs broken lines of contour lines in steep slopes and gullies based on the differential geometric curvature variation law. The output includes the topographic survey results, including standard format digital elevation model raster data, vector contour lines, and topographic curvature distribution maps.

Images