A New Sequential Statistical Test Procedure from Observation Differences for Geodetic Deformation Analysis

1. Introduction

One of today’s major challenges of geodetic data analysis is to detect geometric changes of objects or areas which are subject to displacements and/or deformations, such as man-made structures like dams, dikes, bridges, wind turbines or high towers as well as natural Earth structures like volcanos, mining area, or tectonic plates. The analysis of monitoring measurement can be categorized into four deformation models: congruence model, kinematic model, static model and dynamic model. The congruence model describes the deformations by means of displacement vectors without specifying the time and any factor related to the acting forces, internal and external loads as well. The kinematic models describe the geometric changes in terms of temporal variations (velocities and accelerations), considering that the state of the object is permanently in motion, but also there are no concerns with the causes of the deformations. If there is interest in investigating the functional relationship between causative forces and geometrical reactions of the object, then the static model will be more suitable. In latter, the deformations are described from physical properties of the object (e.g., expansion coefficients, temperature, and lengths), so that the temporal aspects are not explicit in the model. Finally, the dynamic model is a combination of static and kinematic models, i.e., deformations are linked to their influencing factors (causative forces, internal and external loads) and the object’s physical properties [1]. Due to its common usage, we restrict ourselves to the congruence model which only tells us whether the object has moved or not.

In the congruence analysis, the structure under investigation is often monitored by a geodetic network which is measured in at least two epochs in time, and these epochal measurements are then analysed statistically. The geodetic network works as a displacement monitoring system. The statistical test is one of the most widely used approaches to the specification of deformation congruence models [2,3,4,5,6,7]. The robust approach is another one also widely used and has been had important advances in recent years [5,8,9,10], but it is not part of the scope of this work. Typically, the input data is a vector composed of the differences between the coordinates of points estimated by least squares at an initial epoch (say, Epoch I) and a current epoch (say, Epoch 2). The null hypothesis, denoted by ${\mathcal{H}}_0$, is formulated under the condition that all points are stable (points which do have a congruent/rigid geometrical structure at both considered epochs). On the other hand, the alternative hypotheses are stipulated from the assumption that there is at least one unstable point. As it is not known which point or group of points are unstable, a consecutive hypothesis test is often applied to identify one unstable point after the other [11]. Such a test procedure is similar to the outlier screening by iterative data-snooping [12]. However, the iterative consecutive hypothesis tests are non-rigorous because the alternative hypotheses are restricted to only one single unstable point [11]. The point which is flagged as most suspected to be unstable in a given step is not then inspected in the next iteration step. In addition, the result of a misidentification in a given iteration conditions the result of the next iteration [13]. The weakness of the iterative consecutive hypothesis tests for the case where multiple displacements are in play has been reported by several authors [13,14,15,16].

To overcome the problem of iterative procedures, non-iterative combinatorial procedure emerges for the case where all possible combinations of displaced points are considered [11,17,18]. Such procedure consists of comparing all possible candidates for stable points at the same stage, and consequently it is not necessary to consecutively point-by-point specify the congruence model. This method has been applied in some numerical examples and discussed in detail by Velsink [17,18]. Velsink proposed the ratio of the test statistic and its critical value as a decision rule. The point or group of points with the largest ratio is flagged as unstable, if it exceeds the critical value. Another interesting combinatorial method is discussed by Lehmann and Lösler [11]. They use various information criteria, and then select from among all possible candidate models the one which provide the lowest information criterion value as the best model. The idea of ​​using information criteria is outside the scope of the work but will be considered in future works.

Unfortunately, combinatorial-only method is often done from the set of different-dimensional models [13]. The comparison between models of different dimensions is complicated in this case. For example, the more points are modeled as unstable (higher dimension of congruence model – more complex is the model), the larger the occurrence of overfitting (a model is always better fitted to observations with a larger number of parameters). Model complexity can be circumvented by applying penalties. As highlighted by Nowel [13], the goodness of model fit and a penalty term constitutes an identification criterion. However, he warns there are many criteria, and it is not yet clear which of these to adopt. Nowel used the possibilities of combinatorics and generalized likelihood ratio tests performed in an iterative step to overcome the weaknesses related to the both iterative consecutive point-by-point model specification and combinatorial-only method. Although there has been substantial progress in this new combinatorics field, there are still challenges that open up new research perspectives [19,20,21,22,23,24,25,26]. In this contribution, we present an alternative and sophisticated method that integrates combinatorial analysis and likelihood ratio tests in a sequential procedure, which we call Sequential Likelihood Ratio Tests for Unstable Points Identification – SLRTUPI. The procedure is an extension of the sequential method for detecting multiple outliers proposed by Klein et al. [27].

Here, the method makes use of differences of the observations from two epochs instead of estimated coordinates, as proposed by some authors [9]. The idea of ​​using unadjusted (original) observation differences has been around for some time, as can be seen in [28,29]. When adopting the differences of observations as a vector of observations in the Gauss-Markov model, for example, we do not need to concern about the problem of defining the datum of the epochs and applications of the S-transformation [30,31,32,33]. The effect of the network geometry is eliminated [34]. On the other hand, we must guarantee that the campaigns are always carried out with the same occupation of the points in order to be able to compare measurements between epochs – epoch-wise measurements.

Erdogan et al. [32] presented a methodology for identifying unstable points based on this approach of analyzing the differences between measurements taken in two epochs, which they called the univariate approach. As a result, we do not also need to concern about whether it is a linear geodetic network (e.g., levelling or GNSS vector network) or non-linear (e.g., trilateration), since the univariate model is always linear. Linearization of a nonlinear model may reduce the detection power [35]. In addition, the univariate approach has the benefit of reducing the smearing and masking effect of displacements. Smearing means that one unstable point makes another stable point appear as an unstable and masking means that one unstable point prevents another one from being identified.

However, Erdogan et al. [32] do not control for false positive rates (Type I Error rate – detect displacements when in reality there are none). This is because they generally consider that the test involves only a single alternative hypothesis, when in fact they are multiple tests, i.e., it involves multiple alternative hypotheses. Consequently, an approach that allows controlling the Type I Error is needed. Here, for example, we use a Monte Carlo method so that the user-defined Type I Error rate is efficiently controlled [36]. Furthermore, the works cited above do not make it clear how to choose the maximum number of points to be modeled as unstable (displaced). This choice is somewhat arbitrary. To avoid a subjective choice, the maximum number of possible points to be inspected as displaced (unstable) is based on the rank computation of the design matrices constructed for all possible combinations of points modeled as unstable, as well as on the statistical overlap analysis (i.e., when it is not possible to distinguish one model from another, since the computed test statistics present the same values, and therefore identification cannot occur).

The next part of the paper is organized as follows. Section 2 describes the univariate model under null and alternative hypothesis for the case where the point or group of points to be tested is a priori specified. In this section, a trilateration network is presented, which is the object of study throughout the entire paper. In addition, we present the matrix that makes the connection between the network points and the measurements and its conversion to the displacement-design matrix that captures the sign of the differences between the observations from the two epochs in time. This section ends with the test statistic derived from the maximum likelihood ratio test concept for the case where there is only one single alternative hypothesis. In section 3, we describe the proposed SLTRUPI method step-by-step. We provide a Monte Carlo-based procedure for controlling the false positive rates in subsection 3.1. And in the last subsection 3.2., we present the procedure to determine the maximum number of points allowed for identification via SLRTUPI, and its application in some numerical examples. Section 4 is devoted to experiments and results based on computer simulations and real dataset to demonstrate the reliability and performance of SLRTUPI in several displacement pattern scenarios. In that section, we also describe the success and failure rates classes associated with SLRTUPI for the case of having both individually and mutually (simultaneously) unstable points. Section 6 highlights the contributions from the study.

2. Null and Alternative Hypothesis

Let's start by describing the univariate model for identifying unstable points in geodetic networks by the following equations [32,34]:

$$\begin{array}{c}{\mathit{y}}^{\left\{1\right\}}={{\mathit{y}}_{\mathbf{0}}}^{\left\{1\right\}}+{\mathit{e}}^{\left\{1\right\}} \mathrm{for} \mathrm{Epoch} 1\\ {\mathit{y}}^{\left\{2\right\}}={{\mathit{y}}_{\mathbf{0}}}^{\left\{2\right\}}+{\mathit{e}}^{\left\{2\right\}} \mathrm{for} \mathrm{Epoch} 2\end{array}$$

(1)

where${\mathit{y}}^{\left\{.\right\}}\in {ℝ}^{n\times 1}$ are the vectors of geodetic measurements, ${{\mathit{y}}_{\mathbf{0}}}^{\left\{.\right\}}\in {ℝ}^{n\times 1}$ the unknowable true quantity vectors of measurand, and ${\mathit{e}}^{\left\{.\right\}}\in {ℝ}^{n\times 1}$ the unknown vectors of measurement errors (note: the indices ‘1’ and ‘2’ inside the curly braces represent the quantities related to the first and second epoch in time, respectively).

By subtracting the Epoch 2 equation from the Epoch 1 equation, the univariate model can be formulated as a linear Gauss-Markov model, as follows:

$$\begin{array}{c}{\mathit{y}}^{\left\{2\right\}}-{\mathit{y}}^{\left\{1\right\}}={{\mathit{y}}_{\mathbf{0}}}^{\left\{2\right\}}-{{\mathit{y}}_{\mathbf{0}}}^{\left\{1\right\}}+{\mathit{e}}^{\left\{2\right\}}-{\mathit{e}}^{\left\{1\right\}}, \mathrm{making}\\ \begin{array}{c}{\Delta }_{\mathit{y}}={\mathit{y}}^{\left\{2\right\}}-{\mathit{y}}^{\left\{1\right\}},\\ \begin{array}{c}{\mathit{e}}_{{\Delta }_{\mathit{y}}}={\mathit{e}}^{\left\{2\right\}}-{\mathit{e}}^{\left\{1\right\}},\\ \begin{array}{c}x={{\mathit{y}}_{\mathbf{0}}}^{\left\{2\right\}}-{{\mathit{y}}_{\mathbf{0}}}^{\left\{1\right\}}\\ \therefore {\Delta }_{\mathit{y}}=Ax+{\mathit{e}}_{{\Delta }_{\mathit{y}}}\end{array}\end{array}\end{array}\end{array}$$

(2)

with ${\Delta }_{\mathit{y}}\in {ℝ}^{n\times 1}$ being the vector of the two-epoch geodetic observations differences, ${\mathit{e}}_{{\Delta }_{\mathit{y}}}\in {ℝ}^{n\times 1}$ the unknown vector of errors of the two-epoch geodetic observations differences, $\mathit{A}\in {ℝ}^{n\times 1}$ the Jacobian matrix (also called the design matrix) of full rank $u=1$, which in this case is a column vector of ones (i.e., $\mathit{A}={\left[\mathbf{1} \mathbf{1} \mathbf{1} ... \mathbf{1}\right]}^T$) , and $x\in {ℝ}^{1\times 1}$ the estimable unknown parameter, which in this case is a scalar that represents the unknown true difference between the two epochs. It is important to highlight that now our vector of observations corresponds to the differences of the geodetic measurements at two epochs in time. Furthermore, these differences are computed from the identical observations of each epoch.

After having defined the univariate model according to Eq. (2), we now need to resort to hypothesis testing theory to infer whether a subset of point fields measured in two epochs is stable/congruent or not. Typically, in displacement detection, the model under null hypothesis, denoted by ${\mathcal{H}}_0$, is set up for the case where all points to be analyzed are treated as stable points. In other words, the null hypothesis states that the null model in Eq. (2) explains the observations. Assuming normally distributed observation errors with expectation zero, i.e.:

$${\mathit{e}}_{{\Delta }_{\mathit{y}}({\mathscr{H}}_0)} ~N\left(\mathbf{0},{\Sigma }_{{\Delta }_{\mathit{y}}}\right)$$

(3)

the null hypothesis ${\mathcal{H}}_0$ of the standard Gauss–Markov model in the linear form of the Eq. (2) is then given by:

(4)

with $E\left( . \right)$ the expectation operator, $D\left( . \right)$ the dispersion operator, and ${\Sigma }_{{\Delta }_{\mathit{y}}}\in {ℝ}^{n\times n}$ the positive-definite variance matrix of ${\Delta }_{\mathit{y}}$. The variance matrix ${\Sigma }_{{\Delta }_{\mathit{y}}}$ is ​​obtained by applying the variance propagation law to the ${\Delta }_{\mathit{y}}$, which is the result of the sum of the variance matrices of the geodetic observations from Epoch 1 $\left({\Sigma }_{\mathit{y}}^{\left\{1\right\}}\in {ℝ}^{n\times n}\right)$ and Epoch 2 $\left({\Sigma }_{\mathit{y}}^{\left\{2\right\}}\in {ℝ}^{n\times n}\right)$.

When the null model in Eq. (4) is assumed to be true, the scalar parameter $x$ can be estimated by simple least-squares calculus as:

$${\widehat{x}}_{({\mathscr{H}}_0)} ={\left({\mathit{A}}^{\mathrm{T}}\mathit{W}\mathit{A}\right)}^{-1}{\mathit{A}}^{\mathrm{T}}\mathit{W}{\Delta }_{\mathit{y}}$$

(5)

and the estimated observation errors as being:

$${\widehat{\mathit{e}}}_{{\Delta }_{\mathit{y}}({\mathscr{H}}_0)} =\mathit{A}\widehat{x}-{\Delta }_{\mathit{y}}$$

(6)

where $\mathit{W}\in {ℝ}^{n\times n}$ is the known matrix of weights, taken as $\mathit{W}={\mathsf{\sigma }}_0^2{\Sigma }_{{\Delta }_{\mathit{y}}}^{-1}$, where ${\mathsf{\sigma }}_0^2$ is the priori variance factor. The overall degrees of freedom $r$ (redundancy) of the model under ${\mathscr{H}}_0$ is given by:

(7)

$${\Sigma }_{{\widehat{\mathit{e}}}_{{\Delta }_{\mathit{y}}({\mathcal{H}}_0)}}={\mathsf{\sigma }}_0^2{\mathit{W}}^{-1}-{\mathsf{\sigma }}_0^2\mathit{A}{\left({\mathit{A}}^{\mathrm{T}}\mathit{W}\mathit{A}\right)}^{-1}{\mathit{A}}^{\mathrm{T}}$$

(8)

where ${\Sigma }_{{\widehat{\mathit{e}}}_{{\Delta }_{\mathit{y}}({\mathcal{H}}_0)}}\in {ℝ}^{n\times n}$ is the estimated variance matrix of the observation errors.

On the other hand, an alternative model can be proposed when there are doubts about the stability of points in the network. Here, we restrict ourselves – for simplicity of the analyses – measurements not contaminated by outliers. Thus, for the case of univariate approach, the model under the alternative hypothesis, denoted by ${\mathcal{H}}_{\mathrm{A}}$, is to oppose Eq. (4) by an extended model that includes an extra parameter $\nabla \in {ℝ}^{p\times 1}$ which describes the disturbance on the measurements in function of the displacement of a subset “$p$” of network points, as follows:

(9)

with $\mathit{G}\in {ℝ}^{n\times p}$ being the matrix which captures the relationship between the displaced points and the changes in the measurements connected to them (hereafter, this matrix will be called displacement-design matrix). If the alternative model in Eq. (9) holds true, then the unknown parameters will be obtained as:

$${\left(\begin{array}{c}\widehat{x}\\ \widehat{\mathbf{\nabla }}\end{array}\right)}_{({\mathscr{H}}_{\mathrm{A}})} ={\left[{\left(\begin{array}{cc}\mathit{A}& \mathit{G}\end{array}\right)}^{\mathrm{T}}\mathit{W}\left(\begin{array}{cc}\mathit{A}& \mathit{G}\end{array}\right)\right]}^{-\mathbf{1}}{\left(\begin{array}{cc}\mathit{A}& \mathit{G}\end{array}\right)}^{\mathrm{T}}\mathit{W}{\Delta }_{\mathit{y}}$$

(10)

The redundancy of the model under ${\mathcal{H}}_{\mathrm{A}}$ is $n-rank\left(\begin{array}{cc}\mathit{A}& \mathit{G}\end{array}\right)$, with the estimated observation errors ${\widehat{\mathit{e}}}_{{\Delta }_{\mathit{y}}({\mathcal{H}}_{\mathrm{A}})}$and estimated variance matrix of the observation errors ${\Sigma }_{{\widehat{\mathit{e}}}_{{\Delta }_{\mathit{y}}\left({\mathcal{H}}_{\mathrm{A}}\right)}}$ under alternative hypothesis ${\mathcal{H}}_{\mathrm{A}}$ given respectively by:

$${\widehat{\mathit{e}}}_{{\Delta }_{\mathit{y}}({\mathscr{H}}_{\mathrm{A}})} =\left(\begin{array}{cc}\mathit{A}& \mathit{G}\end{array}\right)\left(\begin{array}{c}\widehat{x}\\ \widehat{\nabla }\end{array}\right)-{\Delta }_{\mathit{y}}$$

(11)

$${\Sigma }_{{\widehat{\mathit{e}}}_{{\Delta }_{\mathit{y}}\left({\mathscr{H}}_{\mathrm{A}}\right)}}={\mathsf{\sigma }}_0^2{\mathit{W}}^{-1}-{\mathsf{\sigma }}_0^2\left(\begin{array}{cc}\mathit{A}& \mathit{G}\end{array}\right){\left[{\left(\begin{array}{cc}\mathit{A}& \mathit{G}\end{array}\right)}^{\mathrm{T}}\mathit{W}\left(\begin{array}{cc}\mathit{A}& \mathit{G}\end{array}\right)\right]}^{-\mathbf{1}}{\left(\begin{array}{cc}\mathit{A}& \mathit{G}\end{array}\right)}^{\mathrm{T}}$$

(12)

As a simple example, an alternative hypothesis may then be formulated for the case of having a specific point as displaced ($p=1$), as can be seen in the illustration of a trilateration network in Figure 1. It is presumed that only point F has been moved to a new position F’. Consequently, this causes changes only to the observations connected to it. This means that such distances measured in the second epoch undergo deformations in the sense of stretching and/or compressing in relation to their initial states, whereas the other ones remain stable. The initial state refers to the first epoch the measurements were gathered, so that all network points are assumed to be free of displacements. For point F, for example, the $\left(\begin{array}{cc}\mathit{A}& \mathit{G}\end{array}\right)$ matrix is described as:

$${\left(\begin{array}{cc}\mathit{A}& \mathit{G}\end{array}\right)}^{\mathrm{T}}=\left( \begin{array}{c} 1\\ 0\end{array}\begin{array}{c}1\\ 0\end{array}\begin{array}{c}1\\ -1\end{array}\begin{array}{c} 1\\ 0\end{array}\begin{array}{c} 1\\ 0\end{array}\begin{array}{c} 1\\ -1\end{array}\begin{array}{c}1\\ 0\end{array}\begin{array}{c}1\\ 0\end{array}\begin{array}{c}1\\ 1\end{array} \right)$$

(13)

The sign assigned to the elements of displacement-design matrix $\mathit{G}$ depends on the sign of the two-epoch geodetic observations differences. A positive sign means that at epoch 2 the measure was positively distorted and therefore the parameter acts positively ($+\nabla$). The same occurs for the case in which the measure is reduced in relation to its initial state and, therefore, the parameter acts with a negative sign ($-\nabla$). In this example (Figure 1), the shifted point F to its new position F’ makes the AF' segment smaller than its initial state AF. Consequently, the disturbance parameter on that measurement in function of the displacement of the F point acts with negative sign ($-\nabla$). The same occurs for BF, i.e., BF' segment gets smaller than the BF segment ($-\nabla$), while for the case of CF sign is positive, since CF' segment gets larger than CF ($+\nabla$).

The displacement-design matrix $\mathit{G}$ in expression (10) has been constructed for the case of assuming only one single displaced point, which in our example was point F. In that case, the unknown parameter vector $\mathbf{\nabla }$ has become a scalar $\nabla$, and the displacement-design matrix $\mathit{G}$ has reduced to a unit vector, whose elements are exclusively formed by values of 0, 1 and -1, where 1 or -1 means that the ith parameter of magnitude $+\nabla$ or $-\nabla$, respectively, affects the measurements connected to the point, and 0 means otherwise. However, as given by Eq. (9), we can design an overall displacement-design matrix $\mathit{G}$, so that each column represents the displacement of a given point. For this, we first provide an auxiliary matrix $\mathit{C}\in {ℝ}^{n\times p}$ which describes the general case of the relationship between the measures (matrix rows) and the network points (matrix columns), but without taking into account the sign of the unknown estimable vector $\mathbf{\nabla }$. We call this matrix of Point-to-Measurement Connection Matrix, which for the network in Figure (1) is given by:

(14)

Then, the signs of the coefficients of the displacement-design matrix $\mathit{G}$ can be obtained generically from the application of a signum function as:

(15)

where $sgn\left(.\right)$ is the sign function and $\left|\Delta y_{\left(i\right)}\right|$ gives a non-negativity value (absolute value) for each two-epoch geodetic observation differences. The result of Eq. (15) can be stacked into a vector as $\Delta {\mathit{y}}_{\mathit{s}\mathit{g}\mathit{n}}=sgn\left(\mathsf{\Delta }{\mathit{y}}_{\left(l\right)}\right)$. Thus, the displacement-design matrix for the general case becomes:

$${\mathit{G}}_{\mathit{p}}=diag\left(\Delta {\mathit{y}}_{\mathit{s}\mathit{g}\mathit{n}}\right)\mathit{C}$$

(16)

where $diag\left(.\right)$ is the diagonalize operator which convert a vector to a diagonal matrix. Consequently, the alternative hypothesis is now given by:

(17)

The displacement-design matrix ${\mathit{G}}_{\mathit{p}}$ can be developed in several ways, depending on which point, or group of points will be chosen to be monitored. Taking another example in Figure 2, if $n=10$ and $p=2$, then a possible displacement-design matrix ${\mathit{G}}_{\mathit{p}}$ would be ${\mathit{G}}_{\mathit{p}=\mathbf{2}}={\left(\begin{array}{cc}1& 1\\ 0& 0\end{array}\begin{array}{cc}1& 0\\ 0& -1\end{array}\begin{array}{cc}0& 0\\ -1& -1\end{array}\begin{array}{cc}0& 0\\ 0& 0\end{array}\begin{array}{cc}0& 0\\ 0& 0\end{array}\right)}^{\mathrm{T}}$and ${\mathbf{\nabla }}_{\mathit{p}=\mathbf{2}}=\left(\begin{array}{c}{\nabla }_A\\ {\nabla }_B\end{array}\right)$ for the case where the point A and B were considered simultaneously as displaced points. It is important to note that in this example the first and second columns of the matrix $\mathit{C}$ were selected, which refer to points A and B, respectively.

After postulating the null and alternative hypotheses, the test statistic derived from the concept of the maximum likelihood ratio test are computed as [37]:

$$T_p={\widehat{\mathit{e}}}_{{\Delta }_{\mathit{y}}}^{\mathrm{T}} \mathit{W} {\mathit{G}}_{\mathit{p}}{\left({{\mathit{G}}_{\mathit{p}}}^{\mathrm{T}} \mathit{W} {\Sigma }_{{\widehat{\mathit{e}}}_{{\Delta }_{\mathit{y}}}}\mathit{W} {\mathit{G}}_{\mathit{p}}\right)}^{-1}{{\mathit{G}}_{\mathit{p}}}^{\mathrm{T}} \mathit{W} {\widehat{\mathit{e}}}_{{\Delta }_{\mathit{y}}}$$

(18)

Then, a test decision is then performed as

(19)

Under the null hypothesis ${\mathcal{H}}_0$, the test statistic $T_p$ follows the central chi-squared distribution with $p$ degrees of freedom; under the alternative hypothesis ${\mathcal{H}}_{\mathrm{A}}$, the test statistic $T_p$ follows the non-central chi-squared distribution with $p$ degrees of freedom and non-centrality parameter ${\lambda }_0$. If the test statistic $T_p$ is larger than a critical value $k$ ($T_p>c_0$), then a particular combination of tested points is flagged as displaced, where $c_0$ is the chi-squared critical value at a given significance level ${\mathsf{\alpha }}_0$ with $p$ degrees of freedom, i.e., $c_0={\mathsf{\chi }}_{\left(p,{\mathsf{\alpha }}_0\right)}^2$. The probability level ${\mathsf{\alpha }}_0$ defines the size of a test and is often called the “false alarm probability”. Thus, the null model is assessed by specifying a significance level ${\mathsf{\alpha }}_0$. [Note: the subscript "0" means that the test is performed for the case of having only one single alternative hypothesis. Thus, the non-centrality parameter ${\lambda }_0$ and the significance level ${\mathsf{\alpha }}_0$ with the subscript "0" means that it is a single test involving only one single alternative hypothesis. In that case, rejecting the null hypothesis automatically consists of accepting the alternative hypothesis, and therefore, detection is the same as identification.]

In addition, the hypothesis test-based approach presented here is based on the mean-shift model: the alternative model is formulated under the condition that the displacement of the point (or a subset of network points) acts as a systematic effect by shifting the random error distribution under ${\mathcal{H}}_0$ by its own value. This means that a displacement can cause a shift of the expectation under ${\mathcal{H}}_0$ to a nonzero value. Therefore, hypothesis testing is often employed to check whether the possible shifting of the random error distribution under ${\mathcal{H}}_0$ is, in fact, a systematic effect (bias) coming from a possible deformation or merely a random effect due to measurement uncertainties [38].

The alternative model in Eq. (17) has been formulated under the assumption that a specific group of points has undergone displacement. This group of points can be formulated for the cases where we have a single only one point displaced (Figure 1) or a subset of such points of the geodetic network (Figure 2). In other words, the number of points to be monitored and their locations on the network have been fully specified. This is a special case of testing the null hypothesis against only one single alternative hypothesis, and therefore, the rejection of the null model automatically implies the acceptance of the alternative model and vice versa.

3. Sequential Likelihood Ratio Tests for Unstable Points Identification – SLRTUPI

One of the questions that arises is how to define the number of points to be monitored and their locations. Clearly, geodesists need to work in an interdisciplinary way with other professional/areas so that an alternative model may be properly formulated [18]. However, even with such information, we may still have doubts whether there are points that a priori were assumed to be stable, whereas in fact they may be shifting. From a practical point of view, therefore, we often do not know well how many and which points are or are not stable. Here, all geodetic network points are subjected to be inspected, regardless of whether they are located on structure or not. Thus, a more conservative alternative hypothesis would be to assume that there is at least one point that is moving among all points involved in the geodetic network. Consequently, we would have $n_p$ alternative models for $p=1$ to be tested against the null hypothesis, i.e.:

(20)

Note that ${\mathit{G}}_{\mathit{p}=\mathbf{1}}^{\left(g_i\right)}$ represents the ith column vector of the displacement-design matrix ${\mathit{G}}_{\mathit{p}=\mathbf{1}}$. For the case of the network illustrated in Figure 1 and Figure 2, we would then have 6 alternative hypotheses, since we have a total of $n_p=6$ points, so that $g_1:=\left[A\right];g_2:=\left[B\right];g_3:=\left[C\right];g_4:=\left[D\right];g_5:=\left[E\right];g_6:=\left[F\right].$ This means that we would have to test ${\mathcal{H}}_0$ against ${\mathcal{H}}_{\mathrm{A}}^{\left(A\right)},{\mathcal{H}}_{\mathrm{A}}^{\left(B\right)},{\mathcal{H}}_{\mathrm{A}}^{\left(C\right)},{\mathcal{H}}_{\mathrm{A}}^{\left(D\right)},{\mathcal{H}}_{\mathrm{A}}^{\left(E\right)},{\mathcal{H}}_{\mathrm{A}}^{\left(F\right)}$, i.e.:

(21)

We are now dealing with testing involving multiple hypotheses. Now, we are interested in knowing which of the alternative hypotheses may lead to the rejection of the null hypothesis with a certain probability. To answer this, the test statistic in Eq. (18) is now computed for each of such tests $\left(T_{p=1}^{\left(g_i\right)}\right)$, so that we would have a vector of test statistics, denoted by ${\mathit{T}}_{\mathit{p}=\mathbf{1}}$. In that case, the test statistic coming into effect is the maximum test value, which is computed as

(22)

The decision rule for this case is given by:

$$\mathrm{Accept} {\mathscr{H}}_0 \mathrm{if} \max \left({\mathit{T}}_{\mathit{p}=\mathbf{1}}\right)\le c, \mathrm{Reject} \mathrm{otherwise} \mathrm{in} \mathrm{favour} \mathrm{of} {\mathscr{H}}_{\mathrm{A}}^{\left({\mathrm{g}}_i\right)}$$

(23)

The decision rule in (23) states that if none of the $n_p$ tests get rejected, then we accept the null hypothesis ${\mathcal{H}}_0$. If the $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{1}}\right)$ is larger than some percentile of its probability distribution (i.e., some critical value $c$), then there is evidence that there is a displaced point in the structure. In this case, we can only assume that detection occurred. The identification, however, is not straightforward. The identification (or localization) of the displaced point consists of seeking the point that produced the maximum test statistic $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{1}}\right)$, and whose value is greater than some critical value $c$. Thus, point identification only happens when point detection necessarily exists, i.e., “point identification” only occurs when the null hypothesis ${\mathcal{H}}_0$ is rejected. This means that the correct detection does not necessarily imply correct identification. For instance, we may correctly detect the occurrence of deformation, but wrongly identify a point as displaced, while in fact it is another point that has changed. The description of the success and failure rates will be covered in the next sections.

It is important to highlight that the maximum test statistic, now generally denoted by $\max \left({\mathit{T}}_{\mathit{p}}\right)$, is treated directly as a test statistic. Thus, the distribution of $\max \left({\mathit{T}}_{\mathit{p}}\right)$ cannot be derived from well-known test distributions (e.g., Chi-squared distribution). Therefore, critical values cannot be taken from a statistical table but must be computed numerically. Here, the critical value $c$ is computed by Monte Carlo such that a user-defined Type I decision error rate “${\alpha }_t$” for the proposed procedure is warranted [39]. In the next section, we provide details on how to compute the critical value $c$ for the proposed sequential procedure.

The procedure described so far allows only one single point to be identified. However, we may have multiple points displaced simultaneously. A more appropriate approach would be to apply a sequential test to decide the number and location of shifted points for the case where the null model is rejected. If the maximum test statistic value exceeds the critical value at the significance level adopted, i.e., if $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{1}}\right)>c$, the test statistics proceeds for $p=2$ (two possible unstable points). Since we do not know which group of two points might be shifting, we have to compute the test statistics by Eq. (18) from the possible combinations for two unstable points $\left({\mathit{T}}_{\mathit{p}=\mathbf{2}}\right)$, and then identifying the corresponding candidate group for $p=2$ through the maximum value of the test statistic $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{2}}\right)$. Note that the extreme test statistic $\max \left({\mathit{T}}_{\mathit{p}}\right)$ returns only one group of points or point. In general, the number of possible groups of points, denoted by $K_{n_p}^p$, is given by:

$$K_{n_p}^p=\left(\begin{array}{c}{n}_p\\ p\end{array}\right)=\frac{n_p!}{\left(n_p-p\right)!p!}$$

(24)

For instance, for a geodetic network with $n_p=6$ points and by assuming $p=2$ two-point group to be tested, we would have a total of $K_6^2=15$ possible combinations. That would be the case for the geodetic network presented here, as can be seen either by Figure 1 or Figure 2, then we would have $15$ alternative models, i.e., $\left({\mathcal{H}}_{\mathrm{A}}^{\left(A,B\right)},{\mathcal{H}}_{\mathrm{A}}^{\left(A,C\right)},{\mathcal{H}}_{\mathrm{A}}^{\left(A,D\right)},{\mathcal{H}}_{\mathrm{A}}^{\left(A,E\right)},{\mathcal{H}}_{\mathrm{A}}^{\left(A,F\right)},{\mathcal{H}}_{\mathrm{A}}^{\left(B,C\right)},{\mathcal{H}}_{\mathrm{A}}^{\left(B,D\right)},{\mathcal{H}}_{\mathrm{A}}^{\left(B,E\right)},{\mathcal{H}}_{\mathrm{A}}^{\left(B,F\right)},{\mathcal{H}}_{\mathrm{A}}^{\left(C,D\right)},{\mathcal{H}}_{\mathrm{A}}^{\left(C,E\right)},{\mathcal{H}}_{\mathrm{A}}^{\left(C,F\right)},{\mathcal{H}}_{\mathrm{A}}^{\left(D,E\right)},{\mathcal{H}}_{\mathrm{A}}^{\left(D,F\right)},{\mathcal{H}}_{\mathrm{A}}^{\left(E,F\right)}\right)$. Note that for $p=1$ it is a particular case of having $K_{n_p}^{p=1}=n_p$, as previously presented. Thus, the general cases of the alternative hypotheses can be expressed as:

(25)

If the identified point with $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{1}}\right)$ is also contained in the case of $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{2}}\right)$, i.e. if the one point flagged initially is among one of the pairs flagged subsequently, then the null hypothesis becomes specified according to the model identified at $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{1}}\right)$ and the alternative hypothesis to the model defined at $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{2}}\right)$, as can be seen the example below:

(26)

For example, if $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{1}}\right)>c$ and the point A in the geodetic network illustrated in Figure 1 (or Figure 2) was identified by the $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{1}}\right)$, the model under the null hypothesis would become ${\mathcal{H}}_0^{\left(A\right)}:E\left({\Delta }_{\mathit{y}}^{\left(A\right)}\right)=\mathit{A}x+{\mathit{G}}_{\mathit{p}=\mathbf{1}}^{\left(A\right)}{\nabla }_{p=1}^{\left(A\right)}$. If in the next step for ${\mathit{T}}_{\mathit{p}=\mathbf{2}}$, points A and C were identified by the $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{2}}\right)$ , then the alternative hypothesis would be defined as ${\mathcal{H}}_{\mathrm{A}}^{\left(A,C\right)}:E\left({\Delta }_{\mathit{y}}^{\left(A,C\right)}\right)=\mathit{A}x+{\mathit{G}}_{\mathit{p}=\mathbf{2}}^{\left(A,C\right)}{\mathbf{\nabla }}_{\mathit{p}=\mathbf{2}}^{\left(A,C\right)}$. Note that ${\mathcal{H}}_0^{\left(A\right)}$ is a subset of ${\mathcal{H}}_{\mathrm{A}}^{\left(A,C\right)}$. To decide which of these models to select, we compute the test statistic based on the maximum likelihood ratio between ${\mathcal{H}}_0$ and ${\mathcal{H}}_{\mathrm{A}}$, denoted by ${\mathsf{\Lambda }}_{\mathrm{MLR}}\left({\Delta }_{\mathit{y}}\right)$ [37]:

$${\mathsf{\Lambda }}_{\mathrm{MLR}}\left({\Delta }_{\mathit{y}}\right)={\widehat{\mathit{e}}}_{{\Delta }_{\mathit{y}}({\mathscr{H}}_0)}^{\mathrm{T}}\mathit{W}{\widehat{\mathit{e}}}_{{\Delta }_{\mathit{y}}({\mathscr{H}}_0)}-{\widehat{\mathit{e}}}_{{\Delta }_{\mathit{y}}({\mathscr{H}}_{\mathrm{A}})}^{\mathrm{T}}\mathit{W}{\widehat{\mathit{e}}}_{{\Delta }_{\mathit{y}}({\mathscr{H}}_{\mathrm{A}})}$$

(27)

The result of Eq. (27) comes from the ratio between the maximum of the probability density function of ${\Delta }_{\mathit{y}}$ under ${\mathcal{H}}_0$ and ${\mathcal{H}}_{\mathrm{A}}$, and the fact that the null hypothesis is (and must be) a subset of the alternative hypothesis. Thus, the selection of the most likely hypothesis is based on the following decision:

$$\mathrm{Accept} {\mathscr{H}}_0 {\mathrm{if} \mathsf{\Lambda }}_{\mathrm{MLR}}\left({\Delta }_{\mathit{y}}\right)\le c, \mathrm{Reject} \mathrm{otherwise} \mathrm{in} \mathrm{favour} \mathrm{of} {\mathscr{H}}_{\mathrm{A}}$$

(28)

where ${\widehat{\mathit{e}}}_{{\mathcal{H}}_0}$ and ${\widehat{\mathit{e}}}_{{\mathcal{H}}_A}$ are the least-squares estimated observation errors for the model under the null and alternative hypothesis, respectively. For the example above, we would have ${\widehat{\mathit{e}}}_{{\mathcal{H}}_0} =\mathit{A}\widehat{x}+{\mathit{G}}_{\mathit{p}=\mathbf{1}}^{\left(A\right)}{\widehat{\mathbf{\nabla }}}_{p=1}^{\left(A\right)}-{\Delta }_{\mathit{y}}$ for the null model and ${\widehat{\mathit{e}}}_{{\mathcal{H}}_A} =\mathit{A}\widehat{x}+{\mathit{G}}_{\mathit{p}=\mathbf{2}}^{\left(A,C\right)}{\widehat{\mathbf{\nabla }}}_{\mathit{p}=\mathbf{2}}^{\left(A,C\right)}-{\Delta }_{\mathit{y}}$ for the alternative one.

If the null hypothesis is not rejected, the testing ends and only the point corresponding to $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{1}}\right)$ is flagged as unstable. On the other hand, if the null hypothesis is rejected again, then a new step is started by taking the model under the alternative hypothesis identified in the previous step as the null hypothesis$\left(\mathrm{e}.\mathrm{g}., {\mathcal{H}}_0\equiv {\mathcal{H}}_{\mathrm{A}}^{\left(g_i\right)}\right)$. It proceeds by computing the test statistics for $p=3$ according to Eq. (18) for all possible combinations given by Eq. (24), identifying the corresponding candidate group for $p=3$ through the maximum value of the test statistic $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{3}}\right)$, and then checking if one of the points flagged in the previous step (i.e., for $p=2$) is among those flagged in the current step for $p=3$. If this is the case, then the test is applied through Eq. (26) to decide between the null model defined by the identified model for $p=2$ and the alternative model for $p=3$. This sequential procedure is repeated until any of the following conditions is met:

(i)

The current null hypothesis fails to be rejected.

(ii)

More than one group of geodetic network points is identified by the extreme statistic $\max \left({\mathit{T}}_{\mathit{p}}\right)$ for a given $p$, i.e., for the case where the hypotheses cannot be separable (statistical overlap), and therefore the identification cannot be accomplished.

(iii)

If the group of point(s) flagged by $\mathit{max}\left({\mathit{T}}_{\mathit{p}-\mathbf{1}}\right)$ is not fully contained in the group of points flagged by $\mathit{max}\left({\mathit{T}}_{\mathit{p}}\right)$, then the procedure ends with the null hypothesis given by the group of point(s) flagged by $\mathit{max}\left({\mathit{T}}_{\mathit{p}-\mathbf{1}}\right)$.

(iv)

Iteration reaches the threshold $p_{max}$ (i.e., until the maximum number of points to be evaluated is fully inspected). The definition of the maximum number of points will be detailed in the next section.

In general, therefore, the sequential testing procedure proposed here is based on likelihood ratio, which we now call Sequential Likelihood Ratio Tests for Unstable Points Identification (SLRTUPI). It consists in determining whether the additional subset of displaced point, considered in every new alternative hypothesis, is statistically significant or not, in terms of its impact on the quadratic form of the estimated observation errors, similar to its form for outlier detection [27]. Consequently, we can identify the number and the location of the displaced geodetic network points. Figure 3 provides step-by-step how to run the SLRTUPI procedure.

3.1. Monte Carlo Approach for Controlling the False Detection Rate

There is the probability of committing at least one false discovery, or Type I error when performing multiple hypotheses tests. Thus, Type I Error rate for the case where SLRTUPI is in play corresponds to the probability of incorrectly detecting at least one point as displaced while in fact there is none (i.e., accept at least one alternative hypothesis, when, in fact, the null hypothesis is true). This means that the SLRTUPI Type I decision error does not depend on all subsequent test steps, only on the extreme test statistic computed in its first step, i.e., only for $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{1}}\right)$ in Eq. (22). The risk of rejecting a true ${\mathcal{H}}_0$ is now one-fold: the undesired random event “reject a true ${\mathcal{H}}_0$” can occur in any of the $K_{n_p}^{p=1}=n_p$ tests. Let the probability of rejecting a true ${\mathcal{H}}_0$ in test $i$ be ${\alpha }_i$ (the so-called “experimentwise error rate”) and let ${\alpha }_i\ll 1$. Furthermore, assume the random events “reject a true ${\mathcal{H}}_0$ in test $i$” to be approximately statistically independent. Then the total probability of rejecting a true ${\mathcal{H}}_0$ in the multiple hypothesis test (the so-called “familywise error rate”, denoted here by ${\alpha }_t$) is [11]:

(29)

The classical and well-known procedure to control the ${\alpha }_t$ is to apply the Bonferroni equation by choosing [40]:

$${\alpha }_i:={\alpha }_t/n_p$$

(30)

Unfortunately, the test statistics $T_{p=1}^{\left(g_i\right)}$ and consequently the random events “reject a true ${\mathcal{H}}_0$ in test $i$” are statistically dependent. The extreme statistic $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{1}}\right)$ captures such dependencies, as it is extracted from the $T_{p=1}^{\left(g_i\right)}$. If such dependencies are neglected, then the computed critical values are erroneous, and the test decisions do not have the user-defined familywise error rate ${\alpha }_t$ [34]. Here, the maximum test value $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{1}}\right)$ in Equation (22) is treated directly as a test statistic [39]. Note that when using Equation (22) as a test statistic, the decision rule is based on a one-sided test of the form $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{1}}\right)\le c$, as can be seen in expression (23). However, the distributions of $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{1}}\right)$ cannot be derived from well-known test distributions (e.g., ${\chi }^2$-distribution). Therefore, critical values cannot be taken from a statistical table but should be computed numerically. A rigorous computation of critical values requires a Monte Carlo technique.

The procedure to compute the critical values for $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{1}}\right)$ is given step-by-step as follows:

• Generate a sequence of $m$ random vectors of the measurement errors for both epoch 1 ${\mathit{e}}_k^{\left\{1\right\}}$ and epoch 2 ${\mathit{e}}_k^{\left\{2\right\}},$ k = 1,..., m of the desired distribution. e.g.:

  

  (31)

  where $m$ is known as the number of Monte Carlo experiments. In addition, Matlab's “mvnrnd” command may be used in this step, for example.
• For each pair ${\mathit{e}}_k^{\left\{1\right\}}$ and ${\mathit{e}}_k^{\left\{2\right\}}$, k = 1,..., m compute the differences in errors between the two epochs, i.e.:

  

  (32)
• Apply the least-squares to estimate the observation errors, as follows:

  

  (33)

  where $\mathit{R}\in {ℝ}^{n\times n}$ is known as the redundancy matrix, which is given by:

  $$\mathit{R}=\mathit{I}-\mathit{A}{\left({\mathit{A}}^{\mathrm{T}}\mathit{W}\mathit{A}\right)}^{-1}{\mathit{A}}^{\mathrm{T}}\mathit{W}$$

  (34)

  with $\mathit{I}\in {ℝ}^{n\times n}$ the identity matrix.
• Assemble the displacement-design matrix ${\mathit{G}}_{\mathit{p}=\mathbf{1}}^{\left(g_i\right)}$,i = 1,..., n_p for each Monte Carlo experiment, according to equation (16), but the signs of the coefficients of the displacement-design matrix given by $\mathit{s}\mathit{g}\mathit{n}\left({\Delta }_{{\mathit{e}}_k}\right),$ k = 1,..., m, and compute the test statistic $\max \left({\mathit{T}}_{\mathit{p}={\mathbf{1}}_{\left(k\right)}}\right)$ by (18) and (22). The frequency distribution of $\max \left({\mathit{T}}_{\mathit{p}={\mathbf{1}}_{\left(k\right)}}\right)$ approximates the probability distribution of $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{1}}\right)$.
• Sort in ascending order the $\max \left({\mathit{T}}_{\mathit{p}={\mathbf{1}}_{\left(k\right)}}\right)$, getting a sorted vector $\mathit{T}\mathit{s}$, such that:

  

  (35)

  The sorted values in $\mathit{T}\mathit{s}$provide a discrete representation of the cumulative density function of the maximum test statistic $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{1}}\right)$.
• Determine the critical value $c$:

  

  (36)

  where $\left[·\right]$ denotes rounding down to next integer that indicates the position of the selected elements in the ascending order of $\mathit{T}\mathit{s}$. This position corresponds to a critical value for a stipulated overall false detection probability ${\alpha }_t$. This can be done for a sequence of values ${\alpha }_t$ in parallel [34].

The Matlab function "kslrtupi.m" was elaborated for the computation of critical values. Inputs in the function are $kslrtupi\left(\mathit{C},{\Sigma }_{\mathit{y}}^{\left\{1\right\}},{\Sigma }_{\mathit{y}}^{\left\{2\right\}},{\mathsf{\alpha }}_{\mathit{t}},m\right)$, where “$m$” is the user-defined number of Monte Carlo experiments.

3.2. Determining the maximum possible number of points $p_{max}$

Here, we objectively and universally demonstrate how to obtain the maximum number of points “$p_{max}$” to be inspected by the SLRTUPI sequential procedure. The maximum number of points is also defined sequentially. The procedure to define the maximum number of points “$p_{max}$” consists in finding a regular model (i.e., a matrix $\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{\mathit{p}}\end{array}\right)$ with full column rank) and, of course, that model has enough redundancy for the identification of unstable points. The step-by-step procedure to obtain the maximum possible number of points $p_{max}$ to be inspected by SLRTUPI is given by the flowchart in Figure 5.

If the full displacement-design matrix ${\mathit{G}}_{\mathit{p}}$ is specified by taking the total number of points in the network, i.e., by taking $p=n_p$ according to equation (16), then we fall into a problem in which there is a rank deficiency. Usually, the rank defect is greater than or equal to 1, because $rank\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{\mathit{p}={\mathit{n}}_{\mathit{p}}}\end{array}\right)\le n_p<u=n_p+1$, where $u$is the number of columns of the matrix $\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{\mathit{p}={\mathit{n}}_{\mathit{p}}}\end{array}\right)$. In this situation, the determinant of the normal matrix ${\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{\mathit{p}={\mathit{n}}_{\mathit{p}}}\end{array}\right)}^T\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{\mathit{p}={\mathit{n}}_{\mathit{p}}}\end{array}\right)$ is zero. This means that such matrix cannot be considered invertible – the matrix is then said to be singular. Consequently, it is not possible to compute the test statistics by Equation (18). This demonstrates that $p_{max}<rank\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{\mathit{p}={\mathit{n}}_{\mathit{p}}}\end{array}\right)\le n_p$.

By reducing the number of points by one unit, i.e., $p=n_p-1$, we expect to find regular models (i.e., without rank deficiency). For that, it should be checked again whether there is rank defect or not. But now we have at our disposal not only one single group of points, but a combination of $K_{n_p}^{n_p-1}=n_p$ group of $n_p-1$ points. If at least one of the matrices $\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{\mathit{p}}^{\left(g_i\right)}\end{array}\right), i=1,\dots ,n_p$ presents rank defect, then $p_{max}$ must be reduced again by one unit, i.e., $p=n_p-2$, and the rank computation is then ​​repeated for $K_{n_p}^{n_p-2}=\frac{n_p\left(n_p-1\right)}{2}$ group of $n_p-2$ points. If the rank defect persists, then the rank computation is repeated by decreasing $p$ to $p-1$, otherwise we proceed to evaluate whether the test statistics computed in Equation (18) are different from each other. If there are at least two statistics with equal values, then a statistical overlap will be flagged. In that case, therefore, we should reduce $p$ to $p-1$ and restart the procedure. If the models are found to be regular and there is no statistical overlap, then the process ends with the $p_{max}$ value found.

Important to highlight that the rank of the matrix $\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{\mathit{p}}\end{array}\right)$depends on the behaviour of the signs of the coefficients of its displacement-design submatrix ${\mathit{G}}_{\mathit{p}}$. The signs will also depend on the geometry of the geodetic network, as well as the unknown magnitude and number of points that have been shifted. Remember, however, $\Delta {\mathit{y}}_{\mathit{s}\mathit{g}\mathit{n}}$ captures these quantities. Although we have the limitation of computing the test statistics a priori, the procedure proposed above allows to determine the maximum number of testable points for identification. Generally, the choice of this maximum number is restricted to a subjective choice [11,13]. Here, Matlab function "pmax.m" was developed for the automatic computation of the maximum number of points. Inputs in the function are ${\mathit{y}}^{\left\{1\right\}}, {\mathit{y}}^{\left\{2\right\}},\mathit{C}, {\Sigma }_{\mathit{y}}^{\left\{1\right\}}$and ${\Sigma }_{\mathit{y}}^{\left\{2\right\}}$, i.e., $pmax\left({\mathit{y}}^{\left\{1\right\}}, {\mathit{y}}^{\left\{2\right\}},\mathit{C}, {\Sigma }_{\mathit{y}}^{\left\{1\right\}},{\Sigma }_{\mathit{y}}^{\left\{2\right\}}\right)$. Below are presented two possible scenarios that may occur in obtaining the $p_{max}$ for the geodetic network presented in the scope of this article (Figure 1).

Example 1:

(37)

Note that the number of columns in the matrix (38) is $u=7$​​ but its $rank\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{\mathit{p}=\mathbf{6}}\end{array}\right)=6$. In that case, we would fall into a rank deficiency problem, because $rank\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{\mathit{p}=\mathbf{6}}\end{array}\right)=6<u=7$. If we considered that $p_{max}=p=rank\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{{\mathit{p}}_{\mathit{m}\mathit{a}\mathsf{x}}=\mathbf{6}}\end{array}\right)=n_p=6$, then we would have $K_{n_p=6}^{p_{max}=6}=1$ group by Eq. (24), namely $g_1:\left[A,B,C,D,E,F\right]$. However, in this case, the determinant of the normal matrix would be equal to zero, denoted by$d\mathit{e}t\left[{\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{{\mathit{p}}_{\mathit{m}\mathit{a}\mathsf{x}}=\mathbf{6}}\end{array}\right)}^{\mathrm{T}}\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{{\mathit{p}}_{\mathit{m}\mathit{a}\mathsf{x}}=\mathbf{6}}\end{array}\right)\right]=0$, which would characterize it as a singular matrix, and consequently would not allow the computation of the test statistics by the Eq. (18). By reducing $p_{max}$ to one unit   ($p_{max}=p_{max=6}-1=5$), we would have now a total of $K_{n_p=6}^{ p_{max}=5}=6$ groups of 5 points, which would be the following: $g_1:=\left[A,B,C,D,E\right]$; $g_2:=\left[A,B,C,D,F\right]$; $g_3:=\left[A,B,C,E,F\right]$; $g_4:=\left[A,B,D,E,F\right]$; $g_5:=\left[A,C,D,E,F\right]$; and $g_6:=\left[B,C,D,E,F\right]$. All these groups would now have full column rank with $rank\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{{\mathit{p}}_{\mathit{m}\mathit{a}\mathsf{x}}=\mathbf{5}}\end{array}\right)=\mathrm{u}=6$ and $d\mathit{e}t\left[{\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{{\mathit{p}}_{\mathit{m}\mathit{a}\mathsf{x}}=\mathbf{5}}\end{array}\right)}^{\mathrm{T}}\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{{\mathit{p}}_{\mathit{m}\mathit{a}\mathsf{x}}=\mathbf{5}}\end{array}\right)\right]>0$. Thus, we would end and set $p_{max}=5$.

However, another important analysis is needed. Although the determinants of the matrices constructed for each group of points are non-zero, their values ​​are all equal to $d\mathit{e}t\left[{\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{{\mathit{p}}_{\mathit{m}\mathit{a}\mathsf{x}}=\mathbf{5}}\end{array}\right)}^{\mathrm{T}}\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{{\mathit{p}}_{\mathit{m}\mathit{a}\mathsf{x}}=\mathbf{5}}\end{array}\right)\right]=360$. Mostly, the test statistics computed for each of those groups by Eq. (18) would also have the same values. So, we would not be able to identify the group by the $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{5}}\right)$, since all the statistics are the same, i.e., there would be a statistical overlap. Consequently, we would only be able to get the information about displacement detection (i.e., $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{5}}\right)>k$), but not identification. In order to have the possibility of identification, we would then decrease $p_{max}$ again to one unit (i.e., $p_{max}=p_{max=5}-1=4$), which would lead to a total of $K_{n_p=6}^{p_{max}=4}=15$ groups of 4 points. Now, the test statistics ${\mathit{T}}_{{\mathit{p}}_{\mathit{m}\mathit{a}\mathit{x}}=\mathbf{4}}$ for each of $\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{{\mathit{p}}_{\mathit{m}\mathit{a}\mathit{x}}=\mathbf{4}}\end{array}\right)$ would result of values ​​different from each other, making identification possible. As a result, the maximum number of points to be considered would be $p_{max}=n_p-2$ for identification, which for that example we would get $p_{max}=4$.

Example 2:

(38)

In this scenario, we would fall back into a problem where the matrix $\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{\mathit{p}=\mathbf{6}}\end{array}\right)$ holds rank defect, but now with $rank\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{\mathit{p}=\mathbf{6}}\end{array}\right)=5$. If $p_{max}=n_p-1=5$, we would have $K_{n_p=6}^{ p_{max}=5}=6$ groups of 5 points to be considered, namely: $g_1:=\left[A,B,C,D,E\right]$; $g_2:=\left[A,B,C,D,F\right]$; $g_3:=\left[A,B,C,E,F\right];g_4:=\left[A,B,D,E,F\right];$     $g_5:=\left[A,C,D,E,F\right];$ and $g_6:=\left[B,C,D,E,F\right]$. Unfortunately, the rank of any of these groups would be equal to $rank\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{{\mathit{p}}_{\mathit{m}\mathit{a}\mathit{x}}=\mathbf{5}}\end{array}\right)=5<u=6$, and therefore we would still remain with rank deficiency problem. Moreover, the determinant of any normal matrix of these groups would be zero, i.e., ${\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{{\mathit{p}}_{\mathit{m}\mathit{a}\mathit{x}}=\mathbf{5}}\end{array}\right)}^{\mathrm{T}}\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{{\mathit{p}}_{\mathit{m}\mathit{a}\mathit{x}}=\mathbf{5}}\end{array}\right)=0$. By reducing $p_{max}$ to one unit (i.e., $p_{max}=p_{max=5}-1=4$), we would have a total of $K_{n_p=6}^{ p_{max}=4}=15$ groups of 4 points, but still with rank defect for the following groups: $g_1:=\left[A,B,C,D\right];g_2:=\left[A,B,C,E\right];g_3:=\left[A,B,C,F\right];$ and $g_{15}:=\left[C,D,E,F\right].$ As there would still be a lack of full column rank for some groups, we would again reduce $p_{max}=p_{max=4}-1=3$. Now we would have a larger number of groups than the previous case, but of size equals to $p_{max}=3$, i.e., $K_{n_p=6}^{ p_{max}=3}=20$ groups of 3 points. In the latter case, however, we would still have only one single group of points that would cause a rank defect, namely $g_1:\left[A,B,C\right]$. Again, we would proceed to decrease the maximum number of points followed by computing the rank of the matrices. Finally, we would find the maximum number of points, which in that case would be $p_{max}=p_{max=3}-1=2$. All these groups of 2 points would have their respective matrices $\left(\begin{array}{cc}\mathit{A}& {\mathit{G}}_{{\mathit{p}}_{\mathit{m}\mathit{a}\mathit{x}}=\mathbf{2}}\end{array}\right)$ with full column rank and theirs test statistics ${\mathit{T}}_{{\mathit{p}}_{\mathit{m}\mathit{a}\mathit{x}}=\mathbf{2}}$ would have different values from each other. Thus, $p_{max}=2$ would be the dimension found, so that all groups of 2 points would be testable for identification.

In addition to the examples given above, we may have cases where $p_{max}=1$, which allows us to identify only one single point. In that case, the SLRTUPI procedure would not proceed in its sequential form, but would restrict only in its first step, i.e., only in the identification of a point by means of $\max \left({\mathit{T}}_{\mathit{p}=\mathbf{1}}\right)$. And, if $p_{max}=0$, then in that case we would have a situation in which it is not possible to detect displacements, since the network does not have enough redundancy for this purpose. This latter case may occur if the geodetic network is very poorly designed so that it will not be possible to detect displacements.

4. Results from computational simulation-based approach and real dataset: a trilateration geodetic network

4.1. Simulation setup

For a first performance testing of the SLRTUPI, we use a simulated trilateration network measured in two epochs (Figure 4), which is the same network presented in the previous sections (Figure 1 and Figure 2).

The network points are assumed to be stable at Epoch 1 (Table 1), so measurements connected to them are undisturbed (Table 2). In that case, the “true” distances $d_{0_{ij}}^{\left\{1\right\}}$ for any segment formed by any two points $i$ and $j$ are easily computed as:

(39)

The true quantities in Table 2 are stacked into a vector as ${{\mathit{y}}_{\mathbf{0}}}^{\left\{1\right\}}$. To have simulated observations for Epoch 1 $\left({\mathit{y}}^{\left\{1\right\}}\right)$ – see Equation 1 –, random errors ${\mathit{e}}_{\mathit{y}}^{\left\{1\right\}}$are synthetically generated based on a multivariate normal distribution and added up to the true distances, as follows:

(40)

The observations are assumed to be uncorrelated and with the same known standard deviation of $\sigma =2 mm$, which is a value compatible with real cases. Then the main diagonal elements of the variance matrix for the Epoch 1 are the variances with values equal to ${\sigma }^2=4 m{m}^2$, and its off-diagonal elements are zeros. Here, we use the Mersenne Twister algorithm to generate a sequence of random numbers and Box–Muller to transform it into a normal distribution [41,42]. For instance, Matlab's “mvnrnd” command can be applied to generate the random errors.

For Epoch 2, on the other hand, measurement distortions are simulated from the intentional displacement of the network points. The point displacements are simulated in terms of magnitude, denoted by $\mathbf{\nabla }$, and orientation, denoted by $\theta$, so that they can be propagate to the distance equation. To illustrate this step, consider that the segment formed by any two points$i$ and $j$ has suffered distortion when having point $i$ displaced to its new position $i\prime$(Figure 5). In this case, the true disturbance distance in Epoch 2 is computed as:

(41)

$${\nabla }_{X_i}={\nabla }_{i}sin\theta \mathrm{and} {\nabla }_{{\mathit{y}}_i}={\nabla }_{i}cos\theta$$

(42)

Figure 5. Geometric aspect of simulating the displacement of a point and its effect on distance distortion.

Figure 5. Geometric aspect of simulating the displacement of a point and its effect on distance distortion.

One can generalize the displacement simulation of any point or group of points from the Point-to-Measurement Connection Matrix $\mathit{C}$ (see Equation 14) from the following mathematical operations: first, we store the coordinates of the network points with no displacement into a vector $\mathit{V}\in {ℝ}^{2\times n_p}$. In our example, we have $n_p=6$, which leads to:

$$\mathit{V}=\left(\begin{array}{c}{X}_A\\ Y_A\end{array} \begin{array}{c}{X}_B\\ Y_B\end{array} \begin{array}{c}{X}_C\\ Y_C\end{array} \begin{array}{c}{X}_D\\ Y_D\end{array} \begin{array}{c}{X}_E\\ Y_E\end{array} \begin{array}{c}{X}_F\\ Y_F\end{array}\right)$$

(43)

Next, we use the Point-to-Measurement Connection Matrix in a modified form, denoted by $\overline{\mathit{C}}\in {ℝ}^{n\times n_p}$, with the signs of its coefficients resulting from the orientation of the segment. For example, from point A towards point D, we would have $\mathit{A}\mathit{D}=\mathit{D}-\mathit{A}$ with its corresponding components as ${{\mathsf{\Delta }}_X}_{AD}=X_D-X_A$ and ${{\mathsf{\Delta }}_Y}_{AD}=Y_D-Y_A$, and therefore, the first column of matrix ${\overline{\mathit{C}}}^{ T}$ would be ${\left(\begin{array}{cccccc}\underset{A}{\underbrace{-1}}& \underset{B}{\underbrace{ 0 }}& \underset{C}{\underbrace{ 0 }}& \underset{D}{\underbrace{ 1 }}& \underset{E}{\underbrace{ 0 }}& \underset{F}{\underbrace{ 0 }}\end{array}\right)}^T$. Following this reasoning, we would then have the general matrix $\overline{\mathit{C}}$ in its transposed form for the case of the network in Figure 4 as

(44)

In general, therefore, we can simulate the disturbed distances in Epoch 2 by means the displacement of a point or groups of points from the relationship below:

$${\mathit{y}}_{\mathbf{0}}^{\left\{\mathbf{2}\right\}}={{\Vert \left(\mathit{V}+\left[\begin{array}{c}{{\mathbf{\nabla }}_X}_{g_i}\\ {{\mathbf{\nabla }}_Y}_{g_i}\end{array}\right]\right){\overline{\mathit{C}}}^{ T}\Vert }_{\mathbf{2}}}^T$$

(45)

where ${\Vert ·\Vert }_{\mathbf{2}}$ represents the Euclidean norm (also called the vector magnitude, Euclidean length, or 2-norm) of each column of a given matrix. In other words, this vector-wise 2-norm provides the distances of each segment of the network; ${{\mathbf{\nabla }}_{\mathit{X}}}_{g_i}\in {ℝ}^{1\times n_p}$ and ${{\mathbf{\nabla }}_{\mathit{Y}}}_{g_i}\in {ℝ}^{1\times n_p}$ are the displacement vectors decomposed in the X and Y direction, which consist exclusively of elements with values of ${\nabla }_X$ and ${\nabla }_Y$ for the selected $g_i$ group of points, respectively (see Equation 42), and 0 for points assumed to be free of displacements. The index “$i$” in $g_i$ indicates a possible group, i.e., $\left\{g_i\right| i=1,2,3,\dots ,K_{n_p}^p\}$, as described in the previous section (see Equation 24). Note that the matrix $\mathit{V}$ composes the coordinates of the points with no displacements, whereas the second part of the Equation (45) represents the magnitude of the simulated displacements in the X-direction ${\mathbf{\nabla }}_{\mathit{X}}$ and Y-direction ${\mathbf{\nabla }}_{\mathit{Y}}$ for a group of points $g_i$.

For example, if we want to simulate the displacement of two points (i.e., $p=2$), then for the considered network we would have a total of $K_{n_{p=6}}^{p=2}=15$ groups, i.e., $\left\{g_i\right| i=1,2,3,\dots ,15\}$, namely: $g_1:=\left[A,B\right]$,    $g_2:=\left[A,C\right]$, $g_3:=\left[A,D\right]$, $g_4:=\left[A,E\right]$, $g_5:=\left[A,F\right]$, $g_6:=\left[B,C\right]$, $g_7:=\left[B,D\right]$, $g_8:=\left[B,E\right]$, $g_9:=\left[B,F\right]$,    $g_{10}:=\left[C,D\right]$, $g_{11}:=\left[C,E\right]$, $g_{12}:=\left[C,F\right]$, $g_{13}:=\left[D,E\right]$, $g_{14}:=\left[D,F\right]$ and $g_{15}:=\left[E,F\right]$. To make it clearer, let's take the following specific example of having a simultaneous displacement of points A and D. i.e., for the group $g_3:=\left[A,D\right]$, with the following fixed definitions ${\nabla }_{\mathrm{A}}=2\mathrm{cm}$ @ ${\theta }_A=45°$ and ${\nabla }_{\mathrm{D}}=4\mathrm{cm}$ @ ${\theta }_D=135°$, respectively. In that case, then we would have:

(46)

Thus, the ‘true’ distances in Epoch 2 would be:

$${{\mathit{y}}_{\mathbf{0}}}^{\left\{2\right\}}=\left(\begin{array}{c}129.7700\\ 130.1843\\ 113.5220\\ 129.0901\\ 163.0530\\ 181.8570\\ 91.8161\\ 62.1665\\ 93.7482\end{array}\right)$$

(47)

Next, the observations at the Epoch 2 are generated in the same way as it was done for Epoch 1 (see Equation 1), as follows:

(48)

The variance matrix for the Epoch 2 is ​​taken to be the same as in Epoch 1, i.e., ${\Sigma }_{\mathit{y}}^{\left\{1\right\}}={\Sigma }_{\mathit{y}}^{\left\{2\right\}}$. Here, the value of $\sigma$ is unimportant in this investigation. Now, we have all the necessary elements to apply SLRTUPI, as can be seen in the flowchart displayed in Figure 3. However, the simulation procedure explained so far is for only one single Monte Carlo experiment (i.e., $m=1$). Thus, a sequence of $m$ random observation vectors for both ${\mathit{y}}^{\left\{1\right\}}$ (Epoch 1) and ${\mathit{y}}^{\left\{2\right\}}$ (Epoch 2) are needed to evaluate the success and failure rates of the SLRTUPI. Here, therefore, we set up $m=200,000$ Monte Carlo experiments as suggest by [43], such that:

(49)

Remember that random vectors ${\mathit{y}}^{\left\{2\right\}}$ consist of the disturbed measurements due to the displacements. Displacements can vary in terms of their magnitude ($\mathbf{\nabla }$) and orientation ($\mathsf{\theta }$) as well as the number of points ($p$). Each scenario is simulated individually, i.e., we generate $m=200,000$ Monte Carlo experiments for each combination of ($\mathbf{\nabla },\mathsf{\theta },p$). Here, we evaluate for cases where the number of displaced points is $p=1, p=2$ and $p=3$ for the network in Figure 4.

The critical values were obtained for cases having the following significance levels: ${\alpha }_t=0.001 \left(c=16.75\right)$; ${\alpha }_t=0.01 \left(c=12.27\right)$; ${\alpha }_t=0.05 \left(c=9.06\right)$ and ${\alpha }_t=0.1 \left(c=7.62\right)$, as displayed in Table 3.

These critical values were computed by 2 million of Monte Carlo experiments such that a user-defined Type I decision error ${\alpha }_t$ for SLRTUPI was warranted (see topic 3.1.). To verify the consistency of these critical values found, we ran 200,000 Monte Carlo experiments under the condition of absent of displacements. In this case, the null hypothesis is fixed as true, and the detection of displacements represents the error of SLRTUPI in detecting displacements, when in fact they do not exist (false positives). Table 3 shows the significance level values found for these critical values, denoted by ${\alpha }_t^{\prime }$. The ones found must be as close as the ones setup (see the fourth column of the Table 3). Therefore, it is observed that the control of false positives was efficient by means of Monte Carlo Method, as described in Topic 3.1.

4.2. Numerical example for individually displaced points $p=1$

In the first analysis, we have considered the case of having only one displaced point $p=1$ at a time. In that case, we had then 6 groups individually simulated, namely $\{g_i| i=1,2,3,\dots ,6\}$, with each group consisting of only one point, i.e.: $g_1:=\left[A\right], g_2:=\left[B\right], g_3:=\left[C\right],g_4:=\left[D\right],g_5:=\left[E\right],g_6:=\left[F\right]$. The magnitude of the displacement for each group $\left\{{\nabla }_{g_i}\right|i=1,2,3,4,5,6\}$ was defined within the interval of [1,10]$\sigma$ with increments of 1$\sigma$, where $\sigma$ is the standard deviation of the observation, taken as $\sigma =2mm$ – as mentioned before. The displacements were simulated to act in different orientations. The displacement orientation $\theta$ was simulated for the range of [0°,355°] with 5° increment. As a result, we had ten magnitudes of displacement for each of the 72 orientations, totalling 720 simulated displacements for each point (Figure 6).

200,000 Monte Carlo experiments were generated for each of those 720 cases, totalling 144 million trials for each point. It is important to remember that the simulated displacements of the points were converted to the space of the observations by Equation (45), which resulted in the corresponding disturbances on the measurements. Then, the Epoch 2 observations were generated from Equation (48), considering those true disturbance distances observations.

Finally, the probabilities levels associated with SLRTUPI (denoted by $𝒫$[.]) were computed as the ratio between the occurrence of a particular event – correct detection, wrong identification, overidentifications, statistical overlap or correct identification – and the number of possible cases (i.e., total number of Monte Carlo experiments 𝑚=200,000), which are given respectively by [38]:

$${\mathcal{P}}_{CD}=\frac{n_{CD}}{m}$$

(50)

$${\mathcal{P}}_{WI}=\frac{n_{WI}}{m}$$

(51)

$${\mathcal{P}}_{OI+}=\frac{n_{OI+}}{m}$$

(52)

$${\mathcal{P}}_{OI-}=\frac{n_{OI-}}{m}$$

(53)

$${\mathcal{P}}_{Ol}=\frac{n_{Ol}}{m}$$

(54)

$${\mathcal{P}}_{CI}=\frac{n_{CI}}{m}$$

(55)

where:

• $n_{CD}$ – number of correct detections is the number of experiments in which SLRTUPI procedure correctly detect displacements. Inversely we have the missed detection rate for the case where SLRTUPI does not detect displacements, as follows:

  

  (56)
• $n_{WI}$ – number of wrong identifications is the number of experiments in which SLRTUPI flags only one single point as being unstable while the ‘true’ unstable point remains unidentified.
• $n_{OI+}$– number of overidentifications positive is the number of experiments in which SLRTUPI identifies more than one point as being unstable, and among these points there is one that is correctly identified as unstable.
• $n_{OI-}$– number of overidentifications negative is the number of experiments in which SLRTUPI identifies more than one point as being unstable, whereas the ‘true’ unstable point remains unidentified.
• $n_{Ol}$ – number of statistical overlaps is the number of experiments in which the detector $\max \left({\mathit{T}}_{\mathit{p}}\right)$ flags two (or more) group of points simultaneously during a given iteration of SLRTUPI.
• $n_{CI}$ – number of correct identifications is the number of experiments in which SLRTUPI correctly identifies the displaced point. In that case, we have the following relationship for the success rate of correct identification of the displaced point:

  

  (57)

  $${\mathcal{P}}_{CI}={\mathcal{P}}_{CD}-\left({\mathcal{P}}_{WI}+{\mathcal{P}}_{OI+}+{\mathcal{P}}_{OI-}+{\mathcal{P}}_{Ol}\right)$$

  (58)

For identification to exist, detection must have occurred (see Equation 57). It is important to highlight that “displacement detection” only informs us whether or not there might have been at least one unstable point. For identification to exist, detection must have occurred. Detection does not guarantee correct identification, and it is clearly noted that identification is more difficult than detection (see Equation 58). This is explained in detail by Rofatto et al. [39]. Here it is also important to highlight that on average ~95% of the total experiments had $p_{max}=4$; no case occurred when $p_{max}=3$; ~3% had $p_{max}=2$; and ~2% had $p_{max}=1$, according to the computation given by topic 3.2. Therefore, there were rare cases where $p_{max}<4$.

Figure 7 shows the result of the correct identification. The black circles represent the radial range of magnitude displacements of 4mm, 1cm and 2cm. It is important to note that uncertainty of the two-epoch geodetic observations differences is ${\mathsf{\sigma }}_{{\mathsf{\Delta }}_y}~2.82mm$. Therefore, these magnitudes of 4mm, 1cm and 2cm correspond to approximately $1.4{\sigma }_{{\mathsf{\Delta }}_y}$, $2.8{\sigma }_{{\mathsf{\Delta }}_y}$ and $7{\sigma }_{{\mathsf{\Delta }}_y}$.

Clearly, it is observed that the probability of correct identification depends on the geometry of the connections. Note that the highest success rates occur in the lines of sight. On the other hand, it is more difficult to identify in the directions perpendicular to the lines of sight, mainly for magnitudes close to the uncertainty of the measurement method $\left(\nabla \le 4mm\right)$. The increase in the significance level, and consequently a larger null hypothesis rejection region (smaller critical values) favours the identification of low magnitude displacements. On the other hand, higher significance levels reduce (slightly) the probability of correctly identifying displacements of larger magnitudes.

Figure 8 shows the result of the correct detection. Detection follows the same behaviour as identification in terms of the influence of the geometry of the point connections. There are no significant differences in terms of detection and identification for the case in which a low significance level is adopted $({\alpha }_t=0.001$). The lower the critical value (higher significance level) the higher the detection rate.

Figure 9 provides the result for the wrong identification (also called Type III decision error). This decision error gets more expressive when the level ${\alpha }_t$ is increased. In the worst case, we have ${\mathcal{P}}_{WI}\cong 0.4 \left(40\%\right)$ for the case where point E moves and for ${\alpha }_t=0.01$. In general, the wrong identification occurs for large magnitudes in the region outside the lines of sight.

The overidentification class ${\mathcal{P}}_{OI+}$ is presented in Figure 10. In general, such probability level become more evident when increasing ${\alpha }_t$. The overidentification ${\mathcal{P}}_{OI-}$ was much rarer to occur, with its highest rate value of ${\mathcal{P}}_{OI-}\cong 0.022 \left(2.2\%\right)$for ${\alpha }_t=0.1$. So, it is not shown here. On the other hand, the overidentification positive ${\mathcal{P}}_{OI+}$ seems to occur more frequently than ${\mathcal{P}}_{OI-}$, especially close to the lines of sight and for large magnitudes.

Figure 11 portrays the behavior of the cumulative probability $P\left(x\right)$ from the empirical cumulative distribution function (ECDF) of the identification and detection success when taking ${\alpha }_t=0.1\%,{\alpha }_t=1\%,{\alpha }_t=5\%$ and ${\alpha }_t=10\%$ for small (4mm), medium (1cm) and large (2cm) magnitude of displacement.

Increasing ${\alpha }_t$ improved identification in the case of small (4mm) and medium (1cm) magnitude of displacement. Figure 11(a) reveals that approximately 90% ($P\left(x\right)=0.9$) of the ${\mathcal{P}}_{CI}$ rates for $\nabla =4mm$ achieved $\cong 2\%$ or less for ${\alpha }_t=0.001$ (blue line); $\cong 5\%$ or less for ${\alpha }_t=0.01$ (red line); $\cong 12.5\%$ or less for ${\alpha }_t=0.05$ (yellow line); and $\cong 17.5\%$ or less for ${\alpha }_t=0.1$ (purple line). In the case of medium magnitude of displacement ($\nabla =1cm$), approximately 90% ($P\left(x\right)=0.9$) of the ${\mathcal{P}}_{CI}$ rates achieved $\cong 68\%$ or less for ${\alpha }_t=0.001$ (blue line); $\cong 82\%$ or less for ${\alpha }_t=0.01$ (red line); $\cong 89\%$ or less for ${\alpha }_t=0.05$ (yellow line) and ${\alpha }_t=0.1$ (purple line), as highlighted in Figure 11(b). The impact of the ${\alpha }_t$ choice on the correct identification rate is not as critical for large magnitudes. Figure 11(c) shows that approximately 50% ($P\left(x\right)=0.5$) of the ${\mathcal{P}}_{CI}$ rates for $\nabla =2cm$ achieved $\cong 98\%$ or less for ${\alpha }_t=0.001$ (blue line) and $\cong 82\%$ or less for ${\alpha }_t=0.1$ (purple line). In that latter case, therefore, it is preferable to choose a lower familywise error ate ${\alpha }_t$ to have higher identification success rate, i.e., ${\alpha }_t=0.001$. Detection always improves as the critical value increases, at the cost of having to increase false positive rates ${\alpha }_t$. In the next section, we present the testing performance of the SLRTUPI for the case of having 2 points displaced simultaneously.

4.3. Numerical example for simultaneously displaced points $p=2$

In this second analysis, the scenarios were created for the case of having two points simultaneously unstable. For this, we have simulated all possible cases for $p=2$ at a time, i.e., we had then 15 groups of 2 points individually simulated $\{g_i| i=1,2,3,\dots ,15\}$. The magnitude of the displacement was fixed for each group as being of 10$\sigma =2cm$ $\{{\nabla }_{g_i}=10\sigma |i=1,2,3,\dots ,15\}$. The displacement patterns of the two-point group varied in the same direction and opposite directions, as detailed in Figure 12. Table 4 and Table 5 display the results for correct identification and detection, respectively. [note: the results for each probability level are represented by $\frac{max}{min}$ , where max. and min. are the largest and smallest values].

In terms of detection and identification success rates, SLRTUPI depends on the interaction between the pattern in which the points move and the geometry of the geodetic network (Table 4 and Table 5). One geometry may be best for a given displacement pattern that one wish to monitor. The patterns (d), (e), (g) and (i) are the most difficult to identify for the case of the trilateration network considered in this work. This is due to the fact the displacements occur close to the perpendicular direction to the lines of sight (region more difficult to identify), as well as close to the direction where the Wrong Identification rate is larger, as also can be seen in the previous section for the case $p=1$. On the other hand, displacement patterns with the pair of points displaced in different directions (c, f, j, k, and l in Figure 12) seem to be easier to identify. Furthermore, increasing the significance level does not always improve the identification rate of unstable points, but always improves detection (Table 5). In this scenario, we did not also have the occurrence of statistical overlap. The detection rate was very high (Table 5), so we can infer that it is very rare not to detect the displacement patterns simulated here.

4.4. Numerical example for simultaneously displaced points $p=3$

In this last analysis under simulated conditions, the scenarios were created for the case of having three points simultaneously simulated as unstable. For this, we have simulated all possible cases for $p=3$ at a time, i.e., we had then 20 groups of 3 points individually simulated $\{g_i| i=1,2,3,\dots ,20\}$. The magnitude of the displacement was also fixed for each group as being of 10$\sigma =2cm$ $\{{\nabla }_{g_i}=10\sigma |i=1,2,3,\dots ,20\}$. The displacement patterns of the three-point group varied in the same direction and opposite directions, as detailed in Figure 13. Table 6 and Table 7 display the results for correct identification and detection for the case of having three mutually displaced points $p=3$, respectively.

The identification for the case where we have 3 simultaneous displacements is more difficult than the case of 2 simultaneous displacements (see Table 4 and Table 6), but the correct detection rates are similar (see Table 5 and Table 7). In fact, most of the patterns simulated here for $p=3$ are difficult to identify. The pattern described in (c, k, and l) are the easiest to identify (Table 6).

The numerical examples described here provide a way in which the user may design the geodetic network as best as possible to find the optimum geometry to identify a certain type of displacement pattern. In the next section, we will evaluate the performance of SLRTUPI for simulations of real displacements in the field.

4.4. Real example

In this experiment, the geodetic network described in Figure 4 was materialized in the field. For this, the distances were obtained using the total station FOIF OTS 685, with linear uncertainty of 2mm + 2ppm (manufacturer specifications). Tribrach were used for both the Total Station and the Reflector Prisms. The supplement files are provided for those interested in the dataset for reproducing the experiments or even for testing other procedures. Epoch 1 was recorded as no displacements, while Epoch 2 various patterns were tested. For this, the displacements at the points were intentionally applied radially. Points A, B and C were kept fixed, i.e., they were not subjected to displacement, whereas D, E and F were considered the points to be monitored. The displacements were performed in the field, by moving the reflector 1 cm from the initial position (Epoch 1) to the new positions (Epoch 2), as indicated in Figure 14. Four cases were tested for shifting only one point at a time, as displayed in Figure 14(a,b,c,d); three scenarios for 2 and 3 simultaneous displacements, with the patterns shown in Figure 14(e,f,g) and Figure 14(h,i,j), respectively, totalling 23 experiments (12 for $p=1$; 8 for $p=2$; and 3 $p=3$).

First, SLRTUPI was applied under the condition that all points are unstable, which in that case Point-to-Measurement Connection Matrix $\mathit{C}$ is the one given in (14). Based on the simulated results, we adopted a significance level of ${\mathsf{\alpha }}_t=10\%$, which resulted in a critical value of $c=7.62$. The results are very promising, as can be seen in Table 8. We had a success rate in identification of ~74%, and the correct detection rate was 100% for all scenarios.

The results show that of the 23 experiments, only 6 wrong identifications occurred. Most experiments resulted in $p_{max}=4$, with a few exceptions, such as: [D, E, F] in Figure 15h resulted $p_{max}=1$, consequently the wrong identification was already expected; and [D, E, F] in Figure 15i resulted $p_{max}=3$.

Here, we also apply SLRTUPI based on the a priori knowledge that points A, B and C are known to be stable. In the latter, therefore, the matrix $\mathit{C}$ in (14) was reduced to:

(59)

Because we had a change in the mathematical model in terms of redundancy – 3 points less to be monitored than the previous case, consequently there was better redundancy than monitoring all points –, the critical value found was $c$ = 6.64 for ${\alpha }_t$ = 10%. As a result, the correct identification rates increased from ~74% to ~96% when points A, B and C were fixed (adopted as control points), as can be seen in Table 9. We had only one single case of Type III Error among the 23 experiments. This is due to the fact of increasing redundancy, since we have gone from 6 points to 3 points to be monitored for the same set of observations. Furthermore, $p_{max}$ = 3 for all experiments where A, B, C were taken as fixed.

Finally, SLRTUPI was applied under a scenario in which there are no displacements to verify the efficiency of the Type I Error control (specifically, of the familywise error rate ${\alpha }_t=0.1$). For this, new measurements were performed and stamped to Epoch 2, but now without applying any intentional displacement to the points. As a result, no displacement was detected in the field experiments, which shows that SLRTUPI can guarantee user-defined Type I Error control under real conditions of use.

5. Contributions

In this contribution, we present a statistical method for detecting and identifying unstable points for geodetic deformation analysis, namely Sequential Likelihood Ratio Test for Unstable Points Identification (SLRTUPI). The method makes use of the differences of measurements taken at two epochs in time instead of using the classic difference of estimated coordinates. As an advantage, we will always have a linear model, and therefore the model nonlinearity problem does not affect the test power. In addition, another advantage is to avoid the S-transformation to maintain the same datum between epochs.

The method is not restricted to points subject to monitoring but is applied in such a way that all points participate in its inspection. However, the power of the identification improves considerably if we know in advance which points are stable (see Table 8 and Table 9). The results support that the success of the proposed method depends on the interaction between the displacement pattern and the network geometry. In the case of the trilateration network used as a case study in this contribution, we observe that for the case of having only one single point displaced, the highest success rates occur in the lines of sight (Figure 7 and Figure 8). When more than one point is unstable simultaneously, success rates are generally higher for the scenarios where the displacements occur in opposite directions (see Table 4, Table 5, Table 6 and Table 7). More experiments based on different types of geodetic networks (such as, levelling and GNSS networks) will be conducted in the future works.

The larger the number of points displaced simultaneously, the more difficult the identification. However, success rates depend on the geometry and redundancy of the network. Furthermore, detection is always larger than identification. Although identification is important to localize displacement, detection plays a dominant role for the activities that involves risk assessment, such as the deformation analysis. It is preferable that detection occurs even if identification does not occur in the expected way. It is crucial that a system raises an alarm when displacements are detected, even though their location may be wrong. In this sense, we also observe that increasing the significance level does not always improve the identification rate, but always improves detection.

It is important to highlight here that the false positive rates (Type I errors) are efficiently controlled from the definition of the familywise error rate by the user (${\alpha }_t$). The approach for controlling that decision error is based on Monte Carlo. Of course, this requires some computational cost, but nothing prohibitive. One way to save computational costs would be to develop a methodology based on artificial neural networks, as done for example by [44], or some other alternatives, such as [45]. This will be investigated in future work.

Another important contribution is that now the maximum number of points $p_{max}$ to be inspected by the SLRTUPI procedure is not based on non-objective choices, but rather determined from the rank analysis of the design matrices involved and the occurrence of statistical overlap. Consequently, the method avoids the occurrence of statistical overlap, which is common in multiple test problems [39,45]. More experiments will be needed to understand the effect of determining $p_{max}$ on the detection and identification of unstable points.

The simulation methodology described here can also be applied in the design stage of a geodetic network for deformation analysis purposes. This allows the user to know a priori how the measurement system (in this case, the geodetic network) should be optimally designed so that it is able to detect and/or identify a certain displacement pattern for a given probability. As a result, reliability measures can be easily extracted, such as Minimal Identifiable Displacement (MID) and Minimal Detectable Displacement (MDD) from its corresponding success rates ${\mathcal{P}}_{CI}$ and ${\mathcal{P}}_{CD}$ by a given user-defined significance level ${\mathsf{\alpha }}_t$. More experiments for different geodetic networks will be needed to evaluate the ability of the measurement system to identify possible outliers (internal reliability) and the effect of these undetected errors on the quality of deformation analysis results (external reliability).

In future works we will also evaluate the question of formulating the null hypothesis with parameter-free, so that the differences between the measurements in the two epochs may be taken as being directly the errors. Thus, the test could be simpler, with the benefit that the degrees of freedom of the test will be larger than the model adopted here.

Finally, it is emphasized that the proposed method can be extended to the outlier detection problem. The proposed method can also be extended to the case of Point Clouds from Terrestrial Laser Scanning in case of having to decide between the different models of surface representation in the area-based deformation analyses problems. The dataset and algorithms used in this work are available for those interested in reproducing the results and/or applying them to other methodologies, as follows: https://data.mendeley.com/datasets/msg783rh2y/draft?a=607de244-7c1d-417f-a9c8-81281d8a6056