CN111709981A

CN111709981A - Registration method of laser point cloud and analog image with characteristic line fusion

Info

Publication number: CN111709981A
Application number: CN202010576614.XA
Authority: CN
Inventors: 高小翎; 高宏松
Original assignee: Individual
Current assignee: Individual
Priority date: 2020-06-22
Filing date: 2020-06-22
Publication date: 2020-09-25

Abstract

The invention provides a registration method of laser point cloud and analog image with fused characteristic lines, which aims at the problems that the registration method of the image and the point cloud based on the point cloud distance intensity image or point characteristic in the prior art has a plurality of defects of manual point selection error, point cloud data interpolation error, need of a plurality of observed values and the like, and the precision is very limited. The method for registering the laser point cloud and the simulated image based on the two-dimensional and three-dimensional fusion characteristic lines is provided, a direct conversion model between the image characteristic lines and the point cloud characteristic lines is established, the requirement on the initial value of external orientation elements of the image is not harsh, the robustness on noise is strong, the automatic matching of the two-dimensional and three-dimensional fusion characteristics is realized, the automation degree is improved, the higher real-time performance is achieved, the method has the advantages of no model error, no need of the initial value of the external orientation elements, only a small number of registration elements and the like, the calculation speed is high, and the registration requirement of the simulated image and the ground laser point cloud data is met.

Description

Registration method of laser point cloud and analog image with characteristic line fusion

Technical Field

The invention relates to a registration method of laser point cloud and a simulation image, in particular to a registration method of laser point cloud and a simulation image with characteristic line fusion, and belongs to the technical field of image registration methods.

Background

The photogrammetry mainly obtains the three-dimensional coordinates of the target through a plurality of images with certain overlapping degree and a small number of ground control points, has mature development, is widely applied to the fields of urban surveying and mapping, topographic survey, geographic national conditions census investigation and the like, has rich image texture and semantic information, but needs to obtain the three-dimensional information of the target through homonymy point matching and forward intersection, and is easily influenced by factors such as complexity of ground objects, shielding among the ground objects and the like, so that the high-precision three-dimensional information of the target is difficult to extract only by using the images. The GPS/IMU assisted photogrammetry technology still has the obvious problems of low positioning precision, insufficiently dense homonymous matching points, easily sheltered facade of urban buildings, easily influenced by topographic relief due to image distortion and the like, and simultaneously, because the urban buildings are densely distributed and the urban surveying task is heavy and complex, the actual requirement is difficult to meet only depending on photogrammetry, and the ground survey cannot be completely replaced. The geodetic three-dimensional laser scanning is an active sensor technology for rapidly acquiring three-dimensional spatial data, and actively emits laser pulse signals with certain frequency and wavelength to a target surface through a certain direction angle and a certain height angle, records the time difference of emitted and received signals and the intensity value of echo signals, and acquires the three-dimensional coordinates and the reflection characteristics of a measuring point. The geodetic three-dimensional laser scanning technology integrates spectral information such as three-dimensional coordinate data, reflection intensity, object color and the like, can truly reflect the form, structure and spectral characteristics of a target, and provides abundant data information for target identification and analysis.

Laser scanning measurement has many characteristics and advantages: the method has the advantages that firstly, the method is active and real-time, the detection of the space position and the spectrum information of the target is carried out by receiving the laser pulse echo signal emitted by the method, and the method is not limited by the conditions such as time, weather and the like; secondly, the precision of the coordinate data is ensured due to the fact that the target information is acquired by using a dot matrix or grid which is small in sampling interval and uniform in distribution based on quasi-parallel laser beams and high in density and high in precision; thirdly, digitalization and automation are realized, the processing software directly converts the laser pulse signal into a digital signal, and the full digital characteristic is realized; and fourthly, the system is non-contact and penetrable, a complete non-contact mode is adopted for scanning and measuring the target, the system has strong applicability, laser can penetrate through vegetation which is not dense enough to reach the surface of the target, three-dimensional data of the ground is obtained, and the problem that a simulation image is easily influenced by shielding is avoided. As the laser scanning measurement meets the use requirements and development trends of the surveying and mapping industry, obvious technical advantages and wide application prospects are shown in numerous engineering fields.

However, due to the limitations of the measurement technology itself or the characteristics of the target to be measured, obtaining three-dimensional geodetic information only by ground photogrammetry or geodetic three-dimensional laser scanning cannot meet the actual requirements, and the reasons mainly include: firstly, the precision of the image-based measurement technology is difficult to meet the requirement, and particularly, the automation degree is low, so that simple and rapid visualization is only met; and secondly, the pure point cloud data is only high-precision surface three-dimensional data, lacks texture information of a target surface, and cannot meet the requirements of users on environment perception and immersion experience in scene-type applications such as virtual reality and the like. The laser scanning data and the analog image have respective advantages for describing the target, if the urban modeling completely adopts the camera to carry out data acquisition, only the reconstruction of the concept and the abstract model can be completed, and the ground laser scanning is adopted to carry out data acquisition, so that the efficiency and the precision of the measurement operation and the three-dimensional modeling can be improved. Therefore, the method and the device for realizing the registration of the ground laser point cloud and the simulated image are a precondition for the fusion processing of the two data, and the registration problem of the two data is provided. Because the laser point cloud and the image data have strong complementarity, the integration of the laser point cloud and the image data has wide application in building extraction modeling, urban surveying and mapping and the like. However, due to the installation error of the optical camera and other reasons, the element of the external orientation of the simulated image has a large error, which results in that the simulated image and the laser point cloud cannot be well registered, so that before the two kinds of data are comprehensively applied, the two kinds of data need to be brought into a unified coordinate system, namely, the simulated image and the laser point cloud are registered.

For the registration problem of the laser point cloud and the analog image, the prior art has many limitations, and the following is further summarized as follows: firstly, the method strongly depends on point features, needs a large number of registration primitives and has high calculation complexity; the registration method for generating the point cloud intensity image is based on a two-dimensional registration model, and the registration algorithm for generating the point cloud by image matching is a three-dimensional registration model, so that the defect that the registration model is not strict exists; thirdly, errors are introduced in the processes of point cloud interpolation, projection and image matching; fourthly, manual point selection is needed, and the point characteristics in the point cloud are discrete and difficult to accurately select; fifthly, point characteristics or straight line characteristics need to be manually paired, so that the automation degree is low, and manual errors exist to influence the analysis precision; sixthly, the initial value of the external orientation element of the image is relied on, and the iteration times are more; and seventhly, the method is sensitive to noise and local convergence easily occurs.

In view of the fact that the registration of the laser point cloud and the simulation image in the prior art mainly depends on the characteristics of two-dimensional and three-dimensional fusion points, artificial point selection errors exist, and initial values of positions and postures are needed, the invention realizes the registration method of the laser point cloud and the simulation image based on two-dimensional and three-dimensional fusion characteristic lines under the condition of no initial values of position and posture parameters, judges the registration accuracy, and comprises the following three aspects: firstly, an image characteristic line extraction method, a judgment condition of short straight line segment combination processing and a combination method are adopted; secondly, a high-precision plane cutting and characteristic line extracting method for earth three-dimensional laser point cloud data; and thirdly, a function model based on the position and attitude estimation problem of the two-dimensional and three-dimensional fusion characteristic line, a position and attitude parameter analysis method and an accuracy judgment method under the condition of no initial value.

Disclosure of Invention

Aiming at the defects of the prior art, the registration method of the laser point cloud and the analog image with the fused characteristic lines provided by the invention has the advantages that the constraint effect of registration elements is strong, the required minimum registration element number is small, a strict registration model between two-dimensional and three-dimensional fusion registration elements is established, the requirement on the initial value of external orientation elements of the image is not strict, the calculated time cost is in direct proportion to the observed value number, the robustness on noise is stronger, the automatic matching of two-dimensional and three-dimensional fusion characteristics is realized, the automation degree is improved, the higher real-time performance is achieved, the method has universality and portability, the popularization and application potential is huge, and the wide market application prospect is possessed.

In order to achieve the technical effects, the technical scheme adopted by the invention is as follows:

the registration method of laser point cloud and analog image with fused characteristic line establishes a direct conversion model between the image characteristic line and the point cloud characteristic line, the two-dimensional and three-dimensional fusion registration method based on the characteristic line comprises a camera imaging model and a camera calibration and a position and posture estimation method of a plurality of groups of fused characteristic lines, the position and posture estimation method of the plurality of groups of fused characteristic lines comprises a characteristic line fusion registration function model, a quaternion array expression method of a rotation matrix, an analytic registration parameter and a registration precision judgment method;

extracting laser point cloud characteristic lines by adopting a characteristic line extraction method based on point cloud cutting plane intersection, managing point cloud by adopting KD-Tree, estimating normal vector and curvature by using a covariance analysis method, realizing point cloud plane cutting by adopting a region growing rough cutting method and a random sampling consistency plane fitting fine cutting method, and finally solving intersection of adjacent point cloud planes to extract high-precision straight line characteristics; extracting characteristic lines of the laser point cloud comprises point cloud preprocessing, cutting of a point cloud plane and intersection of the point cloud plane, wherein the point cloud preprocessing comprises point cloud thinning and density homogenization and point cloud denoising, the cutting of the point cloud plane comprises point cloud rough cutting and plane fitting based on region growing, and the intersection of the point cloud plane and the characteristic lines comprises estimation of adjacent straight line parameters, judgment of plane adjacent relation and solving of adjacent straight line segments;

extracting a characteristic line of the simulated image by adopting a characteristic line extraction algorithm based on the gradient direction of the image, searching a similar pixel region based on a gradient tension value and the gradient direction, determining a support region of a straight line segment by adopting a rectangular structure approximation method, extracting parameters of the straight line segment based on the center, the length and the width of a rectangular structure, and providing a short straight line segment merging judgment condition based on a straight line included angle and a distance and a merging processing method comprehensively considering the direction and the relative length of the straight line segment; the method for extracting the characteristic line of the simulated image comprises a straight line segment detection algorithm, straight line parameter estimation, short straight line segment merging processing and characteristic line screening.

The registration method of laser point cloud and simulated image with fused characteristic lines, further, the position and attitude estimation method of several groups of fused characteristic lines respectively extracts characteristic lines from the simulated image and geodetic three-dimensional laser point cloud, analyzes the position and attitude parameters of the simulated image under the scanner coordinate system to implement registration of the two, in the characteristic line fusion registration function model,

the imaging model of an ideal optical camera is:

q_pi＝C_3×3q_ci

q_ci＝B_3×3q_wi+G_3×1

G_3×1＝[g_xg_yg_z]^Tformula 2

The above formula is a functional model of a two-dimensional and three-dimensional fusion registration method based on point features, wherein q is_piCoordinates of the image point in the image coordinate system, C₃x₃Is an internal parameter matrix of the camera, corresponding to internal orientation elements in photogrammetry, which is obtained by calibrating the camera with optical distortion parameters, B₃x₃、G₃x₁The method comprises the following steps that a rotation matrix and a translation matrix from an object space coordinate system to a camera coordinate system are respectively, correspond to external orientation elements in photogrammetry and are parameters to be solved for two-dimensional and three-dimensional fusion registration;

the invention analyzes the registration parameter B by utilizing the geodetic three-dimensional laser point cloud and the characteristic line on the image₃x₃、G₃x₁Establishing a two-dimensional and three-dimensional fusion function model based on the characteristic line; the method for describing the straight line comprises the following steps:

firstly, the parameter of any straight line under the three-dimensional rectangular coordinate system is K ═ q_i,u_iWherein q is_iIs the three-dimensional coordinate of any point on the straight line, v_iIs the unit direction vector of the line, K_iThe coordinate of any point is

Wherein a is any real number;

secondly, the parameter of a certain straight line under the object space coordinate system is K_wi＝{q_wi，u_wi}，q_wi＝[X_iY_iZ_i]^T，u_wi＝[D_iE_iF_i]^TI.e. the equation of the line is D_iX+E_iY+F_iZ+H_i＝0；

Thirdly, the parameter of a certain straight line under the camera coordinate system is K_ci＝{q_ci，u_ci}，q_ci＝[x_iy_if_i]^TI.e. the equation of the line is d_ix+e_iy+j_iz+h_i＝0；

For any group of homonymous two-dimensional and three-dimensional fusion straight lines K_ci、K_wiAnd, object space straight line K_wiWarp B₃x₃、G₃x₁And converting the transformed image into a camera coordinate system to obtain a straight line K_{ci_}＝{q_{ci_}，u_{ci_}Due to K_wiAnd K_ciIs a straight line of the same name, according to the camera imaging principle, K_{ci_}With a straight line K under the camera coordinate system_ciAnd the projection center S of the camera is coplanar, and the plane is recorded as N_i，N_iThe normal vector of (a) is recorded as M_iThe coplanarity condition is a basic condition of two-dimensional and three-dimensional fusion registration based on the characteristic line, and specifically comprises two constraint conditions: one is K_wiDirection vector u of_wiThrough the rotation matrix B₃x₃Direction vector u under camera coordinate system obtained after transformation_{ci_}And M_iVertically; second, K_wiAny point p on_wiThrough B₃x₃、G₃x₁Point p in camera coordinate system obtained after transformation_{ci_}To plane N_iIs 0, i.e. satisfies the plane N_iThe equation of (c); the two constraints are expressed by the equation:

and because:

then the following results are obtained:

the above formula is a function model of a two-dimensional and three-dimensional fusion registration algorithm based on characteristic lines;

note plane N_iNormal vector of (1)

Plane N_iThe equation in the camera coordinate system is:

M_i1x+M_i2y+M_i3z+M_i40-0 formula 6

Plane N_iNormal vector M of_iBased on camera principal distance f and K_ci＝{q_ci，u_ciFind out, because,

and the projection center S of the camera is equal to [0,0 ]]^TIn plane N_iSubstituting into formula 6 to obtain N_i40, then plane N_iThe equation of (a) is:

(y_ij_i-e_if)x+(d_if-x_ij_i)y+(x_ie_i-d_iy_i) z is 0 or 9

Combining formula 4 with formula 8 yields:

namely:

where K is the registration parameter to be found, B₁₁…B₃₃Is 9 elements of a rotation matrix B, which is an orthonormal matrix of units, B₁₁…B₃₃Are not independent of each other.

The invention discloses a registration method of a laser point cloud and an analog image fused with a characteristic line, and further, photogrammetry adopts three independent angle parameters to represent rotation transformation₀+p₁i+p₂j+p₃c, wherein p₀、p₁、p₂、p₃And i, j and c are real numbers and generalized imaginary units, and satisfy the following conditions:

i²＝j²＝c²＝ijc＝-1

ij ═ ji ═ c, jc ═ cj ═ i, Ci ═ iC ═ j formula 12

Then p is called a quaternion, if | p₁1, the method is a unit quaternion array which represents the large inclination angle rotation transformation in the three-dimensional space;

first, the unit quaternion array is converted to a rotation matrix:

let unit quaternion p be p₀+p₁i+p₂j+p₃c, obtaining a rotation matrix:

as long as the quaternion array is a unit quaternion array, the unit orthogonality of the rotation matrix can be ensured without supplementing an additional constraint equation;

second, the rotation matrix is converted to a unit quaternion array:

quaternion groups p and-p represent the same rotation matrix;

thirdly, the euler angle is converted into a unit quaternion array:

fourth, the unit quaternion array is converted to euler angles:

the euler angle formed by the quaternion array is characterized in that the quaternion array is expressed as a rotation matrix firstly, then the rotation matrix is converted into the euler angle, l takes a value from-90 degrees to +90 degrees, m and n take values from-180 degrees to +180 degrees, and the conversion formula from the rotation matrix to the euler angle is as follows:

the registration method of the laser point cloud and the analog image with the characteristic line fused further comprises the steps of analyzing registration parameters and setting a unit quaternion array p as p₀+p₁i+p₂j+p₃And c, expressing the rotation matrix by using a unit quaternion array as follows:

recording:

for s sets of line pairs, 2s constraint equations can be composed:

in the above formula d_ij、e_ijIs q is_wi、u_wiK is a coefficient matrix of 2s rows and 13 columns, K is a column vector consisting of 13 unknowns, | p | ═ 1, and the first 10 elements of K are not independent from each other;

firstly, elements in an unknown number vector K are used as independent parameters, an equation has 13 unknown numbers, the invention adopts a singular value decomposition method to solve a basic solution system, and the formula is as follows:

N＝VHU^T

in the formula 19, H is a singular value of N, and is also N^TN and NN^TThe square root of the non-zero eigenvalues, constituting a diagonal matrix of order 13, V, U being orthogonal matrices of order 2s and 13, respectively, each column of U being N^TN orthogonal unit feature vector, each column of V is NN^TOrthogonal unit feature vector of, N^TThe characteristic corresponding to the minimum characteristic value of NThe vector is the least square solution of 13 parameters of K, elements on the diagonal line of H are arranged from large to small, and the last column of U corresponding to the minimum singular value is the 13 parameters of K;

for equation 18, the general solution is represented by its base solution series as follows:

u_i＝[λ_i1λ_i2… λ_i10]^Tformula 20

m、u_i、w_iRespectively, the number of base solutions, the ith base solution and its corresponding coefficients, if equation 18 has S independent equations, then m is 13-S, since K is₁₃x₁The first 10 elements of (A) are rotation parameters, so K is calculated first₁₃x₁The expression is as follows:

wherein K₁₃x₁(1:10) is K₁₃x₁The first 10 elements of (a), u_i(1:10) As the first 10 elements of each basic solution, u_iSolving 10 coefficients when the calculation is completed; the first 10 elements in the K are not mutually independent, the first 10 elements are correlated, a constraint equation among 10 coefficients is established, the first 10 elements of the K are obtained after 10 coefficients are solved through a singular value decomposition method, namely a quaternion array representation method of a rotation matrix B, the translation vector G is solved by substituting the equation, and the common external orientation elements in the photogrammetry are obtained through conversion.

The invention discloses a registration method of a laser point cloud and a simulation image with fused characteristic lines, and further discloses a registration accuracy qualitative and quantitative evaluation method, which comprises the following steps:

firstly, a point cloud projection effect is achieved, laser point cloud is projected onto an image plane based on camera interior orientation elements and resolved image exterior orientation elements, and the registration effect of the laser point cloud and an image is qualitatively determined by observing the registration effect of point cloud outlines, obvious point, line and surface features and corresponding features on the image;

secondly, projecting the single image characteristic line to a plane based on the inside orientation element and the resolved outside orientation element of the image, calculating the average value of the vertical distances from two end points of the projection straight line segment of the point cloud characteristic line on the image to the corresponding image characteristic line, and calculating the included angle between the projection straight line segment and the corresponding image characteristic line;

and thirdly, based on the stereo measurement of two registration images, after the two images are registered with the ground laser point cloud, the three-dimensional coordinates of a target point are measured through front intersection and compared with the homonymous point in the initial laser point cloud, and the absolute measurement accuracy of the registration method of the invention is calculated by measuring the error in the point location.

The registration method of the laser point cloud and the analog image with the fused characteristic line further comprises the steps of selecting one or a plurality of simple surface patches as seed areas in point cloud rough cutting based on area growth, judging the similarity of point clouds in the neighborhood points and the seed areas on differential geometric attributes according to the local differential geometric attributes of the point clouds, and continuously supplementing the neighborhood points with certain similarity into the area to realize the continuous growth and expansion of the seed areas; in the process, parameters of the fitting plane or curved surface are continuously updated until no point meeting the normal vector included angle or distance critical value exists in the neighborhood range of the plane or curved surface, the cutting effect and the cutting precision mainly depend on the position selection, the similarity measurement and the growth sequence of the seed point or surface, and the method specifically comprises the following steps:

firstly, a local plane formed by fitting seed points and a point set in a certain range of a neighborhood of the seed points is used as a seed surface, a plurality of adjacent points of the seed surface to be tested are determined to determine whether the adjacent points can form a plane or not, if each parameter of a certain fitting plane can meet a preset critical value, the certain fitting plane is set as the seed surface, otherwise, other points are continuously used for searching the seed surface;

secondly, the similarity measure includes that the included angle between the normal vector of the candidate neighborhood point and the normal vector of the seed surface is smaller than a set critical value, and the distance between the candidate neighborhood point and a cut point or a cut plane is smaller than a set critical valueSetting a critical value; the altitude difference and local gradient information of the airborne laser point cloud can assist in cutting; the method specifically comprises the following steps: setting the point to be grown and the region in the range of neighborhood radius H of the growth surface boundary point, and setting the point in the radius as a candidate growth point, if the candidate growth point reaches the initial seed surface S₀If the vertical distance of the point is less than H, the point is brought into a seed point set of the growth surface; if the included angle between the normal vector of the candidate growing point and the normal vector of the growing surface accords with the critical value of the included angle, absorbing the point as a growing point;

thirdly, the growing sequence adopts a gradual cutting mode from local to whole for point cloud data with complex distribution or irregular shape, and after determining seed points or surfaces, local to whole expanding cutting is realized by searching neighborhood points meeting the geometric critical value condition in the neighborhood.

The registration method of the laser point cloud and the analog image with the fused characteristic line further comprises the steps that the point cloud fine cutting of the plane fitting adopts a random sampling consistency algorithm to estimate parameters of a plane model, the point cloud is subjected to fine cutting based on the random sampling consistency algorithm plane fitting, and the process of extracting high-precision plane parameters comprises two aspects of content, namely plane model fitting and point cloud plane parameter estimation based on the random sampling consistency algorithm; the random sampling consistency algorithm carries out multiple random sampling on sample data, the minimum sample is randomly taken out each time to calculate initial model parameters, all data are divided into inner points within a certain error range of the model and outer points outside the error range according to the model parameters, most data are divided into the inner points, after multiple random sampling, the maximum number of outer point sets are determined within an allowed error range, and the optimal model parameters are obtained.

The method for estimating the model parameters of the plane by the random sampling consistency algorithm comprises the following three steps:

step one, randomly selecting s (s >3) points from sample data with M points, and estimating model parameters of an initial fitting plane;

secondly, counting a data point inner point set S which accords with the geometric distance critical value condition and the plane model in the sample according to the geometric distance critical value parameter_iAnd make statistics ofThe ratio m of the inner points conforming to the model is S_iThe geometric distance critical value comprises a normal included angle critical value and a distance critical value, wherein the normal included angle critical value is a coplanarity condition, and the distance critical value is a coplanarity condition;

thirdly, calculating the probability q that the model parameter is correct after c times of sampling as 1- (1-m)^s)^cJudging whether the Q value is larger than a preset probability critical value Q for the correctness of the model parameters; if so, based on the current S_iCalculating final model parameters; if not, returning unmarked points and repeating the cycle until the optimal fitting parameters are obtained.

The registration method of the laser point cloud and the analog image with the fused characteristic line further comprises the following steps of judging the adjacency relation of two point cloud planes, wherein three critical value parameters are mainly adopted:

the first is the minimum distance between the point clouds of the two planes;

the minimum distance from a point in the plane to be analyzed to the reference plane and the minimum distance from a point in the reference plane to the plane to be analyzed;

thirdly, an included angle between the normal vectors of the two planes;

the adjacent line obtained by intersection operation directly based on the equations of two non-parallel adjacent planes is expressed by the coordinates of a certain point on the adjacent straight line and the normalized direction vector of the straight line; the invention accurately extracts the straight line segment characteristics intersected by the house roof surface, also needs to obtain the end point coordinates of two ends of the line characteristics, obtains two groups of collinear point sets by projecting all points of two adjacent surfaces onto the adjacent straight lines, and then solves the intersection of the two groups of point sets, namely the point sets on the adjacent lines of the two planes, thereby obtaining the coordinates of the two end points of the adjacent line segments.

The registration method of the laser point cloud and the analog image with characteristic line fusion further comprises the step of setting the function relation of x and y as d in the estimation of linear parameters₀+d₁x，d₀Represents intercept, d₁Representing the slope, wherein the relation has two undetermined parameters, and at least 2 points are needed for solving, and the method adopts a random sampling consistency algorithm method to estimate the parameters of a straight line; random sampling consistency algorithm method selects a random group in input dataThe sub-set achieves the goal, assuming the selected sub-set as a sample point, the parameters are estimated as follows:

step 1, calculating unknown parameters d by using the selected sample points₀,d₁Obtaining initial straight line parameters;

step 2, testing all data according to the initial model obtained in the step 1, and adding a point into the internal point set if the point is suitable for the estimation model;

step 3, if the number of points in the inner point set is enough, namely the number of points is larger than a critical value, the estimated initial model is considered to be reasonable, and the linear parameters are updated based on the updated inner point set;

step 4, solving optimal parameters by taking the fitting residual errors as evaluation criteria; the execution is repeated until a fixed number of times is reached, and the model parameters obtained each time are discarded with too few points in the sample point set or are approved better than the previous model.

The registration method of the laser point cloud and the analog image with the fused characteristic lines further comprises the steps that a merging method of line segment attributes such as the direction and the relative length of a line segment to be merged is adopted for merging the short line segment, the direction and the position of the line segment in which the merged line segment is located are calculated by taking the direction and the relative length of the line segment to be merged as weight parameters, and a good merging result is obtained; the method comprises the following two steps:

step 1, calculating the position and the direction of a straight line where the merged straight line segment is located; assuming that the point Q is a point on the merged straight line segment, firstly considering the influence of the lengths of the two straight line segments to be merged on the position of the point Q, and defining the point coordinates as follows:

wherein a, b, c and d are two end points of the straight line segments i and j respectively, and k_i、k_jIs the length of two straight line segments (a)_x,a_y) The coordinates of the point a and the coordinates of other 3 points are the same, and then the direction of two straight line segments to be merged is considered to be opposite to the direction of the straight line where the merged straight line segment is positionedThe influence of (2) on the direction of the straight line segment after merging is defined as follows:

wherein t represents the direction of the merged straight line, k_i、t_iRespectively representing the length and the direction of the line segment i;

step 2, determining two end points of the merged straight-line segment; after the position and the direction of the straight line where the merging straight line segment is located are determined, two end points of the straight line segment are determined, a, b, c and d are respectively projected onto a straight line k, and the straight line segment corresponding to the maximum range of the projection points is taken as a merging result.

Compared with the prior art, the invention has the following contributions and innovation points:

the method for registering the laser point cloud and the analog image by fusing the characteristic lines, provided by the invention, aims at the problems that the registration method of the image and the point cloud based on the point cloud distance intensity image or the point characteristic in the prior art has many defects of manual point selection error, point cloud data interpolation error, need of numerous observed values and the like, and the precision is very limited. The method for registering the laser point cloud and the simulated image based on the two-dimensional and three-dimensional fusion characteristic lines is provided, a direct conversion model between the image characteristic lines and the point cloud characteristic lines is established, the method has the advantages of no model error, no need of an external orientation element initial value, only need of a small number of registration elements and the like, the calculation speed is high, and the registration requirement of the simulated image and the ground laser point cloud data is met.

The point cloud characteristic line extraction adopts a characteristic line extraction method based on point cloud cutting plane intersection, KD-Tree management point cloud is adopted, normal vector and curvature are estimated by using a covariance analysis method, a region growing rough cutting method and a random sampling consistency plane fitting fine cutting method are adopted to realize point cloud plane cutting, and finally high-precision straight line characteristics are obtained by solving intersection of adjacent point cloud planes.

The method comprises the steps of extracting the characteristic line of the simulated image, searching the same type of pixel area based on the gradient tension value and the gradient direction by adopting a characteristic line extraction algorithm based on the gradient direction of the image, determining the support area of a straight line section by adopting a rectangular structure approximation method, and finally extracting the parameters of the straight line section based on the center, the length and the width of the rectangular structure. And providing a short straight line section combination judgment condition based on the straight line included angle and the distance and a combination processing method comprehensively considering the direction and the relative length of the straight line section. The method has the advantages of good noise immunity and robustness, high positioning precision of the characteristic line, small error of the registration parameter analyzed by the registration algorithm, high registration precision and capability of meeting the requirements of practical application.

And fourthly, according to the registration method of the laser point cloud and the analog image with the fused characteristic lines, the linear segment detection algorithm has the advantages of real-time performance, accuracy and robustness, high calculation efficiency, no need of setting parameters excessively, linear time complexity, capability of controlling false straight lines and the like. Feature lines are extracted from the aerial image based on a straight-line segment detection algorithm, parameters such as the length, the angle, the distance between endpoints, the vertical distance between line segments and the like are combined for optimization, and finally semi-automatic and high-precision registration of the aerial image of the urban area and the airborne laser radar is achieved.

The registration method of the laser point cloud and the simulated image with the fused characteristic lines has the advantages that the constraint effect of registration elements is strong, the required minimum registration element number is small, a strict registration model between two-dimensional and three-dimensional fusion registration elements is established, the requirement on the initial value of the external orientation element of the image is not strict, the calculated time cost is in direct proportion to the observed value number, the robustness on noise is strong, the automatic matching of two-dimensional and three-dimensional fusion characteristics is realized, the automation degree is improved, and the higher real-time performance is achieved. The method has very important value in the aspects of scientific research and application and social economic development requirements, has small limitation, universality and transportability, has huge popularization and application potential, and has wide market application prospect.

Drawings

FIG. 1 is a schematic flow chart of the registration method of laser point cloud and simulated image with feature line fusion according to the present invention.

Fig. 2 is a schematic diagram of a camera calibration photographing method of the invention.

FIG. 3 is a schematic diagram of extracting feature lines based on plane intersection of point cloud cutting.

FIG. 4 is a flow chart of the point cloud plane cutting of the region growing and RANSAC model fitting of the present invention.

FIG. 5 is a schematic diagram of the invention of region growing based point cloud segmentation.

FIG. 6 is a flow chart of the invention for extracting feature lines by intersecting point cloud planes.

Fig. 7 is a schematic view of the merging method of the present invention considering the direction and relative length of the straight line segment.

Detailed Description

The technical solution of the registration method of the laser point cloud with fused feature lines and the simulated image provided by the present invention is further described below with reference to the accompanying drawings, so that those skilled in the art can better understand the present invention and can implement the present invention.

The registration method of the laser point cloud and the simulated image with the fused characteristic lines, provided by the invention, has the advantages of no model error, no need of an external orientation element initial value, only need of a small amount of registration elements and the like by establishing a direct conversion model between the image characteristic lines and the point cloud characteristic lines, and is high in calculation speed, and the registration requirement of the simulated image and the ground laser point cloud data is met. The two-dimensional and three-dimensional fusion registration method based on the characteristic lines comprises a camera imaging model, camera calibration and a position and posture estimation method of a plurality of groups of fusion characteristic lines, and the position and posture estimation method of the plurality of groups of fusion characteristic lines comprises a characteristic line fusion registration function model, a quaternion array expression method of a rotation matrix, an analytic registration parameter and a registration precision judgment method, as shown in figure 1.

Extracting laser point cloud characteristic lines by adopting a characteristic line extraction method based on point cloud cutting plane intersection, managing point cloud by adopting KD-Tree, estimating normal vector and curvature by using a covariance analysis method, realizing point cloud plane cutting by adopting a region growing rough cutting method and a random sampling consistency plane fitting fine cutting method, and finally solving intersection of adjacent point cloud planes to extract high-precision straight line characteristics, wherein the characteristic lines have higher positioning precision and better noise immunity and robustness; the extraction of the laser point cloud characteristic line comprises point cloud preprocessing, point cloud plane cutting and intersection of the point cloud plane to extract the characteristic line, wherein the point cloud preprocessing comprises point cloud thinning and density homogenization and point cloud denoising, the point cloud plane cutting comprises point cloud rough cutting and plane fitting based on region growing, and the intersection of the point cloud plane to extract the characteristic line comprises estimation of adjacent straight line parameters, judgment of plane adjacent relation and solution of adjacent straight line segments.

Extracting the characteristic line of the simulated image by adopting a characteristic line extraction algorithm based on the gradient direction of the image, searching the similar pixel region based on the gradient tension value and the gradient direction, determining the support region of the straight-line segment by adopting a rectangular structure approximation method, and finally extracting the parameters of the straight-line segment based on the center, the length and the width of the rectangular structure. And providing a short straight line section combination judgment condition based on the straight line included angle and the distance and a combination processing method comprehensively considering the direction and the relative length of the straight line section. Experimental results prove that the method has better noise immunity and robustness, high positioning precision on the characteristic line, small error of the registration parameters analyzed by the registration algorithm, high registration precision and capability of meeting the requirements of practical application. The method for extracting the characteristic line of the simulated image comprises a straight line segment detection algorithm, straight line parameter estimation, short straight line segment merging processing and characteristic line screening.

Feature line-based two-dimensional and three-dimensional fusion registration method

Camera imaging model and camera calibration

The ground laser point cloud and the analog image are registered by using internal orientation elements of a camera, and in order to improve the registration precision, the analog image is subjected to optical distortion correction, a camera imaging model comprises a pinhole imaging model and an optical distortion model, the pinhole imaging model expresses the transformation relation between an object space three-dimensional point and an ideal image point coordinate, and the optical distortion model expresses the difference between the ideal image point coordinate and an actual image point coordinate.

The camera calibration analysis method comprises the following steps of calibrating and analyzing a space conversion relation among a camera coordinate system, an object coordinate system and an image coordinate system of a camera, wherein a single camera calibration parameter comprises the following steps: the first is internal parameter: 3 internal orientation elements, namely a principal distance and a principal point; secondly, distortion parameters: a radial distortion parameter and a tangential distortion parameter; and thirdly, external parameters: 6 external orientation elements, wherein the internal parameters of the solving camera are internal calibration, the external orientation parameters of the solving camera are external calibration, and the internal parameters and the distortion parameters are calibrated mainly.

The camera calibration adopts a camera calibration method of multi-chip space back intersection of a three-dimensional control field, the three-dimensional control field is composed of a plurality of three-dimensional control points with known space coordinates, the coordinates of the three-dimensional control points are accurately measured by angle intersection by using a theodolite or a total station and are obtained by adjustment, and the control points have depth and reliable positions in space distribution and do not deform.

The calibration method adopted by the invention comprises the following steps: selecting a plurality of positions in front of the control field, and acquiring three-dimensional control field images from different angle directions by using a single camera to ensure that the images shot by the same camera at different positions have enough overlapping degree, as shown in FIG. 2; the camera calibration method for multi-slice space back intersection of the three-dimensional control field takes the coordinates of image points as an observed value U according to a collinear condition equation to solve an inner orientation element J of a photo_{Inner part}And distortion addition parameter J_adWhile solving the element J of exterior orientation of the photo_{Outer cover}K is the corresponding mismatch value, and the formula is as follows:

U＝DJ_{outer cover}+EJ_{Inner part}+FJ_ad-K formula 1

D, E, F is coefficient matrix, and the unknown parameters J are obtained by using least square principle and successive iteration_1-outer、J_2-outer、J_3-outer、J_{Inner part}And J_adThe initial value of the optical distortion parameter of the camera is set to be zero, and the initial values of the internal and external orientation elements are obtained by a three-dimensional direct linear transformation method.

According to the invention, if the calibration point is selected on the two images, the three-dimensional coordinate of the calibration point is obtained by carrying out forward intersection on the basis of the calibrated internal parameter and external parameter, and the three-dimensional coordinate is compared with the known coordinate of the calibration point to judge the calibration precision of the camera.

Position and attitude estimation method for (II) multiple groups of fusion characteristic lines

The position and attitude estimation method of the multiple groups of fusion characteristic lines extracts the characteristic lines from the simulated image and the geodetic three-dimensional laser point cloud respectively, analyzes the position and attitude parameters of the simulated image under the scanner coordinate system, and realizes the registration of the two.

1. Feature line fusion registration function model

The imaging model of an ideal optical camera is:

q_pi＝C_3×3q_ci

q_ci＝B_3×3q_wi+G_3×1

G_3×1＝[g_xg_yg_z]^Tformula 2

The above formula is a functional model of a two-dimensional and three-dimensional fusion registration method based on point features, wherein q is_piCoordinates of the image point in the image coordinate system, C₃x₃Is an internal parameter matrix of the camera, corresponding to internal orientation elements in photogrammetry, which is obtained by calibrating the camera with optical distortion parameters, B₃x₃、G₃x₁The three-dimensional coordinate system is a rotation matrix and a translation matrix from an object coordinate system to a camera coordinate system, corresponds to an external orientation element in photogrammetry, and is a parameter to be solved for two-dimensional and three-dimensional fusion registration.

The invention analyzes the registration parameter B by utilizing the geodetic three-dimensional laser point cloud and the characteristic line on the image₃x₃、G₃x₁And establishing a two-dimensional and three-dimensional fusion function model based on the characteristic line.

The method for describing the straight line comprises the following steps:

firstly, the parameter of any straight line under the three-dimensional rectangular coordinate system is K ═ q_i，u_iWherein q is_iIs the three-dimensional coordinate of any point on the straight line, v_iIs the unit direction vector of the line, K_iAll aboveThe coordinates of a point are

Wherein a is any real number;

Thirdly, the parameter of a certain straight line under the camera coordinate system is K_ci＝{q_ci，u_ci}，q_ci＝[x_iy_if_i]^TI.e. the equation of the line is d_ix+e_iy+j_iz+h_i＝0。

For any group of homonymous two-dimensional and three-dimensional fusion straight lines K_ci、K_wiAnd, object space straight line K_wiWarp B₃x₃、G₃x₁And converting the transformed image into a camera coordinate system to obtain a straight line K_{ci_}＝{q_{ci_}，u_{ci_}Due to K_wiAnd K_ciIs a straight line of the same name, according to the camera imaging principle, K_{ci_}With a straight line K under the camera coordinate system_ciAnd the projection center S of the camera is coplanar, and the plane is recorded as N_i，N_iThe normal vector of (a) is recorded as M_iThe coplanarity condition is a basic condition of two-dimensional and three-dimensional fusion registration based on the characteristic line, and specifically comprises two constraint conditions: one is K_wiDirection vector u of_wiThrough the rotation matrix B₃x₃Direction vector u under camera coordinate system obtained after transformation_ciA and M_iVertically; second, K_wiAny point p on_wiThrough B₃x₃、G₃x₁Point p in camera coordinate system obtained after transformation_{ci_}To plane N_iIs 0, i.e. satisfies the plane N_iThe equation of (c); the two constraints are expressed by the equation:

and because:

then the following results are obtained:

the above equation is a function model of a two-dimensional and three-dimensional fusion registration algorithm based on characteristic lines.

Note plane N_iNormal vector of (1)

Plane N_iThe equation in the camera coordinate system is:

M_i1x+M_i2y+M_i3z+M_i40-0 formula 6

Plane N_iNormal vector M of_iBased on camera principal distance f and K_ci＝{q_ci，u_ciGet out due to

(y_ij_i-e_if)x+(d_if-x_ij_i)y+(x_ie_i-d_iy_i) z is 0 or 9

Combining formula 4 with formula 8 yields:

namely:

2. Quaternion array expression method of rotation matrix

The photogrammetry adopts three independent angle parameters to represent rotation transformation, the invention adopts a quaternion array to represent a rotation matrix, and p is equal to p₀+p₁i+p₂j+p₃c, wherein p₀、p₁、p₂、p₃And i, j and c are real numbers and generalized imaginary units, and satisfy the following conditions:

i²＝j²＝c²＝ijc＝-1

ij ═ ji ═ c, jc ═ cj ═ i, ci ═ ic ═ j formula 12

Then p is called a quaternion, if | p₁And 1, the method is a unit quaternion group which represents the large-inclination-angle rotation transformation in the three-dimensional space.

First, the unit quaternion array is converted to a rotation matrix:

let unit quaternion p be p₀+p₁i+p₂j+p₃c, obtaining a rotation matrix:

as long as the quaternion array is a unit quaternion array, the unit orthogonality of the rotation matrix can be ensured, and an additional constraint equation does not need to be supplemented.

Second, the rotation matrix is converted to a unit quaternion array:

the quaternion groups p and-p represent the same rotation matrix.

Thirdly, the euler angle is converted into a unit quaternion array:

fourth, the unit quaternion array is converted to euler angles:

3. resolving registration parameters

Let unit quaternion p be p₀+p₁i+p₂j+p₃And c, expressing the rotation matrix by using a unit quaternion array as follows:

recording:

for s sets of line pairs, 2s constraint equations can be composed:

in the above formula d_ij、e_ijIs q is_wi、u_wiK is a coefficient matrix of 2s rows and 13 columns, K is a column vector consisting of 13 unknowns, | p | ═ 1, and the first 10 elements of K are not independent of each other.

N＝VHU^T

in the formula 19, H is a singular value of N, and is also N^TN and NN^TThe square root of the non-zero eigenvalues, constituting a diagonal matrix of order 13, V, U being orthogonal matrices of order 2s and 13, respectively, each column of U being N^TN orthogonal unit feature vector, each column of V is NN^TOrthogonal unit feature vector of, N^TThe eigenvector corresponding to the minimum eigenvalue of N is the least square solution of 13 parameters of K, the elements on the diagonal line of H are arranged from large to small, and the last column of U corresponding to the minimum singular value is the 13 parameters of K.

u_i＝[λ_i1λ_i2… λ_i10]^Tformula 20

wherein K₁₃x₁(1:10) is K₁₃x₁The first 10 elements of (a), u_i(1:10) As the first 10 elements of each basic solution, u_iThe 10 coefficients have been found, and are solved for this. In KThe first 10 elements are not mutually independent and are correlated, a constraint equation among 10 coefficients is established, the first 10 elements of K, namely a quaternion array representation method of a rotation matrix B, are obtained after 10 coefficients are solved through a singular value decomposition method, the translation vector G is solved by substituting the equation, and the commonly used exterior orientation elements in photogrammetry are obtained through conversion.

4. Registration accuracy evaluation method

After the registration of the geodetic laser point cloud and the simulated image is completed and the exterior orientation element of the simulated image under the laser point cloud coordinate system is obtained, the invention adopts a qualitative and quantitative judgment method for the registration accuracy as follows:

the method comprises the steps of firstly, performing point cloud projection effect, projecting laser point cloud onto an image plane based on camera interior orientation elements and resolved image exterior orientation elements, and qualitatively judging the registration effect of the laser point cloud and an image by observing the registration effect of point cloud outlines, obvious point, line and surface features and corresponding features on the image.

Secondly, projecting the characteristic lines of the point clouds onto a plane based on the inside orientation elements and the resolved outside orientation elements of the image by the projection effect, the distance and the angle errors of the characteristic lines of the single image, and calculating the average value of the vertical distances from two end points of the projection straight line segments of the characteristic lines of the point clouds on the image to the characteristic lines of the corresponding image and the included angle between the projection straight line segments and the characteristic lines of the corresponding image.

The two-dimensional and three-dimensional fusion registration method based on the characteristic line has the characteristics and advantages that: the initial value of the exterior orientation element is not needed; the rotation matrix is expressed based on the quaternion array, so that many defects of expressing the rotation matrix by using the Euler angle are avoided; the large-angle rotation can be represented without linearization and iterative operation; the method has stronger robustness to noise and very high registration precision; the calculation time and the number of the observed values are in a linear relation, and the calculation efficiency is high.

Secondly, extracting laser point cloud characteristic lines

The invention extracts characteristic lines by adopting plane intersection cut based on point cloud, and the contour lines of main ground objects of the city are generally formed by intersecting adjacent planes, such as k shown by a solid line in figure 3₁₂、k₂₃、k₁₃However, based on the discreteness of the scanning point cloud data and the limitation of the scanning density, the scanning points may not be located on all the intersecting lines, and a straight line fitted by a ground object facade point cloud contour point set obtained directly based on a point cloud cutting algorithm is deviated to the interior of a building, as shown by a dotted line in fig. 3, and the outline of the ground object cannot be truly reflected. In order to extract accurate contour lines and other local feature lines of the ground objects, a point cloud cutting algorithm is firstly adopted to extract a plurality of local planes of the surface of the ground objects, such as a plane G in figure 3₁、G₂、G₃Then, the precise contour line, ridge line or local feature line is extracted by the way that the planes intersect each other pairwise, as shown by the broken line K in FIG. 3₁₂. The method has obvious advantages, can extract high-precision characteristic lines, and adopts the method to extract the characteristic lines in the point cloud.

Point cloud preprocessing

When the earth three-dimensional laser point cloud is collected, various random noises exist in the point cloud data due to the influence of instrument measurement errors and scanning environments. The invention mainly carries out preprocessing on point cloud data from two aspects, improves the data quality and improves the precision and efficiency of point cloud cutting and characteristic line extraction.

1. Point cloud thinning and density homogenization

According to the method, the density of the point cloud is homogenized based on a space grid method, the density of the point cloud in each area is enabled to be close to uniform, the global applicability of a critical value parameter of a subsequent cutting algorithm is enhanced, point cloud thinning is achieved, and the calculation efficiency is improved. The space grid method is characterized in that a regular grid structure is constructed in a three-dimensional space based on point cloud data, for each cubic grid, a point closest to the center of the cubic grid is reserved or the gravity center of a point set in the grid is used for approximating all points in the grid, point cloud sparse extraction and density consistency are achieved, key parameters of the method are the size of the cubic grid in X, Y, Z three directions, if the size set value is large, output point clouds are few, the subsequent point cloud cutting time can be shortened, and meanwhile, the characteristic extraction effect and accuracy can be adversely affected.

2. Point cloud denoising

The invention mainly uses two point cloud denoising methods: a Gaussian discrimination denoising method and a sphere filtering method; the Gaussian discrimination denoising method establishes a Gaussian distribution model for point cloud data, takes points with errors larger than a plurality of times of standard deviation as noise and deletes the noise, and the main parameters are as follows: the number N of the neighborhood search range and the standard deviation multiple N are determined, and when N is fixed, the larger N is, the fewer the number of points to be rejected is; when N is constant, the larger N is, the more the number of points to be removed is. The sphere filtering method is to draw a sphere based on a set radius by taking the current point as the center of the sphere, count the number of laser points in the sphere, and if the number is smaller than a set critical value, the point is considered as a free noise point and is deleted. The method comprises the following steps that main parameters comprise a sphere radius r and a critical value M of the number of laser points in a sphere, the critical value of the number of the laser points in the sphere is set based on the space density of point cloud, and when r is fixed, the larger M is, the more the number of points is eliminated; when M is fixed, the larger r is, the fewer the number of points to be removed is.

(II) cutting point cloud plane

Dense point clouds on the surface of a building, which are collected by a ground laser scanner, can reflect many details, but also collect a plurality of point clouds covering other ground objects, such as pedestrians, vehicles, roads, trees and the like, which have adverse effects on the extraction of the building characteristic line. Therefore, the building point cloud data and the point cloud of other ground objects need to be separated, so as to improve the extraction efficiency and the extraction precision of the surface features of the building. Data segmentation is one of the key techniques for data feature extraction.

The invention adopts a point cloud plane fitting method combining region growing and a random sampling consistency algorithm model, takes the normal vector and curvature of a seed point and the Euclidean distance from a neighborhood point to a normal plane as growth constraint conditions, extracts an initial sample point set conforming to the plane model by continuously expanding the seed plane, and then estimates high-precision plane parameters by adopting a steady parameter estimation method based on the random sampling consistency algorithm plane fitting, and the specific flow is as shown in figure 4:

1. rough cutting of point clouds based on region growing

Selecting one or a plurality of simple patches as a seed region, judging the similarity of point clouds in a neighborhood point and the seed region on the differential geometric attribute according to the local differential geometric attribute of the point clouds, and continuously supplementing the neighborhood point with certain similarity into the region to realize the continuous growth and expansion of the seed region; in the process, parameters of the fitting plane or curved surface are continuously updated until no point meeting the normal vector included angle or distance critical value exists in the neighborhood range of the plane or curved surface, as shown in fig. 5, the cutting effect and precision mainly depend on the position selection, similarity measurement and growth sequence of the seed point or surface, specifically:

secondly, similarity measurement comprises that an included angle between a normal vector of a candidate neighborhood point and a normal vector of the seed surface is smaller than a set critical value, and the distance between the candidate neighborhood point and a cut point or a cut plane is smaller than the set critical value; the altitude difference and local gradient information of the airborne laser point cloud can assist in cutting; the method specifically comprises the following steps: setting the point to be grown and the region in the range of neighborhood radius H of the growth surface boundary point, and setting the point in the radius as a candidate growth point, if the candidate growth point reaches the initial seed surface S₀If the vertical distance of the point is less than H, the point is brought into a seed point set of the growth surface; if the included angle between the normal vector of the candidate growing point and the normal vector of the growing surface accords with the critical value of the included angle, absorbing the point as a growing point;

2. Point cloud fine cutting of plane fitting

Under ideal conditions, the point cloud data with the form comparison rule can be directly solved for corresponding model parameters through model fitting, but due to the limitations of parameters of a point cloud denoising algorithm and a rough cutting algorithm, point cloud clusters obtained through rough cutting based on region growing still have wrong cutting points, the optimal model parameters need to be estimated from the point cloud data containing the wrong points, and precise cutting is realized.

The method adopts a random sampling consistency algorithm to estimate parameters of a plane model, performs fine cutting on the point cloud based on random sampling consistency algorithm plane fitting, and extracts high-precision plane parameters, wherein the process comprises two aspects of contents, namely plane model fitting and point cloud plane parameter estimation based on the random sampling consistency algorithm; the random sampling consistency algorithm carries out multiple random sampling on sample data, the minimum sample is randomly taken out each time to calculate initial model parameters, all data are divided into inner points within a certain error range of the model and outer points outside the error range according to the model parameters, most data are divided into the inner points, after multiple random sampling, the maximum number of outer point sets are determined within an allowed error range, and the optimal model parameters are obtained.

secondly, counting a data point inner point set S which accords with the geometric distance critical value condition and the plane model in the sample according to the geometric distance critical value parameter_iAnd counting the ratio m of the inner points conforming to the model as S_iThe geometric distance critical value comprises a normal included angle critical value and a distance critical value, wherein the normal included angle critical value is a coplanarity condition, and the distance critical value is a coplanarity condition;

Regarding the critical value: m is₀The ratio of prior inner points in the data set is obtained, the accurate value of the ratio is unknown, and the estimation is carried out according to the data quality; in the sampling of randomly selecting s points once, the error probability of the analyzed model parameter is

Assuming that sampling is needed at least c times, the probability that a sample containing s points and having no external points is extracted from the three-dimensional point set is ensured; the invention takes the critical value Q of Q as 0.99, and then

And solving a corresponding critical value C of the sampling times C.

(III) extracting characteristic line by intersecting point cloud plane

After the point cloud plane cutting and parameter estimation are realized, the straight line segment features in the point cloud are extracted by adopting the process shown in fig. 6.

1. Estimating adjacent straight line parameters

If two planes C in three-dimensional space₁、C₂Non-parallel, they must intersect and their intersection is a straight line, let the equation for these two planes be:

wherein

Are normal vectors of the two planes respectively, the direction vector of the intersection line of the two planes is

Let the point set on the intersection of the two planes satisfy the following linear equation:

simultaneous two formula yields:

solving to obtain:

substituting equation 23, the equation of the intersecting straight line is obtained:

wherein d, e and g are coefficients, respectively.

2. Determining plane adjacency relation and solving adjacent straight line segment

The adjacent relation of two point cloud planes is judged by mainly adopting three critical value parameters:

the first is the minimum distance between the point clouds of the two planes;

thirdly, an included angle between the normal vectors of the two planes;

the adjacent line obtained by intersection operation directly based on the equations of two non-parallel adjacent planes is expressed by the coordinates of a certain point on the adjacent straight line and the normalized direction vector of the straight line. The invention accurately extracts the straight line segment characteristics intersected by the house roof surface, also needs to obtain the end point coordinates of two ends of the line characteristics, obtains two groups of collinear point sets by projecting all points of two adjacent surfaces onto the adjacent straight lines, and then solves the intersection of the two groups of point sets, namely the point sets on the adjacent lines of the two planes, thereby obtaining the coordinates of the two end points of the adjacent line segments.

Thirdly, extracting the characteristic line of the simulation image

The analog image of the invention is a digital image acquired by a common digital camera, and a large number of characteristic lines exist in artificially constructed structures and main road ground objects.

Straight line segment detection algorithm

The method comprises the steps of carrying out edge extraction by utilizing a gradient direction to replace a gradient tension value, carrying out denoising treatment on an initial image by utilizing a Gaussian model, then calculating the gradient tension value and the gradient direction of each pixel, arranging all the pixels according to the sequence of the gradient tension values from large to small, starting from the pixel with the larger gradient tension value, dividing the pixels with the similarity of the gradient direction into pixel regions with the uniform gradient direction through iterative calculation, adopting a rectangular structure to approximate the pixel regions with the similar gradient direction, and finally extracting the corresponding straight-line segment characteristics of the regions based on the central line of the rectangular structure.

The main principle and steps of the straight-line segment detection algorithm are as follows:

estimating the gradient tension value and direction of a pixel; in order to reduce the influence of noise on the result, a Gaussian model is used for denoising the digital image before the gradient is calculated, the gradient of the digital image is obtained by calculating a derivative of a two-dimensional discrete function, and the differentiation in the x direction and the y direction is completed. In order to make the calculated value close to the true value as much as possible under the condition of random noise and reduce the adverse effect of the neighborhood pixels on the opening gradient value and the gradient direction, the invention calculates the opening gradient value and the gradient direction by differentiating 4 pixels in a 2 multiplied by 2 neighborhood around the pixel, and the mathematical expression of the first order differential of the discrete digital image is corresponding to the difference of two adjacent pixels;

searching areas with similar gradient directions; the digital image pixel gradient is a geometric attribute of the change size and direction of the pixel gray value, wherein the gradient size reflects the intensity of the change of the local pixel gray value of the image, specifically, the edge in the image has a sudden change of the pixel gray value, so that the pixel gradient value is larger, the change of the pixel gray value at the non-edge is basically consistent, and the pixel gradient value is smaller.

Firstly, storing gradient information of each pixel point in a matrix structure, then sequencing all pixels according to the sequence of gradient values of the pixels from large to small, taking the pixel with larger gradient tension as a seed point, searching 8 neighborhood pixels of the seed point through iterative calculation, searching pixel points with similarity in the gradient direction, classifying the pixel points and the seed point into a same type region and marking the same type region as traversed, calculating a gradient direction mean value of the pixel points in the current same type region after traversing all the pixel points in the 8 neighborhood of the seed point and supplementing the gradient similar pixel points to the current region, and continuously increasing the pixel points in the same type region and continuously updating the direction mean value along with the expansion of the search region.

Step three, extracting straight line segments based on rectangular structure approximation

After the search of the areas with similar gradient directions is completed, approximating each area by using a rectangular structure, and determining the position and the direction of a rectangle; specifically, for each region class, the coordinate of the center point of the rectangle is estimated by utilizing the weighted average of the gray value information of each pixel in the region class; the direction of the rectangle is the average value of the directions of the pixels in the region, then the size of the rectangle is adjusted under the condition that the coverage rate of the pixels belonging to the same region in the rectangle reaches 99%, the length and the width of the rectangle are calculated, the approximation of the rectangle is completed, and at the moment, the central line of the rectangle is the straight-line segment characteristic of the region.

(II) estimating straight line parameters

Let the functional relationship between x and y be d₀+d₁x，d₀Represents intercept, d₁The slope is represented, the relation has two undetermined parameters, at least 2 points are needed for solving, and the method adopts a random sampling consistency algorithm method to estimate the parameters of the straight line.

The random sampling consistency algorithm method achieves the goal by selecting a random set of subsets in the input data, assuming the selected subsets as sample points, and estimating the parameters according to the following method:

step 1, calculating unknown parameters by using the selected sample pointsNumber d₀,d₁Obtaining initial straight line parameters;

(III) merging treatment of short straight line segments

The method adopts a merging method of integrating the line segment attributes such as the direction, the relative length and the like of the line segment to be merged, calculates the direction and the position of the line segment after merging by taking the direction and the relative length of the line segment to be merged as weight parameters, and obtains a better merging result; the method comprises the following two steps:

wherein a, b, c and d are two end points of the straight line segments i and j respectively, and k_i、k_jIs the length of two straight line segments (a)_x,a_y) The coordinates of the point a and the coordinates of the other 3 points are the same, then the influence of the directions of the two straight line segments to be merged on the direction of the straight line where the merged straight line segment is located is considered, and the direction of the straight line where the merged straight line segment is located is defined as follows:

wherein t represents the direction of the merged straight line, k_i、t_iRespectively representing the length and direction of the line segment i.

Step 2, determining two end points of the merged straight-line segment; after the position and the direction of the straight line where the merging straight line segment is located are determined, two end points of the straight line segment are determined, as shown in fig. 7, a, b, c and d are respectively projected onto a straight line k, and the straight line segment (e and f corresponding to the straight line segment) corresponding to the maximum range of the projection point is taken as a merging result.

(IV) screening the characteristic line

In order to improve the resolution precision of the registration parameters, the image characteristic lines need to meet the requirements of proper length and uniform distribution on the image, and the quantity of the image characteristic lines needs to be slightly different from the quantity of the three-dimensional characteristic lines participating in registration. Therefore, after the linear features are extracted from the simulated image, the extracted linear features are screened. The specific method for screening the characteristic line comprises the following steps: uniformly dividing an image into a plurality of areas; counting the number and length of straight line segments in each region; setting different quantity and length condition critical values for each area in consideration of image texture and material difference in different areas, wherein the critical value conditions are relatively loose for the area with less quantity of characteristic lines, and the critical value conditions are strict for the area with more quantity of characteristic lines; and respectively eliminating overlong or overlong characteristic lines in each region according to set critical value conditions.

Claims

1. The feature line fused laser point cloud and analog image registration method is characterized in that a direct transformation model between an image feature line and a point cloud feature line is established, a two-dimensional and three-dimensional fusion registration method based on the feature line comprises a camera imaging model and a camera calibration and position and attitude estimation methods of a plurality of groups of fusion feature lines, and the position and attitude estimation methods of the plurality of groups of fusion feature lines comprise a feature line fusion registration function model, a quaternion array expression method of a rotation matrix, an analytic registration parameter and a registration precision judgment method;

2. The method of claim 1, wherein the position and orientation estimation method for multiple sets of fused feature lines extracts feature lines from the simulated image and the ground three-dimensional laser point cloud respectively, analyzes the position and orientation parameters of the simulated image under the scanner coordinate system, realizes the registration of the two, and the feature line fusion in the registration function model,

the imaging model of an ideal optical camera is:

q_pi＝C_3×3q_ci

q_ci＝B_3×3q_wi+G_3×1

G_3×1＝[g_xg_yg_z]^Tformula 2

The above formula is a functional model of a two-dimensional and three-dimensional fusion registration method based on point features, wherein q is_piCoordinates of the image point in the image coordinate system, C_3x3Is an internal parameter matrix of the camera, corresponding to internal orientation elements in photogrammetry, which is obtained by calibrating the camera with optical distortion parameters, B_3x3、G_3x1The method comprises the following steps that a rotation matrix and a translation matrix from an object space coordinate system to a camera coordinate system are respectively, correspond to external orientation elements in photogrammetry and are parameters to be solved for two-dimensional and three-dimensional fusion registration;

the invention analyzes the registration parameter B by utilizing the geodetic three-dimensional laser point cloud and the characteristic line on the image_3x3、G_3x1Establishing a two-dimensional and three-dimensional fusion function model based on the characteristic line; the method for describing the straight line comprises the following steps:

firstly, the parameter of any straight line under the three-dimensional rectangular coordinate system is K ═ q_i,u_iWherein q is_iIs the three-dimensional coordinate of any point on the straight line, v_iIs the unit direction vector of the line, K_iThe coordinate of any point is q'_i＝q_i+au_iWherein a is any real number;

For any group of homonymous two-dimensional and three-dimensional fusion straight lines K_ci、K_wiAnd, object space straight line K_wiWarp B_3x3、G_3x1And converting the transformed image into a camera coordinate system to obtain a straight line K_{ci_}＝{q_{ci_}，u_{ci_}Due to K_wiAnd K_ciIs a straight line of the same name, according to the camera imaging principle, K_{ci_}With a straight line K under the camera coordinate system_ciAnd the projection center S of the camera is coplanar, and the plane is recorded as N_i，N_iThe normal vector of (a) is recorded as M_iThe coplanarity condition is a basic condition of two-dimensional and three-dimensional fusion registration based on the characteristic line, and specifically comprises two constraint conditions: one is K_wiDirection vector u of_wiThrough the rotation matrix B_3x3Direction vector u under camera coordinate system obtained after transformation_{ci_}And M_iVertically; second, K_wiAny point p on_wiThrough B_3x3、G_3x1Point p in camera coordinate system obtained after transformation_{ci_}To plane N_iIs 0, i.e. satisfies the plane N_iThe equation of (c); the two constraints are expressed by the equation:

and because:

then the following results are obtained:

note plane N_iNormal vector of (1)

Plane N_iThe equation in the camera coordinate system is:

M_i1x+M_i2y+M_i3z+M_i40-0 formula 6

Plane N_iNormal vector M of_iBased on phasePrincipal distance f and K_ci＝{q_ci，u_ciFind out, because,

(y_ij_i-e_if)x+(d_if-x_ij_i)y+(x_ie_i-d_iy_i) z is 0 or 9

Combining formula 4 with formula 8 yields:

namely:

3. The method of claim 1, wherein the photogrammetry uses three independent angle parameters to represent the rotation transformation, the invention uses a quaternion array to represent the rotation matrix, and p ═ p represents the rotation matrix₀+p₁i+p₂j+p₃c, wherein p₀、p₁、p₂、p₃And i, j and c are real numbers and generalized imaginary units, and satisfy the following conditions:

i²＝j²＝c²＝ijc＝-1

ij ═ ji ═ c, jc ═ cj ═ i, ci ═ ic ═ j formula 12

first, the unit quaternion array is converted to a rotation matrix:

let unit quaternion p be p₀+p₁i+p₂j+p₃c, obtaining a rotation matrix:

second, the rotation matrix is converted to a unit quaternion array:

quaternion groups p and-p represent the same rotation matrix;

thirdly, the euler angle is converted into a unit quaternion array:

fourth, the unit quaternion array is converted to euler angles:

。

4. the method of claim 1, wherein the registration parameters are analyzed by setting a unit quaternion array p-p₀+p₁i+p₂j+p₃And c, expressing the rotation matrix by using a unit quaternion array as follows:

recording:

for s sets of line pairs, 2s constraint equations can be composed:

N＝VHU^T

in the formula 19, H is a singular value of N, and is also N^TN and NN^TThe square root of the non-zero eigenvalues, constituting a diagonal matrix of order 13, V, U being orthogonal matrices of order 2s and 13, respectively, each column of U being N^TN orthogonal unit feature vector, each column of V is NN^TOrthogonal unit feature vector of, N^TThe least square solution of 13 parameters with K as the characteristic vector corresponding to the minimum characteristic value of N is used for H diagonalThe elements are arranged from large to small, and the last column of U corresponding to the minimum singular value is 13 parameters of K;

u_i＝[λ_i1λ_i2… λ_i10]^Tformula 20

m、u_i、w_iRespectively, the number of base solutions, the ith base solution and its corresponding coefficients, if equation 18 has S independent equations, then m is 13-S, since K is_13x1The first 10 elements of (A) are rotation parameters, so K is calculated first_13x1The expression is as follows:

wherein K_13x1(1:10) is K_13x1The first 10 elements of (a), u_i(1:10) As the first 10 elements of each basic solution, u_iSolving 10 coefficients when the calculation is completed; the first 10 elements in the K are not mutually independent, the first 10 elements are correlated, a constraint equation among 10 coefficients is established, the first 10 elements of the K are obtained after 10 coefficients are solved through a singular value decomposition method, namely a quaternion array representation method of a rotation matrix B, the translation vector G is solved by substituting the equation, and the common external orientation elements in the photogrammetry are obtained through conversion.

5. The registration method of the laser point cloud and the simulation image fused with the characteristic line according to claim 1, which is characterized in that the registration accuracy qualitative and quantitative judgment method adopted by the invention is as follows:

6. The registration method of the laser point cloud and the analog image fused with the characteristic line according to claim 1, characterized in that in the point cloud rough cutting based on the region growth, one or a plurality of simple surface patches are selected as seed regions, according to the local differential geometric attributes of the point cloud, the similarity of the point cloud in the neighborhood and the seed region on the differential geometric attributes is judged, the neighborhood with certain similarity is continuously supplemented into the region, and the continuous growth and expansion of the seed region are realized; in the process, parameters of the fitting plane or curved surface are continuously updated until no point meeting the normal vector included angle or distance critical value exists in the neighborhood range of the plane or curved surface, the cutting effect and the cutting precision mainly depend on the position selection, the similarity measurement and the growth sequence of the seed point or surface, and the method specifically comprises the following steps:

secondly, the similarity measure includes that the included angle between the normal vector of the candidate neighborhood point and the normal vector of the seed surface is smaller than the set critical valueThe distance from the candidate neighborhood point to the cut point or the cut plane is smaller than a set critical value; the altitude difference and local gradient information of the airborne laser point cloud can assist in cutting; the method specifically comprises the following steps: setting the point to be grown and the region in the range of neighborhood radius H of the growth surface boundary point, and setting the point in the radius as a candidate growth point, if the candidate growth point reaches the initial seed surface S₀If the vertical distance of the point is less than H, the point is brought into a seed point set of the growth surface; if the included angle between the normal vector of the candidate growing point and the normal vector of the growing surface accords with the critical value of the included angle, absorbing the point as a growing point;

7. The registration method of the laser point cloud and the analog image fused with the characteristic line according to claim 1, wherein the point cloud fine cutting of the plane fitting adopts a random sampling consistency algorithm to estimate parameters of a plane model, the point cloud is finely cut based on the random sampling consistency algorithm plane fitting, and the process of extracting high-precision plane parameters comprises two aspects of contents, namely plane model fitting and point cloud plane parameter estimation based on the random sampling consistency algorithm; the random sampling consistency algorithm carries out multiple random sampling on sample data, the minimum sample is randomly taken out each time to calculate initial model parameters, all data are divided into inner points within a certain error range of the model and outer points outside the error range according to the model parameters, most data are divided into the inner points, after multiple random sampling, the maximum number of outer point sets are determined within an allowed error range, and the optimal model parameters are obtained;

in the second step, the first step is that,counting a data point inner point set S which accords with the geometric distance critical value condition and the plane model in the sample according to the geometric distance critical value parameter_iAnd counting the ratio m of the inner points conforming to the model as S_iThe geometric distance critical value comprises a normal included angle critical value and a distance critical value, wherein the normal included angle critical value is a coplanarity condition, and the distance critical value is a coplanarity condition;

8. The method of claim 1, wherein the determination of the adjacency between two point cloud planes is mainly performed using three threshold parameters:

the first is the minimum distance between the point clouds of the two planes;

thirdly, an included angle between the normal vectors of the two planes;

9. The method of claim 1, wherein the linear parameter is estimated by setting the function relationship between x and y as y-d₀+d₁x，d₀Represents intercept, d₁Representing the slope, wherein the relation has two undetermined parameters, and at least 2 points are needed for solving, and the method adopts a random sampling consistency algorithm method to estimate the parameters of a straight line; the random sampling consistency algorithm method achieves the goal by selecting a random set of subsets in the input data, assuming the selected subsets as sample points, and estimating the parameters according to the following method:

10. The registration method of the laser point cloud and the simulation image fused with the characteristic line according to claim 1, wherein the merging processing of the short straight line section adopts a merging method of integrating line section attributes such as the direction and the relative length of a straight line section to be merged, and the direction and the position of the straight line section after merging are calculated by taking the direction and the relative length of the line section to be merged as weight parameters, so that a better merging result is obtained; the method comprises the following two steps: