On the impact of correlations on the congruence test: a bootstrap approach

Abstract

The detection of deformation is one of the major tasks in surveying engineering. It is meaningful only if the statistical significance of the distortions is correctly investigated, which often underlies a parametric modelization of the object under consideration. So-called regression B-spline approximation can be performed for point clouds of terrestrial laser scanners, allowing the setting of a specific congruence test based on the B-spline surfaces. Such tests are known to be strongly influenced by the underlying stochastic model chosen for the observation errors. The latter has to be correctly specified, which includes accounting for heteroscedasticity and correlations. In this contribution, we justify and make use of a parametric correlation model called the Matérn model to approximate the variance covariance matrix (VCM) of the residuals by performing their empirical mode decomposition. The VCM obtained is integrated into the computation of the congruence test statistics for a more trustworthy test decision. Using a real case study, we estimate the distribution of the test statistics with a bootstrap approach, where no parametric assumptions are made about the underlying population that generated the random sample. This procedure allows us to assess the impact of neglecting correlations on the critical value of the congruence test, highlighting their importance.

1 Introduction

Statistical tests for deformation provide a way to detect potential risks linked with natural hazards or the collapse of artificial objects, such as bridges. One of their main applications in geodesy is the estimation and minimization of risks. Combined with, for example, the concept of utility theory, they allow their consequences and costs to be estimated (Zhang and Neumann 2014).

Terrestrial laser scanners (TLS) capture a large amount of three-dimensional (3D) points rapidly, with high precision and spatial resolution. In addition to the visualization of 3D topography or buildings, TLS are, thus, increasingly used to detect and quantify displacements or deformations. Software such as CloudCompare (www.danielgm.net/cc), 3DReshaper (Hexagon Metrology) or Geomagic Studio (3DSystems) provide maps of deformation that allow one to perform a first deformation analysis. Such software, however, does not perform congruence tests, which necessitate an adequate parametric approximation of the point clouds (Niemeier 2002). Such an approximation can be obtained by fitting B-spline surfaces to the TLS point clouds with a least-square (LS) adjustment (Koch 2009). Once the mathematical surface has been obtained, an adapted version of the congruence test based on the gridded surface points can be applied (Zhao et al. 2019).

The congruence test is known to be the most powerful test in Gauss-Markov models with normally distributed random deviations and a correctly specified stochastic model (Pelzer 1971). To make use of its full potential for risk assessment, the setting of a correct variance covariance matrix (VCM) of the errors has, thus, to be carried out as correctly as possible. One way is to use an empirical model for computing the variances of the errors of the raw TLS observations. The datasheet from the manufacturers provides a good approximation for the angle variances (Boehler and Marbs 2002; Cosarca et al. 2009). The range variance can be modelled with an empirical intensity model, as proposed in Wujanz et al. (2017) and Wujanz et al. (2018) and simplified in Kermarrec et al. (2018). Mathematical correlations must be accounted for because TLS point clouds approximated with B-spline surfaces have to be transformed into Cartesian co-ordinates. Their impact on the test statistic for deformation between B-spline surface approximations has been studied in previous contributions of the authors: it could be shown that neglecting mathematical correlations affects the congruence test in the case of small deformation and unfavourable scanning geometries (see Zhao et al. 2019 for simulations; Kermarrec et al. 2019a for real data analysis). These results confirmed theoretical derivations (Hoog and Craig 1978) by additionally quantifying the expected impact on the test statistics of neglecting mathematical correlations.

Temporal correlations between TLS raw observations are often disregarded due to the associated computational burden (Kauker and Schwieger 2017). In the same way as mathematical correlations, they have to be taken into consideration to avoid biased test statistics. A usual approach would make use of variance component estimation (Amiri-Simkooei 2007; Searle 1995; Teunissen and Amiri-Simkooei 2008) based on the residuals of the LS adjustment from the B-spline surface fitting. Unfortunately, this procedure may lead to numerical inaccuracies, particularly for small samples: matrices inversion inaccuracies translate into the detection of a non-existing deformation. An alternative is to make use of the residuals by modelling their correlation structure, i.e. estimating the parameters of a known function, such as the Matérn function (Mátern 1960; Stein 1999). Such a statement suggests approximating geometrical and physical effects simultaneously within one function only: we are aware that this strong assumption necessitates deep foundations to identify when and if it holds true. This is the first aim of the present contribution, which provides simple and understandable tools to analyse the to analyse the spectral content of the residuals using an empirical mode decomposition (EMD; Flandrin et al. 2004, 2005; Huang et al. 1998).

Once the temporal correlation of the residuals is modelled, the assumption of the Chi square distributed test statistics of the congruence test can be investigated. This is our second goal: deviations from the theory are likely to happen, particularly when temporal correlations are neglected. The correct critical value under a given confidence level can be empirically derived following Efron and Tibshirani (1993) using a bootstrapping approach. In a first approximation, we propose fitting the distribution of the test statistic with a Gamma distribution, which includes the Chi square distribution as a particular case. We will make use of real TLS observations from a bridge under load to highlight the deviations from the expected theoretical distribution (Kermarrec et al. 2019a).

The paper is organised as follows: in the second section, the principle of B-spline approximation is presented briefly, along with the test statistics for deformation, with a focus on the bootstrapping procedure. The third section is devoted to the case study, where the theoretical results are applied to observations from a bridge under load. The distribution of the test statistics for deformation is derived with a Monte Carlo approach. We conclude this contribution with some recommendations on how to account for correlations for deformation analysis based on B-spline surfaces.

2 Mathematical background

2.1 B-spline surfaces

Scattered and noisy raw point clouds from TLS can be approximated mathematically by means of a linear combination of basic functions, such as B-splines. The mathematical surface in its parametric formulation can be expressed as

$$ {\mathbf{s}}\left( {u,v} \right) = \left( {x_{i,j} ,y_{i,j} ,z_{i,j} } \right) = \sum\limits_{i = 0}^{n} {\sum\limits_{j = 0}^{m} {N_{i,p} } } \left( u \right)N_{j,q} \left( v \right){\mathbf{p}}_{i,j} $$

(1)

where $\left( {u,v} \right) \in \left[ {0,1} \right] \times \left[ {0,1} \right]$ are the parameters in the two directions, so that a B-spline surface maps the unit square to a rectangular surface patch. The choice of the parameters affects the shape of the surface. The uniformly spaced or chord length methods (Bureick et al. 2016; Piegl and Tiller 1997) provide good results to determine $u$ and $v$ for regularly and rectangular-shaped point clouds with rows and w columns. Once a set of parameters is obtained, the knot vector can be obtained over which the basic functions $N_{i,p}$ and $N_{j,q}$ of degree p and q, respectively, are computed. Different methods can be employed; the simplest one being to use a uniformly spaced knot vector or to average the parameters (de Boor 2001). The basic function $N_{i,p} \left( u \right)$ is, a composite curve of degree p polynomials with joining points at knots in the interval $\left[ {u_{i} ,u_{i + p + 1} } \right)$ and can be evaluated by means of a recurrence relationship. ${\mathbf{p}}_{i,j}$ is called the control points vector, with i varying from 0 to n and j from 0 to m. Given the degree of the basic functions, approximating a surface is simply looking to find the control points such that the distance of the data points to the approximated surface is minimized. This step can be performed within the context of a LS approximation. If ${\mathbf{x}}$ of size $\left( {\left( {n + 1} \right) \times \left( {m + 1} \right),3} \right)$ is the matrix containing the co-ordinates of the control points to be estimated and ${\mathbf{A}}$$\left( {3n_{obs} ,3\left( {n + 1} \right)\left( {m + 1} \right)} \right)$ the deterministic design matrix containing the values of the tensor product of the B-spline functions, ${\mathbf{v}}$ is the difference between the Cartesian co-ordinates of the point cloud ${\mathbf{l}}$ of size $\left( {n_{obs} ,3} \right)$ and ${\mathbf{Ax}}$. Interested readers should refer to Zhao et al. (2017) for the description of the design matrix.

We assume $E\left( {\mathbf{l}} \right) = {\mathbf{Ax}}$, where $E( \cdot )$ is the expectation operator. We further call ${\varvec{\Sigma}}_{ll}$, the VCM of the errors in Cartesian co-ordinates with ${\varvec{\Sigma}}_{ll} = \sigma_{0}^{2} {\mathbf{I}}$, where $\sigma_{0}^{2}$ is an a priori variance factor and ${\mathbf{I}}$ the identity matrix of size $\left( {3n_{obs} ,3n_{obs} } \right)$. The estimated co-ordinates of the control points are expressed by the unbiased ordinary LS estimator (Koch 1999)

$$ {\hat{\mathbf{x}}} = \left( {{\mathbf{A}}^{T} {\varvec{\Sigma}}_{ll}^{ - 1} {\mathbf{A}}} \right)^{ - 1} {\mathbf{A}}^{T} {\varvec{\Sigma}}_{ll}^{ - 1} {\mathbf{l}} = \left( {{\mathbf{A}}^{T} {\mathbf{A}}} \right)^{ - 1} {\mathbf{A}}^{T} {\mathbf{l}} $$

(2)

The a priori VCM of the estimated parameters is given by

$$ {\varvec{\Sigma}}_{{{\mathbf{\hat{x}\hat{x}}}}} {\mathbf{ = }}\left( {{\mathbf{A}}^{T} {\varvec{\Sigma}}_{ll}^{ - 1} {\mathbf{A}}} \right)^{ - 1} $$

(3)

The number of control points to estimate can be iteratively adjusted by using information criteria (e.g. Alkhatib et al. 2018 for the specific application of information criteria to B-spline surface approximation and the references inside). Two possibilities exist: (1) the Akaike information criterion (AIC), which minimizes the Kullback–Leibler divergence of the assumed model from the data-generating model, or (2) the Bayesian information criterion (BIC), which assumes that the true model exists and is, thus, more adequate for large samples.

2.2 Stochastic model

The estimated co-ordinates of the control points are not affected by the choice of the estimated ${\hat{\mathbf{\Sigma }}}_{ll}$ instead of ${\varvec{\Sigma}}_{ll}$, due to the unbiasedness of the LS estimator. However, the ordinary LS estimator is no longer the most efficient within the class of linear unbiased estimators when ${\hat{\mathbf{\Sigma }}}_{ll}$ deviates from the true VCM and hypothesis tests, such as the global test, outlier tests or congruence tests, become invalid (Williams et al. 2013). The assumption of homoscedasticity has to be adapted in order to derive realistic test statistics based on Eq. (3). In this contribution, we propose to make use of the residuals from the B-spline surface approximation to approximate ${\hat{\mathbf{\Sigma }}}_{ll}$. The VCM of the Cartesian observation errors will, thus, be considered as the VCM of the residuals ${\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}}}}$, with ${\mathbf{\hat{v} = l - A\hat{x}}}$ the residuals of the LS adjustment. Because the residuals are correlated and heteroscedastic, ${\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}}}}$ is fully populated. The global VCM of the three components is set up by neglecting the correlations between components:

$$ {\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}}}} = \left[ {\begin{array}{*{20}c} {{\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}},x}} } & {\mathbf{0}} & {\mathbf{0}} \\ {\mathbf{0}} & {{\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}},y}} } & {\mathbf{0}} \\ {\mathbf{0}} & {\mathbf{0}} & {{\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}},z}} } \\ \end{array} } \right] $$

(4)

As a powerful alternative to empirical VCE, we will assume that the correlation structure of the residuals can be approximated with a parametric Matérn covariance function (“Appendix 1”). We assume, thus, that only a coloured noise is present in the residuals, which is a strong assumption. The impact of geometrical misspecifications on the residuals has to be investigated prior to a fitting with such a parametric covariance function. One way is to make use of maximum likelihood-based model selection, which, as indicated by its name, is based on predefined models (often a combination of coloured and white noise). To our point of view, such procedures answer only conditional the question whether a parametric covariance function is a relevant choice or not; their result may be biased in case of discrepancies from the noise assumption due to model misspecification. As a less computational demanding alternative, we propose making use of the analysis of the power spectral decomposition of the residuals combined with an EMD (“Appendix 2”) to validate the use of such a model, similarly to denoising strategies (Kopsinis and Mclanglin 2009). We propose for didactical purposes to highlight our methodology with a practical example in Sect. 3.

Three parameters will have to be estimated: the variance, the range and the smoothness of the correlation function from the residuals. A possibility to make a computationally efficient estimation of the parameters in the case of small samples is described in “Appendix 3”. The most important parameter is the smoothness, which it impacts strongly on the first values of the inverse of the VCM and, thus, the value of the test statistics (Kermarrec and Schön 2016).

2.3 Congruence test for deformation analysis

In order to test for deformation between two epochs, we follow Zhao et al. (2019) and make use of an adapted form of the global congruency test derived by Pelzer (1971). No difference between the two vectors of control points co-ordinates (the “parameters”) can be computed because the number of control points estimated at different epochs may vary. Therefore, we test the deformation at the level of the estimated surfaces and not at the parameter level (the control points).

Using gridded observations and assuming a deformation only in the z-direction, we only consider ${\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}},z}}$ from Eq. (4) and determine the surface difference. The surface points at a given epoch are estimated by a corresponding matrix product ${\mathbf{F\hat{x}}}$ from Eq. (1), where the ${\mathbf{F}}$ is evaluated from the tensor product of B-spline functions estimated at the parameters of the gridded points.

The uniformly most powerful invariant test for deformation under the null hypothesis can be shown to be based on the test statistic:

$$ \begin{aligned} T_{apriori} & = {\hat{\mathbf{\Delta }}}^{{\mathbf{T}}} {\varvec{\Sigma}}_{{{\mathbf{\hat{\Delta }\hat{\Delta }}}}}^{{{\mathbf{ - 1}}}} {\hat{\mathbf{\Delta }}} = \frac{1}{{\sigma_{0}^{2} }}{\hat{\mathbf{\Delta }}}^{{\mathbf{T}}} {\mathbf{Q}}_{{{\mathbf{\hat{\Delta }\hat{\Delta }}}}}^{{{\mathbf{ - 1}}}} {\hat{\mathbf{\Delta }}} \sim \chi_{p}^{2} \\ {\varvec{\Sigma}}_{{{\mathbf{\hat{\Delta }\hat{\Delta }}}}}^{{}} & = {\mathbf{H\varSigma }}_{{{\mathbf{\hat{\beta }\hat{\beta }}}}}^{{}} {\mathbf{H}}^{{\mathbf{T}}} \\ \end{aligned} $$

(5)

where ${\varvec{\Sigma}}_{{{\mathbf{\hat{\beta }\hat{\beta }}}}}^{{}}$ is the VCM of the estimated CP for both epochs and ${\hat{\mathbf{\beta }}} = \left[ {\begin{array}{*{20}c} {{\hat{\mathbf{x}}}_{{{\mathbf{epoch1}}}} } \\ {{\hat{\mathbf{x}}}_{{{\mathbf{epoch2}}}} } \\ \end{array} } \right]$ contains the LS estimates of the CP at epoch 1 and 2. The matrix ${\mathbf{H}}$ is split into two parts to determine a difference at the level of the surface points: ${\mathbf{H}} = \left[ {\begin{array}{*{20}c} {{\mathbf{ - F}}_{{{\mathbf{epoch1}}}} } & {{\mathbf{F}}_{{{\mathbf{epoch2}}}} } \\ \end{array} } \right]$ with p rows. The difference between the two surfaces is given by ${\hat{\mathbf{\Delta }}} = {\mathbf{H\hat{\beta }}}$, where we here only consider the difference in the z-direction.

We define the hypotheses for deformation analysis as follows:

$$ H_{0} :{\varvec{\Delta}} = 0\quad {\text{vs}}\quad H_{1} :{\varvec{\Delta}} \ne 0 $$

The null hypothesis states that no deformation happened. The test decision involves the quantile $k_{{1 - {\raise0.7ex\hbox{${\alpha_{test} }$} \!\mathord{\left/ {\vphantom {{\alpha_{test} } 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}}}}^{{\chi_{p}^{2} }}$ based on the $\chi_{p}^{2}$ test distribution and the significance level $\alpha_{test}$. In particular, $H_{0}$ is accepted if $T_{apriori} \le k_{{1 - {\raise0.7ex\hbox{${\alpha_{test} }$} \!\mathord{\left/ {\vphantom {{\alpha_{test} } 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}}}}^{{\chi_{p}^{2} }}$.

The test statistic $T_{apriori}$ is strongly dependent on the VCM of the CP estimates, which is based on ${\varvec{\Sigma}}_{{{\mathbf{\hat{\beta }\hat{\beta }}}}}^{{}}$, see Eq. (2). Consequently, any misspecification will affect the decision step, potentially leading to the erroneous conclusion that no deformation occurs or, contrarily, that a deformation occurs (error of type I and II). The results are particularly sensitive in the presence of correlated errors: they may lead to the detection of deformations which are not the consequence of an existing deformation but only the result of the correlated noise structure of the observations.

We will consider two cases for the VCM in the following:

• (i) ${\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}}}}$ leading to $T_{apriori}^{(i)}$

• (ii) ${\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}}}}^{diag}$ leading to $T_{apriori}^{(ii)}$

The VCM ${\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}},z}}^{diag}$ is determined by multiplying $diag\left( {{\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}},z}}^{{}} } \right)$ component-wise with the identity matrix of corresponding size.

The full methodology is summarized in a flowchart form in Fig. 1.

Fig. 1

Methodology of the congruence test using an estimation of the covariance structure from the residuals of the LS approximation

Additional comment

Our approach is not based on a physical model for the raw (polar) TLS observations. In such cases, mathematical correlations would have to be additionally considered due to the transformation from potentially correlated and heteroscedastic polar raw measurements and Cartesian co-ordinates.

2.4 Validation of the test distribution

The defined a priori test statistic is expected to follow a $\chi_{p}^{2}$ distribution. This result remains theoretical and may not hold when an approximated stochastic model is used or when correlations are present. A way to validate the assumed distribution is to make use of a bootstrapping-based hypothesis test. This elegant procedure allows one to estimate the parameters of a general Gamma distribution, which includes the $\chi_{{}}^{2}$ as a particular case (Kargoll et al. 2018; Appendix A.2).

In the following, we will briefly present the methodology used for validation, which is further summarized in a flowchart in Fig. 2. Our procedure involves three steps, in the same way as Kargoll et al. (2019):

Fig. 2

Flowchart of the computational steps to derive the distribution of the test statistics using a Monte Carlo approach

• 1. Testing: the first step starts with the approximation of scattered (TLS) observations from 2 epochs with B-spline surfaces.

• 2. Generating: in the kth-generating step, we add a noise vector, whose structure corresponds to the positive definite ${\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}}}}$, to the two approximated surfaces. We use a Cholesky decomposition and decompose ${\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}}}}$ as ${\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}}}} = {\mathbf{G}}^{{\mathbf{T}}} {\mathbf{G}}\,$. We further generate a Gaussian random vector ${\mathbf{NO}}_{i,k} ,i = 1,2$ for the two epochs with mean 0 and variance 1 from the Matlab random number generator randn. The correlated noise vector reads: ${\mathbf{N}}_{i,k} = {\mathbf{G}}^{T} {\mathbf{NO}}_{i,k}$. We add this generated noise vector to its corresponding surface and approximate it by a B-spline surface. We compute the a priori test statistic ${\text{T}}_{apriori}$ from the two estimated surfaces for each step. For one iteration k, we call the corresponding test statistic ${\text{T}}_{apriori}^{k}$.

• 3. Evaluation: we perform $k_{total} = 1000$ iterations. The mean of the $k_{total}$ test statistics values in the evaluation step are statistically evaluated using the function fitdist from Matlab. We compute the parameters of a Gamma distribution—the $\chi_{{}}^{2}$ being a particular case corresponding to a scale parameter of 2.

In the next section, we propose applying this theoretical derivation to a case study. We will make use of the same TLS observations as in Kermarrec et al. (2019a). The focus of this contribution is on the impact of neglecting the correlations of the Cartesian residuals on the test statistics. We investigate particularly how they may affect (i) the test statistics and (ii) their corresponding test distribution.

3 Case study

In this section, we use TLS observations from a bridge under load and test for the deformation between two steps of load. We start from a gridded TLS point cloud, which we approximate with B-spline surfaces. Following the methodology described in the previous section, we estimate the correlation structure and the variance from the residuals of the B-spline approximation, which we use to set up the corresponding Toeplitz VCM, assuming stationarity of the variance. The test statistics are computed using both the fully populated VCM of the residuals (case (i)) and the diagonal simplification (case (ii)), i.e. neglecting correlations. After having briefly described the B-splines surface approximation we investigate the impact on the a priori test statistics of the two stochastic models. In a last step, we validate the distribution of the test statistics with the proposed Monte Carlo approach. The justification of the fitting with a parametric covariance function will be didactically explained.

3.1 Description of the data set

We use a real data set of a historic masonry arch bridge over the river Aller near Verden in Germany made of circular brick arches. The aim of the experiment was the combination of numerical models and experimental investigations for model calibration (Paffenholz et al. 2018).

The standard load of 1.0 MN (100 t) was defined and further loadings up to five times the standard load were realised in the scope of the load testing. Thus, a maximum load of approximately 6.0 MN was defined, produced by four hydraulic cylinders mounted on the arch. In the following, we will call the different load steps E00 (reference), E11, E22, E33, E44 and E55. The deformations will always be computed regarding the reference load, which we will denote E11-00, E22-00, E33-00, E44-00, E55-00.

The 3D point cloud acquisition was carried out using TLS sensors of type Zoller + Fröhlich (Z + F) Imager 5006/h in periods of a constant load on the bridge. The 3D point clouds for different load scenarios ranging from 1 up to 6 MN were captured and finally processed. A filtering regarding objects on the arch surface was performed to eliminate interfering objects, i.e. other sensor installations, such as prisms for the laser tracker and strain gauges.

Following Kermarrec et al. (2019a), we selected small patches of observations in order to have an optimal functional model. Using the software CloudCompare, the same surfaces were chosen for each load. The patches chosen were located around two reference Laser Tracker points L8 and L10, as shown in Fig. 3. L10 can be considered as optimally scanned at a short distance in the Up-direction, whereas L8 was scanned with a less favourable geometry. Moreover, L10 has a stronger deformation magnitude than L8 (approximately 10 vs. 4 mm, see Paffenholz et al. 2018).

Fig. 3

A bridge under load: localization of the points L10 and L8

3.2 B-spline surface approximation

The TLS observations were gridded in 15 cells in both directions leading to a total of 225 gridded observations. This number of observations allows for a reliable estimation of the correlation parameters using the Whittle likelihood methods, as presented in “Appendix 3”. The two surfaces are shown in Fig. 4 (left), where it can be seen that they correspond approximately to a plane. All observations falling in one cell were averaged for the gridding. The parameterization was carried out using a uniform method, which is efficient for the uncomplicated geometries under consideration. Similarly, the knot vector was taken as equidistant. The number of CP was determined using the Bayesian information criterion, which gave similar results to the Akaike information criterion. For each dataset, 4 and 3 CP in the u- and v-direction, respectively, were found to be optimal for L8; 5 CP in both directions were considered for L10. The same values were found for all epochs under consideration.

Fig. 4

Top left: the gridded point clouds (yellow) and the approximation (green) for L8. Bottom left: the z-residuals. Top right: scaled Toeplitz VCM ${\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}},z}}$ estimated from the z-component. Bottom right: first ten values of the first row of ${\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}},z}}^{ - 1}$

3.3 Estimation of the VCM

We wish to estimate the VCM of the residuals in the z-direction from the residuals of the B-spline surface approximation and assume that a Matérn model can describe the stochasticity of the residuals.

As mentioned in “Appendix 1”, any Matérn processes will have a typical power spectral decomposition with two different regimes (Lilly et al. 2017):

• at high frequencies, a power law decay governed by the smoothness parameter $\nu$, and

• at low frequencies, a constant value, i.e. a “locally white” behaviour governed by the parameter $\alpha$.

The power spectral decomposition obtained with the Matlab function pwelch is plotted in Fig. 5 (top left). Due to the small sample size under consideration, the identification of the classic Matérn spectral is far from evident, although two different modes can be guessed. In order to justify the use of a Matérn model, we made use of the EMD of the residuals (“Appendix 2”). Since the latter are not smooth, we use ‘pchip’ as an interpolation method for the Matlab function EMD and decompose them in five modes. The corresponding decomposition is shown exemplarily for L8 in Fig. 5 (bottom left).

Fig. 5

Top left: the power spectral decomposition for L8 and L10 residuals (red and green line, respectively) and the decomposition for the IMF2 + 3 of the same residuals as the dashed line. Bottom: the EMD for the L8 residuals. Top right: $\log_{2} \left( {E_{k} } \right)$ versus the IMF $k$

We define $E_{k}$ from the energy-based variance of the IMF, following Flandrin et al. (2004). It was proven from Flandrin et al. (2005) or Wu and Huang (2004) that the mean energy of the different intrinsic model function (IMF) $E_{k}$ from a pure power law noise follows a linear decrease in a semi-log diagram, i.e. $\log_{2} \left( {E_{k} } \right)$ regarding the IMF index $k$. The corresponding plot is shown in Fig. 5 (top right). We observed that after the fourth IMF, the energies diverge significantly from the theoretical model, indicating the presence of significant amounts of “no-coloured noise” signal. This finding is supported by the shape of the corresponding EMD (Fig. 5, bottom left).

Following Flandrin et al. (2005) further, we filtered the noise by excluding the fourth and fifth IMF and interpret these as corresponding to some unknown model misspecifications. Since the ratio of the variances $\frac{{\sigma_{IMF1 + IMF2 + IMF3}^{2} }}{{\sigma_{IMF4 + IMF5}^{2} }} \approx 25$, we, furthermore, chose to neglect these two IMF and concentrated on the third first one, which corresponds to the coloured noise. The first IMF carries primarily the white noise, which we validated with the Kolmogorov–Smirnov test. Consequently, and in order to facilitate the estimation of the Matérn parameters from $IMF2 + IMF3$, we modelled ${\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}},z}}$ as ${\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}},z}} = \sigma_{IMF1}^{2} {\mathbf{I}} + \sigma_{IMF2 + IMF3}^{2} {\hat{\mathbf{\Sigma }}}_{IMF2 + IMF3}$, with ${\mathbf{I}}$ as the identity matrix. We estimated the corresponding variances $\sigma_{IMF1}^{2} ,\sigma_{IMF2 + IMF3}^{2}$ from IMF1 and IMF2 + IMF3, respectively. ${\hat{\mathbf{\Sigma }}}_{IMF2 + IMF3}$ is modelled as a Toeplitz VCM and its corresponding Matérn parameters are estimated using the methodology presented in “Appendix 3”. The power spectral decomposition of IMF2 + IMF3 is shown as a dotted line in Fig. 5 (top left): the two different regimes are now clearly identifiable and support an unbiased estimation of the parameters.

We summarize our methodology in 4 steps:

• compute the z-residuals from the B-spline approximation;

• perform the EMD. Plot the mean energy with respect to the IMF in a log diagram. Identify the power law noise IMF (straight line). If at that step, model misspecifications dominate, no straight line will be identifiable (LS fitting can be used for that purpose). We can still consider the second IMF as being the coloured noise. However, in that case, a power law model only, such as the Matérn covariance function, will be suboptimal, particularly if the variance of the third IMF is still high;

• separate the white noise IMF from the power-law noise and estimate the corresponding Matérn parameters; and

• compute ${\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}},z}}$ assuming stationarity of the variance, which can be checked with IMF1.

The Toeplitz structure of ${\hat{\mathbf{\Sigma }}}_{{{\mathbf{\hat{v}\hat{v}}},z}}$ is shown exemplarily in Fig. 4 (top right) for L8 E00. The small level of correlation is highlighted by the sparse structure of the VCM. Additionally, we show the structure of the inverse of the VCM with the first ten values of one row of the inverse. After a large value in the diagonal (three times higher than the second one), the following nine values oscillate around and close to 0. Compared with results obtained from diagonal matrices, this particular structure has a strong effect on the test statistics, as mentioned in Kermarrec and Schön (2016). The first value is directly linked with the smoothness of the corresponding process: the power spectral decomposition at high frequency should permit a good estimation of this parameter, i.e. a linear decay should be found (Fig. 5, top left, dashed line). The approximation done by neglecting the IMF4 and 5, which corresponds to the low frequency noise, is, thus, fully justified, as well as the splitting between the white noise IMF and the “coloured noise” IMF.

We present the values of the Matérn parameters obtained for L8 and L10 (E00) in Table 1. For the sake of brevity, the other similar values are not shown.

Table 1 $\left[ {\sigma_{IMF1} ,\alpha ,\sigma_{IMF2 + 3} ,\nu } \right]$ for L8 and L10 are computed for the first load step E00. The standard deviation is given in [mm]

Larger values of the smoothness and range parameters are found for L8, which is scanned under an unfavourable geometry, compared to L10, for which the correlations are less pronounced. The standard deviations of the power law noise are, however, similar for both points. The amount of white noise is slightly higher for L10, which we justify by the optimal scanning geometry.

3.4 Impact on the test statistic

Having computed the estimated fully populated VCM from the residuals for each load step following the previous methodology, we highlight their impact on the test statistics in Fig. 6. To that aim, we follow the methodology proposed in Sect. 2.3 and compare the results with the one obtained with a diagonal VCM. We notice that a low level of correlation (L10) combined with a high deformation magnitude leads to a difference between the two stochastic models smaller than for L8. In the latter case, a smaller deformation magnitude together with a higher correlation level makes the difference between $T_{apriori}^{(i)}$ and $T_{apriori}^{(ii)}$ larger. Whereas for L10, $T_{apriori}^{(ii)} \approx 1.2T_{apriori}^{(i)}$, the factor increases to $T_{apriori}^{(ii)} \approx 1.35T_{apriori}^{(i)}$ for L8. We interpret the small difference as due to the relatively low level of correlations. However, the factor for L8 is still slightly higher than for L10, as a consequence of the parameter estimates presented in Table 1.

Fig. 6

Values of the test statistics $T_{apriori}^{(i)}$ and $T_{apriori}^{(ii)}$ for the five deformations E00-11, E00-22, E00-33, E00-44 and E00-55. Top: L8, Bottom L10. The red line corresponds to the diagonal and the blue line to the fully populated VCM, respectively

The main result is the overestimation of the test statistics $T_{apriori}^{(ii)}$ regarding $T_{apriori}^{(i)}$: in all cases, the ratio $\frac{{T_{apriori}^{(ii)} }}{{T_{apriori}^{(i)} }}$ is higher than 1. As mentioned previously, this has to be linked with the particular structure of the inverse of the VCM, which acts as a variance factor: accounting for correlations in the test statistics is similar to having decreased the variance of a corresponding diagonal VCM (Kermarrec et al. 2019b). Even if in our particular case study, the a priori test statistics are in all cases far over the critical value of the Chi square distribution with a confidence level of 0.05, this behaviour should not be underestimated for more subtle cases, where the deformation to detect is much smaller. When correlations are present and a correct VCM is used, the test statistics will become smaller and, thus, more realistic: this will lead to accepting the null hypothesis “no deformation”. This is an important point, as correlations present in the error term should not be considered as true deformation. We, thus, strongly encourage the use of a realistic VCM which can be estimated from the residuals. However, these conclusions assume that the Chi square distribution can be used, which needs to be validated or not.

3.5 Test distribution

Using the methodology described in Sect. 2.4., the ratio about $T_{apriori}^{(i)}$ and $T_{apriori}^{(ii)}$ is computed using a Monte Carlo approach. The corresponding histograms and a fitted gamma distribution are presented in Fig. 7. The parameters of the distribution are estimated using the Matlab function fitdist. We justify the choice of the gamma distribution as it includes the assumed Chi squared distribution as a particular case. Further analysis will test other member of the distribution family, such as the exponential.

Fig. 7

(Left) histogram and corresponding Gamma distribution for case (i) corresponding to an optimal stochastic model including correlations. Right: the same for case (ii), i.e. a diagonal VCM neglecting correlations. The observations are taken from L8, deformation E00-11

For the sake of brevity, the results correspond to the particular case of L8 and the deformation step E00-11. The shape parameters of the Gamma distribution were found to be 3.4 and 3.5 for cases (i) and (ii), respectively. The scale parameters were 575 and 370, respectively. Thus, these parameters are different from those of a Chi square distribution, for which the scale parameter would equal 2. Correspondingly, the critical values for a given significance level do not correspond to the assumed one: they are much higher. Exemplarily, we found a critical value of 5731 with a significance level of 0.05 for (i), whereas a value of 3750 has to be considered for (ii). The ratio of the two critical values is, thus, 1.5 and larger than $\frac{{T_{apriori}^{(ii)} }}{{T_{apriori}^{(i)} }}$. Therefore, the risk of making a wrong testing decision by neglecting correlations is slightly decreased but still exists. Because the spread of the test statistics is, however, much higher in that case, we encourage not disregarding correlations in the testing procedure.

The other cases (deformation steps as well as L10) led to similar conclusions. The parameters of the Gamma distribution, however, changed from case to case and needed to be revaluated for each data set to determine the correct critical value for a given confidence level.

The parameters of the Gamma distribution for L10 were found to be similar for both cases due to the low level of correlations, which are in line with the conclusion of the previous section. The values of the scale parameters were found to be 515 and 585 for case (ii) and (i), respectively, for the deformation step E00-11. We note that they remain between 500 and 600, the same as for L8. A generalization of this finding necessitates further investigations that we propose to carry out in a next step.

To conclude, we briefly resume the proposed methodology to perform a rigorous test for deformation when correlations are suspected:

• 1. from two TLS point clouds, estimate the control points of the two B-spline surfaces with a LS adjustment (Sect. 2.1);

• 2. compute the residuals;

• 3. validate whether a Matérn model is adequate (Sect. 3.3);

• 4. if yes, estimate the corresponding Toeplitz VCM;

• 5. compute the a priori test statistic following Eq. (5); and

• 6. perform a Monte Carlo analysis to derive the correct critical values of the congruence test (Fig. 2).

4 Conclusion

An accurate and realistic knowledge of the stochastic properties of TLS errors is of great importance to avoid untrustworthy test decisions. Standard software allows the visualization of deformation by means of maps of distances, but not for the rigorous statistical testing of deformation. This latter can only be achieved when the point clouds are mathematically modelled with B-spline surfaces so that an adapted congruence test at the level of the estimated surface can be performed. However, test statistics for deformation derived from the surface differences between two epochs are known to be strongly influenced by the VCM of the errors chosen. In this contribution, we proposed a simplified stochastic model with data-driven parameters. We made use of the Matérn covariance function, whose parameters are obtained from the residuals of B-spline approximation. Using the EMD, we justified the use of a power-law model from the linear relationship between the mean energy of the IMF and its IMF index. In addition to being simple, the main advantage of this methodology is in its computational stability and the easy identification of model misspecification. By means of the de-biased Whittle likelihood, the model parameters, including the variance, can be trustworthily estimated, even for small samples.

This procedure allowed us to investigate the impact of correlations on the test statistics and the effect of neglecting them. As a follow-up to a previous contribution dealing with mathematical correlations, we used a real data set from a bridge under load. Two areas corresponding to different scanning geometries and magnitudes of deformation were approximated with B-spline surfaces using optimal fitting parameters. It could be shown that the test statistics obtained by neglecting correlations were overoptimistic: accounting for correlations acts as decreasing the variance in an equivalent diagonal matrix. As expected, the correlation structure was found to depend strongly on the scanning geometry, i.e. an optimized geometry with short range and optimal incidence angle decreased the level of correlations of the LS residuals and, thus, the impact of neglecting correlations in the test statistics for deformation. Correlations should not be neglected for a trustworthy test decision: they provide a more trustworthy test decision for deformation. A Gamma distribution could be fitted from the test statistics using a Monte Carlo approach. Since the critical values of the distribution—as well as the corresponding test statistics—were larger when correlations were neglected, the risk of making a wrong decision that a deformation occurs was decreased. However, the optimal test distribution should be investigated further. This will be the topic of our next contribution.

Appendices

Appendix 1

The Matérn process is a three-parameter continuous Gaussian process defined by an unaliased spectrum $W$:

$$ W\left( \omega \right) = \frac{{\alpha^{2\nu - 1} }}{c}\frac{{\sigma^{2} }}{{\left( {\omega^{2} + \alpha^{2} } \right)^{\nu } }} $$

(6)

where $\omega$ is the frequency. In this formulation, $\sigma^{2}$ is the variance of the process, $c$ a normalizing constant, $\alpha$ the range parameter and $\nu$ the smoothness parameter (Lilly et al. 2017). Thus, whereas the smoothness determines the behaviour of the spectral density as $\omega \to \infty$, the range parameter has a strong impact at low frequency as $\omega \to 0$. The corresponding autocovariance function is found to be (Mátern 1960):

$$ C\left( \tau \right) = \sigma^{2} {\rm M}_{\nu } \left( {\alpha \tau } \right) $$

(7)

where we introduce the Matérn function as

$$ {\rm M}_{\nu } \left( \tau \right) = \frac{2}{{\varGamma \left( {\nu - {\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}}} \right)2^{{\left( {\nu - {\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}}} \right)}} }}\left| \tau \right|^{{\left( {\nu - {\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}}} \right)}} K_{{\left( {\nu - {\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}}} \right)}} \left( {\left| \tau \right|} \right) $$

(8)

with $\varGamma$ the gamma function, $K_{{\left( {\alpha - {\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}}} \right)}}$ the decaying modified Bessel function of the second kind of order $\alpha - {\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}}$. $\tau > 0$ is the time difference, $\alpha$ a range parameter that controls how quickly C_matern decreases as $\tau$ increases, whereas $\nu$ is related to the mean-square differentiability of the process at the origin, as $\tau \to 0$.

Figure 8 (left) displays $C\left( \tau \right)$ scaled to 1 (i.e. the correlation function) for different choices of the shape parameter $\nu$. The smoothness parameter determines the behaviour of the covariance function as $\tau$ approaches 0. It has a great importance since it determines how fully populated the inverse of the VCM will be and how large its impact is.

Fig. 8

Matérn correlation function. Left: variation of the smoothness parameter $\nu$ with $\alpha = 0.05$. Right: variation of the range parameter $\alpha$ by keeping $\nu$ fixed to 2

The following cases are well-known for the one-dimensional field and should be mentioned:

• $\nu = {\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}}$ corresponds to the exponential covariance function, i.e. a short memory process

• $\nu = 1$ to the Markov process of first order, and

• $\nu = \infty$ gives the squared exponential covariance function, which corresponds to a physically less plausible infinitely differentiable Gaussian process at the origin.

Figure 8 (right) shows the different correlation functions obtained by varying the range parameter $\alpha$. As can be seen, $\alpha$ is linked to the decay of the covariance function. Please note that other parametrizations of the Matérn function are possible (Stein 1999)

Appendix 2

Huang et al. (1998) developed the EMD to decompose a given time series into a set of components of different time scales called the Intrinsic Mode Functions (IMF). Combined with the Hilbert Transform, this decomposition provides a time–frequency-energy representation of a signal and is intuitively based on the wave form of the data; it is also called the Hilbert-Huang transform in the literature. The EMD is independent of mother wavelets, such as in wavelet analysis or trigonometric functions, for example, in Fourier transform. Thus, it can be used on non-stationary signals and allows a characterization by direct extraction of the local energy associated with the IMF. The decomposition is, however, not unique and depends exemplarily on the particular interpolation scheme for envelope extraction, the stopping criterion (usually 1e-6) used in the sifting process and the number of IMF to determine.

The EMD algorithm can be summarized as the following loop (Flandrin et al. 2004):

• 1. Identification of the extrema of the signal

• 2. Interpolation between minima (maxima) to get an envelope of the signal

• 3. Computation of the average of the envelope

• 4. Extraction of the detail as difference between the envelope and the original signal

• 5. Iteration on the residual

The IMF should have the following properties: (i) zero mean, (ii) all the maxima and all the minima of h(i)(t) are correspondingly positive and negative (iii) narrowband. The number of zero-crossings can be either N(i) or N(i) − 1. Each IMF has fewer extrema than all the lower order IMFs, leading to fewer and fewer oscillations as the IMF order increases.

The EMD components of coloured noise (i.e. having a spectral decomposition following a power law) have the same rescaled Fourier spectrum. Thus, apart from the first noise-only IMF which carries the main energy and is considered as the white noise component, the energies of the other IMF decrease linearly in a semi-log diagram regarding the IMF index. We understand under “energy” the empirical variance estimate of a given IMF. Any discrepancies from the linear relationship can be used to denoise a signal, see, for example, Kopsinis and McLaughlin (2009). In this contribution, the noise itself is of interest as it is used to estimate the Matérn parameters: we make use of this important property to extract the noise from the other IMF, considered as model misspecifications.

Appendix 3

Assuming Gaussian data, the parameters of the Matérn model—including the variance—can be estimated with Maximum Likelihood, eventually by fixing one parameter to reduce the computational burden (Kaufman et al. 2008). This standard approach has a computational efficiency of the order of $n_{obs}^{2}$ for regularly spaced observations. It produces estimates that are asymptotically efficient, i.e. for smaller samples alternatives such as Restricted Maximum Likelihood have to be employed. Another powerful approach is referred as “quasi-likelihood”, as introduced in Sykulski et al. (2019) using the Whittle likelihood (Percival and Walden 1993). This method for estimating the Matérn parameters of any time series modelled as a second-order stationary stochastic process is employed in the freely available software package www.ulc.ac.uk/statistics/research/spg/software.

Exact maximum likelihood inference can be performed by maximizing the log-likelihood. This strategy, however, involves matrix inversion so that Whittle (1953) proposed approximating the autocovariance matrix in the frequency domain. The “Whittle likelihood” makes use of the periodogram of the observations, which may be a biased measure for finite samples (Percival and Walden 1993) due to blurring and aliasing effects. Sykulski et al. (2019) proposed an alternative called the “de-biased Whittle likelihood”, which is particularly interesting for discrete-time processes of finite length, such as the residuals of the LS adjustment used in this contribution.