FR 3150875 A1 20250110 FR 2307295 A 20230707 Méthode et système de positionnement GNSS Méthode de positionnement différentiel relatif GNSS pour positionner une antenne mobile par rapport à une antenne de référence par le calcul des différences de phase entre les signaux reçus par ces deux antennes et provenant d’une pluralité de satellites, comportant : le traitement du signal numérisé (Src1) reçu par l’antenne de référence dans une première boucle (10) et le traitement du signal numérisé (Src2) reçu par l’antenne mobile dans une deuxième boucle (20), chacune des boucles faisant subir au signal des opérations de démodulation de fréquence, de démultiplexage de code, de corrélation et de calcul de discriminant(s), la première boucle estimant la phase du signal servant à ajuster la fréquence porteuse dans la démodulation de fréquence pour chacune des boucles, la deuxième boucle calculant la différence de phase entre les deux signaux ; et l’estimation de la position de l’antenne mobile à partir de la différence de phase calculée. Figure pour l’abrégé : Fig. 6 Méthode et système de positionnement GNSS La présente invention concerne le domaine du positionnement GNSS (« Global Navigation Satellite System » en anglais), et concerne plus particulièrement le positionnement différentiel relatif GNSS. Un système de positionnement GNSS localise un utilisateur au sol à partir de signaux émis par une constellation de satellites et d’un récepteur qui les reçoit et les traite. Les satellites de navigation qui constituent le système de positionnement GNSS assurent une couverture mondiale à n’importe quel moment. Le positionnement GNSS différentiel consiste à positionner une antenne, dite « rover », par rapport à une antenne de référence, dite « base ». Le récepteur GNSS de l’antenne de référence envoie, par voie hertzienne, des corrections au récepteur de l’antenne à positionner. La précision d’un récepteur GNSS autonome est de l’ordre du mètre, car il calcule seul les positions des satellites et les différentes corrections du signal comme celles liées par exemple à la propagation dans l’ionosphère et la troposphère. Les systèmes de positionnement GNSS différentiel permettent d’améliorer la précision et la robustesse du positionnement. Un récepteur GNSS fournit une pseudo-distance satellite-récepteur à partir des mesures de phase de la porteuse du signal GNSS. Cette mesure de pseudo-distance par la phase est ambiguë et est affectée par différents biais. L’ambigüité de la mesure est la distance satellite- récepteur à l’allumage du récepteur. En effet, le récepteur fournit la variation de la pseudo- distance depuis sa mise en fonctionnement. De plus, le signal est affecté par plusieurs biais qui sont liés à la synchronisation de l’horloge récepteur avec le temps UTC, la synchronisation des horloges des satellites avec le temps UTC, la propagation du signal dans l’atmosphère, et l’erreur de positionnement des satellites. De nombreuses méthodes ont été proposées pour la résolution de l’ambiguïté de phase, et ce sujet est toujours d’actualité surtout pour le positionnement dynamique en temps réel. Parmi ces méthodes, on cite la technique RTK (« Real Time Kinematic » en anglais) ou cinématique en temps réel, décrite par exemple dans les publications Haiyang Li et al. “ A novel partial ambiguity method for multi-GNSS Real-Time Kinematic positioning ” et C.M. Silva et al. “ GNSS ambiguity resolution with ratio and fixed failure ratio tests for long baseline network RTK under ionospheric activity ”. Pour corriger l’erreur de synchronisation des récepteurs avec le temps UTC, la méthode RTK utilise la double différence de phase en temps réel. En effet, dans cette approche, le récepteur de référence et le récepteur à positionner calculent la différence des différences de pseudo-distance obtenues pour deux satellites visibles. La précision obtenue par ce mode de positionnement est de l’ordre du centimètre, mais dépend de nombreux facteurs tels que la ligne de base, le nombre de satellites utilisés et leur configuration géométrique. Comme en temps différé, il faut passer par l’étape d’initialisation des ambiguïtés entières. En temps réel, cette initialisation est réalisée en utilisant la méthode de résolution en vol AROF (« Ambiguity Resolution On the Fly » en anglais) décrite par exemple dans la publication Erwin Frei et al. “ Ambiguity Resolution On the Fly (AROF ); results , facts , limitations ”. Cette méthode repose sur des algorithmes de résolution rapide d’ambiguïtés qui utilisent cinq satellites au minimum. Enfin, la technique RTK peut être utilisée avec des lignes de base allant jusqu’à 10 km, mais au-delà, la résolution des ambiguïtés devient beaucoup plus compliquée. Une autre approche de résolution de l’ambiguïté de phase est la technique PPP (« Precise Point Positioning » en anglais) décrite par exemple dans les publications Zumberge James et al. “ Precise Point Positioning for the Efficient And Robust Analysis of GPS Data from Large Networks ” et Jianghui Geng et al. “ Integer ambiguity resolution in Precise Point Positioning : method comparison ”. Le principe de la technique PPP est de positionner un récepteur qui est connecté à internet et reçoit en continu les positions précises et des corrections d’horloge des satellites. Le récepteur connecté à internet calcule sa position (la position de son antenne), la différence de temps de synchronisation du récepteur avec le temps UTC et l’ambiguïté de phase des mesures de pseudo-distances. À l’origine, la technique PPP est une technique de positionnement réalisé sur des mesures de phase en temps différé. Il a été démontré, à titre d’exemple dans la publication Zumberge James susmentionnée, qu’il est possible d’atteindre des précisions équivalentes à celles obtenues avec une technique différentielle pour un temps d’observation de plusieurs heures. Par la suite, grâce à la mise à disposition, sur internet, par l’IGS (« International GNSS Service » en anglais) des orbites des satellites et des informations d’horloges très précises, la technique PPP est devenue utilisable en mode cinématique et en temps quasi-réel avec une précision au centimètre. Le positionnement par la technique PPP se fait pour un récepteur à double fréquence avec une estimation précise d’orbite et d’horloge satellite issus de l’IGS. La double fréquence sur le récepteur représente un avantage dans la mesure où elle permet de corriger les effets ionosphériques. Cette méthode de positionnement utilise des mesures de phase ambiguës provenant du récepteur au lieu de celles provenant de la station de référence comme dans la technique RTK. Aussi, dans le cas de la technique PPP, aucune différence de mesure n’est effectuée. Ainsi, les erreurs qui s’éliminent par simple ou double différence sont estimées par l’utilisation de modèles et de traitement adaptés. Cela conduit à une résolution flottante des ambiguïtés (fixation des ambiguïtés à des valeurs réelles). Les publications Jianghui Geng susmentionnée, Maorong Ge et al. “ Resolution of GPS carrier-phase ambiguities in Precise Point Positioning (PPP) with daily observations ” et Laurichesse Denis et al. “ Integer ambiguity resolution on undifferenced GPS phase measurements and its application to PPP and satellite precise orbit determination ” proposent diverses méthodes pour la correction des ambiguïtés réelles. Plus récemment, plusieurs techniques de résolution d’ambiguïtés entières intégrées à la technique PPP ont vu le jour. On parle alors du Integer PPP (IPPP), tel que décrit dans la publication Laurichesse Denis susmentionnée, qui améliore la précision du positionnement et diminue le temps de convergence vers la position estimée. Enfin, la technique PPP est utilisable partout dans le monde et ne nécessite pas de station de référence. De ce fait, aucun problème lié par exemple à la stabilité des stations de base ou à leur géo-référencement n’intervient. Cependant, cette technique nécessite une longue période d’initialisation pour l’estimation des ambiguïtés (comme démontré dans la thèse Junbo Shi. Precise Point Positioning integer ambiguity resolution with decoupled clocks ), indispensable pour obtenir un positionnement de l’ordre du centimètre. Dans chacune des approches RTK et PPP, le positionnement horizontal est obtenu avec une précision centimétrique. La hauteur est obtenue avec une précision décimétrique. Pour ces précisions, la plus haute fréquence de mesure annoncée par les constructeurs est de l’ordre de 10 Hz. En pratique, on observe cette précision pour une fréquence de mesure de 1 Hz et au moins 1 heure d’observation du signal. La précision millimétrique est quant à elle obtenue après plusieurs heures de mesures. En effet, les traitements réalisés par le récepteur, qui estime l’ambiguïté des mesures de phase et qui filtre les signaux, amènent une corrélation temporelle des erreurs de positionnement. Cette corrélation augmente le temps d’intégration nécessaire à l’amélioration de la précision du positionnement. Dans les approches différentielles classiques par la phase, ce sont les mesures de doubles différences qui sont utilisées pour calculer la position absolue de l’antenne. La double différence a pour effet d’amplifier la puissance du bruit et de diminuer le nombre d’observations, car pour chaque observation, deux satellites sont nécessaires. Ces deux inconvénients, avec le calcul de l’ambigüité pour chaque couple de satellite, limitent la précision du positionnement des récepteurs actuels. Dans l’article Daniel Macias‐Valadez et al. “ Improving vertical GPS precision with a GPS-over- fiber architecture and real-time relative delay calibration ”, la fibre optique est utilisée pour relier les antennes. Une méthode de calibrage automatique permet de calculer les délais apportés par les fibres. Les signaux des antennes sont traités par un unique récepteur à deux entrées. La méthode classique de la simple différence est utilisée pour calculer la position de l'antenne. Il y a d’abord estimation de l'évolution de la phase déroulée en cycles (résolution d'ambiguïté, à la seconde) pour les deux signaux, puis calcul de la différence entre ces deux signaux. Il existe donc un besoin pour perfectionner les techniques existantes, notamment en termes de précision du positionnement et de durée d’acquisition pour obtenir la précision recherchée. Méthode de positionnement L’invention vise à répondre à cet objectif et a pour objet, selon l’un de ses aspects, une méthode de positionnement différentiel relatif GNSS pour positionner une antenne mobile par rapport à une antenne de référence par le calcul des différences de phase entre les signaux reçus par ces deux antennes respectives et provenant d’une pluralité de satellites GNSS, chacun des signaux ayant un délai de code et une fréquence porteuse, la méthode comportant : - la numérisation de ces signaux ; - le traitement du signal numérisé reçu par l’antenne de référence dans une première boucle de poursuite de code et de phase et le traitement du signal numérisé reçu par l’antenne mobile dans une deuxième boucle de poursuite de code et de phase, chacune des boucles faisant subir au signal des opérations de démodulation de fréquence, de démultiplexage de code, de corrélation et de calcul de discriminant(s), pour estimer le délai de code du signal, la première boucle estimant la phase du signal servant à ajuster la fréquence porteuse dans la démodulation de fréquence pour chacune des boucles, la deuxième boucle calculant la différence de phase Ψ entre les deux signaux à chaque cycle de mesure ; et - l’estimation de la position de l’antenne mobile, notamment par maximum de vraisemblance, à partir de la différence de phase Ψ ainsi calculée. L’antenne de référence n’est pas nécessairement une antenne fixe. Dans le cas où elle serait mobile, sa vitesse de déplacement est suffisamment faible pour ne pas introduire d’erreur inacceptable par rapport à la précision recherchée, et préférentiellement une erreur inférieure ou égale à 1 millimètre. Il faut comprendre par « antenne mobile » une antenne dont la position relativement à l’antenne de référence peut varier, de façon continue ou non. Un cycle de mesure, correspondant à une durée de code, est la durée nécessaire pour obtenir une information permettant de calculer la différence de phase. Il peut être d’une durée aux alentours de 1 ms, selon le type de code. Dans la suite du texte, la notation représente une valeur estimée. De préférence, avant l’estimation de la position de l’antenne mobile, on accumule sur plusieurs cycles de mesure, notamment sur au moins 9 cycles de mesure, des valeurs de différences de phase Ψ pour améliorer la précision de l’estimation de la position de l’antenne. L’intérêt est de pouvoir atteindre la précision recherchée, notamment une précision millimétrique de mesure. L’invention permet d’obtenir une plus grande précision et une plus grande fréquence des mesures (proche de 100 Hz) que les méthodes de positionnement existantes. En effet, grâce à l’invention, on obtient une précision millimétrique pour le positionnement. Un avantage sous-jacent au positionnement précis est l’obtention de mesures à hautes fréquences. En fait, plus la mesure est précise, plus la fréquence de mesure doit être grande pour pouvoir observer l’évolution millimétrique des processus. La précision du positionnement et la fréquence des mesures sont améliorées d’un facteur 10 en comparaison des performances annoncées par les constructeurs. De plus, les performances de la méthode selon l’invention ont été définies théoriquement et les résultats de positionnement obtenus sur les données réelles sont en accord avec les résultats théoriques. La méthode selon l’invention est plus simple que les méthodes différentielles classiques, car elle ne nécessite pas de calculer l’ambigüité de phase et la synchronisation du récepteur avec le temps UTC. Dans la présente invention, on utilise la simple différence de phase ce qui permet de corriger les biais dans la mesure de pseudo-distance qui sont liés à la propagation dans l’atmosphère, à la synchronisation des horloges satellites avec le temps UTC et le biais d’erreur de positionnement des satellites. De plus, l’implémentation matérielle du système selon l’invention permet de synchroniser la numérisation des signaux en provenance des deux antennes. Dans ce cas, la simple différence de phase permet aussi de retirer le biais d’erreur lié à la synchronisation du récepteur avec le temps UTC. Les différences de phase Ψ associées à chaque satellite dépendent préférentiellement des composantes Nord N A2 et Est E A2 de la distance au sol entre les deux antennes dans un repère cartographique de projection Lambert I et de la hauteur h entre celles-ci par la relation d’évolution linéaire circulaire multiple : + k étant un entier relatif, λ étant la longueur d’onde d’un signal GNSS, étant l’angle d’élévation du satellite correspondant vu par l’antenne de référence, étant l’angle d’azimut de ce satellite vu par l’antenne de référence et étant un bruit de mesure, notamment distribué suivant la distribution circulaire normale. Le fait que la phase est observée modulo 2π permet de supprimer l’ambiguïté entière. De préférence, la méthode selon l’invention comporte, préalablement à la numérisation des signaux, l’introduction d’un retard dans le signal reçu par l’une des antennes par rapport au signal reçu par l’autre antenne, via une ligne à délai, le retard étant suffisant pour que les fonctions de corrélation de ces signaux ne soient pas superposées. De préférence, la méthode selon l’invention comporte la somme du signal retardé et du signal non retardé et la numérisation du signal résultant via un numériseur unique, notamment un numériseur monovoie. Ainsi, la méthode selon l’invention est plus économique que les méthodes différentielles classiques, car elle ne nécessite qu’un seul récepteur et pas de liaison hertzienne. Par ailleurs, la méthode selon l’invention est plus générale que les méthodes différentielles classiques, car elle permet d’estimer la position de plusieurs antennes par rapport à l’antenne de référence. En effet, on peut utiliser plusieurs lignes à délai de distances différentes pour séparer les signaux en provenance de plusieurs antennes. Ces signaux peuvent être ensuite additionnés et numérisés par un récepteur mono-entrée. Dans un mode de réalisation, la méthode selon l’invention comporte la conversion des signaux reçus par les antennes en signaux optiques via des photo-modulateurs, la ligne à délai étant en fibre optique. La méthode selon l’invention peut comporter l’acquisition de manière indépendante du signal retardé et du signal non retardé numérisés, en générant pour chacun de ces signaux un signal répliqué local à corréler avec le signal associé pour estimer le délai de code de ce signal et sa fréquence porteuse et pour identifier le satellite correspondant, la corrélation avec le signal répliqué étant notamment calculée pour plusieurs fréquences du signal répliqué. La démodulation de fréquence dans chaque boucle comprend préférentiellement la multiplication du signal par deux fréquences porteuses en quadrature de phase obtenues à partir d’un générateur de fréquence porteuse et d’un déphaseur, ce générateur utilisant initialement la fréquence porteuse estimée à l’étape d’acquisition puis ajustant la fréquence à sa sortie par la phase estimée par la première boucle. Le démultiplexage de code dans chaque boucle comprend préférentiellement la multiplication du signal démodulé par trois codes C/A différents par leurs décalages, à savoir un code synchrone avec le signal démodulé, un code avancé par rapport au signal démodulé d’un décalage de d/2 et un code retardé par rapport au signal démodulé d’un décalage de -d/2, d étant le décalage entre le code avancé et le code retardé, ces codes étant obtenus à partir d’un générateur de code, ce générateur de chaque boucle ajustant les trois codes C/A à sa sortie par le délai de code estimé par la boucle. De préférence, la corrélation dans chaque boucle comprend l’intégration du signal démultiplexé pour produire deux fonctions (encore appelées composantes) de corrélation I et Q. La première boucle peut comprendre un discriminant de code et un discriminant de phase calculant respectivement la différence de délais de code Δ et la différence de phases Δ comme suit : Δ = D ( ), où la fonction D est un discriminant de code cohérent, la phase estimée par la première boucle et le délai de code estimé par la première boucle ; et Δ = . Il est à noter que les grandeurs entre parenthèses dans l’expression de la composante I sont les grandeurs dont dépend cette composante. Un exemple de discriminant cohérent est la fonction différence. Dans un mode de réalisation, la différence de délais de code Δ et la différence de phase Δ en sortie des discriminants de la première boucle sont chacune intégrée pour donner respectivement le délai de code estimé et la phase estimée. La phase et le délai de code estimés par la première boucle peuvent servir à décoder le message de navigation du signal reçu par l’antenne de référence pour dater ce signal en vue d’associer les angles d’élévation et d’azimut des satellites aux calculs de différences de phase Ψ. La deuxième boucle peut comprendre un discriminant de code et un discriminant de phase calculant respectivement la différence de délais de code Δ et la différence de phase Ψ comme suit : Δ = D ( ), où la fonction D est un discriminant de code incohérent, la phase estimée par la première boucle et le délai de code estimé par la deuxième boucle ; Ψ = - , où est la phase du signal reçu par l’antenne mobile, est la phase du signal reçu par l’antenne de référence, bruit résultant du calcul de Ψ et de l’estimation de . Un exemple de discriminant incohérent est Early-Late Power : + ) - + ), E étant le code C/A avancé (Early) et L le code C/A retardé (Late). La différence de délais de code Δ en sortie du discriminant de code de la deuxième boucle peut être intégrée pour donner le délai de code estimé. Système de positionnement L’invention a encore pour objet, selon un autre de ses aspects, un système de positionnement différentiel relatif GNSS comportant un récepteur connecté à une antenne de référence et une antenne mobile, pour positionner l’antenne mobile par rapport à l’antenne de référence par le calcul des différences de phase Ψ entre les signaux reçus par ces deux antennes et provenant d’une pluralité de satellites GNSS, chacun des signaux ayant un délai de code et une fréquence porteuse, le récepteur comprenant : - un numériseur unique, notamment un numériseur monovoie, pour numériser les signaux ; - un module de traitement numérique, comportant : -- deux boucles de poursuite du code et de la phase, une première boucle configurée pour traiter le signal reçu par l’antenne de référence et une deuxième boucle configurée pour traiter le signal reçu par l’antenne mobile, chacune des boucles comportant un démodulateur de fréquence, un démultiplexeur de code, un corrélateur et au moins un calculateur de discriminant, pour estimer le délai de code du signal, la première boucle étant configurée pour estimer la phase du signal servant à ajuster la fréquence porteuse dans le démodulateur de fréquence pour chacune des boucles, la deuxième boucle étant configurée pour calculer la différence de phase Ψ entre les deux signaux à chaque cycle de mesure ; et -- un estimateur de la position de l’antenne mobile, notamment par maximum de vraisemblance, à partir de la différence de phase Ψ ainsi calculée. De préférence, le système est configuré pour accumuler, avant l’estimation de la position de l’antenne mobile, notamment sur au moins 9 cycles de mesure, des valeurs de différences de phase Ψ pour améliorer la précision de l’estimation de la position de l’antenne. Le système comporte préférentiellement en amont du numériseur une ligne à retard configurée pour retarder le signal reçu par l’une des antennes par rapport au signal reçu par l’autre antenne de sorte que les fonctions de corrélation de ces signaux ne soient pas superposées. Dans un mode de réalisation, le système comporte un sous-module d’acquisition du signal retardé et du signal non retardé numérisés, configuré pour générer pour chacun de ces signaux un signal répliqué local à corréler avec le signal associé pour estimer le délai de code de ce signal et sa fréquence porteuse et pour identifier le satellite correspondant. De préférence, le démodulateur de fréquence comprend : - un générateur de fréquence porteuse, ce générateur étant configuré pour utiliser initialement la fréquence porteuse du signal reçu par l’antenne de référence estimée par le sous-module d’acquisition puis ajuster la fréquence à sa sortie par la phase estimée par la première boucle ; - un déphaseur configuré pour introduire un déphasage de à la porteuse générée ; et - deux multiplieurs configurés pour multiplier le signal en entrée par les deux porteuses en quadrature de phase. De préférence, le démultiplexeur comprend : - un générateur de code configuré pour générer trois codes C/A différents par leurs décalages, à savoir un code synchrone avec le signal démodulé, un code avancé par rapport au signal démodulé d’un décalage de d/2 et un code retardé par rapport au signal démodulé d’un décalage de -d/2, d étant le décalage entre le code avancé et le code retardé, le générateur de code étant configuré pour ajuster les trois codes C/A à sa sortie par le délai de code estimé par la boucle ; - deux multiplieurs configurés pour multiplier chacune des composantes en quadrature de phase du signal démodulé par les trois codes en sortie du générateur de code. L’invention pourra être mieux comprise à la lecture de la description détaillée qui va suivre, d’exemples non limitatifs de mise en œuvre de celle-ci, et à l’examen du dessin annexé, sur lequel : La est une vue géométrique 3D d’un système comportant une antenne de référence, une antenne mobile et un satellite GNSS dans différents repères référentiels ; la représente : (a) une vue partielle de la et (b) une vue de dessus de la en deux dimensions ; la est un schéma bloc d’un exemple de module de traitement analogique d’un système selon l’invention ; la est un schéma bloc d’un exemple de module de traitement numérique d’un système selon l’invention ; la représente schématiquement un exemple de sous-module d’acquisition d’un système selon l’invention ; la illustre schématiquement un exemple de deux boucles de poursuite de code et de phase faisant partie d’un système selon l’invention ; la représente schématiquement un exemple de générateur de code faisant partie d’un système selon l’invention ; la représente la courbe d’évolution de la différence de phase entre deux signaux reçus par l’antenne de référence et l’antenne mobile, pour chaque satellite visible, ainsi que la courbe de son modèle théorique non bruité ; la représente la fonction de corrélation d’un signal composite obtenu après introduction d’un retard ; la représente les courbes de la racine carrée de la CRLB (Cramér-Rao Lower Bounds) et de l’écart type d’estimation en fonction de κ, paramètre de concentration du bruit de mesure, (a) dans le cas d’un seul satellite et (b) dans le cas de la fusion des mesures obtenues par deux satellites ; la représente les courbes de la racine carrée de la CRLB pour différentes variations d’évolution de l’azimut et de l’élévation, (a) dans le cas d’une variation de l’azimut, et (b) dans le cas d’une variation de l’élévation ; la illustre schématiquement la position des capteurs (satellites) sur la sphère unité : en trait gras continu V Az = 2π et V θ = π/2, en trait discontinu ---- V Az = 0 et V θ = π/2, en trait discontinu -.-.- V Az = 2π et V θ = 0 ; la représente schématiquement un exemple de dispositif expérimental pour un système selon l’invention, positionné sur un toit ; la représente les courbes d’estimation des trois composantes de la position (NA2, EA2, h) pendant 20 s ; et la illustre schématiquement les écarts d’estimation des paramètres NA2 et EA2. Description détaillée La est une vue géométrique 3D d’un système comportant une antenne de référence A1, une antenne mobile A2 et un satellite GNSS S dans différents repères référentiels. La représente : (a) une vue partielle de la et (b) une vue de dessus de la en deux dimensions. L’antenne A2 est placée dans un plan P2 parallèle au sol représenté par un plan P1. L’antenne A1 est placée à une hauteur h de l’antenne A2. Les deux antennes reçoivent des signaux GNSS Sr1 et Sr2 arrivant du satellite S et sont connectées à un récepteur GNSS (non représenté sur la ) où le traitement du signal est effectué. Les angles d’élévation et d’azimut pour le satellite S sont respectivement notés θ et Az. On représente sur la la différence de chemin δ entre les signaux Sr1 et Sr2. Le point B est le projeté orthogonal de l’antenne A2 sur la droite qui représente la trajectoire du signal allant du satellite S à l’antenne A1. La distance entre A1 et B est la différence de chemin δ. On montre plus loin que cette différence de chemin peut être observée avec la différence de phase entre les signaux Sr1 et Sr2. On considère dans cette représentation trois repères cartographiques. Le repère (O, E, N, z) est le repère de référence qui positionne les différents éléments dans le système de coordonnées cartographiques Lambert I. Le repère (O, x, y, z) est défini tel que le plan vertical (O, x, z) contient la position des deux antennes notées A1 et A2. Le repère (O, X, Y, z) est défini tel que le plan vertical (O, X, z) contient l’antenne A1 et le satellite S. Γ est l’angle de rotation autour de l’axe (O, z) entre les repères (O, x, y, z) et (O, X, Y, z). Γ′ est l’angle de rotation autour de l’axe (O, z) entre les repères (O, E, N, z) et (O, x, y, z) tel que l’on a un angle Az entre les axes (O, X) et (O, N). ∆x est l’abscisse de l’antenne A2 dans le repère (O, x, y, z). On a alors pour la géométrie de la la relation entre les repères (O, x, y, z) et (O, X, Y, Z) qui est définie par la rotation suivante (A.1) : On a : Et (A.2) : Du système d’équations (A.1), on déduit aussi : Soit (A.3) : De même, du système d’équation (A.1), on peut déduire (A.4) : L’angle d’élévation θ étant l’angle entre le plan horizontal contenant l’antenne A1 et la droite (A1B), on a : (A.5) : (A.6) : B étant le projeté orthogonal de l’antenne A2 sur la ligne allant de l’antenne A1 à un satellite en vue considéré, on a (A.7) : Soit (A.8) : (A.9) : Par factorisation, on a (A.10) : Soit (A.11) : Par conséquent (A.12) : On obtient (A.13) : Par simplification, on a (A.14) : Enfin, on a (A.15) : A partir des équations (A.2) et (A.6), on peut écrire : (A.16) : (A.17) : Or d’après la (b) : Donc (A.18) : (A.19) : Par application des formules trigonométriques, on a (A.20) : De la (b), on a les expressions suivantes : (A.21) : (A.22) : En remplaçant les expressions (A.21) et (A.22) dans (A.20), on arrive à (A.23) : Finalement, on montre que l’expression de la différence de chemin δ dans le repère cartographique (O, E, N, z) est donnée par (A.23) où et sont les composantes Nord et Est de ∆x. Dans la présente invention, on propose d’estimer ces composantes à partir des observables de phase des signaux Sr1 et Sr2. La est un schéma bloc d’un exemple de module de traitement analogique 1 (« front-end » en anglais) d’un système selon l’invention. L’objectif de ce module analogique 1 est de numériser les signaux GNSS reçus par les deux antennes. Le module analogique 1 comprend un bloc d’alimentation des antennes 50, deux coupleurs 501, 502 reliés au bloc d’alimentation 50, une ligne à retard 60, un sommateur 70 et un numériseur 80. Dans la mise en œuvre d’un démonstrateur par les inventeurs, l’alimentation des antennes est fournie par deux récepteurs NovAtel qui permettent également de récupérer les éphémérides des satellites et les positions approximatives des antennes. La ligne à retard 60 est de préférence en fibre optique et est placée soit sur la voie du signal reçu par l’antenne A1 soit sur celle du signal reçu par l’antenne A2. La ligne à retard 60 a pour but de décaler dans le temps l’un des signaux. La longueur de la ligne à retard peut être d’au moins 1 km. Sa longueur est égale à 1.5 km dans ce démonstrateur. La longueur de la ligne à délai utilisée dans le travail expérimental des inventeurs introduit un retard d’environ 50 µs. Le sommateur 70 additionne le signal retardé et le signal non retardé. Le numériseur 80 numérise la somme des signaux en sortie du sommateur pour obtenir un signal composite Src. Le numériseur 80 est de préférence un récepteur monovoie, notamment de type 1 voie 1 bit et 1 bande GNSS. La carte de numérisation peut être celle développée par l’entreprise Syntony. Une telle carte enregistre les signaux GNSS sur deux bandes, L1 et L5, sous la forme d’échantillons quantifiés sur un ou deux bits. Dans les expérimentations des inventeurs, une fréquence d’échantillonnage de 25 MHz et une quantification sur 1 bit sont utilisées, et les signaux en sortie d’antenne d’une durée de 20s sont enregistrés. La représente la fonction de corrélation d’un signal composite obtenu après introduction d’un retard. Il s’agit d’une carte à délai Doppler (« Delay Doppler Map » en anglais). Dans cette représentation, on a en abscisse le délai de code en mètres, en ordonnée la fréquence Doppler en Hertz, et la valeur de corrélation définit la troisième dimension, représentée sous forme de nuances de gris. On peut observer sur la deux pics de corrélation séparés et associés au signal de la première antenne et au signal retardé de la deuxième antenne. En pratique, pour différentier les satellites, on utilise différents codes CDMA. Dans le récepteur de la présente invention, pour différentier les signaux en provenance des antennes, on décale dans le temps ces différents codes CDMA. La durée du code CDMA du signal GPS étant de 1 ms (300 km), ce décalage est sans ambiguïté s’il est inférieur à la milliseconde. La est un schéma bloc d’un exemple de module numérique 2 d’un système selon l’invention. Ce module numérique 2 traite le signal composite obtenu en sortie du numériseur 80 pour en extraire les mesures et paramètres utilisés pour calculer la position de l’antenne mobile. Le module numérique 2 comprend un sous-module d’acquisition 30, un sous-module de poursuite de code et de phase 100 et un sous-module d’extraction d’informations 40. Le sous-module d’acquisition 30 a en entrée le signal composite Src. Un exemple de ce sous-module 30 est décrit à la . L’objectif de l’étape d’acquisition est d’obtenir une initialisation des délais de code et fréquences Doppler des signaux satellitaires reçus. L’étape d’acquisition utilise le principe de corrélation pour l’estimation des délais de codes, des fréquences et pour identifier les satellites visibles. Ainsi, à la réception du signal, le récepteur génère un signal répliqué local qu’il corrèle avec le signal reçu. Le calcul de la corrélation peut se faire dans le domaine fréquentiel pour accélérer les calculs. La fonction de corrélation peut être calculée pour plusieurs fréquences du signal répliqué localement par le récepteur. On parle alors d’acquisition fréquentielle. La représente l’acquisition pour une fréquence donnée f i . Le sous-module 30 comprend un générateur 31 de fréquence porteuse f i , un déphaseur configuré pour introduire un déphasage de π/2 à la porteuse générée et deux multiplieurs configurés pour multiplier le signal en entrée par les deux porteuses ainsi obtenues en quadrature de phase. On note sur la FFT la transformée de Fourier Rapide, IFFT son inverse et FFT* le conjugué de la transformée de Fourier. Le calcul de la fonction de corrélation réalisé pour différentes fréquences permet d’obtenir les fonctions d’auto corrélations R(τ, f i ) où τ est le délai du code local généré par le récepteur. Dans la présente invention, l’étape d’acquisition est réalisée indépendamment sur les deux signaux reçus, mais l’étape de poursuite est réalisée conjointement sur les deux signaux. Le sous-module de poursuite de code et de phase 100 comprend deux boucles de poursuite de code et de phase : une première boucle 10 ayant en entrée le signal numérisé Src1 reçu par l’antenne de référence et une deuxième boucle 20 ayant en entrée le signal numérisé Src2 reçu par l’antenne mobile. Ces boucles sont illustrées schématiquement à la et sont décrites en détail plus loin. Le sous-module d’extraction d’informations 40 a en entrée les estimations de la phase et du délai de code du signal reçu par l’antenne de référence en sortie de la première boucle 10. Il peut extraire, à partir de ces estimations, le message de navigation du satellite S. Ce message est utilisé pour dater le signal. La datation du signal consiste à retrouver l’heure d’émission de chaque code CDMA, donc de chaque mesure de différence de phases. Cette heure d’émission est utilisée pour positionner les satellites et donc pour calculer en fonction du temps l’azimut et l’élévation de chaque satellite. Dans l’implémentation expérimentale effectuée par les inventeurs, les éphémérides des satellites sont obtenues à partir des fichiers RINEX fournis par un récepteur NovAtel. Chacune des boucles de poursuite de code et de phase 10, 20 comporte au moins un démodulateur de fréquence 200, un démultiplexeur de code 300, un corrélateur 400 et un calculateur de discriminant 106, 206, pour poursuivre et corriger le délai de code du signal associé. On représente les architectures des boucles 10 et 20 respectivement sur les figures 6(a) (figure du haut) et 6(b) (figure du bas). Le démodulateur de fréquence 200 comprend dans cet exemple un générateur de fréquence porteuse 101, un déphaseur 102 et deux multiplieurs 103. Le générateur 101 est configuré pour utiliser initialement la fréquence porteuse du signal Src1 estimée par le sous-module d’acquisition 30 puis ajuster la fréquence à sa sortie par la phase corrigée du signal Src1, produite en sortie de la première boucle 10. Le déphaseur 102 est configuré pour introduire un déphasage de à la porteuse générée. Les deux multiplieurs 103 sont configurés pour multiplier le signal en entrée par les deux porteuses en quadrature de phase. La première boucle 10 est configurée pour poursuivre et corriger la phase du signal reçu par l’antenne de référence en produisant une phase corrigée . La deuxième boucle est configurée pour utiliser cette phase corrigée afin de calculer la différence de phase Ψ entre les deux signaux Src1 et Src2. Le démultiplexeur 300 comprend dans cet exemple un générateur de code 104 et deux multiplieurs 103. Un exemple d’architecture de générateur de code est illustré à la . Le générateur de code 104 est configuré pour générer trois codes C/A différents par leurs délais de propagation, à savoir un code synchrone (noté P sur la figure 7) avec le signal démodulé, un code avancé par rapport au signal démodulé d’un décalage de d/2 (noté E sur la figure 7) et un code retardé (noté L sur la figure 7) par rapport au signal démodulé d’un décalage de -d/2, d étant le décalage entre le code avancé et le code retardé. Le générateur de code 104 est aussi configuré pour ajuster les trois codes C/A à sa sortie par le délai de code corrigé par la boucle associée 10, 20. Les deux multiplieurs 103 sont configurés pour multiplier chacune des composantes en quadrature de phase du signal démodulé par les trois codes en sortie du générateur de code. Le corrélateur 400 comprend dans cet exemple deux intégrateurs, chaque intégrateur référencé 105 étant configuré pour intégrer l’une des composantes en quadrature de phase du signal démultiplexé et produire une composante de corrélation I ou Q. La première boucle 10 comprend un discriminant de code et un discriminant de phase, schématisés par un même bloc 106, et configurés pour calculer respectivement la différence de délais de code Δ et la différence de phases Δ comme suit : Δ = D ( ), où la fonction D est un discriminant du code cohérent, la phase estimée par la première boucle et le délai de code estimé par la première boucle ; et Δ = . La première boucle 10 comporte un intégrateur final 107 configuré pour intégrer la différence de délais de code Δ et la différence de phases Δ en sortie des discriminants de la première boucle pour donner respectivement le délai de code corrigé qui est utilisé rétroactivement dans le générateur de code 104 de la première boucle, et la phase corrigée qui est utilisée rétroactivement dans le générateur de porteuse 101 de chacune des boucles 10, 20. La deuxième boucle 20 comprend un discriminant de code et un discriminant de phase, schématisés par un même bloc 206, et configurés pour calculer respectivement la différence de délais de code Δ et la différence de phases Ψ comme suit : Δ = D ( ), où la fonction D est un discriminant de code ; Ψ = - , où est la phase du signal reçu par l’antenne mobile, est la phase du signal reçu par l’antenne de référence, bruit résultant du calcul de Ψ et de l’estimation de . Il est à noter que la différence de phases Ψ est obtenue en boucle ouverte. La différence de phase Ψ est proportionnelle à la différence de trajet entre les signaux arrivant sur chaque antenne. En effet, un cycle de phase (2π radians) correspond en distance à la longueur d’onde λ d’un signal GNSS soit 19.04 cm. Dans ce contexte, Ψ est une mesure ambiguë dans [−π, π] de la différence de chemin que l’on peut noter : k étant un entier relatif, λ étant la longueur d’onde d’un signal GNSS, étant l’angle d’élévation du satellite et étant l’angle d’azimut de ce satellite. Cette expression est une régression circulaire linéaire multiple. On représente sur la l’évolution de la différence de phase et de son modèle d’évolution linéaire-circulaire. Les observations de différence de phase sont obtenues pour tous les satellites visibles. On observe sur la que le modèle d’évolution linéaire circulaire a une abscisse définie sur un intervalle [0, 1] et une ordonnée observée par plusieurs satellites positionnés différemment sur l’intervalle [−π, π]. On peut conclure que le modèle d’évolution est linéaire circulaire multiple à trois paramètres EA2, NA2 et h. Ces paramètres sont estimés par maximum de vraisemblance. En effet, si l’on se réfère à l’expression de la différence de phases Ψ établie théoriquement ci-dessus en lien avec la , on note que : La deuxième boucle 20 comporte un intégrateur final 207 configuré pour intégrer la différence de délais de code Δ en sortie du discriminant de code de la deuxième boucle pour donner le délai de code corrigé qui est utilisé rétroactivement dans le générateur de code 104 de la deuxième boucle 20. Ainsi, dans la présente invention, on utilise la simple différence de phases ce qui permet de corriger les biais dans la mesure de pseudo-distance qui sont liés à la propagation dans l’atmosphère, à la synchronisation des horloges satellites avec le temps UTC et le biais d’erreur de positionnement des satellites. De plus, l’implémentation matérielle du système selon l’invention permet de synchroniser la numérisation des signaux en provenance des deux antennes. Dans ce cas, la simple différence de phase permet aussi de retirer le biais d’erreur lié à la synchronisation du récepteur avec le temps UTC. Enfin, la phase est observée modulo 2π ce qui permet de supprimer l’ambiguïté entière. Les figures 10 à 12 concernent des résultats de simulations faites par les inventeurs. Pour évaluer l’estimateur du maximum de vraisemblance proposé dans un contexte simulé, on compare la CRLB (Cramér-Rao Lower Bounds) approchée (obtenue théoriquement) avec la variance de l’estimateur. Dans cette simulation, on considère 1000 réalisations du signal et 90 échantillons de l’angle d’élévation et de l’angle d’azimut. L’élévation varie entre 0 et π/2 et l’azimut entre 0 et 2π. On représente sur la la racine carrée de la CRLB approchée et l’écart type de l’estimation des paramètres NA2, EA2 et h en fonction de κ, paramètre de concentration du bruit de mesure, (a) dans le cas d’un seul satellite et (b) dans le cas de la fusion des mesures obtenues par deux satellites. A partir de la figure 10(a), on montre qu’un estimateur du maximum de vraisemblance défini dans le domaine circulaire peut atteindre la CRLB approchée de l’estimation des paramètres de la régression linéaire-circulaire multiple. Sur la figure 10(a), on considère 1 satellite et on fait varier en abscisse la puissance du bruit (où est approximativement l’inverse de l’écart type de la distribution du bruit sur la mesure de phase). La convergence de notre estimateur est obtenue avec 90 observations et la précision d’estimation obtenue est dans ce cas de l’ordre du millimètre. Sur la (b), on montre l’amélioration de la précision de l’estimation quand on utilise deux satellites. A partir de la (b), on montre qu’un estimateur du maximum de vraisemblance défini dans le domaine circulaire est bien l’opérateur de fusion dans le domaine circulaire, car sa précision atteint la CRLB approchée de l’opérateur de fusion. Cet estimateur combine les mesures capteurs en tenant compte de la précision de chaque capteur (les capteurs étant les satellites GNSS). L’observabilité des paramètres NA2, EA2 et h dépend des prédicteurs de la régression circulaire linéaire multiple qui dépendent eux-mêmes de l’azimut et l’élévation du satellite. On représente sur la la CRLB approchée pour différentes variations d’évolution de l’azimut et de l’élévation. Chaque domaine de variation est observé par 90 mesures. Les abscisses V Az et V θ de la représentent les valeurs maximum des intervalles de variations [0;V Az] et [0;V θ]. Sur la (a), V Az varie entre 0 et 2π. On remarque sur cette figure que l’erreur d’estimation augmente quand V Az tend vers 0. Cet effet est plus remarquable sur la CRLB de EA2. En effet, Az et θ définissent la position des capteurs sur une sphère de rayon unité dans le repère cartographique (O,E,N,z). La illustre schématiquement la position des capteurs (satellites) sur la sphère unité : en trait gras continu V Az = 2π et V θ = π/2, en trait discontinu ---- V Az = 0 et V θ = π/2, en trait discontinu -.-.- V Az = 2π et V θ = 0. Sur la , quand V Az varie entre 0 et 2π, les capteurs sont situés autour de la sphère sur la courbe en trait gras continu. Quand V Az tend vers 0, les capteurs ne parcourent qu’une partie du tour de la sphère. Sur cette même figure, quand V Az = 0, les capteurs parcourent un quart de cercle dans le plan (O,N,z) sur la courbe en trait discontinu ----. Dans ce cas, l’axe E et donc le paramètre EA2 n’est plus observable. Donc, quand V Az tend vers 0, l’identifiabilité du paramètre EA2 diminue et la CRLB de son estimation augmente. Sur la (b), V θ varie entre 0 et π/2. On remarque sur cette figure que l’erreur d’estimation augmente quand V θ tend vers 0. Sur la , quand V θ = 0, les capteurs parcourent un cercle dans le plan (O, E, N) sur la courbe en trait discontinu -.-.-. Dans ce cas, quand V θ tend vers 0, l’identifiabilité du paramètre h diminue et la CRLB de son estimation augmente. Pour évaluer les performances que l’on peut attendre du système de positionnement différentiel que nous proposons, nous évaluons l’estimateur proposé dans un contexte réel. Le contexte réel consiste à utiliser une vraie constellation satellite décrite dans un fichier RINEX. On reporte dans le tableau 1 les élévations et azimuts des satellites qui sont utilisés. Ces informations sont données pour les 5 configurations de cadran étudiées. Constellation de satellites PRN des satellites, élévations et azimuts Côté Ouest 14 (10°, -39°),10 (12°, -94°), 32 (33°, -53°), 25 (34°, -116°), 12 (80°, -121°) Côté Est 6 (9°, 86°), 15 (17°,174°), 17 (22°,41°), 19 (3°, 59°), 24 (69°, 115°) Côté Nord-Est 6 (9°, 86°), 17 (22°, 41°),19 (37°, 59°) Côté Sud-Ouest 10 (12°, -94°), 25 (34°, -116°), 12 (80°, -121°) Quatre cadrans Tous les satellites visibles Le tableau 2 présente l’évaluation du temps d’intégration nécessaire au calcul de l’estimateur du maximum de vraisemblance (MLE : « Maximum Likelihood Estimate » en anglais) défini dans le domaine circulaire par les expérimentateurs pour atteindre les CRLB. Les temps d’intégrations pour l’estimation de NA2, EA2 et h sont donnés en fonction de la constellation de satellites et de C/N0. Dans ce tableau, il est indiqué la CRLB et l’écart type de l’estimateur MLE proposé obtenu en simulation STD( ), STD ) et STD ), de NA2, EA2 et h. Constellation de satellites C/N 0 (dB Hz) Intégration (ms) CRLB ( (cm) STD ( ) (cm) CRLB ) (cm) STD ) (cm) CRLB ( (cm) STD ) (cm) Côté Ouest 40 17 0.27 0.28 0.22 0.22 0.25 0.25 35 25 0.43 0.44 0.35 0.35 0.4 0.4 30 60 0.45 0.49 0.38 0.42 0.43 0.46 Côté Est 40 19 0.21 0.21 0.25 0.25 0.26 0.26 35 45 0.25 0.26 0.31 0.32 0.26 0.26 30 72 0.34 0.37 0.42 0.45 0.44 0.46 Côté Nord-Est 40 100 0.26 0.27 0.13 0.14 0.36 0.36 35 180 0.37 0.37 0.19 0.19 0.52 0.53 30 500 0.37 0.4 0.19 0.2 0.51 0.53 Côté Sud-Ouest 40 140 0.39 0.4 0.1 0.11 0.1 0.1 35 500 0.4 0.4 0.11 0.11 0.1 0.11 30 950 0.48 0.49 0.13 0.14 0.13 0.14 Quatre cadrans 40 9 0.21 0.21 0.16 0.16 0.19 0.19 35 11 0.37 0.37 0.28 0.29 0.33 0.33 30 15 0.53 0.58 0.41 0.43 0.48 0.51 Dans le tableau 2, l’estimateur de régression multiple proposé est supposé avoir convergé lorsque les valeurs de STD( ), STD ) et STD ) sont proches de celles de leur CRLB respective. En effet, on suppose ici que si l’écart-type de l’erreur d’estimation est proche de sa valeur minimum définie par la CRLB, alors l’estimation obtenue est bien celle du maximum de vraisemblance. La deuxième colonne du tableau 2 indique le temps d’intégration minimal N pour atteindre les CRLB. On remarque que la période d’intégration augmente lorsque C/N 0 diminue. On peut observer dans les colonnes 4 à 7 de la première partie du tableau 2 et dans les colonnes 4 et 5 de la deuxième partie du tableau 2 que les meilleures performances, en termes de temps d’intégration, sont obtenues pour les configurations ayant le plus de satellites. La réduction de la période d’intégration est due à l’augmentation du nombre de satellites visibles. Enfin, on remarque que les écarts-types des estimations des paramètres de la régression sont toujours proches des CRLB approchées et que nous obtenons une précision millimétrique sur chacun des paramètres estimés. Par conséquent, on conclut que lorsqu’on combine les observations de plusieurs satellites avec une période d’intégration suffisamment longue, l’estimateur proposé converge. En effet, dans ce cas, les écarts-types des estimations des paramètres sont proches des CRLB approchées. La période d’intégration diminue à mesure que le nombre de satellites augmente. La suite de la description présente l’évaluation du positionnement sur données réelles. Etant donné qu’on souhaite estimer la position relative (NA2, EA2, h) de l’antenne A2 par rapport à l’antenne A1 dans les trois dimensions de l’espace, on a besoin, en plus des élévations des satellites, des azimuts des satellites visibles. Pour l’expérimentation effectuée par les inventeurs, l’ensemble des satellites visibles est utilisé. Les figures 13 à 15 présentent le dispositif expérimental et les résultats des expérimentations effectuées par les inventeurs. La représente schématiquement un exemple de dispositif expérimental pour un système selon l’invention, positionné sur un toit. Il est rapporté dans le tableau 3 les élévations et azimuts des satellites visibles et la puissance du signal GPS pour chaque satellite. La hauteur h entre les deux antennes du dispositif expérimental de la est d’environ 50 cm et la position relative horizontale de l’antenne A2 est d’environ NA2 = 42 cm, EA2 = 139 cm. La numérotation PRN (« Pseudo Random Noise » en anglais) permet d’identifier le satellite et est fournie dans le message de navigation du signal GPS. PRN 1 3 4 6 9 12 17 19 22 31 Az (°) 133.63 65.32 173.26 -55.34 -155.95 -19.99 -99.59 -72.51 81.45 54.17 θ (°) 35.53 71.15 57.78 28.16 25.49 5.3 45.68 42.02 48.59 20.41 C/N 0 (dB Hz) 47.63 53.07 49.67 46.01 46.26 43.32 49.97 49.6 50.33 47.26 On peut observer dans le tableau 3 que la puissance du signal GPS des satellites est supérieure à celle testée en simulation. Ainsi, étant donné la visibilité satellitaire sur 4 cadrans et la puissance des signaux GPS, le temps d’intégration du signal a été fixé à 9 ms. Pour cette durée d’intégration et cette configuration de satellite, les écarts-types théoriques des estimations de NA2, EA2 et h, obtenus à partir de données simulées, sont respectivement de 2.01 mm, 1.92 mm et 1.68 mm. On considère des périodes consécutives de 9 ms de signal pour une durée totale de signal de 20 s. A chaque période de 9 ms est associé un triplé de valeurs estimées , et . L’écart-type sur chaque paramètre estimé nous fournit une précision d’estimation pour chacun d’eux. On représente sur la les 2222 estimations obtenues pour les trois composantes de la position relative de l’antenne A2. Les résultats d’estimation des composantes de la position présentés sur la montrent une variation des valeurs d’ordre millimétrique. De plus, l’évolution des paramètres estimés est constante. Cela montre la stabilité de l’estimateur développé. On reporte dans le tableau 4 ci-dessous la précision estimée pour les trois composantes de la position. Dans le tableau 4, on remarque que les précisions obtenues sont millimétriques. On remarque aussi que la précision obtenue sur données réelles est légèrement meilleure que celle obtenue en simulation. Cela s’explique par le fait que le niveau de bruit utilisé pour les données simulées (C/N 0 = 40 dB Hz) est plus élevé que celui relevé sur les données réelles (C/N 0 = 50 dB Hz). Enfin, la différence de position des antennes obtenue par RTK avec le logiciel GrafNav de novAtel et après 2h d’enregistrement est NA2 = 42.01 cm, EA2 = 138.78 cm et h = 49.69 cm. Ce résultat est consistant avec celui obtenu par la méthode selon l’invention. Paramètre Moyenne (cm) Ecart-type (mm) 42.35 1.77 138.10 1.24 49.45 0.86 On représente sur la les écarts des paramètres NA2 et EA2 estimés par rapport à leurs valeurs moyennes. On représente en abscisse les variations sur en mm et en ordonnée les variations sur en mm. On observe sur la que le nuage de points se répartit de manière équilibrée autour de l’origine, ce qui montre que la méthode selon l’invention est aussi précise dans la direction Nord que dans la direction Est. On remarque aussi à partir des résultats reportés dans le tableau 4 que la méthode selon l’invention est plus précise pour l’estimation de la hauteur que pour la position horizontale. Ce résultat qui confirme les résultats obtenus par simulation est différent des résultats obtenus par les techniques classiques de positionnement différentiel. En effet, par les approches classiques, les erreurs de positionnement dans la direction Nord et dans la direction Est sont différentes et la précision sur la hauteur est nettement inférieure. L’étude la reproductibilité de la mesure est calculée sur des signaux acquis sur plusieurs jours du même mois (septembre 2020) à approximativement la même heure. Il est indiqué dans le tableau 5 ci-dessous l’estimation de NA2, EA2 et h obtenue pour 9 jours du mois de septembre 2020. Les estimations sont obtenues avec un temps d’intégration de 9 ms. On indique aussi dans le tableau 5 les moyennes et écarts-types des estimations. Jour du mois 03 04 07 08 10 11 12 14 17 Temps (UTC) 13h12 13h15 13h14 13h12 13h11 13h12 13h10 13h13 13h15 (cm) 43.17 42.88 43.29 42.41 43.28 42.49 43.24 42.17 43.03 (cm) 138.66 137.88 138.62 137.79 138.56 138.46 138.49 138 138.3 (cm) 49.17 49.1 49.18 49.35 49.28 49.29 49.24 49.28 49.3 Moyenne (cm) 42.89 138.31 49.25 Ecart-type (mm) 4.2 3.3 0.8 On remarque dans le tableau 5 que la précision du positionnement obtenue les 9 jours est millimétrique. En effet, pour une durée d’intégration de 9 ms, les écarts-types de NA2, EA2 et h sont respectivement de 4.2 mm, 3.3 mm et 0.80 mm. Ces résultats sont proches de ceux obtenus en simulation et en théorie. Le faible nombre de mesures ne permet pas de conclure qu’il existe une différence significative entre les résultats obtenus et la théorie. Finalement, on peut conclure que la méthode selon l’invention n’est pas significativement dépendante des conditions de propagation du signal qui changent chaque jour. En effet, la différence entre les positions obtenues chaque jour est de l’ordre du millimètre. L’invention n’est pas limitée à l’exemple de réalisation décrit ci-dessus, et peut s’appliquer au positionnement de plusieurs antennes mobiles par rapport à une antenne de référence. La méthode selon l’invention peut être utilisée pour le positionnement relatif précis en extérieur. Un exemple d’utilisation est le positionnement d’un canon sur un char. En effet, le positionnement du bout du canon est de première importance pour la précision du tir, et compte tenu de la portée des canons, ce positionnement doit être millimétrique. Par l’approche proposée selon l’invention, l’antenne de référence est située au centre de la tourelle et l’antenne à localiser sur le bout du canon. Une autre application est la mesure de déformation d’un avion en vol. En effet, on peut considérer une antenne de référence sur le cockpit et des antennes à localiser sur les ailes et parties qui se déforment pendant le vol. Le système de positionnement permet alors de surveiller ces déformations pendant le vol pour qu’elles restent acceptables par la structure de l’avion. On peut aussi envisager des applications dans le domaine de l’agriculture et des travaux publics qui utilisent des engins avec de multiples outils à positionner. Méthode de positionnement différentiel relatif GNSS pour positionner une antenne mobile (A2) par rapport à une antenne de référence (A1) par le calcul des différences de phase entre les signaux reçus (Sr1, Sr2) par ces deux antennes respectives (A1, A2) et provenant d’une pluralité de satellites GNSS, chacun des signaux ayant un délai de code ( et une fréquence porteuse ( , la méthode comportant : - la numérisation de ces signaux (Sr1, Sr2) ; - le traitement du signal numérisé reçu par l’antenne de référence (Src1) dans une première boucle de poursuite de code et de phase (10) et le traitement du signal numérisé reçu par l’antenne mobile (Src2) dans une deuxième boucle de poursuite de code et de phase (20), chacune des boucles faisant subir au signal des opérations de démodulation de fréquence, de démultiplexage de code, de corrélation et de calcul de discriminant(s), pour estimer le délai de code ( du signal, la première boucle (10) estimant la phase du signal ( ) servant à ajuster la fréquence porteuse dans la démodulation de fréquence pour chacune des boucles, la deuxième boucle (20) calculant la différence de phase Ψ entre les deux signaux à chaque cycle de mesure ; et - l’estimation de la position de l’antenne mobile (A2), notamment par maximum de vraisemblance, à partir de la différence de phase Ψ ainsi calculée. Méthode selon la revendication précédente, dans laquelle, avant l’estimation de la position de l’antenne mobile (A2), on accumule sur plusieurs cycles de mesure, notamment sur au moins 9 cycles de mesure, des valeurs de différences de phase Ψ pour améliorer la précision de l’estimation de la position de l’antenne. Méthode selon la revendication précédente, les différences de phase Ψ associées à chaque satellite (S) dépendant des composantes Nord N A2 et Est E A2 de la distance au sol (Δx) entre les deux antennes (A1, A2) dans un repère cartographique de projection Lambert I (o, E, N, z) et de la hauteur h entre celles-ci par la relation d’évolution linéaire circulaire multiple : + k étant un entier relatif, λ étant la longueur d’onde d’un signal GNSS, étant l’angle d’élévation du satellite correspondant vu par l’antenne de référence (A1), étant l’angle d’azimut de ce satellite vu par l’antenne de référence (A1) et étant un bruit de mesure, notamment distribué suivant la distribution circulaire normale. Méthode selon l’une des trois revendications précédentes, comportant, préalablement à la numérisation des signaux, l’introduction d’un retard dans le signal reçu par l’une des antennes par rapport au signal reçu par l’autre antenne, via une ligne à délai (60), le retard étant suffisant pour que les fonctions de corrélation de ces signaux ne soient pas superposées. Méthode selon la revendication précédente, comportant la somme du signal retardé et du signal non retardé et la numérisation du signal résultant via un numériseur unique (80), notamment un numériseur monovoie. Méthode selon l’une des deux revendications précédentes, comportant la conversion des signaux reçus par les antennes (A1, A2) en signaux optiques via des photo-modulateurs, la ligne à délai (60) étant en fibre optique. Méthode selon l'une quelconque des revendications 4 à 6, comportant l’acquisition de manière indépendante du signal retardé et du signal non retardé numérisés, en générant pour chacun de ces signaux un signal répliqué local à corréler avec le signal associé pour estimer le délai de code ( de ce signal et sa fréquence porteuse ( et pour identifier le satellite correspondant, la corrélation avec le signal répliqué étant notamment calculée pour plusieurs fréquences du signal répliqué. Méthode selon la revendication précédente, la démodulation de fréquence dans chaque boucle (10; 20) comprenant la multiplication du signal par deux fréquences porteuses en quadrature de phase obtenues à partir d’un générateur de fréquence porteuse (101) et d’un déphaseur (102), ce générateur utilisant initialement la fréquence porteuse ( estimée à l’étape d’acquisition puis ajustant la fréquence à sa sortie par la phase estimée ( ) par la première boucle. Méthode selon l'une quelconque des revendications précédentes, le démultiplexage de code dans chaque boucle comprenant la multiplication du signal démodulé par trois codes C/A (E, P, L) différents par leurs décalages, à savoir un code synchrone (P) avec le signal démodulé, un code avancé (E) par rapport au signal démodulé d’un décalage de d/2 et un code retardé (L) par rapport au signal démodulé d’un décalage de -d/2, d étant le décalage entre le code avancé et le code retardé, ces codes étant obtenus à partir d’un générateur de code (104), ce générateur de chaque boucle ajustant les trois codes C/A (E, P, L) à sa sortie par le délai de code estimé ( par la boucle. Méthode selon l'une quelconque des revendications précédentes, la corrélation dans chaque boucle comprenant l’intégration du signal démultiplexé pour produire deux composantes de corrélation I et Q. Méthode selon la revendication précédente, la première boucle (10) comprenant un discriminant de code (106) et un discriminant de phase (106) calculant respectivement la différence de délais de code Δ et la différence de phases Δ comme suit : Δ = D ( ), où la fonction D est un discriminant de code cohérent, la phase estimée par la première boucle et le délai de code estimé par la première boucle ; et Δ = . Méthode selon la revendication précédente, la différence de délais de code Δ et la différence de phase Δ en sortie des discriminants de la première boucle étant chacune intégrée pour donner respectivement le délai de code estimé ( et la phase estimée ( . Méthode selon la revendication 10, la deuxième boucle (20) comprenant un discriminant de code (206) et un discriminant de phase (206) calculant respectivement la différence de délais de code Δ et la différence de phase Ψ comme suit : Δ = D ( ), où la fonction D est un discriminant de code incohérent, la phase estimée par la première boucle et le délai de code estimé par la deuxième boucle ; Ψ = - , où est la phase du signal reçu par l’antenne mobile, est la phase du signal reçu par l’antenne de référence, bruit résultant du calcul de Ψ et de l’estimation de . Méthode selon la revendication précédente, la différence de délais de code Δ en sortie du discriminant de code (206) de la deuxième boucle (20) étant intégrée pour donner le délai de code estimé ( . Système de positionnement différentiel relatif GNSS comportant un récepteur connecté à une antenne de référence (A1) et une antenne mobile (A2), pour positionner l’antenne mobile (A2) par rapport à l’antenne de référence (A1) par le calcul des différences de phase Ψ entre les signaux reçus (Sr1, Sr2) par ces deux antennes (A1, A2) et provenant d’une pluralité de satellites GNSS, chacun des signaux ayant un délai de code ( et une fréquence porteuse ( , le récepteur comprenant : - un numériseur unique (80), notamment un numériseur monovoie, pour numériser les signaux (Sr1, Sr2) ; - un module de traitement numérique (2), comportant : -- deux boucles de poursuite du code et de la phase (10, 20), une première boucle (10) configurée pour traiter le signal reçu par l’antenne de référence (A1) et une deuxième boucle (20) configurée pour traiter le signal reçu par l’antenne mobile (A2), chacune des boucles comportant un démodulateur de fréquence (200), un démultiplexeur de code (300), un corrélateur (400) et au moins un calculateur de discriminant (106 ; 206), pour estimer le délai de code ( du signal, la première boucle (10) étant configurée pour estimer la phase du signal ( ) servant à ajuster la fréquence porteuse dans le démodulateur de fréquence (101) pour chacune des boucles, la deuxième boucle (20) étant configurée pour calculer la différence de phase Ψ entre les deux signaux à chaque cycle de mesure ; et -- un estimateur de la position de l’antenne mobile, notamment par maximum de vraisemblance, à partir de la différence de phase Ψ ainsi calculée. Système selon la revendication précédente, étant configuré pour accumuler, avant l’estimation de la position de l’antenne mobile (A2), notamment sur au moins 9 cycles de mesure, des valeurs de différences de phase Ψ pour améliorer la précision de l’estimation de la position de l’antenne. Système selon la revendication 15 ou 16, comportant en amont du numériseur (80) une ligne à retard (60) configurée pour retarder le signal reçu par l’une des antennes par rapport au signal reçu par l’autre antenne de sorte que les fonctions de corrélation de ces signaux ne soient pas superposées. Système selon la revendication précédente, comportant un sous-module d’acquisition (30) du signal retardé et du signal non retardé numérisés, configuré pour générer pour chacun de ces signaux un signal répliqué local à corréler avec le signal associé pour estimer le délai de code de ce signal et sa fréquence porteuse et pour identifier le satellite correspondant. Système selon la revendication précédente, le démodulateur de fréquence (200) comprenant : - un générateur de fréquence porteuse (101), ce générateur étant configuré pour utiliser initialement la fréquence porteuse ( du signal reçu par l’antenne de référence (A1) estimée par le sous-module d’acquisition (30) puis ajuster la fréquence à sa sortie par la phase estimée ( ) par la première boucle (10) ; - un déphaseur (102) configuré pour introduire un déphasage de à la porteuse générée ; et - deux multiplieurs (103) configurés pour multiplier le signal en entrée par les deux porteuses en quadrature de phase. Système selon l'une quelconque des revendications 15 à 19, le démultiplexeur (300) comprenant : - un générateur de code (104) configuré pour générer trois codes C/A (E, P, L) différents par leurs décalages, à savoir un code synchrone (P) avec le signal démodulé, un code avancé € par rapport au signal démodulé d’un décalage de d/2 et un code retardé (L) par rapport au signal démodulé d’un décalage de -d/2, d étant le décalage entre le code avancé et le code retardé, le générateur de code (104) étant configuré pour ajuster les trois codes C/A à sa sortie par le délai de code estimé ( par la boucle ; - deux multiplieurs (103) configurés pour multiplier chacune des composantes en quadrature de phase du signal démodulé par les trois codes C/A en sortie du générateur de code (104).