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.

Let's start by describing the univariate model for identifying unstable points in geodetic networks by the following equations [32,34]: 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: 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.: the null hypothesis ${\mathcal{H}}_0$ of the standard Gauss–Markov model in the linear form of the Eq. (2) is then given by: 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: and the estimated observation errors as being: 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: 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: 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:

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:

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:

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:

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: 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: where $diag\left(.\right)$ is the diagonalize operator which convert a vector to a diagonal matrix. Consequently, the alternative hypothesis is now given by:

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]:

Then, a test decision is then performed as

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.

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.: 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.:

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

The decision rule for this case is given by:

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: 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: 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: 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]:

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: 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:

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.