GPS, les bases

Structure et histoire du système (quelques mots)

                     
GPS est né à partir de système militaire NAVSTAR, né en 1973. On a commencé à déployer les satellites en 1978, et tout a été pleinement opérationnel entre 1994/95. Ce système demande le contrôle et synchronisation permanente depuis des stations terrestres, stationnaires, dont la principale se trouve à Colorado Springs. Jusqu'à 1983 le secteur civil était affaibli, avec la précision plutôt d'une centaine de mètres qu'une dizaine.

[Utilisable quand même, par ex. pour les navires océaniques].

Le gouvernement Reagan a "libéré" le GPS aux applications civiles à grand échelle, quand en 1983 le vol KAL 007 Sud-Coréen a été abattu par la défense aérienne de l'URSS, tuant 270 passagers. L'explication était que l'"espion" a envahi l'espace soviétique, mais aucun document n'a été fourni à la communauté internationale pendant 10 ans.

En 1986, l'explosion de la navette Challenger a ralenti la mise à niveau du système. Les premiers satellites de la série "Block II" ont été lancés en 1989.

Le système GPS actuel opère avec une trentaine de satellites, dont certains constituant une réserve. La précision militaire est actuellement de quelques centimètres ou mieux, surtout grâce aux techniques différentielles, mais ceci ne sera pas discuté ici.

Les satellites ne sont pas géostationnaires (ceci signifierait des orbites à 36 mille km), mais suivent des orbites presque circulaires de rayon 26 560 km, avec période de 11.967 heures (11h, 58min). Les orbites sont inclinées à 55° par rapport à l'équayeur, et séparées une de l'autre par multiples de 60 degrés.
Au moins 3 satellites doivent être visibles de tout point sur la Terre, mais il faut savoir les identifier. Chacun émet en deux fréquences principales, \(L_1 = 1575.42\) MHz, et \(L_2 = 1227.6\) MHz (ce sont des multiples de la fréquence de base 1.023 MHz). Donc, la longueur d'onde est d'ordre de 20cm, ce sont des micro-ondes. (Si vous avez oublié : pour les fours à micro-ondes, c'est 12cm).

Il y a une troisième fréquence : \(L_5 = 1176.45\) MHz, parfois omise dans les publications.. Les satellites pèsent une tonne chacun, ont 5M de diamètre, et typiquement ne vivent plus de 10 ans. En début des '2000, le producteur principal était Rockwell.

Les antennes des satellites se dirigent en permanence vers la Terre, mais le perte d'énergie reste très importante. Le signal émis ne contient aucune information directionnelle, la seule information présente est le marquage qui identifie le satellite, et deux motifs de bits pseudo-aléatoires (PRN), qui agit comme une balise temporelle, augmentés par des corrections temporelles et atmosphériques. Les "codes" pseudoaléatoires : C/A (Coarse Acquisition) et P (Precise, ou Protected) permettent d'identifier le temps et le satellite en question ; ceci est loin d'être simple, le récepteur reçoit les signaux mélangés, et puisque les fréquences sont les mêmes, on ne peut pas séparer les satellites comme on sépare les stations radio...

Le code C/A, avec la "chip frequency" (+/- la fréq. des bits) de 1.023 MHz se répète toute miliseconde, mais le code précis, 10 fois plus précis, chip frequency de 10.23MHz, est très long, et pour chaque satellite se répète au bout d'une semaine (et le code complet après 269 jours...). Le code C/A passe seulement sur la bande \(L_1\). Le taux de l'information navigationnelle (les éphémérides des satellites) est à peinde de 50 bauds. Les décodeurs civils et militaires sont différents, le code "P" est traité différemment. Mais les récepteurs doivent connaître tous les codes des satellites, on ne peut pas du jour au lendemain de changer ces codes.

Ceci est techniquement critique même avant le décodage. Le signal est très faible, comme les fares d'une voiture distante de 2500km, en dessous du bruit ambiant. Le récepteur, lors de la phase de synchronisation, accumule plusieurs échantillons, qui ont exactement la même forme.

Ainsi, puisque le bruit de fond est vraiment aléatoires, tout se compense, l'amplitude pourra grandir éventuellement comme \(\sqrt{N}\), où \(N\) est le nombre d'échantillons, tandis qu'un signal répétitif sera proportionnel à \(N\).

Faisons cette expérience, avec un simple pattern, et un bruit aléatoire. Le motif qui est 100 fois plus petit que l'amplitude du bruit, est totalement invisible sur ce fond aléatoire. Mais si on est capable de localiser la phase du début du motif, et ajouter, disons, 100 000 échantillons, alors l'image devient très différente.

Ajoutons le motif au bruit aléatoire. On n'observe aucune différence par rapport au bruit seul. Mais si on répète cette superposition du motif, toujours le même, avec des échantillons aléatoires, différentes du bruit, après avoir ajouté 100 000 échantillons, "le miracle arrive". L'amplitude du motif sera d'ordre de 1000, tandis que l'amplitude globale du bruit, de l'ordre de \(\sqrt{100000} \approx 300\). Le temps d'acquisition peut être assez long, car identifier la correcte phase n'est pas facile.

Normalement les étudiants ne s'intéressent pas trop aux détails numériques des équipements réels, et oublient tout. Cependant ceci fait la différence entre un professionnel, et quelqu'un qui a suivi un cours de "GPS pour les nuls", et similaire. En toute indépendance des solutions des équations de trilatération, on voit qu'une période "chip" de l'\(L_1\), qui dure une \(\mu s\), signifie 300 m de distance.
Le traitement du signaux, la reconnaissance des phases dans le pattern émis permet d'obtenir la résolution de 1% de la longueur du "chip" : 3 m. Pour le code "P" : 30 cm.

... et ici commence toute énorme couche d'ajustements, corrections et cryptage (éventuel), ce qui peut avaler 3 heures de ce cours... On a des corrections ionosphériques et troposphériques, relativité, etc.

Distances et temps

Le principe fondamental de radio-mesures des distances est très simple. Si on connaît le temps entre l'émission et la réception de l'onde EM (radio, IR, lumière, gamma...), on connaît la distance. Il y a des variantes : l'effet de Doppler (le changement des fréquences à cause de la vitesse) permet de localiser le moment où le satellite est plus proche de l'observateur, et on peut essayer de reconstruire le reste de la géométrie. Mais le GPS "classique" n'a pas besoin de cela. Voici l'idée de base, qui n'est pas directement applicable ! Il faudra des modfications.

On observe trois ou plus émetteurs, dont les positions dans le temps sont connues, donc on connaît la position des centres des sphères définissant les distances à travers le temps. Deux sphères (deux satellites) spécifient un cercle de leur intersection. Une troisième équation de trilatération donne des points individuels, un ou deux, le redondant peut être parfoit éliminé par des raisonnements de cohérence (par ex, en dessous de la surface terrestre, ou loin d'une mesure antérieure, etc.)

Pour mesurer la vitesse, on peut diviser la distance parcourue par le temps, mais accessoirement inclure les corrections Doppler. On essaie d'exploiter toutes les données disponibles.

Si on ajoute le quatrième, cinquième, etc. satellites, on élimine des ambiguïtés, et on améliore la précision ; les solutions sont toujours un peu approximatives.
Exercice. Étant données les positions des émetteurs : \((x_k,y_k,z_k)\) pour \(k=0,1,2\), et \(R_0, R_1, R_2\) : les distances des émetteurs, écrire un programme qui trouve \(x,y,z\), la position du récepteur (ou les deux possibles, ou aucune).

La solution, par programmation, est conceptuellement facile, techniquement pénible pour les étudiants ici, car dans notre établissement les étudiants n'apprennent vraiment jamais à programmer les constructions géométriques... Les équations seront : $$ (x-x_0)^2 + (y-y_0)^2 + (z-z_0)^2 = R_0^2\,, \\ (x-x_1)^2 + (y-y_1)^2 + (z-z_1)^2 = R_1^2\,, \\ (x-x_2)^2 + (y-y_2)^2 + (z-z_2)^2 = R_2^2\,. $$ et ceci décourage les faibles, plusieurs équations non-linéaires... Il y a plusieurs méthodes de traiter des équations avec trop de variables. Une raisonnable, est d'utiliser dans la mesure du posssible des équations vectorielles, sans la séparation des coordonnées, avec des variables vectorielles, en exploitant les invariances géométriques. Par exemple, si on a l'équation d'un plan et d'une droite, qui coupe la plan, il est avantageux de représenter le plan par une équation "implicite" \(\vec{A} \cdot \vec{x} = D \), ou : \(A \cdot x + B \cdot y +C\cdot z = D\). \(\vec{A}\) est normal (perpendiculaire) au plan.

Alors, si la droite est dans la représentation paramétrique, \(\vec{x} = \vec{x_0} + t \vec{v} \), où \(\vec{v}\) est la direction de la droite (un vecteur normalisé, de longueur 1), et le paramètre \(t\) est la longueur du segment, le déplacement le long de la droite :

\(x = x_0 + t v_x; y=y_0+t v_y; z=z_0+tv_z\), [typiquement : \(|\vec{v}|=1\)],

alors il est inutile, et pénalisant d'utiliser les coordonnées \(x,y,z\) afin de résoudre ces équations. Il faut, si possible, utiliser la notation vectorielle.

Nous aurons : \(\vec{A}\cdot(\vec{x_0} + t \vec{v}) = D \), ce qui donne \(t (\vec{A}\cdot\vec{v}) = D - \vec{A} \cdot \vec{x}_0\). Donc, \(t = (D - \vec{A} \cdot \vec{x}_0)/(\vec{A}\cdot\vec{v})\), et l'équation de la droite nous donne le point comme vecteur. On n'a pas séparé les coordonnées, ce qui rend la programmation plus compacte, et moins susceptible aux erreurs.

De la même façon, on peut obtenir l'intersection d'une droite, paramétrique comme ci-dessus, et d'une sphère, implicite : \((\vec{x} - \vec{x}_c)^2 = R^2\).
$$ (\vec{x}_{0c} + t \vec{v})^2 = R^2\,; \quad t^2 + 2t [\vec{v}\cdot \vec{x}_{0c}] + \vec{x}_{0c}^2 - R^2 = 0\, ; \quad \vec{x}_{0c} = \vec{x}_0 - \vec{x}_c\,. $$ ce qui donne : \(t_{\pm} = -(\vec{v}\cdot \vec{x}_{0c}) \pm \sqrt{(\vec{v}\cdot \vec{x}_{0c})^2 +R^2-\vec{x}_{0c}^2}\).

Je n'envisage pas de vous fatiguer avec du numérique peu inspirant, mais pour votre culture générale il est utile de connaître quelques astuces du métier. Dans ce domaine on aura besoin toujours de nouvelles techniques, si quelqu'un pense qu'il y a des librairies et des paquetages pour tout, il a tort.

Dans des calculs de ce type on doit garder dans la mémoire le fait que le choix du système de coordonnées est arbitraire. Par exemple, on peut décaler tout, et rédéfinir les coordonnées par rapport à un point de référence connu, par ex \(\vec{x}_0\) (Le premier satellite se retrouve à l'origine ; à la fin des calculs on restaure les coordonnées "canoniques"). Ainsi, la première équation se réduit à \((\vec{x}-\vec{x}_0)^2 \to \vec{x}^2 = R_0^2\).

Ensuite, nous pouvons effectuer la rotation des coordonnées, de manière à ce que le satellite 1 possède seulement la coordonnée \(x_1\), avec \(y_1 = z_1 = 0\). La rotation du plan des satellites autour de la droite \(\vec{x}_0 \to \vec{x}_1\) peut "amener" la coordonnée \(z_2\) à zéro, et les équations deviennent plus mangeables.
Ceci est populaire, mais c'est du "bricolage". Avec : $$ x^2+y^2+z^2 = R_0^2\,; \quad (x-x_1)^2 +y^2 + z^2 = R_1^2\,, $$ après la soustraction : \(R_0^2 - R_1^2 = 2 x\cdot x_1 - x_1^2\), ou \(x = (R_0^2 - R_1^2 + x_1^2)/(2 x_1)\). Ceci permet ensuite de trouver \(y\) et \(z\), mais c'est plus difficile, et demande plusieurs transformations accessoires. C'est un bon exercice en programmation...

Conceptuellement, mon approche préférée est géométriquement plus invariante et élégante. Pour l'équation \(i, i=1, 2, 3, \ldots\) on transforme $$ (\vec{x} - \vec{x}_i)^2=R_i^2 = (\vec{x} - \vec{x}_0 - (\vec{x}_i - \vec{x}_0))^2)\,. \text{Ceci donne} \\ (\vec{x} - \vec{x}_0) \cdot \vec{x}_{i0} = \frac{1}{2}\left(R_0^2 + R_{i0}^2 - R_i^2\right)\,. $$ et ceci est une équation linéaire. On obtient une équation (zéro) quadratique et deux linéaires. Avecplus de satellites : plus d'équations linéaires.

Cependant, même sans corrections atmosphériques, et l'augmentation de précision par des équations redondantes, avec 4, 8, etc. satellites, un problème de taille est le suivant : le récepteur GPS n'a pas d'horloge atomique, et sa vraie synchronisation avec le temps-satellite est impossible.

On mesure le temps, mais on calcule des "pseudo-distances", contenant le décalage entre le vrai intervalle temporel du signal, et les données mesurées.

Dans ces circonstances, il n'est pas recommandé d'utiliser des techniques numériques exactes, mais de chercher la meilleure solution combinée de toutes les équations par la méthode de moindres carrées. [Une courte digression si vous n'avez jamais entendu parler de cette méthode. C'est essentiel pour tout scientifique et tout ingénieur. Sérieusement.]

Dans le contexte de nos problèmes, on peut considérer que la partie linéaire du calcul est un ensemble d'équations : \(A \vec{x} \approx \vec{b}\), où $$ A = \pmatrix{x_1 - x_0 & y_1-y_0 & z_1-z_0\cr x_2 - x_0 & y_2-y_0 & z_2-z_0\cr \vdots & \vdots & \vdots \cr x_n - x_0 & y_n-y_0 & z_n-z_0\cr }\,; \quad \vec{x} = \pmatrix{x-x_0 \cr y-y_0 \cr z-z_0}\,; \quad \vec{b} = \pmatrix{b_{10} \cr b_{20} \cr \vdots \cr b_{n0}}\,. $$ Ce système est "sur-conditionné", plus d'équations que de variables, est la solution exacte typiquemen n'existe pas. Cependant, nous pouvons essayer de trouver la meilleure solution. Ceci est un problème-type pour la robotique, pour l'optimisation des trajectoires ou forces, etc. On procède comme suit.

La valeur \(\vec{b}-A\vec{x}\) doit être zéro, mais ne l'est pas. Alors on définit \(S=(\vec{b}-A\vec{x})^T (\vec{b}-A\vec{x})\), et on cherche tel \(\vec{x}\) que \(S\) soit minimalisé.

On doit résoudre l'équation dite "normale" : \(A^T A\vec{x} = A^T \vec{b}\), et si la matrice de ce système est "bien conditionnée" (pas de singularités mathématiques proches, det pas trop proche de zéro), alors on obtient \(\vec{x} = (A^T A)^{-1} A^T \vec{b}\). Sinon, la situation reste gérable, mais je vous épargne les détails. De toute façon, tout n'est qu'un exercice algébrique.

Pour terminer cette section : supposons que l'équipement nous donne les coordonnées canoniques WGS84, \(X,Y,Z\). Mais nous avons besoin des coordonnées géodésiques, \(\varphi, \lambda, h\), avec \(h\) : la hauteur depuis l'ellipsoïde de répère, selon la normale à la surface.

Le problème est que les coordonnées géodésiques et géocentriques sont différentes.
La normale à l'elipsoïde ne passe pas par le centre de la Terre.

Les calculs – conversions entre le Cartésien et le géodésique ne sont pas difficiles, mais posent quelques problèmes de stabilité. Pour votre culture générale, voici comment on le fait.

La construction de \(X,Y,Z\) depuis \(\varphi\), etc., est simple. Avec : \(N(\varphi) = \frac{a}{\sqrt{1 - e^2\sin(\varphi)^2}}\) : $$ X = (N(\varphi)+h)\cos(\phi)\cos(\lambda)\,; \quad Y = (N(\varphi)+h)\cos(\phi)\sin(\lambda)\,;\\ Z = (N(\varphi)(1-e^2)+h)\sin(\phi)\,. $$ ...et ceci ne se renverse pas facilement...

Ce qui est parfois difficile à accepter par les "profanes", est que même si un problème admet une solution "exacte" (pas d'approximations numériques), par ex. la solution de Ferrari d'une équation quartique, ceci est moins avantageux qu'une méthode approximative, itérative, jusqu'à la convergence. La méthode exacte parfois engendre des formules longues et abjectes, et numériquement peu stables. Les itérations bien conçues sont simples et rapides. Ici : $$ r=\sqrt{X^2+Y^2}\,; \quad \varphi = \varphi_0 = \arctan(Z/r)\,, $$ et ensuite on itère quelques fois (2 - 5) : $$ h = \frac{r}{\cos(\varphi)}-N\,; \quad \varphi=\arctan\left(\frac{Z}{r} \left(1 - \frac{e^2 N}{N+h}\right)^{-1} \right)\,. $$ \(N\) comme avant, est aussi recalculé. Ceci converge rapidement pour \(h\) petit, ce qui est presque toujours le cas, et pour \(\varphi\) pas trop proche des pôles (90°) ; pour l'usage typique, c'est aussi le cas. Est-ce qu'un utilisateur lambda a besoin de cela? Si quelqu'un opère avec les données "crus", de bas niveau, fournis par le récepteur GPS, oui. Les logiciels GPS de base le font, et peut-être certains parmi vous construiront un jour quelque chose d'intéressant, qui combine les données GPS avec GLONASS, avec Galileo, avec des corrections différentielles... Pendant quelques années le marché sera assez chaud.