Abstract
Deformation congruence models form the basis for conventional deformation analysis (CDA). In geometrical sense, these models connect an epochal object states—represented by its characteristic points—at stable/congruent points to disclose possible deformations. To this day, the deformation congruence models are usually specified using the global congruence test (GCT) procedure which, however, has a weakness in the case of multiple displacements. More precisely, the GCT procedure is based on consecutive point-by-point specification which may suffer from so-called displacement smearing. To overcome the above weakness, a revolutionary—in the context of GCT—concept (two methods) involving combinatorial possibilities was suggested in recent years. Admittedly, this concept avoids the problem of consecutive point-by-point specification. Nevertheless, it generates another weakness, namely the problem of the comparison of different-dimensional models. This paper makes a step forward in this new combinatorial field and discusses a more sophisticated combinatorial procedure, denoted as CIDIA. It was shown that, thanks to an appropriately used the possibilities of combinatorics and generalized likelihood ratio tests performed in the detection–identification–adaptation (DIA) iterative steps, the above weaknesses can be overcome. In the context of GCT, the suggested procedure has rather evolutionary—than revolutionary—character and the general concepts of both procedures have similar heuristic substantiation. To demonstrate the efficacy of CIDIA against GCT and the two existing combinatorial methods, various deformation scenarios were being randomized independently many times with the use of comprehensive computer simulations and then processed. Generally, the obtained results confirmed the statement that the suggested CIDIA procedure—unlike the existing combinatorial methods—can be substantially more resistant to displacement smearing than the GCT procedure, at no significant costs. The efficacy of CIDIA—unlike the ones of the two existing combinatorial methods—turned out always higher (on average by several percentages) than the one of GCT for all considered deformation scenarios. At the same time, the CIDIA procedure turned out substantially less time-consuming than the other combinatorial methods.
Similar content being viewed by others
1 Introduction
Geodetic absolute or relative control networks—which are measured in two or more epochs of time—form the framework for conventional deformation analysis (CDA), (see, e.g., Welsch and Heunecke 2001). Absolute networks consist of points placed on the deformable object and points placed outside this object. Relative networks are established when the surroundings of the deformable object do not show stability. These networks consist only of points placed on the deformable object. Absolute networks are most often employed in the deformation analysis of engineering structures, while relative networks are most often employed in the deformation analysis of the surface of Earth’s crust (e.g., Chrzanowski and Chen 1990; Caspary 2000).
Practically, CDA consists in the congruence analysis of the state of an object, represented by its characteristic points, at two—or more—measurement epochs (Welsch and Heunecke 2001). In geometrical sense, control networks adjusted in two—or more—epochs are here connected at stable/congruent points to disclose the displacements of possible unstable points. It means that the stable part of network, also called reference base, forms the basis for the deformation estimates of remaining part. There are two widely accepted approaches to the specification of deformation congruence models: the global congruence test (GCT) procedure (Pelzer 1971, 1974; Niemeier 1981, 1985; Kok 1982; Heck 1983; Caspary 2000, pp. 116–130; Neumann and Kutterer 2007) and robust approach (Chen 1983; Caspary and Borutta 1987; Chen et al. 1990; Caspary et al. 1990; Caspary 2000, pp. 130–132; Nowel and Kamiński 2014; Nowel 2015). The congruence analysis with the use of GCT is performed similarly as the outlier screening with the use of iterative data snooping procedure (e.g., of Baarda 1968). GCT uses a so-called mean shift model and employs the least squares (LS) estimation supported by statistical hypothesis tests, in consecutive iteration steps. However, the general concept of the procedure has heuristic substantiation, just as the general concept of iterative data snooping procedure (e.g., Kok 1984, p. 33; Caspary 2000, p. 77). Displaced points are here consecutively one-by-one removed from the congruent part of the model to specify the final model. Whereas the robust approach assumes a so-called variance inflation model and employs a robust M-estimation. Displaced points are not removed from the congruent part of the model, but their participation in this part is gradually suppressed in consecutive iteration steps. A method recently published, which can be regarded as the extension of the robust approach, is presented by Nowel (2019). The idea behind this method is inspired by the theory of Msplit estimation (Wiśniewski 2009, 2010; Wiśniewski and Zienkiewicz 2020; Duchnowski and Wiśniewski 2019, 2020; Wyszkowska and Duchnowski 2019, 2020; Zienkiewicz 2019) and lies in the assumption that control networks may simultaneously realize many competitive congruence models.
The GCT procedure is more widely accepted, long-standing experience has proven its high efficacy and it is recommended in most journal papers and all significant geodetic textbooks comprising deformation analysis (e.g., Caspary 2000, pp. 116–172; Niemeier 2008, pp. 433–457 or Heunecke et al. 2013, pp. 488–521). However, GCT has a serious weakness in the case of multiple displacements (e.g., Hekimoglu et al. 2002, 2010). Namely, alternative hypotheses imply only one single displaced point in a given iteration step and then the alternative hypotheses which are considered in the next steps are conditioned by the alternative hypothesis identified in the previous step. More precisely, the congruent part of model is point-by-point reduced in consecutive iteration steps and the point which is identified as most suspected to be displaced in a given step is not then modeled as stable in the next steps. The point most suspected to be displaced in a given step is the one which returns the greatest contribution to the global test statistic. The problem is that the contributions are expressed by the quadratic forms of estimated displacements which may be substantially burdened in the case of multiple displacements. Firstly, and definitely most important, the possible displacements of points modeled as stable (in congruent part) smear—more or less—on the estimated displacements of the points modeled as displaced and, in consequence, on their contributions. Secondly, since estimated displacements are interrelated through a network, the possible displacements of points modeled as displaced also slightly smear on the estimated displacements of other points modeled as displaced and, in consequence, on their contributions (see, e.g., Hekimoglu et al. 2002, 2010; Erdogan and Hekimoglu 2014; Durdag et al. 2018). The above smearing may lead to the identification of a stable point as most suspected to be displaced in a given iteration step. If it is the case, the models which are specified in the next steps are—from now—conditioned by the result of the erroneous identification from the previous step. It means, all these models are—from now—conditioned by the erroneous information. The computational load of the GCT procedure is relatively low but it is done at the above cost. Note that the above weakness of GCT under multiple displacements has similar nature as the weakness of iterative data snooping procedure under multiple outliers (e.g., Kok 1984, pp. 35–36; Ding and Coleman 1996).
Recently, a revolutionary LS concept—which overcomes the above weakness—was presented by Velsink (2015) and Lehmann and Lösler (2017). Since the advent of powerful computers, non-iterative combinatorial concept—where all possible combinations of displaced points come into considerations—may be today performed in near-real time. It is not necessary here to consecutively point-by-point specify the congruence model, because all candidate congruence models can be compared at the same stage. Additionally, this concept allows imposing some geometrical conditions to displaced points (extended congruence analysis), which is not possible in GCT or robust approach. For instance, it is the case when some displaced points can be moved in the same way. Such a combinatorial method is discussed in detail by Velsink (2015). Generally, the author considers many congruence models and proposes the ratio of test statistic and its critical value as a decision rule. This means that among all candidate models, the one associated with the largest ratio pinpoints the best approximating model. Another combinatorial method is discussed by Lehmann and Lösler (2017). The authors have made a step forward and consider the use of various information criteria as measures of model approximation. Among all candidate models, the one with the lowest value is selected as the best. The same idea but dedicated to the screening of multiple outliers in geodetic observations is presented by Lehmann and Lösler (2016) and Nowel et al. (2020). It should be noted that there are also other applications of combinatorics in geodesy and surveying, e.g., Baselga (2011), Biagi and Caldera (2013) or Wujanz et al. (2016).
Admittedly, the above combinatorial concept of congruence model specification overcomes the weakness of GCT, i.e., the problem of the consecutive point-by-point model specification which may suffer from displacement smearing. Nevertheless, it has other serious weakness which is absent in GCT. The best approximating model is identified here from the set of different-dimensional models and the crucial question is: How can one identify the best approximating model from such set? Of course, the goodness of model fit—represented by the model likelihood—cannot itself pinpoint the best approximating model. The higher dimension of congruence model (more points modeled as displaced) always realizes a better fit to observations, regardless of whether all the points modeled as displaced are displaced or not. To obtain some balance between the goodness of fit and the parsimony of model, the authors apply some penalties for the increase of model dimension. This means that a function of the goodness of model fit and a penalty term constitutes an identification criterion. However, there are many such criteria (ratio, several information criteria or p value) which may be used. The question is: Which one is the best and why? The answer is not unique (e.g., Burnham and Anderson 2002). The numerical experiments which were performed in the cited papers did not say much about this weakness and, generally, about the efficacy of this concept. Those experiments referred to particular scenarios of deformations and measurement errors and the obtained results may only lead to the conclusion that the suggested concept somehow works and gives results similar to GCT in those particular cases.
This paper discusses a more sophisticated combinatorial procedure to address the weaknesses associated both with the GCT procedure and the existing combinatorial methods. It is shown that, thanks to an appropriately used the possibilities of combinatorics and generalized likelihood ratio (GLR) tests (e.g., Teunissen 2006b, pp. 53–123; Koch 2010, pp. 276–295; Kargoll 2012, pp. 29–33) performed in the detection–identification–adaptation (DIA) steps (e.g., Baarda 1968; Teunissen 1990, 2018; Zaminpardaz and Teunissen 2019), both the problem of the consecutive point-by-point model specification which may suffer from displacement smearing and the problem of the comparison of different-dimensional models can be overcome, at no significant cost. At the same time, the suggested procedure allows to extended congruence analysis just as the existing combinatorial methods. In the context of GCT, the suggested procedure has rather evolutionary—than revolutionary—character and the general concepts of both procedures have similar heuristic substantiation.
The next part of the paper is organized as follows. Section 2 describes the GCT procedure and the two existing combinatorial methods. Particular attention is paid to the forms of functional congruence models and the criteria used for the identification of the best approximating model. Section 3 presents, in relation to Sect. 2, the new combinatorial procedure. A research motivation is presented in Sect. 3.1. The next sections focus on the form of dedicated functional congruence model, its parameter estimation and, most of all, the concept of the suggested combinatorial procedure. Section 4 is devoted to numerical experiments based on comprehensive computer simulations to objectively demonstrate the suggested method against conventional ones. Section 5 gives the conclusions from the study.
2 Congruence analysis using the LS approach
Appropriate matching of two (or more) control networks at mutually stable points (i.e., at the point group which has a congruent/rigid geometrical structure at all considered epochs) plays a crucial role in the congruence analysis. However, the stable points may not be a priori known. Geodetic control networks are free networks suffering from datum defect and can realize infinitely many possible matchings (Fig. 1). This is why the congruence analysis is a quite difficult task, especially when there are many displaced points (see, e.g., Hekimoglu et al. 2002, 2010; Nowel 2019).
Control networks at two epochs (inspired from Niemeier 2008, p. 437)
The suggested—in this paper—procedure of the specification of deformation congruence models is addressed to the weaknesses of the LS approach, i.e., the GCT procedure and two existing combinatorial methods. Therefore, Sect. 2 is briefly devoted to the congruence analysis using the LS approach. For the sake of simplicity and without a loss of generality, the following and further discussions are restricted to two epochs, however, nothing stands in the way to extend them to more epochs.
2.1 Global congruence test procedure
The GCT procedure can be realized based on two different approaches; both approaches, however, provide identical results. The mathematical proof of this equivalence can be found in Caspary (2000, pp. 63–64, 123–124) and some empirical results in Nowel (2018). The first approach is often named the explicit hypothesis method [two variants exist: hypotheses can be formulated before adjustment (e.g., Kok 1982; Heck 1983, pp. 158–162; Lehmann and Lösler 2017) or after adjustment (e.g., Pelzer 1971, 1974; Niemeier 1981, 2008, pp. 436–440)] and the second is named the implicit hypothesis method (e.g., Niemeier 1985; Heunecke et al. 2013, pp. 498–508). Which approach suits best depends on the software available. The first approach realized by hypotheses formulated after adjustment is best grounded in the literature and, hence, it is very briefly demonstrated below.
After screening observations for outliers, the—minimally constrained—LS adjustment of control network (e.g., using the Lagrange method) is performed for both epochs. The linear(ized) functional model of one-epoch observations has the following form:
where i is the designation of considered epoch (here i = 1, 2), \({\mathbf{A}}_{i} (n_{i} \times u_{i} )\) is the design matrix with rank deficiency \(d_{i} = u_{i} - {\text{rank}}({\mathbf{A}}_{i} )\), \({\mathbf{A}}_{i,a} = [{\mathbf{A}}_{i,s} \, {\mathbf{A}}_{i,o} ]\), \({\mathbf{x}}_{i,a} = [{\mathbf{x}}_{i,s} \, {\mathbf{x}}_{i,o} ]\) is the vector of the coordinates (or increments) of all analyzed points, i.e., all points which are common in both epochs, \({\mathbf{x}}_{i,s}\) is the vector of the coordinates (or increments) of analyzed points modeled as stable, \({\mathbf{x}}_{i,o}\) the vector of the coordinates (or increments) of remaining analyzed points, \({\mathbf{x}}_{i,n}\) is the vector of the coordinates (or increments) of non-common points (if such points exist) or nuisance parameters, \({\mathbf{y}}_{i}\) is the vector of observations (or observed minus computed values of observations), \({\mathbf{e}}_{i}\) is the vector of measurement errors and \({\mathbf{D}}_{i}^{{\text{T}}} (d_{i} \times u_{i} )\) is a datum matrix such that \({\text{rank}}({\mathbf{D}}_{i} ) = d_{i}\) and \({\text{rank}}([{\mathbf{A}}_{i}^{{\text{T}}} \, {\mathbf{D}}_{i} ]^{{\text{T}}} ) = u_{i}\). This matrix may read \({\mathbf{D}}_{i}^{{\text{T}}} = [{\mathbf{H}}_{i,a(s)}^{{\text{T}}} \;\left| {\;{\mathbf{0}}_{i,n} } \right.]\) where \({\mathbf{H}}_{i,a(s)}\) is the Helmert matrix of S-transformation filled with the Helmert components (e.g., Caspary 2000, pp. 31–33) for all the analyzed points modeled as stable and with zeros for the remaining analyzed points, and \({\mathbf{0}}_{i,n}\) is the null matrix. The above functional model is supplemented with the Gaussian probabilistic model: \({\mathbf{y}}_{i} \sim N({\mathbf{A}}_{i} {\mathbf{x}}_{i} ,\sigma_{0i}^{2} {\mathbf{Q}}_{{y_{i} }} )\) where N means the Gaussian distribution, \(\sigma_{0i}^{2}\) is the variance factor and \({\mathbf{Q}}_{{y_{i} }}\) is the cofactor matrix of observations; the probabilistic model is needed in the further part for statistical testing (3), (4). The appropriate congruence analysis is performed after the LS estimation, for the results of at least two epochs.
Since GCT resembles the iterative data snooping procedure, it can be demonstrated using the DIA framework which is often used to demonstrate that procedure (e.g., Teunissen 2006b, p. 132). Note that the DIA framework is also used for the new combinatorial procedure discussed in Sect. 3.
2.1.1 Detection
The objective is to detect whether or not the point group modeled as stable/congruent is stable/congruent at both epochs. This group may consist of all network points common to the epochs, of the points establishing the reference net or those belonging to a specified block of the network. For this point group, the null hypothesis/model: “all considered points are stable” against the most relaxed alternative hypothesis/model (no restrictions are placed on considered points): “all analyzed points are individually displaced”:
is tested using the global test of congruence (Pelzer 1971, 1974):
where \(t_{\Delta } = {{\varvec{\Delta}}}_{{\hat{x},s}}^{{\text{T}}} {\mathbf{Q}}_{{{{\varvec{\Delta}}}_{{\hat{x},s}} }}^{ + } {{\varvec{\Delta}}}_{{\hat{x},s}}\), \({{\varvec{\Delta}}}_{{\hat{x},s}}^{{}} = {\hat{\mathbf{x}}}_{2,s}^{{}} - {\hat{\mathbf{x}}}_{1,s}^{{}}\), \({\hat{\mathbf{x}}}_{i,s}\) is the vector of the estimated coordinates (or increments) of analyzed points modeled as stable at epoch i, \({\mathbf{Q}}_{{{{\varvec{\Delta}}}_{{\hat{x},s}} }}\) is the cofactor matrix of \({{\varvec{\Delta}}}_{{\hat{x},s}}^{{}}\), \(( \cdot )^{ + }\) means pseudo-inverse, \(h = {\text{rank}}({\mathbf{Q}}_{{{{\varvec{\Delta}}}_{{\hat{x},s}} }} )\), \(\sigma_{0}^{2}\) is the known variance factor (pooled for epoch 1 and 2), \(f = \infty\) and \(F_{1 - \alpha } ( \cdot )\) is the critical value from the central Fisher’s cumulative distribution function (CDF) for a chosen significance level \(\alpha\). For the estimated variance factor it goes analogously where f is the sum of the degrees of freedom from both epochs. Note the test (4) may be treated as the likelihood ratio (LR) test in today’s language (e.g., Teunissen 2006b, p. 53; Koch 2010, p. 276; Kargoll 2012, p. 29). The result of the global test of congruence (4) is either a rejection of H0, leading to an identification step, or a non-rejection in which case it may be possible to skip the identification step.
2.1.2 Identification
The rejection of H0 (3) allows to suspect the existence of displacement(s) in the point group modeled as stable and leads to the identification of the point most suspected to be displaced, i.e., the point with the greatest contribution, \(t_{\Delta }^{\max }\), to global quadratic form, \(t_{\Delta }\), (4), (e.g., with the use of Gauss‐elimination, S-transformation or implicit hypothesis method), (e.g., Caspary 2000, pp. 125–129; Niemeier 2008, pp. 446–449; Heunecke et al. 2013, pp. 500–508). One should also note that—unlike the iterative data snooping procedure—no local test is performed for these contributions in the primary version of GCT, however, nothing stands in the way of additionally doing it as some authors suggest.
2.1.3 Adaptation
The point identified as the most suspected to be displaced is treated as displaced and is transferred from the group of stable points to the group of remaining analyzed points; the vectors/matrices of model (1) are reordered. Now, the new detection step is again performed.
The above DIA testing procedure is repeated until the complete separation of stable and displaced points is accomplished, i.e., H0 (3) passes. More information about the above realization of GCT and the other realizations can be found in the papers cited in the first paragraph of Sect. 2.1.
2.2 Two existing combinatorial methods
Firstly, it is worth noting that both presented methods do not divide the network into reference and object points. However, nothing stands in the way of using these methods for the standard case, when the sub-group of potentially stable points is a priori known and only this group is modeled as stable/congruent.
The first existing combinatorial method involves statistical hypothesis tests and was presented by Velsink (2015). The deformation congruence model is here inspired by the theory of the connection of two epoch adjustment results (Teunissen 1985). The author starts from the vectors of coordinates and their covariance matrices, which may be adjusted in different datums at both epochs [e.g., different Di (2)]. To overcome the problem of different datums, the author builds the model of connection adjustment. Thanks to this, the S-transformation of the two-epoch adjustment results (the vectors of estimated coordinates and their covariance matrices) to a common datum—what is usually done to overcome the problem of different datums—can be omitted here. For all network points, the following H0: “all analyzed points are stable” against the most relaxed HA: “all analyzed points are individually displaced”:
(supplemented with the Gaussian probabilistic model) is tested by the generalized LR (GLR) test (e.g., Teunissen 2006b), where \({{\varvec{\Delta}}}_{{\hat{x},a}} = {\hat{\mathbf{x}}}_{2,a} - {\hat{\mathbf{x}}}_{1,a}\), \({\hat{\mathbf{x}}}_{1,a}\) and \({\hat{\mathbf{x}}}_{2,a}\) are the vectors of the estimated coordinates (or increments) of all analyzed points in epoch 1 and 2 [e.g., from (1)], \({\mathbf{H}}(u \times t)\) is the Helmert matrix of S-transformation (e.g., Caspary 2000, pp. 31–33; Aydin 2017), \({\mathbf{t}}(t \times 1)\) is the vector of the S-transformation parameters, t is the number of degrees of freedom between both epochs, Ht gives the vector of point displacements resulting only from different datums at both epochs (rigid body movements of adjusted networks), \({\mathbf{B}}_{I} (u \times q)\) is the specification matrix of full rank such that the redundancy of the model under HA equals zero, q is the number of the redundancy of null model and \(\nabla (q \times 1)\) is the unknown vector of the displacements of all analyzed points. The test of H0 (5), as that of H0 (3), are the global tests of congruence; no restrictions are placed on considered points under HA (here, BI is identity matrix). Thus, the first step is the same as in GCT and it may be also treated as a detection step.
The rejection of H0 (5) states that deformation is detected but not yet localized. Now, the best approximating congruence model is identified. To identify such a model the author considers by “brute force” a large amount, m, of different-dimensional congruence models; each defined/conditioned by its specification matrix of hypothetical deformations, \({\mathbf{B}}_{I,i} (u \times q_{i} )\). It means, one point may be displaced, but also two, three, or even more points may be displaced (either in the same or different way). Each of such considered cases is expressed by the following alternative hypothesis:
where \({\mathbf{B}}_{I,i}\) and \(\nabla_{i}\) refer to a candidate congruence model. The H0 (5) against the alternative hypothesis (6) leads to the GLR test statistic of a candidate model:
where \(p_{N}\) is the Gaussian probability density function (PDF) of \({{\varvec{\Delta}}}_{{\hat{x},a}}\). Finally, the ratio of the appropriate test statistic \((T_{i} /q_{i} ) \sim F\left( {q_{i} ,\infty ,0} \right)\) and their critical value:
gives an indication of the approximation of candidate model, where \(F_{{1 - \alpha_{i} }} ( \cdot )\) means the critical value from the central Fisher’s CDF, \(q_{i} = {\text{rank}}({\mathbf{B}}_{I,i} )\) and \(\alpha_{i}\) is the significance level which realizes the same non-centrality parameter and power as the global test of H0 (5), see also Baarda (1968). The model with the largest value of this ratio is identified as the best approximating congruence model. More information, e.g., the examples of specification matrices, \({\mathbf{B}}_{I,i}\), can be found in Velsink (2015). Note also that the general problem of the penalty for the complexity of a candidate model is discussed in detail by Burnham and Anderson (2002).
The second existing combinatorial method, presented by Lehmann and Lösler (2017), is based on information theory and uses some relative estimators of the Kullback–Leibler divergence as measures of the approximation of a candidate congruence model. The functional model is formed as an ordinary model of observations (not of adjusted coordinates) with the same datum at all considered epochs and the problem of different datums is not considered here. The authors do not carry out the detection of deformation (global test of congruence) and start the identification of the best approximating congruence model straight away. Similar to Velsink (2015), they consider by “brute force” many different-dimensional candidate models; each defined/conditioned by its specification matrix of hypothetical deformations, \({\mathbf{B}}_{II,i}\). The functional linear(ized) model employed here has the following form (two epochs):
where \(i = 1,...,m + 1\), \({\mathbf{x}}_{1,a}\) and \( \, {\mathbf{x}}_{2,a}\) are the vectors of coordinates (or increments) of all analyzed points in epoch 1 and 2, \({\mathbf{x}}_{1,n}\) and \({\mathbf{x}}_{2,n}\) are the vectors of coordinates (or increments) of non-common points (if such points exist) in epoch 1 and 2, \({\mathbf{B}}_{II,i}^{{\text{T}}} (n_{c} \times u) = \left[ {{\mathbf{B}}_{{a_{i} }}^{{\text{T}}} \;\left| {\begin{array}{*{20}c} {\;{\mathbf{0}}_{1,n} } & {{\mathbf{0}}_{2,n} } \\ \end{array} } \right.} \right]\), \(n_{c}\) is the number of constraints and \({\mathbf{B}}_{{a_{i} }}\) is a specification matrix referring to all analyzed points. The other denotations are the same as earlier. For example, if a leveling control network including twelve analyzed points (without non-common points), (e.g., the network which will be depicted in Fig. 2) was processed where points 1–5 were be modeled as displaced—in consequence, the remaining analyzed points (6–12) as stable—the specification matrix (10) would read:
where I is identity matrix. Additionally, if points 1, 3 were modeled as shifted with a certain joint value and points 2, 4 as shifted with another joint value, this matrix would read:
and if, for example, all conditioned displaced points (1–4) were modeled as shifted with a certain joint value, it would be enough to introduce an additional condition into matrix (13), e.g., it could be the condition:
which would model points 1, 2 as shifted with a certain joint value. The matrix \({\mathbf{D}}_{i}^{{\text{T}}} (d \times u)\) is a datum matrix such that \({\text{rank}}({\mathbf{D}}_{i} ) = d\), \(d = u - {\text{rank}}([{\mathbf{A}}^{{\text{T}}} \, {\mathbf{B}}_{II,i} ]^{{\text{T}}} )\) and \({\text{rank}}([{\mathbf{A}}^{{\text{T}}} \, {\mathbf{B}}_{II,i} \, {\mathbf{D}}_{i} ]^{{\text{T}}} ) = u\). This matrix is not unique and it may read:
where \({\mathbf{H}}_{a(s)}\) is filled with the Helmert components only for the analyzed points modeled as stable and with zeros for the remaining analyzed points. Note that the probabilistic model, e.g., the Gaussian model, is also necessary for the presented method. After the maximum likelihood estimation, which is reduced here to the LS estimation under the Gaussian probabilistic model, the authors use—unlike Velsink (2015)—information criteria to identify the best approximating congruence model from the set of all candidate models (9), (10), i.e., from the set containing null model and all alternative models which together form m + 1 candidate models (9), (10). For instance, the Akaike’s information criterion with correction (AICc)—suggested by Lehmann and Lösler (2017)—for a candidate model (9), (10) reads (see also Burnham and Anderson 2002, pp. 66, 374–380):
where \(k = u - r - d\), \(r = {\text{rank}}({\mathbf{B}}_{II,i}^{{\text{T}}} {\mathbf{Q}}_{{\hat{x}}} {\mathbf{B}}_{II,i} )\) and \({\mathbf{Q}}_{{\hat{x}}}\) is the cofactor matrix of externally estimated coordinates. If the variance factor is unknown (estimated), k is increased by one because the variance factor also counts as a parameter. Among all candidate models, the one with the lowest value of information criterion [e.g., AICc (16)] is identified as the one best approximating. One can see the described information criterion includes a trade-off between the goodness of fit and the dimension of model. Thus, such a measure (16) plays the same role as the ratio used by Velsink (2015), [see Eq. (8)]. More information, e.g., the ready-made form of AICc which is used in congruence analysis can be found in Lehmann and Lösler (2017).
Configuration of the leveling control network (based on Lehmann and Lösler 2017)
Summing up, the concepts of the above combinatorial methods are similar. They differ mainly in the criterion (8), (16) using to the identification of the best approximating model.
3 Combinatorial iterative DIA testing procedure
This section is devoted to a new LS combinatorial method which—unlike the existing ones—involves a combinatorial iterative DIA (CIDIA) testing procedure. The tests involving in the suggested procedure are based on residuals and, as such, are invariant under a change of S-system (reference system in which parameters are defined), as it is postulated in the literature (e.g., Koch 1985; Velsink 2015).
3.1 Motivation
As mentioned in Sect. 1, the GCT procedure has a serious weakness in the case of multiple displacements. It involves the consecutive point-by-point model specification which is conditioned by the results of previous iterations. It means, models considered in the second iteration are conditioned by the result of the first iteration (i.e., the point which is identified as most suspected to be displaced in the first iteration is not modeled as stable in the second iterations). Then, the models considered in the third iteration are conditioned by the results of the first and the second iterations (i.e., the point which is identified as most suspected to be displaced in the first iteration and the point which is identified as most suspected to be displaced in the second iteration are not modeled as stable in the third iterations), and so on. The problem is that the consecutive identifications of single-point displacements may be burdened by the smearing of displacements which mainly results from the possible displacements of points modeled as stable (in the congruent part of network—reference base) and slightly from network interrelations (e.g., Hekimoglu et al. 2002, 2010; Erdogan and Hekimoglu 2014). In consequence, a stable point may be identified as most suspected to be displaced in a given iteration. If it is the case, the models which are specified in the next iterations are—from now—conditioned by the result of erroneous identification from the previous iteration.
To avoid the above problem, Velsink (2015) and Lehmann and Lösler (2017) discuss a revolutionary—in the context of GCT—combinatorial concept. The authors suggest the identification of the best approximating model directly from the set of all possible candidate models. Thanks to this, the above problem of the consecutive specification which may suffer from displacement smearing can be avoided. Unfortunately, another problem arises, i.e., the problem of the comparison of different-dimensional models. To identify the best approximating model, the authors use the penalized goodness of model fit and the model with an optimal value is identified as the one best approximating. The goodness of model fit itself can be clearly expressed; it is defined as log-likelihood function and is characterized by the function of residuals in practice. However, a serious problem is related to the term which penalizes an increase in the dimension (number of parameters) of model. It means the identified model should not be either underfitting to the observations (lower dimension of model than it is) or overfitted to the observations (higher dimension of model than it is), (see also Burnham and Anderson 2002). Since the candidate congruence models are different-dimensional, it is not exactly clear how one should properly penalize the goodness of fit for an increase in the dimension of model. To penalize it, Velsink (2015) divides the goodness of fit by corresponding critical value, which is the function of the dimension of candidate congruence model (8). This idea of penalizing comes from the paper de Heus et al. (1994, p. 266) and it has an only intuitive justification. Indeed, the penalty terms suggested by Lehmann and Lösler (2017) are additionally justified by information theory, nonetheless, this is a quite arbitrary area. It means, there are suggested various penalty terms, but it is not exactly clear which one provides the best trade-off between the goodness of model fit and its parameter number. In the author’s opinion, the above problem of the comparison of different-dimensional models may raise some fears—especially in the context of surveying and geodetic observation models—and it still deserves deeper theoretical and empirical research.
Furthermore, although the power of computers is growing, the existing combinatorial concept can still suffer from a computational load. For instance, Velsink (2015, p. 1086) states: “With an adjustment model that covers more than two epochs, it is possible to specify alternative hypotheses that describe deformation processes during more epochs. The search for the best alternative hypothesis will be even more complicated. It is therefore important to find ways to reduce the amount of alternative hypotheses that have to be tested.” What is worse, this problem can also concern the two-epoch analysis of large and very large control networks, where hundreds of thousands—or even more—of time-consuming models would have to be considered and it also would not allow obtaining calculation results in near-real time, not to mention different simulation studies in scientific research, where one performs the thousands of calculation repetitions for a given method. It becomes clear that the computational load can also be a certain problem concerning the existing combinatorial methods.
The further part of Sect. 3 is devoted to the new combinatorial procedure which is addressed to the above problems. In the context of GCT, the new combinatorial procedure has rather evolutionary—than revolutionary, as the existing combinatorial methods—character.
3.2 General form of mathematical model
As it is known, the congruence models either in explicit [e.g., (9), (10)] or implicit forms may be equivalently used in congruence analysis. For the suggested procedure, the congruence model will be designed in implicit form, however, a little different than the conventional implicit one (e.g., Niemeier 1985; Caspary 2000, pp. 123–125; Heunecke et al. 2013, pp. 498–508). The idea behind the suggested form of congruence model comes from Teunissen (2006b, pp. 71–123), where the general problem of misspecifications (e.g., observation blunders) in an underlying functional model of random sample is discussed. Since point displacements may be interpreted as misspecifications in underlying congruence model (i.e., in the model which specifies all analyzed points as congruent/stable) and the results of deformation measurements constitute a random sample, the general model discussed by Teunissen (2006b, p. 72) may be implemented here as an elementary deformation congruence model. Furthermore, such implicit congruence model [first published in Nowel (2018)] is more transparent and easier to implement (ibid), than the implicit model presented in the literature, hence it motivates the transfer of it to the suggested CIDIA procedure.
For the sake of simplicity and without a loss of generality, the discussed congruence model—as the congruence models discussed earlier (Sect. 2)—does not divide the network into reference and object points, however, nothing prevents using it.
3.2.1 Elementary congruence model
In relation to the general form of functional model discussed in Teunissen (2006b, pp. 71–123), one can formulate the deformation congruence functional model in the following implicit form (two epochs):
where \({\mathbf{x}}_{a}\) is the vector of coordinates (or increments) of all analyzed points (the same values for both epochs), \({\mathbf{E}}_{i} (n \times q_{i} )\) may be interpreted as an elementary specification matrix, \({\mathbf{A}}_{2,o}^{i}\) is the sub-matrix of \({\mathbf{A}}_{2,a}\) referring to analyzed points modeled as displaced and \(\nabla_{i} (q_{i} \times 1)\) is the unknown vector of point displacements. The other denotations are the same as earlier. Note that the general form of model (17), i.e., \({\mathbf{y}} = {\mathbf{Ax}} + {\mathbf{E}}_{i} \nabla_{i} + {\mathbf{e}}\), comes from Teunissen (2006b, p. 72). The part \({\mathbf{y}} = {\mathbf{Ax}} + {\mathbf{e}}\) refers here to the underlying congruence model (without point displacements) and the possible displaced points in this congruent part are specified by \({\mathbf{E}}_{i}\). The matrix \({\mathbf{D}}_{i}^{{\text{T}}} \left( {d \times (u + q_{i} )} \right)\) is a datum matrix such that \({\text{rank}}({\mathbf{D}}_{i} ) = d\), \(d = (u + q_{i} ) - {\text{rank}}([{\mathbf{A}} \, {\mathbf{E}}_{i} ])\) and \({\text{rank([}}[{\mathbf{A}} \, {\mathbf{E}}_{i} ]^{{\text{T}}} \, {\mathbf{D}}_{i} ]^{{\text{T}}} ) = u + q_{i}\). This matrix is not unique and it may read:
where \({\mathbf{H}}_{a(s)}\) is filled with Helmert components only for the analyzed points modeled as stable and with zeros for the remaining analyzed points. For example, if a leveling control network including twelve analyzed points was processed where points 6–12 was modeled as stable (and, in consequence, points 1–5 as displaced), the Helmert matrix would read: \({\mathbf{H}}_{a(s)}^{{\text{T}}}\) = [0 0 0 0 0 1 1 1 1 1 1 1]. Note that for the most relaxed model, i.e., where all analyzed points are specified as individually displaced and, in consequence, no stable points exist, \({\mathbf{E}}(n \times q) = [{\mathbf{0}} \, {\mathbf{A}}_{2,a}^{{\text{T}}} ]^{{\text{T}}}\), the datum matrix may read:
where \(d_{E} = (u + q) - {\text{rank}}([{\mathbf{A}} \, {\mathbf{E}}])\) and \({\mathbf{H}}_{a}\) is filled with the Helmert components for all analyzed points.
One can see that—unlike the explicit model (9), (10)—no constraints specifying stable and displaced points are needed because they are already implied in implicit model (17). For example, if a leveling control network was processed where points 1–5 were modeled as displaced, the sub-matrix \({\mathbf{A}}_{2,o}^{i}\) would only include the columns of matrix \({\mathbf{A}}_{2,a}\) which refer to the points 1–5. Note that the suggested implicit congruence model (17) is equivalent to the conventional implicit congruence model and, simultaneously, to the conventional explicit congruence model (9), (10). Experimental confirmation of this equivalence can be found in Nowel (2018) or in the further part of this paper (Sect. 4.3, Table 5—the last two column).
Let us assume that the probabilistic model is Gaussian, as earlier.
3.2.2 Extending part
Admittedly, the implicit model (17) makes it possible to set up displaced points. Nevertheless, these points can only be modeled there as displaced in different ways (individual displacements) and this model may be only used for such elementary deformation scenarios. However, this can be easily developed in extended scenarios which are considered in the existing combinatorial methods. The possible extended deformation scenarios, e.g., shifting away with the same value in the same direction, can be here easily specified using the constraints in the following general form:
where \({\mathbf{F}}\left( {n_{c} \times (u + q_{i} )} \right)\) may be interpreted as extended specification matrix, \({\mathbf{F}}_{\nabla } \left( {n_{c} \times q_{i} } \right)\) is the sub-matrix which constitutes some—linearly independent—constraints on the point displacements and \(n_{c}\) is the number of constraints. For example, if a leveling control network including twelve analyzed points (without non-common points) was processed where points 1–5 were modeled as displaced and, additionally, points 1, 3 as shifted with a certain joint value (\(\nabla_{1} - \nabla_{3} = 0\)) and points 2, 4 as shifted with another joint value (\(\nabla_{2} - \nabla_{4} = 0\)), this matrix would read:
and if, for instance, all conditioned displaced points (1–4) were additionally specified as shifted with joint value, it would be enough to introduce an additional condition, e.g., the condition: \(\nabla_{1} - \nabla_{2} = 0\) which would generate the third row:
into matrix (22). The reader is also referred to the corresponding sub-matrices (13) and (14) of the matrix \({\mathbf{B}}_{II,i}^{T}\) which specify there the same—as above—extended deformation scenarios. One can see that the constraints (22), (23) are easier to formulate and more transparent.
3.3 Estimation
Let us assume that the deformation constraints (21) are not applied and the functional model has the elementary form (17). Keep in mind that the probabilistic model of (17) is Gaussian, the maximum likelihood optimization problem leads to the LS optimization problem. However, since the estimator of parameter vector has to satisfy also the datum constraints (18), the ordinary LS optimization criterion has to be replaced by the following Lagrange optimization criterion:
where \({\mathbf{e}}_{{}} = {\mathbf{y}} - {\mathbf{A}}_{{E_{i} }} {\mathbf{x}}_{{\nabla_{i} }}\), \({\mathbf{A}}_{{E_{i} }} = [{\mathbf{A}} \, {\mathbf{E}}_{i} ]\) and \({\mathbf{k}}_{i}\) is the vector of Lagrange multipliers. The necessary condition for the existence of a minimum of this function has the form: \(\varphi_{L} \left( {{\mathbf{x}}_{{\nabla_{i} }} } \right)/\partial {\mathbf{x}}_{{\nabla_{i} }} = {\mathbf{0}}|_{{{\mathbf{x}}_{{\nabla_{i} }} = {\hat{\mathbf{x}}}_{{\nabla_{i} }} }}\) (e.g., Koch 2010, pp. 170–173) and it leads to the estimator:
where \({\mathbf{Q}}_{{{\hat{\mathbf{x}}}_{{\nabla_{i} }} }}^{(0)} = {\mathbf{N}}^{ - }\), \({\mathbf{N}} = {\mathbf{A}}_{{E{}_{i}}}^{{\text{T}}} {\mathbf{Q}}_{y}^{ - 1} {\mathbf{A}}_{{E_{i} }}\) and \(( \cdot )^{ - }\) means g-inverse. For instance, an easy to use g-inverse of N may be pseudo-inverse \({\mathbf{N}}^{ - } \to {\mathbf{N}}^{ + } = ({\mathbf{N}} + {\mathbf{GG}}^{{\text{T}}} ) - {\mathbf{GG}}^{{\text{T}}}\), where \({\mathbf{G}}((u + q_{i} ) \times d)\) is the orthonormal matrix which includes d eigenvectors from the modal matrix of N which correspond to the d zero-eigenvalues from the spectral matrix of N and \({\mathbf{NG}} = {\mathbf{0}}\), \({\mathbf{AG}} = {\mathbf{0}}\) (e.g., Niemeier 2008, pp. 236–237, 262–264). Note the modal and spectral matrices can be easily obtained, e.g., by the MATLAB software with the use of function eig. Since estimator (25) has to satisfy also datum constraints (18), one may additionally write:
Finally, after the substitution of \({\mathbf{k}}_{i}\) (26) into (25), the parameter estimator of elementary congruence model (17), (18) and its cofactor matrix take the form:
where \({\mathbf{I}}\left( {(u + q_{i} ) \times (u + q_{i} )} \right)\) is identity matrix, \({\hat{\mathbf{x}}}_{{\nabla_{i} }}^{(0)}\) is the ordinary LS estimator of former—without datum constraints (18)—model (17) and \({\mathbf{Q}}_{{{\hat{\mathbf{x}}}_{{\nabla_{i} }} }}^{(0)}\) is its cofactor matrix. It can be easily proved that the corresponding residuals and their cofactor matrix take here the form:
One can see that the estimator (27) is found from updating the estimator of former model (17), \({\hat{\mathbf{x}}}_{{\nabla_{i} }}^{(0)}\), in a formulation based on latter constraints (18). From the Lagrange solution (27) it also follows that a parameter estimator updating—in a formulation based on some linearly independent constraints—requires only the former parameter estimator, its cofactor matrix (here: \({\hat{\mathbf{x}}}_{{\nabla_{i} }}^{(0)}\), \({\mathbf{Q}}_{{{\hat{\mathbf{x}}}_{{\nabla_{i} }} }}^{(0)}\)) and the latter coefficient matrix of constraints (here: \({\mathbf{D}}_{i}\)).
Now, let us assume that the deformation conditions (21) are applied and the functional model has the extended form (17), (18), (21). According to the updating principle disclosed above in the Lagrange solution, the parameter estimator of the extended model and its cofactor matrix can be derived with the use of the former—now from elementary model (17), (18)—parameter estimator, \({\hat{\mathbf{x}}}_{{\nabla_{i} }}^{(1)}\) (27), its cofactor matrix, \({\mathbf{Q}}_{{{\hat{\mathbf{x}}}_{{\nabla_{i} }} }}^{(1)}\) (28), and the coefficient matrix of latter constraints, F (21), and they take the form:
The corresponding residuals and their cofactor matrix take the form:
Note that all the above LS estimators can be also derived based on the theory of recursive LS estimation; in fact, they are particular forms of recursive LS estimators (e.g., Teunissen 2006a, pp. 32–36).
3.4 Testing procedure
Most generally, the core idea behind the suggested procedure will be to determine the final congruence model using consecutive increasing the dimension of deformations where: Iteration I will refer to all possible—to occur—candidate models with \(q_{I}\)-dimensional deformations and \(q_{I}\) is the lowest possible dimension of deformations (the dimension of single-point displacement vector: 1, 2 or 3 for the 1D, 2D or 3D network), Iteration II will refer to all possible candidate models with \(2q_{I}\)-dimensional deformations, etc. What is most important, unlike GCT, the deformation scenarios specified in subsequent iteration steps will not be conditioned by the results of previous iteration steps. Thanks to this, if a stable point was—erroneously—specified as displaced in a given iteration step, it would be still able to be specified as stable in the next iteration step(s). The only information which comes from the previous iteration is: “the dimension of deformations considered in the previous iteration is too small”. In the context of the GCT procedure, the general concept of the suggested procedure is based on consecutive dimension-by-dimension specification, instead of consecutive point-by-point.
After screening observations for outliers (in each epoch), e.g., by an iterative data snooping procedure, the suggested CIDIA procedure goes as follows. Note that, however, the sets of screened observations are not fully Gaussian. In consequence, the global test statistics involving in the suggested procedure—as those involving in GCT—are not fully the realizations of the Fisher PDF and imposed testing parameters are not strictly realized.
Detection I Since it first should be known whether one can have any confidence in the following null hypothesis: “all analyzed points are stable”, the first step may consist of a check on the global validity of the H0. This implies that the H0 has to be opposed to the most relaxed—in deformation context—alternative hypothesis possible. The most relaxed HA is the one that leaves the point displacements completely free and it states that “all analyzed points are individually displaced”. For the suggested functional model (17) there is the situation:
where \({\mathbf{E}}(n \times q) = [{\mathbf{0}} \, {\mathbf{A}}_{2,a}^{{\text{T}}} ]^{{\text{T}}}\) is the most relaxed specification matrix and \(\nabla (q \times 1) \in {\mathbb{R}}^{q} /\{ 0\}\) is the unknown vector of the displacements of all analyzed points. According to Teunissen (2006b, p. 53), the appropriate GLR test can be computed from the PDF of y under H0 and HA. However, since the same Gaussian PDF, \(p_{N}\), is assigned for the null and alternative model (35), only the function proportional to the Gaussian PDF is sufficient here and the constant term, \(c_{N} = (2\pi )^{ - n/2} |\sigma_{0}^{2} {\mathbf{Q}}_{y} |^{ - 1/2}\), may be omitted for the GLR test. The function proportional to the Gaussian PDF of y under H0 reads:
and under HA it reads analogously, where \(\propto\) means ‘is proportional to’, \({\mathbf{e}}_{0} = {\mathbf{y}} - {\mathbf{Ax}}\) and \({\mathbf{e}}_{A} = {\mathbf{y}} - {\mathbf{Ax}} - {\mathbf{E}}\nabla\). In consequence, the GLR test statistic takes the form:
where \({\hat{\mathbf{e}}}_{0} , \, {\hat{\mathbf{e}}}_{A}\) are the vectors of LS residuals under \(H_{0}\), \(H_{A}\) and \(h_{I} = {\text{rank}}([{\mathbf{A}} \, {\mathbf{E}}]) - {\text{rank}}({\mathbf{A}})\) is the degrees of freedom of the nominator of test statistic. Finally, the appropriate GLR test for testing H0 against HA (35) can be written as:
for the known variance factor, \(\sigma_{0}^{2}\), and
for the unknown variance factor, \(\hat{\sigma }_{0}^{2} = {\hat{\mathbf{e}}}_{A}^{{\text{T}}} {\mathbf{Q}}_{y}^{ - 1} {\hat{\mathbf{e}}}_{A} /f\); where \(f = n - {\text{rank}}\left( {[{\mathbf{A}} \, {\mathbf{E}}]} \right)\) and \(F_{{1 - \alpha_{I} }} ( \cdot )\) is the critical value from the central Fisher’s CDF for a chosen significance level, \(\alpha_{I} = \alpha\). For the known variance factor the significance level means:
and analogously for the unknown variance factor (\(\infty : = f\)), where \(c_{F}\) means the Fisher’s CDF. A useful discussion about the choice of \(\alpha\) in the context of geodetic deformation analysis is given by Koch (2010, p. 281). Note that the GLR test (38) [or (39)] is equivalent to the global test used in the GCT procedure [see Eq. (4)]. The global test (4) is associated with the explicit hypothesis of congruence (3) and the global test (38) [or (39)] is associated with implicit hypothesis of congruence (35) and, this is why, the forms of both test are different.
If the global test (38) [or (39)] fails, this is a reason to suspect that the H0 is wrong and it may be rejected.
Identification I The failure of the global test allows to suspect the existence of at least \(q_{I}\)-dimensional deformations and it may lead to the identification of the best fitting model from all possible models with \(q_{I}\)-dimensional deformation, where \(q_{I}\) means the lowest possible dimension of deformations (recall that \(q_{I}\) = 1, 2 or 3 for the 1D, 2D or 3D network, respectively, is imposed in this work). Two scenarios of such deformations may be considered in geodesy and surveying. The first—elementary—scenario of \(q_{I}\)-dimensional deformations includes all possible models with one displaced point and the candidate models can be specified as follows:
where \({\mathbf{E}}_{I,i} (n \times q_{I} )\) is the elementary specification matrix referring to a candidate model i and \(w\) is the number of all analyzed points. For each model (41), one may assign the following log-likelihood function:
Since the constant terms \(\ln c_{N}\), \(- \,0.5\) and \(\sigma_{0}^{2}\) are the same for each candidate model, these terms may be omitted and, finally, the following quadratic form:
may be assigned for each candidate model (41) as a measure of the goodness of model fit [see also (29)]. One can see that the lower the value of the quadratic form (43) is, the higher the value of log-likelihood function (42) and, simultaneously, the better fit of the i model to observations is [see also (8), (16)]. The second—extended—scenario of \(q_{I}\)-dimensional deformations can include all possible models with two, three, …, \(s_{\max }\) points displaced with a joint value (block shift) and the models can be specified as follows:
where elementary specification matrix \({\mathbf{E}}_{I,j}\) refers to the candidate model of j deformation scenario, \({\mathbf{F}}_{I,j}\) constitutes constraints (21) which specify all the point displacements of j deformation scenario as shifted with a joint value [the constraints of type (22), (23)] and
where \(s_{\max }\) is a chosen maximal number of points modeled as displaced. According to (42), j candidate model (44) produces the following quadratic form:
which corresponds to the goodness of j model (44) fit [see also (33)]. In consequence of the above, the set of all possible models with \(q_{I}\)-dimensional deformations, i.e., (41), (44), with their quadratic forms, i.e., \(t_{I,i}\), \(t_{I,j}\) can be obtained. Recall that the lower value of the quadratic form of a candidate model is, the better fit of the model to observations is. Thus, the candidate model which returns the minimum quadratic form, i.e., \(t_{I}^{{\min}} = \mathop {\min }\limits_{i,j} \{ t_{I,i} ,t_{I,j} \}\), may be treated as the best fitting model with \(q_{I}\)-dimensional deformations, \(H_{I}^{{{\text{best}}}}\). Note that the constant-dimensional (not different-dimensional) models are compared and no penalty term is needed.
Remark The suggested CIDIA procedure is mainly addressed to elementary deformation scenarios and the extending part (21) should be treated only as optional. The extended model (17), (18), (21) is worth to be, additionally, considered only for some rigid movement phenomena expected with a high probability, e.g., the expected rigid movement of Earth’s crust, glacier or shelf ice. However, for the deformation phenomena where rigid movements are not expected as, e.g., slope creep, ground subsidence or various deformations of man-made structures, only elementary congruence models (17), (18) should be considered. It is because it may happen that an extended congruence model may turn out the best fitting model in the case of elementary deformations. If it was the case, it would be of course disinformation.
Adaptation I The best fitting model with \(q_{I}\)-dimensional deformations may be treated as currently the best model and it may replace the current H0. It means, the identified model becomes the new null model:\(H_{0} : = H_{I}^{{{\text{best}}}}\). Once the adaptation step is completed, one still has to check if the newly specified model is valid. This implies a repetition of the detection step.
Note if the above Iteration I considered only the elementary deformation scenarios (41), the same new H0 as in the first iteration in GCT would be obtained. If extended deformation scenarios (44) were also considered, it could not be the case.
Detection II The null hypothesis has been just actualized and the global test again may be performed to validate it. This implies the following H0: “the best fitting \(q_{I}\)-dimensional deformation model is valid” against the most relaxed HA possible: “all analyzed points are individually displaced”:
The appropriate GLR test for testing H0 against HA (47) can be written as:
for the known variance factor and analogously \(T_{{II,\hat{\sigma }}} \le F_{{1 - \alpha_{II} }} (h_{II} ,f,0)\) for the unknown variance factor; where \(h_{II} = h_{I} - q_{I}\) and \(\alpha_{II}\) is the current significance level which realizes the same non-centrality parameter and power as the initial global test (38) [or (39)], (Baarda 1968). For the known variance factor, this coupling may be symbolically written as:
and analogously for the unknown variance factor (\(\infty : = f\)), where \(\gamma_{I}\) is a given power of initial global test. The choice of equal values for the non-centrality parameter \(\lambda_{I} = \lambda_{II}\) and power \(\gamma_{I} = \gamma_{II}\) in both tests implies that certain model misspecification(s) [here, point displacement(s)] can be detect with the same probability by both tests. In practice, for a given power, the procedure (49) may involve computing the non-centrality parameter, \(\lambda_{I}\), from the relation (e.g., Koch 2010, pp. 279–280):
and then \(\alpha_{II}\) (associated with \(F_{{1 - \alpha_{II} }}\)) from the analogous relation:
For the unknown variance factor, the above procedure goes analogously (\(\infty : = f\)). The algorithms presented by Gaida and Koch (1985), Aydin and Demirel (2005) or the ready-made tables given by Aydin (2012) can be used to compute \(\lambda_{I}\) and then \(\alpha_{II}\). It can be also easily done with the use of function ncfcdf (corresponding to \(c_{F}\)) from the MATLAB Statistics Toolbox. As can be seen, the suggested concept of the coupling of different-dimensional tests (49)—using the same non-centrality parameters for a given test power—comes from the B-method (Baarda 1968). Note that, however, the consecutive detection tests—in different iteration steps—are coupled here, while the detection and identification tests—in the same iteration steps—are coupled in the B-method. A more detailed discussion devoted to the concepts of the coupling of different-dimensional tests can be found, e.g., in Arnold (1981, e.g., pp. 128–137, 182–187, 195–196).
Identification II The failure of the global test allows to suspect the existence of at least \(q_{II}\)-dimensional deformations where \(q_{II} = 2q_{I}\) and it may lead to the identification of the best fitting model from all possible models with \(q_{II}\)-dimensional deformations. The elementary scenario of \(q_{II}\)-dimensional deformations includes all possible models with two displaced points and the models can be specified as follows:
where \({\mathbf{E}}_{II,i} (n \times q_{II} )\) is elementary specification matrix referring to a candidate model i and
is the number of all possible models with \(q_{II}\)-dimensional elementary deformations; where s (currently, \(s = 2\)) is the number of displaced points. Note that—unlike GCT—the result of Iteration I does not in any way restrict the current set of candidate models (52). It means, if a stable point was—erroneously—specified as displaced in Iteration I, it would be still able to be specified as stable. According to (42), i candidate model (52) produces the following quadratic form:
which may be assigned for each candidate model (52) as a measure of the goodness of model fit. The candidate model which returns the minimum quadratic form, i.e., \(t_{II}^{{\min}} = \mathop {\min }\limits_{i} \{ t_{II,i} \}\), may be treated as the best fitting model with \(q_{II}\)-dimensional elementary deformations, \(H_{II}^{{{\text{best}}}}\).
Remark Theoretically, the extended deformation scenarios can be also considered for the \(q_{II}\)-dimensional deformations, i.e., some points can be specified as shifted with a certain joint value and some other point(s) as shifted with another (joint) value. For higher-dimensional deformations, more such extended scenarios can be specified. Nevertheless, the number of all possible models violently increases for such higher-dimensional extended scenarios, in relation to \(q_{I}\)-dimensional extended scenario. In consequence, they are not considered in this paper.
Adaptation II As before, the best fitting model with \(q_{II}\)-dimensional deformations may be treated as currently the best model and, hence, it may replace the current null hypothesis: \(H_{0} : = H_{II}^{{{\text{best}}}}\). It may be tested once again whether the newly specified H0 is valid. This, of course, implies the repetition of the detection step.
Detection, identification, adaptation III, …: If the currently valid H0 is again rejected by the global test (now, \(\alpha_{III}\) comes from \(\lambda_{I} (\alpha_{I} ,h_{I} ,\gamma_{I} ) = \lambda_{III} (\alpha_{III} ,h_{III} ,\gamma_{III} = \gamma_{I} )\) and \(h_{III} = h_{I} - q_{II}\)), one may suspect that at least \(q_{III}\)-dimensional deformations may exist in the group of all, w, analyzed points, where \(q_{III} = 3q_{I}\). This may lead to the new identification of the best fitting model from the set of all possible \(c_{III}\) (53) elementary models with three displaced points (currently, \(s = 3\)). Note that—unlike GCT—the results of Iterations I, II do not in any way restrict the current set of candidate models. Each of such elementary models can be expressed by the alternative hypothesis of type (52) through the new elementary matrices \({\mathbf{E}}_{III,i} (n \times q_{III} )\) specifying three displaced points. As before, the obtained \(c_{III}\) alternative models constitute the set of quadratic forms \(t_{III,i}\) of type (54) and the best fitting \(q_{III}\)-dimensional elementary model may be identified, \(H_{III}^{{{\text{best}}}}\) (the model which returns \(t_{III}^{{\min}} = \mathop {\min }\limits_{i} \{ t_{III,i} \}\)). The identified model replaces the currently valid null hypothesis: \(H_{0} : = H_{III}^{{{\text{best}}}}\) and, then, the detection step may be repeated.
The above combinatorial iterative DIA testing procedure should be carried out until a null hypothesis is accepted in the detection step. At this moment, one may assume that the deformation congruence model is finally specified.
3.5 Quality analysis
To analytically describe the quality of the global test in a given detection step, the reliability can be measured by the power of the test for a user-defined displacement vector, \(\nabla\). The necessary non-centrality parameter of the distribution of both \(T_{\sigma } \sim F(h,\infty ,\lambda )\) and \(T_{{\hat{\sigma }}} \sim F(h,f,\lambda )\) under an alternative hypothesis reads (see, e.g., Teunissen 2006b, pp. 74–77):
where \({\mathbf{Q}}_{{\hat{e}_{0} }}\) is the cofactor matrix of residuals under a null hypothesis [see (30), (34)]. The non-centrality parameter (55) is valid either for elementary or extended deformations; a possible extended information would be included in \(\nabla\). The corresponding power of analyzed test (the probability that the test will detect a user-defined displacement vector) can be calculated based on relation (50), e.g., with the use of function ncfcdf (corresponding to \(c_{F}\)) from the MATLAB Statistics Toolbox. Note that if outliers were screened, the actual test power would be higher than the one calculated from (50).
Admittedly, the above reliability analysis can give the power of a given global test (or a lower bound of this power) and it can be somehow useful. Nevertheless, it does not solve the problem of the quality analysis of estimated parameters in the finally specified model. The analytical analysis of such quality would require a more advanced study. It would have to refer to the probabilistic properties of DIA estimators (given by Teunissen 2018) which, however, is outside the scope of this paper.
Nevertheless, a useful empirical quality analysis can be performed for a user-defined displacement vector, \(\nabla\). For example, the efficacy of correct model specification can be measured by the mean success rate (MSR) with the use of stochastic computer simulations and such empirical quality analysis, if necessary, is recommended in this paper. Note that this simulation technique has been widely applied in geodesy and surveying since the pioneering work of Hekimoglu and Koch (1999, 2000).
4 Numerical experiments
To demonstrate the suggested CIDIA procedure against the existing methods, four experiments are presented where the control networks considered by Velsink (2015) and Lehmann and Lösler (2017) are used. Both the deformation scenarios given by the authors and generated using comprehensive computer simulations are analyzed. Outliers in computer-simulated observations are screened using the iterative data snooping procedure without detection steps. The identification steps come from the Baarda’s data snooping and the significance level of the w-tests is set to \(\alpha_{1} = 0.001\). Deformation congruence models are specified with the use of: the GCT procedure, the method presented by Velsink (2015) denoted as the V method, the method presented by Lehmann and Lösler (2017) denoted as the AICc method and the suggested CIDIA procedure.
4.1 Experiment 1: arbitrary scenarios 1D
Experiment 1 is transparent and can be easily repeated by the reader. A simple leveling control network discussed by Lehmann and Lösler (2017) is engaged in it (Fig. 2). The network consists of 12 analyzed points (without non-common points) and is “measured” in two epochs with the same known standard deviation \(\sigma = 1\) mm; the values of the true heights of the network points are unimportant. The authors simulated the measurement error vectors of 17 height differences (black lines in Fig. 2) and they imposed the following elementary deformation scenario: the height of point 6 was increased by 2.5 mm and the height of point 7 was decreased by 2.5 mm, in epoch 2. To simulate the error vectors, pseudorandom numbers were generated 10 times from www.random.org/gaussian-distributions: 34 numbers from the Gaussian distribution in two columns with mean 0.0, standard deviation 1.0 and using 10 significant digits, based on persistent identifiers: 1, 2, …, 10. However, since the observation scheme designed by the authors is quite weak in the sense of network quality control, another—stronger—scheme is investigated in this paper. More precisely, the measures of local internal reliability of the authors’ scheme read: \(r_{jj} \in \left\langle 0.29 - 0.43 \right\rangle\), where \({\mathbf{R}} = {\mathbf{I}} - {\mathbf{A}}_{i} ({\mathbf{A}}_{i}^{{\text{T}}} {\mathbf{Q}}_{{y_{i} }}^{ - 1} {\mathbf{A}}_{i} )^{ - } {\mathbf{A}}_{i}^{{\text{T}}} {\mathbf{Q}}_{{y_{i} }}^{ - 1}\), \(j = 1,...,17\), \(i = 1,2\) and, hence, this scheme does not realize the criterion recommended for geodetic networks: \(r_{jj} > 0.5\) (Prószyński 1994; Prószyński and Kwaśniak 2018). Thus, it is worth extending the original scheme to obtain more realistic results. To do it, 12 additional height differences (grey lines in Fig. 2) are added and \(r_{jj} \in \left\langle 0.55 - 0.73 \right\rangle\) are obtained, where \(i = 1,...,29\). The error vectors of 29 height differences, \({\mathbf{e}}_{1}\), \({\mathbf{e}}_{2}\), are simulated 10 times with the use of the given above web page under the same assumptions as the cited authors (the obtained error vectors are also accessible at: https://www.pracownicy.uwm.edu.pl/krzysztof.nowel/errors.txt) and the same deformation scenario is imposed. In consequence, the following vectors of “measured” height differences are obtained 10 times: \({\mathbf{y}}_{1} (29 \times 1) = {\mathbf{e}}_{1}\) in epoch 1 and \({\mathbf{y}}_{2} (29 \times 1) = {\mathbf{e}}_{2} + {\mathbf{A}}_{2} \nabla\) in epoch 2, where \(\nabla\) = [0 0 0 0 0 2.5 − 2.5 0 0 0 0 0]T mm. These data constitute the first variant. The testing parameters are set to: \(\alpha = 0.10\) and \(\gamma = 0.50\), as it is recommended by Velsink (2015). It is not a priori assumed that a block (a part of the network) shift can be expected, hence only elementary models are engaged in the combinatorial methods. The results of congruence model specification are presented in Table 1. The second variant is based on the same errors as the first variant but differs in deformation scenario. Namely, the following extended deformation scenario is imposed: the heights of points 5, 6, 7, 8 are decreased by 2.5 mm (block shift/subsidence), in epoch 2; \(\nabla\) = [0 0 0 0 − 2.5 − 2.5 − 2.5 − 2.5 0 0 0 0]T mm. It is a priori assumed that a block can be moved, hence both elementary and extended congruence models are engaged in the combinatorial methods. The results are presented in Table 2.
One can see that the CIDIA and V method specify the correct model most frequently (80%, 90%), for both variants. Whereas AICc turns out worse for the elementary deformation scenario (60%) and GCT turns out worse for the extended deformation scenario (50%). Furthermore, GCT and CIDIA specify the same congruence models for the elementary deformation scenario—except run 7. However, only 2/12 points are displaced in this scenario and such similarity can be expected. A more detailed study is left to the reader.
To gain further insight into the working of the CIDIA procedure—especially in the context of GCT—a third variant is discussed. This variant is very easy to repeat and the author suggests such repetition for the reader to check the understanding of CIDIA. This time, only the sub-network consisting of the points: 1, 2, 3, 5, 6, 7 (Fig. 2, dotted area) is analyzed. The elementary deformation scenario is imposed where points denoted as 2, 5, 6 are individually displaced: \(\nabla\) = [0 2.5 0 3.0 − 2.5 0]T mm and only the errors from the run 4 are taken into consideration: \({\mathbf{e}}_{1}\) = [− 0.05 0.24 1.59 0.37 0.27 − 0.38 0.34 2.24 0.28 0.97 0.53]T mm, \({\mathbf{e}}_{2}\) = [0.32 0.14 1.44 − 2.91 0.14 − 1.67 1.46 1.64 1.35 0.54 − 1.24]T mm. The detailed results of congruence model specification (in fact, the points specified as displaced) are presented in Table 3. It is not a priori assumed that a block shift can be expected, hence only elementary models are engaged in the combinatorial methods.
Generally, the results obtained in Experiment 1 show that the suggested CIDIA procedure works well. However, the experiments do not say much about the efficacy of the considered methods. These experiments concern only some arbitrary deformation scenarios and measurement errors, from many possible. Hence, one should not draw any decisive conclusions here. To measure the efficacy of the suggested procedure against the efficacy of the existing ones, more experiments must be carried out, using comprehensive computer simulations. Such experiments are quite complicated to carry out and very time-consuming for combinatorial methods, and they have not yet been performed/published for any combinatorial method.
4.2 Experiment 2: comprehensive scenarios 1D
For Experiment 2, the same leveling control network is used (Fig. 2). However, this time Monte Carlo (MC) simulations (e.g., Koch 2018) are engaged to generate both deformations and measurement errors. The errors are randomized from the Gaussian distribution under the same stochastic model as before, whereas point displacements are randomized from the uniform interval of \(\nabla_{i} \in \left\langle {2 - 25} \right\rangle\) mm, in two variants: displacements with different signs (variant 1) and with the same sign (variant 2; more often found in geodetic and surveying practice). For each variant, the nine sub-variants of displaced points are randomized: one, two, …, eight and \(\left\langle {{\text{zero}} - {\text{eight}}} \right\rangle\) points are randomly selected to be displaced. The vectors of \(2 \times 29\) height difference errors are simulated independently 1000 times for each sub-variant, using MATLAB software. It is not a priori assumed that a block shift can be expected, hence only elementary congruence models are engaged in the combinatorial methods. Table 4 and Fig. 3 present the efficacy of correct model specification with the use of individual methods. As before, the testing parameters are set to: \(\alpha = 0.10\) and \(\gamma = 0.50\).
Efficacy of correct model specification (%)
One can note that the V method is much less efficacious than the well-established GCT procedure. Indeed, the efficacy is very high and similar for both methods when none or one point is displaced. Nevertheless, V is substantially less efficacious than GCT in the other cases. The general efficacy of AICc is also substantially higher than the efficacy of V. Indeed, the efficacy of AICc is much lower and highly unsatisfactory when none, one or two point are displaced. Nevertheless, AICc is substantially more efficacious in the other cases. The relationship between AICc and GCT turns out quite ambiguous. AICc begins to be more efficacious than GCT only when the number of displaced points is greater than five or six. Earlier, the efficacy of AICc is lower than of GCT. CIDIA turns out to be the most efficacious method and its advantage over GCT is especially visible in variant 2. Indeed, the efficacy of CIDIA is similar to the efficacy of AICc when seven and eight points are displaced. Nevertheless, CIDIA is substantially more efficacious in the other cases. In the general case, i.e., when the number of displaced points is randomly selected from interval \(\left\langle {{0} - {8}} \right\rangle\), the order is as follows: CIDIA (89.8%, 88.6%), GCT (84.3%, 75.6%), AICc (68.9%, 64.5%), V (36.7%, 39.7%). In the author’s opinion, the obtained advantage of CIDIA can be explained by the reasons which already have been discussed in Sect. 1 and Sect. 3.1. For instance, one can see that the efficacy of GCT severely decreases as the number of displaced points increases, especially in variant 2 where displacement smearing can be particularly severe. At the same time, the efficacy of CIDIA does not substantially decrease. Hence, one may say that the displacement smearing can burden the results of GCT much more than the results of CIDIA.
The obtained results show that the efficacy of correct model specification substantially differs between the considered methods. However, this does not necessarily mean that, for example, the accuracy of estimated deformations in the specified models would be substantially different. Admittedly, this paper is devoted to the specification of congruence models. Nevertheless, it seems also worth investigating the effect on estimate accuracy which wrong specifications of individual models have. Therefore, the displacements of individual points are additionally estimated—based on earlier specified models—and the root mean squares (RMSs) of the true errors of estimated displacements are calculated, independently for each of 1000 runs. Only the general deformation scenario when the number of displaced points is randomly selected from interval \(\left\langle {{0} - {8}} \right\rangle\) and displacements have different signs is analyzed here. Figure 4 depicts the empirical distributions of the individual RMSs of these errors.
Distributions of true error RMS of estimated deformations (mm)
The results appear quite interesting and AICc especially positively surprises. This means, although the MSRs of correct model specification by AICc have been substantially lower than by CIDIA for the general case (Table 4), the accuracies of estimated deformations are very similar (Fig. 4), i.e., mean RMSs of estimate errors are in the range: \(\left\langle {{0}{\text{.72}} - {0}{\text{.73}}} \right\rangle\) mm. A deeper analysis allows to explain it. It turns out that the majority of wrong models obtained by AICc specify more displaced points than necessary. More precisely, all displaced point(s) and, additionally, some stable point(s) are often jointly modeled as displaced (model overfitting). The results in Table 1 also clearly show it. Note that this kind of model misspecifications is less dangerous for the estimate accuracy than underfitting where not all displaced point(s) are modeled as displaced. Despite overfitting, the congruent part of deformation model includes only stable points and, as a consequence, the datum for estimates is properly defined. Indeed, there are less datum points than there should be, although the most important is that the reference base does not include any displaced point(s). This is why the estimated deformations based on models specified by AICc have such high accuracies. The reverse kind of misspecification interferes with GCT. This means, although the efficacy of GCT has been substantially higher than of AICc, the accuracies of estimated deformations are lower, i.e., the mean RMS of estimate errors is equal to 0.91 mm. It is due to underfitting. A deeper analysis shows that the majority of wrong models obtained by GCT specify less displaced points than necessary and, in consequence, some displaced point(s) also constitute a reference base.
4.3 Experiment 3: arbitrary scenario 2D
Although the conclusions drawn from Experiment 2 may be convincing, however, they refer only to some leveling control network and, hence, another network should also be taken into consideration. This time, the very popular geodynamic ‘Delft’ network (e.g., Welsch and Zhang 1983, pp. 329–336; Caspary 2000, pp. 154–157; Niemeier 2008, pp. 450–456) designed by the FIG working group to check and compare different approaches into deformation analysis is the subject of investigation. Note that this network was also used either by Velsink (2015) or Lehmann and Lösler (2017) to demonstrate their combinatorial methods.
The ‘Delft’ network is a 2D scheme with 14 common control points, where directions and distances are measured in two epochs (106 and 96 observations in epoch 1 and 2, respectively). Four—free of outliers—datasets for deformation analysis were created (e.g., Welsch and Zhang 1983, pp. 329–336). Only the dataset denoted as 1-2B in Welsch and Zhang (1983, p. 331, 333), where all the common points of the eastern part are moved as a block in the X-direction of 20.0 cm and in the Y-direction of 12.0 cm (see also Niemeier 2008, pp. 450–456), is engaged in this experiment (Fig. 5; the green arrows). Note that the dataset was also used by Velsink (2015) and Lehmann and Lösler (2017). To be comparable to the authors’ studies, the same standard deviations, i.e., \(\sigma_{r} = 0.50{\text{ mgon}}\), \(\sigma_{s} = 1.0{\text{ cm}}\) for epoch 1 and \(\sigma_{r} = 0.12{\text{ mgon}}\), \(\sigma_{s} = 10.5{\text{ cm}}\) for epoch 2 (2B), and the same testing parameters: \(\alpha = 0.10\) and \(\gamma = 0.50\) are used. Admittedly, the stochastic model is given by Niemeier (2008, pp. 452–454). Nevertheless, it is quite inadequate to the simulated observation sets and the variance factor—unlike in Velsink (2015) and Lehmann and Lösler (2017)—is estimated in this experiment (\(\hat{\sigma }_{0}^{2} = 0.46\)); the adequate stochastic model for this data is given by Welsch and Zhang (1983, p. 329) or Caspary (2000, p. 157). The estimated standard deviations of estimated displacement components are in the ranges: \(\hat{\sigma }_{{\hat{\nabla }_{{X_{i} }} }} \in \left\langle {2.7 - 7.9} \right\rangle\) cm and \(\hat{\sigma }_{{\hat{\nabla }_{{Y_{i} }} }} \in \left\langle {2.5 - 5.8} \right\rangle\) cm.
‘Delft’ network together with the a priori 95% confidence ellipses of estimated displacements (based on Niemeier 2008, p. 451)
Following Lehmann and Lösler (2017), specification matrices (scenarios) with the maximum possible number of 12 displaced points are considered in combinatorial methods; two points are needed to define the joint network datum. Note that the fault zone is a priori unknown in this experiment. The results of congruence model specification together with estimated displacements are presented in Table 5. For cognitive purposes, either elementary or elementary and extended congruence models are engaged in the combinatorial methods.
To supplement the above results, it is also worth showing how GCT and CIDIA have been specified displaced points in subsequent iteration steps. The results have been going as follows: I (41), II (39), III (11), IV (21), V (3), VI (5), VII (15) in GCT and I (41), II (39, 41), III (11, 39, 41), IV (3, 5, 11, 39), V(3, 5, 11, 39, 41) in CIDIA. As one can see, GCT has given erroneous results: in identification IV where stable point 21 (\(t_{\Delta }^{\max } = 7.56\)) has been identified as displaced instead of displaced point 3 (\(t_{\Delta } = 5.65\)) or displaced point 5 (\(t_{\Delta } = 2.63\)) and, then, in detection VII where \(T = 1.59\) under \(F_{0.90} (13,114,0) = 1.58\).
Generally, from the obtained results one can see that CIDIA works well. Admittedly, CIDIA—unlike V—is not able to specify the displaced points as a block shift in the second case. Nevertheless, it is due to the unfavorable choice of significance level (\(\alpha = 0.10\)). For lower values, e.g. \(\alpha = 0.05\), the results also become proper for CIDIA (not published here). GCT—unlike CIDIA—specifies overfitted model. The reason may be the displacement smearing which leads to the erroneous identification IV and too high significance level which leads to unnecessary positive detection VII (significance level does not decrease in GCT, in consecutive iterations). Note also that the V method pinpoints two equivalent solutions using extended models. Nevertheless, it is not a paradox, but a justifiable result. It is not known—based only on the terrestrial geodetic observations available here—which block is moved and which remained at rest, or if both moved differently in case of relative deformation network.
4.4 Experiment 4: comprehensive scenarios 2D
For Experiment 4, the ‘Delft’ control network is also used, but this time the MC simulations are engaged to generate both deformations and measurement errors (as in Sect. 4.2). The errors are randomized from the Gaussian distribution under the same stochastic model as before, whereas the lengths of point displacements are randomized from the uniform interval: \(d_{i} \in \left\langle {l_{i} - gl_{i} } \right\rangle\), where \(l_{i} = \sqrt {a_{i} \cdot b_{i} }\) are the radiuses of the circles whose areas are equal to the areas of a priori confidence ellipses of estimated displacements (here, \(l_{i} \in \left\langle {8.8 - 18.0} \right\rangle\) cm), \(a_{i} , \, b_{i}\) are the semi-axes of these ellipses, g > 1 and i indicates a given point. The confidence level of ellipses is coupled with the assumed global confidence level, \(1 - \alpha = 0.90\), by the B-method (\(\gamma = 0.50\)) and it reads \(1 - \alpha_{i} = 0.9934\). The general idea of such simulation method comes from the paper Hekimoglu et al. (2010). Two scenarios of deformations are considered: points move as a block and points move individually. The sets of 106 (epoch 1) and 96 (epoch 2) observations are simulated independently 1,000 times for each considered variant, using MATLAB software. Since standard deviations are a priori known and they should be quite adequate for the simulated observations here, it is assumed that the variance factor is known.
4.4.1 The eastern block is homogeneously displaced
The deformations have the same nature as those previously analyzed in Experiment 3, i.e., points 3, 5, 11, 39, 41 move as a block in the X and Y-direction of some positive values, but, this time, the values of the movement lengths, \(d_{i}\), are randomly selected 1,000 times from the interval: \(\left\langle {l_{i} - 3l_{i} } \right\rangle\). The term ‘positive values’ means that the azimuths/directions of movements, \(A_{i}\), are randomized from the uniform interval: \(\left\langle {0 - 100} \right\rangle\) gon. Finally, the simulated displacement vectors of displaced—as a block—points 3, 5, 11, 39, 41 takes the form:
Of course, the displacement vectors with the same values and directions have to be assigned for all five displaced points here, hence the confidence ellipse referred only to one (arbitrary indicated) point should be used to calculate the lengths of these vectors and point 3 was indicated for it (\(l_{3} = 15.2\) cm). The areas where the eastern sub-network can now be moved in the second epoch are depicted using the green dotted lines in Fig. 5. It is a priori assumed that a block shift can be expected here, hence both elementary and extended models are considered in the combinatorial methods. The obtained efficacy (MSR) of correct model specification (i.e., the model which specifies points 3, 5, 11, 39, 41 as a block shift) with the use of GCT, V, AICc and CIDIA turns out equal to: 37.2%, 93.6%, 93.6%, 87.8%, respectively. It is clearly shown that either CIDIA or the other combinatorial methods provide relatively high efficacy, substantially higher than the well-established GCT procedure. This can be explained by the possibility of rigid body movement specification and, hence, the combinatorial methods, especially V and AICc may be recommended for such geodynamical studies.
4.4.2 The entire network is variously deformed
This experiment refers to the elementary scenario of deformations which is most often found in geodetic and surveying practice. Some randomly selected points—from the group of all common, 14 points (Fig. 5)—move individually (as in Sect. 4.2). The lengths, \(d_{i}\), of simulated displacement vectors (56) are randomized independently 1000 times from the interval: \(\left\langle {l_{i} - 5l_{i} } \right\rangle\) but, this time, the individual vectors can have different values and directions for selecting points. Two deformation variants are considered: point displacements with different signs (variant 1): \(A_{i} \in \left\langle {0 - 400} \right\rangle\) gon and with the same sign (variant 2; more often found in geodetic and surveying practice): \(A_{i} \in \left\langle {0 - 100} \right\rangle\) or \(A_{i} \in \left\langle {100 - 200} \right\rangle\) or \(A_{i} \in \left\langle {200 - 300} \right\rangle\) or \(A_{i} \in \left\langle {300 - 400} \right\rangle\) gon. For each variant, the three sub-variants of displaced points are randomized: four, eight and \(\left\langle {{\text{zero}} - {\text{eight}}} \right\rangle\) points are randomly selected to be displaced. It is not a priori assumed that a block shift can be expected here, hence only elementary models are engaged in the combinatorial methods. Table 6 shows the efficacy of correct model specification with the use of individual methods. For cognitive purposes, three sets of testing parameters are considered in the methods involving hypothesis tests.
Generally, the results are fairly consistent with those obtained for the 1D network (Table 4). It should be emphasized, however, that AICc works here better than before. In the general case, the efficacy of AICc is similar or even higher (variant 2) than of GCT. It can be explained by the better working of penalty term. More precisely, the subsequent congruence models differ in two explanatory parameters in a 2D network and, hence, the AICc values can better separate the groups of constant-dimensional models. Despite the above, the efficacy of AICc is still lower than the one of CIDIA, especially in variant 2 which is more often found in practice. However, the differences are not so large as for the 1D network. In the general case, the efficacy of CIDIA (under \(\alpha = 0.10\), \(\gamma = 0.50\)) read: 77.5%, 77.2% and of AICc read: 72.0%, 69.1%.
The debatable testing parameters do not have such importance as it could be expected. To clearly show this, Fig. 6 illustrates the efficacy under the three considered sets of testing parameters, for the general case in variant 1 (full bars) and variant 2 (empty bars). One can see that the efficacies are similar for the considered testing parameters and it is hard to draw any decisive conclusions. One may only conclude that the testing parameters recommended by the cited authors, i.e., \(\alpha = 0.10\) and \(\gamma = 0.50\), work well and they may be also recommended for the CIDIA procedure.
Efficacy of individual methods under various testing parameters for the general case in variant 1 (full bars) and variant 2 (empty bars)
To supplement all the above results, it is also worth investigating the nature of incorrect models, i.e., how often the incorrect models are underfitted and overfitted. Table 7 presents such mean failure rates (MFRs) for variant 1 and under testing parameters: \(\alpha = 0.10\), \(\gamma = 0.50\).
It can be seen that the results somehow correspond to those obtained for the 1D network. Recall that AICc tended to overfitting and GCT, CIDIA to underfitting of incorrect models. In consequence, AICc produced higher—than expected—estimate accuracies (Fig. 4). One can see that the same tendency is visible here. AICc—unlike GCT and CIDIA—specifies more overfitted than underfitted models (Table 7). Note also that the V method strongly tends to specify underfitted models and this diagnosis is in line with Velsink (2018, p. 10).
Finally, since the combinatorial approach has—more or less—“brute force” nature, it is also worth presenting the computational load of the considered methods. Table 8 shows how long it takes to specify a model in one run, for variant 1 and testing parameters: \(\alpha = 0.10\), \(\gamma = 0.50\).
Computation times refer to a PC using Core i7-7700 CPU at 3.60 GHz.
The high computational costs of V, AICc and CIDIA must not be surprising in the light of their combinatorial nature. Nevertheless, thanks to its iterative nature, the CIDIA procedure is substantially less time-consuming than the other combinatorial methods. (The smaller the number of displaced points is, the faster the CIDIA is.) Note that the sub-groups of control network points considered as congruent mainly consist of stable points (or even all these points are stable) in geodetic and surveying practice. Hence, the computational advantage of CIDIA over the other combinatorial methods can be substantial in practice. However, one may conclude that the combinatorial methods are not substantially limited here by computational time. All methods specify congruence models in near real-time and, as such, may be successfully used in most geodetic and surveying deformation analyses. Furthermore, computer technology is still making progress, which will enable the working of these methods increasingly faster in the future.
5 Summary and conclusions
The global congruence test (GCT) procedure is the most commonly used method for the specification of deformation congruence models in geodesy and surveying. Nevertheless, this procedure has a serious weakness in the case of multiple displacements. Namely, the congruent part of model is here point-by-point reduced in consecutive iteration steps; the point which is identified as most suspected to be displaced in a given step is not then modeled as stable in the next steps. The problem is that the consecutive point-by-point identifications may be burdened by displacement smearing which can lead to the identification of a stable point as most suspected to be displaced. If it is the case in a given iteration step, the models which are specified in the next steps are—from now—conditioned by the erroneous information. In recent years, to overcome the above weakness, a revolutionary—in the context of GCT—concept involving combinatorial possibilities was suggested. Since the best approximating model is identified there directly from the set of all possible candidate models, the problem of the consecutive point-by-point specification which may suffer from displacement smearing is avoided. Unfortunately, another serious problem arises, i.e., the problem of the comparison of different-dimensional models. It is not exactly clear what criterion—from many available—should be used to identify the best approximating congruence model from the set of different-dimensional models.
This paper makes a step forward in this new field and discusses a more sophisticated combinatorial procedure—designated as CIDIA—to address the weaknesses associated both with the GCT procedure and the existing combinatorial concept. Thanks to the appropriately used the possibilities of combinatorics and GLR tests performed in the DIA steps, both the problem of the consecutive point-by-point specification which may suffer from displacement smearing and the problem of the comparison of different-dimensional models have been overcome. In the context of GCT, the CIDIA procedure has rather evolutionary—than revolutionary—character and the general concepts of both procedures have similar heuristic substantiation.
To objectively measure the efficacy of CIDIA against GCT and two existing combinatorial methods, various deformation scenarios have been randomized independently many times with the use of comprehensive computer simulations and then processed. Thanks to this, the objective mean success rates of correct model specification could be obtained for the considered methods under individual deformation scenarios. It is worth noting that such studies are quite complicated to carry out, very time-consuming and were not earlier performed/published for the combinatorial concept. Most generally, the obtained results have shown that the best (most efficacious) existing combinatorial method—supported by one of the estimators of Kullback–Leibler divergence, designated as AICc—has not turned out unambiguously more efficacious than GCT. Admittedly, it has given a higher efficacy when the number of displaced points has been relatively high (higher than 50%). Nevertheless, it has not been the case when the number of displaced points has been lower. One may conclude that the problem of the comparison of different-dimensional models has somehow come out here. Whereas, unlike AICc, CIDIA has turned out unambiguously more efficacious than GCT. As expected, the efficacy of CIDIA has been higher than the one of GCT for all the deformation scenarios with multiple displacements. The more the displaced points have been, the higher the advantage of CIDIA over GCT has been. Note also that the advantage of CIDIA has been more visible when displacements have had the same sign. Since the problem of displacement smearing occurs especially under such deformation scenarios, i.e., when the number of displaced points is relatively high and/or the displacements have the same sign, one may conclude that CIDIA has turned out substantially more resistant to this problem than GCT.
Finally, the discussed CIDIA procedure can be employed for other issues in geodesy and surveying. This is due to the fact that the procedure deals with the universal problem of discrepancies between a random sample and its underlying functional model [the functional model has the universal form discussed by Teunissen (2006b, pp. 71–123)]. Thanks to this, CIDIA can be easily implemented, e.g., for the specification of various observation models with outliers. It can be expected that the procedure should be more efficacious for the cases of multiple outliers than the well-established iterative data snooping procedure.
6 Data availability statement
Part of the input data generated during this research are included in this article. The rest of input data can be found at: www.random.org/gaussian-distributions or https://www.pracownicy.uwm.edu.pl/krzysztof.nowel/errors.txt—Sect. 4.1 and in the given references: Welsch and Zhang (1983, pp. 329–336), Caspary (2000, pp. 154–157) or Niemeier (2008, pp. 450–456)—Sect. 4.3. Furthermore, both all the input data and the MATLAB codes are available from the author on reasonable request.
References
Arnold SF (1981) The theory of linear models and multivariate analysis. Wiley, New York
Aydin C (2012) Power of global test in deformation analysis. J Surv Eng 138(2):51–55
Aydin C (2017) Effects of displaced reference points on deformation analysis. J Surv Eng 143(3):04017001
Aydin C, Demirel H (2005) Computation of Baarda’s lower bound of the non-centrality parameter. J Geod 78(7–8):437–441
Baarda W (1968) A testing procedure for use in geodetic networks. Netherlands Geodetic Commission, Publ. on Geodesy, New Series 2(5)
Baselga S (2011) Exhaustive search procedure for multiple outlier detection. Acta Geod Geophys Hung 46(4):401–416
Biagi L, Caldera S (2013) An efficient leave one block out approach to identify outliers. J Appl Geod 7(1):11–19
Burnham KP, Anderson DR (2002) Model selection and multimodel inference: a practical information-theoretic approach, 2nd edn. Springer, Berlin
Caspary WF (2000) Concepts of network and deformation analysis. The University of New South Wales, Kensington
Caspary WF, Borutta H (1987) Robust estimation in deformation models. Surv Rev 29(223):29–45
Caspary WF, Haen W, Borutta H (1990) Deformation analysis by statistical methods. Technometrics 39(1):49–57
Chen YQ (1983) Analysis of deformation surveys—a generalized method. Technical Report No. 94. University of New Brunswick. Fredericton
Chen YG, Chrzanowski A, Secord JM (1990) A strategy for the analysis of the stability of reference points in deformation surveys. CISM J ACSGG 44(2):141–149
Chrzanowski A, Chen YQ (1990) Deformation monitoring, analysis and prediction-status report. FIG XIX International Congress. Helsinki 6(604.1):83–97
de Heus HM, Martens MHF, Verhoef HME (1994) Stability-analysis as part of the strategy for the analysis of the Groningen gas field levellings. In: Proceedings of the perlmutter workshop on dynamic deformation models, Haifa. pp 259–272
Ding X, Coleman R (1996) Multiple outlier detection by evaluating redundancy contributions of observations. J Geod 70:489–498
Duchnowski R, Wiśniewski Z (2019) Robustness of squared Msplit(q) estimation: a theoretical approach. Stud Geophys Geod 63:390–417
Duchnowski R, Wiśniewski Z (2020) Robustness of squared Msplit(q) estimation: empirical analyses. Stud Geophys Geod 64:153–171
Durdag UM, Hekimoglu S, Erdogan B (2018) Reliability of models in kinematic deformation analysis. J Surv Eng 144(3):04018004
Erdogan B, Hekimoglu S (2014) Effect of subnetwork configuration design on deformation analysis. Surv Rev 46(335):142–148
Gaida W, Koch KR (1985) Solving the cumulative distribution function of the noncentral F-distribution for the noncentrality parameter. Sci Bull Stanisl Staszic Univ Min Metall Geod 90(1024):35–43
Heck B (1983) Das analyseverfahren des geodätischen instituts der Universität Karlsruhe. HSBW. Heft 9:153–172
Hekimoglu S, Koch KR (1999) How can reliability of the robust methods be measured? In: Altan MO, Gründig L (eds) Third Turkish–German joint geodetic days, Istanbul. pp 179–196
Hekimoglu S, Koch KR (2000) How can reliability of the test for outliers be measured? Allg Vermess Nachr 107(7):247–253
Hekimoglu S, Demirel H, Aydin C (2002) Reliability of conventional deformation analysis methods for vertical networks. In: FIG XXII international congress. Washington, pp 1–13
Hekimoglu S, Erdogan B, Butterworth S (2010) Increasing the efficacy of the conventional deformation analysis methods: alternative strategy. J Surv Eng 136(2):53–62
Heunecke O, Kuhlmann H, Welsch WM, Eichhorn A, Neuner H (2013) Handbuch Ingenieurgeodäsie: auswertung geodätischer Überwachungsmessungen, 2nd edn. Wichmann, Heidelberg
Kargoll B (2012) On the theory and application of model misspecification tests in geodesy. Deutsche Geodätsche Kommission Reihe C, Nr. 674, München
Koch KR (1985) Ein statistisches auswerteverfahren für deformationsmessungen. Allg Vermess Nachr 92(3):97–108
Koch KR (2010) Parameter estimation and hypothesis testing in linear models. Springer, Berlin
Koch KR (2018) Monte Carlo methods. GEM Int J Geomath 9(1):117–143
Kok JJ (1982) Statistical analysis of deformation problems using Baarda’s testing procedures. Published in forty years of thought Delft university of technology. Delft 2:470–488
Kok JJ (1984) On data snooping and multiple outlier testing. NOAA technical report NOS NGS 30, national geodetic survey, Rockville, Md
Lehmann R, Lösler M (2016) Multiple outlier detection: hypothesis tests versus model selection by information criteria. J Surv Eng 142(4):04016017
Lehmann R, Lösler M (2017) Congruence analysis of geodetic networks—hypothesis tests versus model selection by information criteria. J Appl Geod 11(4):271–283
Neumann I, Kutterer H (2007) Congruence tests and outlier detection in deformation analysis with respect to observation imprecision. J Appl Geod 1(1):1–7
Niemeier W (1981) Statistical tests for detecting movements in repeatedly measured geodetic networks. Tectonophysics 71:335–351
Niemeier W (1985) Deformationsanalyse. Geodaetische Netze in Landes-und Ingenieurvermessung II. In: Pelzer (Hrsg), Wittwer, Stuttgart, 559–624
Niemeier W (2008) Ausgleichungsrechnung, statistische auswertemethoden, 2nd edn. de Gruyter, Berlin
Nowel K (2015) Robust M-estimation in analysis of control network deformations: classical and new method. J Surv Eng 141(4):04015002
Nowel K (2018) A rigorous and user-friendly algorithm for global congruency test in geodetic deformation analysis. In: Baltic geodetic congress (BGC Geomatics), Olsztyn, Poland
Nowel K (2019) Squared Msplit(q) S-transformation of control network deformations. J Geod 93:1025–1044
Nowel K, Kamiński W (2014) Robust estimation of deformation from observation differences for free control networks. J Geod 88(8):749–764
Nowel K, Cellmer S, Fischer A (2020) Validation of GNSS baseline observation models using information criteria. Surv Rev. https://doi.org/10.1080/00396265.2020.1790168
Pelzer H (1971) Zur analyse geodätischer deformationsmessung. Deutsche Geodätische Kommission, Reihe C, Heft 164
Pelzer H (1974) Neuere ergebnisse bei der statistischen analyse von deformationsmessungen. In: FIG XIV international congress. Washington 608.3
Prószyński W (1994) Criteria for internal reliability of linear least squares models. Bull Geod 68:162–167
Prószyński W, Kwaśniak M (2018) Analytic tools for investigating the structure of network reliability measures with regard to observation correlations. J Geod 92(3):321–332
Teunissen PJG (1985) The geometry of geodetic inverse linear mapping and non-linear adjustment. Publications on Geodesy New Series, Netherlands Geodetic Commission 8(1)
Teunissen PJG (1990) An integrity and quality control procedure for use in multi sensor integration. In: Proceedings ION GPS, Colorado Spring. pp 513–522 (republished in ION Red Book Series, vol 7, 2010)
Teunissen PJG (2006a) Network Quality Control. Delft University Press, Delft
Teunissen PJG (2006b) Testing theory: an introduction. Delft University Press, Delft
Teunissen PJG (2018) Distributional theory for the DIA method. J Geod 92(1):59–80
Velsink H (2015) On the deformation analysis of point fields. J Geod 89(11):1071–1087
Velsink H (2018) Testing methods for adjustment models with constraints. J Surv Eng 144(4):04018009
Welsch MW, Heunecke O (2001) Models and terminology for the analysis of geodetic monitoring observations. In: Official report of the ad-hoc committee of FIG working group 6.1. X FIG international symposium on deformation measurements. Orange, USA, pp 390–412
Welsch WM, Zhang Y (1983) Einige methoden zur untersuchung kongruenter und affiner beziehungen in geodätischen überwachungsnetzen zur ermittlung von deformationen. HSBW. Heft 9:299–336
Wiśniewski Z (2009) Estimation of parameters in a split functional model of geodetic observations (Msplit estimation). J Geod 83(2):105–120
Wiśniewski Z (2010) Msplit(q) estimation: estimation of parameters in a multi split functional model of geodetic observations. J Geod 84(6):355–372
Wiśniewski Z, Zienkiewicz MH (2020) Estimators of covariance matrices in Msplit(q) estimation. Surv Rev. https://doi.org/10.1080/00396265.2020.1733817
Wujanz D, Krueger D, Neitzel F (2016) Identification of stable areas in unreferenced laser scans for deformation measurement. Photogram Rec 31(155):261–280
Wyszkowska P, Duchnowski R (2019) Msplit(q) estimation based on L1 norm condition. J Surv Eng 145(3):04019006
Wyszkowska P, Duchnowski R (2020) Iterative process of Msplit(q) estimation. J Surv Eng 146(3):06020002
Zaminpardaz S, Teunissen PJG (2019) DIA-datasnooping and identifiability. J Geod 93(1):85–101
Zienkiewicz MH (2019) Deformation analysis of geodetic networks by applying msplit estimation with conditions binding the competitive parameters. J Surv Eng 145(2):04019001
Acknowledgements
The author thanks the four anonymous reviewers for their time and constructive comments which have improved the quality of the paper. The author also feels indebted to Prof. Rüdiger Lehmann for his time and some explanations.
Author information
Authors and Affiliations
Contributions
KN designed and performed the research and wrote the manuscript.
Corresponding author
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Nowel, K. Specification of deformation congruence models using combinatorial iterative DIA testing procedure. J Geod 94, 118 (2020). https://doi.org/10.1007/s00190-020-01446-9
Received:
Accepted:
Published:
DOI: https://doi.org/10.1007/s00190-020-01446-9