Calibrating Range Measurements of Lidars Using Fixed Landmarks in Unknown Positions

We consider the problem of calibrating range measurements of a Light Detection and Ranging (lidar) sensor that is dealing with the sensor nonlinearity and heteroskedastic, range-dependent, measurement error. We solved the calibration problem without using additional hardware, but rather exploiting assumptions on the environment surrounding the sensor during the calibration procedure. More specifically we consider the assumption of calibrating the sensor by placing it in an environment so that its measurements lie in a 2D plane that is parallel to the ground. Then, its measurements come from fixed objects that develop orthogonally w.r.t. the ground, so that they may be considered as fixed points in an inertial reference frame. Moreover, we consider the intuition that moving the distance sensor within this environment implies that its measurements should be such that the relative distances and angles among the fixed points above remain the same. We thus exploit this intuition to cast the sensor calibration problem as making its measurements comply with this assumption that “fixed features shall have fixed relative distances and angles”. The resulting calibration procedure does thus not need to use additional (typically expensive) equipment, nor deploy special hardware. As for the proposed estimation strategies, from a mathematical perspective we consider models that lead to analytically solvable equations, so to enable deployment in embedded systems. Besides proposing the estimators we moreover analyze their statistical performance both in simulation and with field tests. We report the dependency of the MSE performance of the calibration procedure as a function of the sensor noise levels, and observe that in field tests the approach can lead to a tenfold improvement in the accuracy of the raw measurements.


Introduction
Localization is essential for applications where robots shall move precisely in the surroundings, and is typically performed by leveraging on measurements acquired through distance sensors. Calibrating these distance measurement sensors so that their readings are as accurate as possible is thus a preliminary task that is essential for achieving good navigation performance at a later stage.
To make a practical and typical (but not exhaustive) example, calibrating a distance sensor may mean considering a measurement model of the type r = f bias (d) + f st.dev (d)e (1) where r is the noisy sensor reading, d is the true distance, f bias (·) includes the true distance plus an unknown bias term, and f st.dev (·) is a factor modulating the measurement noise whose stochasticity is induced by a random variable e typically assumed standard Gaussian and independent and identically distributed (iid). Such a model is motivated by practical real-life situations, such as the one depicted in Figure 1, and the uncertainty constructed according to the work in [1]. A strategy for calibrating a model such as (1) would in this case correspond to estimating f bias and f st.dev so that these factors can be accounted for when postprocessing new measurements from the sensor.
To make a practical and typical (but not exhaustive) example, calibrating a distance sensor may mean considering a measurement model of the type where r is the noisy sensor reading, d is the true distance, f bias (·) includes the true distance plus an 27 unknown bias term, and f st.dev (·) is a factor modulating the measurement noise whose stochasticity is 28 induced by a random variable e typically assumed standard Gaussian and independent and identically 29 distributed (iid). Such a model is motivated by practical real-life situations, such as the one depicted 30 in Figure 1, and the uncertainty constructed according to [1]. A strategy for calibrating a model such 31 as (1) would in this case correspond to estimating f bias and f st.dev so that these factors can be accounted 32 for when post-processing new measurements from the sensor. The operation of calibrating a sensor is typically performed by comparing the raw measurements 34 of the sensor against readings from an external and sufficiently more accurate system (e.g., a motion 35 capture system) that is considered as ground truth (e.g., as in [2]). Acquiring this ground truth (and thus 36 this external and sufficiently more accurate system), however, may be expensive and time-consuming. 37 It may thus be beneficial to find strategies that substitute this information-acquisition step with more 38 easily implementable and cheaper approaches. 39 For example, this substitution can be performed as follows: assume that in structured 40 environments certain structures do not move (e.g., walls, doors and corners in the built environment; 41 trunks of trees in a forest, etc.). Assume moreover that these structures may produce specific and easily 42 recognizable signatures in the readings (e.g., trunks of trees in a forest produce somehow circular 43 shaped features, as in Figure 2). 44 Figure 1. An example of a series of raw measurements obtained using a noncalibrated distance sensor (in this case the triangulation lidar described in Section 4.3).
The operation of calibrating a sensor is typically performed by comparing the raw measurements of the sensor against readings from an external and sufficiently more accurate system (e.g., a motion capture system) that is considered as ground truth (e.g., as in [2]). Acquiring this ground truth (and thus this external and sufficiently more accurate system), however, may be expensive and time-consuming. It may thus be beneficial to find strategies that substitute this information acquisition step with more easily implementable and cheaper approaches.
For example, this substitution can be performed as follows: assume that in structured environments certain structures do not move (e.g., walls, doors, and corners in the built environment trunks of trees in a forest, etc.). Assume moreover that these structures may produce specific and easily recognizable signatures in the readings (e.g., trunks of trees in a forest produce somehow circular shaped features, as in Figure 2). Version December 18, 2020 submitted to Sensors 3 of 23 1m 1m As soon as the structures do not move, these signatures may be considered as fixed points 45 in an inertial reference frame. If a distance sensor is moved within this environment, then the 46 measurements from the sensor referring to these fixed points should be such that the relative distances 47 and angles among these fixed points remain the same. The calibration process may then be cast, from 48 an intuitive perspective, as finding models like (1) for which the measurement process complies with 49 the assumption that "fixed features shall have fixed relative distances and angles". 50 The goal of this paper is thus to understand how to leverage these assumptions on the structure of the 51 environment surrounding the distance sensors for the purpose of building statistically accurate distance sensors 52 calibration strategies. As soon as the structures do not move, these signatures may be considered as fixed points in an inertial reference frame. If a distance sensor is moved within this environment, then the measurements from the sensor referring to these fixed points should be such that the relative distances and angles among these fixed points remain the same. The calibration process may then be cast, from an intuitive perspective, as finding models like (1) for which the measurement process complies with the assumption that "fixed features shall have fixed relative distances and angles".
The goal of this paper is thus to understand how to leverage these assumptions on the structure of the environment surrounding the distance sensors for the purpose of building statistically accurate distance sensors calibration strategies.
To do so we will thus consider using a simple strategy: (a) place some artificial landmarks (i.e., some poles) in random positions in space; (b) calibrate the sensor by making its measurements comply with the fixed-world assumption above.

Literature Review
The strategy described in the previous subsection relates to the existing literature as in the following. First of all, distance sensors are often used for reconstructing environmental maps used by robots to move without colliding with obstacles, as, for example, in [3] (We also note that map generation is not the only application; for example, forestry applications use lidars to measure and monitor the growth of forests, compute trunk's diameters, and calculate the density of trunks or canopies [4,5]). Several strategies have been proposed to improve distance sensor performance and accuracy through statistical manipulation of their measurement processes. For example, the statistical sensor model for ultrasound sensors was presented in [6] with calibration algorithms in [7] and a good review on odometer calibration is presented in [8]. As a generic definition, statistical sensor calibration is the process of improving sensor accuracy and/or precision through transforming the measurements into something closer to ground truth via combining information about the same sensed quantity that is obtained using different sources of information like ground truth data. Sensor calibration for installation error has been well studied in the literature and, in general, is solved using nonlinear optimization, see for example in [9]. However, in this research, we are not considering the extrinsic calibration of the sensor but we focus on the intrinsic calibration that is dealing with the sensor nonlinearity and heteroskedastic, range-dependent, measurement error. Unfortunately, ground truth is not always available and in some cases, if exists, it is very expensive. Therefore, it is usually substituted with one of the following strategies.
• Certain assumptions about the sensor movement and about the surrounding environment, in which the calibration process is shaped as joint parameters and state estimation, for example, lidar calibration from linear motion [10]. • Another strategy for substituting ground truth information with some other information is to implement appropriate sensor fusion strategies, i.e., to combine redundant information from independent distance sensors. Such a strategy has been used in [10,11], where approximated Expectation Maximization (EM) procedures (in the former) and Markov chain Monte Carlo (MCMC) techniques under Bayesian frameworks (in the later) are used for joint parameter and state estimation combining information from lidars, odometry, and ultrasound sensors. Calibrating the intrinsic parameters of one beam based on other beams of rotating multi-beam lidar attracted large amount of research, for example, in [12][13][14][15]. We note that sensor fusion is a vast topic and there are many publications on calibrating other sensors, for example, magnetometer calibration using inertial sensors in [16], camera and IMU calibration in [17], and lidar and camera calibration in [18]. However, here, we are interested only in calibration that is related to lidars. • The last strategy is to use assumptions on the environment, for example, odometer calibration with localization [19]. Another example is to use the planar feature in the environment to calibrate lidars. Originally plane-based calibration was presented for calibrating airborne lidars in [20,21]. Then, the authors of [22]  both the geometric and temporal parameters based on Rényi Quadratic Entropy to formulate an optimization problem that maximizes the quality of the point cloud.
As said above, here we specifically investigate how to substitute ground truth information with assumptions on the environment. Our strategy will intrinsically require localizing (in a sense to be specified later) the position of the sensor within the surrounding environment. This means that our paper relates to the existing literature on localizing sensors in space using noisy measurements of distances from landmarks or beacons. To the best of our knowledge, it is possible to do so using three different approaches: • triangulation, where the position is determined through measuring the angles between the sensing device and the known landmarks (see, e.g., in [24]); • trilateration, where the position is determined through measuring the distances between the device and the landmarks (see, e.g., in [25,26]); • triangulateration, a strategy that combines both of the above (see, e.g., in [27,28]).
Generally, most algorithms use either triangulation or trilateration alone, as they require less information from the sensor (measuring both distances and angles, indeed, normally requires more expensive hardware). For this reason several studies on how to localize the position of a sensor using landmarks or beacons mostly use triangulation or trilateration approaches. For example, the authors of [29][30][31] all propose different techniques for selflocalization using landmarks or beacons and triangulation concepts, while the authors of [25,26,32,33] all use trilateration.
In the literature mentioned above the solutions are based on the assumption that he landmark positions are known which is not the case in our setup. Extensive research has been done to solve mobile device localization given several known mobile base stations [34]. However, in this paper, we relax this assumption into a more general case where the landmark positions are assumed to be completely unknown. Instead we assume to know imprecise information about the sensor's new position with respect to the previous one. This kind of information is actually the control commands to the robot which is always available for robot moving in its environment. Furthermore, to make the calibration process independent on the robot dynamical model, we assume to take measurements only when the robot is not moving (stand still).
Scan-matching techniques like ICP [35] and NDT [36,37] are also commonly used to localize a sensor in space. This class of methods determines the relative transformation between two lidar scans by minimizing surface-to-surface distances using all points in the scans, as opposed to a sparse set of extracted landmarks as is the case with triangulation and trilateration. Scan matching could potentially be used instead of, or in addition to, the control commands used as input to estimate the sensor pose in this work. However, this pose estimation is not the focus of the present paper, but rather the calibration of the range-dependent sensor noise.
Finally, we note that our strategy is specifically designed for cheap distance sensors: generally, the more accurate and precise a sensor is, the more useful (and, at the same time, likely expensive) it is. Our focus is on enabling software-based improvements of cheap sensors so that, by adding a bit of statistical processing, we extend their applicability. For this reason, we use triangulation lidar sensors as a practical and motivating example. Therefore, our paper relates also with the literature around the calibration of these sensors, and thus to the analyses on the effect of color of the target on the measured distance [38]; the works in [39,40], that build two partially different statistical models (the former homoskedastic, the latter heteroskedastic) and thus two slightly different calibration procedures based on ground truth information, using Weighted Least Squares (WLS) for parameter estimation and Akaike Information Criterion (AIC) for model selection; and the work in [41], that extends upon the work in [40] by including the effects of beam angles in the calibration process.

Statement of Contributions
Summarizing, we propose and validate a strategy that uses triangulateration concepts for calibrating distance sensors that return 2D measurements (i.e., both angles and distances). The algorithm leverages on placing the distance sensor inaccurately in equally distant positions along straight paths and making the measurements from the sensor comply with the assumption that the landmarks do not move, plus some other practical assumptions listed exhaustively in Section 2. The strategy is intended to be applicable at least (but not only) for the very specific situation where vacuum cleaning robots move within an apartment, and are able-by moving around and detecting obstacles-to self-calibrate their distance sensors.
More specifically, our strategy works as follows. First, we assume the pre-existence of a procedure that correctly identifies and distinguishes landmarks in the 2D measurements stream. Then, laddering on this knowledge, we perform two steps: (1) use the measured angles to the identified landmarks and the knowledge of the sensor movement to obtain an unbiased estimate of the landmarks positions in the fixed frame, and (2) calibrate the distance measurement model using these estimated landmarks positions.
We moreover validate the strategy using two approaches: (a) via simulated datasets, to analyze the limitations of the proposed procedure in a Monte Carlo (MC) fashion, and (b) to quantitatively assess the performance of the proposed procedure in real-life scenarios via field datasets recorded in a lab equipped with high-fidelity motion capture system.

Organization of the Manuscript
Section 2 formulates the calibration problem from a mathematical perspective. Section 3 describes the proposed algorithm, while Section 4 quantitatively assesses its performance. Section 5 concludes by listing the most important discoveries and research questions opened in the process.

Problem Formulation
We pose the following assumptions.
(A1) The environment from which we collect the measurements to be used for the calibration process has particular structures that produce easily recognizable features in the sensor readings. For example, the situation is as in Figure 3, where corners and poles produce clear features in the 2D plane of the measurements. Note that this means that our strategy cannot work in environments that miss these easily recognizable structures (such as natural places like deserts, or flat areas without trees). However, generally we consider applications where robots shall move precisely in the surroundings, and this calls for objects to be avoided. If there are no such obstacles/structures then the need for precise calibration becomes feeble. Given this, without loss of generality we require static and detectable landmarks; in this paper, we will use cylinders with known radius, but it could be anything as long as we have a detector for it. (A2) The sensor measurements lie in a 2D plane that is parallel to the ground. Moreover, the objects that produce the above mentioned features develop orthogonally w.r.t. the ground. This implies that the distances measurements are not affected by tilt effects. This requirement may not hold in generic situations; however, our envisioned calibration strategy is to be carried out within a building, where the conditions above hold. The problem of removing these assumptions is considered as a potential future extension. (A3) The statistical model underlying the distance readings contains heteroskedastic noise (for which the variance of the noise increases with the actual distance that shall be measured) and a bias whose amplitude also increases with the distance above. More specifically, we will focus on the situation where there exist l = 1, . . . , L objects in the environment, and k = 1, . . . , K places where the sensor can be placed. We then let [x l , y l ] and [ x k , y k ] be, respectively, the Cartesian coordinates of the L objects and of the K sensor positions. Accordingly, the actual distance between the sensor position k and the object placement l is We then assume that the distance readings are distributed as the polynomial model with e l,k ∼ N (0, 1) iid. The model parameters are thus {α i }, {β i }, with n being the corresponding model order (hereafter assumed for simplicity equal for both the bias and noise terms, i.e., n = n b = n h ). Note that in the following we may also use a simplified distance model that, for the sake of numerical tractability, neglects the heteroskedastic term in (3) so that the model reduces to We will refer to this model as to the "simplified distance model". (A4) Finally, we also assume that the angular readings θ l,k are noisy measurements of the actual angles θ l,k from which the object l is seen by the sensor at position k with respect to the reference frame of the horizontal axis. More precisely, we assume where the measurement noise is ν l,k ∼ N 0, σ 2 θ iid. Note that in practice this is a simplificative assumption that we use for analytical tractability reasons and that, a posteriori, is motivated by the numerical results we got during our experiments (For the sake of precision, it would be more formally correct to model the angle measurement noise through a Von Mises distribution with circular mean and non-circular concentration parameter. However, such a distribution converges to a normal one as the concentration parameter grows larger. In our case thus the approximation is justified in practice). Note, moreover, that this assumption implies that σ 2 θ is an unknown parameter of the model. We also assume that the error characteristics of (3)-(5) are time-invariant and do not depend on the absolute positions of the landmark (while they obviously depend on the relative distances "sensor vs. obstacle"). The angle θ l,k is thus the sum of the angle from which the object l is seen by the sensor with respect to the robot reference frame, plus the robot heading angle plus the rotation angle of the lidar's internal coordinate system with the robot coordinate system which is assumed to be a known constant. Note also that the measurement noise ν l,k in (5) incorporates imprecisions in the knowledge of the robot heading and rotation angle.
Summarizing, the calibration procedure shall return a reasonable model order n and an estimated parameter vector Θ = α 1 , . . . , α n , β 1 , . . . , β n , σ 2 θ . The problems are thus (P1) design a statistically optimal or near-optimal (in the Mean Squared Error (MSE) sense) algorithm that can be computed using closed-form expressions, and that can simultaneously estimate: the sensor coordinates x k , y k for each sampling position k, the position of the objects x l , y l for each object l, the model order n and the model parameters vector Θ above; (P2) quantitatively characterize the statistical performance of these estimators using appropriate mathematical analysis and field tests.   For illustration purposes and to make the paper self-contained, we now present a simple 225 landmarks position estimation algorithm and show that the method still works, even though better 226 landmark detectors may be available. We also note that the focus of the paper is on sensor calibration.

227
However, the same method developed for the sensor calibration part can be used to build a simple

A Triangulateration Strategy for Calibrating Distance Sensors
Optimally solving the problem 2 above requires jointly solving a nonlinear parameter estimation and a nonlinear state estimation problem. The solution is in general not available in closed form, and a viable numerical strategy could be using Monte Carlo techniques. However, this would require extensively long simulations and high processing power, which we assume is not available or usable. Recall indeed our initial idea: our goal is to develop strategies that can be used to endow cheap embedded systems (such as vacuum cleaning robots) to autonomously self-calibrate their distance sensors when desired and without the need to connect to external computing infrastructure. Therefore, we proceed to solve the problem using an ad hoc strategy that is easily implementable on normal embedded systems at the cost of sacrificing optimality of the estimates in the MSE sense.
More precisely, we propose to construct an estimator that computes solutions performing the following steps: 1.
Assume to know that there exist L landmarks, and to be able to identify and label them at each time instant from the raw measurements stream; 2.
place the sensor in a finite number of ideally equally spaced positions along an ideally straight line (say s k where k = 1, . . . , K); 3.
collect noisy measurements of the angles θ l,k and distances d l,k between the sensor and the various landmarks l = 1, . . . , L at each position s k (we recall that the stochastic models for these processes are (3) and (5)); 4.
estimate the 2D positions of the L landmarks in the inertial frame based on the sensor angle measurements θ l,k only, using the strategy defined in Section 3.1 below; and 5.
given the estimated landmark positions above, and the measured distances d l,k , estimate the model order and model parameters (i.e., do the actual sensor calibration step) with the strategy proposed in Section 3.2 below.

Estimating the 2D Positions of Circular Landmarks
For illustration purposes and to make the paper self-contained, we now present a simple landmark position estimation algorithm and show that the method still works, even though better landmark detectors may be available. We also note that the focus of the paper is on sensor calibration. However, the same method developed for the sensor calibration part can be used to build a simple landmark position estimation algorithm. Recall then that we assume the landmarks l = 1, . . . , L to be circular geometrical features in the measurement stream that are induced by distinct circular objects in the sensor environment as in Figure 2. To contextualize this assumption, consider Figures 3 and 4 and their captions, showing the sensor mounted on a robot moving in between the various circular landmarks. We used a lidar sensor from Neato. (It measure ranges from 0.2 m to 6 m and cover full 360 degree planar scan with angular resolution of 1 degree, [42]). (Neato Robotics www.neatorobotics.com).
Assume then, as in step 2 above, to place the sensor in k = 1, . . . , K ideally equally spaced positions along an ideally straight line, and to collect the corresponding raw measurements from the sensor. Note that it is possible to do the calibration without moving in a straight line, but it was chosen so to prevent the error from increasing largely during rotation and also to keep motion model (8) as simple as possible by excluding the robot heading angle. We selected equally spaced positions to minimize the error propagation from rover controller to the calibration process, as different step sizes might have different error levels in the controller. The next step is to compute, starting from this raw data, an estimate of the center of each circular landmark l using the information obtained at each sensor position k. Intuitively, we estimate the center of the landmark l as the point that minimizes the sum of its distances with the lines obtained at each k pointing towards the landmark, as shown in Figure 5 and described more formally in its caption. Given that in this paper we assume that the sensor sees round landmarks, this means that we 246 implicitly assume the presence of an offset in the distance measurements that is equal to the landmark 247 radius. For brevity, here we assume to know this parameter. Otherwise, circular landmarks lead to 248 raw measurements like the ones in Figure 2, from which it is not difficult to obtain practically accurate 249 estimates of such radii.

250
To summarize, thus, we assume that the robot moves along a straight line and that the sensor Given that in this paper we assume that the sensor sees round landmarks, this means that we implicitly assume the presence of an offset in the distance measurements that is equal to the landmark radius. For brevity, here we assume to know this parameter. Otherwise, circular landmarks lead to raw measurements like the ones in Figure 2, from which it is not difficult to obtain practically accurate estimates of such radii.
Therefore, to summarize, we assume that the robot moves along a straight line and that the sensor takes measurements in equally spaced positions along this line. For each landmark we can compute the straight lines that aim from all the various sensor positions to the (individually estimated) landmark centers. If there was no error, we could find the center of each landmark simply by intersecting the various lines referring to the same landmark. As we are in practice always far from this ideal condition (i.e., do not have accurate information neither about the sensor positions nor the measurement angles), the proposal is to find an estimate of each landmark center, say x l , y l , with a least-squares solution that minimizes the sum of perpendicular distances from the unique solution point to all these lines.
The remainder of this section is then dedicated to deriving the analytical structure of such estimator. To help readability, the notation is general so that l means "landmark", s k means "sensor position", x and y relate to the Cartesian reference frame, and refer to respectively the "ideal" and "noisy" versions of the same quantity. For example, the relation θ l,k = θ l,k + ν l,k indicates that the measured angle θ l,k of landmark l w.r.t. the position s k is a noisy version of the actual angle θ l,k . To find the estimated landmark position x l , y l in a closed form, consider then that in its k-th position the sensor is located at s k = [ x k y k ] T (importantly, a quantity that is unknown to the system). Ideally, the (k + 1)-th position of the sensor should be along a straight line (i.e., a fixed heading angle) and at a fixed distance, i.e., be where δx k and δy k are either zeros when the robot is not moving or constants determined by the step size when the robot is moving. In practice, though, the ideal conditions are not satisfied. We thus model the actual sensor position to be where e x,k and e y,k are zero mean Gaussian iid with the same variance σ 2 s . Note that this Gaussianity assumption is once again instrumental for the purpose of being able to devise computationally efficient schemes that can be implemented in embedded systems. Note, moreover, that we are recording the sensor measurements after the robot reached its new position (i.e., we are not considering measurements recorded during transients). As the actual sensor positions s k are not available, the best we have is only the reference (noiseless) sensor positions [x k y k ] T which can be determined using (7).
Moreover, consider the actual sensor position s k and the measured angle θ l,k of landmark l w.r.t. the position s k (line whose slope is then tan( θ l,k ), as the dotted lines in Figure 5). The equation of this line is then sin( θ l,k )x l − cos( θ l,k )y l − sin( θ l,k ) x k + cos( θ l,k ) y k = 0.
Moving the stochastic terms to the right hand side of the equation, we then obtain sin(θ l,k )x l − cos(θ l,k )y l − sin(θ l,k ) x k + cos(θ l,k ) y k = g k tan(−ν l,k ) where g k := sin(θ l,k )(y l − y k ) + cos(θ l,k )(x l − x k ).
Now substituting the geometrical identities sin(θ l,k ) = (y l − y k ) d l,k and cos(θ l,k ) = (x l − x k ) d l,k (immediately proved upon inspecting Figure 6) into (15), and then simplifying leads to n December 18, 2020 submitted to Sensors re g k := sin(θ l,k ) (y l − y k ) + cos(θ l,k ) (x l − x k ) . Now substituting the geometrical identities sin(θ l,k ) = (y l − y k ) d l,k and cos(θ l,k ) = ediately proved upon inspecting Figure 6) into (15), and then simplifying leads to ll that d l,k is the ground truth for the sensor-to-landmark distances -a ground truth t lable. Consider then that our main goal is to calibrate the sensor without using such gro Recall that d l,k is the ground truth for the sensor-to-landmark distances-a ground truth that is not available. Consider then that our main goal is to calibrate the sensor without using such ground truth information. Therefore, the best we can do instead is to plug in the measured distances d l,k as estimates (or the sample mean if more than one measurement is available at the same sensor position).
For small values of ν l,k (i.e., for the case where the dotted lines in Figure 5 aim decently at their target) we can simplify (14) using the approximations sin(ν l,k ) ν l,k cos(ν l,k ) 1.
Recall then that in our assumptions the measurement noise e l,k in (3) is statistically independent of the measurement noise ν l,k in (5). For that reason, the residual error in (17) will be heteroskedastic since the variance is σ 2 θ multiplied by the heteroskedastic variance of d l,k (that, we recall, is different in each sensor position). In the ideal case of no measurement noise, the line should pass through the center of the landmark l. In this special case, the point of intersection between these lines will be the position of the landmark center. As mentioned above, the presence of the measurement noise will however make the lines drift. The K lines corresponding to the K sensor positions in general will not intersect in a unique point, but in pairs. In this case, we may then solve the problem in a least-squares sense: the idea is to minimize the weighted sum of the squared distances (the solid red lines in Figure 5), i.e., to solve where Consequently, solving this system according to Aitken's generalized least square method [43] gives the following Best Linear Unbiased Estimator (BLUE) of the landmark centers, The computations above assume the full knowledge of the sensor positions s k = [ x k y k ]. As this information is not available, the best we can do is to plug in, instead, the expected sensor positions [x k y k ] in (7). Replacing [ x k y k ] with [x k y k ] in the least squares problem (18) thus leads to the problem x l y l = arg min where W * l := W l (as explaind below Equation (25) which in turns gives the weighted least squares estimator whose weights matrix W * is motivated below and defined in (25). This estimator is solvable, in the sense that the embedded system has all the information necessary to compute and minimize the cost. In other words, solving (21) is computationally feasible as all the required information is available while solving (18) is not. However, solving (21) leads to solving an approximate model of the intersection problem which will indeed result in a biased estimator, as it corresponds to solving the system of equations that is obtained after substituting x k and y k from (8) into (17), i.e., sin(θ l,k )x l − cos(θ l,k )y l − sin(θ l,k )x k + cos(θ l,k )y k sin(θ l,k )e x,k − cos(θ l,k )e y,k + d l,k ν l,k .
To characterize the error of this estimator, we notice that the residual error includes two different terms: the first is homoskedastic (specifically, corresponding to the first two terms in (24)) and with a variance of var sin(θ l,k )e x,k − cos(θ l,k )e y,k = sin(θ l,k ) 2 σ 2 s + cos(θ 2 l,k )σ 2 The second term is heteroskedastic with variance σ 2 θ var d l,k (specifically, corresponding to the last terms in (24)). Consequently, we suggest to set the weighting matrix W * l as W * l := σ −2 s I + σ −2 θ W l . (25) Note that in our assumptions both σ 2 s and σ 2 θ are assumed unknown constants. However, minimizing the sum of the squared residual errors in (24) is equivalent to minimizing the sum of the squared residual errors d l,k ν l,k , as the transformation between these errors is affine. This means that replacing W * l with W l will give exactly the BLUE for the parameters in (24).

Calibrating the Sensor
Once all the L landmark positions are estimated as x l , y l as in the previous section, we can easily estimate the various landmark-sensor position distances d l,k simply through computing the distance between each landmark with all the sensor positions and subtracting the radius of the landmark (which, as we said above, is either assumed to be known or assumed to be inferrable from the raw data). For notational compactness, define the (KL × 1)-dimensional distance measurement vector and rewriting (4) through a Vandermonde matrix, i.e., where the Vandermonde matrix Φ is of size KL × n + 1 and the parameter vector of size n + 1 × 1. Given this notation, the calibration procedure consists of three phases: • phase#1: model parameters estimation. After obtaining the estimates of the distances between the sensor and landmarks, estimate the parameters α casting the problem as a linear regression on (26) and the measurement vector d for model orders n = 0, 1, 2, . . . , n max , where n max is a user-defined parameter. This means solving for each potential n the problem which has the closed-form solution Note that, once again, the estimator α is unbiased; however, due to the simplification of the noise term in (3) (i.e., ignoring the heteroskedastic part of the noise), α will not be efficient. • phase#2: model order selection. We note that there exist various alternatives for selecting the optimal model order n ∈ {0, 1, 2, . . . , n max }: fitting opportune test sets, using crossvalidation, or also using model order selection criteria, for example, AIC. In the setups we considered for this paper we actually found that the model order selection problem has quite clear solutions, implying that all the various alternatives clearly indicated the very same number (see Section 4), implying in its turn that for our specific case all the various approaches tend to give equivalent results. It may, however, be that in other cases different strategies lead to different results; • phase#3: filtering new measurements. Once the model order selection and the model parameters estimation problems are solved, this means rewriting the "object distance vs. sensor reading" measurement model (3) as where d is the raw measurement, and d is the actual distance. To estimate d from d and the trained model one should thus invert (29). This inversion is not immediate; for example, one may solve the Least Squares (LS)-type optimization problem which requires finding the roots of a polynomial of order 2 n. Thus, despite its apparent simplicity, the problem of finding polynomial roots requires numerical methods for polynomial orders greater than 3.

Numerical Results
In this section, we verify in an empirical way the performances of the proposed calibration procedure, first with simulations using Matlab® and then with laboratory experiments using real sensors and landmarks in an environment endowed with a localization infrastructure that can return accurate assumingly ground truth information.
In general terms we thus consider d j , j = 1, . . . , J raw measurements from a noncalibrated sensor. To each raw distance measurement d j there corresponds also a true distance d j , and a filtered distance d j , i.e., the corresponding filtered version of these raw data.
As for the statistical performance index, the goal is to assess if and how much the calibration algorithm is actually leading to improved estimates of the distances, i.e., whether d j is statistically closer to d j than d j , and if so how much. To do so we use the MSE, i.e.,

MSE(a, b)
More precisely we will compute the ratio between the MSE computed with the raw data (distances measured by the sensor) and the MSE computed with the estimated data (distances estimated by the algorithm), i.e., use the Thus, in order to get an improvement with the estimation, this ratio has to be greater than 1.

Analyzing the Statistical Properties of the Landmark Position Estimator through Simulation Results
Unless otherwise stated, for all our simulations' plots, each point is the average of 1000 simulations of five landmarks and 20 sensor positions.
We start with statistically characterizing the landmark position estimator described in Section 3.1 by simulating a measurement model of the type (3) characterized by n = 2, α = [0.0525 , 0.8838 , 0.0584] and β = 0.05α, values that seem representing typical distance measurement systems mounted in modern autonomous vacuum systems.
We then investigate the MSE of the landmark position estimation procedure by analyzing how its bias and variance depend on three specific quantities: the standard deviation σ θ associated to the uncertainty of the sensor-to-landmark angle measurement in (6), 2.
the standard deviation σ s associated to the uncertainty in the sensor position evolution in (8), and 3.
the total number of landmarks L present in the scene.
The results are summarized in Figure 7, plotting the dependencies on σ s for a set of given σ θ and L, and in Figure 8, plotting the dependencies on L given σ s and σ θ .
In words, the results shown in Figures 7 and 8 confirm the obvious intuition that the smaller the noises, the better the estimator. However, we also note that, from numerical standpoints, it seems that guaranteeing σ s < 10 cm is important, and that guaranteeing σ s < 1 cm is instead not a necessity. This result is of practical importance, because the assumption that the sensor will be placed in a perfectly straight line will always be violated. However it seems that, at least for standard "domestic" cases like autonomous vacuum cleaners, violating this assumption will not disrupt the final results. We also note that Figure 8 suggests to set L to be around 5, i.e., an environment that is sufficiently rich while not being cluttered. data-filtering algorithm. We analyse how the MSE ratio (31) depends on the standard deviation σ θ of 331 the sensor-to-landmark angle measurement error in (6), and the standard deviation σ s of the sensor 332 position evolution uncertainty in (8).

333
The results are summarized in Figure 9, and show that the overall approach seems to be robust:

334
increasing gradually the standard deviations σ θ and σ s does not lead to abrupt decays of the overall 335 statistical performance. Moreover for values of σ θ and σ s that are meaningful in autonomous vacuum 336 cleaners situations we note MSE ratios that may reach 100. 337 We also remark that the overall strategy seems robust in its model order selection step. More 338 precisely, in all our simulations we selected a model order n = 2, which is the value we obtained in our 339 previous work [40] while calibrating the same lidar sensor from field data (a value that is numerically 340 convenient also because n = 2 leads to closed-form solutions for (30)). In the simulations considered 341 in this section, the overall estimation approach leads to estimating the model order n as the correct one, i.e., 2, when the standard deviations σ θ and σ s are reasonably low. However, as the noises increase, we 343 noted that the order selection process tends to become more and more conservative, and select the 344 simpler model n = 1 (see graphically Figure 10).

Analyzing the Statistical Properties of the Sensor Calibration Procedure through Simulation Results
We then pass to the second part of the estimation procedure, i.e., calibrating the model parameter that we presented in Section 3.2. Recall that the calibration algorithm is based on the estimated landmarks position, i.e., there is the need to estimate the landmarks positions first, to then proceed to the calibration step. We then define as main performance index the MSE ratio (31), i.e., a measure of how much worse the raw data measured by the sensor is w.r.t. the distances estimated by the data filtering algorithm. We analyze how the MSE ratio (31) depends on the standard deviation σ θ of the sensor-to-landmark angle measurement error in (6), and the standard deviation σ s of the sensor position evolution uncertainty in (8).
The results are summarized in Figure 9, and they show that the overall approach seems to be robust: increasing gradually the standard deviations σ θ and σ s does not lead to abrupt decays of the overall statistical performance. Moreover, for values of σ θ and σ s that are meaningful in autonomous vacuum cleaners situations, we note MSE ratios that may reach 100.  measurements are for all the practical purposes noiseless and considerable as ground truth. 355 We then apply both the landmarks position estimator in Section 3.1 and the consequent parameters 356 calibration procedure in Section 3.2 to calibrate the triangulation lidar shown in Figure 3. We recall that 357 this type of lidar is not really accurate, as its measurements are affected by both a systematic bias and 358 a heteroskedastic variance (see Figure 1) that lead to increasing measurement errors as the measured 359 distance increases.   We also remark that the overall strategy seems robust in its model order selection step. More precisely, in all our simulations we selected a model order n = 2, which is the value we obtained in our previous work [40] while calibrating the same lidar sensor from field data (a value that is numerically convenient also because n = 2 leads to closed-form solutions for (30)). In the simulations considered in this section, the overall estimation approach leads to estimating the model order n as the correct one, i.e., 2, when the standard deviations σ θ and σ s are reasonably low. However, as the noises increase, we noted that the order selection process tends to become more and more conservative, and select the simpler model n = 1 (see graphically Figure 10). Finally, we again investigate the effect of using different numbers of landmarks on the whole calibration process. The results, summarized in Figure 11, show again that increasing the number of landmarks from one to three leads to noticeable improvements in the MSE ratio. However, increasing the number of landmarks further will not lead to further improvements while, at the same time, increasing the computation complexity. proper calibration of the sensor we need to have a "sufficient" calibration dataset. In general, for range 370 Figure 11. Dependency of the MSE ratio on the number of landmarks L as a function of the sensor position standard deviation σ s for the case σ θ = 2. Decreasing σ s is, as expected, always beneficial, while increasing L is not so important.

Field Experiments
We consider field experiments in a laboratory provided with a Vicon motion capture system that uses triangulation to compute the position of the objects inside the laboratory. Such a Vicon system is very accurate, compared to the sensors we aim at calibrating. For this reason, we assume that the Vicon measurements are for all the practical purposes noiseless and considerable as ground truth.
We then apply both the landmarks position estimator in Section 3.1 and the consequent parameters calibration procedure in Section 3.2 to calibrate the triangulation lidar shown in Figure 4. We recall that this type of lidar is not really accurate, as its measurements are affected by both a systematic bias and a heteroskedastic variance (see Figure 1) that lead to increasing measurement errors as the measured distance increases.
For practical purposes we placed the lidar sensor on top of a Pioneer 3AT mobile robot, as shown in Figure 4, controlled through a computer using Robot Operating System (ROS). We also consider five hand-made cylindrical landmarks with a radius of 12 cm, scattered within the field of view of the Vicon system and in a way that all of them are always visible and distinguishable by the sensor from all the positions [x k , y k ] from where it will take measurements. As a practical indication, because of the intrinsic limits of the considered sensor, each landmark has to be placed not farther than 5 meters and closer than 20 cm from all the sensor positions. We then programmed the robot to move on a straight line path, and oriented it so to not hit the landmarks while moving. More precisely, we programmed the robot to move and stop 10 times, doing each time an incremental step of 30 cm. In order to have proper calibration of the sensor we need to have a "sufficient" calibration dataset. In general, for range sensors, a sufficient dataset should cover all the sensor ranges of interest. During our field experiments we moved the sensor in a straight line for about 3 m to ensure the richness of the recorded dataset. Figure 3 shows a photo of one of the experiments. We repeated this type of experiment with three different placements, so to take three datasets of Vicon measurements of distances and angles (i.e., ground truth) of the five landmarks for all the sensor positions. In other words, thanks to the Vicon system we were able to compute all the actual distances and angles between all the various landmarks and the sensor in its various positions.
We then split each dataset in three parts: the first two to be used as a training set (the first third to estimate the sensor parameters in (3) and the second third to choose the model order n) and the third part to be used as a test set. The field results presented for landmarks number other than 5 is basically obtained with the same recorded dataset of five landmarks after removing the receded data of extra land marks. For example, the 3-landmark datasets, is the 5-landmark datasets after removing the recorded data associated with the last two landmarks, and the 4-landmark datasets is the same as the 5-landmark datasets after removing the data of the fifth landmark, and so on. The obtained results, summarized in Figure 12, show again that the improvements change as more landmarks are involved in the calibration process. In other words, the field results are in good agreement with the simulated ones. However, a few outliers still exist which might be due to increased noise variance in one of the unconsidered processes like the landmark association problem. Moreover, in all the calibrations we performed on the different datasets, we obtain a selected model order n = 1, which indicates, based on our simulations, that there may be a high noise variance associated to the noise in measuring the angle with the landmarks.
Finally, Figure 13 shows how the proposed calibration procedure can help in real-life situations by reporting the measurement process relative to the third placement considered in Figure 12. Here, the landmarks' borders are plotted as gray circles, the series of actual positions of the sensor as a stripe of blue dots. Moreover, the red crosses plot the raw measurements d l,k obtained by the sensor when observing the various landmarks, while the green dots plot the filtered measurements d l,k obtained by applying the proposed calibration and filtering algorithms. One may note how the d l,k 's capture the actual positions of the various landmarks in a qualitatively much more precise way.  Through field experiments we saw that the overall proposed calibration approach may be quite 407 robust: even if one does not get results that are as good as the ones achievable using external ground 408 truth systems, our algorithm has been able to lead to a reduction of the norm of the measurement 409 errors between the pre-calibration raw data and the post-calibration ones by a factor 10; in comparison, 410 using ground truth calibration as in [40] led to a reduction factor of 17 (thus better but not of orders of 411 magnitude, and at the cost of having to buy, set up and use a ground truth collection system). 412 We thus remark that this factor 10 achievement is through using just software logic and 413 assumptions on the landmarks being fixed, and no additional hardware nor special conditions. In  Figure 13. Example of the effects of the proposed calibration procedure on a field experiment. The actual sensor positions are plotted with a series of practically aligned blue dots, the true cylindrical landmarks in gray circles, the raw measurements taken by the sensor with red crosses, and the filtered distances, computed using the proposed strategy, in green circles. Ideally, the measurements should lie on the borders of the landmarks. It is immediately noticeable how the calibrated measurements lie much closer to such borders than the non-calibrated ones.

Author Contributions:
The authors are equally contributed in conceptualization, methodology, validation, 420 writing-review and editing.  Figure 13. Example of the effects of the proposed calibration procedure on a field experiment. The actual sensor positions are plotted with a series of practically aligned blue dots, the true cylindrical landmarks in gray circles, the raw measurements taken by the sensor with red crosses, and the filtered distances, computed using the proposed strategy, in green circles. Ideally, the measurements should lie on the borders of the landmarks. It is immediately noticeable how the calibrated measurements lie much closer to such borders than the non-calibrated ones.

Conclusions
The nonlinear and heteroskedastic model of distance sensors can be calibrated by exploiting just the structure of a fixed environment. In other words, if the environment presents some particular features that may be used as generic landmarks, then one may use the fact that the landmarks do not move to infer the own movement. This means the possibility of estimating the landmarks' positions minimizing opportune cost functions, and in this way obtain information useful to learn the characteristics of a distance sensor without the need for external distance measuring devices to be used as providers of ground truth information.
Through field experiments we saw that the overall proposed calibration approach may be quite robust: even if one does not get results that are as good as the ones achievable using external ground truth systems, our algorithm has been able to lead to a reduction of the norm of the measurement errors between the precalibration raw data and the postcalibration ones by a factor 10; in comparison, using ground truth calibration as in [40] led to a reduction factor of 17 (thus better but not of orders of magnitude, and at the cost of having to buy, set up and use a ground truth collection system).
We thus remark that this factor 10 achievement is through using just software logic and assumptions on the landmarks being fixed, and no additional hardware nor special conditions. In conclusion, the here proposed calibration procedure is expected to lessen the time to prepare the calibration setup, and is expected to be implementable well beyond laboratory setups.
We though recall that one standing assumption we exploited is that sensor measurements lie in a 2D plane that is parallel to the ground. As this requirement may not hold in some practical situations, we devise this as the most important future research direction spanned by the current work.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available at the current time.