CN116645348A - High-speed loop linear fitting method based on point cloud random sampling - Google Patents
High-speed loop linear fitting method based on point cloud random sampling Download PDFInfo
- Publication number
- CN116645348A CN116645348A CN202310611822.2A CN202310611822A CN116645348A CN 116645348 A CN116645348 A CN 116645348A CN 202310611822 A CN202310611822 A CN 202310611822A CN 116645348 A CN116645348 A CN 116645348A
- Authority
- CN
- China
- Prior art keywords
- point
- representing
- curve
- points
- point cloud
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 61
- 238000005070 sampling Methods 0.000 title claims abstract description 28
- 238000001514 detection method Methods 0.000 claims abstract description 4
- 238000007781 pre-processing Methods 0.000 claims abstract description 4
- 238000004364 calculation method Methods 0.000 claims description 30
- 238000001914 filtration Methods 0.000 claims description 19
- 239000011159 matrix material Substances 0.000 claims description 18
- 238000005096 rolling process Methods 0.000 claims description 14
- 238000013461 design Methods 0.000 claims description 13
- 230000010354 integration Effects 0.000 claims description 9
- 238000012937 correction Methods 0.000 claims description 4
- 230000001133 acceleration Effects 0.000 claims description 3
- 238000004422 calculation algorithm Methods 0.000 description 7
- 238000000605 extraction Methods 0.000 description 7
- 238000005259 measurement Methods 0.000 description 5
- 230000011218 segmentation Effects 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 4
- 238000010276 construction Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 230000005484 gravity Effects 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 230000007547 defect Effects 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000018109 developmental process Effects 0.000 description 2
- 201000010099 disease Diseases 0.000 description 2
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000005192 partition Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 210000000988 bone and bone Anatomy 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 238000012850 discrimination method Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000013178 mathematical model Methods 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 230000002040 relaxant effect Effects 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 238000004901 spalling Methods 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 239000013589 supplement Substances 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/13—Edge detection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/70—Determining position or orientation of objects or cameras
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10028—Range image; Depth image; 3D point clouds
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
- G06T2207/10044—Radar image
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Quality & Reliability (AREA)
- Image Analysis (AREA)
Abstract
The invention discloses a high-speed loop linear fitting method based on point cloud random sampling, and relates to the technical field of loop repair. The method comprises the following steps: s1, preprocessing an original point cloud to obtain a Lu Miandian cloud; s2, carrying out boundary detection on the road point cloud by utilizing the alpha shapes section to obtain a boundary point cloud; s3, carrying out circular curve boundary fitting and circle center determination on the boundary point cloud based on random sampling consistency, and carrying out straight line segment boundary fitting on the boundary point cloud based on random sampling consistency to determine the size of a field boundary; s4, reversely pushing the length of the field relief curve and the length of the field straight line according to the field boundary size and the position of the circular curve. The invention realizes that the measuring precision and efficiency can be greatly improved and the degree of automation can be improved by extracting the loop line shape through the LiDAR point cloud data of the laser radar.
Description
Technical Field
The invention relates to the technical field of loop repair, in particular to a high-speed loop linear fitting method based on point cloud random sampling.
Background
From the construction of the high-speed loop in the eighties of the last century, the construction of the domestic high-speed loop has undergone a rapid development period, and the early high-speed loop has reached the life cycle successively at present, wherein the high-speed loop typically comprises a plurality of test sites which are used for more than 30 years, such as Hainan, constant distance, xiangyang and the like. The high-speed loop is affected by the use strength, pavement diseases occur to different degrees, and main disease types comprise crushing plates, plate corner fracture, corner spalling, cracks, exposed bones and the like, wherein the highest occurrence of the high-speed loop is achieved by the medium-speed lane. Under the premise that the middle and small repair cannot meet the use requirement of the loop, the development of the loop overhaul engineering is necessary.
Different from a newly built loop, the overhaul of the established loop is influenced by various factors including initial design indexes, construction process, post-construction settlement, road surface use conditions and the like, the design method is based on the basic design theory of the loop, the actual condition of the existing loop is combined, and the engineering design of the loop overhaul is developed by fully utilizing the new technology and the new method. By carrying out analysis on the total life cycle cost of the loop, a reference basis can be provided for future loop creation or reconstruction.
Before the overhaul engineering is developed, the prior line shape of the high-speed loop is required to be subjected to early engineering measurement, and the measurement requirement of the high-speed loop cannot be met by the traditional total station and level measuring method limited by the measurement precision and the curved surface operation condition.
In chinese patent CN2021112910934, an automatic system and method for repairing unevenness of highway road are provided, and the integrated design of motor vehicle is combined to improve repairing efficiency, so as to reduce the blockage or safety accident of highway caused by road surface repairing; obtaining two road surface models through two processing modes and solving root mean square values of the road surface models, and more accurately calculating the volume of the to-be-filled road surface to be repaired; the pavement evenness can be higher by filling the concave space for multiple times, the quality is better, meanwhile, the same concave space can be filled in a layered mode according to different concave degrees, and the repairing effect is further improved. However, the prior art cannot provide assistance for reconstructing a road surface with a large area, and the camera is arranged on a running motor vehicle and cannot provide powerful assistance for high-precision measurement.
Therefore, providing a high-speed loop line fitting method is a problem that needs to be solved by those skilled in the art. The three-dimensional laser scanning technology can effectively solve the problem, and the measuring precision and efficiency can be greatly improved by extracting the loop line shape through the laser radar LiDAR point cloud data in the topographic survey, and the degree of automation is improved.
Disclosure of Invention
In view of the above, the invention provides a high-speed loop linear fitting method based on point cloud random sampling, so as to achieve the purpose of measuring the road surface condition with high precision.
In order to achieve the above purpose, the present invention adopts the following technical scheme:
the high-speed loop linear fitting method based on the point cloud random sampling comprises the following specific steps:
s1, preprocessing an original point cloud to obtain a Lu Miandian cloud; extracting linear point clouds of a high-speed loop through a laser radar; the road point cloud comprises a straight line segment road point cloud and a curve segment Lu Miandian cloud;
s2, carrying out boundary detection on the road point cloud by utilizing the Alpha shapes section to obtain a boundary point cloud; the boundary point cloud comprises a straight line segment boundary point cloud and a curve segment boundary point cloud;
s3, carrying out circular curve boundary fitting and circle center determination on the boundary point cloud based on random sampling consistency, and carrying out straight line segment boundary fitting on the boundary point cloud based on random sampling consistency to determine the size of a field boundary;
S4, reversely pushing the length of the field relief curve and the length of the field straight line according to the field boundary size and the position of the circular curve.
Optionally, the specific step S1 includes:
s1.1, filtering each point in the original point cloud and generating a denoising point cloud; the method for filtering the original point cloud adopts statistical filtering;
s1.2, dividing and denoising point cloud based on combination of pavement normal and elevation mean square error, and generating a pavement point cloud.
Optionally, the specific step of S1.1 includes:
s1.1.1 denoising the original point cloud once, selecting any point in the original point cloud as a first point in each denoising, and adding the first point and N nearest to the first point 1 -1 point constitutes a set of points; wherein, after denoising all points in the original point cloud, denoising point cloud is generated, N 1 Is an integer greater than 2; the point set formula:
wherein: x represents a point set, X 1 Represents a first point nearest to the first point, x 2 Representing a second point closest to the first point,indicating the third point nearest to the first point, < ->Representing a first point; wherein the point set accords with normal distribution;
s1.1.2 for each ofPoint, calculate +.>Point-to-point integration divide +.>The specific formula of the distance of any point except the point is as follows:
wherein:representation- >Point-to-point integration divide +.>Distance of any point other than point, +.>Representing the +.>Any point other than;
s1.1.3 atPoint to except->Selecting the minimum kmean distances from the distances of any point to generate a distance set, wherein the distance set formula is as follows:
wherein:representing distance set,/->Representation->Except->Distance of nearest first point outside the point, +.>Representation->Except->Distance of nearest second point outside the point, +.>Representation->Except->The distance of the closest kmean point outside the point; wherein kmean is an integer greater than 1;
s1.1.4 the mean of the distance setThe average value calculation formula of the characteristic value of the distance set is as follows:
wherein:representing a distance set mean;
s1.1.5 judging outliers based on a distance mean formula and a variance mean formula, removing outliers from original cloud points to reserve effective points, and generating denoising point clouds; the distance average formula is:
wherein: mean represents the mean of the distance set means of all points in the point set;
the distance variance formula is:
wherein: stddev represents the variance of the distance set mean of all points in the point set and the distance set mean of all points in the point set;
If d i1 >When mean+C stddev, then x j Is an outlier; if it isWhen then x j Is the effective point, where the C value is the outlier of the point.
Optionally, the specific step S1.2 includes:
s1.2.1 optionally selecting one point in the point cloud as a query pointBy inquiring about the point->Nearest k 2 Solving k by adjacent points 2 A covariance matrix A created by each neighboring point is fitted with a feature vector and a feature value of the covariance matrix by a least square method to obtain a point cloud normal; each point in the pair point set +.>Its corresponding covariance matrix can be expressed as:
wherein: a represents covariance matrix, k 2 Representation pointsAdjacent points number, P i Represents a query point->Representing the three-dimensional centroid of a neighboring point, +.>J-th representing covariance matrix 2 A feature vector lambda j2 J-th representing covariance matrix 2 A characteristic value;
s1.2.2 when lambda j2 Feature vector corresponding to minimum feature valueThe normal vector of the calculation point is:
wherein:represents the point normal +_>Representing the normal x-direction component, +.>Representing the normal y-direction component, +.>Representing the normal z-direction component, T being the transposed symbol;
s1.2.3 calculates the gradient value of each point, and the gradient value calculation formula is as follows:
wherein: tan alpha 1 Representing a grade value;
s1.2.4 comparing the gradient value with a gradient threshold value through gradient, and determining a dividing point and a pavement point of a high-speed loop line; wherein, if tan alpha 1 ≤tanα 2 Then it is the dividing point, if tan alpha 1 >tanα 2 Then is the road surface point, tan alpha 2 Representing a grade threshold;
s1.2.5 a local cylinder is set in the denoising point cloud, and the heights of each point in the range of the local cylinder are respectively as follows:
wherein: z is Z G Represent high Cheng Ge set, z G1 Representing the elevation, z, of a first point within the partial cylinder G2 Representing the elevation of the second point within the partial cylinder,representing the Nth in a partial cylinder 2 Elevation of point N 2 Representing the number of points within the partial cylinder; wherein N is 2 Is an integer greater than 1;
s1.2.6 calculating the elevation mean value and the elevation standard deviation of each point in the range of the local cylinder, and judging the local cylinder as a partition point when the elevation standard deviation exceeds the elevation standard deviation limit value; the calculation formula of the elevation mean value is as follows:
wherein:represent elevation mean, i 3 Representing the ith in a partial cylinder 3 A dot;
the calculation formula of the elevation standard deviation is as follows:
wherein:representing the standard deviation of elevation; wherein if->Then each point in the range of the local cylinder is a foreign point, < >>Representing a preset elevation standard deviation limit value;
s1.2.7 the cloud of denoising points is divided by dividing points with the road surface points as road surface areas, and Lu Miandian cloud is obtained.
Optionally, the specific step S2 includes:
s2.1 selecting any point in the plane projection point cloud as a point P to be judged e (x e ，y e ) According to the preset rolling circle radius alpha e Searching distance P from plane projection point cloud e The point distance is less than 2 alpha e All points in the set are marked as a point set Q; the generation method of the plane projection point cloud comprises the steps of enabling all points z of the road surface point cloud to be e The coordinates are all set to the same constant;
s2.2 selecting any point in the point set Q as a point P e1 (x e1 ，y e1 ) According to point P e Point P e1 And rolling circle radius alpha e Calculating the center coordinates O of the first circle 1 (x o1 ，y o1 ) And the center coordinates O of the second circle 2 (x o2 ，y o2 ) The method comprises the steps of carrying out a first treatment on the surface of the Wherein the first circle intersects with the second circle and has two points of intersection, the radius alpha of the rolling circle e Is the center of the first circle and two intersecting pointsIs equal to the rolling circle radius alpha e Is the distance between the center of the second circle and two intersecting points, P e 、P e1 Two intersecting points of the first circle and the second circle respectively;
s2.3 dividing P for Point set Q e1 Each point outside the circle is calculated to the center O 1 And to O 2 When all points are at the distance of O 1 And to O 2 Are all greater than alpha e When the point P is marked e Terminating the point P as the contour point e Is calculated;
s2.4 dividing P in the Point set Q e1 Any point outside to the centre of a circle O 1 And to O 2 When the distance of (2) is less than or equal to alpha, then traversing all the points in the point set Q to rotate as P e1 Point, recalculate the center coordinates and repeat step S2.3 to point P e Judging; if there is one point is point P e If the mark is a contour point, the point P can be judged e Is the contour point, and the contour point set obtained is the boundary point set S, otherwise, the point P is judged e Is a non-contour point.
Optionally, the specific step S3 includes:
s3.1, fitting a point cloud boundary contour by adopting RANSAC;
s3.2, an output model obtained by fitting a circular curve of the point cloud boundary profile through S3.1 is as follows:
(x h -x ho ) 2 +(y h -y ho ) 2 ＝R h 2
wherein: (x) h ,y h ) Representing the coordinates of points on a circle curve, (x) ho ,y ho ) Represents the center coordinates of a circular curve, R h Representing the radius of the circular curve.
S3.3, an output model obtained by fitting the S3.1 to the straight line segment of the point cloud boundary contour is as follows:
k h (x h -x h1 )＝(y h -y h1 )
wherein: (x) h ,y h ) Representing point coordinates on a straight line, (x) h1 ,y h1 ) Representing the coordinates of a known point on a straight line, y h1 Representation, k h The slope of the straight line is shown.
Optionally, the specific step S3.1 includes:
s3.1.1 randomly selecting a subset S from the boundary point set S, taking the randomly selected subset S as assumed local points, and fitting the initial model by the local points S;
estimating model parameters according to the subset s, calculating the coordinates of the assumed local points as (x) by using the circular curve parameters β1 ,y β1 ),(x β2 ,y β2 ),(x β3 ,y β3 ) With S3.2 circular curve model
(x β1 -x ho ) 2 +(y β1 -y ho ) 2 ＝R h 2
(x β2 -x ho ) 2 +(y β2 -y ho ) 2 ＝R h 2
(x β3 -x ho ) 2 +(y β3 -y ho ) 2 ＝R h 2
Solving the equation set to obtain
R h 2 ＝(x β1 -x ho ) 2 +(y β1 -y ho ) 2
Straight line parameter calculation assuming that the local point coordinates are (x β1 ,y β1 ),(x β2 ,y β2 ) With S3.3 straight line model
k h (x β1 -x h1 )＝(y β1 -y h1 )
k h (x β2 -x h1 )＝(y β2 -y h1 )
Solving the equation set to obtain
k h ＝(y β2 -y β1 )/(x β2 -x β1 )
x h1 ＝x β1
y h1 ＝y β1
S3.1.2 traversing all data except subset S in boundary point set S and calculating radius R of circular curve h Or slope k of straight line h If the error is within the error limit e, marking as an inner point, otherwise marking as an outer point;
s3.1.3 forming a consistent set by all the interior points obtained in S3.1.2, re-estimating model parameters by using all the interior points in the consistent set if the number of the consistent set points meets a given threshold T, and outputting a result; if the number of the inner points in the consistent set is less than the threshold value T, a new subset s is reselected, and the steps S3.1.1-S3.1.2 are repeated;
s3.1.4 through multiple iterations of S3.1.1-S3.1.3, selecting a consistent set with the largest number of interior points, re-estimating model parameters by using all interior points in the consistent set, and outputting a result; wherein the number of iterations is an integer number of iterations greater than 1.
Optionally, the specific step S4 includes:
s4.1, calculating the overall field width of the high-speed loop and the overall field length of the high-speed loop through a drift distance formula, and then respectively obtaining a relaxation curve longitudinal drift distance and a relaxation curve transverse drift distance; the formula of the relation between the width and the length of the field and the offset distance of the relaxation curve is as follows:
wherein: w represents the whole field width of the high-speed loop, and k 1 The slope of the straight line is shown,the abscissa representing the midpoint of the center of the two sides, +.>An ordinate representing the center point of two sides, y 1 Representing the ordinate, x, of a known point on the first fitted line 1 An abscissa, y, representing a known point on the first fitted line 2 Representing the ordinate, x, of a known point on the second fitting line 2 Represents the abscissa of the known point on the second fitting straight line, M represents the lane width adjustment value between the straight line segment center line and the fitting straight line, beta represents the moderation curve angle, R represents the radius of the circular curve, L y Represents the longitudinal offset distance of a field relaxation curve, H represents the whole field length of a high-speed loop, and X o1 X represents the circular center abscissa of the first end point of the high-speed loop o2 Represents the abscissa of the center of a circle of a second endpoint of the high-speed loop, Y o1 Represents the ordinate of the circular center of the first end point of the high-speed loop, Y o2 Represents the ordinate of the center of a circle of a second endpoint of the high-speed loop, L represents the length of a straight line of a field, and L x Representing the ground relief curve lateral offset.
S4.2, calculating a moderation curve angle by a curvature quadratic integration method; the relaxation curve angle β is calculated using the following formula:
wherein: ls (Ls) 0 Representing the assumed length of the relaxation curve, k (l) representing the function of the curvature of any point on the relaxation curve as a function of the length of the relaxation curve;
the moderation curve lateral offset and moderation curve longitudinal offset are calculated by adopting the following steps:
S4.3 after linear fitting and correction, the length Ls of the assumed relaxation curve is 0 And calculating errors of the corresponding field width and the actual field width W, obtaining the final relaxation curve length Ls through an approximation method, and calculating the field straight line length L through the field length H.
Optionally, the relaxing curve curvature function k (l) in S4.2 adopts any one of a Bloss curve curvature function, a McConnell curve curvature function, a Grabowski curve curvature function, an Auberlen curve curvature function and a Klein curve curvature function;
the Bloss curve curvature function is:
wherein: k (l) Bloss Representing a Bloss curve curvature function;
the McConnell curve curvature function is:
wherein: k (l) McConnell Representing a McConnell curve curvature function, wherein ψ represents a maximum superhigh angle, v represents a design speed, and g represents a gravitational acceleration;
the Grabowski curve curvature function is:
wherein: k (l) Grabowski Representing a Grabowski curve curvature function;
the Auberlen curve curvature function is:
wherein: k (l) Auberlen Representing the Auberlen curve curvature function;
the Klein curve curvature function is:
wherein: k (l) Klein Representing the Klein curve curvature function.
Compared with the prior art, the invention discloses a high-speed loop linear fitting method based on point cloud random sampling, thereby obtaining the following beneficial effects:
1. The loop line shape can be extracted through the LiDAR point cloud data of the laser radar, so that the measurement precision and efficiency can be greatly improved, and the degree of automation is improved;
2. the random sampling consistency is utilized to fit the segmented road surface boundary, so that the accidental calculation of the circle center result by the single dynamic circle radius is eliminated, and the self-adaption capability of road data fitting is improved.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings that are required to be used in the embodiments or the description of the prior art will be briefly described below, and it is obvious that the drawings in the following description are only embodiments of the present invention, and that other drawings can be obtained according to the provided drawings without inventive effort for a person skilled in the art.
FIG. 1 is a flow chart of a high-speed loop line fitting method based on point cloud random sampling;
FIG. 2 is a schematic drawing of the extraction of Alpha shapes contour points according to the present invention;
FIG. 3 is a schematic view of the Alpha shapes contour point judgment of the present invention;
FIG. 4 shows the effect of Alpha values of the Alpha shapes method of the present invention.
Detailed Description
The following description of the embodiments of the present invention will be made clearly and completely with reference to the accompanying drawings, in which it is apparent that the embodiments described are only some embodiments of the present invention, but not all embodiments. All other embodiments, which can be made by those skilled in the art based on the embodiments of the invention without making any inventive effort, are intended to be within the scope of the invention.
As shown in fig. 1, the embodiment of the invention discloses a high-speed loop linear fitting method based on point cloud random sampling, which comprises the following specific steps:
s1, preprocessing an original point cloud to obtain a Lu Miandian cloud; extracting linear point clouds of a high-speed loop through a laser radar; the road point cloud comprises a straight line segment road point cloud and a curve segment Lu Miandian cloud;
s2, carrying out boundary detection on the road point cloud by utilizing the Alpha shapes section to obtain a boundary point cloud; the boundary point cloud comprises a straight line segment boundary point cloud and a curve segment boundary point cloud;
s3, carrying out circular curve boundary fitting and circle center determination on the boundary point cloud based on random sampling consistency, and carrying out straight line segment boundary fitting on the boundary point cloud based on random sampling consistency to determine the size of a field boundary;
s4, reversely pushing the length of the field relief curve and the length of the field straight line according to the field boundary size and the position of the circular curve.
Further, the specific step of S1 includes:
s1.1, filtering an original point cloud and generating a denoising point cloud; the method for filtering the original point cloud adopts statistical filtering; when point cloud data is acquired through LiDAR, noise points inevitably appear in the point cloud data due to the influences of equipment precision, operator experience, environmental factors and the like, as well as electromagnetic wave diffraction characteristics, surface property changes of an object to be measured and the data stitching registration operation process. Common point cloud filtering methods include direct filtering, voxel filtering, statistical filtering, gaussian filtering and the like, and in various filtering methods, the high-speed loop point cloud filtering is suitable for adopting the statistical filtering in consideration of the characteristic that the high-speed loop is taken as the whole point cloud.
S1.2, dividing and denoising point cloud based on combination of pavement normal and elevation mean square error, and generating a pavement point cloud.
Further, the specific steps of S1.1 include:
s1.1.1 denoising the original point cloud once, selecting any point in the original point cloud as a first point in each denoising, and adding the first point and N nearest to the first point 1 -1 point constitutes a set of points; wherein, after denoising all points in the original point cloud, denoising point cloud is generated, N 1 Is an integer greater than 2; the point set formula:
wherein: x represents a point set, X 1 Represents a first point nearest to the first point, x 2 Representing a second point closest to the first point,indicating the third point nearest to the first point, < ->Representing a first point; wherein the point set accords with normal distribution;
s1.1.2 for each ofPoint, calculate +.>Point-to-point integration divide +.>The specific formula of the distance of any point except the point is as follows:
wherein:representation->Point-to-point integration divide +.>Distance of any point other than point, +.>Representing the +.>Any point other than;
s1.1.3 atPoint to except->Selecting a distance group with the smallest kmean number from the distances of any point as a distance set, wherein the distance set formula is as follows:
wherein:representing distance set,/- >Representation->Except->Distance of nearest first point outside the point, +.>Representation->Except->Distance of nearest second point outside the point, +.>Representation->Except->The distance of the closest kmean point outside the point; wherein kmean is an integer greater than 1;
s1.1.4 the mean of the distance setThe average value calculation formula of the characteristic value of the distance set is as follows:
wherein:representing a distance set mean;
s1.1.5 judging outliers based on a distance mean formula and a variance mean formula, removing outliers from original cloud points to reserve effective points, and generating denoising point clouds; the distance average formula is:
wherein: mean represents the mean of the distance set means of all points in the point set;
the distance variance formula is:
wherein: stddev represents the variance of the distance set mean of all points in the point set and the distance set mean of all points in the point set;
if it isWhen then x j Is an outlier; if->When then x j Is the effective point, where the C value is the outlier of the point.
Furthermore, based on the analysis of the characteristics of the high-speed loop, as the road surface point cloud of the high-speed loop is dense and the noise points on two sides are more obvious, the N value can be properly increased to remove more non-road surface points, and the filtering test shows that when the number of points N is more than 100 and the outlier C is more than 1, most of the noise points can be removed. The point cloud after the outliers are removed is smoother, so that the method is suitable for the next step of point cloud clustering segmentation processing and is easier to converge.
Besides the noise point, other road auxiliary facilities such as marks, guardrails, side ditches, greening and pavement workers, equipment and the like can also influence the calculation and analysis of the point cloud pavement. The method is limited by the limitation of a filtering algorithm, the characteristics of the point cloud are analyzed by further combining the characteristics of a high-speed loop in the processing of non-road point cloud, and for the high-speed loop, obvious characteristic differences exist between a straight line segment and a curve segment, various facilities are distributed on the two sides of the straight line segment with small transverse slopes without ultrahigh, and boundary characteristics are not obvious; the curve part has large transverse slope, and the two sides are only provided with guardrails, so that the boundary characteristics are obvious.
Further, the specific steps of S1.2 include:
s1.2.1 optionally selecting one point in the point cloud as a query pointBy inquiring about the point->Nearest k 2 Obtaining a plane point set to be fitted as +.> The plane equation to be constructed is:
wherein:representing the normal vector of the fitting plane, D representing the distance from the origin to the plane;
any pointDistance to plane->Is that
The best fit plane needs to satisfy the conditional function
Taking a function
Extreme point with f to D offset of 0
Derived from
Wherein:derived->
Wherein:and f is->Extreme point with deviation obtaining guide of 0
Converting the above formula to obtain
In the middle of
Solving for k 2 A covariance matrix A created by each neighboring point is fitted with a feature vector and a feature value of the covariance matrix by a least square method to obtain a point cloud normal; each point in the pair point setIts corresponding covariance matrix can be expressed as:
wherein: a represents covariance matrix, k 2 Representation pointsAdjacent points number, P i Represents a query point->Representing the three-dimensional centroid of a neighboring point, +.>J-th representing covariance matrix 2 A feature vector lambda j2 J-th representing covariance matrix 2 A characteristic value;
s1.2.2 when it isFeature vector corresponding to the minimum feature value is obtained>The normal vector of the calculation point is: />
Wherein:represents the point normal +_>Representing the normal x-direction component, +.>Representing the normal y-direction component, +.>Representing the normal z-direction component, T being the transposed symbol;
s1.2.3 calculates the gradient value of each point, and the gradient value calculation formula is as follows:
wherein: tan alpha 1 Representing a grade value;
s1.2.4 comparing the gradient value with a gradient threshold value through gradient, and determining a dividing point and a pavement point of a high-speed loop line;wherein, if tan alpha 1 ≤tanα 2 Then it is the dividing point, if tan alpha 1 >tanα 2 Then is the road surface point, tan alpha 2 Representing a grade threshold;
according to the characteristics of the high-speed loop, the gradient threshold value of the road surface is the outer edge point of the circular curve section, and the gradient threshold value is the maximum superhigh of the high-speed loop, so that the road surface can be obtained through the history data of the high-speed loop; in the absence of high speed loop history, the grade threshold may be determined by local statistical analysis.
The S1.2.5 gradient classification method can realize the segmentation of most road surface point clouds, but the single gradient separation method has inherent defects due to overlarge composite gradient of the high-speed loop road surface, and the characteristics of the high-speed loop road surface point clouds also need to be further analyzed. The method is characterized by integrating the characteristics of the high-speed loop road surface, and can be used as a supplement method for the segmentation of the high-speed loop point cloud road surface by adopting the standard deviation based on the local point cloud elevation. Setting a local cylinder in the denoising point cloud, wherein the radius of the cylinder is preferably 0.05-0.1 meter, and the heights of each point in the range of the local cylinder are respectively as follows:
wherein: z is Z G Represent high Cheng Ge set, z G1 Representing the elevation, z, of a first point within the partial cylinder G2 Representing the elevation of the second point within the partial cylinder,representing the Nth in a partial cylinder 2 Elevation of point N 2 Representing the number of points within the partial cylinder; wherein N is 2 Is an integer greater than 1;
s1.2.6 calculating the elevation mean value and the elevation standard deviation of each point in the range of the local cylinder, and judging the local cylinder as a partition point when the elevation standard deviation exceeds the elevation standard deviation limit value; the calculation formula of the elevation mean value is as follows:
wherein:represent elevation mean, i 3 Representing the ith in a partial cylinder 3 A dot;
the calculation formula of the elevation standard deviation is as follows:
wherein:representing the standard deviation of elevation; wherein if- >Then each point in the range of the local cylinder is a foreign point, < >>Representing a preset elevation standard deviation limit value;
similar to the gradient discrimination method, obvious inflection points can be found on the sequential point cloud elevation standard deviation change chart to serve as standard deviation discrimination basis. The road surface point cloud after the point cloud pretreatment can realize ideal segmentation effect, and the road surface point cloud segmentation effect can reach more than 99%.
S1.2.7 the cloud of denoising points is divided by dividing points with the road surface points as road surface areas, and Lu Miandian cloud is obtained.
Further, the specific step of S2 includes:
the preprocessed high-loop surface point cloud needs to be further subjected to boundary extraction to determine the outline dimension of the high-speed loop. The common point cloud boundary extraction algorithm includes a warp and weft scanning method, a grid dividing method, a normal line estimation method, an Alpha shapes algorithm and the like. Because the high-ring point cloud plane extraction is mainly based on two-dimensional contour points, the Alpha shapes algorithm in the method is a simpler and effective algorithm for rapidly extracting boundary points, and the method overcomes the defect of influence of the shape of the point cloud boundary points and can rapidly and accurately extract the boundary points. The principle is that for plane point cloud with arbitrary shape, if a circle with radius alpha rolls around it, the rolling track forms a point which is the contour point. When extracting the outline point of the point cloud, the point needs to be projected on a two-dimensional plane, typically a projection plane, and at this time, the z coordinates of all points of the road surface point cloud are set to be the same constant, namely, the z=0 plane.
S2.1 selecting any point in the plane projection point cloud as a point P to be judged e (x e ，y e ) According to the preset rolling circle radius alpha e Searching distance P from plane projection point cloud e The point distance is less than 2 alpha e All points in the set are marked as a point set Q; the generation method of the plane projection point cloud comprises the steps of enabling all points z of the road surface point cloud to be e The coordinates are all set to the same constant;
s2.2 As shown in FIG. 2, any point in the point set Q is selected as a point P e1 (x e1 ，y e1 ) According to point P e Point P e1 And rolling circle radius alpha e Calculating the center coordinates O of the first circle 1 (x o1 ，y o1 ) And the center coordinates O of the second circle 2 (x o2 ，y o2 ) The method comprises the steps of carrying out a first treatment on the surface of the As shown in fig. 3, wherein the first circle intersects the second circle and has two points of intersection, the rolling circle radius α e The radius alpha of the rolling circle is the distance between the center of the first circle and two intersection points e Is the distance between the center of the second circle and two intersecting points, P e 、P e1 Two intersecting points of the first circle and the second circle respectively;
s2.3 dividing P for Point set Q e1 Each point outside the circle is calculated to the center O 1 And to O 2 When all points are at the distance of O 1 And to O 2 Are all greater than alpha e When the point P is marked e Terminating the point P as the contour point e Is calculated;
s2.4 dividing P in the Point set Q e1 Any point outside to the centre of a circle O 1 And to O 2 When the distance of (2) is less than or equal to alpha, then traversing all the points in the point set Q to rotate as P e1 Point, recalculate the center coordinates and repeat step S2.3 to point P e Judging; if there is one point is point P e If the mark is a contour point, the point P can be judged e As contour points, otherwise, judging point P e Is a non-contour point.
According to the above calculation steps, a reasonable rolling circle radius α needs to be determined to complete extraction of the high-ring point cloud plane contour, and the point cloud contour calculation is shown in fig. 4, because the dynamic circle radius α has a certain influence on the point cloud contour, the value rationality of the dynamic circle radius α has a larger influence on the linear fitting. When the alpha value is smaller than the distance at the point cloud boundary, the point cloud boundary cannot be extracted; when the alpha value is larger than the point cloud boundary interval and smaller than the point cloud average interval, extracting a large number of holes contained in the point cloud, and influencing the fitting result; when the alpha value is larger than the average distance of the point clouds, the boundary point clouds are gradually sparse along with the increase of the alpha value, and the influence on subsequent fitting is larger.
Further, the specific step of S3 includes:
s3.1, fitting a point cloud boundary contour by adopting RANSAC;
s3.2, an output model obtained by fitting a circular curve of the point cloud boundary profile through S3.1 is as follows:
(x h -x ho ) 2 +(y h -y ho ) 2 ＝R h 2
wherein: (x) h ,y h ) Representing the coordinates of points on a circle curve, (x) ho ,y ho ) Represents the center coordinates of a circular curve, R h Representing the radius of the circular curve.
S3.3, an output model obtained by fitting the S3.1 to the straight line segment of the point cloud boundary contour is as follows:
k h (x h -x h1 )＝(y h -y h1 )
wherein: (x) h ,y h ) Representing point coordinates on a straight line, (x) h1 ,y h1 ) Representing the coordinates of a known point on a straight line, y h1 Representation, k h The slope of the straight line is shown.
For the whole model, the result output by RANSAC is influenced by the boundary extraction result and the fitting threshold value, and a certain accidental exists in the result of a single parameter. In order to reduce the influence of accidental results, the influence of parameter variation on the fitting result needs to be reduced, and the standard deviation of calculation results of the dynamic circle radius of different boundaries under the same fitting threshold condition is calculated by analyzing the combined parameters of the dynamic circle radius alpha and the fitting threshold value, so that the fitting result has lower discreteness and better fitting result when the fitting threshold value is between 0.004 and 0.006.
In S3.2, the circle radius α is moved with different boundary contours for the circle curve b1 ，α b2 ，…，α bn Calculation will obtain different circle centers (x bo1 ，y bo1 )，(x bo2 ，y bo2 )，…，(x bon ，y bon ) To eliminate the accidental calculation of the center result by a single radius of a movable circle, the center of gravity of a set of different calculated centers is preferably selectedFor the discrete center point set, the calculation method of the barycentric coordinates is the arithmetic average value of the coordinates of each point, namely, the calculation formula of the barycentric coordinates is as follows:
Wherein:represents the abscissa of gravity, x boi Represents the abscissa of discrete circle center,/->Represents the abscissa of gravity, y boi Representing the ordinate of the discrete center of circle.
In S3.2, for straight lines, because the straight line of the high-speed loop is longer, when the RANSAC is adopted for linear fitting, when the fitting threshold value is smaller, the fitting straight line precision is high, but the number of boundary sub-set points is small, the maximum distance of fitting points is smaller, and the integrity of the fitted straight line is poor; on the contrary, when the fitting threshold is larger, the number of boundary subset points is larger, the maximum distance of the fitting points is larger, but the RANSAC output model accuracy is lower. For a complete high-speed loop, 4 boundary lines exist in two sections of straight lines, the maximum lengths of fitting point clouds are different due to the influence of point cloud acquisition precision and entering and exiting section edge line widening, and when the maximum length of each straight line fitting point reaches stability, the minimum threshold value can be used as a straight line fitting result.
Further, the specific step S3.1 includes:
after the high-speed loop point cloud boundary contour extraction is completed, corresponding fitting is required to be carried out twice on the straight line and the circular curve, and the linear index can be obtained. In the step S3, a random sampling consistency RANSAC method is adopted to respectively develop and fit the circular curve and the straight line segment, and the method can estimate the parameters of the mathematical model in an iterative mode from a group of observation data sets containing 'local outer points'. This is an uncertain algorithm that can derive a reasonable result from a certain probability, and the number of iterations must be increased in order to increase the probability.
The basic assumption of RANSAC is: (1) the data consists of "intra-office points", such as: the distribution of the data may be interpreted with some model parameters; (2) the "outlier" cannot fit the data of the model; the data other than (3) belongs to noise. The reasons for the generation of the off-site points are: extreme values of noise; an erroneous measurement method; false assumptions about data.
The inputs to the RANSAC algorithm are a set of boundary data, a parameterized model that can interpret or adapt to the observed data, and some parameters. RANSAC achieves this goal by iteratively selecting a random subset of the data. The selected subset is assumed to be an intra-office point and verified by:
s3.1.1 randomly selecting a subset S from the boundary point set S, taking the randomly selected subset S as assumed local points, and fitting the initial model by the local points S;
estimating the modulus from the subset sThe calculation of the parameters of the circular curve assumes that the coordinates of the local points are (x β1 ,y β1 ),(x β2 ,y β2 ),(x β3 ,y β3 ) With S3.2 circular curve model
(x β1 -x ho ) 2 +(y β1 -y ho ) 2 ＝R h 2
(x β2 -x ho ) 2 +(y β2 -y ho ) 2 ＝R h 2
(x β3 -x ho ) 2 +(y β3 -y ho ) 2 ＝R h 2
Solving the equation set to obtain
R h 2 ＝(x β1 -x ho ) 2 +9y β1 -y ho ) 2
Straight line parameter calculation assuming that the local point coordinates are (x β1 ,y β1 ),(x β2 ,y β2 ) With S3.3 straight line model
k h (x β1 -x h1 )＝(y β1 -y h1 )
k h (x β2 -x h1 )＝(y β2 -y h1 )
Solving the equation set to obtain
k h ＝(y β2 -y β1 )/(x β2 -x β1 )
x h1 ＝x β1
y h1 ＝y β1
S3.1.2 traversing all but subset S of the set of boundary points SFrom this, the radius R of the circular curve is calculated h Or slope k of straight line h If the error is within the error limit e, marking as an inner point, otherwise marking as an outer point;
s3.1.3 forming a consistent set by all the interior points obtained in S3.1.2, re-estimating model parameters by using all the interior points in the consistent set if the number of the consistent set points meets a given threshold T, and outputting a result; if the number of the inner points in the consistent set is less than the threshold value T, a new subset s is reselected, and the steps S3.1.1-S3.1.2 are repeated;
s3.1.4 through multiple iterations of S3.1.1-S3.1.3, selecting a consistent set with the largest number of interior points, re-estimating model parameters by using all interior points in the consistent set, and outputting a result; wherein the number of iterations is an integer number of iterations greater than 1.
Further, the specific step of S4 includes:
s4.1, calculating the overall field width of the high-speed loop and the overall field length of the high-speed loop through a drift distance formula, and then respectively obtaining a relaxation curve longitudinal drift distance and a relaxation curve transverse drift distance; the formula of the relation between the width and the length of the field and the offset distance of the relaxation curve is as follows:
wherein: w represents the whole field width of the high-speed loop, and k 1 The slope of the straight line is shown,the abscissa representing the midpoint of the center of the two sides, +.>An ordinate representing the center point of two sides, y 1 Representing the ordinate, x, of a known point on the first fitted line 1 An abscissa, y, representing a known point on the first fitted line 2 Representing the ordinate, x, of a known point on the second fitting line 2 Represents the abscissa of the known point on the second fitting straight line, M represents the lane width adjustment value between the straight line segment center line and the fitting straight line, beta represents the moderation curve angle, R represents the radius of the circular curve, L y Represents the longitudinal offset distance of a field relaxation curve, H represents the whole field length of a high-speed loop, and X o1 X represents the circular center abscissa of the first end point of the high-speed loop o2 Represents the abscissa of the center of a circle of a second endpoint of the high-speed loop, Y o1 Represents the ordinate of the circular center of the first end point of the high-speed loop, Y o2 Represents the ordinate of the center of a circle of a second endpoint of the high-speed loop, L represents the length of a straight line of a field, and L x Representing the ground relief curve lateral offset.
Further, the result of the model output in S3 is based on the calculation and fitting of the point cloud, but there is still a difference between the real point cloud and the design index. The design index fitted by the circular curve part cannot be completely matched with the design index fitted by the straight line, and the slope of the connecting line between the two straight lines and the circle center is unequal as shown in the following figure, namely V 1 /U 1 ≠V 2 /U 2 ≠(y o2 -y o1 )/(x o2 -x o1 ). For this reason, the high loop shape requires further fitting analysis and correction.
The high-speed loop fitting line shape needs to be defined, and the indexes comprise center coordinates (Xo 1, yo 1), (Xo 2, yo 2) at two ends, a field straight line length L, a field relaxation curve length Ls and a circular curve radius R; the circle center coordinates are obtained through RANSAC fitting, the radius of a circle curve can be calculated through the boundary fitting radius and the offset value of a design line, and the circle center coordinates are comprehensively determined after being compared with the original design indexes, and the length of a straight line and the length of a moderation curve are unknown parameters; according to the characteristics of the moderation curve, the known R, ls can calculate moderation curve offset distances x and y and moderation curve angle beta, so as to further derive the overall field width W and length H of the high-speed loop, when the calculation is performed, the W adopts the fitting straight line spacing at two sides, the lane width adjustment is performed according to the center line position of the straight line segment, and the H adopts the center coordinate distance of two ends to determine.
S4.2, calculating a moderation curve angle by a curvature quadratic integration method; the relaxation curve angle β is calculated using the following formula:
wherein: ls (Ls) 0 Representing the assumed length of the relaxation curve, k (l) representing the function of the curvature of any point on the relaxation curve as a function of the length of the relaxation curve;
the moderation curve lateral offset and moderation curve longitudinal offset are calculated by adopting the following steps:
s4.3 after linear fitting and correction, the length Ls of the assumed relaxation curve is 0 And calculating errors of the corresponding field width and the actual field width W, obtaining the final relaxation curve length Ls through an approximation method, and calculating the field straight line length L through the field length H to obtain the index, thereby completing the loop plane line shape high-precision fitting.
Further, the relaxation curve curvature function k (l) in S4.2 is any one of a Bloss curve curvature function, a McConnell curve curvature function, a Grabowski curve curvature function, an Auberlen curve curvature function, and a Klein curve curvature function;
the Bloss curve curvature function is:
wherein: k (l) Bloss Representing a Bloss curve curvature function;
the McConnell curve curvature function is:
wherein: k (l) McConnell Representing a McConnell curve curvature function, wherein ψ represents a maximum superhigh angle, v represents a design speed, and g represents a gravitational acceleration;
the Grabowski curve curvature function is:
wherein: k (l) Grabowski Representing a Grabowski curve curvature function;
the Auberlen curve curvature function is:
wherein: k (l) Auberlen Representing the Auberlen curve curvature function;
the Klein curve curvature function is:
wherein: k (l) Klein Representing the Klein curve curvature function.
In the present specification, each embodiment is described in a progressive manner, and each embodiment is mainly described in a different point from other embodiments, and identical and similar parts between the embodiments are all enough to refer to each other. For the device disclosed in the embodiment, since it corresponds to the method disclosed in the embodiment, the description is relatively simple, and the relevant points refer to the description of the method section.
The previous description of the disclosed embodiments is provided to enable any person skilled in the art to make or use the present invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the generic principles defined herein may be applied to other embodiments without departing from the spirit or scope of the invention. Thus, the present invention is not intended to be limited to the embodiments shown herein but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.
Claims (9)
1. The high-speed loop linear fitting method based on the point cloud random sampling is characterized by comprising the following specific steps of:
s1, preprocessing an original point cloud to obtain a Lu Miandian cloud; extracting linear point clouds of a high-speed loop through a laser radar; the road point cloud comprises a straight line segment road point cloud and a curve segment Lu Miandian cloud;
s2, carrying out boundary detection on the road point cloud by utilizing the Alpha shapes section to obtain a boundary point cloud; the boundary point cloud comprises a straight line segment boundary point cloud and a curve segment boundary point cloud;
s3, carrying out circular curve boundary fitting and circle center determination on the boundary point cloud based on random sampling consistency, and carrying out straight line segment boundary fitting on the boundary point cloud based on random sampling consistency to determine the size of a field boundary;
S4, reversely pushing the length of the field relief curve and the length of the field straight line according to the field boundary size and the position of the circular curve.
2. The high-speed loop line fitting method based on point cloud random sampling as claimed in claim 1, wherein the specific steps of S1 include:
s1.1, filtering an original point cloud and generating a denoising point cloud; the method for filtering the original point cloud adopts statistical filtering;
s1.2, dividing and denoising point cloud based on combination of pavement normal and elevation mean square error, and generating a pavement point cloud.
3. The high-speed loop line fitting method based on point cloud random sampling according to claim 2, wherein the specific step of S1.1 comprises the following steps:
s1.1.1 denoising the original point cloud once, selecting any point in the original point cloud as a first point in each denoising, and adding the first point and N nearest to the first point 1 -1 point constitutes a set of points; wherein, after denoising all points in the original point cloud, denoising point cloud is generated, N 1 Is an integer greater than 2; the point set formula:
wherein: x represents a point set, X 1 Represents a first point nearest to the first point, x 2 Representing a second point closest to the first point,indicating the third point nearest to the first point, < ->Representing a first point; wherein the point set accords with normal distribution;
S1.1.2 for each ofPoint, calculate +.>Point-to-point integration divide +.>The specific formula of the distance of any point except the point is as follows:
wherein:representation->Point-to-point integration divide +.>Distance of any point other than point, +.>Representing the +.>Any point other than;
s1.1.3 atPoint to except->Selecting the minimum kmean distances from the distances of any point to generate a distance set, wherein the distance set formula is as follows:
wherein:representing distance set,/->Representation->Except->The distance of the nearest first point outside the point,representation->Except->Distance of nearest second point outside the point, +.>Representation->Except->The distance of the closest kmean point outside the point; wherein kmean is an integer greater than 1;
s1.1.4 the mean of the distance setThe average value calculation formula of the characteristic value of the distance set is as follows:
wherein:representing a distance set mean;
s1.1.5 judging outliers based on a distance mean formula and a variance mean formula, removing outliers from original cloud points to reserve effective points, and generating denoising point clouds; the distance average formula is:
wherein: mean represents the mean of the distance set means of all points in the point set;
the distance variance formula is:
Wherein: stddev represents the variance of the distance set mean of all points in the point set and the distance set mean of all points in the point set;
if it isWhen then x j Is an outlier; if->When then x j Is the effective point, where the C value is the outlier of the point.
4. The high-speed loop line fitting method based on point cloud random sampling as claimed in claim 3, wherein the specific steps of S1.2 include:
s1.2.1 optionally selecting one point in the point cloud as a query pointBy inquiring about the point->Nearest k 2 Solving k by adjacent points 2 A covariance matrix A created by each neighboring point is fitted with a feature vector and a feature value of the covariance matrix by a least square method to obtain a point cloud normal; each point in the pair point set +.>Its corresponding covariance matrix can be expressed as:
wherein: a represents covariance matrix, k 2 Representation pointsAdjacent points number, P i Represents a query point->Representing the three-dimensional centroid of a neighboring point, +.>J-th representing covariance matrix 2 Individual feature vectors->J-th representing covariance matrix 2 A characteristic value;
s1.2.2 when it isFeature vector corresponding to the minimum feature value is obtained>The normal vector of the calculation point is:
wherein:represents the point normal +_>Representing the normal x-direction component, +.>Representing the normal y-direction component, +. >Representing the normal z-direction component, T being the transposed symbol;
s1.2.3 calculates the gradient value of each point, and the gradient value calculation formula is as follows:
wherein: tan alpha 1 Representing a grade value;
s1.2.4 comparing the gradient value with a gradient threshold value through gradient, and determining a dividing point and a pavement point of a high-speed loop line; wherein, if tan alpha 1 ≤tanα 2 Then it is the dividing point, if tan alpha 1 >tanα 2 Then is the road surface point, tan alpha 2 Representing a grade threshold;
s1.2.5 a local cylinder is set in the denoising point cloud, and the heights of each point in the range of the local cylinder are respectively as follows:
wherein: z is Z G Represent high Cheng Ge set, z G1 Representing the elevation, z, of a first point within the partial cylinder G2 Representing the elevation of the second point within the partial cylinder,representing the Nth in a partial cylinder 2 Elevation of point N 2 Representing the number of points within the partial cylinder; wherein N is 2 Is an integer greater than 1;
s1.2.6 calculating the elevation mean value and the elevation standard deviation of each point in the range of the local cylinder, and judging the local cylinder as a division point when the elevation standard deviation exceeds the elevation standard deviation limit value; the calculation formula of the elevation mean value is as follows:
wherein:represent elevation mean, i 3 Representing the ith in a partial cylinder 3 A dot;
the calculation formula of the elevation standard deviation is as follows:
wherein:representing the standard deviation of elevation; wherein if- >The points within the local cylinder are outlier points,representing a preset elevation standard deviation limit value;
s1.2.7 the cloud of denoising points is divided by dividing points with the road surface points as road surface areas, and Lu Miandian cloud is obtained.
5. The high-speed loop line fitting method based on point cloud random sampling as claimed in claim 4, wherein the specific step of S2 comprises:
s2.1 selecting any point in the plane projection point cloud as a point P to be judged e (x e ，y e ) According to the preset rolling circle radius alpha e Searching distance P from plane projection point cloud e The point distance is less than 2 alpha e All points in the set are marked as a point set Q; the generation method of the plane projection point cloud comprises the steps of enabling all points z of the road surface point cloud to be e The coordinates are all set to the same constant;
s2.2 selecting any point in the point set Q as a point P e1 (x e1 ，y e1 ) According to point P e Point P e1 And rolling circle radius alpha e Calculating the center coordinates O of the first circle 1 (x o1 ，y o1 ) And the center coordinates O of the second circle 2 (x o2 ，y o2 ) The method comprises the steps of carrying out a first treatment on the surface of the Wherein the first circle intersects with the second circle and has two points of intersection, the radius alpha of the rolling circle e The radius alpha of the rolling circle is the distance between the center of the first circle and two intersection points e Is the distance between the center of the second circle and two intersecting points, P e 、P e1 Two intersecting points of the first circle and the second circle respectively;
S2.3 dividing P for Point set Q e1 Each point outside the circle is calculated to the center O 1 And to O 2 When all points are at the distance of O 1 And to O 2 Are all greater than alpha e When the point P is marked e Terminating the point P as the contour point e Is calculated;
s2.4 in the Point set QP removal e1 Any point outside to the centre of a circle O 1 And to O 2 When the distance of (2) is less than or equal to alpha, then traversing all the points in the point set Q to rotate as P e1 Point, recalculate the center coordinates and repeat step S2.3 to point P e Judging; if there is one point is point P e If the mark is a contour point, the point P can be judged e Is the contour point, and the contour point set obtained is the boundary point set S, otherwise, the point P is judged e Is a non-contour point.
6. The high-speed loop line fitting method based on point cloud random sampling as claimed in claim 5, wherein the specific step of S3 comprises:
s3.1, fitting a point cloud boundary contour by adopting RANSAC;
s3.2, an output model obtained by fitting a circular curve of the point cloud boundary profile through S3.1 is as follows:
(x h -x ho ) 2 +(y h -y ho ) 2 ＝R h 2
wherein: (x) h ,y h ) Representing the coordinates of points on a circle curve, (x) ho ,y ho ) Represents the center coordinates of a circular curve, R h Representing the radius of a circular curve;
s3.3, an output model obtained by fitting the S3.1 to the straight line segment of the point cloud boundary contour is as follows:
k h (x h -x h1 )＝(y h -y h1 )
wherein: (x) h ,y h ) Representing point coordinates on a straight line, (x) h1 ,y h1 ) Representing the coordinates of a known point on a straight line, y h1 Representation, k h The slope of the straight line is shown.
7. The high-speed loop line fitting method based on point cloud random sampling as claimed in claim 6, wherein the specific step of S3.1 comprises the following steps:
s3.1.1 randomly selecting a subset S from the boundary point set S, taking the randomly selected subset S as assumed local points, and fitting the initial model by the local points S;
estimating model parameters according to the subset s, calculating the coordinates of the assumed local points as (x) by using the circular curve parameters β1 ,y β1 ),(x β2 ,y β2 ),(x β3 ,y β3 ) With S3.2 circular curve model
(x β1 -x ho ) 2 +(y β1 -y ho ) 2 ＝R h 2
(x β2 -x ho ) 2 +(y β2 -y ho ) 2 ＝R h 2
(x β3 -x ho ) 2 +(y β3 -y ho ) 2 ＝R h 2
Solving the equation set to obtain
R h 2 ＝(x β1 -x ho ) 2 +(y β1 -y ho ) 2
Straight line parameter calculation assuming that the local point coordinates are (x β1 ,y β1 ),(x β2 ,y β2 ) With S3.3 straight line model
k h (x β1 -x h1 )＝(y β1 -y h1 )
k h (x β2 -x h1 )＝(y β2 -y h1 )
Solving the equation set to obtain
k h ＝(y β2 -y β1 )/(x β2 -x β1 )
x h1 ＝x β1
y h1 ＝y β1
S31.2 traversing all data except subset S in boundary point set S, calculating radius R of circular curve h Or slope k of straight line h If the error is within the error limit e, marking as an inner point, otherwise marking as an outer point;
s3.1.3 forming a consistent set by all the interior points obtained in S3.1.2, re-estimating model parameters by using all the interior points in the consistent set if the number of the consistent set points meets a given threshold T, and outputting a result; if the number of the inner points in the consistent set is less than the threshold value T, a new subset s is reselected, and the steps S3.1.1-S3.1.2 are repeated;
S3.1.4 through multiple iterations of S3.1.1-S3.1.3, selecting a consistent set with the largest number of interior points, re-estimating model parameters by using all interior points in the consistent set, and outputting a result; wherein the number of iterations is an integer number of iterations greater than 1.
8. The high-speed loop line fitting method based on point cloud random sampling as claimed in claim 6, wherein the specific step of S4 comprises:
s4.1, calculating the overall field width of the high-speed loop and the overall field length of the high-speed loop through a drift distance formula, and then respectively obtaining a relaxation curve longitudinal drift distance and a relaxation curve transverse drift distance; the formula of the relation between the width and the length of the field and the offset distance of the relaxation curve is as follows:
wherein: w represents the whole field width of the high-speed loop, and k 1 The slope of the straight line is shown,the abscissa representing the center point of the circle at the two sides,an ordinate representing the center point of two sides, y 1 Representing the ordinate, x, of a known point on the first fitted line 1 An abscissa, y, representing a known point on the first fitted line 2 Representing the ordinate, x, of a known point on the second fitting line 2 Represents the abscissa of the known point on the second fitting straight line, M represents the lane width adjustment value between the straight line segment center line and the fitting straight line, beta represents the moderation curve angle, R represents the radius of the circular curve, L y Represents the longitudinal offset distance of a field relaxation curve, H represents the whole field length of a high-speed loop, and X o1 X represents the circular center abscissa of the first end point of the high-speed loop o2 Represents the abscissa of the center of a circle of a second endpoint of the high-speed loop, Y o1 Represents the ordinate of the circular center of the first end point of the high-speed loop, Y o2 Represents the ordinate of the center of a circle of a second endpoint of the high-speed loop, L represents the length of a straight line of a field, and L x Representing the field relaxation curve lateral offset;
s4.2, calculating a moderation curve angle by a curvature quadratic integration method; the relaxation curve angle β is calculated using the following formula:
wherein: ls (Ls) 0 Representing the assumed length of the relaxation curve, k (l) representing the function of the curvature of any point on the relaxation curve as a function of the length of the relaxation curve;
the moderation curve lateral offset and moderation curve longitudinal offset are calculated by adopting the following steps:
s4.3 after linear fitting and correction, the length Ls of the assumed relaxation curve is 0 And calculating errors of the corresponding field width and the actual field width W, obtaining the final relaxation curve length Ls through an approximation method, and calculating the field straight line length L through the field length H.
9. The high-speed loop linear fitting method based on point cloud random sampling according to claim 8, wherein the moderation curve curvature function k (l) in S4.2 adopts any one of a Bloss curve curvature function, a McConnell curve curvature function, a Grabowski curve curvature function, an Auberlen curve curvature function and a Klein curve curvature function;
The Bloss curve curvature function is:
wherein: k (l) Bloss Representing a Bloss curve curvature function;
the McConnell curve curvature function is:
wherein: k (l) McConnell Representing a McConnell curve curvature function, wherein ψ represents a maximum superhigh angle, v represents a design speed, and g represents a gravitational acceleration;
the Grabowski curve curvature function is:
wherein: k (l) Grabowski Representing a Grabowski curve curvature function;
the Auberlen curve curvature function is:
wherein: k (l) Auberlen Representing the Auberlen curve curvature function;
the Klein curve curvature function is:
wherein: k (l) Klein Representing the Klein curve curvature function.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310611822.2A CN116645348A (en) | 2023-05-29 | 2023-05-29 | High-speed loop linear fitting method based on point cloud random sampling |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310611822.2A CN116645348A (en) | 2023-05-29 | 2023-05-29 | High-speed loop linear fitting method based on point cloud random sampling |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116645348A true CN116645348A (en) | 2023-08-25 |
Family
ID=87624116
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310611822.2A Pending CN116645348A (en) | 2023-05-29 | 2023-05-29 | High-speed loop linear fitting method based on point cloud random sampling |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116645348A (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117553686A (en) * | 2024-01-12 | 2024-02-13 | 成都航空职业技术学院 | Laser radar point cloud-based carriage bulk cargo overrun detection method |
-
2023
- 2023-05-29 CN CN202310611822.2A patent/CN116645348A/en active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117553686A (en) * | 2024-01-12 | 2024-02-13 | 成都航空职业技术学院 | Laser radar point cloud-based carriage bulk cargo overrun detection method |
| CN117553686B (en) * | 2024-01-12 | 2024-05-07 | 成都航空职业技术学院 | Laser radar point cloud-based carriage bulk cargo overrun detection method |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN115424232B (en) | Method for identifying and evaluating pavement pit, electronic equipment and storage medium | |
| CN105783779B (en) | The real-time form identification of rail profile and distortion calibration method based on three layers of matching | |
| CN109887015A (en) | A kind of point cloud autoegistration method based on local surface feature histogram | |
| CN107301648B (en) | Redundant point cloud removing method based on overlapping area boundary angle | |
| Mizoguchi et al. | Quantitative scaling evaluation of concrete structures based on terrestrial laser scanning | |
| CN107610223A (en) | Power tower three-dimensional rebuilding method based on LiDAR point cloud | |
| CN109658431B (en) | Rock mass point cloud plane extraction method based on region growth | |
| CN110807781B (en) | Point cloud simplifying method for retaining details and boundary characteristics | |
| CN116645348A (en) | High-speed loop linear fitting method based on point cloud random sampling | |
| CN103295232A (en) | SAR (specific absorption rate) image registration method based on straight lines and area | |
| CN106529593B (en) | Pavement disease detection method and system | |
| CN102750449B (en) | Point cloud linear feature extraction method based on substep three-dimensional space and feature dimension mapping | |
| CN105387801A (en) | Subway tunnel segment dislocation quantity detection method | |
| CN105069395A (en) | Road marking automatic identification method based on terrestrial three-dimensional laser scanning technology | |
| CN111354083A (en) | Progressive building extraction method based on original laser point cloud | |
| Dong et al. | Intelligent segmentation and measurement model for asphalt road cracks based on modified mask R-CNN algorithm | |
| CN108931206A (en) | Method for distinguishing is known for rail profile outlier detection and effective profile | |
| CN113779672B (en) | Rail profile abrasion calculation method, device, equipment and storage medium | |
| CN102201060B (en) | Method for tracking and evaluating nonparametric outline based on shape semanteme | |
| CN117162225B (en) | Demoulding and forming method and system for concrete prefabricated part | |
| CN102141385A (en) | Method for testing curved surface morphology of bituminous pavement | |
| CN104899592A (en) | Road semi-automatic extraction method and system based on circular template | |
| CN116701969A (en) | Processing and axis estimation method for air film Kong Dianyun | |
| Ma et al. | Surface roughness detection based on image analysis | |
| CN115661164A (en) | Disease three-dimensional information extraction method based on precise three-dimensional pavement |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |