CN106780459A - A kind of three dimensional point cloud autoegistration method - Google Patents

Abstract

The invention discloses a method for automatic registration of three-dimensional point cloud data, comprising the following steps: sampling two point clouds to be registered to obtain feature points, calculating the rotation invariant feature of these feature points, The rotation-invariant feature sub-characteristics of the feature points are matched and searched to obtain the initial correspondence between the feature points; then the random sampling consensus algorithm is used to judge and remove the false matching points existing in the initial matching point set, and the optimized feature point correspondence is obtained. relationship, calculate the roughly rigid transformation relationship between the two point clouds, and realize rough registration; and provide a rigid transformation consistency detection algorithm, use the local coordinate system transformation relationship between the matching feature points to perform rough registration results Constraint detection, to complete the verification of the correctness of the rough registration results; and use the ICP algorithm to optimize the rigid transformation relationship between point cloud data, and finally realize the automatic and accurate registration of point clouds.

Description

A method for automatic registration of 3D point cloud data

technical field

The invention belongs to the field of three-dimensional measurement, and more specifically relates to an automatic registration method for three-dimensional point cloud data.

Background technique

The 3D model is of great significance in the fields of industrial inspection, cultural relics protection, biomedicine and so on. With the development of 3D measurement technology, the processing technology of point cloud model has become a research hotspot in recent years. An important step in point cloud processing is to unify the point cloud data obtained from different viewpoints in the same scene into the same coordinate system, that is, point cloud registration.

The point cloud registration methods are mainly divided into the following categories: (1) Manual registration: Manually select the corresponding points between two point clouds to achieve merging. This method requires manual intervention, and the effect of merging often depends on the experience of the operator; ( 2) Registration method based on external assistance: using a turntable or pasting mark points on the object to achieve auxiliary registration, this type of method has high environmental requirements, and the registration accuracy is easily affected by external assistance; (3 ) registration method based on its own topographic features: use the topographic features of the measured object to achieve matching between point clouds and then realize automatic registration of point clouds, without using additional equipment or marker points, and can meet different measurement occasions However, the stability and accuracy of the algorithm need to be strengthened compared with the assisted point cloud registration technology.

Aiming at the registration method based on its own shape features, some research units have carried out related research and achieved certain results; Rusu used the normal vector angle between the feature point and the surrounding neighbor points as a feature to construct a fast point feature histogram Fast Point Feature Histograms (FPFH) are used to search and match feature points. After completing the initial registration to obtain an ideal initial position, the precise registration algorithm is used to improve the registration accuracy; Chen and Besl proposed an iterative closest point algorithm ( Iterative closest point (ICP) is used to perform accurate registration of two point cloud data. This method iteratively finds the closest point in the corresponding point set as the corresponding point, continuously optimizes the rigid transformation matrix, and finally obtains accurate registration. However, this method has high requirements for the initial positions of the two point clouds, and the initial relative positions between the point clouds need to be roughly similar.

Existing point cloud registration methods usually use a combination of coarse registration and precise registration. Among them, coarse registration generally finds corresponding points by calculating the characteristics of points, and estimates the initial position relationship between the two point clouds; The registration criterion is to use ICP and its improved algorithm to further optimize the rough registration results to achieve precise registration of point cloud data; And interference factors such as point cloud density changes are relatively sensitive, and a large number of mismatching points are prone to appear, which reduces the accuracy and stability of registration.

Contents of the invention

The purpose of the present invention is to provide an automatic registration method for three-dimensional point cloud data in view of the existing technical problems, which has high registration accuracy and strong stability, and completes automatic judgment of the correctness of registration results.

In order to achieve the above object, the technical solution of the present invention is an automatic point cloud registration method based on rotation-invariant feature descriptors. First, the two point clouds to be registered are sampled to obtain feature points, and the rotation invariance of these feature points is calculated respectively. Variable feature sub, to search for the rotation invariant features of the feature points in the two point clouds, and obtain the initial correspondence between the feature points; then use the random sampling consensus algorithm (Random Sample Consensus, RANSAC) Judging and removing the mis-matched points, obtaining the optimized corresponding relationship of feature points, calculating the roughly rigid transformation relationship between the two point clouds, and realizing rough registration; and providing a rigid transformation consistency detection algorithm, using the matching The local coordinate system transformation relationship between the feature points is constrained to detect the rough registration results, and the correctness of the rough registration results is verified; finally, the iterative closest point algorithm is used to optimize the rigid transformation relationship between the point cloud data, and finally the point cloud is realized. Automatic precise registration.

Specifically include the following steps:

Step 1: Read in the collected source point cloud P and target point cloud Q to be registered;

Step 2: Calculate the density of the source point cloud P and the target point cloud Q respectively, and randomly select several points from the source point cloud P to form the source feature point set S ₁ , and randomly select several points from the target point cloud Q Constitute the target feature point set S ₂ ;

Step 3: Calculate the local rotation-translation invariant coordinate system of each feature point in the source feature point set S ₁ and the target feature point set S ₂ according to the density of the source point cloud and the target point cloud;

Step 4: Calculate the high-dimensional feature description of each feature point according to the local rotation and translation invariant coordinate system of the feature point; perform feature point matching on the target point cloud Q and the source point cloud P, and obtain the initial matching point set C;

Step 5: Use the random sampling consensus algorithm to remove the mismatching in the initial matching point set obtained in the step 4, use the rigidity change estimation algorithm based on the singular value decomposition method to calculate and obtain the rotation matrix R and the translation matrix T, and obtain the a coarse registration result between the target point cloud and said source point cloud;

Step 6: Use a registration error detection algorithm based on rigid transformation consistency to determine whether the rough registration result obtained in step 5 is correct, and if the result is correct, go to step 7, and if the result is incorrect, return the result of registration failure;

Step 7: Iteratively optimize the rotation matrix R and translation matrix T obtained by the above estimation by using the iterative closest point algorithm, so as to realize the automatic and accurate registration of the target point cloud and the source point cloud.

Generally speaking, compared with the prior art, the above technical solutions conceived by the present invention can achieve the following beneficial effects:

(1) The automatic registration method of 3D point cloud data provided by the present invention, by selecting a feature point set from the point cloud data, constructs a local rotation-translation invariant coordinate system for each feature point, and calculates the coordinate system on the feature point High-dimensional feature description realizes the replacement of traditional artificial markers, and can be used in application scenarios such as digital scanning of cultural relics, which can reduce the workload of manual pasting and removal of markers;

(2) The automatic registration method of 3D point cloud data provided by the present invention can effectively realize the search and matching of corresponding points between different point cloud data by using the obtained feature point set and high-dimensional feature description, and realize Automatic registration of point cloud data;

(3) The method for automatic registration of three-dimensional point cloud data provided by the present invention compares the result of the local rotation matrix calculated by using the local rotation-invariant coordinate system with the result of the global rotation-translation matrix obtained by using the random sampling consistency algorithm, Judging whether there is rotation-translation consistency between them can further effectively judge whether the rough registration result is successful, and realize precise registration.

Description of drawings

Fig. 1 is a schematic diagram of the overall flow of the method for automatic registration of 3D point cloud data according to the present invention.

detailed description

In order to make the object, technical solution and advantages of the present invention more clear, the present invention will be further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described here are only used to explain the present invention, not to limit the present invention. In addition, the technical features involved in the various embodiments of the present invention described below may be combined with each other as long as they do not constitute a conflict with each other.

The method for automatic registration of three-dimensional point cloud data provided by the embodiment of the present invention has a process as shown in FIG. 1 and includes the following steps:

Step 1: Read in the target point cloud and source point cloud collected by the 3D measuring equipment;

Step 2: Randomly select several points in the source point cloud and the target point cloud to form two sets of feature points, specifically:

According to the preset feature point sampling ratios h ₁ and h ₂ , randomly sample from the source point cloud P and the target point cloud Q respectively to obtain the corresponding feature point sets

Among them, n ₁ refers to the number of feature points in the feature point set of the source point cloud P, n ₂ refers to the number of feature points in the feature point set of the target point cloud Q, n ₁ =h ₁ ·N ₁ , n ₂ =h ₂ • N ₂ .

Step 3: Calculate the local rotation and translation invariant coordinate system of each feature point for the sampled feature point set; specifically include the following sub-steps:

Step 3.1: For each point p _i in the source point cloud and each point q _i in the target point cloud Q, use the nearest neighbor search method to calculate the nearest neighbor point p _i ′ of point p _i and the nearest point q _i neighboring point q _i ′

Obtain the point cloud density d _p of point cloud p _i , and the point cloud density d _q of point cloud q _i ;

Among them, ||p _i -p _i ′|| is the Euclidean distance between point p _i and p _i ′, ||q _i -q _i ′|| is the distance between point q _i and point q _i ′ Euclidean distance; i=1,2,3,...,N, N is the number of points in the point cloud data; the source point cloud P contains N ₁ points, for the source point cloud P, N=N ₁ ; Q in the target point cloud contains N ₂ points, and for the target point cloud Q, N=N ₂ ;

Step 3.2: Randomly select a point p _i in the set of feature points, and construct the covariance matrix of point p _i :

Among them, w _i ＝1/||p _j -p _i ||, The radius r _loc is the calculated radius of the neighborhood of p _i ; in the embodiment, r _loc =15d _p ;

Step 3.3: use the following formula to solve the eigenvalues and eigenvectors of the covariance matrix;

COV(p _i )V=EV

where E is given by The 3×3 diagonal matrix formed by

is the eigenvalue Corresponding eigenvectors, l=1,2,3;

Step 3.4: Put Set to x ⁺ , y ⁺ , z ⁺ respectively, The reverse vector of is set to x ^- , y ^- , z ^- ;

The direction of the z-axis in the local coordinate system is determined by the following formula:

Step 3.5: Use the method of step 3.4 to determine the direction of the x-axis in the local coordinate system; the direction of the y-axis in the local coordinate system is determined by z×x, thus establishing a local rotation-translation invariant coordinate system F _i = with p _i as the origin {x, y, x},; use the same method as the above steps 3.2 to 3.4 to establish a local rotation and translation invariant coordinate system with the point q _i as the origin.

Step 4: According to the local rotation and translation invariant coordinate system of the calculated feature points, calculate the high-dimensional feature description of the feature points, match the feature points of the two point clouds, and obtain the initial matching point set; specifically include the following sub-steps:

Step 4.1: For any feature point p _i in the source point cloud feature point set, its high-dimensional features are described in the form of high-dimensional vectors In the high-dimensional space formed by the feature geometry corresponding to the feature point set of the target point cloud, search and nearest feature and next closest feature

In the present invention, the search of the nearest neighbor feature is accelerated calculation by fast approximate nearest neighbor search method (FastApproximate Nearest Neighbor Search);

Step 4.2: Calculate features to the nearest feature next closest feature Euclidean distance of with calculate and Ratio, and judge the size of the ratio and the preset distance ratio threshold τ;

Specifically, it is determined according to the following formula and Is there a correct correspondence between

like then it shows and The match is successful, otherwise the match fails; if and If the matching is successful, the corresponding feature points in the source point cloud and the target point cloud are also correctly corresponding, and the point pair constitutes a matching point pair;

The above processing is also performed on other feature points of the source point cloud, and the initial matching point set C is formed by all matching point pairs.

Step 5: Use Random Sample Consensus (RANSAC) to remove false matches in the initial matching point set, and use Singular Value Decomposition (SVD) to calculate the rotation matrix R and translation matrix T, and obtain the distance between the source point cloud and the target point cloud. Coarse registration relationship among them, specifically including the following sub-steps:

Step 5.1: Set the number of random sampling S _num and the initial residual error E _rr , randomly sample from the initial matching point set C, and select 3 pairs of matching points as the initial point;

If there are multiple matching points of a selected feature point, one of them is randomly selected from the corresponding point set as the corresponding point of the search point;

Step 5.2: Use the SVD algorithm to solve the rotation matrix R and translation matrix T, as follows:

Step 5.21: Assuming that the matching point set obtained in step 5.1 is P={p ₁ ,p ₂ ,p ₃ ,}, P′={p′ ₁ ,p′ ₂ ,p′ ₃ }, use the following formula to calculate the point set Centroid;

Wherein, p _k and p′ _k are the three-dimensional coordinates of any pair of matching points;

Step 5.22: Use the following formula to translate the point sets P and P′ relative to their respective centroids to obtain new point sets Q={q ₁ ,q ₂ ,q ₃ ,} and Q′={q′ ₁ ,q′ ₂ , q′ ₃ };

q _k =p _k -g,q' _k =p' _k -g',(k=1,2,3);

Step 5.23: Calculate the 3×3 matrix M using the following formula:

Step 5.24: performing singular value decomposition M=UΛV ^T on the M matrix;

Among them, the superscript T is the transposition of the matrix, U and V are 3×3 unitary matrices, Λ is a 3×3 diagonal matrix, and the 3×3 diagonal matrix A is defined as:

Step 5.25: Use the following formula to calculate the 3×3 rotation matrix R and the 3×1 translation matrix T:

R = UAV ^T , T = g'-Rg;

Step 5.3: For all matching point pairs in the initial matching point set, calculate the distance error d _err after rotation and translation according to the following formula;

d _err ＝||p _τ -(R·q _τ +T)|| ₂ ;

If the distance error d _err is smaller than the given initial residual error E _rr , it is determined that the matching point is an interior point; all interior points are calculated according to the above method and the number m′ of interior points is counted;

In the embodiment, the residual error is the distance threshold; (p _τ , q _τ ) is a matching point pair in the initial matching point set, τ=1,2,...r, r is the number of matching point pairs in the initial matching point set ;

Step 5.4: If the number of interior points m' is greater than the given number of interior points m, use the points in the interior point data set to re-solve the rotation matrix R and translation matrix T using the SVD algorithm;

Calculate the residual error E′ _rr according to the following formula;

If the residual error is smaller than the initial residual error, then use the rotation-translation matrix calculated at this time as the best estimated target model parameter to update the rotation-translation matrix R, T and the initial residual error E′ _rr ; otherwise, go to step 5.5; if the interior point is less than the number of given interior points, then go directly to step 5.5;

Step 5.5: Repeat random sampling S _num times, repeat steps 5.2 to 5.4 to obtain the number of inliers corresponding to the group of samples, sort the number of inliers in all samples, and select the sampling result with the largest number of inliers as the best sampling;

Using the interior point data set C obtained under the optimal sampling, solve the rotation matrix R and translation matrix T according to step 5.2 as the optimal rotation and translation matrices R _ran and T _ran ; where C ₁ contains S pairs of matching points.

Step 6: Use the rigid transformation consistency detection algorithm to check whether the rough registration relationship is correct;

Step 6.1: Use the following formula to calculate the rotation-translation matrix (R _loc , T _loc ) between matching point pairs (p _loc , q _loc ) in the initial matching point set C:

Among them, R _loc is the local rotation matrix between two matching points, and T _loc is the local translation matrix between two matching points;

F _ploc and F _qloc are local coordinate systems established at points p _loc and q _loc respectively, loc=1,2,...S, S is the number of matching point pairs in matching point set C;

Step 6.2: Convert the obtained rotation matrix R _loc into Euler angle representation, namely

R _loc →(α _loc ,β _loc ,γ _loc ); convert the rotation matrix R _ran obtained in step 5.6 into Euler angle representation, that is, R _ran →(α _ran ,β _ran ,γ _ran );

Step 6.3: Calculate the angle difference d _a between the Euler angles according to the following formula:

Among them, Δ(η ₁ , η ₂ ) ² = (η ₁ -η ₂ ) ² , η ₁ , η ₂ are Euler angles;

Step 6.4: Calculate according to the following formula The distance difference d _t between T _ran = (t _x , t _y , t _z ) ^T :

Step 6.5: According to whether d _a and d _t are both smaller than the given thresholds σ _a and σ _t , determine whether the rotation and translation matrices obtained by the two methods are consistent;

In the embodiment, σ _a takes 0.5236, σ _t takes 15D _den , and D _den is the point cloud density, that is, when the rotation matrix difference is less than 30°, and when the translation matrix distance is less than 15 times the point cloud density, it indicates that the rotation calculated in this step The translation matrix is consistent with the rotation-translation matrix obtained in step 5;

Step 6.6: Perform steps 6.1 to 6.5 for other matching point pairs in the inlier data set C obtained under the best sampling, and obtain the local rotation-translation matrix between all matching point pairs and the above-mentioned optimal rotation-translation matrix (R _ran , T _ran ) consistency relationship between;

The final statistics obtain the number s of all matching point pairs satisfying the consistency relationship in the interior point data set C obtained under the best sampling;

Step 6.7: Calculate the consistency ratio λ=s/S of the number s of matching point pairs satisfying the consistency relationship and the number S of matching point pairs in the matching point set C;

If λ≥threshold τ _λ , it indicates that the rotation matrix obtained by using RANSAC solution is consistent with the rotation matrix calculated by using the local rotation-invariant coordinate system, and it is determined that the registration is successful; otherwise, it is determined that the registration has failed; in the embodiment , the threshold τ _λ =0.7.

Step 7: Use the improved ICP algorithm to optimize the rigid transformation relationship between point clouds to realize automatic and accurate registration of point clouds; the rigid transformation relationship includes the rotation and translation matrix; the specific steps are as follows:

Step 7.1: Set the distance threshold ω as the condition for iteration termination; where, ω>0; the distance threshold ω is based on the point cloud density _dp of the source point cloud P; in the embodiment, ω= _5dp ;

Step 7.2: Randomly select several points in the source point cloud as points to be matched;

Step 7.3: use the reverse projection method to find the corresponding point of the point to be matched in the source point cloud in the target point cloud;

Step 7.4: Using the point-to-plane distance measurement as the objective function to be solved by the ICP algorithm, iteratively calculate the rigid transformation relationship from the source point cloud to the target point cloud;

Step 7.5: When the objective function value is less than the distance threshold ω, stop the iteration; and use the rigid transformation relationship obtained at this time as the final result to complete the automatic precise registration of the point cloud.

Those skilled in the art can easily understand that the above descriptions are only preferred embodiments of the present invention, and are not intended to limit the present invention. Any modifications, equivalent replacements and improvements made within the spirit and principles of the present invention, All should be included within the protection scope of the present invention.

Claims (7)

1. A kind of three-dimensional point cloud data automatic registration method, is characterized in that, comprises the steps: Step 1: Read in the collected source point cloud P and target point cloud Q to be registered; Step 2: Calculate the density of the source point cloud P and the target point cloud Q respectively, and randomly select several points from the source point cloud P to form the source feature point set S ₁ , and randomly select some points from the target point cloud Q points constitute the target feature point set S ₂ ; Step ₃ : Calculate the local rotation-translation invariant coordinate system of each feature point in the source feature point set S1 and the target feature point set S2 according to the density _; Step 4: Calculate the high-dimensional feature description of each feature point according to the local rotation-translation invariant coordinate system; perform feature point matching on the target point cloud Q and the source point cloud P, and obtain the initial matching point set C; Step 5: Use the random sampling consensus algorithm to remove the false match of C in the initial matching point set, calculate and obtain the rotation matrix R and translation matrix T by using the singular value decomposition method, and obtain the rough match between the target point cloud and the source point cloud quasi-results; Step 6: Judging whether the rough registration result satisfies the rigid transformation consistency condition, if so, proceed to step 7; if not, judge that the registration fails; Step 7: Iteratively optimize the rotation matrix R and translation matrix T by using the iterative closest point algorithm to obtain the precise registration result of the target point cloud Q and the source point cloud P. 2. The three-dimensional point cloud data automatic registration method as claimed in claim 1, is characterized in that, described step 2 comprises following sub-steps: Step 2.1: Use the nearest neighbor search method to calculate the nearest neighbor point of each point in the source point cloud P and the target point cloud Q; And calculate the point cloud density of each point cloud Among them, ||calculation point-nearest neighbor point|| refers to the Euclidean distance between the calculation point and the nearest neighbor point; the source point cloud P contains N ₁ points, and for the source point cloud P, N=N ₁ ; Q in the target point cloud contains N ₂ points, and for the target point cloud Q, N=N ₂ ; Step 2.2: Set the feature point sampling ratio h ₁ and h ₂ , and randomly sample from the source point cloud P and the target point cloud Q respectively to obtain the feature point set of the source point cloud P and the feature point set of the target point cloud Q Among them, i is the serial number of the feature point in the feature point set; n ₁ refers to the number of feature points in the feature point set of the source point cloud P, n ₂ refers to the number of feature points in the feature point set of the target point cloud Q, n ₁ =h ₁ ·N ₁ , n ₂ =h ₂ ·N ₂ . 3. the three-dimensional point cloud data automatic registration method as claimed in claim 2, is characterized in that, comprises following sub-step in the described step 3: Step 3.1: Randomly select a point p _i in the feature point set S ₁ of the source point cloud P, obtain the neighborhood calculation radius of p _i according to the point cloud density of p _i , and use the following formula to construct the covariance matrix of point p _i : $\mathrm{<span\; class="notranslate">\; C</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; V</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\underset{<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; <</span>{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; c</span>}}}{<span\; class="notranslate">\; \&Sigma;</span>}{\mathrm{<span\; class="notranslate">\; w</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; )</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; /</span>\underset{<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; <</span>{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; c</span>}}}{<span\; class="notranslate">\; \&Sigma;</span>}{\mathrm{<span\; class="notranslate">\; w</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}$ Among them, w _i =1/||p _j -p _i ||, The radius r _loc is calculated for the neighborhood of p _i ; For a point q _i in the feature point set S ₂ of the target point cloud Q, also perform the above processing to obtain the covariance matrix COV(q _i ) of q _i ; Step 3.2: Solve the eigenvalues and eigenvectors of the covariance matrix; COV(p _i )V=EV where E is given by The 3×3 diagonal matrix formed by is the eigenvalue Corresponding eigenvectors, l=1,2,3, Step 3.3: According to the feature vector Construct the x, y, z coordinate axes of the local rotation and translation invariant coordinate system, and establish the local rotation and translation invariant coordinate system F _i ={x, y, x} with the point p _i as the origin; Step 3.4: According to the covariance matrix of the point q _i in the feature point set S ₂ of the target point cloud Q, adopt the method described in steps 3.2 to 3.3 to establish a local rotation-translation invariant coordinate system with the point q _i as the origin. 4. three-dimensional point cloud data automatic registration method as claimed in claim 3, is characterized in that, comprises following sub-step in the described step 4: Step 4.1: For any feature point p _i in the feature point set of the source point cloud, its high-dimensional features are described in the form of high-dimensional vectors In the high-dimensional space formed by the feature geometry corresponding to the feature point set of the target point cloud, search and nearest feature and next closest feature Step 4.2: Calculate features to the nearest feature next closest feature Euclidean distance of with Judgment according to the following formula and Is there a correct correspondence between like then it shows and The match is successful, otherwise the match fails; if and If the matching is successful, the corresponding feature points in the source point cloud and the target point cloud are also correctly corresponding, and this point pair constitutes a matching point pair; all matching point pairs form an initial matching point set. 5. three-dimensional point cloud data automatic registration method as claimed in claim 4, is characterized in that, comprises following sub-step in the described step 5: Step 5.1: set the number of times of sampling S _num and the distance threshold μ; each time sampling, randomly select t pairs of matching points from the set of initial matching points as the initial point; Step 5.2: Solve the rotation matrix R and translation matrix T by using the singular value decomposition method; Step 5.3: For all matching points in the initial matching point set, calculate the distance error d _err after rotation and translation; if the distance error d _err is smaller than the distance threshold μ, then determine that the matching point is an inlier, and find all inliers and Count the number m' of interior points; Among them, d _err ＝||p _i -(R·q _i +T)|| ₂ , (p _i ,q _i ) is a pair of matching points in the initial matching point set, i=1,2,…r,r is the number of matching point pairs in the initial matching point set; Step 5.4: Repeat random sampling S _num times, and for each sampling result, obtain the number of inliers corresponding to this sampling according to steps 5.2 to 5.3; Sort the number of interior points in all samples, and select the sampling result with the largest number of interior points as the best sampling; use the interior point data set C ₁ obtained under the best sampling, and obtain the rotation matrix R and translation matrix T according to the method in step 5.2, as The optimal rotation-translation transformation matrices are denoted as R _ran and T _ran . 6. The three-dimensional point cloud data automatic registration method as claimed in claim 5, is characterized in that, described step 5.2 comprises following sub-steps: Step 5.21: Calculate the centroid of the initial matching point set according to the following formula; ${\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; t</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; t</span>}}{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; t</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; t</span>}}{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; ;</span>$ Among them, the initial matching point set is P={p ₁ ,p ₂ ,…,p _t }, P′={q ₁ ,q ₂ ,…,q _t }; Step 5.22: Translate the initial matching point sets P and P' relative to the centroid obtained in step 5.21, and obtain new point sets P'={p' ₁ ,p' ₂ ,...,p' _t } and Q' ={q′ ₁ ,q′ ₂ ,…,q′ _t }; p' _k ＝p _k -p ₀ , q' _k ＝q _k -q ₀ Step 5.23: Compute the t×t matrix M: $\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; t</span>}}{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}^{<span\; class="notranslate">\; \&prime;</span>}{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}^{<span\; class="notranslate">\; \&prime;</span>\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; ;</span>$ Step 5.24: Carry out singular value decomposition on the M matrix M=UΛV ^T , where the superscript T is the transposition of the matrix, U and V are 3×3 unitary matrices, Λ is a 3×3 diagonal matrix, define 3×3 The diagonal matrix A of : $\mathrm{<span\; class="notranslate">\; A</span>}<span\; class="notranslate">\; =</span>\left\{\begin{array}{cc}{\mathrm{<span\; class="notranslate">\; I</span>}}^{\mathrm{<span\; class="notranslate">\; 3</span>}}\mathrm{<span\; class="notranslate">\; det</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; det</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; V</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&Greater\; Equal;</span>\mathrm{<span\; class="notranslate">\; 0</span>}\\ \mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; a</span>}\mathrm{<span\; class="notranslate">\; g</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; det</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; )</span>\mathrm{<span\; class="notranslate">\; det</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; V</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; <</span>\mathrm{<span\; class="notranslate">\; 0</span>}\end{array}\right.$ Step 5.25: Compute the t×t rotation matrix R and the t×1 translation matrix T: R=UAV ^T , T=q ₀ −R·p ₀ . 7. The three-dimensional point cloud data automatic registration method as claimed in claim 5, is characterized in that, comprises following sub-step in the described step 6: Step 6.1: Perform the following processing on each matching point pair in the interior point data set C ₁ obtained under the optimal sampling: For a matching point pair (p _loc , q _loc ) in the interior point data set C ₁ , the local rotation matrix R _loc and the local translation matrix T _loc between the two matching points are solved: ${\mathrm{<span\; class="notranslate">\; R</span>}}_{\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; c</span>}}<span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; f</span>}}_{{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; c</span>}}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; T</span>}}{\mathrm{<span\; class="notranslate">\; f</span>}}_{{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; c</span>}}}$ T _loc =p _loc -q _loc R _loc in are the local coordinate systems established at points p _loc and q _loc respectively, loc=1, 2,...S, S is the number of matching point pairs in the interior point data set C ₁ ; Step 6.2: Calculate the angular distance between the rotation matrix R _loc obtained in this step and the rotation matrix R _ran obtained in step 5.4: ${\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; a</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; a</span>}\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; (</span>\frac{\mathrm{<span\; class="notranslate">\; t</span>}\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; a</span>}\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; e</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; R</span>}}_{\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; c</span>}}{\mathrm{<span\; class="notranslate">\; R</span>}}_{\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; a</span>}\mathrm{<span\; class="notranslate">\; no</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ;</span>$ Step 6.3: Calculate the translation distance between the translation matrix T _loc obtained in this step and the translation matrix T _ran obtained in step 5.4: ${\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; t</span>}}<span\; class="notranslate">\; =</span>\frac{<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; c</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; a</span>}\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>}{{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}}<span\; class="notranslate">\; ;</span>$ Step 6.4: If the above d _a and d _t are both smaller than the given thresholds σ _a and σ _t , it means that the rotation-translation matrix obtained in step 5 and step 6 is consistent; Step 6.5: Perform steps 6.1 to 6.4 on other matching point pairs in the interior point data set C ₁ to obtain the local rotation-translation matrix between all matching point pairs and the rotation-translation matrix obtained in step 5 (R _ran ,T _ran ) the consistency relation between; Count the number s of all matching point pairs satisfying the consistency relationship in the interior point data set _C1 ; Step 6.6: Calculate the consistency ratio λ=s/S between s and the number S of matching point pairs in the interior point data set _C1 ; If λ≥threshold τ _λ , it indicates that the local rotation matrix calculated using the local rotation-invariant coordinate system is consistent with the global rotation-translation matrix obtained by using the random sampling consensus algorithm in step 5, and it is determined that the registration is successful; otherwise , which is judged as registration failure.