Laser radar inertial odometer method based on block-by-block updating

Abstract

The invention discloses a laser radar inertial odometer method based on block-by-block updating, which comprises the following steps: constructing a kinematics equation and an observation equation of a laser radar inertial odometer system based on a block-by-block updating scheme, and performing state propagation on point cloud block laser radar data and inertial measurement data to obtain a prior prediction state; compensating the pose of each point cloud point in the point cloud block, and aligning the point cloud point to the pose of the point cloud block at the last moment; registering the aligned point cloud blocks with a local map, executing nearest neighbor plane search, and then calculating a point-surface residual error; combining the prior state prediction with the posterior state residual error, and performing state iteration updating and iteration convergence to obtain a final result; and then local map updating is carried out. According to the method, the continuous laser radar point cloud data is divided into a plurality of small blocks for batch processing according to the extremely short time interval, so that the time overhead of single update can be remarkably shortened.

Description

Technical Field

[0001] This invention relates to the field of odometer technology, and specifically to a method for lidar inertial odometer based on block-by-block updates. Background Technology

[0002] Existing frame-by-frame LiDAR inertial odometry (IO) schemes have many well-known projects, such as LINS, FAST-LIO, FAST-LIO2, and FASTER-LIO, which have made outstanding contributions to the development of LiDAR IO. LiDAR IO utilizes high-resolution environmental geometry information provided by LiDAR, combined with high-frequency motion prediction from an IMU, to achieve high-precision pose estimation. It boasts strong real-time performance with low latency in frame-by-frame estimation; good environmental adaptability, independent of visual features, and stable operation even in environments with varying lighting and weak textures (such as tunnels and nighttime); and it requires no prior map, enabling mapless localization based on scan-to-scan or scan-to-map matching, making it suitable for exploring unknown environments.

[0003] However, these lidar inertial odometry systems based on frame-by-frame state estimation schemes have many inherent shortcomings:

[0004] (1) It consumes a lot of computing resources. The number of points in the whole frame is large. Even with downsampling, there are still thousands of point cloud objects. Furthermore, subsequent point cloud processing is required, including feature extraction and feature matching, which requires a lot of computing resources and is not suitable for edge devices such as airborne computers.

[0005] (2) Severe motion distortion. The frame-by-frame scanning of the 3D world by lidar itself introduces motion distortion, especially during high-angular-velocity rotation or rapid translation. Existing methods to improve motion distortion include using high-frequency IMU measurements for compensation or assuming uniform velocity. However, these compensation methods are still insufficient for low-frequency radar data or high-speed moving acquisition equipment. Moreover, due to errors in the state-space motion model, motion distortion cannot be effectively compensated, making it difficult to obtain accurate point cloud data to provide reliable support for subsequent inter-frame registration, state propagation, and pose updates. These rapidly accumulating errors lead to inaccurate positioning. Figure 1 Poor consistency or dispersion.

[0006] (3) Low odometry status output and map update frequency. The odometry output frequency estimated frame by frame depends on the point cloud frame release frequency of the acquisition device. The point cloud frame release frequency of general lidar is generally between 10Hz and 100Hz. Therefore, the odometry output frequency estimated frame by frame is also the corresponding frequency. This low output frequency is unfavorable for downstream tasks that require fast response (trajectory planning, motion output and decision analysis, etc.). In fact, downstream tasks may need to wait for pose output for a long time, causing system delay or limited bandwidth for task processing.

[0007] (4) State estimation for high-speed or high-angular-velocity motion is prone to divergence. Frame-by-frame estimation schemes must wait for a frame of point cloud to be packaged before processing the point cloud. However, at this time, due to the large intra-frame motion distortion caused by high-speed motion, the severe impact of the accumulation of distortion error and model error renders various compensation strategies ineffective. Therefore, the odometry will fall into divergence after a short period of time. Therefore, timely and rapid updating of point cloud data is essential. Summary of the Invention

[0008] The purpose of this invention is to provide a method for lidar inertial odometry based on block-by-block updates.

[0009] The technical solution to achieve the objective of this invention is: a method for lidar inertial odometry based on block-by-block updates, comprising the following steps:

[0010] S1. Based on the block-by-block update scheme, construct the kinematic equations and observation equations of the lidar inertial odometry system, perform state propagation on the point cloud block lidar data and inertial measurement data, and obtain the prior predicted state.

[0011] S2. Compensate the pose of each point in the point cloud block and align it to the pose of the point cloud block at the last moment; register the aligned point cloud block with the local map, perform nearest neighbor plane search, and then calculate the point-to-surface residual.

[0012] S3. Combine the prior state prediction and the posterior state residual to perform iterative state updates. The final result can be obtained by iterative convergence. Then, perform local map updates.

[0013] Furthermore, step S1 specifically includes:

[0014] The IMU coordinate system is denoted as I, the global coordinate system is denoted as G, and the first IMU frame is used as the global coordinate system; the radar coordinate system is denoted as L, and the motion model is:

[0015]

[0016] in G p I , G R I It refers to the position and orientation of the IMU in the global coordinate system. G v I This refers to the IMU speed in the global coordinate system. G g is the gravity vector in the global coordinate system, a m ω m It refers to the acceleration and angular velocity measured by the IMU, n a n ωIt is the noise from the IMU's acceleration measurement and angular velocity measurement, b a b ω It is by IMU bias determined by random walk noise, sign This represents the skew-symmetric cross product matrix of vector a;

[0017] Define two operators to parameterize a manifold of dimension n. State error on:

[0018]

[0019] Wherein the formula It is an exponential mapping over SO(3), and Log(·) is its inverse mapping. R and r are the rotation angle and rotation perturbation, respectively, and a and b are n-dimensional vectors. n For Euclidean space;

[0020] For complex manifolds It is the Cartesian product of its submanifold components, and has the following operation:

[0021]

[0022] Let i be the index of the IMU measurement value, and j be the sampling index of the point cloud point, as defined above by the operators. Discretize at a single point cloud sampling period Δt, i.e., Δt = ρ j+1 -ρ j This transforms the continuous model into a discrete model; and each point cloud point sampling timestamp is included within the IMU sampling interval, i.e., ρ j ∈[τ i , τ i+1 ], ρ j τ represents the actual measurement timestamp of a single point cloud point. i This represents the actual measurement timestamp of a single IMU, assuming the IMU measurement value remains constant throughout the entire IMU measurement interval, expressed by the following formula:

[0023]

[0024] Where, x j and x j+1 These represent the timestamps ρ and ρ, respectively. j and ρ j+1 The definitions of the individual states of point clouds, function f, state x, input u, and noise w are as follows:

[0025]

[0026] The index k represents the index of the LiDAR point cloud block, {Lp jLet {j = 1, 2, ..., m} be the points in the k-th point cloud block. L p j Also written as The form indicates that the coordinate system of the three-dimensional point cloud point is the coordinate system of the lidar during the point cloud sampling time ρ. j The local coordinate system, L is omitted when there is no ambiguity. j The subscript j; due to lidar measurement noise, each measurement point L p j Noise consisting of ranging and beam-directing noise L n j Pollution; removing this noise allows us to obtain the true point positions in the local lidar coordinate system.

[0027]

[0028] Use the corresponding lidar pose G T I =( G R I , G p I ) and external parameters I T L After being projected onto the global frame, this ground truth point should be accurately located within the nearest small planar region on the local map, i.e., the observation model should be...

[0029]

[0030] in G u j It is the normal vector of the corresponding plane. G q j It is a point located on that plane; G T I It is contained in the state vector x; therefore, the j-th measurement point L p j The measurement of contribution can be summarized from the above formula into a more compact form, as shown below:

[0031]

[0032] It defines an implicit observation model for the state vector x;

[0033] Once the system receives the point cloud data input from the lidar point cloud block, it begins the forward propagation step, during which process noise w is reduced. j Setting it to 0, expressed by the formula as follows:

[0034]

[0035] Substituting the prior estimated parameters, we obtain the following formula:

[0036]

[0037] Where Δt=ρ j+1 -ρ j That is, the difference in sampling time between two adjacent point cloud points in a LiDAR point cloud block scan;

[0038] To transmit covariance, an error state model is established using the following results:

[0039]

[0040] and F w From arrive Jacobian matrix;

[0041] Let the covariance of the white noise w be denoted as The covariance of the forward propagation can then be calculated iteratively using the following formula:

[0042]

[0043] Furthermore, step 2 specifically involves:

[0044] Forward propagation will continue until the end time t of scanning a point cloud block is reached. k The propagation state and covariance at this point are respectively expressed as: and Or written as and m represents the number of points in a point cloud block, and the sequence number m is the sequence number of the estimated state of the points at the end of the point cloud block; where Represents the true state x k Estimated state with forward propagation The covariance of the errors between them and They are completely equivalent, without distinction, and can be written directly without causing ambiguity.

[0045] When the state of the point cloud block propagates point by point to the end time t of the point cloud block k At this point, a set of prior estimated states is obtained. That is, the set of point cloud blocks { L p1, L p2, ..., L p m} at their respective real sampling time sets {ρ1, ρ2, ..., ρ mThe prior state estimation of the true timestamp, where the time comparison relationship is t. k-1 <ρ1<ρ2<...<ρ m <t k ;

[0046] The state update time is selected as the end time of the point cloud block, i.e., t. k In order to determine the time ρ j to t k The relative motion and motion distortion caused by the time difference are compensated using the following formula:

[0047]

[0048] Equivalent to:

[0049]

[0050] Equivalent to:

[0051]

[0052] Therefore, at time ρ j and point cloud block end time t k This generates a relative pose, i.e., a compensated pose:

[0053]

[0054] After obtaining the motion distortion compensated pose, ρ j Momentary point cloud measurement Projected onto the point cloud block at time t k The point cloud measurement after projection is obtained.

[0055] After the data propagation process, the residual can be constructed; assuming the estimated value of the iterative state x being updated is... The state undergoes its first iteration. Then, each observed LiDAR point L p j Mapped to the global coordinate system, it looks like this:

[0056]

[0057] Next, in the map represented by the voxel structure, find its 5 closest points, and use the found closest neighboring points to fit a normal vector. G u j The center of mass is G q j Local planes, which will be used in the observation model; the observation model in The first-order Tate expansion approximates the following:

[0058]

[0059] in Equivalent to It is the observation model h j (x k L n j )about A Jacobian matrix that takes values ​​at 0; Noise is measured by the original LiDAR. L n j The total measurement noise caused is expected to have a covariance of 0. Will Set to a constant; It is the residual, expressed as follows:

[0060]

[0061] Furthermore, step S3 specifically includes:

[0062] The residuals can be calculated using the constructed residual formula, followed by an iterative update step; the propagation state is then obtained. Covariance For unknown state x k Apply a prior Gaussian distribution; The covariance of the following error states is represented as:

[0063]

[0064] J κ yes about The partial derivative taking the value at 0:

[0065]

[0066] For the first iteration of the Kalman filter, we have J κ =I;

[0067] In addition to the prior distribution, the state distribution can also be obtained from the observation model:

[0068]

[0069] By combining the prior distribution with the measurement model, we obtain the result from... The equivalent representation of state x k The posterior distribution and its maximum a posteriori probability estimate:

[0070]

[0071] Where m represents the number of measurement points, for any invertible matrix M and vector a, we have This maximum a posteriori probability estimation problem can be solved using an iterative Kalman filter;

[0072] make but:

[0073]

[0074] The above process will be repeated until the convergence condition is met; the optimal state and the estimated covariance after convergence are:

[0075]

[0076] Use the updated post-hoc state For each lidar point in the k-th point cloud block L p j Transform to the global coordinate system using the following formula:

[0077]

[0078] Converted lidar points It is inserted into the map represented by the ikdtree structure to update the map.

[0079] A computer device includes a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the steps of the method described above.

[0080] A computer-readable storage medium having a computer program stored thereon that, when executed by a processor, implements the steps of the above-described method.

[0081] A computer program product includes a computer program that, when executed by a processor, implements the steps of the above-described method.

[0082] Compared with existing technologies, the significant advantages of this invention are as follows: This invention proposes a lidar inertial odometry method based on a block-by-block update scheme. By dividing continuous lidar point cloud data into several small blocks for batch processing at extremely short time intervals, the time overhead of a single update can be significantly shortened. Simultaneously, the reduction in point cloud processing volume directly reduces the system's computational load, while the more frequent update rhythm effectively shortens the error accumulation time window, thereby mitigating point cloud distortion caused by the dynamic characteristics of the sensor during motion. Furthermore, the high-frequency odometry status output not only improves the system's response speed to environmental changes but also enhances the ability to track states in high-speed motion or high angular velocity scenarios. Real-time error correction effectively avoids state estimation divergence, providing more reliable technical support for accurate positioning and navigation in complex dynamic environments. Attached Figure Description

[0083] Figure 1 Flowchart for the lidar inertial odometry calculation method for block-by-block updates.

[0084] Figure 2 These are schematic diagrams of the scene near Dingxin Building. The top image is an overall overview of the scene near Dingxin Building; the left side of the middle image shows the mapping effect created by FAST-LIO2; the right side shows the mapping effect created by Block-LIO-input; the left side of the bottom image shows the mapping effect created by FAST-LIO2; the right side shows the mapping effect created by Block-LIO-input.

[0085] Figure 3 The top image shows an overall overview of the fast-shaking scene; the middle image shows the FAST-LIO2 mapping effect on the left and the Block-LIO-input mapping effect on the right; the bottom image shows the FAST-LIO2 mapping effect on the left and the Block-LIO-input mapping effect on the right. Detailed Implementation

[0086] This invention proposes a lidar inertial odometry (IO) based on a block-by-block state estimation scheme to address the problems mentioned in frame-by-frame estimation schemes. First, after receiving a point cloud frame, the block-by-block estimation scheme sorts all points in the frame by time. Then, it slices the point cloud frame into blocks for estimation based on time regions. Within each block, the inertial measurement unit's (IMU) measurements (angular velocity and linear acceleration) are still used as compensation for motion distortion. Therefore, this lidar IO calculation method is referred to as "Block-LIO-input," where "input" indicates that the IMU's measurements are used as system input. This scheme shortens the state update interval, ensuring timely state updates even when accumulated errors are not yet significant, thus avoiding the long-term accumulation of motion distortion and state motion model errors. Timely updates allow for state correction even before the point cloud has accumulated significantly, as the correction update phase requires inverting a portion of the gain matrix, which is very time-consuming if the point cloud is large. It can output odometer status at high frequency to improve system bandwidth and provide a fast and accurate response for downstream tasks. Mitigation of accumulated errors can improve positioning accuracy and system establishment. Figure 1 The accuracy is greatly improved, avoiding divergence during localization and mapping. It can also track high-speed or high-angular-velocity motion, enabling the odometry to handle more complex tasks, such as high-speed autonomous driving, drone racing, and quadruped robot mountain inspection.

[0087] The operation of lidar inertial odometry falls under the category of state estimation theory. A relatively elegant and intuitive way to represent states and disturbances in state estimation theory is using Lie groups and Lie algebras. Lie groups represent the system state, while Lie algebras represent small disturbances, i.e., infinitesimal transformations. This representation also has a good counterpart in filtering theory: the error-state extended Kalman filter. This filter decomposes the true state into a nominal state (known estimated state) and an error state (unknown estimated state, which is a random distribution containing the true state). The correspondence is that Lie groups represent the true state, nominal state, and a priori estimated state, while Lie algebras represent the true error state and the a priori error state.

[0088] In state estimation, the use of Lie groups and Lie algebras offers significant advantages. As manifolds satisfying specific mathematical properties, Lie groups naturally satisfy physical constraints. For example, in attitude estimation, the rotation matrix is ​​a Lie group, representing the set of rotation operations, thus avoiding the occurrence of singularities. Furthermore, Lie groups and Lie algebras have clear geometric meanings, describing the local properties of the system through exponential mappings, making state estimation more physically intuitive. Preserving higher-order terms in Lie algebras reduces linearization errors, thereby improving estimation accuracy. Simultaneously, the combination of Lie groups and Lie algebras leads to faster convergence to the optimal solution, enhances filter stability, and effectively suppresses the effects of noise and outliers. Moreover, this combination reduces computational complexity and improves computational efficiency, making state estimation more efficient.

[0089] In addition to these properties, Lie groups and Lie algebras possess many excellent characteristics. The group operations (multiplication and inverses) of Lie groups are smooth mappings. This means that the group operations are continuous and differentiable on the manifold, facilitating the computation of the Jacobian matrix of the state-space equations. Every Lie group is associated with its corresponding Lie algebra. The Lie algebra is a vector space describing the local properties of the Lie group. The two can be correlated through exponential and logarithmic mappings, and the Lie group satisfies the closure property of multiplication, ensuring arbitrary switching between state and error states. The group operations of Lie groups satisfy the adjoint property, allowing mappings between the covariance matrices of randomly distributed states of different Lie groups; this property is crucial for maintaining consistency.

[0090] In state estimation, the Error-State Extended Kalman Filter (ESEKF) offers significant advantages. By representing the state as an error form, ESEKF avoids the biases introduced when directly linearizing nonlinear systems, thus improving estimation accuracy. Through multiple Taylor expansions to preserve higher-order terms, ESEKF further reduces linearization errors, enhancing estimation accuracy. Furthermore, ESEKF performs local linearization near the reference point, avoiding instabilities caused by global nonlinearity, thereby accelerating convergence and improving filter stability. When facing noise and outliers, ESEKF handles noise effects better and automatically suppresses the influence of outliers through geometric constraints, enhancing robustness. Simultaneously, by updating only local states rather than the entire state space, ESEKF reduces computational complexity and improves computational efficiency. Moreover, the combination of ESEKF and Lie groups and Lie algebras naturally satisfies physical constraints; for example, in attitude estimation, the error state representation directly reflects the true physical motion, making the results easier to interpret.

[0091] To fully leverage the performance advantages of the error-state Kalman filter, iterative optimization, similar to the Gauss-Newton method, is added. The Iterated Error-State Extended Kalman Filter (IESEKF) has significant advantages over the Error-State Extended Kalman Filter (ESEKF). IESEKF updates the state estimate through multiple iterations, continuously correcting the error state, thus significantly improving the estimation accuracy. In contrast, ESEKF performs only one linearization and update, easily leading to the accumulation of linearization errors. IESEKF performs multiple Taylor expansions on the nonlinear function, retaining more higher-order terms, further reducing linearization errors and enhancing estimation accuracy. Furthermore, IESEKF uses the previous correction result as new input in each iteration, progressively optimizing the state estimate, thereby improving convergence speed and filter stability. ESEKF, on the other hand, is prone to divergence in highly nonlinear systems. In terms of robustness, IESEKF handles noise better by adaptively adjusting the state estimate through multiple iterations, effectively suppressing the impact of outliers, while ESEKF is prone to large errors under strong noise. Furthermore, IESEKF updates only the local state in each iteration, significantly reducing computational complexity and improving computational efficiency.

[0092] This invention employs a lidar inertial odometry calculation method that combines Lie group Lie algebra and iterative error state Kalman filtering. This algorithm is computationally efficient, intuitive, and robust.

[0093] The framework diagram of the lidar inertial odometry calculation method based on iterative error-varying Kalman filtering is as follows: Figure 1As shown. The raw point cloud data block of the LiDAR first undergoes state propagation (the prediction step in filtering) in conjunction with the IMU, estimating a reasonable state value for each point in the point cloud at its actual sampling time. After the prediction step of the iterative error state-extended Kalman filter, compensation values ​​are calculated for all points in the point cloud block with the reference point cloud (generally, the point cloud point at the intermediate time is taken as the reference point cloud point). After compensation, all points in the point cloud block are unified to the radar coordinate system of the reference point cloud point. Then, each point in the point cloud block is registered to the local map. This local map is maintained by the ikdtree dynamic data space structure, and the point-to-surface distance residual between all point cloud points and the corresponding plane of the local map is calculated. Then, iterative updates are performed until convergence. The optimized pose is updated, the point cloud block is transformed again, registered, and inserted into the local map. During this process, the ikdtree dynamic data structure will perform addition and deletion operations. In addition, during the local map maintenance process, the local map has a certain maximum boundary value. The local map is maintained in a large cube. Once the current point cloud block position exceeds the boundary of the local map, the local map will delete the previously out-of-bounds point cloud data from the ikdtree dynamic data space structure, that is, the point cloud data farthest from the current point cloud block, to ensure that the local map can always be tracked by the current point cloud block.

[0094] In the LiDAR inertial odometry framework, directly registering the original point cloud to the map without feature extraction and registration modules offers significant advantages. This method simplifies the process, reduces preprocessing steps such as feature extraction and matching, thereby saving computation time and resources and improving the algorithm's real-time performance. Simultaneously, direct registration of the original point cloud enhances robustness, reducing localization failures due to feature loss or interference, especially in complex environments (such as tunnels or dynamic scenes), where the original point cloud contains rich geometric information and better adapts to environmental changes. Furthermore, direct registration of the original point cloud fully utilizes the information within the point cloud, improving positioning accuracy and reducing accumulated errors, particularly in high-speed driving or dynamic environments. This method also has broader applicability, suitable for various environmental conditions, including scenarios with sparse or dynamically changing features.

[0095] The motion model of the laser inertial odometry is described below. The IMU coordinate system is denoted as I, the global coordinate system as G, and the first IMU frame is used as the global coordinate system; the radar coordinate system is denoted as L, and the extrinsic parameters between the radar and IMU are represented as follows:

[0096] I T L =( I R L , I t L The motion model is:

[0097]

[0098] in G p I , G R I It refers to the position and orientation of the IMU in the global coordinate system. G v I This refers to the IMU speed in the global coordinate system. G g is the gravity vector in the global coordinate system, a m ω m It refers to the acceleration and angular velocity measured by the IMU, n a n ω It is the noise from the IMU's acceleration measurement and angular velocity measurement, b a b ω It is by IMU bias determined by random walk noise, sign Let represent the skew-symmetric cross product matrix of vector a.

[0099] In addition, two operators are defined to parameterize a manifold of dimension n. State error on:

[0100]

[0101] Wherein the formula It is an exponential mapping over SO(3), and Log(·) is its inverse mapping.

[0102] For complex manifolds It is the Cartesian product of its submanifold components, and has the following operation:

[0103]

[0104] Let i be the index of the IMU measurement value, and j be the sampling index of the point cloud point, as defined above by the operators. Discretization can be performed at a single point cloud sampling period Δt, i.e., Δt = ρ j+1 -ρ j This transforms the continuous model into a discrete model. Furthermore, the sampling timestamp of each point cloud point is included within the IMU sampling interval, i.e., ρ. j ∈[τ i , τ i+1 ], ρ j τ represents the actual measurement timestamp of a single point cloud point. i This represents the actual measurement timestamp of a single IMU, assuming the IMU measurement value remains constant throughout the entire IMU measurement interval, expressed by the following formula:

[0105]

[0106] Where, x j and x j+1 These represent the timestamps ρ and ρ, respectively. j and ρ j+1 The individual states of the point cloud points, here the function f, state x, input u, and noise w are defined as follows:

[0107]

[0108] An iterative error state-expanded extended Kalman filter is employed. The estimated covariance is characterized in the tangent space (Lye algebra space) of the state estimation. It is assumed that the previous lidar point cloud patch was in state t... k-1 The optimal state estimate at time t is Its covariance matrix is but The covariance of a random error state vector is defined as follows:

[0109]

[0110] in The attitude error is the only one; the rest are standard additive errors (i.e., state variables). The estimated value The error is Intuitively, attitude error This describes the (small) deviation between the true pose and the estimated pose. A major advantage of this error definition is that it can be expressed using a 3×3 covariance matrix. This represents attitude uncertainty. Since attitude has three degrees of freedom (DOF), this is the least parameterized representation. Furthermore, the covariance of the Lie group states... Essentially, it is transmitted through error states, that is... Lie groups themselves lack the traditional concept of "covariance" (because they are not vector spaces); their uncertainty must be represented by the covariance of the tangent space (Lie algebra). The main reason is that Lie groups are manifolds, which are nonlinear, and their elements do not satisfy the linear superposition property of vector spaces, making it impossible to directly define the covariance. By using the error state of the tangent space (Lie algebra), the nonlinear problem can be transformed into a linear space suitable for Gaussian distributions. In error-state Kalman filtering, both prediction and updating occur in the tangent space (error state), and the covariance matrix... Iteratively updated. The covariance matrix of the main state x. Always with Synchronization. In Lie group state estimation, the covariance matrix of the state is the same as that of the error state because the uncertainty of a Lie group can only be described by the tangent space, and the covariance of the error state is the only reasonable mathematical representation. The error state defines the local deviation between the true state and the estimated state, and its distribution directly reflects the confidence level of the estimated state (the covariance matrix of the state describes the uncertainty or confidence level of the estimated state). The design consistency of filtering algorithms requires that the covariances of both be updated synchronously (e.g., ESKF).

[0111] LiDAR point cloud scanning is performed sequentially. Therefore, when the LiDAR undergoes continuous high-speed movement, it samples the point cloud in the 3D world space at different poses, but the state of the point cloud block is only updated at a specific, simultaneous moment. To eliminate this artificially introduced motion distortion or abnormality during point cloud block scanning, a forward sequential propagation method for point cloud points is employed. Based on this, IMU observation data is used to assist in estimating the pose of each point after the LiDAR point cloud block scanning is completed; that is, IMU measurements are used to describe the motion state. Using the estimated positions, these points can be mapped to the end time of a point cloud block scan based on the accurate sampling time of each point in a single scan. Therefore, the points in the point cloud block can be considered as all points sampled simultaneously at the end time of the point cloud block.

[0112] The index k represents the index of the LiDAR point cloud block, { L p j Let {j = 1, 2, ..., m} be the points in the k-th point cloud block. L p j It can also be written as The form indicates that the coordinate system of the three-dimensional point cloud point is the coordinate system of the lidar during the point cloud sampling time ρ. j The local coordinate system, L is omitted when there is no ambiguity. j The subscript j. Due to lidar measurement noise, each measurement point... L p j Noise typically consists of ranging and beamforming noise. L n j Pollution. Removing this noise allows us to obtain the true point positions in the local lidar coordinate system.

[0113]

[0114] Use the corresponding lidar pose G T I =( G R I , G p I ) and external parameters I T L(The fixed extrinsic parameters between the lidar and the inertial measurement unit) After being projected onto the global frame, this real point should be accurately located within the nearest small planar region on the local map, i.e., the observation model is...

[0115]

[0116] in G u j It is the normal vector of the corresponding plane. G q j It is a point located on this plane. It should be noted that... G T I It is contained in the state vector x. Therefore, the j-th measurement point L p j The measurement of contribution can be summarized from the above formula into a more compact form, as shown below:

[0117]

[0118] It defines an implicit observation model for the state vector x.

[0119] Once the system receives the point cloud data input from the lidar point cloud block, it begins the forward propagation step, during which process noise w is reduced. j Setting it to 0, expressed by the formula as follows:

[0120]

[0121] Substituting the prior estimated parameters, we obtain the following formula:

[0122]

[0123] Where Δt=ρ j+1 -ρ j It is the difference in sampling time between two adjacent point cloud points in a LiDAR point cloud block scan.

[0124] To transmit covariance, an error state model can be established using the following results:

[0125] and F w From arrive The Jacobian matrix. Without losing generality, the matrix is... and F w All subscripts are removed from the table to facilitate partial differential calculations. To calculate them, according to the chain rule for partial differential equations, the matrix in equation (1.16) is... and F w The calculation is as follows:

[0126]

[0127] All partial differentials in the equation are calculated as follows:

[0128]

[0129]

[0130] in A(u) -1 The definition is as follows:

[0131]

[0132] Let the covariance of the white noise w be denoted as The covariance of the forward propagation can then be calculated iteratively using the following formula:

[0133]

[0134] Because in the block-by-block update scheme, a point cloud block is aligned to the coordinate system at the end of the last point cloud block, i.e., ρ. m =t k The coordinate system for time is used, so in the algorithm, the end time t from the previous cloud block can be calculated in one step. k-1 up to the end time t of the current point cloud block k The algorithm first propagates the prior state and then compensates for other points in the current point cloud block in turn. However, compared to other algorithms, propagating the state point by point before motion distortion compensation makes the algorithm flow more natural and concise. Compared to FAST-LIO and FAST-LIO2 algorithms, it avoids the backpropagation process and is more concise and easier to understand. Therefore, this invention uses the method of propagating the state point by point.

[0135] Forward propagation will continue until the end time t of scanning a point cloud block is reached. k The propagation state and covariance at this point are respectively expressed as: and It can also be written as and m represents the number of points in a point cloud block, and the sequence number m is the sequence number of the estimated state of the points at the end of the point cloud block. Represents the real state x k Estimated state with forward propagation The covariance of the errors between them, here and As described above, they are completely equivalent, so no distinction is made, and they can be written directly without causing ambiguity. When the state of the point cloud block propagates point by point to the end time t of the point cloud block k At this point, a set of prior estimated states is obtained. That is, the set of point cloud blocks { L p1, L p2, ..., L p m} at their respective real sampling time sets {ρ1, ρ2, ..., ρ m The prior state estimation of the true timestamps, and this neat correspondence, provide the basis for subsequent motion distortion compensation, where the time comparison relationship is t. k-1 <ρ1<ρ2<...<ρ m <t k Because of this temporal inconsistency, and because state updates must occur at a specific moment, this inconsistency leads to inconsistencies in the lidar coordinate system, resulting in motion distortion. Compensating for this artificially introduced error is necessary. Generally, the state update time is chosen as the end time of the point cloud block, i.e., t. k In order to determine the time ρ j to t k The relative motion and motion distortion caused by the time difference can be compensated using the following formula:

[0136]

[0137] Equivalent to:

[0138]

[0139] Equivalent to:

[0140]

[0141] Therefore, at time ρ j and point cloud block end time t k This generates a relative pose, i.e., a compensated pose:

[0142]

[0143] After obtaining this motion distortion compensation pose, as described in the formula, ρ can be... j Momentary point cloud measurement Projected onto the point cloud block at time t k The point cloud measurement after projection is obtained. The next step is to construct the residual.

[0144] After the data propagation process described above, the residual can be constructed. Assume the estimated value of the iterative state x being updated is... The predicted state originates from the propagation in formula (1.15), and the state undergoes its first iteration. Then, each observed LiDAR point L pj Mapped to the global coordinate system, as follows:

[0145]

[0146] Next, the five closest points are found in the map represented by the voxel structure. Then, a normal vector is fitted using the found closest neighboring points. G u j The center of mass is G q j These are local planes that will be used in the observation model. The observation model will be placed in... The first-order Tate expansion approximates the following:

[0147]

[0148] in (equivalent to) ), It is the observation model h j (x k L n j )about A Jacobian matrix that takes values ​​at 0. Noise is measured by the original LiDAR. L n j The total measurement noise caused is expected to have a covariance of 0. The algorithm of this invention will Set to a constant. It is the residual, expressed as follows:

[0149]

[0150] The residuals can be calculated using the constructed residual formula, followed by an iterative update step. The propagation state obtained earlier... Covariance For unknown state x k A prior Gaussian distribution was applied. Specifically, The covariance of the following error states is represented by:

[0151]

[0152] J κ yes about Partial derivative taking the value at 0:

[0153]

[0154] For the first iteration of the Kalman filter, we have J κ=I.

[0155] In addition to the prior distribution, the state distribution can also be obtained from the observation model (1.32):

[0156]

[0157] Combining the prior distribution in equation (1.34) with the measurement model in equation (1.36), we can obtain the result from... The equivalent representation of state x k The posterior distribution and its maximum a posteriori probability estimate:

[0158]

[0159] Where m represents the number of measurement points, for any invertible matrix M and a vector a of appropriate dimension, we have: The maximum a posteriori probability estimation problem can be solved by using an iterative Kalman filter.

[0160] For ease of calculation, let but:

[0161]

[0162] It is worth noting that the Kalman gain K is not calculated using an observation dimension matrix, but rather by finding the inverse of a state dimension matrix.

[0163] The above process will be repeated until the harvest conditions are met, for example... The optimal state and the estimated covariance after convergence are:

[0164]

[0165] Use the updated post-hoc state For each lidar point in the k-th point cloud block L p j Transform to the global coordinate system using the following formula:

[0166]

[0167] Converted lidar points It is inserted into the map represented by the ikdtree structure to update the map.

[0168] The following describes in detail the point cloud compression method based on deep learning and utilizing the implicit representation of the distance field according to the present invention, with reference to the embodiments and accompanying drawings.

[0169] Example

[0170] The proposed block-by-block update algorithm is implemented in C++ and the Robot Operating System (ROS). The iterative error extended Kalman filter is implemented based on the IKFoM toolset. An incremental kd-tree, or ikd-Tree, developed in FAST-LIO2, is used as the dynamic spatial data structure with its default parameters: local map size 2000m, spatial downsampling resolution 0.25m, ikd-Tree rebalancing thresholds of 0.6 and 0.5, and a subtree size threshold of 1500 for parallel reconstruction (in the second thread).

[0171] Although the system of this invention is designed to perform state estimation after each received LiDAR point cloud block, in practice, due to limitations imposed by the available drivers provided by the LiDAR manufacturer, LiDAR points are packaged after accumulating complete scan frames and then sent to the LIO system. To accommodate this practical limitation, this invention sorts the LiDAR points and IMU data contained in the received data packet according to their respective timestamps. The sorted data is then processed block by block.

[0172] In all evaluations, the method of this invention was compared with the state-of-the-art frame-based odometry FAST-LIO2. All experimental results for FAST-LIO2 were collected using the public version of FAST-LIO2 and its default parameter values. Since FAST-LIO2 spatially downsamples each received LiDAR scan at a resolution of 0.3, this invention also performs such spatial downsampling before data sorting and block-by-block updates for a fair comparison.

[0173] This invention develops a LiDAR inertial sensor suite for collecting point cloud data of 3D scenes in real-world environments. It consists of a solid-state 3D LiDAR, namely LivoxAvia. The LivoxAvia LiDAR features a conical field of view of 70.4 degrees horizontally and 77.2 degrees vertically, and an unconventional non-repeating scan mode, generating point cloud measurements at 240 kHz, as well as a built-in IMU (model BMI088) capable of generating 200 Hz IMU data. The point cloud data and IMU data are packaged at frequencies adjustable between 10 Hz and 100 Hz. The onboard computer uses an NVIDIA Jetson Xavier NX edge computing module, featuring a 6-core NVIDIA Carmel ARM@v8.2 64-bit processor with a single-core maximum speed of 1.9 GHz, 6MB L2 + 4MB L3 cache, and other hardware performance. It comes pre-installed with a customized version of Ubuntu 18.04 LTS (Bionic Beaver) operating system, with a default kernel of 4.9.253-tegra (version 32.7.3), making it very user-friendly for Linux users. Furthermore, this invention manually installed the Robot Operating System (ROS) as the core framework for robot development, using the ROS SmartCelodicMorenia version, which is the officially recommended version of Ubuntu 18.04 (BionicBeaver). Finally, a 5-inch display screen with a resolution of 800×480 was incorporated for real-time display of localization and mapping effects.

[0174] The current mainstream LiDAR inertial odometry calculation method based on Iterative Error State Extended Kalman Filter (IESEKF) uses a frame-by-frame estimation approach, causing the odometry output frequency to depend on the device's output frame rate. The full-state output frequency of mainstream algorithms typically ranges from tens of hertz to hundreds of hertz. Compared to the real-time full-state odometry output frequency of hundreds of hertz achieved by the method of this invention, the fast response capability of these algorithms is severely inadequate. Furthermore, if a frame has a long time span (e.g., 100ms), it can lead to severe motion distortion, significantly affecting positioning accuracy and system establishment. Figure 1 To address these issues, this invention proposes an algorithm based on block-by-block state estimation, namely Block-LIO-input. This scheme shortens the state estimation time, with the point cloud block time span on the order of milliseconds, compared to the tens to hundreds of milliseconds typically found in mainstream algorithms. This section will verify the effectiveness of this scheme.

[0175] In the lidar inertial odometry framework proposed in this invention, the algorithm updates the state odometry in real time upon receiving each lidar point cloud block, thereby achieving block-by-block state estimation. It can be seen that the block-by-block update scheme used in this invention differs significantly from the odometry and mapping algorithms of the existing mainstream frame-by-frame update scheme. Block-by-block updates shorten the time interval between state updates, thus enabling extremely high-frequency odometry state output and mapping updates. Ideally, the odometry output frequency of the block-by-block update scheme can reach thousands of hertz. This high-rate odometry state update prevents the point cloud from accumulating too large a quantity, which would overload the computational load during the correction and update phase. Furthermore, it allows for timely correction and updating of the estimated state before the odometry model error and estimation error become excessive during the prediction propagation step, thereby increasing the odometry bandwidth and ensuring that the algorithm effectively tracks state changes during high-speed or high-angular-velocity motion.

[0176] Another advantage of the block-by-block update scheme proposed in this invention is that it greatly alleviates the problem of intra-frame motion distortion or motion loss. This characteristic can improve the accuracy and consistency of odometry in subsequent positioning and mapping processes. Currently, the mainstream framework for frame-by-frame state estimation accumulates all point clouds within a single frame's sampling interval, assuming that all point clouds within the frame were sampled at the same time for odometry state update correction. However, in reality, point clouds are sampled sequentially, each with its own actual sampling time. This blind accumulation inevitably leads to intra-frame point cloud motion distortion or loss. Moreover, to counteract this motion distortion, some compensation measures must be adopted, such as inertial measurement unit integration, uniform velocity assumption, or B-spline analysis of continuous-time trajectories.

[0177] In contrast, the block-by-block state update scheme of the present invention shortens the point cloud accumulation time, and with the addition of IMU compensation, it greatly improves this motion distortion or aberration.

[0178] To verify the effectiveness of the proposed block-by-block state update scheme in resolving motion distortion, we used... Figure 1 The handheld mapping device shown collected several data packet sequences. Two different scenarios were tested: an outdoor semi-structured scenario near Dingxin Building at Nanjing University of Science and Technology, named the "Dingxin Building Nearby Scene," and an indoor structured scenario involving rapid handheld device rotation and shaking, named the "Rapid Shaking Scene." In all sequences, the actual motion trajectory returned to the starting point, enabling drift calculations. The LiDAR data packet frequency was 10 Hz.

[0179] The mapping results for the scene near Dingxin Building and the fast-shaking scene are as follows: Figure 2 and Figure 3As shown in the image above. In each set of images, a global view of the final mapping result is first presented (i.e., the top image of the set), and then the focus is on a specific local area containing large planes (such as walls) (the middle and bottom images of the set). The consistency of points on the walls is studied (the middle and bottom images of the set), from which the mapping accuracy and intra-frame distortion of FAST-LIO2 (i.e., the left image of the middle image, the left image of the bottom image) and Block-LIO-input (i.e., the right image of the middle image, the right image of the bottom image) can be compared.

[0180] A comparison of the middle and lower images in Figure 2 shows that the Block-LIO-input module... Figure 1 The consistency is better than FAST-LIO2. In the lower middle image, the FAST-LIO2 (left) portion, compared to the Block-LIO-input (right) portion, shows blurred and divergent walls, distortion, poor consistency, and numerous outliers near windows. However, the Block-LIO-input mapping portion (middle and lower right images) maintains high consistency, with neat wall outlines without outliers, sharp and clear edges, flat walls, and a gentle and reasonable point cloud distribution. From Figure 3 A side-by-side comparison of the middle and lower images reveals that the FAST-LIO2 mapping (left side of the middle and lower images) exhibits severe ghosting and misalignment, resulting in mapping failure. In contrast, the Block-LIO-input mapping remains effective and exhibits high consistency, with clear wall outlines, no ghosting or misalignment, and evenly distributed point clouds converging on the wall surface.

[0181] FAST-LIO2 uses backpropagation based on IMU measurements to project all points from a single scan onto the pose at the end of the scan. Since a single frame of point cloud spans 100 milliseconds, the time interval is too long. Even with appropriate motion distortion compensation, it's difficult to eliminate the severe impact of this long-term accumulated error. Therefore, intra-frame motion distortion compensation for such long time intervals is futile under high-speed motion. Furthermore, long-time motion distortion compensation is easily affected by IMU measurement noise, changes in IMU angular velocity and acceleration offsets, and the limited IMU sampling rate. Significant changes in motion state and offset during IMU sampling can lead to large IMU propagation errors. Assuming that angular velocity, acceleration, and their respective offsets remain constant within an IMU sampling interval, this motion distortion compensation scheme using IMU measurements introduces even more uncertainty. Therefore, as the experiments show, while the mapping process can be maintained at low speeds, localization and mapping fail under high-speed motion conditions.

[0182] The block-by-block update scheme proposed in this invention not only shortens the point cloud processing interval and allows for timely state updates, but also avoids the accumulation of large system and model errors during forward propagation. For motion distortion compensation, a shorter time frame yields better results, as it minimizes the accumulation of IMU noise and model errors, thus providing reliable data support for subsequent state updates. For mapping updates, a shorter time frame minimizes the impact of various artificially introduced noises. Figure 1 The higher the consistency, the better. Coupling between different system modules (propagation, update, and mapping, etc.) also affects the overall effect. For example, the smaller the cumulative motion distortion error, the more accurate other modules become, and these system modules will then affect the subsequent motion distortion compensation effect. In summary, the experimental mapping results show that the block-by-block update scheme is significantly better than the frame-by-frame update scheme in terms of motion distortion compensation capability, reducing the accumulation and impact of cumulative errors.

[0183] The drift characteristics of FAST-LIO2 and Block-LIO-input are summarized in Table 1. To ensure accurate comparison of drift magnitude, the starting positions (including position and orientation) of all point cloud data sequences were marked to ensure that the displacement distance between the starting and ending positions is exactly zero, and the starting and ending orientations coincide. It can be seen that for scene sequences near Dingxin Building, the drift characteristics of FAST-LIO2 and Block-LIO-input are basically equivalent. For rapidly jittering scene sequences, because FAST-LIO2 failed to construct the map, the drift is much greater than that of Block-LIO-input.

[0184] Table 1. Odometer Drift Comparison (m)

[0185] FAST-LIO2 Block-LIO-input Scenes near Dingxin Building: 0.098 0.056 Rapid shaking scenario: 0.431 0.078

[0186] The above embodiments verify that the proposed solution of the present invention has advantages in odometer output frequency, construction, and other aspects. Figure 1 Improvement in sexual function.

Claims

1. A method for lidar inertial odometry based on block-by-block updates, characterized in that, Includes the following steps: S1. Based on the block-by-block update scheme, construct the kinematic equations and observation equations of the lidar inertial odometry system, perform state propagation on the point cloud block lidar data and inertial measurement data, and obtain the prior predicted state. S2. Compensate the pose of each point in the point cloud block and align it to the pose of the point cloud block at the last moment; register the aligned point cloud block with the local map, perform nearest neighbor plane search, and then calculate the point-to-surface residual. S3. Combine the prior state prediction and the posterior state residual to perform iterative state updates. The final result can be obtained by iterative convergence. Then, perform local map updates.

2. The lidar inertial odometry method based on block-by-block update according to claim 1, characterized in that, Step S1 is as follows: The IMU coordinate system is denoted as I, the global coordinate system is denoted as G, and the first IMU frame is used as the global coordinate system; the radar coordinate system is denoted as L, and the motion model is: Among them G pI G RI The position and orientation of the IMU in the global coordinate system, G vI The IMU velocity in the global coordinate system, G g The gravity vector in the global coordinate system, a m ω m It refers to the acceleration and angular velocity measured by the IMU, n a n ω It is the noise from the IMU's acceleration measurement and angular velocity measurement, b a b ω It is by IMU bias determined by random walk noise, sign This represents the skew-symmetric cross product matrix of vector a; Define two operators to parameterize a manifold of dimension n. State error on: Wherein the formula It is an exponential mapping over SO(3), and Log(·) is its inverse mapping. R and r are the rotation angle and rotation perturbation, respectively, and a and b are n-dimensional vectors. n For Euclidean space; For complex manifolds It is the Cartesian product of its submanifold components, and has the following operation: Let i be the index of the IMU measurement value, and j be the sampling index of the point cloud point, as defined above by the operators. Discretize at a single point cloud sampling period Δt, i.e., Δt = ρ j+1 -ρ j The continuous model is converted into a discrete model; and the sampling timestamp of each point cloud point is included in the IMU sampling interval, i.e., ρ j ∈[τ i , τ i+1 , ρ j τ represents the actual measurement timestamp of a single point cloud point. i This represents the actual measurement timestamp of a single IMU, assuming the IMU measurement value remains constant throughout the entire IMU measurement interval, expressed by the following formula: Where, x j and x j+1 These represent the timestamps ρ and ρ, respectively. j and ρ j+1 The definitions of the individual states of point clouds, function f, state x, input u, and noise w are as follows: The index k represents the index of the LiDAR point cloud block, { L p j Let {j = 1, 2, ..., m} be the points in the k-th point cloud block. L p j Also written as The form indicates that the coordinate system of the three-dimensional point cloud point is the coordinate system of the lidar during the point cloud sampling time ρ. j The local coordinate system, L is omitted when there is no ambiguity. j The subscript j; due to lidar measurement noise, each measurement point L p j Noise consisting of ranging and beam-directing noise L n j Pollution; removing this noise allows us to obtain the true point positions in the local lidar coordinate system. Use the corresponding lidar pose G T I =( G R I , G p I ) and external parameters I T L After being projected onto the global frame, this ground truth point should be accurately located within the nearest small planar region on the local map, i.e., the observation model should be... in G u j It is the normal vector of the corresponding plane. G q j It is a point located on that plane; G T I It is contained in the state vector x; therefore, the j-th measurement point L p j The measurement of contribution can be summarized from the above formula into a more compact form, as shown below: It defines an implicit observation model for the state vector x; Once the system receives the point cloud data input from the lidar point cloud block, it begins the forward propagation step, during which process noise w is reduced. j Setting it to 0, expressed by the formula as follows: Substituting the prior estimated parameters, we obtain the following formula: Where Δt=ρ j+1 -ρ j That is, the difference in sampling time between two adjacent point cloud points in a LiDAR point cloud block scan; To transmit covariance, an error state model is established using the following results: and F w From arrive The Jacobian matrix; Let the covariance of the white noise w be denoted as The covariance of the forward propagation is then calculated iteratively using the following formula:

3. The lidar inertial odometry method based on block-by-block update according to claim 2, characterized in that, Step 2 is as follows: Forward propagation will continue until the end time t of scanning a point cloud block is reached. k The propagation state and covariance at this point are respectively expressed as: and Or written as and m represents the number of points in a point cloud block, and the sequence number m is the sequence number of the estimated state of the points at the end of the point cloud block; where Represents the true state x k Estimated state with forward propagation The covariance of the errors between them and They are completely equivalent, without distinction, and can be written directly without causing ambiguity. When the state of the point cloud block propagates point by point to the end time t of the point cloud block k At this point, a set of prior estimated states is obtained. That is, the set of point cloud blocks { L p1, L p2, ..., L p m } at their respective real sampling time sets {ρ1, ρ2, ..., ρ m The prior state estimation of the true timestamp, where the time comparison relationship is t. k-1 <ρ1<ρ2<...<ρ m <t k ; The state update time is selected as the end time of the point cloud block, i.e., t. k In order to determine the time ρ j to t k The relative motion and motion distortion caused by the time difference are compensated using the following formula: Equivalent to: Equivalent to: At time ρ j and point cloud block end time t k A relative pose is generated between them, namely the compensated pose: After obtaining the motion distortion compensated pose, ρ j Momentary point cloud measurement Projected onto the point cloud block at time t k The point cloud measurement after projection is obtained. After the data propagation process, the residual can be constructed; assuming the estimated value of the iterative state x being updated is... The state undergoes its first iteration. Then, each observed LiDAR point L p j Mapped to the global coordinate system, it looks like this: Next, in the map represented by the voxel structure, find its 5 closest points, and use the found closest neighboring points to fit a normal vector. G u j The center of mass is G q j The local planes, which will be used in the observation model; the observation model in The first-order Tate expansion approximates the following: in Equivalent to It is the observation model h j (x k L n j )about A Jacobian matrix that takes values ​​at 0; Noise is measured by the original LiDAR. L n j The total measurement noise caused is expected to have a covariance of 0. Will Set to a constant; It is the residual, expressed as follows:

4. The lidar inertial odometry method based on block-by-block update according to claim 3, characterized in that, Step S3 is as follows: The residuals are calculated using the pre-constructed residual formula, followed by an iterative update step; the propagation state is then obtained. Covariance For unknown state x k Apply a prior Gaussian distribution; The covariance of the following error states is represented as: J κ yes about The partial derivative taking the value at 0: For the first iteration of the Kalman filter, we have J κ =I; In addition to the prior distribution, the state distribution can also be obtained from the observation model: By combining the prior distribution with the measurement model, we obtain the result from... The equivalent representation of state x k The posterior distribution and its maximum a posteriori probability estimate: Where m represents the number of measurement points, for any invertible matrix M and vector a, we have This maximum a posteriori probability estimation problem can be solved using an iterative Kalman filter; make but: The above process will be repeated until the convergence condition is met; the optimal state and the estimated covariance after convergence are: Use the updated post-hoc state For each lidar point in the k-th point cloud block L p j Transform to the global coordinate system using the following formula: Converted lidar points It is inserted into the map represented by the ikdtree structure to update the map.

5. A computer device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the steps of the method described in any one of claims 1-4.

6. A computer-readable storage medium having a computer program stored thereon, characterized in that, When executed by a processor, the program implements the steps of the method described in any of claims 1-4.

7. A computer program product, comprising a computer program, characterized in that, When executed by a processor, the computer program implements the steps of the method described in any one of claims 1-4.

Images