A comparison of different solutions of the Bursa–Wolf
model and of the 3D, 7-parameter datum transformation
Jo
´
zsef Za
´
voti
1
•
Ja
´
nos Kalma
´
r
1
Received: 18 March 2015 / Accepted: 12 June 2015 / Published online: 11 July 2015
Ó Akade
´
miai Kiado
´
2015
Abstract The present work deals with an important theoretical problem of geodesy: we are
looking for a mathematical dependency between two spatial coordinate systems utilizing
common pairs of points whose coordinates are given in both systems. In geodesy and pho-
togrammetry the most often used procedure to move from one coordinate system to the other is
the 3D, 7 parameter (Helmert) transformation. Up to recent times this task was solved either by
iteration, or by applying the Bursa–Wolf model. Producers of GPS/GNSS receivers install
these algorithms in their systems to achieve a quick processing of data. But nowadays algebraic
methods of mathematics give closed form solutions of this problem, which require high level
computer technology background. In everyday usage, the closed form solutions are much more
simple and have a higher precision than earlier procedures and thus it can be predicted that these
new solutions will find their place in the practice. The paper discusses various methods for
calculating the scale factor and it also compares solutions based on quaternion with those that
are based on rotation matrix making use of skew-symmetric matrix.
Keywords Quaternion-algebra Bursa–Wolf model Rotation matrix Scale factor
3D or 7-parameter datum (Helmert) transformation
1 Introduction
The conventional treatment of the 3D, 7-parameter datum transformation is given in
Grafarend and Krumm (1995), in Grafarend and Kampmann (1996), and in Grafarend and
Shan (1997). Subsequently Awange et al. (2004) has added extensions to the solutions.
& Jo
´
zsef Za
´
voti
zavoti@ggki.hu;
https://www.ggki.hu
Ja
´
nos Kalma
´
r
kalmar.janos@csfk.mta.hu
1
Research Centre for Astronomy and Earth Sciences, Hungarian Academy of Sciences,
Csatkai u. 6-8, Sopron 9400, Hungary
123
Acta Geod Geophys (2016) 51:245–256
DOI 10.1007/s40328-015-0124-6
Za
´
voti (1999) has given a solution with limited conditions in L1 norm. In the application of
computer algebraic systems to datum transformations, Awange and Grafarend (2002,
2003a, b, c) have taken new directions. In Hungary, Za
´
voti (2005) gives the first algebraic
solution to the problem, and his solution also proposes correction to the mathematical
model. Za
´
voti and Jancso
´
(2006) give the basic idea for bringing the problem to linear
form, and this concept is described in more detail in Za
´
voti (2012). Battha and Za
´
voti
(2009a, b) have applied computerized algebra to the problem of the intersection problem.
Za
´
voti and Fritsch (2011) have given a totally new solution for the outer orientation
problem of photogrammetry. Horn (1987) is one of the earliest works to give a solution to
the absolute orientation problem, but his solution is different from Za
´
voti (2012). Za
´
voti
and Kalma
´
r (2014) give a good summary of the differences between the two solutions.
2 The model of the new solution for the 3D, 7-parameter similarity
transformation
Suppose that we have two distinct coordinate systems with n common points given by their
coordinates.
The 3D, 7-parameter (Helmert) overdetermined similarity transformation is given by
the matrix equation:
s
i
¼ t þ kRp
i
; i ¼ 1; 2; ...; n; ð1Þ
where s
i
¼ X
i
; Y
i
; Z
i
½
T
are the target coordinates, t ¼ X
0
; Y
0
; Z
0
½
T
is the unknown shift, k is
the unknown scale-factor, R(a, b, c) is the rotation matrix, p
i
¼ x
i
; y
i
; z
i
½
T
are the coor-
dinates of the object points.
The R rotation matrix is parametrized using the three independent and unknown a, b
and c Cardan-angles (Awange 2002), which belong to rotations about axes z-y-x by the
angle c, b and a consecutively.
R ¼ R
1
aðÞR
2
bðÞR
3
cðÞ: ð2Þ
Obviously, changing the order of the three rotations or the direction of the angles leads to
rotation matrix.
The rotation angles can be obtained from the elements of the rotation matrix:
a ¼arctan
r
23
r
33
; b ¼ arcsinðr
13
Þ; c ¼arctan
r
12
r
11
; ð3Þ
where r
ij
is the j-th entry in the i-th row of the matrix R. Thus our goal is the determination
of the rotation matrix.
Awange and Grafarend (2002) have introduced the (5) skew-symmetric matrix C
0
that
has the property:
R ¼ I
3
C
0
1
I
3
þ C
0
; ð4Þ
where I
3
is the 3-dimesional identity matrix, and C
0
is given by
C
0
¼
0 cb
c 0 a
ba 0
2
4
3
5
; ð5Þ
with parameters a, b and c.
246 Acta Geod Geophys (2016) 51:245–256
123
Multiplying from the left Eq. (1) with I
3
C
0
, and applying Eq. (4), we obtain:
1 c b
c 1 a
b a 1
2
6
4
3
7
5
X
i
Y
i
Z
i
2
6
4
3
7
5
¼
1 c b
c 1 a
b a 1
2
6
4
3
7
5
X
0
Y
0
Z
0
2
6
4
3
7
5
þ k
1 cb
c 1 a
ba 1
2
6
4
3
7
5
x
i
y
i
z
i
2
6
4
3
7
5
;
i ¼ 1; 2; ...; n:
ð6Þ
These equations are the basis of the algebraic solution of the 3D, 7-parameter (Helmert)
transformation.
3 Determination of the scale-factor of the 3D, 7-parameter similarity
transformation
Za
´
voti (2012) has eliminated the shift parameters by reducing the coordinates to the center
of gravity. He also has shown that during the solution of the overdetermined system of
equations parameters a, b and c are eliminated, and the following overdetermined system
of equations, quadratic in the unknown k parameter, are obtained:
k
2
x
2
is
þ y
2
is
þ z
2
is
¼ X
2
is
þ Y
2
is
þ Z
2
is
; i ¼ 1; 2; ...; n; ð7Þ
where X
is
¼ X
i
X
s
; Y
is
¼ Y
i
Y
s
; Z
is
¼ Z
i
Z
s
i ¼ 1; 2; ...; n; x
is
¼ x
i
x
s
; y
is
¼
y
i
y
s
; z
is
¼ z
i
z
s
; i ¼ 1; 2; ...; n; X
s
; Y
s
; Z
s
ðÞ; x
s
; y
s
; z
s
ðÞare coordinates the centre of
gravity.
The system of Eq. (7) is overdetermined, with several solutions according to a chosen
error function.
3.1 Solution I
The (7) system of equations is written as a product:
k
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
x
2
is
þ y
2
is
þ z
2
is
q
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
X
2
is
þ Y
2
is
þ Z
2
is
q
k
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
x
2
is
þ y
2
is
þ z
2
is
q
þ
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
X
2
is
þ Y
2
is
þ Z
2
is
q
¼ 0;
i ¼ 1; 2; ...; n: ð8Þ
Let us consider the first factor in the above product. The system of equations to be
solved is:
k
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
x
2
is
þ y
2
is
þ z
2
is
q
¼
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
X
2
is
þ Y
2
is
þ Z
2
is
q
; i ¼ 1; 2; ...; n: ð9Þ
Add up all these equations! The solution for k of the overdetermined system of equa-
tions—taking into account that only the positive root has a physical meaning for us—is
given, according to Za
´
voti (2012), by the following well known equation:
k
1
¼
P
n
i¼1
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
X
2
is
þ Y
2
is
þ Z
2
is
p
P
n
i¼1
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
x
2
is
þ y
2
is
þ z
2
is
p
: ð10Þ
Acta Geod Geophys (2016) 51:245–256 247
123
Albertz and Kreiling (1975) have shown that the k scale-factor can be obtained as the
quotient of the sums of the coordinates of the points, taken in the coordinate system with
origin in the center of gravity. Thus we have transformed the quadratic Eq. (7) to linear
equations, as opposed to the procedure given in Awange and Grafarend (2002), which
requires the cumbersome separation of the roots of a polynomial of degree 4.
3.2 Solution II
Adding up all equations in (7) gives
k
2
X
n
i¼1
x
2
is
þ y
2
is
þ z
2
is
¼
X
n
i¼1
X
2
is
þ Y
2
is
þ Z
2
is
: ð11Þ
There is a simple solution of this equation in non-negative real numbers without trans-
forming it to a product. After taking into account that a root with physical meaning must be
positive, we obtain for the lambda scale-factor the following equation, given in Horn
(1987) using quaternions, which is also the solution of the Bursa–Wolf model:
k
2
¼
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
P
n
i¼1
X
2
is
þ Y
2
is
þ Z
2
is
ðÞ
P
n
i¼1
x
2
is
þ y
2
is
þ z
2
is
ðÞ
v
u
u
u
u
u
t
: ð12Þ
Thus we can obtain the unique k scale-factor from the quadratic equations as opposed to
the complicated separation procedure of the roots of a polynomial of degree 4 given by
Awange and Grafarend (2002).
3.3 Solution III
Our starting point is again the system of Eq. (7). We want to obtain a least squares solution
for k using intermediary equations. Elementary steps, given in full detail in (23)–(26),
yield:
k
3
¼
P
n
i¼1
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
x
2
is
þ y
2
is
þ z
2
is
ðÞX
2
is
þ Y
2
is
þ Z
2
is
ðÞ
p
P
n
i¼1
x
2
is
þ y
2
is
þ z
2
is
ðÞ
: ð13Þ
Thus there are three different solution procedures (estimations) for the k scale factor of the
3D, 7-parameter (Helmert) transformation.
4 Determination of the linear and shift parameters
After having determined the scale-factor, the problem can be written in linear form, and the
adjustment model for the linear problem can be given. This procedure makes it possible to
include arbitrarily many equations (common points), and give a solution for the parameters
a, b, and c.
Za
´
voti (2013) has determined the normal matrix of the problem:
248 Acta Geod Geophys (2016) 51:245–256
123
P
n
i¼1
ky
is
þ Y
is
ðÞ
2
þ kz
is
þ Z
is
ðÞ
2
hi
P
n
i¼1
kx
is
þ X
is
ðÞky
is
þY
is
ðÞ
P
n
i¼1
kx
is
þX
is
ðÞkz
is
þ Z
is
ðÞ
P
n
i¼1
kx
is
þ X
is
ðÞ
2
þ kz
is
þ Z
is
ðÞ
2
hi
P
n
i¼1
ky
is
þ Y
is
ðÞkz
is
þ Z
is
ðÞ
P
n
i¼1
kx
is
þ X
is
ðÞ
2
þ ky
is
þ Y
is
ðÞ
2
hi
2
6
6
6
6
6
6
4
3
7
7
7
7
7
7
5
:
ð14Þ
(The symmetric entries of the normal matrix are not listed).
The normal vector is obtained in a similar way:
2k
P
n
i¼1
ðy
is
Z
is
z
is
Y
is
Þ
P
n
i¼1
ðz
is
X
is
x
is
Z
is
Þ
P
n
i¼1
x
is
Y
is
y
is
X
is
ðÞ
2
6
6
6
6
6
6
4
3
7
7
7
7
7
7
5
: ð15Þ
There are several procedures to obtain the parameters a, b and c from the 3 9 3 normal
system of equations, we have chosen the singular value decomposition because of its
stability. Utilizing special properties of the normal matrix the rotation parameters in (3) can
also be determined.
The shift parameters X
0
, Y
0
and Z
0
can be determined after Eq. (1) is rewritten to a
system with the center of gravity as origin
X
0
Y
0
Z
0
2
4
3
5
¼
X
s
Y
s
Z
s
2
4
3
5
kR
x
s
y
s
z
s
2
4
3
5
: ð16Þ
Parameters for the precision, variance and covariance are computed in the conventional
way.
5 The Bursa–Wolf model of the datum-transformation
Equation (13) can also be obtained in the following way (
s and
p are the centers of gravity
in the two systems):
Ds
i
¼ s
i
s ) s
i
¼ Ds
i
þ
s ;
Dp
i
¼ p
i
p ) p
i
¼ Dp
i
þ
p :
ð17Þ
Substituting this into the transformation Eq. (1) gives:
Ds
i
þ
s ¼ t þ k R Dp
i
þ
p
ðÞ
; i ¼ 1; 2; ...; n: ð18Þ
Reordering terms gives:
Ds
i
þ
s ¼ t þk R
p þ k R Dp
i
; i ¼ 1; 2; ...; n: ð19Þ
Since (1) holds for the centers of gravity
s and
p too, so the term in the middle of (19)
can be omitted and thus we are left with:
Ds
i
¼ k R Dp
i
; i ¼ 1; 2; ...; n: ð20Þ
Acta Geod Geophys (2016) 51:245–256 249
123