CN114821522A - Urban road cross slope and super height value calculation method based on vehicle-mounted laser point cloud data - Google Patents
Urban road cross slope and super height value calculation method based on vehicle-mounted laser point cloud data Download PDFInfo
- Publication number
- CN114821522A CN114821522A CN202210321423.8A CN202210321423A CN114821522A CN 114821522 A CN114821522 A CN 114821522A CN 202210321423 A CN202210321423 A CN 202210321423A CN 114821522 A CN114821522 A CN 114821522A
- Authority
- CN
- China
- Prior art keywords
- point cloud
- point
- road
- points
- value
- 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
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
- G01S17/88—Lidar systems specially adapted for specific applications
- G01S17/89—Lidar systems specially adapted for specific applications for mapping or imaging
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Computer Networks & Wireless Communication (AREA)
- Electromagnetism (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Length Measuring Devices By Optical Means (AREA)
Abstract
The invention discloses an urban road cross slope and super-elevation calculation method based on vehicle-mounted laser point cloud data, which comprises the following steps: reconstructing the original point cloud data into an aligned scanning format grid; providing a layered road extraction algorithm to separate non-ground points including surrounding buildings, street lamps and signs, wherein the non-ground points include a height histogram, plane division and Euclidean clustering; establishing a point cloud space index, and removing ground object points including vehicles and vegetation running on a road by adopting a two-step self-adaptive height threshold method; filling the point cloud hole by adopting an interpolation algorithm based on a virtual grid; and extracting the cross section of the road according to a preset interval, and performing regression on the elevation and the transverse numerical value by adopting a random consistency sampling algorithm to obtain a transverse slope value. The method can completely extract the road part from the high-density laser point cloud data under the condition that a large amount of ground object noise exists, calculate the road cross slope and the ultrahigh value, and check the road section index safety.
Description
Technical Field
The invention belongs to the technical field of traffic safety analysis, and particularly relates to an urban road cross slope and ultrahigh value calculation method based on vehicle-mounted laser point cloud data.
Background
The cross slope is an important element of road geometric design. It is generally measured on a plane perpendicular to the road direction of travel, from the highest road surface center to both edges (on straight lines) or from the outer edge to the inner edge (on curves). The cross-slope measurements are typically made manually. These manual field measurements require engineers to place equipment on the road to acquire the cross slope. Figure 2 shows that three chinese engineers are using digital gradiometers to measure the cross slope of an established road. Furthermore, if traffic cones are placed around the measurement points, traffic will be affected. The field measurement is relatively time-consuming and labor-intensive, and the measurement of the large-scale road cross slope is difficult to be completed. Therefore, there is a need to propose a more automated and highly accurate method. Currently, mobile laser scanning systems are an emerging, promising surveying technology that integrates laser scanners, navigation sensors (global navigation satellite system (GNSS) and Inertial Measurement Unit (IMU)), and image data acquisition sensors on mobile platforms (panoramic and digital cameras). By continuous laser scanning, dense point clouds are collected on the road surface and roadside facilities as a vehicle carrying the point clouds travels along a given road. The point cloud data has higher accuracy and rich information inclusion, and is widely applied to detection and extraction of various road objects such as road surfaces, road marks, driving lines, road cracks, road manholes and the like.
However, although some studies have been done to evaluate the cross slope using point cloud data, there are some areas that could be improved. The noise discussed herein is primarily non-infrastructure noise, such as vehicles, pedestrians, and cyclists, rather than noise caused by airborne dust. Ideally, when the laser emitted by the laser radar is reflected from the road surface, a point cloud of the road surface can be obtained. And the existence of non-infrastructure noise enables the emitted laser to be reflected back on the surface of the vehicle halfway, the obtained point cloud is the non-infrastructure point cloud, and meanwhile, holes appear in the ground point cloud due to the shielding of the vehicle. These non-infrastructure noise and missing point regions caused by vehicle noise can adversely affect the fit of the cross-section. Collecting data at low traffic hours or combining multiple data sets may be effective in mitigating noise problems. However, both measures may increase the cost. Therefore, assessing the cross slope from noisy point cloud data is a more versatile way to reduce both the need to collect point cloud data at a particular time and the number of data collections.
The invention with patent number CN114170149A provides a method for extracting road geometric information based on laser point cloud, which includes: performing radius filtering and grid downsampling on the point cloud to simplify the point cloud; considering that the elevation distribution of the road surface points is concentrated and the road surface is smoother, extracting the elevation features and the local normal vector features of the point cloud to distinguish the ground points from the non-ground points; connecting the road surface points on the ground points by using a region growing method, and obtaining complete road surface point cloud after point deleting recovery by mistake; and finally, calculating a track vector and cutting a road cross section according to the track information of the collected vehicle, and obtaining the geometric parameters of the road by using a least square method. The method and the device have the advantages that the extraction precision and the program operation efficiency are considered simultaneously, and the automatic extraction of the road geometric information under different road environments can be considered comprehensively. However, the method needs a road point cloud extraction method by means of vehicle tracks and does not relate to the condition that a large amount of ground object noise exists.
Disclosure of Invention
The technical problem to be solved is as follows: aiming at the defects in the prior art, the invention provides the urban road cross slope and ultrahigh value calculation method based on the vehicle-mounted laser point cloud data, which can completely extract the road part from the high-density laser point cloud data in the presence of a large amount of ground noise, calculate the road cross slope and ultrahigh value and check the safety of road section indexes.
The technical scheme is as follows:
a method for calculating the cross slope and the ultrahigh value of an urban road based on vehicle-mounted laser point cloud data comprises the following steps:
s10, preprocessing the point cloud data, and reconstructing the original point cloud data into an aligned scanning format grid;
s20, segmenting the road point cloud by adopting a filtering technology based on a height histogram, a k-nearest neighbor and a clustering technology based on Euclidean distance, and separating non-ground points including surrounding buildings, street lamps and signs to obtain complete road point cloud;
s30, establishing a point cloud space index, and removing ground feature points including vehicles and vegetation running on the road by adopting a two-step self-adaptive height threshold method on the basis;
s40, filling the point cloud holes by adopting an interpolation algorithm based on a virtual grid to obtain complete road point cloud;
s50, extracting road cross sections according to a preset interval, and for each extracted cross section, performing regression on the elevation and the transverse numerical value by adopting a random consistency sampling algorithm to obtain a transverse slope value; and comparing the calculated cross slope value with the standard cross slope value, and determining the positions of the non-standard cross slope and the ultrahigh value.
Further, in step S10, the process of preprocessing the point cloud data and reconstructing the original point cloud data into an aligned scan format grid includes the following steps:
s11, performing coordinate conversion on the original laser point cloud data, and reconstructing the point cloud into a scanning format grid according to the corresponding relation between the original laser point cloud data and the track data;
s12, smoothing a plurality of parameters in the point cloud coordinate conversion by using a robust local weighted regression algorithm, and straightening curve segments in the road point cloud, wherein all the road point clouds are in a similar elevation range.
Further, in step S11, the process of reconstructing the point cloud into the scan format grid includes the following steps:
s111, dividing all laser points into scanning lines by using the time stamp of each laser point, wherein each scanning line has a corresponding track point; let T { T 1 ，T 2 ...T k ...T n |1≤k＜n，k，n∈N + Is the set of trace points, S { S } 1 ，S 2 ...S k ...S n |1≤k＜n，k，n∈N + LiDAR point sets for scan lines corresponding to the trajectory points;
s112, for the k scanning line, let T k Being the start of the trajectory vector, T k+1 As end point of the trajectory vector, forward vector
Is shown as
By T k As an origin, in
Establishing a local three-dimensional coordinate system for an X 'axis, wherein the Y' axis is orthogonal to the X 'axis along the horizontal direction, and a Z' axis is vertical to an upward vector of an X '-Y' plane; transforming geodesic points into local points in X ' Y ' Z ' space using a coordinate transformation matrix:
in the formula:
is a track point T k The coordinates of (a); (x) k y k z k ) T Is described by a geodetic coordinate system S k Point coordinates of (a); (x' k y′ k z′ k ) T Is described by a local coordinate system S k Point coordinates of (a); beta is a k Is the angle between the X axis and the ground plane; gamma ray k Is a vector
An angle to the X axis;
converting the original point cloud into a scanning mode grid, and starting an interval d; beta is a k And gamma k And calculating by using the coordinates of the adjacent path points to obtain:
further, in step S20, the process of obtaining the complete road point cloud includes the following steps:
s21: filtering ground objects higher than the road surface by adopting a height histogram method, equally dividing the elevation value range of all points according to the interval of 0.2m when reconstructing point clouds to obtain the number, index and vertical coordinate of each range of points as the histogram of the point number, finding the elevation value range corresponding to the highest strip section, and reserving the points with the elevation value near the elevation value range;
s22: filtering facilities on two sides of a road in a geometric range defining mode, reconstructing point clouds, and removing outliers on the sides of the region road in a region interest defining mode to obtain preliminarily defined road point clouds;
s23: establishing K neighbor point cloud indexes for the preliminarily defined road point clouds, and calculating the average distance of all the point clouds; traversing all the point clouds, searching K nearest neighbor points near each query point, and calculating the average distance between the query point and the K nearest neighbor points; if the average distance of the neighbor points of any query point is greater than the average distance of all point clouds, marking the query point as an outlier; after all points are traversed, all outliers in the road point cloud are removed;
s24: for the road point clouds with outliers removed, the point clouds with Euclidean distances of all the points smaller than a given threshold value are combined into a cluster, and the threshold value is empirically determined to be twice of the average distance of all the points; and after clustering is finished, all points in the point cloud have a cluster label, the cluster with the largest number of points is defined as the road point cloud, and the road point cloud is divided according to the label.
Further, in step S30, the multi-step adaptive height threshold feature filtering process includes the following steps:
s31, establishing a virtual grid, dividing a plane area where the discrete three-dimensional point cloud is located into a plurality of virtual grids with the same size by using the grids, wherein each virtual grid is equivalent to a subspace container of the point cloud space, and each laser point must fall into one of the grids;
s32, traversing all grids, and searching the lowest value z 'of the elevation in the grids' min Calculating average value z 'of point cloud elevation in grid network' mean Calculating the elevation fluctuation parameter delta z′ (ii) a Traversing point clouds in the grid network, and removing the elevation value larger than z' min +δ z′ The laser spot of (2);
s33, traversing all grids, and searching the lowest value z 'of the elevations of all laser points in 8 neighborhood grids around the grids' min2 Calculating the elevation fluctuation parameter omega at the moment z′ (ii) a Traversing point clouds in the grid network, and removing the elevation value larger than z' min2 +ω z′ Calculating the standard deviation of the elevation value of each point;
and S34, repeating the step S33 until the standard deviation change of two adjacent times is less than the set critical value, and stopping filtering.
Further, in step S31, the process of creating a virtual grid, and dividing the planar region where the discrete three-dimensional point cloud is located into a plurality of virtual grids with the same size by using the grid includes the following steps:
s311, establishing a virtual grid in Matlab environment, determining the size epsilon of the virtual grid, and creating a virtual grid of y' l -y′ r )/ε+1]×[x′ max /ε+1]The blank cell array of (1);
s312, traversing all the original laser point cloud data, solving the coordinate range of the original laser point cloud data in the XOZ plane, and putting the original laser point cloud data into a corresponding cellular array; any point (x' i ，y′ i ，z′ i ) The number of rows and columns of the virtual grid is as follows:
w＝[x′ i -x′ min /ε]+1
l＝[(y′ i -y′ min )/ε]+1
wherein w and l represent the row and column number, x 'of the virtual grid where the point is located' min 、y′ min Represents the minimum coordinate in the set of points, [.]And expressing rounding, wherein epsilon is the size of the virtual grid.
Further, for any point (x ') in the point cloud' i ，y′ i ，z′ i ) Parameters in adaptive altitude thresholding
And
the calculation method of (2) is as follows:
wherein epsilon is the size of the virtual grid,
is the height threshold for the initial filtering,
is the height threshold for the second filtering,
is the second filtering influence range.
Further, in step S40, the process of filling the point cloud hole by using the interpolation algorithm based on the virtual grid includes the following steps:
s41, extracting a road boundary point cloud;
s42, smoothing the road boundary by adopting a steady weighted local weighted regression algorithm to obtain complete road boundary point cloud;
s43, traversing all virtual grids, and if a certain grid is in the road boundary and has no point inside, regarding the region as a cavity region; searching a cavity area in a road range, and uniformly generating a plane two-dimensional coordinate of a point to be interpolated;
s44, calculating the vertical coordinates of the points to be interpolated, searching 8 non-blank grids around the grid where the points to be interpolated are located, determining the vertical coordinates of the points to be interpolated by adopting Delaunay triangular space interpolation, and completing filling of the holes after the vertical coordinates of the interpolation points in all the hole areas are calculated.
Further, in step S42, the process of obtaining a complete road boundary point cloud by smoothing the road boundary with a robust weighted local weighted regression algorithm includes the following steps:
s421, defining the window width of the filter, wherein the window width represents the proportion of the data points for calculating the smooth value to the total number of all the data points;
s422, traversing all the points, finding all the points in the given point window width, and calculating the weight values of all the adjacent points of the point, wherein the calculation method comprises the following steps:
wherein x represents the point to be smoothed, x i Representing the proximate point within the span x, dis being the horizontal distance from x to the farthest predicted value within the window width;
s423, obtaining a temporary smooth value x of x t Using calculated weights w i Performing weighted linear regression calculation on X to enable X to be t Is a temporary smoothed point set;
s424, calculating the point set X t Point-to-point difference R { R ] from X interior points 1 ，r 2 ，r 3 ...r k ...r n |1≤k＜n，k，n∈N + And on the basis of point-by-point difference, calculating the weight values of all the adjacent points in the window width again, wherein the calculation method comprises the following steps:
in the formula: r is i The residual error of the ith point is delta, and the median absolute deviation of the residual error is delta; when the residual size is 6 delta and the robustness weight is 0, relevant abnormal values are eliminated in the calculation process;
s425, using w i And
smoothing road boundaries。
Further, in step S50, the step of obtaining the cross slope value includes the steps of:
s51, along the road advancing direction, in the reconstruction space, the direction is consistent with the Y-axis direction, and sections with certain thickness are extracted at given intervals to be used as cross sections;
s52, for each extracted cross section, fitting the elevation and the transverse numerical value of the cross section point set by using a random sampling consistency algorithm, wherein the slope of the model is a transverse slope value and is represented as an ultrahigh value in a curve section;
and S53, referring to the current urban road design specification, comparing the specification value with the calculated value, and marking the station where the non-specification cross slope value is located.
Has the beneficial effects that:
the urban road cross slope and super-high value calculation method based on the vehicle-mounted laser point cloud data can accurately calculate the cross slope from dense MLS data under the condition that considerable vehicle noise and data difference exist. The method can be used to assess the cross-sectional condition of an established road. Is generally suitable for pavement maintenance and road reconstruction. By determining the position of the side slope which is out of the standard, the relevant road mechanism can be repaired in time, and accidents caused by water damage of the road surface are avoided. Therefore, the research result can also provide a basis for road widening in the case of considering road settlement or lack of design data. The method can accurately acquire the road cross section information. With advances in building information modeling and lidar technology, such road inventory can be imported into design and modeling platforms for full-life design and maintenance.
Drawings
Fig. 1 is a flow chart of an urban road cross slope and super-elevation calculation method based on vehicle-mounted laser point cloud data according to an embodiment of the invention.
Fig. 2 is a schematic diagram illustrating an effect of reconstructing an original point cloud into a scan format grid according to an embodiment of the present invention.
FIG. 3 is a schematic diagram illustrating the effect of extracting roads by height histogram and geometric partitioning according to an embodiment of the present invention; wherein, fig. 3(a) is a schematic diagram of reconstructing a spatial point cloud; fig. 3(b) is a schematic diagram of the road point cloud after division.
Fig. 4 is a schematic diagram illustrating a final effect of the road extraction method according to the embodiment of the present invention.
Fig. 5 is a schematic structural diagram of the Delaunay triangular space interpolation according to the embodiment of the present invention.
FIG. 6 is a schematic diagram illustrating a road void filling effect according to an embodiment of the present invention; wherein, fig. 6(a) is a schematic view of a road surface before filling a void; fig. 6(b) is a schematic diagram of the road after hole filling.
Fig. 7 is a schematic diagram of cross-sectional effects of RANSAC fitting roads according to an embodiment of the present invention.
Detailed Description
The following examples are presented to enable one of ordinary skill in the art to more fully understand the present invention and are not intended to limit the invention in any way.
FIG. 1 is a schematic structural diagram of an urban road cross slope and super-elevation calculation method based on vehicle-mounted laser point cloud data. The embodiment is suitable for the urban road cross slope and super-elevation calculation method based on vehicle laser point cloud data, and as shown in reference to fig. 1, the calculation method comprises the following steps:
and S10, preprocessing the point cloud data, and reconstructing the original point cloud data into an aligned scanning format grid.
And S20, segmenting the road point cloud by adopting a filtering technology based on a height histogram and a k-nearest neighbor and a clustering technology based on Euclidean distance, and separating most non-ground points (surrounding buildings, street lamps, signs and the like) to obtain complete road point cloud.
And S30, establishing a point cloud space index, and removing feature points (vehicles and vegetation running on the road) by adopting a two-step self-adaptive height threshold method on the basis.
And S40, filling the point cloud holes by adopting an interpolation algorithm based on a virtual grid to obtain complete road point cloud.
And S50, extracting road cross sections according to a given distance, and for each extracted cross section, performing regression on the elevation and the transverse numerical value by adopting a random consistency sampling algorithm to obtain a transverse slope value. And comparing the calculated cross slope value with the standard, and determining the positions of the non-standard cross slope and the ultrahigh value.
FIG. 2 is a schematic diagram of the effect of reconstructing an original point cloud into a scanning format grid according to the present invention. In one embodiment, the collected point cloud data preprocessing process in step S10 includes the following steps:
and S11, performing coordinate transformation on the original laser point cloud data, and reconstructing the point cloud into a scanning format grid according to the corresponding relation between the original laser point cloud data and the track data.
First, the points are divided into scan lines using the time stamp of each laser point. Each scan line has a corresponding trace point. Let T { T 1 ，T 2 ...T k ...T n |1≤k＜n，k，n∈N + Is the set of trace points. Let S { S 1 ，S 2 ...S k ...S n |1≤k＜n，k，n∈N + The LiDAR point set for the scan line that the trace point corresponds to. For the k-th scan line, let T k Being the start of the trajectory vector, T k+1 Is the end point of the trajectory vector. Thus, a forward vector
Is shown as
By T k As an origin, in
For the X' axis, a local three-dimensional coordinate system is established. The Y 'axis is orthogonal to the X' axis in the horizontal direction. The Z ' axis is then the vector oriented perpendicular to the X ' -Y ' plane. Transforming geodesic points into local points in X ' Y ' Z ' space using a coordinate transformation matrix:
the meaning of each parameter in the formula is as follows:
(x k y k z k ) T ＝S k The coordinates of points described by a geodetic coordinate system.
(x′ k y′ k z′ k ) T ＝S k Point coordinates described by a local coordinate system.
β k The angle between the X axis and the ground plane.
γ k Is a vector
And the X axis.
The original point cloud is converted to a scanning pattern grid, starting at interval d. Beta is a k And gamma k Is a key parameter in the meter and can be calculated using the coordinates of the neighboring path points:
s12, smoothing important parameters in the point cloud coordinate conversion by using a steady local weighted regression algorithm to ensure the reconstruction effect, straightening curve segments of the reconstructed road point cloud, and basically keeping all the road point clouds in a similar elevation range.
FIG. 3 is a schematic diagram illustrating the effect of extracting roads by using height histogram and geometric partitioning according to the present invention; wherein, fig. 3(a) is a schematic diagram of reconstructing a spatial point cloud; fig. 3(b) is a schematic diagram of the road point cloud after division. Fig. 4 is a schematic diagram of the final effect of the road extraction method of the present invention.
In one embodiment, in step S20, the process of performing road surface identification segmentation processing on the point cloud data includes the following steps:
s21: and (3) filtering ground objects higher than the road surface by adopting a height histogram method: when point clouds are reconstructed, the elevation value ranges of all the points are equally divided according to the interval of 0.2m, the number and the index of the points in each range can be obtained, the vertical coordinate is a histogram of the number of the points, the elevation value range corresponding to the highest strip section is found, and the points with the elevation values near the range are reserved.
S22: most facilities on both sides of the road are filtered out in a geometric range defining mode: after the point cloud is reconstructed, the fluctuation of the road point cloud in the transverse range (namely the Y value) is not large, so that the outliers on the road side can be removed in a mode of dividing the region of interest.
S23: performing K-neighbor based filtering: and for the preliminarily defined road point cloud, establishing a K neighbor point cloud index, calculating the average distance of all the point clouds, traversing all the point clouds, searching K nearest neighbor points nearby the point clouds, calculating the average distance between a query point and K neighbors, marking the point as an outlier if the average distance of the neighbor point of the point is greater than the average distance of all the point clouds, and removing all the outliers after the traversal of all the points is finished.
S24: performing Euclidean distance clustering-based ground point cloud segmentation: and (3) for the road point clouds with the outliers removed, forming a cluster by the point clouds with Euclidean distances of all the points smaller than a given threshold value, wherein the threshold value is empirically determined to be twice the average distance of all the points, after the clustering is completed, all the points in the point clouds have a cluster label, and the cluster with the largest number of points is the road point cloud and is segmented according to the label.
Fig. 4 is a schematic diagram of the final effect of the road extraction according to the present invention. In step S24 of this embodiment, the processing procedure of the euclidean clustering method includes the following steps:
s241: find a point p in space 11 Finding n points nearest to the KD-Tree by using KD-Tree, and judging the n points to p 11 Euclidean distance of. Points p with a distance less than a threshold r 12 ，p 13 ，p 14 .., placed in class Q.
S242: at Q (p) 11 ) Find a point p 12 Step S241 is repeated.
S243: at Q (p) 11 ，p 12 ) Find a point, repeat step S241, find p 22 ，p 23 ，p 24 .., all put into Q.
S244: when Q can not add new points any more, the search is completed.
In one embodiment, in step S30, the adaptive height threshold feature filtering process includes the following steps:
s31: establishing a virtual grid: dividing a plane area where the scattered three-dimensional point cloud is located into a plurality of virtual grids with the same size by using grids, wherein each virtual grid is equivalent to a subspace container of the point cloud space, and each laser point must fall into one of the grids.
Specifically, to build the virtual mesh in a Matlab environment, the virtual mesh size ε is first determined and a value of y 'is created' l -y′ r )/ε+1]×[x′ max /ε+1]And then traversing all the original laser point cloud data, solving the coordinate range of the original laser point cloud data in the XOZ plane, and putting the original laser point cloud data into the corresponding cell array. Any point (x' i ，y′ i ，z′ i ) The number of rows and columns of the virtual grid is as follows:
w＝[x′ i -x′ min /ε]+1
l＝[(y′ i -y′ min )/ε]+1
the meaning of each parameter in the formula is as follows: w and l represent the row and column number, x 'of the virtual grid where the point is located' min 、y′ min Represents the minimum coordinate in the set of points, [.]Indicating rounding.
S32: performing a first step of height threshold terrain filtering: traversing all grids, and searching the lowest value z 'of the elevation in the grids' min Calculating average value z 'of point cloud elevation in grid network' mean Calculating the elevation fluctuation parameter delta z′ . Traversing point clouds in the grid network, and removing the elevation value larger than z' min +δ z′ The laser spot of (1).
And (3) performing height threshold terrain filtering in the second step: traversing all grids, and searching the lowest value z 'of the elevations of all laser points in 8 neighborhood grids around the grids' min Calculating the elevation fluctuation parameter omega at the moment z′ . Traversing point clouds in the grid network, and removing the elevation value larger than z' min +ω z′ And calculating the standard deviation of the elevation value of each point position, and comparing the standard deviation with the last execution result. When the standard deviation variation is below criticalAt value (the empirical value in this study was set to 0.2), the filtering was stopped.
The design of the algorithm takes into account three factors: (1) the road point cloud has smoothness, and the fluctuation of the elevation of the road point cloud is very small in a certain range; furthermore, the road points are continuous in cross section, and their elevation is not significantly increased or decreased except for intermediate obstacles and roadside obstacles. (2) The road alignment is converted into a straight line at S20, eliminating the adverse effect of the curve segment in the mesh generation process. For any point (x ') in the point cloud' i ，y′ i ，z′ i ) Parameters in adaptive altitude thresholding
And
the calculation method of (2) is as follows:
in one embodiment, S40, the filling the point cloud hole by the virtual grid-based interpolation algorithm includes the following steps:
and S41, extracting the road boundary point cloud.
The searching and filling of the hollow holes mainly aim at the area in the road boundary, and the hollow holes outside the road are not in the subsequent hollow hole filling and cross section calculation range. Therefore, extraction of the road boundary is required. And dividing the road surface in the reconstructed scene into strip-shaped units, wherein the divided size belongs to the E. And searching to obtain the leftmost point and the rightmost point of each strip-shaped unit respectively. Considering that the calculations within the bar cells at this stage are independent and similar to each other, parallel calculations can be used to improve the practical operating efficiency.
And S42, smoothing the road boundary by using a robust local weighted linear regression algorithm (RLWLR) to obtain a complete road boundary point cloud.
After the processing of step S41, the left and right boundaries of the road are basically extracted. However, because the holes which are partially and incompletely filled exist on the road surface, deviation points exist in the extracted road boundary point cloud. Therefore, it is necessary to obtain a boundary curve point that is closer to the road surface boundary by using a smoothing algorithm.
In step S42 of this embodiment, the processing procedure of the robust local weighted linear regression algorithm includes the following steps:
s421, a window width of the filter is defined, which represents a ratio of the number of data points used for calculating the smoothing value to the total number of data points.
S422, traversing all the points, finding all the points in the given point window width, and calculating the weight values of all the adjacent points of the point, wherein the calculation method comprises the following steps:
in the formula: x denotes the point to be smoothed, x i Representing the near point within the span x, dis is the horizontal distance from x to the farthest predicted value within the window width.
S423, obtaining a temporary smooth value x of x t Weight w calculated using the above formula i Performing weighted linear regression calculation on X to enable X to be t Is a temporary set of smoothed points.
S424, calculating the point set X t Point-to-point difference (i.e., residual) R { R } from points within X 1 ，r 2 ，r 3 ...r k ...r n |1≤k＜n，k，n∈N + Based on the difference, the weight values of all the adjacent points in the window width are calculated again, and the calculation method is as follows:
in the formula: r is i Is the residual error of the ith point, and is the median absolute deviation of the residual error.
When the residual is greater than 6 δ, the robustness weight is 0, and the related outliers will be eliminated in the calculation process.
S424, smoothed data usage w i And
to ensure the smoothing effect, steps S422 and S423 are repeated a plurality of times. Empirically, a window width of 20 was used.
And S43, searching the hollow area in the road range, and uniformly generating the plane two-dimensional coordinates of the points to be interpolated.
And traversing all the virtual grids, and if a certain grid is in the road boundary and has no point inside, regarding the region as a cavity region. And after the cavity area is determined, uniformly generating two-dimensional points in the virtual grid, namely the plane two-dimensional coordinates of the points to be interpolated.
S44, calculating the vertical coordinates of the points to be interpolated, searching the points in 8 non-space grids around the grid where the points to be interpolated are located, determining the vertical coordinates of the points to be interpolated by using Delaunay triangular space interpolation with the points in the point set as reference points, and completing filling of the cavity after the vertical coordinates of the interpolation points in all the cavity areas are calculated.
FIG. 5 is a schematic diagram of the Delaunay triangle spatial interpolation structure of the present invention. FIG. 6 is a schematic diagram of the road surface void filling effect of the present invention. Wherein, fig. 6(a) is a schematic view of a road surface before filling a void; fig. 6(b) is a schematic diagram of the road after hole filling. The invention adopts a linear interpolation method to complete the filling work of the cavity. The size of the cavity area is related to the size of vehicles on the road and the angle of the laser beam, and according to the empirical evaluation of different laser point cloud data sets, the size of the cavity is generally approximate to a rectangular area with the length of 2-5 m and the width of 7-15 m. Therefore, the range of the reference point should be larger than the rectangular area to ensure that the center point of the hole can also find the reference point. As shown in FIG. 6, let [ x' j ，y′ j ]Is the grid coordinate of the jth query point. Let [ x' j1 ，y′ j1 ，z′ j1 ]，[x′ j2 ，y′ j2 ，z′ j2 ]Is [ x' j3 ，y′ j3 ，z′ j3 ]And performing linear interpolation fitting on the points to be interpolated by using a Delaunay triangular space interpolation algorithm to obtain the elevation values of the points to be interpolated. The calculation process is as follows:
in the formula: z' j Is the elevation value of the point to be interpolated.
In one embodiment, S50, the cross section extraction and the cross slope value calculation include the following steps:
s51, extracting cross sections at given intervals, along the advancing direction of the road, in the reconstruction space, consistent with the Y-axis direction, and extracting sections with certain thickness at given intervals, namely the cross sections;
in the step S20 point cloud preprocessing, the complex road line shape is converted into a straight line in a point cloud reconstruction mode. At this time, the road advancing direction is the same as the x' axis direction, and it is simpler to acquire a section perpendicular to the road advancing direction than the known method. Theoretically all abscissas are x' i The plane formed by the points of (a) is the cross section there. However, considering that the distribution of dots in the dot cloud is not continuous, a region [ x 'of a certain width may be used' i -d，x′ i +d]And (5) extracting the cross section. It is contemplated that in field cross-sectional measurements are typically made at 10-30m intervals. Therefore, the cross sections are taken at intervals of 10m, and the width of the cross section is 0.4 m. Since the width value extracted from the cross section may affect the range of the data interval for performing the cross slope calculation, the selection of the width value may affect the accuracy of the cross slope calculation. To determine the optimal value for this value, the present invention has performed an allergy analysis on a plurality of measured road segments.
And S52, calculating a cross slope value and an ultrahigh value, and fitting the elevation and the transverse value of the cross section point set by using a random sampling consistency algorithm for each extracted cross section, wherein the slope of the model is the cross slope value (the ultrahigh value in a curve segment).
And according to the definition of the transverse slope, performing linear regression on the elevation and the transverse numerical value of the laser point in each cross-sectional point set to obtain the transverse slope value. When the road surface is extracted, all the ground object point clouds cannot be completely segmented, so a small amount of noise points still exist in the extracted cross section. Least Squares (LS) is a simple fitting algorithm whose principle is model solution based on mean square error minimization. In linear regression, the least squares approach is to try to find a straight line, minimizing the sum of the euclidean distances from all the sample points to the straight line. However, this algorithm adapts to all points as much as possible, resulting in a poor fitting accuracy in the presence of a large number of noise points.
The invention adopts a Random Sample Consensus (RANSAC) method to calculate the corresponding cross slope value of each cross section. The random sample consensus algorithm solves the model from a set of data containing outliers by using an iterative method. The idea of the algorithm is as follows: 1) selecting SS points as sample points randomly, and fitting the model; 2) finding out points in the distance fitting straight line tolerance range MD, and counting the number of the points; 3) then, randomly selecting MD points, and repeating the steps 2) and 3) until the iteration is finished; 4) and finding out the condition that the data points are the most after certain fitting, namely the model for solving.
At this stage, the user needs to manually specify the minimum sample size SS and the tolerance range MD. Since the model in this example is a straight line, it is sufficient to set SS to 2. Based on empirical evaluation of two different laser point cloud data sets, any point having a distance of more than 0.2m from the fitting model is regarded as an outlier (MD ═ 0.2m, distance metric is the square of the euclidean distance, and the maximum number of iterations is 1000). Fig. 7 gives 4 examples of fitting a straight line with RANSAC. As can be seen from fig. 7, the fitted straight line is very close to the fitting of the point set, and the RANSAC method has better fitting effect on the points than the LS method.
And S53, evaluating the safety of the cross slope value, comparing the standard value with the calculated value according to the current urban road design standard, and marking the station where the non-standard cross slope value is located.
The above is only a preferred embodiment of the present invention, and the protection scope of the present invention is not limited to the above-mentioned embodiments, and all technical solutions belonging to the idea of the present invention belong to the protection scope of the present invention. It should be noted that modifications and embellishments within the scope of the invention may be made by those skilled in the art without departing from the principle of the invention.
Claims (10)
1. A method for calculating the cross slope and the ultrahigh value of an urban road based on vehicle-mounted laser point cloud data is characterized by comprising the following steps:
s10, preprocessing the point cloud data, and reconstructing the original point cloud data into an aligned scanning format grid;
s20, segmenting the road point cloud by adopting a filtering technology based on a height histogram, a k-nearest neighbor and a clustering technology based on Euclidean distance, and separating non-ground points including surrounding buildings, street lamps and signs to obtain complete road point cloud;
s30, establishing a point cloud space index, and removing ground feature points including vehicles and vegetation running on the road by adopting a two-step self-adaptive height threshold method on the basis;
s40, filling the point cloud holes by adopting an interpolation algorithm based on a virtual grid to obtain complete road point cloud;
s50, extracting road cross sections according to a preset interval, and for each extracted cross section, performing regression on the elevation and the transverse numerical value by adopting a random consistency sampling algorithm to obtain a transverse slope value; and comparing the calculated cross slope value with the standard cross slope value, and determining the positions of the non-standard cross slope and the ultrahigh value.
2. The urban road cross slope and ultra-high value calculation method based on vehicle-mounted laser point cloud data as claimed in claim 1, wherein in step S10, the point cloud data is preprocessed, and the process of reconstructing the original point cloud data into the aligned scan format grid comprises the following steps:
s11, carrying out coordinate conversion on the original laser point cloud data, and reconstructing the point cloud into a scanning format grid according to the corresponding relation between the original laser point cloud data and the track data;
s12, smoothing a plurality of parameters in the point cloud coordinate conversion by using a robust local weighted regression algorithm, and straightening curve segments in the road point cloud, wherein all the road point clouds are in a similar elevation range.
3. The method for calculating the urban road cross slope and the super elevation based on the vehicle-mounted laser point cloud data according to claim 2, wherein the step S11 of reconstructing the point cloud into the scan format grid comprises the following steps:
s111, dividing all laser points into scanning lines by using the time stamp of each laser point, wherein each scanning line has a corresponding track point; let T { T 1 ，T 2 ...T k ...T n |1≤k＜n，k，n∈N + Is the set of trace points, S { S } 1 ，S 2 ...S k ...S n |1≤k＜n，k，n∈N + LiDAR point sets for scan lines corresponding to the trajectory points;
s112, for the k scanning line, let T k Being the start of the trajectory vector, T k+1 As end point of the trajectory vector, forward vector
Is shown as
By T k As an origin, in
Establishing a local three-dimensional coordinate system for an X 'axis, wherein the Y' axis is orthogonal to the X 'axis along the horizontal direction, and a Z' axis is vertical to an upward vector of an X '-Y' plane; transforming geodesic points into local points in X ' Y ' Z ' space using a coordinate transformation matrix:
in the formula:
is a track point T k The coordinates of (a); (x) k y k z k ) T Is described by a geodetic coordinate system S k Point coordinates of (a); (x' k y′ k z′ k ) T Is described by a local coordinate system S k Point coordinates of (a); beta is a k Is the angle between the X axis and the ground plane; gamma ray k Is a vector
An angle to the X axis;
converting the original point cloud into a scanning mode grid, and starting an interval d; beta is a k And gamma k And calculating by using the coordinates of the adjacent path points to obtain:
4. the method for calculating the urban road cross slope and the ultra-high value based on the vehicle-mounted laser point cloud data as claimed in claim 1, wherein in the step S20, the process of obtaining the complete road point cloud comprises the following steps:
s21: filtering ground objects higher than the road surface by adopting a height histogram method, equally dividing the elevation value range of all points according to the interval of 0.2m when reconstructing point clouds to obtain the number, index and vertical coordinate of each range of points as the histogram of the point number, finding the elevation value range corresponding to the highest strip section, and reserving the points with the elevation value near the elevation value range;
s22: filtering facilities on two sides of a road in a geometric range defining mode, reconstructing point clouds, and removing outliers on the sides of the region road in a region interest defining mode to obtain preliminarily defined road point clouds;
s23: establishing K neighbor point cloud indexes for the preliminarily defined road point clouds, and calculating the average distance of all the point clouds; traversing all the point clouds, searching K nearest neighbor points near each query point, and calculating the average distance between the query point and the K nearest neighbor points; if the average distance of the neighbor points of any query point is greater than the average distance of all point clouds, marking the query point as an outlier; after all points are traversed, all outliers in the road point cloud are removed;
s24: for the road point clouds with outliers removed, the point clouds with Euclidean distances of all the points smaller than a given threshold value are combined into a cluster, and the threshold value is empirically determined to be twice of the average distance of all the points; and after clustering is finished, all points in the point cloud have a cluster label, the cluster with the largest number of points is defined as the road point cloud, and the road point cloud is divided according to the label.
5. The method for calculating the cross slope and the ultra-high value of the urban road based on the vehicle-mounted laser point cloud data as claimed in claim 1, wherein in the step S30, the multi-step adaptive height threshold feature filtering process comprises the following steps:
s31, establishing a virtual grid, dividing a plane area where the discrete three-dimensional point cloud is located into a plurality of virtual grids with the same size by using the grids, wherein each virtual grid is equivalent to a subspace container of the point cloud space, and each laser point must fall into one of the grids;
s32, traversing all grids, and searching the lowest value z 'of the elevation in the grids' min Calculating the average value z 'of the point cloud elevation in the grid network' mean Calculating the elevation fluctuation parameter delta z′ (ii) a Traversing point clouds in the grid network, and removing the elevation value larger than z' min +δ z′ The laser spot of (2);
s33, traversing all grids, and searching the lowest value z 'of the elevations of all laser points in 8 neighborhood grids around the grids' min2 Calculating the elevation fluctuation parameter omega at the moment z′ (ii) a Traversing point clouds in the grid network, and removing the elevation value larger than z' min2 +ω z′ Calculating the standard deviation of the elevation value of each point;
and S34, repeating the step S33 until the standard deviation change of two adjacent times is less than the set critical value, and stopping filtering.
6. The method for calculating the urban road cross slope and the ultrahigh value based on the vehicle-mounted laser point cloud data as claimed in claim 5, wherein in the step S31, a virtual grid is established, and the process of dividing the plane area where the dispersed three-dimensional point cloud is located into a plurality of virtual grids with consistent sizes by using the grid comprises the following steps:
s311, establishing a virtual grid in Matlab environment, determining the size epsilon of the virtual grid, and creating a virtual grid of y' l -y′ r )/ε+1]×[x′ max /ε+1]The blank cell array of (1);
s312, traversing all the original laser point cloud data, solving the coordinate range of the original laser point cloud data in the XOZ plane, and putting the original laser point cloud data into a corresponding cellular array; any point (x' i ，y′ i ，z′ i ) The number of rows and columns of the virtual grid is as follows:
w＝[x′ i -x′ min /ε]+1
l＝[(y′ i -y′ min )/ε]+1
wherein w and l represent the row and column number, x 'of the virtual grid where the point is located' min 、y′ min Represents the minimum coordinate in the set of points, [.]And expressing rounding, wherein epsilon is the size of the virtual grid.
7. The method for calculating the urban road cross slope and the ultrahigh value based on the vehicle-mounted laser point cloud data as claimed in claim 5, wherein the method is applied to any point (x' i ，y′ i ，z′ i ) Parameters in adaptive altitude thresholding
And
the calculation method of (2) is as follows:
wherein epsilon is the size of the virtual grid,
is the height threshold for the initial filtering,
is the height threshold for the second filtering,
is the second filtering influence range.
8. The method for calculating the urban road cross slope and the ultrahigh value based on the vehicle-mounted laser point cloud data according to claim 1, wherein in the step S40, the process of filling the point cloud hole by using the interpolation algorithm based on the virtual grid comprises the following steps:
s41, extracting a road boundary point cloud;
s42, smoothing the road boundary by adopting a steady weighted local weighted regression algorithm to obtain complete road boundary point cloud;
s43, traversing all virtual grids, and if a certain grid is in the road boundary and has no point inside, regarding the region as a cavity region; searching a cavity area in a road range, and uniformly generating a plane two-dimensional coordinate of a point to be interpolated;
s44, calculating the vertical coordinates of the points to be interpolated, searching 8 non-blank grids around the grid where the points to be interpolated are located, determining the vertical coordinates of the points to be interpolated by adopting Delaunay triangular space interpolation, and completing filling of the holes after the vertical coordinates of the interpolation points in all the hole areas are calculated.
9. The method for calculating the urban road cross slope and the ultrahigh value based on the vehicle-mounted laser point cloud data as claimed in claim 8, wherein in the step S42, the road boundary is smoothed by adopting a robust weighted local weighted regression algorithm, and the process of obtaining the complete road boundary point cloud comprises the following steps:
s421, defining the window width of the filter, wherein the window width represents the proportion of the data points for calculating the smooth value to the total number of all the data points;
s422, traversing all the points, finding all the points in the given point window width, and calculating the weight values of all the adjacent points of the point, wherein the calculation method comprises the following steps:
wherein x represents a point to be smoothed, x i Representing the proximate point within the span x, dis being the horizontal distance from x to the farthest predicted value within the window width;
s423, obtaining a temporary smooth value x of x t Using calculated weights w i Performing weighted linear regression calculation on X to enable X to be t Is a temporary smoothed point set;
s424, calculating the point set X t Point-to-point difference R { R ] from X interior points 1 ，r 2 ，r 3 ...r k ...r n |1≤k＜n，k，n∈N + And on the basis of point-by-point difference, calculating the weight values of all the adjacent points in the window width again, wherein the calculation method comprises the following steps:
in the formula: r is i The residual error of the ith point is delta, and the median absolute deviation of the residual error is delta; when the residual size is 6 delta and the robustness weight is 0, relevant abnormal values are eliminated in the calculation process;
s425, using w i And
and smoothing the road boundary.
10. The method for calculating the cross slope and the ultra-high value of the urban road based on the vehicle-mounted laser point cloud data as claimed in claim 1, wherein in the step S50, the step of obtaining the cross slope value comprises the following steps:
s51, along the road advancing direction, in the reconstruction space, the direction is consistent with the Y-axis direction, and sections with certain thickness are extracted at given intervals to be used as cross sections;
s52, for each extracted cross section, fitting the elevation and the transverse numerical value of the cross section point set by using a random sampling consistency algorithm, wherein the slope of the model is a transverse slope value and is represented as an ultrahigh value in a curve section;
and S53, referring to the current urban road design specification, comparing the specification value with the calculated value, and marking the station where the non-specification cross slope value is located.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210321423.8A CN114821522A (en) | 2022-03-29 | 2022-03-29 | Urban road cross slope and super height value calculation method based on vehicle-mounted laser point cloud data |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210321423.8A CN114821522A (en) | 2022-03-29 | 2022-03-29 | Urban road cross slope and super height value calculation method based on vehicle-mounted laser point cloud data |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN114821522A true CN114821522A (en) | 2022-07-29 |
Family
ID=82533125
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210321423.8A Pending CN114821522A (en) | 2022-03-29 | 2022-03-29 | Urban road cross slope and super height value calculation method based on vehicle-mounted laser point cloud data |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114821522A (en) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116895022A (en) * | 2023-09-11 | 2023-10-17 | 广州蓝图地理信息技术有限公司 | Building boundary extraction method based on point cloud data processing |
| CN117392632A (en) * | 2023-12-11 | 2024-01-12 | 中交第二公路勘察设计研究院有限公司 | Road element change monitoring method and device |
-
2022
- 2022-03-29 CN CN202210321423.8A patent/CN114821522A/en active Pending
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116895022A (en) * | 2023-09-11 | 2023-10-17 | 广州蓝图地理信息技术有限公司 | Building boundary extraction method based on point cloud data processing |
| CN116895022B (en) * | 2023-09-11 | 2023-12-01 | 广州蓝图地理信息技术有限公司 | Building boundary extraction method based on point cloud data processing |
| CN117392632A (en) * | 2023-12-11 | 2024-01-12 | 中交第二公路勘察设计研究院有限公司 | Road element change monitoring method and device |
| CN117392632B (en) * | 2023-12-11 | 2024-03-15 | 中交第二公路勘察设计研究院有限公司 | Road element change monitoring method and device |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110570428B (en) | Method and system for dividing building roof sheet from large-scale image dense matching point cloud | |
| Yan et al. | Automated extraction of structural elements in steel girder bridges from laser point clouds | |
| Rodriguez-Cuenca et al. | An approach to detect and delineate street curbs from MLS 3D point cloud data | |
| CN109633674A (en) | Three-dimensional Track automatic planning is maked an inspection tour in transmission of electricity based on laser point cloud data | |
| CN106780524A (en) | A kind of three-dimensional point cloud road boundary extraction method | |
| CN114821522A (en) | Urban road cross slope and super height value calculation method based on vehicle-mounted laser point cloud data | |
| CN106401643A (en) | Tunnel back-break detection method based on laser-point cloud | |
| CN111325138B (en) | Road boundary real-time detection method based on point cloud local concave-convex characteristics | |
| CN103017739A (en) | Manufacturing method of true digital ortho map (TDOM) based on light detection and ranging (LiDAR) point cloud and aerial image | |
| CN104376595A (en) | Method for generating three-dimensional roads on basis of airborne LiDAR (light detection and ranging) and GIS (geographic information system) cooperation | |
| CN114549879B (en) | Target identification and central point extraction method for tunnel vehicle-mounted scanning point cloud | |
| Cheng et al. | Three-dimensional reconstruction of large multilayer interchange bridge using airborne LiDAR data | |
| CN113009453B (en) | Mine road edge detection and mapping method and device | |
| CN111652241B (en) | Building contour extraction method integrating image features and densely matched point cloud features | |
| CN112561944A (en) | Lane line extraction method based on vehicle-mounted laser point cloud | |
| Zhang et al. | A patch-based method for the evaluation of dense image matching quality | |
| Zhou et al. | Highway alignments extraction and 3D modeling from airborne laser scanning point clouds | |
| CN114842450A (en) | Driving region detection method, device and equipment | |
| CN115267815B (en) | Road side laser radar group optimization layout method based on point cloud modeling | |
| CN112184725A (en) | Structured light strip center extraction method for asphalt pavement image | |
| Liu et al. | Detection and reconstruction of static vehicle-related ground occlusions in point clouds from mobile laser scanning | |
| CN114842166A (en) | Negative obstacle detection method, system, medium, and apparatus applied to structured road | |
| Xu et al. | A method of 3d building boundary extraction from airborne lidar points cloud | |
| CN117036326A (en) | Defect detection method based on multi-mode fusion | |
| CN109697418B (en) | Post-processing method for extracting image for remote sensing image road network for scene restoration |
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 |