Target reference based highway traffic structure three-dimensional geometry multi-layer calibration method and system

Abstract

The present application belongs to the field of highway traffic geometry measurement, and discloses a highway traffic structure three-dimensional geometry multi-layer calibration method and system based on target reference, comprising: collecting point cloud data of highway infrastructure structure by a three-dimensional laser scanner and carrying out denoising pretreatment; based on the denoising pretreated point cloud data, calibrating the plane error of the three-dimensional laser scanning device by a standard plane target; based on the denoising pretreated point cloud data, calibrating the three-dimensional error of the three-dimensional laser scanning device by a standard spherical target; based on the calibration results of the plane error and the three-dimensional error of the three-dimensional laser scanning device, converting different station coordinates by a four-parameter method, and registering the same name area by the maximum geometric error. The present application improves the reliability, accuracy and adaptability of the highway infrastructure scanning device.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of highway traffic geometric quantity measurement, and particularly relates to a multi-layer calibration method and system for three-dimensional geometric quantity of highway traffic structure based on target reference. BACKGROUND

[0002] In the field of geometric quantity measurement, the existing technology often adopts an error calibration method based on a target, but still has the following key problems:

[0003] 1. Serious noise interference:

[0004] The traditional denoising method (such as mean filtering) has limited ability to remove global random noise, and cannot be combined with the spatial characteristics of three-dimensional point cloud in depth, resulting in high effective data deletion rate;

[0005] 2. Insufficient plane calibration accuracy:

[0006] The traditional plane fitting algorithm such as least squares is sensitive to point cloud density distribution, and in the non-uniform sampling scene, the geometric center recognition error cannot meet the high-precision target calibration demand;

[0007] 3. Limitations of stereo parameter optimization:

[0008] The existing spherical calibration method (such as least squares) only optimizes single-direction residual error, without fusing multi-dimensional gradient information, resulting in significant measurement deviation of the sphere center and radius (typical error ≥3mm);

[0009] 4. Algorithm efficiency bottleneck:

[0010] The traditional iterative algorithm needs multiple random sampling, and has high computational complexity, which is difficult to support real-time processing of large-scale point cloud.

[0011] The above problems restrict the development of high-precision geometric quantity measurement technology, and an calibration system that takes into account theoretical rigor and engineering efficiency is urgently needed. Therefore, it is of great significance to study the target multi-layer calibration technology for three-dimensional geometric quantity measurement of highway traffic structure, to improve the reliability, accuracy and adaptability of highway infrastructure scanning equipment. SUMMARY

[0012] To solve the problems in the prior art, the application provides a multi-layer calibration method and system for three-dimensional geometric quantity of highway traffic structure based on target reference, to solve the problem of lack of calibration technology that takes into account theoretical rigor and engineering efficiency in the industry, and to improve the reliability, accuracy and adaptability of highway infrastructure scanning equipment.

[0013] To achieve the above purpose, the application provides the following solutions:

[0014] The multi-layer calibration method for three-dimensional geometric quantity of highway traffic structure based on target reference comprises:

[0015] Collecting point cloud data of highway infrastructure structure by three-dimensional laser scanner and carrying out denoising preprocessing;

[0016] Calibrating the plane error of the three-dimensional laser scanning device based on the denoising preprocessed point cloud data by using a standard plane target;

[0017] Calibrating the three-dimensional error of the three-dimensional laser scanning device based on the denoising preprocessed point cloud data by using a standard spherical target;

[0018] Based on the calibration results of the plane error and the three-dimensional error of the three-dimensional laser scanning device, converting different station coordinates by four-parameter method, and using maximum geometric error for homonymous area registration.

[0019] Preferably, collecting point cloud data of highway infrastructure structure by three-dimensional laser scanner and carrying out denoising preprocessing, comprising:

[0020] Converting three-dimensional point cloud coordinates to the Euclidean distance of the origin, and constructing a one-dimensional data set;

[0021] Calculating the mean value and Z-score of the one-dimensional data set, i.e. Z i ;

[0022] According to the calculation result of Z-score, setting a dynamic threshold;

[0023] Eliminating points satisfying ∣Z i ∣> threshold.

[0024] Preferably, calibrating the plane error of the three-dimensional laser scanning device based on the denoising preprocessed point cloud data by using a standard plane target, comprising:

[0025] Obtaining the area of the to-be-measured region;

[0026] Calculating the point cloud density distribution, and calculating the neighborhood radius according to the area of the to-be-measured region and the point cloud density distribution;

[0027] According to the neighborhood radius, calculating the initial value of the geometric center;

[0028] Based on the initial value of the geometric center, iteratively optimizing the geometric center coordinates by weighted mean algorithm to complete the calibration of the plane error of the three-dimensional laser scanning device;

[0029] Wherein, obtaining the area of the to-be-measured region, comprising:

[0030]

[0031] Calculating the point cloud density distribution, comprising:

[0032]

[0033] In the formula, N r is the actual sampling point number of the to-be-measured area, which is determined according to field measurement requirements; N tot is a theoretical sampling point number; λ is an empirical parameter of point cloud density; r is a neighborhood radius; and D is a point cloud density distribution function;

[0034] According to the area of the to-be-measured area and the point cloud density distribution, a neighborhood radius is calculated, including:

[0035]

[0036] According to the neighborhood radius, an initial value of a geometric center is calculated, including:

[0037] A neighborhood is constructed with the neighborhood radius r, containing N point cloud data, and an initial geometric center coordinate C0=(x0, y0) is calculated:

[0038]

[0039] In the formula, N i is a sampling point number of an i-th sub-area in the neighborhood, and n is a number of sub-areas;

[0040] Based on the initial value of the geometric center, a geometric center coordinate is iteratively optimized by a weighted mean value algorithm, to complete planar error calibration of the point cloud data, including:

[0041] Based on a distance d i of the point cloud data to a temporary geometric center, a weight function ω i is defined:

[0042]

[0043] In the formula, is a temporary geometric center coordinate of the k-th iteration;

[0044] The geometric center coordinate C=(x, y) is iteratively optimized by the weighted mean value algorithm:

[0045]

[0046] An iteration termination condition is that an absolute value of a difference between two adjacent times is less than a threshold ∈.

[0047] Preferably, based on the point cloud data after denoising preprocessing, a three-dimensional laser scanning device is calibrated for a stereoscopic error by using a standard spherical target, including:

[0048] A spherical model is established;

[0049] Based on the spherical model, a loss function is constructed;

[0050] based on the loss function, a gradient and a Hessian matrix are calculated;

[0051] based on the gradient and the Hessian matrix, a parameter iterative update is performed to complete the calibration of the three-dimensional laser scanning device;

[0052] wherein the spherical model is established, comprising:

[0053] The preset spherical equation is:

[0054] (x-a) 2 +(y-b) 2 +(z-c) 2 =r 2 ;

[0055] wherein (x, y, z) is the coordinate of a point on the preset sphere, (a, b, c) is the three-dimensional coordinate of the center of the sphere, and r is the radius;

[0056] based on the spherical model, a loss function is constructed, comprising:

[0057] minimizing the sum of the square of the distance of all point clouds to the fitted sphere, and the loss function is:

[0058]

[0059] based on the loss function, a gradient and a Hessian matrix are calculated, comprising:

[0060] the gradient vector is:

[0061]

[0062] the Hessian matrix is:

[0063]

[0064] based on the gradient and the Hessian matrix, a parameter iterative update is performed, comprising:

[0065] the update formula of the parameter vector is:

[0066]

[0067] wherein, is the updated parameter vector, is the parameter vector before updating, H is the Hessian matrix, is the gradient vector;

[0068] the iteration termination condition is: the radius error Δr≤0.003m or the maximum number of iterations is reached.

[0069] Preferably, based on the calibration results of the planar error and the stereoscopic error of the three-dimensional laser scanning device, the different station coordinates are converted by a four-parameter method, and the same name area registration is performed by using the maximum geometric error, including:

[0070] Calibrating the four-parameter coordinates of the point cloud data;

[0071] Optimizing the calibration results by using the NDT method;

[0072] Based on the optimized structure, the same name area of the point cloud data is fixed;

[0073] Wherein, the four-parameter coordinates of the point cloud data are calibrated, including:

[0074] The same name point pair of the source point cloud and the target point cloud is preset as: the source point cloud coordinate: P t,i =(x i ,y i ), the target point cloud coordinate: P s,i =(x s,i ,y s,i ); the planar coordinate calibration is performed by using the four-parameter method, and the calculation formula is:

[0075]

[0076] Wherein, a is a scale transformation parameter, θ is a rotation angle, t x and t y are translation amounts;

[0077] The calibration results are optimized by using the NDT method, including:

[0078] The target point cloud is divided into a fixed-size cubic grid;

[0079] The mean value μ and the covariance matrix Γ of the point cloud in each grid are calculated, and a probability density model is constructed;

[0080]

[0081] Based on the probability density model, the target is optimized by using the maximum overall likelihood function;

[0082]

[0083] The transformation matrix is updated by using the Levenberg-Marquardt algorithm until the difference between the adjacent two times of transformation is less than a threshold value:

[0084]

[0085] In the formula, M k+1 and M ktransform parameters of the kth and (k+1)th iteration, a sum of gradients of the log-likelihood function of all grid cells or data points, a sum of second-order derivatives of the log-likelihood function, λ is a damping factor of the algorithm, used to balance the gradient descent and the Gauss-Newton method, when λ is larger, the algorithm is close to the gradient descent, when λ is smaller, the algorithm is close to the Gauss-Newton method, I is a bit matrix, combined with λ to ensure that the matrix is invertible, avoiding singular problems;

[0086] Based on the optimized structure, the same name area of the point cloud data is fixed, including:

[0087] Calculate the maximum Euclidean distance between the source point cloud and the target point cloud:

[0088]

[0089] Threshold setting: filter points that meet d i ≤ε, wherein the threshold

[0090] Match the point cloud by using the KDTree structure of the Open3D library.

[0091] The application also provides a target reference-based multi-layer calibration system for three-dimensional geometric quantities of a highway traffic structure, which is used to implement the method and comprises a preprocessing module, a plane error calibration module, a stereo error calibration module and a region matching module.

[0092] The preprocessing module is used to collect point cloud data of a highway infrastructure structure by using a three-dimensional laser scanner and to perform denoising preprocessing.

[0093] The plane error calibration module is used to calibrate the plane error of the three-dimensional laser scanning device by using a standard plane target based on the denoising-preprocessed point cloud data.

[0094] The stereo error calibration module is used to calibrate the stereo error of the three-dimensional laser scanning device by using a standard spherical target based on the denoising-preprocessed point cloud data.

[0095] The region matching module is used to convert different station coordinates by using a four-parameter method and to perform same-name region registration by using the maximum geometric error based on the calibration results of the plane error and the stereo error of the three-dimensional laser scanning device.

[0096] Preferably, the point cloud data of the highway infrastructure structure is collected by using a three-dimensional laser scanner and is subjected to denoising preprocessing, including:

[0097] Convert the three-dimensional point cloud coordinates to the Euclidean distance of the origin, and construct a one-dimensional data set;

[0098] The mean value, Z-score, i.e. Z i ;

[0099] According to the calculation result of the Z-score, a dynamic threshold is set;

[0100] The points satisfying ∣Z i ∣> threshold are removed.

[0101] Preferably, based on the point cloud data after denoising preprocessing, a standard plane target is used to calibrate the plane error of the three-dimensional laser scanning device, including:

[0102] Obtaining the area of the to-be-measured region;

[0103] Calculating the point cloud density distribution, and calculating the neighborhood radius according to the area of the to-be-measured region and the point cloud density distribution;

[0104] According to the neighborhood radius, calculating the initial value of the geometric center;

[0105] Based on the initial value of the geometric center, the geometric center coordinates are iteratively optimized by a weighted mean value algorithm to complete the calibration of the plane error of the three-dimensional laser scanning device;

[0106] Wherein, obtaining the area of the to-be-measured region, including:

[0107]

[0108] Calculating the point cloud density distribution, including:

[0109]

[0110] In the formula, N r is the actual sampling point number of the to-be-measured region, which is determined according to the field measurement requirement; N tot is the theoretical sampling point number; λ is an empirical parameter of the point cloud density; r is the neighborhood radius; and D is the point cloud density distribution function;

[0111] According to the area of the to-be-measured region and the point cloud density distribution, calculating the neighborhood radius, including:

[0112]

[0113] According to the neighborhood radius, calculating the initial value of the geometric center, including:

[0114] A neighborhood is constructed with the neighborhood radius r, containing N point cloud data, and the initial geometric center coordinates C0=(x0, y0) are calculated:

[0115]

[0116] In the formula, N iis the number of sampling points in the i-th sub-region in the neighborhood, and n is the number of sub-regions;

[0117] Based on the initial value of the geometric center, the geometric center coordinates are iteratively optimized by a weighted mean algorithm to complete the plane error calibration of the point cloud data, including:

[0118] Based on the distance d of the point cloud data to the temporary geometric center i , a weight function ω i is defined:

[0119]

[0120] In the formula, is the temporary geometric center coordinate of the k-th iteration;

[0121] The geometric center coordinates C=(x,y) are iteratively optimized by a weighted mean algorithm:

[0122]

[0123] The iteration termination condition is that the absolute value of the difference between the coordinates of the adjacent two times is less than a threshold ∈.

[0124] Preferably, based on the point cloud data after denoising preprocessing, the stereoscopic error of the three-dimensional laser scanning device is calibrated by using a standard spherical target, including:

[0125] A spherical model is established;

[0126] Based on the spherical model, a loss function is constructed;

[0127] Based on the loss function, the gradient and Hessian matrix are calculated;

[0128] Based on the gradient and Hessian matrix, parameter iterative updating is performed to complete the calibration of the stereoscopic error of the three-dimensional laser scanning device;

[0129] Wherein, the spherical model is established, including:

[0130] The preset spherical equation is:

[0131] (x-a) 2 +(y-b) 2 +(z-c) 2 =r 2 ;

[0132] Wherein, (x,y,z) is the coordinate of a point on the preset sphere, (a,b,c) is the three-dimensional coordinate of the sphere center, and r is the radius;

[0133] Based on the spherical model, a loss function is constructed, including:

[0134] Minimizing the sum of squared distances of all point clouds to the fitted sphere, the loss function is:

[0135]

[0136] Based on the loss function, the gradient and Hessian matrix are calculated, including:

[0137] The gradient vector is:

[0138]

[0139] The Hessian matrix is:

[0140]

[0141] Based on the gradient and Hessian matrix, the parameter iterative update is performed, including:

[0142] The update formula of the parameter vector is:

[0143]

[0144] wherein, is the updated parameter vector, is the parameter vector before updating, H is the Hessian matrix, is the gradient vector;

[0145] The iteration termination condition is: the radius error Δr≤0.003m or the maximum iteration number is reached.

[0146] Preferably, based on the calibration results of the plane error and the stereo error of the three-dimensional laser scanning device, the four-parameter method is used to convert the coordinates of different stations, and the maximum geometric error is used for homonymous area registration, including:

[0147] Calibrating the four-parameter coordinates of the point cloud data;

[0148] Using the NDT method to optimize the calibration results;

[0149] Based on the optimized structure, the homonymous areas of the point cloud data are fixed;

[0150] Wherein, the calibration of the four-parameter coordinates of the point cloud data includes:

[0151] The homonymous point pairs of the source point cloud and the target point cloud are preset as: the source point cloud coordinates: P t,i =(x i ,y i ), the target point cloud coordinates: P s,i =(x s,i ,y s,i); the plane coordinate calibration is performed by a four-parameter method, and the calculation formula is:

[0152]

[0153] wherein a is a scale transformation parameter, θ is a rotation angle, t x and t y are translation amounts;

[0154] The NDT method is used to optimize the calibration result, including:

[0155] The target point cloud is divided into a fixed-size cubic grid;

[0156] The mean μ and covariance matrix Γ of the point cloud in each grid are calculated to construct a probability density model;

[0157]

[0158] Based on the probability density model, the target is optimized by maximizing the overall likelihood function;

[0159]

[0160] The transformation matrix is updated by the Levenberg-Marquardt algorithm until the difference between the adjacent two transformations is less than a threshold:

[0161]

[0162] wherein M k+1 and M k represent the transformation parameters of the kth and k+1th iterations, represents the sum of the gradients of the log-likelihood function of all grid units or data points, is the sum of the second derivatives of the log-likelihood function, λ is the damping factor of the algorithm, used to balance the gradient descent and the Gauss-Newton method, when λ is larger, the algorithm is close to the gradient descent, when λ is smaller, the algorithm is close to the Gauss-Newton method, I is an identity matrix, combined with λ to ensure that the matrix is invertible, avoiding singular problems;

[0163] Based on the optimized structure, the same name area of the point cloud data is fixed, including:

[0164] The maximum Euclidean distance between the source point cloud and the target point cloud is calculated:

[0165]

[0166] Threshold setting: screening points satisfying d i ≤ε, wherein the threshold

[0167] The KDTree structure of the Open3D library is used to match the point cloud.

[0168] Compared with the prior art, the application has the following beneficial effects:

[0169] 1. Efficiently remove global noise:

[0170] Based on the Z-score dynamic threshold method, combined with the statistical characteristics of three-dimensional Euclidean distance, the accuracy of abnormal point removal is improved to 95%, while the integrity of the main distribution point cloud is preserved;

[0171] The time complexity is effectively optimized to support million-level point cloud processing per second.

[0172] 2. High-precision plane calibration:

[0173] The regional adaptive weighted mean algorithm dynamically adjusts the neighborhood radius according to the target area, and the plane geometric center error is stable within 0.1mm;

[0174] The Newton iteration method fuses the Hessian matrix to optimize the spherical parameter, the spherical center error is ≤0.5mm, the radius error is ≤0.05mm, and the accuracy is improved by 80% compared with the traditional method.

[0175] 3. Multi-dimensional space theory innovation:

[0176] A multi-level calibration framework of four-parameter method+NDT is proposed, which realizes the joint optimization of coordinate system conversion and point cloud registration, and reduces the global geometric error by 73.6%;

[0177] The dynamic weight function (ω=1 / (1+d 2 )) is introduced to suppress the interference of outliers on the geometric center calculation.

[0178] 4. Standardization and universality:

[0179] The target design conforms to the specifications of the highway transportation industry, and is suitable for multiple types of geometric models such as spherical and planar;

[0180] The algorithm is compatible with the geometric measurement of infrastructure such as bridges and tunnels, and has strong expansibility. BRIEF DESCRIPTION OF DRAWINGS

[0181] In order to more clearly illustrate the technical solutions of the present application, the following briefly introduces the drawings needed in the embodiments. Obviously, the drawings described in the following are only some embodiments of the present application, and other drawings can also be obtained by those skilled in the art without creative labor.

[0182] Figure 1 The flowchart of the highway traffic structure three-dimensional geometric quantity multi-layer calibration method based on target reference of the embodiment of the present application;

[0183] Figure 2 Figure 1 is a schematic diagram of a multi-layer calibration system structure for target reference-based three-dimensional geometric quantity of a highway traffic structure according to an embodiment of the present application. DETAILED DESCRIPTION

[0184] The technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of protection of the present application.

[0185] In order to make the above-mentioned purposes, features and advantages of the present application more obvious and easy to understand, the present application will be further described in detail below with reference to the drawings and specific embodiments.

[0186] Embodiment one

[0187] As shown in the figure, the present embodiment provides a target multi-layer calibration method for three-dimensional geometric quantity measurement of a highway traffic structure, including the following steps: Figure 1 Step one: point cloud data preprocessing. A point cloud data denoising method for three-dimensional laser scanner acquisition is proposed, and the data quality is improved by Z-score method;

[0188] Step two: point cloud data plane error calibration. Based on the denoised point cloud data, the plane error of the three-dimensional laser scanning device is calibrated by using the standard plane target, the geometric center recognition accuracy is improved and the point cloud data plane error is reduced by using the region adaptive improved weighted mean algorithm;

[0189] Step three: point cloud data three-dimensional error calibration. In order to solve the problem of insufficient three-dimensional geometric quantity measurement accuracy of the scanning device, the three-dimensional error of the three-dimensional laser scanning device is calibrated by using the standard spherical target, and the data accuracy of the three-dimensional point cloud of the device is improved by using Newton iteration algorithm to iterate the parameter vector multiple times.

[0190] Step four: multi-station acquisition device point cloud data calibration. In order to solve the problem of data mismatch and low precision of multi-station scanning device, the coordinates of different stations are converted by four-parameter method, the same name area registration is carried out by using the maximum geometric error, the point cloud data geometric error is calibrated by using NDT optimization method, and the point cloud data precision is improved by algorithm compensation.

[0191] In the present embodiment, the point cloud data preprocessing: global point cloud random noise preprocessing based on Z-score: a point cloud data denoising method for three-dimensional laser scanner acquisition is proposed, and the abnormal points and noise points are eliminated by Z-score method to improve the data quality.

[0192] ​

[0193] 1 Three-dimensional point cloud coordinate conversion

[0194] Convert the three-dimensional point cloud coordinates to the Euclidean distance to the origin, and construct a one-dimensional data set:

[0195]

[0196] 2 Statistical parameter calculation

[0197] Mean:

[0198] Standard deviation:

[0199] Z-score calculation:

[0200] 3 Dynamic threshold setting

[0201] Initial threshold: According to the 3σ principle of normal distribution, set the initial threshold to 3.

[0202] Dynamic adjustment: Gradually reduce the threshold according to the actual removal effect (considering the distribution characteristics of the data set and the actual application status of the highway transportation industry, generally take 2.5, and dynamically adjust downward according to the removal effect of data outliers), until the proportion of outliers is stable.

[0203] 4 Outlier removal

[0204] Remove points that satisfy ∣Z i ∣> threshold, and keep the main distribution data, where Z i is the Z-score value of each point.

[0205] 5 Example

[0206] Step 1: Input the simulated point cloud data containing noise

[0207] The coordinate values of the simulated 6 points are as follows:

[0208] {(1,2,3),(2,3,3),(3,4,3),(4,5,3),(5,6,3),(15,15,15)}

[0209] Step 2: Calculate the Euclidean distance of each point to the origin

[0210] d i = {3.74, 4.69, 5.83, 7.07, 8.37, 25.98}

[0211] Step 3: Calculate the mean and standard deviation

[0212] μ = 9.28, σ = 7.61

[0213] Step 4: Calculate Z-score and remove outliers

[0214] Z i = {-0.73, -0.60, -0.45, -0.29, -0.12, +2.19}

[0215] Set threshold 2.0, point 6 (Z6 = 2.14) is determined as an outlier and removed.

[0216] 6 Technical effects

[0217] Global denoising: suitable for random noise (such as sensor speckle), retaining the main distribution point cloud.

[0218] Dynamic adaptability: threshold can be dynamically optimized according to data distribution, avoiding excessive deletion of effective data.

[0219] Computational efficiency: low time complexity, supporting large-scale point cloud real-time processing.

[0220] In this embodiment, the plane error calibration of point cloud data: based on the region adaptive improved weighted mean plane error calibration method: based on the denoising pretreated point cloud data, the plane error of three-dimensional laser scanning equipment is calibrated by using the standard plane target, the region adaptive improved weighted mean algorithm is used to improve the geometric center recognition accuracy and reduce the plane error of point cloud data.

[0221] 1 Region classification

[0222] According to the actual application situation of common structure to be detected in highway industry, reference is made to "JGJ8 Building Deformation Measurement Specification",

[0223] "DB42 / T 1496-2019 Highway Slope Monitoring Technical Specification", "JT / T 1037-2022 Highway Bridge Structure Monitoring Technical Specification", "JTG F80 Highway Engineering Quality Inspection and Evaluation Standard", "GB / T 39410-2020 Low-Orbit Spaceborne GNSS Measurement Receiver General Specification", etc., the geometric size of structure monitoring refers to JTGD60, JTG3362, JTG / TD65-05, JTG / T3365-01, JTG / TD65-06, JTG / T3360-01, JTG / TD65-05, JTG5120-2021, JTG / TH21, etc. The area S of the region to be detected is divided into four types of small, medium, large and extra-large, and the specific classification standards are as follows: u

[0224]

[0225] The region classification provides a basis for the adaptive adjustment of the subsequent neighborhood radius and point cloud density.​

[0226] 2Point cloud density distribution and weighted neighborhood radius calculation

[0227] Define the point cloud density distribution function D as the judgment basis of neighborhood radius r, then:

[0228]

[0229] In the formula, N r is the actual number of sampling points in the region to be measured, which is determined according to the field measurement requirements;

[0230] N tot is the theoretical number of sampling points;

[0231] λ is an empirical parameter of point cloud density, and for three-dimensional scanning of highway infrastructure structures, the recommended value range is 0.79-0.83.

[0232] Combined with regional classification, the calculation formula of neighborhood radius is:

[0233]

[0234] 3Calculation of initial value of geometric center

[0235] A neighborhood is constructed with neighborhood radius r, containing N point cloud data, and the initial geometric center coordinates C0=(x0, y0) are calculated:

[0236]

[0237] In the formula, N i is the number of sampling points in the i-th sub-region in the neighborhood, and n is the number of sub-regions.

[0238] 4Weight calculation and iterative optimization

[0239] Based on the distance d i of point cloud data to the temporary geometric center, the weight function ω i is defined:

[0240]

[0241] In the formula, is the temporary geometric center coordinates of the k-th iteration.

[0242] The geometric center coordinates C=(x, y) are iteratively optimized by the weighted mean algorithm:

[0243]

[0244] The iteration termination condition is that the absolute value of the difference between the coordinates of the adjacent two times is less than the threshold value ∈ (recommended ∈=0.1mm).

[0245] 5 Embodiment

[0246] Step 1: Design of standard target

[0247] A circular planar target with a standard area of 300 square centimeters and a standard geometric center coordinate of (0 mm, 0 mm) is designed, and the target surface reflectivity needs to meet the specification requirements of “JT / T 1037-2022”.

[0248] When only using the traditional scanning method without plane error calibration, the measured target geometric center positioning error is:

[0249]

[0250] The coordinate positioning error is:

[0251]

[0252] Step 2: Neighborhood radius calculation

[0253] Area type: medium (S u = 300 cm 2 )

[0254] Parameter selection: λ = 0.81

[0255] Point cloud density distribution function:

[0256] Optimal neighborhood radius:

[0257] Step 3: Calculation of initial value of geometric center

[0258]

[0259] Step 4: Iterative optimization calibration of plane geometric center error

[0260] The weighted mean algorithm is used to calibrate the plane geometric center error, and the iteration limit is set to 20 and the iteration threshold is set to 0.1 mm. The core code segment is as follows:

[0261]

[0262]

[0263] In this embodiment, the final iteration number is 5, and the detailed iteration process is as follows:

[0264] Iteration 1: X = 0.2512 mm, Y = 0.2476 mm,

[0265] Iteration 2: X = 0.1589 mm, Y = 0.1534 mm

[0266] Iteration 3: X = 0.1135 mm, Y = 0.1128 mm

[0267] Iteration 4: X = 0.1041 mm, Y = 0.1023 mm,

[0268] Iteration 5: X = 0.1023 mm, Y = 0.0987 mm

[0269] Final coordinates: (0.1023 mm, 0.0987 mm)

[0270] Final coordinate error: 0.1421 mm, achieving sub-millimeter level precision calibration, meeting the precision requirements of the highway transportation industry.

[0271] 6 Technical effects

[0272] Region adaptation: dynamically adjust the neighborhood radius to adapt to the point cloud distribution characteristics of different area regions;

[0273] Standardized implementation: design targets combined with industry specification requirements to ensure the repeatability and engineering applicability of calibration results;

[0274] Efficiency optimization: fast convergence speed, precision can reach millimeter level.

[0275] In this embodiment, the point cloud data stereo error calibration: the laser scanning device stereo error calibration method based on Newton iteration algorithm: aiming at the problem of insufficient three-dimensional feature detection accuracy of laser scanning device, the stereo error of the device is calibrated by using standard spherical target, and the parameter vector is iterated multiple times by Newton iteration algorithm to improve the data precision of three-dimensional point cloud of the device. Before actual application, the stereo error of the three-dimensional scanning device is calibrated by using the spherical target to improve the device precision.

[0276] 1. Spherical model establishment

[0277] Define the spherical equation as:

[0278] (x-a) 2 +(y-b) 2 +(z-c) 2 =r 2

[0279] Where (x, y, z) is the coordinate of a point on the defined sphere, (a, b, c) is the three-dimensional coordinate of the sphere center, and r is the radius.

[0280] 2. Loss function construction

[0281] Minimize the sum of the squared distances of all point clouds to the fitted sphere, and the loss function is defined as:

[0282]

[0283] 3 Gradient and Hessian matrix calculation

[0284] Gradient vector:

[0285]

[0286] Hessian matrix:

[0287]

[0288] Take for example

[0289]

[0290] The remaining elements are calculated according to the second-order partial derivative rule to construct a 4x4 symmetric matrix.

[0291] 4 Parameter iterative update

[0292] Parameter vector The update formula is:

[0293]

[0294] In the formula, is the updated parameter vector, is the parameter vector before updating, H is the Hessian matrix, is the gradient vector.

[0295] The iteration termination condition is: radius error Δr≤0.003m or the maximum number of iterations is reached.

[0296] 5 Example

[0297] Step 1: Experimental setup

[0298] Target specifications: as standard spherical target, theoretical radius r=75mm.

[0299] Data acquisition: use a three-dimensional laser scanner to obtain target point cloud data.

[0300] Initial parameters: initial values are calculated by taking the mean of point cloud coordinates, and the median of the distance from the point to the centroid is selected to improve robustness.

[0301] Step 2: Algorithm implementation

[0302] The Python core code is as follows:

[0303]

[0304]

[0305] Step 3: Optimization parameters and error calculation

[0306] After iterative calculation optimization,

[0307] Geometric center and error:

[0308] P0(x0, y0, z0) = (0.4134, 0.1856, 0.4432) mm = φ(P0)

[0309] Radius and error:

[0310] r = 74.9953 mm; φ(r) = 0.0047 mm

[0311] Traditional error is as follows:

[0312] P0(x0, y0, z0) = (6.6385, 3.3251, 2.2295) mm = φ(P0)

[0313] r = 78.5553 mm; φ(r) = 3.5553 mm

[0314] 6 Technical effects

[0315] High-precision error compensation: By using Newton's iterative method to optimize the geometric center and radius of the spherical target, the final geometric center error is reduced to (0.4134, 0.1856, 0.4432), and the radius error is 0.0047 mm, which is a significant improvement in precision compared to the traditional least squares method (error 3.5553 mm).

[0316] Fast convergence: Newton's method uses second-order derivative information (Hessian matrix) to accelerate convergence, and the average error threshold can be reached within 100 iterations, significantly improving calibration efficiency.

[0317] Enhanced anti-interference capability: By fusing the spatial distribution characteristics of three-dimensional point cloud, the robustness of the algorithm to environmental and edge point interference is significantly better than traditional methods, ensuring the reliability of three-dimensional feature measurement such as road crack depth and unevenness.

[0318] Practical application advantages: After optimization, the accuracy of device point cloud data is improved, which can directly reduce the calculation deviation of three-dimensional geometric quantities, providing high-reliability calibration data support for road industry geometric measurement.

[0319] In this embodiment, multi-station acquisition device point cloud data calibration: To solve the problem of mismatching of point cloud data and low accuracy of point cloud data of three-dimensional scanning equipment in different stations, four-parameter method is used to convert different station coordinates, maximum geometric error is used for same name area registration, NDT optimization method is used to calibrate point cloud data geometric error, and point cloud data accuracy is improved through algorithm compensation.

[0320] 1 Four-parameter coordinate calibration of point cloud data

[0321] In multi-station scanning, point clouds in different coordinate systems need to be unified to a global coordinate system. Suppose the homonymous point pairs of the source point cloud and the target point cloud are: source point cloud coordinates: P t,i =(x i ,y i ), target point cloud coordinates: P s,i =(x s,i ,y s,i ). The plane coordinate calibration is performed by the four-parameter method, and the formula is as follows:

[0322]

[0323] Wherein, a is the scale transformation parameter, θ is the rotation angle, t x and t y are the translation amounts.

[0324] Four-parameter calibration steps:

[0325] ①Centroid decentralization: calculate the centroids of the source point cloud and the target point cloud and decentralize the processing.

[0326] ②Calculate intermediate variables

[0327] C1=∑(x t x s +y t y s ); C2=∑(y t x s -x t y s )

[0328] C3=∑x s ; C4=∑y s ; C5=∑x t ; C6=∑y t

[0329] ③Solve the rotation angle:

[0330] θ=arctan2(C2,C1)

[0331] ④Calculate the translation amount:

[0332] t x =(C5-C3cosθ+C4sinθ) / N

[0333] t y =(C6-C3sinθ-C4cosθ) / N

[0334] ⑤Longitudinal axis calibration:

[0335] By Compensate for elevation error.

[0336] 2 NDT calibration optimization of point cloud data

[0337] To improve calibration accuracy and efficiency, the NDT (Normal Distributions Transform) method is used to optimize the initial calibration results:

[0338] Mesh partitioning: divide the target point cloud into fixed-size cubic meshes

[0339] Probability density modeling: calculate the mean μ and covariance matrix Γ of the point cloud in each mesh:

[0340]

[0341] Optimization goal: maximize the overall likelihood function

[0342]

[0343] Iterative optimization: update the transformation matrix through the Levenberg-Marquardt algorithm until the difference between adjacent two transformations is less than the threshold

[0344]

[0345] In the formula, M k+1 and M k represent the transformation parameters of the kth and k+1th iterations, represents the sum of the gradients of the log-likelihood function of all mesh units or data points, is the sum of the second derivatives of the log-likelihood function, λ is the damping factor of the algorithm, used to balance the gradient descent and Gauss-Newton method, when λ is larger, the algorithm is close to gradient descent, when λ is smaller, the algorithm is close to Gauss-Newton method, I is the identity matrix, combined with λ to ensure that the matrix is invertible, avoiding singular problems.

[0346] 3 Fixing of the same name region of point cloud data

[0347] To reduce manual intervention, an automatic calibration strategy for the same name region based on the maximum geometric error is proposed:

[0348] Distance calculation: calculate the maximum Euclidean distance between the source point cloud and the target point cloud

[0349]

[0350] Threshold setting: filter points that satisfy d i ≤ε, where the threshold

[0351] KDTree acceleration: Use the KDTree structure of the Open3D library to quickly match point clouds and improve computing efficiency.

[0352] 4 Engineering implementation site application

[0353] Select a standard road section of a test field in Beijing for multi-site scanning test:

[0354] Device configuration: FARO Focus X Series long-range laser scanner, laser wavelength 1550nm, ambient temperature 16.8°C.

[0355] Data acquisition: Lay out 4 scanning stations to cover a test space of 10m x 10m x 3m and obtain raw point cloud data.

[0356] Preprocessing: Remove noise points by Z-score method.

[0357] Plane calibration: Combine with standard plane targets and use regional adaptive weighted mean algorithm to calibrate the plane error of the device.

[0358] Stereoscopic calibration: Combine with standard spherical targets and use Newton iteration method to calibrate the stereoscopic error of the device.

[0359] Multi-site device calibration: Calibrate the coordinate system by four-parameter method for multi-site device data coarse registration; optimize point cloud alignment by NDT for multi-site device data fine registration.

[0360] 5 Examples

[0361] Step 1: Test setup

[0362] Test object: Use laser scanning equipment to scan a crack disease on the test road section

[0363] Parameter comparison: Use traditional method to measure the geometric parameters of the crack before calibration, use high-level standard instrument (Renishaw XL-80 laser interferometer, China Institute of Metrology Research Institute calibration CDlx2024-01568) to calculate the measurement error of the traditional method, the area calculation error is 0.0181m 2 , the perimeter calculation error is 18.1805mm, and the depth calculation error is 5.3902mm.

[0364] NDT parameter setting: Voxel size is set to 0.01, maximum iteration number is set to 100, Open3D parallel acceleration is enabled.

[0365] Step 2: Algorithm implementation

[0366] Python core code as follows:

[0367]

[0368] Step 3: Comparison of measured parameters

[0369] After data preprocessing, plane error calibration, stereo error calibration and multi-station error calibration, the pit slot compensation topography is obtained.

[0370] The area calculation error is 0.0031m 2 (Compared with the optimization before calibration, 82.91%), the perimeter calculation error is 1.2326mm (Compared with the optimization before calibration, 93.22%), and the depth calculation error is 0.7725mm (Compared with the optimization before calibration, 85.67%).

[0371] 6Technical effects

[0372] Precision verification: After calibration, the data meet the detection precision requirements of the geometric quantity of the highway transportation industry, and the error of feature geometric quantity detection is reduced.

[0373] Efficiency improvement: The multi-station registration time (about 8min) is shortened to 60% of the traditional method.

[0374] In summary, 1. Multi-level error calibration system based on target:

[0375] The plane calibration adopts the region adaptive improved weighted mean algorithm; the stereo calibration adopts the Newton iteration method; the multi-station calibration fuses the four-parameter method (coarse registration) and NDT (fine registration), to realize global point cloud alignment.

[0376] 2. Automatic data processing flow:

[0377] Noise removal adopts Z-score method; the same name area is automatically screened through dynamic threshold .

[0378] 3. Standardized engineering implementation scheme:

[0379] The target design conforms to the specifications of the highway transportation industry, and is suitable for geometric quantity calibration of various types of infrastructure structures and multi-type geometric models such as spherical surface and plane.

[0380] Example two

[0381] As Figure 2 shown, the application also provides a target multi-layer calibration system for highway transportation structure three-dimensional geometric quantity measurement, which is used to realize the method, and the system comprises a preprocessing module, a plane error calibration module, a stereo error calibration module and a region matching module.

[0382] The preprocessing module is used to collect point cloud data of highway infrastructure structures using a 3D laser scanner and perform noise reduction preprocessing.

[0383] The plane error calibration module is used to calibrate the plane error of the 3D laser scanning equipment based on the denoised preprocessed point cloud data and using a standard planar target.

[0384] The stereo error calibration module is used to calibrate the stereo error of the 3D laser scanning equipment based on the denoised preprocessed point cloud data and using a standard spherical target.

[0385] The region matching module is used to perform registration of the same region based on the calibration results of the planar and stereo errors of the 3D laser scanning equipment. It transforms the coordinates of different stations using the four-parameter method and uses the maximum geometric error.

[0386] In this embodiment, point cloud data of highway infrastructure structures is acquired using a 3D laser scanner and subjected to noise reduction preprocessing, including:

[0387] Transform the 3D point cloud coordinates to the Euclidean distance from the origin to construct a one-dimensional dataset;

[0388] Calculate the mean and Z-score of a one-dimensional dataset. i ;

[0389] Based on the Z-score calculation results, set a dynamic threshold;

[0390] Eliminate those that satisfy |Z i |> Threshold point.

[0391] In this embodiment, based on the denoised point cloud data, the planar error of the 3D laser scanning equipment is calibrated using a standard planar target, including:

[0392] Obtain the area of ​​the region to be measured;

[0393] Calculate the point cloud density distribution, and calculate the neighborhood radius based on the area of ​​the region to be measured and the point cloud density distribution;

[0394] Calculate the initial value of the geometric center based on the neighborhood radius;

[0395] Based on the initial value of the geometric center, the coordinates of the geometric center are iteratively optimized through a weighted average algorithm to complete the calibration of the planar error of the 3D laser scanning equipment;

[0396] Obtaining the area of ​​the region to be measured includes:

[0397]

[0398] Calculating point cloud density distribution includes:

[0399]

[0400] In the formula, N r is the actual number of sampling points of the to-be-measured area, which is determined according to field measurement requirements; N tot is the theoretical number of sampling points; λ is an empirical parameter of point cloud density; r is a neighborhood radius; and D is a point cloud density distribution function;

[0401] According to the area of the to-be-measured area and the point cloud density distribution, a neighborhood radius is calculated, including:

[0402]

[0403] According to the neighborhood radius, an initial value of a geometric center is calculated, including:

[0404] A neighborhood is constructed with the neighborhood radius r, containing N point cloud data, and an initial geometric center coordinate C0=(x0, y0) is calculated:

[0405]

[0406] In the formula, N i is the number of sampling points of the i-th sub-area in the neighborhood, and n is the number of sub-areas;

[0407] Based on the initial value of the geometric center, the geometric center coordinate is iteratively optimized by a weighted mean algorithm to complete the planar error calibration of the point cloud data, including:

[0408] Based on the distance d i of the point cloud data to the temporary geometric center, a weight function ω i is defined:

[0409]

[0410] In the formula, is the temporary geometric center coordinate of the k-th iteration;

[0411] The geometric center coordinate C=(x, y) is iteratively optimized by a weighted mean algorithm:

[0412]

[0413] The iteration termination condition is that the absolute value of the difference between two adjacent coordinates is less than a threshold ∈.

[0414] In this embodiment, based on the point cloud data after denoising preprocessing, the stereoscopic error of a three-dimensional laser scanning device is calibrated by using a standard spherical target, including:

[0415] A spherical model is established;

[0416] Based on the spherical model, a loss function is constructed;

[0417] Based on the loss function, a gradient and a Hessian matrix are calculated;

[0418] Based on the gradient and the Hessian matrix, a parameter iterative update is performed to complete the calibration of the three-dimensional laser scanning device;

[0419] Wherein, the spherical model is established, including:

[0420] The preset spherical equation is:

[0421] (x-a) 2 +(y-b) 2 +(z-c) 2 =r 2 ;

[0422] Wherein, (x, y, z) is the coordinate of the point on the preset sphere, (a, b, c) is the three-dimensional coordinate of the sphere center, and r is the radius;

[0423] Based on the spherical model, a loss function is constructed, including:

[0424] The loss function is:

[0425]

[0426] Based on the loss function, a gradient and a Hessian matrix are calculated, including:

[0427] The gradient vector is:

[0428]

[0429] The Hessian matrix is:

[0430]

[0431] Based on the gradient and the Hessian matrix, a parameter iterative update is performed, including:

[0432] The update formula of the parameter vector is:

[0433]

[0434] In the formula, is the updated parameter vector, is the parameter vector before updating, H is the Hessian matrix, and is the gradient vector.

[0435] The iteration termination condition is that the radius error Δr≤0.003m or the maximum number of iterations is reached.

[0436] In the embodiment, based on the calibration results of the planar error and the stereo error of the three-dimensional laser scanning device, the different station coordinates are converted by the four-parameter method, and the same name area registration is performed by using the maximum geometric error, including:

[0437] Calibrating the four-parameter coordinates of the point cloud data;

[0438] Optimizing the calibration results by using the NDT method;

[0439] Based on the optimized structure, the same name area of the point cloud data is fixed;

[0440] Wherein, the calibration of the four-parameter coordinates of the point cloud data includes:

[0441] The same point pairs of the source point cloud and the target point cloud are preset as: the source point cloud coordinates: P t,i =(x i ,y i ), the target point cloud coordinates: P s,i =(x s,i ,y s,i ); the planar coordinate calibration is performed by using the four-parameter method, and the calculation formula is:

[0442]

[0443] Wherein, a is a scale transformation parameter, θ is a rotation angle, t x and t y are translation amounts;

[0444] The NDT method is used to optimize the calibration results, including:

[0445] The target point cloud is divided into a fixed-size cubic grid;

[0446] The mean value μ and the covariance matrix Γ of the point cloud in each grid are calculated, and a probability density model is constructed;

[0447]

[0448] Based on the probability density model, the target is optimized by using the maximum overall likelihood function;

[0449]

[0450] The transformation matrix is updated by using the Levenberg-Marquardt algorithm until the difference between the adjacent two times is less than a threshold value:

[0451]

[0452] In the formula, M k+1 and M kdenote the transformation parameters of the kth and k+1th iteration, denote the sum of the gradient of the log-likelihood function of all grid cells or data points, denote the sum of the second derivative of the log-likelihood function, λ is the damping factor of the algorithm, used to balance the gradient descent and Gauss-Newton method, when λ is larger, the algorithm is close to the gradient descent, when λ is smaller, the algorithm is close to the Gauss-Newton method, I is the identity matrix, combined with λ to ensure that the matrix is invertible, avoiding singular problems.

[0453] Based on the optimized structure, the same name area of the point cloud data is fixed, including:

[0454] Calculate the maximum Euclidean distance between the source point cloud and the target point cloud:

[0455]

[0456] Threshold setting: filter the points that meet d i ≤ ε, where the threshold

[0457] Match the point cloud by using the KDTree structure of the Open3D library.

[0458] The above-described embodiments are only descriptions of the preferred modes of the present application and do not limit the scope of the present application. Without departing from the design spirit of the present application, various modifications and improvements to the technical solutions of the present application made by those skilled in the art shall fall within the protection scope determined by the claims of the present application.

Claims

1. A multi-layer calibration method for three-dimensional geometry of a highway traffic structure based on target reference, characterized in that, The method comprises: acquiring point cloud data of a highway infrastructure structure by using a three-dimensional laser scanner and performing denoising preprocessing; calibrating a plane error of a three-dimensional laser scanning device based on the denoising preprocessed point cloud data by using a standard plane target; calibrating a three-dimensional error of the three-dimensional laser scanning device based on the denoising preprocessed point cloud data by using a standard spherical target; based on the calibration results of the plane error and the three-dimensional error of the three-dimensional laser scanning device, converting different station coordinates by using a four-parameter method, and performing homonymic area registration by using a maximum geometric error; calibrating the plane error of the three-dimensional laser scanning device based on the denoising preprocessed point cloud data by using the standard plane target, comprising: acquiring an area of a to-be-measured region; calculating a point cloud density distribution, and calculating a neighborhood radius according to the area of the to-be-measured region and the point cloud density distribution; calculating an initial value of a geometric center according to the neighborhood radius; based on the initial value of the geometric center, iteratively optimizing a geometric center coordinate by using a weighted mean value algorithm to complete the calibration of the plane error of the three-dimensional laser scanning device; wherein acquiring the area of the to-be-measured region comprises: ； calculating the point cloud density distribution comprises: ； In the formula, is the actual number of sampling points in the region to be measured, which is determined according to the field measurement requirement; is the theoretical number of sampling points; is an empirical parameter of point cloud density; r is the neighborhood radius; D is the point cloud density distribution function; calculating the neighborhood radius according to the area of the to-be-measured region and the point cloud density distribution comprises: ； calculating the initial value of the geometric center according to the neighborhood radius comprises: with a neighborhood radius r constructing a neighborhood, comprising N point cloud data, calculating an initial geometric center coordinate C0= (x0, y0): ； In the formula, N i is the number of sampling points in the i-th sub-region in the neighborhood, i is the number of sampling points in the i-th sub-region in the neighborhood, is the number of sub-regions. based on the initial value of the geometric center, iteratively optimizing the geometric center coordinate by using the weighted mean value algorithm to complete the plane error calibration of the point cloud data comprises: Distance from point cloud data to temporary geometric center d i , defining a weight function ω i : ； In the formula, ( , ) is the first k The temporary geometric center coordinates for the next iteration; Iterative optimization of geometric center coordinates by weighted mean algorithm : ； The iteration termination condition is that the absolute value of the difference between two adjacent coordinates is less than a threshold value .

2. The method of claim 1, wherein, acquiring point cloud data of a highway infrastructure structure by using a three-dimensional laser scanner and performing denoising preprocessing, comprising: converting three-dimensional point cloud coordinates to the Euclidean distance of the origin, and constructing a one-dimensional data set; Computing the mean, Z-score, i.e. Z i ; setting a dynamic threshold according to the calculation result of the Z-score; Z i points satisfying |Z > threshold value.

3. The method of claim 1, wherein, calibrating a three-dimensional error of the three-dimensional laser scanning device based on the denoising preprocessed point cloud data by using a standard spherical target, comprising: establishing a spherical model; based on the spherical model, constructing a loss function; based on the loss function, calculating a gradient and a Hessian matrix; based on the gradient and the Hessian matrix, performing parameter iterative updating to complete the calibration of the three-dimensional error of the three-dimensional laser scanning device; wherein establishing the spherical model comprises: presetting a spherical equation as: ； wherein, is the coordinate of the point on the preset sphere, is the three-dimensional coordinate of the sphere center, r is the radius; based on the spherical model, constructing a loss function, comprising: minimizing the sum of squares of distances from all point clouds to the fitted sphere, and the loss function is: ； based on the loss function, calculating a gradient and a Hessian matrix, comprising: a gradient vector: ； a Hessian matrix: ； based on the gradient and the Hessian matrix, performing parameter iterative updating, comprising: Parameter vector The update formula for is ； wherein is the updated parameter vector, is the updated parameter vector, H is the Hessian matrix, is the gradient vector; The iteration termination condition is a radius error or a maximum number of iterations is reached.

4. The method of claim 3, wherein, based on the calibration results of the plane error and the three-dimensional error of the three-dimensional laser scanning device, converting different station coordinates by using a four-parameter method, and performing homonymic area registration by using a maximum geometric error, comprising: calibrating four-parameter coordinates of the point cloud data; optimizing the calibration result by using an NDT method; fixing the homonymic area of the point cloud data based on the optimized structure; wherein calibrating the four-parameter coordinates of the point cloud data comprises: The homonymic point pair of the preset source point cloud and the target point cloud is: source point cloud coordinate: , target point cloud coordinate: ; plane coordinate calibration is performed through a four-parameter method, and the calculation formula is: ； where a is a scale transformation parameter, θ is a rotation angle, t x and t y are translation amounts; optimizing the calibration result by using the NDT method, comprising: dividing the target point cloud into a fixed-size cubic grid; Compute the mean of the point cloud within each grid and covariance matrix , construct a probability density model; ； Based on the probability density model, the target is optimized by maximizing the overall likelihood function; ； The transformation matrix is updated by the Levenberg-Marquardt algorithm until the difference between adjacent two transformations is less than a threshold value: ； wherein and denote the transformation parameters of the kth and k+1th iteration, denotes the sum of the gradients of the log-likelihood function for all grid cells or data points, is the sum of the second derivatives of the log-likelihood function, λ is a damping factor of the algorithm, used to balance the gradient descent and Gauss-Newton method, when λ is larger, the algorithm is closer to the gradient descent, when λ is smaller, the algorithm is closer to the Gauss-Newton method, I is an identity matrix, combined with λ, ensures that the matrix is invertible, avoiding singular problems; Based on the optimized structure, the same name area of the point cloud data is fixed, including: Calculate the maximum Euclidean distance between the source point cloud and the target point cloud: ； Threshold setting: screening points that meet threshold value ; Match the point cloud by using the KDTree structure of the Open3D library.

5. A multi-layer calibration system for three-dimensional geometry of a road traffic structure based on target reference, the system being for implementing the method according to any one of claims 1 to 4, characterized in that, The system comprises a preprocessing module, a plane error calibration module, a stereo error calibration module, and a region matching module. The preprocessing module is configured to collect point cloud data of a highway infrastructure structure by using a three-dimensional laser scanner and perform denoising preprocessing. The plane error calibration module is configured to calibrate the plane error of the three-dimensional laser scanning device based on the denoised and preprocessed point cloud data by using a standard plane target. The stereo error calibration module is configured to calibrate the stereo error of the three-dimensional laser scanning device based on the denoised and preprocessed point cloud data by using a standard spherical target. The region matching module is configured to convert different site coordinates by a four-parameter method based on the calibration results of the plane error and the stereo error of the three-dimensional laser scanning device, and perform same-name region registration by using the maximum geometric error.

6. The system of claim 5, wherein, The point cloud data of the highway infrastructure structure is collected by using a three-dimensional laser scanner and denoising preprocessing, including: Convert the three-dimensional point cloud coordinates to the Euclidean distance of the origin, and construct a one-dimensional data set; Computing the mean, Z-score, i.e. Z i ; According to the calculation result of the Z-score, a dynamic threshold is set; Z i points satisfying |Z i threshold.

7. The system of claim 5, wherein, Based on the denoised and preprocessed point cloud data, the plane error of the three-dimensional laser scanning device is calibrated by using a standard plane target, including: Obtain the area of the to-be-tested region; Calculate the point cloud density distribution, and calculate the neighborhood radius based on the area of the to-be-tested region and the point cloud density distribution; According to the neighborhood radius, calculate the initial value of the geometric center; Based on the initial value of the geometric center, the geometric center coordinates are iteratively optimized by a weighted mean algorithm to complete the calibration of the plane error of the three-dimensional laser scanning device; The area of the to-be-tested region is obtained, including: ； The point cloud density distribution is calculated, including: ； In the formula, is the actual number of sampling points of the region to be measured, which is determined according to the field measurement requirement; is the theoretical number of sampling points; is an empirical parameter of the point cloud density; r is the neighborhood radius; D is the point cloud density distribution function; The neighborhood radius is calculated based on the area of the to-be-tested region and the point cloud density distribution, including: ； The initial value of the geometric center is calculated based on the neighborhood radius, including: with a neighborhood radius r constructing a neighborhood, comprising N point cloud data, calculating an initial geometric center coordinate C0= (x0, y0): ； In the formula, N i For the neighboring region i Number of sampling points in each sub-region The number of sub-regions; Based on the initial value of the geometric center, the geometric center coordinates are iteratively optimized by a weighted mean algorithm to complete the plane error calibration of the point cloud data, including: Distance from point cloud data to temporary geometric center d i , defining a weight function ω i : ； In the formula, ( , ) is the first k The temporary geometric center coordinates for the next iteration; Iterative optimization of geometric center coordinates by weighted mean algorithm : ； The iteration termination condition is that the absolute value of the difference between two adjacent coordinates is less than a threshold value .

8. The system of claim 7, wherein, Based on the denoised and preprocessed point cloud data, the stereo error of the three-dimensional laser scanning device is calibrated by using a standard spherical target, including: Establish a spherical model; Based on the spherical model, a loss function is constructed; Based on the loss function, a gradient and a Hessian matrix are calculated; Based on the gradient and the Hessian matrix, parameter iterative updating is performed to complete the calibration of the stereo error of the three-dimensional laser scanning device; The spherical equation is preset as: Based on the spherical model, a loss function is constructed, including: ； wherein, is the coordinate of the point on the preset sphere, is the three-dimensional coordinate of the sphere center, r is the radius; Minimize the sum of the squared distances of all point clouds to the fitted sphere, and the loss function is: Based on the loss function, a gradient and a Hessian matrix are calculated, including: ； Gradient vector: ​ ； Hessian matrix: ； Based on the gradient and Hessian matrix, the parameter iterative update includes: Parameter vector The update formula for is ； wherein is the updated parameter vector, is the updated parameter vector, is the gradient vector; The iteration termination condition is a radius error or a maximum number of iterations is reached.

9. The system of claim 8, wherein, Based on the calibration results of the planar error and the stereo error of the three-dimensional laser scanning device, the coordinates of different stations are converted by the four-parameter method, and the same name area registration is performed by using the maximum geometric error, including: Calibration of four-parameter coordinates of point cloud data; Optimize the calibration results by NDT method; Based on the optimization structure, fix the same name area of the point cloud data; Wherein, the calibration of the four-parameter coordinates of the point cloud data includes: The homonymic point pair of the preset source point cloud and the target point cloud is: source point cloud coordinate: , target point cloud coordinate: ; plane coordinate calibration is performed through a four-parameter method, and the calculation formula is: ； where a is a scale transform parameter, θ is a rotation angle, t x and t y are translation amounts; Optimize the calibration results by NDT method, including: Divide the target point cloud into fixed-size cubic grids; computing the mean of the point cloud within each grid and the covariance matrix , constructing a probability density model; ； Based on the probability density model, optimize the target by maximizing the overall likelihood function; ； Update the transformation matrix by Levenberg-Marquardt algorithm until the difference between the adjacent two transformations is less than the threshold: ； wherein and denote the transformation parameters of the kth and k+1th iteration, denotes the sum of the gradients of the log-likelihood function for all grid cells or data points, is the sum of the second derivatives of the log-likelihood function, λ is a damping factor of the algorithm, used to balance the gradient descent and Gauss-Newton method, when λ is larger, the algorithm is closer to the gradient descent, when λ is smaller, the algorithm is closer to the Gauss-Newton method, I is an identity matrix, combined with λ, ensures the matrix is invertible, avoiding singular problem; Based on the optimization structure, fix the same name area of the point cloud data, including: Calculate the maximum Euclidean distance between the source point cloud and the target point cloud: ； Threshold setting: screening points that meet , wherein the threshold ; Match the point cloud by using the KDTree structure of the Open3D library.

Images