Minimum-variance-based outlier detection method using forward-search model error in geodetic networks

Geodetic observations are crucial for monitoring landslides, crustal movements, and volcanic activity. They are often integrated with data from interdisciplinary studies, including paleo-seismological, geological, and interferometric synthetic aperture radar observations, to analyze earthquake hazards. However, outliers in geodetic observations can significantly impact the accuracy of estimation results if not reliably identified. Therefore, assessing the outlier detection model's reliability is imperative to ensure accurate interpretations. Conventional and robust methods are based on the additive bias model, which may cause type-I and type-II errors. However, outliers can be regarded as additional unknown parameters in the Gauss–Markov model. It is based on modeling the outliers as unknown parameters, considering as many combinations as possible of outliers selected from the observation set. In addition, this method is expected to be more effective than conventional methods as it is based on the principle of minimal variance and eliminates the interdependence of decisions made in iterations. The primary purpose of this study is to seek an efficient outlier detection model in the geodetic networks. The efficiency of the proposed model was measured and compared with the robust and conventional methods by the mean success rate (MSR) indicator of different types and magnitudes of outliers. Thereby, this model enhances the MSR by almost 40 %–45 % compared to the Baarda and Danish (with the variance unknown case) method for multiple outliers. Besides, the proposed model is 20 %–30 % more successful than the others in the low-controllability observations of the leveling network.

1 Introduction

Conventional tests for outliers and robust M estimation are based on the least squares estimation (LSE). If an observation contains an outlier, the LSE method ceases to be the optimal estimation method in terms of a minimum-variance unbiased estimator. Once outliers are detected and isolated, the LSE can be called an efficient estimation. Otherwise, an undetected outlier has a slight deviation from the normality assumption that may cause a smearing effect on all estimated parameters regardless of whether using LSE directly or indirectly which may be named the local influence function (IF) of the LSE (Gao et al., 1992; Hekimoglu et al., 2010; Nowel, 2020). For different bias intervals, the smearing effect of the LSE that behaves systematically as a function of partial redundancy has been proven by Durdag et al. (2022). Normalized residuals, which would be exposed to the smearing effect, are investigated to identify and isolate outliers by conventional tests for outliers and some robust methods such as M estimation (Zienkiewicz and Dąbrowski, 2023; Wang et al., 2021; Batilović et al., 2021). Thereby the falsified test result may induce type-I errors. In addition to the low efficiency of LSE results, the low success of the F test shown by Hekimoglu led researchers to seek a more reliable and effective method such as the univariate method and original observations (Hekimoglu, 1999; Erdogan, 2014; Hekimoglu et al., 2014). Although these methods boost the reliabilities of the conventional methods, the identification of outliers in these models is based on the same procedure as the conventional and robust methods. If the normalized residual exceeds 3 times its standard deviation (SD), also called the 3σ rule, an observation is flagged as an outlier (Lehmann, 2013). However, tests for outliers can be dealt with sufficiently with a single outlier since the LSE has an unbounded IF (Duchnowski, 2011; Maronna et al., 2019; Huber, 1981; Durdag et al., 2022). Studies show that the reliability of these techniques, established with the additive bias model, decreases significantly as the number of outliers increases. In the decision stage, the outliers that mask or swamp other observations can produce a type-I error (false negative) and type-II error (false positive). Multiple outliers can be identified at most as the number of possible outliers ($m_{\max }\le \frac{n-u}{\mathrm{2}}$) by repetitive test procedures (Hekimoglu, 2005). However, the efficiency of conventional tests is rather small when the outlier value is close to the critical value, i.e., when small outliers lie between 3–6σ.

If the rate of successful detection of an outlier using conventional and robust methods is 50 % and one outlier is determined incorrectly, the probability of correctly determining two outliers remains below 50 %. This condition is based on the interdependence of each iteration. Incorrect determination at each step also reduces the possibility of identifying more than one outlier in the next step. Therefore, besides modeling the outliers as unknown, the proposed method is based on two essential factors: the principle of the slightest variance and the assumption of looking at all points with suspicion in each iteration. It has been proven by Hekimoglu et al. (2015) for linear regression that the method with which the outlier is modeled as an additional unknown gives more successful results than the robust methods. The method suggests carrying outlier detection until all possible combinations are investigated. In the $C_k^n$ combination, observation(s) is (are) included as an additional unknown parameter(s) in the proposed model. Then, observations are viewed with suspicion considering combinations of k elements (groups of two, three, and so forth) selected from a set of n elements ($C_k^n$), where n is the number of observations and k denotes the number of outliers. The observation with the smallest variance among $C_{\mathrm{1}}^n$ combinations is determined. Considering the $C_{\mathrm{2}}^n$ combinations, the pair of observations was regarded as a model error and the two observations with the smallest variance were flagged as candidates. All possible combinations will be regarded until the maximum number of burdened observations that would occur up to one-half of the degrees of freedom ($\cong f/\mathrm{2}$) for the geodetic network. The potential observations are clustered separately and compared with the specified critical values for each combination step. The model errors of the potential outlier(s) exceeding the critical value are flagged as suspicious for each combination step. The test values of all potential outliers must exceed the critical value for each combination step, and, if not, the previous candidates are detected as outliers.

The primary purpose of this study is to apply the proposed outlier detection method to geodetic networks and to establish its efficiency. The suggested model was compared with the robust methods by the mean success rate (MSR) indicator of different types and magnitudes of outliers. As in the classic models, the number of outliers is inversely proportional to the success of the presented method. When outliers have various magnitudes (e.g., small, large, gross, and extreme outliers) and specific observations are not available in the network (observations with low controllability), it has been found that the proposed method is quite successful compared to the conventional and robust methods.

2 Gauss–Markov model

Let A_nxu be a design matrix which has full column rank, i.e., rank (A)=u; P a weight matrix of the observations; x_ux1 a vector of the unknown parameter; l_nx1 an observation vector; ${{\mathbf{C}}_{ll}}_{nxn}$ an a priori covariance matrix of observations; ${{\mathbf{Q}}_{ll}}_{nxn}$ a weighted coefficient matrix of observations; and ${\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}}$ an a priori variance factor, where n and u are the number of observation and number of unknowns, respectively. By adding v_nx1 as a residual vector, one can get $\widehat{\mathit{x}}$, an estimated vector of unknown parameters presented in the following Gauss–Markov model (Koch, 1999):

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathit{l}+\mathit{v}=\mathbf{A}\widehat{\mathit{x}},{\mathbf{C}}_{ll}={\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}}{\mathbf{P}}^{-\mathrm{1}}={\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}}{\mathbf{Q}}_{ll},\text{(2)}& {\displaystyle }\widehat{\mathit{x}}=({\mathbf{A}}^{\mathrm{T}}\mathbf{PA}{)}^{+}{\mathbf{A}}^{\mathrm{T}}\mathbf{P}\mathit{l},\text{(3)}& {\displaystyle }{\mathbf{Q}}_{xx}=({\mathbf{A}}^{\mathrm{T}}\mathbf{PA}{)}^{+},\text{(4)}& {\displaystyle }{\mathbf{Q}}_{vv}={\mathbf{P}}^{-\mathrm{1}}-{\mathbf{AQ}}_{xx}{\mathbf{A}}^{\mathrm{T}},\end{array}$$

where Q_xx denotes a cofactor matrix of the unknown parameters and Q_vv implies the cofactor matrix of the residuals.

2.1 Test for outliers

In Geodesy, procedures for the outlier detection were developed by Baarda (1968) and Pope (1976). If the observations come from the normal distribution, they are called good observations, whereas the burdened observations that contain outliers originate from another distribution. Let l_i be a burdened observation that has δl_i as an outlier; the null hypothesis,

$$\begin{array}{}\text{(5)}& H_{\mathrm{0}}:\mathit{\delta }{l}_i=\mathrm{0}\text{against}{H}_{\mathrm{1}}:\mathit{\delta }{l}_i\ne \mathrm{0},\end{array}$$

is tested. If the observations are uncorrelated and the variance ${\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}}$ is known, the normalized residuals can be written as

$$\begin{array}{}\text{(6)}& w_i={\displaystyle \frac{\left|v_i\right|}{{\mathit{\sigma }}_{\mathrm{0}}\sqrt{q_{v_i{v}_i}}}},\end{array}$$

where w_i is the test value and q_vv is the cofactor of the residual for i=1…n. This is known as Baarda's method (i.e., a data-snooping test). A posteriori variance ($m_{\mathrm{0}}^{\mathrm{2}}$) is calculated in Pope's method given by

$$\begin{array}{}\text{(7)}& {\mathit{\tau }}_i={\displaystyle \frac{\left|v_i\right|}{m_{\mathrm{0}}\sqrt{q_{v_i{v}_i}}}},\end{array}$$

where τ_i is the test value. The observation with the biggest normalized or studentized residual is tested in one loop of the iterations. The test for outliers is used iteratively if the observations contain more than one outlier. The flagged observation is removed when H₀ is rejected. The remaining observations are adjusted once more. Until no more outliers are found, this process is repeated. However the multiple outliers cause swamping or masking effects that make it impossible to distinguish the burdened observations from the good ones. In the following sections, the robust and the proposed methods will be demonstrated to prevent the smearing effect of the LSE.

2.2 Robust methods

M estimation (Huber, 1964) is a generalized form of maximum likelihood estimation. In this paper M estimation from the Huber and Danish methods, commonly chosen to handle outliers in robust statistics, was used to compare the results of the proposed method.

2.2.1 M estimation

The re-weighted LSE is applied iteratively to the non-linear normal equation of the M estimation as follows:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\widehat{\mathit{x}}}^r=({\mathbf{A}}^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathbf{W}}}^r\mathbf{A}{)}^{+}{\mathbf{A}}^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathbf{W}}}^r\mathit{l},\text{(9)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathbf{W}}}^r={\mathbf{PW}}_{\left({\stackrel{\mathrm{‾}}{\mathit{v}}}^{r-\mathrm{1}}\right)},\text{(10)}& {\displaystyle }{\mathbf{W}}_{\left({\stackrel{\mathrm{‾}}{\mathit{v}}}^{\mathrm{0}}\right)}=\mathbf{E},\text{(11)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{v}}}^r=\mathbf{A}{\widehat{\mathit{x}}}^r-\mathit{l},\text{(12)}& {\displaystyle }\mathbf{W}\left(\stackrel{\mathrm{‾}}{\mathit{v}}\right)=\mathrm{diag}\left(\right(W\left({\stackrel{\mathrm{‾}}{v}}_{\mathrm{1}}\right),W\left({\stackrel{\mathrm{‾}}{v}}_{\mathrm{2}}\right),\mathrm{\dots },W\left({\stackrel{\mathrm{‾}}{v}}_n\right),\end{array}$$

where ${\widehat{\mathit{x}}}^r$ equals the $\widehat{\mathit{x}}$ from Eq. (2) for the first iteration and E stands for a unit matrix. r implies a number of iterations and is chosen as 5 in this paper. The weight function of Huber's M estimation is given as follows:

$$\begin{array}{}\text{(13)}& W\left({\stackrel{\mathrm{‾}}{v}}_i^r\right)=\left\{\begin{array}{ll}\mathrm{1}& \left|{\stackrel{\mathrm{‾}}{v}}_i^r\right|\le c\\ & \\ \frac{c}{\left|{\stackrel{\mathrm{‾}}{v}}_i^r\right|}& \left|{\stackrel{\mathrm{‾}}{v}}_i^r\right|>c\end{array}\right.;i=\mathrm{1}\mathrm{\dots }n,\end{array}$$

and the weight function of the Danish method is given by

$$\begin{array}{}\text{(14)}& W\left({\stackrel{\mathrm{‾}}{v}}_i^r\right)=\left\{\begin{array}{ll}\mathrm{1}& \left|{\stackrel{\mathrm{‾}}{v}}_i^r\right|<c\\ \exp \left(-\frac{\left|{\stackrel{\mathrm{‾}}{v}}_i^r\right|}{c}\right)& \left|{\stackrel{\mathrm{‾}}{v}}_i^r\right|\ge c\end{array}\right.;i=\mathrm{1}\mathrm{\dots }n,\end{array}$$

where ${\stackrel{\mathrm{‾}}{v}}_i$ is the residual and c is taken as 1.5σ₀. After the diagonal elements of the $\stackrel{\mathrm{‾}}{\mathbf{W}}$ weight matrix are determined, ${\stackrel{\mathrm{‾}}{\mathit{v}}}^r$ and ${\widehat{\mathit{x}}}^r$ are recalculated for each iteration. The residual that is computed in the final iteration is detected as an outlier if it exceeds 3σ₀.

3 Forward search of model error

The Gauss–Markov model (Eq. 1) is now expanded by the u×1 vector ϵ of additional unknown parameters, also with the n×u design matrix M:

$$\begin{array}{}\text{(15)}& \mathit{l}+\mathit{v}=\left[\mathbf{A}\mathbf{M}\right]\left[\begin{array}{l}\widehat{\mathit{x}}\\ \widehat{\mathit{ϵ}}\end{array}\right],{\mathbf{C}}_{ll}={\mathit{\sigma }}^{\mathrm{2}}{\mathbf{P}}^{-\mathrm{1}}={\mathit{\sigma }}^{\mathrm{2}}{\mathbf{Q}}_{ll},\end{array}$$

where the variance σ² stands for the unit weight of the augmented model and the vector ϵ contains the outliers which are subtracted from the observations. If only the outlier Δ_j is present in the observation l_j, then one should define ϵ=Δ_j and M=e_j, where ${\mathit{e}}_j=\left[\mathrm{0},\mathrm{\dots },\mathrm{0},\mathrm{1},\mathrm{0},\mathrm{\dots },\mathrm{0}\right]$ for j=1…n. The jth component of e_j gets the value 1. For the jth observation with $\mathbf{A}=[{\mathbf{A}}_{\mathrm{1}},\mathrm{\dots },{\mathbf{A}}_j,\mathrm{\dots }{]}^{\mathrm{T}}$, the observation equation is given as

$$\begin{array}{}\text{(16)}& {\mathit{l}}_j+{\mathit{v}}_j={\mathbf{A}}_j^{\mathrm{T}}\widehat{\mathit{x}}+{\widehat{\mathrm{\Delta }}}_j,\end{array}$$

where ${\mathbf{A}}_j^{\mathrm{T}}$ is the jth row vector of A and for the remainder of the observations ${\mathit{l}}_k+{\mathit{v}}_k={\mathbf{A}}_k^{\mathrm{T}}\widehat{\mathit{x}}(k=\mathrm{1},\mathrm{2},\mathrm{\dots },n)$, k≠j. If the outliers exist in the observations, ϵ and M are rewritten as follows:

$$\begin{array}{}\text{(17)}& \mathit{ϵ}={\left|{\widehat{\mathrm{\Delta }}}_j,{\widehat{\mathrm{\Delta }}}_{j+\mathrm{1}}\mathrm{\dots },{\widehat{\mathrm{\Delta }}}_t\right|}^{\mathrm{T}}\text{and}\mathbf{M}={\left[{\mathit{e}}_j,{\mathit{e}}_{j+\mathrm{1}},\mathrm{\dots },{\mathit{e}}_t\right]}^{\mathrm{T}}.\end{array}$$

The estimate of unknown parameters of the augmented model can, therefore, be expressed as follows (Koch, 1999):

$$\begin{array}{}\text{(18)}& \left[\begin{array}{l}\widehat{\mathit{x}}\\ \widehat{\mathit{ϵ}}\end{array}\right]={\left[\begin{array}{ll}{\mathbf{A}}^{\mathrm{T}}\mathbf{PA}& {\mathbf{A}}^{\mathrm{T}}\mathbf{PM}\\ {\mathbf{M}}^{\mathrm{T}}\mathbf{PA}& {\mathbf{M}}^{\mathrm{T}}\mathbf{PM}\end{array}\right]}^{-\mathrm{1}}\left[\begin{array}{l}{\mathbf{A}}^{\mathrm{T}}\mathbf{P}\mathit{l}\\ {\mathbf{M}}^{\mathrm{T}}\mathbf{P}\mathit{l}\end{array}\right],\end{array}$$

where

$$\begin{array}{}\text{(19)}& {\displaystyle }\begin{array}{rl}& {\left[\begin{array}{ll}{\mathbf{A}}^{\mathrm{T}}\mathbf{PA}& {\mathbf{A}}^{\mathrm{T}}\mathbf{PM}\\ {\mathbf{M}}^{\mathrm{T}}\mathbf{PA}& {\mathbf{M}}^{\mathrm{T}}\mathbf{PM}\end{array}\right]}^{-\mathrm{1}}=\\ & \left[\begin{array}{ll}\left({\mathbf{A}}^{\mathrm{T}}\mathbf{PA}{)}^{-\mathrm{1}}\right(\mathbf{E}+{\mathbf{A}}^{\mathrm{T}}{\mathbf{PMSM}}^{\mathrm{T}}\mathbf{PA}({\mathbf{A}}^{\mathrm{T}}\mathbf{PA}{)}^{-\mathrm{1}}& -({\mathbf{A}}^{\mathrm{T}}\mathbf{PA}{)}^{-\mathrm{1}}{\mathbf{A}}^{\mathrm{T}}\mathbf{PMS}\\ -{\mathbf{SM}}^{\mathrm{T}}\mathbf{PA}({\mathbf{A}}^{\mathrm{T}}\mathbf{PA}{)}^{-\mathrm{1}}& \mathbf{S}\end{array}\right],\end{array}\text{(20)}& {\displaystyle }\mathbf{S}=[{\mathbf{M}}^{\mathrm{T}}\left(\mathbf{P}-\mathbf{PA}({\mathbf{A}}^{\mathrm{T}}\mathbf{PA}{)}^{-\mathrm{1}}{\mathbf{A}}^{\mathrm{T}}\mathbf{P}\right)\mathbf{M}{]}^{-\mathrm{1}}=({\mathbf{M}}^{\mathrm{T}}{\mathbf{PQ}}_{vv}\mathbf{PM}{)}^{-\mathrm{1}},\text{(21)}& {\displaystyle }\widehat{\mathit{ϵ}}={\mathbf{SM}}^{\mathrm{T}}\mathbf{P}(\mathbf{E}-\mathbf{A}({\mathbf{A}}^{\mathrm{T}}\mathbf{PA}{)}^{-\mathrm{1}}{\mathbf{A}}^{\mathrm{T}}\mathbf{P})\mathit{l},\end{array}$$

The residuals are expressed for the Gauss–Markov model in Eq. (1) by

$$\begin{array}{}\text{(22)}& \mathit{v}=\mathbf{A}\widehat{\mathit{x}}-\mathit{l}=\mathbf{A}({\mathbf{A}}^{\mathrm{T}}\mathbf{PA}{)}^{+}{\mathbf{A}}^{\mathrm{T}}\mathbf{P}\mathit{l}-\mathit{l}=\left(\mathbf{E}-\mathbf{A}({\mathbf{A}}^{\mathrm{T}}\mathbf{PA}{)}^{+}{\mathbf{A}}^{\mathrm{T}}\mathbf{P}\right)(-\mathit{l}),\end{array}$$

whose right-hand side can be replaced in Eq. (21) as follows:

$$\begin{array}{}\text{(23)}& \widehat{\mathit{ϵ}}={\mathbf{SM}}^{\mathrm{T}}\mathbf{P}\left(\mathbf{E}-\mathbf{A}({\mathbf{A}}^{\mathrm{T}}\mathbf{PA}{)}^{-\mathrm{1}}{\mathbf{A}}^{\mathrm{T}}\mathbf{P}\right)\mathit{l}=-{\mathbf{SM}}^{\mathrm{T}}\mathbf{P}\mathit{v},\end{array}$$

and considering Eq. (20) the following equation yields

$$\begin{array}{}\text{(24)}& \widehat{\mathit{ϵ}}=-({\mathbf{M}}^{\mathrm{T}}{\mathbf{PQ}}_{vv}\mathbf{PM}{)}^{-\mathrm{1}}{\mathbf{M}}^{\mathrm{T}}\mathbf{P}\mathit{v}.\end{array}$$

3.1 Testing procedure

The alternative hypothesis, in the case of the presence of outliers, takes the form against the null hypothesis as follows:

$$\begin{array}{}\text{(25a)}& {\displaystyle }{H}_{\mathrm{0}}:E\left\{l\right\}=\mathbf{A}\widehat{\mathit{x}},\text{(25b)}& {\displaystyle }{H}_A:E\left\{l\right\}=\left[\mathbf{A}\mathbf{M}\right]\left[\begin{array}{l}\widehat{\mathit{x}}\\ \widehat{\mathit{ϵ}}\end{array}\right].\end{array}$$

One should consider all possible combinations of potentially burdened observations for the correct specification of the alternative hypothesis (Teunissen, 2006). All potential alternative hypotheses $C_b^n$, where n is the number of observations and b is the number of potential outliers, are considered in the detection step. Firstly, the observations are assumed to be unknown one by one in the model. The additional unknowns of the model $\widehat{\mathit{ϵ}}$ are calculated by rewriting the relevant rows for each observation in the coefficient matrix iteratively. The design matrix can be rewritten as follows by including a dimension in the model as an unknown:

$$\begin{array}{}\text{(26)}& \begin{array}{rl}{\mathbf{A}}_{C_{\mathrm{1}}^n}^{\mathrm{1},\mathrm{1}}& =\left[\mathbf{A}\mathbf{M}\right]=\left[\begin{array}{lllllllll}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& |& \mathrm{1}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& |& \mathrm{0}\\ \mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& |& \mathrm{0}\\ \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& |& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& |& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{1}& |& \mathrm{0}\end{array}\right],\\ & {\mathbf{A}}_{C_{\mathrm{1}}^n}^{\mathrm{1},\mathrm{2}}=\left[\begin{array}{lllllllll}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& |& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& |& \mathrm{1}\\ \mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& |& \mathrm{0}\\ \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& |& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& |& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{1}& |& \mathrm{0}\end{array}\right],\end{array}\end{array}$$

where A^b,i denotes the matrix of coefficients for $b=\mathrm{1},\mathrm{\dots },f/\mathrm{2}$ and $i=\mathrm{1},\mathrm{\dots },n$.

Figure 1 Flowchart of the forward search of model error.

3.1.1 Calculation steps for model error

The rows of the additional column vector are rewritten iteratively for each observation, and the corresponding one is modeled as an unknown using calculation steps given below.

1. After calculating the cofactor matrix, the unknowns matrix is obtained:

   $$\begin{array}{}\text{(27)}& {\displaystyle }{\mathbf{Q}}_{xx}^b=({{\mathbf{A}}^b}^{\mathrm{T}}{\mathbf{PA}}^b{)}^{+},\text{(28)}& {\displaystyle }{\widehat{\mathit{x}}}^b={\left({{\mathbf{A}}^b}^{\mathrm{T}}{\mathbf{PA}}^b\right)}^{+}{{\mathbf{A}}^b}^{\mathrm{T}}\mathbf{P}\mathit{l}.\end{array}$$

2. To determine the observation that gives the smallest-variance value, the step of calculating the residuals is given by

   $$\begin{array}{}\text{(29)}& {\mathit{v}}^b={\mathbf{A}}^b{\widehat{\mathit{x}}}^b-\mathit{l}.\end{array}$$

3. The posteriori variance is calculated as

   $$\begin{array}{}\text{(30)}& ({\mathrm{s}}^b{)}^{\mathrm{2}}={{\mathit{v}}^b}^{\mathrm{T}}\mathbf{P}{\mathit{v}}^b/{\mathit{f}}^b.\end{array}$$

4. Determining the observation with minimum variance,

   $$\begin{array}{}\text{(31)}& j=min(s^b{)}^{\mathrm{2}}.\end{array}$$

5. After the relevant observation is determined, the test value is calculated as given by

   $$\begin{array}{}\text{(32)}& T={\widehat{\mathrm{\Delta }}}_j/\left(s_{\mathrm{0}}\sqrt{{\mathit{q}}_{jj}}\right).\end{array}$$

Thus, the unit-weighted posteriori variances for each additional unknown parameter are calculated by

$$\begin{array}{}\text{(33)}& {\widehat{s}}^{\mathrm{2}}={\displaystyle \frac{{\mathit{v}}^{\mathrm{T}}\mathbf{P}\mathit{v}}{n-u_k}},i=\mathrm{1}\mathrm{\dots }n,\end{array}$$

where $u_k=u+\mathrm{1}$ represents the number of the unknowns calculated for the model given in Eq. (15). The number of elements in the set of the posterior variances calculated for each observation appears as $C_{\mathrm{1}}^n$. After the acceptance or rejection of the H₀ hypothesis is evaluated in the identification phase mentioned below, the decision is made to rewrite the model, where the unknowns are expanded for the observations two by two for the $C_{\mathrm{2}}^n$ combination. It is identified to which observation the smallest-variance value $min\mathit{\left\{}{\widehat{s}}_i^{\mathrm{2}}\mathit{\right\}}$ belongs, and the unknown of the relevant observation is compared with the critical value. When $min\left\{{\widehat{s}}_i^{\mathrm{2}}\right\}={\widehat{s}}_k^{\mathrm{2}}$, the absolute value of T is compared with the t test. If $\left|T\right|\ge t_{f-\mathrm{1},\mathrm{1}-\mathit{\alpha }}$, H₀ is rejected and the kth observation is flagged as an outlier. If the null hypothesis is accepted, the process ends. The model is expanded for another alternative hypothesis which assume two potential blunders in case H₀ is rejected. The coefficient matrix is rewritten for each combination of $C_{\mathrm{2}}^n$ as

$$\begin{array}{}\text{(34)}& \begin{array}{rl}{\mathbf{A}}_{C_{\mathrm{2}}^n}^{\mathrm{2},\mathrm{1}}& =\left[\begin{array}{llllllllll}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& |& \mathrm{1}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& |& \mathrm{0}& \mathrm{1}\\ \mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& |& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& |& \mathrm{0}& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& |& \mathrm{⋮}& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{1}& |& \mathrm{0}& \mathrm{0}\end{array}\right],\\ & {\mathbf{A}}_{C_{\mathrm{2}}^n}^{\mathrm{2},\mathrm{2}}=\left[\begin{array}{llllllllll}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& |& \mathrm{1}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& |& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& |& \mathrm{0}& \mathrm{1}\\ \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& |& \mathrm{0}& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& |& \mathrm{⋮}& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{1}& |& \mathrm{0}& \mathrm{0}\end{array}\right].\end{array}\end{array}$$

An important point to be emphasized here is that all combinations are taken into account independently of the previous result (i.e., regardless of the biased observation flagged in the previous step). For example, all potential $C_{\mathrm{2}}^n$ combinations are considered, neglecting the previous result where the kth observation was flagged. The test values of model errors (Δ_i, Δ_j) for i=1…n and j=1…n and i≠j, which have the smallest variance, are compared with the $t_{f-\mathrm{1},\mathrm{1}-\mathit{\alpha }}$ threshold value where α=0.05, whether the model errors of the observations that give the smallest-variance value are higher than the critical value or not. If both are greater than the critical value, the relevant observations are flagged as outliers. Possible $C_{\mathrm{3}}^n$ combinations are sought, and the coefficient matrix is rewritten as follows:

$$\begin{array}{}\text{(35)}& \begin{array}{rl}{\mathbf{A}}_{C_{\mathrm{3}}^n}^{\mathrm{3},\mathrm{1}}& =\left[\begin{array}{lllllllllll}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& |& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& |& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& |& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& |& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& |& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{1}& |& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right],\\ & {\mathbf{A}}_{C_{\mathrm{3}}^n}^{\mathrm{3},\mathrm{2}}=\left[\begin{array}{lllllllllll}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& |& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& |& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& |& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& |& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& |& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{1}& |& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right].\end{array}\end{array}$$

If all three values of unknowns exceed the critical value, they are flagged as outliers. This process is repeated for four or more combinations until all the combinations of potentially burdened observations have been considered. The $\widehat{\mathit{ϵ}}$ vector of the observations corresponding to the minimum-variance value calculated for each combination step is compared with the critical value. If at least one of the relevant unknowns of the observations does not exceed the critical value, the H₀ hypothesis is accepted and the observations flagged in the previous step (i.e., the latest rejected H₀) are approved as outliers. The flowchart of the FSME (forward search of model error) model is presented in Fig. 1.

Table 1 MSR of models (small outliers).

Two cases which the variance is known and unknown were considered for robust methods as follows: ^* The a priori variance is known. Here, the c was taken to be 1.5, where c=1.5σ₀ and σ₀=1. When the residual from the robust techniques exceeded the 3σ₀ threshold value, it was regarded as an outlier. In the case where the a priori variance is known, it can be seen in Table 1 that the MSRs of the robust methods are higher by 0.001 than the Baarda test with α considered. Pope's test had a lower MSR than Baarda's test. However, the MSRs of the FSME are higher than the robust methods in both cases where the a priori variance is known and unknown. ${}^{**}$ The a priori variance is unknown. The standard deviation from the first iteration (LSE) was obtained for robust methods. So the c was taken from 1.5m₀. α was chosen as 0.05 for the Pope's test, which had a lower MSR than Baarda's. Except for the Danish* method, all other models identified an excellent observation as an outlier with a risk ranging from 0.01 % to 5 % if there was no outlier in the observations. The a posteriori variance negatively impacted the robust method's results, and the outlier's spoilt variance significantly contributed to the false detection. The a priori variance significantly affects how reliable the procedures are.

Table 2 Local MSRs (small outliers).

4 Leveling network

In statistics, there are different indicators to measure the reliability of tests and estimators. Hekimoglu and Koch (2000) showed that a finite-sample breakdown point determined the global reliability of an estimator and a test procedure. Using the power function of the global test, a capacity in deformation networks is explored as suggested by Niemeier (1985). Also, it has been shown that the MSR results of the two testing procedures (χ² and F test) are identical to their respective test powers known beforehand (Aydin, 2012). MSR depends on the number of outliers, the magnitude of an outlier, the number of unknowns, the number of observations, and the type of outliers. Since it considers these different cases, MSR is more reliable, whereas the power of the test is the same for all disparate conditions. Also, Erdogan et al. (2019) has proven that the MSR is the empirical estimation of the power of the test in outlier detection. In this study, therefore, MSR is used to specify the ability of the conventional, robust, and proposed models. By this purpose three different leveling networks have been simulated. The random errors ε_i for i=1…n were generated using a normal distribution N(0, σ²) with a mean of zero and a variance of σ². Also, the good and contaminated observations were acquired by simulation technique as described in detail by Hekimoglu and Erenoglu (2007). Since the outliers are produced through simulation, it is easy to determine whether an observation is burdened before analyzing. The method is deemed successful if the observation recognized as an outlier matches the really burdened observation. The process is considered unsuccessful if it fails. When the simulated observation is chosen randomly, the successful rate indicates that the global MSR and the local MSR can be computed for each particular observation in the leveling network for 10 000 samples. The same samples were subjected to conventional, robust, and proposed methods to compare their MSRs with different scenarios. This study simulated outliers randomly chosen from small- and large-magnitudes outliers (variously described as gross and influential outliers) for three leveling networks. An influential outlier is a situation that, either independently or when combined with other biased observations, adversely affects the outcomes of an analysis. Even a single influential outlier may ruin the estimation parameters. A leveling network used for the simulation has 7 points and 15 observations as seen in Fig. 2. The precision is considered to be ${\mathit{\sigma }}_i^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{0}}/\sqrt{S}$, where S is the length of the leveling line in kilometers and σ₀=1 mm ${\sqrt{\mathrm{1}\mathrm{km}}}^{-\mathrm{1}}$. MSRs for 10 000 samples were calculated for each method when there were different magnitudes and different numbers of outliers in the network. The small and large outliers were generated in the intervals of [3–6σ] and [6–12σ], respectively.

Figure 2 Leveling network.

Figure 3 Leveling networks 2a and 2b with low redundancy.

As Table 1 shows, even if the number of outliers changes, the MSRs of the proposed method increase significantly compared to the conventional and robust methods. In cases where there is no outlier (e.g., m=0), the results in which the H₀ hypothesis is rejected are also seen in Table 1. The proposed method generated type-II errors at the rate of 5 %, where the significance level was at 0.05.

Additionally, the a posteriori variance is easily influenced by outliers in the data set, which harms the abilities of methods that use the a posteriori variance. The a posteriori variance from the LSE is typically utilized as a threshold value instead of the a priori variance if the a priori variance is unknown. Therefore, the MSRs of the robust techniques of the former case are higher than the latter. As a result of these findings, only the case where the variance is known, which is less affected by an outlier, is taken into account in the results shown in the tables (Tables 2–8) to compare with the FSME hereafter.

5 Results

Extensive experiments have been done to compare the proposed method with robust methods, such as the Danish and Huber methods, besides the conventional outlier detection procedures (i.e., Baarda and Pope). The redundancies are an important indicator to recognize the observations most vulnerable to bias (Durdag, 2022). The redundancy matrix is calculated from $\mathbf{R}=\mathbf{H}-\mathbf{E}$, where $\mathbf{H}=\mathbf{A}({\mathbf{A}}^{\mathrm{T}}\mathbf{PA}{)}^{-\mathrm{1}}{\mathbf{A}}^{\mathrm{T}}\mathbf{P}$ is a hat matrix. The local MSRs have been calculated for the specific observations with the highest and lowest redundancy in the leveling network. Among the observations, those with the two largest redundancies are h₁₃ and h₉, and the three lowest are h₁, h₇, and h₈. As can be seen from the Table 4 below, MSRs increase as the redundancy does.

Table 3 MSR of models (large outliers).

Table 4 The effect of large and low partial redundancies on MSRs for a pair of observations.

Table 5 MSRs for gross and influential outliers.

It is apparent from Table 2 that the highest MSRs for the biased observation h₁ amongst the conventional and robust methods are from the Danish method (54 %). In addition, the MSR of the proposed method is higher than the Danish one by 30 %. The highest MSR has been obtained by the FSME as 94 % for the observation with the highest redundancy h₁₃. As the redundancy gets smaller, the difference in MSR between the proposed method and other methods increases.

As shown in Table 3, the highest MSRs are obtained by the FSME in contrast with other techniques for different numbers of outliers. When Tables 1 and 3 are compared, the MSRs increase with an enlargement in the magnitude of outlier.

The smearing effect of the LSE, almost equivalent to its SC (sensitivity curve), behaves systematically as a function of the partial redundancy (Durdag, 2022). For this reason, the MSRs have been calculated for the pair of observations with the lowest and largest partial redundancy with small outliers in Table 4. The neighboring observations, especially the point that has three leveling lines, are some of the most vulnerable to bias (e.g., h₆h₇ and h₈h₇) in the leveling network (Fig. 2). The local MSRs are lower than the global MSRs, in the case m=2 in Table 1, for observations most vulnerable to bias with both the lowest and highest redundancies as shown in Table 4.

The results, as shown in Table 4, indicate that the MSRs of observation h₆h₇ with the highest partial redundancies increases compared to h₁₃h₁ with the lowest ones by almost 30 % for the Baarda and Danish models and 50 % for the FSME model.

It is apparent from Table 5 that Baarda and Danish are the two most successful methods against gross and influential outliers among classical and robust methods. In the case of small outliers with gross or influential outliers, the robustness of the models has been tested. Different types of outliers have been generated to evaluate the MSRs of the models for various scenarios as follows: (I) gross outlier (50σ), (II) influential outliers (1000σ), (III) a small outlier and a gross outlier, (IV) a small outlier and an influential outlier, (V) two small outliers and a gross outlier, and (VI) two small outliers and an influential outlier.

Table 6 MSRs for small outliers with gross or influential outliers.

Comparing Table 5 with Table 6, it is observable how MSRs of these two methods were affected in case of one or two small outliers. If a small outlier occurs, the MSRs drop dramatically by about 40 %. The MSRs drop to 20 % with two small outliers in the network. Furthermore, this loss is around 35 %–40 % for the proposed method. The FSME, however, stands out as the model with an MSR of 60 %–70 % in scenarios involving small outliers.

When the redundancies of the observations decrease in the leveling network, difficulties arise in determining the outliers due to the swamping and masking effects. Two different leveling networks are considered to obtain the MSR of the methods in such cases. In the first of these, the MSR of the models has been compared by excluding an observation of the network. As seen in Table 7, the MSR decreased by 30 % in all models when m=2, compared with the case m=1. Although the number of small outliers changes, the highest MSRs have been obtained by the FSME for the leveling network 2a. The network is further weakened, so only two lines of the corner point P.5 remain in the leveling network 2b.

Table 7 MSR of models (small outliers) for leveling network 2a.

The results, as shown in Table 8, indicate that the FSME is the model with the highest MSR for m=1. When m>1 is compared with the case m=1 in Table 8, MSRs of the conventional and robust methods show a more dramatic decrease than the FSME model. Comparing the estimated results for the networks 2a and 2b reveals an approximately 15 % drop in MSR values when m=1. MSRs decrease as the controllability of the observations in the network decreases.

Table 8 MSR of models (small outliers) for leveling network 2b.

6 Conclusion

Since geodetic observations are utilized in studies requiring high accuracy for determining deformations, detecting and identifying outliers become increasingly critical. Researchers commonly favor conventional and robust methods based on the additive bias model. However, this study contributes to our understanding by advancing the modeling of outliers as an additional unknown parameter within the Gauss–Markov model. The aim of this study was to evaluate the suitability of the FSME method within geodetic networks. To achieve this objective, the FSME method was applied to a leveling network. The design of the FSME method is based on identifying the minimum variance from all possible combinations that assume observations as model errors in the Gauss–Markov model. Although only a leveling network has been simulated, the functional and stochastic models of the FSME methods can be applied to all type of geodetic networks. This model yields more reliable results by preventing the swamping and masking effect. The MSRs of the suggested method were obtained for various outliers in three different leveling networks. The results of this investigation show that the FSME is a more efficient model than the robust and conventional methods. Specifically, the proposed method enhanced the MSR by almost 40 %–45 % compared to the Baarda and Danish (with the variance unknown case) methods for multiple outliers (i.e., 1<m<4). Moreover, in scenarios where specific observations were absent at corner points in leveling network 1, the proposed method exhibited 20 %–30 % greater success than alternative methods. Despite the proposed model demonstrating a higher MSR than other methods, the FSME method may encounter numerous combinations depending on the presence of observations and outliers. To address this challenge, particularly in real-world applications, the MSS (maximum subsample method) proposed by Neitzel (2004) and Ebeling (2014) offers a promising approach to reduce the number of combinations in outlier detection procedures. As demonstrated by Ebeling in deformation monitoring, MSS holds potential as a valuable tool to enhance the applicability of the proposed method, particularly within extensive geodetic networks.