Positioning and Navigation Using the Russian Satellite System

R
− y
S
)
2
+ (z
R
− z
S
)
2
(6.3.2)
The position vectors of receiver x
R
and satellite x
S
in Eq. (6.3.2) must be given in the same coordinate
frame.
In the given case of direct determination of transformation parameters between PZ-90 and WGS84
from range measurements to GLONASS satellites, the known coordinates of the observation sites are
given in WGS84, whereas the coordinates of the GLONASS satellites are determined from ephemeris
data and are given in PZ-90. Since both WGS84 and PZ-90 are Earth-centered Earth-fixed (ECEF)
coordinate frames, they rotate along with the Earth. Their orientation therefore is a function of time.
More precise, thus, the coordinates of the observation sites are given in WGS84, as valid at the time
of signal reception t
RX
, and the coordinates of the satellites are given in PZ-90, as valid at the time of
signal transmission t
T X
.
x
R
= x
R,W GS(t
RX
)
x
S
= x
S
P Z(t
T X
)
To obtain the actual geometrical distance from receiver to satellite, one of these two sets of coordinates
must be transformed to the coordinate frame of the other set. Here it is chosen to transform the satellite
coordinates to the coordinate frame of the observer. This requires two steps of transformation:
•
Transformation from PZ-90 to WGS84: x
S
P Z(t
T X
)
→ x
S
W GS(t
T X
)
•
Correction of Earth rotation: x
S
W GS(t
T X
)
→ x
S
W GS(t
RX
)
Correction of Earth rotation
While the satellite signal is travelling towards the observer, the Earth – and along with it the Earth-
fixed coordinate frame – keeps rotating. During this signal travel time, it rotates by an angle of α =
S
R
/c · ω
E
. This is a positive rotation around the z-axis. Thus, the satellite coordinates transform by
x
S
P Z(t
RX
)
=



1 α 0
−α
1 0
0
0 1


 · x
S
P Z(t
T X
)
(6.3.3)
for small angles α.
The geometrical distance therefore becomes
S
R
=
x
R,P Z(t
RX
)
−x
S
P Z(t
RX
)
2
+ y
R,P Z(t
RX
)
−y
S
P Z(t
RX
)
2
+ z
R,P Z(t
RX
)
−z
S
P Z(t
RX
)
2
(6.3.4)

64
6 DETERMINATION OF TRANSFORMATION PARAMETERS
or
S
R
2
=
x
R,P Z(t
RX
)
−x
S
P Z(t
T X
)
−
ω
E
c
S
R
· y
S
P Z(t
T X
)
2
+
(6.3.5)
y
R,P Z(t
RX
)
−y
S
P Z(t
T X
)
+
ω
E
c
S
R
· x
S
P Z(t
T X
)
2
+ z
R,P Z(t
RX
)
−z
S
P Z(t
T X
)
2
=
x
R,P Z(t
RX
)
−x
S
P Z(t
T X
)
2
−2 x
R,P Z(t
RX
)
−x
S
P Z(t
T X
)
ω
E
c
S
R
y
S
P Z(t
T X
)
+
ω
2
E
c
2
S
R
2
y
S
P Z(t
T X
)
2
+
y
R,P Z(t
RX
)
−y
S
P Z(t
T X
)
2
+2 y
R,P Z(t
RX
)
−y
S
P Z(t
T X
)
ω
E
c
S
R
x
S
P Z(t
T X
)
+
ω
2
E
c
2
S
R
2
x
S
P Z(t
T X
)
2
+
z
R,P Z(t
RX
)
−z
S
P Z(t
T X
)
2
(6.3.6)
=
x
R,P Z(t
RX
)
−x
S
P Z(t
T X
)
2
+ y
R,P Z(t
RX
)
−y
S
P Z(t
T X
)
2
+ z
R,P Z(t
RX
)
−z
S
P Z(t
T X
)
2
+
(6.3.7)
ω
2
E
c
2
S
R
2
· x
S
P Z(t
T X
)
2
+y
S
P Z(t
T X
)
2
−2
ω
E
c
S
R
· x
R,P Z(t
RX
)
y
S
P Z(t
T X
)
−y
R,P Z(t
RX
)
x
S
P Z(t
T X
)
This transforms into
0 =
S
R
2
· 1 −
ω
E
c
2
x
S
P Z(t
T X
)
2
+ y
S
P Z(t
T X
)
2
+
2
ω
E
c
S
R
· x
R,P Z(t
RX
)
y
S
P Z(t
T X
)
− y
R,P Z(t
RX
)
x
S
P Z(t
T X
)
−
(6.3.8)
x
R,P Z(t
RX
)
− x
S
P Z(t
T X
)
2
− y
R,P Z(t
RX
)
− y
S
P Z(t
T X
)
2
− z
R,P Z(t
RX
)
− z
S
P Z(t
T X
)
2
The term
ω
E
c
2
· x
S
P Z(t
T X
)
2
+ y
S
P Z(t
T X
)
2
in Eq. (6.3.8) is of the order of magnitude 10
−12
and can
be neglected with respect to the 1:
0 =
S
R
2
+ 2
ω
E
c
S
R
· x
R,P Z(t
RX
)
y
S
P Z(t
T X
)
− y
R,P Z(t
RX
)
x
S
P Z(t
T X
)
−
(6.3.9)
x
R,P Z(t
RX
)
− x
S
P Z(t
T X
)
2
− y
R,P Z(t
RX
)
− y
S
P Z(t
T X
)
2
− z
R,P Z(t
RX
)
− z
S
P Z(t
T X
)
2
This quadratic equation has two possible solutions:
S
R1,2
= −
ω
E
c
· x
R,P Z(t
RX
)
y
S
P Z(t
T X
)
− y
R,P Z(t
RX
)
x
S
P Z(t
T X
)
±
ω
E
c
2
· x
R,P Z(t
RX
)
y
S
P Z(t
T X
)
−y
R,P Z(t
RX
)
x
S
P Z(t
T X
)
2
+
(6.3.10)
x
R,P Z(t
RX
)
−x
S
P Z(t
T X
)
2
+ y
R,P Z(t
RX
)
−y
S
P Z(t
T X
)
2
+ z
R,P Z(t
RX
)
−z
S
P Z(t
T X
)
2 1/2
Since the square root is always positive and of larger magnitude than the term in front of it, the
negative sign in front of the square root does not provide a physically meaningful solution. Therefore the
distance between satellite and observer, corrected for Earth rotation during the signal travel time, can
be written as:
S
R
= −
ω
E
c
· x
R,P Z(t
RX
)
y
S
P Z(t
T X
)
−y
R,P Z(t
RX
)
x
S
P Z(t
T X
)
+
ω
E
c
2
· x
R,P Z(t
RX
)
y
S
P Z(t
T X
)
−y
R,P Z(t
RX
)
x
S
P Z(t
T X
)
2
+
(6.3.11)
x
R,P Z(t
RX
)
−x
S
P Z(t
T X
)
2
+ y
R,P Z(t
RX
)
−y
S
P Z(t
T X
)
2
+ z
R,P Z(t
RX
)
−z
S
P Z(t
T X
)
2 1/2
For deriving Eq. (6.3.11), the PZ-90 coordinate frame was used. Of course, Eq. (6.3.11) is also valid in
the WGS84 coordinate frame:
S
R
= −
ω
E
c
· x
R,W GS(t
RX
)
y
S
W GS(t
T X
)
−y
R,W GS(t
RX
)
x
S
W GS(t
T X
)
+
ω
E
c
2
· x
R,W GS(t
RX
)
y
S
W GS(t
T X
)
−y
R,W GS(t
RX
)
x
S
W GS(t
T X
)
2
+
(6.3.12)
x
R,W GS(t
RX
)
−x
S
W GS(t
T X
)
2
+ y
R,W GS(t
RX
)
−y
S
W GS(t
T X
)
2
+ z
R,W GS(t
RX
)
−z
S
W GS(t
T X
)
2 1/2

6.3 Direct Estimation of Transformation Parameters
65
Transformation from PZ-90 to WGS84
The transformation of a position vector given in PZ-90 to WGS84 can be written as (cf. Section 5.5):



x
y
z



W GS
=



∆x
∆y
∆z


 + (1 + δs) ·



1
δω
−δψ
−δω
1
δε
δψ
−δε
1


 ·



x
y
z



P Z
(6.3.13)
or in individual coordinates, valid at time t
T X
:
x
S
W GS(t
T X
)
= ∆x + (1 + δs) · x
S
P Z(t
T X
)
+ δω · y
S
P Z(t
T X
)
− δψ · z
S
P Z(t
T X
)
y
S
W GS(t
T X
)
= ∆y + (1 + δs) · −δω · x
S
P Z(t
T X
)
+ y
S
P Z(t
T X
)
+ δε · z
S
P Z(t
T X
)
(6.3.14)
z
S
W GS(t
T X
)
= ∆z + (1 + δs) · δψ · x
S
P Z(t
T X
)
− δε · y
S
P Z(t
T X
)
+ z
S
P Z(t
T X
)
Inserting this into Eq. (6.3.12) yields
S
R
= −
ω
E
c
x
R,W GS(t
RX
)
∆y + (1+δs) · −δω · x
S
P Z(t
T X
)
+ y
S
P Z(t
T X
)
+ δε · z
S
P Z(t
T X
)
−
y
R,W GS(t
RX
)
∆x + (1+δs) · x
S
P Z(t
T X
)
+ δω · y
S
P Z(t
T X
)
− δψ · z
S
P Z(t
T X
)
+
ω
E
c
2
x
R,W GS(t
RX
)
∆y + (1+δs) · −δω · x
S
P Z(t
T X
)
+ y
S
P Z(t
T X
)
+ δε · z
S
P Z(t
T X
)
−
y
R,W GS(t
RX
)
∆x + (1+δs) · x
S
P Z(t
T X
)
+ δω · y
S
P Z(t
T X
)
− δψ · z
S
P Z(t
T X
)
2
+
(6.3.15)
x
R,W GS(t
RX
)
− ∆x + (1+δs) · x
S
P Z(t
T X
)
+ δω · y
S
P Z(t
T X
)
− δψ · z
S
P Z(t
T X
)
2
+
y
R,W GS(t
RX
)
− ∆y + (1+δs) · −δω · x
S
P Z(t
T X
)
+ y
S
P Z(t
T X
)
+ δε · z
S
P Z(t
T X
)
2
+
z
R,W GS(t
RX
)
− ∆z + (1+δs) · δψ · x
S
P Z(t
T X
)
− δε · y
S
P Z(t
T X
)
+ z
S
P Z(t
T X
)
2 1/2
Here, the coordinates x
R,W GS(t
RX
)
are the known station coordinates, given in WGS84, valid at
the time of signal reception. x
S
P Z(t
T X
)
are the satellite coordinates at the time of signal transmission as
obtained from ephemeris data. They are given in the PZ-90 frame, valid at the time of signal transmission.
The transformation parameters ∆x, ∆y, ∆z, δs, δε, δψ, δω are unknown in this geometrical range in
different coordinate frames Eq. (6.3.15). Inserting this geometrical range into the observation equation
Eq. (6.3.1), these unknowns can be solved for, provided there is a sufficient number of observations.
However, the receiver clock error δt
R
in Eq. (6.3.1) is also unknown. The satellite clock error δt
S
can be determined from ephemeris data, whereas tropospheric and ionospheric path delays δt
S,T rop
R
and
δt
S,Iono
R
can be modelled. Since GLONASS offers free and unobstructed access to the second frequency,
dual-frequency ionospheric corrections can be applied alternatively.
This leaves seven unknown transformation parameters and one unknown receiver clock error to solve
for. Thus, with range measurements to eight GLONASS satellites at one observation site, a complete
set of transformation parameters could be determined. Besides having eight GLONASS satellites in view
is rather unlikely in times of depleted GLONASS constellation, this one observation site will provide a
poor geometry to separate origin and orientation parameters. As discussed in Section 5.6.1, this will lead
to a set of parameters that are only valid at a small area around the observation site. More stations
will add more strength to the geometry. However, each additional observation site does also mean an
additional receiver and thus one more receiver clock error as a further unknown. For simultaneous
observations at five stations, as was the case during the IfEN/IfAG/DLR measurement campaign, this
totals in twelve unknowns. Thus, with two to three observations at each station, it was possible to
determine the transformation parameters directly from pseudorange observations. With seven or more
observation sites, which can easily be reached with data from the IGEX campaign, only two satellites
in view per site are required. This is another bonus of this method for determining the transformation

66
6 DETERMINATION OF TRANSFORMATION PARAMETERS
parameters: For the conventional method of determining point coordinates in PZ-90 from GLONASS
satellite observations and then comparing these coordinates to known coordinates in WGS84, requires at
least four satellites visible at a station to calculate station coordinates. This approach of direct estimation
of the transformation parameters may work with as little as two observations per site. Depending on the
number of stations involved, at some sites only one observation may mathematically be sufficient to get
a solution, but this one measurement contributes only to the station clock error.
The geometrical range Eq. (6.3.15) and with it the observation equation Eq. (6.3.1) is non-linear in
the unknown transformation parameters. Before trying to solve a system of observation equations, the
Eq. (6.3.15) has to be linearized. Therefore, this equation is expanded in a Taylor series around the
approximate values ∆x
0
, ∆y
0
, ∆z
0
, δs
0
, δε
0
, δψ
0
, δω
0
:
S
R
(∆x, ∆y, ∆z, δs, δε, δψ, δω) =
S
R,0
+
∂
S
R
∂∆x
0
· (∆x − ∆x
0
) +
∂
S
R
∂∆y
0
· (∆y − ∆y
0
)+
∂
S
R
∂∆z
0
· (∆z − ∆z
0
) +
∂
S
R
∂δs
0
· (δs − δs
0
) +
∂
S
R
∂δε
0
· (δε − δε
0
) +
(6.3.16)
∂
S
R
∂δψ
0
· (δψ − δψ
0
) +
∂
S
R
∂δω
0
· (δω − δω
0
)
with
S
R,0
=
S
R
(∆x
0
, ∆y
0
, ∆z
0
, δs
0
, δε
0
, δψ
0
, δω
0
)
∂
S
R
∂d
0
=
∂
S
R
∂d
(∆x
0
, ∆y
0
, ∆z
0
, δs
0
, δε
0
, δψ
0
, δω
0
),
d ∈ {∆x, ∆y, ∆z, δs, δε, δψ, δω}
For reasons of simplicity, Eq. (6.3.15) is not differentiated directly. Rather it is preferred to differen-
tiate Eq. (6.3.14) and then insert the results into the differentiated Eq. (6.3.12).
Differentiating Eq. (6.3.12) by any of the unknown transformation parameters d, d ∈ {∆x, ∆y, ∆z,
δs, δε, δψ, δω} yields:
∂
S
R
∂d
0
= −
ω
E
c
· x
R,W GS(t
RX
)
·
∂y
S
W GS(t
T X
)
∂d
0
− y
R,W GS(t
RX
)
·
∂x
S
W GS(t
T X
)
∂d
0
+
1
S
0
·
ω
E
c
2
s
0
· x
R,W GS(t
RX
)
·
∂y
S
W GS(t
T X
)
∂d
0
− y
R,W GS(t
RX
)
·
∂x
S
W GS(t
T X
)
∂d
0
−
(6.3.17)
ξ
0
·
∂x
S
W GS(t
T X
)
∂d
0
− υ
0
·
∂y
S
W GS(t
T X
)
∂d
0
− ζ
0
·
∂z
S
W GS(t
T X
)
∂d
0
with
S
0
=
ω
E
c
2
· s
2
0
+ ξ
2
0
+ υ
2
0
+ ζ
2
0
s
0
= x
R,W GS(t
RX
)
· y
S
W GS(t
T X
),0
− y
R,W GS(t
RX
)
· x
S
W GS(t
T X
),0
= x
R,W GS(t
RX
)
· ∆y
0
+ (1 + δs
0
) · −δω
0
· x
S
P Z(t
T X
)
+ y
S
P Z(t
T X
)
+ δε
0
· z
S
P Z(t
T X
)
−
y
R,W GS(t
RX
)
· ∆x
0
+ (1 + δs
0
) · x
S
P Z(t
T X
)
+ δω
0
· y
S
P Z(t
T X
)
− δψ
0
· z
S
P Z(t
T X
)
ξ
0
= x
R,W GS(t
RX
)
− x
S
W GS(t
T X
),0
= x
R,W GS(t
RX
)
− ∆x
0
− (1 + δs
0
) · x
S
P Z(t
T X
)
+ δω
0
· y
S
P Z(t
T X
)
− δψ
0
· z
S
P Z(t
T X
)
υ
0
= y
R,W GS(t
RX
)
− y
S
W GS(t
T X
),0
= y
R,W GS(t
RX
)
− ∆y
0
− (1 + δs
0
) · −δω
0
· x
S
P Z(t
T X
)
+ y
S
P Z(t
T X
)
+ δε
0
· z
S
P Z(t
T X
)
ζ
0
= z
R,W GS(t
RX
)
− z
S
W GS(t
T X
),0
= z
R,W GS(t
RX
)
− ∆z
0
− (1 + δs
0
) · δψ
0
· x
S
P Z(t
T X
)
− δε
0
· y
S
P Z(t
T X
)
+ z
S
P Z(t
T X
)

6.3 Direct Estimation of Transformation Parameters
67
x
S
W GS(t
T X
),0
= x
S
W GS(t
T X
)
(∆x
0
, ∆y
0
, ∆z
0
, δs
0
, δε
0
, δψ
0
, δω
0
)
y
S
W GS(t
T X
),0
= y
S
W GS(t
T X
)
(∆x
0
, ∆y
0
, ∆z
0
, δs
0
, δε
0
, δψ
0
, δω
0
)
z
S
W GS(t
T X
),0
= z
S
W GS(t
T X
)
(∆x
0
, ∆y
0
, ∆z
0
, δs
0
, δε
0
, δψ
0
, δω
0
)
∂x
S
W GS(t
T X
)
∂d
0
=
∂x
S
W GS(t
T X
)
∂d
(∆x
0
, ∆y
0
, ∆z
0
, δs
0
, δε
0
, δψ
0
, δω
0
)
∂y
S
W GS(t
T X
)
∂d
0
=
∂y
S
W GS(t
T X
)
∂d
(∆x
0
, ∆y
0
, ∆z
0
, δs
0
, δε
0
, δψ
0
, δω
0
)
∂z
S
W GS(t
T X
)
∂d
0
=
∂z
S
W GS(t
T X
)
∂d
(∆x
0
, ∆y
0
, ∆z
0
, δs
0
, δε
0
, δψ
0
, δω
0
)
From Eq. (6.3.14) one obtains the following partial derivatives:
∂x
S
W GS(t
T X
)
∂∆x
(∆x
0
, ∆y
0
, ∆z
0
, δs
0
, δε
0
, δψ
0
, δω
0
) = 1
∂x
S
W GS(t
T X
)
∂∆y
(∆x
0
, ∆y
0
, ∆z
0
, δs
0
, δε
0
, δψ
0
, δω
0
) = 0
∂x
S
W GS(t
T X
)
∂∆z
(∆x
0
, ∆y
0
, ∆z
0
, δs
0
, δε
0
, δψ
0
, δω
0
) = 0
(6.3.18)
∂x
S
W GS(t
T X
)
∂δs
(∆x
0
, ∆y
0
, ∆z
0
, δs
0
, δε
0
, δψ
0
, δω
0
) = x
S
P Z(t
T X
)
+ δω
0
· y
S
P Z(t
T X
)
− δψ
0
· z
S
P Z(t
T X
)
∂x
S
W GS(t
T X
)
∂δε
(∆x
0
, ∆y
0
, ∆z
0
, δs
0
, δε
0
, δψ
0
, δω
0
) = 0
∂x
S
W GS(t
T X
)
∂δψ
(∆x
0
, ∆y
0
, ∆z
0
, δs
0
, δε
0
, δψ
0
, δω
0
) = −(1 + δs
0
) · z
S
P Z(t
T X
)
∂x
S
W GS(t
T X
)
∂δω
(∆x
0
, ∆y
0
, ∆z
0
, δs
0
, δε
0
, δψ
0
, δω
0
) = (1 + δs
0
) · y
S
P Z(t
T X
)
∂y
S
W GS(t
T X
)
∂∆x
(∆x
0
, ∆y
0
, ∆z
0
, δs
0
, δε
0
, δψ
0
, δω
0
) = 0
∂y
S
W GS(t
T X
)
∂∆y
(∆x
0
, ∆y
0
, ∆z