Joint Temperature-Lasing Mode Compensation for Time-of-Flight LiDAR Sensors

We propose an expectation maximization (EM) strategy for improving the precision of time of flight (ToF) light detection and ranging (LiDAR) scanners. The novel algorithm statistically accounts not only for the bias induced by temperature changes in the laser diode, but also for the multi-modality of the measurement noises that is induced by mode-hopping effects. Instrumental to the proposed EM algorithm, we also describe a general thermal dynamics model that can be learned either from just input-output data or from a combination of simple temperature experiments and information from the laser’s datasheet. We test the strategy on a SICK LMS 200 device and improve its average absolute error by a factor of three.


Introduction
ToF LiDARs estimate distances by emitting short bursts of laser light and by measuring the time it takes for the reflected photons to arrive back to the device [1]. Despite being based on a very simple principle, they are very much both accurate and precise devices [2]: for example, precisions can reach 10 mm of standard error when the object is 10 m away. Due to these favorable properties, they are commonly used in critical industrial applications where there is the need for high quality measurements.
It is well known that these devices need temperature compensation mechanisms, since changing their temperature leads to changes in the statistics of the returned measurements. The effect of temperature may be huge: experiments by [3] on an amplitude-modulated continuous-wave laser radar pointing at a target six meters away from the sensor showed that measurements at 21°C and measurements at 45°C were differing by 40 cm. Since thermal stabilization of a laser scanner may take up to 30 min [4], it is clear that these sensors are affected by a warm-up-induced time drift that must be compensated. Manufacturers of ToF devices thus usually embed opportune algorithms in their products that implement this temperature compensation mechanism.
Unfortunately, temperature is not the only physical factor that deserves compensation: as described in detail in Section 2, lasers can suddenly change their lasing mode. This property, called the mode-hopping effect, has a substantial impact on the measure returned by ToF devices, since changing lasing mode means to change the spectral content of the laser burst, i.e., change its time of flight. Remarkably, to the best of our knowledge, the existing literature does not focus on managing this effect, but rather, considers only temperature compensation mechanisms.
We would like to mention here that there are also other methods with high temperature compensation at the nano-scale measurement, such as [5]. The method reduces the offset, the temperature characteristic of the main sensing element, the temperature drift and the noise by the switching method.

Statement of Contributions
We propose an expectation maximization (EM) algorithm that compensates mode-hopping effects by modeling the induced measurement noise as a Gaussian mixture. Thus, from the mathematical perspective, we introduce some latent variables (namely, from which Gaussian the noise comes from) as additional estimands. This EM algorithm is also coupled to a temperature compensation filter that is built on a physics-based linear model for the thermal dynamics of the laser scanner. Summarizing, thus, our contributions are, with respect to the aforementioned literature:

1.
A thorough motivation for why it is meaningful to consider mode-hopping effects in laser scanners, arising from a physical description of the lasing mechanism in laser diodes; 2.
A thermodynamical model for the thermal dynamics of a whole laser scanner, needed by the proposed strategy to account for temperature effects; 3.
A statistical model describing the measurement process that decouples the effects of the mode-hopping from temperature effects; 4.
A numerically-efficient EM strategy based on the statistical model above; 5.
A validation of the proposed compensation strategy on real devices.
With the validation, we also show that it is possible to improve the absolute error of a SICK LMS 200 device by a factor of three.

Organization of the Manuscript
Section 2 analyzes the effects of the laser temperature on the measured distance. Section 3 proposes a general model for the thermal dynamics of a pulsed ToF LiDAR. Section 4 presents a general measurement model accounting for both temperature and mode-hopping effects. Section 5 describes how to train the EM algorithm, while Section 6 describes how to use the same algorithm for testing purposes. Section 7 then proposes a likelihood ratio test for calibrating the hyperparameters of the EM algorithm. Section 8 presents some numerical results on commercial devices. Finally, Section 9 draws some conclusions and future research directions.

Effects of the Laser Temperature on the Measured Distance
This section lays down interpretations motivating the structure of the novel compensation procedures. We thus here describe the functioning principle of ToF scanners, explain why the measured distance depends on the temperature of the device and motivate why the measurement noise of a LiDAR is intrinsically multi-modal.
Consider then the basic operation of pulsed ToF LiDARs in Figure 1 and in its caption. The measurement of the distance derives from ideal considerations: if the temporal width of the pulse is null, then the distance d between the sensor and the object should satisfy: where c is the speed of light and τ is the measured ToF between when the laser pulse is emitted and when it is received. Assume ideally that the laser pulse contains photons with a unique nominal wavelength λ 0 . Since: with f the light frequency and λ 0 the nominal wavelength, and since the light frequency f remains the same through different media, to know c, it is sufficient to know λ 0 . Thus, from knowing λ 0 , one can compute d, since τ is measured. Graphical description of the operating principle for pulsed ToF LiDARs. A pulsed infra-red laser beam is first emitted from the transmitter. The case of the transmitter, in dark gray, encloses a laser junction and a laser cavity. The emitted laser beam is then deflected by a rotating mirror (resulting in a fan-shaped scan pattern) and, finally, reflected from the object surface. The time of flight τ between the transmission and the reception of the laser beam is then used to estimate the distance d between the scanner and the object.
We can already now notice the first effect of the temperature of the device on the measurement: according to its datasheet, the nominal wavelength of the laser diode SPL PL90 from OSRAM [6] is 905 nm at 25°C, and 907.8 nm at 35°C. Assuming the target to be at a one-meter distance, this 2.8 nm variation in the wavelength λ then results in a 10.3-ps variation in the ToF τ, i.e., a variation in the measured distance of approximately 3.1 mm. The previous consideration is nonetheless simplistic. We can indeed notice another three distinct effects:

1.
In general, lasers do not emit at a unique frequency λ. Indeed, the average spectral distribution of the laser pulses follows a "comb"-like density, like the one in Figure 2.

2.
Lasers are affected by the so-called mode-hopping effect [8] and, indeed, oscillate between different lasing modes (the teeth of the "comb" of Figure 2), for which two different pulses generated under the same temperature and external conditions may have different λ's (e.g., referring to the same figure, the first pulse may contain only photons with wavelengths λ 1 , while the second pulse may contain only photons with wavelengths λ 2 ). In other words, the actual distribution of one specific pulse may contain only a subset of the teeth of the average spectral distribution. Using naively Equations (1) and (2) to estimate d without being aware of the mode hopping, i.e., assuming a certain λ 0 without actually knowing that the average λ jumps between different lasing modes, reflects thus in a multimodal measurement of d, as clearly shown in Figure 3.

3.
The average spectral distribution of the laser pulses is not fixed, but rather depends on the temperature of the transmitter [9]. More precisely, the positions and amplitudes of the modes in Figure 2 depend on both the current flowing through the laser junction and the geometry of the laser cavity, but eventually, these two effects are inter-combined: the current flow produces heat that will modify the geometry of the cavity. Eventually, thus, the temperature affects the position and amplitude of the modes of the average spectral distribution. This temperature effect can be clearly seen in Figure 4: even if the device is nominally already compensated in temperature, one can clearly see two different lasing modes shifting in temperature.    . Dependency of the distance measurements on the device temperature for a SICK 200 device pointing to a fixed object in a controlled environment. Despite the true distance and other parameters potentially affecting the measurements being constant in time, the distributions of the measurements are temperature varying. We can notice how the device compensates for the temperature change by adding a temperature-varying negative bias, but that it does not compensate the mode-hopping effect.
To summarize, the actual λ of a laser pulse is in general different from the nominal λ 0 because of two effects: first, the laser may oscillate between different modes; second, these positions of these modes vary with temperature.

A General Model for the Thermal Dynamics of a Pulsed ToF LiDAR
Instrumental for compensating the effects of temperature offsets highlighted in the previous section, we here discuss a general thermodynamic model describing the temperature dynamics of the whole laser scanner that will lead us to being able to represent the temperature changes of the laser cavity.

Physical Modeling
Just like any power electronic device, LiDARs generate heat that is then exchanged with the environment, so that the temperature of a scanner depends on the temperature of the environment. The main sources of heat inside the LiDAR are thus the laser diode, the motor and the electronic components of the system. The heat generated inside the scanner is then transferred to the surrounding environment through the case. Since our experience indicates that motors and other electronic components induce negligible thermal effects, we consider only the heat produced by the laser diode.
We thus represent the thermal model of a generic scanner as the equivalent electrical circuit shown in Figure 5, interpretable as follows: when the laser is turned on, the heat generated by the laser junction is dissipated in the surrounding environment through first the transmitter case and then, second, through the laser scanner case. Figure 5. The proposed thermodynamic model for a generic pulsed ToF LiDAR.
Considering the notation: P Heat power generated by the junction (equal to zero when the device is off) t j Temperature of the junction t c Temperature of the transmitter case t a Temperature of the external ambient t s Noisy measurement of the temperature of the transmitter case C 1 Thermal inertia of the transmitter case C 2 Thermal inertia of the laser scanner case R 1 Thermal resistance between the transmitter case and the laser scanner case R 2 Thermal resistance between the laser scanner case and the ambient It follows that the dynamics for the temperatures of the laser junction and scanner case when the laser is on are thus: Discretizing the previous dynamics with a discretization step of T s and letting: where t s k is a noisy temperature measurement of the case at time k, ν k and µ k are independent process and measurement noises, and where: Notice that through Equation (5), we introduce t s k , i.e., a noisy measurement of the temperature of the transmitter case. This correspond to the practical assumption that perfect knowledge of the actual case temperature t c k is in general unavailable, since temperature sensors attached to the case of the scanner will never give noiseless recordings.

Identifying Model Equation (5)
If one has a dataset of recorded temperatures t j k , t c k and t a k , then one may identify Equation (5) using standard system identification approaches, e.g., a prediction error method (PEM) as in [10]. For the common case where it may be difficult to obtain direct measurements of the quantities, we propose to resort to the following general strategy that uses the datasheet of the laser scanner in conjunction with noisy case temperature measurements t s k .
Algorithm 1 Identification of model Equation (5) starting from the datasheet of a laser scanner.
1: From the datasheet of the laser scanner infer: • its thermal resistance R 1 (directly from the datasheet); • its thermal capacity C 1 (by estimating the volume of the laser diode from the datasheet and multiplying it for the heat capacity of the material, also indicated in the datasheet); 2: From measuring the weight and material of the case, infer its thermal capacity C 2 ; 3: Estimate the generated heat power P by measuring the electrical power absorbed by the device and multiplying this quantity by 0.5 (for the estimated efficiency of generic laser diodes [11]); 4: From situations where the scanner is in thermal equilibrium, calculate R 2 by calculating the difference between the measured case temperature t s k and the ambient temperature t a k divided by the estimated generated heat power P.

Estimating t j from t s
Assume that model Equation (5) has been identified, either from measured data or using Algorithm 1, and observe that the thermal model is observable, reachable and has stable dynamics. Due to their favorable theoretical and numerical properties, we thus devise to estimate t j k from noisy measurements of the case temperature t s k via Kalman smoothers/filters, as in [12], i.e., to let: for an opportune (and potentially time varying) matrix W. As before, notice that the filtering starts from noisy measurements of the case temperature t s k , rather than from perfect measurements t c k .

A General Measurement Model Accounting for Temperature and Mode Hopping Effects
We recall that our aim is to model the effect of the temperature of the laser junction t j k on the measured distance d and understand how noisy case temperature measurements t s k help improving the accuracy on the final estimate of d. To this aim, we propose the following measurement model at the generic time instant k: where: • y k is the distance returned by the sensor; • d is the true distance from the object (assumed deterministic); • t j k is the temperature of the laser cavity at time k; • The two modes w 1 k ∼ N µ 1 , σ 2 1 and w 2 k ∼ N µ 2 , σ 2 2 account for a bimodal Gaussian and white additive measurement noise. The Bernoulli random variable (r.v.) ∆ k ∼ B (π) selects the active mode at time k, so that π reflects the relative importance of the modes. Intuitively, ∆ k represent which lasing mode has been active during measurement y k . We notice that here, we consider bimodal noises (i.e., only two lasing modes) just for notational simplicity. It is nonetheless immediate to generalize the subsequent findings for an M-modal case; • H (·) θ is a non-linear transformation of the temperature of the laser junction t j k in a measurement bias. In the following Examples 1 and 2, we show how different H (·)'s and θ's express different maps from the junction temperature t j k to the measurement bias.
Notice that, induced by our experience, we let the measurement noise modes w 1 k and w 2 k be independent from the laser junction temperature t j k . Our experiments indeed indicated that the moments of the noises are not affected by changing temperatures (at least for a range between 0°C and 40°C). Notice moreover that model Equation (8) is linear in θ; this restrictive assumption is nonetheless essential for building distance estimation algorithms that are numerically fast.
i.e., a measurement model where the temperature plays the role of an N-th order polynomial bias (see Figure 6). y k N = 2 N = 5   (11) and θ := θ 1 , θ 1 , . . . , θ N , θ N , where the fundamental frequency f 0 is assumed to be known. Then, the generic model Equation (8) specializes into: i.e., a measurement model where the temperature plays the role of a bias that is periodic with frequency f 0 (see Figure 7).  We then use model Equation (8) for practical purposes following the classical scheme: • Design the model, i.e., decide the structure for H (·) (e.g., between Example 1, Example 2 or also different ones depending on the need) starting from data collected in a controlled environment (see Section 7); • Train the model, i.e., estimate θ and the statistics of w 1 k , w 2 k and ∆ k from data collected in a controlled environment (see Section 5); • Test the model, i.e., use the previous estimated quantities during the normal operation of the laser scanner, so as to improve the estimation of d from data collected in a non-controlled environment (see Section 6).
The following sections are dedicated in detail to how to implement the previous three points. Notice that Section 7, on the design of the model, is presented after Sections 5 and 6 for notational convenience.

Training Model Equation (8)
In this section, we devise a numerical algorithm for learning θ (i.e., the coefficients multiplying H(·)), µ 1 , µ 2 , σ 2 1 , σ 2 2 , π (i.e., the statistics of the noises w 1 k , w 2 k and of the random variables ∆ k ) and the values of the mode selection variables ∆ k starting from a dataset containing: • y k for k ∈ {1, . . . , K}; • d, i.e., the real distance; • t s k , i.e., measurements of the temperature of the laser scanner case (to be transformed into estimates of t j k through Equation (7)).
Let then: Given the frequentist assumptions on the unknowns, we would like to perform maximum likelihood (ML) estimation for θ, i.e., seek for: Since the laser junction temperature t j is unavailable, Equation (14) cannot be solved. We thus resort to solving the approximated problem: where t j is the estimate of t j given t s as in Equation (7). We also notice that, instead of considering the joint density of y and t j , it is sufficient to consider the conditional density of y given t j : It is important to notice that the ML problem in Equation (16) contains the latent variables ∆; to estimate them, we thus resort to a tailored EM approach [13]. To this aim, define then the auxiliary variables: with 1 being a vector of K ones. The computation of θ, ∆ is performed through the iteration up to convergence (stopping criteria for EM algorithms are usually based on relative or absolute changes in the parameter estimate or in the value of the log likelihood; see, e.g., [14]; in our implementations, we used the absolute changes in the parameter estimate π < ; where is a small number 1 × 10 −5 in the implementation) of the following two steps:

M-step
In our numerical experiments, we empirically found it convenient to use the following initial conditions: As for the convergence of the EM to the true ML estimate, we notice that EM algorithms are not in general ensured to have convergence properties. A sufficient condition for convergence is in [15], where the authors show that EM algorithms are convergent if the maximizer of the M-step is unique (a condition that is almost always satisfied in practice). In our case, the M-step maximizer is unique as long as in the update for θ in Equation (19), the matrix H T C −1 H admits the inverse. In general, e.g., in both the polynomial case of Example 1 and in the Fourier expansion case of Example 2, this translates into the need for at least N samples associated with N different laser junction temperatures t j k .

Testing Model Equation (8)
We now devise a numerical algorithm for estimating d and the values of the lasing mode selection variables ∆ k starting from the model trained in Section 5 (i.e., an estimated vector θ := θ T , µ 1 , µ 2 , σ 2 1 , σ 2 2 , π T and the statistics of w 1 k , w 2 k and ∆ k ) and a set of measurements y k and t s k for k = 1, . . . , K.
Assuming once again to transform the temperature sensor measurements t s into estimated laser junction temperatures t j = Wt s through Equation (7), the problem of estimating d and the ∆ k s can be cast as: As before, we compute this ML estimate through an EM approach:

E-step
The same values of µ 1 , µ 2 , σ 2 1 and σ 2 2 at the end of the training step will be used in the testing step. There is no need to recompute them again. In our experiments, we found it beneficial to start from the initial conditions d = 1 T y K and ∆ k = 0.5. The convergence properties of this EM procedure are then very similar to the EM in Section 5.

Designing Model Equation (8)
This section is divided into two parts: 1.
Section 7.1, suggesting some hints for designing different structures for H (·) (e.g., choosing the order for Example 1, for Example 2 or also designing different functional structures depending on the collected information); 2.
Section 7.2, suggesting a numerical algorithm for discriminating between different competing structures for H starting from data collected in a controlled environment.

Designing H (·)
The proposed EM algorithms have the numerically-favorable property of having both the E and the M steps solvable in closed form. It is important to notice that this is induced by the fact that model Equation (8) is linear in θ, i.e., Thus, with this structure, the designer can model the effect of the temperature t j k on the measurement y k as the sum of N independent effects, each one represented as an opportune generic function of t j k (the weight of which is actually assumed unknown before the training phase). As shown in Figures 6 and 7, the structures proposed in Examples 1 and 2 have quite general generalization capabilities. Nonetheless, the designer can tailor H so that it resembles other structures; our suggestion is to start from raw measured data spanning different temperatures, check visually how the macroscopic temperature trend behaves and then decompose this trend as the sum of different functions that will become the various H n (·) in model Equation (8).

Determining the Best H (·) among Different Competing Potential Structures
The process described in Section 7.1 may lead to different competing structures for H (·). In other words, the designer may propose different structures H (1) , H (2) , etc., and would like to choose the "best" H (i) given a dataset containing y, t s and the true distance d.
We propose to use the classical approach of discriminating the various H (i) 's considering their goodness of fit, i.e., to use generalized likelihood ratios (GLRs) [16], for which we first estimate the best estimates θ (i) and ∆ (i) for each H (i) given the dataset and then select the best hypothesis considering their resulting log-likelihoods. More formally, the suggested procedure is as in Algorithm 2. Compute θ (i) and ∆ (i) as in Section 5; 3: where Σ

Experiments
We now validate the proposed thermal model and estimation strategy on real data, aiming to show their effectiveness. We thus consider the simple scenario where a target is fixed in front of a sensor, under constant light and electrical conditions, so that the actual distance between the sensor and the target d is fixed. We consider a SICK LMS 200 LiDAR, one of the most widely-used LiDARs in industry and robotics applications.

Training and Validation of the Thermal Model Equation (5)
To train and validate the thermal model Equation (5), we conducted the experiment summarized in Figure 8: in a thermally-controlled room at 22°, we performed several on-off cycles of the device and measured the corresponding case temperature t s k . We then used the data represented as "training data" in the figure, the data-sheet of the SICK LMS 200 LiDAR and Algorithm 1 to train model Equation (5). After that, we drove the trained model with the u k in the "test data" as inputs and obtained a predicted temperature y s k . The goodness of fit of the predicted temperatures, computed as: is then 92.79%. This indicates a very good fit, i.e., a good approximation capability of our proposed thermal model (and associated learning algorithm).

Selection of the Optimal H(·) for the SICK 200 LiDAR
We used Algorithm 2 on real data from a SICK LMS 200 to discriminate between different H (i) 's in Examples 1 and 2, i.e., the hypotheses: against the null hypothesis: so that the log likelihood ratio test between two hypotheses can be defined as: where we explicitly mention the dependence of the log likelihood from the various parameters of the considered model. Figures 9 and 10 show the likelihood ratios obtained for different candidates for polynomial and Fourier expansions, respectively. For the polynomial case, we obtained the highest likelihood for the polynomial order N = 16, while for the Fourier case, we obtained the optimal value for N = 180. We motivate this latter order, much higher than that of the polynomial case, to be due to the periodic nature of the Fourier model.

Assessment of the Performance Improvement for the SICK 200 LiDAR
Assume considering the best model order as in Section 8.2 and applying the main EM estimation algorithm for testing purposes. Figure 11 then shows an example of the outcome of this test on real data: first, the figure plots for each raw measurement the associated lasing modes detected by the EM testing procedure. The same plot shows also the bias induced by the temperature, i.e., H θ.
Once the raw measurements have been associated with the corresponding lasing modes, it is then possible to subtract from them the corresponding noise bias µ 1 or µ 2 ; in other words, it is possible to remove a certain lasing mode-dependent bias from each of the various measurements. This leads to a new distribution of the measurements, as shown in Figure 12. This figure shows the cumulative distribution of the measurements that we plotted in Figure 11 before and after subtracting the biases induced by the lasing modes. It is clearly visible that after this removal, the cumulative distribution becomes sharper, indicating that the novel dataset has a smaller variance. CDF after the EM algorithm CDF before the EM algorithm actual distance Figure 12. Effects of compensating some raw data with the biases corresponding to the lasing mode of each sample.

Convergence Properties for the EM Algorithms
We now report some numerical results for the convergence properties of the EM training and test algorithms. More precisely, Figure 13 shows the evolution of the estimates during the training phase, while Figure 14 shows the same evolution for the testing phase. We also plot in Figure 15 the computational burden of performing a fixed number of EM test steps for different measurement vector lengths and show empirically how the computational efforts for testing some points are quadratic with the number of samples. Importantly, this computational effort does not depend heavily on the order of the model.
We show in Figure 16 how the precision and accuracy of the estimate d changes with the number of samples in the test set. More precisely, we show the following cumulative distribution: given the set of 20,000 raw measurements, extract every subset of 200, 400 or 800 consecutive measurements and, for each of these subsets, compute an estimate d. The result is intuitive: estimating from more measurements leads to a better estimator (i.e., a cumulative distribution that steepens in correspondence of the true distance).  We show in Figure 17 the realizations of the errors committed by three different estimators: d (our estimator), y (the sample average of a window of measurements) and y − H θ (the sample average of a window of measurements corrected by the temperature compensation term H θ). We notice that d commits the smallest errors almost everywhere. To integrate the information, we also plot the temperature of the device in the lower plot. This helps with recognizing that d has better performance, especially where the effect of the temperature is higher.  We also plot in Figures 18 and 19 how the statistical moments of the empirical error and of the empirical absolute error of these estimators behave with the test sample size. Each point in the plot represents the mean or the variance of an error (or absolute error) sequence similar to those shown in Figure 17. We see that the performance of d improves monotonically, while for the other estimators, this does not happen. Similar results happen also in the plots for the variances, where the proposed estimator outperforms the other two.

Conclusions
Physical considerations on the mode-hopping effect lead to the consideration that the measurement noise of time of flight (ToF) laser scanners is intrinsically multi-modal. In its turn, this implies that estimating the actual distance between scanners and the surrounding objects should be performed using latent variable-based statistical models, where the latent variables correspond to the lasing modes of the laser. Since no literature seems to account for this multi-modality, we aimed at closing this gap.
We thus proposed an expectation maximization (EM) algorithm on top of a model that captures biases on the measured distance induced by temperature changes, plus mode-hopping effects through an opportune Gaussian mixture on the measurement noise. Importantly, thanks to a separable model, the EM iterations can be performed analytically. The computational advantages are clear: non-separable models may indeed need to perform the EM iterations numerically, and this would lead to a computational burden hindering the usability of mode-hopping correction procedures in on-line settings.
The proposed strategy incorporates an accurate model of the temperature dynamics of the laser diode. Moreover, to account for the fact that it may be difficult to collect temperature data on the laser, we proposed a strategy for the identification of the parameters of this model that exploit the datasheet of the laser device and some very simple experiments.
Overall, the proposed temperature compensation strategy led to diminishing the spread of the distribution of real measurements from a SICK LMS 200 around the true value d, as shown in Figure 12, with corresponding decays of the variance of the absolute error from 2.0 mm 2 to 0.68 mm 2 and shift of the expected absolute error from 3 mm to 1 mm, as shown in Figure 19.
We eventually notice that the proposed strategy is still in its infancy: indeed, we considered a frequentist case for which the actual distance d is a deterministic and fixed quantity, aiming at showing that it is possible to improve the overall precision of a laser scanner through accounting for mode-hopping effects. Nonetheless, in real scenarios, d will vary; we thus devise future efforts focusing on strategies for which d is a stochastic process that varies according to some given a priori distribution.