A visual evaluation method for urban road intersections based on point cloud data

Abstract

The present invention discloses a visualization evaluation method for urban road intersections based on point cloud data. The method collects point cloud data of road intersections, performs road surface recognition and segmentation processing on the point cloud data, separates the road surface part from other discrete plane points to obtain complete road surface point cloud data, establishes a plane grid dot matrix, uses the recognized road surface points as benchmark interpolation references, uses a neighboring point interpolation algorithm to obtain height information of the plane grid points, constructs a digital elevation model, constructs a coordinate system with the starting point of the driver's line of sight as the origin, converts the three-dimensional coordinate system into a two-dimensional coordinate system, and then converts the two-dimensional rectangular coordinate system into a polar coordinate system to determine the position of obstacles blocking the driver's line of sight. The present invention uses a cylindrical perspective projection model to realize the conversion of a three-dimensional problem into a two-dimensional problem, constructs a visual domain model, analyzes and determines the position of obstacles that affect the driver's visibility, and analyzes the visualization of the intersection.

Description

Technical Field

[0001] The present invention relates to the technical field of road three-dimensional sight distance detection, and in particular to a method for visualizing and evaluating urban road intersections based on point cloud data. Background Art

[0002] Urban road intersections are frequent sections of traffic accidents. In the road design stage, the sight triangle of the intersection is used as an important visibility evaluation index to ensure good intersection visibility. However, compared with the roads in the design stage, the three-dimensional and real-time road intersection environment in reality is more complex. There are many types of potential sight obstacles that are not considered or insufficiently considered in the design stage, such as vegetation growth, new facilities, and debris accumulation. The current means of obtaining road infrastructure data, such as surveying-level laser radar, drones, and communication satellites, have the advantages of high data accuracy and large data volume, and have obvious advantages for macro-scale traffic road data acquisition. However, it has certain limitations for micro-scale traffic infrastructure such as intersections, mainly affecting traffic travel and high equipment loss costs. Small laser radars are low-cost and easy to operate, but their application in engineering practice is also limited due to their small viewing angle of 38.4°. Three-dimensional sight distance calculation is the basis for visual evaluation of urban road intersections. Domestic road three-dimensional sight distance calculation methods have fundamentally simplified the three-dimensional environment of the road to varying degrees, and their application is greatly limited in road intersections with complex traffic flow composition and many roadside obstacles. Due to modeling limitations, many foreign road three-dimensional sight distance methods can only be called 2.5D, and cannot handle objects in the same position but at different heights. A few calculation methods based on 3D models also have the limitations of complex algorithms and long calculation time, and there is room for further improvement. The existing point cloud data acquisition methods and intersection visualization evaluation programs still have major defects in economy and efficiency, and the evaluation consideration of actual driver conditions and road environment conditions is not comprehensive enough. Summary of the invention

[0003] In view of the deficiencies in the prior art, the present invention provides a method for visualizing and evaluating urban road intersections based on point cloud data, which takes into account different parameters such as the driver's line of sight height, range, and viewing angle position.

[0004] To achieve the above object, the present invention adopts the following technical solutions:

[0005] In a first aspect, an embodiment of the present invention provides a method for visual evaluation of urban road intersections based on point cloud data, wherein the method for calculating three-dimensional sight distance of a road intersection comprises the following steps:

[0006] S1, collects point cloud data of road intersections;

[0007] S2, performing road surface recognition and segmentation processing on the point cloud data, distinguishing the road surface part from other discrete plane points to obtain complete road surface point cloud data;

[0008] S3, establishing a plane grid dot matrix, using the identified road surface points as the base interpolation reference, using the neighboring point interpolation algorithm to obtain the height information of the plane grid points, and constructing a digital elevation model;

[0009] S4, quantifies the concept of visual field into a fan-shaped shape composed of countless rays emitted from the origin of the driver's line of sight. The angle of the fan-shaped shape is related to the driver's visual range. Through the cylindrical perspective projection model, the three-dimensional coordinate points in space are converted into two-dimensional coordinate points, and the three-dimensional problem is converted into a two-dimensional problem.

[0010] A coordinate system is constructed with the starting point of the driver's line of sight as the origin, the three-dimensional coordinate system is converted into a two-dimensional coordinate system, and then the two-dimensional rectangular coordinate system is converted into a polar coordinate system to determine the position of the obstacle blocking the driver's line of sight;

[0011] S5, analyze the variable factors that may affect the visibility of urban road intersections.

[0012] Furthermore, in step S1, the process of collecting point cloud data of the road intersection includes the following steps:

[0013] S11, fix the small laser radar on the rotating gimbal, place it at the center of the urban road intersection, and keep the laser radar mounting bracket fixed; collect data multiple times through the rotating gimbal to obtain the sub-file of the overall point cloud data of the intersection, and the overlap degree of two adjacent point cloud files is 1 / 3;

[0014] S12, performing data splicing on the obtained point cloud files through corresponding feature points to obtain complete point cloud data of urban road intersections.

[0015] Furthermore, in step S2, the process of performing road surface recognition and segmentation processing on the point cloud data includes the following steps:

[0016] S21: Rasterize the original 3D point cloud data. Each grid point of the rasterized point cloud data has a corresponding linear index value. The data points at the same position in the same matrix will have the same linear index value. When a data point crosses a matrix, the linear index value ζ of the data point will have a step.

[0017] S22: Using the height threshold h t The points with lower elevation in the columnar unit are screened out, and the remaining data points are analyzed using the principal component analysis method to obtain the main direction of the columnar unit;

[0018] S23: Use the height histogram method to sort the detected road surface points and extract features from the obtained plane point cloud data; the feature variables include the number of point clouds per unit area, and the standard deviation and median of the difference between adjacent items along the x and y directions used to describe the distribution characteristics of the point cloud on the xy plane, and obtain an observation value array consistent with the number of cylindrical units; use the unsupervised classification K-Means clustering algorithm to divide the cylindrical units into three categories: uniform plane points, non-uniform plane points and pseudo-plane points; standardize all variables to eliminate their differences so that the intervals of the variables fall within [0,1]; eliminate non-uniform plane points and pseudo-plane points, and merge the uniform plane points into the cylindrical units after preliminary filtering to obtain optimized plane point cloud data;

[0019] S24: Use the supervoxel clustering method to distinguish the road surface part from other discrete plane points to obtain the complete road surface point cloud data.

[0020] Furthermore, in step S21, the process of rasterizing the original three-dimensional point cloud data includes the following steps:

[0021] S211, assuming that the original 3D point cloud data space coordinates are (X, Y, Z) T , create a temporary array (x,y,z) T :

[0022]

[0023] S212, assuming (min x, max y) T and (max x,min y) T They are the starting point and end point of the linear index calculation. After the point cloud data is rasterized, the horizontal distance D between the starting point and the end point is calculated by the formula x The vertical distance D y , where D x With D y All are positive integers, and the calculation formula is as follows:

[0024]

[0025] S213, establish a size of (1+D y )×(1+D x )’s zero value matrix Ψ and the empty set cell matrix Ω; calculate the horizontal distance D from any grid point to the starting point of the data index x and vertical distance D y ; Use line number 1+d y and column number 1+d x Represents the elements of the zero-value matrix Ψ and the empty set cell matrix Ω; by the formula ζ = d x ·(1+D y )+dy +1 converts row and column numbers into linear indices represented by ζ;

[0026] When the point cloud data is divided into blocks, the point cloud data is arranged in ascending order according to the value of the linear index ζ. The difference Δζ between adjacent data points is calculated to determine whether the point is a step point. If Δζ < 1, the point is not a step point, otherwise it is a step point. The step point ζ is j With the previous step point ζ j-1 The point cloud data between is stored in the cell Ω matrix ζ j-1 In each unit, give the ζth j-1 The element is assigned the value 1.

[0027] Further, in step S23, the processing process of the height rectangular diagram sorting method includes the following steps:

[0028] S231, the height δ of a single stripe in the height histogram h Set it to 1m~4m, create a frequency histogram along the Z direction, and create an empty cell array Φ;

[0029] S232, sort the frequencies in descending order, set F{F1,F2...F i ...F n |1≤i≤n,i,n∈N +} is the height value corresponding to each stripe after arrangement;

[0030] S233, establish a loop i=2, k=1, and set F i-1 The corresponding point cloud data are merged and stored in Φ{k};

[0031] S234, such as F i -F i-1 ≤δ h , then F i With F i-1 The corresponding point cloud data is merged and stored in Φ{k}, otherwise k=k+1, and F i The corresponding point cloud data is stored in Φ{k};

[0032] S235, let i=i+1, if i≤N h , then return to step S234, otherwise the loop ends.

[0033] Furthermore, in step S24, the processing process of the supervoxel clustering method is as follows:

[0034] Using the index position information, the corresponding three-dimensional point cloud data in the cell matrix is ​​retrieved, and the three-dimensional point cloud data is spliced ​​to obtain the complete road surface point cloud data in the final road surface supervoxel clustering result, which distinguishes the road surface part from other discrete plane points.

[0035] Furthermore, in step S4, the three-dimensional coordinate points in space are converted into two-dimensional coordinate points through the cylindrical perspective projection model, and the processing process of converting the three-dimensional problem into a two-dimensional problem is as follows:

[0036] Set the viewpoint starting point to the origin and the coordinates to (x m ,y m ,z m ) T , the vehicle's forward direction is set to the Y' axis, the X' axis is parallel to the XOY plane and perpendicular to the Y axis, and the Z' axis is set perpendicular to the X'OY' plane to construct a three-dimensional space coordinate system;

[0037] Construct a corresponding two-dimensional coordinate system u-υ, the origin of which is set on the Y' axis and at a distance R from the origin of the sight line, the υ axis is set parallel to the Z' axis along the vector, and the u axis is set as a clockwise arc segment, which is perpendicular to both the υ axis and the Y' axis;

[0038] In the two spatial coordinate systems, the coordinate calculation formula of the two-dimensional point on the cylindrical surface is as follows:

[0039]

[0040]

[0041]

[0042] in:

[0043] (x, y, z) T represents the coordinates of points around the driver in geodetic space;

[0044] (x m ,y m , z m ) T Represents the coordinates of the sight origin;

[0045] (x′, y′, z′) T Represents the coordinates converted to the X′Y′Z′ coordinate system;

[0046] λ represents the scale factor between the two reference coordinate systems and its value is 1.0;

[0047] θ represents the rotation matrix of the line of sight rigid transformation;

[0048] θ x ,θ y ,θ z Represents the rotation matrix around the X, Y and Z axes;

[0049] α mRepresents the rotation angle around the Z axis;

[0050] γ m Represents the rotation angle around the X axis;

[0051] R represents the radius of the cylindrical surface, which is 3.0;

[0052] (u,υ) T represents the coordinates of the cylindrical surface;

[0053] α m and γ m They represent the azimuth and vertical angle of the origin of the line of sight respectively. The direction of the Y′ axis is consistent with the forward direction of the vehicle. Therefore, in this case, the direction of the Y′ axis is consistent with the forward direction of the vehicle.

[0054] In the locally constructed local three-dimensional space, the coordinate calculation formula of all sight end points is as follows:

[0055]

[0056] η=θ / θ res +1

[0057] [ρ]={ρ1,ρ2...ρ j ... η-1 , ρ η |ρ i =D(1≤j≤η), η=θ / θ res +1}

[0058] [y′ e , x′ e ]=Polar To Cartesian([θ,ρ])

[0059] [z′ e ]={z1,z2...z j ...z η-1 , z η |z i =D·tanθ v (1≤j≤η), η=θ / θ res +1}

[0060] in:

[0061] [x′ e ,y′e,z′ e ] T Represents the 3D coordinates of the end point of the line of sight in the local 3D coordinate system;

[0062] θ represents the driver’s perspective;

[0063] θ resrepresents the angular distance between sight lines;

[0064] j represents the sight index, when 1≤j≤η, then η is the sight index;

[0065] The area where the projection point coordinates are located is regarded as having a side length of W i The projection point is located at the center of the square area, and the coordinates are set to Set up collection Store the two-dimensional coordinates of all sight endpoints, and set the target point coordinates to Ob i (x i ,y i ,z i ) T , through the KD tree algorithm, with parameter W i Chebyshev sequence search for / 2 Points in the neighborhood;

[0066] Set Collection i Store the two-dimensional coordinates of all points in the square area where the projection point is located, the set Ψ i All two-dimensional points have corresponding three-dimensional coordinates; let κ i represents the three-dimensional space distance between the sight origin and the target point, and any three-dimensional space distance from the sight origin exceeding κ i The points are excluded from the set Ψ i The beneficial effects of the present invention are:

[0067] 1. The present invention proposes to fix a small laser radar on a gimbal, rotate it multiple times to collect point cloud data and perform data stitching, thereby solving the problem that the small laser radar has a small viewing angle and incomplete and inaccurate data collection.

[0068] 2. The point cloud data sub-unit division method based on linear indexing proposed in the present invention improves the efficiency of data segmentation and indexing.

[0069] 3. The present invention proposes the use of preliminary filtering method, height histogram sorting method and supervoxel clustering method to achieve fast and efficient segmentation and extraction processing of road surface points in point cloud data.

[0070] 4. The present invention proposes to use a digital elevation model to convert the spatial coordinate system and the local plane coordinate system according to a cylindrical perspective projection model, construct a visual domain model, analyze and determine the location of obstacles that affect the driver's visibility, and analyze the visualization of the intersection.

[0071] 5. The present invention improves the processing and calculation efficiency of point cloud data through a three-dimensional line of sight calculation method, and solves the problem of processing massive point cloud data.

[0072] 6. The present invention realizes the transformation of three-dimensional coordinate points in space into two-dimensional coordinate points through a cylindrical perspective projection model, thereby converting a three-dimensional problem into a two-dimensional problem. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] Figure 1 It is a structural schematic diagram of the urban road intersection visualization evaluation method based on point cloud data of the present invention.

[0074] Figure 2 The present invention is a schematic diagram of the structure of the method for constructing a digital elevation model.

[0075] Figure 3 Schematic diagram of the effect of the preliminary filtering method of the present invention.

[0076] Figure 4 It is a structural schematic diagram of the height histogram method of the present invention.

[0077] Figure 5 It is a structural schematic diagram of the cylindrical perspective projection calculation of the present invention.

[0078] Figure 6 This is a structural diagram of the visual field blind area of ​​the present invention.

[0079] Figure 7 It is a corresponding schematic diagram of the visual range and viewing distance curve diagram of the present invention. DETAILED DESCRIPTION

[0080] The present invention will now be described in further detail with reference to the accompanying drawings.

[0081] It should be noted that the terms such as "upper", "lower", "left", "right", "front", "back", etc. cited in the invention are only for the convenience of description and are not used to limit the scope of implementation of the present invention. Changes or adjustments in their relative relationships should be regarded as the scope of implementation of the present invention without substantially changing the technical content.

[0082] Embodiment 1

[0083] Figure 1 It is a structural schematic diagram of the urban road intersection visualization evaluation method based on point cloud data of the present invention. Figure 2 This is a schematic diagram of the structure of the method for constructing a digital elevation model of the present invention. This embodiment is applicable to a method for calculating the three-dimensional sight distance of an urban road intersection based on point cloud data, and the method for calculating the three-dimensional sight distance of an urban road intersection includes the following steps:

[0084] S1, collect point cloud data of road intersections.

[0085] S2, performing road surface recognition and segmentation processing on the point cloud data, separating the road surface part from other discrete plane points to obtain complete road surface point cloud data.

[0086] S3, establish a plane grid dot matrix, use the identified road surface points as the benchmark interpolation reference, use the neighboring point interpolation algorithm to obtain the height information of the plane grid points, and build a digital elevation model.

[0087] S4, quantifies the concept of field of view into a fan-shaped shape composed of countless rays emitted from the origin of the driver's line of sight. The angle of the fan-shaped shape is related to the driver's viewing angle range. Through the cylindrical perspective projection model, the three-dimensional coordinate points in space are converted into two-dimensional coordinate points, and the three-dimensional problem is converted into a two-dimensional problem.

[0088] A coordinate system is constructed with the starting point of the driver's line of sight as the origin, the three-dimensional coordinate system is converted into a two-dimensional coordinate system, and then the two-dimensional rectangular coordinate system is converted into a polar coordinate system to determine the position of the obstacle blocking the driver's line of sight.

[0089] S5, analyze the variable factors that may affect the visibility of urban road intersections.

[0090] 1. Collecting point cloud data of road intersections

[0091] In the above step S1, the process of collecting point cloud data of the road intersection includes the following sub-steps:

[0092] S11, fix the small laser radar on the rotating gimbal, place it at the center of the urban road intersection, and keep the laser radar mounting bracket fixed; collect data multiple times through the rotating gimbal to obtain the sub-file of the overall point cloud data of the intersection, and the overlap degree of two adjacent point cloud files is 1 / 3.

[0093] S12, performing data splicing on the obtained point cloud files through corresponding feature points to obtain complete point cloud data of urban road intersections.

[0094] 2. Road surface recognition and segmentation processing of point cloud data

[0095] In step S2, the process of performing road surface recognition and segmentation processing on the point cloud data includes the following sub-steps:

[0096] S21: Rasterize the original 3D point cloud data. Each grid point of the rasterized point cloud data has a corresponding linear index value. The data points at the same position in the same matrix will have the same linear index value. When a data point crosses a matrix, the linear index value ζ of the data point will jump. The specific rasterization processing steps are as follows:

[0097] First, the original 3D point cloud data must be rasterized. Assume that the spatial coordinates of the original 3D point cloud data are (X, Y, Z) T , create a temporary array (x,y,z) T , (x,y,z) TIt can be calculated by the following formula:

[0098]

[0099] Where: [.] represents the integer symbol, ε x , ε y , ε z Represents the grid interval along the X-axis, Y-axis and Z-axis respectively. The smaller the interval, the smaller the size of the block, and vice versa. The block division process in different dimensions is similar, so this article takes the division of two-dimensional cylindrical units as an example to explain. Assume (min x, max y) T and (min x, min y) T They are the starting point and end point of linear index calculation. After the point cloud data is rasterized, the formula The horizontal distance D between the starting and ending points can be obtained by calculation x The vertical distance D y , D x With D y are all positive integers. At the same time, the size is (1+D y )×(1+D x ) of the zero value matrix Ψ and the empty set cell matrix Ω. The horizontal distance d from any grid point to the starting point of the data index x and vertical distance d y All are calculated in the same way. Therefore, in the computer storage system, the elements in the zero-value matrix Ψ and the empty set cell matrix Ω can be considered to use row numbers 1+d y and column number 1+d x Indicates (0<d x ≤D x ,0<d y ≤D y ). In order to facilitate data indexing, the formula ζ=d x ·(1+D y )+d y +1 converts the row and column numbers into linear indices represented by ζ. Each grid point of the point cloud data after rasterization has a corresponding linear index value. The data points at the same position in the same matrix will have the same linear index value, and each time the data point crosses a matrix, the ζ of the data point will have a step, which we call a step point, and ζ is used to represent the linear index. j Definition. When dividing the point cloud data into blocks, the point cloud data can be first sorted in ascending order according to the value of the linear index ζ, and then the difference Δζ between adjacent data points is calculated to determine whether the point is a step point. If Δζ < 1, the point is not a step point, otherwise it is a step point. For all step points ζ j The same process is performed, that is, the distance between it and the previous step point is ζ j-1 The point cloud data is stored in the ζth cell Ω matrixj-1 In each unit, give the ζth j-1 The value of each element is 1, thereby achieving the block division of the point cloud data.

[0100] S22: Using the height threshold h t The points with lower elevations in the cylindrical unit are screened out, and the remaining data points are analyzed using the principal component analysis method to obtain the main direction of the cylindrical unit.

[0101] Figure 3 This is a schematic diagram of the effect of the preliminary filtering method of the present invention. In this embodiment, the unit grid is divided into 0.5m*0.5m. In each cylindrical unit, the height range of the road surface point cloud data is extremely limited and approximately equal, so the height threshold h can be used. t First, select the points with lower elevation in the columnar unit. Then use the principal component analysis method to analyze the remaining data points to obtain their main directions. Principal component analysis is a commonly used statistical method, that is, a set of variables that may be correlated is converted into linearly uncorrelated variables, namely, principal components, through matrix orthogonal transformation. Specifically in this embodiment, it is assumed that (x, y, z) T are the coordinates of any point in the cylindrical unit, (λ,u,υ) T is the converted variable. For a point cloud with good plane features, its main direction v is consistent with the plane normal vector it constitutes, that is, they are perpendicular to each other. Therefore, compared with non-feature points, the average absolute error in the υ direction is smaller. The average absolute error is simply the average of the absolute values ​​of the deviations of all individual data from their overall arithmetic mean. For a set of m data, the average absolute error is calculated as calculate, Since the elevation of the ground point is at a relatively low position in all point cloud data, the height threshold is set to 0.2 m in this embodiment, that is, the size of each cylindrical unit is 0.5 m*0.5 m*0.2 m.

[0102] S23: Use the height histogram method to sort the detected road surface points and extract features from the obtained plane point cloud data; the feature variables include the number of point clouds per unit area, as well as the standard deviation and median of the difference between adjacent items along the x and y directions used to describe the distribution characteristics of the point cloud on the xy plane, and obtain an observation value array consistent with the number of cylindrical units; use the unsupervised classification K-Means clustering algorithm to divide the cylindrical units into three categories: uniform plane points, uneven plane points and pseudo-plane points; standardize all variables to eliminate their differences so that the intervals of the variables fall within [0,1]; eliminate uneven plane points and pseudo-plane points, merge the uniform plane points into the cylindrical units after preliminary filtering, and obtain optimized plane point cloud data.

[0103] Figure 4Schematic diagram of the structure of the height histogram method of the present invention. In step S23 of this embodiment, the processing process of the height rectangular diagram sorting method includes the following steps:

[0104] S231, the height of a single stripe in the height histogram δ h Set it to 1m~4m, create a frequency histogram along the Z direction, and create an empty cell array Φ.

[0105] S232, sort the frequencies in descending order, set F{F1,F2...F i ...F n |1≤i≤n,i,n∈N +} is the height value corresponding to each stripe after arrangement.

[0106] S233, establish a loop i=2, k=1, and set F i-1 The corresponding point cloud data are merged and stored in Φ{k}.

[0107] S234, such as F i -F i-1 ≤δ h , then F i With F i-1 The corresponding point cloud data is merged and stored in Φ{k}, otherwise k=k+1, and F i The corresponding point cloud data is stored in Φ{k}.

[0108] S235, let i=i+1, if i≤N h , then return to step S234, otherwise the loop ends.

[0109] Considering that the algorithm process in each cylindrical unit is similar and independent of each other, this embodiment makes full use of the parallel computing advantages of multi-threaded processors to improve the efficiency of automatic detection of plane points in each cylindrical unit. Considering the differences in data points in different cylindrical units, this embodiment also extracts features from the plane point cloud data obtained by the above algorithm. The extracted feature variables are mainly the standard deviation and median of the difference between adjacent items along the x and y directions. When calculating, they are first arranged in ascending order and then the standard deviation and median are calculated. The feature variable mainly describes the distribution characteristics of the point cloud on the xy plane. Because the distribution of road surface points is continuous and regular, while the distribution of non-road surface points is uneven. In addition, considering that the road surface points in some cylindrical units are not completely evenly distributed, the point cloud density of the number of point clouds per unit area is also used as one of the feature variables. The point cloud density calculation of the number of point clouds per unit area adopts the method of block calculation and then taking the average value. Each cylindrical unit is evenly divided into 100 sub-areas of 10*10. The point cloud density is equal to the ratio of the area of ​​the sub-area containing the point to the total number of points.

[0110] After the above-mentioned feature vector extraction, an array of observation values ​​consistent with the number of cylindrical units is obtained. Then, the cylindrical units can be divided into three categories: uniform plane points, non-uniform plane points, and pseudo-plane points using the unsupervised classification K-Means clustering algorithm. Among them, uniform plane points are distributed regularly in the unit, and non-uniform plane points are distributed in a disorderly manner; pseudo-plane points are points that are mistakenly classified as plane points due to a small height range in the Z direction, so they present a linear structure. First, all variables need to be standardized to eliminate their differences so that the intervals of the variables all fall within [0,1]. Then, the non-uniform plane points and pseudo-plane points are eliminated, and the uniform plane points are incorporated into the cylindrical units after preliminary filtering to obtain the optimized plane point cloud data.

[0111] S24: Use the supervoxel clustering method to distinguish the road surface from other discrete plane points to obtain complete road surface point cloud data, use the index position information to retrieve the corresponding three-dimensional point cloud data in the cell matrix, and combine them to obtain the final road surface supervoxel clustering result. The supervoxel-based clustering method can effectively distinguish the road surface from other discrete plane points and obtain complete road surface point cloud data.

[0112] 3. Building a Digital Elevation Model

[0113] In this embodiment, the construction of the digital elevation model includes the following sub-steps:

[0114] On the plane with (minx, maxy) T and (maxx,miny) T As the corner point, a plane grid matrix is ​​established with an interval of 0.2m;

[0115] The identified road surface points are used as the reference for interpolation, and the height information of the plane grid points is obtained using the neighboring point interpolation algorithm. When interpolating under this digital elevation model, the neighboring grid points can be located based on the plane position of the observation point, and the height information of the observation point position can be quickly obtained to facilitate subsequent visibility analysis.

[0116] 4. Convert 3D Problems to 2D Problems

[0117] Figure 5 It is a structural schematic diagram of the cylindrical perspective projection calculation of the present invention. Figure 6 This is a structural diagram of the blind area of ​​the visual field of the present invention. In step S4 of this embodiment, the three-dimensional coordinate points in space are converted into two-dimensional coordinate points through a cylindrical perspective projection model. The processing process of converting a three-dimensional problem into a two-dimensional problem is as follows:

[0118] In a two-dimensional plane, the concept of vision is quantified as a fan composed of countless rays emitted from the origin of the driver's line of sight. The distance between the starting and ending points of the line of sight is the fan radius D, and the driver's viewing angle is θ. Both D and θ are set as variable parameters for a more comprehensive evaluation. In this fan-shaped area, when the line of sight is blocked, the obstacle will be detected, and θ is set to 120°.

[0119] Through the cylindrical perspective projection model, the transformation of three-dimensional coordinate points in space to two-dimensional coordinate points is realized, and the three-dimensional problem is converted into a two-dimensional problem. First, two local coordinate systems are established to obtain the coordinate data of the cylindrical surface points. Now set the viewpoint starting point to the origin and the coordinates to (x m ,y m ,z m ) T , the vehicle's forward direction is set to the Y' axis, the X' axis is parallel to the XOY plane and perpendicular to the Y axis, and the Z' axis is set perpendicular to the X'OY' plane. On the basis of this three-dimensional space, another two-dimensional coordinate system u-υ can be constructed. The origin of the two-dimensional coordinate system is set on the Y' axis and the distance from the origin of the line of sight is R. The υ axis is set parallel to the Z' axis along the vector, and the u axis is set to a clockwise arc segment, which is perpendicular to both the υ axis and the Y' axis. In the two spatial coordinate systems, the coordinates of the two-dimensional point on the cylindrical surface can be calculated by the following formula:

[0120]

[0121]

[0122] Where: (x, y, z) T represents the coordinates of points around the driver in geodetic space; (x m ,y m ,z m ) T Represents the coordinates of the origin of the line of sight; (x′, y′, z′) T represents the coordinates converted to the X'Y'Z' coordinate system; λ represents the scale factor between the two reference coordinate systems, which is 1.0; θ represents the rotation matrix of the rigid transformation of the line of sight; θ x ,θ y ,θ z Represents the rotation matrix around the X, Y, and Z axes; α m Represents the rotation angle around the Z axis; γ m represents the rotation angle around the X axis; R represents the radius of the cylindrical surface, which is 3.0; (u,υ) T represents the coordinates of the cylindrical surface; α m and γ mRespectively represent the azimuth and vertical angle of the origin of the sight line, and the direction of the Y' axis is consistent with the vehicle's forward direction, so in this case, the direction of the Y' axis is consistent with the vehicle's forward direction. Therefore, even if the alignment of the road has horizontal or vertical curves, the perspective projection can still work. Subsequently, in the local 3D space constructed locally, the coordinates of all sight end points are calculated using the following formula.

[0123]

[0124] η=θ / θ res +1

[0125] [ρ]={ρ1,ρ2…ρ j …ρ η-1 ,ρ η |ρ i =D(1≤j≤η),η=θ / θ res +1}

[0126] [y′ e ,x′ e ]=Polar To Cartesian([θ,ρ])

[0127] [z′ e ]={z1,z2…z j …z η-1 ,z η |z i =D·tanθ v (1≤j≤η),η=θ / θ res +1}

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

[0129] [x′ e ,y′ e ,z′ e ] T represents the three-dimensional coordinates of the end point of the line of sight in the local three-dimensional coordinate system; θ represents the driver's perspective; θ res represents the angular distance between lines of sight; j represents the line of sight index, when 1≤j≤η, then η is the line of sight index.

[0130] Through the above calculations, the two-dimensional coordinates of the points representing the road environment and all the sight end points can be obtained.

[0131] The next step is to obtain Figure 5 The projection point shown is the coordinate of the projection point of the target point on the cylindrical surface. The area where the projection point coordinates are located can be regarded as a side length of W. i The projection point is located at the center of the square area, and the coordinates are set to Set up collection Store the two-dimensional coordinates of all sight endpoints, and set the target point coordinates to Ob i (x i ,y i ,z i ) T , and then through the KD tree algorithm, with parameter W i Chebyshev sequence search for / 2 Points in the neighborhood. Set the set Ψ i Store the two-dimensional coordinates of all points in the square area where the projection point is located, the set Ψ i All two-dimensional points have corresponding three-dimensional coordinates. Let κ i Represents the three-dimensional space distance between the sight origin and the target point. Any three-dimensional space distance from the sight origin exceeds κ i The points will be excluded from the set Ψ i In addition, this can greatly improve the efficiency of searching along the line of sight.

[0132] In this embodiment, the process of visual domain modeling mainly includes: constructing a coordinate system with the starting point of the driver's line of sight as the origin, converting the three-dimensional coordinate system to the two-dimensional coordinate system, and converting the two-dimensional rectangular coordinate system to the polar coordinate system. The position of the obstacle blocking the driver's line of sight is determined through two coordinate conversions. Since the coordinate system is converted based on a specific matrix equation, when the coordinates in the polar coordinate system are mapped to the real three-dimensional space coordinate system, the construction of the three-dimensional visual domain model in the coordinate system can be realized.

[0133] Figure 7 The corresponding schematic diagram of the visual range and the sight distance curve of the present invention. According to the corresponding schematic diagram of the visual range modeling and the sight distance curve obtained in the embodiment, the position of obstacles can be quickly located and the sight distance can be quantitatively calculated through the present invention, so as to obtain the visibility evaluation of the urban road intersection.

[0134] 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 visual evaluation method for urban road intersections based on point cloud data, characterized in that: The method for calculating the three-dimensional sight distance at a road intersection comprises the following steps: S1, collects point cloud data of road intersections; S2, performing road surface recognition and segmentation processing on the point cloud data, distinguishing the road surface part from other discrete plane points to obtain complete road surface point cloud data; S3, establishing a plane grid dot matrix, using the identified road surface points as the base interpolation reference, using the neighboring point interpolation algorithm to obtain the height information of the plane grid points, and constructing a digital elevation model; S4, quantifies the concept of visual field into a fan-shaped shape composed of countless rays emitted from the origin of the driver's line of sight. The angle of the fan-shaped shape is related to the driver's visual range. Through the cylindrical perspective projection model, the three-dimensional coordinate points in space are converted into two-dimensional coordinate points, and the three-dimensional problem is converted into a two-dimensional problem. A coordinate system is constructed with the starting point of the driver's line of sight as the origin, the three-dimensional coordinate system is converted into a two-dimensional coordinate system, and then the two-dimensional rectangular coordinate system is converted into a polar coordinate system to determine the position of the obstacle blocking the driver's line of sight; S5, analyze the variable factors that may affect the visibility of urban road intersections; In step S4, the three-dimensional coordinate points in space are converted into two-dimensional coordinate points through the cylindrical perspective projection model, and the processing process of converting the three-dimensional problem into a two-dimensional problem is as follows: Set the viewpoint starting point to the origin and the coordinates to (x m ,y m , z m ) T , the vehicle's forward direction is set to the Y′ axis, the X′ axis is parallel to the XOY plane and perpendicular to the Y axis, and the Z′ axis is set perpendicular to the X′OY′ plane to construct a three-dimensional space coordinate system; Construct a corresponding two-dimensional coordinate system u-υ, the origin of which is set on the Y′ axis and at a distance R from the origin of the line of sight, the v axis is set parallel to the Z′ axis along the vector, and the u axis is set as a clockwise arc segment, which is perpendicular to both the v axis and the Y′ axis; In the two spatial coordinate systems, the coordinate calculation formula of the two-dimensional point on the cylindrical surface is as follows: Where: (x, y, z) T represents the coordinates of points around the driver in geodetic space; (x m ,y m , z m ) T Represents the coordinates of the origin of the line of sight; (x′, y′, z′) T represents the coordinates converted to the X′Y′Z′ coordinate system; λ represents the scale factor between the two reference coordinate systems, which is 1.0; θ represents the rotation matrix of the line of sight rigid transformation; θ x ,θ y ,θ z represents the rotation matrix around the X, Y and Z axes; α represents the rotation angle around the Z axis; γ m represents the rotation angle around the X axis; R represents the radius of the cylindrical surface, which is 3.0; (u, v) T represents the coordinates of the cylindrical surface; α m and γ m They represent the azimuth and vertical angle of the origin of the line of sight respectively. The direction of the Y′ axis is consistent with the forward direction of the vehicle. Therefore, in this case, the direction of the Y′ axis is consistent with the forward direction of the vehicle. In the locally constructed local three-dimensional space, the coordinate calculation formula of all sight end points is as follows: n = θ / θ res +1 [ρ]={ρ1,ρ2...ρ j ...r η-1 ,r η |r i =D,η=θ / θ res +1} [y′ e ,x′ e ]=Polar To Cartesian([θ,ρ]) [With' e ]={z1,z2...z j ...With η-1 ,With η |from i =D tanθ v ,η=θ / θ res +1} Where: [x′ e , y′ e , z′ e ] T represents the three-dimensional coordinates of the end point of the line of sight in the local three-dimensional coordinate system; θ represents the driver's perspective; θ res represents the angular spacing between sight lines; ρ represents the distance parameter in polar coordinates after coordinate conversion; D represents the sector radius, i.e., the driver’s sight distance; j represents the sight index, when 1≤j≤η, then η is the sight index; The area where the projection point coordinates are located is regarded as having a side length of W i The projection point is located at the center of the square area, and the coordinates are set to Setting up the collection Store the two-dimensional coordinates of all sight endpoints, and set the target point coordinates to Ob i (x i ,y i , z i ) T , through the KD tree algorithm, with parameter W i Chebyshev sequence search for / 2 Points in the neighborhood; Set Collection i Store the two-dimensional coordinates of all points in the square area where the projection point is located, the set Ψ i All two-dimensional points have corresponding three-dimensional coordinates; let κ i represents the three-dimensional space distance between the sight origin and the target point, and any three-dimensional space distance from the sight origin exceeding κ i The points are excluded from the set Ψ i Besides.

2. The method for visual evaluation of urban road intersections based on point cloud data according to claim 1, characterized in that: In step S1, the process of collecting point cloud data of a road intersection includes the following steps: S11, fix the small laser radar on the rotating gimbal, place it at the center of the urban road intersection, and keep the laser radar mounting bracket fixed; collect data multiple times through the rotating gimbal to obtain the sub-file of the overall point cloud data of the intersection, and the overlap degree of two adjacent point cloud files is 1 / 3; S12, performing data splicing on the obtained point cloud files through corresponding feature points to obtain complete point cloud data of urban road intersections.

3. The method for visual evaluation of urban road intersections based on point cloud data according to claim 1, characterized in that: In step S2, the process of performing road surface recognition and segmentation processing on the point cloud data includes the following steps: S21: Rasterize the original 3D point cloud data. Each grid point of the rasterized point cloud data has a corresponding linear index value. The data points at the same position in the same matrix will have the same linear index value. When a data point crosses a matrix, the linear index value ζ of the data point will jump. S22: Using the height threshold h t The points with lower elevation in the columnar unit are screened out, and the remaining data points are analyzed using the principal component analysis method to obtain the main direction of the columnar unit; S23: Use the height histogram method to sort the detected road surface points, and extract features from the obtained plane point cloud data; the feature variables include the number of point clouds per unit area, and the standard deviation and median of the difference between adjacent items along the x and y directions used to describe the distribution characteristics of the point cloud on the xy plane, and obtain an observation value array consistent with the number of cylindrical units; use the unsupervised classification K-Means clustering algorithm to divide the cylindrical units into three categories: uniform plane points, non-uniform plane points, and pseudo-plane points; standardize all variables to eliminate their differences, so that the intervals of the variables all fall within [0, 1]; eliminate non-uniform plane points and pseudo-plane points, and merge the uniform plane points into the cylindrical units after preliminary filtering to obtain optimized plane point cloud data; S24: Use the supervoxel clustering method to distinguish the road surface part from other discrete plane points to obtain the complete road surface point cloud data.

4. The method for visual evaluation of urban road intersections based on point cloud data according to claim 3 is characterized in that: In step S21, the process of rasterizing the original three-dimensional point cloud data includes the following steps: S211, assuming that the original 3D point cloud data space coordinates are (X, Y, Z) T , with ε x , ε y , ε z Represents the grid intervals along the X-axis, Y-axis, and Z-axis respectively, creating a temporary array (x, y, z) T : S212, assuming (minx, maxy) T and (maxx, miny) T They are the starting point and end point of the linear index calculation. After the point cloud data is rasterized, the horizontal distance D between the starting point and the end point is calculated by the formula x The vertical distance D y , where D x With D y All are positive integers, and the calculation formula is as follows: S213, establish a size of (1+D y )×(1+D x )’s zero value matrix Ψ and the empty set cell matrix Ω; calculate the horizontal distance D from any grid point to the starting point of the data index x and vertical distance D y ; Use line number 1+d y and column number 1+d x Represents the elements of the zero-value matrix Ψ and the empty set cell matrix Ω; by the formula ζ = d x ·(1+D y )+d y +1 converts row and column numbers into linear indices represented by ζ; When the point cloud data is divided into blocks, the point cloud data is arranged in ascending order according to the value of the linear index ζ. The difference Δζ between adjacent data points is calculated to determine whether the point is a step point. If Δζ<1, the point is not a step point, otherwise it is a step point. The step point ζ is j With the previous step point ζ j-1 The point cloud data between is stored in the cell Ω matrix ζ j-1 In each unit, give the ζth unit of the Ψ matrix j-1 The element is assigned the value 1.

5. The method for visual evaluation of urban road intersections based on point cloud data according to claim 2 is characterized in that: In step S23, the processing process of the height rectangular diagram sorting method includes the following steps: S231, the height δ of a single stripe in the height histogram h Set it to 1m~4m, create a frequency histogram along the Z direction, and create an empty cell array Φ; S232, sort the frequencies in descending order, set F{F1, F2...F i ...F n |1≤i≤n,i,n∈N + } is the height value corresponding to each stripe after arrangement; S233, establish a loop i=2, k=1, and set F i-1 The corresponding point cloud data are merged and stored in Φ{k}; S234, such as F i -F i-1 ≤δ h , then F i With F i-1 The corresponding point cloud data are merged and stored in Φ{k}, otherwise k=k+1, and F i The corresponding point cloud data is stored in Φ{k}; S235, let i=i+1, if i≤N h , then return to step S234, otherwise the loop ends.

6. The method for visual evaluation of urban road intersections based on point cloud data according to claim 3 is characterized in that: In step S24, the processing process of the supervoxel clustering method is as follows: Using the index position information, the corresponding three-dimensional point cloud data in the cell matrix is ​​retrieved, and the three-dimensional point cloud data is spliced ​​to obtain the complete road surface point cloud data in the final road surface supervoxel clustering result, which distinguishes the road surface part from other discrete plane points.

Images