CN117152215A - Automatic registration method for fast vehicle-mounted laser point cloud in complex urban environment - Google Patents

Abstract

The invention relates to the technical field of autonomous navigation of robots, in particular to a rapid vehicle-mounted laser point cloud automatic registration method in a complex urban environment. Global fine registration was then performed using ICP plane-to-plane algorithm modified by the face-to-face and Levenberg-Marquardt (LM) algorithm. Not only reduces the requirement on the initial position, but also accelerates the convergence speed of registration. The improved Iterative Closest Point (ICP) algorithm and normal vector constraint take the rod object and edge characteristic points as registration primitives, reject error point pairs by using SAC-IA algorithm, and use the characteristic point relation corresponding to the two-way KD tree rapidly, so as to accelerate the registration speed and improve the accuracy and enhance the robustness.

Description

Automatic registration method for fast vehicle-mounted laser point cloud in complex urban environment

Technical Field

The invention relates to the technical field of autonomous navigation of robots, in particular to a rapid vehicle-mounted laser point cloud automatic registration method in a complex urban environment.

Background

The three-dimensional laser point cloud scanning technology has the advantages of high effective rate, wide range and accurate precision in the process of three-dimensional modeling. In recent years, unmanned and mobile robots have developed rapidly, and laser synchronous positioning and drawing (SLAM, simultaneousLocalizationAndMapping) are widely applied to various scenes. The point cloud registration technology is an important step of three-dimensional modeling, and is to scan point cloud data from different angles to enable the point cloud data to be in the same coordinate system. The method has great effect on the development of the fields of reverse engineering, VR scene, unmanned mapping, cultural relic restoration, medical treatment and the like, and the research significance of the point cloud registration technology also becomes more important.

The traditional point cloud registration method is a manual or semi-automatic mode, for example, three or more targets are distributed in a scene, and the common target positions of the cloud data of two adjacent sites are manually or semi-automatically identified, so that a transformation matrix is obtained to complete registration. In practical application, due to the limitations of scenes, cost and the like, the target layout cannot be performed, and the registration process is complicated and the efficiency is low by using the method. Therefore, it is necessary to research a target-free, automatic, rapid and efficient point cloud registration method, and a good registration result is formed, so that subsequent three-dimensional reconstruction, target recognition segmentation and other works of the point cloud are facilitated.

Algorithms based on global search strategies and algorithms based on statistical probability are currently mainly used. The algorithm of the global search strategy is to search for a point pair meeting a certain limiting condition in the global range of the point cloud as a corresponding point pair by setting the certain limiting condition, and select an optimal transformation matrix through iteration. The most classical algorithm is an iterative nearest point algorithm (ICP, iterativeClosestPoint) proposed by Besl and Mckay (1992), points with the nearest distance in another point cloud are continuously and iteratively searched in a global range to serve as corresponding points until the sum of the distances of all the points is minimum, the algorithm requires that two point clouds have good initial positions and are in a containing relation, otherwise, the two point clouds are prone to being in local optimum, and a good effect cannot be achieved. Therefore, many researchers improve the ICP algorithm and only use it for accurate registration after obtaining good initial positions, the efficiency is lower, and the robustness is high. And the subsequent proposal replaces the point-to-point distance of ICP by calculating the point-to-plane distance, thereby reducing the iteration times. In addition, the four-point consistent set algorithm (4 PCS, 4-PointsCongrentSets) proposed by Aiger et al (2008) is also an algorithm with higher efficiency based on a global search strategy, and the algorithm uses a RANSAC algorithm framework to construct four coplanar points by utilizing affine invariance constraint and searches the best matched congruent four-point set in a target point cloud, so that a transformation matrix is calculated, but the efficiency is lower and the calculation complexity is higher.

The algorithm based on statistical probability refers to estimating the distribution of the point clouds by using a probability density function, and obtaining a transformation matrix between two pieces of point clouds by using an optimization technology method. Compared with the traditional ICP algorithm, the registration effect of the algorithm is more excellent. The Biber et al (2003) uses a normal distribution transformation algorithm (NDT, normalDistributionsTransform) for two-dimensional point cloud data registration for the first time, and achieves good effects. Magnusson et al (2007) improves on the basis of a 2D-NDT algorithm, and proposes a 3D-NDT algorithm for point cloud registration, which converts the registration problem into a nonlinear optimization problem by using the point cloud as a gaussian distribution function, and requires an initial estimation before registration, otherwise, a good effect is not obtained. In addition, jian et al (2011) proposes a new probability framework, namely, the point cloud is regarded as a gaussian mixture model (GMM, gaussianMixtureModel), and conversion between two point clouds is estimated through a GMM, but the algorithm has high computational complexity, long time consumption and sensitivity to the point clouds with noise or uneven density distribution.

The vehicle-mounted laser scanning system for acquiring the urban scene data has the problems of complex ground features, multiple targets, unstructured description, multiple discrete points, low reflection distance errors, low precision and the like, and the method still has the limitations of sensitivity to noise points, higher complexity, low registration efficiency, weak robustness and the like.

Disclosure of Invention

The invention aims to provide a rapid vehicle-mounted laser point cloud automatic registration method in a complex urban environment. The method particularly relates to a laser synchronous positioning and drawing (SLAM, simultaneousLocationandMapping) technology, and aims to realize three-dimensional reconstruction by registering and fusing upper and lower frame scanning data. According to the long-term excellent stability and extraction consistency of the rod in the urban road scene, the rod is firstly extracted to serve as a target feature, and the point cloud with the ground points filtered is converted into the depth map to extract edge feature points, so that the requirement can be met under the condition that the rod is restrained less. Global fine registration was then performed using ICP plane-to-plane algorithm modified by the face-to-face and Levenberg-Marquardt (LM) algorithm. Not only reduces the requirement on the initial position, but also accelerates the convergence speed of registration. The improved Iterative Closest Point (ICP) algorithm and normal vector constraint take the rod object and edge characteristic points as registration primitives, reject error point pairs by using SAC-IA algorithm, and use the characteristic point relation corresponding to the two-way KD tree rapidly, so as to accelerate the registration speed and improve the accuracy and enhance the robustness.

In order to achieve the technical purpose and the technical effect, the invention is realized by the following technical scheme:

A rapid vehicle-mounted laser point cloud automatic registration method in a complex urban environment comprises the following steps:

s1: in-vehicle laser scanning typically generates point cloud datasets with different point densities, performs point cloud filtering and voxel downsampling;

s2: filtering a cloth model, and setting a distance threshold value to quickly separate ground points;

s3: extracting the rod-shaped object as a registration element;

s4: extracting edge key points by using a Normal Aligned Radial Feature (NARF) algorithm;

s5: global fine registration was performed using ICP plane-to-plane algorithm modified by the face-to-face and Levenberg-Marquardt (LM) algorithm.

Further, the point cloud filtering and voxel downsampling specifically include: in-vehicle laser scanning typically generates point cloud data sets with different point densities, creating outliers for which some of the information is unusable. The statistical filter is mainly used for filtering outliers in the point cloud, calculating the average distance from each point to the nearest k points, and assuming that the obtained result is a gaussian distribution, the shape of the gaussian distribution is determined by the mean value and the standard deviation, the points with the average distance outside the standard range can be defined as outliers and removed from the data. The probability density function of the Gaussian random variable n is

Wherein: u represents the distance from which the object is to be moved,mean of all u, σ standard deviation.

Distance threshold d _max Can be expressed as

d _max ＝u+α·σ (2)

α is a constant, which may be referred to as a scaling factor, which depends on the number of neighbor points; finally, traversing the point cloud again, and eliminating that the average distance between the point cloud and k neighbor points is greater than d _max Is a point of (2).

Voxel size calculation is the key point of voxelized sampling, and a reasonable voxel size can effectively simplify the original point cloud. Firstly, determining a proper L value, wherein the size of the L value can directly influence the downsampling rate of the point cloud, and when L is overlarge, the points falling into voxels become more, the sampling rate is low, so that the local characteristics of the point cloud are lost; when the L value is small, the number of points in the voxel is small, the sampling rate becomes high, and although the local characteristics can be well reserved, the sampling result is not obvious. To solve this problem, the combination of the downsampling rate ρ and the point cloud density estimates the voxel side length L, i.e

Wherein: l' represents the estimated voxel side length; s is the average point cloud number within a unit voxel; gamma is the point cloud density. s and gamma are expressed as

Wherein: n is the number of point cloud points. Substituting the formulas (3) and (4) into the formula (2) to obtain

Obtaining a proper voxel side length L to divide a bounding box of the point cloud P, calculating the centers of gravity of all points in the voxel for better retaining the feature of the point cloud, and using the calculated center of gravity point P _k ' as downsampling points. P (P) _k ' is denoted as

Wherein:to calculate the number of k-th voxel inliers, k ε [1,2, …, m]. A new point cloud P ', P' can be obtained from equation (6), which is the final downsampled point cloud.

Further, the step S2 specifically includes: and inverting the point cloud, and covering the inverted curved surface with cloth from above. And determining the final shape of the cloth model by analyzing the interaction between the nodes and the corresponding points of the cloth model, and dividing the ground and the non-ground according to the shape of the final shape.

Further, the extracting rod specifically includes: firstly, non-ground point clouds are classified in an unsupervised mode by using a K mean value, and then, the road trees, the street lamps, the identification plates and the telegraph poles of the scene point clouds are marked after classification, and trunks, the street lamps, the identification plates and the telegraph poles of the road trees are extracted to be rod-shaped targets. In order to ensure the extraction effect of the rod, firstly, a box packing method is used for establishing a space grid for the classified point clouds, and the superfluous point clouds higher than the rod are removed by utilizing an elevation threshold; traversing the grid aerial view by using a random sampling consistency algorithm (RANSAC), setting a cylindrical radius by taking a trunk as a center, and searching an arc-shaped point set for the radius neighborhood; and finally, judging the grid density, the minimum number of support points and the angle threshold value perpendicular to the vertical face, and regarding the grid as rod-shaped point cloud meeting the requirements. The lamp cap part of the urban street lamp has a horizontal extension area, the number of point clouds of the lamp cap part and the signboard part is small, and the rod-shaped point clouds can be accurately extracted by setting different densities and minimum supporting point cloud number thresholds.

According to the coordinate extreme value in the point cloud, an m multiplied by n grid with the side length L is constructed, the point cloud data are divided into a regular two-dimensional grid, and the m and n calculation modes are as follows:

wherein: x is X _max 、X _min 、Y _max And Y _min For the extreme value of the point cloud value, L is the grid width, and the invention is set to be 1 meter according to the setting of the point cloud density.

And searching a fitting circle by adopting a RANSAC algorithm to extract the cylindrical point cloud, and setting the range of a radius buffer area to be 0.05-0.5 meter to extract the cylindrical point cloud.

Wherein:for the coordinates of each point, +.>For the center of each fitting circle, R _i Radius d _i E for the distance to each cylinder _i Representing the degree of difference of the grid blocks from the fitted circle.

Further, the extracting the edge key point specifically includes: edge keypoints are extracted using a Normal Aligned Radial Features (NARF) algorithm. NARF is a 3D feature point detection and description algorithm, and a distance image generated by using point cloud is used for extracting a boundary from a foreground to a background, and key points are proposed by identifying objects from the distance image, so that the operation efficiency is higher than that of directly operating on the point cloud. The extraction process considers the edge and object surface change information, the position of the key point can be stably and repeatedly detected at different visual angles, a stable supporting area is needed, and a descriptor and a unique normal vector can be calculated.

The edge key point extraction method specifically comprises the following substeps:

s4.1: and traversing each depth image point, and searching a position with depth mutation in a neighboring area for edge detection. For a point p in each range image _i Selecting a square window field point in the same plane, and marking asWhere s is the window width. Calculation of p _i And->The 3D distance (euclidean distance) of each point in the interior, resulting inFor->Arranged in ascending order to obtain ∈ ->Let N be _pi At least M points and p _i On the same surface, select δ=d _M ' as p _i At p _i -distances of non-coplanar points within s-neighborhood. For p _i Four confidence levels are calculated, up, down, left, right, representing the likelihood of an edge being present, the following being taken as an example on the right:

firstly, calculating the average distance of the right point:wherein m is _p To consider the number of right dots, p _x,y Representing points at (x, y) positions in the image. Calculating 3D distance D _right ＝||p _x,y -p _right Sum of s _right ∈[0,1)：

Wherein: s is(s) _right The larger the value, the more marked point p _i The more pronounced the right distance change, the greater the probability that an edge will be present.

For the obtained s _right Smoothing the graph to obtain continuous edges while improving noise robustness, determining p _x,y The type of point is according to p _x,y And p is as follows _right P _x,y A smaller distance indicates that an object boundary is possible, and vice versa, a shadow boundary is obtained, and the window width s corresponding to the distance between the shadow boundary is calculated _right ′＝max(0.9,(1-s _shadow ) ³ )×s _right If s _right ' greater than the threshold and is (p) _x-1,y ,p _x+1,y ) Maximum in the field, p is considered to be _x,y For object boundaries, shadow points and shadow boundaries are found according to the formula, and the shadow points, namely the extracted interest feature points, are located between the shadow points.

S4.2: the coefficient of measuring the surface variation, and the principal direction of variation, is determined from the surface variation of the neighboring region. For each point p _i The method belongs to a range image, and a normal vector of the range image is estimated by using a 3D position of a point with Euclidean distance within 2 delta in the 2D field and a Principal Component Analysis (PCA) algorithm. And calculating a main direction v and an amplitude lambda, wherein v is an edge direction for boundary points, v is a main curvature direction for other points, and lambda is the amplitude of the direction. For each v, a weight ω∈0, 1) is assigned, for the boundary point ω=1; for other points ω=1- (1- λ) ³ . One parameter sigma is specified to represent the supportable range size. For each point p in the image, find it so that the Euclidean distance is atPoints in which there is no edge in between (two points on the same surface), denoted as N _p ＝{p ₀ ,p ₁ ,…,p _N }。

For a pair ofWill v _pi Projected onto a plane P to obtain an angle alpha _pi For alpha _pi The following adjustments were made:

so that alpha' e (-180 deg., 180 deg.).

S4.3: and smoothing and filtering the interest value. Smoothing ω and α using bounded gaussian kernels to obtain a point of interest I for point p _p Points exceeding the threshold I are taken as interest points:

I(P)＝I ₁ (P)·I ₂ (P) (14)

wherein: i ₁ (P) high weight point values, i.e. points where the plane changes drastically, can be extracted; i ₂ (P) points in the vicinity where adjacent points exist and the main direction changes drastically can be added.

S4.4: and carrying out maximum value-free compression to find a final key point, namely the NARF key point.

Further, the step S5 specifically includes:

assume that two sets of point clouds a= { a ₀ ,a ₁ ,…,a _N Sum b= { B ₀ ,b ₁ ,…,b _N Performing nearest neighbor distance lookup according to their corresponding index relationship such thatThis holds true, where T is the transform stiffness matrix and E (T) is the error loss functionThe number d is the Euclidean distance threshold. In the probabilistic model, it is assumed that there are two point sets +.>Andaccording to->And->Yielding a and B which, in this case,and->Is the covariance matrix associated with the measurement points. Assume that T is ^* To correctly transform the rigid matrix, so as to obtain perfect corresponding relation

For any rigid matrix T, defineAnd deriving +.>Is given by a _i And b _i Derived from an independent Gaussian distribution, by the above formula

Iterative calculation of T using maximum likelihood estimation

Simplifying the above steps into

Wherein: and t opt is an updated rigid transformation matrix, and t=r·t is a rigid matrix obtained by LM decomposition.

To improve performance with respect to point-to-plane and increase symmetry of the model, point-to-plane is used to excite the probabilistic model. Essentially, each measurement point provides only constraints along its surface normal. To construct such a model, each sample point is considered to have a high covariance distribution along its local plane and a very low covariance in the direction of the surface normal. If the surface normal of the point is n _i The covariance matrix becomes

Wherein: alpha represents a small constant of covariance along the normal. This corresponds to a very high confidence and knows the position along the normal, but does not determine its position in the plane.

Wherein: s is(s) _i ，t _i Is a as _i And b _i Normal vector of corresponding point, R _x Is a rotation matrix for transforming the basis vectors.

And substituting (20) and (21) into the above formula (18) to calculate the transformation rigid matrix T opt. Namely, the construction of the ICP algorithm from plane to plane is completed;

to reduce the effects of the point cloud pairing process noise and outliers, the Cauchy loss function is utilized as a weighting function to resist noise and outliers. The modified ICP algorithm is as follows:

wherein: ρ _Cauchy For Cauchy loss function weights, k is the loss parameter, and typically k is 1.345 and r is the residual at the time of optimization.

R obtained by decomposing LM is substituted into the above formula to obtainAnd->Then substituting to obtain updated transformation matrix T opt and obtaining +.>By this iteration +.>Less than the set Euclidean distance threshold τ or the transform update matrix is made less than the set thresholdSetting the iteration times to 200 times, and stopping iteration when one of the conditions is met.

The euclidean distance RMSE of the point pairs is used as the registration error accuracy for evaluation. RMSE obtainable from equation (22) is as follows:

wherein: s is(s) _i ＝(s _ix ,s _iy ,s _iz ) Coordinates of cloud points of source points, t _i ＝(t _ix ,t _iy ,t _iz ) And n is the number of registered point pairs for the coordinates of the point cloud to be registered.

Further, the step S5 further includes: the method comprises the steps of using an improved Iterative Closest Point (ICP) algorithm and normal vector constraint, using rod objects and edge feature points as registration primitives, using a sampling geometric consistency algorithm to remove error point pairs, and using a two-way KD tree to quickly correspond to the relation of the feature points, so that the registration speed is increased and the precision is improved.

The invention has the beneficial effects that:

1. fast registration: the method realizes the rapid automatic registration of the vehicle-mounted laser point cloud through an improved ICP plane-to-plane algorithm and a Levenberg-Marquardt (LM) algorithm. The improved ICP algorithm uses an iterative closest point matching strategy to achieve fast and accurate registration by matching planes in the point cloud. Compared with the traditional ICP algorithm, the method can be used for converging to the optimal solution more efficiently while maintaining the matching effect, and the registration speed is greatly improved.

2. High-precision three-dimensional reconstruction: the method utilizes the extracted rod-shaped objects and the edge characteristic points as registration primitives, and realizes accurate three-dimensional reconstruction through high-precision registration. First, the extracted shaft has long term excellent stability and extraction consistency, which can be a reliable registration feature. And secondly, the registration process is more accurate by adopting normal vector constraint, so that the accuracy of the reconstruction result is improved. In this way, the method can realize highly accurate three-dimensional reconstruction in a complex urban environment.

3. The robustness is strong: the method improves the robustness of registration through an improved ICP algorithm and the use of a two-way KD tree. The improved ICP algorithm uses face-to-face and LM algorithms for global registration to reduce registration errors. The corresponding relation of the feature points is quickly established in the registration process of the two-way KD tree, so that the registration speed is increased and the registration accuracy is improved. In addition, error point pairs are removed by using SAC-IA algorithm, so that the robustness of registration is further enhanced, and reliable registration results are ensured to be obtained in complex urban environment.

4. The requirement for the initial position is low: according to the method, the shaft and the edge feature points are used as registration primitives, and the requirement on the initial position is reduced in the registration process. The shaft acts as a reliable registration feature, providing stable constraint information, reducing dependence on initial position. This means that the method can still obtain high quality registration results even if the initial position is not completely accurate. The reduction of the requirement on the initial position improves the robustness and the practicability of the method, and is particularly suitable for scenes in which the accurate initial position is difficult to acquire in practical application.

In conclusion, the method has the characteristics of quick registration, high-precision three-dimensional reconstruction, strong robustness and low requirement on the initial position, and provides an effective solution for automatic registration of the vehicle-mounted laser point cloud in a complex urban environment.

Of course, it is not necessary for any one product to practice the invention to achieve all of the advantages set forth above at the same time.

Drawings

In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings that are needed for the description of the embodiments will be briefly described below, and it is obvious that the drawings in the following description are only some embodiments of the present invention, and that other drawings may be obtained according to these drawings without inventive effort for a person skilled in the art.

FIG. 1 is a schematic overall flow chart of the present invention;

fig. 2 is a schematic diagram of an initial position of a two-point cloud according to an embodiment of the present invention;

fig. 3 is a diagram of ICP registration results according to an embodiment of the invention;

FIG. 4 is a graph of point-to-face ICP registration results according to an embodiment of the present invention;

FIG. 5 is a graph showing the results of the method according to the present invention

Fig. 6 is a schematic diagram of a cloth model filtering algorithm according to an embodiment of the present invention;

FIG. 7 is a schematic view of a point cloud color chart according to an embodiment of the present invention;

FIG. 8 is a schematic view of a separated ground point cloud according to an embodiment of the present invention;

FIG. 9 is a schematic view of a non-ground point cloud according to an embodiment of the present invention;

FIG. 10 is a schematic diagram of a non-ground point cloud classification result according to an embodiment of the present invention;

FIG. 11 is a schematic view of the ground of a shaft according to an embodiment of the present invention;

fig. 12 is a schematic view of projection extraction of a tree trunk and a street lamp according to an embodiment of the present invention; (a) light pole extraction; (b) trunk extraction;

FIG. 13 is a schematic view of the result of extraction of a shaft and partial details according to an embodiment of the present invention;

FIG. 14 is a schematic view of NARF extracted edge feature points according to an embodiment of the present invention;

FIG. 15 is a schematic view of a plane-to-plane matching according to an embodiment of the present invention;

Detailed Description

The following description of the embodiments of the present invention will be made clearly and completely with reference to the accompanying drawings, in which it is apparent that the embodiments described are only some embodiments of the present invention, but not all embodiments. All other embodiments, which can be made by those skilled in the art based on the embodiments of the invention without making any inventive effort, are intended to be within the scope of the invention.

Example 1

In order to verify the feasibility of the proposed invention, a large city outdoor point cloud data set Toronto-3D was selected. The point cloud in the data set is obtained through a vehicle-mounted laser scanning system TeledyneOptechMaverick 2. The system consists of a 32-line laser radar sensor, an ladybug5 panoramic camera, a GNSS system and an SLAM system. The invention selects the regional data (21567172 points) set, and the coverage range is about 250 meters. Compared with the ICP and the point-to-surface ICP in the prior art, the method disclosed by the invention is high in efficiency, high in precision and better in robustness, and can be used for realizing efficient and accurate registration of the vehicle-mounted laser point cloud of the complex urban scene.

A rapid vehicle-mounted laser point cloud automatic registration method in a complex urban environment comprises the following steps:

s1: in-vehicle laser scanning typically generates point cloud datasets with different point densities, performs point cloud filtering and voxel downsampling;

s2: filtering a cloth model, and setting a distance threshold value to quickly separate ground points;

s3: extracting the rod-shaped object as a registration element;

s4: extracting edge key points by using a Normal Aligned Radial Feature (NARF) algorithm;

s5: global fine registration was performed using ICP plane-to-plane algorithm modified by the face-to-face and Levenberg-Marquardt (LM) algorithm.

Example 2

When the number of the point clouds is large, the traditional ICP algorithm has low registration efficiency and poor fusion effect, the algorithm operation time is long and the convergence speed is low by using a search mode of adjacent point pairs, the iteration threshold number is reached, the method is sensitive to noise, the robustness is poor, and obvious position deviation can be seen from fig. 2-5; the registration result error of the point-to-surface ICP algorithm is reduced greatly, the registration efficiency is improved by 24 seconds although the convergence is faster than that of the traditional ICP algorithm, the iteration threshold times are reached, the effect is not good, and obvious position deviation still exists. Under the complex urban environment with lower initial overlapping degree and noise, the method provided by the invention can be seen to be completely fused, is less in time consumption and rapid in convergence, and can reach the set error threshold after 20 times of iteration, so that compared with the traditional registration method, the efficiency, registration precision and robustness are obviously improved, and the light-weighted rapid and accurate point cloud registration is realized.

Table 1 registration accuracy and runtime of algorithms

Example 3

Point cloud filtering and voxel downsampling

In-vehicle laser scanning typically generates point cloud data sets with different point densities, creating outliers for which some of the information is unusable. The statistical filter is mainly used for filtering outliers in the point cloud, calculating the average distance from each point to the nearest k points, and assuming that the obtained result is a gaussian distribution, the shape of the gaussian distribution is determined by the mean value and the standard deviation, the points with the average distance outside the standard range can be defined as outliers and removed from the data. The probability density function of the Gaussian random variable n is

Wherein: u represents the distance from which the object is to be moved,mean of all u, σ standard deviation.

Distance threshold d _max Can be expressed as

d _max ＝u+α·σ (2)

α is a constant, which may be referred to as a scaling factor, which depends on the number of neighbor points; finally, traversing the point cloud again, and eliminating that the average distance between the point cloud and k neighbor points is greater than d _max Is a point of (2).

Voxel size calculation is the key point of voxelized sampling, and a reasonable voxel size can effectively simplify the original point cloud. Firstly, determining a proper L value, wherein the size of the L value can directly influence the downsampling rate of the point cloud, and when L is overlarge, the points falling into voxels become more, the sampling rate is low, so that the local characteristics of the point cloud are lost; when the L value is small, the number of points in the voxel is small, the sampling rate becomes high, and although the local characteristics can be well reserved, the sampling result is not obvious. To solve this problem, the combination of the downsampling rate ρ and the point cloud density estimates the voxel side length L, i.e

Wherein: l' represents the estimated voxel side length; s is the average point cloud number within a unit voxel; gamma is the point cloud density. s and gamma are expressed as

Wherein: n is the number of point cloud points. Substituting the formulas (3) and (4) into the formula (2) to obtain

Obtaining a proper voxel side length L to divide a bounding box of the point cloud P, calculating the centers of gravity of all points in the voxel for better retaining the feature of the point cloud, and using the calculated center of gravity point P _k ' as downsampling points. P (P) _k ' is denoted as

Wherein:to calculate the number of k-th voxel inliers, k ε [1,2, …, m]. A new point cloud P ', P' can be obtained from equation (6), which is the final downsampled point cloud.

Example 4

Cloth model filtering (CFS) ground separation

As shown in fig. 6-9, the ground point information is less and is not suitable as a characteristic point, and removing the ground point not only reduces the redundancy of the point cloud data, but also accelerates the registration and classification extraction accuracy. Generally, the gradient fluctuation of the urban road ground point cloud is small, so that the filtering of a distribution model is selected, and the ground points are quickly separated by setting a distance threshold. An overview of the algorithm is shown in fig. 6. Firstly, the point cloud is inverted, and then the cloth model is covered on the inverted curved surface from above. The final shape of the cloth is determined by analyzing the interaction between the nodes and corresponding points of the cloth, and the ground and non-ground are divided according to the shape thereof.

Example 5

Extracting the rod-like substance

As shown in fig. 10-13, the shaft is rich in urban ground object types, relatively stable and regular, has more remarkable object characteristics, and is suitable for being used as a registering primitive. In order to avoid the situation that the vertical face of a building is extracted as a rod object by mistake to cause the wrong separation of the rod object when the rod object is extracted, the invention firstly uses K-means unsupervised classification for non-ground point clouds, and tags the scene point cloud street trees, the street lamps, the identification plates and the telegraph poles after classification, and the trunks, the street lamps, the identification plates and the telegraph poles of the street trees are extracted as rod objects. In order to ensure the extraction effect of the rod, firstly, a box packing method is used for establishing a space grid for the classified point clouds, and the superfluous point clouds higher than the rod are removed by utilizing an elevation threshold; traversing the grid aerial view by using a random sampling consistency algorithm (RANSAC), setting a cylindrical radius by taking a trunk as a center, and searching an arc-shaped point set for the radius neighborhood; and finally, judging the grid density, the minimum number of support points and the angle threshold value perpendicular to the vertical face, and regarding the grid as rod-shaped point cloud meeting the requirements. The lamp cap part of the urban street lamp has a horizontal extension area, the number of point clouds of the lamp cap part and the signboard part is small, and the rod-shaped point clouds can be accurately extracted by setting different densities and minimum supporting point cloud number thresholds.

According to the coordinate extreme value in the point cloud, an m multiplied by n grid with the side length L is constructed, the point cloud data are divided into a regular two-dimensional grid, and the m and n calculation modes are as follows:

wherein: x is X _max 、X _min 、Y _max And Y _min For the extreme value of the point cloud value, L is the grid width, and the invention is set to be 1 meter according to the setting of the point cloud density.

And searching a fitting circle by adopting a RANSAC algorithm to extract the cylindrical point cloud, and setting the range of a radius buffer area to be 0.05-0.5 meter to extract the cylindrical point cloud.

Wherein:for the coordinates of each point, +.>For the center of each fitting circle, R _i Radius d _i E for the distance to each cylinder _i Representing the degree of difference of the grid blocks from the fitted circle.

Example 6

Extracting edge feature points

As shown in fig. 14, the edge points have strong stability and repeatability, particularly in large scene data, have good effects on feature detection and matching, and have high extraction speed. The invention takes into account that fewer rods affect the quality of point cloud registration, so a Normal Aligned Radial Feature (NARF) algorithm is used to extract edge keypoints. NARF is a 3D feature point detection and description algorithm, and a distance image generated by using point cloud is used for extracting a boundary from a foreground to a background, and key points are proposed by identifying objects from the distance image, so that the operation efficiency is higher than that of directly operating on the point cloud. The extraction process considers the edge and object surface change information, the position of the key point can be stably and repeatedly detected at different visual angles, a stable supporting area is needed, and a descriptor and a unique normal vector can be calculated. In order to meet the above requirements, the following detection steps are proposed to perform the key point extraction:

(1) And traversing each depth image point, and searching a position with depth mutation in a neighboring area for edge detection. For a point p in each range image _i Selecting a square window field point in the same plane, and marking asWhere s is the window width. Calculation of p _i And->The 3D distance (Euclidean distance) of each point in the interior, get +.>For->Arranged in ascending order to obtain ∈ ->Let N be _pi At least M points and p _i On the same surface, select δ=d _M ' as p _i At p _i -distances of non-coplanar points within s-neighborhood. For p _i Four confidence levels are calculated, up, down, left, right, representing the likelihood of an edge being present, the following being taken as an example on the right:

firstly, calculating the average distance of the right point:wherein m is _p To consider the number of right dots, p _x,y Representing points at (x, y) positions in the image. Calculating 3D distance D _right ＝||p _x,y -p _right Sum of s _right ∈[0,1)：

Wherein: s is(s) _right The larger the value, the more marked point p _i The more pronounced the right distance change, the greater the probability that an edge will be present.

For the obtained s _right Smoothing the graph to obtain continuous edges while improving noise robustness, determining p _x,y The type of point is according to p _x,y And p is as follows _right P _x,y A smaller distance indicates that an object boundary is possible, and vice versa, a shadow boundary is obtained, and the window width s corresponding to the distance between the shadow boundary is calculated _right ′＝max(0.9,(1-s _shadow ) ³ )×s _right If s _right ' greater than the threshold and is (p) _x-1,y ,p _x+1,y ) Maximum in the field, p is considered to be _x,y For object boundaries, shadow points and shadow boundaries are found according to the formula, and the shadow points, namely the extracted interest feature points, are located between the shadow points.

(2) The coefficient of measuring the surface variation, and the principal direction of variation, is determined from the surface variation of the neighboring region. For each point p _i The method belongs to a range image, and a normal vector of the range image is estimated by using a 3D position of a point with Euclidean distance within 2 delta in the 2D field and a Principal Component Analysis (PCA) algorithm. And calculating a main direction v and an amplitude lambda, wherein v is an edge direction for boundary points, v is a main curvature direction for other points, and lambda is the amplitude of the direction. For each v, a weight ω∈0, 1) is assigned, for the boundary point ω=1; for other points ω=1- (1- λ) ³ . One parameter sigma is specified to represent the supportable range size. For each point p in the image, find it so that the Euclidean distance is atPoints in which there is no edge in between (two points on the same surface), denoted as N _p ＝{p ₀ ,p ₁ ,…,p _N }。

For a pair ofWill v _pi Projecting onto a plane P, obtaining an angle +.>For->The following adjustments were made:

so that alpha' e (-180 deg., 180 deg.).

(3) To interest value The rows are smoothed and filtered. Smoothing ω and α using bounded gaussian kernels to obtain a point of interest I for point p _p Points exceeding the threshold I are taken as interest points:

I(P)＝I ₁ (P)·I ₂ (P)(14)

wherein: i ₁ (P) high weight point values, i.e. points where the plane changes drastically, can be extracted; i ₂ (P) points in the vicinity where adjacent points exist and the main direction changes drastically can be added.

(4) And carrying out maximum value-free compression to find a final key point, namely the NARF key point.

Example 7

Improved ICP registration

The traditional ICP algorithm adopts the Euclidean distance nearest point as a corresponding point, so that noise interference is easy to occur, a plurality of error point pairs can be introduced, and a solution obtained by direct least square optimization calculation can be influenced to a certain extent; the eigenvalues decomposed by a Singular Value Decomposition (SVD) method are sensitive to noise, and the influence on the accuracy of registration leads to large errors. Thus, the present invention utilizes a probabilistic model to construct plane-to-plane distances to search for corresponding points, and the improved ICP is based on the probabilistic model's additional simplified steps. The rest of the algorithm is kept unchanged, the complexity of calculation is reduced, the registration speed is increased, and the Euclidean distance is still used for calculating the corresponding relation instead of the probability measurement. This is done to allow the KD-tree to be used when finding the closest point and thus to maintain the main advantages of ICP over other full probability techniques, speed and ease of application. The minimization of the probability function can be carried out by nonlinear optimization by using a steepest descent method, a Newton method and the like, but the methods have large calculated amount and slow convergence. To this end, the present invention uses the improved ICP plane-to-plane algorithm of the Levenberg-Marquardt (LM) algorithm for global fine registration. Not only reduces the requirement on the initial position, but also accelerates the convergence speed of registration.

Assume that two sets of point clouds a= { a ₀ ,a ₁ ,…,a _N Sum b= { B ₀ ,b ₁ ,…,b _N Performing nearest neighbor distance lookup according to their corresponding index relationship such thatWhere T is the transform stiffness matrix, E (T) is the error loss function, and d is the Euclidean distance threshold. In the probabilistic model, it is assumed that there are two point sets +.>Andaccording to->And->Yielding a and B which, in this case,and->Is the covariance matrix associated with the measurement points. Assume that T is ^* To correctly transform the rigid matrix, so as to obtain perfect corresponding relation

For any rigid matrix T, defineAnd deriving +.>Is given by a _i And b _i Derived from an independent Gaussian distribution, by the above formula

Iterative calculation of T using maximum likelihood estimation

Simplifying the above steps into

Wherein: and t opt is an updated rigid transformation matrix, and t=r·t is a rigid matrix obtained by LM decomposition.

To improve performance with respect to point-to-plane and increase symmetry of the model, point-to-plane is used to excite the probabilistic model. Essentially, each measurement point provides only constraints along its surface normal. To construct such a model, each sample point is considered to have a high covariance distribution along its local plane and a very low covariance in the direction of the surface normal. If the surface normal of the point is n _i The covariance matrix becomes

Wherein: alpha represents a small constant of covariance along the normal. This corresponds to a very high confidence and knows the position along the normal, but does not determine its position in the plane.

Wherein: s is(s) _i ，t _i Is a as _i And b _i Normal vector of corresponding point, R _x Is a rotation matrix for transforming the basis vectors.

And substituting (20) and (21) into the above formula (18) to calculate the transformation rigid matrix T opt. I.e. the construction of the plane-to-plane ICP algorithm is completed, the schematic diagram is shown in fig. 15.

Fig. 15 provides an illustration of the effect of the algorithm in an extreme case. In this case, all points along the vertical portion of the green scan are incorrectly associated with a single point in the red scan. Because of the non-uniform surface orientation, the plane-to-plane will automatically ignore these matches; each corresponding final summed covariance matrix will be isotropic and will form a very small contribution to the objective function relative to the precisely defined corresponding covariance matrix. Another view of this behavior is as each corresponding soft constraint. The non-uniform matching allows the red scan point to move along the horizontal axis while the green scan point is free to move along the vertical direction. Thus, incorrect correspondence forms a very weak and informationless constraint on the overall alignment.

To reduce the effects of the point cloud pairing process noise and outliers, the Cauchy loss function is utilized as a weighting function to resist noise and outliers. The modified ICP algorithm is as follows:

wherein: ρ _Cauchy For Cauchy loss function weights, k is the loss parameter, and typically k is 1.345 and r is the residual at the time of optimization.

R obtained by decomposing LM is substituted into the above formula to obtainAnd->Then substituting to obtain updated transformation matrix T opt and obtaining +.>By this iteration +.>Less than the set Euclidean distance threshold τ or the transform update matrix is made less than the set thresholdSetting the iteration times to 200 times, and stopping iteration when one of the conditions is met.

The euclidean distance RMSE of the point pairs is used as the registration error accuracy for evaluation. RMSE obtainable from equation (22) is as follows:

wherein: s is(s) _i ＝(s _ix ,s _iy ,s _iz ) Coordinates of cloud points of source points, t _i ＝(t _ix ,t _iy ,t _iz ) And n is the number of registered point pairs for the coordinates of the point cloud to be registered.

The preferred embodiments of the invention disclosed above are intended only to assist in the explanation of the invention. The preferred embodiments are not exhaustive or to limit the invention to the precise form disclosed. Obviously, many modifications and variations are possible in light of the above teaching. The embodiments were chosen and described in order to best explain the principles of the invention and the practical application, to thereby enable others skilled in the art to best understand and utilize the invention. The invention is limited only by the claims and the full scope and equivalents thereof.

Claims (7)

1. The automatic registration method of the fast vehicle-mounted laser point cloud in the complex urban environment is characterized by comprising the following steps of:

s1: in-vehicle laser scanning typically generates point cloud datasets with different point densities, performs point cloud filtering and voxel downsampling;

s2: filtering a cloth model, and setting a distance threshold value to quickly separate ground points;

s3: extracting the rod-shaped object as a registration element;

s4: extracting edge key points by using a Normal Aligned Radial Feature (NARF) algorithm;

s5: global fine registration was performed using ICP plane-to-plane algorithm modified by the face-to-face and Levenberg-Marquardt (LM) algorithm.

2. The method for automatically registering the fast vehicle-mounted laser point cloud in the complex urban environment according to claim 1, wherein the method comprises the following steps: the point cloud filtering and voxel downsampling specifically comprises the following steps: vehicle-mounted laser scanning generally generates point cloud data sets with different point densities, and generates outliers, a part of which cannot be used; the statistical filter is mainly used for filtering outliers in the point cloud, calculating the average distance from each point to k nearest points, and assuming that the obtained result is a Gaussian distribution, the shape of the Gaussian distribution is determined by the mean value and the standard deviation, and then the points with the average distance being out of the standard range can be defined as outliers and removed from the data; the probability density function of the Gaussian random variable n is

Wherein: u represents the distance from which the object is to be moved,mean value of all u, sigma standard deviation;

distance threshold d _max Can be expressed as

d _max ＝u+α·σ (2)

α is a constant, which may be referred to as a scaling factor, which depends on the number of neighbor points; finally, traversing the point cloud again, and eliminating that the average distance between the point cloud and k neighbor points is greater than d _max Is a point of (2);

voxel size calculation is the key point of voxel sampling, and a reasonable voxel size can effectively simplify the original point cloud; firstly, determining a proper L value, wherein the size of the L value can directly influence the downsampling rate of the point cloud, and when L is overlarge, the points falling into voxels become more, the sampling rate is low, so that the local characteristics of the point cloud are lost; when the L value is small, the number of points in the voxel is small, the sampling rate is high, and although the local characteristics can be well reserved, the sampling result is not obvious; to solve this problem, the combination of the downsampling rate ρ and the point cloud density estimates the voxel side length L, i.e

Wherein: l' represents the estimated voxel side length; s is the average point cloud number within a unit voxel; gamma is the point cloud density; s and gamma are expressed as

Wherein: n is the number of point cloud points; substituting the formulas (3) and (4) into the formula (2) to obtain

Obtaining a proper voxel side length L to divide a bounding box of the point cloud P, so as to better preserve the feature of the point cloud and calculate the inside of the voxel The center of gravity of all points is calculated by the calculated center of gravity point P _k ' as downsampling points; p (P) _k ' is denoted as

Wherein:to calculate the number of k-th voxel inliers, k ε [1,2, …, m]The method comprises the steps of carrying out a first treatment on the surface of the A new point cloud P ', P' can be obtained from equation (6), which is the final downsampled point cloud.

3. The method for automatically registering the fast vehicle-mounted laser point cloud in the complex urban environment according to claim 1, wherein the method comprises the following steps: the step S2 specifically comprises the following steps: inverting the point cloud, and covering the inverted curved surface with cloth from above; and determining the final shape of the cloth model by analyzing the interaction between the nodes and the corresponding points of the cloth model, and dividing the ground and the non-ground according to the shape of the final shape.

4. The method for automatically registering the fast vehicle-mounted laser point cloud in the complex urban environment according to claim 1, wherein the method comprises the following steps: the extraction shaft specifically comprises: firstly, performing non-supervision classification on non-ground point clouds by using a K mean value, marking scene point cloud street trees, street lamps, identification plates and telegraph poles after classification, and extracting trunks, street lamps, identification plates and telegraph poles of the street trees as rod-shaped targets; in order to ensure the extraction effect of the rod, firstly, a box packing method is used for establishing a space grid for the classified point clouds, and the superfluous point clouds higher than the rod are removed by utilizing an elevation threshold; traversing the grid aerial view by using a random sampling consistency algorithm (RANSAC), setting a cylindrical radius by taking a trunk as a center, and searching an arc-shaped point set for the radius neighborhood; finally, judging the grid density, the minimum number of support points and the angle threshold value vertical to the vertical face, and regarding the grid as rod-shaped point cloud meeting the requirements; the lamp cap part of the urban street lamp has a horizontal extension area, the number of point clouds of the lamp cap part and the signboard part is small, and the rod-shaped point clouds can be accurately extracted by setting different densities and minimum supporting point cloud number thresholds;

According to the coordinate extreme value in the point cloud, an m multiplied by n grid with the side length L is constructed, the point cloud data are divided into a regular two-dimensional grid, and the m and n calculation modes are as follows:

wherein: x is X _max 、X _min 、Y _max And Y _min The value extremum of the point cloud is set to be 1 meter according to the setting of the point cloud density, wherein L is the grid width;

searching a fitting circle by adopting a RANSAC algorithm to extract a cylindrical point cloud, and setting the range of a radius buffer area to be 0.05-0.5 m to extract the cylindrical point cloud;

wherein:for the coordinates of each point, +.>For the center of each fitting circle, R _i Radius d _i E for the distance to each cylinder _i Representing the degree of difference of the grid blocks from the fitted circle.

5. The method for automatically registering the fast vehicle-mounted laser point cloud in the complex urban environment according to claim 1, wherein the method comprises the following steps: the extracting edge key points specifically comprises: extracting edge key points by using a Normal Aligned Radial Feature (NARF) algorithm; NARF is a 3D characteristic point detection and description algorithm, the distance image generated by using the point cloud is used for extracting the boundary from the foreground to the background, and the key point is provided by identifying an object from the distance image, so that the operation efficiency is higher than that of directly operating on the point cloud; the extraction process considers the edge and object surface change information, the position of the key point can be stably and repeatedly detected at different visual angles, a stable supporting area is required, and a descriptor and a unique normal vector can be calculated;

The edge key point extraction method specifically comprises the following substeps:

s4.1: traversing each depth image point, and searching a position with depth mutation in a neighboring area for edge detection; for a point p in each range image _i Selecting a square window field point in the same plane, and marking asWherein s is the window width; calculating pi and +.>The 3D distance (Euclidean distance) of each point in the interior, get +.>For->Arranged in ascending order to obtain ∈ ->Let->At least M points are located on the same surface as pi, and delta=d is selected _M ' as p _i At p _i -distances of non-coplanar points within s-neighborhood; for p _i Four confidence levels are calculated, up, down, left, right, representing the likelihood of an edge being present, the following being taken as an example on the right:

firstly, calculating the average distance of the right point:wherein m is _p To consider the number of right dots, p _x,y Representing a point at the (x, y) position in the image; calculating 3D distance D _right ＝||p _x,y -p _right Sum of s _right ∈[0,1)：

Wherein: s is(s) _right The larger the value, the more marked point p _i The more obvious the right distance change, the greater the probability of an edge being present;

for the obtained s _right Smoothing the graph to obtain continuous edges while improving noise robustness, determining p _x,y The type of point is according to p _x,y And p is as follows _right P _x,y A smaller distance indicates that an object boundary is possible, and vice versa, a shadow boundary is obtained, and the window width s corresponding to the distance between the shadow boundary is calculated _right ′＝max(0.9,(1-s _shadow ) ³ )×s _right If s _right ' greater than the threshold and is (p) _x-1,y ,p _x+1,y ) Maximum in the field, p is considered to be _x,y Finding shadow points and shadow boundaries for object boundaries according to a formula, wherein the shadow points are occlusion points, namely the extracted interest feature points, located between the shadow points and the shadow boundaries;

s4.2: determining a coefficient for measuring the surface change and a main direction of the change according to the surface change of the adjacent area; for each point p _i The method comprises the steps of (1) estimating a normal vector of a distance image by utilizing a 3D position of a point with Euclidean distance within 2 delta in a 2D field and a Principal Component Analysis (PCA) algorithm; calculating a main direction v and an amplitude lambda, wherein v is an edge direction for boundary points, v is a main curvature direction for other points, and lambda is the amplitude of the direction; for each v, a weight ω∈0, 1) is assigned, for the boundary point ω=1; for other points ω=1- (1- λ) ³ The method comprises the steps of carrying out a first treatment on the surface of the Specifying a parameter sigma representing the supportable range size; for each point p in the image, find it so that the Euclidean distance is atPoints in which there is no edge in between (two points on the same surface), denoted as N _p ＝{p ₀ ,p ₁ ,…,p _N }；

For a pair ofWill->Projecting onto a plane P, obtaining an angle +.>For->The following adjustments were made:

so that alpha' is E (-180 degrees, 180 degrees);

s4.3: smoothing and filtering the interest value; smoothing ω and α using bounded gaussian kernels to obtain a point of interest I for point p _p Points exceeding the threshold I are taken as interest points:

I(P)＝I ₁ (P)·I ₂ (P) (14)

wherein: i ₁ (P) high weight point values, i.e. points where the plane changes drastically, can be extracted; i ₂ (P) points in the vicinity where there are adjacent points and the main direction changes drastically can be added;

s4.4: and carrying out maximum value-free compression to find a final key point, namely the NARF key point.

6. The method for automatically registering the fast vehicle-mounted laser point cloud in the complex urban environment according to claim 1, wherein the method comprises the following steps: the step S5 specifically includes:

assume that two sets of point clouds a= { a ₀ ,a ₁ ,…,a _N Sum b= { B ₀ ,b ₁ ,…,b _N Performing nearest neighbor distance lookup according to their corresponding index relationship such thatThe method is established, wherein T is a transformation rigid matrix, E (T) is an error loss function, and d is a Euclidean distance threshold; in the probabilistic model, it is assumed that there are two point sets +.>Andaccording to->And->Yielding a and B which, in this case,and->Is associated withMeasuring covariance matrix associated with the points; assume that T is ^* To correctly transform the rigid matrix, so as to obtain perfect corresponding relation

For any rigid matrix T, defineAnd deriving +.>Is given by a _i And b _i Derived from an independent Gaussian distribution, by the above formula

Iterative calculation of T using maximum likelihood estimation

Simplifying the above steps into

Wherein: t opt is an updated rigid transformation matrix, and t=r·t is a rigid matrix obtained by LM decomposition;

to improve performance with respect to point-to-plane and increase symmetry of the model, a point-to-plane excitation probability model is used; essentially, each measurement point provides only constraints along its surface normal; to construct such a model, each sample point is considered to have a high covariance distribution along its local plane and to have non-normal to the surfaceA very low covariance; if the surface normal of the point is n _i The covariance matrix becomes

Wherein: alpha represents a small constant of covariance along the normal; this corresponds to a very high confidence and knows the position along the normal, but does not determine its position in the plane;

wherein: s is(s) _i ，t _i Is a as _i And b _i Normal vector of corresponding point, R _x A rotation matrix for transforming the basis vectors;

substituting (20) and (21) into the formula (18) to calculate a transformation rigid matrix T opt; namely, the construction of the ICP algorithm from plane to plane is completed;

to reduce the influence of noise and outliers in the point cloud pairing process, the Cauchy loss function is used as a weighting function to resist the noise and the outliers; the modified ICP algorithm is as follows:

Wherein: ρ _Cauchy For Cauchy loss function weight, k is a loss parameter, and under normal conditions, k is 1.345, and r is a residual error during optimization;

dividing LM intoR obtained by the solution is substituted into the above formula to obtainAnd->Then substituting to obtain updated transformation matrix T opt and obtainBy this iteration +.>Less than the set Euclidean distance threshold τ or the transform update matrix is made less than the set threshold +.>Setting the iteration times to 200 times, and stopping iteration when one of the conditions is met;

estimating by using the Euclidean distance RMSE of the point pairs as registration error precision; RMSE obtainable from equation (22) is as follows:

wherein: s is(s) _i ＝(s _ix ,s _iy ,s _iz ) Coordinates of cloud points of source points, t _i ＝(t _ix ,t _iy ,t _iz ) And n is the number of registered point pairs for the coordinates of the point cloud to be registered.

7. The method for automatically registering the fast vehicle-mounted laser point cloud in the complex urban environment according to claim 6, wherein the method comprises the following steps: the step S5 further includes: the method comprises the steps of using an improved Iterative Closest Point (ICP) algorithm and normal vector constraint, using rod objects and edge feature points as registration primitives, using a sampling geometric consistency algorithm to remove error point pairs, and using a two-way KD tree to quickly correspond to the relation of the feature points, so that the registration speed is increased and the precision is improved.