An Autonomous Reference Frame for Relativistic GNSS

0.67 MB·

by  Uroš Kostić

ANAUTONOMOUSREFERENCEFRAMEFORRELATIVISTICGNSS
UrosˇKostic´1,2,MartinHorvat1,SanteCarloni3,PacoˆmeDelva4,andAndrejaGomboc1,2
1FacultyofMathematicsandPhysics,UniversityofLjubljana,Jadranska19,1000Ljubljana,Slovenia
2Sloveniancentreofexcellenceforspacesciencesandtechnologies,Asˇkercˇeva12,1000Ljubljana,Slovenia
3InstituteofTheoreticalPhysics,FacultyofMathematicsandPhysics,CharlesUniversity,VHolesovickach2,18000
Praha8,CzechRepublic
4LNE-SYRTE,ObservatoiredeParis,CNRSetUPMC,61avenuedelObservatoire,75014,Paris,France
ABSTRACT well.
Our previous work shows that it is possible to construct
4
1 Current GNSS systems rely on global reference frames such a system and do the positioning within it: the or-
0 which are fixed to the Earth (via the ground stations) bital parameters of the GNSS satellites can be deter-
2 so their precision and stability in time are limited by mined and checked internally by the GNSS system it-
ourknowledgeoftheEarthdynamics. Thesedrawbacks self through inter-satellite links (Cˇadezˇ et al. 2011). In
n
a could be avoided by giving to the constellation of satel- thisway, wecanconstructanAutonomousBasisofCo-
J lites the possibility of constituting by itself a primary ordinates(ABC)whichisindependentofanyEarthbased
4 and autonomous positioning system, without any a pri- coordinatesystem. Thissystemconstructsitselfaseach
2 orirealizationofaterrestrialreferenceframe. Ourwork satellite receives proper times of all other satellites, and
shows that it is possible to construct such a system, an by comparison with its own proper time determines the

] Autonomous Basis of Coordinates, via emission coordi- parametersoftheABCsystemwithgreataccuracy.
c
nates. HerewepresenttheideaoftheAutonomousBasis
q
of Coordinates and its implementation in the perturbed Here we present the results of our recent work, where
-
r space-timeofEarth,wherethemotionofsatellites,light wehavefurtherdevelopedrelativisticpositioningandthe
g
propagation,andgravitationalperturbationsaretreatedin ABCbyincludingallrelevantgravitationalperturbations,
[
theformalismofgeneralrelativity. suchasEarthmultipoles(uptothe6th), Earthsolidand

2 ocean tides, the Sun, the Moon, Jupiter, Venus, and the
v Kerreffect.
Keywords: RelativisticGNSS.
6
0
7
2. MOTIONOFSATELLITESINAPERTURBED
4
1. INTRODUCTION
SPACE-TIME
.
1
0
4 TheclassicalconceptofpositioningsystemforaGlobal The perturbed satellite orbits were calculated from the
1 navigation satellite system (GNSS) would work ideally perturbedHamiltonian(Goldstein1980)
: ifallsatellitesandthereceiverwereatrestinaninertial
v
i referenceframe. Butatthelevelofprecisionneededbya H = 1gµνp p , (1)
X GNSS,onehastoconsidercurvatureandrelativisticiner- 2 µ ν
r tialeffectsofspacetime,whicharefarfrombeingnegligi-
a ble.Theseeffectsaremosteasilyandelegantlydealtwith where the perturbed metric gµν is a sum of the unper-
inarelativisticpositioningsystembasedonemissionco- turbed(i.e. Schwarzschild)metricgµ(0ν) andtheperturba-
ordinates (Coll & Morales 1991; Rovelli 2002; Blago- tivemetrichµν
jevic´ et al. 2002; Cˇadezˇ et al. 2010; Delva et al. 2011).
Theydependonthesetoffoursatellitesandtheirdynam- gµν =gµ(0ν)+hµν , |hµν|(cid:28)gµ(0ν). (2)
ics, and can be linked to a terrestrial reference system.
The Hamiltonian (Eq.1) is rewritten with the canonical
Consequently,thedifficultynolongerliesintheconcep-
momentap =g x˙ν as
tion of the primary reference frame but in its link with µ µν
terrestrialreferenceframes(Cˇadezˇ etal.2010). Thisal- 1 1
H = g(0)µνp p − hµνp p =H(0)+∆H , (3)
lowstocontrolmuchmorepreciselyalltheperturbations 2 µ ν 2 µ ν
thatlimittheaccuracyandthestabilityoftheprimaryref-
erenceframe,ifthedynamicsoftheGNSSsatellites,de- whereH(0) and∆H aretheunperturbedandthepertur-
scribedbytheirorbitalparameters,isknownsufficiently bativepartoftheHamiltonian,respectively.
Theunperturbedorbitisthengivenbytheequations
10000 400
∂H(0) ∂H(0) 300
x˙µ = ∂p(0) and p˙(µ0) =− ∂xµ , (4) 8000 200
µ
m 6000(cid:68) 100
where the coordinates and the momenta k (cid:64) 0
for Schwarzschild case are (xµ,pµ) = (cid:68)L 4000 0 2 4 6 8 10
(t,r,θ,φ,p ,p ,p ,p ). The unperturbed Hamilto-
t r θ ϕ
nian H(0) admits 6 constants of motion (Qk,P ): the 2000
k
orbital energy P = E, the magnitude of the angular
1
momentumP = l,thez-componentoftheangularmo- 0
2
mentum P = l , the time of the first apoapsis passage 0 50 100 150 200 250 300 350
3 z
Q1 = −t , the longitude of the apoapsis Q2 = ω, and t day
a
the longitude of the ascending node Q3 = Ω. These
constantsareusedtofindtheanalyticalsolutionsforthe Figure1. Thedifferencesintheposition∆Lofthesatel-
orbitsintheform: lite due to sum of all perturb(cid:64)atio(cid:68)ns. The gray graph
shows the long-term changes in one year, while the red
t=t(λ|Qk,Pk) r =r(λ|Qk,Pk) inset shows the short time-scale changes within the first
(5)
θ =θ(λ|Qk,P ) φ=φ(λ|Qk,P ), 10days. The(Schwarzschild)timeonx-axiscountsdays
k k
from1January2012at12:00noon.Axesonredandgray
whereλisthetrueanomaly(Kostic´ 2012). plotshavethesameunits. Theinitialvaluesofparame-
ters are: t = 12 h, ω = 0◦, Ω = 0◦, a = 29602 km,
a
If the metric is perturbed, (Qk,P ) are no longer con- ε=0.007,ι=56◦,correspondingtotheGalileosystem.
k
stants of motion – they are slowly changing functions
of time. Labelling these new canonical variables as
(Qˆk,Pˆ ),weobtaintheirtimeevolutionfromthefollow- oneyear. Adetailedstudyofeffectsofeachperturbation
k
ingexpressions(Goldstein1980): oneachorbitalparametershowsthatthesechangesarise
mostlyfromtheprecessionsofΩandω;thechangesdue
Qˆ˙k (cid:39) ∂∆H(cid:12)(cid:12)(cid:12) =−1∂hµνpµpν (6) to other orbital parameters are 2-3 orders of magnitude
∂P (cid:12) 2 ∂P smallerandoscillating.
k Qk,Pk k
Pˆ˙k (cid:39)− ∂∂∆QHk (cid:12)(cid:12)(cid:12)(cid:12) = 12∂h∂µνQpkµpν . (7)
Qk,Pk 3. RELATIVISTICPOSITIONINGSYSTEM
TheEqs.6and7arenumericallyintegratedtoobtainthe
solutionsforQˆk(λ|Qk,P )andPˆ (λ|Qk,P ),whichare
k k k We simulate a constellation of four satellites moving
then used to replace (Qk,P ) in the analytical expres-
k alongtheirtime-likegeodesics(Cˇadezˇ etal.2010;Delva
sionsforunperturbedorbits(Eq.5). Inthisway,theper-
et al. 2011; Kostic´ 2012). The initial orbital parameters
turbedorbitisdescribedastime-evolvingunperturbedor-
ofthegeodesicsareknownandtheirevolutionduetoper-
bit:
turbationsiscalculatedasshowninSection2.
t=t(λ|Qˆk(λ),Pˆ (λ)) r =r(λ|Qˆk(λ),Pˆ (λ))
k k
Ateverytime-stepofthesimulation,eachsatelliteemits
θ =θ(λ|Qˆk(λ),Pˆk(λ)) φ=φ(λ|Qˆk(λ),Pˆk(λ)). a signal and a user on Earth receives signals from all
(8) satellites – the signals are the proper times of satellites
attheiremissioneventsandconstitutetheemissioncoor-
dinates of the user. The emission coordinates determine
Before the Eqs.6 and 7 are solved numerically, we cal-
the user’s “position” in this particular relativistic refer-
culatethemetrich forKerreffect, Jupiter, Venus, the
µν ence frame defined by the four satellites and allow him
Moon, the Sun, Earth multipoles, solid, and ocean tides
to calculate his position and time in the more custom-
using relativistic multipole expansion (Gomboc et al.
ary Schwarzschild coordinates. So, to simulate the po-
2013;Horvatetal.2014).
sitioningsystem,twomainalgorithmshavetobeimple-
mented: (1) determination of the emission coordinates,
We calculated the evolution of orbital parameters (with
and(2)calculationoftheSchwarzschildcoordinates.
initialvaluest =12h,ω =0◦,Ω=0◦,a=29602km,
a
ε=0.007,ι=56◦,correspondingtotheGalileosystem)
foreachperturbationindividuallyaswellasforthesum
of all perturbations and thus obtained the corresponding Determinationoftheemissioncoordinates Thesatel-
orbits.Becausetheperturbedorbitsdonotdifferfromthe lites’trajectoriesareparametrizedbytheirtrueanomaly
unperturbedonessignificantly,inFig.1weshowonlythe λ. The event P = (t ,x ,y ,z ) marks user’s
o o o o o
differencesbetweentheperturbedandunperturbedposi- Schwarzschildcoordinatesatthemomentofreceptionof
tionsforallperturbationsincluded: thedifferenceinpo- the signals from four satellites. Each satellite emitted a
sitions∆Lare∼300kmin10daysand∼10000kmin signal at event P = (t ,x ,y ,z ), corresponding to λ
i i i i i i
(i=1,...,4). EmissioncoordinatesoftheuseratP are,
o Table 1. Orbital parameters for 4 satellites: longitude
therefore, the proper times τ (λ ) of the satellites at P .
i i i ofascendingnode(Ω),longitudeofperigee(ω),inclina-
Taking into account that the events P and P are con-
o i tion(ι),majorsemi-axis(a),eccentricity(ε),andtimeof
nectedwithalight-likegeodesic,1 wecalculateλ atthe
i apogeepassage(t ).
emissionpointP usingtheequation a
i
t −t (λ |Qˆk(λ ),Pˆ (λ ))= # Ω[◦] ω[◦] ι[◦] a[km] ε ta[s]
o i i i k i
(9) 1 0 270 45 30000 0.007 0
T (R(cid:126) (λ |Qˆk(λ ),Pˆ (λ )),R(cid:126) ),
f i i i k i o 2 0 315 45 30000 0.007 0
where R(cid:126) = (x ,y ,z ) and R(cid:126) = (x ,y ,z ) are the 3 0 275 135 30000 0.007 0
i i i i o o o o
spatialvectorsofthesatellitesandtheuser,respectively. 4 0 320 135 30000 0.007 0
The function T calculates the time-of-flight of photons
f
betweenPo andPi asshownbyCˇadezˇ &Kostic´ (2005) 1 2
andCˇadezˇetal.(2010). TheEq.9isactuallyasystemof τ
e [k+1]
fourequationsforfourunknownλ –oncethevaluesof m
λiaredetermined,itisstraightforwiardtocalculateτifor e ti
eachsatelliteandthusobtainuser’semissioncoordinates at
n τ
atPo =(τ1,τ2,τ3,τ4). di [k+1]
r
o
o
c τ
[k]
Calculation of the Schwarzschild coordinates
Here we solve the inverse problem of calculating
τ
Schwarzschild coordinates of the event Po from proper [k]
times (τ ,τ ,τ ,τ ) sent by the four satellites. We
1 2 3 4
do this in the following way: For each satellite, we
numericallysolvetheequation
space coordinate
τ(λ |Qˆk(λ ),Pˆ (λ ))=τ , (10)
i i k i i
Figure 2. A pair of satellites exchanging their proper
toobtainλ ,whereτ(λ|Qˆk(λ),Pˆ (λ))isaknownfunc- times. At every time-step k, the satellite 1 sends the
i k
proper time of emission τ[k] to the satellite 2, which
tionforpropertimeontime-likegeodesics(Kostic´2012).
receives it at the time of reception τ[k]. The emission
The Schwarzschild coordinates of the satellites are then
andreceptioneventpairsareconnectedwithalight-like
calculated from λ using Eq.8. With the satellites’
i
geodesic.
coordinates known, we can take the geometrical ap-
proach presented by Cˇadezˇ et al. (2010) to calculate the
Schwarzschild coordinates of the user. The final step in
the true values of orbital parameters would be transmit-
thismethodrequiresusagaintosolveEq.9,however,this
tedtotheusertogetherwiththeemissioncoordinates,so
timeitistreatedasasystemof4equationsfor4unknown
to account for this in our simulations, we calculated the
usercoordinates,i.e.,solvingit,gives(t ,x ,y ,z ).
o o o o evolution of parameters from their initial values before
starting the positioning, and (2) the position of the user
Theaccuracyofthisalgorithmhasbeentestedforsatel-
iscompletelyunknown,i.e.,wedonotstartfromthelast
lites on orbits with initial parameters given in Table1
and a user at coordinates r = 6371 km, θ = 43.97◦, known position. If we did, the times for calculating the
o o
φ =14.5◦. Theuser’scoordinatesremainconstantdur- positionwouldbeevenshorter.
o
ingthesimulation. Therelativeerrors,definedas
t −te R(cid:126) −R(cid:126)e 3.1. AutonomousBasisofCoordinates
(cid:15) = o o, (cid:15) = o o , (11)
t to x,y,z R(cid:126)o
Toconstructanautonomouscoordinatesystem,weapply
are of the order 10−32 − 10−30 for coordinate t, and
theideaoftheAutonomousBasisofCoordinates(ABC)
10−28 −10−26 for x, y, and z; here teo and R(cid:126)oe are user presented in Cˇadezˇ et al. (2011) to a perturbed satellite
timeandcoordinatesascalculatedfromtheemissionco- system,i.e.,wesimulatethemotionofapairofsatellites
ordinates. Usingalaptop2forcalculations,theuser’spo- alongtheirperturbedorbits.
sition(withsucherrors)wasdeterminedin0.04s,where
weassumedthat(1)inrealapplicationsofthepositioning At each time-step of the simulation, both satellites ex-
change emission coordinates as shown in Fig.2, where,
1Thelight-likegeodesicsarecalculatedinSchwarzschildspace-time
forclarity, onlycommunicationfromsatellite1tosatel-
withoutperturbations,becausetheeffectsofperturbationsonlightprop-
lite2isplotted.Theseeventsofemissionatpropertimeτ
agationarenegligible.
2Withthefollowingconfiguration: Intel(R)Core(TM)i7-3610QM ofthefirstsatelliteandreceptionatτ ofthesecondsatel-
[email protected],8GBRAM,IntelC/C++/Fortrancompiler13.0.1. liteareconnectedwithalight-likegeodesic,i.e.,thedif-
ference between the coordinate times of emission t (τ)
1
andreceptiont2(τ)mustbeequaltothetimeofflightof 1017
aphotonbetweenthetwosatellites(cf. Eq.9)
109
T =t (τ)−t (τ). (12)
f 2 1 n
ctio 10
However,thisisonlytrueifweknowtheexactvaluesof A
the initial orbital parameters of each satellite, as well as 10(cid:45)7
theirevolution.
10(cid:45)15
When constructing the relativistic positioning system, it 10(cid:45)23
isreasonabletoassumethattheinitialorbitalparameters 0 2000 4000 6000 8000 10000
(Qˆi(0),Pˆ(0))arenotknownveryprecisely. Toimprove
i step
their values, we sum the differences between RHS and
LHS of Eq.12 for all communication events into an ac- Figure 3. The action S(Qˆi(0),Pˆ(0)) during the mini-
i
tion mizationprocess. Thefirststagetakes1765stepsandthe
S(Qˆi(0),Pˆ(0))=(cid:88)(cid:16)t (τ[k]|Qˆi(τ[k]),Pˆ(τ[k]))− secondonetakes8513steps.
i 1 i
k
t (τ[k]|Qˆi(τ[k]),Pˆ(τ[k]))− Byrepeatingtheminimizationprocedureforallpossible
2 i pairs of satellites, we can reconstruct the orbital param-
T (R(cid:126) (τ[k]|Qˆi(τ[k]),Pˆ(τ[k])), etersofeverysatelliteinthesystemwithouttrackingthe
f 1 i
(cid:17)2 satellites from the Earth and thus obtain an autonomous
R(cid:126)2(τ[k]|Qˆi(τ[k]),Pˆi(τ[k]))) , coordinatesystem.
(13)
which has a minimum value (close to zero) for the true 4. SUMMARY
initialvaluesoforbitalparameters–forthe2×6orbital
parametersthatwehave(Qˆi(0),Pˆ(0)forbothsatellites),
i
this becomes a problem of finding a minimum of a 12D In this contribution we have shown how to construct an
function. Becausetheorbitalparametersdependontime, AutonomousBasisofCoordinates(ABC)forrelativistic
theirtimeevolutionhastoberecalculated(aspresentedin GNSSinaperturbedSchwarzschildspace-time.
Section2)ateverystepoftheminimization,whichmakes
The ABC concept establishes a local inertial frame,
theminimizationprocessveryslow.
which is based solely on dynamics of GNSS satellites
The minimization was done in two stages. In the first andisthuscompletelyindependentofaterrestrialrefer-
stage, we use the PRAXIS minimization method (Brent ence. Generalrelativitywasusedtocalculatethedynam-
1973)implementedintheNLOPTlibrary(Johnson2013) ics of the GNSS satellites and the gravitational pertur-
todeterminetheparameterswithindoubleprecision.The bations affecting the dynamics (Earth multipoles, Earth
resultingvaluesarethenusedasinitialvaluesforthesec- solidandoceantides,theSun,theMoon,Jupiter,Venus,
ondstage,whereweusethesimplexmethodto“polish” andtheKerreffect),aswellastheinter-satellitecommu-
the parameters within 128-bit quad precision.3 In Fig.3 nication,whichactuallymakespossibleforconstruction
weplotthevaluesoftheactionduringtheminimization of the whole system. The numerical codes for position-
process;thefirststagetakes1765steps,whilethesecond ing determine all four coordinates in 40 ms with 25-30
onetakes8513steps. digitaccuracy,showingthatgeneralrelativistictreatment
presentsnotechnicalobstaclesforcurrentGPSdevices.
The number of time-steps along the orbits was suffi-
Because the system does not rely on Earth based refer-
cientlylarge(k = 1...433)tocoverapproximatelytwo
enceframesandisconstructedonlythroughpropertime
orbital times. The initial values of the orbital param-
exchange between satellites, it offers unprecedented ac-
eters used as starting point in the minimization differ
curacyandstability. Infact,presenttechnology,planned
fromthetruevaluesbyanamountwhichinducestheer-
to be used in the Galileo system, should be able to rou-
ror of ∼ 2 − 3 km in the satellites’ positions. At the
tinelyreachmillimetreaccuracywithrespecttoanabso-
beginning of the minimization, the value of the action
is S ≈ 1024(r /c)2, at the end of the first stage it is lute local inertial frame defined independently of Earth
g
2 × 1010(r /c)2, and at the end of the second stage it based coordinates. At this level of accuracy it seems
g necessary to decouple the local inertial frame from the
drops to 8 × 10−24(r /c)2. The relative errors of the
g geodetic Earth frame, and allow the comparison of the
orbitalparameters(Qˆi(0),Pˆ(0))aftertheminimization
i two, to tell us fine details about Earth rotation, gravita-
areoftheorderof10−22.
tionalpotentialanddynamicsofEarthcrust. Lastbutnot
least,trackingofsatelliteswithgroundstationsisneces-
3Quadprecisionisrequirediftheresultingparametersareusedin
sary only to link the relativistic positioning system to a
Eq.8,wherecancellationeffectsbecomesignificantincaseofquasi-
circularorbits. terrestrial frame, although this link can also be obtained
byplacingseveralreceiversattheknownterrestrialposi-
tions.
ACKNOWLEDGMENTS
U.K.,M.H,andA.G.acknowledgefinancialsupportfrom
ESA PECS project Relativistic Global Navigation Sys-
tem.
REFERENCES
Blagojevic´, M., Garecki, J., Hehl, F. W., & Obukhov,
Y.N.2002,Phys.Rev.D,65,044018
Brent, R. 1973, Algorithms for Minimization With-
outDerivatives,DoverBooksonMathematics(Dover
Publications)
Cˇadezˇ,A.&Kostic´,U.2005,Phys.Rev.D,72,104024
Cˇadezˇ, A., Kostic´, U., & Delva, P. 2010, Mapping the
Spacetime Metric with a Global Navigation Satellite
System, EuropeanSpaceAgency, theAdvancedCon-
cepts Team, Ariadna Final Report (09/1301), Tech.
rep.,EuropeanSpaceAgency
Cˇadezˇ, A., Kostic´, U., & Delva, P. 2011, Mapping the
Spacetime Metric with a Global Navigation Satellite
System-extensionstudy: Recoveryoforbitalconstants
using intersatellite links, European Space Agency,
the Advanced Concepts Team, Ariadna Final Report
(09/1301CCN),Tech.rep.,EuropeanSpaceAgency
Coll,B.&Morales,J.A.1991,JournalofMathematical
Physics,32,2450
Delva, P., Kostic´, U., & Cˇadezˇ, A. 2011, Advances in
SpaceResearch,47,370
Goldstein, H. 1980, Classical Mechanics, 2nd edn.
(Reading,MA:Addison-WesleyPublishingCompany)
Gomboc, A., Kostic´, U., Horvat, M., Carloni, S., &
Delva, P. 2013, in Acta Futura: Proceedings of the
workshop Relativistic Positioning Systems and their
ScientificApplications,ESA
Horvat,M.,Kostic´,U.,Carloni,S.,&Gomboc,A.2014,
submitted
Johnson, S. 2013, The NLopt nonlinear-optimization
package,http://ab-initio.mit.edu/nlopt
Kostic´, U.2012, GeneralRelativity and Gravitation, 44,
1057
Rovelli,C.2002,Phys.Rev.D,65,044017