A method for precise registration of laser point clouds in a mobile measurement system

Abstract

This invention belongs to the field of point cloud registration processing technology, specifically relating to a method for accurate registration of laser point clouds in a mobile measurement system. The steps include: detecting repeated scanning regions in the point clouds acquired by various airborne and vehicle-mounted laser scanning systems; performing non-rigid registration of the point clouds in the repeated scanning regions using a piecewise cubic polynomial method; adaptively segmenting the non-rigidly registered vehicle-mounted laser point clouds, including coarse segmentation based on angle changes and fine segmentation based on velocity changes; performing rigid registration of the segmented vehicle-mounted laser point clouds with the airborne laser point clouds based on multi-feature constraints, registering the segmented vehicle-mounted laser point clouds onto the airborne laser point clouds; after registration, smoothing the edges of the vehicle-mounted laser point clouds, and globally optimizing the smoothed vehicle-mounted and airborne laser point cloud data to complete the final registration. This invention effectively solves the non-rigid error problem of laser point clouds in mobile measurement systems and improves registration accuracy.

Description

Technical Field

[0001] This invention belongs to the field of point cloud registration processing technology, specifically relating to a method for accurate laser point cloud registration in a mobile measurement system. Background Technology

[0002] With the rapid development of autonomous driving and 3D mapping technologies, LiDAR has become a crucial sensing device for data acquisition and scene reconstruction. However, in practical applications, the acquisition of multi-source, multi-platform data is often affected by various factors, such as the unstable attitude of the mobile acquisition platform, attitude sensor errors, GPS signal loss, or trajectory drift of the mobile platform. These factors cause non-rigid deformation of the mobile platform when acquiring point cloud data, making traditional registration methods difficult to handle effectively, thus affecting registration accuracy or leading to local registration errors. The problem of non-rigid registration is widespread and critical, and it is a pressing issue that needs to be addressed.

[0003] Furthermore, airborne and vehicle-mounted laser point clouds differ significantly in resolution, data volume, and acquisition environment, posing a significant challenge to directly fusing these two types of data. Therefore, researching an efficient non-rigid registration method for airborne and vehicle-mounted laser point clouds is of great importance. Summary of the Invention

[0004] In view of the shortcomings of the prior art, the purpose of this invention is to provide a method for accurate registration of laser point clouds in a mobile measurement system, which effectively solves the problem of non-rigid error of laser point clouds in mobile measurement systems and improves the registration accuracy.

[0005] To achieve the above objectives, this invention provides a method for precise registration of laser point clouds in a mobile measurement system, comprising the following steps:

[0006] S1. Detect duplicate scanning areas in the point clouds collected by various airborne and vehicle-mounted laser scanning systems;

[0007] S2. Use the piecewise cubic polynomial method to perform non-rigid registration on the point cloud in the repeated scanning area to eliminate the layering or ghosting of the point cloud in the repeated scanning area.

[0008] S3. Adaptive segmentation is performed on the non-rigidly registered vehicle-mounted laser point cloud, including coarse segmentation based on angle changes and fine segmentation based on velocity changes.

[0009] S4. Based on multi-feature constraints, perform rigid registration of segmented vehicle-mounted laser point clouds and airborne laser point clouds, and register the segmented vehicle-mounted laser point clouds onto the airborne laser point cloud.

[0010] S5. After registration, perform edge smoothing on the vehicle-mounted laser point cloud and then perform global optimization on the smoothed vehicle-mounted laser point cloud data and airborne laser point cloud data to complete the final registration.

[0011] As a preferred embodiment of the present invention, the method for detecting the repeated scanning region in S1 is as follows:

[0012] S1.1 Trajectory Preprocessing: Obtain trajectory data containing position and pose from various airborne laser scanning systems (ALS) and vehicle-mounted laser scanning systems (MLS), and preprocess the trajectory data.

[0013] S1.2 Neighborhood Search and Density Analysis: Set a neighborhood search radius r. For each point in the trajectory, perform a neighborhood search, count the number of point clouds in its neighborhood, and calculate the point cloud density of the point by dividing the number of point clouds in the neighborhood by the neighborhood volume. Set a density threshold I to determine whether a point belongs to a repeated scan area.

[0014] S1.3 Repeated Region Detection: Compare the point cloud density of each point in the trajectory with the set density threshold I. If the point cloud density of a certain point exceeds the density threshold I, it is determined that the point is located in the repeated scanning region and is marked as a repeated point. Obtain each repeated point and form the boundary of the repeated region by connecting adjacent repeated points, thereby identifying the complete repeated scanning region.

[0015] S1.4 Parameter Optimization and Validation: Repeat steps S1.2 and S1.3, and test with different combinations of neighborhood search radius r and density threshold I to evaluate the effect of duplicate region detection under different parameter combinations, including accuracy, completeness and computational efficiency, and select the optimal parameter combination;

[0016] S1.5 Determination of Repeated Scan Region: Based on the optimal parameter combination, the repeated scan region is obtained.

[0017] As a preferred embodiment of the present invention, the non-rigid registration process in S2 includes:

[0018] S2.1 Voxel division: The detected repetitive region point cloud is divided into multiple small cubes, i.e., voxels;

[0019] S2.2 Model Construction: A piecewise cubic polynomial model is constructed within each voxel. Assuming p is a point within the target point cloud voxel, after non-rigid registration, point p is transformed to... The non-rigid deformation of a point within the voxel can be described as follows:

[0020] (1);

[0021] in, This represents the change in position of each point in the point cloud during non-rigid registration:

[0022] (2);

[0023] In the formula, , , The coefficients are those of a cubic polynomial. , , To translate point p to its coordinates in a coordinate system with the lower left corner of the voxel as the origin; i, j, and k are the exponential variables corresponding to the three dimensions in the cubic polynomial model, respectively;

[0024] S2.3, Minimizing the error to solve for the coefficients of the cubic polynomial: Solving for the coefficients of the cubic polynomial by minimizing the objective function. , , The expression is:

[0025] (3);

[0026] In the formula, Q is a voxel in the source point cloud; P is a voxel in the target point cloud; and q is the corresponding 3D point in the source point cloud voxel.

[0027] S2.4 Non-rigid transformation: Within each voxel, based on the piecewise cubic polynomial model, the points within each voxel are transformed according to the values ​​of the cubic polynomial coefficients calculated in S2.3 to achieve non-rigid registration.

[0028] As a preferred embodiment of the present invention, in S3, the coarse segmentation based on angle change specifically involves considering both the horizontal and vertical planes of the vehicle's motion path, determining the segmentation points on the horizontal and vertical planes using the changes in heading and pitch angles, calculating the rate of change of heading and pitch angles, and comparing it with a set angular velocity change threshold to determine whether segmentation is necessary at that point. The angular velocity calculation formula is:

[0029] (4);

[0030] (5);

[0031] In the formula, This is the vertical angular velocity, i.e., the rate of change of the pitch angle; This is the horizontal angular velocity, i.e., the rate of change of the heading angle; , This represents the pitch angle between two adjacent data points; , This represents the heading angle between two adjacent data points; , This represents the timestamp of two adjacent data points when the vehicle-mounted laser scanning system collects point cloud data.

[0032] set up , These represent the threshold values ​​for the angular velocity changes in pitch and yaw angles, respectively; for two adjacent data points, when or At that time, the point cloud data between two adjacent data points is divided into different segments.

[0033] As a preferred embodiment of the present invention, in S3, the subdivision based on velocity change specifically involves, based on the coarse segmentation based on angle change, calculating the velocity change rate between two adjacent data points and comparing it with a set velocity change rate threshold. When the rate of change of velocity is greater than or equal to When this occurs, add a segmentation point at the latter of the two adjacent data points;

[0034] The formula for calculating the rate of change of velocity is:

[0035] (6);

[0036] In the formula, a is the rate of change of velocity; , This indicates the velocity between two adjacent data points.

[0037] As a preferred embodiment of the present invention, in S3, the geometric complexity of the segmented region is dynamically adjusted. , and To avoid over-segmentation or under-segmentation, the dynamic threshold adjustment process includes:

[0038] Step 1: Quantization of geometrical complexity of point cloud segments:

[0039] For each point cloud segment, principal component analysis is used to calculate the local curvature. Taking the current point as the center, a set of points within a neighborhood search radius r is taken, and the covariance matrix is ​​calculated to obtain the eigenvalues. , , Then the local curvature value at that point The calculation formula is:

[0040] (7);

[0041] In the formula, The largest eigenvalue, For intermediate eigenvalues, It is the smallest eigenvalue. ;

[0042] Let M be the number of points within a point cloud segment. The local curvature values ​​of all points within the point cloud segment are calculated. The mean curvature is obtained. :

[0043] (8);

[0044] In the formula, m is the index of a point within a point cloud segment. This represents the local curvature value at the m-th point, where m = 1, 2, ..., M;

[0045] Then the standard deviation of the curvature of the point cloud segment The calculation formula is:

[0046] (9);

[0047] Divide the point cloud segment into cubic voxels with side length L. Assume the total number of cubic voxels after the point cloud segment is G, and the point density of the g-th cube is... The calculation formula is:

[0048] (10);

[0049] In the formula, It is the number of points inside the g-th cube;

[0050] Average point density of point cloud segments The calculation formula is:

[0051] (11);

[0052] Standard deviation of point density in statistical point cloud segments :

[0053] (12);

[0054] Normalizing the standard deviation of curvature and the standard deviation of point density to [0, 1] yields the normalized geometric complexity C:

[0055] (13);

[0056] Where U represents the total number of all point cloud segments, This represents the standard deviation of curvature corresponding to the u-th point cloud segment; This represents the standard deviation of the point density corresponding to the u-th point cloud segment;

[0057] Step 2: Dynamically adjust the threshold based on geometric complexity C , and :

[0058] when At this time, the threshold is lowered to increase segmentation sensitivity. The threshold adjustment calculation formula is as follows:

[0059] (14);

[0060] In the formula, , , These are the weighting coefficients; , , For the updated , and ;

[0061] when At this time, the threshold is not adjusted;

[0062] when At this time, the threshold is increased to reduce redundant segments. The threshold adjustment calculation formula is:

[0063] (15).

[0064] As a preferred embodiment of the present invention, the process of registering the segmented vehicle-mounted laser point cloud to the airborne laser point cloud coordinate system in step S4 is as follows:

[0065] S4.1 Feature extraction, including:

[0066] Point feature extraction: In the segmented vehicle-mounted laser point cloud and airborne laser point cloud, a curvature-based feature point extraction method is adopted. By calculating the geometric curvature of each point in the point cloud in its neighborhood, points whose curvature changes exceed a set threshold are identified as obvious feature points.

[0067] Line feature extraction: The Hough transform algorithm is used to extract straight line segments or curve segments as line features from two types of point clouds.

[0068] Surface feature extraction: The RANSAC algorithm, a plane fitting algorithm, is used to extract planar regions as surface features from two types of point clouds.

[0069] S4.2 Constructing constraint relationships, including:

[0070] Point constraints: For the extracted point features, the corresponding point features in the vehicle-mounted and airborne laser point clouds are found through nearest neighbor search or feature matching methods, establishing the correspondence between point pairs, and constructing a registration model optimization criterion with point features as registration elements.

[0071] (16);

[0072] In the formula, and These are the corresponding point features in vehicle-mounted and airborne laser point clouds, respectively. It is a rotation matrix; It is a translation vector; Represents the L2 norm;

[0073] Line constraints: For the extracted line features, the corresponding line features are found through the similarity measure of line segments, the correspondence between line pairs is established, and the registration model optimization criterion with line features as registration elements is constructed.

[0074] (17);

[0075] (18);

[0076] In the formula, and These are the corresponding line features in vehicle-mounted and airborne laser point clouds, respectively. yes The direction vector; yes The direction vector; , They are , The point on;

[0077] For surface constraints, for the extracted surface features, the corresponding surface features are found through the plane's normal vector, area, and position. The correspondence between surfaces is established, and a registration model optimization criterion is constructed using surface features as registration elements.

[0078] (19);

[0079] (20);

[0080] In the formula, and These are the corresponding surface features in vehicle-mounted and airborne laser point clouds, respectively. , They are , The normal vector; , They are , The point on;

[0081] S4.3 Construct a multi-feature joint registration optimization model and solve for the optimal rigid transformation matrix. The optimization model is expressed as:

[0082] (twenty one);

[0083] In the formula, E, F, and W represent the number of point, line, and surface features, respectively. , These are weighting coefficients used to balance the contributions of point, line, and surface features;

[0084] Equation (21) is solved using the nonlinear least squares method. Based on the optimal rotation matrix and translation vector obtained from the solution, the rigid transformation matrix T is expressed as:

[0085] (twenty two);

[0086] This is the optimal rotation matrix, which represents the rotation matrix that transforms one segment of the vehicle-mounted laser point cloud into the airborne laser point cloud coordinate system. This is the optimal translation vector, which represents the translation vector that transforms one segment of the vehicle-mounted laser point cloud into the airborne laser point cloud coordinate system.

[0087] As a preferred embodiment of the present invention, in step S4, the obtained rigid transformation matrices are applied to the corresponding segmented vehicle-mounted laser point cloud, transforming it to the coordinate system of the airborne laser point cloud. A spatial error color distribution map is then constructed on the registered point cloud using the Viridis color mapping scheme for visualization inspection and verification. The steps are as follows:

[0088] Two sets of registered point cloud data are obtained, where the vehicle-mounted point cloud is the source point cloud and the airborne point cloud is the target point cloud. Using the nearest neighbor search method, for each point in the source point cloud, the point closest to it in the target point cloud is found as its "corresponding point", thus constructing a one-to-one correspondence between the source point cloud and the target point cloud.

[0089] For each pair of matching points, calculate its spatial residual and curvature difference, and normalize them to convert them into standard values ​​between 0 and 1. Calculate the joint error factor for subsequent color space mapping.

[0090] Using the Viridis color mapping scheme, the source point cloud is color-mapped based on the joint error factor to construct a spatial error color distribution map, which intuitively represents the distribution differences of registration errors in different regions.

[0091] As a preferred embodiment of the present invention, in step S5, edge smoothing of the vehicle-mounted laser point cloud is performed to eliminate the jagged effect and irregular noise of the registered vehicle-mounted point cloud edges. The smoothing process specifically includes:

[0092] S5.1 Point cloud edge detection based on local geometric features identifies edges by analyzing the geometric attributes of the point's neighborhood. First, the normal vector of the neighborhood of each point is calculated or a local surface is fitted. The degree of geometric change is quantified by the angle between the normal vectors, the curvature value, or the ratio of the eigenvalues ​​after PCA decomposition. If the difference in the direction of the normal vector of a point's neighborhood exceeds a set difference threshold, or the curvature / eigenvalue ratio exceeds a set ratio threshold, then the point is determined to be an edge point.

[0093] S5.2 Adaptive smoothing model construction: Bilateral filtering is used to balance smoothing effect and edge preservation. The filtering formula is expressed as follows:

[0094] (twenty three);

[0095] In the formula, It is an edge point; express The set of neighborhood points, express The s-th neighboring point; s is the index in the set of neighboring points, used to traverse and weight the neighboring points to calculate the smoothed point position; yes The weight of the s-th neighboring point; Indicates to Smoothed point position;

[0096] S5.3 Execute a layered smoothing strategy, perform layered processing based on edge strength, and define... , These represent the high threshold and low threshold, respectively. express The edge strength, which quantifies the degree to which it belongs to the edge, is calculated by the angle between normal vectors, curvature value, or PCA eigenvalue ratio.

[0097] when At this time, the original coordinates are preserved and no smoothing is performed;

[0098] when When necessary, apply bilateral filtering for appropriate smoothing;

[0099] when At that time, Gaussian filtering was used for noise reduction.

[0100] As a preferred embodiment of the present invention, in S5, global optimization is performed by using the established multi-feature joint registration optimization model to perform global optimization calculation based on the corresponding point, line, and surface features obtained by searching the global scope of vehicle-mounted and airborne laser point clouds, and obtaining the optimal rigid transformation matrix for registering the vehicle-mounted laser point cloud to the airborne laser point cloud in the whole scene, thereby achieving overall accurate registration of vehicle-mounted and airborne laser point clouds.

[0101] The algorithm involved in this invention can be executed by an electronic device, which includes a memory, a processor, and a computer program stored in the memory and capable of running on the processor. The processor executes the software to implement the above-mentioned algorithm calculation.

[0102] The beneficial effects of this invention are:

[0103] This invention, through a step-by-step process, first detects repetitive regions and performs non-rigid registration, then adaptive segmentation and rigid registration, and finally smooths and optimizes, effectively solving the non-rigid error problem of laser point clouds in mobile measurement systems and improving registration accuracy. The advantages of this invention are mainly reflected in:

[0104] Efficiently solves non-rigid errors: For non-rigid deformations caused by factors such as attitude instability in mobile measurement systems, techniques such as piecewise cubic polynomial models are used to accurately eliminate errors such as point cloud layering and ghosting, significantly improving registration accuracy and laying a good foundation for subsequent data processing and analysis.

[0105] Multi-feature constraint optimization registration: The optimization model is constructed by comprehensively using point, line and surface multi-feature constraints. It fully considers the geometric attributes and correspondence of different types of features, and more comprehensively describes the spatial information between point cloud data. This makes the registration results of vehicle-mounted laser point cloud and airborne laser point cloud more accurate and reliable, and effectively avoids local registration errors.

[0106] Adaptive segmentation strategy: Based on the angle and speed changes of the vehicle-mounted laser scanning system, the segmentation can be adaptively segmented. It can dynamically adjust the segmentation method according to the characteristics of the data itself, which not only improves the processing efficiency, but also enhances the adaptability to complex scenarios and different data characteristics, ensuring that each segment of data can be accurately registered.

[0107] Global optimization and smoothing: After registration, the edges of adjacent segments are smoothed and optimized globally to eliminate jagged edges and irregular noise, making the final fused point cloud data smoother and more natural overall. This improves the integrity and consistency of the data and provides higher-quality point cloud data for subsequent applications such as 3D modeling and map building. Attached Figure Description

[0108] Figure 1 This is a flowchart illustrating the principle of this invention. Detailed Implementation

[0109] The embodiments of the present invention will be further described below with reference to the accompanying drawings:

[0110] like Figure 1 As shown, a method for accurate registration of laser point clouds in a mobile measurement system includes the following steps:

[0111] S1. Detect duplicate scanning areas in the point clouds collected by various airborne and vehicle-mounted laser scanning systems;

[0112] S2. Use the piecewise cubic polynomial method to perform non-rigid registration on the point cloud in the repeated scanning area to eliminate the layering or ghosting of the point cloud in the repeated scanning area.

[0113] S3. Adaptive segmentation is performed on the non-rigidly registered vehicle-mounted laser point cloud, including coarse segmentation based on angle changes and fine segmentation based on velocity changes.

[0114] S4. Based on multi-feature constraints, perform rigid registration of segmented vehicle-mounted laser point clouds and airborne laser point clouds, and register the segmented vehicle-mounted laser point clouds onto the airborne laser point cloud.

[0115] S5. After registration, perform edge smoothing on the vehicle-mounted laser point cloud and then perform global optimization on the smoothed vehicle-mounted laser point cloud data and airborne laser point cloud data to complete the final registration.

[0116] In S1, the method for detecting repeated scan regions is as follows:

[0117] S1.1 Trajectory Preprocessing: Obtain trajectory data containing position and pose from various airborne laser scanning systems (ALS) and vehicle-mounted laser scanning systems (MLS), and preprocess the trajectory data (including filtering and noise reduction).

[0118] S1.2 Neighborhood Search and Density Analysis: Set a neighborhood search radius r (which can be adjusted according to the scanning characteristics and flight altitude of the mobile measurement system). For each point in the trajectory, perform a neighborhood search, count the number of point clouds in its neighborhood, and calculate the point cloud density of the point by dividing the number of point clouds in the neighborhood by the neighborhood volume. Set a density threshold I to determine whether a point belongs to a repeated scanning area.

[0119] S1.3 Repeated Region Detection: Compare the point cloud density of each point in the trajectory with the set density threshold I. If the point cloud density of a certain point exceeds the density threshold I, it is determined that the point is located in the repeated scanning region and is marked as a repeated point. Obtain each repeated point and form the boundary of the repeated region by connecting adjacent repeated points, thereby identifying the complete repeated scanning region.

[0120] S1.4 Parameter Optimization and Verification: Repeat steps S1.2 and S1.3, and use different combinations of neighborhood search radius r and density threshold I to test and evaluate the effect of repeated region detection under different parameter combinations, including accuracy, completeness and computational efficiency. Select the optimal parameter combination to ensure that the best fit is achieved between detection accuracy and computational efficiency.

[0121] S1.5 Determination of Repeated Scan Region: Based on the optimal parameter combination, the repeated scan region is obtained.

[0122] In S2, the non-rigid registration process includes:

[0123] S2.1 Voxel partitioning: The detected repetitive point cloud regions are divided into multiple small cubes, or voxels, to segment the point cloud data into smaller, more manageable units, allowing for the construction of a piecewise cubic polynomial model within each unit. The size of the voxel partition depends on the density of the point cloud and the required registration accuracy;

[0124] S2.2 Model Construction: A piecewise cubic polynomial model is constructed within each voxel. Assuming p is a point within the target point cloud voxel, after non-rigid registration, point p is transformed to... The non-rigid deformation of a point within the voxel can be described as follows:

[0125] (1);

[0126] in, This represents the change in position of each point in the point cloud during non-rigid registration:

[0127] (2);

[0128] In the formula, , , The coefficients are those of a cubic polynomial. , , To translate point p to its coordinates in a coordinate system with the lower left corner of the voxel as the origin; i, j, and k are the exponential variables corresponding to the three dimensions in the cubic polynomial model, respectively;

[0129] S2.3, Minimizing the error to solve for the coefficients of the cubic polynomial: Solving for the coefficients of the cubic polynomial by minimizing the objective function. , , The expression is:

[0130] (3);

[0131] In the formula, Q is a voxel in the source point cloud; P is a voxel in the target point cloud; p is a 3D point in the target point cloud voxel; and q is the corresponding 3D point in the source point cloud voxel.

[0132] S2.4 Non-rigid transformation: Within each voxel, based on the piecewise cubic polynomial model, the points within each voxel are transformed according to the values ​​of the cubic polynomial coefficients calculated in S2.3 to achieve non-rigid registration.

[0133] In S3, the coarse segmentation based on angle changes specifically considers both the horizontal and vertical planes of the vehicle's motion path. It uses the changes in heading and pitch angles to determine the segmentation points on the horizontal and vertical planes, calculates the rates of change of heading and pitch angles, and compares them with a set threshold for angular velocity change to determine whether segmentation is necessary at that point. The formula for calculating angular velocity is:

[0134] (4);

[0135] (5);

[0136] In the formula, This is the vertical angular velocity, i.e., the rate of change of the pitch angle; This is the horizontal angular velocity, i.e., the rate of change of the heading angle; , This represents the pitch angle between two adjacent data points; , This represents the heading angle between two adjacent data points; , The timestamp represents the timestamp of two adjacent data points (in a vehicle-mounted LiDAR system, a data point refers to each independent measurement point in a scan frame) when the vehicle-mounted LiDAR scanning system collects point cloud data. (The set of point cloud data obtained by the vehicle-mounted LiDAR system in a single scan at a certain point in time contains the three-dimensional information of the surrounding environment perceived by the LiDAR system at that moment.)

[0137] set up , These represent the threshold values ​​for the angular velocity changes in pitch and yaw angles, respectively; for two adjacent data points, when or At that time, the point cloud data between two adjacent data points is divided into different segments.

[0138] To segment the point cloud data more precisely, the speed variation of the onboard equipment also needs to be considered. Changes in speed may alter the density and distribution of the point cloud data, thus affecting the subsequent registration accuracy.

[0139] The sub-segmentation based on velocity change specifically involves, based on the coarse segmentation based on angle change, calculating the rate of velocity change between two adjacent data points and comparing it to a set velocity change rate threshold. When the rate of change of velocity is greater than or equal to When this occurs, add a segmentation point at the latter of the two adjacent data points;

[0140] The formula for calculating the rate of change of velocity is:

[0141] (6);

[0142] In the formula, a is the rate of change of velocity; , This indicates the velocity between two adjacent data points.

[0143] Dynamically adjust based on the geometric complexity of the segmented region , and To avoid over-segmentation or under-segmentation, the dynamic threshold adjustment process includes:

[0144] Step 1: Quantization of geometrical complexity of point cloud segments:

[0145] For each point cloud segment, principal component analysis is used to calculate the local curvature. Taking the current point as the center, a set of points within a neighborhood search radius r is taken, and the covariance matrix is ​​calculated to obtain the eigenvalues. , , Then the local curvature value at that point The calculation formula is:

[0146] (7);

[0147] In the formula, The largest eigenvalue, For intermediate eigenvalues, It is the smallest eigenvalue. ;

[0148] Let M be the number of points within a point cloud segment. The local curvature values ​​of all points within the point cloud segment are calculated. The mean curvature is obtained. :

[0149] (8);

[0150] In the formula, m is the index of a point within a point cloud segment. This represents the local curvature value at the m-th point, where m = 1, 2, ..., M;

[0151] Then the standard deviation of the curvature of the point cloud segment The calculation formula is:

[0152] (9);

[0153] Divide the point cloud segment into cubic voxels with side length L. Assume the total number of cubic voxels after the point cloud segment is G, and the point density of the g-th cube is... The calculation formula is:

[0154] (10);

[0155] In the formula, It is the number of points inside the g-th cube;

[0156] Average point density of point cloud segments The calculation formula is:

[0157] (11);

[0158] Standard deviation of point density in statistical point cloud segments :

[0159] (12);

[0160] Normalizing the standard deviation of curvature and the standard deviation of point density to [0, 1] yields the normalized geometric complexity C:

[0161] (13);

[0162] Where U represents the total number of all point cloud segments, This represents the standard deviation of curvature corresponding to the u-th point cloud segment; This represents the standard deviation of the point density corresponding to the u-th point cloud segment;

[0163] Step 2: Dynamically adjust the threshold based on geometric complexity C , and :

[0164] when At this time, the threshold is lowered to increase segmentation sensitivity. The threshold adjustment calculation formula is as follows:

[0165] (14);

[0166] In the formula, , , These are the weighting coefficients; , , For the updated , and ;

[0167] when At this time, the threshold is not adjusted;

[0168] when At this time, the threshold is increased to reduce redundant segments. The threshold adjustment calculation formula is:

[0169] (15).

[0170] In S4, the process of registering the segmented vehicle-mounted laser point cloud to the airborne laser point cloud coordinate system is as follows:

[0171] S4.1 Feature extraction, including:

[0172] Point feature extraction: In the segmented vehicle-mounted laser point cloud and airborne laser point cloud, a curvature-based feature point extraction method is adopted. By calculating the geometric curvature of each point in the point cloud in its neighborhood, points whose curvature changes exceed a set threshold are identified as obvious feature points, such as corner points and edge points. These point features should have high stability and distinguishability.

[0173] Line feature extraction utilizes the Hough transform algorithm to extract straight or curved segments as line features from two types of point clouds; these line features should reflect the geometric structure and edge information of the point cloud.

[0174] Surface feature extraction involves using the RANSAC (Random Sample Consensus) algorithm, a plane fitting algorithm, to extract planar regions as surface features from two types of point clouds. These surface features should reflect the planar structure and main planar regions of the point cloud.

[0175] S4.2 Constructing constraint relationships, including:

[0176] Point constraints: For the extracted point features, the corresponding point features in the vehicle-mounted and airborne laser point clouds are found through nearest neighbor search or feature matching methods, establishing the correspondence between point pairs, and constructing a registration model optimization criterion with point features as registration elements.

[0177] (16);

[0178] In the formula, and These are the corresponding point features in vehicle-mounted and airborne laser point clouds, respectively. It is a rotation matrix; It is a translation vector; Represents the L2 norm;

[0179] Line constraints: For the extracted line features, the corresponding line features are found through the similarity measure of line segments, the correspondence between line pairs is established, and the registration model optimization criterion with line features as registration elements is constructed.

[0180] (17);

[0181] (18);

[0182] In the formula, and These are the corresponding line features in vehicle-mounted and airborne laser point clouds, respectively. yes The direction vector; yes The direction vector; , They are , The point on;

[0183] For surface constraints, for the extracted surface features, the corresponding surface features are found through the plane's normal vector, area, and position. The correspondence between surfaces is established, and a registration model optimization criterion is constructed using surface features as registration elements.

[0184] (19);

[0185] (20);

[0186] In the formula, and These are the corresponding surface features in vehicle-mounted and airborne laser point clouds, respectively. , They are , The normal vector; , They are , The point on;

[0187] S4.3 Construct a multi-feature joint registration optimization model and solve for the optimal rigid transformation matrix. The optimization model is expressed as:

[0188] (twenty one);

[0189] In the formula, E, F, and W represent the number of point, line, and surface features, respectively. , These are weighting coefficients used to balance the contributions of point, line, and surface features;

[0190] Equation (21) is solved using the nonlinear least squares method. Based on the optimal rotation matrix and translation vector obtained from the solution, the rigid transformation matrix T is expressed as:

[0191] (twenty two);

[0192] This is the optimal rotation matrix, which represents the rotation matrix that transforms one segment of the vehicle-mounted laser point cloud into the airborne laser point cloud coordinate system. This is the optimal translation vector, which represents the translation vector that transforms one segment of the vehicle-mounted laser point cloud into the airborne laser point cloud coordinate system.

[0193] The obtained rigid transformation matrices are applied to the corresponding segmented vehicle-mounted laser point cloud to transform it into the coordinate system of the airborne laser point cloud. A spatial error color distribution map is then constructed using the Viridis color mapping scheme for visualization inspection and verification. The steps are as follows:

[0194] Two sets of registered point cloud data are obtained, where the vehicle-mounted point cloud is the source point cloud and the airborne point cloud is the target point cloud. Using the nearest neighbor search method, for each point in the source point cloud, the point closest to it in the target point cloud is found as its "corresponding point", thus constructing a one-to-one correspondence between the source point cloud and the target point cloud.

[0195] For each pair of matching points, calculate its spatial residual and curvature difference, and uniformly normalize them to convert them into standard values ​​between 0 and 1. Calculate the joint error factor (based on spatial residual and curvature difference, combined with weight parameters) for subsequent color space mapping.

[0196] Using the Viridis color mapping scheme, the source point cloud is color-mapped based on the joint error factor to construct a spatial error color distribution map, which intuitively represents the distribution differences of registration errors in different regions.

[0197] After completing the adaptive segmentation of the vehicle-mounted laser point cloud and registering each segment of the vehicle-mounted laser point cloud data to the airborne laser point cloud data, the edges between adjacent segments may appear uneven due to possible errors or data discontinuities during the registration process. To improve this issue, edge smoothing technology is employed.

[0198] In S5, edge smoothing of the vehicle-mounted laser point cloud is performed to eliminate the jagged edges and irregular noise of the registered vehicle-mounted point cloud edges, while preserving key geometric features (such as road boundaries and building outlines). The smoothing process is as follows:

[0199] S5.1 Point cloud edge detection based on local geometric features identifies edges by analyzing the geometric attributes of the point's neighborhood (such as normal vectors, curvature, or eigenvalues ​​of the covariance matrix). First, for each point, the normal vector of its neighborhood (such as K-nearest neighbors or radius search) or a fitted local surface is calculated. The degree of geometric change is quantified by the angle between normal vectors, curvature values, or the ratio of eigenvalues ​​after PCA decomposition (such as linearity / flatness). If the difference in the direction of the normal vectors of a point's neighborhood exceeds a set difference threshold, or the curvature / eigenvalue ratio exceeds a set ratio threshold, then the point is determined to be an edge point.

[0200] S5.2 Adaptive smoothing model construction: Bilateral filtering is used to balance smoothing effect and edge preservation. The filtering formula is expressed as follows:

[0201] (twenty three);

[0202] In the formula, It is an edge point; express The set of neighborhood points, express The s-th neighboring point; s is the index in the set of neighboring points, used to traverse and weight the neighboring points to calculate the smoothed point position; yes The weight of the s-th neighboring point is determined by spatial distance (closest to the nearest neighbor). The weight of points with similar normal vectors / curvatures is determined by both the point weight and feature similarity (points with similar normal vectors / curvatures have higher weights). Indicates to Smoothed point position;

[0203] S5.3 Execute a layered smoothing strategy, performing layered processing based on edge strength to avoid over-smoothing key features, and define... , These represent high threshold (used to distinguish between strong and weak edges) and low threshold (used to distinguish between weak edges and noise), respectively. express The edge strength, which quantifies the degree to which it belongs to the edge, is calculated by the angle between normal vectors, curvature value, or PCA eigenvalue ratio.

[0204] when At this time, the original coordinates are preserved and no smoothing is performed;

[0205] when When necessary, apply bilateral filtering for appropriate smoothing;

[0206] when At that time, Gaussian filtering was used for noise reduction.

[0207] The global optimization involves using the corresponding point, line, and surface features obtained from a global search of vehicle-mounted and airborne laser point clouds, and then employing an established multi-feature joint registration optimization model to perform a global optimization solution, thereby obtaining the optimal rigid transformation matrix for registering vehicle-mounted laser point clouds to airborne laser point clouds across the entire scene. This enables precise overall registration of vehicle-mounted and airborne laser point clouds.

Claims

1. A method for precise registration of laser point clouds in a mobile measurement system, characterized in that... Includes the following steps: S1. Detect duplicate scanning areas in the point clouds collected by various airborne and vehicle-mounted laser scanning systems; S2. Use the piecewise cubic polynomial method to perform non-rigid registration on the point cloud in the repeated scanning area to eliminate the layering or ghosting of the point cloud in the repeated scanning area. S3. Adaptive segmentation is performed on the non-rigidly registered vehicle-mounted laser point cloud, including coarse segmentation based on angle change and angular velocity change thresholds, and fine segmentation based on velocity change and velocity change rate thresholds, where the angles include pitch angle and yaw angle. The threshold for angular velocity change of the pitch angle is dynamically adjusted based on the geometric complexity C of the segmented region. Threshold for angular velocity change of heading angle and the threshold of rate of change of velocity :when When, lower the threshold and increase segmentation sensitivity; when When the threshold is not adjusted, when... At the same time, increase the threshold and reduce redundant segments; S4. Based on multi-feature constraints, perform rigid registration of segmented vehicle-mounted laser point clouds and airborne laser point clouds, and register the segmented vehicle-mounted laser point clouds onto the airborne laser point cloud. S5. After registration, perform edge smoothing on the vehicle-mounted laser point cloud and then perform global optimization on the smoothed vehicle-mounted laser point cloud data and airborne laser point cloud data to complete the final registration.

2. The method for precise registration of laser point clouds in a mobile measurement system according to claim 1, characterized in that: In S1, the method for detecting repeated scan regions is as follows: S1.1 Trajectory Preprocessing: Obtain trajectory data containing position and pose from various airborne laser scanning systems (ALS) and vehicle-mounted laser scanning systems (MLS), and preprocess the trajectory data. S1.2 Neighborhood Search and Density Analysis: Set a neighborhood search radius r. For each point in the trajectory, perform a neighborhood search, count the number of point clouds in its neighborhood, and calculate the point cloud density of the point by dividing the number of point clouds in the neighborhood by the neighborhood volume. Set a density threshold I to determine whether a point belongs to a repeated scan area. S1.3 Repeated Region Detection: Compare the point cloud density of each point in the trajectory with the set density threshold I. If the point cloud density of a certain point exceeds the density threshold I, it is determined that the point is located in the repeated scanning region and is marked as a repeated point. Obtain each repeated point and form the boundary of the repeated region by connecting adjacent repeated points, thereby identifying the complete repeated scanning region. S1.4 Parameter Optimization and Validation: Repeat steps S1.2 and S1.3, and test with different combinations of neighborhood search radius r and density threshold I to evaluate the effect of duplicate region detection under different parameter combinations, including accuracy, completeness and computational efficiency, and select the optimal parameter combination; S1.5 Determination of Repeated Scan Region: Based on the optimal parameter combination, the repeated scan region is obtained.

3. The method for precise registration of laser point clouds in a mobile measurement system according to claim 1, characterized in that: In S2, the non-rigid registration process includes: S2.1 Voxel division: The detected repetitive region point cloud is divided into multiple small cubes, i.e., voxels; S2.2 Model Construction: A piecewise cubic polynomial model is constructed within each voxel. Assuming p is a point within the target point cloud voxel, after non-rigid registration, point p is transformed to... The non-rigid deformation of a point within the voxel can be described as follows: （1）； in, This represents the change in position of each point in the point cloud during non-rigid registration: （2）； In the formula, , , The coefficients are those of a cubic polynomial. , , To translate point p to its coordinates in a coordinate system with the lower left corner of the voxel as the origin; i, j, and k are the exponential variables corresponding to the three dimensions in the cubic polynomial model, respectively; S2.3, Minimizing the error to solve for the coefficients of the cubic polynomial: Solving for the coefficients of the cubic polynomial by minimizing the objective function. , , The expression is: （3）； In the formula, Q is a voxel in the source point cloud; P is a voxel in the target point cloud; and q is the corresponding three points in the source point cloud voxel. S2.4 Non-rigid transformation: Within each voxel, based on the piecewise cubic polynomial model, the points within each voxel are transformed according to the values ​​of the cubic polynomial coefficients calculated in S2.3 to achieve non-rigid registration.

4. The method for precise registration of laser point clouds in a mobile measurement system according to claim 1, characterized in that: In S3, the coarse segmentation based on angle changes specifically involves considering both the horizontal and vertical planes of the vehicle's motion path. The segmentation points on the horizontal and vertical planes are determined using the changes in the heading and pitch angles. The rates of change of the heading and pitch angles are calculated and compared with a set threshold for angular velocity change to determine whether segmentation is necessary at that point. The formula for calculating angular velocity is: （4）； （5）； In the formula, This is the vertical angular velocity, i.e., the rate of change of the pitch angle; This is the horizontal angular velocity, i.e., the rate of change of the heading angle; , This represents the pitch angle between two adjacent data points; , This represents the heading angle between two adjacent data points; , This represents the timestamp of two adjacent data points when the vehicle-mounted laser scanning system collects point cloud data. set up , These represent the threshold values ​​for the angular velocity changes in pitch and yaw angles, respectively; for two adjacent data points, when or At that time, the point cloud data between two adjacent data points is divided into different segments.

5. The method for precise registration of laser point clouds in a mobile measurement system according to claim 4, characterized in that: In S3, the subdivision based on velocity change specifically involves, based on the coarse segmentation based on angle change, calculating the rate of velocity change between two adjacent data points and comparing it with a set velocity change rate threshold. When the rate of change of velocity is greater than or equal to When this occurs, add a segmentation point at the latter of the two adjacent data points; The formula for calculating the rate of change of velocity is: （6）； In the formula, a is the rate of change of velocity; , This indicates the velocity between two adjacent data points.

6. The method for precise registration of laser point clouds in a mobile measurement system according to claim 5, characterized in that: In S3, the geometric complexity of the segmented region is dynamically adjusted. , and To avoid over-segmentation or under-segmentation, the dynamic threshold adjustment process includes: Step 1: Quantization of geometrical complexity of point cloud segments: For each point cloud segment, principal component analysis is used to calculate the local curvature. Taking the current point as the center, a set of points within a neighborhood search radius r is taken, and the covariance matrix is ​​calculated to obtain the eigenvalues. , , Then the local curvature value at that point The calculation formula is: （7）； In the formula, The largest eigenvalue, For intermediate eigenvalues, It is the smallest eigenvalue. ; Let M be the number of points within a point cloud segment. The local curvature values ​​of all points within the point cloud segment are calculated. The mean curvature is obtained. : （8）； In the formula, m is the index of a point within a point cloud segment. This represents the local curvature value at the m-th point, where m = 1, 2, ..., M; Then the standard deviation of the curvature of the point cloud segment The calculation formula is: （9）； Divide the point cloud segment into cubic voxels with side length L. Assume the total number of cubic voxels after the point cloud segment is G, and the point density of the g-th cube is... The calculation formula is: （10）； In the formula, It is the number of points inside the g-th cube; Average point density of point cloud segments The calculation formula is: （11）； Standard deviation of point density in statistical point cloud segments : （12）； Normalizing the standard deviation of curvature and the standard deviation of point density to [0, 1] yields the normalized geometric complexity C: （13）； Where U represents the total number of all point cloud segments, This represents the standard deviation of curvature corresponding to the u-th point cloud segment; This represents the standard deviation of the point density corresponding to the u-th point cloud segment; Step 2: Dynamically adjust the threshold based on geometric complexity C , and : when At this time, the threshold is lowered to increase segmentation sensitivity. The threshold adjustment calculation formula is as follows: （14）； In the formula, , , These are the weighting coefficients; , , For the updated , and ; when At this time, the threshold is not adjusted; when At this time, the threshold is increased to reduce redundant segments. The threshold adjustment calculation formula is: （15）。 7. The method for precise registration of laser point clouds in a mobile measurement system according to claim 1, characterized in that: In S4, the process of registering the segmented vehicle-mounted laser point cloud to the airborne laser point cloud coordinate system is as follows: S4.1 Feature extraction, including: Point feature extraction: In the segmented vehicle-mounted laser point cloud and airborne laser point cloud, a curvature-based feature point extraction method is adopted. By calculating the geometric curvature of each point in the point cloud in its neighborhood, points whose curvature changes exceed a set threshold are identified as obvious feature points. Line feature extraction: The Hough transform algorithm is used to extract straight line segments or curve segments as line features from two types of point clouds. Surface feature extraction: The RANSAC algorithm, a plane fitting algorithm, is used to extract planar regions as surface features from two types of point clouds. S4.2 Constructing constraint relationships, including: Point constraints: For the extracted point features, the corresponding point features in the vehicle-mounted and airborne laser point clouds are found through nearest neighbor search or feature matching methods, establishing the correspondence between point pairs, and constructing a registration model optimization criterion with point features as registration elements. （16）； In the formula, and These are the corresponding point features in vehicle-mounted and airborne laser point clouds, respectively. It is a rotation matrix; It is a translation vector; Represents the L2 norm; Line constraints are applied. For the extracted line features, the corresponding line features are found through the similarity measure of line segments, the correspondence between line pairs is established, and the registration model optimization criterion with line features as registration elements is constructed. （17）； （18）； In the formula, and These are the corresponding line features in vehicle-mounted and airborne laser point clouds, respectively. yes The direction vector; yes The direction vector; , They are , The point on; For surface constraints, for the extracted surface features, the corresponding surface features are found through the plane's normal vector, area, and position. The correspondence between surfaces is established, and a registration model optimization criterion is constructed using surface features as registration elements. （19）； （20）； In the formula, and These are the corresponding surface features in vehicle-mounted and airborne laser point clouds, respectively. , They are , The normal vector; , They are , The point on; S4.3 Construct a multi-feature joint registration optimization model and solve for the optimal rigid transformation matrix. The optimization model is expressed as: （21）； In the formula, E, F, and W represent the number of point, line, and surface features, respectively. , These are weighting coefficients used to balance the contributions of point, line, and surface features; Equation (21) is solved using the nonlinear least squares method. Based on the optimal rotation matrix and translation vector obtained from the solution, the rigid transformation matrix T is expressed as: （22）； This is the optimal rotation matrix, which represents the rotation matrix that transforms one segment of the vehicle-mounted laser point cloud into the airborne laser point cloud coordinate system. This is the optimal translation vector, which represents the translation vector that transforms one segment of the vehicle-mounted laser point cloud into the airborne laser point cloud coordinate system.

8. The method for precise registration of laser point clouds in a mobile measurement system according to claim 7, characterized in that: In step S4, the obtained rigid transformation matrices are applied to the corresponding segmented vehicle-mounted laser point cloud, transforming it to the coordinate system of the airborne laser point cloud. A spatial error color distribution map is then constructed using the Viridis color mapping scheme for visualization inspection and verification. The steps are as follows: Two sets of registered point cloud data are obtained, where the vehicle-mounted point cloud is the source point cloud and the airborne point cloud is the target point cloud. Using the nearest neighbor search method, for each point in the source point cloud, the point closest to it in the target point cloud is found as its "corresponding point", thus constructing a one-to-one correspondence between the source point cloud and the target point cloud. For each pair of matching points, calculate its spatial residual and curvature difference, and normalize them to convert them into standard values ​​between 0 and 1. Calculate the joint error factor for subsequent color space mapping. Using the Viridis color mapping scheme, the source point cloud is color-mapped based on the joint error factor to construct a spatial error color distribution map, which intuitively represents the distribution differences of registration errors in different regions.

9. The method for precise registration of laser point clouds in a mobile measurement system according to claim 1, characterized in that: In step S5, edge smoothing of the vehicle-mounted laser point cloud is performed to eliminate the jagged edges and irregular noise of the registered vehicle-mounted point cloud edges. The smoothing process is as follows: S5.1 Point cloud edge detection based on local geometric features identifies edges by analyzing the geometric attributes of the point's neighborhood. First, the normal vector of the neighborhood of each point is calculated or a local surface is fitted. The degree of geometric change is quantified by the angle between the normal vectors, the curvature value, or the ratio of the eigenvalues ​​after PCA decomposition. If the difference in the direction of the normal vector of a point's neighborhood exceeds a set difference threshold, or the curvature / eigenvalue ratio exceeds a set ratio threshold, then the point is determined to be an edge point. S5.2 Adaptive smoothing model construction: Bilateral filtering is used to balance smoothing effect and edge preservation. The filtering formula is expressed as follows: （23）； In the formula, It is an edge point; express The set of neighborhood points, express The s-th neighboring point; s is the index in the set of neighboring points, used to traverse and weight the neighboring points to calculate the smoothed point position; yes The weight of the s-th neighboring point; Indicates to Smoothed point position; S5.3 Execute a layered smoothing strategy, perform layered processing based on edge strength, and define... , These represent the high threshold and low threshold, respectively. express The edge strength, which quantifies the degree to which it belongs to the edge, is calculated by the angle between normal vectors, curvature value, or PCA eigenvalue ratio. when At this time, the original coordinates are preserved and no smoothing is performed; when When necessary, apply bilateral filtering for appropriate smoothing; when At that time, Gaussian filtering was used for noise reduction.

10. The method for precise registration of laser point clouds in a mobile measurement system according to claim 7, characterized in that: In S5, global optimization is performed by using the established multi-feature joint registration optimization model to perform global optimization calculation based on the corresponding point, line, and surface features obtained by searching the global scope of vehicle-mounted and airborne laser point clouds. This yields the optimal rigid transformation matrix for registering the vehicle-mounted laser point cloud to the airborne laser point cloud across the entire scene, thereby achieving overall accurate registration of the vehicle-mounted and airborne laser point clouds.

Images