Automatic digital modeling method, system and equipment for urban road

Abstract

The invention provides an automatic digital modeling method, system and equipment for an urban road. The method comprises the following steps: synchronously acquiring a road laser point cloud and a panoramic image; ground point extraction is carried out on the laser point cloud, and space-time registration is carried out on the ground point cloud and the panoramic image; road element extraction is carried out based on the registered ground point cloud, road elements are distributed to a vector layer according to types, and semantic attributes and unique codes are configured; calling a matched model template from the parameterized template library according to the semantic attributes of the road elements, and generating a road element three-dimensional model based on the model template; mapping the image texture of the image source matched with each model surface to the corresponding model surface according to the visibility of the model surface in the panoramic image to obtain a three-dimensional model fused with texture information; and carrying out superposition display on the three-dimensional model fused with the texture information, the original laser point cloud and the panoramic image in the same three-dimensional scene so as to carry out model verification, correcting an unqualified part, and obtaining a three-dimensional model of the target road.

Description

Technical Field

[0001] This invention relates to the fields of computer vision and laser point cloud processing technology, and in particular to an automated digital modeling method, system and device for urban roads. Background Technology

[0002] With the deepening of digital twin city and new infrastructure construction, the demand for high-precision, high-efficiency, and full-element digital modeling of urban road infrastructure is becoming increasingly urgent. Traditional 3D road modeling methods mainly rely on manual field surveys and manual 3D modeling, which have drawbacks such as low efficiency, high cost, long construction period, and difficulty in unifying accuracy.

[0003] Currently, although some studies use airborne lidar or oblique photography for large-scale terrain and building modeling, the ability to model roads and their ancillary facilities (such as road markings, streetlights, signs, guardrails, etc.) is insufficient. Airborne data has blind spots in acquiring facade and bottom information, making it difficult to meet the accuracy requirements of road component-level modeling. Another approach is to use a vehicle-mounted mobile measurement system (MMS) to acquire high-precision point cloud and image data. Existing technologies usually use point clouds to extract road frameworks and images to textures, but the two have not achieved deep coupling, resulting in the following prominent problems: (1) Low degree of automation: the identification and extraction of road elements (such as lane lines and signs) still require a lot of manual intervention; (2) Lack of semantic information: the generated models are mostly "white models" or "oblique photography models", which cannot distinguish component types, lack attribute information, and are difficult to integrate with business management systems (BIM / GIS); (3) Fragmented modeling process: vector extraction, model construction, texture mapping, and other links are independent of each other and have not formed an integrated automatic process, making it difficult to promote and apply on a large scale. Summary of the Invention

[0004] In view of the above problems, the present invention is proposed to provide an automated digital modeling method, system and device for urban roads that overcomes or at least partially solves the above problems.

[0005] One aspect of the present invention provides an automated digital modeling method for urban roads, the method comprising: Acquire the original laser point cloud and panoramic image of the target road simultaneously; Ground points are extracted from the original laser point cloud, and the extracted ground point cloud is spatiotemporally registered with the panoramic image. Based on the spatiotemporally registered ground point cloud, extract road features such as lane lines, road center lines, directional arrows, and curb stones. Based on the non-ground point cloud in the original laser point cloud, guardrail road elements and road ancillary facilities road elements are extracted, and the location and geometric attributes of the road ancillary facilities road elements are identified. The lane surface road features are constructed based on the extracted lane line road features, curb road features, and guardrail road features. The extracted road features are assigned to the corresponding vector layers according to their type and configured with structured semantic attributes and globally unique codes. The semantic attributes include information related to the road segment, feature type information, location information, and / or geometric parameter information. Based on the semantic attributes of the road elements, a matching model template is called from the parametric template library, and a 3D model of the corresponding road elements is generated based on the corresponding model template. Based on the visibility of the surface of the 3D model of each road element in the panoramic image, an image source matching each model surface is selected, and the image texture of the image source is mapped to the corresponding model surface to obtain a 3D model with fused texture information. The 3D model with fused texture information is overlaid with the original laser point cloud and panoramic image in the same 3D scene to verify the geometric accuracy, semantic attribute consistency and texture quality of the 3D model, and to correct the unqualified parts of the verification results, so as to obtain the 3D model of the target road.

[0006] Another aspect of the present invention provides an automated digital modeling system for urban roads, the system comprising: The multi-source data acquisition module is used to acquire the original laser point cloud and panoramic image of the target road collected simultaneously. The collaborative preprocessing module is used to extract ground points from the original laser point cloud and perform spatiotemporal registration between the extracted ground point cloud and the panoramic image. The first element extraction module is used to extract lane line road elements, road center line road elements, directional arrow road elements, and curbstone road elements from the ground point cloud after spatiotemporal registration. The second feature extraction module is used to extract guardrail road features and road ancillary facility road features based on the non-ground point cloud in the original laser point cloud, and to identify the location and geometric attributes of the road ancillary facility road features. The feature construction module is used to construct lane surface road features based on the extracted lane line road features, curb road features, and guardrail road features; The semantic annotation module is used to assign the extracted road features to the corresponding vector layers according to their type and configure structured semantic attributes and globally unique codes. The semantic attributes include information related to the road segment, feature type information, location information and / or geometric parameter information. The road element modeling module is used to call the matching model template from the parametric template library according to the semantic attributes of the road elements, and generate the corresponding 3D model of the road elements based on the corresponding model template. The model texture fusion module is used to select image sources that match the surface of each model based on the visibility of the surface of the 3D model of each road element in the panoramic image, and map the image texture of the image source to the corresponding model surface to obtain a 3D model with fused texture information. The model optimization module is used to overlay the 3D model with fused texture information with the original laser point cloud and panoramic image in the same 3D scene to verify the geometric accuracy, semantic attribute consistency and texture quality of the 3D model, and correct the unqualified parts of the verification results to obtain the 3D model of the target road.

[0007] Another aspect of the present invention provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor; when the computer program is executed by the processor, it implements the steps of the automated digital modeling method for urban roads as described above.

[0008] Another aspect of the present invention provides a computer program product storing a computer program that, when executed by a processor, implements the steps of the automated digital modeling method for urban roads as described above.

[0009] The automated digital modeling method, system, and equipment for urban roads provided in this invention, through deep fusion and high-precision registration of vehicle-mounted laser point clouds and panoramic images, can combine the geometric accuracy of vehicle-mounted point clouds with the visual detail advantages of panoramic images, ensuring the realism and accuracy of the model's geometry and texture. This invention realizes a fully automated construction technology chain for road elements, from identification and extraction to semantic modeling, significantly reducing manual intervention and improving modeling efficiency. The 3D road element model generated by this invention is not merely a geometric object, but carries rich structured semantic attributes, making the final constructed 3D urban road model semantically complete and easy to integrate and apply with business management systems such as BIM / GIS.

[0010] The above description is merely an overview of the technical solution of the present invention. In order to better understand the technical means of the present invention and to implement it in accordance with the contents of the specification, and to make the above and other objects, features and advantages of the present invention more apparent and understandable, specific embodiments of the present invention are described below. Attached Figure Description

[0011] Various other advantages and benefits will become apparent to those skilled in the art upon reading the following detailed description of preferred embodiments. The accompanying drawings are for illustrative purposes only and are not intended to limit the invention. In the drawings: Figure 1 A flowchart illustrating an automated digital modeling method for urban roads provided in an embodiment of the present invention; Figure 2 This is a diagram showing the spatiotemporal registration effect of point cloud and image of the automated digital modeling method for urban roads proposed in this embodiment of the invention. Figure 3 A schematic diagram illustrating the morphological characteristics of different arrow marker lines in the automated digital modeling method for urban roads provided in this embodiment of the invention; Figure 4 A schematic diagram of corner distances for an automated digital modeling method for urban roads provided in an embodiment of the present invention; Figure 5 A curbstone road element structure diagram for an automated digital modeling method for urban roads provided in an embodiment of the present invention; Figure 6 This is a texture mapping effect diagram of the automated digital modeling method for urban roads proposed in an embodiment of the present invention; Figure 7 This is a structural block diagram of the automated digital modeling system for urban roads proposed in an embodiment of the present invention. Detailed Implementation

[0012] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

[0013] Those skilled in the art will understand that, unless specifically stated otherwise, the singular forms “a,” “an,” “the,” and “the” used herein may also include the plural forms. It should be further understood that the term “comprising” as used in this specification means the presence of the stated features, integers, steps, operations, elements, and / or components, but does not exclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and / or groups thereof.

[0014] It will be understood by those skilled in the art that, unless otherwise defined, all terms used herein (including technical and scientific terms) have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. It should also be understood that terms such as those defined in general dictionaries should be understood to have the same meaning as in the context of the prior art and should not be interpreted in an idealized or overly formal sense unless specifically defined.

[0015] Figure 1 A flowchart illustrating an embodiment of the automated digital modeling method for urban roads according to the present invention is shown. (Refer to...) Figure 1The automated digital modeling method for urban roads in this embodiment of the invention specifically includes the following steps: S11. Acquire the original laser point cloud and panoramic image of the target road acquired simultaneously.

[0016] In this embodiment, a vehicle-mounted mobile measurement system (MMS) is used to simultaneously acquire road laser ranging data, panoramic image data, and POS data. The POS data is refined and processed to generate a high-precision trackline file with continuous time, accurate position, and attitude. This file contains the precise latitude, longitude, elevation, roll, pitch, and heading of the acquisition vehicle at each timestamp (laser point cloud and image synchronization). Based on the processed trackline file, the laser ranging data is converted into the original laser point cloud in an absolute coordinate system. Specifically, the sensors on the vehicle-mounted mobile measurement system (MMS) include a laser emitter, a panoramic camera system, an inertial navigation system, a GNSS system, a main control system, and wheel encoders. The MMS simultaneously acquires road laser ranging data, panoramic image data, and high-precision positioning and attitude data (POS data). All sensors are synchronized via a precision clock source (PPS signal), with a synchronization error ≤1ms, ensuring that each frame of laser point cloud and each image is strictly aligned with the corresponding POS pose data in time.

[0017] S12. Extract ground points from the original laser point cloud and perform spatiotemporal registration between the extracted ground point cloud and the panoramic image.

[0018] S13. Extract lane line road elements, road centerline road elements, directional arrow road elements, and curbstone road elements from the ground point cloud after spatiotemporal registration.

[0019] S14. Extract guardrail road features and road ancillary facility road features from the non-ground point cloud in the original laser point cloud, and identify the location and geometric attributes of the road ancillary facility road features.

[0020] S15. Construct lane surface road elements based on the extracted lane line road elements, curb road elements, and guardrail road elements.

[0021] Specifically, based on the extracted lane line vector lines, curb 3D lines, and guardrail 3D lines, a closed 3D vector surface is constructed through topological relationships, and the 3D vector lane surfaces of the main road and the intersection are divided by the intersection.

[0022] S16. Assign the extracted road features to the corresponding vector layers according to their type and configure structured semantic attributes and globally unique codes. Semantic attributes include information related to the road segment, feature type information, location information, and / or geometric parameter information. Specifically, information related to the road segment includes the road segment name code, road segment code, administrative region code, road segment classification code, road segment sequence code, driving direction code, facility classification code, facility sequence code, and road feature code. The globally unique code is the same as the road feature code, and its format is: road segment code - administrative region code - road segment classification code - road segment sequence code - driving direction code - facility classification code - facility sequence code.

[0023] S17. The semantic attributes of road elements are retrieved from the parametric template library to match the model template, and a 3D model of the corresponding road elements is generated based on the corresponding model template.

[0024] S18. Based on the visibility of the surface of the 3D model of each road element in the panoramic image, select the image source that matches the surface of each model, and map the image texture of the image source to the corresponding model surface to obtain a 3D model with fused texture information.

[0025] S19. The 3D model with fused texture information is overlaid with the original laser point cloud and panoramic image in the same 3D scene to verify the geometric accuracy, semantic attribute consistency and texture quality of the 3D model, and the unqualified parts of the verification results are corrected to obtain the 3D model of the target road.

[0026] The automated digital modeling method for urban roads provided in this invention, through deep fusion and high-precision registration of vehicle-mounted laser point clouds and panoramic images, combines the geometric accuracy of vehicle-mounted point clouds with the visual detail advantages of panoramic images, ensuring the realism and accuracy of the model's geometry and texture. This invention realizes a fully automated construction technology chain for road elements, from identification and extraction to semantic modeling, significantly reducing manual intervention and improving modeling efficiency. The 3D road element model generated by this invention is not merely a geometric object, but carries rich structured semantic attributes, making the final constructed 3D urban road model semantically complete and facilitating integration and application with BIM / GIS and other business management systems.

[0027] In this embodiment of the invention, step S12, which involves extracting ground points from the original laser point cloud, includes the following specific steps: S121. Based on the acquisition time of each point in the original laser point cloud, find two adjacent track points for each point in the trackline, and obtain a trackline reference point that is synchronized with the time of each point by linear interpolation of the two adjacent track points. This invention first removes noise from the original laser point cloud data through "statistics + radius" filtering, then performs coarse screening based on the 3D spatial distance and azimuth of the laser point cloud and the flight path, and finally uses a parameter-optimized cloth simulation filtering algorithm to extract ground points in a refined manner.

[0028] Specifically, before determining the reference point for the flight path, the average distance between each point in the original laser point cloud and its k nearest neighbors is first calculated. and its standard deviation Set distance threshold ( = +N* Remove points whose distance exceeds a distance threshold, then count the number of points in the neighborhood of each point within a set radius threshold. When the number of points is lower than the set number threshold, the point is determined to be a small-scale dense noise and removed to obtain the denoised laser point cloud.

[0029] For each point in the laser point cloud, a neighboring track point that satisfies the acquisition time is found in the trackline. A trackline reference point that is strictly time-synchronized with the point is obtained through linear interpolation. Then, the 3D spatial distance and azimuth between the point and the trackline reference point are calculated. Thresholds for the 3D distance and azimuth are set to obtain the laser point cloud of the road surface and the components on both sides. Finally, a cloth simulation filtering algorithm is used to separate the ground points and non-ground points from the extracted laser point cloud.

[0030] Let the original vehicle-mounted laser point cloud be... ,in: For point Three-dimensional coordinates (unit: m); For point Data acquisition time (unit: s); This represents the total number of points in the original laser point cloud. Let the track recorded by the vehicle-mounted POS system be... ,in: for waypoints Timestamp (unit: seconds); for Corresponding IMU center three-dimensional coordinates (unit: m); for The heading angle at the location (0° is due north, and positive is clockwise, unit: °); This represents the total number of waypoints.

[0031] Target Collection time Find adjacent waypoints on waypoint T and (satisfy ≤ ≤ ), obtained through linear interpolation Corresponding track reference point The formula is as follows:

[0032] in, As a time interpolation weight, the reference point is strictly synchronized with the laser point cloud acquisition time.

[0033] S122. Calculate the 3D spatial distance and azimuth angle between each point and the track reference point. Remove points from the original laser point cloud whose 3D spatial distance is greater than a preset spatial distance threshold and whose azimuth angle is greater than a preset azimuth angle threshold to obtain a subset of candidate ground points. The method for calculating 3D spatial distance is as follows: calculate the points in the laser point cloud. To the trackline reference point Euclidean 3D distance, setting a 3D spatial distance threshold Remove laser point clouds that are far from the trackline. The formula for the distance between points in the laser point cloud and the trackline is:

[0034] in, The unit is meters. This is a spatial distance threshold, which can be set according to the road width, such as 10-15m for urban roads and 7-12.5m for highways.

[0035] The azimuth angle is calculated as follows: the azimuth angle is a point... Relative to the direction of travel of the track ( The horizontal angle is calculated by removing laterally irrelevant points on the track. First, the horizontal vector is calculated: Then calculate the midpoint of the laser point cloud. Horizontal azimuth (0° is due north, clockwise is positive): Next, calculate the azimuth difference (take the absolute value to avoid the influence of direction): ,in, The unit is °. This is the azimuth threshold.

[0036] Points from the laser point cloud that satisfy both "3D distance ≤ 3D distance threshold" and "azimuth angle ≤ azimuth angle threshold" are retained to form a candidate ground point subset.

[0037] .

[0038] S123. Use the cloth simulation filtering algorithm to separate ground points and non-ground points from the candidate ground point subset. Specifically, this includes: constructing a two-dimensional particle mesh covering the range of a subset of candidate ground points; driving particle settling under gravity; correcting the updated particle positions based on inter-particle elastic constraints and the particle settling velocity based on damping force during particle position iteration; performing collision detection on each particle and the point cloud in the subset of candidate ground points during particle position iteration; extracting the maximum elevation value from the collision points when collisions exist; determining whether the Z-coordinate value of the current particle is less than or equal to the sum of the maximum elevation value and a preset collision threshold; if the Z-coordinate value of the current particle is less than or equal to the sum of the maximum elevation value and the preset collision threshold, the current particle stops settling and the maximum elevation value is taken as the final elevation value of the current particle; after the particles stop settling, for each point in the subset of candidate ground points, finding the cloth particle closest to the current candidate point; if the difference between the elevation value of the current candidate point and the final elevation value of the corresponding particle is less than or equal to the preset elevation threshold, the current candidate point is determined to be a ground point.

[0039] Specifically, build coverage Two-dimensional particle mesh of range Initial position of each particle: ;

[0040] in: :according to The XOY plane is uniformly divided, and the interparticle spacing is... ( for The average density of the laser point cloud (unit: points / square meter); m is the row index of the grid (from 1 to M'); n is the column index of the grid (from 1 to N'); M' is the total number of rows in the grid; N' is the total number of columns in the grid.

[0041] , (Initial fabric height).

[0042] In this embodiment, particle position updates consider gravity, elastic constraints, and damping forces, and the iteration step is set to kD (kD=50). The initial displacement under gravity (z direction only, x and y directions are fixed).

[0043] , in: The velocity of the particle at step kD-1 (unit: m / s); (Selection of time step); (Acceleration due to gravity).

[0044] In this embodiment, to avoid excessive stretching of the fabric, inter-particle elastic constraints are used for correction: for each particle... Its distance from adjacent particles (up, down, left, and right) is Initial spacing between adjacent particles (Same as particle mesh spacing), current spacing: ; in , For particles Coordinates on the horizontal plane; , represents the horizontal coordinates of adjacent particles. This represents the elevation value of the adjacent particle after the gravitational action at step kD.

[0045] when > (Stretching) fixes x and y coordinates, then corrects the z coordinate: ; In this embodiment, damping force is used to reduce particle velocity and avoid oscillation: ; ; in, The damping coefficient can be selected as 0.15, which is optimized to adapt to the smooth characteristics of the road laser point cloud.

[0046] In this embodiment, the collision detection implementation steps are as follows: For the Step particles ,exist Search horizontal distance ≤ point set :

[0047] Optionally, collision radius It can be optimized to 1.5 times the interparticle spacing, that is... .like (A collision point exists), take the maximum elevation of the collision point:

[0048] If the particle's z-coordinate Then the particles stop settling, and the final z-coordinate is: .

[0049] right Each point in Find the nearest cloth particles If point Elevation difference between the final elevation of the particle If it is, then it is determined to be a ground point: ; in, The elevation threshold can be set to 0.03m. This is an optimized ground detection threshold that adapts to the slight undulations of urban roads.

[0050] In this embodiment of the invention, step S12, which involves spatiotemporal registration of the extracted ground point cloud with the panoramic image, specifically includes: S124. After distortion correction processing, the panoramic image is projected onto a uniform projection surface. The projection surface is then resampled and fused to generate a spherical panoramic image.

[0051] First, panoramic camera calibration is performed. The image plane coordinates of marker points are measured with high precision on the panoramic image, and the absolute coordinates of the corresponding marker points are simultaneously obtained from the laser point cloud. Using the image plane coordinates and absolute coordinates of corresponding image points, a rigorous collinearity equation is established between the image point coordinates, object point coordinates, and the instantaneous pose of the camera. The formula is as follows. The overall adjustment is then performed to solve for the rotation matrix R0 and translation vector T0 between the camera and the IMU, ultimately obtaining the fixed mounting parameters of the panoramic camera relative to the IMU.

[0052]

[0053] in: , , Let these be the coordinates of the point in the camera coordinate system. , These are the coordinates in the IMU coordinate system.

[0054] Then, based on the focal length (fx, fy), principal point (cx, cy), radial distortion coefficient (k1, k2, k3), and tangential distortion coefficient (p1, p2) of each camera, the original fisheye image of each camera is subjected to distortion correction to correct lens distortion. Then, based on the camera's extrinsic parameters (R0, T0), the pixels of each camera image after distortion correction are projected onto a uniform spherical or cubic surface. Pixel resampling is performed on the uniform projection surface, and the images are fused to generate a seamless 360°×180° spherical panoramic image.

[0055] S125. Calculate the color histogram of the overlapping area of ​​adjacent images for the spherical panoramic image, use the histogram matching method to perform color matching, and fuse the overlapping area in the gradient domain to obtain a color-balanced panoramic image.

[0056] In this step, the color histogram of the overlapping area of ​​adjacent images is calculated for the stitched panoramic image, and color matching is performed to initially eliminate color differences; the overlapping areas are then fused in the gradient domain (rather than the pixel domain), using the following formula:

[0057] in: A function to guide the panoramic image stitching and fusion process; It is the panoramic image of the target to be obtained; This is the image after the initial stitching; yes The gradient; It is the expected gradient field (usually derived from) However, a smooth transition is performed at the overlapping boundaries. E It is a parameter that balances the weights of the two items.

[0058] S126. Acquire synchronously collected POS data, construct a joint adjustment model based on POS data, original laser point cloud and camera parameters, and calculate the optimal system parameters using the least squares adjustment method.

[0059] S127. Based on the optimal system parameters, the extracted ground point cloud is mapped onto the panoramic image to obtain the image coordinates of the point cloud in the panoramic image coordinate system.

[0060] In this embodiment, systematic errors between different source data are eliminated by using observations such as POS data, laser point clouds, and camera parameters to obtain joint adjustment optimization system parameters. The optimized parameters are then used in the coordinate transformation chain to batch map massive amounts of laser point clouds onto panoramic imagery. Figure 2 The image shows the spatiotemporal registration effect of point cloud and image of the automated digital modeling method for urban roads proposed in this embodiment of the invention.

[0061] Specifically, POS data, laser scanning angle, camera parameters, and other observations are incorporated into a unified mathematical model to solve for the optimal system parameters and point coordinates. The POS observations ( ), laser point cloud distance observations ( ), image matching point observations ( This is incorporated into the least squares model. The model is shown below: ; ; in: : Residual vector; Unknown parameter vector (including POS system error) Camera intrinsic parameter correction amount Correction amount of lidar installation parameters A: Design matrix (composed of the partial derivatives of each observation with respect to the unknown parameters); L: Observation vector ( (splicing); Weight matrix (POS observation weights) Image matching point weights (where σ is the observation precision).

[0062] Then, based on the optimal system parameters, a series of coordinate transformations are performed to convert the laser coordinates in the lidar coordinate system to... Mapped to the pixel coordinate system of the panoramic image The specific coordinate transformation chain includes: laser scanner coordinate system to vehicle IMU coordinate system, vehicle IMU coordinate system → ECEF geocentric coordinate system, ECEF geocentric coordinate system → camera center coordinate system, camera center coordinate system → normalized coordinate system, normalized planar coordinate system → image coordinate system.

[0063] Specifically, the coordinate system of the laser scanner is converted to the coordinate system of the vehicle-mounted IMU: The relationship between the laser scanner coordinate system and the vehicle-mounted IMU coordinate system is established as shown below.

[0064] ; in: The coordinates of the point in the IMU coordinate system; The corrected rotation matrix of the lidar relative to the IMU is determined by the installation angle. Obtained through Euler angle transformation . This refers to the original mounting roll angle of the lidar relative to the IMU; The original installation elevation angle of the lidar relative to the IMU; This is the original mounting heading angle of the lidar relative to the IMU; This is the correction amount for the roll angle; This is the correction amount for the pitch angle; This is the correction amount for the heading angle.

[0065] The single-axis rotation matrix (Roll around X-axis, Pitch around Y-axis, Heading around Z-axis) is: ; ; ; , The unit is radians.

[0066] : The corrected translation vector of the lidar relative to the IMU. , . This represents the original installation translation (X direction) of the lidar relative to the IMU. This represents the original installation translation (Y direction) of the lidar relative to the IMU. This represents the original installation translation (Z direction) of the lidar relative to the IMU. This represents the corrected translation (X-direction) of the lidar relative to the IMU. This represents the corrected translation (Y direction) of the lidar relative to the IMU. This represents the corrected translation (Z direction) of the lidar relative to the IMU.

[0067] It is a fundamental matrix describing rotation in three-dimensional space, corresponding to rotations around the Z-axis, Y-axis, and X-axis, respectively, and is used to describe the attitude rotation relationship between the lidar coordinate system and the vehicle-mounted IMU coordinate system.

[0068] α, β, and γ are the rotation angles corresponding to these three rotation axes.

[0069] Specifically, the coordinate system of the vehicle-mounted IMU changes to the geocentric coordinate system of ECEF: The relationship between the vehicle-mounted IMU coordinate system and the ECEF geocentric coordinate system is established as follows.

[0070] ; in: Coordinates of a point in the ECEF coordinate system (unit: m) The corrected rotation matrix of the IMU relative to ECEF, derived from the POS attitude angle. )calculate; The raw roll angle measured by the IMU; The raw pitch angle measured by the IMU; The original heading angle measured by the IMU; This is the correction amount for the roll angle; This is the correction amount for the pitch angle; This is the correction amount for the heading angle.

[0071] The IMU is unknown after correction in the ECEF coordinate system, and is originally unknown from the POS coordinate system. Add correction amount to obtain . The original X coordinates of the IMU in the ECEF coordinate system; The original Y coordinate of the IMU in the ECEF coordinate system; The original Z coordinate of the IMU in the ECEF coordinate system; This is the correction amount for the X coordinate; This is the correction amount for the Y-coordinate; This is the correction amount for the Z-coordinate.

[0072] Specifically, ECEF geocentric coordinate system → camera center coordinate system: By using the corrected camera placement parameters, the relationship between the ECEF geocentric coordinate system and the camera center coordinate system is established, as shown below.

[0073] ; in: : The coordinates of the point in the camera center coordinate system (the camera center is the origin, the Z-axis is forward, the X-axis is right, and the Y-axis is downward); The rotation matrix of ECEF relative to the camera is obtained by combining the camera's mounting angle relative to the IMU and the rotation matrix of the IMU relative to ECEF. ,in This is the corrected rotation matrix of the camera relative to the IMU, therefore .

[0074] The position of the camera center in the ECEF coordinate system is obtained by combining the IMU position and the translation vector of the camera relative to the IMU. ,in This is the original translation vector of the camera relative to the IMU. This is the translation correction amount.

[0075] Specifically, camera center coordinate system → normalized coordinate system: Based on the pinhole camera model, projecting 3D points onto a 2D normalized plane (without distortion), the relationship between the camera center coordinate system and the normalized coordinate system is established using the following formula: ; in: , : The coordinates of the point in the normalized plane coordinate system (unit: mm); The corrected camera focal length is obtained by adding an intrinsic parameter correction amount to the original focal length. , To calibrate the camera's original focal length, This is the focal length correction amount. If it is a multi-camera system, the focal length of the corresponding camera needs to be taken.

[0076] The negative sign means that the camera coordinate system's Z-axis is forward, while the normalized plane is behind the camera, therefore the X and Y directions need to be reversed.

[0077] Specifically, normalized planar coordinate system → image coordinate system: The distortion of the camera lens is corrected and intrinsic parameters are mapped to establish the relationship between the corrected normalized coordinate system and the image coordinate system.

[0078] 1) Through the corrected distortion coefficients The focal length (fx, fy), principal point (cx, cy), radial distortion coefficients (k1, k2, k3), and tangential distortion coefficients (p1, p2) of each camera are used to correct for distortion in the camera lens. The formula is as follows: First, calculate the radial distance in the normalized planar coordinates: ; Then radial distortion correction:

[0079]

[0080] in The radial distortion coefficient is corrected; then the tangential distortion is corrected.

[0081] ; ; in This is the corrected tangential distortion coefficient.

[0082] 2) Intrinsic parameter mapping The corrected normalized coordinates are mapped to image coordinates using the corrected camera intrinsic parameters, using the following formula: ; ; in: u,v: Pixel coordinates of the point on the panoramic image (u is the column number, v is the row number, both are integers); Corrected pixel focal length (dx, dy are camera pixel sizes, unit: mm / pixel); Corrected principal point pixel coordinates.

[0083] Furthermore, after obtaining the image coordinates of the point cloud in the panoramic image coordinate system, the method also includes registration accuracy verification. Specifically, in the panoramic image and the laser point cloud, clearly corresponding feature points (such as road marking corner points, sign corner points, etc.) are selected as independent checkpoints. The pixel coordinates of the selected independent checkpoints in the panoramic image are ( The corresponding points in the laser point cloud are projected to the image pixel coordinates through a coordinate transformation chain to obtain the calculated values. The pixel error at each checkpoint is then... , The Euclidean distance of the pixel position error is And calculate the average pixel error for all checkpoints. and root mean square error .

[0084] The laser point cloud coordinates of the independent checkpoint are compared with the three-dimensional coordinates obtained through image backprojection. Let the three-dimensional coordinates of the checkpoint in the point cloud be... The three-dimensional coordinates obtained through image back projection (using camera parameters and POS data) are: Then the three-dimensional spatial error is , , The Euclidean distance of the three-dimensional error is Then calculate the average three-dimensional error of all laser point clouds corresponding to the independent checkpoints. and root mean square error Based on the application requirements of digital city road modeling, the average pixel error... Less than or equal to 1.5 pixels and average 3D spatial error The accuracy must be within 0.05 meters. If the accuracy does not meet the requirements, the joint adjustment parameters or data quality must be rechecked, and iterative optimization must be performed until the accuracy standard is met.

[0085] In one embodiment of the present invention, step S13, which involves extracting lane line road features, road centerline road features, directional arrow road features, and curbstone road features from the spatiotemporally registered ground point cloud, specifically includes the following steps: S131. Construct a sparse triangular network based on the ground point cloud, and optimize the sparse triangular network to form a road reference surface. Specifically, this includes: selecting a very sparse initial seed point set from the ground point cloud, and constructing a sparse triangular network based on the initial seed point set; setting a point cloud interval threshold Li and a height threshold Hi, which can be selected as Li=50cm (Hi=0.15cm); traversing every ground point in the ground point cloud except for the sparse triangular network, obtaining the distance and height difference between the ground point and the current sparse triangular network, as well as the distance between the ground point and the three vertices of the current sparse triangular network. When the distance to the three vertices is greater than or equal to the preset interval threshold, or the height difference is greater than or equal to the height difference threshold, the ground point is added as a feature point to the current sparse triangular network to obtain an updated sparse triangular network. Otherwise, the point is considered a "redundant point" and will be ignored and discarded. Repeat the current step to update the sparse triangular network until no ground point meets the conditions, and finally form a road reference surface based on the feature points.

[0086] S132. Assign color features to the point cloud data in the road reference plane based on the panoramic image to obtain a reconstructed 3D point cloud with color features.

[0087] In this step, color-balanced panoramic images are used to assign color feature values ​​to the reconstructed point cloud, so that the reconstructed 3D point cloud has accurate color features.

[0088] S133. Extract lane line vectors, road centerline vectors, and directional arrow vectors from reconstructed 3D point clouds; S134. Extract 3D lines of curb stones based on ground point clouds.

[0089] In this embodiment, before assigning color feature values ​​to the point cloud data in the road reference plane based on the panoramic image, the method further includes the steps of point cloud intensity correction and multimodal super-resolution reconstruction. Specifically, this includes: enhancing and correcting the intensity values ​​of the point cloud in the road reference plane to obtain the effective intensity values ​​of the corresponding point cloud; generating a depth map of the road reference plane based on the intensity information of the enhanced and corrected point cloud; and performing an inverse mapping from the image pixel coordinate system to the point cloud coordinate system on the depth map to achieve point cloud reconstruction. Further, assigning color feature values ​​to the point cloud data in the road reference plane based on the panoramic image specifically involves: assigning color feature values ​​to the reconstructed point cloud of the road reference plane based on the panoramic image to obtain a reconstructed 3D point cloud with both color and grayscale features.

[0090] In this embodiment, the intensity value of the point cloud in the road reference plane is enhanced and corrected, specifically including: Calculate the perpendicular distance from each point cloud in the road reference plane to the data acquisition trajectory line. Based on the extracted road reference plane, using point cloud intensity information, employ a "distance constraint + intensity correction" method. For coarse ground point clouds, strictly synchronize the vehicle-mounted LiDAR point cloud P with the trajectory line T based on GNSS timestamps, and calculate the perpendicular distance from each point cloud Pi to the trajectory line T. .

[0091] The original intensity values ​​of the point cloud within a preset vertical threshold range are statistically analyzed. A preset distance attenuation correction model is used to perform nonlinear mapping on the original intensity values ​​to enhance the point cloud intensity, significantly improving the reflection contrast between the road markings and the road surface. The distance attenuation correction model is as follows: ; in, The vertical distance from point cloud Pi in the road reference plane to the data acquisition trajectory line is denoted as . Vertical distance The corrected effective intensity value of the point cloud. Vertical distance The original intensity value of the point cloud. This is a standard distance, a reference distance set for normalizing the intensity values ​​of the laser point cloud. The standard distance for scanning urban roads is 20-30 meters. This is the preset atmospheric attenuation coefficient.

[0092] Preset vertical threshold range D G The value can be 5-10 meters. The above operation can be used to exclude point cloud data with strong near-end contrast and areas with insufficient point cloud density at the far end.

[0093] In this embodiment, a depth map of the road reference surface is generated based on the intensity information of the enhanced and corrected point cloud, specifically including: Step 1: Grayscale DOM (Digital Orthophoto Map) Creation. A sparse depth map is generated by selecting point cloud data from the enhanced and corrected point cloud data that are within a preset height threshold range from the ground. Specifically, this includes: selecting point cloud data from the enhanced and corrected point cloud data that are within a preset height threshold range from the ground to obtain target point cloud data; determining the coverage area of ​​the digital orthophoto map to be generated based on the two-dimensional planar range of the road reference plane, and creating a grayscale image grid based on the coverage area and a preset ground resolution, with each pixel in the grayscale image grid having an initial value of 0 or NULL; vertically projecting the target point cloud data onto the grayscale image grid, and determining the effective pixel value of each grid in the grayscale image grid based on the effective intensity value of the target point cloud data, resulting in an image matrix; and linearly contrast-stretching the grayscale values ​​of the image matrix to convert them to an integer range of 0-255 to obtain the sparse depth map.

[0094] The specific implementation method is as follows: Determine the range and resolution: Based on the two-dimensional planar range (X,Y) of the classified ground point cloud and the range of 1 meter above the ground points, obtain the point cloud for generating the grayscale DOM, and determine the coverage range and ground resolution of the grayscale DOM.

[0095] Create a grid: Create an empty grayscale image grid based on the parameters above, with each pixel initialized to 0 or NULL.

[0096] Point cloud projection and assignment: The ground point cloud is vertically projected onto an XOY plane grid. For each grid pixel, the maximum or average value of all point cloud intensity values ​​within that pixel is assigned to that pixel.

[0097] When grid pixels (x) P ,y P (This refers to a function that retrieves the maximum intensity value of a point cloud within a grid.) .

[0098] When grid pixels (x) P ,y PIf there are no point clouds (invalid pixels) within a certain area, search for the K0 nearest valid pixel values ​​within that area, calculate their interpolation weights, and assign them. The formula is:

[0099] in, The intensity values ​​are the K0 nearest neighboring valid pixels, optionally K0=4. Let be the distance between pixel (x, y) and the j-th valid pixel.

[0100] Intensity normalization: The gray values ​​of the interpolated image matrix are linearly contrast stretched to convert them to the integer range of 0-255, resulting in a sparse depth map D1.

[0101] when hour,

[0102] when When the intensity is uniform across the entire region, (in grayscale).

[0103] The normalized pixel grayscale value, ranging from 0 to 255, is the pixel value ultimately used to generate the sparse depth map. It is the minimum gray value of all pixels in the original image matrix. This represents the maximum grayscale value of all pixels in the original image matrix.

[0104] Step 2: Upsample the sparse depth map using bilinear interpolation based on the panoramic image resolution to obtain an initial depth map with the same resolution as the panoramic image. Specifically, upsample D1 using bilinear interpolation to make its resolution consistent with the panoramic image resolution, resulting in the initial depth map D2; for the hole pixel p in D2, select non-empty depth pixels with grayscale similarity within a 5×5 neighborhood, and take the average value as the pre-fill value.

[0105] Step 3: Perform pixel-by-pixel local structure analysis on the panoramic image (i.e., the color-balanced panoramic image) to generate filtering parameters to guide subsequent filtering. Then, based on the panoramic image, apply the generated filtering parameters to the initial depth map for adaptive joint bilateral filtering to obtain a dense depth map. The pixel-by-pixel local structure analysis to generate filtering parameters includes: calculating the structure tensor of the panoramic image in the neighborhood of each pixel; determining the local structure type of the neighborhood centered on the current pixel based on the feature values ​​of the structure tensor (local structure types include flat regions, edge regions, or textured regions); and calculating the kernel function parameters for the corresponding pixel in real time based on the local structure type of each pixel.

[0106] In a specific example, the structure tensor of the panoramic image within the neighborhood of each pixel p0 is calculated, and the region is classified as a flat region, an edge region, or a textured region based on the eigenvalue relationship. The formula is as follows:

[0107] in: , These are the Sobel gradients of pixel q in the x and y directions, respectively; right Eigenvalue decomposition yields two eigenvalues. Determine the region type based on the eigenvalue relationship: when When the gradient is small and there is no obvious structure, the region type is a flat region, and a standard Gaussian kernel is selected; when When the gradient is large in one direction and there is a clear edge, the region type is an edge region, and the Cauchy kernel with a heavy-tailed distribution is selected as the range kernel; when When there are large gradients in multiple directions and complex textures, the region type is a textured region, and an anisotropic Gaussian kernel is selected.

[0108] Based on local statistical properties, the kernel function parameters (spatial parameters) of each pixel p are calculated in real time. , range parameter ).

[0109] ;

[0110] ;

[0111] in: This represents the proportion of non-empty pixels in the neighborhood to the total number of pixels. , To preset basic parameters, It represents the variance of the image intensity (color) within the neighborhood of pixel P, reflecting the texture complexity of that region.

[0112] Based on the panoramic image I, adaptive joint bilateral filtering is performed on D2 to obtain the dense depth map D3. The formula is:

[0113]

[0114] in: To output the depth value of the dense depth map at pixel p; Normalization factor; It represents all pixels q within a neighborhood window Ω centered at pixel p; Let be the spatial Euclidean distance between pixel P and its neighboring pixel q; The standard deviation parameter of the spatial kernel (in pixels) controls the decay rate of the spatial weights; This represents the depth value of the initial depth map D2 at pixel q. The range kernel standard deviation controls the decay rate of color weights.

[0115] It is a Gaussian function.

[0116] In this embodiment of the invention, the implementation of lane line vector line extraction based on reconstructed 3D point cloud in step S133 specifically includes lane feature seed point extraction and lane line fitting, wherein: Lane feature seed point extraction: The Otsu's method is used to extract lane line feature seed points from the reconstructed 3D point cloud. Specifically, this includes: traversing the reconstructed 3D point cloud and assigning point cloud intensity... Points with intensity values ​​greater than a preset global intensity threshold are identified as candidate lane line points. (Global intensity threshold) To maximize the value of the intensity threshold T of the objective function, that is, the variance between the background (road surface) and the target (marking) classes. The objective function for finding the maximum value of T is:

[0117] The variance between the point clouds of road surface and lane lines is given. These represent the proportions of point cloud numbers below and above the threshold T, respectively. These represent the average intensity of point clouds below and above the threshold T, respectively. The candidate lane line points are divided into several local grids based on their spatial location. For each grid, the intra-class variance of the intensity values ​​is calculated. For grids with intra-class variance exceeding a set variance threshold, the points within those grids are sorted by intensity value to remove points whose intensity values ​​deviate from the mean. The remaining candidate lane line points are used as lane line feature seed points. Specifically, based on the candidate lane line points, several local grids are divided according to their spatial location. The grid size is adaptively set according to the point cloud density, typically 0.5m × 0.5m. For the point cloud within each grid, the intra-class variance of its intensity values ​​is calculated. The formula is as follows:

[0118] Where: NI is the number of points in the grid. This represents the average intensity within the grid. When the intensity of a certain grid... If the intensity value exceeds a set threshold (indicating a discrete intensity distribution and potential noise), then the points within the grid are sorted by intensity, and those with intensity values ​​deviating from the mean by more than a certain amount are removed. * point, You can choose 2-3. The standard deviation of the intensity within the grid; the points ultimately retained are the lane line feature seed points, forming the refined seed point set. .

[0119] Lane line fitting: Based on a preset index distance, continuous points of lane line feature seed points that satisfy azimuth constraints are searched, and vectorized lane lines are fitted based on the searched continuous points. Specifically, this includes: obtaining a set of candidate feature seed points adjacent to each lane line feature seed point based on the index distance; calculating the geometric center point of each candidate feature seed point set and using it as the current seed point; obtaining the azimuth angle of the track line where each geometric center point is located; using each geometric center point as the origin and the corresponding azimuth angle as the initial direction, continuing to search for the next level of candidate feature seed point set according to the index distance; calculating its geometric center point and using it as the subsequent candidate seed point of the current seed point; calculating the azimuth angle of the line connecting the current seed point to the previous seed point and the subsequent candidate seed point, respectively, and calculating the azimuth angle change of the azimuth angle of the line connecting the current seed point to the previous seed point and the subsequent candidate seed point; if the azimuth angle change is greater than a preset azimuth angle change threshold, the subsequent candidate seed point is removed, and the current index distance is shortened. The search continues until the next set of candidate feature seed points is searched again until the preset termination conditions are met. The termination conditions include: if the updated index distance is less than the preset iteration distance threshold (optional 0.3 meters), or the number of iterations exceeds the preset iteration number threshold (optional 5 times), or the number of candidate points is insufficient to form a valid geometric relationship, or the azimuth angle and azimuth angle change of the line connecting the current seed point to the previous seed point and the subsequent candidate seed point do not meet the requirements, the search is terminated; if the azimuth angle change is less than or equal to the preset azimuth angle change threshold, the subsequent candidate seed point is confirmed as a seed point and the search continues until the next seed point cannot be found; the finally obtained ordered seed point sequence is fitted into a smooth vector curve using the least squares method to obtain the lane line vector line.

[0120] Specifically, based on the lane line feature seed point set extracted above... Set the index distance Dseed, which can be 2 meters, and obtain the set of nearby next-level candidate feature seed points according to the index distance. The density clustering algorithm was used to obtain the candidate clusters for the marking line. Each set of feature seed points corresponds to a candidate set of feature seed points, i.e., a candidate cluster for the datum line. For each cluster... Calculate its geometric center point .

[0121] ; Obtain the azimuth of the current flight path. From any center point Departure, based on its current position and azimuth. Using the initial direction, search for the next set of feature seed points according to the index distance. And calculate its geometric center point. ; Calculate the azimuth angle of the seed point and the line connecting the previous and next seed points. The formula is: ; Set the azimuth angle abrupt change threshold (Usually 15°), calculate the change in azimuth angle of the seed point before and after. ,when Greater than When the index is empty, remove the next seed point information, shorten the current index distance, and re-search the feature seed point set; when Less than When a point is found, it is accepted, and the search continues until no further seed point can be extracted. The final ordered sequence of center points is then obtained. The least squares method is used to fit a smooth vector curve. Finally, the fitted lane lines are overlaid with the grayscale DOM of the point cloud, the original vehicle point cloud, and the panoramic image, and local discontinuities or deviations are corrected. The formula is: ;in: For p-th degree B-spline basis functions; These are control points.

[0122] In this embodiment of the invention, the implementation of extracting the road centerline vector line based on the reconstructed 3D point cloud in step S133 specifically includes: Bidirectional inner lane line extraction and matching: Identifying bidirectional inner lane lines based on the type attributes and location distribution characteristics of lane line vector lines. Specifically, this includes: obtaining the vectorized lane line set LC, and identifying and extracting two inner lane lines representing bidirectional traffic flow based on the lane line type attributes (long markings) and spatial distribution characteristics. and .

[0123] Two-way lane line resampling: Resampling of one inner lane line yields an ordered set of dense sampling points. The unit normal vector at each sampling point is calculated. A perpendicular line is constructed passing through each sampling point and along the direction of the corresponding unit normal vector. The geometric intersection of this perpendicular line with the other inner lane line and the corresponding sampling point are then used to construct a perpendicular pair. Specifically, this includes: selecting one inner lane line... Resampling is performed at fixed arc length intervals ΔS to obtain an ordered set of dense sampling points. traversal adjacent points ( Calculate the tangent direction vector at that point. Normalize the tangent vector to obtain the unit tangent vector. Rotating the unit tangent vector by 90 degrees will give you the unit normal vector of the perpendicular line at that point. ,in, and It is the unit tangent vector The two components. Passing through the reference point. Along its unit normal vector Construct a perpendicular line in the direction of the direction. The formula is: Where t is a parameter representing the distance moved along the normal direction. (Obtain) The perpendicular line and the other lane line geometric intersection That is, relative to the reference point Points that have a strict geometric perpendicular correspondence.

[0124] Centerline point calculation and generation: The midpoint of the line connecting each pair of perpendicularly corresponding points is taken as the road centerline point of the current road surface. All the obtained road centerlines are arranged sequentially to form the road centerline point set. Specifically, this includes: for each pair of matched sampling points ( , ), calculate the midpoint of the line connecting the two points. As the center point of the road at this cross-section:

[0125] Arrange all the center points of the cross sections in order of station number to form the road centerline point set: .

[0126] Centerline smoothing and optimization: Based on the road centerline point set, B-spline curve fitting is used to smooth the road centerline vector line and assign semantic attributes. Specifically, this includes: using a B-spline smoothing curve fitting algorithm on the centerline point set... A fitting process is performed. The formula is: ;

[0127] in: The basis functions are cubic B-spline functions. For control points, from the original centerline point set The unknown parameters are obtained by inverse calculation; The parameters of the curve are ([0-1]). The road centerline C obtained by fitting is entered into the road coding attribute, and the start and end points of the centerline are determined. If the direction of the start and end points is opposite to the station number order, the direction of the centerline is reversed.

[0128] In this embodiment of the invention, the implementation of extracting directional arrow vector files based on reconstructed 3D point clouds in step S133 specifically includes: Initial screening of arrow candidate point clouds: Lane line vector lines are buffered with a preset width to form a buffer region covering the lane lines; point clouds within the buffer region are removed from the reconstructed 3D point cloud, and the remaining point clouds are vertically projected onto the XOY plane to generate a binary digital orthophoto image, where the pixel value of the marker line region is 1, and the pixel value of the road background region is 0; arrow candidate regions are selected based on the geometric attributes of each marker line in the binary image. Specifically, this includes: a refined set of lane line feature seed points. Based on the reconstructed point cloud using intensity correction and multimodal super-resolution, the fitted lane line vector curve is processed with a width... The buffer, optionally set to 0.2-0.3 meters, generates a continuous buffer area covering all lane lines. Points within the buffer are removed from the multimodal super-resolution point cloud and vertically projected onto the XOY plane to generate a binarized digital orthophoto image, where the arrow region pixel value is 1 and the road surface background pixel value is 0. Morphological opening operations are performed on the binary image to eliminate fine noise; morphological closing operations are performed to fill small holes and breaks inside the arrows. Connectivity component analysis is then performed on the closed image to obtain the geometric properties of each region, such as area A and the length L of the minimum bounding rectangle. R and width W R Elongation E L =L R / W R Set the threshold for the area range of the arrow. } and elongation threshold , when A AREA satisfy{ And E L Less than If the condition is met, the region is determined to be a candidate region for an arrow; otherwise, the region is determined not to be a candidate region for an arrow.

[0129] Arrow category selection: Traverse the candidate arrow region to extract the smallest bounding rectangle of the marker line. Perform preliminary arrow category classification for each marker line based on the number of marker line corner points, the length of the bounding rectangle, the distance between the widest point of the marker line and the bottom edge of the bounding rectangle, the positional relationship between the starting point of the marker line and the midpoint of the bottom edge of the bounding rectangle, and / or the area of ​​the bounding rectangle. Match each marker line with the binary image of each arrow template in the preset arrow template library, and select the most relevant arrow template as the best template for the corresponding marker line.

[0130] In a specific example, the candidate arrow regions are traversed, and the smallest bounding rectangle of the arrow is taken as a unit. The various types of marker lines are initially classified based on constraints such as the length of the bounding rectangle, the distance between the widest part of the marker line and the bottom edge of the rectangle, the positional relationship between the starting point of the marker line and the midpoint of the bottom edge of the rectangle, and the area of ​​the rectangle. Then, all arrows are oriented and rotated to a uniform horizontal direction. Finally, they are matched with binary images in the arrow template library, and the correlation coefficient is calculated to obtain the best template. Figure 3 This is a schematic diagram showing the morphological characteristics of different arrow markers; where L is the length of the circumscribed rectangle, D is the distance between the widest point of the marker and the bottom edge of the rectangle, Q is the starting point of the marker, and Z is the midpoint of the bottom edge of the rectangle.

[0131] For each candidate arrow region, first extract the smallest circumscribed rectangle contour of the arrow, obtain the coordinates of the four vertices and the center coordinates of the rectangle contour, and use Hough lines to detect the longer side of the circumscribed rectangle. (Formula) ,in Let be the slope, where It is the ordinate value of any point on the line (usually the row coordinate in the image coordinate system). It is the x-coordinate of any point on the line (usually the column coordinate in an image coordinate system). Let be the slope of the line. It is the intercept of the line.

[0132] Using the center point of the rectangle as the origin, calculate the rotation angle required to rotate the longer side to horizontal. Construct the rotation matrix M. R Rotation transformation is performed on the coordinates of all pixels within the region to obtain the orientation-corrected binary image R. R This is done to uniformly rotate arrows from different directions to a standard horizontal orientation, ensuring accurate matching and recognition with the standard template library. The formula is shown below: ; ; Rotation angle (unit: radians). The angle required to rotate the longer side of the rectangle circumscribed by the arrow to the horizontal direction. D : The original X coordinate of the point to be rotated (e.g., a pixel within the arrow area in the image). y D Ox: The original Y-coordinate of the point to be rotated. Ox: The X-coordinate of the rotation center, typically the X-coordinate of the center point of the smallest bounding rectangle of the arrow candidate region. Oy: The Y-coordinate of the rotation center, typically the Y-coordinate of the center point of the smallest bounding rectangle of the arrow candidate region. x′: The new X-coordinate after rotation. y′: The new Y-coordinate after rotation.

[0133] Then, the binary image R is extracted. RThe minimum bounding rectangle of the 255 intensity region is used as a seed region to be identified. The Harris corner detection algorithm is used to detect corners in the seed region to be identified, and the number of corners N is counted. R The number of diagonal points determines the matching template for the seed region to be identified. N R ≤10, call the straight arrow and turn arrow templates; when N R >10, call the straight-ahead turn and straight-ahead bidirectional turn arrow templates.

[0134] Arrow vector file extraction: The vectorized data of the best template is retrieved from the arrow template library and overlaid on the corresponding marker line positions to form a directional arrow vector file. In a specific example, the arrow template is moved to the seed area to be identified until the center point of the rectangle, the midpoints of the top and bottom edges of the rectangle coincide, as shown below. Figure 4 As shown in the diagram, the corner distance is illustrated. Block C3 represents a corner point of the template library identifier line, recording its horizontal and vertical coordinates. Block A2 represents the corner point corresponding to the marker line in the binary image, recording the horizontal and vertical coordinates. First, calculate the Euclidean distance between each pair of corresponding corner points. Then calculate the average distance of all corresponding corner points. Finally, calculate the standard deviation of the distance between the corner points of the two figures. , Figure 4 A schematic diagram of the corner distance for the automated digital modeling method for urban roads provided in this embodiment of the invention is shown below, with the formula as follows: ; ; ; ; Set the standard deviation threshold When the standard deviation < When the matching of the marker line is successful, vectorized data is retrieved from the template library and applied to the position of the marker line to form a vector file.

[0135] In this embodiment of the invention, the implementation of the three-dimensional line extraction of curb stones based on ground point clouds in step S134 specifically includes: S1341. Vertically project the ground point cloud onto the XOY plane, divide the projected image into a grid, and filter out points at a preset curb height threshold based on the point cloud elevation. The grid corresponding to the point cloud region within the specified range is used as the candidate grid region for curb stones. This invention uses a height difference threshold ΔZthreshold to initially filter grid regions that may contain curb stones from ground points. The calculation formula is as follows: ,in These represent the average elevations of adjacent grid cells.

[0136] S1342. Cut cross sections perpendicular to the road direction at preset first cutting distance intervals along the road line direction, and identify candidate curbstone points from the points in each cross section based on elevation change information.

[0137] In this embodiment, along the road line direction, the point cloud of the candidate grid region for the curbstone is spaced at intervals (D). L Cut a cross-section perpendicular to the road direction (5 meters). Sort the points within each cross-section according to their distance from the centerline of the cross-section (or their X-coordinate). For the sorted point sequence, calculate the distances between adjacent points. The ratio of elevation difference to horizontal distance difference ,when The absolute value is greater than the set threshold. ,but These are considered candidate curb points. Abrupt changes in positive and negative gradients will occur at the curb point. The formula is:

[0138] Then, select three consecutive points from the candidate curbstone points. Points can be calculated The approximate curvature at a given point is used to obtain the point where the curvature is maximized (corresponding to a corner or edge point). The formula is:

[0139] Set the height difference between the curb point and its adjacent road point. When the point of maximum curvature has an elevation difference H with the road surface Z Greater than or less than At that time, noise points were eliminated.

[0140]

[0141] in: Elevation of the curbstone; This refers to the road surface elevation.

[0142] S1343. Cluster the candidate curbstone points identified from each cross section.

[0143] In this embodiment, calculations are performed from the point cloud of candidate curb stones in each cross-section. distance , Less than the threshold The points are grouped into the same category. The formula is: ;

[0144] Lower limit: Must be greater than the maximum spacing of the point cloud. (For example, if the average point cloud density is 0.05m, then...) It should be set to at least around 0.1m. Upper limit: Must be less than the minimum distance between different curbstone facets.

[0145] S1344. Randomly select several points from the curbstone cluster points and fit a straight line model. Calculate the distance from other points in the curbstone cluster points to the straight line model. If the distance is less than a preset point determination threshold, determine that the current point is an interior point of the straight line model; otherwise, it is an exterior point. Count the number of interior points in the current straight line model and compare the number of interior points with a preset interior point number threshold. When the number of interior points is greater than the interior point number threshold, update the interior point number threshold to the current number of interior points and record the current straight line model and the corresponding interior point set. Repeat this step to re-obtain the straight line model and count the number of interior points in the current straight line model according to the updated interior point number threshold. Finally, select the straight line model with the most interior points to fit the three-dimensional curbstone line.

[0146] In this embodiment, a small number of points are randomly selected from the curbstone cluster points to obtain a straight line model. Then, the distance from all other points in the curbstone cluster points to the newly obtained straight line model d is calculated. L The distance.

[0147]

[0148] B1 and B1 are coefficients that determine the direction of the line; C1 is a constant term that is related to the position (intercept) of the line; X6 and y6 are the planar coordinates of the point (the coordinates of the curbstone cluster points projected onto the XY plane). When d L Less than a set distance threshold D L If the value is true, then the point is an interior point of the current model; otherwise, it is an exterior point. Count the number of interior points I. L .

[0149] Current model I L The number of interior points and the number of interior points in the current model Compare. If I L > Then update =I L The current model and its corresponding set of interior points are recorded. The model and interior points are then iteratively acquired again following the steps above. Finally, the model with the most interior points is selected to fit the 3D line of the curb. The formula is: , : Direction vector (principal eigenvector); : Centroid coordinates.

[0150] In one embodiment of the present invention, step S14, which involves extracting guardrail road elements from non-ground point clouds in the original laser point cloud, specifically includes: selecting point clouds from the non-ground point cloud whose vertical distance from the road reference plane is within a preset guardrail height threshold range to obtain a candidate guardrail point cloud set; cutting cross sections perpendicular to the road direction along the road line direction at preset second cutting distance intervals; extracting candidate guardrail point clouds in the buffer zone of each cross section; sorting the point clouds in each cross section according to their distance from the center line of the cross section; dividing the point clouds in the cross section into different guardrail point cloud clusters according to the time sequence of the flight path; selecting the highest point of each guardrail point cloud cluster according to the elevation information; and fitting and generating a three-dimensional guardrail line based on the highest point of each guardrail point cloud cluster.

[0151] In this embodiment, based on the extracted road reference plane, non-ground points are first screened. The vertical distance from the point cloud to the road reference plane is calculated, and point clouds within a typical guardrail height range are obtained according to a distance threshold (1-1.5 meters). Upper and lower limits for intensity are set. Highly reflective objects (such as traffic signs) and low-reflective objects (such as dark vegetation) are removed, and then candidate guardrail crossbar points are obtained by combining manual inspection and modification.

[0152] Then, along the road trajectory, cross-sections perpendicular to the trajectory line are generated at fixed intervals ΔS (e.g., 0.5m). For each cross-section, candidate guardrail point clouds are extracted from its buffer zone. The point clouds in the current cross-section are sorted according to their distance d from the cross-section centerline. Based on the time sequence of the flight track, the point clouds in the profile are divided into different guardrail point cloud clusters, each cluster representing the point cloud set of the guardrail profile. The highest point is selected according to the elevation information of the point cloud set, and the highest point of each profile is fitted to generate a 3D line. The formula is as follows: ; ; ; Parameter values ​​(define the "order" and "position" of points on the curve, thus mapping scattered three-dimensional points onto a one-dimensional, ordered parametric axis); , , for The corresponding three-dimensional coordinate values ​​represent the continuous path of the guardrail's three-dimensional lines. , ... The polynomial fitting coefficients are the x-coordinates. , ... The polynomial fitting coefficients are for the y-coordinate. , ... The coefficients are the polynomial fitting coefficients for the z-coordinate.

[0153] The extracted 3D vector lines of the guardrail (including horizontal bars and posts) are overlaid with the extraction results of S2.2 and S2.3, point cloud intensity DOM, panoramic images, etc., and a human-computer interaction interface is provided for operators to check, edit and confirm. Areas that are difficult for the algorithm to handle, such as occlusion and severe damage, are manually measured or corrected.

[0154] In one embodiment of the present invention, step S14, which involves extracting road features of road ancillary facilities from the non-ground point cloud in the original laser point cloud and identifying the positional and geometric attributes of these road features, specifically includes: selecting cross-sectional point sets from the non-ground point cloud whose vertical distance from the road reference plane is within a preset road facility height threshold range to obtain candidate road facility point clouds; searching for point sets exhibiting an arc-shaped distribution within the cross-sectional point set and clustering adjacent isolated arc point sets to form a candidate set of rod-shaped column parts; performing region growing on the candidate set of rod-shaped column parts based on a grid index to search for target point clouds of rod-shaped facilities; fitting rods based on the target point clouds of rod-shaped facilities to identify the positional, geometric, and azimuth attributes of the target rod-shaped facilities and configuring them into the attribute information; and for traffic signs and non-rod-shaped facilities, extracting target point cloud regions from the candidate road facility point cloud according to user instructions, identifying the spatial positional attributes of traffic signs and non-rod-shaped facilities, and extracting geometric and azimuth attributes.

[0155] Specifically, based on the road reference plane, an elevation threshold (1-1.2 meters) is set, and cross-sectional point sets within the threshold range are extracted from non-ground points. Then, points exhibiting an arc-shaped distribution are searched within the cross-sectional point set, and adjacent isolated arc-shaped point sets are clustered to form candidate sets for the rod-shaped column portion. Furthermore, based on the statistical regularities of the three-dimensional features of the arc-shaped point sets, residual noise points in the candidate sets are further eliminated. Finally, region growing is performed based on a grid index to search for and identify complete rod-shaped facility targets. This process automatically identifies the type of rod (e.g., T-shaped, F-shaped), the length and radius of horizontal and vertical rods, and other attributes, directly inputting them into the three-dimensional point attributes. Simultaneously, the algorithm fits the rod body based on the local display of the rod-shaped point cloud, accurately calculates the rod height, and confirms the rod's attribute information through necessary human-computer interaction. For traffic signs and other non-rod-shaped facilities, target point cloud regions can be manually selected from the road facility point cloud to identify their spatial location and extract key structural dimension information (e.g., the length and width of the sign). Furthermore, the point cloud clusters of pole-shaped and sign-shaped facilities such as streetlights, signs, cameras, and guardrails are projected onto the synchronized panoramic image (the exterior orientation elements of the facility point cloud clusters and the flight path). The pre-trained target detection model (such as YOLO and Faster R-CNN) is used to accurately identify the facility type (such as "speed limit 60 sign" and "streetlight - double arm"), and visual size information is extracted from it to assist in the fine contour extraction and azimuth calculation of the point cloud data.

[0156] In one embodiment of the present invention, step S15, which involves constructing lane surface road features based on the extracted lane line road features, curb road features, and guardrail road features, specifically includes: In a grayscale DOM, stop lines, median strips, green belts, pedestrian boundary lines, intersections, entrances / exits, and channelized islands are manually drawn. Elevation values ​​are assigned to the vector data, and the elevations are pressed. The extracted lane markings ("long markings"), curb stones, guardrails, and manually drawn edges are matched with the road reference plane for elevation matching; these are then summarized in the "Road Boundary Line" layer. The closure of the boundary lines is manually checked, and the 3D vector surface of the road surface (including motor vehicle, non-motor vehicle, pedestrian, median strip, and green belt surfaces) is obtained through topological relationships. Then, using the intersection stop line as the boundary line of the crossroads, line-surface segmentation is performed with the 3D vector surface of the road surface to obtain two types of 3D vector surfaces: main road and crossroads. The layer names are marked in the attribute table, and then exported to various layers based on the layer names.

[0157] In one embodiment of the present invention, step S16, which involves assigning the extracted road features to the corresponding vector layers according to their type and configuring structured semantic attributes and globally unique codes, includes: The extracted road elements (curbstones, road surface polygons, road marking vectors, facility locations, facility surfaces) are automatically assigned to different vector layers according to their type. Each element is then assigned semantic attributes such as road segment name, road segment code, administrative region code, and unique code, obtained according to preset coding rules, forming an attributed vector base map. The complete coding system is shown in Table 1 below, including the following key fields and unique codes.

[0158] Table 1. Encoding System of Semantic Attributes

[0159] In one embodiment of the present invention, before calling the matching model template from the parametric template library based on the semantic attributes of the road elements in step S17, the method further includes the step of constructing a parametric template library. Specifically, based on the extraction results of road elements of road ancillary facilities, facility point cloud clusters corresponding to various types of road ancillary facilities are obtained. The size parameters of the facility point cloud clusters are calculated and the facility point cloud clusters are projected onto a synchronized panoramic image to capture image images of the facility point cloud clusters from different perspectives. In the created 3D scene, the captured image images from different perspectives are set as reference images. Based on the calculated size parameters, basic geometry corresponding to the corresponding road ancillary facilities is created. Based on the component shapes in the background image, the basic geometry is refined to obtain a preliminary model. A parametric controller is created for the key driving parameters of the preliminary model to obtain a parametric model, which is then stored in a preset parametric template library.

[0160] In a specific example, a facility point cloud cluster is selected from the feature extraction results. Using point cloud processing tools, the dimensional parameters of each component of the facility are precisely measured and recorded. Then, based on the established high-precision registration relationship, the facility point cloud is projected onto a synchronized panoramic image, and clear image images of the facility from multiple perspectives are captured. A new scene is created in 3D modeling software, and the captured panoramic image component images are set as the background reference image for the orthogonal view. Strictly adhering to the precise dimensions measured from the point cloud, basic geometries are created in the software by inputting numerical parameters; for example, cylinders are used for rods, and cuboids for sign panels. Referring to the component shapes and details in the background image, the basic geometries are refined using relevant tools in the 3D modeling software, including subdivision, extrusion, and chamfering. The different parts of the refined model are then rationally grouped and linked. Subsequently, parametric controllers are created for the model's key driving parameters (such as the height of the support rods, the extension radius of the crossarms, and the width and height of the sign). Finally, the completed parametric model is saved according to the preset naming rules (station number + facility classification code) and archived in the system's parametric model template library.

[0161] Furthermore, step S17, which involves calling a matching model template from the parametric template library based on the semantic attributes of the road features and generating a 3D model of the corresponding road features based on the corresponding model template, specifically includes: For directional road facilities, their azimuth angles are calculated and recorded as semantic attribute fields in the corresponding vector data. Specifically, for directional facilities such as streetlights and signs, the azimuth angle is calculated. This azimuth angle is defined as the angle between the perpendicular line connecting the center point of the facility and the centerline of the road, and true north. The coordinates of the facility's center point are... Its perpendicular coordinates to the centerline of the road are The calculated azimuth angle will be automatically filled into the "azimuth" attribute field of the vector data. The formula is: ; .

[0162] Parameter files corresponding to different road feature types are extracted based on a preset mapping table. The attribute fields in these parameter files correspond one-to-one with the key driving parameters of the corresponding road feature parameterized model. The mapping table includes the correspondence between road feature types and the attribute fields contained in the parameter files. Table 2 shows a mapping table in one embodiment. In a specific example, encoded vector data is read, and a parameter file (TXT format) for modeling is generated. During the export process, the offset values ​​of the planar coordinates (X, Y) need to be set, ensuring that the coordinate values ​​of all vector data after offset are less than 10000. Then, based on the different feature types, the corresponding attribute fields are automatically identified and extracted, and the exported parameters correspond one-to-one with the parameters required by the modeling template. Finally, based on the above mapping relationship, the vector geometric information and semantic attributes are automatically written into a text file to generate structured parameter data.

[0163] Table 2 Examples of mapping relationships

[0164] Obtain the parameter file of the lane surface road elements. The parameter file of the lane surface road elements includes the coordinate position and elevation of each vertex in the ordered vertex sequence that constitutes the closed vector boundary line of the lane surface. The constrained Delaunay triangulation algorithm is used to construct the three-dimensional polygon boundary formed by the vertex coordinates of the lane surface road elements and to perform interpolation elevation assignment based on the road reference surface elevation to generate the three-dimensional surface model of the lane surface. Obtain the parameter file of the guardrail road elements. The parameter file of the guardrail road elements includes the bottom coordinate elevation of the guardrail 3D line node, the guardrail height and the road element encoding information. According to the preset guardrail sampling interval threshold, the guardrail vector line is sampled with equal arc length. The guardrail direction vector at each sampling point is calculated. Then, according to the guardrail element type information, the corresponding parametric guardrail template is called from the parametric template library. The parametric guardrail template is instantiated and placed in batches along the vector line path at the sampling points according to the guardrail direction vector to generate a continuous three-dimensional guardrail model. Obtain the parameter file for the curbstone road features. This file contains the start-point coordinates, end-point coordinates, node coordinates, elevation, and road feature encoding information for each curbstone 3D line. Based on preset curbstone height and profile width thresholds, perform equally spaced profile scans along the curbstone 3D line path. For adjacent sampling points of the curbstone road features... and At this point, each profile contains four vertices, and every two consecutive vertices on the profile... and The vertex corresponding to the next sampling point and A quadrilateral is formed, which is then divided into two triangles to form a continuous triangular mesh. The vertices on the profile at the same sampling point are connected in sequence to construct the triangular mesh of the top and bottom surfaces. The profiles corresponding to the start and end points of the path are triangulated to form closed end faces. Finally, the collection of all triangular facets constitutes a three-dimensional solid model, resulting in a three-dimensional curbstone solid model. Figure 5 The curbstone road element structure diagram is provided for the automated digital modeling method for urban roads in this embodiment of the invention.

[0165] Obtain the parameter file of pole-shaped facilities (such as streetlights and signs). The parameter file of pole-shaped facilities contains facility type, location information, geometric parameter information, azimuth angle and road element coding information. According to the facility type, call the corresponding parametric pole-shaped facility model from the parametric template library. Its geometry is determined by geometric parameter information, such as pole radius, crossarm length, height, etc., and drive the instantiation and placement of the 3D model of pole-shaped facilities through azimuth angle and location attributes. Obtain the parameter files of point facilities (such as warning posts, kilometer markers, and crash barriers). The parameter files of point facilities include the coordinate position, elevation, azimuth, facility type, and other road element coding information of the point facilities. According to the facility type, call the corresponding parametric point facility model from the parametric template library, and drive the instantiation and placement of the parametric point facility model based on the coordinate position, elevation, and azimuth of the point facilities.

[0166] In one embodiment of the present invention, step S18, which involves selecting image sources matching each model surface based on the visibility of the 3D model surface of each road element in the panoramic image, and mapping the image texture of the image source to the corresponding model surface, specifically includes: S181. For each surface of the 3D model of the road element, calculate the angle between its normal vector and the camera's line of sight in the panoramic image. If the angle is less than 90°, determine that the current model surface is visible in the panoramic image. Figure 6 This is a texture mapping effect diagram of the automated digital modeling method for urban roads proposed in an embodiment of the present invention.

[0167] In this embodiment, for each model surface, the angle between its normal vector and the camera's line of sight is calculated, and visible surfaces are selected. Let the normal vector of the model surface be... The camera's line-of-sight vector is This calculates whether the surface is visible in the current image. The formula is: ;when When g < 90°, the surface is determined to be visible in the current image.

[0168] S182. For each vertex on the model surface, project the vertex of the model surface onto the panoramic image pixel coordinate system, calculate the texture coordinates of the projected point, and map the texture coordinates of the projected point to the corresponding vertex.

[0169] In this embodiment, for the vertices of the model surface The model surface vertices are projected to the panoramic image pixel coordinate system through a coordinate transformation chain to calculate the texture coordinates (uj, vj).

[0170] S183. For each model surface, based on the resolution of the panoramic image and the viewing quality of the panoramic image relative to the model surface, match an image source for each model surface and map the image texture of the image source to the corresponding model surface.

[0171] In this embodiment, for the same model surface, an optimization strategy based on viewpoint quality and resolution is adopted to select the image with the highest score as the main texture of that surface. The formula is: ; in: The angle between the surface normal and the line of sight; Distance from the camera to the surface; Texture resolution rating; , Weighting coefficients, satisfying + =1.

[0172] For regular geometries (such as guardrails, lampposts, road markings, and asphalt), parametric UV mapping is used; for complex curved surfaces (such as curbs and irregularly shaped signs), the image with the highest score is used for projected texture mapping.

[0173] In this embodiment of the invention, step S19 involves overlaying the 3D model with fused texture information with the original laser point cloud and panoramic image in the same 3D scene to visually check the automatic modeling results and verify the geometric accuracy, semantic attribute consistency, and texture quality of the 3D model.

[0174] The specific implementation method for geometric accuracy verification includes: in a 3D scene, overlaying the road model with the original vehicle point cloud and panoramic image; selecting feature locations on the model (such as lane corner points, curb vertices, facility center points, etc.) as check points, and calculating the 3D spatial error between these check points and their corresponding points in the point cloud. S. Its formula is shown below.

[0175] in: These are the coordinates of the model's feature points; These are the coordinates of the corresponding point in the point cloud.

[0176] And calculate the average error of all checkpoints. and root mean square error ,when ≥0.3 meters or When the distance is ≥0.3 meters, backtrack to the corresponding modeling stage for parameter optimization or local reconstruction.

[0177] The specific implementation of semantic attribute consistency verification includes: based on a preset semantic coding system, verifying whether the attribute fields of model elements (such as road segment codes, facility classification codes, azimuth angles, etc.) are complete and accurate, and consistent with the extracted vector data. It also automatically checks whether key attributes such as facility type, size, and orientation match the actual objects in the point cloud and imagery.

[0178] The specific implementation methods for texture mapping quality verification include: checking the clarity, color consistency, and seam treatment effect of the model surface texture based on the texture mapping results. For areas with blurred, misaligned, or missing textures, re-projection mapping or manual repair is performed using multi-view panoramic images.

[0179] Furthermore, this invention also includes correcting deviations in the position or orientation of road ancillary facilities (such as streetlights and signs) by dragging, rotating, or inputting precise parameters, based on the results of geometric accuracy verification, semantic attribute consistency checks, and texture mapping quality assessments. For facilities with type recognition errors or mismatched geometric shapes, appropriate models are called from the parametric template library for replacement, or manual reconstruction is performed based on point clouds and images. The semantic attributes of model elements are directly modified to align with actual business rules (such as maintenance codes and traffic management classifications). Finally, for areas not fully covered by automatic modeling (such as intersections and complex green belts), manual modeling is performed using point cloud and image data to ensure model continuity and integrity. Finally, the quality-checked and corrected 3D road model is converted into a standardized format (such as OSGB, OBJ, IFC, 3DTiles) so that each model element carries its globally unique code.

[0180] The automated digital modeling method for urban roads proposed in this invention has the following beneficial effects: High degree of automation: It realizes an automated pipeline operation from data to model, which significantly reduces human intervention and improves modeling efficiency.

[0181] High modeling accuracy: Based on high-precision registered point clouds and images, the model's geometry and texture are made realistic and accurate.

[0182] Complete semantic information: The generated model is not a simple geometric object, but carries rich structured semantic attributes, which facilitates integration and application with business management systems such as BIM / GIS.

[0183] For the sake of simplicity, the method embodiments are described as a series of actions. However, those skilled in the art should understand that the embodiments of the present invention are not limited to the described order of actions, because according to the embodiments of the present invention, some steps can be performed in other orders or simultaneously. Furthermore, those skilled in the art should also understand that the embodiments described in the specification are preferred embodiments, and the actions involved are not necessarily essential to the embodiments of the present invention.

[0184] Another embodiment of the present invention provides an automated digital modeling system for urban roads, the system including functional modules for implementing the automated digital modeling method for urban roads as described in any of the preceding claims. Figure 7 A schematic diagram illustrating the structure of an automated digital modeling system for urban roads according to another embodiment of the present invention is shown. (Refer to...) Figure 7 The automated digital modeling system for urban roads in this embodiment specifically includes: The multi-source data acquisition module 201 is used to acquire the original laser point cloud and panoramic image of the target road collected simultaneously. The collaborative preprocessing module 202 is used to extract ground points from the original laser point cloud and perform spatiotemporal registration between the extracted ground point cloud and the panoramic image. The first element extraction module 203 is used to extract lane line road elements, road center line road elements, directional arrow road elements, and curbstone road elements based on the ground point cloud after spatiotemporal registration. The second element extraction module 204 is used to extract guardrail road elements and road ancillary facility road elements based on the non-ground point cloud in the original laser point cloud, and to identify the location and geometric attributes of the road ancillary facility road elements. The feature construction module 205 is used to construct lane surface road features based on the extracted lane line road features, curb road features, and guardrail road features; The semantic annotation module 206 is used to assign the extracted road features to the corresponding vector layers according to their type and configure structured semantic attributes and globally unique codes. The semantic attributes include information related to the road segment, feature type information, location information and / or geometric parameter information. The road element modeling module 207 is used to call the matching model template from the parametric template library according to the semantic attributes of the road elements, and generate the corresponding three-dimensional model of the road elements based on the corresponding model template. The model texture fusion module 208 is used to select image sources that match the surface of each model based on the visibility of the surface of the 3D model of each road element in the panoramic image, and map the image texture of the image source to the corresponding model surface to obtain a 3D model with fused texture information. The model optimization module 209 is used to overlay the 3D model with fused texture information with the original laser point cloud and panoramic image in the same 3D scene to verify the geometric accuracy, semantic attribute consistency and texture quality of the 3D model, and to correct the unqualified parts of the verification results to obtain the 3D model of the target road.

[0185] As the system implementation is basically similar to the method implementation, the description is relatively simple, and relevant parts can be found in the description of the method implementation.

[0186] Furthermore, another embodiment of the present invention provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor; when the computer program is executed by the processor, it implements the steps of the automated digital modeling method for urban roads as described above.

[0187] In addition, another embodiment of the present invention provides a computer program product on which a computer program is stored, and when the computer program is executed by a processor, it implements the steps of the automated digital modeling method for urban roads as described above.

[0188] Those skilled in the art will understand that although some embodiments herein include certain features included in other embodiments but not others, combinations of features from different embodiments are intended to be within the scope of the invention and form different embodiments. For example, any of the claimed embodiments can be used in any combination.

[0189] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. An automated digital modeling method for urban roads, characterized in that, The method includes: Acquire the original laser point cloud and panoramic image of the target road simultaneously; Ground points are extracted from the original laser point cloud, and the extracted ground point cloud is spatiotemporally registered with the panoramic image. Based on the spatiotemporally registered ground point cloud, extract road features such as lane lines, road center lines, directional arrows, and curb stones. Based on the non-ground point cloud in the original laser point cloud, guardrail road elements and road ancillary facilities road elements are extracted, and the location and geometric attributes of the road ancillary facilities road elements are identified. The lane surface road features are constructed based on the extracted lane line road features, curb road features, and guardrail road features. The extracted road features are assigned to the corresponding vector layers according to their type and configured with structured semantic attributes and globally unique codes. The semantic attributes include information related to the road segment, feature type information, location information, and / or geometric parameter information. Based on the semantic attributes of the road elements, a matching model template is called from the parametric template library, and a 3D model of the corresponding road elements is generated based on the corresponding model template. Based on the visibility of the surface of the 3D model of each road element in the panoramic image, an image source matching each model surface is selected, and the image texture of the image source is mapped to the corresponding model surface to obtain a 3D model with fused texture information. The 3D model with fused texture information is overlaid with the original laser point cloud and panoramic image in the same 3D scene to verify the geometric accuracy, semantic attribute consistency and texture quality of the 3D model, and to correct the unqualified parts of the verification results, so as to obtain the 3D model of the target road.

2. The method according to claim 1, characterized in that, Ground point extraction is performed on the original laser point cloud, including: Based on the acquisition time of each point in the original laser point cloud, find two adjacent track points for each point in the track line, and obtain a track line reference point that is synchronized with the time of each point by linear interpolation of the two adjacent track points. Calculate the 3D spatial distance and azimuth of each point relative to the track reference point, and remove points from the original laser point cloud whose 3D spatial distance is greater than a preset spatial distance threshold and whose azimuth is greater than a preset azimuth threshold to obtain a subset of candidate ground points; A cloth-based simulated filtering algorithm is used to separate ground points from non-ground points in a subset of candidate ground points.

3. The method according to claim 2, characterized in that, The method of using a cloth-simulation filtering algorithm to separate ground points and non-ground points from a subset of candidate ground points includes: Construct a two-dimensional particle mesh covering the range of a subset of candidate ground points; Under the action of gravity, the particles are driven to settle. When the particle position is updated iteratively, the updated particle position is corrected based on the elastic constraints between particles, and the particle settling velocity is corrected based on the damping force. During the particle position iteration process, collision detection is performed between each particle and the point cloud in the candidate ground point subset. When a collision point exists, the maximum elevation value in the collision point is extracted, and it is determined whether the Z coordinate value of the current particle is less than or equal to the sum of the maximum elevation value and the preset collision threshold. If the Z coordinate value of the current particle is less than or equal to the sum of the maximum elevation value and the preset collision threshold, the current particle stops settling and the maximum elevation value is taken as the final elevation value of the current particle. After the particles stop settling, for each point in the candidate ground point subset, find the cloth particle closest to the current candidate point. If the difference between the elevation value of the current candidate point and the final elevation value of the corresponding particle is less than or equal to a preset elevation threshold, then the current candidate point is determined to be a ground point.

4. The method according to claim 1, characterized in that, The extracted ground point cloud was spatiotemporally registered with the panoramic image, including: The panoramic image is calibrated for distortion and then projected onto a uniform projection surface. The projection surface is then resampled and fused to generate a spherical panoramic image. The color histogram of the overlapping region of adjacent images is calculated for the spherical panoramic image. The histogram matching method is used for color matching and the overlapping region is fused in the gradient domain to obtain a color-balanced panoramic image. Acquire synchronously collected POS data, construct a joint adjustment model based on POS data, original laser point cloud and camera parameters, and calculate the optimal system parameters using the least squares adjustment method; Based on the optimal system parameters, the extracted ground point cloud is mapped onto the panoramic image to obtain the image coordinates of the point cloud in the panoramic image coordinate system.

5. The method according to claim 1, characterized in that, Based on the spatiotemporally registered ground point cloud, lane line road features, road centerline road features, directional arrow road features, and curbstone road features are extracted, including: A sparse triangular network is constructed based on ground point clouds, and the sparse triangular network is optimized to form a road reference surface; Based on the panoramic image, color features are assigned to the point cloud data in the road reference plane to obtain a reconstructed 3D point cloud with color features. Extraction of lane line vectors, road centerline vectors, and directional arrow vectors based on reconstructed 3D point clouds; 3D line extraction of curb stones based on ground point cloud.

6. The method according to claim 1, characterized in that, Before calling a matching model template from the parametric template library based on the semantic attributes of the road features, the method further includes: Based on the extraction results of road elements of road ancillary facilities, obtain facility point cloud clusters corresponding to various road ancillary facilities, calculate the size parameters of the facility point cloud clusters, and project the facility point cloud clusters onto the synchronized panoramic image to capture image images of the facility point cloud clusters from different perspectives. In the created 3D scene, set the captured image images from different perspectives as reference images, create the basic geometry corresponding to the road ancillary facilities according to the calculated size parameters, refine the basic geometry according to the component shapes in the background image to obtain a preliminary model, create a parametric controller for the key driving parameters of the preliminary model to obtain a parametric model, and store it in the preset parametric template library.

7. The method according to claim 6, characterized in that, Based on the semantic attributes of the road features, a matching model template is retrieved from the parametric template library. A 3D model of the corresponding road features is then generated based on the corresponding model template, including: For directional road facilities among road ancillary facilities, calculate their azimuth angle and record it as a semantic attribute field of the corresponding vector data; The parameter files corresponding to different road feature types are extracted according to the preset mapping relationship table. The attribute fields included in the parameter files correspond one-to-one with the key driving parameters of the corresponding road feature parameterization model. The mapping relationship table includes the correspondence between road feature types and attribute fields contained in the parameter files. Obtain the parameter file of the lane surface road elements. The parameter file of the lane surface road elements includes the coordinate position and elevation of each vertex in the ordered vertex sequence that constitutes the closed vector boundary line of the lane surface. The constrained Delaunay triangulation algorithm is used to construct the three-dimensional polygon boundary formed by the vertex coordinates of the lane surface road elements and to perform interpolation elevation assignment based on the road reference surface elevation to generate the three-dimensional surface model of the lane surface. Obtain the parameter file of the guardrail road elements. The parameter file of the guardrail road elements includes the bottom coordinate elevation of the guardrail 3D line node, the guardrail height and the element type information. According to the preset guardrail sampling interval threshold, the guardrail vector line is sampled with equal arc length. The guardrail direction vector at each sampling point is calculated. Then, according to the guardrail element type information, the corresponding parametric guardrail template is called from the parametric template library. The parametric guardrail template is instantiated and placed in batches along the vector line path at the sampling points according to the guardrail direction vector to generate a continuous three-dimensional guardrail model. Obtain the parameter file for the curbstone road features. This file contains the start coordinates, end coordinates, node coordinates, and elevation of each curbstone's 3D line. Based on preset curbstone height and profile width thresholds, perform equally spaced profile scans along the curbstone's 3D line path. For adjacent sampling points of the curbstone road features... and At this point, each profile contains four vertices, and every two consecutive vertices on the profile... and The vertex corresponding to the next sampling point and A quadrilateral is formed, which is then divided into two triangles to form a continuous triangular mesh. The vertices on the profile at the same sampling point are connected in sequence to construct the triangular mesh of the top and bottom surfaces. The profiles corresponding to the start and end points of the path are triangulated to form closed end faces. Finally, the collection of all triangular facets constitutes a three-dimensional solid model, resulting in a three-dimensional curbstone solid model. Obtain the parameter file of the pole-shaped facility. The parameter file of the pole-shaped facility contains the facility type, location information, geometric parameter information, azimuth angle and road element coding information. According to the facility type, call the corresponding parametric pole-shaped facility model from the parametric template library. Its geometric shape is determined by the geometric parameter information, and the 3D model of the pole-shaped facility is instantiated and placed by driving the azimuth angle and location attributes. Obtain the parameter file of the point facility. The parameter file of the point facility includes the coordinate position, elevation, azimuth, facility type and road element coding information of the point facility. According to the facility type, call the corresponding parametric point facility model from the parametric template library. Drive the instantiation and placement of the parametric point facility model based on the coordinate position, elevation and azimuth of the point facility.

8. The method according to claim 1, characterized in that, Based on the visibility of the 3D model surfaces of each road element in the panoramic image, image sources matching each model surface are selected, and the image textures of the image sources are mapped to the corresponding model surfaces, including: For each surface of the 3D model of the road element, calculate the angle between its normal vector and the camera's line of sight in the panoramic image. When the angle is less than 90°, it is determined that the current model surface is visible in the panoramic image. For each vertex on the model surface, project the vertex of the model surface onto the panoramic image pixel coordinate system, calculate the texture coordinates of the projected point, and map the texture coordinates of the projected point to the corresponding vertex; For each model surface, an image source is matched for each model surface based on the resolution of the panoramic image and the spectral quality of the panoramic image relative to the model surface, and the image texture of the image source is mapped to the corresponding model surface.

9. The method according to claim 1, characterized in that, The road segment related information in the structured semantic attributes includes road segment name code, road segment code, administrative region code, road segment classification code, road segment sequence code, driving direction code, facility classification code, facility sequence code, and road element code; The globally unique code is the same as the road element code, and the coding format is: road segment code - administrative region code - road segment classification code - road segment sequence code - driving direction code - facility classification code - facility sequence code.

10. An automated digital modeling system for urban roads, characterized in that, The system includes: The multi-source data acquisition module is used to acquire the original laser point cloud and panoramic image of the target road collected simultaneously. The collaborative preprocessing module is used to extract ground points from the original laser point cloud and perform spatiotemporal registration between the extracted ground point cloud and the panoramic image. The first element extraction module is used to extract lane line road elements, road center line road elements, directional arrow road elements, and curbstone road elements from the ground point cloud after spatiotemporal registration. The second feature extraction module is used to extract guardrail road features and road ancillary facility road features based on the non-ground point cloud in the original laser point cloud, and to identify the location and geometric attributes of the road ancillary facility road features. The feature construction module is used to construct lane surface road features based on the extracted lane line road features, curb road features, and guardrail road features; The semantic annotation module is used to assign the extracted road features to the corresponding vector layers according to their type and configure structured semantic attributes and globally unique codes. The semantic attributes include information related to the road segment, feature type information, location information and / or geometric parameter information. The road element modeling module is used to call the matching model template from the parametric template library according to the semantic attributes of the road elements, and generate the corresponding 3D model of the road elements based on the corresponding model template. The model texture fusion module is used to select image sources that match the surface of each model based on the visibility of the surface of the 3D model of each road element in the panoramic image, and map the image texture of the image source to the corresponding model surface to obtain a 3D model with fused texture information. The model optimization module is used to overlay the 3D model with fused texture information with the original laser point cloud and panoramic image in the same 3D scene to verify the geometric accuracy, semantic attribute consistency and texture quality of the 3D model, and correct the unqualified parts of the verification results to obtain the 3D model of the target road.

11. A computer device, characterized in that, It includes a memory, a processor, and a computer program stored in the memory and executable on the processor; when executed by the processor, the computer program implements the steps of the method as described in any one of claims 1-9.

12. A computer program product, characterized in that, The computer program product stores a computer program that, when executed by a processor, implements the steps of the method as described in any one of claims 1-9.

Images