Point cloud registration method based on neighborhood grid gaussian weighted top view

Abstract

The application discloses a point cloud registration method based on a neighborhood grid Gaussian weighted overhead view, which comprises four steps of point cloud overhead view generation, image matching, point cloud registration and global optimization, and the point cloud is registered by projecting point cloud data along a z axis into an image and obtaining corresponding feature matching on the image. The improved steps in the registration mainly obtain better images and correct matching of endpoints through neighborhood weighting. The neighborhood Gaussian weighted point cloud projection is that after the point cloud is rasterized, the Gaussian weighted density of points in the neighborhood of each grid to the center of the grid is calculated, and finally normalized to the pixel range of a gray image 0-255. The image matching is feature searching and matching on two gray images, and the registration efficiency is improved through dimension reduction.

Description

Technical Field

[0001] This invention belongs to the field of 3D reconstruction and computer graphics, specifically relating to a method for rapid registration of point clouds based on Gaussian weighted top view. Background Technology

[0002] Ground-based lidar-based point cloud acquisition provides an effective means of obtaining intuitive information about the three-dimensional world. Currently, the need and application of 3D description of the real world and the digitization of scenes are becoming increasingly widespread. For example, with the rapid development of sensor technology and the emergence of numerous software processing tools, point cloud data has become widely available. Point cloud data, through a large set of 3D points, digitizes scenes, densely and accurately representing the 3D geometry of objects in the environment.

[0003] In real life, the following two scenarios often occur:

[0004] Scenario 1: When designing the interior of a building, it is necessary to obtain a model of the building and its precise dimensions in advance. Since manual measurement is time-consuming and costly, it is essential to obtain a high-precision 3D building model quickly.

[0005] Scenario 2: With the development of robotics and autonomous driving technology, autonomous vehicles and other vehicles require greater environmental perception capabilities. High-precision maps and rapid environmental perception are needed when vehicles are driving autonomously.

[0006] To quickly obtain high-precision 3D building models and outdoor 3D maps, researchers have proposed many effective methods. Most current mainstream point cloud registration methods are based on 3D spatial normal vectors, geometric relationships, or density transformations to obtain corresponding matching points or point sets. Then, they iteratively optimize the obtained matching sets to obtain the optimal transformation. The advantages of existing methods are: feature points or point sets can represent the overall point cloud information with less data; registration efficiency is improved by finding matching feature points or point sets; and feature points are obtained by performing keypoint detection on the original point cloud and designing a robust feature descriptor. The corresponding point pairs are then found in the feature point cloud to calculate the transformation matrix, thus completing the registration.

[0007] Existing methods still have shortcomings, regardless of whether they are based on geometric features, normal vector features, or other features in 3D space, there is room for improvement. Some point cloud registration algorithms have proposed corresponding improvement strategies to address issues such as efficiency, rotation invariance, and noise, but these usually lead to increased computational costs. Some researchers have also combined point cloud registration with image matching, proposing image-based point cloud registration methods, which provides a new direction for the registration of terrestrial laser point cloud data with low overlap rates and large data volumes. However, with the rapid development of new technologies such as laser SLAM, robot navigation, and augmented reality, as well as the social demand for digital 3D reconstruction of smart cities, higher requirements have been placed on the accuracy and efficiency of point cloud registration. Therefore, it is particularly necessary to develop a point cloud registration method that meets the requirements in terms of accuracy, efficiency, and computational cost. Summary of the Invention

[0008] The purpose of this invention is to provide a point cloud registration method based on a neighborhood grid Gaussian weighted top view, in order to solve the problem that existing point cloud registration methods mentioned in the background art cannot simultaneously achieve the desired registration accuracy, efficiency, and computational cost. The method finds feature points in the top view and obtains the transformation of the three-dimensional point cloud through image feature matching, thereby improving registration efficiency by using the idea of ​​dimensionality reduction.

[0009] To achieve the above objectives, the present invention provides the following technical solution: a point cloud registration method based on a neighborhood raster Gaussian weighted top view, comprising the following steps:

[0010] Step 1: Generating the point cloud top view; specifically including the following steps:

[0011] A) Traverse the point cloud, find the minimum point value (minX, minY, minZ) and maximum point value (maxX, maxY, maxZ) in the X, Y, and Z directions, and create [(maxX-minX) / VoxelSize]×[(maxY-minY) / VoxelSize]×[(maxZ-minZ) / VoxelSize] grids given the grid size VoxelSize;

[0012] B) Traverse the grids on the x and y planes. That is, all grids with the same x and y but different z are considered as one grid in the top view. Take the weighted density of the center point of each grid as the corresponding pixel value in the image.

[0013] Step 2: Image Matching; After generating the grayscale point cloud image, feature point search, feature description, and feature matching are performed on the image; specifically, the following steps are included:

[0014] C) First, SIFT feature detection, description and matching in OpenCV are used to obtain corresponding points. Then, Lowe's algorithm is used to select excellent matching corresponding points. Finally, the RANSAC framework and least squares method are used to obtain image transformation.

[0015] D) Correct matches need to be selected; Hough transform is used to detect line segments on the image, and then vectors are used to determine line segments with similar angles: approximately parallel line segments are divided into one category, and then the line segments of each category are sorted according to their quantity. Line segments with fewer quantity are definitely false detections; for parallel line segments, the distance between the line segments is calculated, and the endpoints of the closest lines are selected for line segment fitting to solve the problem of line segment overlap; based on the detected line segments, their endpoints are obtained, and their intersections are used to obtain the corner points; red endpoints are detected using I_target, and green endpoints are detected using I_source. The results calculated in step B) are applied to the endpoints of I_source, and then the distance dist between the red endpoint and the transformed green endpoint is calculated. If dist < two pixels, it is considered a matching endpoint; if at least two endpoints are found to be matched, it is considered a correct transformation;

[0016] Step 3: Point Cloud Registration; The point cloud is projected onto the image and scaled. The image transformation is converted to a point cloud transformation. A connection between the point cloud and the image needs to be established. The point cloud transformation is obtained by combining α, β, and displacement parameters. Let P1:(x1,y1,z1) and P2:(x2,y2,z2) be a pair of corresponding points in the two point clouds, satisfying the following relationship:

[0017] P1 = R_cloud·P2 + T_cloud; the corresponding points of P1 and P2 on the projected grayscale image are I1:(u1,v1) and I2:(u2,v2), respectively, and they satisfy the following relationship:

[0018] I1 = R_Image·I2 + T_Image;

[0019] ((u1@v1))=((r11&r12@r21&r22))·((u2@v2))+((t1@t2));

[0020] R_Image and T_Image can be determined from the calculation results in step two;

[0021] Based on the coordinate relationship between the point cloud and the image, I can obtain their scaling relationship and displacement. Let the scaling on X be α and the displacement be tx, and the scaling on Y be β and the displacement be ty. The 3×3 rotation and 3×1 translation transformation matrices in three dimensions can be written as 2×2 and 2×1 planar transformation matrices; the transformation formula from P1 to I1 is:

[0022] ((u1@v1))=((α1&0@0&β1))((x1@y1))+((tx1@ty1));

[0023] The transformation formula from P2 to I2 is:

[0024] ((u2@v2))=((α2&0@0&β2))((x2@y2))+((tx2@ty2));

[0025] Substitute the two equations into the formula.

[0026] From ((u1@v1))=((r11&r12@r21&r22))·((u2@v2))+((t1@t2)), we get:

[0027] ((α1x1+tx1@β1y1+ty1))=((r11&r12@r21&r22))((α2x2+tx2@β2y2+ty2))+((t1@t2))

[0028] Simplifying, we get:

[0029] ((x1@y1))=((α2r11 / α1&β2r12 / α1@α2r21 / β1&β2r22 / β1))((x2@y2))+(((tx2r11+ty2r12+t1-tx1) / α1@(tx2r21+ty2r22+t2-ty1) / β1))

[0030] Therefore, the obtained image transformation is converted into a point cloud transformation;

[0031] Step 4: Global Optimization; Global registration optimization is performed through pose optimization to reduce overall error; this includes the following steps:

[0032] E) First, the transformation relationship between each pair of multi-site clouds is obtained through steps one to three;

[0033] F) Obtain the transformation relationship between each pair of multi-site clouds through steps one to three;

[0034] G) Construct a graph optimization problem for all the point clouds that can be registered. The vertices of the graph are the poses of the point cloud at each station, and the edges are the transformation matrices between any two point clouds that can be registered.

[0035] H) The error is minimized through overall optimization, and finally the ICP point cloud fine registration algorithm is used to further reduce the error.

[0036] Preferably, step C) specifically includes the following steps:

[0037] (1) For the source image I_target and the target image I_source, the SIFT algorithm is used to extract n pairs of corresponding points. If n is less than 2, a match cannot be found and the process is terminated directly.

[0038] (2) In n pairs of points, I take only two sets of corresponding points each time. First, I calculate the rotation angle θ through these two sets of corresponding points, and then rotate I_source by θ to get I_(source-θ);

[0039] (3) Then, calculate the distance Δd between I_target and I_(source-θ) through the midpoint of the two sets of corresponding points, and translate I_(source-θ) by Δd to obtain I_(source-θ-d);

[0040] (4) Using the angles and displacements obtained from steps (2) and (3), calculate the transformation T_(i,j) according to the following formula.

[0041] T_(i,j)=((cosθ&-sinθ&Δd_x@sinθ&cosθ&Δd_y))

[0042] (5) Apply the transformation T_(i,j) to the n feature points of I_source, calculate the distance between these n points and their corresponding n points of I_target, and the points with a distance less than the threshold are interior points. Count the number of interior points m.

[0043] (6) Repeat steps (2)-(5) to obtain the transformation with the most interior points m_max, which is considered the best transformation. Then use these m_max interior points together to calculate the final image transformation T_final.

[0044] Preferably, in step four, the Levenberg-Marquadt method from the g2o library is used for global registration optimization.

[0045] Compared with the prior art, the beneficial effects of the present invention are:

[0046] This method registers point clouds by projecting them along the z-axis into an image and obtaining corresponding feature matches on the image. Point cloud projection involves: point cloud thinning, rasterization, and neighborhood Gaussian weighting. Image matching involves: feature matching, filtering for better matches based on RANSAC principles, and obtaining correct matches using endpoints. Point cloud registration involves: establishing the scaling relationship between the point cloud and the image, pairwise point cloud registration, and multi-point point cloud registration. The improved registration steps primarily involve obtaining better images through neighborhood weighting and obtaining correct matches using endpoints. Neighborhood Gaussian weighted point cloud projection involves calculating the Gaussian weighted density from each grid cell and its neighborhood points to the grid center after rasterization, and finally normalizing it to the 0-255 pixel range of a grayscale image. Image matching involves feature lookup and matching on two grayscale images, using stable endpoints such as corners or the ends of lines to filter image transformations, and improving registration efficiency through dimensionality reduction. Attached Figure Description

[0047] Figure 1 The overall registration method of this invention uses multi-station input point cloud maps;

[0048] Figure 2 This is a top-down view of the cloud generated for every two sites in the overall registration method of this invention;

[0049] Figure 3 This is a feature matching diagram of the overall registration method of the present invention;

[0050] Figure 4 It is the image matching result;

[0051] Figure 5 This is a schematic diagram of point cloud registration between two sites obtained through point cloud changes;

[0052] Figure 6 It is a topology graph obtained by registering every two sites in a multi-site cloud.

[0053] Figure 7 This is a schematic diagram for global optimization;

[0054] Figure 8 To output the final multi-station registration results;

[0055] Figure 9 It is the original point cloud that generates a grayscale image based on neighborhood weighted density from point cloud data;

[0056] Figure 10 This is a schematic diagram of obtaining the neighborhood weighted density;

[0057] Figure 11 This is a top-down view of the generated point cloud;

[0058] Figure 12It is the process of rasterizing the point cloud after rotation, based on the neighborhood-weighted feature point rotation invariance analysis.

[0059] Figure 13 This is a grayscale image generated after neighborhood weighting;

[0060] Figure 14 This is an unweighted grayscale image. Unweighted images detect more feature points after point cloud rotation, exhibiting poor rotation invariance. Weighted images show more stable feature points.

[0061] Figure 15 It is the process of obtaining the image transformation matrix by using two sets of corresponding feature points to obtain the rotation angle and displacement transformation.

[0062] Figure 16 This is a top-down view generated from a point cloud by obtaining image endpoints.

[0063] Figure 17 To detect a straight line using the Hough transform;

[0064] Figure 18 To detect erroneous lines and overlapping lines;

[0065] Figure 19 The correct line is selected after analyzing the slope and the distance between the lines.

[0066] Figure 20 After detecting the endpoints, the endpoint verification image is transformed. This involves transforming the images containing the detected endpoints from two different images.

[0067] Figure 21 These are endpoints that are close to each other after image transformation.

[0068] Figure 22 It represents the transformation relationship between the point cloud coordinate system and the image coordinate system.

[0069] Figure 23 It is a multi-site point cloud registration, which is the transformation relationship between pairs of point clouds. The green lines represent point clouds that can be registered, and the red lines represent point clouds that cannot be registered.

[0070] Figure 24 This paper compares various point cloud projection methods, including methods for generating top views based on 1×1, 3×3, and 5×5 neighborhood weighting, as well as methods for directly calculating point density using binary maps and unweighted methods. In my method, point cloud registration uses 3×3 neighborhood weighting.

[0071] Figure 25 A comparison of the time and accuracy of pairwise registration of point clouds, including 4PCS, super-4PCS, wk-4PCS, and our own registration method.

[0072] Figure 26 The experimental results across multiple datasets are presented as an unstitched point cloud top view.

[0073] Figure 27 To stitch together the top view of the point cloud;

[0074] Figure 28 For a spliced ​​multi-site cloud. Detailed Implementation

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

[0076] Please see Figures 1 to 8 The present invention provides a point cloud registration method based on a domain-weighted top view, which includes four steps: point cloud top view generation, image matching, point cloud registration, and global optimization.

[0077] The collected point cloud data were all acquired from a LiDAR equipped with an automatic leveling device, meaning the z-axis of the point clouds was parallel. These point clouds were projected along the z-axis onto the xy-plane to generate a grayscale image normalized based on the point cloud density. Here, a neighborhood Gaussian weighting method was used to smooth the image. Next, feature matching was performed. Due to the presence of many incorrect corresponding points in the features, endpoint verification was used to obtain the correct matching image. After obtaining the image transformation matrix, the transformation matrix of the point cloud was obtained through the scaling and translation relationship between the point cloud and the projected image, thus registering the point clouds. After performing the above coarse registration steps between multiple point clouds, the poses of the point clouds were placed in a registration relationship graph. Nodes represent point clouds, and edges represent the transformation relationships between point clouds. Transformation optimization was performed to reduce the overall error. Finally, the ICP algorithm was used for further refinement to obtain the complete point cloud.

[0078] Step 1: Generate a top-down view of the point cloud density based on neighborhood Gaussian weighted average. Specifically, this includes:

[0079] Step 1.1: Point cloud data is often very large. If each point in the point cloud is directly projected into pixel coordinates, i.e., the z-value of each point is set to 0, the computational load is very high. Therefore, a point cloud rasterization method is first used to make the point cloud projection more efficient. Traverse the point cloud, find the minimum point values ​​(minX, minY, minZ) and maximum values ​​(maxX, maxY, maxZ) in the X, Y, and Z directions, and give a raster size VoxelSize to establish... A grid.

[0080] Step 1.2: Traverse the grids on the x and y planes. All grids with the same x and y coordinates but different z coordinates are considered a single grid in the top view. Take the weighted density of the center point of each grid as the corresponding pixel value in the image. For each grid center coordinate (x... c ,y c Consider all points (x) in its neighborhood raster. i ,y i Here we use the 8-neighborhood as an example. Figure 2 The density value is calculated as shown in formula (1):

[0081]

[0082]

[0083] Where σ is the standard deviation in the two-dimensional Gaussian function. The larger the value of σ, the larger the Gaussian weighting range and the smoother the image. In order to smooth the image while better preserving the point cloud features, I only consider the influence of points in the neighborhood and take the value of σ as the length of the raster / 2. It is the distance from each point in the grid to the center point. Here, we only consider the distances on the x and y planes as shown in formula (2).

[0084] After obtaining the density value of each grid cell, I need to normalize all the values ​​to a range of 0 to 255 so that a grayscale image can be generated. The normalization formula is shown in formula (3):

[0085]

[0086] Where g max and g min These are the maximum and minimum density values, respectively, where g(x,y) is the density value for each grid cell, and I... G (x,y) represents the grayscale value of the corresponding pixel grid in the image.

[0087] Finally, we obtained the image from Figures 9 to 11 The effect is that the dot cloud on the wall is denser and has a higher gray value, while the dots on the ceiling are sparser and lighter in color compared to the wall.

[0088] For the theory of weighted averages in the point cloud domain, please refer to... Figures 15 to 17 Rotating a point cloud representing a straight line by an angle, due to the consideration of neighborhood weighting, results in a point cloud that remains a straight line after rotation. Both detected feature points are endpoints of this straight line, demonstrating good rotation invariance. In contrast, the image generated without considering neighborhood transforms a single straight line into two after rotation, leading to the detection of incorrect feature points. These extra feature points are unstable and interfere with subsequent image matching.

[0089] Step 2: After generating the grayscale point cloud image, feature point lookup, feature description, and feature matching are performed on the image. Specifically, this includes:

[0090] Step 2.1: First, feature detection, description, and matching in OpenCV are used to obtain corresponding points. Then, Lowe's algorithm is used to select excellent matching points. Finally, the RANSAC framework and least squares method are used to obtain the image transformation. The specific process is as follows:

[0091] Step 1: For the source image I target and target image I source The algorithm extracts n pairs of corresponding points. If n is less than 2, no match can be found and the process exits directly.

[0092] Step 2: In n pairs of points, I only take two sets of corresponding points each time. First, I calculate the rotation angle θ using these two sets of corresponding points, such as... Figure 15 Show, and will I source Rotate θ to obtain I source-θ .

[0093] Step 3: Then calculate I using the midpoint of the two sets of corresponding points. target and I source-θ The distance Δd between them, and I source-θ Translate Δd to obtain I source-θ-d .

[0094] Step 4: Calculate the transformation T from the angles and displacements obtained in Step 2 and Step 3. i,j As shown in formula (4).

[0095]

[0096] Step 5: Transform T i,j Using I source For n feature points, calculate the sum of these n points and their corresponding I. target The distance between n corresponding points is considered as 'interior points' if the distance is less than a certain threshold. The number of 'interior points' is counted, m.

[0097] Step 6: Repeat Steps 2-5 to obtain the maximum number of 'interior points', m. max The transformation is considered the optimal transformation, and then this m max The in-place points are used together to calculate the final image transformation T. final .

[0098] Step 2.2: The image transformation result may not be correct; correct matches need to be filtered. Hough transform is used to detect line segments on the image. Hough transform is a line segment detection method based on a voting strategy. Next, vectors are used to determine line segments with similar angles. Approximately parallel line segments are grouped into one category, and then the line segments in each category are sorted according to their quantity. Line segments with fewer counts are considered incorrect detections, such as... Figure 20 The circled section on the right. For parallel line segments, I calculated the distances between them, selected the closest endpoints of the lines, and performed line segment fitting to resolve the overlapping issue, such as... Figure 20 The circled part on the left.

[0099] Based on the detected line segments, their endpoints are obtained, and their intersections are used to locate the corner points. For example... Figure 22 I target The red endpoint was detected, I source Detect the green endpoints and use the T calculated in step 1.2. final Application to I source The red endpoint is then used to calculate the distance *dist* between it and the transformed green endpoint. If *dist* is smaller, it is considered a matching endpoint; typically, *dist* < two pixels. If at least two endpoints are found to be matching, the transformation is considered correct.

[0100] Step 3: Projecting a point cloud onto an image involves scaling. While the top-down view of the point cloud doesn't consider the z-axis, it does have transformation scales α and β, as well as displacement, on the x and y axes. The image transformation needs to be converted to a point cloud transformation. Establishing the connection between the point cloud and the image requires combining the α, β, and displacement parameters. This yields the point cloud transformation.

[0101] Let P1:(x1,y1,z1) and P2:(x2,y2,z2) be a pair of corresponding points in two site clouds, and let them satisfy the following relationship:

[0102] P1 = R cloud ·P2+T cloud (5)

[0103] Where R cloud T cloud This is the result I need.

[0104] The corresponding points of P1 and P2 on the projected grayscale image are I1:(u1,v1) and I2:(u2,v2), respectively, and they satisfy the following relationship:

[0105] I1 = R Image ·I2+T Image (6)

[0106]

[0107] Where R Image and T Image T can be obtained final To determine.

[0108] Based on the coordinate relationship between the point cloud and the image, I can obtain their scaling relationship and displacement. I assume that the scaling on X is α and the displacement is tx, and the scaling on Y is β and the displacement is ty. Since the transformation on the z-axis is not considered, the scaling and displacement on the z-axis are not used. Therefore, the 3×3 rotation and 3×1 translation transformation matrices in three dimensions can be written as 2×2 and 2×1 planar transformation matrices. The transformations from P1 to I1 and from P2 to I2 are as shown in formulas (8) and (9) respectively:

[0109]

[0110]

[0111] Substituting formulas (8) and (9) into formula (7) yields:

[0112]

[0113] Simplifying, we get:

[0114]

[0115] The resulting image transformation is then converted into a point cloud transformation.

[0116] Step 4: The main idea of ​​multi-site cloud registration is to transform all point clouds so that they can be unified into a single coordinate system. Since the registration error between any two point clouds is different, this error will inevitably accumulate in the next registration, resulting in a larger error in the multi-site cloud registration.

[0117] To reduce the overall error, global registration optimization is performed using pose optimization. First, the transformation relationships between each pair of multi-site clouds are obtained through the steps described above, such as... Figure 25 As shown, Accept indicates accepted registration, and Reject indicates registration failure. To speed up registration, all point clouds are first projected into binary images, and then all images are read into memory, instead of reading both point clouds each time before projection. Then, the point clouds that can be registered are added to the pose graph for optimization. In the graph structure, nodes represent the pose information of the point clouds, and edges represent the transformation relationships between pairs of points. Let there be two point clouds, i and j, with the pose of point cloud i being T. i The pose of point cloud j is T j Then the transformation T from i to j i,j It can be written in the following form:

[0118] T i,j =T i -1 T j (12)

[0119] We can move the left side of formula (12) to the right side to obtain an error value, which is the deviation between the observed value and the actual value:

[0120]

[0121] Sum up the errors e of all the point clouds that can be registered:

[0122]

[0123] The goal is to minimize equation (14). This can be solved using gradient descent or other methods for finding extrema. I added a trust region to the gradient descent method to control its effectiveness, which can avoid the singularity problem of the coefficient matrix of the thread equation to a certain extent and provide a more stable incremental value.

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

Claims

1. A point cloud registration method based on a neighborhood grid Gaussian-weighted top view, characterized in that, Comprising the following steps: Step one: generation of point cloud overhead view; specifically comprising the following steps: A) traversing the point cloud to find the minimum point values (minX, minY, minZ) and maximum values (maxX, maxY, maxZ) in the X, Y, Z directions of the point cloud, and given a grid size VoxelSize, establish [(maxX-minX) / VoxleSize]x[(maxY-minY) / VoxleSize]x[(maxZ-minZ) / VoxleSize] grids; B) Traverse the grid on the x,y plane, i.e. all the grids with the same x,y but different z are considered as one grid in the top view, take the Gaussian weighted density of the center point of each grid as the pixel value corresponding to the image; for each grid center coordinate (x c ,y c ), consider all the points (x i ,y i ) in its neighborhood grid, calculate its density value as follows: where σ is the standard deviation in a two-dimensional Gaussian function, and the value of σ is taken as the length of the grid / 2; is the distance of each point in the grid to the center point; Step two: image matching; after generating the point cloud grayscale image, feature point searching, feature description and feature matching are performed on the image; specifically comprising the following steps: C) first, SIFT feature detection, description and matching in openCV are used to obtain the same name points, then Lowe's algorithm is used to screen out the optimal matching same name points, and then RANSAC framework and least squares method are used to obtain image transformation; D) Selecting the correct match; using Hough transform on the image to detect lines, from the detected lines get their endpoints, and use their intersections to get the corner points; using I target detecting the endpoints of the red color, I source detecting the endpoints of the green color, applying the results calculated in step C) to the endpoints of I source and then calculating the distance dist of the red endpoints and the green endpoints after the transformation, if dist < 2 pixels then the endpoints are considered to be matching, if at least two endpoints are found to be matching then the transformation is considered to be correct; Step three: point cloud registration; point cloud projection to image scale, image transformation to point cloud transformation, to establish the relationship between point cloud and image, there are transformation scales α and β and displacement in x and y axes, and the point cloud transformation is obtained by combining α, β and displacement parameters; let a set of corresponding points of two station point clouds be P1:(x1,y1,z1), P2:(x2,y2,z2), which satisfy the relationship: P1 = R cloud • P2 + T cloud Wherein, the corresponding points of P1 and P2 on the projected grayscale image are I1:(u1,v1), I2:(u2,v2), which satisfy the relationship: I1 = R Image • I2 + T Image ; wherein R Image and T Image determined from the calculation in step two; The formula for converting image change to point source transformation is: Step four: global optimization; global registration optimization is performed by means of pose optimization to reduce the overall error; specifically comprising the following steps: E) first, the transformation relationship between the point clouds of each station is obtained by steps one to three; F) all the point clouds that can be registered are constructed into a graph optimization problem, the vertices of the graph are the poses of each station point cloud, and the edges are the transformation matrices between each two station point clouds that can be registered; G) the errors of all the point clouds that can be registered are accumulated, and the overall optimization is performed to minimize the error, and finally the ICP point cloud fine registration algorithm is used to further reduce the error.

2. The method of claim 1, wherein, The step C) specifically comprises the following steps: (1) For source image I source and target image I target , n groups of corresponding point pairs are extracted by using SIFT algorithm, if n is less than 2, it is impossible to obtain matching, and it is directly exited; (2) In n pairs of points, only two groups of corresponding points are taken each time. First, the rotation angle θ is calculated through the two groups of corresponding points, and I source is obtained by rotating θ (source-θ) ; (3) Then calculate I through the midpoint of the two sets of corresponding points. target and I (source-θ) The distance Δd between them, and I (source-θ) Translate Δd to obtain I (source-θ-d) ; (4) From the angles and displacements obtained in steps (2) and (3), compute the transformation T according to i，j : (5) Transform T i，j is applied to n feature points of I source , the distance between the n points and their corresponding points of I target is calculated, and the points with a distance less than a threshold are inliers, and the number of inliers m is counted. (6) repeat steps (2)-(5) to get m max inner points, and the transformation is considered as the best transformation, then use the m max inner points together to compute the final image transformation T final .

3. The method of claim 1, wherein, The Levenberg-Marquadt method in the g2o library is used for global registration optimization in the step four.

Images