CN104134216B - The laser point cloud autoegistration method described based on 16 dimensional features and system - Google Patents

Abstract

The invention discloses a laser point cloud automatic registration method and system based on 16-dimensional feature description, comprising the following steps: step 1, reorienting the unit normal vector of each laser point in the laser point cloud data; step 2, constructing laser point neighbors The local coordinate system between any two adjacent laser points in the domain; step 3, obtain the feature vector of the laser point in the local coordinate system; step 4, extract the feature point from the point cloud data based on the feature vector of the laser point, and obtain the feature point The optimal scale of ; step 5, based on the feature vectors of the feature points at the optimal scale, the laser point clouds of the two stations are registered. The invention can improve the automation degree and matching accuracy of laser point cloud registration.

Description

Laser point cloud automatic registration method and system based on 16-dimensional feature description

technical field

The invention belongs to the application field of laser point cloud data processing, and in particular relates to a laser point cloud automatic registration method and system based on 16-dimensional feature description.

Background technique

Terrestrial Laser Scanning Technology (Terrestrial Laser Scanning Technology) is a new type of three-dimensional measurement and scanning technology developed in recent years. Many domestic surveying and mapping manufacturers have launched lidar products with independent intellectual property rights. At present, in China, ground laser scanning technology has been applied in various fields such as earthwork calculation, traffic accident handling, urban planning, resource detection, emergency relief, cultural relics protection, etc., but the application ratio of domestic ground laser scanning technology in various industries is still at the low level. One of the main reasons is that the laser point cloud data processing software matched with domestic hardware equipment still has many shortcomings and deficiencies.

The registration of laser point cloud data is the first step in laser point cloud data processing, and it is also the basis for post-processing such as laser point cloud segmentation, classification, and modeling. It is very important in laser point cloud data processing. The registration of laser point cloud data is generally carried out by placing targets and identifying them or by manually selecting points with the same name, but the above methods have great limitations. Therefore, the research on the target-free laser point cloud registration method also highlights its necessity and importance. Target-free laser point cloud registration is mainly based on feature extraction and matching of laser point cloud, but this laser point cloud registration method is difficult to apply to all situations, because the scenes corresponding to laser point cloud data are often more complex, Many algorithms can only perform registration for some of the scenes. Therefore, through continuous improvement and perfection, it is important to find a laser point cloud registration method with good scene adaptability, strong anti-noise ability, and high registration efficiency for the application of ground laser scanning equipment and laser point cloud data in actual production. value.

At present, laser point cloud feature extraction methods mainly focus on geometric feature extraction. This type of feature extraction method calculates more advanced and stable point features by fitting the normal vector and curvature of each laser point. For example, the three-dimensional integral descriptor (calculate the volume of the space formed by the spherical neighborhood of the laser point and the fitting surface passing through the laser point by integral), the sine value of the angle between the normal vector and the direction of the radius of curvature, 3D-SITF features, invariant moments, and spherical harmonics Invariants and other point features. In addition to point features, many methods also use multi-dimensional features such as line features, surface features, ring features, and ball features to describe and extract laser point clouds. After extracting point cloud features, at present, the nearest neighbor search in the feature space is mainly used to determine the same-named points in the laser point cloud for registration, but this method often has many mismatching points; The extracted features are often different.

The following related documents are involved in this paper:

[1] Gelfand N, Niloy J M, Leonidas J G, et al. Robust Global Registration. SGP’05: Proceedings of the third Eurographics symposium on Geometry processing, 2005.197–206.

[2] Liu R, Hirzinger G.. Marker-free automatic matching of range data. Proceedings of In: R. Reulke and U. Knauer (eds), Panoramic Photogrammetry Workshop, Proceedings of the ISPRS working group V/5, 2005.

[3]Sharp G C, Lee S W, Wehe D K.ICP registration using invariant features[J].Pattern Analysis and Machine Intelligence,IEEE Transactions on,2002,24(1):90-102.

[4]Bae K H, Lichti D D.Automated registration of unorganized point clouds from terrestrial laser scanners[M].Curtin University of Technology.2006.

[5] Sadjadi F A, Hall E L. Three-dimensional moment invariants [J]. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 1980(2):127-136.

[6] Burel G, Hénocq H. Three-dimensional invariants and their application to object recognition [J]. Signal Processing, 1995, 45(1): 1-22.

[7 I Stamos, M Leordeanu, Automated Feature-based Range Registration of UrbanScenes of Large Scale, in IEEE Conference on Computer Vision and Pattern Recognition (2003)

[8] J Yao, MR Ruggeri, P Taddei, V Sequeira, Automatic scan registration using 3D linear and planar features. 3D Res. 1(3), 1–18(2010)

[9] C Chao, I Stamos, Semi-automatic Range to Range Registration: a Feature-based Method, in International Conference on 3-D Digital Imaging and Modeling (3DIM) (2005)

[10]C Dold,C Brenner,Registration of Terrestrial Laser Scanning Data Using PlanarPatches and Image Data,in International Society for Photogrammetry and Remote Sensing(2006)

[11] C Chao, I Stamos, Range Image Registration Based on Circular Features, in Proceedings of International Symposium on 3D Data Processing Visualization and Transmission (3DPVT) (2006)

[12]M Franaszek,GS Cheok,C Witzgall,Fast automatic registration of range images from 3D imaging systems using sphere targets.Autom Constr.18(3),265–274(2009).doi:10.1016/j.autcon.2008.08. 003

Contents of the invention

Aiming at the problems existing in the prior art, the present invention provides a laser point cloud automatic registration method and system based on 16-dimensional feature description with more accurate matching.

In order to solve the problems of the technologies described above, the present invention adopts the following technical solutions:

1. A laser point cloud registration method based on 16-dimensional feature description, including steps:

Step 1. Reorient the unit normal vector of each laser point in the laser point cloud data, that is, select the viewpoint, if the angle between the vector vector starting from the viewpoint and ending at the laser point and the unit normal vector of the laser point is greater than 90 degrees , then the unit normal vector of the laser point is reversed; otherwise, the direction of the unit normal vector of the laser point remains unchanged;

Step 2. Construct the local coordinate system between any two laser points in the laser point neighborhood, specifically:

Obtain the acute angle between the unit normal vector of the two neighboring laser points and the line connecting the two neighboring laser points respectively, take the neighboring laser point corresponding to the smaller acute angle as the origin, and the other neighboring laser point as the target point; take the origin The unit normal vector is the u axis, the cross product result of the vector vector with the starting point as the origin and the end point as the target point and the origin unit normal vector is the v axis, and the cross product result of the u axis and the v axis direction vector is the w axis;

Step 3, obtain the feature vector of the laser point in the local coordinate system, specifically:

3.1 Under the local coordinate system, calculate the point multiplication relationship f1 between the unit normal vector of the target point and the direction vector of the v-axis, the distance f2 between the origin and the target point, and the clip between the vector vector with the starting point as the origin and the end point as the target point and the u-axis The angle f3 and the unit normal vector of the target point form the arcsine value f4 of the projection on the plane formed by the u-axis and the v-axis;

3.2 Compare the size of fi and threshold t _i , if fi>t _i , then s(t _i ,fi)=1; otherwise s(t _i ,fi)=0; i=1, 2, 3, 4, t ₁ and t ₂ take values in the range of [-1,1], t ₄ takes values in the range of [-π/2, π/2], and t ₃ represents the scale;

3.3 Obtain the eigenvalues of any two laser points in the laser point neighborhood Count the frequencies of the integers in [0, 15] of the eigenvalues of any two neighboring laser points to form the 16-dimensional eigenvector of the laser point;

Step 4, extract feature points from point cloud data based on the feature vector of the laser point, and obtain the probability combination of the dimensional characteristics of feature points at different scales, and take the scale with the smallest Shannon entropy of the probability combination as the optimal scale of feature points;

Step 5, based on the eigenvectors of the feature points at the optimal scale, the laser point clouds of the two stations are registered.

t ₁ , t ₂ and t ₄ described in sub-step 3.2 are all set to zero.

The feature vector based on the laser point described in step 4 extracts feature points from point cloud data specifically as follows:

According to the laser point feature vectors under different scales, the average feature vectors under each scale are respectively obtained, and the average feature vector is the mean value of all laser point feature vectors in the point cloud data;

At different scales, measure the distance between the laser point feature vector and the average feature vector, and select the initial feature point at the current scale according to the distance between the laser point feature vector and the average feature vector;

The laser point that is the initial feature point on two continuous scales is the final feature point.

The above-mentioned initial feature points at the current scale are selected according to the distance between the laser point feature vector and the average feature vector, specifically:

The laser point whose distance from the average feature vector is greater than the standard deviation σ is selected as the initial feature point, and the standard deviation σ is the standard deviation of the distance between all laser point feature vectors and the average feature vector in the point cloud data.

The distance between the above-mentioned laser point eigenvector and the average eigenvector is measured by KL distance: Among them, D _KL represents the KL distance between the laser point feature vector and the average feature vector, Represents the i-th dimension element of the laser point feature vector, μ _i is the i-th dimension element of the average feature vector.

The probability combination of obtaining the dimensional characteristics of feature points at different scales described in step 4 is specifically:

For the neighborhood laser point set (X ₁ ,...,X _i ,...,X _n ) of the feature point, obtain the matrix and M=B ^T B, where,

Calculate the eigenvalues of the matrix M, and sort the eigenvalues from large to small, and the eigenvalues after sorting are λ ₁ ≥ λ ₂ ≥ λ ₃ ;

According to the eigenvalues of the matrix M, the probability values of the feature point dimensional characteristics are obtained: a ₁ =(λ ₁ -λ ₂ )/λ ₁ , a ₂ =(λ ₂ -λ ₃ )/λ ₁ and a ₃ =λ ₃ / λ ₁ , obtain the probability combination (a ₁ , a ₂ , a ₃ );

At the same time, the probabilistic combination Shannon entropy described in step 4 ${\mathrm{E.}}_{f_{\mathrm{Vp}}}^r=-a_1*{\mathrm{ina}}_1-a_{12}*{\mathrm{ina}}_2-a_3*{\mathrm{ina}}_3.$

Step 6 further includes:

Based on the eigenvectors of the feature points at the optimal scale, the feature points in the laser point clouds of the two stations are roughly registered to obtain the initial point pairs with the same name;

Based on the root mean square error of the distance between the same-name point pairs, the hierarchical greedy method is used to screen the initial same-name point pairs to obtain the screened same-name point pairs;

According to the screened points with the same name, the laser point clouds of the two stations are registered.

Based on the root mean square error of the distance between the same-name point pairs, the layered greedy method is used to screen the initial same-name point pairs to obtain the screened same-name point pairs, specifically:

Merge any two pairs of initial homonymous point pairs whose distance root mean square error is less than the threshold r _Threshold in the initial homonymous point set, and add the second-order point object set;

For any object e _ki in the k-order point object set, search for an object e _kj that has no duplicate node with the object e _ki in the _k -order point object set _. If the root mean square error of the distance is less than the threshold r _Threshold , the object pair (e _ki , e _kj ) is merged into the 2k-order point object set; at the same time, the k-order point object set has the same node as the object pair (e _ki , e _kj ) Among them, k takes 2, 4, and 8 in sequence, and finally obtains a set of 16-order point objects; the threshold r _Threshold is set according to the point cloud density of the two laser point clouds;

Add the objects in the 16th-order point object set to the same-name point pair set, and obtain the conversion parameters (R, t) according to the same-name point pair in the same-name point pair set, and the remaining same-name point pairs (p _i ',q _i ') calculate |R*p _i '+t|-q _i ', and add the remaining point pairs with the same name that |R*p _i '+t|-q _i ' is less than the preset threshold to the set of point pairs with the same name , the preset threshold is set according to the point cloud density of the two laser point clouds.

According to the above-mentioned points with the same name after screening, the laser point clouds of the two stations are registered, specifically:

According to the point pairs with the same name, the conversion parameters between the laser point clouds of the two stations are obtained, and the conversion parameters are used to register the laser point clouds of the two stations.

2. A laser point cloud registration system based on 16-dimensional feature description, including:

(1) The unit normal vector orientation module is used to reorient the unit normal vector of each laser point in the laser point cloud data, that is, to select a viewpoint, if the vector vector with the viewpoint as the starting point and the laser point as the end point is the same as the laser point unit If the included angle of the normal vector is greater than 90 degrees, the unit normal vector of the laser point is reversed; otherwise, the direction of the unit normal vector of the laser point remains unchanged;

(2) The local coordinate system construction module is used to construct the local coordinate system between any two adjacent laser points in the laser point neighborhood. This module further includes sub-modules:

The origin determination module is used to obtain the acute angle included between the unit normal vector of the two adjacent laser points and the line connecting the two adjacent laser points, and the adjacent laser point corresponding to the smaller acute angle is taken as the origin, and the other adjacent laser point as the target point;

The coordinate axis determination module is used to take the origin unit normal vector as the u axis, the cross product result of the vector vector with the origin as the origin and the end point as the target point and the origin unit normal vector as the v axis, and the cross of the u axis and the v axis direction vector The multiplication result is the w axis;

(3) The eigenvector building block is used to obtain the eigenvector of the laser point under the local coordinate system. This module further includes sub-modules:

The geometric feature calculation module is used to calculate the point multiplication relationship f1 between the unit normal vector of the target point and the v-axis direction vector, the distance f2 between the origin and the target point, and the vector with the starting point as the origin and the end point as the target point in the local coordinate system The angle f3 between the vector and the u-axis and the arcsine value f4 of the unit normal vector of the target point projected on the plane formed by the u-axis and the v-axis;

The comparison module is used to compare the size of fi and the threshold t _i , if fi>t _i , then s(t _i ,fi)=1; otherwise s(t _i ,fi)=0; i=1, 2, 3, 4. t ₁ and t ₂ take values in the range of [-1,1], t ₄ takes values in the range of [-π/2, π/2], and t ₃ represents the scale;

The eigenvector obtaining module is used to obtain the eigenvalues of any two adjacent laser points in the laser point neighborhood Count the frequencies of the integers in [0, 15] of the eigenvalues of any two neighboring laser points to form the 16-dimensional eigenvector of the laser point;

(4) The optimal scale acquisition module is used to extract feature points from the point cloud data based on the feature vector of the laser point, and obtain the probability combination of the dimension characteristics of the feature points at different scales, and the minimum scale of the probability combination Shannon entropy is taken as The optimal scale of feature points;

(5) The registration module is used to register the laser point clouds of two stations based on the feature vectors of the feature points at the optimal scale.

Compared with the prior art, the present invention has the following characteristics and beneficial effects:

Combined with Shannon entropy to analyze the optimal scale of feature vectors, obtain the optimal scale of feature points, and extract features and match based on the optimal scale; combined with the invariant distance characteristic in rigid body transformation, use greedy thinking to further screen the initial point pairs with the same name, and obtain possible Set of identically named point pairs to improve matching accuracy.

The invention can improve the automation degree and matching accuracy of laser point cloud registration.

Description of drawings

Fig. 1 is a schematic diagram of the determination process of the origin and the target point;

Fig. 2 is the concrete flowchart of the inventive method;

Fig. 3 is a specific flowchart of obtaining the unit normal vector of the laser point in the laser point cloud.

detailed description

The technical solutions of the present invention will be further described below in conjunction with the drawings and specific embodiments.

The present invention is based on the laser point cloud automatic registration method described by 16-dimensional features, and the specific steps are as follows:

Step 1. Obtain the unit normal vector of each laser point in the laser point cloud data.

The matrix A is obtained according to the laser point p neighbor point Xi = ( _xi , y _{, zi )} :

A＝(X ₁ ,...,X _i ,...,X _n ) ^T (1)

Among them, ( _xi , y _i , _zi ) represent the coordinates of neighborhood point Xi _, and n is the number of laser points in the neighborhood of laser point p.

List the error equation V=AX+L according to the principle of least squares to obtain the normal vector of the laser point p:

X＝( ^AT PA) ^{-1 AT} PL (2)

Among them, V represents the error matrix of size n×1; all elements in the matrix L are -1, and the size is n×1; P represents the weighting matrix of size n×n, and in general, the weighting matrix P is the unit matrix; X Represents a normal vector matrix with a size of 3×1, X=(a,b,c) ^T , that is, the laser point p normal vector n' _p =(a,b,c).

make The unit normal vector n _p =n' _p /d of the laser point p is obtained. The unit normal vector of each laser point in the laser point cloud data is obtained by the above method.

Step 2. Reorient the unit normal vector of the laser point.

Reorient the direction of the unit normal vector of the laser point p to unify the direction of the unit normal vector of the laser point in the point cloud, specifically: let O be the viewpoint, and the coordinates of point O are generally taken as (0, 0, 0). If the laser The unit normal vector n _p of point p and the vector The included angle is greater than 90 degrees, so that the unit normal vector n _p of the laser point p is reversed, as follows:

$\mathrm{if}\frac{<o-{p,\mathrm{no}}_p>}{\left|\right|o-p\left|\right|}<0,$ Let _np = _-np (3)

In formula (3), <Op, n _p > represents the vector Point multiplication operation with the unit normal vector n _p of the laser point p; ||Op|| represents the distance between the viewpoint O and the laser point p.

Step 3, construct the local coordinate system within the neighborhood of the laser point.

For any two neighborhood laser points (X _i , X _j ) in the neighborhood of laser point p, the unit normal vectors of neighborhood laser points Xi and X _j are _ni and _{n j respectively} , and the unit normal vectors _ni and The acute angle between n _j and the line connecting the laser points X _i and X _j in the neighborhood, the laser point in the neighborhood corresponding to the unit normal vector with the smaller angle is taken as the origin X _s , and the other laser point in the neighborhood is the target point X _t .

The method of obtaining the origin point X _s and the target point X _t is as follows:

$\mathrm{if}\frac{<o-{p,\mathrm{no}}_p>}{\left|\right|o-p\left|\right|}<0,$ Let n _p =-n _p (n _i ,n _j )

if<n _i ,X _i -X _j ＞≤<n _i ,X _j -X _i ＞

X _s =X _j ,X _t =X _i (4)

else

X _t =X _j ,X _s =X _i

Use the origin X _s and the target point X _t to construct the local coordinate system (u, v, w), let n _s and n _t be the unit normal vectors of the origin X _s and the target point X _t respectively, and the local coordinate system is defined as follows:

u=n _s

$\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; =</span>\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; t</span>}}}<span\; class="notranslate">\; \&times;</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 5</span>}<span\; class="notranslate">\; )</span>$

w=u×v

The process of determining the origin and the target point will be illustrated below, as shown in Figure 1. The unit normal vector of the neighborhood laser point X _s is n _s , the unit normal vector of the neighborhood laser point X _t is n _t , and the vector The acute angle included with the unit normal vector n _s of the neighborhood laser point X _s is α, and the vector The acute angle included with the unit normal vector n _t of the neighborhood laser point X _t is β, since the acute angle α is smaller than the acute angle β, the neighborhood laser point X _s is taken as the origin, and the neighborhood laser point X _t is the target point . Take the unit normal vector n _s of the origin X _s as the u axis, the unit normal vector n _s and The result of the cross product of is the v axis, and the result of the cross product of the u axis and the direction vector of the v axis is the w axis, thus constructing the local coordinate system X _s -uvw.

Step 4, obtain the 16-dimensional feature description.

For any two neighboring laser points (X _i , X _j ) in the neighborhood of laser point p, construct the local coordinate system X _s -uvw between the two neighboring laser points. In the local coordinate system X _s -uvw, calculate (1) the point product relationship f1=<v,n _t > between the unit normal vector n _t of the target point X _t and the v-axis direction vector, (2) the origin X _s and the target point The distance between X _t f2=||X _t -X _s ||, (3) vector The included angle f3=< _u , X _t -X>/f2 with the _u axis and (4) the _arcsine value f4=atan(<w ,n _t >,<u,n _t >).

Let t _i be the threshold of fi, both t ₁ and t ₂ can take values in the range of [-1,1], and t ₄ can take values in the range of [-π/2, π/2], as the preferred t ₁ , The values of t ₂ and t ₄ are both set to 0; t ₃ is set as the neighborhood radius of the laser point p. If fi>t _i , then s(t _i , fi)=1; otherwise s(t _i ,fi)=0. Based on the value of s(t _i , fi), obtain the statistical characteristic value fx of the relationship between two points (X _i , X _j ):

$\mathrm{<span\; class="notranslate">\; fx</span>}<span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; 4</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}<span\; class="notranslate">\; [</span>{\mathrm{<span\; class="notranslate">\; 2</span>}}^{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; fi</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ]</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 6</span>}<span\; class="notranslate">\; )</span>$

fx is an integer in [0, 15], calculate the statistical feature value fx of any point pair in the neighborhood of laser point p, count the frequency of 16 integers in [0, 15], and constitute the 16 statistics based on the laser point p dimension feature vector.

Step 5, extract feature points based on the 16-dimensional feature vector of the laser point p.

At different scales (i.e., neighborhood radius), calculate the mean μ of the 16-dimensional feature vectors of all laser points in the point cloud data, and measure the distance between the laser point p feature vector and the average feature vector. Specifically, KL distance or Euclidean distance can be used The distance is measured, where the KL distance formula is as follows:

${\mathrm{<span\; class="notranslate">\; D.</span>}}_{\mathrm{<span\; class="notranslate">\; KL</span>}}<span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; 16</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}<span\; class="notranslate">\; [</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}^{\mathrm{<span\; class="notranslate">\; f</span>}}{<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; \&mu;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; In</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}^{\mathrm{<span\; class="notranslate">\; f</span>}}<span\; class="notranslate">\; /</span>{\mathrm{<span\; class="notranslate">\; \&mu;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ]</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 7</span>}<span\; class="notranslate">\; )</span>$

Among them, D _KL represents the KL distance between the laser point p eigenvector and the average eigenvector, Represents the i-th dimension element of the laser point p feature vector, μ _i is the i-th dimension element of the average feature vector.

According to the distance between the laser point p eigenvector and the average eigenvector, select the laser point whose distance from the average eigenvector is greater than the standard deviation σ as the initial feature point, and the standard deviation σ is the eigenvector and the average eigenvector of all laser points in the point cloud data The standard deviation of the distance.

Feature points are extracted at different scales. Only laser points that are initial feature points on two consecutive scales are marked as final feature points. Let P _fi be the set of feature points at r _i scale, then the final set of feature points Among them, P _fi and P _fi+1 are feature point sets under adjacent scales r _i and r _i+1 , and i is the scale number.

Step 6, analyze the optimal scale of feature points, that is, the optimal neighborhood radius.

For each laser point in the point cloud, the radius of the neighborhood used to calculate the feature is different, and the characteristics of the laser point will also be different, even very different. For each laser point in the point cloud, theoretically there is a neighborhood radius r, so that the information in the neighborhood can best describe the characteristics of the laser point.

For laser point p neighborhood laser point set (X ₁ ,...,X _i ,...,X _n ), X _i =( _xi ,y _i , _zi ), calculate its center of gravity According to center of gravity get matrix and M=B ^T B, the M matrix is a symmetric matrix with a size of 3×3. Calculate the eigenvalues of the matrix M, sort the eigenvalues from large to small, and the eigenvalues after sorting are λ ₁ ≥λ ₂ ≥λ ₃ . According to the eigenvalues of the matrix M, the probability values of the p-dimensional characteristic performance of the laser point are obtained: a ₁ =(λ ₁ -λ ₂ )/λ ₁ , a ₂ =(λ ₂ -λ ₃ )/λ ₁ and a ₃ =λ ₃ / λ ₁ , a ₁ +a ₂ +a ₃ =1, a ₁ , a ₂ , and a ₃ are three probability values after unitization.

For any laser point p in the point cloud, the preset neighborhood radius range is [r _min , r _max ]. Take the neighborhood radius with the smallest Shannon entropy of probability combination (a ₁ , a ₂ , a ₃ ) as the best neighborhood radius of laser point p. Shannon entropy Calculated as follows:

${\mathrm{<span\; class="notranslate">\; E.</span>}}_{{\mathrm{<span\; class="notranslate">\; f</span>}}_{\mathrm{<span\; class="notranslate">\; Vp</span>}}}^{\mathrm{<span\; class="notranslate">\; r</span>}}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; a</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; *</span>{\mathrm{<span\; class="notranslate">\; ina</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; a</span>}}_{\mathrm{<span\; class="notranslate">\; 12</span>}}<span\; class="notranslate">\; *</span>{\mathrm{<span\; class="notranslate">\; ina</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; a</span>}}_{\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; *</span>{\mathrm{<span\; class="notranslate">\; ina</span>}}_{\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 8</span>}<span\; class="notranslate">\; )</span>$

The optimal scale of the feature points is obtained by the above method.

Shannon entropy is used to solve the problem of quantitative measurement of information. If you want to understand something unclear, you need a lot of information. If you want to know something you already know a little bit about, you don't need much information. Therefore, the amount of uncertainty in one thing can represent the amount of information in the calculation.

For the neighborhood of laser point p, a ₁ , a ₂ , a ₃ represent the possibility that the point cloud belongs to lines, surfaces and scattered points. In the optimal scale analysis, the thing to be understood is whether the point cloud is finally presented as a line, a plane or a scattered point at the laser point p. Therefore, the present invention utilizes the Shannon entropy to find the scale within [r _min ,r _max ] that minimizes the Shannon entropy of events occurring at the laser point p.

Step 7, the same-name point search and mismatch screening under the optimal scale.

At the optimal scale, the 16-dimensional feature vector of the feature point is calculated. For the two site clouds to be registered, use the KD tree (k-dimension tree) to perform the nearest search in the 16-dimensional feature space to obtain the coarse registration point, that is, the initial point pair with the same name.

The transformation between different laser point cloud stations is a rigid body transformation, and the distance between the corresponding laser points does not change. Therefore, this feature can be used to filter the initial point pairs with the same name. In the present invention, the root mean square error of the distance between the initial point of the same name is used to screen the mismatching points, and the root mean square error is calculated as follows:

${\mathrm{<span\; class="notranslate">\; wxya</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; p</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; q</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{{\mathrm{<span\; class="notranslate">\; no</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; no</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}\underset{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; no</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{<span\; class="notranslate">\; (</span><span\; class="notranslate">\; |</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">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 9</span>}<span\; class="notranslate">\; )</span>$

Among them, the point pairs (p _i , q _i ) and (p _j , q _j ) are the i-th pair and the j-th pair of initial vertex pairs with the same name, that is, any two pairs of vertex pairs with the same name.

Based on the root mean square error of the distance between pairs of points with the same name, a layered greedy algorithm is used to screen the initial set of points with the same name. The specific steps are as follows:

7.1 Use formula (9) to calculate the root mean square error dRMS(p,q) of the distance between any two pairs of initial points with the same name (p _i ,q _i ) and (p _j ,q _j ), at this time, formula (9) where n is 2; if dRMS(p,q) is less than the threshold r _Threshold , the initial point pairs (p _i ,q _i ) and (p _j ,q _j ) with the same name will be merged into structure e ₂ , and the structure e ₂ Join the second-order point object set E ₂ . The threshold r _Threshold is set according to the point cloud density of the two laser point clouds, and generally takes 1 to 10 times the smaller point cloud density of the two laser point cloud densities.

Perform the above operations on all initial point pairs with the same name in different point cloud stations, and sort the objects in the second-order point object set E ₂ according to dRMS(p,q) from small to large.

7.2 Merge the objects in the _second -order point object set E2.

According to sorting, for the object e _2i in the second-order point object set E ₂ , the operations are as follows:

Search for an object e _2j that has no duplicate nodes with object e _2i in the second-order point object set E ₂ , and for the object pair (e _2i , e _2j ), calculate the dRMS(p,q ), at this time, n in formula (9) is 4; if dRMS(p,q) is less than the threshold r _Threshold , the object pair (e _2i , e _2j ) is merged into the fourth _- order point object set E ₄ , and at the same time Remove the object that has the same node as the object pair (e _2i , e _2j ). Then repeat the above operation for the remaining objects in _E2 .

After the fourth-order point object set E ₄ is generated, sort the objects in it according to dRMS(p,q) from small to large.

7.3 According to sub-step 7.2, continuously merge the point object sets of each order, that is, generate the point object set E _2k of order 2k according to the point object set E _k of order k until the point object set E ₁₆ of order 16 is obtained.

7.4 Add the objects in the 16th-order point object set E ₁₆ to the same-name point pair set, calculate the conversion parameters (R, t) according to the same-name point pairs in the same-name point pair set, and apply the conversion parameters (R, t) to the initial same-name point pair set. Add the same-named point pair (p _i ',q _i ') in the 16-order point object set E ₁₆ , if |R*p _i '+t|-q _i ' is less than the preset threshold, then the same-named point pair (p _i ',q _i ') are also added to the point pair set with the same name. The preset threshold is set according to the point cloud densities of the laser point clouds of the two stations, and generally takes 1 to 10 times the smaller point cloud density of the laser point cloud densities of the two stations.

Step 8: Obtain the conversion parameters between the clouds of the two sites according to the point pairs with the same name in the same-name point pairs, that is, the optimal conversion parameters, and use the optimal conversion parameters to register the laser point clouds of the two sites.

Claims (8)

1. The laser point cloud registration method based on 16-dimensional feature description, is characterized in that, comprises steps: Step 1, obtain the unit normal vector of each laser point in the laser point cloud data: The matrix A is obtained according to the laser point p neighbor point Xi = ( _xi , y _{, zi )} : A＝(X ₁ ,...,X _i ,...,X _n ) ^T Among them, ( _xi , y _i , _zi ) represent the coordinates of neighborhood point Xi _, and n is the number of laser points in the neighborhood of laser point p; List the error equation V=AX+L according to the principle of least squares to obtain the normal vector of the laser point p: X＝( ^AT PA) ^{-1 AT} PL Among them, V represents the error matrix of size n×1; all elements in the matrix L are -1, and the size is n×1; P represents the weighting matrix of size n×n, and in general, the weighting matrix P is the unit matrix; X Represents a normal vector matrix with a size of 3×1, X=(a,b,c) ^T , that is, the normal vector n' _p =(a,b,c) of the laser point p; make Obtain the unit normal vector n _p =n' _p /d of the laser point p; Step 2, re-orientate the unit normal vector of each laser point in the laser point cloud data, that is, select the viewpoint, if the angle between the vector vector starting from the viewpoint and ending at the laser point and the unit normal vector of the laser point is greater than 90 degrees , then the unit normal vector of the laser point is reversed; otherwise, the direction of the unit normal vector of the laser point remains unchanged; Step 3, construct the local coordinate system between any two laser points in the laser point neighborhood, specifically: Obtain the acute angle between the unit normal vector of the two neighboring laser points and the line connecting the two neighboring laser points respectively, take the neighboring laser point corresponding to the smaller acute angle as the origin, and the other neighboring laser point as the target point; take the origin The unit normal vector is the u axis, the cross product result of the vector vector with the starting point as the origin and the end point as the target point and the origin unit normal vector is the v axis, and the cross product result of the u axis and the v axis direction vector is the w axis; Step 4, obtain the feature vector of the laser point in the local coordinate system, specifically: 4.1. In the local coordinate system, calculate the point product relationship f1 between the unit normal vector of the target point and the v-axis direction vector, the distance f2 between the origin and the target point, and the clip between the vector vector with the starting point as the origin and the end point as the target point and the u-axis The angle f3 and the unit normal vector of the target point form the arcsine value f4 of the projection on the plane formed by the u-axis and the v-axis; 4.2, compare the size of fi and the threshold t _i , if fi>t _i , then s(t _i ,fi)=1; otherwise s(t _i ,fi)=0; i=1, 2, 3, 4, t ₁ and t ₂ take values in the range of [-1,1], t ₄ takes values in the range of [-π/2, π/2], and t ₃ represents the scale; 4.3. Obtain the eigenvalues of any two laser points in the laser point neighborhood Count the frequencies of the integers in [0, 15] of the eigenvalues of any two neighboring laser points to form the 16-dimensional eigenvector of the laser point; Step 5, extract feature points from point cloud data based on the feature vector of the laser point, and obtain the probability combination of feature point dimensional characteristics at different scales, so that the scale with the smallest Shannon entropy of the probability combination is the best scale of feature points; The probability combination of obtaining the dimensional characteristics of feature points at different scales is specifically: For the neighborhood laser point set (X ₁ ,...,X _i ,...,X _n ) of the feature point, obtain the matrix and M=B ^T B, where, Calculate the eigenvalues of the matrix M, and sort the eigenvalues from large to small, and the eigenvalues after sorting are λ ₁ ≥ λ ₂ ≥ λ ₃ ; According to the eigenvalues of the matrix M, the probability values of the feature point dimensional characteristics are obtained: a ₁ =(λ ₁ -λ ₂ )/λ ₁ , a ₂ =(λ ₂ -λ ₃ )/λ ₁ and a ₃ =λ ₃ / λ ₁ , obtain the probability combination (a ₁ , a ₂ , a ₃ ); At the same time, the probability combination Shannon entropy Step 6, based on the eigenvectors of the feature points at the optimal scale, the laser point clouds of the two stations are registered. 2. the laser point cloud registration method described based on 16 dimension features as claimed in claim 1, is characterized in that: t ₁ , t ₂ and t ₄ described in sub-step 4.2 are all set to zero. 3. the laser point cloud registration method described based on 16 dimension features as claimed in claim 1, is characterized in that: The feature vector based on the laser point described in step 5 extracts feature points from point cloud data specifically as follows: According to the laser point feature vectors under different scales, the average feature vectors under each scale are respectively obtained, and the average feature vector is the mean value of all laser point feature vectors in the point cloud data; At different scales, measure the distance between the laser point feature vector and the average feature vector, and select the initial feature point at the current scale according to the distance between the laser point feature vector and the average feature vector; The laser point that is the initial feature point on two continuous scales is the final feature point. 4. the laser point cloud registration method described based on 16 dimension features as claimed in claim 3, is characterized in that: The initial feature point under the current scale is selected according to the distance between the laser point feature vector and the average feature vector, specifically: Select the laser point whose distance from the average feature vector is greater than the standard deviation σ as the initial feature point, and the standard deviation σ is the standard deviation of the distance between all laser point feature vectors and the average feature vector in the point cloud data. 5. the laser point cloud registration method described based on 16 dimension features as claimed in claim 3, is characterized in that: The distance between the measured laser point feature vector and the average feature vector is measured by KL distance: Among them, D _KL represents the KL distance between the laser point feature vector and the average feature vector, Represents the i-th dimension element of the laser point feature vector, μ _i is the i-th dimension element of the average feature vector. 6. the laser point cloud registration method described based on 16 dimension features as claimed in claim 1, is characterized in that: Step 6 further includes: Based on the eigenvectors of the feature points at the optimal scale, the feature points in the laser point clouds of the two stations are roughly registered to obtain the initial point pairs with the same name; Based on the root mean square error of the distance between the same-name point pairs, the hierarchical greedy method is used to screen the initial same-name point pairs to obtain the screened same-name point pairs; According to the filtered point pairs with the same name, the laser point clouds of the two stations are registered. 7. the laser point cloud registration method described based on 16 dimension features as claimed in claim 6, is characterized in that: The registration of the two laser point clouds according to the screened point pairs with the same name is specifically: The conversion parameters between the two laser point clouds are obtained according to the same-name point pairs in the same-name point pair collection, and the conversion parameters are used to register the two laser point clouds. 8. A laser point cloud registration system based on 16-dimensional feature description, characterized in that it comprises: (1) The unit normal vector acquisition module is used to obtain the unit normal vector of each laser point in the laser point cloud data: The matrix A is obtained according to the laser point p neighbor point Xi = ( _xi , y _{, zi )} : A＝(X ₁ ,...,X _i ,...,X _n ) ^T Among them, ( _xi , y _i , _zi ) represent the coordinates of neighborhood point Xi _, and n is the number of laser points in the neighborhood of laser point p; List the error equation V=AX+L according to the principle of least squares to obtain the normal vector of the laser point p: X＝( ^AT PA) ^{-1 AT} PL Among them, V represents the error matrix of size n×1; all elements in the matrix L are -1, and the size is n×1; P represents the weighted matrix of size n×n; X represents the normal vector matrix of size 3×1, X =(a,b,c) ^T , i.e. the normal vector n' _p =(a,b,c) of the laser point p; make Obtain the unit normal vector n _p =n' _p /d of the laser point p; (2) The unit normal vector orientation module is used to reorient the unit normal vector of each laser point in the laser point cloud data, that is, to select a viewpoint, if the vector vector with the viewpoint as the starting point and the laser point as the end point is the same as the laser point unit If the included angle of the normal vector is greater than 90 degrees, the unit normal vector of the laser point is reversed; otherwise, the direction of the unit normal vector of the laser point remains unchanged; (3) The local coordinate system construction module is used to construct the local coordinate system between any two adjacent laser points in the laser point neighborhood. This module further includes sub-modules: The origin determination module is used to obtain the acute angle included between the unit normal vector of the two adjacent laser points and the line connecting the two adjacent laser points, and the adjacent laser point corresponding to the smaller acute angle is taken as the origin, and the other adjacent laser point as the target point; The coordinate axis determination module is used to take the origin unit normal vector as the u axis, the cross product result of the vector vector with the origin as the origin and the end point as the target point and the origin unit normal vector as the v axis, and the cross of the u axis and the v axis direction vector The multiplication result is the w axis; (4) eigenvector building block, used to obtain the eigenvector of the laser point under the local coordinate system, this module further includes sub-modules: The geometric feature calculation module is used to calculate the point multiplication relationship f1 between the unit normal vector of the target point and the v-axis direction vector, the distance f2 between the origin and the target point, and the vector vector with the starting point as the origin and the end point as the target point in the local coordinate system The angle f3 with the u-axis and the arcsine value f4 of the projection of the unit normal vector of the target point on the u-axis and the v-axis to form a plane; The comparison module is used to compare the size of fi and the threshold t _i , if fi>t _i , then s(t _i ,fi)=1; otherwise s(t _i ,fi)=0; i=1, 2, 3, 4. t ₁ and t ₂ take values in the range of [-1,1], t ₄ takes values in the range of [-π/2, π/2], and t ₃ represents the scale; The eigenvector obtaining module is used to obtain the eigenvalues of any two adjacent laser points in the laser point neighborhood Count the frequencies of the integers in [0, 15] of the eigenvalues of any two neighboring laser points to form the 16-dimensional eigenvector of the laser point; (5) The optimal scale acquisition module is used to extract feature points from point cloud data based on the feature vector of the laser point, and obtain the probability combination of feature point dimension performance at different scales, so that the scale with the smallest Shannon entropy of the probability combination is the feature The optimal scale of points; The probability combination of obtaining the dimensional characteristics of feature points at different scales is specifically: For the neighborhood laser point set (X ₁ ,...,X _i ,...,X _n ) of the feature point, obtain the matrix and M=B ^T B, where, Calculate the eigenvalues of the matrix M, and sort the eigenvalues from large to small, and the eigenvalues after sorting are λ ₁ ≥ λ ₂ ≥ λ ₃ ; According to the eigenvalues of the matrix M, the probability values of the feature point dimensional characteristics are obtained: a ₁ =(λ ₁ -λ ₂ )/λ ₁ , a ₂ =(λ ₂ -λ ₃ )/λ ₁ and a ₃ =λ ₃ / λ ₁ , obtain the probability combination (a ₁ , a ₂ , a ₃ ); At the same time, the probability combination Shannon entropy (6) The registration module is used to register the laser point clouds of two stations based on the feature vectors of the feature points at the optimal scale.