Using Least-Squares Residuals to Assess the Stochasticity of Measurements—Example: Terrestrial Laser Scanner and Surface Modeling

: Terrestrial laser scanners (TLS) capture a large number of 3D points rapidly, with high precision and spatial resolution. These scanners are used for applications as diverse as modeling architectural or engineering structures, but also high-resolution mapping of terrain. The noise of the observations cannot be assumed to be strictly corresponding to white noise: besides being heteroscedastic, correlations between observations are likely to appear due to the high scanning rate. Unfortunately, if the variance can sometimes be modeled based on physical or empirical considerations, the latter are more often neglected. Trustworthy knowledge is, however, mandatory to avoid the overestimation of the precision of the point cloud and, potentially, the non-detection of deformation between scans recorded at different epochs using statistical testing strategies. The TLS point clouds can be approximated with parametric surfaces, such as planes, using the Gauss– Helmert model, or the newly introduced T-splines surfaces. In both cases, the goal is to minimize the squared distance between the observations and the approximated surfaces in order to estimate parameters, such as normal vector or control points. In this contribution, we will show how the residuals of the surface approximation can be used to derive the correlation structure of the noise of the observations. We will estimate the correlation parameters using the Whittle maximum likelihood and use comparable simulations and real data to validate our methodology. Using the least-squares adjustment as a “ﬁlter of the geometry” paves the way for the determination of a correlation model for many sensors recording 3D point clouds.


Introduction
The widespread use of three-dimensional (3D) laser-scanning technology offers various possibilities to digitize real-world objects (see, e.g., [1]).For instance, the latest generation of terrestrial laser scanners (TLS) is able to record millions of points during a short time period, allowing various applications from standard deformation monitoring [2] to agricultural uses [3].The noise of the observations is often characterized as being normally distributed; temporal correlations are neglected.This disregard can affect the computation of distances between point clouds recorded at different epochs, as well as test statistics for deformation [4].Furthermore, the deeper study of correlations can provide precious information about the sensor noise with the aim to decrease its level, and innovative applications such as an analysis of the turbulent atmosphere traveled by the electromagnetic signals (see [5] for an application with GPS observations).
The study of correlations is linked with the filtering of noise from the point cloud, which has been the topic of various publications (see, e.g., [6] for a review).Most of the Eng.Proc.2021, 5, 59 2 of 8 strategies cannot ensure that the filtered counterpart-also called residuals-will reflect the statistical property of the observation noise.For example, the wavelet chosen will affect the point cloud mode decomposition and, potentially, the correlation structure of the residuals.This is not the case for the least-squares (LS) approximation, provided that the model linking the observations with some parameters to estimate is optimal.In this contribution, we will make use of this property and filter the geometry of a point cloud with the aim of studying the stochasticity of the underlying sensor noise.The surface fitting will be performed with T-splines [7].This method uses iterative local refinement, as in [8]; it is a powerful and computationally efficient strategy to approximate a point cloud [9], as it is not restricted to predefined forms such as circles or ellipses [10].Both simulated and real data will highlight the feasibility and potential of our methodology to derive the correlation structure of the underlying noise from the LS residuals.We will focus on correlations modeled as a fractional Gaussian noise (fGn, [11]).
The remainder of the paper is as follows: In the second section, we will present the mathematical background of LS and surface fitting.The third section presents the results of simulations and real data.The fourth section concludes this contribution.

Least-Squares
In the following section, we assume that a linear or linearized functional model describes our observations: where I is the n × 1 observation vector, A is the non-stochastic n × u design matrix with full column rank, x the u × 1 parameter vector to be estimated, and v the n × 1 observation noise vector.The error term has zero mean and is normally distributed with E vv T = Σ vv , where Σ vv is the n × n positive definite fully populated variance covariance matrix of the error term.E(.) denotes the mathematical expectation.
The system is overdetermined.A solution can be given by minimizing the sum of the squares of the residuals.In that case, the generalized LS estimator x reads: We call v = y−A x the estimated n × 1 residual vector of the LS adjustment and Σ vv its variance-covariance matrix.We have [12]: In many applications, as, for example, in surface approximation, Σ vv is often replaced with the identity matrix, assuming equal variances for the observations.From (3) we have access to the variance and covariance of the noise observations from the residuals of the approximation, if we assume that the terms contained in Σ vv_est are much smaller than those of Σ vv .The knowledge of Σ vv allows us to avoid an overestimation of the precision of the LS estimator and potentially biased statistical tests.It can be used to describe and quantify the sources of correlations: the cables, the atmosphere, or the processing at the receiver level.Such information provides precious data for stochastic modeling.

General Principle of Surface Approximation
In this paper, we focus on the residuals of surface approximation from scattered data using a mesh optimization approach.In the following, we only provide a short introduction on that topic; interested readers should refer to other publications [7,13].
From now on, we consider the observations to be the Cartesian coordinates of a 3D point cloud, recorded by, for example, a TLS.We assume the point cloud to have been parametrized in advance and the observations to be temporally sorted.We use the T-splines surface approximation to estimate a net of control points-sometimes called vertices.The denser the net, the closest the approximation will be to the noisy measurements.We minimize the squared distance between the noisy points and the mathematical surface and focus on the stochasticity of the residuals.

T-Splines Surface
A T-splines surface is defined as i,Ub (s) as the cubic blending basis functions defined by a recurrence relationship and associated with a knot quintuple Ub = [u i,0 , u i,1 , u i,2 , u i,3 , u i,4 ]; N 3 j,Vb (s) is defined similarly [7].We call x = {d i } the vector matrix, containing the nt control points to be estimated iteratively.The main element of the T-splines surface is the T-mesh, which consists of control points connected by several straight lines.Figure 1 (right) shows an example for the surface under consideration in this paper, depicted in Figure 1 (left).

General Principle of Surface Approximation
In this paper, we focus on the residuals of surface approximation from scattered data using a mesh optimization approach.In the following, we only provide a short introduction on that topic; interested readers should refer to other publications [7,13].
From now on, we consider the observations to be the Cartesian coordinates of a 3D point cloud, recorded by, for example, a TLS.We assume the point cloud to have been parametrized in advance and the observations to be temporally sorted.We use the Tsplines surface approximation to estimate a net of control points-sometimes called vertices.The denser the net, the closest the approximation will be to the noisy measurements.We minimize the squared distance between the noisy points and the mathematical surface and focus on the stochasticity of the residuals.

T-Splines Surface
A T-splines surface is defined as as the cubic blending basis functions defined by a recurrence relationship and associated with a knot quintuple , , , , the vector matrix, containing the nt control points to be estimated iteratively.The main element of the T-splines surface is the T-mesh, which consists of control points connected by several straight lines.Figure 1 (right) shows an example for the surface under consideration in this paper, depicted in Figure 1 (left).The surface approximation is performed iteratively via LS optimization, and usually starts with a basic rectangular and regular mesh structure.The parameters to be estimated are the nt control points from the T-mesh.The optimization problem can be written in matrix form as Ax = B , where A is an ( ) , nt nt matrix that contains the estimation of blending functions at the parameter location: Matrix A depends on the underlying T-mesh via the two local knot vectors . With this notation, we have a unique solution 1 ˆ− x = A B which corresponds to (2) by taking vv Σ as the identity matrix.The residuals of the approximation can be computed for each parametrized observation as . We defined TH as being a user-defined threshold.The surface approximation is performed iteratively via LS optimization, and usually starts with a basic rectangular and regular mesh structure.The parameters to be estimated are the nt control points from the T-mesh.The optimization problem can be written in matrix form as Ax=B, where A is an (nt, nt) matrix that contains the estimation of blending functions at the parameter location: a gi = n ∑ j=1 N 3 g,Ub s j N 3 g,Vb t j , g = 1 . . .nt. Matrix A depends on the underlying T-mesh via the two local knot vectors Ub, Vb.We With this notation, we have a unique solution x = A −1 B which corresponds to (2) by taking Σ vv as the identity matrix.The residuals of the approximation can be computed for each parametrized observation as e i = ŜT (s i , t i ) − l(s i , t i ) 2 , i = 1 . . .n.We defined TH as being a user-defined threshold.If e i > TH, the cells containing the corresponding data points are refined following the algorithm proposed in [14].The local refinement allows an economic surface approximation in terms of the parameters to estimate.A wise choice of TH further reduced the risk of overfitting, i.e., fitting the noise of the observations rather than the underlying surface.We propose to set TH ≈ 2σ Z , with σ Z the variance in the Z-direction, to prevent such unnecessary refinement.The refinement is ended if the errors do not exceed TH or after a given number of iterations.

Residual Analysis
Once the surface approximation is performed, the residuals can be further analyzed.In this paper, we concentrate on the Z-direction only, as the vertical component is known to be noisier than the horizontal ones [10].The observations are sorted temporally so that the residuals can be seen as a time series (1D) rather than as a surface (2D).
From physical considerations, it is justifiable to assume the noise to be an fGn, or a combination of them: the high rate of measurements of most sensors induces a long dependency between the observations.More specifically, we are predisposed to think that flicker and white noise (WN, coming from the electronic component and the processing of the raw observations), or atmospheric noise (from the propagation of the signal) will be present.Although they will be found in different bandwidths, a global correlation parameter can be estimated using the WhiE as proposed by [15].The fGn is fully described by its bounded variance and Hurst exponent H, which is related to the slope of the psd.Following [16], an unbiased likelihood is given by: with Ω the set of discrete Fourier frequencies, f (ω, H) the continuous-time process spectral density and X H,j e −ijω 2 .The Hurst exponent can be estimated as Ĥ = argmax(l W (H)).The slope of the psd β is given by H = β + 1/2.In this paper, we use the WhiE as implemented in MATLAB by [17]; this function has the main advantage of allowing the estimation of the Hurst parameter when the psd saturates at low frequencies (the so-called Matérn covariance model [18]).
The methodology used in this contribution to derive the analysis of the observation noise from the LS residuals is summarized in a flowchart form in Figure 2.
direction, to prevent such unnecessary refinement.The refinement is ended if the errors do not exceed TH or after a given number of iterations.

Residual Analysis
Once the surface approximation is performed, the residuals can be further analyzed.In this paper, we concentrate on the Z-direction only, as the vertical component is known to be noisier than the horizontal ones [10].The observations are sorted temporally so that the residuals can be seen as a time series (1D) rather than as a surface (2D).
From physical considerations, it is justifiable to assume the noise to be an fGn, or a combination of them: the high rate of measurements of most sensors induces a long dependency between the observations.More specifically, we are predisposed to think that flicker and white noise (WN, coming from the electronic component and the processing of the raw observations), or atmospheric noise (from the propagation of the signal) will be present.Although they will be found in different bandwidths, a global correlation parameter can be estimated using the WhiE as proposed by [15].The fGn is fully described by its bounded variance and Hurst exponent H , which is related to the slope of the psd.
Following [16], an unbiased likelihood is given by: ( ) with Ω the set of discrete Fourier frequencies, ( ) .In this paper, we use the WhiE as implemented in MATLAB by [17]; this function has the main advantage of allowing the estimation of the Hurst parameter when the psd saturates at low frequencies (the so-called Matérn covariance model [18]).
The methodology used in this contribution to derive the analysis of the observation noise from the LS residuals is summarized in a flowchart form in Figure 2.

Data Analysis
In this section, we will validate the methodology presented in Section 2 to derive the correlation structure of the measurement noise from the residuals of the LS surface

Data Analysis
In this section, we will validate the methodology presented in Section 2 to derive the correlation structure of the measurement noise from the residuals of the LS surface approximation.Comparable simulated and real data will be used to highlight the potential of the method.

Reference Point Cloud
We simulated 4000 Terrestrial Laser Scanner observations in Cartesian coordinates from a surface corresponding to Figure 1.This reference surface contains different shapes: it mimics a dam in the middle, a mountain in its side, and a wave structure on the top part of it.

Noise Generation
The reference surface is noised by adding:

•
to the X-and Y-components: a Gaussian noise with a standard deviation of 1 × 10 −4 m-generated with the Matlab function randn; • to the Z-component: an fGn or a combination of fGn and WN.We use the MATLAB function ffgn [19].
Two noise vectors are generated: a pure flicker noise (case 1) with H ≈ 0.9, and a combination of 30% white and 70% flicker noise (case 2).The standard deviation of both is set to σ Z = 1 × 10 −3 m, with the aim of approximating the reality of the noise of a TLS sensor [10].

Surface Approximation
The approximation is performed with T-splines using the concept presented in Section 2. We use TH = 2σ Z .6 iterations were performed until the error term did not exceed the threshold.The final T-mesh is shown in Figure 1 (right).

Residual Analysis
The residuals of the approximation are shown for case 1 as an example in Figure 3 (top, red line), together with their power spectral density (psd, Figure 3, bottom).A few outliers caused by the mesh approximation are visible when compared with the reference noise vector (Figure 3, blue line).Fortunately, the shape of the original noise vector is still preserved, which highlights the goodness of the surface approximation.To validate that the correlation structure of the residuals is the same as the one of the measurement noise, we propose to justify empirically that Σ vv_est can be discounted (see [3]).To that end, we firstly assume an equal variance of 1 for the error term in the surface approximation.In Figure 3, we plot the diagonal values and the first 100 lines of Σ vv_est : the diagonal values are more than 5 times smaller than 1; the sparsity of the matrix is evident.The mean value of Σ vv_est is found to be of the order of 1 × 10 −3 << 1: this justifies discounting Σ vv_est as its impact will be small.The psd of the residuals (Figure 3, right bottom) confirms the validity of this assumption; the slopes of the reference noise psd and the residuals are similar (see Table 1 for the corresponding estimates using the WhiE).This result is valid for cases 1 and 2. We mention that the additional WN component acts to decrease the estimated global slope (0.88 instead of 0.9); This is expected, since the WN was not filtered out.We further point out that low frequencies are deleted in the residuals, which can be interpreted intuitively by considering the T-mesh as a high-pass filter (see Figure 1).Fortunately, the WhiE is not affected by the loss of power at low frequencies, as shown in Table 1.The variance estimated from the residuals is slightly underestimated, which may be due to the aforementioned effect.Table 1.Results of the residual analysis for case 1 (no WN) and case 2 (30% WN): slope of the psd and standard deviation (abbreviated as std).

Slope/std [m]
Case 1 Case 2 Original noise 0.9/1 × 10 −3 0.88/1 × 10 The simulated surface generated in Section 3.1 was created with a high-quality 3D Table 1.Results of the residual analysis for case 1 (no WN) and case 2 (30% WN): slope of the psd and standard deviation (abbreviated as std).

Slope/std [m] Case 1 Case 2
Original noise 0.9/1 × 10 The simulated surface generated in Section 3.1 was created with a high-quality 3D printer (see Figure 4a).To that end, the reference surface was firstly approximated with T-splines to serve as a basis for the printing.To generate a 3D-printable plane from a point cloud, this latter was converted into a solid.A surface model was first created from the points using Delaunay triangulation, and a solid was finally generated by uniformly thickening this model.The final export to a 3D printer-compatible file, such as the .stlformat, makes the surface generated from a point cloud 3D-printable.This latter was printed using plastic and stabilized.Grey paint was used in order to avoid strong reflections during the scanning process with a TLS.The additional WN generated by the scanner to generate the surface is considered to be below 1 mm from manufacturer specifications.Consecutively, the printed surface is not exactly the reference one but is superimposed with WN.As the residuals are computed with respect to the generated surface and not the true one, this will have no impact on our conclusions about the correlation structure.

Scanning
The 3D-printed surface was scanned using a phase TLS Zoller+Fröhlich Imager 5016.The scanning configuration was optimal (no angle of incidence, and a surface aligned centrally at the height of the tilting axis of the TLS).The measurements took place indoors at the measuring laboratory of the Geodetic Institute in Hannover (see Figure 4a).
The panel was observed from a total of four standpoints, with distances between 2 and 6 m.No over-radiation occurred.In this contribution, we selected as an example the distance of 2 m, scanned with the angle resolution "high" and the scanning duration "high quality", resulting in a total of 143,490 points.The point cloud was parametrized and approximated with a T-splines surface with TH = 0.001, i.e., we assumed a standard deviation of the noise of about 0.5 mm, following the manufacturer specification.The obtained residuals and their psd are depicted in Figure 4 (r, top) and Figure 5, respectively.

Results
The psd of the whole residuals of the approximation is physically barely interpretable (Figure 5, blue line), i.e., it depicts effects coming from the functional model or T-mesh which acts by adding high frequencies-or identically down-weighting low frequencies, see Section 3.1-yielding a serrated psd.For a better understanding of the correlation structure, we thus selected 5000 epochs of the residuals, which are shown in Figure 4 (right, bottom).Other parts were selected: the same results were comparable and are not presented here for the sake of brevity.We eliminated outliers with the mean absolution deviation method [20] and performed a low-pass Butterworth filter of the first order with a normalized cutoff frequency of 0.05, to further eliminate the WN component following the methodology proposed in [21] (see Figures 4 and 5, red line).As expected from [22], we can identify different noises present in different bandwidths: a strong WN component at high frequencies, an fGn with a slope close to −8/3, followed by an fGn with a slope close to −2/3 and a saturation at low frequencies.The lack of power at low frequencies is to be interpreted following the results of the simulation, and comes from the T-splines surface approximation.The global slope estimated with the WhiE for the non-filtered residuals was found to be −0.8 due to the WN component and close to −2.8 for the filtered residuals.Atmospheric turbulence affects optical signals traveling through a random medium close to the earth's surface and acts to correlate the measurements.This slope is close to the one that is expected from the turbulence theory [21].The variance of the non-filtered residuals after outlier elimination was 5.8 × 10 −4 , and is close to the expected value from the manufacturer specification.
surface approximation.The global slope estimated with the WhiE for the non-filtered residuals was found to be −0.8 due to the WN component and close to −2.8 for the filtered residuals.Atmospheric turbulence affects optical signals traveling through a random medium close to the earth's surface and acts to correlate the measurements.This slope is close to the one that is expected from the turbulence theory [21].The variance of the nonfiltered residuals after outlier elimination was 5.8 × 10 −4 , and is close to the expected value from the manufacturer specification.

Conclusions
With the continuous increase of the data rate, the measurement noise of sensors is more likely to become strongly temporally correlated.In this contribution, we have demonstrated how the residuals of the LS approximation can be used to analyze the temporal correlation structure of measurement noise from a TLS.As an example, we surface approximation.The global slope estimated with the WhiE for the non-filtered residuals was found to be −0.8 due to the WN component and close to −2.8 for the filtered residuals.Atmospheric turbulence affects optical signals traveling through a random medium close to the earth's surface and acts to correlate the measurements.This slope is close to the one that is expected from the turbulence theory [21].The variance of the nonfiltered residuals after outlier elimination was 5.8 × 10 −4 , and is close to the expected value from the manufacturer specification.

Conclusions
With the continuous increase of the data rate, the measurement noise of sensors is more likely to become strongly temporally correlated.In this contribution, we have demonstrated how the residuals of the LS approximation can be used to analyze the temporal correlation structure of measurement noise from a TLS.As an example, we

Conclusions
With the continuous increase of the data rate, the measurement noise of sensors is more likely to become strongly temporally correlated.In this contribution, we have demonstrated how the residuals of the LS approximation can be used to analyze the temporal correlation structure of measurement noise from a TLS.As an example, we performed the extraction of the geometry of a TLS point cloud using an LS approximation with a T-splines surface.This latter has the main advantage of being computationally efficient; it provides an optimal functional model for the further analysis of the correlation structure of the residuals of the approximation.Simulations were set up to validate the methodology with a noised reference surface.We used the unbiased WhiE to estimate the slope parameter of the psd.The same surface was printed in 3D and scanned with a TLS.The residuals could be successfully interpreted based on the results of the simulations.The psd of the real data analysis showed noises present in different bandwidths.Besides a WN component, some of them were identified as coming from the propagation of optical signals through a turbulent medium.This analysis is a validation of the potential of the LS surface approximation to quantify and analyze the correlation structure of sensor noise.Further investigations will focus on its modelization, as well as the impact of functional mismodeling on the spectral decomposition of the residuals.

Figure 1 .
Figure 1.Reference surface under consideration and the corresponding T-mesh after the 6th iteration; threshold (TH) 0.002.

Figure 1 .
Figure 1.Reference surface under consideration and the corresponding T-mesh after the 6th iteration; threshold (TH) 0.002.

.
The slope of the psd β is given by 1 2H β = +

Figure 2 .
Figure 2. Flowchart explaining the methodology to extract the geometry of a point cloud to analyze the correlation structure of the residuals.

Figure 2 .
Figure 2. Flowchart explaining the methodology to extract the geometry of a point cloud to analyze the correlation structure of the residuals.

Figure 4 .
Figure 4. (a) 3D print of the reference surface used for the simulations.(b) The residuals of the approximation.

Figure 5 .
Figure 5.The psd of the whole residuals (blue line), of the selection section (yellow line), and of the filtered section (red line).

Figure 4 .
Figure 4. (a) 3D print of the reference surface used for the simulations.(b) The residuals of the approximation.

Figure 4 .
Figure 4. (a) 3D print of the reference surface used for the simulations.(b) The residuals of the approximation.

Figure 5 .
Figure 5.The psd of the whole residuals (blue line), of the selection section (yellow line), and of the filtered section (red line).

Figure 5 .
Figure 5.The psd of the whole residuals (blue line), of the selection section (yellow line), and of the filtered section (red line).