- Le calcul de position GPS consiste à résoudre simultanément les pseudo-distances (pseudorange) jusqu’aux satellites, les éphémérides et l’erreur d’horloge du récepteur ; l’article traite des données brutes dans Matlab jusqu’à l’implémentation d’une estimation de position réelle
- Pour les calculs, le repère WGS 84 ECEF est plus adapté que la latitude et la longitude, et le repère local ENU est utilisé en complément pour l’azimut et l’angle d’élévation vus par l’utilisateur, ainsi que pour l’analyse des erreurs
- La position des satellites est obtenue avec la procédure de la GPS Interface Specification et les paramètres orbitaux ; comme la Terre tourne pendant la propagation du signal, les coordonnées au moment de l’émission doivent être corrigées pour correspondre au repère ECEF au moment de la réception
- La position de l’utilisateur et le biais d’horloge du récepteur sont estimés par moindres carrés itératifs à partir des pseudo-distances corrigées d’au moins 4 satellites ; dans l’analyse d’exemple, les délais ionosphérique et troposphérique sont exclus
- Dans une expérience avec un récepteur fixe utilisant u-blox NEO-M8T/6T, RTKLib STRSVR, RTCM 1002/1019 et goGPS, les écarts types de position étaient East 14,00 m, North 39,88 m, Up 47,35 m, et le biais d’horloge dérivait linéairement à 4.27e-7sec/sec
Le problème de base du calcul de position GPS
- La tâche centrale du GPS est de calculer la position de l’utilisateur
- La latitude, la longitude et l’altitude sont familières pour représenter une position à la surface de la Terre, mais elles sont peu pratiques pour les calculs mathématiques, car une différence d’un degré de latitude ou de longitude ne correspond pas toujours à la même distance physique
- La distance correspondant à un degré de longitude est maximale à l’équateur et devient proche de 0 aux pôles
- Les calculs nécessitent un repère cartésien dans lequel une différence unitaire de coordonnées représente une distance physique constante
- Le GPS détermine la position de l’utilisateur à partir des distances vers plusieurs satellites et de la position de ces satellites
- Il faut d’abord calculer la distance jusqu’à chaque satellite ainsi que la position de chacun
Repères : ECEF, WGS 84, ENU
- Le repère cartésien fixé à la Terre et tournant avec elle est appelé ECEF (Earth Centered, Earth Fixed)
- Il convient à la représentation de la position utilisateur, car les coordonnées d’un utilisateur immobile à la surface restent constantes dans le temps
- Le repère ECEF le plus courant est le WGS 84, développé par le département de la Défense des États-Unis
- Son origine est le centre de masse de la Terre
- L’axe z passe par le CTP (Conventional Terrestrial Pole)
- Le CTP est la position moyenne du pôle terrestre entre 1900 et 1905 ; la position réelle du pôle se déplace dans un cercle d’environ 15 m de rayon
- L’axe x passe par l’intersection du plan équatorial du CTP et du méridien de référence, le Mean Greenwich Meridian
- Il est naturel de traiter le mouvement des satellites dans un repère inertiel, conformément à la mécanique newtonienne, mais la GPS Interface Specification fournit une procédure pour calculer la position d’un satellite à un instant donné dans le repère ECEF
- Pour les applications locales, le repère ENU (East-North-Up), centré sur la position de l’utilisateur, est pratique
- Les coordonnées ECEF peuvent être converties en coordonnées ENU par multiplication matricielle à partir de la latitude et de la longitude de l’utilisateur
- La transformation ENU est utilisée pour calculer l’azimut et l’angle d’élévation des satellites
Définition de l’altitude : ellipsoïde de référence et géoïde
- Pour l’altitude, il faut d’abord définir « par rapport à quoi » elle est mesurée
- L’ellipsoïde de référence (reference ellipsoid) est un modèle qui abstrait la Terre sous forme d’ellipsoïde aplati
- Il est centré sur la Terre, et son axe de rotation coïncide avec l’axe z ECEF
- Le rayon de 6 371 km couramment utilisé lorsque l’on modélise la Terre comme une sphère se situe entre le demi-grand axe et le demi-petit axe
- Un point réel de la surface terrestre peut se trouver au-dessus ou au-dessous de l’ellipsoïde de référence
- Le géoïde (geoid) est l’ensemble des points ayant la même valeur de potentiel gravitationnel ; c’est une surface ayant une signification physique
- L’altitude par rapport au géoïde est appelée altitude orthométrique ou hauteur au-dessus du niveau moyen de la mer (MSL)
- Le géoïde est généralement spécifié comme un ensemble de valeurs de hauteur au-dessus de l’ellipsoïde de référence
- La latitude, la longitude et l’altitude sont définies comme des coordonnées ellipsoïdales
- La latitude géodésique est l’angle entre le plan équatorial et la droite perpendiculaire à la surface de l’ellipsoïde au point P
- L’angle de la droite reliant le centre de la Terre au point P est la latitude géocentrique ; si la Terre était une sphère parfaite, elle coïnciderait avec la latitude géodésique
- La conversion des coordonnées ellipsoïdales vers les coordonnées cartésiennes peut se faire en une étape, mais la conversion d’ECEF vers les coordonnées ellipsoïdales nécessite une procédure itérative à convergence rapide
Calcul de la position des satellites
- L’orbite idéale d’un satellite est une orbite elliptique décrite par 6 éléments orbitaux képlériens
- 5 éléments déterminent la taille et la forme de l’ellipse ainsi que l’orientation du plan orbital
- Le 6e élément détermine la position du satellite à une epoch donnée
- Les orbites réelles des satellites ne sont pas des ellipses parfaites, en raison des inhomogénéités de la composition terrestre et des effets gravitationnels du Soleil et de la Lune
- Le GPS diffuse 16 paramètres orbitaux pour corriger ces perturbations
- Le tableau 20-IV de la GPS Interface Specification contient la procédure de calcul de la position des satellites, corrections orbitales incluses
- La position de l’utilisateur est calculée à l’instant de réception t, mais le signal GPS est parti du satellite à l’instant t-τ
- La position du satellite est calculée à l’instant d’émission du signal, t-τ
- Comme la Terre tourne pendant la propagation du signal, durant τ, le vecteur de position du satellite est tourné de la quantité de rotation terrestre correspondante afin de l’aligner avec le repère ECEF utilisateur à l’instant de réception t
- Ce n’est pas équivalent à simplement calculer la position du satellite à l’instant t
Pseudo-distance et biais d’horloge
- Un récepteur GPS compare l’horodatage d’émission inclus dans le signal du satellite avec l’heure du récepteur, puis multiplie la différence par la vitesse de la lumière pour calculer approximativement la distance au satellite
- Cette mesure est la pseudo-distance (pseudorange)
- Si l’horloge du satellite et celle du récepteur étaient parfaitement synchronisées, et si le signal se propageait en ligne droite dans le vide à la vitesse de la lumière, elle serait égale à la distance réelle
- En pratique, elle diffère de la distance réelle à cause des décalages d’horloge et des délais atmosphériques
- Le biais d’horloge du satellite peut entraîner des erreurs de position de plusieurs milliers de mètres et doit impérativement être corrigé
- Il est calculé avec un polynôme utilisant les coefficients du message d’éphéméride GPS et avec un terme relativiste
- Le polynôme fournit l’essentiel de la correction, tandis que l’effet relativiste contribue à hauteur d’environ 1 à 10 m selon la position du satellite
- Le biais d’horloge du récepteur est une inconnue qui doit être estimée avec la position de l’utilisateur
- Dans l’algorithme, le biais d’horloge est traité en unité de distance après multiplication par la vitesse de la lumière
- Le délai atmosphérique se divise en composantes ionosphérique et troposphérique
- Le délai ionosphérique provoque généralement une erreur de position d’environ 25 m
- Le délai troposphérique provoque généralement une erreur de position d’environ 2 m
- L’analyse expérimentale de l’article ignore ces délais
Algorithme d’estimation de la position utilisateur et du biais d’horloge
- Les mesures de pseudo-distance corrigées s’expriment comme la somme de la distance réelle utilisateur-satellite, du biais d’horloge du récepteur et des erreurs non modélisées
- La position utilisateur et le biais d’horloge sont recherchés comme les valeurs minimisant l’écart entre les pseudo-distances mesurées et prédites
- La solution est une procédure de moindres carrés itératifs
- La valeur initiale de la position utilisateur est
[0 0 0] - La valeur initiale du biais d’horloge utilisateur est
0 - À chaque itération, la matrice G est construite en empilant les vecteurs unitaires de direction des satellites depuis l’estimation courante de la position
- On résout les corrections de position et de biais d’horloge, puis on itère jusqu’à ce que la variation soit inférieure au seuil
- La valeur initiale de la position utilisateur est
- S’il y a exactement 4 satellites et une configuration non dégénérée, une solution directe peut être obtenue
- Lorsque le ciel n’est pas masqué, davantage de satellites sont visibles, et on utilise généralement une solution par moindres carrés
- L’implémentation suit le flux suivant
- Prendre en entrée les pseudo-distances brutes et les éphémérides des satellites
- Calculer le biais d’horloge de chaque satellite et corriger les pseudo-distances
- Appliquer les corrections ionosphériques et troposphériques lorsque c’est possible
- Corriger les pseudo-distances avec le biais d’horloge courant du récepteur
- Diviser la pseudo-distance par la vitesse de la lumière pour obtenir le temps de propagation du signal τ
- Calculer la position du satellite à l’instant t-τ
- Tenir compte de la rotation de la Terre pendant τ pour aligner la position du satellite avec le repère ECEF utilisateur
- Construire la matrice G et les écarts de pseudo-distance, puis calculer les corrections de position et de biais d’horloge
Détails d’implémentation dans Matlab
- La plupart du code Matlab consiste à évaluer en une seule fois les inconnues de gauche à partir des valeurs connues de droite
- Certains calculs n’ont pas de solution analytique fermée et nécessitent un solver
- Un exemple est l’étape, dans le calcul de la position du satellite, qui consiste à obtenir l’anomalie excentrique E à partir de l’anomalie moyenne M
- La relation
E - e*sin(E) == Mne peut pas être résolue sous forme fermée ;vpasolveest donc utilisé
- Le code en annexe comprend les fonctions suivantes
- Calcul de la position utilisateur et du biais d’horloge
- Calcul de la position des satellites
- Calcul de la solution par moindres carrés pour la position utilisateur et le biais d’horloge
- Calcul du biais d’horloge des satellites
- Conversion des coordonnées ECEF WGS84 en coordonnées ellipsoïdales
- Conversion du format des données d’éphémérides
Configuration de collecte de données GPS brutes
- Pour obtenir des données GPS brutes, il faut un récepteur qui produise des informations de timing telles que les pseudo-distances brutes et les éphémérides des satellites, et non un appareil GPS ordinaire qui calcule seulement une position en interne avant de l’exporter
- Les puces u-blox NEO-M8T et 6T répondent à ces exigences
- Un assemblage matériel comprenant une unité GPS, une antenne et un port de sortie série peut être acheté sur Amazon pour environ 40 dollars
- L’utilitaire STRSVR de RTKLib est utilisé pour recevoir et enregistrer les signaux GPS bruts
- RTKLib est un ensemble de programmes open source prenant en charge le positionnement standard et précis pour des GNSS comme GPS, Glonass, Galileo ou Baidu
- STRSVR convertit la sortie au format personnalisé du récepteur u-blox vers le format standard RTCM
- Les informations nécessaires se trouvent dans les messages RTCM 1002 et 1019
- 1002 contient les informations de pseudo-distance brute
- 1019 contient les informations d’éphémérides des satellites
- STRSVR est configuré pour recevoir les données via le port série à 9600 Baud et les enregistrer dans un fichier au format RTCM 3
- La collecte des données est effectuée sur le toit d’un immeuble d’habitation
- Le récepteur GPS est placé à un endroit où le ciel n’est pas masqué
- Le logiciel u-blox u-center est utilisé pour vérifier qu’un nombre suffisant de satellites est visible et qu’un bon position fix est possible
- Environ 1 heure de données GPS brutes est collectée
Traitement RTCM et utilisation de goGPS
- STRSVR enregistre les données GPS brutes au format binaire RTCM3
- Pour les traiter dans Matlab, les données RTCM3 doivent être décodées en structures de données Matlab
- Au lieu d’écrire directement un décodeur RTCM, la fonction
load_streamde la bibliothèque Matlab goGPS est utilisée- Elle lit les fichiers au format RTCM et extrait les messages RTCM
- Les données extraites sont enregistrées dans un fichier
.matet utilisées comme entrée de l’algorithme de calcul de position
- Le fichier
rtcm_dataest également fourni- En raison des restrictions de sécurité de WordPress, il est fourni avec l’extension
.txtau lieu de.mat - Après téléchargement, il faut le renommer en
.mat
- En raison des restrictions de sécurité de WordPress, il est fourni avec l’extension
Résultats expérimentaux : variation de position et dérive d’horloge
- Le récepteur étant immobile pendant la collecte, la variation temporelle de la position calculée montre les performances réelles de l’algorithme de calcul de position
- Dans le repère ENU centré sur l’utilisateur, les écarts types des composantes de position sont les suivants
- East : 14,00 m
- North : 39,88 m
- Up : 47,35 m
- La variation de position est d’environ 30 m dans les directions East et North, et plus importante dans la direction Up
- Le biais d’horloge du récepteur n’est pas constant : il dérive linéairement dans le temps
- Dans l’algorithme, le biais d’horloge est traité en unité de distance
- Dans le graphique de résultat, il est converti en unité de temps en le divisant par la vitesse de la lumière
- La quantité de dérive est de 4.27e-7sec/sec
Calcul de l’azimut et de l’angle d’élévation des satellites
- L’azimut et l’angle d’élévation des satellites étant définis du point de vue de l’utilisateur, ils sont calculés dans le repère ENU centré sur l’utilisateur
- La procédure de calcul est la suivante
- Calculer, dans le repère ECEF, le vecteur de position allant de l’utilisateur vers le satellite
- Convertir la position de l’utilisateur en coordonnées ellipsoïdales, c’est-à-dire latitude et longitude
- Faire tourner ce vecteur de position vers le repère ENU centré sur l’utilisateur
- Calculer l’azimut et l’angle d’élévation à partir des coordonnées ENU
- Pour l’epoch d’exemple, les angles d’élévation calculés des 8 satellites sont tous positifs
- Les azimuts peuvent être aussi bien positifs que négatifs
- Comme l’utilisateur ne peut pas voir les satellites sous l’horizon, il est naturel que les angles d’élévation soient positifs
- En calculant la position des satellites avec la même procédure sur plusieurs epochs, on peut produire un satellite track chart comme ceux affichés par les logiciels de traitement GPS
DOP : facteur géométrique de la qualité d’estimation de position
- Le DOP (Dilution of Precision) est un indicateur évaluant la qualité d’une estimation de position
- L’erreur de position dépend non seulement du bruit de mesure, mais aussi de la géométrie utilisateur-satellites
- Plus les mesures de pseudo-distance et de position des satellites sont bruitées, plus l’erreur de position augmente
- Plus les satellites sont largement répartis en azimut et en élévation, plus la géométrie est favorable et plus le DOP est faible
- La covariance des erreurs de position et de biais d’horloge se décompose comme une fonction de l’erreur de distance utilisateur et de la matrice G
- La matrice G est constituée des vecteurs unitaires allant de l’utilisateur vers les satellites
- La matrice G dans le repère ECEF est tournée vers le repère ENU afin de faciliter le calcul du DOP
- Les composantes du DOP sont définies dans les directions East, North et Up
- Le HDOP est le DOP horizontal combinant les composantes East et North
- Le VDOP est le DOP vertical de la composante Up
- Dans les données réelles, le HDOP et le VDOP sont généralement inférieurs à 2,5
- Cette valeur est considérée comme suffisante
- Le VDOP est supérieur au HDOP
- Le VDOP augmente parce qu’un utilisateur à la surface ne peut pas observer les satellites sous l’horizon, et que les signaux des satellites situés à moins de 10 degrés au-dessus de l’horizon sont généralement trop bruités pour être utilisés
L’échelle de l’infrastructure GPS
- La mise en place de la constellation GPS a coûté environ 30 milliards de dollars, et le gouvernement américain dépense environ 1 milliard de dollars par an pour sa maintenance
- La valeur d’Uber, qui n’aurait pas pu exister sans le GPS, est mentionnée comme dépassant 70 milliards de dollars
- Si l’on inclut les nombreuses applications rendues possibles par le GPS, l’investissement public dans le GPS peut être vu comme un cas ayant généré des retombées économiques et technologiques très importantes
1 commentaires
Commentaires sur Hacker News
Android donne accès depuis assez longtemps déjà à l’approche par phase de porteuse (carrier phase), qui permet d’obtenir la position relative d’appareils situés à peu près dans le même quartier avec une précision telle qu’il faut se demander où l’antenne GNSS est cachée dans l’appareil.
À lui seul, ce n’est pas extraordinaire, mais cela devient intéressant si l’on y combine les accéléromètres et gyroscopes de chaque appareil.
Comme les variations des mesures de pseudodistance GNSS sont prévisibles même lorsque l’on n’est pas à l’arrêt, cela peut fonctionner en temps réel avec une faible dégradation des performances.
Par exemple, on peut faire atterrir automatiquement un avion radiocommandé sans roues dans la benne d’un camion, tout en évitant les rayures ou la dépendance à une piste en herbe.
Si la consommation électrique n’est pas un critère absolument critique, il n’y a pas besoin de rendre un assez bon récepteur GNSS coûteux ; je ne comprends pas pourquoi on ne peut pas simplement en acheter une paire pour 100 dollars.
Voici un récepteur GNSS à 30 dollars. Il est livré avec quelques babioles inutiles en plus.
Pour ceux qui veulent fabriquer eux-mêmes un récepteur GPS, il existe un projet entièrement open source qui explique aussi beaucoup de théorie : http://www.aholme.co.uk/GPS/Main.htm
https://lea.hamradio.si/~s53mv/navsats/theory.html
Un article qui revient souvent dans les commentaires sur le GPS, mais à juste titre : https://ciechanow.ski/gps/
Une autre explication, probablement plus interactive :
https://ciechanow.ski/gps/
Il existe aussi une autre bonne implémentation open source :
https://m.youtube.com/watch?v=dVD1Yws__v0
J’ai vu des chercheurs collecter des données GPS d’animaux aquatiques qui ne remontent à la surface qu’occasionnellement, et seulement très brièvement.
En enregistrant les données brutes puis en les post-traitant, on réduit fortement à la fois la consommation électrique et la durée minimale d’exposition nécessaire au signal des satellites, jusqu’à moins d’une seconde.
En lisant le passage « la figure ci-dessous montre comment la géométrie utilisateur-source influe sur l’incertitude de la position de l’utilisateur », je me suis dit qu’il serait bien d’avoir, dans les applis de cartographie sur téléphone, un réglage remplaçant le cercle d’incertitude de position par ce genre de forme d’intersection d’arcs.
J’ai entendu dire que le GPS était l’une des rares applications du quotidien où il faut tenir compte des effets relativistes. Est-ce que cela veut dire que ces effets relativistes ont déjà été retirés des données générées ?
Si tu parles de la sortie d’un appareil GPS grand public, alors oui. Le traitement post-acquisition qui produit la sortie corrige toutes sortes d’effets sources d’erreurs.
Cet article traite des données GPS brutes diffusées par plusieurs satellites ; pour en tirer des valeurs de sortie, il faut les traiter, et l’on ajoute souvent des entrées supplémentaires, comme des corrections provenant de stations au sol ou maritimes, afin d’améliorer la précision.
Les différents fournisseurs d’équipements GPS font globalement des choses similaires, mais les détails sont essentiels.
Le lien https://ciechanow.ski/gps/ cité dans d’autres commentaires est aussi une bonne lecture.
Le temps s’écoule aussi plus lentement pour les objets qui se déplacent plus vite que l’observateur, et les satellites se déplacent assez vite.
Le GPS exige que le temps soit synchronisé entre l’observateur et les satellites ; l’écoulement du temps est donc corrigé pour tenir compte des effets de la relativité restreinte et de la relativité générale.
La relativité est donc importante, mais il n’est pas nécessaire d’en savoir beaucoup sur elle pour résoudre sa propre position.
Cela dit, plusieurs formes de RTK à longue base ou en réseau peuvent nécessiter une modélisation plus sophistiquée.
L’étape suivante, c’est le PPP ou le RTK. Le GNSS est un terrier de possibilités très amusant.
Exercice pour les platistes : expliquez comment les cartes GPS de votre téléphone fonctionnent sans satellites orbitant autour d’une Terre sphérique. Montrez votre raisonnement.
Il faut d’abord établir un raisonnement assez complexe montrant qu’un téléphone ne peut pas faire cela sans satellites.
Un exemple un peu plus simple, et plus difficile à contourner pour un platiste, est que l’ISS est presque visible à l’œil nu et clairement visible avec un télescope de jardin. Les satellites Starlink aussi, de façon similaire.
Le platisme n’est pas une position à laquelle on parvient par la raison ; il naît presque toujours de la confusion ou comme conséquence inévitable d’une croyance centrale inébranlable. Le plus souvent, il vient d’une lecture extrêmement littérale de la Bible, ou d’un délire paranoïaque selon lequel « tout ce qui est officiel est faux ».