A spatio-temporal deformation model for laser scanning point clouds

The establishment of the terrestrial laser scanner changed the analysis strategies in engineering geodesy from point-wise approaches to areal ones. During recent years, a multitude of developments regarding a laser scanner-based geometric state description were made. However, the areal deformation analysis still represents a challenge. In this paper, a spatio-temporal deformation model is developed, combining the estimation of B-spline surfaces with the stochastic modelling of deformations. The approach’s main idea is to model the acquired measuring object by means of three parts, similar to a least squares collocation: a deterministic trend, representing the undistorted object, a stochastic signal, describing a locally homogeneous deformation process, and the measuring noise, accounting for uncertainties caused by the measuring process. Due to the stochastic modelling of the deformations in the form of distance-depending variograms, the challenge of defining identical points within two measuring epochs is overcome. Based on the geodetic datum defined by the initial trend surface, a point-to-surface- and a point-to-point-comparison of the acquired data sets is possible, resulting in interpretable and meaningful deformation metrics. Furthermore, following the basic ideas of a least squares collocation, the deformation model allows a time-related space-continuous description as well as a space- and time-continuous prediction of the deformation. The developed approach is validated using simulated data sets, and the respective results are analysed and compared with respect to nominal surfaces.


Introduction
Deformation analysis has always been part of a large range of application fields: the monitoring of an object's change over time is of high interest in gas and oil production, civil and mechanical engineering, hydrology or environmental sciences (Velsink 2015). Classical geodetic approaches like levelling, GNSS or total station measurements allow a point-based determination of deformations by repeatedly measuring representative points of an object. The evaluation of the resulting coordinate differences represents the object's deformation (Wunderlich et al. 2016). Although there exists a variety of sophisticated analysis strategies for such pointbased approaches, they also entail some drawbacks: the B Corinna Harmening corinna.harmening@tuwien.ac.at Hans Neuner hans.neuner@tuwien.ac.at 1 Department of Geodesy and Geoinformation, TU Wien, Vienna, Austria appropriate choice of representative points requires some prior knowledge of the expected deformations and the needed signalization of those points makes such approaches only applicable for accessible measuring objects. Furthermore, the resulting information is always sparse. Finally, point-based approaches can be very expensive as well as time-and labourintensive, especially for large measuring objects (Paffenholz et al. 2017;Shamshiri et al. 2014;Li et al. 2015).
With the development of the terrestrial laser scanner (TLS), a measuring instrument which allows a fast, efficient and contactless data acquisition even of inaccessible measuring objects moved into focus of engineering geodesy. The acquired data are of high resolution giving a quasi-continuous description of the measuring object (Paffenholz et al. 2017). All in all, TLS provides the best conditions for an areal deformation analysis, overcoming the drawbacks of classical approaches.
However, when performing a laser scanner-based deformation analysis, new challenges occur (cf. Wunderlich et al. 2016;Mukupa et al. 2016;Holst and Kuhlmann 2016). The impossibility of reproducing measured points in different epochs and the resulting question of how to compare two point clouds in an appropriate way have to be mentioned in this context. Furthermore, the lack of an appropriate error model for laser scanner may result in a superimposition of actual deformations and systematics caused by the measuring process. These challenges are the reason why-although laser scanning has established as a measuring techniquethere is still a lack of appropriate analysis strategies and of significance tests regarding an areal deformation analysis (Wunderlich et al. 2016).
A first step for meeting these challenges was made by defining three possibilities to compare different point clouds (Mukupa et al. 2016): the point-to-point-based comparison, the point-to-surface-based comparison as well as the surfaceto-surface-based comparison. Approaches belonging to all of these three classes can be found in the literature: in Little (2006), the points of two different point clouds are directly compared in order to determine the deformations of a slope. However, as it is impossible to measure the same point twice (Zamecnikova and Neuner 2017), advanced approaches to define point correspondences are needed. In Paffenholz et al. (2017), the Multiscale Model to Model Cloud Comparison (M3C2) algorithm (Lague et al. 2013) is used to determine deformations of a historic bridge. A stochastic approach to define point correspondences can be found in Wujanz (2018).
When using the second strategy, one of the point clouds is approximated by a surface (either a mesh or a fitted analytical function) and the deformation is reflected by the distances of the second point cloud to this surface (Mukupa et al. 2016). Examples can be found in Ioannidis and Valani (2006),  and Erdélyi et al. (2017). In these three publications, the point clouds are compared to analytical functions like non-uniform rational B-splines (NURBS), paraboloids or planes, whereas in Serantoni and Wieser (2016) a mesh is used as a reference surface.
One common way to use the surface-to-surface-based approach is the comparison of the estimated parameters characterizing the respective analytical functions over time. Examples can be found in Vezočnik et al. (2009) or Lindenbergh and Pfeifer (2005).
In this contribution, an approach to an areal deformation analysis is developed which allows a point-to-surface-based comparison of laser scanning point clouds with respect to an initial and undistorted reference surface as well as a pointto-point-based comparison between two distorted states of the measuring object. Furthermore, the approach allows a prediction of the deformation within a measured epoch as well as into a non-measured epoch. Thus, unlike the deformation models introduced so far, the developed model allows a time-continuous areal description of deformations. A first step towards a time-continuous areal deformation analysis was made in Kutterer et al. (2009) by applying time series analysis to laser scanning profiles.
The basis of the approach developed in this paper is formed by an initial approximation of the point cloud of the first measuring epoch by means of a B-spline surface. This surface serves as a reference surface which is assumed to represent the undistorted measuring object. The deformations with respect to this reference surface are modelled stochastically similar to the signal in a least squares collocation.
The paper is structured as follows: Sect. 2 provides the methodological basics for the development of the presented approach. In Sect. 3 the data sets, on which the developed approach is applied, are introduced. As the simulation process motivates the developed analysis approach, Sect. 3 provides the basis for Sect. 4, which is the main part of this contribution. It deals with the development of a space-and time-continuous areal deformation model and the application to simulated data sets. The results of different data sets are analysed, evaluated and compared in Sect. 5. In Sect. 6, the prediction of the deformations within a measured epoch as well as into a non-measured epoch is presented and applied to simulated data sets. Finally, a conclusion is drawn and the limitations of the developed approach as well as future investigations are discussed in Sect. 7.

Estimation of B-spline surfaces
A B-spline surface of degree p and q is defined by: (1) According to Eq. (1), a surface pointŜ (u, v) is expressed as the weighted average of the (n P + 1) × (m P + 1) control points P i j (Piegl and Tiller 1995, p. 100 ff.). The corresponding weights are defined by the functional values of the B-spline basis functions N i, p (u) and N j,q (v), which can be recursively computed (Cox 1972;de Boor 1972). Two knot vectors, one in direction of the surface parameter u (U = [u 0 , . . . , u r ]) and one in direction of the surface parameter v (V = [v 0 , . . . , v s ]), split the B-spline's domain into knot spans. As a consequence, the shifting of one control point changes the surface only locally. Usually, only the location of the control points is estimated in a linear Gauß-Markov model when estimating a best-fitting B-spline surface. The choice of the optimal number of control points to be estimated (n P + 1) and (m P + 1), respectively, is a model selection problem and can be solved by classical model selection criteria or by structural risk minimization Neuner 2016a, 2017). In order to obtain a linear relationship between the 3n l observations l k = S(u k , v k ) = [x k , y k , z k ] T with (k = 1, . . . n l ) and the unknown control points P i j , the B-spline's knots as well as its degrees are usually specified a priori. In this study, the use of cubic B-splines with p = q = 3 is applied as it covers the geometric continuity properties of a large amount of real-world man-made structures. Methods for determining appropriate knot vectors can be found in Schmitt and Neuner (2015) or Bureick et al. (2016). Furthermore, convenient surface parameters u and v, locating the observations on the surface to be estimated, have to be a-priori allocated to the observations (Harmening and Neuner 2015).

Least squares collocation
Least squares collocation (LSC) was originally developed for purposes of physical geodesy (cf. Borre and Krarup 2006;Moritz 1989), but it is also applied in other fields of geodesy and geostatistics nowadays (Mysen 2014;. The extension of the functional model of the classical least squares adjustment leads to the functional model of LSC (Heunecke et al. 2013, p. 204 ff.): The observed phenomenon is modelled as the sum of a deterministic trend Ax, characterized by the vector of unknowns x, and the component Rs. The stochastic signal s carries information about the phenomenon in the form of stochastic relationships and is assumed to be normally distributed with expectation 0 and covariance matrix Σ ss : The discrepancy of the measurements l and the phenomenon is described by the stochastic measurement noise , which is also assumed to be normally distributed: Correlations between signal and noise are excluded: The aim of LSC is threefold (Höpcke 1980, p. 210 f.): in an adjustment, the parameter vector x is estimated with respect to an optimality criterion. The filtering reduces the noise in the measured points. Finally, the prediction aims to determine trend and signal in unobserved locations. In order to combine those three tasks, Eq.
(2) is complemented by the signal values to be predicted s .
Equation (6) corresponds to a Gauß-Helmert model resulting in the following formulas for the filtering and the prediction: with the Lagrange multiplierŝ and the estimated trend The respective cofactor matrices arise from variance covariance propagation to For a detailed derivation of these formulas, we refer to Heunecke et al. (2013, p. 205 ff.).

Spatio-temporal stochastic processes
The observations of a LSC are interpreted to be a realization of a stochastic process Z . Such a process can either be a pure function of time Z (t), a function of the phenomenon's location in space Z (x) with the three-dimensional space vector x or a spatio-temporal function Z (x, t).
Time series analysis provides computational tools for analysing data dependent on time, whereas geostatistics deals with the analysis of stochastic processes dependent on position and their extension by the time domain (Cressie and Wikle 2015;Matheron 1963). Spatio-temporal kriging is one of the geostatistical main tools. It is used in hydrological applications to model rainfall (Bargaoui and Chebbi 2009) or to combine data from different altimetry missions (Boergens et al. 2017), in environmental applications to forecast irradiance (Aryaputera et al. 2015) as well as in soil science in order to predict soil water content (Snepvangers et al. 2003). A geodetic application of kriging to determine a regional ionospheric model can be found in Abdelazeem et al. (2018).

Properties of stochastic processes
Classical methods to estimate statistical moments of the observed process require several realizations of the stochastic process. However, in practice only one single realization is available, making additional assumptions necessary (Schlittgen and Streitberg 2013, p. 100 f.): a stochastic process is defined to be stationary if its statistical moments are constant over time and if its joint statistical moments depend only on the time lag τ between two observations. Special cases among stationary processes are ergodic processes: different realizations of an ergodic process result in identical mean values and autocorrelation functions (Bendat and Piersol 2010, p. 12).
The transfer of these definitions to the space domain leads to the definition of homogeneity (location-invariant statistical moments and distant-depending joint statistical moments). Furthermore, a spatial stochastic process is called isotropic if the joint statistical moments are direction-independent.
In some literature, stochastic processes with locationinvariant statistical moments are also defined to be stationary. However, this paper follows the definition of Bendat and Piersol (2010), allowing a clear distinction between space and time domain.

Modelling dependencies of stationary/homogeneous stochastic processes
Dependencies of a detrended (z = 0) and equidistantly acquired time series z(t j ) with j = 1, . . . , n l are analysed by means of autocovariance functions. Assuming stationarity, this function can be biasedly estimated (Koch et al. 2010): with τ = t j+κ − t j . In order to ensure a stable estimation, the parameter κ is chosen to be κ = 0, 1, . . . , n l /10.
A covariogram C s (d) is the geostatistical analogon to autocovariance functions with d being a spatial distance (Matheron 1963). However, the majority of the geostatistical literature put emphasis on the calculation and interpretation of variograms (e.g. Aryaputera et al. 2015;Snepvangers et al. 2003;Tapoglou et al. 2014): averaging the squared differences over all point pairs whose distances d i j = ||x i − x j || are contained in the interval N k (k = 1, . . . , n N ) (Cressie and Wikle 2015, p. 131). The sub-division of the range of separation distances d i j into n N consecutive intervals N k is necessary as equidistant data cannot be taken for granted in the spatial domain. As a consequence, the variogram is a function of the mean separation distance d k of all point pairs belonging to interval N k . Contrary to temporal data, spatial data may be anisotropic (Matheron 1963). In these cases, the above-described grouping of the point pairs has to be realized according to the absolute distance d i j , but also according to the orientation θ i j of the separation vector d i j = x i − x j , resulting in directional variograms.
As the estimation of variograms is more stable than the estimation of covariograms (cf. Smith 2016, p. 4-29 ff.), solely variograms are estimated to describe spatial relationships in the following. Afterwards, the estimated variograms are transformed into covariograms: with σ 2 being the variance of the process. Standardizing the autocovariance (eitherĈ t (τ ) orĈ s (d k )) by the data sets' variance σ 2 =Ĉ t/s (0) gives the estimator of the autocorrelation function/correlogram: Equation (13) can be generalized to the crosscovariance function, describing the similarity of two different detrended time series (Heunecke et al. 2013, p. 348): κ = −n l /10, . . . , n l /10. Standardizing Eq. (17) results in the crosscorrelation function: The respective transfer to the space domain is straightforward and results in cross(co)variograms and crosscorrelograms: When analysing the estimated empirical covariances, it has to be taken into account that, usually, the observed phenomenon The first one is represented by the measuring noise and, in this study, is initially assumed to be a white noise process (top left picture in Fig. 1). The second process causes the correlated random variables of the signal s. A typical covariogram of such a process can be seen in the bottom left picture in Fig. 1. The estimation of empirical correlations of such a combined process results in a covariogram which is the sum of both individual covariograms (blue dots in the right picture in Fig. 1) (Smith 2016, p. 4-6 ff.). One possibility to account for this superimposition is the estimation of an analytical covariance function based on the empirical values with the exception of the first one. The resulting analytical function directly separates the two processes as can be seen in Fig. 1. Alternatively, the Dirac function can be used to model the superimposition of the two models (cf. Koch et al. 2010).

Locally stationary stochastic processes
In reality, stationary/homogeneous processes occur quite rarely and the assumptions of stationarity and homogeneity are only approximations, allowing for the use of a wide range of computation tools (Bendat and Piersol 2010, p. 344). Among the variety of non-stationary processes, there exists the special case of locally stationary processes (Silverman 1957): considering two zero-mean non-stationary random processes Z 1 (t) and Z 2 (t), the autocorrelation functions as well as the crosscorrelation function are time-dependent (Bendat and Piersol 2010, p. 358 ff.). The transformations lead to the crosscorrelation function: If this function can be split up into a product the process is said to be locally stationary. In equation (25), ρ μ (t) is a slowly varying scale factor, whereas ρ Δ (τ ) is a stationary correlation function. The class of locally stationary processes is of great importance as these processes are usually a good approximation to non-stationary processes.

Confidence bands for spatial dependent data
Usually, bootstrap methods are used to estimate confidence bands for correlograms. However, classical bootstrap methods are based on strong assumptions such as independently and identically distributed data and, therefore, are not suitable to estimate confidence bands in case of highly correlated data. An overview of alternative methods can be found in Clark and Allingham (2011). In this article the parametric spatial bootstrap introduced in Tang et al. (2006) is used: the basis for this bootstrapping approach is formed by an uncorrelated data set * ∼ N (0, 1). In each of B bootstrapping steps, a bootstrap sample * B is drawn with replacement. The initial independent data are correlated according to the estimated correlation structure of the actual data set, represented by the estimated covariance matrixΣ using its Cholesky decompo-sitionΣ =L ·L T : Each of these B correlated data sets is afterwards used to estimate a correlogram. The resulting variance over these estimated replicates is used to establish confidence bands. Even a misspecified covariogram model for setting upΣ leads to appropriate confidence bands (Tang et al. 2006).

Data simulation
The deformation model developed in Sect. 4 is applied to two types of data sets: the stochastically simulated data sets help to understand the general procedure of the developed approach, as it is a "backwards" application of the simulation process, whereas the functionally simulated data sets are used to demonstrate the applicability of a stochastic approach to deterministically deformed data sets.  The basis for both simulation processes is provided by a B-spline surface with 63 control points ((n p + 1) = 9, (m p + 1) = 7) and with dimensions of approximately 40 cm x 40 cm x 18 cm (cf. Fig. 2).

Stochastically simulated data sets
The basic idea of the developed approach is that the observations consist of three parts as defined in Eq. (2): suring object and, thus, is identical in all measuring epochs. In this study, three measuring epochs are simulated and the B-spline surface in Fig. 2 serves as the trend surface. In order to simulate the TLS measuring process, this surface is sampled with a spatial resolution of approximately 6 mm. -The signal component Rs is of particular importance as it captures the deformation. Because of the signal's great significance, it will be treated in detail subsequent to this listing.
-The measuring noise presents the uncertainty which is caused by the measuring process. Due to missing or incomplete realistic models to describe the stochastic behaviour of a laser scanner's measuring process, white noise with a standard deviation of σ = 1 mm is used in the following. This simplifying assumption is only a first step in the development of a spatio-temporal deformation analysis and will be adapted in future investigations.
The main idea how to provoke a deformation by means of the signal component Rs and, thus, purely caused by stochastic relationships between observed points can be seen in Figs. 3, 4 and 5. The basis of this example is formed by a one-dimensional time series of normally distributed and uncorrelated random variables Z ∼ N (0, I) (blue dots in all three subplots of Fig. 3). Three positive definite exponential correlation functions of the type which differ in their correlation length due to the choice of the parameter b i (see Fig. 4), are used to set up covariance matrices Σ i with i = 1, 2, 3. Computing the Cholesky decompositions Σ i = G T i G i , the random variable Z can be transformed to a normally distributed and correlated random variable (Koch 1997, p. 167): The resulting three correlated time series can also be seen in Fig. 3 (red dots). Comparing those three time series, it becomes apparent that the stronger the correlations, the less pronounced the stochastic behaviour within the time series and the slower the fluctuation around the expectation value of zero. In case a section of the time series, whose length is considerably smaller than the correlation length, is observed, strong correlations can lead to a situation in which the correlated time series seems to be shifted away from the respective expectation value. This situation occurs in the lower picture in Fig. 3 when observing the time interval between 200 and 350 s. While in this simulation study the expectation value of zero is guaranteed by Eq. (28), another indicator for an expectation value of zero is the (possibly very slow) convergence of the square root of the correlation function towards the random variable's expectation value μ. This convergence has not necessarily to happen within the observed interval, but only for τ → ∞ (Bendat and Piersol 2010, p. 20): In order to stochastically convert this pure translation into a more general shape of a deformation, the principle of locally stationary stochastic processes (see Sect. 2.3.3) is utilized, which are characterized by a stationary correlation function and a slowly varying variance (cf. Eq. 25). Thus, a timedependent variance is assigned to the original time series values, while the respective stationary correlation structures given by Fig. 4 are maintained over the entire time series.
The results of this modified simulation process can be seen in Fig. 5, depicting only the 150 s-time interval between 200 and 350 s. The variance level was chosen to be a tenth of the original variance level for the main part of the presented time interval. In these parts, the red points in Fig. 5 scatter minimally around the expectation value of zero. In the interval 240-260 s the variance level is linearly increased with the maximum value at 250 s and afterwards decreased again to the minimal variance level. This section is clearly noticeable in all three time series by the increased variation range of the resulting values. However, only in case of a very slowly decaying correlation function an effect occurs which can be interpreted to be a typical deformation. The example in Fig. 5 (bottom) can be interpreted to represent the bending of a beam, observed at one single point over time. After 240 s an increasing load is starting to act on the beam, resulting in a deflection of the beam which is represented by the red curve. When the load decreases after 250 s, the beam slowly returns to the initial state. The type of deformation (elastic, plastic, linear, periodic, etc.) is to a large extent controlled by means of the variance level's behaviour.
In order to stochastically simulate areal deformations, this principle is extended to the spatio-temporal domain, interpreting the entirety of observations in all measuring epochs as a single realization of a stochastic process: with the exception of the point cloud of the first measuring epoch, the sampled trend surface is divided into a distorted and a non-distorted part for each measuring epoch. The former is subsequently subdivided into areas with constant variance, whereas the latter is no longer considered in the simulation of the stochastic deformation. For reasons of simplicity and due to the definition of the coordinate system, the deformation is realized solely in direction of the z-coordinate in this study. Thus, autocorrelograms for the z-coordinate of the second and third measuring epochs as well as crosscorrelograms between the z-coordinates of the second and third measuring epochs are set up. Based on these correlograms and the a-priori defined variance levels, a normally distributed and uncorrelated signal can be transformed into a normally distributed and correlated signal describing a spatio-temporal deformation.
Adding this deformation to the noisy trend surface leads to the data set which can be seen in Fig. 6. The blue point cloud presents the undistorted measuring object, which is superimposed by the measuring noise. A clearly visible deformation which is simulated solely by means of stochastic relationships occurs in the other two measuring epochs (red and yellow point clouds).

Functionally simulated data sets
In order to functionally simulate deformations, the shape of the B-spline surface is directly changed by linearly moving the control point P 4,5 (encircled in red in Fig. 2) in the direction of the z-coordinate, following the equation of motion The resulting surfaces at defined points of times in arbitrary time units (t.u.) t 1 = 0, t 2 = 1, t 3 = 1.5, t 4 = 1.75, t 5 = 2 represent various states of the deformed measuring object, with the surface at t 1 = 0 representing the undistorted object. In Fig. 7, two of these states can be seen: due to the local support of B-spline functions, the movement of a single control point affects the surface only locally. We assume that four of the states mentioned above are acquired by means of a terrestrial laser scanner. The scanning is realized by sampling the respective surfaces and by subsequently adding uncorrelated noise with σ = 1 mm to the sampled point clouds (PC). The sampling resolution was chosen to be approximately 6 mm, resulting in 2500 points for each measuring epoch. The non-measured surface at t 4 = 1.75 is only used as a verification of the prediction step.
In the following, particular emphasis is put on the stochastically simulated data set PC s , whereas the functionally simulated data set PC f is used to demonstrate the applicability of a stochastic modelling to functionally simulated deformations. Thus, unless otherwise stated, the discussions  on the empirical issues always refer to data set PC s . Both data sets are summarized in Table 1. mentation issues and by the results obtained for the simulated data sets, which were introduced in the previous section.

Definition of the approach's framework
The approach developed in the following is placed into the last step of the process chain of a laser scanner-based deformation monitoring (cf. Wujanz 2018; Holst and Kuhlmann 2016) and, therefore, deals solely with the quantification of deformations. All foregoing steps are assumed to be solved in a satisfying manner so that no systematics due to different scanning conditions, the registration process or an incorrect error model occur. The point clouds of the different measuring epochs are assumed to be available in the same geodetic datum, making a comparison possible. When using the term "deformation", it has to be distinguished between rigid body movements (the object underlies solely rotations and translations) and distortions (the object changes in shape) (Heunecke et al. 2013, p. 92 ff.). In this contribution, the focus lies on distortions. Thus, unless otherwise stated, the term "deformations" always refers to a change of shape in this study.

Basic ideas of the developed approach
The basic ideas of the developed approach can be summarized as follows: -The major assumption of this approach is that the deterministic trend describes the initial geometry of the object, captured in the first measuring epoch. This trend solely refers to the space domain and is the same for each epoch. Consequently, the flexibility and approximation power of B-spline surfaces can be used to model the object's initial geometry. -The deformation process is interpreted to be a Gaussian multivariate spatio-temporal stochastic process (Genton and Kleiber 2015) with: Hence, the deformation process is exclusively modelled by means of a stochastic signal similar to a least squares collocation. In principle, the signal's definition in the space of the surface parameters u and v is imaginable. However, as the deformation itself has to be expressed in the Euclidean space to be meaningful and interpretable, in this study, the signal is defined in the Euclidean space.
-It has to be noted that a functional modelling of the deformations by approximating the point clouds with different B-spline surfaces and by expressing the deformations in terms of the B-spline surfaces' parameter groups' changes is possible in principle. However, the proposed approach of a stochastic modelling of the deformation offers two advantages: on the one hand, the challenge of a consistent surface parametrization, which is necessary to compare different B-spline surfaces (cf. Harmening and Neuner 2017), is circumvented. On the other hand, only in a unified B-spline-based framework for handling rigid body movements and distortions, the former can be a-priori eliminated according to Harmening and Neuner (2016b). This elimination requires trend surfaces with fixed number of control points in each measuring epoch.
In this interpretation of the deformation process, the B-spline surfaces with fixed number of control points correspond to the reference potential of the gravitational field in the classical interpretation of a least squares collocation established in physical geodesy (cf. Moritz 1989, p. 99;Koch 1997, p. 241).
Remaining systematics in the point clouds when having subtracted the trend surfaces are deliberately accepted in order to exploit these advantages. Although the deformation process is per se a deterministic phenomenon, a treatment within a stochastic framework can be justified: in systems theory, a distinction is made between causal models and descriptive ones. Up to now, the developments concerning the former models, which are based on a physical model of the processes causing the deformations, are still far from being applicable to a space-continuous deformation analysis. The latter ones usually require point correspondences in order to functionally describe the deformation. An alternative to the functional description is the treatment within a stochastic framework which does not require explicit point correspondences and can be seen as a good approximation of the deterministic phenomenon. Following these considerations, the deformation process is interpreted to be mean-homogeneous with This assumption is valid in the elastic deformation domain, where most deformation activities take place. It implies that in expectation no deformations occur. This complies with the well-established null hypothesis of the congruence model, stating that no deformations occur (cf. Heunecke et al. 2013, p. 252). Furthermore, it covers the lack of knowledge regarding the direction of the deformation. In case of the stochastically simulated data sets, this assumption corresponds to a mean value of s = 0 with respect to the entirety of realizations, whereas one single realization expresses the deformation (cf. Moritz 1989, p. 99 ff. for a similar interpretation of the signal's expectation value). Consequently, the deformations, expressed as residuals with respect to the estimated trend surface, are characterized by the (co)variances of the spatio-temporal stochastic process. As the magnitude of these residuals may strongly change over the whole distorted area and over the entire measuring period, the process is furthermore interpreted as variance-inhomogeneous. However, when excluding discontinuous deformations, a locally homogeneous model is a reasonable choice. Thus, in this model the correlation structure is identical over the entire surface, whereas the variance is a slowly varying function of the location, resulting in the spatio-temporal extension of Eq. (25), following (Genton and Kleiber 2015): with:

Derivation of a spatio-temporal deformation model
Starting point of the deformation model is the functional model of a least squares collocation [Eq.
(2)], which is extended by k epochs. Choosing the matrix R to be the identity matrix leads to: . . .
According to the idea of a stochastically modelled deformation, the trend can be estimated once and is afterwards excluded from the functional model: ⎡ . . .
Hence, the actual observations of the newly developed approach are the measurements' residuals with respect to the trend: ⎡ . . .
The extension by the predicted signal gives: . . .
which can be written more compactly: Although the starting point of this approach is equal to that of a LSC, the strict separation of trend (equals the undistorted object) and signal (equals the deformation) results in an approach which is only similar to a LSC. Nevertheless, the common indications like "trend", "signal", "filtering" and "prediction" used in LSC are maintained in the following.
The formulas above indicate that the application of the developed deformation model can be divided into three main parts: the modelling of the trend, resulting in the estimated parameter vectorx (1) , the modelling of the signal, which is fully described by means of its covariance matrix Σ ss as well as the modelling of the noise which is represented by the covariance matrix Σ .
The determination of those parts is a backward application of the simulation process described in Sect. 3.1. It is summarized by means of the flowchart in Fig. 8 and is described and illustrated in detail in the following sections.

Modelling of the trend
Based on the ideas presented in Sect. 4.2 and the methodological basics given in Sect. 2.1, the point cloud of the first measuring epoch is used to estimate the B-spline surface's control points P (1) = x (1) in a linear Gauß-Markov model. For the sake of simplicity, the remaining B-spline parameter groups (degrees of the B-spline basis functions, number of control points, surface parameters, knot vectors) are chosen to be the nominal ones, which are known due to the simulation process. Based on the surface parameters u i and v i (i = 1, . . . , n l ), which locate the observations The control points can be estimated: Afterwards, these control points are used to estimate the observations describing the undistorted object in each measuring epoch: For this step, it is assumed that the surface parameters (u, v) endure during all measuring epochs.
The residuals of the trend estimation, introduced in equation (38), are then composed of the measurement noise and the object's distortion. In Fig. 9, the estimated point cloudPC s can be seen. The colouring of the estimated observations corresponds to As can be seen, the object is divided into a distorted part in the middle of the object with dimensions of approximately 10 cm x 15 cm and with solely negative residuals as well as a non-distorted part with randomly scattering residuals.
Although the simulation of the deformation is realized by provoking the signal solely in the z-direction, the residuals of the B-spline estimation also reveal a distortion in x-and y-direction.

Detection of distorted regions
As the signal solely occurs in distorted regions, its modelling requires a distinction between distorted and non-distorted regions. Especially the inspection of the residuals of the functionally simulated data sets reveals that the extents of the distorted regions are not identical for the three coordinate directions. This is caused by the coordinate-wise estimation of the B-spline surface. Thus, for reasons of consistency, the detection of the distorted regions is also performed coordinate-wise. For the sake of simplicity, in the following the equations are specified solely for the z-coordinate, which is indicated by the index z. Naturally, the respective computations are performed for x and y, too.
The detection is based on the knowledge that there do not exist any distortions in the first measuring epoch. Consequently, the point cloud of this epoch contains information about the magnitude of the measurement noise. In the following, the variance of the first measuring epoch is used as a threshold to detect distortions: Assuming that regions, in which the discrepancy between observations and estimated trend is by a certain amount larger than the measurement noise, are distorted, each coordinate whose residual fulfils is marked to belong to the distorted region. The choice of this heuristic threshold has a serious impact on the filtering results. It is chosen according to equation (49) in order to keep the type II error low. A detailed justification can be found in Sect. 4.6.5. Due to the choice of the relatively small threshold, there is a relatively high probability of type I errors. Such points are automatically detected and allocated to the non-distorted area in a post-processing step.
This threshold consideration results in a rough distinction between distorted and non-distorted regions for each coordinate direction, which can be seen in Fig. 10 for the z-direction of point cloud PC s . It has to be noted that the method is only applicable if the measuring process does not affect the measuring noise that is to say that a changed measuring configuration does not lead to an increase of .
An alternative way to determine the noise level by setting up synthetic covariance matrices can be found in Kauker et al. .

Modelling of the signal
The complete modelling of the object's distorted parts requires an estimation of the signal, which includes the determination of the locally homogeneous variances and the modelling of the homogeneous correlation structure. This is a multi-stage procedure, which is introduced in the following subsections.

Creating locally homogeneous areas
In order to account for the signal's local homogeneity, in a first step the distorted part is subdivided into areas in which the corresponding signal can be considered homogeneous. This subdivision is achieved by a k-means clustering (Lloyd 1982). The only parameter of the algorithm is the number of clusters n c and has to be chosen strategically in dependence on the size of the distorted area. In Fig. 11, the result of the clustering can be seen for the z-coordinate of the point cloud PC (3) s , when the number of clusters is chosen to be n c = 12. The choice of the number of clusters will be justified in Sect. 4.6.5. In these areas, the first factor in Eq. (33), the variance of the mean-homogeneous process, is assumed to be constant. In compliance with Fig. 5, the discrepancies with respect to the trend are solely determined by the locally homogeneous standard deviations which can be computed for each cluster c j in accordance with the 3-sigma rule:

Establishing global homogeneity
In order to determine the homogeneous correlation structure, the residuals of epoch (i) with respect to the trend are normalized for each cluster c j by means of the linear transformation: Due to the homogeneous magnitude of these residuals, they are used in the following to compute empirical correlograms according to Sect. 2.3.2.

Estimation of empirical correlograms
After having standardized the residuals according to equation (51), a variety of empirical correlograms can be estimated: The . . , k) of each measuring epoch reflect the stochastic relationships between the same coordinate types within this epoch. In order to maintain consistency with the Bspline estimation, these computations are also performed coordinate-wise according to Eqs. (14)-(16). In Fig. 12 (top), these empirical correlograms can be seen for the point cloud PC (3) s . As could be expected from the simulation process in Sect. 3.1 (cf. Fig. 4), the autocorrelograms show a very slow decay from full correlation and especially the autocorrelograms in the direction of y-and z-coordinate are far from reaching the expectation value of zero within the observed area. Comparing the curves of ρ x x , ρ yy and ρ zz , the different lengths of the curves attract attention. This behaviour is Due to the coordinate-wise approach, spatial crosscor- (21)), reflecting the stochastic dependencies between the three coordinate directions within each measuring epoch. Figure 12 (bottom) shows the respective results for point cloud PC (3) s . All three crosscorrelograms reveal strong correlations between the residuals of the different coordinate directions.
Following the ideas of Sect. 4.2, crosscorrelograms are also used to model the correlations between two different epochs. Although the time dimension is only indirectly taken into account, those crosscorrelogramsρ x (i) x ( j) ,ρ y (i) y ( j) and ρ z (i) z ( j) (i, j = 2, . . . , k, i = j) are denoted as temporal crosscorrelograms from now on. In Fig. 13, those temporal crosscorrelations between the point clouds PC can be seen: The three curves behave very similarly to the autocorrelations in Fig. 12 (top), with the exception that they do not start at ρ = 1.
An important property of spatial data is its directionality. The correlation structure causing the deformation is assumed to be isotropic in this study (cf. Genton and Kleiber 2015). Nevertheless, anisotropies can be taken into account by using directional correlograms (cf. Sect. 2.3.2).

Setting up the stochastic model of the signal
In order to set up the stochastic model, the empirical correlations are modelled by means of analytical functions (cf. Fig. 1). The respective selection was accomplished according to the requirement of their positive semi-definiteness (Koch et al. 2010). For the available data sets, simple functions have been proven to be sufficient. The analytical functions used are listed in Table 2. Due to the coordinate-wise definition of the signal, the one-dimensional forms of these functions are used. In addition to their positive semi-definiteness, all functions satisfy equation (29) with μ = 0. Thus, the use of these functions guarantees the modelling of a signal with an expectation value of zero.  It has to be noted that the modelling of crosscorrelations is an extensively studied field. (see, for example, Genton and Kleiber 2015 and the references herein.) As the corresponding in-detail analysis is far beyond the scope of this paper, in this first step towards a spatio-temporal deformation analysis the relatively simple models in Table 2 are used. The validity is proven by checking the resulting variance covariance matrix for positive definiteness.
In Fig. 14, the empirical autocorrelations and the empirical temporal crosscorrelations of the z-coordinate of data set PC s can be seen (crosses).
The solid lines represent the estimated analytical functions. For all three illustrated correlograms, the Gauß function was chosen, which, obviously, is an appropriate choice.
By means of these analytical correlation functions, the correlation matrix R ss , consisting of k x k submatrices, can be set up: The submatrices on the main diagonal reflect the correlations within one measuring epoch, whereas the remaining submatrices model the temporal correlations between two measuring epochs. Each of these submatrices of dimensions 3n l × 3n l is structured according to the following schema: For reasons of clarity, the number of observations n l is not equipped with the superscript defining the epoch. Nevertheless, the number of observations may vary between the measuring epochs. The correlations ρ = ρ(d) are computed by means of the respective analytical correlation function, the threedimensional Euclidean distance d between the respective points being the only input of the correlation function.
The resulting submatrices in Eq. (53) and, in further consequence, the correlation matrix in Eq. (52) have an epochal block-wise structure and are sparse, as a large part of the matrix covers the stochastic relationships between two observations of the non-distorted parts, which are modelled by zero-correlations.
Having set up the correlation matrix R ss , it has to be converted into the covariance matrix Σ ss by taking into account the locally stationary variances computed by means of Eq. (50) (see Heunecke et al. 2013, p. 144).
Due to the "regularisation" by means of Σ before inverting Σ ss [cf. Eq. (43)], the slowly decaying correlation functions do not lead to numerical instabilities.

In-depth considerations of the assumptions adopted for the deformation modelling
In the subsections, above three essential assumptions were made. This inserted section justifies these assumptions and, therefore, is an important building block for the development of this approach.
-The normal distribution of the signal [Eq. (3)] is a basic assumption of the general LSC-approach, allowing for the complete modelling of the signal by means of the first two statistical moments, the mean vector μ and the variance-covariance matrix Σ ss . When interpreting deformations as a stochastic signal, this assumption is very restrictive and not necessarily met. However, the simulation process in Sect. 3.1 has illustrated that a normal distributed signal may cause typical deformation patterns. Reversely, a normal distributed signal can be adopted as an approximation of typical deformation patterns. -In Eq. (49), the threshold for distinguishing between distorted and non-distorted parts of the object was heuristically chosen. In order to justify this choice, the The correlogram of the z-coordinate in the lower picture shows a similar behaviour, indicating that those regions which show distortions in the z-coordinate were successfully detected. Contrary, the correlograms of the x-and y-coordinate do not immediately decrease to ≈ 0, but to ρ x (0.005) ≈ 0.4 and to ρ y (0.005) ≈ 0.3. With increasing distance d, these correlations gradually decrease and reach their minimum at d ≈ 0.13, showing a periodic behaviour. Although the magnitude of these correlations does not significantly differ from zero as indicated by the 97,5%-confidence intervals, these correlograms show the existence of remaining correlations in those areas which were classified to be non-distorted. These systematics are caused by the relatively simple approach to detect the distorted parts of the object: the threshold consideration has difficulties to find a strict distinction between distorted and non-distorted areas as there is a smooth transition between those two areas. When choosing a larger threshold, a larger part of the dis-torted areas is allocated to the non-distorted part. This is reflected by a slower decrease in the correlograms. However, when decreasing the threshold, the non-distorted area is more and more polluted by small areas which are classified to belong to the distorted area. For this reason, the choice of Eq. (49) is an appropriate compromise between these two cases: there do not exist any pollutions in the non-distorted parts, and the remaining systematics in the non-distorted areas are small. Nevertheless, the limits of this procedure are revealed by this example: in case of small deformations, which can hardly be distinguished from the noise, a simple threshold consideration might not work sufficiently well. However, for this study, which has to be considered as a first step towards an areal and time-continuous deformation analysis, the threshold considerations will be maintained.
-The purpose of the standardization using Eq. (51) is to represent the homogeneous correlation structure by the resulting homogeneous pseudo-residuals. In order to prove its feasibility, the correlogram is repeatedly calculated, while in each calculation one of the clusters is excluded. Afterwards, an assessment is made whether the resulting correlograms differ significantly. This strategy is chosen, because in case of homogeneity the correlation structure does not change when individual sub-areas of the homogenized distorted region are excluded from the correlogram's estimation. The correlograms' similarity analysis is based on the confidence band of the correlogram of the whole distorted area which is computed according to the parametric spatial bootstrap introduced in Sect. 2.3.4. In case the correlograms of all sub-areas lie within this confidence band, the effectiveness of this approach is validated.
In the upper picture of Fig. 16, the respective results can be seen for the z-coordinate of point cloud PC (3) s in case the number of clusters is chosen to be n c = 12. Apart from small variations, all 12 correlograms are almost identical and lie within the estimated 90%-confidence bands. Thus, the standardized residuals can indeed be assumed to represent the homogeneous correlation structure. The validation of homogeneity is closely linked to the choice of an appropriate number of clusters. This is a critical issue as the number of clusters has a strong influence on the filtering results: in case the number of clusters is chosen too small, the homogeneity assumption is not fulfilled as can be seen in the lower picture of Fig. 16 which shows the results when only two clusters are created. Both correlograms differ significantly. In case the number of clusters is chosen too large, the cluster size becomes so small that a meaningful standardization is not longer possible. The smallest possible number of clusters leading to homogeneous standardized residuals is therefore chosen to be appropriate.

Modelling of the noise
The remaining quantity to be determined before being able to filter the observations is the covariance matrix of the noise Σ . As this study is the first step towards an areal deformation analysis, the noise is initially assumed to be noncorrelating and, thus, Σ is modelled to be a diagonal matrix. The variances of the non-distorted parts result from the trend estimation, whereas for the computation of the variances of the distorted parts the considerations depicted in Fig. 1 have to be taken into account. Theoretically, the wanted variances of the noise would directly result from the estimation of the analytical covariance functions as the difference of the analytical variance and the empirical one. However, due to the normalization of the residuals, no covariance functions can be computed when assuming locally stationary processes. Nevertheless, knowing the ratio of the analytical correlations and the empirical ones for d(0) (cf. Fig. 14) and the local variances of the signal, the respective variances of the noise can be computed.

Filtering results
With the establishment of the stochastic models, the approach introduced in this section can be used to filter the measured point clouds and, therefore, to describe the deformations.

Filtering results for the stochastically simulated data sets
In order to evaluate the results, the residuals of the filtering, which, optimally, solely contain the measuring noise, are analysed in a first step. The residuals' spatial distribution for the z-coordinate of point cloud PC s can be seen in Fig.  17.
Evidently, the magnitude of the residuals randomly varies from approximately −3 mm to 3 mm with few exceptions in the transition zone between the distorted and the nondistorted part of the object (dark blue points). Thus, apart from those outliers, the residuals lie within the 3-σ range of the measuring noise used for the data simulation.
For a more detailed analysis, the residuals of the filtering of point cloud PC (3) s are plotted in the form of histograms (see Fig. 18). It has to be noted that both the residuals of the distorted and those of the non-distorted areas are included in the histograms, whereas the three largest residuals (the outliers coloured dark blue in Fig. 17) are ignored in this representation. The respective residuals in direction of the z-coordinate can be approximated very well by means of a normal distribution, whereas in case of the other two coordinate directions the kurtosis is much too large compared to that of a normal distribution.
A closer look at the statistical moments of the filtering's residuals over time in Table 3 reveals that in case of the non-distorted object (PC (1) s ), the residuals closely follow the normal distribution which was used to generate the noisy data. For the remaining two data sets, this statement is approximately valid for the residuals in the direction of the z-coordinate as the respective standard deviations are only slightly smaller than 1 mm. Skewness and kurtosis give an idea about how close the respective distributions are to a normal distribution (skewness = 0, kurtosis = 3). As already indicated by Fig. 18, kurtosis and skewness confirm a very good approximation of the residuals by means of a normal  distribution in the z-direction, whereas this does not hold for the residuals in x-and y-direction. The residuals' behaviour is caused by the definition of the coordinate system (see for example Fig. 7): as x-and y-direction correspond to the two principal components of the surface, the data sets are very insensitive for determining deformations in these directions.
Comparing the statistical moments of the residuals in the z-direction over time, a decrease in similarity with respect to a normal distribution can be observed when the deformation increases.
The use of simulated data allows a comparison of the filtered data with respect to nominal surfaces, aiming to highlight systematic residual errors. A graphic representation can be seen in Fig. 19 showing the discrepancies between the fil- s . Due to the different treatments of the two areas, a distinction between the distorted and the non-distorted part of the point cloud by means of the values' magnitude is possible: the precision of the measurements is increased by the trend estimation, whereas it is maintained by the filtering. However, the discrepancies behave randomly within these two parts and no pattern corresponding to the deformation can be recognized. It is noticeable, though, that the majority of the points showing the largest discrepancies accumulate in the transition zone between distorted and non-distorted areas. Apart from that, the systematics caused by the deformation process are compensated by means of the stochastic modelling in a satisfying way.
In analogy to the residuals' analysis, the discrepancies to the nominal surfaces are plotted in the form of histograms (Fig. 20). As can be seen, the distribution of the discrepancies can be approximated very well by means of a normal distribution in x-and y-direction, whereas the histogram of the residuals in the z-direction shows a relatively large kurtosis. As, optimally, the discrepancies to the nominal surface equal zero, a large kurtosis is an indicator for a successful modelling of the deformation.
The empirical statistical moments characterizing the distributions of the discrepancies are summarized in Table 4. For all measuring epochs and for all coordinate directions, the distribution of the discrepancies is centred close to zero. Thus, there does not exist a significant bias with respect to the nominal surface.
The standard deviations of the discrepancies give some insight into the discrepancies' variation. The respective values of the first measuring epoch are significantly smaller than the standard deviation of the measuring noise as, due to the estimation of the B-spline surface, the precision is increased. As already indicated by Fig. 19, the precision is maintained rather than increased in the distorted region. That is the reason  Kurtosis (d z ) 6 .8 6 .8 6 .5 why the standard deviations of the second and third measuring epochs are considerably larger than those of the first epoch. In compliance with Table 3, the standard deviations in x-and y-direction are larger than the standard deviation in z-direction, indicating that a detection of the deformation is only successful in the direction perpendicular to the object's surface. This statement is supported by the values of the kurtosis, which are considerably larger for the residuals in direction of the z-coordinate. Finally, the comparison of the statistical moments of the three distorted epochs does not reveal a relationship between the parameters' magnitude and the size of the deformation. For a more meaningful analysis, the filtering computations for data set PC s are repeated 500 times using different realizations of this data set each time.
The results in terms of the mean values and the standard deviations of the first two statistical moments over the 500 repetitions are summarized in Table 5.
The averaged mean values are almost identical with the respective values of point cloud PC s in Table 4 with the averaged mean values in the direction of the z-coordinate being slightly larger in case of the Monte Carlo simulation. However, as the deviation clearly falls within the range defined by the threefold standard deviation, the mean values in Table  4 can be seen to be representative. The small standard deviations of the means being almost equal to the averaged mean values indicate an unbiased reproducibility of the results.
The mean values of the standard deviations take similar values to the standard deviations of data set PC s in Table 4, too. With the mean values of the standard deviation being more than an order of magnitude larger than the standard deviation σ , the results of Table 4 can be seen to be representative.
A deeper analysis reveals that the variation in the results is not only caused by the different data sets, but also by the random choice of the initial cluster centres of the k-means algorithm: even when using the same data set, this random choice leads to a variation in the respective results.

Filtering of the functionally deformed data sets
The computations are repeated using data set PC f . In analogy to the results in the previous section, the results of this computation are summarized by means of Tables 6 and 7 as well as Figs. 21 and 22.  The residuals of the filtering show a similar behaviour as in case of the stochastically simulated data set. The statistical parameters in Table 6 reveal a centring of the respective distributions close to zero and standard deviations which, in case of the residuals in the direction of the z-coordinate, are very close to that one which was used to generate the data. Kurtosis and skewness indicate a good approximation of the residuals' distribution by means of a normal distribution in the direction of the z-coordinate, whereas the distributions corresponding to x-and y-coordinate are leptokurtic.
The residuals' spatial distribution is presented in Fig.  21. As in case of the stochastically simulated data set, the residuals show a random pattern both in the distorted and in the non-distorted parts of the surface and attain maximal values of approximately ±4 mm. However, the transition area between distorted and non-distorted zone is more pro- nounced. It is recognizable in terms of a yellow coloured outline of the distorted region which contains the majority of residuals being larger than 3 mm in absolute value.
This outline is even more apparent in the discrepancies to the nominal surface represented in Fig. 22. This systematics demonstrate the limits of the threshold-based detection of distorted regions. Apart from that, the discrepancies vary randomly and attain smaller values in the non-distorted region and larger values, which are, however, within the scope of the measuring accuracy, in the distorted area.
The resulting statistical parameters are listed in Table 7 showing the same characteristics as the corresponding values for the stochastically simulated data set: on average, the discrepancies are close to zero. The relatively small standard deviation as well as the large kurtosis demonstrates the applicability of the stochastic deformation model to functionally simulated deformations.

Prediction
Usually, not only the filtering of the data is of importance, but also the prediction of values at locations where no measurements were performed, allowing a space-and time-continuous description of the deformation as well as the construction of identical points. The prediction into the non-distorted parts is straightforward by determining the respective position on the B-spline surface. Thus, in the following only the prediction into the distorted parts is treated.

General procedure
Due to the situation's set-up, two cases have to be distinguished: the prediction of the signal in a measured epoch and the prediction of the signal into a non-measured epoch.
The challenge in both cases is that the spatial location on the deformed surface, where the deformation shall be predicted, cannot be specified (x, y, z in equation (31)). However, this information is required for deriving the stochastic relationships which are a function of the spatial distance in case of the correlations and a function of the point's location in case of the variance.
-In case of a prediction into a measured epoch i, the measured point cloud gives information about the deformation. The relationship between the measured points' Cartesian coordinates (x, y, z) (i) and their known position on the trend surface (u, v) (i) is used to determine an approximate position of the point to be predicted. For this purpose, the assumption of small deformations with respect to the surface's dimension is relevant. Starting point is the location on the trend surface defined by the surface parameters u p and v p (lower red point in Fig.  23). Based on these surface parameters, the four nearest neighbours in terms of the measured and projected points are determined (black points in Fig. 23). The relative position of the point to be predicted with respect to these four neighbours is defined by an inverse bilinear interpolation (lower green lines). Applying the same bilinear interpolation to the measured correspondences of the four nearest neighbours defining the deformed surface of epoch i (grey points and upper green lines) yields the initial position of the point to be predicted (upper red point). It has to be noted that this strategy does not yield an orthogonal projection with respect to the trend surface.
Based on this position, on the one hand, the point to be predicted is allocated to the cluster of the nearest measured neighbour. This allocation simultaneously defines the point's locally homogeneous variance defined by equation (50). On the other hand, the correlations with respect to all the measured points can be computed by means of the analytical auto-and crosscorrelograms defined in Sect. 4.6.4. Afterwards, the covariance matrix Σ s s is set up, and the signal is predicted using equation (42). Adding the values of the predicted signal to the initial positions on the B-spline surface results in the final prediction into a measured epoch. -In case of a prediction into a non-measured epoch i, no direct information regarding the object's deformation in this epoch is available. However, as the deformation process is assumed to be continuous, the deformation of the neighbouring epochs is used to derive information about the deformation in the non-measured epoch. For this reason, in a first step, the auto-and crosscorrelograms as well as the locally stationary variances which are needed to predict into the two neighbouring and measured epochs i − 1 and i + 1 are determined. Afterwards, they are used to derive the stochastic relationships of the unmeasured epoch. Figure 24 sketches the case of three measured epochs 1, 2 and 4: two crosscorrelograms are given representing the stochastic dependencies between epoch 1 and epoch 2 (blue curve) as well as between epoch 1 and epoch 4 (green curve). The continuous nature of these curves allows the identification of functional values corresponding to the same d on these curves and, thus, the analysis of how these correspondences change over time (indicated by the black lines). Therefore, in this example, the required crosscorrelogram between epoch 1 and the unmeasured epoch 3 results as the linear interpolant between the two given correlograms (red curve). This strategy is applied to the entirety of auto-and crosscorrelograms, which are needed to describe the signal, as well as on the locally stationary variance information. More complex relationships than the linear interpolation are imaginable in case a larger number of measuring epochs is available. In these cases, the change of the resulting correlograms over time has to be analysed.

Prediction results
Due to the better controllability of functionally simulated deformations, data set PC f is used in the following to evaluate the prediction process. In order to illustrate the prediction results, 2500 locations on the estimated B-spline surface are randomly chosen. The prediction of those points into three measured epochs (t = 1, t = 1.5 as well as t = 2) and into the non-measured epoch (t = 1.75) is realized according to Sect. 6.1. In Fig. 25, the discrepancies of the prediction into the last measuring epoch with respect to the nominal surface can be seen for the z-coordinate. As can be expected, the prediction into the non-distorted part of the object is very precise as this part is described very accurately by means of the trend. The discrep- ancies in the distorted part of the object range between −1.5 and 3 mm. Similar to the filtering results in Fig. 22, a slightly systematic behaviour of the discrepancies is noticeable in the transition zone between distorted and non-distorted regions. However, over the remaining part of the measuring object, the discrepancies behave randomly.
The results for all epochs and all coordinate directions are summarized in Table 8.
Comparing the parameters of the non-measured epoch PC (4) f with those of the measured ones, it can be seen that the former is seamlessly integrated into the measuring sequence: mean values and standard deviations are almost identical for all epochs. As the standard deviation is computed over the entirety of the object and as the distorted parts of the object vary in size, this parameter has only limited relevance. Thus, the table is extended by the minimum and maximum values of the discrepancies, which always refer to the distorted part of the object. These parameters indicate an asymmetric distribution, especially in the direction of the z-coordinate, which is also supported by the values of the skewness. These parameters reflect the systematics which are already indicated in Fig. 25. In spite of these systematics, the prediction results are very promising: the large kurtosis and the standard deviation being significantly smaller than the noise's standard deviation reveal a very satisfying approximation of the nominal surface.
However, it has to be noted that the repetition of the combined filtering and prediction with the same data set happens to yield the non-satisfying results with discrepancies to the nominal surface in the magnitude of 10-15 mm. Due to the random choice of the cluster centres during the k-means clustering, a clustering which does not reflect the variance's locality in an appropriate way is possible. A deterministic approach which selects the cluster centres using the variance information would provide a remedy.

Summary
In this contribution, an areal deformation model for laser scanning data was developed which allows the filtering of measured data as well as the prediction of non-measured data.
The approach's basis is provided by the estimation of a B-spline surface which represents the non-distorted object. Similar to a least squares collocation, the deformation is modelled solely stochastically, allowing for the a-priori determination of rigid body movements based on the B-spline's control points. Due to the introduction of the non-deforming B-spline surface, the optimization problem is solved in the conditional adjustment model, which constitutes a major difference with regard to a least squares collocation.
The deformation is interpreted to be a locally homogeneous process with homogeneous mean values and correlation structures as well as location-dependent variances. Its description requires a multi-step procedure consisting of the distinction between distorted and non-distorted regions, the determination of local variances as well as the estimation of a homogeneous correlation structure.
Due to the stochastic interpretation of the deformation, the challenge of finding identical points in different measuring epochs is circumvented. The developed model allows a pointto-surface-based comparison with respect to the undistorted object as well as a point-to-point-based comparison between two distorted states of the measuring objects. The resulting deformation measures in the form of distances are interpretable and meaningful. Furthermore, the combination of filtering and prediction allows a space-and time-continuous description of the deformed surface.
The approach was applied to different simulated data sets in order to filter them and in order to predict unmeasured values. The comparison of the results with nominal surfaces are promising: even when, as in these simulated cases, the deformation is relatively large compared to the object's dimension, the deformation is modelled with an accuracy in the order of the measuring noise both, in case of the filtering and in case of the prediction.

Limitations of the approach and future investigations
In spite of the promising results, there do exist some limitations which have to be investigated in future. First of all, it is assumed that the non-distorted measuring object has a continuous surface which can be approximated sufficiently well by means of a B-spline surface.
Furthermore, the approach is based on the strong assumption that the deformation process can be described by means of a locally homogeneous stochastic process in a satisfying manner. However, there do exist some types of deformations, such as fractures and other discontinuous deformations, which cannot be assumed to meet this assumption. In those cases, the approach is only applicable with limits, for example, by segmenting the distorted area into sub-areas which can be assumed to represent locally homogeneous deformation processes. The same applies for continuous deformation processes, for which the assumption of a homogeneous correlation structure is too restrictive.
Based on these thoughts, the relationship between different types of deformations (periodic, linear, fractures, etc.) and the resulting correlograms will be investigated in future. These investigations are closely connected to the assessment of the correlograms' influence on the respective filtering and prediction results.
In order to keep the deformation interpretable, the signal is defined in the Euclidean space in this study. Nevertheless, further investigations are required whether the values of the covariance function need to be computed in the domain of the surface parameters.
Apart from these general issues, the approach's implementation posed some concrete problems, which will be solved in future and which should further reduce the current discrepancies between filtered/predicted values and the nominal surfaces: -The simple threshold consideration for detecting distorted regions is applicable only in limited situations. Thus, a more sophisticated approach for distinguishing between distorted and non-distorted areas will be developed in future. -The randomness of the k-means clustering causes the non-reproducible-and in very few cases unsatisfactoryresults. A deterministic approach for defining the locally stationary areas promises more stable results and, as a consequence, more meaningful Monte Carlo simulations. -The computation and definition of suitable crosscorrelation functions are wide topics for themselves, which require an in-detail analysis. -In this study, only a brief glance was taken at the prediction and the applicability was demonstrated solely for a linear deformation process. In future, the influence of a linear treatment of nonlinear deformation processes will be investigated. These investigations will also devote the question whether it is possible to use the timely change of the correlograms in order to derive information about the nonlinearity of the deformation process.
Furthermore, for these initial investigations, the framework was restricted by using simulated data. Thus, in further studies, the approach will be tested on real measuring data which poses new challenges: -The simulated data sets allow the use of nominal values for the B-spline's surface parameters, knot vectors and degrees. When using real measuring data, those parameter groups have to be determined in an appropriate way. -Up to now, no attention was given to the stochastic modelling of the measuring noise: in the simulation process, white noise was generated so that the identity matrix was an appropriate choice. However, in reality the measuring process itself induces systematics into the point clouds which have to be modelled in an appropriate way by means of Σ , as otherwise the signal absorbs these systematics and the modelling of the deformation is distorted.
Finally, the observation of large measuring objects over long periods of time results in a huge amount of data. A sequential variant of the proposed approach is one possibility to deal with such a huge amount of data. right holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.