Calculation method of cross slope and superelevation value of urban roads based on vehicle-borne laser point cloud data

Abstract

The present invention discloses a method for calculating the cross slope and superelevation of urban roads based on vehicle-mounted laser point cloud data, including: reconstructing the original point cloud data into an aligned scanning format grid; proposing a hierarchical road extraction algorithm to separate non-ground points including surrounding buildings, street lights and signs, including height histogram, plane division and Euclidean clustering; establishing a point cloud spatial index, using a two-step adaptive height threshold method to remove ground object points including vehicles and vegetation on the road; using an interpolation algorithm based on a virtual grid to fill point cloud holes; extracting road cross sections at preset intervals, using a random consistency sampling algorithm to regress elevation and lateral values, and obtaining cross slope values. The present invention can completely extract a road section from high-density laser point cloud data in the presence of a large amount of ground object noise, and calculate the road cross slope and superelevation value, and verify the safety of road section indicators.

Description

Technical Field

[0001] The invention belongs to the technical field of traffic safety analysis, and in particular relates to a method for calculating cross slope and superelevation value of urban roads based on vehicle-mounted laser point cloud data. Background Art

[0002] Cross slope is an important element of road geometry design. It is usually measured in a plane perpendicular to the direction of travel of the road, from the highest center of the pavement to both edges (on a straight section) or from the outer edge to the inner edge (on a curve). Cross slope measurements are usually made manually. These manual field measurements require engineers to place equipment on the pavement to obtain the cross slope. Figure 2 The video shows three Chinese engineers using a digital level to measure the cross slope of a built road. In addition, if traffic cones are set around the measurement point, it will affect traffic. On-site measurement is relatively time-consuming and labor-intensive, and it is difficult to complete large-scale road cross slope measurements. Therefore, a more automated and high-precision method is needed. At present, the mobile laser scanning system is an emerging and promising measurement technology that integrates laser scanners, navigation sensors (global navigation satellite system (GNSS) and inertial measurement unit (IMU)), and image data acquisition sensors (panoramic and digital cameras) on mobile platforms. Through continuous laser scanning, dense point clouds on the road surface and roadside facilities are collected while the vehicle loaded with the point cloud travels along a given road. Due to the high accuracy and rich information inclusiveness of point cloud data, it is widely used in the detection and extraction of various road objects such as road surfaces, road markings, lanes, road cracks, and road manholes.

[0003] However, although some studies have used point cloud data to evaluate cross slope, there are still some areas for improvement. The noise discussed in this paper is mainly non-infrastructure noise, such as vehicles, pedestrians and cyclists, rather than noise caused by airborne dust. Ideally, when the laser emitted by the lidar is reflected from the road surface, a point cloud of the road surface can be obtained. However, the presence of non-infrastructure noise causes the emitted laser to reflect back from the vehicle surface halfway, and the point cloud obtained is a non-infrastructure point cloud. At the same time, holes appear in the ground point cloud due to vehicle occlusion. These non-infrastructure noises and missing point areas caused by vehicle noise will have an adverse effect on the fitting of the cross section. Collecting data during low traffic hours or merging multiple data sets can effectively alleviate the noise problem. However, these two measures may increase costs. Therefore, evaluating cross slope from noisy point cloud data is a more widely applicable way, which can reduce the need to collect point cloud data at a specific time and reduce the amount of data collection.

[0004] The invention with patent number CN114170149A mentions a method for extracting road geometry information based on laser point cloud, including: performing radius filtering and grid downsampling on the point cloud to simplify the point cloud; taking into account that the elevation distribution of road surface points is relatively concentrated and the road surface is smoother, the elevation features and local normal vector features of the point cloud are extracted to distinguish ground points from non-ground points; the road surface points are connected using the regional growing method on the ground points, and a complete road surface point cloud is obtained after the points are restored after being deleted by mistake; finally, based on the collected vehicle trajectory information, the trajectory vector is calculated and the road cross section is cut, and the geometric parameters of the road are obtained using the least squares method. The present invention takes into account both extraction accuracy and program operation efficiency, and can comprehensively consider the automated extraction of road geometry information under different road environments. However, the invention requires the use of vehicle trajectories, and does not involve a road point cloud extraction method in the presence of a large amount of ground noise. Summary of the invention

[0005] Technical problem to be solved: In view of the deficiencies in the prior art, the present invention provides a method for calculating the cross slope and superelevation of urban roads based on vehicle-mounted laser point cloud data. The method can completely extract a road section from high-density laser point cloud data in the presence of a large amount of ground noise, calculate the cross slope and superelevation of the road, and verify the safety of road section indicators.

[0006] Technical solution:

[0007] A method for calculating the cross slope and superelevation value of an urban road based on vehicle-mounted laser point cloud data, the method comprising the following steps:

[0008] S10, preprocessing the point cloud data to reconstruct the original point cloud data into an aligned scan format grid;

[0009] S20, uses height histogram-based, k-nearest-neighbor-based filtering and Euclidean distance-based clustering technology to segment the road point cloud, separates non-ground points including surrounding buildings, street lights and signs, and obtains a complete road point cloud;

[0010] S30, establishing a point cloud spatial index, based on which a two-step adaptive height threshold method is used to remove ground object points including vehicles on the road and vegetation;

[0011] S40, uses a virtual grid-based interpolation algorithm to fill the point cloud holes to obtain a complete road point cloud;

[0012] S50, extracting road cross sections at preset intervals, and for each extracted cross section, regressing the elevation and lateral values ​​using a random consistency sampling algorithm to obtain a cross slope value; comparing the calculated cross slope value with the standard cross slope value to determine the location of the cross slope that does not meet the standard and the superelevation value.

[0013] Furthermore, 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:

[0014] S11, coordinate conversion is performed on the original laser point cloud data, and the point cloud is reconstructed into a scanning format grid according to the corresponding relationship between the original laser point cloud data and the trajectory data;

[0015] S12, using a robust local weighted regression algorithm to smooth multiple parameters in the point cloud coordinate transformation, straightening the curve segments in the road point cloud, and all road point clouds are in a similar elevation range.

[0016] Furthermore, in step S11, the process of reconstructing the point cloud into a scan format grid includes the following steps:

[0017] S111, using the timestamp of each laser point, divide all laser points into scan lines, each scan line has a corresponding trajectory point; let T{T1, T2...T k ...T n |1≤k<n,k,n∈N +} is a set of trajectory points, S{S1, S2...S k ...S n |1≤k<n,k,n∈N +} is the LiDAR point set of the scan line corresponding to the trajectory point;

[0018] S112, for the kth scan line, let T k is the starting point of the trajectory vector, T k+1 is the end point of the trajectory vector, the positive vector Expressed as T k As the origin, The local three-dimensional coordinate system is established with the X' axis as the X' axis, the Y' axis is orthogonal to the X' axis in the horizontal direction, and the Z' axis is perpendicular to the upward vector of the X'-Y' plane; the geodesic points are transformed into local points in the X'Y'Z' space using the coordinate transformation matrix:

[0019]

[0020] Where: is the trajectory point T k The coordinates of (x k y k z k ) T S is described by the geodetic coordinate system k The coordinates of the point; (x′ k y′ k z′ k )T S is described by the local coordinate system k The point coordinates of β k is the angle between the X-axis and the geodetic plane; γ k is a vector The angle between the X axis and the

[0021] The original point cloud is converted into a scanning pattern grid with a starting interval of d; β k and γ k Calculated using the coordinates of adjacent path points:

[0022]

[0023] Furthermore, in step S20, the process of obtaining a complete road point cloud includes the following steps:

[0024] S21: Use the height histogram method to filter out the objects above the road surface. When reconstructing the point cloud, divide the elevation value range of all points into equal parts at intervals of 0.2 m, obtain the number of points in each range, the index, and the histogram with the ordinate as the number of points, find the elevation value range corresponding to the highest segment, and retain the points with elevation values ​​near this elevation value range;

[0025] S22: Filter out the facilities on both sides of the road by geometric range demarcation, reconstruct the point cloud, and remove the outliers on the roadside of the area by demarcating the area of ​​interest to obtain a preliminary demarcated road point cloud;

[0026] S23: For the initially delineated road point cloud, establish a K nearest neighbor point cloud index and calculate the average distance of all point clouds; traverse all point clouds, find the K nearest neighbor points near each query point, and calculate the average distance between the query point and its K nearest neighbors; if the average distance of the neighbor points of any query point is greater than the average distance of all point clouds, mark the query point as an outlier; after all points are traversed, remove all outliers in the road point cloud;

[0027] S24: For the road point cloud from which outliers have been removed, the point cloud whose Euclidean distance of all points is less than a given threshold is grouped into a cluster. The threshold is empirically determined as twice the average distance of all points. After clustering, 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, which is segmented according to the label.

[0028] Furthermore, in step S30, the multi-step adaptive height threshold object filtering process includes the following steps:

[0029] S31, establishing a virtual grid, dividing the plane area where the discrete three-dimensional point cloud is located into a plurality of virtual grids of the same size using a grid, each virtual grid is equivalent to a subspace container of the point cloud space, and each laser point will fall into one of the grids;

[0030] S32, traverse all grids and find the lowest value z′ of the elevation in the grid min , calculate the average value z′ of the point cloud elevation in the grid mean , calculate the elevation fluctuation parameter δ z′ ; Traverse the point cloud in the grid and remove the points with elevation values ​​greater than z′ min +δ z′ The laser point;

[0031] S33, traverse all grids and find the lowest value z′ of the elevation of all laser points in the 8 neighboring grids around the grid min2 , calculate the height fluctuation parameter ω at this time z′ ; Traverse the point cloud in the grid and remove the points with elevation values ​​greater than z′ min2 +ω z′ The laser is used to calculate the standard deviation of the elevation values ​​at each point;

[0032] S34, repeat step S33 until the standard deviation change between two adjacent times is less than the set critical value, and stop filtering.

[0033] Furthermore, in step S31, a virtual grid is established, and the process of dividing the plane area where the discrete three-dimensional point cloud is located into a plurality of virtual grids of the same size by using the grid includes the following steps:

[0034] S311, establish a virtual grid in the Matlab environment, determine the virtual grid size ε, create a [(y′ l -y′ r ) / ε+1]×[x′ max / ε+1];

[0035] S312, traverse all the original laser point cloud data, solve its coordinate range in the XOZ plane, and put it into the corresponding cell array; any point (x′ i , y′ i , z′ i ) is located in the following virtual grid rows and columns:

[0036] w=[x′ i -x′ min / ε]+1

[0037] l=[(y′ i -y′ min ) / ε]+1

[0038] In the formula, w and l represent the row and column numbers of the virtual grid where the point is located, and x′ min , y′ min represents the minimum coordinate in the point set, [.] represents rounding, and ε is the virtual grid size.

[0039] Furthermore, for any point (x′ i , y′ i , z′ i ), parameters in the adaptive height threshold method and The calculation method is as follows:

[0040]

[0041]

[0042] Where ε is the virtual grid size, is the height threshold of the initial filtering, is the height threshold of the secondary filtering, is the influence range of the secondary filter.

[0043] Furthermore, in step S40, the process of filling the point cloud holes using the virtual grid-based interpolation algorithm includes the following steps:

[0044] S41, extracting road boundary point cloud;

[0045] S42, using a robust weighted local weighted regression algorithm to smooth the road boundary and obtain a complete road boundary point cloud;

[0046] S43, traversing all virtual grids, if a grid satisfies the condition that it is within the road boundary and has no points inside, then the area is regarded as a hole area; searching for the hole area within the road range, and uniformly generating the plane two-dimensional coordinates of the points to be interpolated;

[0047] S44, calculate the vertical coordinates of the point to be interpolated, find 8 non-empty grids around the grid where the point to be interpolated is located, use Delaunay triangulation space interpolation to determine the vertical coordinates of the point to be interpolated, and when the vertical coordinates of the interpolation points in all the hole areas are calculated, the hole filling is completed.

[0048] Furthermore, in step S42, a robust weighted local weighted regression algorithm is used to smooth the road boundary, and the process of obtaining a complete road boundary point cloud includes the following steps:

[0049] S421, defining the window width of the filter, where the window width represents the ratio of the number of data points used to calculate the smoothing value to the total number of all data points;

[0050] S422, traverse all points, find all points within the window width of a given point, and calculate the weight values ​​of all neighboring points of the point. The calculation method is as follows:

[0051]

[0052] In the formula, x represents the point to be smoothed, x i represents the nearby point within the span of x, and dis is the horizontal distance from x to the farthest predicted value within the window width;

[0053] S423, obtaining a temporary smoothing value x of x t , using the calculated weight w i Perform weighted linear regression on x, let X t is a temporary smoothed point set;

[0054] S424, calculate the point set X t The point-by-point differences R{r1, r2, r3...r k ...r n |1≤k<n,k,n∈N +}, based on the point-by-point difference, the weight values ​​of all adjacent points within the window width are calculated again. The calculation method is as follows:

[0055]

[0056] Where: r i is the residual of the ith point, δ is the median absolute deviation of the residual; when the residual is greater than 6δ and the robustness weight is 0, the related outliers will be eliminated during the calculation process;

[0057] S425, using w i and Smooth the road boundaries.

[0058] Further, in step S50, the process of obtaining the cross slope value includes the following steps:

[0059] S51, along the direction of the road, consistent with the Y-axis direction in the reconstruction space, extract a section of a certain thickness at a given interval as a cross section;

[0060] S52, for each extracted cross section, the elevation and lateral values ​​of the cross section point set are fitted using a random sampling consensus algorithm. The slope of the model is the cross slope value, which is expressed as a superelevation value in the curve segment;

[0061] S53, referring to the current urban road design specifications, compare the specification values ​​with the calculated values, and mark the sites where the cross slope values ​​do not meet the specifications.

[0062] Beneficial effects:

[0063] The method for calculating the cross slope and superelevation of urban roads based on vehicle-mounted laser point cloud data mentioned in the present invention can accurately calculate the cross slope from dense MLS data in the presence of considerable vehicle noise and data gaps. This method can be used to evaluate the cross-sectional conditions of built roads. It is generally applicable to pavement maintenance and road reconstruction. By determining the location of substandard slopes, relevant road agencies can make repairs in a timely manner to avoid accidents caused by water damage to the pavement. Therefore, the research results can also provide a basis for road widening when considering road settlement or lack of design data. The method can accurately obtain road cross-sectional information. With the advancement of building information modeling and lidar technology, this road inventory can be imported into the design and modeling platform for full life cycle design and maintenance. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] Figure 1 The present invention is a flowchart of a method for calculating the cross slope and superelevation value of an urban road based on vehicle-mounted laser point cloud data according to an embodiment of the present invention.

[0065] Figure 2 It is a schematic diagram of the effect of reconstructing the original point cloud into a scanned format grid according to an embodiment of the present invention.

[0066] Figure 3 This is a schematic diagram of the effect of extracting roads by using a height histogram and geometric division according to an embodiment of the present invention; wherein: Figure 3 (a) is a schematic diagram of the reconstructed spatial point cloud; Figure 3 (b) is a schematic diagram of the road point cloud after division.

[0067] Figure 4 It is a schematic diagram of the final effect of the road extraction method according to an embodiment of the present invention.

[0068] Figure 5 Schematic diagram of the structure of Delaunay triangulated space interpolation according to an embodiment of the present invention.

[0069] Figure 6 This is a schematic diagram of the road hole filling effect of an embodiment of the present invention; wherein: Figure 6 (a) is a schematic diagram of the road surface before the hole is filled; Figure 6 (b) is a schematic diagram of the road surface after the void is filled.

[0070] Figure 7 This is a schematic diagram of the effect of RANSAC fitting a road cross section according to an embodiment of the present invention. DETAILED DESCRIPTION

[0071] The following examples will enable those skilled in the art to more fully understand the present invention, but are not intended to limit the present invention in any way.

[0072] Figure 1This is a structural schematic diagram of a method for calculating the cross slope and superelevation value of urban roads based on vehicle-mounted laser point cloud data. This embodiment is applicable to a method for calculating the cross slope and superelevation value of urban roads based on vehicle-mounted laser point cloud data. Figure 1 As shown, the calculation method includes the following steps:

[0073] S10, preprocessing the point cloud data to reconstruct the original point cloud data into an aligned scan format grid.

[0074] S20,adopts the height histogram-based, k-nearest-neighbor-based filtering and Euclidean distance-based clustering technology to segment the road point cloud,separates most of the non-ground points (surrounding buildings, street lights, signs, etc.) and obtains a complete road point cloud.

[0075] S30, establishing a point cloud spatial index, based on which a two-step adaptive height threshold method is used to remove ground object points (vehicles on the road, vegetation).

[0076] S40, uses a virtual grid-based interpolation algorithm to fill the point cloud holes to obtain a complete road point cloud.

[0077] S50, extracting road cross sections at a given interval, and for each extracted cross section, regressing the elevation and lateral values ​​using a random consistency sampling algorithm to obtain a cross slope value. Compare the calculated cross slope value with the specification to determine the location of the cross slope and superelevation value that do not meet the specification.

[0078] Figure 2 The schematic diagram of the effect of reconstructing the original point cloud into a scanned format grid is shown in the present invention. In one embodiment, in step S10, the point cloud data preprocessing process includes the following steps:

[0079] S11, coordinate conversion is performed on the original laser point cloud data, and the point cloud is reconstructed into a scanning format grid according to the corresponding relationship between the original laser point cloud data and the trajectory data.

[0080] First, use the timestamp of each laser point to divide these points into scan lines. Each scan line has a corresponding trajectory point. Let T{T1, T2...T k ...T n |1≤k<n,k,n∈N +} is a set of trajectory points. Let S{S1, S2...S k ...S n |1≤k<n,k,n∈N +} is the LiDAR point set of the scan line corresponding to the trajectory point. For the kth scan line, let T k is the starting point of the trajectory vector, T k+1 is the end point of the trajectory vector. Therefore, the forward vector Expressed as T k As the origin, The local three-dimensional coordinate system is established with the X' axis as the X' axis. The Y' axis is orthogonal to the X' axis in the horizontal direction. Then the Z' axis is the upward vector perpendicular to the X'-Y' plane. Using the coordinate transformation matrix, the geodesic point is transformed into a local point in the X'Y'Z' space:

[0081]

[0082] The meanings of the parameters in the formula are as follows:

[0083] Track point coordinates.

[0084] (x k y k z k ) T =S k The coordinates of a point described in the geodetic coordinate system.

[0085] (x′ k y′ k z′ k ) T =S k The coordinates of a point described in a local coordinate system.

[0086] β k = the angle between the X-axis and the geodetic plane.

[0087] γ k = Vector The angle between the X axis and the

[0088] The original point cloud is converted into a scan pattern grid with a starting interval of d. k and γ k is the key parameter in the calculation, which can be calculated using the coordinates of adjacent path points:

[0089]

[0090] S12, using a robust local weighted regression algorithm to smooth important parameters in the point cloud coordinate transformation to ensure the reconstruction effect. The curve segments of the reconstructed road point cloud will be straightened, and all road point clouds will be basically in a similar elevation range.

[0091] Figure 3 This is a schematic diagram of the effect of extracting roads by using a height histogram and geometric division in the present invention; wherein: Figure 3 (a) is a schematic diagram of the reconstructed spatial point cloud; Figure 3 (b) is a schematic diagram of the road point cloud after division. Figure 4 Schematic diagram of the final effect of the road extraction method of the present invention.

[0092] In one embodiment, in step S20, the process of performing road surface recognition and segmentation processing on the point cloud data includes the following steps:

[0093] S21: Use the height histogram method to filter out objects above the road surface: When reconstructing the point cloud, divide the elevation value range of all points into equal parts at intervals of 0.2m. You can get the number and index of points in each range, and a histogram with the ordinate as the number of points. Find the elevation value range corresponding to the highest segment, and retain points with elevation values ​​near this range.

[0094] S22: Filter out most of the facilities on both sides of the road by defining the geometric range: After reconstructing the point cloud, the lateral range of the road point cloud (i.e., the Y value) does not fluctuate much. Therefore, the outliers on the road side can be removed by defining the area of ​​interest.

[0095] S23: Perform K-nearest-neighbor based filtering: For the initially delineated road point cloud, establish a K-nearest-neighbor point cloud index, calculate the average distance of all point clouds, then traverse all point clouds, find the K nearest neighbor points nearby, and calculate the average distance between the query point and its K nearest neighbors. If the average distance between the neighbor points of the point is greater than the average distance of all point clouds, mark the point as an outlier. After all points are traversed, all outliers are removed.

[0096] S24: Perform ground point cloud segmentation based on Euclidean distance clustering: for the road point cloud after removing outliers, group all points whose Euclidean distance is less than a given threshold into a cluster. The threshold is empirically determined to be twice the average distance of all points. After clustering, all points in the point cloud have a cluster label. The cluster with the largest number of points is the road point cloud, which can be segmented according to the label.

[0097] Figure 4 The final effect diagram of the road extraction of the present invention is shown in FIG. 1 . In step S24 of this embodiment, the processing process of the Euclidean clustering method includes the following steps:

[0098] S241: Find a point p in space 11 , use KD-Tree to find the n points closest to it, and determine the distance from these n points to p 11 The Euclidean distance of the point p whose distance is less than the threshold r 12 , p 13 , p 14 ...and put it in class Q.

[0099] S242: In Q(p 11 ) to find a little p 12 , repeat step S241.

[0100] S243: In Q(p11 , p 12 ) find a point, repeat step S241, find p 22 , p 23 , p 24 ...put it all in Q.

[0101] S244: When no new points can be added to Q, the search is completed.

[0102] In one embodiment, in step S30, the process of filtering out objects using the adaptive height threshold comprises the following steps:

[0103] S31: Establish a virtual grid: Use a grid to divide the plane area where the discrete three-dimensional point cloud is located into a number of virtual grids of the same size. Each virtual grid is equivalent to a subspace container of the point cloud space, and each laser point will fall into one of the grids.

[0104] Specifically, to establish a virtual grid in the Matlab environment, first determine the virtual grid size ε, create a [(y′ l -y′ r ) / ε+1]×[x′ max / ε+1], and then traverse all the original laser point cloud data, solve its coordinate range in the XOZ plane, and put it into the corresponding cell array. i , y′ i , z′ i ) is located in the following virtual grid rows and columns:

[0105] w=[x′ i -x′ min / ε]+1

[0106] l=[(y′ i -y′ min ) / ε]+1

[0107] The meanings of the parameters in the formula are as follows: w, l represent the row and column numbers of the virtual grid where the point is located, x′ min , y′ min Indicates the minimum coordinate in the point set, and [.] indicates rounding.

[0108] S32: Perform the first step of height threshold feature filtering: traverse all grids and find the lowest value z′ of the elevation in the grid min , calculate the average value z′ of the point cloud elevation in the grid mean , calculate the elevation fluctuation parameter δ z′ Traverse the point cloud in the grid and remove the points with elevation values ​​greater than z′ min +δ z′ of laser dots.

[0109] Perform the second step of height threshold feature filtering: traverse all grids and find the lowest value z′ of the elevation of all laser points in the 8 neighboring grids around the grid min , calculate the height fluctuation parameter ω at this time z′ Traverse the point cloud in the grid and remove the points with elevation values ​​greater than z′ min +ω z′ The standard deviation of the elevation values ​​of each point is calculated and compared with the last execution result. When the standard deviation changes below the critical value (the critical value is set to 0.2 in this study), the filtering is stopped.

[0110] The design of this algorithm has three considerations: (1) The road point cloud is flat and its elevation fluctuates very little within a certain range; in addition, the road points are continuous along the cross section, and their elevations do not increase or decrease significantly except for the intermediate obstacles and roadside obstacles. (2) The road line shape is converted into a straight line in S20, eliminating the adverse effects of the curve segment in the mesh generation process. For any point (x′ i , y′ i , z′ i ), parameters in the adaptive height threshold method and The calculation method is as follows:

[0111]

[0112]

[0113] In one embodiment, S40, the filling of point cloud holes by the virtual grid-based interpolation algorithm comprises the following steps:

[0114] S41, extracting road boundary point cloud.

[0115] Holes are mainly found and filled in the area within the road boundary. Holes outside the road are not within the scope of subsequent hole filling and cross-section calculation. Therefore, it is necessary to extract the road boundary. The road surface in the reconstructed scene is divided into strip units with a segmentation size of ∈. For each strip unit, search for its leftmost point and rightmost point respectively. Considering that the calculations in each strip unit are independent and similar at this stage, parallel computing can be used to improve the actual operation efficiency.

[0116] S42, a robust local weighted linear regression algorithm (RLWLR) is used to smooth the road boundary to obtain a complete road boundary point cloud.

[0117] After processing in step S41, the left and right boundaries of the road are basically extracted. However, due to the presence of some incompletely filled holes on the road surface, there are deviation points in the extracted road boundary point cloud. Therefore, it is necessary to use a smoothing algorithm to obtain boundary curve points that are more closely aligned with the road surface boundary.

[0118] In step S42 of this embodiment, the processing process of the robust local weighted linear regression algorithm includes the following steps:

[0119] S421, defining the window width of the filter, that is, the ratio of the number of data points used to calculate the smoothing value to the total number of all data points.

[0120] S422, traverse all points, find all points within the window width of a given point, and calculate the weight values ​​of all neighboring points of the point. The calculation method is as follows:

[0121]

[0122] Where: x represents the point to be smoothed, x i Represents the adjacent point within the span of x, and dis is the horizontal distance from x to the farthest predicted value within the window width.

[0123] S423, obtaining a temporary smoothing value x of x t , the weight w calculated using the above formula i Perform weighted linear regression on x, let X t is a temporary smoothed point set.

[0124] S424, calculate the point set X t The point-by-point differences (i.e., residuals) R{r1, r2, r3...r k ...r n |1≤k<n,k,n∈N +}, based on this difference, the weight values ​​of all adjacent points within the window width are calculated again, and the calculation method is as follows:

[0125]

[0126] Where: r i is the residual of the ith point, and δ is the median absolute deviation of the residual.

[0127] When the residual is larger than 6δ, the robustness weight is 0 and the related outliers will be eliminated during the calculation process.

[0128] S424, smoothed data using w i and

[0129] To ensure the smoothing effect, step S422 and step S423 need to be repeated multiple times. According to experience, the window width is 20.

[0130] S43, searching for a hole area within the road range, and uniformly generating the plane two-dimensional coordinates of the points to be interpolated.

[0131] Traverse all virtual grids, if a grid satisfies the road boundary and has no points inside, then the area is considered as a hole area. After the hole area is determined, two-dimensional points are uniformly generated in the virtual grid, which are the plane two-dimensional coordinates of the points to be interpolated.

[0132] S44, calculate the vertical coordinates of the point to be interpolated, find the points in the 8 non-empty grids around the grid where the point to be interpolated is located, use the points in the point set as reference points, use Delaunay triangulation space interpolation to determine the vertical coordinates of the point to be interpolated, and when the vertical coordinates of the interpolation points in all the hole areas are calculated, the hole filling is completed.

[0133] Figure 5 A schematic diagram of the structure of the Delaunay triangulated space interpolation of the present invention. Figure 6 Schematic diagram of the road hole filling effect of the present invention. Figure 6 (a) is a schematic diagram of the road surface before the hole is filled; Figure 6 (b) is a schematic diagram of the road surface after the hole is filled. The present invention uses a linear interpolation method to complete the hole filling work. The size of the hole area is related to the size of the vehicle on the road and the angle of the laser beam. According to the empirical evaluation of different laser point cloud data sets, the size of the hole is usually approximately a rectangular area with a length of 2 to 5 meters and a width of 7 to 15 meters. Therefore, the range of the reference point should be larger than the rectangular area to ensure that the reference point can also be found at the center of the hole. Figure 6 As shown, 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 ] are three reference points. The Delaunay triangular space interpolation algorithm is used to perform linear interpolation fitting on the interpolation points to obtain the elevation values ​​of the interpolation points. The calculation process is as follows:

[0134]

[0135] Where: z′ j The elevation value of the point to be interpolated.

[0136] In one embodiment, S50, the cross section extraction and cross slope value calculation includes the following steps:

[0137] S51, extracting cross sections at given intervals, along the direction of the road, consistent with the Y-axis direction in the reconstruction space, and extracting sections of a certain thickness at given intervals, which are cross sections;

[0138] In the point cloud preprocessing step S20, the present invention converts the complex highway line shape into a straight line by means of point cloud reconstruction. At this time, the road advance direction is the same as the x' axis direction, and obtaining the section perpendicular to the road advance direction is simpler than the known method. Theoretically, all horizontal coordinates are x' i The plane formed by the points is the cross section at that location. However, considering that the distribution of points in the point cloud is not continuous, we can take a certain width of area [x′ i -d, x′ i +d] to extract the cross section. Considering that when conducting on-site cross-section measurement, it is generally measured at intervals of 10-30m. Therefore, this paper extracts the cross section at intervals of 10m, and the width of the cross section extraction is 0.4m. Since the width value of the cross section extraction will affect the data interval range for cross slope calculation, the selection of this value may affect the accuracy of cross slope calculation. In order to determine the optimal value of this value, the present invention has conducted sensitivity analysis on multiple measured sections.

[0139] S52, calculate the cross slope value and superelevation value. For each extracted cross section, use the random sampling consistency algorithm to fit the elevation and lateral values ​​of the cross section point set. The slope of the model is the cross slope value (superelevation value in the curve section).

[0140] According to the definition of cross slope, the cross slope value can be obtained by linear regression of the elevation and lateral values ​​of the laser points in each cross-section point set. Since it is impossible to completely segment all the ground point clouds when extracting the road surface, there are still a small number of noise points in the extracted cross section. The least squares method (LS) is a simple fitting algorithm, and its principle is to solve the model based on minimizing the mean square error. In linear regression, the least squares method tries to find a straight line that minimizes the sum of the Euclidean distances of all sample points to the straight line. However, this algorithm will try to adapt to all points as much as possible, resulting in low fitting accuracy when there are a large number of noise points.

[0141] The present invention uses the Random Sample Consensus (RANSAC) algorithm to calculate the cross slope value corresponding to each cross section. The random sampling consensus algorithm uses an iterative method to solve the model from a set of data containing outliers. The idea of ​​the algorithm is as follows: 1) Randomly select SS points as sample points and fit the model; 2) Find the points within the tolerance range MD of the fitting line and count the number of points; 3) Randomly select MD points again, repeat steps 2) and 3) until the iteration ends; 4) Find the case with the most data points after a certain fitting, which is the solved model.

[0142] 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 the empirical evaluation of two different laser point cloud datasets, any point that is more than 0.2m away from the fitted model is considered an outlier (MD = 0.2m, the distance metric is the square of the Euclidean distance, and the maximum number of iterations is 1000). Figure 7 Four examples of using RANSAC to fit a straight line are given. Figure 7 It can be seen that the fitted straight line is very close to the point set, and the RANSAC method has a better fitting effect on the points than the LS method.

[0143] S53, safety evaluation of cross slope value, refers to the current urban road design specifications, compares the specification value with the calculated value, and marks the sites where the cross slope value does not meet the specifications.

[0144] The above are only preferred embodiments of the present invention, and the protection scope of the present invention is not limited to the above embodiments. All technical solutions under the concept of the present invention belong to the protection scope of the present invention. It should be pointed out that for ordinary technicians in this technical field, some improvements and modifications without departing from the principle of the present invention should be regarded as the protection scope of the present invention.

Claims

1. A method for calculating the cross slope and superelevation value of urban roads based on vehicle-mounted laser point cloud data, characterized in that: The road cross slope and superelevation value calculation method comprises the following steps: S10, preprocessing the point cloud data to reconstruct the original point cloud data into an aligned scan format grid; S20, uses height histogram-based, k-nearest-neighbor-based filtering and Euclidean distance-based clustering technology to segment the road point cloud, separates non-ground points including surrounding buildings, street lights and signs, and obtains a complete road point cloud; S30, establishing a point cloud spatial index, based on which a two-step adaptive height threshold method is used to remove ground object points including vehicles on the road and vegetation; S40, uses a virtual grid-based interpolation algorithm to fill the point cloud holes to obtain a complete road point cloud; S50, extracting road cross sections at preset intervals, and for each extracted cross section, regressing the elevation and lateral values ​​using a random consistency sampling algorithm to obtain a cross slope value; comparing the calculated cross slope value with a standard cross slope value to determine the location of a cross slope that does not meet the standard and a superelevation value; In step S30, the process of removing ground object points including vehicles and vegetation on the road by using the two-step adaptive height threshold method includes the following steps: S31, establishing a virtual grid, dividing the plane area where the discrete three-dimensional point cloud is located into a plurality of virtual grids of the same size using a grid, each virtual grid is equivalent to a subspace container of the point cloud space, and each laser point will fall into one of the grids; S32, traverse all grids and find the lowest value z′ of the elevation in the grid min , calculate the average value z′ of the point cloud elevation in the grid mean , calculate the elevation fluctuation parameter δ z′ ; Traverse the point cloud in the grid and remove the points with elevation values ​​greater than z′ min +δ z′ The laser point; S33, traverse all grids and find the lowest value z′ of the elevation of all laser points in the 8 neighboring grids around the grid min2 , calculate the height fluctuation parameter ω at this time z′ ; Traverse the point cloud in the grid and remove the points with elevation values ​​greater than z′ min2 +ω z′ The laser is used to calculate the standard deviation of the elevation values ​​at each point; S34, repeat step S33 until the standard deviation change between two adjacent times is less than the set critical value, and stop filtering.

2. The method for calculating the cross slope and superelevation value of urban roads based on vehicle-mounted laser point cloud data according to claim 1 is characterized in that: In step S10, the point cloud data is preprocessed to reconstruct the original point cloud data into an aligned scan format grid, including the following steps: S11, coordinate conversion is performed on the original laser point cloud data, and the point cloud is reconstructed into a scanning format grid according to the corresponding relationship between the original laser point cloud data and the trajectory data; S12, using a robust local weighted regression algorithm to smooth multiple parameters in the point cloud coordinate transformation, straightening the curve segments in the road point cloud, and all road point clouds are in a similar elevation range.

3. The method for calculating the cross slope and superelevation value of urban roads based on vehicle-mounted laser point cloud data according to claim 2 is characterized in that: In step S11, the process of reconstructing the point cloud into a scan format grid includes the following steps: S111, using the timestamp of each laser point, divide all laser points into scan lines, each scan line has a corresponding trajectory point; let T{T1,T2…T k …T n |1≤k <n,k,n∈N + } is a set of trajectory points, S{S1,S2…S k …S n |1≤k <n,k,n∈N + } is the LiDAR point set of the scan line corresponding to the trajectory point; S112, for the kth scan line, let T k is the starting point of the trajectory vector, T k+1 is the end point of the trajectory vector, the positive vector Expressed as T k As the origin, For X ′ Axis, establish a local three-dimensional coordinate system, the Y' axis is orthogonal to the X' axis in the horizontal direction, and the Z' axis is perpendicular to the upward vector of the X'-Y' plane; use the coordinate transformation matrix to transform the geodesic point into a local point in the X'Y'Z' space: Where: is the trajectory point T k The coordinates of (x k y k z k ) T S is described by the geodetic coordinate system k The coordinates of the point; (x ′ k y ′ k z ′ k ) T is described by the local coordinate system S k The point coordinates of β k is the angle between the X-axis and the geodetic plane; γ k is a vector The angle between the X axis and the The original point cloud is converted into a scanning pattern grid with a starting interval of d; β k and γ k Calculated using the coordinates of adjacent path points:

4. The method for calculating the cross slope and superelevation value of urban roads based on vehicle-mounted laser point cloud data according to claim 1 is characterized in that: In step S20, the process of obtaining a complete road point cloud includes the following steps: S21: Use the height histogram method to filter out the objects above the road surface. When reconstructing the point cloud, divide the elevation value range of all points into equal parts at intervals of 0.2 m, obtain the number of points in each range, the index, and the histogram with the ordinate as the number of points, find the elevation value range corresponding to the highest segment, and retain the points with elevation values ​​near this elevation value range; S22: Filter out the facilities on both sides of the road by geometric range demarcation, reconstruct the point cloud, and remove the outliers on the roadside of the area by demarcating the area of ​​interest to obtain a preliminary demarcated road point cloud; S23: For the initially delineated road point cloud, establish a K nearest neighbor point cloud index and calculate the average distance of all point clouds; traverse all point clouds, find the K nearest neighbor points near each query point, and calculate the average distance between the query point and its K nearest neighbors; if the average distance of the neighbor points of any query point is greater than the average distance of all point clouds, mark the query point as an outlier; after all points are traversed, remove all outliers in the road point cloud; S24: For the road point cloud from which outliers have been removed, the point cloud whose Euclidean distance of all points is less than a given threshold is grouped into a cluster. The threshold is empirically determined as twice the average distance of all points. After clustering, 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, which is segmented according to the label.

5. The method for calculating the cross slope and superelevation value of urban roads based on vehicle-mounted laser point cloud data according to claim 1 is characterized in that: In step S31, a virtual grid is established, and the process of dividing the plane area where the discrete three-dimensional point cloud is located into a plurality of virtual grids of the same size by using the grid includes the following steps: S311, establish a virtual grid in the Matlab environment, determine the virtual grid size ε, and create a [(y′ l -y′ r ) / ε+1]×[x′ max / ε+1]; S312, traverse all the original laser point cloud data, solve its coordinate range in the XOZ plane, and put it into the corresponding cell array; any point (x′ i ,y′ i ,z′ i ) is located in the following virtual grid rows and columns: w=[x′ i -x′ min / ε]+1 l=[(and′ i -and' min ) / ε]+1 In the formula, w and l represent the row and column numbers of the virtual grid where the point is located, and x′ min , y′ min represents the minimum coordinate in the point set, [.] represents rounding, and ε is the virtual grid size.

6. The method for calculating the cross slope and superelevation value of urban roads based on vehicle-mounted laser point cloud data according to claim 1, characterized in that: For any point (x′ i ,y′ i ,z′ i ), parameters in the adaptive height threshold method and The calculation method is as follows: Where ε is the virtual grid size, is the height threshold of the initial filtering, is the height threshold of the secondary filtering, is the influence range of the secondary filter.

7. The method for calculating the cross slope and superelevation value of urban roads based on vehicle-mounted laser point cloud data according to claim 1 is characterized in that: In step S40, the process of filling the point cloud holes using the interpolation algorithm based on the virtual grid includes the following steps: S41, extracting road boundary point cloud; S42, using a robust weighted local weighted regression algorithm to smooth the road boundary and obtain a complete road boundary point cloud; S43, traversing all virtual grids, if a grid satisfies the condition that it is within the road boundary and has no points inside, then the area is regarded as a hole area; searching for the hole area within the road range, and uniformly generating the plane two-dimensional coordinates of the points to be interpolated; S44, calculate the vertical coordinates of the point to be interpolated, find 8 non-empty grids around the grid where the point to be interpolated is located, use Delaunay triangulation space interpolation to determine the vertical coordinates of the point to be interpolated, and when the vertical coordinates of the interpolation points in all the hole areas are calculated, the hole filling is completed.

8. The method for calculating the cross slope and superelevation value of urban roads based on vehicle-mounted laser point cloud data according to claim 7 is characterized in that: In step S42, a robust weighted local weighted regression algorithm is used to smooth the road boundary, and the process of obtaining a complete road boundary point cloud includes the following steps: S421, defining the window width of the filter, where the window width represents the ratio of the number of data points used to calculate the smoothing value to the total number of all data points; S422, traverse all points, find all points within the window width of a given point, and calculate the weight values ​​of all neighboring points of the point. The calculation method is as follows: In the formula, x represents the point to be smoothed, x i represents the nearby point within the span of x, and dis is the horizontal distance from x to the farthest predicted value within the window width; S423, obtaining a temporary smoothing value x of x t , using the calculated weight w i Perform weighted linear regression on x, let X t is a temporary smoothed point set; S424, calculate the point set X t The point-by-point differences R{r1,r2,r3…r k …r n |1≤k <n,k,n∈N + }, based on the point-by-point difference, the weight values ​​of all adjacent points within the window width are calculated again. The calculation method is as follows: Where: r i is the residual of the ith point, δ is the median absolute deviation of the residual; when the residual is greater than 6δ and the robustness weight is 0, the related outliers will be eliminated during the calculation process; S425, using w i and Smooth the road boundaries.

9. The method for calculating the cross slope and superelevation value of urban roads based on vehicle-mounted laser point cloud data according to claim 1, characterized in that: In step S50, the process of obtaining the cross slope value includes the following steps: S51, along the direction of the road, consistent with the Y-axis direction in the reconstruction space, extract a section of a certain thickness at a given interval as a cross section; S52, for each extracted cross section, the elevation and lateral values ​​of the cross section point set are fitted using a random sampling consensus algorithm. The slope of the model is the cross slope value, which is expressed as a superelevation value in the curve segment; S53, referring to the current urban road design specifications, compare the specification values ​​with the calculated values, and mark the sites where the cross slope values ​​do not meet the specifications.

Images