A laser radar point cloud DTM reconstruction method for geological disaster surveying

Abstract

The application discloses a kind of laser radar point cloud DTM reconstruction methods for geological disaster survey, to solve the problem of microrelief distortion and feature recognition accuracy caused by vegetation obstruction.The method comprises the following steps: data acquisition, obtaining laser radar point cloud data and auxiliary data such as orthophoto, slope vegetation data; data preprocessing, original data solving, point cloud denoising, retaining intensity and normal vector information, and exporting.las format; point cloud elevation flip constructs "inverted mold", extracts ground feature points through mirror transformation and rasterization, and generates flip DSM; DTM reconstruction and optimization, flip DSM restores true geographic reference, and high-precision DTM is obtained after optimization.The present application can effectively solve the problem of important microrelief distortion and feature recognition accuracy caused by vegetation obstruction when establishing terrain model by using the present point cloud classification technology.

Description

Technical Field

[0001] This invention relates to the field of railway exploration technology, specifically to a lidar point cloud (DTM) reconstruction method for geological disaster surveying. Background Technology

[0002] In geological hazard surveying in mountainous areas, the accuracy and reliability of the survey are often constrained by surface vegetation and complex terrain. Traditional survey methods are difficult to penetrate dense forest canopies and cannot effectively obtain true surface information of vegetated areas, resulting in the omission of key micro-topographic features such as unstable rocks and fissures, which brings great uncertainty to geological hazard exploration.

[0003] LiDAR technology provides a novel approach to geological hazard surveying in complex environments. Its active remote sensing capabilities allow it to penetrate dense vegetation canopies to acquire high-precision ground point cloud data, providing centimeter-level or even millimeter-level surface morphology information. This significantly improves the ability to detect potential geological hazards in vegetated areas, overcoming the limitations of traditional optical remote sensing in identifying concealed geological hazards. Currently, the general process for generating digital elevation models using airborne LiDAR includes data acquisition, noise reduction, point cloud classification, and terrain modeling.

[0004] However, existing point cloud classification methods have significant shortcomings when dealing with steep slopes, protruding rocks, and abrupt changes in micro-topography. Traditional filtering algorithms often misclassify these point clouds, which represent geological hazards, as non-ground points and discard them, resulting in overly smoothed DTM terrain and severe distortion of important micro-topographic features. This problem seriously weakens the application effect of lidar technology in the detailed surveying of geological hazards such as small unstable rock masses, shallow cracks, and creeping slopes.

[0005] Therefore, in response to the urgent need for refined geological hazard surveys, a refined DTM reconstruction method that can be applied to densely vegetated areas and preserve micro-topographic features to the maximum extent is designed. This method has significant engineering implications and promotional value for improving the accuracy of geological hazard surveys and achieving early identification of geological hazards. Summary of the Invention

[0006] This invention aims to at least solve one of the technical problems existing in the prior art. To this end, one objective of this invention is to propose a lidar point cloud DTM reconstruction method for geological disaster surveying, in order to solve the problems of distortion of important micro-topography and insufficient feature recognition accuracy caused by vegetation occlusion when building terrain models in existing point cloud classification technologies.

[0007] To address the above problems, this invention provides a lidar point cloud (DTM) reconstruction method for geological disaster surveying, comprising the following steps:

[0008] Step S1, data acquisition: acquire lidar point cloud data and auxiliary data for the target area; the auxiliary data includes orthophotos, slope vegetation data and geological data;

[0009] Step S2, data preprocessing, includes solving the raw lidar data, denoising the point cloud, preserving the intensity and normal vector information, and exporting in the .las standard format;

[0010] Step S3: Flip the point cloud elevation and construct a "reverse model". Perform a mirror symmetric transformation on the preprocessed point cloud elevation value "z" relative to the xoy plane to obtain "-z". Rasterize the flipped point cloud and extract the highest elevation points with rock and soil characteristics and reasonable normal vector distribution in each grid as surface feature points by fusing intensity and normal vector information. Based on the extracted feature point set, use interpolation methods to generate the flipped digital surface model (DSM).

[0011] Step S4, DTM reconstruction and optimization: The digital surface model obtained by "molding" is mirrored relative to the xoy plane and restored to the real-world geographic reference system to obtain the initial DTM; the initial DTM is optimized and smoothed to remove abnormal surface protrusions, and the terrain features are enhanced based on slope, aspect and curvature information to obtain the final high-precision DTM for geological disaster survey.

[0012] Preferably, in step S1, lidar data is acquired by a drone platform carrying a lidar sensor flying along a designed route; orthophotos are acquired synchronously with lidar data; vegetation data is acquired by collecting vegetation data and conducting on-site surveys; and geological data is acquired by collecting geological information and conducting on-site investigations.

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

[0014] Step S201, Data processing: By fusing the raw observation data collected by the lidar, the absolute position data provided by the GNSS receiver, and the attitude data provided by the IMU, the point cloud is accurately mapped from the sensor coordinate system to the geodetic coordinate system, thereby generating a structured point cloud result with absolute geospatial position for terrain modeling and geological hazard feature extraction. At the same time, the normal vector of each point is calculated based on the principal component analysis (PCA) method.

[0015] Step S202, Point Cloud Denoising: A statistical outlier removal algorithm is used, and a dynamic threshold is set according to the point cloud intensity information to remove noise points that are significantly off the ground or outliers that are completely separated from the main point cloud due to equipment errors or environmental interference.

[0016] Preferably, in step S3, when flipping the point cloud elevation: first, a point cloud expansion matrix P is established based on the preprocessed point cloud.

[0017] ;

[0018] The extended matrix includes point coordinates (x, y, z) + intensity (I) + normal vector (n). x n y n z ); (x N ,y N ,z N Let I be the spatial coordinates of the Nth point. N Let n be the reflection intensity value at the Nth point. xN ,n yN ,n zN () is the normal vector of the Nth point.

[0019] Then, the elevation value z of each point is extracted from the extended matrix, and a transformation of "z = -z" is performed on each z value to obtain the transformed matrix P'. The transformed point cloud normal vector is then transformed into (n x n y -n z );

[0020] ;

[0021] Then, surface feature point extraction is performed in units of graticules. The graticule size is set to 0.1m×0.1m~0.5m×0.5m. The graticule size is adjusted according to the point cloud density and terrain complexity. The highest elevation point with rock and soil reflection characteristics and reasonable normal vector distribution in each graticule is extracted to form a feature point set S. The rock and soil characteristics include (1) reflection intensity characteristics: the intensity value of the point is within the preset rock and soil intensity threshold range, preferably 80~150; (2) geometric characteristics: the angle between the normal vector of the point and the normal vector of the neighboring points is less than the set threshold, generally <30°, which conforms to the continuous or gradually changing characteristics of the surface. The feature point set S is represented as follows.

[0022] ;

[0023] The continuity of the normal vector is constrained by calculating the angle between the normal vectors of a point and its neighboring points:

[0024] ;

[0025] Where n j With n k These are the unit normal vectors of the current point and its neighboring points, respectively. ， The threshold for the angle between the normal vectors is taken from 10° to 45° based on the terrain features.

[0026] Finally, when constructing the "reverse mold", the extracted feature point set S is used to generate a DSM that is mirror-flipped relative to the xoy plane by using Kriging interpolation or progressive triangulation encryption method, so as to obtain the "reverse mold".

[0027] Preferably, in step S4, the initial DTM optimization includes needle removal and feature enhancement, comprising the following steps:

[0028] Step S401, Removal of needle-like structures: A 3×3 window Gaussian filter is used to smooth the DTM. The filter parameters are dynamically adjusted based on the intensity information to eliminate small protrusions and noise caused by vegetation. The Gaussian kernel function is:

[0029] ;

[0030] in, Let be the coordinate offset of any point in the core relative to the center point (0, 0). The standard deviation controls the "width" of the Gaussian distribution, thus determining the strength of the smoothing.

[0031] The needle-like protrusions caused by the vegetation residue of weeds were removed using morphological operations;

[0032] Step S402, terrain feature enhancement: To highlight the features of dangerous rocks or steep slopes, the terrain feature lines are enhanced based on slope, aspect and curvature. High curvature areas correspond to the locations of micro-topographical abrupt changes.

[0033] The advantages of this invention compared to the prior art are:

[0034] The present invention provides a lidar point cloud DTM reconstruction method for geological disaster survey. By deeply integrating point cloud intensity and normal vector information, it significantly improves the accuracy of terrain reconstruction in vegetated areas and can effectively identify micro-topographic features that are difficult to detect by traditional methods.

[0035] In addition, the elevation mirror flipping and feature enhancement strategy proposed in this invention can preserve the geological disaster features such as steep rock walls and protruding rocks to the maximum extent, and solve the problem of excessive smoothing in traditional filtering algorithms; it greatly improves the efficiency and accuracy of geological disaster investigation and provides reliable technical support for geological disaster prevention and control. Attached Figure Description

[0036] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0037] Figure 1 This is a flowchart of the fine three-dimensional terrain reconstruction method based on lidar in this invention;

[0038] Figure 2 This is a schematic diagram of the UAV-based data acquisition route design in this invention;

[0039] Figure 3 This is a point cloud feature map of an isolated rock on a slope in this invention;

[0040] Figure 4 The image shows a comparison between the DTM constructed by the traditional point cloud classification method and the "reverse molding" method in this invention, where (a) is the traditional point cloud classification method and (b) is the DTM constructed by the "reverse molding" method in this invention. Detailed Implementation

[0041] The embodiments of this application are described in detail below. Examples of the embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain this application, and should not be construed as limiting this application.

[0042] In the description of this application, it should be noted that, unless otherwise expressly specified and limited, the terms "installation," "connection," and "joining" should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral connection; they can refer to a mechanical connection or an electrical connection; they can refer to a direct connection or an indirect connection through an intermediate medium; they can refer to the internal communication of two components or the interaction between two components. Those skilled in the art can understand the specific meaning of the above terms in this application according to the specific circumstances.

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

[0044] like Figure 1 As shown in the figure, this is a technical flowchart of a high-precision DTM reconstruction method for lidar point clouds used in geological disaster surveying, according to this embodiment of the present application, including the following steps:

[0045] Step S1, Data Acquisition: Conduct a rockfall hazard survey on a railway slope in a mountainous area in southern China. Use a DJI M300RTK drone equipped with a Riegl VUX-160 lidar system to collect lidar data on the target area, obtain high-density point cloud data, simultaneously collect orthophotos of the target area, and collect geological conditions, vegetation and other data as supplementary information.

[0046] Specifically, in this embodiment, a drone flight platform equipped with a lidar sensor is used to collect laser data of the target area, thereby obtaining lidar point cloud data of the target area and simultaneously obtaining orthophoto data.

[0047] To ensure point cloud acquisition density and data quality, such as Figure 2 As shown, data was collected using an orthogonal flight path design and a terrain-following flight mode. The flight altitude was 100 meters, the overlap rate was 70%, and the obtained point cloud density was 200 points per m. 2 ~400 points per m 2 .

[0048] Specifically, geological data of the target area is collected to obtain geological condition data, and vegetation cover data is obtained through collection and field investigation;

[0049] Step S2, data preprocessing, including data calculation, point cloud denoising, and exporting the data to the .las standard format;

[0050] Specifically, in step S201, data processing: by fusing the raw observation data collected by the lidar, the absolute position data provided by the GNSS receiver, and the attitude data provided by the IMU, the registration of multi-site cloud data is achieved through the ICP algorithm to ensure the consistency of spatial information, and the point cloud is accurately mapped from the sensor coordinate system to the geodetic coordinate system, thereby generating a structured point cloud result with absolute geographic spatial position that can be used for terrain modeling and geological disaster feature extraction;

[0051] Specifically, in step S202, point cloud denoising: a statistical discrete group removal method is used, setting the neighborhood point data to 25, the standard deviation multiplier to 1.0, and the intensity threshold to 80~150, to remove discrete noise points and outliers that are significantly deviated from the ground. The processed point cloud is exported in a .las standardized format. Figure 3 As shown, the preprocessed point cloud data can reflect information such as vegetation and landform. By extracting profiles from the point cloud, abrupt changes in local topography under vegetation can be clearly observed.

[0052] Step S3: Mirror and flip the point cloud elevation to construct a "reverse mold";

[0053] Specifically, in step S301, point cloud elevation mirroring and flipping: First, a point cloud extension matrix P is established based on the preprocessed point cloud.

[0054] ;

[0055] The extended matrix includes point coordinates (x, y, z) + intensity (I) + normal vector (n). x ,n y ,n z ); (x N,y N ,z N Let I be the spatial coordinates of the Nth point. N Let n be the reflection intensity value at the Nth point. xN ,n yN ,n zN () is the normal vector of the Nth point.

[0056] Then, the elevation value z of each point is extracted from the extended matrix, and a transformation of "z = -z" is performed on each z value to obtain the transformed matrix P'. The transformed point cloud normal vector is then transformed into (n x ,n y ,-n z );

[0057] ;

[0058] Specifically, in step S302, "reverse modeling" is constructed: surface feature point extraction is performed in units of grids. The grid size is generally set to 0.2m×0.2m. The grid size can be adjusted according to the point cloud density and terrain complexity. The highest elevation point with soil and rock mass characteristics in each grid is extracted to form a feature point set S. The soil and rock mass characteristics include (1) reflection intensity characteristics: the intensity value of the point is within the preset soil and rock mass intensity threshold range, preferably 80~150; (2) geometric characteristics: the angle between the normal vector of the point and the normal vector of the neighboring points is less than the set threshold, generally <30°, which conforms to the continuous or gradually changing characteristics of the surface. The feature point set S is represented as follows:

[0059]

[0060] Specifically, in step S303, the "molding" construction step is to use the extracted feature point set S, which can be encrypted using a progressive triangulation method, to generate a DSM that is mirror-flipped relative to the xoy plane.

[0061] Step S4, DTM reconstruction and optimization: The digital surface model obtained by "molding" is mirrored relative to the xoy plane and restored to the real-world geographic reference system to obtain the initial DTM; the initial DTM is optimized and smoothed to remove abnormal surface protrusions, and the terrain features are enhanced based on slope and aspect information to obtain the final DTM used for geological hazard survey.

[0062] Specifically, in step S401, Python code " The data is rewritten into a new file to achieve the overall flipping of the digital surface model, restore it to the real-world geographic reference system, and generate the initial DTM under the real reference system.

[0063] Specifically, in step S402, needle-like structures are removed: a 3×3 window Gaussian filter is used to smooth the DSM, eliminating small protrusions and noise. The Gaussian kernel function is:

[0064] ;

[0065] in Using a value of 0.8, morphological operations were employed to remove needle-like protrusions caused by residual vegetation such as weeds.

[0066] Specifically, in step S403, terrain feature enhancement: based on slope and aspect, the terrain feature lines are enhanced to ensure the reliability and realism of the terrain morphology and obtain a high-precision DTM for geological disaster survey.

[0067] like Figure 4 As shown in (a), the optimized DTM can clearly display isolated rock masses protruding from the slope surface, and its depiction of terrain details is significantly better than that of terrain surface models generated by conventional point cloud classification methods, such as... Figure 4 (b) Using the final high-precision DTM, six unstable rock masses were successfully identified on the slope, with a maximum volume of approximately 15 m³. 3 .

[0068] The method described in this invention effectively solves the problem of significant defects in traditional point cloud segmentation methods for generating 3D terrain models, and can meet the requirements of refined geomorphological analysis of geological hazards. The method is applicable to large-scale point cloud data processing, especially for areas with complex terrain and dense vegetation cover. The high-precision 3D surface model established by this invention can restore ground morphological features to the greatest extent possible and can be applied to the identification and analysis of hidden geological hazards in mountainous areas.

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

Claims

1. A lidar point cloud DTM reconstruction method for geological disaster surveying, characterized in that, Includes the following steps: Step S1, data acquisition: acquire lidar point cloud data and auxiliary data for the target area; the auxiliary data includes orthophotos, slope vegetation data and geological data; Step S2, data preprocessing, includes solving the raw lidar data, denoising the point cloud, preserving the intensity and normal vector information, and exporting in the .las standard format; Step S3: Flip the point cloud elevation and construct a "reverse model". Perform a mirror symmetric transformation on the preprocessed point cloud elevation value "z" relative to the xoy plane to obtain "-z". Rasterize the flipped point cloud and extract the highest elevation points with rock and soil characteristics and reasonable normal vector distribution in each grid as surface feature points by fusing intensity and normal vector information. Based on the extracted feature point set, use interpolation methods to generate the flipped digital surface model (DSM). Step S4, DTM reconstruction and optimization: The digital surface model obtained by "molding" is mirrored relative to the xoy plane and restored to the real-world geographic reference system to obtain the initial DTM; The initial DTM is optimized and smoothed to remove abnormal surface protrusions, and the terrain features are enhanced based on slope, aspect and curvature information to obtain a high-precision DTM for geological hazard survey.

2. The lidar point cloud DTM reconstruction method for geological disaster surveying according to claim 1, characterized in that: In step S1, lidar data is acquired by a drone platform carrying a lidar sensor flying along a designed route; orthophotos are acquired synchronously with lidar data; vegetation data is acquired by collecting vegetation data and conducting on-site surveys; and geological data is acquired by collecting geological information and conducting on-site investigations.

3. The lidar point cloud DTM reconstruction method for geological disaster surveying according to claim 1, characterized in that: In step S2, data preprocessing includes the following steps: Step S201, Data processing: By fusing the raw observation data collected by the lidar, the absolute position data provided by the GNSS receiver, and the attitude data provided by the IMU, the point cloud is accurately mapped from the sensor coordinate system to the geodetic coordinate system, thereby generating a structured point cloud result with absolute geospatial position for terrain modeling and geological hazard feature extraction. At the same time, the normal vector of each point is calculated based on the principal component analysis (PCA) method. Step S202, Point Cloud Denoising: A statistical outlier removal algorithm is used, and a dynamic threshold is set according to the point cloud intensity information to remove noise points that are significantly off the ground or outliers that are completely separated from the main point cloud due to equipment errors or environmental interference.

4. The lidar point cloud DTM reconstruction method for geological disaster surveying according to claim 1, characterized in that: In step S3, when flipping the point cloud elevation: firstly, a point cloud expansion matrix P is established based on the preprocessed point cloud. ； The extended matrix includes point coordinates (x, y, z) + intensity (I) + normal vector (n). x n y n z ); (x N ,y N ,z N Let I be the spatial coordinates of the Nth point. N Let n be the reflection intensity value at the Nth point. xN ,n yN ,n zN () is the normal vector of the Nth point; Then, the elevation value z of each point is extracted from the extended matrix, and a transformation "z=-z" is performed on each z value to obtain the transformed matrix P'. The transformed point cloud normal vector is then transformed into (n x n y -n z ); P’= ； Subsequently, surface feature point extraction is performed on a grid-by-grid basis, with grid sizes ranging from 0.1m×0.1m to 0.5m×0.5m. The grid size is adjusted based on the point cloud density and terrain complexity. The highest elevation points within each grid, possessing both soil and rock mass characteristics and a reasonable normal vector distribution, are extracted to form a feature point set S. These soil and rock mass characteristics include: Reflection intensity characteristics: The intensity value of the point is within the preset rock and soil strength threshold range; Geometric features: The angle between the normal vector of a point and the normal vectors of its neighboring points is less than a set threshold, conforming to the characteristics of continuous or gradually changing land surface. The feature point set S is defined as follows. ； The continuity of the normal vector is constrained by calculating the angle between the normal vectors of a point and its neighboring points: ； Where n j With n k These are the unit normal vectors of the current point and its neighboring points, respectively. ， The threshold for the angle between the normal vectors is taken from 10° to 45° based on the terrain features. Finally, when constructing the "reverse mold", the extracted feature point set S is used to generate a DSM that is mirror-flipped relative to the xoy plane by using Kriging interpolation or progressive triangulation encryption method, so as to obtain the "reverse mold".

5. The lidar point cloud DTM reconstruction method for geological disaster surveying according to claim 1, characterized in that: In step S4, the initial DTM optimization includes needle removal and feature enhancement, comprising the following steps: Step S401, Removal of needle-like structures: A 3×3 window Gaussian filter is used to smooth the DTM. The filter parameters are dynamically adjusted based on the intensity information to eliminate small protrusions and noise caused by vegetation. The Gaussian kernel function is: ； in, Let be the coordinate offset of any point in the core relative to the center point (0, 0). The standard deviation controls the "width" of the Gaussian distribution, thus determining the strength of the smoothing. The needle-like protrusions caused by the vegetation residue of weeds were removed using morphological operations; Step S402, terrain feature enhancement: To highlight the features of dangerous rocks or steep slopes, the terrain feature lines are enhanced based on slope, aspect and curvature. High curvature areas correspond to the locations of micro-topographical abrupt changes.

Images