# Robust Prediction Algorithm Based on a General EIV Model for Multiframe Transformation

Authors: Lihui Yao; Peng Lin; Jingxiang Gao; Chao Liu
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

In modern geodesy, there are cases in which the target frame is unique and there is more than one source frame. Helmert transformations, which are extensively used to solve transformation parameters, can be separately solved between the target frame and one of the source frames. However, this is not globally optimal, even though each transformation is locally optimal on its own. Additionally, this also generates the problem of multiple solutions in the noncommon station of the target frame. Moreover, least squares solutions can cause estimation value distortion, with a gross error existing in observations. Thus, in this paper, Helmert transformations among three frames, that is, one target frame and two source frames, are studied as an example. A robust prediction algorithm based on the general errors-in-variables prediction algorithm and the robust estimation is derived in detail and is applied to achieve multiframe total transformation. Furthermore, simulation experiments were conducted and the results validated the superiority of the proposed total transformation method over classical separate approaches.


---

## Body

## 1. Introduction

The precision and reliability of terrestrial reference frames have improved with the rapid development of the Global Navigation Satellite System (GNSS) [1]. At present, common terrestrial reference frames include WGS-84, PE-90, CGCS-2000, and ITRF2008. The three-dimensional (3D) space coordinates of ground control points in the source coordinate system are obtained using GNSS observations, whereas local coordinate systems are often used for geodetic surveying. Therefore, it is necessary to use coordinate transformation to unify the data. Similarity datum transformation is widely used in 3D coordinate transformation. Common point coordinates are typically used to solve the transformation parameters, and then the coordinates of noncommon points in the old coordinate system are converted to the new coordinate system according to the transformation parameter solution.Seven parameters of 3D datum transformations can be used to relate two reference frames: three for translation, one for scale, and three for rotation. Using two sets of coordinates measured in two frames, these parameters can be solved using the Gauss–Markov (GM) model of least squares (LS) [2]. The assumption of the GM model is that there are no random errors in the coefficient matrix, but the observed coordinates of the reference station in the source coordinate system of the coefficient matrix have random errors; thus, the GM model cannot be used to describe the coordinate transformation problem accurately. For the case in which the coefficient matrix contains random errors, the error-in-variables (EIV) model can be used to describe the GM model, and the LS solution of the EIV model is called the total least squares (TLS).In recent years, TLS has become a research topic of great interest in the field of geodesy and many algorithms have been developed and examined [3–6]. When the small rotation angle and scale are close to one, the datum transformation model can adopt the Bursa–Wolf model. Compared with classical LS, the accuracy of estimation results using TLS is greatly improved [7–10]. In the case of a large rotation angle and arbitrary scale, Felus et al. proposed an algorithm for a multiple EIV model [11]. Fang transformed the EIV model into a nonlinear Gauss–Markov model and proposed the 3D datum transformation algorithm under any rotation angle by taking the coordinates of common points in the source coordinate system as virtual observation values [12]. Chang and Li et al. proposed the analytical and closed-form solutions of 3D datum transformation, respectively [13, 14]. Because the rotation matrix of the datum transformation is a completely orthogonal matrix, there are strict orthogonal constraints. Fang proposed the universal constraint model of weighted TLS [5]. Lin et al. proposed the constraint solution of 3D datum transformation based on the nonlinear Gauss–Helmert model [15].However, most TLS algorithms for datum transformation that focus on how to solve transformation parameters do not combine the two independent processes of transformation parameter estimation and noncommon station transformation; that is, noncommon station transformation is performed at the same time as the estimation of transformation parameters. In practice, there is a statistically significant correlation between the common and noncommon stations [16]; for example, the coordinates of common and noncommon stations come from the same control network adjustment results (e.g., GPS baseline network), but the existing model algorithm does not take into consideration the correlation between reference stations, so the estimation results cannot obtain the optimal solution in theory. To overcome this problem, Li et al. [17, 18] proposed the concept of seamless frame datum transformation, which combines the estimation parameters and coordinate transformation of noncommon stations. Compared with the traditional model, this transformation effectively improves the accuracy of the transformation coordinates. Wang et al. proposed the estimation model of 3D datum transformation under a large rotation angle and arbitrary scale conditions based on the generalized EIV model [19]. Lin et al. proposed the prediction model of 3D datum transformation of source coordinates in the form of geodetic coordinates based on the nonlinear Gauss–Helmert model [20].Most of the current literature mainly concentrates on the datum transformation between two frames. When multiple (two or more) frames are involved, the transformation between the multiple source frames and a target frame, as mentioned earlier, clearly can be used to convert each source frame in turn with the target frame. The advantage is that we can draw lessons from the existing research achievements of Helmert transformation between the two frames directly. However, the existing methods cannot obtain the optimal solution. In recent years, a few researchers have studied the problem of simultaneous Helmert transformation with multiple frames. For example, Han et al. and Aktug proposed the closed solutions of this problem from the perspectives of geometry and algebra, respectively [21, 22], but the solutions all lack statistical optimality; in other words, the solutions are not LS.The parameter estimation results of TLS are seriously affected when the observation of source frames contains gross errors. Mahboub et al. proposed a method of robust TLS based on selecting weight iteration for the EIV model [23]. Lu et al. and Tao et al. studied the robust TLS (RTLS) algorithm for 3D datum transformation and GPS height fitting, respectively, based on the Gauss–Helmert model [24, 25]. Jin et al. and Zhao et al. studied the total minimum L1 norm estimation method, taking into account the error of all variables [26, 27]. Wang et al. proposed the antidifference solution of weighted TLS (WTLS) by standardizing the residual to construct the equal-weight function and using the median method [6]. To solve the problem of processing outliers in universal 3D symmetrical Helmert transformation (probably with large rotation angles), Chang et al. derived the M-estimator by introducing the reweighted LS scheme [28]; in addition, the data snooping procedure was extended to the generalized EIV model by Wang et al. [29], consequently solving this problem by taking another approach.Continuing from the analysis given above, we propose a robust predicted estimation for multiframe transformation, which constructs the weight factor function with standardized residuals and uses the median method to obtain the unit weight mean error. We need to consider the random errors of all observations and their correlation to finally achieve multiframe total transformation based on the generalized EIV model.


## 2. Extrapolation Estimation Method Based on the Generalized EIV Model

The function model for the extrapolation estimation method based on the generalized EIV model [19, 20] is as follows:(1)L1-eL1=fa1-ea1,ξ(2)L-2=fa2-ea2,ξ,where the n × 1 vector L1 denotes measurements and eL1 is the corresponding measurement error; the m × 1 vector L-2 denotes estimated values; the n × 1 vector a1 and m × 1 vector a2 are the measurements of the common and noncommon points, respectively, in the old frame; ea1 and ea2 are the corresponding measurement errors; ξ is the t × 1 vector of estimated parameters; and f is an abstract function.The corresponding stochastic model can be described as(3)e~N0,σ02Q(4)e=eL1ea1ea2,Q=QL1L1QL1a1QL1a2Qa1L1Qa1a1Qa1a2Qa2L1Qa2a1Qa2a2,where QL1L1,Qa1a1, and Qa2a2 are the cofactor matrices of L1, a1, and a2, respectively; Qa1L1, Qa2L1, and Qa2a1 are the intercofactor matrices of L1, a1, and a2, respectively.Taylor series expansion is applied to the right-hand side of (1) and (2) at the approximate values ξ0, ea10, and ea20 of ξ, ea1, and ea2, respectively. Thus, (1) and (2) are now expressed as(5)L1-eL1=fa1-ea10,ξ0+A1dξ+B1ea1-ea10(6)L-20+dL-2=fa2-ea20,ξ0+A2dξ+B2ea2-ea20,where dξ denotes the correction vector of parameters, dξ=ξ-ξ0,(7)A1=∂f∂ξTea10,ξ0,B1=∂f∂ea1Tea10,ξ0,A2=∂f∂ξTea20,ξ0,B2=∂f∂ea2Tea20,ξ0.Then the standard LS form of (5) and (6) is described as follows: (8)l-Ge=Adθwhere l=L1-f(a1-ea10,ξ0)+B1ea10L-20-f(a2-ea20,ξ0)+B2ea20, G=I3B1000B2, A=A10A2-Im, dθ=dξdL-2.On the basis of the weighted LS error criterion, the LS solutions of the parameters and the estimation of the error vector can be calculated using (8):(9)dθ^=ATQG-1A-1ATQG-1l(10)e^=QGTQG-1l-Adθ^,where QG=GQGT.


## 3. Robust Estimation Method Based on the Standardized Residuals and Median

The accurate estimated parameters are difficult to obtain using (8) if there are gross errors in the observation vector or coefficient matrix A. For this, we propose a robust prediction algorithm based on the generalized EIV model. We use the robust estimation introduced by Wang et al. [6]; it is based on the standardized residuals and median robust model. The robust prediction model is more reasonable in theory than the existing TLS robust methods and has two advantages: the weight factor function constructed by the standardized residuals can take into account the robustness of both the observation space and structure space; the estimation of the unit weight mean error obtained by the median method has greater robustness in the iterative process. The estimation criterion of the robust prediction algorithm is described as follows:(11)eTQ--1e=min,where (12)Q-=Q-L1L1Q-L1a1Q-L1a2Q-a1L1Q-a1a1Q-a1a2Q-a2L1Q-a2a1Q-a2a2=P--1,where Q-L1L1, Q-a1a1, and Q-a2a2 are the equivalent cofactor matrices of observation vectors L1, a1, and a2, respectively; Q-a1L1, Q-a2L1, and Q-a2a1 are the equivalent intercofactor matrices of observation vectors L1, a1, and a2, respectively.Then, the solutions of the parameters and the estimation of the error vector are as follows:(13)dθ^=ATQ-G-1A-1ATQ-G-1l(14)e^=Q-GTQ-G-1l-Adθ^.Because the correlation between the observations is taken into account in the stochastic model of the extrapolation estimation method, it is necessary to construct the equivalent cofactor matrices in the form of double factors as follows [6, 30]:(15a)q-i=qiRii(15b)q-ij=qijRiiRjj,where the subscripts i and j represent the serial numbers of the observed values; qi and qij denote the elements in cofactor matrix Q; q-ij and q-ij denote the elements in equivalent cofactor matrix Q-; and Rii and Rjj are the reciprocals of the weight factors, which are called the cofactor factors [30]. Additionally, because the generalized EIV model is adopted in this paper, the prior information is directly the weight matrix/cofactor matrix of the observation vector, and there is no weight reduction problem of the standard EIV for the fixed elements of the coefficient matrix.The corresponding cofactor factor function can be obtained by inverting the IGGIII weight factor function as follows.(16)Rii=1.0e~i≤k0e~ik0k1-k0k1-e~i2k0<e~i≤k11030k1<e~iIt should be noted that the theoretical value of the corresponding cofactor should be infinite because the weight factor of the third line in (16) is zero. However, to facilitate calculation, a large number is used to satisfy the requirement in the numerical calculation, and the abnormal observation information has little effect on parameter estimation.Because of the accuracy difference between variables, the cofactors are calculated with standardized residuale~i to guarantee the antierror effect. The formula of the standardized residual e~i is as follows:(17)e~i=e^iσ0Qei,where e^i is the ith element of residual vector e^, Qei is the cofactor of e^i, and σ0 is the standard error of unit weight, which is expressed as follows:(18)σ0=1.4826×medi=1neiQei.According to the law of cofactor propagation, we can easily obtain cofactor matrices Qe of residual errors; see Wang et al. [6] for the detailed derivation process. Therefore, cofactor matrices Qe can be estimated as(19)Qe=MQRMT,where QR=QG-AATQG-1A-1AT and M=QGTQG-1.It should be pointed out that if the number of redundancies is low, outliers would significantly affect the residuals of normal observations and the observation with larger residuals may not contain outliers. In other words, almost all methods that process outliers will break down when there are only a few redundant observations.


## 4. Robust Prediction Algorithms for Multiframe Transformation

Without loss of generality, in this paper, we consider the transformation between two source frames and one target frame. The two source frames are called the first and second source frames. Assume the following four observation datasets: (i) the common stations of three framesΠ1 = [1,2,…,m1]; (ii) the common stations of the first and second source frames Π2 = [1,2,…,m2]; (iii) the noncommon stations of the first source frame Π3 = [1,2,…,m3]; and (iv) the noncommon stations of the second source frame Π4 = [1,2,…,m4]. The terms X, Y, Z, and W are used to represent the 3D coordinate vectors of the common stations in the above four observation datasets, respectively. Considering the first dataset as an example, X, X′, and X′′ denote the 3D coordinate vectors of the target frame, first source frame, and second source frame, respectively. The observation or estimation equations with respect to the ith element in each dataset can be expressed as(20)Xi-eXi=t1+s1R1Xi′-eXi′Xi-eXi=t2+s2R2Xi′′-eXi′′(21)Y-i=t1+s1R1Yi′-eYi′Y-i=t2+s2R2Yi′′-eYi′′Z-i=t1+s1R1Zi′-eZi′W-i=t2+s2R2Wi′′-eWi′′where e denotes the error vector that corresponds to the coordinate observation vector; t1 and t2 are the 3 × 1 translation parameters between the first source frame or second source frame and the target frame, respectively; s1 and s2 are the scale parameters between the first source frame or second source frame and the target frame, respectively; and R1 and R2 are the 3 × 3 rotation matrices between the first source frame or second source frame and the target frame, respectively, and can be described as follows:(22)R1=Rγ1Rβ1Rα1,R2=Rγ2Rβ2Rα2(23)Rα=1000cos⁡αsin⁡α0-sin⁡αcos⁡α;Rβ=cos⁡β0-sin⁡β010sin⁡β0cos⁡β;Rγ=cos⁡γsin⁡γ0-sin⁡γcos⁡γ0001,where α1, β1, and γ1 are the rotation parameters between the first source frame and the target frame; α2, β2, and γ2 are the rotation parameters between the second source frame and the target frame.The corresponding random model is as follows:(24)e~N0,σ02Q(25)e=eXTeX′TeX′′TeY′TeY′′TeZ′TeW′′TT(26)Q=QXX0000000QX′X′QX′Y′QX′Z′0000QY′X′QY′Y′QY′Z′0000QZ′X′QZ′Y′QZ′Z′0000000QX′′X′′QX′′Y′′QX′′W′′0000QY′′X′′QY′′Y′′QY′′W′′0000QW′′X′′QW′′Y′′QW′′W′′.In a similar manner, (20) and (21) are expanded using the Taylor series expansion method. Then,(27)L1i-eXi+s10R10eXi′=A1,jdξ1L2i-eXi+s20R20eXi′′=A2,jdξ2L3i+s10R10eYi′=A3,jdξ1-dY-iL4i+s20R20eYi′′=A4,jdξ2-dY-iL5i+s10R10eZi′=A5,jdξ1-dZ-iL6i+s20R20eWi′′=A6,jdξ2-dW-iwhere(28)L1i=Xi-t10-s10R10Xi′,L2i=Xi-t20-s20R20Xi′′,L3i=Yi0-t10-s10R10Yi′,L4i=Yi0-t20-s20R20Yi′′,L5i=Zi0-t10-s10R10Zi′,L6i=Wi0-t20-s20R20Wi′′A1,j=I3R10Xi′-eXi′0s10δR1δα1Xi′-eXi′0s10δR1δβ1Xi′-eXi′0s10δR1δγ1Xi′-eXi′0,A2,j=I3R20Xi′′-eXi′′0s20δR2δα2Xi′′-eXi′′0s20δR2δβ2Xi′′-eXi′′0s20δR2δγ2Xi′′-eXi′′0,δRδα=RγRβδRαδα,δRδα=RγδRβδβRα,δRδα=δRγδγRβRα,δRαδα=0000-sin⁡αcos⁡α0-cos⁡α-sin⁡α,δRβδβ=-sin⁡β0-cos⁡β000cos⁡β0-sin⁡β,δRγδγ=-sin⁡γcos⁡γ0-cos⁡γ-sin⁡γ0000,where the format of A3,j and A5,j is the same as that of A1,j; the format of A4,j and A6,j is the same as that of A2,j.All observations can be linked to obtain(29)L1-eX+B1eXi′=A1dξ1L2-eX+B2eXi′′=A2dξ2L3+B3eY′=A3dξ1-dY-L4+B4eY′′=A4dξ2-dY-L5+B5eZ′=A5dξ1-dZ-L6+B6eW′′=A6dξ2-dW-where(30)L1=L1,1TL1,2T⋯L1,m1TT,L2=L2,1TL2,2T⋯L2,m1TT,L3=L3,1TL3,2T⋯L3,m2TT,L4=L4,1TL4,2T⋯L4,m2TT,L5=L5,1TL5,2T⋯L5,m3TT,L6=L6,1TL6,2T⋯L6,m4TT,B1=s10R10⊗Im1,B2=s20R20⊗Im1,B3=s10R10⊗Im2,B4=s20R20⊗Im2,B5=s10R10⊗Im3,B6=s20R20⊗Im4.Then the standard LS form of (29) is described as follows: (31)l-Ge=Adθ,where(32)A=A100000A2000A30-Im2000A4-Im200A400-Im300A600-Im4,G=Im1-B100000Im10-B20000000-B30000000-B40000000-B50000000-B6,l=L1TL2TL3TL4TL5TL6TT,dθ=dξ1Tdξ2TdY-dZ-dW-.Following (26) and (31), the robust LS solution of (31) can be obtained using the robust prediction estimation method in the previous section.


## 5. Experimental Analysis

Assume that (i) the transformation parameters between the first source frame and the target frame aret1 = [100-99200]T,s1 = 0.83, α1= 0.45, β1= 0.37, γ1= 0.54; (ii) the transformation parameters between the second source frame and the target frame are t2 = [200120-170]T, s2 = 0.79, α2= 0.23, β2= 0.48, γ2= 0.33. The four dimensions of the corresponding observation datasets are m1= 10, m2 = 6, m3 = 5, and m4 =4, respectively.Considering the WGS-84 frame as the target frame, the geodetic coordinates of the target frame, including the latitude, longitude, and geodetic height, which are randomly generated in the form of a uniform distribution within the ranges-π/2,π/2, [-π,π], and [10,3000], are transformed into the corresponding 3D rectangular coordinates. The observation station coordinates of the first and second source frames are obtained from the true value of the given transformation parameters. The 3D distribution of the simulation data is shown in Figure 1.Figure 1
3D distribution of observation stations.The detailed procedures for generating covariance matrices of 10 observation stations in the target frame are as follows:( 1 ) Generate 30 diagonal matrices D, that is, the standard error, using uniform distribution U [0.03m,0.04m].( 2 ) Generate 30 × 30 random matrix M using uniform distribution U [-0.80.8], M=T1T2 using QR decomposition forM.( 3 ) Generate 10 covariance matrices of 10 observation stations in the target frame: Q1=T1DT1T.The variance-covariance matrix generation process for the first source frame and second source frame is as above, with the uniform distributions being U[0.05m,0.7m] and U[0.08m,0.1m], respectively. The correlation coefficients are generated using uniform distribution U [-0.8,0.8].The truth values of the parameters and coordinates above are constant when the Monte Carlo experiment is run 1,000 times, but the measurement error of each experiment is generated independently using the variance-covariance matrices mentioned above with zero mean Gaussian distribution. Additionally, assume that the number of gross errors is three, the location of gross errors are randomly generated at the common stations, and the sizes are in the range (−30, −10) and (10, 30) times the prior standard error. Observation data with random errors are added each time.Set the prior unit weight mean error to 0.01, and the calculation scheme is as follows:Before the gross errors are added, three schemes are proposed:(1) the prediction estimation algorithm of the independent transformation between the first source frame and the target frame [19]; (2) the prediction estimation algorithm between the second source frame and the target frame; and (3) the prediction estimation multiframe total transformation algorithm proposed in this paper.After the gross errors are added, four additional schemes are proposed:(4) the prediction estimation algorithm of the total transformation of the multiple frames proposed in this paper; (5) the robust prediction estimation algorithm between the first source frame and the target frame; (6) the robust prediction estimation algorithm between the second source frame and the target frame; and (7) the robust prediction estimation multiframe total transformation algorithm proposed in this paper.According to the results of different schemes, the seven transformation parameters and the prediction points in the target frame are calculated with Root Mean Square (RMS) as the accuracy evaluation index, and the specific calculation formula is described as follows:(33)σ2=∑inδ^-δn,where δ^ is the estimated value, δ is the true value, and n is the number of estimations.From Tables1 and 2 and Figures 1, 2, and 3, we can see the following:Table 1
Transformation parameter accuracies in terms of RMSE between the first source frame and the target frame.
Scheme Δ X / m Δ Y / m Δ Z / m s / p p m ε x / m a s ε y / m a s ε z / m a s ( 1 ) 0.0262 0.0254 0.0242 0.0045 360.1 968.5 516.8 ( 3 ) 0.0244 0.0226 0.0217 0.0033 113.5 101.6 182.6 ( 4 ) 0.1343 0.1252 0.1038 0.1048 1265.3 1742.5 1683.7 ( 5 ) 0.0353 0.0398 0.0344 0.0078 432.4 1023.8 653.3 ( 7 ) 0.0268 0.0287 0.0259 0.0057 328.6 562.4 423.5Table 2
Transformation parameter accuracies in terms of the RMSE between the second source frame and the target frame.
scheme Δ X/m Δ Y/m Δ Z/m s/ppm ε x/mas ε y/mas ε z/mas ( 2 ) 0.0455 0.0370 0.0381 0.0096 370.1 814.9 988.7 ( 3 ) 0.0365 0.0324 0.0323 0.0073 344.4 396.1 391.6 ( 4 ) 0.1562 0.1437 0.1517 0.1359 1472.2 1553.6 1842.4 ( 6 ) 0.0524 0.0418 0.0437 0.0103 417.6 931.2 1024.3 ( 7 ) 0.0393 0.0361 0.0376 0.0087 382.3 448.7 512.8Figure 2
RMSE of the coordinate estimation accuracies in the target frame for the stations in setΠ2.Figure 3
RMSE of the coordinate estimation accuracies in the target frame for the stations in setΠ3.( 1 ) Before the gross errors are added, comparing the results from schemes 1, 2, and 3, the result of scheme 3 is optimal, and it lies not only in the transformation parameter estimation but also particularly in common stations of the two source frames (the observation set Π2). There are two transformation solutions for the same reference station in the target frame using schemes 1 and 2, as shown in Figure 2, so it is difficult to decide which one is the optimal solution. The multiframe unified transformation model avoids the problem of solution selection; in addition, unlike an independent transformation, it also takes into full account the existing information of two source frames; thus, the accuracy of the solution is improved.( 2 ) Figures 3 and 4 show the transformation results of the noncommon stations of a single source frame. Better transformation results are obtained because the multiframe transformation can improve the transformation accuracy.Figure 4
RMSE of the coordinate estimation accuracies in the target frame for the stations in setΠ4.( 3 ) After the gross errors are added, a distortion of the estimated parameters often appears when the robust algorithm is not adopted, as in scheme 4. Compared with scheme 4, the robust estimation adopted in schemes 5, 6, and 7 results in different degrees of improvement. The total transformation model of multiple frames takes more existing information into account; that is, it has better spatial geometric distribution. Therefore, under the premise of adopting the same robust method, the robust prediction estimation algorithm (scheme 7), which adopts the total transformation of multiple frames, has the best performance and can obtain more accurate solutions than independent transformation.( 4 ) There are multiple advantages in taking the coordinates of the points to be converted in the target frame as parameters: it not only extends the coverage of the new frame site but also obtains a unified and complete network frame, that is, one that includes coordinate valuation and precision information.


## 6. Conclusions

In this paper, Helmert transformations among three frames, that is, one target frame and two source frames, were studied as an example. A robust prediction algorithm based on the general EIV prediction algorithm and robust estimation was derived in detail and was applied to achieve multiframe total transformation. Experiments demonstrated that the algorithm proposed in this paper not only effectively avoids the problem of multiple solutions for reference station coordinates by using existing observations to obtain an optimal solution for the total transformation of multiple frames, but also effectively resists interference caused by gross errors. Additionally, we assumed that gross errors can appear in a common station; if a gross error appears in a common station, it can only be eliminated and not solved. Then it will be impossible to solve the coordinates of noncommon stations with gross errors in the target frame; this latter problem will need to be considered in a further in-depth study.



---

*Source: 5173956-2019-02-11.xml*