Map generation method based on laser radar point cloud density analysis

Abstract

The invention discloses a map generation method based on laser radar point cloud density analysis. The map generation method comprises a modeling stage and a mapping stage. In the modeling stage, a map model and a radar density model are established, the map model comprises geometric information and spatial density information of a terrain, and the radar density model is used for representing a spatial density normalization coefficient of a point cloud. In the mapping stage, firstly, laser radar point cloud data are preprocessed; secondly, calculating voxel point cloud occupancy density and obtaining spatial density representation in combination with a density normalization coefficient, calculating a height difference between adjacent voxel columns to obtain a step height feature, and performing plane fitting on each region by adopting a local terrain region division strategy to extract gradient and roughness attributes; and finally, integrating the elements to obtain the elevation-like map. According to the method, multi-dimensional description of spatial density information and topographic features in a non-structural environment can be effectively expressed, and the path planning efficiency and safety of a ground intelligent mobile platform are improved.

Description

Technical Field

[0001] The present invention belongs to the field of ground intelligent mobile platform positioning and navigation, and in particular relates to a map generation method based on laser radar point cloud density analysis. Background Art

[0002] Real-time mapping is one of the key technologies in the ground intelligent mobile platform system. The complete representation of spatial information helps to improve the navigation performance and efficiency of the intelligent mobile platform. As a major environmental perception sensor, LiDAR plays an important role in the intelligent mobile platform as a sensor for real-time positioning and mapping methods. In unstructured environments, the terrain is complex and changeable, with a large number of obstacles and inaccessible areas; maps for structured scenes cannot fully express the spatial information of unstructured environments; therefore, a new solution is needed, an efficient and accurate map creation method that can improve the performance of autonomous navigation and exploration of ground intelligent mobile platforms in unstructured environments.

[0003] Chinese patent application: A method for constructing a navigation cost map for an unmanned vehicle in a complex unknown environment (application number: CN202210930007.8) discloses a method for constructing a navigation cost map for an unmanned vehicle in a complex unknown environment. However, this method does not analyze the density information of spatial objects, and the environmental representation is incomplete.

[0004] Chinese patent application: Method for constructing a local semantic traversable probability grid map for an autonomous mobile robot (application number: CN202211523195.9) discloses a method for generating a local semantic traversable probability grid map based on semantic segmentation and lidar. However, this method only constructs a traversable probability map based on semantic segmentation and lidar, and lacks spatial density analysis. Summary of the invention

[0005] The present invention aims to solve the above problems of the prior art. A map generation method based on laser radar point cloud density analysis is proposed. The technical solution of the present invention is as follows:

[0006] A map generation method based on LiDAR point cloud density analysis includes a modeling stage and a mapping stage, wherein:

[0007] The modeling stage is used to establish a map model and a radar density model; the map model contains the geometric information and spatial density information of the terrain, and is composed of a three-dimensional voxel grid. Each voxel grid contains attributes such as coordinate information, voxel density, step height, terrain slope and roughness; the radar density model is used to characterize the spatial density normalization coefficient of the point cloud, and the emission direction of each laser beam is calculated by setting the horizontal resolution and vertical resolution of the laser radar; voxelization is performed at the radar coordinate origin, and the beam flux of each voxel grid is calculated by ray tracing, and the inverse of the beam flux is used as the voxel density normalization coefficient of the voxel grid;

[0008] In the mapping stage, the laser radar is used to collect environmental data in real time, and the collected point cloud is preprocessed. The point cloud of the area of ​​interest is obtained after straight-through filtering and radius filtering, and voxelized, and the point cloud data is mapped to the corresponding voxel grid; then the map attributes of each voxel grid are calculated to obtain attributes including voxel density, terrain slope and roughness, and the above information is combined to construct a quasi-elevation map.

[0009] Furthermore, the map model is:

[0010] The map G is stored in a three-dimensional voxel grid, G = {T(i,j,k)|i∈[-N x ,N x ],j∈[-N y ,N y ],z∈[-H z ,H z ]}, T(i,j,k)=[i,j,k,roughness,step_height,density,slope]; where N x , N y , H z They represent the boundary range of the x, y, and z axes of the map respectively, T represents the voxel grid, i, j, and k represent the coordinate positions of the three-dimensional voxel grid, roughness represents the roughness, step_height represents the step height, density represents the normalized density, and slope represents the terrain slope.

[0011] Furthermore, the establishment of the radar density model specifically includes the following steps:

[0012] The radar density model is stored in a three-dimensional voxel grid.

[0013] G={T(i,j,k)|i∈[-N x ,N x ],j∈[-N y ,N y ],k∈[-H z ,Hz ]}, T(i,j,k)=[i,j,k,ρ];

[0014] Where ρ is the voxel density normalization coefficient;

[0015] 3.1 First, set the beam distribution mode, and the vertical resolution is a set of different angles, which represent the vertical angle between each beam and the center line of the LiDAR: vertical distribution angle is a set of preset values ​​{θ1,θ2,...,θ i} represents the vertical angle between each laser beam and the center line, where i represents the number of laser beams; and the horizontal distribution angle θ k , then through θ k = k*Δθ calculation, where k is the index of the horizontal angle, and its value range is Δθ represents the angular interval between adjacent scanning points in the horizontal direction, and then the beam direction is calculated

[0016]

[0017] 3.2 The origin of the laser radar coordinate system is at a point outside the space (o x ,o y ,o z ), and the direction of the light is given by the vector d = (d x ,d y ,d z ) is given; any point on the ray (x t ,y t ,z t ) can be expressed as

[0018]

[0019] Where t≥0 is the side length of the voxel grid;

[0020] For each laser beam, starting from the laser radar position, the voxels passed along the direction of the light are calculated in sequence, and the number of times each voxel is passed through is recorded; the voxel index (i t ,j t ,k t ) is calculated as

[0021]

[0022] The voxel grid beam flux is T(i,j,k):

[0023]

[0024] Among them, δ is an indicator function. When the condition is true, the δ value is 1, otherwise it is 0; is the density normalization coefficient. ray represents each laser beam, which is determined by the previous d k,i The beam direction is determined.

[0025] Furthermore, the construction of the elevation map comprises the following steps:

[0026] 4.1 Use the LiDAR device to collect point cloud data in the environment, and pre-process the data using straight-through filtering and radius filtering to retain the point cloud in the area of ​​interest;

[0027] 4.2 Perform voxel processing, map the point cloud of the area of ​​interest to the corresponding voxel grid, and calculate the number of point clouds in each voxel grid, multiply it by the density normalization coefficient of the corresponding voxel grid in the radar density model to obtain the normalized density of the voxel;

[0028] 4.3 Taking the origin of the laser radar coordinate system, differentiate each column of voxel columns on the xoy plane to obtain the step height between adjacent voxel columns;

[0029] 4.4 Classify each voxel into ground points and non-ground points, and then divide the entire map into multiple L×L local areas; perform plane fitting on each local area, the normal vector of the plane fitting represents the slope of the area, and the residual after fitting represents the roughness of the terrain;

[0030] 4.5 The calculated voxel normalized density, step height, terrain slope and roughness are combined to generate a quasi-elevation map.

[0031] Furthermore, the step 4.2 of voxel grid normalized density calculation includes the following steps:

[0032] 5.1 Calculation of actual voxel occupancy density

[0033] The process of voxelizing point clouds is to map the collected point cloud data into a discrete voxel grid; point cloud (x i ,y i ,z i ) and voxel size v, the voxel index (i,j,k) is calculated as:

[0034]

[0035] Define a voxel set V, where each voxel contains O(i,j,k) points.

[0036]

[0037] Where P represents the point cloud dataset, δ is the indicator function, and when the condition is true, the δ value is 1, otherwise it is 0;

[0038] 5.2 Density Normalization

[0039] Normalized density D(i,j,k), which is the actual occupancy density O(i,j,k) of voxel (i,j,k) and the density normalization coefficient The product of

[0040]

[0041] Furthermore, the step 4.4 classifies each voxel into ground points and non-ground points, and then divides the entire map into multiple L×L local areas; plane fitting is performed on each local area, and the normal vector of the plane fitting represents the slope of the area, and the residual after fitting represents the roughness of the terrain, which specifically includes the following steps:

[0042] 6.1 Ground point selection

[0043] In the map G, for the voxel column P(i,j)={T(i,j,k)|k∈[-H z ,H z ]}, when the voxel grid density is below a threshold: when density(T(i,j,k))<ρ low , T(i,j,k) is a non-ground point; when the voxel height exceeds the threshold: Then P(i,j) is marked as a non-ground point, and the remaining voxels are marked as ground points; where ρ low is the density threshold, Hth is the height threshold, density(T) represents the normalized density of voxel T, z(T) represents the height value of voxel T; H z Indicates the height area of ​​the map.

[0044] 6.2 Local terrain meshing

[0045] The basic unit of the terrain grid is B p,q :B p,q ={T(i,j,k)∈G|i∈[pL,(p+1)L],j∈[qL,(q+1)L]}, where Represents the grid index in the x direction, Represents the grid index in the y direction;

[0046] 6.3 Ground Point Downsampling

[0047] Downsample the local terrain grid: B′ p,q =S(B p,q ,Δs), where Δs is the downsampling interval, in B p,q Select the set of downsampling points from: B′ p,q ={T(i,j,k)∈B p,q |i=pL+mΔs,j=qL+nΔs}, where is the row sampling index, is the column sampling index; B′ p,q Represents the local terrain after sampling.

[0048] 6.4 Plane Fitting

[0049] For each terrain grid B′ p,q The ground point set in Let the plane equation be:

[0050] ax+by+cz+d=0#(7)

[0051] Use the least square method to perform plane fitting and obtain the plane coefficients a, b and d;

[0052] 6.5 Residual Calculation

[0053] The above calculations yield the plane coefficients a, b, and d; then the residual error R′ of the fit p,q That is

[0054]

[0055] Indicates the roughness of the ground.

[0056] An electronic device comprises a memory, a processor and a computer program stored in the memory and executable on the processor, wherein the processor implements the map generation method based on laser radar point cloud density analysis as claimed in any one of claims 1 to 6 when executing the program.

[0057] A non-transitory computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements a map generation method based on laser radar point cloud density analysis as described in any one of the items.

[0058] A computer program product comprises a computer program, wherein when the computer program is executed by a processor, the map generation method based on laser radar point cloud density analysis as described in any one of the items is implemented.

[0059] The advantages and beneficial effects of the present invention are as follows:

[0060] The present invention proposes a map generation method based on laser radar point cloud density analysis, which aims to solve the limitations of existing maps in characterizing complex environment spatial information, so as to improve the effective path planning ability of ground intelligent mobile platforms. The core innovation of the present invention lies in the accurate characterization of spatial density information, which is achieved by establishing a laser radar sensing model and density normalization processing. The specific implementation process includes: first, constructing a laser radar sensing model, and parameterizing the beam configuration, horizontal and vertical resolution; second, using a voxel processing method to spatially discretize the point cloud data, combining a ray tracing algorithm to calculate the beam flux of each voxel grid, and establishing a radar density model; then, the point cloud data acquired in real time is mapped to the corresponding voxel, and the voxel occupancy density is calculated and multiplied by the normalized density coefficient to obtain a standardized spatial density representation. The map model of the present invention not only contains normalized density information, but also integrates key terrain attributes such as height difference, slope, and roughness. This method improves the path planning efficiency and safety performance of intelligent mobile platforms in complex environments, and optimizes the accuracy of traffic area identification through accurate spatial density information representation. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] Figure 1 This is a general architecture diagram of a map generation method based on laser radar point cloud density analysis according to a preferred embodiment of the present invention;

[0062] Figure 2 It is a schematic diagram of the map voxel grid division;

[0063] Figure 3 is a flow chart of voxel density normalization calculation according to the present invention;

[0064] Figure 4 It is a schematic diagram of local terrain grid;

[0065] Figure 5 is the local terrain grid B p,q Schematic diagram of downsampling of ground points. DETAILED DESCRIPTION

[0066] The following will describe the technical solutions in the embodiments of the present invention in detail in conjunction with the accompanying drawings in the embodiments of the present invention. The described embodiments are only part of the embodiments of the present invention.

[0067] The technical solution of the present invention to solve the above technical problems is:

[0068] A map generation method based on LiDAR point cloud density analysis first establishes a map model and a LiDAR density model, and then respectively calculates the attributes representing the spatial information in real time.

[0069] Figure 1The figure shows a specific implementation process of a map generation method based on LiDAR point cloud density analysis, which includes the following steps:

[0070] (1) Map model

[0071] The map G is stored in a three-dimensional voxel grid, G = {T(i,j,k)|i∈[-N x ,N x ],j∈[-N y ,N y ],z∈[-H z ,H z ]}, T(i,j,k)=[i,j,k,roughness,step_height,density,slope]; where N x , N y , H z Respectively represent the boundary range of the x, y, and z axes of the map. T represents the voxel grid, where i, j, and k represent the coordinate positions of the three-dimensional voxel grid, roughness is the roughness, step_height is the step height, density is the normalized density, and slope is the terrain slope.

[0072] (2) Radar density model

[0073] The horizontal resolution of the laser radar is Δθ, and the vertical resolution is The direction of the laser radar light can be calculated, thus representing the laser radar of any beam. With the origin of the laser radar coordinate system as the center, the point cloud data is voxelized, and the density normalization coefficient of each voxel grid is calculated using the ray tracing method.

[0074] (2.1) Beam distribution mode setting

[0075] The vertical resolution is a set of different angles representing the perpendicular angle of each beam to the centerline of the lidar:

[0076]

[0077] The horizontal resolution Δθ represents the horizontal angle per increment. k = k*Δθ, where k is the index of the horizontal angle, and its value range is

[0078] (2.2) Calculation of beam direction:

[0079]

[0080] in Right now

[0081]

[0082] The above method can be used to represent a laser radar with any beam and resolution.

[0083] (2.3) Theoretical density calculation

[0084] Assume that the origin of the laser radar coordinate system is at a point outside the space (o x ,o y ,o z ), and the direction of the light is given by the vector d = (d x ,d y ,d z ). Any point on the ray (x t ,t t ,z t ) can be expressed as

[0085]

[0086] Here, t≥0 is the side length of the voxel grid.

[0087] For each ray, starting from the LiDAR position, the voxels passed along the ray direction are calculated one by one, and the number of times each voxel is passed through is recorded. t ,j t ,k t ) is calculated as

[0088]

[0089] The voxel grid beam flux is T(i,j,k):

[0090]

[0091] Where δ is an indicator function. When the condition is true, the δ value is 1, otherwise it is 0. is the density normalization coefficient.

[0092] (3) Data preprocessing: The point cloud data in the environment is collected by the laser radar device. The acquired raw point cloud data is filtered directly to retain the point cloud data of the area of ​​interest. With the laser radar as the center, its space is divided into N×N×H voxels. The collected point cloud data will be mapped to the corresponding voxel grid, and empty voxels without point clouds will be filtered out.

[0093] (4) Map attribute calculation: The preprocessed point cloud data is used to obtain the actual occupancy density by counting the number of point clouds in each voxel grid. The actual occupancy density is multiplied by the corresponding voxel grid density normalization coefficient to obtain the voxel normalized density. At the same time, the height difference between adjacent voxel columns is calculated to obtain the step height; the slope and roughness of the local terrain are calculated by plane fitting.

[0094] (4.1) Calculate the normalized voxel density

[0095] First, the laser radar point cloud data is mapped to the corresponding voxel. A point in the point cloud data (x i ,y i ,z i ) and voxel size v, the voxel index (i,j,k) is calculated as:

[0096]

[0097] Define a voxel set V, where each voxel contains O(i,j,k) points.

[0098]

[0099] Where P represents the point cloud dataset, δ is the indicator function, and when the condition is true, the δ value is 1, otherwise it is 0. O(i,j,k) represents the actual occupancy density of the voxel grid. Then the corresponding density normalization coefficient is obtained from the radar density model The normalized density of the voxel is then express.

[0100] (4.2) Calculate the voxel grid step height

[0101] In each column of voxel columns on the xy plane, the voxel grid calculates the height difference between it and the adjacent voxel column. For each voxel grid that may contain multiple point clouds, its height is defined as the average z value of all point clouds in the grid. Specifically, for each voxel grid point cloud set V e , whose grid height is for:

[0102]

[0103] Where n is the number of point clouds in the voxel grid. The height of the grid is represented by taking the average z coordinate of the point cloud. For each voxel column, the height difference attribute Δh between it and the adjacent voxel column is defined:

[0104] Δh=z(T)-z(T)′

[0105] Where z(T) and z(T)′ are the heights of the current voxel column and the adjacent voxel column, respectively. If Δh>0, it means the height difference is upward, and Δh<0, it means the height difference is downward.

[0106] (4.3) Calculate local terrain roughness and slope

[0107] First, we need to select ground points. In the map G, for the voxel column P(i,j)={T(i,j,k)|k∈[-H z ,H z ]}, when the voxel grid density is below a threshold: when density(T(i,j,k))<ρ low , T(i,j,k) is a non-ground point. When the voxel height exceeds the threshold:

[0108] Then P(i,j) is marked as a non-ground point, and the remaining voxels are marked as ground points. low is the density threshold, Hth is the height threshold, density(T) represents the normalized density of voxel T, and z(T) represents the height value of voxel T. Then the local terrain grid is divided, and the basic unit of the terrain grid is B p,q :B p,q ={T(i,j,k)∈G|i∈[pL,(p+1)L],j∈[qL,(q+1)L]}, where Represents the grid index in the x direction, Represents the grid index in the y direction. Then perform ground point downsampling and downsample the local terrain grid: B′ p,q =S(B p,q ,Δs), where Δs is the downsampling interval, in B p,q Select the set of downsampling points from: B′ p,q ={T(i,j,k)∈B p,q |i=pL+mΔs,j=qL+nΔs}, where is the row sampling index, is the column sampling index. Finally, plane fitting is performed on each local terrain. For each terrain grid B′ p,q The ground point set in Let the plane equation be:

[0109] ax+by+cz+d=0#(7)

[0110] Use the least squares method to fit the plane and get the plane coefficients a, b and d. The above calculations give the plane coefficients a, b and d. Then the residual error R′ of the fit is p,q That is

[0111]

[0112] Indicates the roughness of the ground.

[0113] (5) Generation of elevation map: The calculated terrain attributes and spatial density features are merged and filled into the map model, and the filled map is sent to the downstream module.

[0114] Figure 2 is a schematic diagram of the voxel grid division of the map described in the present invention, which shows the xoz plane projection of the map centered on the laser radar. Under this map model, the beam flux T(i, j, k) passing through each grid is calculated by the ray tracing method of equations (1) and (2) to obtain the normalized density coefficient

[0115] Figure 3 Schematic diagram of the method for normalizing voxel density of the present invention. First, the actual occupancy density of each voxel is calculated by segmenting, filtering and voxelizing the laser radar point cloud data. Then, the actual density is compared with the normalized density coefficient in the radar density model. Multiply them together to get the density value of the normalized density of the voxel grid.

[0116] (1) Point cloud data filtering

[0117] It is necessary to filter the point cloud data collected by the LiDAR. First, use a straight-through filter to retain the points within a certain coordinate axis range. The straight-through filter is as follows: Each point p i In three-dimensional space, there are coordinates (x i ,y i ,z i ); Use a pass filter to keep the point cloud within the x-axis range.

[0118] P filtered ={p i ∈P|x min ≤x i ≤x max}

[0119] That is, only keep min , x max ] i . For the y and z axes, the formula is similar. Then a radius filter is used to remove outliers, which determines whether to keep each point by checking the number of neighbors within a certain radius r. Assume that for a point p i , we define the radius r and the neighbor number threshold k. For each point p i ∈P, define its neighbor set as

[0120] N i ={pj ∈P|||p i -p j ||≤r}

[0121] where ||p i -p j || means p i With p j The Euclidean distance between them. If |N i | is less than the threshold k, then p i is an outlier, so it is deleted. Therefore, the result after filtering is

[0122] P filtered ={p i ∈P||N i |≥k}

[0123] (2) Voxel density calculation

[0124] First, the laser radar point cloud data is mapped to the corresponding voxel.

[0125] (x i ,y i ,z i ) and voxel size v, the voxel index (i,j,k) is calculated as:

[0126]

[0127] Define a voxel set V, where each voxel contains O(i,j,k) points.

[0128]

[0129] Where P represents the point cloud dataset, δ is the indicator function, and when the condition is true, δ is 1, otherwise it is 0. O(i,j,k) represents the actual occupancy density of the voxel grid. Then the corresponding density normalization coefficient T(i,j,k) is obtained from the radar density model, and the normalized density of the voxel is used express.

[0130] Figure 4 It is a schematic diagram of the local terrain grid division of the map model G. This figure shows the three-dimensional voxel grid structure G = {T(i,j,k)|i∈[-N x ,N x ],j∈[-N y ,N y ],z∈[-H z ,H z ]}, the grid division of the xy plane is realized. In the figure, the xy plane is divided into multiple local terrain grids B with equal spacing L×L. p,q, where p and q represent the grid indices in the x and y directions, respectively. Each local grid cell B p,q It contains all the original voxel data T(i,j,k) in the area and its attribute information, including roughness, step height (step_height), normalized density (density) and terrain slope (slope).

[0131] Figure 5 For the local terrain grid B p,q Schematic diagram of ground point downsampling in Figure 1. The left figure shows the distribution of ground points in the original voxel grid, which has a high spatial sampling density; the right figure shows the result after adopting the equal-interval downsampling strategy, where the sampling interval Δs ensures that the representativeness of terrain features is maintained while reducing the amount of data. Point set B after downsampling p,q The key features of the original terrain are retained, providing a data basis for subsequent plane fitting and roughness calculation. This downsampling method effectively reduces data redundancy and improves processing efficiency.

[0132] The systems, devices, modules or units described in the above embodiments may be implemented by computer chips or entities, or by products with certain functions.

[0133] Computer readable media include permanent and non-permanent, removable and non-removable media that can be implemented by any method or technology to store information. Information can be computer readable instructions, data structures, program modules or other data. Examples of computer storage media include, but are not limited to, phase change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technology, compact disk read-only memory (CD-ROM), digital versatile disk (DVD) or other optical storage, magnetic cassettes, magnetic tape magnetic disk storage or other magnetic storage devices or any other non-transmission media that can be used to store information that can be accessed by a computing device. As defined herein, computer readable media does not include temporary computer readable media (transitory media), such as modulated data signals and carrier waves.

[0134] It should also be noted that the terms "include", "comprises" or any other variations thereof are intended to cover non-exclusive inclusion, so that a process, method, commodity or device including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, commodity or device. In the absence of more restrictions, the elements defined by the sentence "comprises a ..." do not exclude the existence of other identical elements in the process, method, commodity or device including the elements.

[0135] The above embodiments should be understood to be only used to illustrate the present invention and not to limit the protection scope of the present invention. After reading the contents of the present invention, technicians can make various changes or modifications to the present invention, and these equivalent changes and modifications also fall within the scope defined by the claims of the present invention.

Claims

1. A map generation method based on laser radar point cloud density analysis, characterized in that: It includes the modeling stage and the mapping stage, in which: The modeling stage is used to establish a map model and a radar density model; the map model contains the geometric information and spatial density information of the terrain, and is composed of a three-dimensional voxel grid. Each voxel grid contains attributes such as coordinate information, voxel density, step height, terrain slope and roughness; the radar density model is used to characterize the spatial density normalization coefficient of the point cloud, and the emission direction of each laser beam is calculated by setting the horizontal resolution and vertical resolution of the laser radar; voxelization is performed at the radar coordinate origin, and the beam flux of each voxel grid is calculated by ray tracing, and the inverse of the beam flux is used as the voxel density normalization coefficient of the voxel grid; In the mapping stage, the laser radar is used to collect environmental data in real time, and the collected point cloud is preprocessed. The point cloud of the area of ​​interest is obtained after straight-through filtering and radius filtering, and voxelized, and the point cloud data is mapped to the corresponding voxel grid; then the map attributes of each voxel grid are calculated to obtain attributes including voxel density, terrain slope and roughness, and the above information is combined to construct a quasi-elevation map.

2. The map generation method based on laser radar point cloud density analysis according to claim 1 is characterized in that: The map model is: The map G is stored in a three-dimensional voxel grid, G = {T(i, j, k)|i∈[-N x , N x ],j∈[-N y , N y ],z∈[-H z , H z ]}, T(i, j, k) = [i, j, k, roughness, step_height, density, slope]; where N x , N y , H z They represent the boundary range of the x, y, and z axes of the map respectively, T represents the voxel grid, i, j, and k represent the coordinate positions of the three-dimensional voxel grid, roughness represents the roughness, step_height represents the step height, density represents the normalized density, and slope represents the terrain slope.

3. The map generation method based on laser radar point cloud density analysis according to claim 2 is characterized in that: The establishment of the radar density model specifically comprises the following steps: The radar density model is stored in a three-dimensional voxel grid. G={T(i,j,k)|i∈[-N x ,N x ],j∈[-N y ,N y ],k∈[-H z ,H z ]},T(i,j,k)=[i,j,k,ρ]; Where ρ is the voxel density normalization coefficient; 3.1 First, set the beam distribution mode, and the vertical resolution is a set of different angles, which represent the vertical angle between each beam and the center line of the LiDAR: vertical distribution angle is a set of preset values ​​{θ1, θ2, ..., θ i } represents the vertical angle between each laser beam and the center line, where i represents the number of laser beams; and the horizontal distribution angle θ k , then through θ k = k*Δθ calculation, where k is the index of the horizontal angle, and its value range is Δθ represents the angular interval between adjacent scanning points in the horizontal direction, and then the beam direction is calculated 3.2 The origin of the laser radar coordinate system is at a point outside the space (o x , o y , o z ), and the direction of the light is given by the vector d = (d x , d y , d z ) is given; any point on the ray (x t ,y t , z t ) can be expressed as Where t≥0 is the side length of the voxel grid; For each laser beam, starting from the laser radar position, the voxels passed along the direction of the light are calculated in sequence, and the number of times each voxel is passed through is recorded; the voxel index (i t , j t , k t ) is calculated as The voxel grid beam flux is T(i, j, k): Among them, δ is the indicator function. When the condition is true, the δ value is 1, otherwise it is 0; is the density normalization coefficient. ray represents each laser beam, which is determined by the previous d k，i The beam direction is determined.

4. The map generation method based on laser radar point cloud density analysis according to claim 1, characterized in that: The constructed class elevation map, The following steps are involved: 4.1 Use the LiDAR device to collect point cloud data in the environment, and pre-process the data using straight-through filtering and radius filtering to retain the point cloud in the area of ​​interest; 4.2 Perform voxel processing, map the point cloud of the area of ​​interest to the corresponding voxel grid, and calculate the number of point clouds in each voxel grid, multiply it by the density normalization coefficient of the corresponding voxel grid in the radar density model to obtain the normalized density of the voxel; 4.3 Taking the origin of the laser radar coordinate system, differentiate each column of voxel columns on the xoy plane to obtain the step height between adjacent voxel columns; 4.4 Classify each voxel into ground points and non-ground points, and then divide the entire map into multiple L×L local areas; perform plane fitting on each local area, the normal vector of the plane fitting represents the slope of the area, and the residual after fitting represents the roughness of the terrain; 4.5 The calculated voxel normalized density, step height, terrain slope and roughness are combined to generate a quasi-elevation map.

5. The map generation method based on laser radar point cloud density analysis according to claim 4 is characterized in that: The step 4.2 voxel grid normalized density calculation comprises the following steps: 5.1 Calculation of actual voxel occupancy density The process of voxelizing point clouds is to map the collected point cloud data into a discrete voxel grid; point cloud (x i ,y i , z i ) and voxel size v, the voxel index (i, j, k) is calculated as: Define a voxel set V, where each voxel contains O(i, j, k) points. Where P represents the point cloud dataset, δ is the indicator function, and when the condition is true, the δ value is 1, otherwise it is 0; 5.2 Density Normalization Normalized density D(i, j, k), which is the actual occupancy density O(i, j, k) of voxel (i, j, k) and the density normalization coefficient The product of 6. According to the map generation method based on laser radar point cloud density analysis of claim 4, it is characterized in that the step 4.4 classifies each voxel into ground points and non-ground points, and then divides the entire map into multiple L×L local areas; plane fitting is performed on each local area respectively, the normal vector of the plane fitting represents the slope of the area, and the residual after fitting represents the roughness of the terrain, specifically comprising the following steps: 6.1 Ground point selection In the map G, for a voxel column P(i, j) = {T(i, j, k) | k∈[-Hz, Hz]}, when the voxel grid density is below a threshold: when density(T(i, j, k)) < ρ low , T(i, j, k) is a non-ground point; when the voxel height exceeds the threshold: Then P(i, j) is marked as a non-ground point, and the remaining voxels are marked as ground points; where ρ low is the density threshold, Hth is the height threshold, density(T) represents the normalized density of voxel T, z(T) represents the height value of voxel T; H z Indicates the height area of ​​the map. 6.2 Local terrain meshing The basic unit of the terrain grid is B p，q :B p，q ={T(i, j, k)∈G|i∈[pL, (p+1)L], j∈[qL, (q+1)L]}, where Represents the grid index in the x direction, Indicates the grid index in the y direction; L represents the basic unit B of the terrain grid p，q length and width. 6.3 Ground Point Downsampling Downsample the local terrain grid: B′ p，q =S(B p，q , Δs), where Δs is the downsampling interval, at B p，q Select the set of downsampling points from: B′ p，q ={T(i, j, k)∈B p，q |i=pL+mΔs,j=qL+nΔs}, where is the row sampling index, is the column sampling index; B′ p，q Represents the local terrain after sampling. 6.4 Plane Fitting For each terrain grid B′ p，q The ground point set in Let the plane equation be: ax+by+cz+d=0#(7) Use the least square method to fit the plane and obtain the plane coefficients a, b and d; 6.5 Residual Calculation The above calculations yield the plane coefficients a, b, and d; then the residual R′ of the fit p，q That is Indicates the roughness of the ground.

7. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the program, the map generation method based on lidar point cloud density analysis as described in any one of claims 1 to 6 is implemented.

8. A non-transitory computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the map generation method based on laser radar point cloud density analysis as described in any one of claims 1 to 6 is implemented.

9. A computer program product, comprising a computer program, characterized in that When the computer program is executed by a processor, the map generation method based on laser radar point cloud density analysis as described in any one of claims 1 to 6 is implemented.

Images