An autonomous calibration method for terrestrial laser scanner based on intensity characteristics

Abstract

This invention discloses an intensity-based autonomous calibration method for terrestrial laser scanners. This method requires no additional deployment costs and can simultaneously calibrate both bifacially sensitive and non-bifacially sensitive parameters in the calibration model, achieving millimeter / submillimeter and arc-second accuracy for offset and angle parameters, respectively. The method comprises Step S1: acquiring point cloud data at two measuring stations and roughly calculating the relative pose between the measuring stations using a point cloud registration method; Step S2: uniformly selecting a number of core points from the point cloud and generating an intensity image based on their nearby points; Step S3: extracting and matching feature points from the intensity image; and Step S4: constructing a system of linear equations to solve the calibration model parameters.

Description

Technical Field

[0001] The present invention relates to the technical field of terrestrial laser scanning, and in particular to an autonomous calibration method for a terrestrial laser scanner based on intensity characteristics. Background Art

[0002] Terrestrial laser scanners (TLS) are efficient and highly accurate tools for 3D data acquisition and modeling. They are widely used in fields such as geodesy, engineering, architecture, and archaeology. Due to the inaccuracy of the scanner's internal optical and mechanical components, manufacturers typically pre-calibrate them. However, over extended periods of operation, the relative positions of the scanner's internal components can shift, introducing systematic errors into the measurement results, exceeding the scanner's nominal accuracy. This presents a serious problem for users with high precision requirements.

[0003] In recent years, researchers have been working on user-side scanner calibration. Installing well-spaced targets at a specific calibration site is a common approach in scientific research. By leveraging the equivalence of target centers observed from different scans, and if scans are taken from multiple stations, calibration model parameters and scanner pose parameters can be estimated simultaneously. Thanks to accurate target center estimation (TCE), target-based scanner calibration can achieve the measurement accuracy specified in the instrument manual.

[0004] In addition, various geometric features can also be used for scanner registration and calibration, such as planes, cylinders, paraboloids, and point features in the surrounding environment. In point feature-based methods, some classic feature operators, such as the Scale-Invariant Feature Transform (SIFT) operator, can be used for feature point matching. These feature operators are usually used for coarse registration of point clouds, and then combined with the Iterative Closest Point (ICP) method to achieve fine registration. In addition, in recent years, many advanced key point extraction algorithms, including some learning-based 3D feature point detection, have also been successfully used for point cloud registration.

[0005] However, target-based methods have high requirements for experimental sites and labor costs, and are relatively cumbersome to implement. The ideal scanner calibration solution is to calibrate the scanner in situ during routine scanning projects without the need for additional manual work. Among the methods based on environmental features, although cylinders and paraboloids have the potential to achieve high-precision scanner calibration, these features are not common in general scanning projects. In addition, some plane-based self-calibration methods can achieve similar accuracy to target-based methods. However, the current general plane-based calibration algorithms require manual extraction of planes or manual layout of planar targets and are not completely autonomous. In addition, although planar features are widely present in urban environments, we cannot rely on planar features when operating in mountainous areas in geographic monitoring applications. Some scholars have studied TLS correction algorithms based on point cloud intensity features, but the current algorithms can only perform feature matching on point clouds from different observation surfaces of the same station, and thus can only calibrate the sensitive parameters of TLS on both sides. Summary of the Invention

[0006] To overcome the shortcomings of the existing technology, the present invention provides an intensity-based autonomous calibration method for terrestrial laser scanners. This method requires no additional deployment costs and can simultaneously calibrate both bifacially sensitive and non-bifacially sensitive parameters in the calibration model. It achieves millimeter / submillimeter and arc-second accuracy for offset and angle parameters, respectively.

[0007] In order to achieve the above object, the present invention provides the following technical solutions:

[0008] A method for autonomously calibrating a terrestrial laser scanner based on intensity characteristics comprises the following steps:

[0009] Step S1: Obtain point cloud data at two measuring stations respectively, and roughly calculate the relative pose between the measuring stations using the point cloud registration method;

[0010] Step S2: uniformly select several core points from the point cloud and generate intensity images based on the points near them;

[0011] Step S3: extracting and matching feature points in the intensity image;

[0012] Step S4: constructing a system of linear equations to solve the calibration model parameters;

[0013] The step S1 comprises:

[0014] Step S11: Two measuring stations are set up at a certain distance from each other. The collected point cloud data are recorded as point cloud 1 and point cloud 2 respectively. The overlap between the two point clouds should be greater than 50%;

[0015] Step S12: using an iterative closest point method or other point cloud registration algorithm to obtain the initial relative pose between point cloud 1 and point cloud 2;

[0016] The step S2 comprises:

[0017] Step S21: Select core points;

[0018] In order to reduce the computational burden while ensuring that all different features can be identified, point cloud 1 is divided into voxels, and the point closest to the center of each voxel will be selected as the core point;

[0019] Step S22: generating an intensity image;

[0020] For each core point, the corresponding intensity image is generated using the points within a certain radius around it; the radius is determined by the following empirical formula:

[0021]

[0022] Among them, ρ1 and ρ2 are the distances from the core point to the two measuring stations, and the value of R is limited to 0.3m to 4.0m. If the number of points in the neighborhood of a core point is less than 100, it will not participate in the subsequent processing. On the contrary, if the number of points in the neighborhood of a core point exceeds 1,000,000, the points around it will be downsampled. p The points in the neighborhood are analyzed using principal component analysis (PCA), and a local coordinate system O-XYZ is constructed with the three main directions as coordinate axes and the core point as the origin. l ; All neighboring points will be projected onto the image plane O-XY l The intensity values ​​are interpolated by the griddata function in Matlab to generate a grid with a size of m. I ×m I The pixel size of the intensity image can be calculated by the following formula:

[0023]

[0024] In this invention, for the observation data with a resolution of 3mm@10m, we take m I = 1001 to obtain sub-pixel matching accuracy (with the increase of resolution, m I can be increased accordingly); in order to obtain the image pairs for subsequent matching, the initial relative pose is used to align the core point C p Perform coordinate transformation to obtain C ′ p , similarly, C ′ p Neighboring points in point cloud 2 will be projected onto the same image plane O-XY l and generate the corresponding intensity image;

[0025] The step S3 comprises:

[0026] Step S31: Utilize The operator locates the salient features in the image, and then extracts and matches the feature points in the image pair based on the ORB (OrientedFAST and Rotated BRIEF) algorithm;

[0027] Step S32: Detect and remove mismatched feature point pairs using an estimated sample consistency algorithm;

[0028] Step S33: Calculate the three-dimensional coordinates corresponding to the two-dimensional feature points;

[0029] According to the pixel coordinates and pixel size of the feature point, its XY l Coordinate, Z l The component is obtained by Z l The coordinates are obtained by linear interpolation; then, the above three-dimensional coordinates are converted to the coordinate system of each station and expressed in polar coordinates. In addition, to further eliminate potential mismatches, the three-dimensional coordinates of the matched feature points are used to estimate the rigid body transformation between point clouds. Feature point pairs with residuals exceeding a threshold (residuals between the same station do not exceed 10 mm, and residuals between different stations do not exceed 200 mm) will be discarded.

[0030] The step S4 comprises:

[0031] Step S41: constructing a conditional equation;

[0032] The following equation can be constructed based on the three-dimensional coordinates of the feature points obtained in step S33:

[0033]

[0034] Where i = 1…I, n = 1…N, I and N are the number of characteristic points and conditional equations respectively. ′ is the station number, R l and T l is the rotation and translation between the station coordinate system and the global coordinate system. and is the coordinate of the feature point after calibration; [f x,n f y,n f z,n ] T is the coordinate difference between a pair of feature points in the global coordinate system. By minimizing it, the calibration parameters and posture parameters of the scanner can be obtained.

[0035] The calibration of the observed values ​​can be performed using the following formula:

[0036]

[0037]

[0038] Where ξ is the calibration parameter, [r i.k θ i,k α i,k ] T is the polar coordinate of the i-th feature point at station k. The calibration model can use the simplified NIST model;

[0039] Step S42: solving the conditional equations;

[0040] First, linearize each equation, and then use the Gauss-Helmert model to estimate the parameters. Assuming that the coordinate components of the observations are uncorrelated, define their covariance matrix as:

[0041]

[0042] The covariance matrix of all observations is then a block diagonal matrix:

[0043] Q f =diag([q1…q u …q U ]) U×U #(7)

[0044] Among them, U is the number of feature points. The initial values ​​are 3.5 mm, 20″, and 20″ respectively. During the iterative solution process, the variance of the observations is updated through variance component estimation (VCE) and F-test, and the Danish method is used to reweight the residuals. At the same time, observations with residuals exceeding three times the mean error will be eliminated.

[0045] The above technical solution is only a feasible technical solution of the present invention. The protection scope of the present invention is not limited thereto. Those skilled in the art can reasonably adjust the specific design according to actual needs.

[0046] The above invention has the following advantages or beneficial effects:

[0047] Existing high-precision TLS calibration methods often rely on artificial targets, which increases the manpower and material costs of calibration. The calibration method based on the existing features of the current environment faces the problems of unavailable geometric features and low calibration accuracy. The present invention makes full use of the intensity features in the environment and can be applied to a variety of scenarios. It does not require additional layout costs and can achieve high-precision calibration parameter calibration. In addition, the present invention can simultaneously calibrate the double-sided sensitive and non-double-sided sensitive parameters in the calibration model. Compared with the existing TLS calibration based on intensity features, this method can achieve more comprehensive TLS parameter calibration. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 It is a schematic diagram of the process flow of the autonomous calibration method of the terrestrial laser scanner of the present invention;

[0049] Figure 2 It is a schematic diagram of feature point extraction and matching of the present invention. DETAILED DESCRIPTION

[0050] The structure of the present invention is further described below with reference to the accompanying drawings and specific embodiments, but is not intended to limit the present invention.

[0051] To address the high-precision calibration of the Leica RTC360, we conducted experiments using data from two stations collected at the 75×33×9m calibration field of the University of Bonn.

[0052] Step S1: Obtain point cloud data at two measuring stations respectively, and roughly calculate the relative pose between the measuring stations using general point cloud registration methods;

[0053] Specifically, step S1 includes:

[0054] Step S11: Two measuring stations are set up at a certain distance from each other, and the collected point cloud data are recorded as point cloud 1 and point cloud 2 respectively, and sufficient overlap should be ensured between the two;

[0055] Step S12: using an iterative closest point method or other point cloud registration algorithm to obtain the initial relative pose between point cloud 1 and point cloud 2;

[0056] Step S2: uniformly select several core points from the point cloud and generate intensity images based on the points near them;

[0057] Specifically, step S2 includes:

[0058] Step S21: Select core points;

[0059] In order to reduce the computational burden while ensuring that all different features can be identified, point cloud 1 is divided into voxels, and the point closest to the center of each voxel will be selected as the core point.

[0060] Step S22: generating an intensity image;

[0061] For each core point, the corresponding intensity image is generated using the points within a certain radius around it. The radius R is determined by the following empirical formula:

[0062]

[0063] Where ρ1 and ρ2 are the distances from the core point to the two stations, and the value of R is limited to 0.3m to 4.0m. If the number of points in the neighborhood of a core point is less than 100, it will not participate in the subsequent processing; on the contrary, if the number of points in the neighborhood of a core point exceeds 1,000,000, the points around it will be downsampled. p The points in the neighborhood are analyzed using principal component analysis (PCA), and a local coordinate system O-XYZ is constructed with the three main directions as coordinate axes and the core point as the origin. l All neighboring points will be projected onto the image plane O-XY l The intensity values ​​are interpolated by the griddata function in Matlab to generate a grid with a size of m. I ×m I The pixel size of the intensity image can be calculated by the following formula:

[0064]

[0065] In this embodiment, m I = 1001. In order to obtain the image pair for subsequent matching, the initial relative pose is used to align the core point C p Perform coordinate transformation to obtain C′ p , similarly, C′ p Neighboring points in point cloud 2 will be projected onto the same image plane O-XY l and generate the corresponding intensity image;

[0066] Step S3: extracting and matching feature points in the intensity image;

[0067] Specifically, step S3 includes:

[0068] Step S31: Utilize The operator locates the salient features in the image, and then extracts and matches the feature points in the image pair based on the ORB (OrientedFAST and Rotated BRIEF) algorithm;

[0069] Step S32: using the estimated sample consistency algorithm (MSAC) algorithm to detect and remove mismatched feature point pairs;

[0070] Step S33: Calculate the three-dimensional coordinates corresponding to the two-dimensional feature points;

[0071] According to the pixel coordinates and pixel size of the feature point, its XY l Coordinate, Z l The component is obtained by Z lThe coordinates are obtained by linear interpolation. These 3D coordinates are then converted to their respective station coordinate systems and expressed in polar coordinates. Furthermore, to further eliminate potential mismatches, the 3D coordinates of the matched feature points are used to estimate the rigid body transformation between the point clouds. Feature point pairs whose residuals exceed a certain threshold are discarded.

[0072] Step S4: constructing a system of linear equations to solve the calibration model parameters;

[0073] Specifically, step S4 includes:

[0074] Step S41: constructing a conditional equation;

[0075] The following equation can be constructed based on the three-dimensional coordinates of the feature points obtained in step S33:

[0076]

[0077] Where i = 1…I, n = 1…N, I and N are the number of characteristic points and conditional equations respectively. ′ is the station number, R l and T l is the rotation and translation between the station coordinate system and the global coordinate system. and is the coordinate of the feature point after calibration. [f x,n f y,n f z,n ] T is the coordinate difference between a pair of feature points in the global coordinate system. By minimizing it, the calibration parameters and posture parameters of the scanner can be obtained.

[0078] The calibration of the observed values ​​can be performed using the following formula:

[0079]

[0080]

[0081] Where ξ is the calibration parameter, [r i.k θ i,k α i,k ] T is the polar coordinate of the i-th feature point at station k. The calibration model can use the simplified NIST model.

[0082] Step S42: solving the conditional equations;

[0083] First, linearize each equation, and then use the Gauss-Helmert model to estimate the parameters. Assuming that the coordinate components of the observations are uncorrelated, define their covariance matrix as:

[0084]

[0085] The covariance matrix of all observations is then a block diagonal matrix:

[0086] Q f =diag([q1…q u …q U ]) U×U #(7)

[0087] Among them, U is the number of feature points. The initial values ​​are 3.5 mm, 20 inches, and 20 inches, respectively. During the iterative solution process, the variance of the observations is updated using variance component estimation (VCE) and F-tests, and the residuals are reweighted using the Danish method. Observations with residuals exceeding three times the mean error are eliminated.

[0088] In summary, the autonomous calibration method of the terrestrial laser scanner based on intensity features of the present invention can achieve high-precision calibration of TLS calibration parameters by extracting intensity feature points from the point cloud and constructing a set of conditional equations according to their correspondence.

[0089] The effect of the present invention was verified using TLS point cloud data collected at the University of Bonn. For a simplified NIST model containing 10 parameters, the calibration accuracy of each parameter is shown in the following table:

[0090]

[0091] As can be seen from the above table, the present invention can achieve sub-millimeter and arc-second level calibration accuracy for the offset and angle parameters in the calibration model, respectively.

[0092] Those skilled in the art should understand that they can implement variations by combining the prior art with the above embodiments, which will not be described in detail here. Such variations do not affect the essence of the present invention and will not be described in detail here.

[0093] The above describes the preferred embodiments of the present invention. It should be understood that the present invention is not limited to the above-mentioned specific embodiments, and the devices and structures that are not described in detail should be understood to be implemented in a common manner in the art; any technician familiar with the art can use the above-mentioned disclosed methods and technical contents to make many possible changes and modifications to the technical solutions of the present invention without departing from the scope of the technical solutions of the present invention, or modify them into equivalent embodiments of equivalent changes, which does not affect the essential content of the present invention. Therefore, any simple modifications, equivalent changes and modifications made to the above embodiments based on the technical essence of the present invention that do not depart from the content of the technical solutions of the present invention are still within the scope of protection of the technical solutions of the present invention.

Claims

1. A method for autonomous calibration of a terrestrial laser scanner based on intensity characteristics, characterized by: The steps include: Step S1: Obtain point cloud data at two measuring stations respectively, and roughly calculate the relative pose between the measuring stations using the point cloud registration method; Step S2: uniformly select several core points from the point cloud and generate intensity images based on the points near them; Step S3: extracting and matching feature points in the intensity image; Step S4: constructing a system of linear equations to solve the calibration model parameters; The step S3 comprises: Step S31: Utilize The operator locates the salient features in the image, and then extracts and matches the feature points in the image pairs based on the ORB algorithm; Step S32: Detect and remove mismatched feature point pairs using an estimated sample consistency algorithm; Step S33: Calculate the three-dimensional coordinates corresponding to the two-dimensional feature points; the XY coordinates of the feature points can be directly calculated based on their pixel coordinates and pixel size. l Coordinate, Z l The component is obtained by Z l The coordinates are obtained by linear interpolation; then, the above three-dimensional coordinates are converted to the coordinate system of each station and expressed in polar coordinates; the three-dimensional coordinates of the matching feature points are used to estimate the rigid body transformation between the point clouds, and the feature point pairs whose residuals exceed the threshold will be discarded; The step S4 comprises: Step S41: Construct a conditional equation; construct the following equation based on the three-dimensional coordinates of the feature points obtained in step S33: ， Where i = 1…I, n = 1…N, I and N are the number of characteristic points and conditional equations respectively; k and k′ are the numbers of the measuring stations, R l and T l is the rotation and translation between the station coordinate system and the global coordinate system; and is the coordinate of the feature point after calibration; [fx,n fy,n fz,n] T The coordinate difference between a pair of feature points in the global coordinate system is minimized to obtain the calibration parameters and posture parameters of the scanner; The calibration of the observations is performed using the following formula: ， ， Where ξ is the calibration parameter, [r i,k θi,k α i,k ] T is the polar coordinate of the i-th feature point at station k; the calibration model uses the simplified NIST model; Step S42: Solve the conditional equations; [r i,k θ i,k α i,k ] T is the polar coordinate of the i-th feature point at the measuring station k; First, each equation is linearized, and then the Gauss-Helmert model is used to estimate the parameters. Assuming that the coordinate components of the observations are uncorrelated, the covariance matrix is ​​defined as: ， The covariance matrix of all observations is then a block diagonal matrix: Qf=diag([q1…q u …q U ]) U×U # (7) Among them, U is the number of feature points; The initial values ​​are approximate values ​​in the range of 1-10 mm, 10″-50″, and 10″-50″ according to the lidar scanning accuracy and registration accuracy; during the iterative solution process, the variance of the observations is updated through variance component estimation and F test, and the Danish method is used to reweight the residuals; at the same time, observations with residuals exceeding three times the mean error will be eliminated.

2. The method for autonomous calibration of a terrestrial laser scanner based on intensity characteristics according to claim 1, characterized in that: The step S1 comprises: Step S11: Two measuring stations are set up at a certain distance from each other. The collected point cloud data are recorded as point cloud 1 and point cloud 2 respectively. The overlap between the two point clouds must be guaranteed. Step S12: Obtain the initial relative pose between point cloud 1 and point cloud 2 using the iterative closest point method or other point cloud registration algorithms.

3. The method for autonomous calibration of a terrestrial laser scanner based on intensity characteristics according to claim 1, characterized in that: The step S2 comprises: Step S21: Select core points. To reduce the computational burden while ensuring that all different features can be identified, the point cloud 1 is divided into voxels, and the point closest to the center of each voxel is selected as the core point. Step S22: Generate an intensity image; for each core point, use the points within the radius around it to generate a corresponding intensity image.

4. The method for autonomous calibration of a terrestrial laser scanner based on intensity characteristics according to claim 1, characterized in that: The threshold value in step S33 is: the residual error between the same stations does not exceed 10 mm, and the residual error between different stations does not exceed 200 mm.

5. The method for autonomous calibration of a terrestrial laser scanner based on intensity characteristics according to claim 2, characterized in that: The degree of overlap in step S11 is greater than 50%.

6. The method for autonomous calibration of a terrestrial laser scanner based on intensity characteristics according to claim 3, characterized in that: The radius in step S22 is determined by the following empirical formula: ， Among them, ρ1 and ρ2 are the distances from the core point to the two measuring stations, and the value of R is limited to 0.3m to 4.0m. If the number of points in the neighborhood of a core point is less than 100, it will not participate in the subsequent processing. On the contrary, if the number of points in the neighborhood of a core point exceeds 1,000,000, the points around it will be downsampled. p The points in the neighborhood are analyzed using principal components, and a local coordinate system O-XYZl is constructed with the three main directions as coordinate axes and the core point as the origin; all neighboring points will be projected onto the image plane O-XY l The intensity values ​​are interpolated by the griddata function in Matlab to generate an intensity image with a size of m1×m1. The pixel size of the image can be calculated by the following formula: ， Take m1≥1001; in order to obtain the image pairs for subsequent matching, use the initial relative pose to align the core point C p Perform coordinate transformation to obtain C′ p , similarly, C′ p Neighboring points in point cloud 2 will be projected onto the same image plane O-XY l and generate the corresponding intensity image.

Images