Next Article in Journal
A Complete Automatic Target Recognition System of Low Altitude, Small RCS and Slow Speed (LSS) Targets Based on Multi-Dimensional Feature Fusion
Next Article in Special Issue
Multi-Temporal InSAR Analysis for Monitoring Ground Deformation in Amorgos Island, Greece
Previous Article in Journal
Body Dimension Measurements of Qinchuan Cattle with Transfer Learning from LiDAR Sensing
Previous Article in Special Issue
Detection of Structural Vibration with High-Rate Precise Point Positioning: Case Study Results Based on 100 Hz Multi-GNSS Observables and Shake-Table Simulation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Efficacy of Msplit Estimation in Displacement Analysis

by
Zbigniew Wiśniewski
,
Robert Duchnowski
* and
Andrzej Dumalski
Institute of Geodesy, University of Warmia and Mazury in Olsztyn, 1 Oczapowskiego St., 10-957 Olsztyn, Poland
*
Author to whom correspondence should be addressed.
Sensors 2019, 19(22), 5047; https://doi.org/10.3390/s19225047
Submission received: 18 October 2019 / Revised: 14 November 2019 / Accepted: 18 November 2019 / Published: 19 November 2019

Abstract

:
Sets of geodetic observations often contain groups of observations that differ from each other in the functional model (or at least in the values of its parameters). Sets of observations obtained at various measurement epochs is a practical example in such a context. From the conventional point of view, for example, in the least squares estimation, subsets in question should be separated before the parameter estimation. Another option would be application of Msplit estimation, which is based on a fundamental assumption that each observation is related to several competitive functional models. The optimal assignment of every observation to the respective functional model is automatic during the estimation process. Considering deformation analysis, each observation is assigned to several functional models, each of which is related to one measurement epoch. This paper focuses on the efficacy of the method in detecting point displacements. The research is based on example observation sets and the application of Monte Carlo simulations. The results were compared with the classical deformation analysis, which shows that the Msplit estimation seems to be an interesting alternative for conventional methods. The most promising are results obtained for disordered observation sets where the Msplit estimation reveals its natural advantage over the conventional approach.

1. Introduction and Motivation

Consider the classical functional model of geodetic observations, which is given for l = 1 , , q different measurement epochs, namely
y = A X + v     y l = A l X l + v l
where y l = [ y 1 , l , , y n l , l ] T are the observation vectors whose elements belong to the respective sets Φ l = { y 1 , l , , y n l , l } ; X l = [ X 1 , l , , X r , l ] T are the parameter vectors; v l = [ v 1 , l , , v n l , l ] T are vectors of random errors; and A l R n l , r are known coefficient matrices. Such models are the basis of deformation analysis, namely for determining the shifts Δ X ( k , l ) = X l X k between the epochs l and k (for example, the changes of the point coordinates between such epochs).
The vectors Δ X ( k , l ) can be estimated by applying different methods or strategies (e.g., [1,2,3]). The least squares method (LS-method) is still the most popular approach in such an analysis, note that LS-estimates are often supplemented with respective statistical tests (e.g., [4,5,6]). However, some unconventional methods are also in use, for example, robust M-estimation [7,8] or R-estimation [9,10,11,12,13,14]. In the case of relative networks, one can also apply methods of free adjustment (e.g., [15,16,17,18]). Some methods as well as their properties are well known, other methods are still being researched.
The Msplit estimation surely belongs in the latter group. The Msplit estimation was proposed by Wiśniewski [19,20] and has been applied to some practical problems in which each observation could be assigned to several different functional models. For example, it was used in remote sensing (terrestrial laser scanning or ALS data) for data modeling [21] and in some geodetic problems, for example, in deformation analysis [22,23,24,25], and robust estimation (e.g., [26]). Automatic assignment of each observation to the best fitted model is one of the most important features of Msplit estimation. It is also very useful in deformation analysis, when the observation set might include observations from all measurement epochs (the set is an unrecognized mixture of such observations). Note that there is usually no problem with separating observations from different epochs and hence with separate analyses. However, there are some cases when the application of the Msplit estimation is advisable. For example, when a point is displaced during an observation session, thus, one should consider two pseudo-epochs, and the Msplit estimation allows us to estimate the parameters of the functional models for such pseudo-epochs. Such models can also be applied when an observation set is disturbed by outliers [23,24]. Note that the method under investigation can be applied in all observation sets that are an unrecognized (and/or unordered) mixture of observation aggregations. Such data can result from different sources or instrumentations. In fact, the source of data does not matter here. It can be, for example, geodetic instruments: total stations, GNSS receivers, etc., or in remote sensing such as terrestrial or airborne laser scanners.
The main properties of the Msplit estimation are discussed in the papers cited above; that study focused on the efficacy of the method in estimating parameters of the competitive functional models, hence also in estimating point displacements. The analyses were based on simulations of the crude Monte Carlo method and the application of elementary functional models or models of a leveling network. The results were compared with the results of the LS-method.

2. Theoretical Foundations

Without loss of generality, we can assume two measurement epochs, thus in the model of Equation (1), we have q = 2 . Then, the optimization criterion of the LS-method and its solution can be written in the following way ( l = 1 , 2 )
φ L S ( X l ) = i = 1 n v i , l 2 p i , l = v l T P l v l = min X ^ L S , l = D L S , l y l
where v i , l = y i , l a i , l X l , D L S , l = ( A l T P l A l ) 1 A l T P l , P l are respective weight matrices, ( a i , l ith row of matrix A ( l ) ). The difference Δ X ^ L S ( 1 , 2 ) = X ^ L S , 2 X ^ L S , 1 is a LS-estimate of the shift Δ X ( 1 , 2 ) .
In the case of the Msplit estimation, we assumed that each observation belonged to either of two sets Φ 1 or Φ 2 ; however, there is one observation set Φ = Φ 1 Φ 2 and one observation vector y = [ y 1 , , y n ] T , n = n 1 + n 2 . There are two competitive functional models
y = A X ( 1 ) + v ( 1 ) , y = A X ( 2 ) + v ( 2 )
with two competitive versions of the parameter X , namely X ( 1 ) and X ( 2 ) ( A R n , r , r a n k ( A ) = r ). The vectors v ( 1 ) , v ( 2 ) R n are two competitive versions of the observation errors related to all elements of the vector y.
The theoretical basis of the Msplit estimation is an assumption that every observation y i can be assigned to either of two density function f ( y i ; X ( 1 ) ) or f ( y i ; X ( 2 ) ) . If y i occurs, it brings the f-information I f ( y i ; X ( 1 ) ) = ln f ( y i ; X ( 1 ) ) or the f-information I f ( y i ; X ( 2 ) ) = ln f ( y i ; X ( 2 ) ) , which are competitive to each other. Msplit estimates of the parameters X ( 1 ) and X ( 2 ) , namely X ^ ( 1 ) and X ^ ( 2 ) , minimize the following global information that is brought by all elements of the vector y [19]
I f ( y ; X ( 1 ) , X ( 2 ) ) = i = 1 n I f ( y i ; X ( 1 ) ) I f ( y i ; X ( 2 ) ) = i = 1 n [ ln f ( y i ; X ( 1 ) ) ] [ ln f ( y i ; X ( 2 ) ) ]
In other words, the estimators in question are the solutions of the following optimization problem:
min I f ( y ; X ( 1 ) , X ( 2 ) ) = min I f ( y ; X ^ ( 1 ) , X ^ ( 2 ) )
For such solutions, the occurrence of the particular observation vector is the most probable. If X ( 1 ) = X ( 2 ) = X , then I f ( y ; X ( 1 ) , X ( 2 ) ) = φ M L ( X ) = i = 1 n [ ln f ( y i ; X ) ] , which is the objective function of the maximum likelihood method (ML-method). In such a context, the Msplit estimation is a special development of the ML-method. Huber [27,28] generalized the ML-method to M-estimation by introducing φ M ( X ) = i = 1 n ρ ( y i ; X ) , where ρ ( y i ; X ) is an arbitrary function for which estimators obtain the desired properties (for example, they are robust against outliers). A similar generalization was also proposed for the Msplit estimation [19,20]. The objective function of Equation (4) is replaced by the following function
φ ρ ( X ( 1 ) , X ( 2 ) ) = i = 1 n ρ ( 1 ) ( y i ; X ( 1 ) ) ρ ( 2 ) ( y i ; X ( 2 ) )
Of course, Msplit estimation is also a development of classical M-estimation, if only X ( 1 ) = X ( 2 ) = X , and hence φ ρ ( X ( 1 ) , X ( 2 ) ) = φ M ( X ) .
There are several variants of Msplit estimation that differ from one another in the objective function or assumed parameters [19,22,29]. So far, the most popular is the squared Msplit estimation for which ρ ( 1 ) ( y i ; X ( 1 ) ) = p i v i ( 1 ) 2 and ρ ( 2 ) ( y i ; X ( 2 ) ) = p i v i ( 2 ) 2 . Hence, one can write the following optimization problem of such a method [19,30] as
φ s q ( X ( 1 ) , X ( 2 ) ) = i = 1 n p i 2 v i ( 1 ) 2 v i ( 2 ) 2 = ( v ( 1 ) v ( 1 ) ) T P 2 ( v ( 2 ) v ( 2 ) ) = min
where P = Diag ( p 1 , , p n ) is a diagonal weight matrix of the observations y ( – the Hadamard product). It is obvious that if X ( 1 ) = X ( 2 ) = X , then φ s q ( X ( 1 ) , X ( 2 ) ) = i = 1 n p i v i 2 , which means that the squared Msplit estimation is a development of the LS method. Considering such a relationship and the range of the practical applications of the Msplit estimation, we will only discuss the squared Msplit estimation. To compute Msplit estimates, one can use the sufficient conditions for the minimum of the objective function. Considering the optimization problem (7), one can write the following equations
g ( 1 ) ( X ( 1 ) , X ( 2 ) ) = ( φ ( X ( 1 ) , X ( 2 ) ) / X ( 1 ) ) X ( 1 ) = X ^ ( 1 ) X ( 2 ) = X ^ ( 2 ) T = 0 g ( 2 ) ( X ( 1 ) , X ( 2 ) ) = ( φ ( X ( 1 ) , X ( 2 ) ) / X ( 2 ) ) X ( 1 ) = X ^ ( 1 ) X ( 2 ) = X ^ ( 2 ) T = 0
A T w ( 1 ) ( v ^ ( 2 ) ) v ^ ( 1 ) = A T w ( 1 ) ( v ^ ( 2 ) ) ( y A X ^ ( 1 ) ) = 0 A T w ( 2 ) ( v ^ ( 1 ) ) v ^ ( 2 ) = A T w ( 2 ) ( v ^ ( 1 ) ) ( y A X ^ ( 2 ) ) = 0
where g ( 1 ) ( X ( 1 ) , X ( 2 ) ) and g ( 2 ) ( X ( 1 ) , X ( 2 ) ) are the gradients of the function φ s q ( X ( 1 ) , X ( 2 ) ) . The following matrices
w ( 1 ) ( v ( 2 ) ) = Diag ( , w ( 1 ) ( v i ( 2 ) ) , ) ,   w ( 2 ) ( v ( 1 ) ) = Diag ( , w ( 2 ) ( v i ( 1 ) ) , )
are diagonal weight matrices that are based on the cross-weighting functions [20,31]
w ( 1 ) ( v i ( 2 ) ) = p i 2 v i ( 1 ) 2 v i ( 2 ) 2 2 v i ( 1 ) v i ( 1 ) = p i 2 v i ( 2 ) 2 ,   w ( 2 ) ( v i ( 1 ) ) = p i 2 v i ( 1 ) 2 v i ( 2 ) 2 2 v i ( 2 ) v i ( 2 ) = p i 2 v i ( 1 ) 2
The solutions of Equation (8) are the following Msplit estimators
X ^ ( 1 ) = D ( 1 ) ( v ^ ( 2 ) ) y , X ^ ( 2 ) = D ( 2 ) ( v ^ ( 1 ) ) y
where
D ( 1 ) ( v ^ ( 2 ) ) = [ A T w ( 1 ) ( v ^ ( 2 ) ) A ] 1 A T w ( 1 ) ( v ^ ( 2 ) ) ,   D ( 2 ) ( v ^ ( 1 ) ) = [ A T w ( 2 ) ( v ^ ( 1 ) ) A ] 1 A T w ( 2 ) ( v ^ ( 1 ) )
Thus, X ^ ( 1 ) is a function of v ^ ( 2 ) = y A X ^ ( 2 ) , whereas X ^ ( 2 ) is a function of v ^ ( 1 ) = y A X ^ ( 1 ) . For this reason, this solution has an asymptotic character. The following iterative procedure can be applied to compute the sought estimates ( j = 1 , , m )
X ( 1 ) j + 1 = D ( 1 ) ( v ( 2 ) j ) y , v ( 1 ) j + 1 = y A X ( 1 ) j + 1 X ( 2 ) j + 1 = D ( 2 ) ( v ( 1 ) j + 1 ) y , v ( 2 ) j + 1 = y A X ( 2 ) j + 1
(for the given starting point, for example, v ( 2 ) 0 = y A X ^ L S ). The process stops when for each l = 1 , 2 , it holds that g ( l ) ( X ^ ( 1 ) , X ^ ( 2 ) ) = 0 and hence X ^ ( l ) = X ( l ) m = X ( l ) m 1 . Note that other iterative processes that use both the gradients and the Hessians of φ ( X ( 1 ) , X ( 2 ) ) , namely Newton’s method, can be found in [19,20,29].
Now, the elementary property the of Msplit estimates is shown. Here, we consider a basic example that precedes the more detail analysis presented in the next section. Let us assume the functional model y i = X + v i , i = 1 , , 7 , and the observation set Φ as a following vector y = [ 1.2 0.9 1.8 1.3 2.2 1.1 1.9 ] T (Figure 1a). Then, X ^ L S = 1.49 . For the sake of comparison, let the robust M-estimate be computed. By applying the Huber method [27,28], where the weight function is w ( v ) = min { 1 ,   k / | v | } and k = 3 , one can obtain X ^ M = 1.27 (Figure 1b). Both estimates in question are not satisfactory, and do not reflect the nature of the observation set. The robust estimate X ^ M lies closer to the “bigger” aggregation of observations. Next, the question of how to treat the observations that are furthest from that estimate arises. In the classical approach, such observations are regarded as outliers (for example, affected by gross errors), and we are no longer interested in such observations. Different conclusions follow the Msplit estimation where X ^ ( 1 ) = 1.10 and X ^ ( 2 ) = 2.00 (Figure 1c).
The Msplit estimates show that set Φ consists of two subsets Φ 1 and Φ 2 (Figure 1d), whose elements can be regarded as realizations of two different random variables that differ from each other in location parameters X 1 and X 2 , respectively. Similar assumptions can also be found in other estimation problems, for example, cluster analysis (e.g., [32,33]); or in a mixed model estimation applied in geosciences (e.g., [34,35]). Such approaches can be regarded as alternatives; however, we should have some understanding that they differ significantly in their general ideas.
Assigning each observation to the model that is the most suitable for it is a natural process in Msplit estimation. This property can be applied in the analysis of network deformation where there are two functional models: y 1 = A 1 X 1 + v 1 and y 2 = A 2 X 2 + v 2 for two measurement epochs, respectively. Thus, one can create one common observation vector y = [ y 1 T , y 2 T ] T , the common weight matrix P = Diag ( P ( 1 ) , P ( 2 ) ) , and the coefficient matrix A = [ A 1 T , A 2 T ] T . It is noteworthy that the order of the observation within vector y can be arbitrary. The actual order of the observations must coincide with the order of the rows within matrix A and order of the weights in weight matrix P . Here, the shift Δ X ( 1 , 2 ) can be estimated by Δ X ^ ( 1 , 2 ) = X ^ ( 2 ) X ^ ( 1 ) . It is worth noting that Δ X ( 1 , 2 ) can also be estimated directly by applying the Shift-Msplit estimation proposed by Duchnowski and Wiśniewski [22].

3. Empirical Analyses

3.1. Elementary Tests

The elementary analysis was based on the univariate models and simulations of observations related to such models. Thus,
y i , 1 = X 1 + v i , 1 ,   i = 1 , , n 1 y i , 2 = X 2 + v i , 2 ,   i = 1 , , n 2 y 1 = 1 n 1 X 1 + v 1 y 2 = 1 n 2 X 2 + v 2
where 1 n l = [ 1 1 , , 1 n l ] T ; X 1 and X 2 are parameters that differ from each other in the shift Δ X ( 1 , 2 ) = X 2 X 1 . The measurements, namely the elements of vectors y 1 and y 2 , were simulated by using the Gaussian generator r a n d n ( n , 1 ) of MATLAB. We assumed that σ = 1 , and the following theoretical values of the parameters: X 1 t = 0 and hence X 2 t = X 1 t + Δ X ( 1 , 2 ) = Δ X ( 1 , 2 ) . Considering the LS-estimation of X 1 and X 2 we can apply the model of Equation (14) or Equation (1) where A 1 = 1 n 1 and A 2 = 1 n 2 . In the case of Msplit estimation, we assumed the model of Equation (3), taking y = [ y 1 T , y 2 T ] T R n , n = n 1 + n 2 , and A = [ 1 n 1 T , 1 n 2 T ] T = 1 n . We also applied the iterative procedure of Equation (13) by taking LS-estimates as the starting point (note that the starting point can usually be arbitrary).
Let us now consider an example of observation simulation for which Δ X ( 1 , 2 ) = 5 σ = 5 and n 1 = 50 , n 2 = 10 . The parameter estimates, together with the respective residuals, are presented in Figure 2.
Now, let us consider more simulated observation sets. By applying the crude Monte Carlo method (MC) for N simulations, one can compute the MC estimates by applying the formula
θ ^ M C = 1 N i = 1 N θ ^ i
where θ ^ i are the estimates obtained for the ith simulation. The location of the MC estimates for N = 5000 and Δ X ( 1 , 2 ) = 5 or Δ X ( 1 , 2 ) = 20 is presented in Figure 3.
This shows that the MC estimates that were obtained for both estimation methods were close to the respective theoretical values (considering the simulated standard deviation). Generally, the LS estimates seemed more satisfactory. Please note that the results obtained for different values of shift Δ X ( 1 , 2 ) indicate that Msplit estimation is more satisfactory for bigger shifts than for smaller ones. Thus, let us examine how efficient the Msplit estimation is for different shifts.
Let the measure of efficacy be defined in relation to the LS estimates, thus
λ ( l ) ( X ^ ( l ) , X ^ L S , l ) = a b s ( X ^ ( l ) X l t ) a b s ( X ^ L S , l X l t )
Note that when λ ( l ) ( X ^ ( l ) , X ^ L S , l ) < 0 , then the Msplit estimate is closer to the theoretical value than the LS estimate. Now, we can define the following function of an elementary success of Msplit estimation
s ( l ) ( X ^ ( l ) , X ^ L S , l ) = { 1 f o r   λ ( X ^ ( l ) , X ^ L S , l ) < 0 0 f o r   λ ( X ^ ( l ) , X ^ L S , l ) > 0
The application of MC simulations allowed us to present the success rate (SR), which can be computed for different values of the shift Δ X ( 1 , 2 )
γ ( l ) ( X ^ ( l ) , X ^ L S , l ; Δ X ( 1 , 2 ) ) = 1 N i = 1 N s ( l ) i ( X ^ ( l ) , X ^ L S , l )
where s ( l ) i ( X ^ ( l ) , X ^ L S , l ) is the value of Equation (17) at the ith simulation.
Note that such a SR is defined in a very similar way to the mean success rate (MSR) given by Hekimoglu and Koch [36]. SRs for different Δ X ( 1 , 2 ) and for N = 5000 simulations are presented in Figure 4.

3.2. Vertical Displacement Analysis

Let us now consider the efficacy of Msplit estimates in a more practical example, namely the analysis of vertical displacements within the leveling network, which is presented in Figure 5. Such a network has already been under investigation in previous papers [24,25].
The network consists of four reference points R 1 , , R 4 with the known heights H R 1 = = H R 4 = 0 m and five object points of P 1 , , P 5 . We assumed that each of the height differences h 1 , , h 16 was measured twice at each of two measurement epochs, and that σ = 2 mm was the known standard deviation of all measurements. We also assumed that at the first epoch X 1 t = [ H 1 , 1 = 0 , , H 5 , 1 = 0 ] T = 0 , where H i , 1 is the height of the ith object point at the first epoch. The shift of the object points between the measurement epochs is given by Δ X ( 1 , 2 ) = [ Δ H 1 ( 1 , 2 ) , , Δ H 5 ( 1 , 2 ) ] T , where Δ H i ( 1 , 2 ) = H i , 2 H i , 1 . In the classical approach to the estimation of the point displacements, we used the functional model of Equation (1). Since all height differences were measured twice at two measurement epochs, namely, we had two series of measurements at each epoch, then we should assume that y l R 32 , X l = [ H 1 , l , , H 5 , l ] T , and A = A 1 2 R 32 , 5 where A R 16 , 5 is a known coefficient matrix related to one series of measurements, 1 2 = [ 1 , 1 ] T , and is the Kronecker product. On the other hand, in the case of Msplit estimation, we should apply the functional model of Equation (3) for which y = [ y 1 T , y 2 T ] T R 64 , A = [ A T , A T ] T R 64 , 5 , and X ( 1 ) , X ( 2 ) R 5 are the competitive versions of the parameter vector, hence v ( 1 ) , v ( 2 ) R 64 are the respective competitive versions of the measurement errors.
When analyzing the efficacy of Msplit estimation, we can use two measures, namely the local measure of the distance between the LS and Msplit estimates
λ ( l ) j ( [ X ^ ( l ) ] j , [ X ^ L S , l ] j ) = = a b s ( [ X ^ ( l ) ] j [ X l t ] j ) a b s ( [ X ^ L S , l ] j [ X l t ] j )
as well as the global one
λ ( l ) ( X ^ ( l ) , X ^ L S , l ) = X ^ ( l ) X l t X ^ L S , l X l t
where [ ] j is jth element of the vector and is the Euclidean norm. The local distance, which is just another form of Equation (16), is related to a particular parameter, for example, the height of a displacing point. The global distance describes the whole parameter vector. Thus, we can define the local and global success rates in the following way
γ ( l ) , j ( [ X ^ ( l ) ] j , [ X ^ L S , l ] j ; Δ X ( 1 , 2 ) ) = 1 N i = 1 N s ( l ) , j i ( [ X ^ ( l ) ] j , [ X ^ L S , l ] j ) γ ( l ) ( X ^ ( l ) , X ^ L S , l ; Δ X ( 1 , 2 ) ) = 1 N i = 1 N s ( l ) i ( X ^ ( l ) , X ^ L S , l )
where s ( l ) j i ( [ X ^ ( l ) ] j , [ X ^ L S , l ] j ) and s ( l ) i ( X ^ ( l ) , X ^ L S , l ) are functions of an elementary success from Equation (17) and indexed with the respective arguments.
The empirical analysis, which was based on the MC method for N = 5000 simulations, was carried out for several variants of the point displacements. First, we assumed that only point P 5 was displaced. The respective MC estimates obtained for the LS and Msplit estimations and Δ H 5 ( 1 , 2 ) = 50 , Δ H 5 ( 1 , 2 ) = 100 , or Δ H 5 ( 1 , 2 ) = 200 mm are presented in Table 1, which also presents the local and global SRs.
The MC estimates were similar for both estimation methods and the stable points. The SRs indicate that the LS estimates were closer to the theoretical values in the vast majority of the simulations. Note that the local SRs obtained for point P 5 were much higher than the global ones. All estimates of the point heights obtained in the MC simulations (for the variant Δ H 5 ( 1 , 2 ) = 50 mm) are presented in Figure 6.
In the second variant, we assumed that there were two unstable points, namely P 5 and P 4 . The results, which were obtained for the different point shifts, are presented in Table 2. Here, the MC estimates obtained for both methods were also similar. Figure 7 presents the LS and Msplit estimates that were obtained for all of the MC simulations. Generally, this confirmed the correctness of both estimation methods; however, differences between these two estimation methods were also apparent. The main difference was the dispersion, which was larger in the case of the Msplit estimation, especially for the stable points, which suggests that the accuracy of the Msplit estimation was worse than LS estimation. It is also worth noting that the SRs of the Msplit estimation achieved bigger values in this variant. In the case of point P 5 , the results of the Msplit estimation were better than the results of the classical approach in almost one third of the simulations.
The results, which are presented here, show that both methods, namely LS and Msplit estimation, yielded satisfactory solutions. However, such a conclusion was valid for the ordered observation sets, namely when each observation was properly assigned to its measurement epoch. If such a condition is not met, then the observation from another epoch will usually be regarded as an outlier. Since LS estimation as well as Msplit estimation are not robust against outliers, they both break down (please note that Msplit estimation is generally not robust unless we introduce an additional virtual model for outliers). Note that in the context addressed here, the outliers result from the assignment of an observation to the wrong measurement epoch, but not from gross errors. The natural feature of Msplit estimation is the automatic assignment of each observation to the proper epoch. Thus, we can suppose that this estimation method will not break down if such outliers occur. To illustrate this feature of Msplit estimation, we simulated that point P 5 was displaced and that Δ H 5 ( 1 , 2 ) = 50 mm. Now, let us consider the following variants of the observation sets: variant A, where both observation sets were correct (all observations were assigned to their epochs properly); variant B, where the observation h 16 at the second epoch was equal to h 16 at the first one, namely h 16 2 = h 16 1 ; and variant C, where h 16 2 = h 16 1 , but also h 15 2 = h 15 1 . Thus, in variants B and C, we simulated that some observations that were assigned to the second measurement epoch should be related to the first one. The results obtained for all variants are presented in Table 3. In the case of variant A, the results were very close to the respective results presented in Table 1. If the observation sets are not ordered correctly, then the local SRs at the second epoch are close to 1, which means that almost always, the height of point P 5 at the second measurement epoch is better assessed by the Msplit estimation than by LS estimation. Additionally, the global SRs were very high at the second epoch, hence one can say that the heights of all network points were better estimated by the application of Msplit estimation.

4. Conclusions

The paper showed that Msplit estimation can be successfully applied in deformation analysis. The results were generally similar to the results of the more conventional LS estimation; however, the latter method usually yielded slightly better outcomes. The elementary tests showed that the efficacy of the Msplit estimation grew with an increasing shift between the observation sets. In the case of geodetic networks, where a parameter vector usually consists of several point coordinates, the shift of one or two such coordinates between measurement epochs does not influence the efficacy of the Msplit estimation in a significant way. The real advantage of the Msplit estimation was revealed for the disordered observation sets, for example, when the observations from at least two measurement epochs were mixed for some reason. Note that the LS estimates break down in such cases, in contrast with the Msplit estimation, for which the ordering of all observations within the combined observation set can be arbitrary and does not influence the final results of the method as well as its iterative process. Such a feature results directly from the theoretical foundations of the method, which are based on the concept of the split potential. In short, each observation “chooses” the functional model that fits it best. In this context, Msplit estimates are robust against some kind of “outliers”, namely observations that come from other observation sets. Referring to the presented example, there were four height differences regarding the height of network point P 5 . If one of them does not fit the other, then the method tries to fit such an “outlying” observation into another epoch. If it works, then the whole estimation process succeeds. However, if such an observation is in fact affected by a gross error, then it does not fit any epoch, and the estimation must break down. The introduction of a virtual epoch, which is not related to any real measurements, is one solution to this problem. One can say that such an epoch can collect all “loners” that do not fit any real measurement epochs. Generally speaking, one can say that the Msplit estimation is not robust against outliers, which results from the occurrence of gross errors. However, if one assumes an additional competitive functional model (dedicated to outliers), then the Msplit estimation can estimate the location parameters for “good” observation aggregations as well as outlier(s). Increasing the number of competitive functional models protects the estimation of location parameters of good observations from the bad influence of outliers. Note that in this context, outliers are no longer “outlying”, and become regular observations of the third (or more generally next) aggregation. This concept, which is out of the scope of this paper, was discussed in [23,24,30].

Author Contributions

Conceptualization, Z.W. and R.D.; Methodology, Z.W.; Software, Z.W and R.D. Validation Z.W. and R.D.; Formal analysis, Z.W.; Investigation, R.D. and A.D.; Writing—original draft preparation, Z.W. and R.D.; Writing—review and editing, R.D. and A.D.; Visualization, A.D.; Supervision, Z.W.

Funding

This research was funded by the Institute of Geodesy, University of Warmia and Mazury in Olsztyn, statutory research no. 28.610.002-300.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Pelzer, H. Zur Analyse Geodätischer Deformationsmessungen; Deutsche Geodätische Kommission: Munich, Germany, 1971. [Google Scholar]
  2. Caspary, W.F.; Haen, W.; Borutta, H. Deformation analysis by statistical methods. Technometrics 1990, 32, 49–57. [Google Scholar] [CrossRef]
  3. Hekimoglu, S.; Erdogan, B.; Butterworth, S. Increasing the efficacy of the conventional deformation analysis methods: Alternative strategy. J. Surv. Eng. 2010, 136, 1–8. [Google Scholar] [CrossRef]
  4. Niemeier, W. Statistical tests for detecting movements in repeatedly measured geodetic networks. In Developments in Geotectonics; Elsevier: Amsterdam, The Netherlands, 1981; Volume 71, pp. 335–351. [Google Scholar]
  5. Setan, H.; Singh, R. Deformation analysis of a geodetic monitoring network. Geomatica 2001, 55, 333–346. [Google Scholar]
  6. Denli, H.H.; Deniz, R. Global congruency test methods for GPS networks. J. Surv. Eng. 2003, 129, 95–98. [Google Scholar] [CrossRef]
  7. Chen, Y.Q. Analysis of Deformation Surveys—A Generalized Method; Technical Report; UNB Geodesy and Geomatics Engineering, University of New Brunswick: Fredericton, NB, Canada, 1983. [Google Scholar]
  8. Caspary, W.F.; Borutta, H. Robust estimation in deformation models. Surv. Rev. 1987, 29, 29–45. [Google Scholar] [CrossRef]
  9. Duchnowski, R. Median-based estimates and their application in controlling reference mark stability. J. Surv. Eng. 2010, 136, 47–52. [Google Scholar] [CrossRef]
  10. Duchnowski, R. Hodges–Lehmann estimates in deformation analyses. J. Geod. 2013, 87, 873–884. [Google Scholar] [CrossRef]
  11. Duchnowski, R.; Wiśniewski, Z. Comparison of two unconventional methods of estimation applied to determine network point displacement. Surv. Rev. 2014, 46, 401–405. [Google Scholar] [CrossRef]
  12. Duchnowski, R.; Wiśniewski, Z. Accuracy of the Hodges-Lehmann estimates computed by applying Monte Carlo simulations. Acta Geod. Geophys. 2017, 52, 511–525. [Google Scholar] [CrossRef]
  13. Duchnowski, R.; Wiśniewski, Z. Msplit and Mp estimation. A wider range of robustness. In Proceedings of the International Conference on Environmental Engineering, Vilnius, Lithuania, 27–28 April 2017; pp. 1–6. [Google Scholar]
  14. Wyszkowska, P.; Duchnowski, R. Subjective breakdown points of R-estimators applied in deformation analysis. In Proceedings of the International Conference on Environmental Engineering, Vilnius, Lithuania, 27–28 April 2017; pp. 1–6. [Google Scholar]
  15. Erdogan, B.; Hekimoglu, S. Effect of subnetwork configuration design on deformation analysis. Surv. Rev. 2014, 46, 142–148. [Google Scholar] [CrossRef]
  16. Nowel, K.; Kamiński, W. Robust estimation of deformation from observation differences for free control networks. J. Geod. 2014, 88, 749–764. [Google Scholar] [CrossRef]
  17. Nowel, K. Robust M-Estimation in analysis of control network deformations: classical and new method. J. Surv. Eng. 2015, 141, 04015002. [Google Scholar] [CrossRef]
  18. Amiri-Simkooei, A.R.; Alaei-Tabatabaei, S.M.; Zangeneh-Nejad, F.; Voosoghi, B. Stability analysis of deformation-monitoring network points using simultaneous observation adjustment of two epochs. J. Surv. Eng. 2017, 143, 04016020. [Google Scholar] [CrossRef]
  19. Wiśniewski, Z. Estimation of parameters in a split functional model of geodetic observations (Msplit estimation). J. Geod. 2009, 83, 105–120. [Google Scholar] [CrossRef]
  20. Wiśniewski, Z. Msplit(q) estimation: Estimation of parameters in a multi split functional model of geodetic observations. J. Geod. 2010, 84, 355–372. [Google Scholar] [CrossRef]
  21. Janowski, A.; Rapiński, J. M–Split Estimation in Laser Scanning Data Modeling. J Indian Soc. Remote Sens. 2013, 41, 15–19. [Google Scholar] [CrossRef]
  22. Duchnowski, R.; Wiśniewski, Z. Estimation of the shift between parameters of functional models of geodetic observations by applying Msplit estimation. J. Surv. Eng. 2011, 138, 1–8. [Google Scholar] [CrossRef]
  23. Zienkiewicz, M.H. Application of Msplit estimation to determine control points displacements in networks with unstable reference system. Surv. Rev. 2015, 47, 174–180. [Google Scholar] [CrossRef]
  24. Wiśniewski, Z.; Zienkiewicz, M.H. Shift- M split estimation in deformation analyses. J. Surv. Eng. 2016, 142, 04016015. [Google Scholar] [CrossRef]
  25. Velsink, H. Testing methods for adjustment models with constraints. J. Surv. Eng. 2018, 144, 04018009. [Google Scholar] [CrossRef]
  26. Li, J.; Wang, A.; Wang, X. Msplit estimate the relationship between LS and its application in gross error detection. Mine Surv. (China) 2013, 2, 57–59. [Google Scholar]
  27. Huber, P.J. Robust estimation of location parameter. In Breakthroughs in Statistics; Springer: Berlin/Heidelberg, Germany, 1992; pp. 492–518. [Google Scholar]
  28. Huber, P.J. Robust Statistics; Springer: Berlin/Heidelberg, Germany, 2011. [Google Scholar]
  29. Wyszkowska, P.; Duchnowski, R. Msplit estimation based on L1 norm condition. J. Surv. Eng. 2019, 145, 04019006. [Google Scholar] [CrossRef]
  30. Zienkiewicz, M.H. Determination of an adequate number of competitive functional models in the square Msplit(q) estimation with the use of a modified Baarda’s approach. Surv. Rev. 2018, 1–11. [Google Scholar] [CrossRef]
  31. Duchnowski, R.; Wiśniewski, Z. Robustness of Msplit(q) estimation: A theoretical approach. Stud. Geophys. Geod. 2019, 63, 390–417. [Google Scholar] [CrossRef]
  32. Sebert, D.M.; Montgomery, D.C.; Rollier, D.A. A clustering algorithm for identifying multiple outliers in linear regression. Comput. Stat. Data Anal. 1998, 27, 461–484. [Google Scholar] [CrossRef]
  33. Soto, J.; Vigo Aguiar, M.I.; Flores-Sintas, A. A fuzzy clustering application to precise orbit determination. J. Comput. Appl. Math. 2007, 204, 137–143. [Google Scholar] [CrossRef]
  34. Spurr, B.D. On estimating the parameters in mixtures of circular normal distributions. J. Int. Assoc. Math. Geol. 1981, 13, 163–173. [Google Scholar] [CrossRef]
  35. Hsu, J.S.; Walker, J.J.; Orgen, D.E. A stepwise method for determining the number of component distributions in a mixture. Math. Geol. 1986, 18, 153–160. [Google Scholar] [CrossRef]
  36. Hekimoglu, S.; Koch, K.R. How can reliability of the test for outliers be measured? Allg. Vermes. Nachr. 2000, 7, 247–253. [Google Scholar]
Figure 1. Least squares, robust M-estimate and both Msplit estimates within a sample observation set (a) observation set; (b) classical estimates; (c) Msplit estimates; (d) observation subsets.
Figure 1. Least squares, robust M-estimate and both Msplit estimates within a sample observation set (a) observation set; (b) classical estimates; (c) Msplit estimates; (d) observation subsets.
Sensors 19 05047 g001
Figure 2. LS estimates, Msplit estimates, and respective residuals (elementary functional models for Δ X ( 1 , 2 ) = 5 ).
Figure 2. LS estimates, Msplit estimates, and respective residuals (elementary functional models for Δ X ( 1 , 2 ) = 5 ).
Sensors 19 05047 g002
Figure 3. Location of the Monte Carlo estimates for Δ X ( 1 , 2 ) = 5 or Δ X ( 1 , 2 ) = 20 (for N = 5000 ).
Figure 3. Location of the Monte Carlo estimates for Δ X ( 1 , 2 ) = 5 or Δ X ( 1 , 2 ) = 20 (for N = 5000 ).
Sensors 19 05047 g003
Figure 4. The success rates of Msplit estimates X ^ ( 1 ) and X ^ ( 2 ) for the growing value of Δ X ( 1 , 2 ) .
Figure 4. The success rates of Msplit estimates X ^ ( 1 ) and X ^ ( 2 ) for the growing value of Δ X ( 1 , 2 ) .
Sensors 19 05047 g004
Figure 5. Tested leveling network.
Figure 5. Tested leveling network.
Sensors 19 05047 g005
Figure 6. The LS and Msplit estimates of the MC simulations ( Δ H 5 ( 1 , 2 ) = 50 mm).
Figure 6. The LS and Msplit estimates of the MC simulations ( Δ H 5 ( 1 , 2 ) = 50 mm).
Sensors 19 05047 g006
Figure 7. The LS and Msplit estimates of the MC simulations ( Δ H 4 ( 1 , 2 ) = 50 and Δ H 5 ( 1 , 2 ) = 100 mm).
Figure 7. The LS and Msplit estimates of the MC simulations ( Δ H 4 ( 1 , 2 ) = 50 and Δ H 5 ( 1 , 2 ) = 100 mm).
Sensors 19 05047 g007
Table 1. The Monte Carlo estimates of the point heights and success rates for one unstable point.
Table 1. The Monte Carlo estimates of the point heights and success rates for one unstable point.
Δ H 5 ( 1 , 2 ) = 50 Δ H 5 ( 1 , 2 ) = 100 Δ H 5 ( 1 , 2 ) = 200
X ^ L S , 1 X ^ ( 1 ) X ^ L S , 2 X ^ ( 2 ) X ^ L S , 1 X ^ ( 1 ) X ^ L S , 2 X ^ ( 2 ) X ^ L S , 1 X ^ ( 1 ) X ^ L S , 2 X ^ ( 2 )
0.2−3.10.42.9−0.5−1.1−0.60.60.9−1.8−0.7−0.4
1.4−1.2−1.00.9−0.4−1.30.72.90.50.70.5−0.8
2.1−0.6−0.6−0.60.1−0.60.51.3−0.3−2.8−0.11.1
-0.8−3.6−0.60.61.1−0.9−1.01.20.0−0.4−1.51.5
0.8−1.9−50.4−49.10.3−1.2−99.8−98.70.8−0.8−200.1−199.7
γ ( 1 ) = 0.018
γ ( 1 ) , 5 = 0.172
γ ( 2 ) = 0.019
γ ( 2 ) , 5 = 0.182
γ ( 1 ) = 0.020
γ ( 1 ) , 5 = 0.177
γ ( 2 ) = 0.017
γ ( 2 ) , 5 = 0.187
γ ( 1 ) = 0.025
γ ( 1 ) , 5 = 0.171
γ ( 2 ) = 0.024
γ ( 2 ) , 5 = 0.196
Table 2. The MC estimates of the point heights and SRs for two unstable points.
Table 2. The MC estimates of the point heights and SRs for two unstable points.
Δ H 5 ( 1 , 2 ) = 50 ;   Δ H 4 ( 1 , 2 ) = 50 Δ H 5 ( 1 , 2 ) = 100 ;   Δ H 4 ( 1 , 2 ) = 50 Δ H 5 ( 1 , 2 ) = 200 ;   Δ H 4 ( 1 , 2 ) = 50
X ^ L S , 1 X ^ ( 1 ) X ^ L S , 2 X ^ ( 2 ) X ^ L S , 1 X ^ ( 1 ) X ^ L S , 2 X ^ ( 2 ) X ^ L S , 1 X ^ ( 1 ) X ^ L S , 2 X ^ ( 2 )
−0.2−0.5−0.30.6−0.50.1−0.40.30.0−1.3−0.3−1.1
−0.4−2.0−0.1−0.10.20.5−0.2−0.40.0−0.2−0.10.9
−0.4−0.4−0.90.20.1−0.30.2−0.3−0.2−0.2−0.3−0.4
−0.1−0.3−50.5−50.10.40.5−49.9−49.6−0.1−0.4−50.0−50.1
−0.5−1.4−50.1−50.2−0.6−0.4−100.1−99.8−0.5−0.8−200.3−200.2
γ ( 1 ) = 0.070
γ ( 1 ) , 5 = 0.272
γ ( 2 ) = 0.070
γ ( 2 ) , 5 = 0.268
γ ( 1 ) = 0.080
γ ( 1 ) , 5 = 0.281
γ ( 2 ) = 0.080
γ ( 2 ) , 5 = 0.288
γ ( 1 ) = 0.103
γ ( 1 ) , 5 = 0.314
γ ( 2 ) = 0.105
γ ( 2 ) , 5 = 0.312
Table 3. The MC estimates of the point heights and SRs for the disturbed observation sets.
Table 3. The MC estimates of the point heights and SRs for the disturbed observation sets.
Variant A: Correct Order Variant   B :   h 16 2 = h 16 1 Variant   C :   h 15 2 = h 15 1 ,   h 16 2 = h 16 1
X ^ L S , 1 X ^ ( 1 ) X ^ L S , 2 X ^ ( 2 ) X ^ L S , 1 X ^ ( 1 ) X ^ L S , 2 X ^ ( 2 ) X ^ L S , 1 X ^ ( 1 ) X ^ L S , 2 X ^ ( 2 )
0.02.20.3−1.10.4−1.5−6.80.3−0.80.4−4.5−5.2
0.4−0.11.10.4−0.5−1.52.11.8−0.2−0.8−5.3−7.7
0.60.80.3−1.5−0.6−3.63.41.6−0.1−1.04.97.4
−0.7−0.90.01.00.3−1.42.02.4−1.3−0.65.27.1
−0.20.5−49.8−50.30.4−1.5−36.2−46.5−2.0−1.025.3−42.6
γ ( 1 ) = 0.018
γ ( 1 ) , 5 = 0.172
γ ( 2 ) = 0.019
γ ( 2 ) , 5 = 0.210
γ ( 1 ) = 0.127
γ ( 1 ) , 5 = 0.321
γ ( 2 ) = 0.875
γ ( 2 ) , 5 = 0.986
γ ( 1 ) = 0.263
γ ( 1 ) , 5 = 0.474
γ ( 2 ) = 0.887
γ ( 2 ) , 5 = 0.998

Share and Cite

MDPI and ACS Style

Wiśniewski, Z.; Duchnowski, R.; Dumalski, A. Efficacy of Msplit Estimation in Displacement Analysis. Sensors 2019, 19, 5047. https://doi.org/10.3390/s19225047

AMA Style

Wiśniewski Z, Duchnowski R, Dumalski A. Efficacy of Msplit Estimation in Displacement Analysis. Sensors. 2019; 19(22):5047. https://doi.org/10.3390/s19225047

Chicago/Turabian Style

Wiśniewski, Zbigniew, Robert Duchnowski, and Andrzej Dumalski. 2019. "Efficacy of Msplit Estimation in Displacement Analysis" Sensors 19, no. 22: 5047. https://doi.org/10.3390/s19225047

APA Style

Wiśniewski, Z., Duchnowski, R., & Dumalski, A. (2019). Efficacy of Msplit Estimation in Displacement Analysis. Sensors, 19(22), 5047. https://doi.org/10.3390/s19225047

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop