- Email: [email protected]
Helmert transformation with mixed geodetic and Cartesian coordinates
Available online at www.sciencedirect.com
ScienceDirect Advances in Space Research xxx (2017) xxx–xxx www.elsevier.com/locate/asr
Helmert transformation with mixed geodetic and Cartesian coordinates Peng Lin a,b, Guobin Chang a,b, Jingxiang Gao a,b,⇑, Qianxin Wang a,b, Hefang Bian a,b a
NASG Key Laboratory of Land Environment and Disaster Monitoring, China University of Mining & Technology, Xuzhou 221116, Jiangsu, China b School of Environment Science and Spatial Informatics, China University of Mining & Technology, Xuzhou 221116, Jiangsu, China Received 27 March 2017; received in revised form 18 September 2017; accepted 24 November 2017
Abstract It is common in geodetic practice that station geodetic coordinates in an old frame are measured/estimated by adjusting conventional leveling and vertical networks (together with a vertical datum), whereas their Cartesian coordinates in a new frame are measured/estimated by adjusting modern 3D Global Navigation Satellite System (GNSS) networks. To estimate Helmert transformation parameters between the two frames, one can simply convert the geodetic coordinates to their Cartesian counterparts in the same frame, and then perform a least-squares transformation with Cartesian coordinates in both frames. This stepwise approach is not optimal because of nonlinearity in the first step. In this work, a direct transformation is conducted. A functional model with mixed geodetic and Cartesian measurements is followed, and a realistic stochastic model considering measurement errors in both frames is adopted. Correlations between common and non-common stations are also taken into account. Weighted least-squares estimates of the transformation parameters and transformed/adjusted coordinates in the new frame are derived in detail, the latter representing the enlarged network of the new frame. Simulations are conducted, and the results validate the superiority of the proposed method compared with the stepwise method. Ó 2017 COSPAR. Published by Elsevier Ltd. All rights reserved.
Keywords: Helmert transformation; Geodetic coordinates; Cartesian coordinates; Weighted least-squares
1. Introduction It is one of the fundamental tasks of the geodetic community to progressively improve the quality of a terrestrial reference frame in terms of coverage, density, and accuracy (Lu et al., 2014). With development of the Global Navigation Satellite System (GNSS), it has become an indispensable technology in the field of geodesy. Station coordinates surveyed by GNSS are defined in the corresponding coordinate reference frame, such as WGS-84, CGCS2000, PE-90 and ITRF08. Whereas new stations are surveyed using new space techniques, old stations surveyed using conventional techniques should not be ⇑ Corresponding author at: School of Environment Science and Spatial Informatics, China University of Mining & Technology, Xuzhou 221116, Jiangsu, China. E-mail address: [email protected] (J. Gao).
discarded. To incorporate the old stations surveyed in the old/initial frame into a new/target frame, a coordinate frame transformation should be done. Among others, the similarity model, also called the Helmert model, is widely used (Lehmann, 2014). Using coordinates in both frames for common stations, the transformation parameters are estimated. Then some non-common stations whose coordinates are available only in the old frames can be incorporated into the new frame via coordinate transformation, using the estimated parameters to generate an enlarged network that includes more stations than the original one. Measurements in the new frame often have higher accuracies than those in the old one. However, measurement errors in both frames should fully considered in order to be statistically correct. The resulting error model may be called the Gauss–Helmert or error-in-variable model (Neitzel, 2010; Xu et al., 2012). This model is precisely a special case of the general mixed model (Leick et al.,
https://doi.org/10.1016/j.asr.2017.11.029 0273-1177/Ó 2017 COSPAR. Published by Elsevier Ltd. All rights reserved.
Please cite this article in press as: Lin, P., et al. Helmert transformation with mixed geodetic and Cartesian coordinates. Adv. Space Res. (2017), https://doi.org/10.1016/j.asr.2017.11.029
2
P. Lin et al. / Advances in Space Research xxx (2017) xxx–xxx
2015). The least-squares solution for this model may be called total least-squares (Schaffrin and Felus, 2008; Fang, 2013). Based on realistic considerations, the following is assumed: all available horizontal geodetic coordinates in the old frame, i.e., geodetic latitudes and longitudes, are obtained from a single leveling network adjustment; similarly, all available geodetic heights are from a single vertical network adjustment; all available Cartesian coordinates are from a single 3D network. Thus, in the error model, for coordinates from the same network, full correlations should be permitted among a pair of any two variables. The case of more than one network can be safely treated as a special case of the above assumption; we simply let the corresponding covariances be zero. The method developed in this paper can be easily extended to even more general error models, e.g., those with crosscorrelations among different networks. However, we do not believe that this extension has any practical significance. The final aim of the Helmert transformation problem is not merely transformation parameter estimation but transformation of the non-common station coordinates surveyed in the old frame to those in the new frame. It is very likely that the common and non-common stations are in the same network. Then, owing to non-zero correlations between the two types of stations (in the same network), the non-common stations should also be considered. This fact has been recently revealed by several researchers. Li et al. (2012, 2013) proposed the 3D seamless similarity and affine transformation model respectively. Kotsakis et al. (2014) indicated that correlation need to be considered not only between common and noncommon stations in the old frame but between the employed common station and the other reference station in new frame. Wang et al. (2017) proposed a general solution of 3D Helmert transformation, but neglected to consider correlation between common and non-common stations in the new frame. Consequently, besides parameter estimation and transformation of non-common stations surveyed in the old frame, the adjustment of coordinates of the common stations and non-common stations in the new frame is necessary. In special cases, Cartesian coordinates in a new frame can be directly determined using the modern space techniques (Seeber, 2003). However, geodetic coordinates may be first obtained in the old frame using conventional leveling and vertical network surveys. Generally, one can easily convert the geodetic coordinates to their Cartesian counterparts (Grafarend and Okeke, 1998; Vanı´cˇek and Steeves, 1996). Of course, together with the coordinate conversion itself, covariance propagation should also be carried out. Then, with the Cartesian coordinates converted in the old frame and in the new frame, the Helmert transformation can be performed (Fang, 2015). This two-step approach cannot be optimal in the leastsquares sense, though each step is optimal on its own. The loss of optimality comes from information loss in the
first step. The mean and covariance propagation from the nonlinear geodetic-to-Cartesian mapping is often firstorder approximate. Thus, the Cartesian coordinates with the covariance matrix fail to be equivalent to the original geodetic coordinates. In statistical terms, the Cartesian coordinates with an assumed Gaussian distribution fail to have sufficient statistics (Kay, 2013). Even with the mean and covariance corrected by considering the nonlinearity in the first step, (see, e.g., Xue et al., 2016a,b), the nonequivalence statement above is still valid. This is because even with a more accurate mean and covariance, the nonlinearity-introduced non-Gaussian properties are still not taken into account. Given the above analysis, in the special case, we must consider correlation between the common and noncommon stations in the same frame and the mixed geodetic and Cartesian coordinates are directly used to perform a globally optimal Helmert transformation. We define three index sets, i.e., U0 = [1, 2, . . ., m0], U1 = [1, 2, . . ., m1], and U2 = [1, 2, . . ., m2] to represent the common stations, non-common stations surveyed in the old frame, and non-common stations surveyed in the new frame, respectively. Of course, in common geodetic practice, along with the parameter/coordinate estimation, the variance-covariance matrix should be calculated as a measure of estimation accuracy. The weighted least-squares criterion is used in the entire estimation. The estimated or adjusted coordinates in the new frame, together with their variances and covariances, represent the enlarged or updated frame network. This task, as the main aim of the current study, is illustrated in Fig. 1. The remainder of the paper is organized as follows. In Section 2, the main theoretical development is presented. Section 3 describes the simulation to compare the
Input
Output
Fig. 1. Sketch of coordinate frame transformation. Blue circles represent stations whose coordinates are surveyed in new frame, blue stars those surveyed in old frame, and red circles stations whose coordinates are estimated or adjusted in new frame. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Please cite this article in press as: Lin, P., et al. Helmert transformation with mixed geodetic and Cartesian coordinates. Adv. Space Res. (2017), https://doi.org/10.1016/j.asr.2017.11.029
P. Lin et al. / Advances in Space Research xxx (2017) xxx–xxx
developed and stepwise methods. Concluding remarks are made in Section 4. Some derivations and supplementary materials are presented in the appendixes. 2. Mixed geodetic and Cartesian coordinate direct approach For common stations, the 3D Helmert model of a station can be expressed as " ~j ex~j ¼ t þ sRf x
#
!
~ 1;j eu~ 1;j u ;~ h1;j e~h1;j ; ~1;j e~ k
ð1Þ
k1;j
~j and where subscript ‘‘j” denotes the jth common station; x ex~j are Cartesian coordinates in the new frame and its mea~ 1;j , ~ surement errors, respectively; u k1;j and ~ h1;j represent longitude, latitude and geodetic height in the old frame, and eu~ 1;j , e~k1;j and e~h1;j denote their corresponding measurement errors, respectively; t, s and R stand for the translation, scale and rotation of Helmert transformation parameters, respectively. The expression for mapping f that converts the geodetic coordinates to their Cartesian counterparts is given in Appendix A. For non-common stations in the old frame, the function equation is " yj ¼ t þ sRf
# ! ~ 2;j eu~ 2;j u ;~ h2;j e~h2;j ; ~ k2;j e~
ð2Þ
k2;j
where yj denotes the Cartesian coordinate in the new frame that must be estimated; other variables are similar to the above definition. For non-common stations in the new frame, the Cartesian coordinates, which are the same as those of common stations in the new frame, are also of interest and should be treated as parameters to be estimated (Kotsakis et al., 2014). The equation can be expressed as ~j ex~j ¼ xj ; x
ð3Þ
~zj e~zj ¼ zj ;
ð4Þ
where x and z denote the Cartesian coordinates in the new frame that must be estimated; other variables are similar to the above definition. Eq. (3) is modified as follows to estimate these coordinates directly. However, this dimension increase does not change the degree of freedom, which is an apparent result of the equivalence between the original and modified. The corresponding stochastic model should obey the Gaussian distribution with zero mean: e N ð0; Qee Þ; where
ð5Þ
3 T
e ¼ ½ eT~u eT~h1 e~Tv eT~h2 eTx~ e~Tz ; 2 0 Qu~~v 0 0 Qu~~u 6 0 Q~ ~ 0 Q~h1 ~h2 0 6 h1 h1 6 6 Qv~~u 0 Qv~~v 0 0 Qee ¼ 6 6 0 Q~ ~ 0 Q~h2 ~h2 0 6 h2 h1 6 4 0 0 0 0 Qx~ ~x 0 e~uj ¼ ½ eu~ 1;j
0
0
0
T
e~k1;j ; e~vj ¼ ½ eu~ 2;j
Qz~~x
0
3
7 7 7 7 7; 0 7 7 7 Qx~ ~z 5 0 0
Qz~~z
T
e~k2;j :
Qee is a positive definition cofactor matrix of total measurement errors e. Preferably, to accord with the realistic signification, it is assumed that cross-correlation among different frames should be a zero matrix. Eqs. (1) and (2) are clearly nonlinear models. To derive the iterative solution, the Gauss-Newton approach of nonlinear least-squares is used. One may use the following process to solve the problem: (1) Construct the loss function as the weighted sum of the estimated errors (i.e., residuals); (2) construct the Lagrangian to incorporate Eqs. (1)–(4) into the loss function; (3) construct the equation systems to be solved according to the first-order necessary condition of minimizing the Lagrangian, i.e., letting the first-order derivative of the Lagrangian with respect to all variables be estimated as zero (these variables should include the transformation parameters, coordinates in the new frame, estimated errors, and the Lagrange multipliers); (4) solve the above nonlinear equation system using numerical algorithms (Bjorck, 1996); (5) perform error or covariance propagation by linearizing the equations obtained in step (3) with the estimated values. However, an alternative approach, namely Pope iteration, is followed. This alternately executes the following two until convergence: linearizing the nonlinear equations and standard linear estimation (Pope, 1972). The Pope approach, essentially a Gauss-Newton iterative procedure, has a merit as a by-product, i.e., the covariance matrix of the estimate can be readily retrieved during the solution. The linearization is done with known approximate values of the variables involved. Letting subscript or superscript ‘‘0” denote approximate values, we define the following increments as ‘‘d” values: t ¼ t0 dt; s ¼ s0 ds;
ð6Þ 1
R ¼ ðI 3 ½dhÞðI 3 þ ½dhÞ R0 ðI 3 2½dhÞR0 ;
ð7Þ
with the cross-product matrix of the 3 1 vector dh defined as 2 3 0 dh3 dh2 6 7 ½dh ¼ 4 dh3 0 dh1 5: ð8Þ dh2 dh1 0 The rotation matrix is the approximate expression used in the linearization. However, in correcting the approximate values using the estimated increments, the exact expression should be used instead, because it can always
Please cite this article in press as: Lin, P., et al. Helmert transformation with mixed geodetic and Cartesian coordinates. Adv. Space Res. (2017), https://doi.org/10.1016/j.asr.2017.11.029
4
P. Lin et al. / Advances in Space Research xxx (2017) xxx–xxx
ensure orthogonality of R given an orthogonal R0. In this exact expression, dh is just the Gibbs vector or Rodrigues parameters of the delta/increment rotation; here it is not the additive error of the Gibbs vector of the full rotation. For more detail, see e.g., Wang et al. (2016). Besides the above variables whose increments are to be estimated, approximate values of the following variables are needed in the linearization: e0~u ,e0~h ,e~0v and e0~h . 1 2 The approximate values used here should be the most accurate available before the subsequent estimation, in order to perform the best possible linearization. For unknowns, i.e., the transformation parameters and coordinates in the new frame, they should be estimates obtained from the immediately previous iteration. For the geodetic coordinates whose corresponding measurements are available, they should be the adjusted ones, i.e., the original measurements minus measurement error estimates obtained in the immediately previous iteration. Either of two candidates can be used to initialize the iteration. The first, based on practical considerations, one can simply let t0 = 0, s0 = 1, and R0 = I3. In the second, one can apply analytical methods from the literature to the common stations to attain a set of estimates of the parameters (e.g., Chang, 2015, 2016). These methods are used only to obtain parameter estimates, and one does not bother with a corresponding error analysis (Chang, 2016; Chang et al., 2017a,b). The first method is simple but the second
Eqs. (1) and (2) are used to perform the following linearization: L1;j ex~j ¼ B1j n s0 R0 G 1j el~j s0 R0 F 1j e~hj L2;j ¼ B 2j n yj s0 R0 G 2j e~vj s0 R0 F 2j e~hj 8 ~ t0 s0 R0 f ð~uj e0~uj ; ~h1;j e0~h Þ L1;j ¼ x > > 1;j > > > j 0 j 0 > < s0 R0 G 1 e~uj s0 R0 F 1 e~h 1;j ; with 0 ~ 0 > > > L2;j ¼ t0 s0 R0 f ð~vj e~vj ; h2;j e~h2;j Þ > > > : s R G j e0 s R F j e0 0
0
2 ~vj
0
0
ð9Þ ð10Þ
2 ~ h2;j
T T where u~j ¼ ½ u ~ 1;j ~k1;j ; ~vj ¼ ½ u ~ 2;j ~k2;j ; and 8 @f ð~ uj e0~u ;~ h1;j e0~ Þ > h1;j j j > > G1 ¼ > @eT > ~j l > e~uj ¼e0~u ;e~h ¼e0~ > h1;j > 1;j j > > > 0 ;~ 0 Þ > @f ð~ v e h e j > j ~vj 2;j ~ h2;j > > G2 ¼ > @e~T > vj > < e~vj ¼e~0v ;e~h ¼e0~ h2;j 2;j j ; 0 ;h 0 Þ ~ @f ð~ u e e > j ~ ~ uj 1;j > h1;j j > > F1 ¼ > @e~h > j > e~uj ¼e0~u ;e~h ¼e0~ > h1;j 1;j j > > > > 0 ;~ 0 Þ > @f ð~ v e h e j > ~vj 2;j ~ h2;j > > F 2j ¼ > @e~h > : j 0 0 e~vj ¼e~v ;e~h j
2;j
¼e~
h2;j
whose detailed expression is presented in Appendix A,
should be more accurate and can be used in cases of nonsmall-angle-rotation. In addition, the approximate values of measurement errors can be given the zero matrix of T ð6 m0 þ 3 m1 þ 3 m2 Þ 1. Letting n ¼ ½ dt ds dhT ,
For all stations, we can define the following variables:
Please cite this article in press as: Lin, P., et al. Helmert transformation with mixed geodetic and Cartesian coordinates. Adv. Space Res. (2017), https://doi.org/10.1016/j.asr.2017.11.029
P. Lin et al. / Advances in Space Research xxx (2017) xxx–xxx
2
2 6 6 6 F1 ¼ 6 6 4
3
s0 R0 G 11
6 6 6 G1 ¼ 6 6 4
s0 R0 G 21 ..
. s0 R0 G m1 0
6 6 6 G2 ¼ 6 6 4
s0 R0 F 21 ..
. s0 R0 F m1 0
6 6 6 F2 ¼ 6 6 4
1 ^ ^e ¼ Qee C T ðCQee C T Þ ðL AbÞ
3
s0 R0 G 22 ..
. s0 R0 G m2 1
s0 R0 F 12 s0 R0 F 22 ..
. s0 R0 F m2 1
7 7 7 7; 7 5 3 7 7 7 7: 7 5
Then we have 8 L1 þ G 1 e~u þ F 1 e~h1 ¼ B1 n > > > < L þ G e þ F e~ ¼ B n y 2 2 ~v 2 h2 1 > ~ z e ¼ z ~z > > : ~ ex~ ¼ x x
3. Monte Carlo simulation
ð11Þ
# ¼ xT
yT
z
2
T T
; b ¼ nT
#
T T
; L ¼ LT1
LT2
~T x
~zT
T
;
3
0 0 B1 0 6B 0 I 0 7 m2 6 2 7 A¼6 7; 4 0 I m1 0 0 5 0 0 0 I m3 2 0 0 G 1 F 1 6 0 0 G 2 F 2 6 C ¼6 4 0 0 0 0 0 0 0 0
I m1 0 I m1 0
0
3
0 7 7 7: 0 5 I m3
Finally, we have the overall measurement equation: L C e ¼ Ab
ð12Þ
According to the least-squares theory, the Lagrange loss function is expressed as: Uðe; b; KÞ ¼ eT Pe þ 2K T ðL C e AbÞ
The true transformation parameters from the new frame to the old were chosen as follows: ttrue = [100 5 200]T, strue = 0.89, and htrue = [0.035 0.029 0.033]T. The true value of the rotation matrix is Rtrue ¼ ðI 3 ½htrue ÞðI 3 þ ½htrue Þ
We define the following variables:
ð16Þ
The top part of (14) and the first four parts of (16) are 0 the A, C and L1 matrices in the iterative procedure. The bottom part of (14) and bottom-right part of (15) constitute the updated network of the new frame. As stated previously, the linearization-estimation pair should be performed iteratively; and in the intermediate iterations, only the parameter estimation (14) is necessary, while the variance-covariance matrix estimation (15) can be performed and stored in the final iteration. ^ represents coordinates of the new Because parameter # frame, whereas ^n symbolizes increments of Helmert transformation parameters, the object calculating the conver^ which ^ but ^n associated with D#, gence condition is not b ^ in the adjacent iteration. In addidenote the difference of # tion, in the iterative procedure, the scale parameter should be non-positive. If it is less than zero, it should be assigned to1again to avoid iterative divergence.
7 7 7 7; 7 5
s0 R0 G 12
2
ð15Þ
7 7 7 7; 7 5
3
s0 R0 F 11
2
5
ð13Þ
where K denotes the ð6 m0 þ 3 m1 þ 3 m2 Þ 1 vector of ‘‘Lagrange multipliers”. Then, the solution is readily expressed as: ð14Þ
1
Assuming m0 = 15, m1 = 20 and m2 = 20. All 45 stations are selected randomly, each of whose true geodetic coordinates in the old are generated as follows: (1) Generate longitude ktrue from the uniform distribution U[p p); (2) and latitude utrue from U[0.45p 0.45p]; 3) generate height htrue from U[10 1000]. Then, we calculate the true Cartesian coordinates in the new frame of the 15 common stations and 20 non-common stations surveyed only in the new frame: utrue xtrue ¼ ttrue þ strue Rtrue f ; htrue ktrue The distribution of all 45 stations in the new frame is illustrated in Fig. 2. The 105 105 matrix for all Cartesian coordinate measurements in the new frame, namely Q1, is generated as follows: (1) Designate the diagonal matrix with its diagonal elements the same as those of Q1 as D1, and generate the square-root of each of its 105 diagonal elements from U [0.05 0.15]; (2) generate an arbitrary 105 105 matrix X1 and perform the QR decomoposition as X1 = Q1R1; 3) let Q1 = Q1D1QT1 . The 60 60 matrix for all horizontal geodetic coordinates, namely Q2, and the 30 30 matrix for all vertical coordinate measurements, namely Q3, are generated in the same way; the only differences are such that
Please cite this article in press as: Lin, P., et al. Helmert transformation with mixed geodetic and Cartesian coordinates. Adv. Space Res. (2017), https://doi.org/10.1016/j.asr.2017.11.029
6
P. Lin et al. / Advances in Space Research xxx (2017) xxx–xxx
Fig. 2. 3D distribution of stations involved in simulation. Coordinates are Cartesian in the new frame. Red, green, and blue spheres denote common stations, non-common stations surveyed in the new frame, and non-common stations surveyed in the old frame, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
0.10
0.40 0.35
12
RMS (ppm)
0.25 0.20 0.15 0.10
RMS (mas)
0.08
0.30
RMS (m)
14
0.06 0.04
8 6 4
0.02
2
0.05 0.00
10
0.00
0
translations
scale
rotation angles
Fig. 3. Transformation parameter estimation accuracies in terms of RMSE. Blue and red bars represent stepwise and direct approaches, respectively. Unit length along the vertical axis represents 1 m, 1 ppm (part per million) and 1 mas (milli-arc-second), for the translation, scale and rotation parameters, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
0.10 0.08 0.06
Z component
Y component
X component
0.12
0.04
0.16
0.16
0.14
0.14 0.12 0.10 0.08
4 6
8 10 12 14 16
common station
0.12 0.10 0.08 0.06
0.06
0 2
0.14
0
2
4
6
8 10 12 14 16
common station
0
2
4
6
8 10 12 14 16
common station
Fig. 4. Cartesian coordinate estimation accuracies in terms of RMSE for the 15 common stations. Blue, red and yellow bars represent estimates by the stepwise approach, by the direct approach, and the original measurements, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
for the former and latter, the distributions from which the square-root of their diagonal elements are generated are U [1.5 3]107 radians and U[1 5] meters, respectively. The results are " # Q~h1 ~h1 Q~h1 ~h2 Qx~ Qu~~u Qu~~v ~x Q x~ ~z Q1 ¼ ;Q2 ¼ ;Q3 ¼ ; Q~h2 ~h1 Q~h2 ~h2 Qz~~x Qz~~z Qv~~u Qv~~v Then, the total variance covariance Qee can be constructed using Q1, Q2 and Q3.
Overall, 100 Monte Carlo experiments were conducted. In all experiments, the above true parameters and coordinates were fixed. Measurement errors were generated independently across different experiments from zero-mean Gaussian distributions with the above variance covariance matrices (Q1, Q2 and Q3). Two approaches were evaluated, the direct approach proposed in this paper and stepwise approach. To be selfcontained, the stepwise approach is presented in Appendix B. The following estimation errors for both approaches
Please cite this article in press as: Lin, P., et al. Helmert transformation with mixed geodetic and Cartesian coordinates. Adv. Space Res. (2017), https://doi.org/10.1016/j.asr.2017.11.029
7
3.0 3.0
2.5
2.5
Z component
3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0
Y componet
X componet
P. Lin et al. / Advances in Space Research xxx (2017) xxx–xxx
2.0 1.5 1.0 0.5 0.0
0 1 2 3 4 5 6 7 8 9 10 11
non-common station (old frame)
2.0 1.5 1.0 0.5 0.0
0 1 2 3 4 5 6 7 8 9 10 11
0 1 2 3 4 5 6 7 8 9 10 11
non-common station (old frame)
non-common station (old frame)
0 2 4 6 8 10 12 14 16 18 20
non-common station(new frame)
0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00
Z component
0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00
Y component
X component
Fig. 5. Cartesian coordinates estimation accuracies in terms of RMSE for the 10 non-common stations surveyed only in the old frames. Blue, red and yellow bars represent estimates by stepwise approach, direct approach, and transforming measured geodetic coordinates. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
0 2 4 6 8 10 12 14 16 18 20
non-common station(new frame)
0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00
0 2 4 6 8 10 12 14 16 18 20
non-common station(new frame)
Fig. 6. Cartesian coordinate estimation accuracies in terms of RMSE for the 20 non-common stations surveyed only in the new frame. Blue, red and yellow bars represent estimates by stepwise approach, direct approach, and original measurements, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
were checked: translation parameter estimation errors ^t ttrue ; scale parameter estimation errors ^s strue ; rotation parameter estimation errors (in radians), whose cross
T ^ R ^ Rtrue ^ Rtrue R ^T þ R product is 12 R ; finally, coordinate estimation errors for all 45 stations x xtrue . To weaken random effects, root mean squares of all error terms above were calculated across the 100 Monte Carlo runs. These root mean squared errors (RMSE) are exactly the indices representing accuracies of the corresponding estimates. Bar charts were used to explicitly compare the performance of the two approaches. The results are presented in Figs. 3–6. These correspond to the transformation parameter, coordinates of the 15 common stations, coordinates of the 10 non-common stations surveyed only in the old frame, and coordinates of the 20 non-common stations surveyed only in the new frame, respectively. In all figures, blue bars represent RMSEs of the stepwise approach, and red bars represent those of the direct approach. In Figs. 4 and 6, RMSEs for the originally measured coordinates are depicted by yellow bars. In Fig. 5, RMSEs for coordinates obtained by transforming the measured geodetic coordinates using the estimated parameters are also shown, by yellow bars. In Fig. 3, the unit length along the vertical axis represents 1 m, 1 ppm (part per million) and 1 mas (milli-arc-second), for the translation, scale and rotation parameters, respectively. Horizontal axes in Figs. 4–6 represent indices of the stations.
Fig. 3 shows that although performances of the two approaches were comparable in estimating the rotation parameters, the direct approach performed better than the stepwise approach in estimating the translation and scale parameters. Figs. 4–6 reveal the following. First, for all stations with their coordinates directly measured in the new frame (i.e., the stations in Figs. 4 and 6), the adjusted/estimated coordinates using either of the two approaches are more accurate than the measured coordinates. Second, for the non-common stations surveyed in the old frame, (those in Fig. 5), the adjusted coordinates are more accurate than those transformed using the estimated transformation parameters. These two observations clearly show the necessity of incorporating the coordinates as parameters in addition to the transformation parameters in the Helmert transformation problem. This again validates the statements made in Kotsakis et al. (2014). Third, for most of the coordinates, the direct approach did better than the stepwise one, though the former could not outperform the latter for all stations or coordinates. Of course, there were still some stations or coordinate components for which the stepwise approach outperformed the direct approach. This appears somewhat irregular. Even more irregular is that there were some coordinate components for which the original measurement was more accurate than estimates by both approaches. It is believed that the randomness effect was compensated by the Monte Carlo experiments. Therefore, this ‘‘irregular” phenomenon
Please cite this article in press as: Lin, P., et al. Helmert transformation with mixed geodetic and Cartesian coordinates. Adv. Space Res. (2017), https://doi.org/10.1016/j.asr.2017.11.029
8
P. Lin et al. / Advances in Space Research xxx (2017) xxx–xxx
cannot be attributed to the randomness effect. Further studies are needed to interpret the observed phenomenon and to possibly further optimize the proposed approach. 4. Concluding remarks Given geodetic coordinates in one frame and Cartesian coordinates in the other, a potential approach to performing the Helmert transformation is to first convert geodetic coordinates to their Cartesian counterparts within the same frame and then do a Helmert transformation with Cartesian coordinates in both frames. However, in this work, we clearly show that this stepwise approach is not globally optimal. A direct approach is proposed instead, using the geodetic coordinates directly. By recalling that one of the tasks in the Helmert transformation problem, possibly the most important one, is to estimate/adjust coordinates of all relevant stations in the new frame, and that measured coordinates of different stations may be correlated, the augmented problem suggested by Kotsakis et al. (2014) is solved. This was accomplished by treating the coordinates of all stations in the new frame as parameters, in addition to the transformation parameters. The least-squares criterion was used to solve the problem, and the Pope iteration used to obtain the solution numerically. Monte Carlo simulations were conducted, and the results validate our theoretical statements; overall performance of the proposed direct method is superior to that of the stepwise one. Acknowledgments This work was supported by the following funds: National Natural Science Foundation of China (No. 41774005 and 41674008) and Research Innovation Program for College Graduates of Jiangsu Province (KYZZ16_0217, SJCX17_0522). Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j. asr.2017.11.029. References Bjorck, A., 1996. Numerical Methods for Least Squares Problems. SIAM, Philadelphia. Chang, G., 2015. On least-squares solution to 3D similarity transformation problem under Gauss-Helmert model. J. Geodesy 89, 573–576. Chang, G., 2016. Closed form least-squares solution to 3D symmetric Helmert transformation with rotational invariant covariance structure. Acta Geodaetica et Geophysica 51, 237–244.
Chang, G., Xu, T., Wang, Q., 2017a. Error analysis of the 3D similarity coordinate transformation. GPS Solut. 21 (3), 963–971. Chang, G., Xu, T., Wang, Q., Zhang, S., et al., 2017b. A generalization of the analytical least-squares solution to the 3D symmetric Helmert coordinate transformation problem with an approximate error analysis. Adv. Space Res. 59 (10), 2600–2610. Fang, X., 2013. Weighted total least squares: necessary and sufficient conditions, fixed and random parameters. J. Geodesy 87, 733–749. Fang, X., 2015. Weighted total least-squares with constraints: a universal formula for geodetic symmetrical transformations. J. Geodesy 89, 459– 469. Grafarend, E., Okeke, F., 1998. Transformation of conformai coordinates of type Mercator from a global datum (WGS 84) to a local datum (Regional, national). Mar. Geodesy 21 (3), 169–180. Kay, S.M., 2013. Fundamentals of Statistical Signal Processing, Volume I: Estimation Theory. Pearson Education, New York. Kotsakis, C., Vatalis, A., Sanso, F., 2014. On the importance of intraframe and inter-frame covariances in frame transformation theory. J. Geodesy 88, 1187–1201. Lehmann, R., 2014. Transformation model selection by multiple hypotheses testing. J. Geodesy 88, 1117–1130. Leick, A., Rapoport, L., Tatarnikov, D., 2015. GPS Satellite Survey, fourth ed. Wiley, New Jersey. Li, B., Shen, Y., Li, W., 2012. The seamless model for three-dimensional datum transformation. Sci. China Earth Sci. 55, 2099–2108. Li, B., Shen, Y., Zhang, X., Li, C., Lou, L., 2013. Seamless multivariate affine error-in-variables transformation and its application to map rectification. Int. J. Geogr. Inform. Sci. 27, 1572–1592. Lu, Z., Qu, Y., Qiao, S., 2014. Geodesy: Introduction to Geodetic Datum and Geodetic Systems. Springer, Berlin. Neitzel, F., 2010. Generalization of total least-squares on example of unweighted and weighted 2D similarity transformation. J. Geodesy 84, 751–762. Pope, A.J., 1972. Some pitfalls to be avoided in the iterative adjustment of nonlinear problems. In: Proceedings of the 38th Annual Meeting. American Society of Photogrammetry, Washington DC, pp. 449–477. Schaffrin, B., Felus, Y.A., 2008. On the multivariate total least-squares approach to empirical coordinate transformations. Three algorithms. J. Geodesy 82, 373–383. Seeber, G., 2003. Satellite Geodesy. Walter de Gruyter, Berlin. Vanı´cˇek, P., Steeves, R.R., 1996. Transformation of coordinates between two horizontal geodetic datums. J. Geodesy 70 (11), 740–745. Wang, B., Li, J., Liu, C., Yu, J., 2017. Generalized total least squares prediction algorithm for universal 3D similarity transformation. Adv. Space Res. 59, 815–823. Wang, Q., Chang, G., Xu, T., Zou, Y., 2016. Representation of the rotation parameter estimation errors in the Helmert transformation model. Surv. Rev., 1–13 https://doi.org/10.1080/ 00396265.2016.1234806. Xu, P., Liu, J., Shi, C., 2012. Total least squares adjustment in partial errors-in-variables models: algorithm and statistical analysis. J. Geodesy 86, 661–675. Xue, S., Dang, Y., Liu, J., Mi, J., Dong, C., Cheng, Y., Wang, X., Wan, J., 2016a. Bias estimation and correction for triangle-based surface area calculations. Int. J. Geogr. Inform. Sci. 30, 2155–2170. Xue, S., Yang, Y., Dang, Y., 2016b. Formulas for precisely and efficiently estimating the bias and variance of the length measurements. J. Geogr. Syst. 18, 399–415.
Please cite this article in press as: Lin, P., et al. Helmert transformation with mixed geodetic and Cartesian coordinates. Adv. Space Res. (2017), https://doi.org/10.1016/j.asr.2017.11.029
ScienceDirect Advances in Space Research xxx (2017) xxx–xxx www.elsevier.com/locate/asr
Helmert transformation with mixed geodetic and Cartesian coordinates Peng Lin a,b, Guobin Chang a,b, Jingxiang Gao a,b,⇑, Qianxin Wang a,b, Hefang Bian a,b a
NASG Key Laboratory of Land Environment and Disaster Monitoring, China University of Mining & Technology, Xuzhou 221116, Jiangsu, China b School of Environment Science and Spatial Informatics, China University of Mining & Technology, Xuzhou 221116, Jiangsu, China Received 27 March 2017; received in revised form 18 September 2017; accepted 24 November 2017
Abstract It is common in geodetic practice that station geodetic coordinates in an old frame are measured/estimated by adjusting conventional leveling and vertical networks (together with a vertical datum), whereas their Cartesian coordinates in a new frame are measured/estimated by adjusting modern 3D Global Navigation Satellite System (GNSS) networks. To estimate Helmert transformation parameters between the two frames, one can simply convert the geodetic coordinates to their Cartesian counterparts in the same frame, and then perform a least-squares transformation with Cartesian coordinates in both frames. This stepwise approach is not optimal because of nonlinearity in the first step. In this work, a direct transformation is conducted. A functional model with mixed geodetic and Cartesian measurements is followed, and a realistic stochastic model considering measurement errors in both frames is adopted. Correlations between common and non-common stations are also taken into account. Weighted least-squares estimates of the transformation parameters and transformed/adjusted coordinates in the new frame are derived in detail, the latter representing the enlarged network of the new frame. Simulations are conducted, and the results validate the superiority of the proposed method compared with the stepwise method. Ó 2017 COSPAR. Published by Elsevier Ltd. All rights reserved.
Keywords: Helmert transformation; Geodetic coordinates; Cartesian coordinates; Weighted least-squares
1. Introduction It is one of the fundamental tasks of the geodetic community to progressively improve the quality of a terrestrial reference frame in terms of coverage, density, and accuracy (Lu et al., 2014). With development of the Global Navigation Satellite System (GNSS), it has become an indispensable technology in the field of geodesy. Station coordinates surveyed by GNSS are defined in the corresponding coordinate reference frame, such as WGS-84, CGCS2000, PE-90 and ITRF08. Whereas new stations are surveyed using new space techniques, old stations surveyed using conventional techniques should not be ⇑ Corresponding author at: School of Environment Science and Spatial Informatics, China University of Mining & Technology, Xuzhou 221116, Jiangsu, China. E-mail address: [email protected] (J. Gao).
discarded. To incorporate the old stations surveyed in the old/initial frame into a new/target frame, a coordinate frame transformation should be done. Among others, the similarity model, also called the Helmert model, is widely used (Lehmann, 2014). Using coordinates in both frames for common stations, the transformation parameters are estimated. Then some non-common stations whose coordinates are available only in the old frames can be incorporated into the new frame via coordinate transformation, using the estimated parameters to generate an enlarged network that includes more stations than the original one. Measurements in the new frame often have higher accuracies than those in the old one. However, measurement errors in both frames should fully considered in order to be statistically correct. The resulting error model may be called the Gauss–Helmert or error-in-variable model (Neitzel, 2010; Xu et al., 2012). This model is precisely a special case of the general mixed model (Leick et al.,
https://doi.org/10.1016/j.asr.2017.11.029 0273-1177/Ó 2017 COSPAR. Published by Elsevier Ltd. All rights reserved.
Please cite this article in press as: Lin, P., et al. Helmert transformation with mixed geodetic and Cartesian coordinates. Adv. Space Res. (2017), https://doi.org/10.1016/j.asr.2017.11.029
2
P. Lin et al. / Advances in Space Research xxx (2017) xxx–xxx
2015). The least-squares solution for this model may be called total least-squares (Schaffrin and Felus, 2008; Fang, 2013). Based on realistic considerations, the following is assumed: all available horizontal geodetic coordinates in the old frame, i.e., geodetic latitudes and longitudes, are obtained from a single leveling network adjustment; similarly, all available geodetic heights are from a single vertical network adjustment; all available Cartesian coordinates are from a single 3D network. Thus, in the error model, for coordinates from the same network, full correlations should be permitted among a pair of any two variables. The case of more than one network can be safely treated as a special case of the above assumption; we simply let the corresponding covariances be zero. The method developed in this paper can be easily extended to even more general error models, e.g., those with crosscorrelations among different networks. However, we do not believe that this extension has any practical significance. The final aim of the Helmert transformation problem is not merely transformation parameter estimation but transformation of the non-common station coordinates surveyed in the old frame to those in the new frame. It is very likely that the common and non-common stations are in the same network. Then, owing to non-zero correlations between the two types of stations (in the same network), the non-common stations should also be considered. This fact has been recently revealed by several researchers. Li et al. (2012, 2013) proposed the 3D seamless similarity and affine transformation model respectively. Kotsakis et al. (2014) indicated that correlation need to be considered not only between common and noncommon stations in the old frame but between the employed common station and the other reference station in new frame. Wang et al. (2017) proposed a general solution of 3D Helmert transformation, but neglected to consider correlation between common and non-common stations in the new frame. Consequently, besides parameter estimation and transformation of non-common stations surveyed in the old frame, the adjustment of coordinates of the common stations and non-common stations in the new frame is necessary. In special cases, Cartesian coordinates in a new frame can be directly determined using the modern space techniques (Seeber, 2003). However, geodetic coordinates may be first obtained in the old frame using conventional leveling and vertical network surveys. Generally, one can easily convert the geodetic coordinates to their Cartesian counterparts (Grafarend and Okeke, 1998; Vanı´cˇek and Steeves, 1996). Of course, together with the coordinate conversion itself, covariance propagation should also be carried out. Then, with the Cartesian coordinates converted in the old frame and in the new frame, the Helmert transformation can be performed (Fang, 2015). This two-step approach cannot be optimal in the leastsquares sense, though each step is optimal on its own. The loss of optimality comes from information loss in the
first step. The mean and covariance propagation from the nonlinear geodetic-to-Cartesian mapping is often firstorder approximate. Thus, the Cartesian coordinates with the covariance matrix fail to be equivalent to the original geodetic coordinates. In statistical terms, the Cartesian coordinates with an assumed Gaussian distribution fail to have sufficient statistics (Kay, 2013). Even with the mean and covariance corrected by considering the nonlinearity in the first step, (see, e.g., Xue et al., 2016a,b), the nonequivalence statement above is still valid. This is because even with a more accurate mean and covariance, the nonlinearity-introduced non-Gaussian properties are still not taken into account. Given the above analysis, in the special case, we must consider correlation between the common and noncommon stations in the same frame and the mixed geodetic and Cartesian coordinates are directly used to perform a globally optimal Helmert transformation. We define three index sets, i.e., U0 = [1, 2, . . ., m0], U1 = [1, 2, . . ., m1], and U2 = [1, 2, . . ., m2] to represent the common stations, non-common stations surveyed in the old frame, and non-common stations surveyed in the new frame, respectively. Of course, in common geodetic practice, along with the parameter/coordinate estimation, the variance-covariance matrix should be calculated as a measure of estimation accuracy. The weighted least-squares criterion is used in the entire estimation. The estimated or adjusted coordinates in the new frame, together with their variances and covariances, represent the enlarged or updated frame network. This task, as the main aim of the current study, is illustrated in Fig. 1. The remainder of the paper is organized as follows. In Section 2, the main theoretical development is presented. Section 3 describes the simulation to compare the
Input
Output
Fig. 1. Sketch of coordinate frame transformation. Blue circles represent stations whose coordinates are surveyed in new frame, blue stars those surveyed in old frame, and red circles stations whose coordinates are estimated or adjusted in new frame. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Please cite this article in press as: Lin, P., et al. Helmert transformation with mixed geodetic and Cartesian coordinates. Adv. Space Res. (2017), https://doi.org/10.1016/j.asr.2017.11.029
P. Lin et al. / Advances in Space Research xxx (2017) xxx–xxx
developed and stepwise methods. Concluding remarks are made in Section 4. Some derivations and supplementary materials are presented in the appendixes. 2. Mixed geodetic and Cartesian coordinate direct approach For common stations, the 3D Helmert model of a station can be expressed as " ~j ex~j ¼ t þ sRf x
#
!
~ 1;j eu~ 1;j u ;~ h1;j e~h1;j ; ~1;j e~ k
ð1Þ
k1;j
~j and where subscript ‘‘j” denotes the jth common station; x ex~j are Cartesian coordinates in the new frame and its mea~ 1;j , ~ surement errors, respectively; u k1;j and ~ h1;j represent longitude, latitude and geodetic height in the old frame, and eu~ 1;j , e~k1;j and e~h1;j denote their corresponding measurement errors, respectively; t, s and R stand for the translation, scale and rotation of Helmert transformation parameters, respectively. The expression for mapping f that converts the geodetic coordinates to their Cartesian counterparts is given in Appendix A. For non-common stations in the old frame, the function equation is " yj ¼ t þ sRf
# ! ~ 2;j eu~ 2;j u ;~ h2;j e~h2;j ; ~ k2;j e~
ð2Þ
k2;j
where yj denotes the Cartesian coordinate in the new frame that must be estimated; other variables are similar to the above definition. For non-common stations in the new frame, the Cartesian coordinates, which are the same as those of common stations in the new frame, are also of interest and should be treated as parameters to be estimated (Kotsakis et al., 2014). The equation can be expressed as ~j ex~j ¼ xj ; x
ð3Þ
~zj e~zj ¼ zj ;
ð4Þ
where x and z denote the Cartesian coordinates in the new frame that must be estimated; other variables are similar to the above definition. Eq. (3) is modified as follows to estimate these coordinates directly. However, this dimension increase does not change the degree of freedom, which is an apparent result of the equivalence between the original and modified. The corresponding stochastic model should obey the Gaussian distribution with zero mean: e N ð0; Qee Þ; where
ð5Þ
3 T
e ¼ ½ eT~u eT~h1 e~Tv eT~h2 eTx~ e~Tz ; 2 0 Qu~~v 0 0 Qu~~u 6 0 Q~ ~ 0 Q~h1 ~h2 0 6 h1 h1 6 6 Qv~~u 0 Qv~~v 0 0 Qee ¼ 6 6 0 Q~ ~ 0 Q~h2 ~h2 0 6 h2 h1 6 4 0 0 0 0 Qx~ ~x 0 e~uj ¼ ½ eu~ 1;j
0
0
0
T
e~k1;j ; e~vj ¼ ½ eu~ 2;j
Qz~~x
0
3
7 7 7 7 7; 0 7 7 7 Qx~ ~z 5 0 0
Qz~~z
T
e~k2;j :
Qee is a positive definition cofactor matrix of total measurement errors e. Preferably, to accord with the realistic signification, it is assumed that cross-correlation among different frames should be a zero matrix. Eqs. (1) and (2) are clearly nonlinear models. To derive the iterative solution, the Gauss-Newton approach of nonlinear least-squares is used. One may use the following process to solve the problem: (1) Construct the loss function as the weighted sum of the estimated errors (i.e., residuals); (2) construct the Lagrangian to incorporate Eqs. (1)–(4) into the loss function; (3) construct the equation systems to be solved according to the first-order necessary condition of minimizing the Lagrangian, i.e., letting the first-order derivative of the Lagrangian with respect to all variables be estimated as zero (these variables should include the transformation parameters, coordinates in the new frame, estimated errors, and the Lagrange multipliers); (4) solve the above nonlinear equation system using numerical algorithms (Bjorck, 1996); (5) perform error or covariance propagation by linearizing the equations obtained in step (3) with the estimated values. However, an alternative approach, namely Pope iteration, is followed. This alternately executes the following two until convergence: linearizing the nonlinear equations and standard linear estimation (Pope, 1972). The Pope approach, essentially a Gauss-Newton iterative procedure, has a merit as a by-product, i.e., the covariance matrix of the estimate can be readily retrieved during the solution. The linearization is done with known approximate values of the variables involved. Letting subscript or superscript ‘‘0” denote approximate values, we define the following increments as ‘‘d” values: t ¼ t0 dt; s ¼ s0 ds;
ð6Þ 1
R ¼ ðI 3 ½dhÞðI 3 þ ½dhÞ R0 ðI 3 2½dhÞR0 ;
ð7Þ
with the cross-product matrix of the 3 1 vector dh defined as 2 3 0 dh3 dh2 6 7 ½dh ¼ 4 dh3 0 dh1 5: ð8Þ dh2 dh1 0 The rotation matrix is the approximate expression used in the linearization. However, in correcting the approximate values using the estimated increments, the exact expression should be used instead, because it can always
Please cite this article in press as: Lin, P., et al. Helmert transformation with mixed geodetic and Cartesian coordinates. Adv. Space Res. (2017), https://doi.org/10.1016/j.asr.2017.11.029
4
P. Lin et al. / Advances in Space Research xxx (2017) xxx–xxx
ensure orthogonality of R given an orthogonal R0. In this exact expression, dh is just the Gibbs vector or Rodrigues parameters of the delta/increment rotation; here it is not the additive error of the Gibbs vector of the full rotation. For more detail, see e.g., Wang et al. (2016). Besides the above variables whose increments are to be estimated, approximate values of the following variables are needed in the linearization: e0~u ,e0~h ,e~0v and e0~h . 1 2 The approximate values used here should be the most accurate available before the subsequent estimation, in order to perform the best possible linearization. For unknowns, i.e., the transformation parameters and coordinates in the new frame, they should be estimates obtained from the immediately previous iteration. For the geodetic coordinates whose corresponding measurements are available, they should be the adjusted ones, i.e., the original measurements minus measurement error estimates obtained in the immediately previous iteration. Either of two candidates can be used to initialize the iteration. The first, based on practical considerations, one can simply let t0 = 0, s0 = 1, and R0 = I3. In the second, one can apply analytical methods from the literature to the common stations to attain a set of estimates of the parameters (e.g., Chang, 2015, 2016). These methods are used only to obtain parameter estimates, and one does not bother with a corresponding error analysis (Chang, 2016; Chang et al., 2017a,b). The first method is simple but the second
Eqs. (1) and (2) are used to perform the following linearization: L1;j ex~j ¼ B1j n s0 R0 G 1j el~j s0 R0 F 1j e~hj L2;j ¼ B 2j n yj s0 R0 G 2j e~vj s0 R0 F 2j e~hj 8 ~ t0 s0 R0 f ð~uj e0~uj ; ~h1;j e0~h Þ L1;j ¼ x > > 1;j > > > j 0 j 0 > < s0 R0 G 1 e~uj s0 R0 F 1 e~h 1;j ; with 0 ~ 0 > > > L2;j ¼ t0 s0 R0 f ð~vj e~vj ; h2;j e~h2;j Þ > > > : s R G j e0 s R F j e0 0
0
2 ~vj
0
0
ð9Þ ð10Þ
2 ~ h2;j
T T where u~j ¼ ½ u ~ 1;j ~k1;j ; ~vj ¼ ½ u ~ 2;j ~k2;j ; and 8 @f ð~ uj e0~u ;~ h1;j e0~ Þ > h1;j j j > > G1 ¼ > @eT > ~j l > e~uj ¼e0~u ;e~h ¼e0~ > h1;j > 1;j j > > > 0 ;~ 0 Þ > @f ð~ v e h e j > j ~vj 2;j ~ h2;j > > G2 ¼ > @e~T > vj > < e~vj ¼e~0v ;e~h ¼e0~ h2;j 2;j j ; 0 ;h 0 Þ ~ @f ð~ u e e > j ~ ~ uj 1;j > h1;j j > > F1 ¼ > @e~h > j > e~uj ¼e0~u ;e~h ¼e0~ > h1;j 1;j j > > > > 0 ;~ 0 Þ > @f ð~ v e h e j > ~vj 2;j ~ h2;j > > F 2j ¼ > @e~h > : j 0 0 e~vj ¼e~v ;e~h j
2;j
¼e~
h2;j
whose detailed expression is presented in Appendix A,
should be more accurate and can be used in cases of nonsmall-angle-rotation. In addition, the approximate values of measurement errors can be given the zero matrix of T ð6 m0 þ 3 m1 þ 3 m2 Þ 1. Letting n ¼ ½ dt ds dhT ,
For all stations, we can define the following variables:
Please cite this article in press as: Lin, P., et al. Helmert transformation with mixed geodetic and Cartesian coordinates. Adv. Space Res. (2017), https://doi.org/10.1016/j.asr.2017.11.029
P. Lin et al. / Advances in Space Research xxx (2017) xxx–xxx
2
2 6 6 6 F1 ¼ 6 6 4
3
s0 R0 G 11
6 6 6 G1 ¼ 6 6 4
s0 R0 G 21 ..
. s0 R0 G m1 0
6 6 6 G2 ¼ 6 6 4
s0 R0 F 21 ..
. s0 R0 F m1 0
6 6 6 F2 ¼ 6 6 4
1 ^ ^e ¼ Qee C T ðCQee C T Þ ðL AbÞ
3
s0 R0 G 22 ..
. s0 R0 G m2 1
s0 R0 F 12 s0 R0 F 22 ..
. s0 R0 F m2 1
7 7 7 7; 7 5 3 7 7 7 7: 7 5
Then we have 8 L1 þ G 1 e~u þ F 1 e~h1 ¼ B1 n > > > < L þ G e þ F e~ ¼ B n y 2 2 ~v 2 h2 1 > ~ z e ¼ z ~z > > : ~ ex~ ¼ x x
3. Monte Carlo simulation
ð11Þ
# ¼ xT
yT
z
2
T T
; b ¼ nT
#
T T
; L ¼ LT1
LT2
~T x
~zT
T
;
3
0 0 B1 0 6B 0 I 0 7 m2 6 2 7 A¼6 7; 4 0 I m1 0 0 5 0 0 0 I m3 2 0 0 G 1 F 1 6 0 0 G 2 F 2 6 C ¼6 4 0 0 0 0 0 0 0 0
I m1 0 I m1 0
0
3
0 7 7 7: 0 5 I m3
Finally, we have the overall measurement equation: L C e ¼ Ab
ð12Þ
According to the least-squares theory, the Lagrange loss function is expressed as: Uðe; b; KÞ ¼ eT Pe þ 2K T ðL C e AbÞ
The true transformation parameters from the new frame to the old were chosen as follows: ttrue = [100 5 200]T, strue = 0.89, and htrue = [0.035 0.029 0.033]T. The true value of the rotation matrix is Rtrue ¼ ðI 3 ½htrue ÞðI 3 þ ½htrue Þ
We define the following variables:
ð16Þ
The top part of (14) and the first four parts of (16) are 0 the A, C and L1 matrices in the iterative procedure. The bottom part of (14) and bottom-right part of (15) constitute the updated network of the new frame. As stated previously, the linearization-estimation pair should be performed iteratively; and in the intermediate iterations, only the parameter estimation (14) is necessary, while the variance-covariance matrix estimation (15) can be performed and stored in the final iteration. ^ represents coordinates of the new Because parameter # frame, whereas ^n symbolizes increments of Helmert transformation parameters, the object calculating the conver^ which ^ but ^n associated with D#, gence condition is not b ^ in the adjacent iteration. In addidenote the difference of # tion, in the iterative procedure, the scale parameter should be non-positive. If it is less than zero, it should be assigned to1again to avoid iterative divergence.
7 7 7 7; 7 5
s0 R0 G 12
2
ð15Þ
7 7 7 7; 7 5
3
s0 R0 F 11
2
5
ð13Þ
where K denotes the ð6 m0 þ 3 m1 þ 3 m2 Þ 1 vector of ‘‘Lagrange multipliers”. Then, the solution is readily expressed as: ð14Þ
1
Assuming m0 = 15, m1 = 20 and m2 = 20. All 45 stations are selected randomly, each of whose true geodetic coordinates in the old are generated as follows: (1) Generate longitude ktrue from the uniform distribution U[p p); (2) and latitude utrue from U[0.45p 0.45p]; 3) generate height htrue from U[10 1000]. Then, we calculate the true Cartesian coordinates in the new frame of the 15 common stations and 20 non-common stations surveyed only in the new frame: utrue xtrue ¼ ttrue þ strue Rtrue f ; htrue ktrue The distribution of all 45 stations in the new frame is illustrated in Fig. 2. The 105 105 matrix for all Cartesian coordinate measurements in the new frame, namely Q1, is generated as follows: (1) Designate the diagonal matrix with its diagonal elements the same as those of Q1 as D1, and generate the square-root of each of its 105 diagonal elements from U [0.05 0.15]; (2) generate an arbitrary 105 105 matrix X1 and perform the QR decomoposition as X1 = Q1R1; 3) let Q1 = Q1D1QT1 . The 60 60 matrix for all horizontal geodetic coordinates, namely Q2, and the 30 30 matrix for all vertical coordinate measurements, namely Q3, are generated in the same way; the only differences are such that
Please cite this article in press as: Lin, P., et al. Helmert transformation with mixed geodetic and Cartesian coordinates. Adv. Space Res. (2017), https://doi.org/10.1016/j.asr.2017.11.029
6
P. Lin et al. / Advances in Space Research xxx (2017) xxx–xxx
Fig. 2. 3D distribution of stations involved in simulation. Coordinates are Cartesian in the new frame. Red, green, and blue spheres denote common stations, non-common stations surveyed in the new frame, and non-common stations surveyed in the old frame, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
0.10
0.40 0.35
12
RMS (ppm)
0.25 0.20 0.15 0.10
RMS (mas)
0.08
0.30
RMS (m)
14
0.06 0.04
8 6 4
0.02
2
0.05 0.00
10
0.00
0
translations
scale
rotation angles
Fig. 3. Transformation parameter estimation accuracies in terms of RMSE. Blue and red bars represent stepwise and direct approaches, respectively. Unit length along the vertical axis represents 1 m, 1 ppm (part per million) and 1 mas (milli-arc-second), for the translation, scale and rotation parameters, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
0.10 0.08 0.06
Z component
Y component
X component
0.12
0.04
0.16
0.16
0.14
0.14 0.12 0.10 0.08
4 6
8 10 12 14 16
common station
0.12 0.10 0.08 0.06
0.06
0 2
0.14
0
2
4
6
8 10 12 14 16
common station
0
2
4
6
8 10 12 14 16
common station
Fig. 4. Cartesian coordinate estimation accuracies in terms of RMSE for the 15 common stations. Blue, red and yellow bars represent estimates by the stepwise approach, by the direct approach, and the original measurements, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
for the former and latter, the distributions from which the square-root of their diagonal elements are generated are U [1.5 3]107 radians and U[1 5] meters, respectively. The results are " # Q~h1 ~h1 Q~h1 ~h2 Qx~ Qu~~u Qu~~v ~x Q x~ ~z Q1 ¼ ;Q2 ¼ ;Q3 ¼ ; Q~h2 ~h1 Q~h2 ~h2 Qz~~x Qz~~z Qv~~u Qv~~v Then, the total variance covariance Qee can be constructed using Q1, Q2 and Q3.
Overall, 100 Monte Carlo experiments were conducted. In all experiments, the above true parameters and coordinates were fixed. Measurement errors were generated independently across different experiments from zero-mean Gaussian distributions with the above variance covariance matrices (Q1, Q2 and Q3). Two approaches were evaluated, the direct approach proposed in this paper and stepwise approach. To be selfcontained, the stepwise approach is presented in Appendix B. The following estimation errors for both approaches
Please cite this article in press as: Lin, P., et al. Helmert transformation with mixed geodetic and Cartesian coordinates. Adv. Space Res. (2017), https://doi.org/10.1016/j.asr.2017.11.029
7
3.0 3.0
2.5
2.5
Z component
3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0
Y componet
X componet
P. Lin et al. / Advances in Space Research xxx (2017) xxx–xxx
2.0 1.5 1.0 0.5 0.0
0 1 2 3 4 5 6 7 8 9 10 11
non-common station (old frame)
2.0 1.5 1.0 0.5 0.0
0 1 2 3 4 5 6 7 8 9 10 11
0 1 2 3 4 5 6 7 8 9 10 11
non-common station (old frame)
non-common station (old frame)
0 2 4 6 8 10 12 14 16 18 20
non-common station(new frame)
0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00
Z component
0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00
Y component
X component
Fig. 5. Cartesian coordinates estimation accuracies in terms of RMSE for the 10 non-common stations surveyed only in the old frames. Blue, red and yellow bars represent estimates by stepwise approach, direct approach, and transforming measured geodetic coordinates. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
0 2 4 6 8 10 12 14 16 18 20
non-common station(new frame)
0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00
0 2 4 6 8 10 12 14 16 18 20
non-common station(new frame)
Fig. 6. Cartesian coordinate estimation accuracies in terms of RMSE for the 20 non-common stations surveyed only in the new frame. Blue, red and yellow bars represent estimates by stepwise approach, direct approach, and original measurements, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
were checked: translation parameter estimation errors ^t ttrue ; scale parameter estimation errors ^s strue ; rotation parameter estimation errors (in radians), whose cross
T ^ R ^ Rtrue ^ Rtrue R ^T þ R product is 12 R ; finally, coordinate estimation errors for all 45 stations x xtrue . To weaken random effects, root mean squares of all error terms above were calculated across the 100 Monte Carlo runs. These root mean squared errors (RMSE) are exactly the indices representing accuracies of the corresponding estimates. Bar charts were used to explicitly compare the performance of the two approaches. The results are presented in Figs. 3–6. These correspond to the transformation parameter, coordinates of the 15 common stations, coordinates of the 10 non-common stations surveyed only in the old frame, and coordinates of the 20 non-common stations surveyed only in the new frame, respectively. In all figures, blue bars represent RMSEs of the stepwise approach, and red bars represent those of the direct approach. In Figs. 4 and 6, RMSEs for the originally measured coordinates are depicted by yellow bars. In Fig. 5, RMSEs for coordinates obtained by transforming the measured geodetic coordinates using the estimated parameters are also shown, by yellow bars. In Fig. 3, the unit length along the vertical axis represents 1 m, 1 ppm (part per million) and 1 mas (milli-arc-second), for the translation, scale and rotation parameters, respectively. Horizontal axes in Figs. 4–6 represent indices of the stations.
Fig. 3 shows that although performances of the two approaches were comparable in estimating the rotation parameters, the direct approach performed better than the stepwise approach in estimating the translation and scale parameters. Figs. 4–6 reveal the following. First, for all stations with their coordinates directly measured in the new frame (i.e., the stations in Figs. 4 and 6), the adjusted/estimated coordinates using either of the two approaches are more accurate than the measured coordinates. Second, for the non-common stations surveyed in the old frame, (those in Fig. 5), the adjusted coordinates are more accurate than those transformed using the estimated transformation parameters. These two observations clearly show the necessity of incorporating the coordinates as parameters in addition to the transformation parameters in the Helmert transformation problem. This again validates the statements made in Kotsakis et al. (2014). Third, for most of the coordinates, the direct approach did better than the stepwise one, though the former could not outperform the latter for all stations or coordinates. Of course, there were still some stations or coordinate components for which the stepwise approach outperformed the direct approach. This appears somewhat irregular. Even more irregular is that there were some coordinate components for which the original measurement was more accurate than estimates by both approaches. It is believed that the randomness effect was compensated by the Monte Carlo experiments. Therefore, this ‘‘irregular” phenomenon
Please cite this article in press as: Lin, P., et al. Helmert transformation with mixed geodetic and Cartesian coordinates. Adv. Space Res. (2017), https://doi.org/10.1016/j.asr.2017.11.029
8
P. Lin et al. / Advances in Space Research xxx (2017) xxx–xxx
cannot be attributed to the randomness effect. Further studies are needed to interpret the observed phenomenon and to possibly further optimize the proposed approach. 4. Concluding remarks Given geodetic coordinates in one frame and Cartesian coordinates in the other, a potential approach to performing the Helmert transformation is to first convert geodetic coordinates to their Cartesian counterparts within the same frame and then do a Helmert transformation with Cartesian coordinates in both frames. However, in this work, we clearly show that this stepwise approach is not globally optimal. A direct approach is proposed instead, using the geodetic coordinates directly. By recalling that one of the tasks in the Helmert transformation problem, possibly the most important one, is to estimate/adjust coordinates of all relevant stations in the new frame, and that measured coordinates of different stations may be correlated, the augmented problem suggested by Kotsakis et al. (2014) is solved. This was accomplished by treating the coordinates of all stations in the new frame as parameters, in addition to the transformation parameters. The least-squares criterion was used to solve the problem, and the Pope iteration used to obtain the solution numerically. Monte Carlo simulations were conducted, and the results validate our theoretical statements; overall performance of the proposed direct method is superior to that of the stepwise one. Acknowledgments This work was supported by the following funds: National Natural Science Foundation of China (No. 41774005 and 41674008) and Research Innovation Program for College Graduates of Jiangsu Province (KYZZ16_0217, SJCX17_0522). Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j. asr.2017.11.029. References Bjorck, A., 1996. Numerical Methods for Least Squares Problems. SIAM, Philadelphia. Chang, G., 2015. On least-squares solution to 3D similarity transformation problem under Gauss-Helmert model. J. Geodesy 89, 573–576. Chang, G., 2016. Closed form least-squares solution to 3D symmetric Helmert transformation with rotational invariant covariance structure. Acta Geodaetica et Geophysica 51, 237–244.
Chang, G., Xu, T., Wang, Q., 2017a. Error analysis of the 3D similarity coordinate transformation. GPS Solut. 21 (3), 963–971. Chang, G., Xu, T., Wang, Q., Zhang, S., et al., 2017b. A generalization of the analytical least-squares solution to the 3D symmetric Helmert coordinate transformation problem with an approximate error analysis. Adv. Space Res. 59 (10), 2600–2610. Fang, X., 2013. Weighted total least squares: necessary and sufficient conditions, fixed and random parameters. J. Geodesy 87, 733–749. Fang, X., 2015. Weighted total least-squares with constraints: a universal formula for geodetic symmetrical transformations. J. Geodesy 89, 459– 469. Grafarend, E., Okeke, F., 1998. Transformation of conformai coordinates of type Mercator from a global datum (WGS 84) to a local datum (Regional, national). Mar. Geodesy 21 (3), 169–180. Kay, S.M., 2013. Fundamentals of Statistical Signal Processing, Volume I: Estimation Theory. Pearson Education, New York. Kotsakis, C., Vatalis, A., Sanso, F., 2014. On the importance of intraframe and inter-frame covariances in frame transformation theory. J. Geodesy 88, 1187–1201. Lehmann, R., 2014. Transformation model selection by multiple hypotheses testing. J. Geodesy 88, 1117–1130. Leick, A., Rapoport, L., Tatarnikov, D., 2015. GPS Satellite Survey, fourth ed. Wiley, New Jersey. Li, B., Shen, Y., Li, W., 2012. The seamless model for three-dimensional datum transformation. Sci. China Earth Sci. 55, 2099–2108. Li, B., Shen, Y., Zhang, X., Li, C., Lou, L., 2013. Seamless multivariate affine error-in-variables transformation and its application to map rectification. Int. J. Geogr. Inform. Sci. 27, 1572–1592. Lu, Z., Qu, Y., Qiao, S., 2014. Geodesy: Introduction to Geodetic Datum and Geodetic Systems. Springer, Berlin. Neitzel, F., 2010. Generalization of total least-squares on example of unweighted and weighted 2D similarity transformation. J. Geodesy 84, 751–762. Pope, A.J., 1972. Some pitfalls to be avoided in the iterative adjustment of nonlinear problems. In: Proceedings of the 38th Annual Meeting. American Society of Photogrammetry, Washington DC, pp. 449–477. Schaffrin, B., Felus, Y.A., 2008. On the multivariate total least-squares approach to empirical coordinate transformations. Three algorithms. J. Geodesy 82, 373–383. Seeber, G., 2003. Satellite Geodesy. Walter de Gruyter, Berlin. Vanı´cˇek, P., Steeves, R.R., 1996. Transformation of coordinates between two horizontal geodetic datums. J. Geodesy 70 (11), 740–745. Wang, B., Li, J., Liu, C., Yu, J., 2017. Generalized total least squares prediction algorithm for universal 3D similarity transformation. Adv. Space Res. 59, 815–823. Wang, Q., Chang, G., Xu, T., Zou, Y., 2016. Representation of the rotation parameter estimation errors in the Helmert transformation model. Surv. Rev., 1–13 https://doi.org/10.1080/ 00396265.2016.1234806. Xu, P., Liu, J., Shi, C., 2012. Total least squares adjustment in partial errors-in-variables models: algorithm and statistical analysis. J. Geodesy 86, 661–675. Xue, S., Dang, Y., Liu, J., Mi, J., Dong, C., Cheng, Y., Wang, X., Wan, J., 2016a. Bias estimation and correction for triangle-based surface area calculations. Int. J. Geogr. Inform. Sci. 30, 2155–2170. Xue, S., Yang, Y., Dang, Y., 2016b. Formulas for precisely and efficiently estimating the bias and variance of the length measurements. J. Geogr. Syst. 18, 399–415.
Please cite this article in press as: Lin, P., et al. Helmert transformation with mixed geodetic and Cartesian coordinates. Adv. Space Res. (2017), https://doi.org/10.1016/j.asr.2017.11.029