CN1312633C - Automatic registration method for large-scale three-dimensional scene multi-view laser scanning data - Google Patents

Abstract

The invention relates to an automatic registration method for large-scale three-dimensional scene multi-viewpoint laser scanning data, which comprises the following steps: acquiring data, placing a laser range finder in front of a detected scene, and adjusting the laser range finder to enable a Z axis of the laser range finder to be vertical to the ground, wherein X and Y axes are parallel to the ground; scanning the detected scene row by row and column by column to obtain three-dimensional data of each viewpoint; the scanning data of adjacent viewpoints are required to keep 10-20% of overlap; extracting the structure of the measured object and extracting the characteristics of the data of the measured object; calculating virtual characteristics and constructing characteristic units, performing coarse registration and fine registration and global registration, and finally establishing a registration model, wherein the method can solve the problem of automatic registration of large-scale three-dimensional scene multi-view laser scanning data without any external cooperation equipment and manual participation; the method has universal applicability, and is suitable for scenes of modern structures and complex great-fun cultural relics and ancient buildings.

Description

Extensive many viewpoints of three-dimensional scenic laser scanning data autoegistration method

Technical field

The invention belongs to the three-dimensional data method for registering in computer information processing field, particularly a kind of extensive many viewpoints of three-dimensional scenic laser scanning data autoegistration method.

Background technology

The three-dimensional data registration is the gordian technique of digitizing and reverse engineering.Particularly along with the appearance of medium and long distance three-dimensional laser scanning technique; simplified the three-dimensional data acquisition process; promoted the development of association areas such as digital city, digital navigation, virtual reality, city planning, cultural relic digitalization protection, enjoyed attention based on the 3 D laser scanning data on ground.Describe in order to obtain extensive three-dimensional scenic full geometry, need gather the laser three-dimensional scanning data, so existence is with the laser three-dimensional scanning data-switching of the different points of view registration problems at the same coordinate system at tens even a hundreds of viewpoint position.

The research spininess of traditional many viewpoints three-dimensional data registration is to single object, most representative is the closest approach iterative algorithm, be that this method of ICP algorithm is sought corresponding point by iteration in two viewpoint laser scanning point sets, the position of calculating two viewpoints concerns that weak point need to be the artificial initial conversion of setting to estimate, easily be absorbed in local minimum, iteration time is long, be not suitable for extensive registration and have two kinds for large scale scene method for registering commonly used: target is put in purchased equipment method target reference mark method requirements such as target reference mark method and GPS in the adjacent viewpoint overlapping region, utilize laser that the special reflection characteristic of target is sought the target reference mark automatically, as long as the reference mark of the same name more than three is arranged in the laser scanning data of adjacent viewpoint, if can obtain and have the material that has identical reflection characteristic with target in the coordinate transformation relation scene of adjacent viewpoint, can mistake know target, registration failure target reference mark method needs a large amount of manual interventions, and it mainly is to obtain the pose parameter of laser scanner in each viewpoint by purchased equipment that some environment can't be placed purchased equipment method such as target GPS, equipment price costliness not only, and the pose parameters precision is not high, therefore registration results and actual conditions variant, be badly in need of a kind ofly directly the seamless amalgamation of 3 D laser scanning data of each viewpoint not being analyzed by literature search at automatically extensive many viewpoints of three-dimensional scenic three-dimensional data method for registering of the same coordinate system by any purchased equipment, find that I.Stamos has proposed employing feature straight line and realized adjacent viewpoint three-dimensional laser data automatic registration method in the article of delivering in 2003 that is called " based on the city large scale scene range data autoregistration (AutomatedFeature-Based Range Registration of Urban Scenes of Large Scale) of feature ", and claim that linear feature has reduced geometric complexity, being fit to this method weak point of extensive registration is only to adopt linear feature not have ubiquity, usually fail in addition for the three-dimensional scenic registrations such as large-scale historical relic ancient architecture that linear feature is less, failing to solve because of practical factor such as blocking influences the registration problems that the overlapping region lacks effective linear feature

Said method is intended to solve two viewpoint three-dimensional data registrations of overlapping region, as for tens even a hundreds of viewpoint, must take the effective registering strategy its three-dimensional data can not to be had the cutting edge of a knife or a sword amalgamation after the same coordinate system registration strategies commonly used has sequence registration and synchronous registration sequence registration to refer to a pair of viewpoint registration, another begins registration to viewpoint, and must comprise last to a viewpoint in the viewpoint, all linked with one another till all viewpoints of traversal the advantage of sequence registration strategies be that the same time has only a pair of viewpoint to participate in registration, committed memory is few; Shortcoming is to exist bigger cumulative errors, closed-loop shaped scene particularly, occurring the synchronous registration in slit between the first and last viewpoint then is all viewpoints registrations simultaneously, cumulative errors can not appear, the registration accuracy height, but calculated amount is big, requires high sequence registration and synchronous registration all to be not suitable for large scale scene to computing equipment

Summary of the invention

The purpose of this invention is to provide a kind of not by any equipment for public use directly with the method for registering of extensive many viewpoints of three-dimensional scenic laser scanning data automatic seamless amalgamation at the same coordinate system

Technical scheme of the present invention is as follows:

Extensive many viewpoints of three-dimensional scenic laser scanning data autoegistration method provided by the invention comprises the steps:

1) data are obtained

Because large scale scene is outdoor scene especially, not only comprise main building structure, and contain shelters such as trees Pedestrians and vehicles, and laser passes at random some shelter of transparent substance formation between three-dimensional laser scanner and main building structure, and laser passes transparent substance and forms point at random and be positioned at after the main building structure when measuring, and places laser range finder before tested scene, adjusts three-dimensional laser scanner, make its Z axle vertical ground, X and Y-axis are parallel to the ground; Tested scene is scanned line by line, obtain the three-dimensional data of each viewpoint; The adjacent viewpoint scan-data keep 10%～20% overlapping;

2) extract main building structure

Import the laser three-dimensional scanning data of each viewpoint correspondence, calculate the horizontal range of each vertical scan line up-sampling point and laser range finder, and, utilize hough transform detection of vertical line segment according to horizontal range; Calculate the length of perpendicular line with its contained counting, and to select the horizontal range of the longest vertical line segment correspondence be benchmark, δ before the benchmark ₁δ after rice and the benchmark ₂Data point in the rice scope is main building structure; Described δ ₁Be 1 meter, δ ₂It is 0.5 meter;

Because laser range finder is to sample line by line to scene, for each vertical scan line, the point that belongs to same vertical plane on it is equal to the horizontal range of laser range finder, and main building structure is generally perpendicular to the ground.The present invention analyzes each vertical scan line one by one, and according to the horizontal range of sampled point to laser range finder, the method for employing statistical distribution is determined main building structure to be extracted the regional extent of main building structure;

3) data point of the main building structure that keeps is carried out feature extraction

Three-dimensional feature has reflected the architectural feature of scene inherence, has simplified the geometry complexity of scene, is effective primitive of many viewpoints three-dimensional data registration.The present invention at first adopts the local curve fitting of weighting processing horizontal sweep trace and vertical scan line respectively, calculates a subdifferential of sweep trace each point; Then according to the differential sign of the 2 neighborhood points up and down of point and point with its about the distance detecting scene marginal point of 2 neighborhood points; Utilize the chain code following technology to form different boundary chains then, and the point of same boundary chain is fitted to end to end one group of line segment, be called the feature line segment.The intersection point of adjacent feature line segment is a unique point;

A) each point differential on the sweep trace

Carry out the local secondary match for the sampled point on every sweep trace with (1) formula, and with this point along scan-line direction up and down k neighbor point be fit interval [K, K], K is 3-5;

$\hat{d}\left(i\right)=a_0+a_i+a_2{i}^2---\left(1\right)$

In the formula: the ordering of i-sweep trace up-sampling point, i=0,1 ..., n;

-Di i point is to the calculated level distance of laser range finder;

a ₀, a ₁, a ₂The fitting coefficient of-fit equation (1) formula;

Adopt adaptive weighted approximating method, ask for fitting coefficient a ₀, a ₁, a ₂, and the weighted fitting error function is

$E^{\text{'}}\left(i\right)=\mathrm{\Σw}\left(i\right){\{d\left(i\right)-\hat{d}\left(i\right)\}}^2---\left(2\right)$

In the formula: d (i)-i point is to the real standard distance of laser range finder;

The weights that i is ordered on w (i)-sweep trace, the reliability of representative point; Adopt iteration weights revised law, progressively eliminate the influence of partial points, make only seemingly to close and in reliable point, carry out; Concrete grammar is as follows: suppose that coefficient is a in known the m time iteration ₀ ^m, a ₁ ^m, a ₂ ^m, ask that coefficient is a in the m+1 time iteration ₀ ^M+1, a ₁ ^M+1, a ₂ ^M+1At first count

Calculate the match residual volume of the m time iteration and the weights of the m time iteration;

$e^k\left(i\right)={(d\left(i\right)-{\hat{d}}^m\left(i\right))}^2---\left(3\right)$

$w^{m+1}\left(i\right)=\left\{\begin{array}{cc}1& e^m\left(i\right)\≤b_s\\ \frac{b_s}{e^m\left(i\right)}& e^m\left(i\right)>b_s\end{array}\right.---\left(4\right)$

In the formula: b _s-error allows threshold value;

Initial weight w ⁰(i)=1, constantly revise fitting coefficient and weights with (3), (4) formula, until matched curve satisfy error requirements and comprise count at most till; Mark participates in the point of match;

Then, (1) formula is asked a subdifferential, and calculate a subdifferential of underlined each point on the sweep trace;

$\frac{\∂\hat{d}\left(i\right)}{\∂i}=a_1+2a_{2}i---\left(5\right)$

B) detected edge points

For markd each point on the sweep trace, investigate it along a scan-line direction underlined differential sign up and down at 2, selecting the point of contrary sign is the normal direction point of discontinuity, and is labeled as 0;

For markd each point on the sweep trace, calculate itself and underlined 2 distances about scan-line direction, investigate wherein ultimate range, selecting ultimate range is degree of depth point of discontinuity greater than threshold epsilon, and is labeled as 1;

C) unique point line drawing

At first adopt the chain code following method that marginal point is distinguished into different boundary chains, and with the clockwise ordering of pressing on each boundary chain, the method that merges by division fits to one group of end to end feature line segment with each boundary chain then, concrete grammar is as follows:

Initialization

Boundary chain is divided into isometric plurality of sections L ₀, L ₁..., L _n

Division

Utilize least square method to current each section L _iFit to straight line.Selection from this fitting a straight line farthest and least residual E _iGreater than what allow threshold value T be benchmark a bit, and this section boundary chain is broken as two sections, and continuous division is till satisfy error of fitting; Handle other each section successively;

Merge

Select the match residual error to be merged into one section and to be the feature line segment less than adjacent two sections of threshold value T, the end points of feature line segment is a unique point;

4) calculate imaginary characteristics and structural attitude unit

The marginal information of the corresponding scene of characteristic curve, the angle point or the cusp of the corresponding scene of unique point, yet they have constituted scene basic structure framework, because real scene exists inevitably and blocks phenomenon, if serious shielding, in the overlay region of two viewpoints, lack validity feature, cause registration failure the present invention to propose the imaginary characteristics notion, the disjoint characteristic curve elongated segment of reality is intersected formation imaginary characteristics line segment and imaginary characteristics point, and add in the feature line segment and unique point of said extracted, imaginary characteristics has illustrated the geometry site between the feature;

The present invention is nested into feature unit PU{P|L with unique point together with its two adjacent feature line segments ₁, L ₂:＜n ₁, n ₂, n ₃, p〉}.Feature unit comprises abundant geological information, i.e. three direction vectors and a position vector: two feature line segment L ₁And L ₂Direction vector (n ₁, n ₂), direction vector points to line segment by unique point, and the 3rd vector is the resultant vector of preceding two vectors, reflected that adjacent two feature line segments constitute the normal direction n on plane ₃=n ₁* n ₂, * for the vector multiplication cross position vector be unique point p;

5) registration

The adjacent viewpoint that local registration refers to have the overlapping region is registration in twos.The present invention adopts by thick registration mode to essence, and two viewpoint data are merged fully;

Carry out thick registration earlier: utilize feature unit that the automatic amalgamation of the three-dimensional data of two viewpoints just can be constructed space conversion relation at a pair of feature unit of the same coordinate system: from the direction of characteristic of correspondence line and synthetic method thereof to calculating rotation matrix; From characteristic of correspondence point to calculating the given a pair of viewpoint S that belongs to respectively of translation vector ₁With viewpoint S ₂Feature unit PU{P|L ₁, L ₂:＜n ₁, n ₂, n ₃, P〉} ∈ S ₁And PU ' P ' | L ₁', L ₂':＜n ₁', n ₂', n ₃', p '〉∈ S ₂, rotation matrix is by n ₁And n ₁Calculate ' (i=1,2,3), though following formula error function minimum, and adopt the plain method of quaternary to obtain rotation matrix;

$\mathrm{Err}=\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{\left|\right|n_1^{\′}-R{n}_i\left|\right|}^2---\left(6\right)$

If p and p ' at the coordinate of local coordinate system separately be (x, y, z) with (x ', y ', z '), compute vector;

$t=R\left[\begin{array}{c}x\\ y\\ z\end{array}\right]-\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\\ z^{\′}\end{array}\right]---\left(7\right)$

Travel through N candidate conversion of all feature units structure and estimate, it is as follows with the constraint condition of the feature unit number calculated characteristics units match that obtains mating to adopt matching degree to determine optimal registration conversion estimation matching degree:

Position constraint

P ' is viewpoint S ₂In with viewpoint S ₁In the nearest unique point of p point, and p and between's p ' distance is less than the threshold value Δ ₁

Direction constrain

L _iWith L _iDistance between ' two line segment place straight lines and angular separation are respectively less than Δ ₂And Δ ₃

Carry out smart registration again: through after the thick registration, the basic amalgamation of three-dimensional data of two viewpoints that has the overlapping region is in the same coordinate system, only poorer slightly in some detail section quality of registration, need further smart registration concrete grammar: 1) on the basis of thick registration, determine overlay region 2) overlay region is divided into the N piece, every contains u * v point 3) calculate the mean difference of every corresponding point, finding out the M piece that differs greatly utilizes the ICP algorithm to carry out smart registration in adaptive district as adaptive district, can significantly reduce operation time 4) utilize the laser reflectivity of laser sampling point and put position information, seek corresponding point in conjunction with mahalanobis distance

Above the many viewpoints three-dimensional scenic more than 2, described step also comprises global registration for viewpoint:

Global registration refers to all viewpoint three-dimensional data amalgamations to calculate matching degree at the strategy of the same coordinate system all viewpoints registration in twos at first; Being node with the viewpoint again, is to link power with the matching degree, generates undirected weighted graph; Utilize the minimum spanning tree principle to set up the global registration illustraton of model then and determine a fixed view S by the global registration illustraton of model _a, and the data that make other viewpoint to transform the concrete grammar of wherein setting up the global registration model to it as follows:

If the transition matrix that i calculates two viewpoints that the overlay region is arranged is T _l=[R _i, t _i] and corresponding local matching degree be g (T _i); With the viewpoint is node, with g (T _i) be the node connection weight, generate one and comprise all viewpoints at interior undirected cum rights connected graph G=＜V, E, V is a node, E is the limit, and e _i=g (T _i) ∈ E, and adopt progressively the short circuit method to obtain minimum spanning tree, set up the global registration model

Technique effect of the present invention is:

The present invention is not by any purchased equipment and the artificial method that participates in, having solved extensive many viewpoints of three-dimensional scenic laser scanning data autoregistration it is advantageous that: 1) structural attitude unit, relation between the different characteristic has been described, simplified geometric complexity, be more suitable for extensive registration 2) introduce imaginary characteristics solved because of factor affecting adjacent viewpoint overlay region such as blocking lack the registration problems 3 of validity feature) employing minimum spanning tree principle set up the global registration illustraton of model, registration accuracy reliable 4) method provided by the invention has general applicability, not only is applicable to the scene of modern structure but also be applicable to complicated large-scale historical relic and ancient architecture

Description of drawings

Fig. 1 carries out laser scanning for the present invention to the scene measured object synoptic diagram;

Fig. 2 (a) is a vertical scanning face inner laser sampling point distributions situation;

Fig. 2 (b) is the statistical conditions of vertical scanning face inner laser sampled point by the horizontal projection value;

Fig. 3 (a) and 3 (b) are the example that main building structure is extracted;

Fig. 4 (a), 4 (b), 4 (c) 4 (d) and 4 (e) are two viewpoint three-dimensional data registration examples under the serious shielding situation;

Fig. 5 is a global registration model generative process;

Fig. 6 (a), 6 (b) are the global registration example;

Embodiment

Below in conjunction with accompanying drawing and specific embodiment technical scheme of the present invention is further described:

The present invention requires to import is the laser scanning data that the three-dimensional laser stadimeter of medium and long distance obtains.The operating distance of three-dimensional laser stadimeter is that 5cm/200m concrete implementation step of the present invention is as follows generally greater than 200m along the beam direction measuring accuracy:

Embodiment 1

1) data are obtained

Laser range finder is placed on the preceding 20-50 rice of tested scene, regulate stadimeter and make Z axle vertical ground, X and Y-axis are parallel to ground, tested scene is scanned (as Fig. 1) line by line, the three-dimensional data that obtains corresponding viewpoint in addition, the adjacent viewpoint scan-data keep 10%～20% overlapping so that make the smooth registration of each viewpoint three-dimensional data

2) extract main building structure

Import the laser scanning data of each viewpoint correspondence, calculate the horizontal range of each vertical scanning face each point and laser range finder, and according to horizontal range, utilize hough transform detection of vertical line segment [as Fig. 2 (a) and Fig. 2 (b)] to calculate the length of perpendicular line with its contained counting, selecting the horizontal range of the longest vertical line segment correspondence is benchmark d _Mi, δ before the benchmark ₁δ after rice and the benchmark ₂Data in the rice scope are the distribution plan that main building structure keeps a certain vertical scanning face of Fig. 2 (a) laser sampling point, and as seen from the figure, keeping in the section is main building structure, and it keeps before the section and afterwards scenery and point at random are left out; Fig. 2 (b) is this vertical scanning face laser sampling point horizontal range statistical distribution synoptic diagram among Fig. 2 (a), and as seen from the figure, the point on the main building wall is on same perpendicular line, and promptly horizontal range equates

Present embodiment δ ₁=1.0m, δ ₂=0.5m can extract main building structure, Fig. 3 (a) and 3 (b) extract the panorama measured drawing of instance graph 3 (a) for measured object is carried out laser scanning line by line for main building structure, wherein A is main building structure (metope), and B is ground, and C is the shelter of main building structure; Fig. 3 (b) extracts the figure as a result of main building structure for using said method, and A is the main building structure (metope) that extracts among the figure

3) feature extraction

A) subdifferential of each point on the calculating sweep trace

Carry out the local secondary match for every sweep trace up-sampling point with (1) formula, and with this point along scan-line direction up and down k neighbor point be fit interval [k, k], K is 3-5; (K also can be 4 or 5 to present embodiment k=3 certainly; Can decide according to analyzing spot density,, belong to those skilled in the art of the present technique's required technical knowledge and skills) as for should what be got actually;

$\hat{d}\left(i\right)=a_0+a_{1}i+a_2{i}^2---\left(1\right)$

In the formula: the ordering of i-sweep trace up-sampling point, i=0,1 ..., n; -Di i point is to the calculated level distance of laser range finder;

a ₀, a ₁, a ₂The fitting coefficient of-fit equation (1) formula

Traditional method utilizes least square method to calculate fitting coefficient a according to the error of fitting function E (i) of (2) formula ₀, a ₁, a ₂

$E\left(i\right)=\mathrm{\Σ}{\{d\left(i\right)-\hat{d}(i\left)\right\}}^2---\left(2\right)$

In the formula: d (i)-i point is to the real standard distance of laser range finder;

But least square method is insensitive to point of discontinuity drawn game exterior point, so adopt adaptive weighted curve-fitting method to ask for fitting coefficient a ₀, a ₁, a ₂, formula (3) is in the weighted fitting error function formula: the weights that i is ordered on w (i)-sweep trace, the reliability of representative point.Adopt iteration weights revised law, progressively eliminate the influence of partial points, match is only carried out in reliable point.Concrete grammar is as follows:

Suppose that fitting coefficient is a in known the m time iteration ₀ ^m, a ₁ ^m, a ₂ ^m, ask fitting coefficient a in the m+1 time iteration ₀ ^M+1, a ₁ ^M+1, a ₂ ^M+1At first calculate the match residual volume of the m time iteration and the weights of the m time iteration

$e^k\left(i\right)={(d\left(i\right)-{\hat{d}}^m\left(i\right))}^2---\left(4\right)$

$w^{m+1}\left(i\right)=\{\begin{array}{cc}1& e^m\left(i\right)\≤b_s\\ \frac{b_s}{e^m\left(i\right)}& e^m\left(i\right)>b_s\end{array}---\left(5\right)$

In the formula: b _s-error allows threshold value.

Initial weight w ⁰(i)=1, utilize (4), (5) formula constantly to revise fitting coefficient and weights, to matched curve satisfy error requirements and comprise count at most till.Mark participates in the point of match.

Then, (1) formula is asked a subdifferential, and calculate a subdifferential of underlined each point on the sweep trace.

$\frac{\∂\hat{d}\left(i\right)}{\∂i}=a_1+2a_{2}i---\left(6\right)$

B) detected edge points

For any one markd some p on the sweep trace, and markd up and down neighborhood point (p ₁And p ₂), detect neighborhood point (p ₁And p ₂) symbol of a subdifferential, if contrary sign then the p point be that normal direction is discontinuous, and be labeled as 0.

Calculate p and p ₁And p ₂Distance, if dis=max (dis (p, p ₁), dis (p, p ₂) ε (ε is a threshold value), then the p point is that the degree of depth is discontinuous, is labeled as 1.

C) unique point line drawing

At first adopt the chain code following method that marginal point is distinguished into different boundary chains, and with the clockwise ordering of pressing on each boundary chain, the method that merges by division fits to one group of end to end feature line segment with each boundary chain then, concrete grammar is as follows:

Initialization

Boundary chain is divided into isometric plurality of sections L ₀, L ₁..., L _n

Division

Division

Utilize least square method to current each section L _iFitting to straight line selects from this fitting a straight line farthest and least residual E _iGreater than what allow threshold value T is benchmark a bit, this section boundary chain is broken as two sections, and constantly division is handled other each section successively till satisfying error of fitting

Merge

Select the match residual error to be merged into one section and to be the feature line segment less than adjacent two sections of threshold value T, the end points of feature line segment is a unique point

4) calculate imaginary characteristics and structural attitude unit

Select actual disjoint characteristic curve elongated segment intersect to form in the feature line segment of imaginary characteristics line segment and imaginary characteristics point adding said extracted and the unique point unique point is nested into feature unit PU{P|L together with its two adjacent feature line segments ₁, L ₂:＜n ₁, n ₂, n ₃, p〉} n ₁And n ₂Be feature line segment L ₁And L ₂Direction vector, n ₃=n ₁* n ₂Be the 3rd vector that preceding two vectors are synthetic, * for vector multiplication cross p be unique point

5) registration

A) thick registration: the given a pair of two viewpoint S that belong to respectively ₁And S ₂Feature unit PU{P|L ₁, L ₂:＜n ₁, n ₂, n ₃, p〉} ∈ S ₁And PU ' P ' | L ₁', L ₂':＜n ₁', n ₂', n ₃', p '〉} ∈ S ₂, rotation matrix can be by n so _iAnd n _i' (i=1,2,3) obtain, even following formula error function minimum, and adopt the plain method of quaternary to obtain rotation matrix

$\mathrm{Err}=\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{\left|\right|n_1^{\′}-R{n}_i\left|\right|}^2---\left(7\right)$

If p and p ' at the coordinate of local coordinate system separately be (x, y, z) with (x ', y ', z '), compute vector

$t=R\left[\begin{array}{c}x\\ y\\ z\end{array}\right]-\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\\ z^{\′}\end{array}\right]---\left(8\right)$

The every pair of character pair unit obtains the registration conversion of a correspondence to be estimated, select the matching degree maximum to estimate forbidden word feature unit matching constraints as follows for the optimal registration conversion:

Position constraint

P ' is viewpoint S ₂In with viewpoint S ₁In the nearest unique point of p point, and p and between's p ' distance is less than the threshold value Δ ₁

Direction constrain

L _iWith L _iDistance between ' two line segment place straight lines and angular separation are respectively less than Δ ₂And Δ ₃B) smart registration: estimate by the conversion that thick registration obtains, with the three-dimensional data amalgamation of two viewpoints at the same coordinate system, determine the overlapping region overlay region of two viewpoint correspondences is divided into 10 respectively, calculating 3 of mean difference selection differences maximum of every interior corresponding point utilizes the ICP algorithm to carry out smart registration for increasing the robustness of corresponding point search in adaptive district as adaptive district, utilize the reflectivity and the three-dimensional information of laser sampling point, with the difference of 2 P ∈ S1 of mahalanobis distance measurement and P ' ∈ S2, promptly

$d^2(p,p^{\′})=\underset{p^{\′}\∈S_2}{\min }{(p-p^{\′})}^T{C}_{p^{\′}}^{-1}(p-p^{\′})---\left(9\right)$

In the formula: d ²(p ' p ')-p and p ' mahalanobis distance;

C _{P '}The covariance of-P ' and k * k neighborhood point, present embodiment k=3

Fig. 4 (a), 4 (b) 4 (c) 4 (d) and 4 (e) are hostel two viewpoint three-dimensional data registration results; The laser sampling point cloud of Fig. 4 (a) for obtaining in the S1 viewpoint; The laser sampling point cloud of Fig. 4 (b) for obtaining in the S2 viewpoint; Fig. 4 (c) and Fig. 4 (d) are for utilizing the thick registration results of said method, and wherein Fig. 4 (c) does not adopt imaginary characteristics, and Fig. 4 (d) has adopted imaginary characteristics; Fig. 4 (e) is for to utilize said method to carry out the result of smart registration

Embodiment 2

For using method of the present invention the laser scanning data autoregistration of 8 viewpoints of extensive three-dimensional scenic, its step 1)-5 are carried out in the student dining room) with embodiment 1, different is: because viewpoint outnumbers 2, need carry out global registration;

At first all (overlay region is arranged) viewpoints are carried out registration in twos, utilize the minimum spanning tree principle to set up the global registration model, then according to viewpoint S of global registration Model Selection _aFixing, the laser scanning data conversion of other viewpoint is at S _aIn the coordinate system at place

Specifically, each all can calculate transition matrix T to two viewpoints that the overlay region is arranged _i=[R _i, t _i] and local matching degree g (T _i) (to obtain counting of registration) be node with each viewpoint, with g (T _i) be the connection weight of node, then can generate one and comprise all viewpoints at interior undirected cum rights connected graph G=＜V, E, V is a node, E is limit (weight), and e _i=g (T _i) ∈ E.According to local matching degree g (T _i), at undirected connection weighted graph G=＜V, E〉in, adopt progressively the short circuit method can ask minimum spanning tree, set up the global registration model

Fig. 5 is its global registration model generative process synoptic diagram, and the node 12345 and 6 among the figure has been represented 6 viewpoints; In the 1st step, calculate all viewpoints local matching degree of registration in twos; The 2nd step is for to carry out short circuit with maximized two viewpoints of matching degree; The 3-5 step is carried out short circuit for two viewpoints that matching degree is taken second place, and finishes until all viewpoint short circuits; The 6th step was fixed as Sa for the viewpoint 4 that branch is maximum

Fig. 6 (a) is 8 viewpoint global registrations of present embodiment " student dining room " figure as a result; Fig. 6 (b) is the partial enlarged drawing at window position among Fig. 6 (a).

Claims (2)

1, a kind of extensive many viewpoints of three-dimensional scenic laser scanning data autoegistration method comprises the steps:

1) obtains data

Place laser range finder before tested scene, regulate laser range finder and make its Z axle vertical ground, X and Y-axis are parallel to ground, and tested scene is scanned line by line, obtain the three-dimensional data of each viewpoint; The adjacent viewpoint scan-data keep 10%～20% overlapping;

2) extract main building structure

Import the laser three-dimensional scanning data of each viewpoint correspondence, calculate the horizontal range of each vertical scan line up-sampling point and laser range finder, and, utilize hough transform detection of vertical line segment according to horizontal range; Calculate the length of perpendicular line with its contained counting, and to select the horizontal range of the longest vertical line segment correspondence be benchmark, δ before the benchmark ₁δ after rice and the benchmark ₂Data point in the rice scope is main building structure; Described δ ₁Be 1 meter, δ ₂It is 0.5 meter;

3) data point of the main building structure that keeps is carried out feature extraction

A) each point differential on the sweep trace

Carry out the local secondary match for the sampled point on every sweep trace with (1) formula, and with this point along scan-line direction up and down k neighbor point be fit interval, its fit interval [K, K] is [3,3];

$\hat{d}\left(i\right)={\mathrm{\α}}_0+{\mathrm{\α}}_{1}i+{\mathrm{\α}}_2{i}^2---\left(1\right)$

In the formula: the ordering of i-sweep trace up-sampling point, i=0,1 ..., n;

-Di i point is to the calculated level distance of laser range finder;

α ₀, α ₁, α ₂The fitting coefficient of-fit equation (1) formula;

Adopt adaptive weighted approximating method, ask for fitting coefficient α ₀, α ₁, α ₂, and the weighted fitting error function is

$E^{\′}\left(i\right)=\mathrm{\Σw}\left(i\right)\{d\left(i\right)-\hat{d}\left(i\right){\}}^2---\left(2\right)$

In the formula: d (i)-i point is to the real standard distance of laser range finder;

The weights that i is ordered on w (i)-sweep trace, the reliability of representative point; Adopt iteration weights revised law, progressively eliminate the influence of partial points, match is only carried out in reliable point; Concrete grammar is as follows: suppose that coefficient is α in known the m time iteration ₀ ^m, α ₁ ^m, α ₂ ^m, ask that coefficient is α in the m+1 time iteration ₀ ^M+1, α ₁ ^M+1, α ₂ ^M+1At first calculate the match residual volume of the m time iteration and the weights of the m time iteration;

$e^k\left(i\right)={(d\left(i\right)-\widehat{d^m}\left(i\right))}^2---\left(3\right)$

$W^{m+1}\left(i\right)=\left\{\begin{array}{cc}1& e^m\left(i\right)\≤b_s\\ \frac{b_s}{e^m\left(i\right)}& e^m\left(i\right)>b_s\end{array}\right.---\left(4\right)$

In the formula: b _s-error allows threshold value;

Initial weight w ⁰(i)=1, constantly revise fitting coefficient and weights with (3), (4) formula, until matched curve satisfy error requirements and comprise count at most till; Mark participates in the point of match;

Then, (1) formula is asked a subdifferential, and calculate a subdifferential of underlined each point on the sweep trace;

$\frac{\∂\widehat{d\left(i\right)}}{\∂i}={\mathrm{\α}}_1+2{\mathrm{\α}}_{2}i---\left(5\right)$

B) detected edge points

For markd each point on the sweep trace, investigate it along a scan-line direction underlined differential sign up and down at 2, selecting the point of contrary sign is the normal direction point of discontinuity, and is labeled as 0;

For markd each point on the sweep trace, calculate itself and underlined 2 distances about scan-line direction, investigate wherein ultimate range, selecting ultimate range is degree of depth point of discontinuity greater than threshold epsilon, and is labeled as 1;

C) unique point line drawing

At first adopt the chain code following method that marginal point is distinguished into different boundary chains, and with the clockwise ordering of pressing on each boundary chain, the method that merges by division fits to one group of end to end feature line segment with each boundary chain then, concrete grammar is as follows:

Initialization

Boundary chain is divided into isometric plurality of sections L ₀, L ₁..., L _n

Division

Utilize least square method to current each section L _iFit to straight line, select farthest and least residual E from this fitting a straight line _iGreater than what allow threshold value T be benchmark a bit, and this section boundary chain is broken as two sections, and continuous division is till satisfy error of fitting; Handle other each section successively;

Merge

Select the match residual error to be merged into one section and to be the feature line segment less than adjacent two sections of threshold value T, the end points of feature line segment is a unique point;

4) calculate imaginary characteristics and structural attitude unit

Select actual disjoint characteristic curve elongated segment to intersect and form in the feature line segment and unique point of imaginary characteristics line segment and imaginary characteristics point adding said extracted; Unique point is nested into feature unit PU{P|L together with its two adjacent feature line segments ₁, L ₂:＜n ₁, n ₂, n ₃, p〉}; n ₁And n ₂Be feature line segment L ₁And L ₂Direction vector, n ₃=n ₁* n ₂Be the 3rd vector that preceding two vectors are synthetic, * for vector multiplication cross .p be unique point;

5) registration

A) thick registration

The given a pair of two viewpoint S that belong to respectively ₁And S ₂Feature unit PU{P|L ₁, L ₂:＜n ₁, n ₂, n ₃, p〉} ∈ S ₁And PU ' P ' | L ₁', L ₂':＜n ₁', n ₂', n ₃', p '〉} ∈ S ₂, rotation matrix can be by n so ₁And n ₁' (i=1,2,3) obtain, and make following formula error function minimum that is:, and adopt the plain method of quaternary to obtain rotation matrix;

$\mathrm{Err}=\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{\left|\right|n}_i^{\′}-R{n}_i{\left|\right|}^2---\left(6\right)$

If p and p ' at the coordinate of local coordinate system separately be (x, y, z) with (x ', y ', z '), compute vector;

$t=R\begin{array}{c}\left[\begin{array}{c}x\\ y\\ z\end{array}\right]-\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\\ z^{\′}\end{array}\right]---\left(7\right)\end{array}$

Every pair of character pair unit obtains the registration conversion estimation of a correspondence, and what select the matching degree maximum is that the optimal registration conversion estimates that the feature unit matching constraints is as follows:

Position constraint

P ' is viewpoint S ₂In with viewpoint S ₁In the nearest unique point of p point, and p and between's p ' distance is less than the threshold value Δ ₁

Direction constrain

L ₁With L ₁Distance between ' two line segment place straight lines and angular separation are respectively less than Δ ₂And Δ ₃

B) smart registration

Estimate by the conversion that thick registration obtains, the three-dimensional data amalgamation of two viewpoints at the same coordinate system, is determined the overlapping region; The overlay region of two viewpoint correspondences is divided into the N piece respectively, calculates the mean difference of every interior corresponding point; The M piece of selection differences maximum utilizes the ICP algorithm to carry out smart registration in adaptive district as adaptive district, for increasing the robustness of corresponding point search, utilizes the reflectivity and the three-dimensional information of laser sampling point, with the difference of 2 p ∈ S1 of mahalanobis distance measurement and p ' ∈ S2, promptly

$d^2(p,p^{\′})=\underset{p\′\∈S_2}{\min }{(p-p^{\′})}^T{C}_{P^{\′}}^{-1}(p-p^{\′})---\left(8\right)$

In the formula: d ²(p, p ')-p and p ' 's mahalanobis distance;

C _pThe covariance of '-p ' and k * k neighborhood point, k=3.

By described extensive many viewpoints of the three-dimensional scenic laser scanning data autoegistration method of claim 1, it is characterized in that 2, above the many viewpoints three-dimensional scenic more than 2, described step also comprises global registration for viewpoint:

At first all have the viewpoint of overlay region to carry out registration in twos, utilize the minimum spanning tree principle to set up the global registration model, then according to viewpoint S of global registration Model Selection _aFixing, the laser scanning data conversion of other viewpoint is at S _aIn the coordinate system at place, the concrete grammar of wherein setting up the global registration model is as follows:

If the transition matrix that i calculates two viewpoints that the overlay region is arranged is T ₁=[R ₁, t ₁] and corresponding local matching degree be g (T _i); With the viewpoint is node, with g (T _i) be the node connection weight, generate one and comprise all viewpoints at interior undirected cum rights connected graph G=＜V, E, V is a node, E is the limit, and e ₁=g (T ₁) ∈ E, and adopt progressively the short circuit method to obtain minimum spanning tree, set up the global registration model.

Images