Décodage des signaux de satellites GPS reçus par récepteur de télévision numérique terrestre DVB-T

Nous avons déjà décrit plusieurs fois dans ces pages l’utilisation d’un récepteur de télévision numérique terrestre pour recevoir des signaux radiofréquences, traités par GNURadio pour en extraire les informations contenues dans la modulation de la porteuse [1]. Nous sommes allés jusqu’à démontrer la réception de signaux émis par des satellites en orbite basse, à quelques centaines de kilomètres d’altitude, encodés selon des modulations classiques de fréquence et d’amplitude. Les signaux GPS avaient été reçus au moyen de l’utilisation de gnss-sdr¹, un outil libre, mais trop complexe pour appréhender les subtilités du traitement de signaux issus de la constellation de satellites de navigation par la lecture des codes sources. Dans la prose qui va suivre, le montage expérimental pour acquérir les signaux se réduit à un récepteur de télévision numérique terrestre compatible avec rtl-sdr² - à savoir un récepteur radiofréquence de type E4000 ou plus probablement R820T(2) suivi d’un convertisseur analogique-numérique Realtek RTL2832U - une antenne GPS active telle que par exemple celle commercialisée par Lextronic³ pour 12 euros, et un T de polarisation formé d’un condensateur pour protéger l’entrée radiofréquence du récepteur (une centaine de picofarads) et une inductance (quelques microhenrys) pour empêcher la remontée du signal radiofréquence dans l’alimentation de l’antenne active qui sera fournie par un port USB. Notez que gnss-sdr ne peut fonctionner qu’avec le récepteur Elonics E4000, désormais obsolète, et que l’application au Rafael Micro R820T est nouvelle, d’autant plus que nous verrons que ce dernier est affecté d’un biais de fréquence énorme, encore plus important que le E4000 dont l’exactitude était déjà médiocre. Ainsi, l’investissement financier pour reproduire ces expériences est de l’ordre d’une trentaine d’euros.

Nous allons voir que le décodage du signal GPS se traduit par l’enregistrement du signal radiofréquence (quelques secondes suffiront), et l’intercorrélation par des copies locales des codes représentant tous les identifiants de satellites possibles. Puisque la démodulation d’un signal modulé en phase n’est possible que si la porteuse a été annulée par l’oscillateur local, la recherche de la présence d’un signal GPS dans une acquisition radiofréquence passe de plus par un calcul relativement lourd dans lequel nous balayons, pour chaque code identifiant chaque satellite, tous les décalages de fréquence Doppler possibles.

1. Introduction

Nous désirons pousser au cran supérieur la démonstration de la puissance du traitement numérique de signaux radiofréquence et la flexibilité du récepteur de télévision numérique terrestre (DVB-T) par l’analyse de signaux dans la bande radiofréquence occupée par le GPS [2, 3]. Le GPS - Global Postionning System - est une des constellations de satellites permettant la localisation du récepteur par triangulation, et fait partie des objectifs plus généraux du GNSS - Global Navigation Satellite Service - aussi mis en œuvre par la Russie avec GLONASS (selon un protocole de communication très différent qui ne sera pas abordé ici) et sera peut être un jour mis en œuvre par l’Europe avec Galileo, la Chine avec Beidou ou, moins avancé, par l’Inde avec IRNSS. Le positionnement du récepteur résulte d’une triangulation précise des signaux émis par plusieurs satellites dont la position dans l’espace est supposée très bien connue, et dont le temps de vol des signaux entre émetteur et récepteur est estimé avec précision. Rappelons qu’un signal électromagnétique se propage à environ 300 m par microseconde, et que le positionnement avec une résolution de 3 m nécessite d’estimer le temps de vol du signal à environ 10 ns près après un temps de vol de plus de 70 ms. La capacité à retrouver une information temporelle aussi précise permet de remonter à des informations bien plus intéressantes et riches que la position [4] : mentionnons par exemple la densité d’électrons dans l’ionosphère, la concentration d’eau dans l’atmosphère, voire les ondes propagées dans l’air suite à un séisme. Les signaux reçus au sol sont aussi exploitables pour des applications de RADAR⁴ dans lesquelles les propriétés du réflecteur sont déduites de l’analyse du signal direct reçu du satellite et du signal réfléchi par la cible : taux d’humidité (ou coefficient de réflexion) ou distance de la cible.

Finalement, les perspectives de reproduire au sol une horloge combinant les propriétés de stabilité des diverses sources atomiques embarquées dans les satellites ouvrent la perspective de diffuser un temps précis en de nombreux sites géographiquement distincts [5, 6], un objectif qu’il est nécessaire d’atteindre pour toutes les applications de localisation d’événements par interférométrie (radiotélescope, positionnement d’impacts de foudre par exemple⁵). Bien que nous appelions communément GPS la constellation de satellites émettant des signaux de fréquence très stables dans le temps, les horloges embarquées ont connu plusieurs évolutions pour en améliorer les performances et la robustesse : initialement équipé de trois horloges au rubidium, référence secondaire de fréquence, une quatrième horloge au césium (référence primaire puisque la seconde est définie à partir de la fréquence de transition d’un atome de césium [7]) a été ajoutée une fois une version spatialisée fonctionnelle, pour finalement être retirée lorsque les gains escomptés n’ont pas été observés [8]. Dans l’application qui va suivre, nous allons nous intéresser à une solution de traitement exclusivement logicielle, offrant ainsi la souplesse du prototypage rapide et la flexibilité de s’adapter aux diverses applications mentionnées ci-dessus : l’approche du traitement de signaux radiofréquence défini par logiciel (SDR) démontrera une fois de plus toute sa puissance.

Le traitement de signaux issus de la constellation de satellites GPS nécessite de répondre à plusieurs questions :

1. les satellites sont en mouvement puisque situés sur une orbite entre 20000 et 21000 km d’altitude. Bien que les horloges atomiques embarquées dans les satellites soient extrêmement stables, le mouvement du satellite induit un décalage de fréquence de la porteuse par rapport à sa valeur nominale qu’il faut compenser au niveau du récepteur. Par ailleurs, l’oscillateur cadençant le récepteur de DVB-T est de piètre qualité : il faut que l’identification de la fréquence de la porteuse soit en mesure de compenser les fluctuations de fréquence de la copie locale de la fréquence de porteuse du fait des fluctuations de température (variations aléatoires) et l’écart à la fréquence nominale du fait de la médiocrité de l’assemblage du circuit d’oscillation sur le récepteur DVB-T (biais) ;
2. tous les satellites communiquent sur la même fréquence de porteuse : il faut implémenter un mécanisme pour identifier quel satellite a envoyé quel message, et ce malgré les autres messages transmis par les autres satellites sur la même porteuse. Le partage de la bande spectrale par un encodage différent par chaque satellite se nomme CDMA - Code Division Multiple Access - un concept très répandu en communication numérique en général ;
3. ayant identifié quelle information est émise par quel satellite et sur quel écart de la porteuse à la valeur nominale, il nous restera à extraire la charge informative encodée dans la phase du signal (Fig. 1). Le mélangeur en émission induit (bit 1) ou pas (bit 0) une rotation de phase de 180° de la porteuse : il s’agit d’une modulation BPSK [1] (Binary Phase Shift Keying) dans laquelle la phase du signal porte un OU-exclusif des bits d’identification (code) et de navigation (message utile).

De façon générale, pour N variables transmises par liaison radiofréquence, il faut N boucles d’asservissement pour retrouver chacun des paramètres de l’émetteur. En liaison radiofréquence, le choix de travailler sur une porteuse - donc de décaler le spectre par rapport à la bande de base (centrée sur 0) - répond aux soucis de partage du spectre radiofréquence et de miniaturisation des antennes permettant la conversion du signal électrique en onde électromagnétique. Nous avons donc au moins deux paramètres à retrouver, la fréquence de la porteuse et la fréquence cadençant le signal transmis, afin de retrouver les informations contenues dans la modulation de la porteuse. En effet, nous n’avons aucune raison de croire que les deux oscillateurs aux deux extrémités de la liaison se comportent de la même façon, et l’identification de la porteuse est une façon de resynchroniser les deux oscillateurs de l’émetteur et du récepteur. Dans le cas du GPS, nous aurons N=3 paramètres libres : la fréquence de la porteuse (entre le signal émis par le satellite et le récepteur au sol), le code identifiant le satellite émetteur, et la fréquence (ou la phase, puisque la fréquence est la dérivée de la phase) à laquelle ce code est transmis. Nous verrons plus bas que la principale subtilité tiendra au mode de modulation qui induit une exploitation de la phase pour deux fonctions - l’asservissement de l’oscillateur local sur la porteuse du signal reçu et l’extraction du signal de navigation émis par le satellite.

En conclusion, les problèmes qu’il nous faudra aborder sont résumés ci-dessous :

1. tous les satellites de la constellation émettent sur la même fréquence nominale (contrairement à la constellation GLONASS russe). Le récepteur radiofréquence enregistre donc un signal représentatif de la somme des contributions de tous les satellites visibles par l’antenne, dont il faut extraire la composante de chaque satellite individuel. Chaque satellite est identifié par un code unique : il s’agit d’une stratégie de CDMA (Fig. 2) ;
2. le signal reçu au sol est sous le bruit thermique, compte tenu de la bande passante d’acquisition nécessaire à traiter les signaux GPS. En effet, tout composant électronique est affecté d’un bruit intrinsèque lié à l’agitation thermique des électrons qui portent l’information électrique dans le conducteur : ce bruit de très grand encombrement spectral est intégré dans la bande de mesure. Il est donc d’autant plus important que la température est élevée et que la bande passante de mesure est importante. Comme il est courant de le retrouver en thermodynamique, la relation entre la température et la densité spectrale de bruit s’effectue au travers de la constante de Boltzman, k_B=1,38 x 10⁻²³ J/K, qui se traduit en passant à une unité logarithmique en 10 x log₁₀ (k_B· T)=−204 dBW (=-174 dBm) à l’ambiante (T=293 K). L’exploitation d’un signal sous le plancher de bruit thermique est quelque peu déroutante au premier abord, et signifie qu’il n’est pas possible d’observer « simplement » un signal GPS sur un analyseur de spectre sans un traitement spécifique. Nous verrons ainsi comment l’analyse sur une longue durée (1 ms) d’un code bien spécifique permettra d’accumuler de façon cohérente la puissance du signal émis par le satellite afin d’augmenter la puissance du signal et de diminuer le bruit thermique jusqu’à obtenir un rapport signal à bruit positif ;
3. enfin, tous ces traitements nécessitent d’extraire l’information contenue dans le signal reçu du satellite et pour ce faire, de s’affranchir de la porteuse radiofréquence. Or un satellite a une fâcheuse tendance à être en mouvement le long de son orbite et donc d’induire un décalage de fréquence variable au cours du temps : nous devons identifier ce décalage de fréquence et l’annuler pour extraire l’information contenue dans la modulation en phase (BPSK).

Lors de nos mesures de sensibilité visant à établir la puissance incidente nécessaire d’une porteuse autour de 1,5 GHz modulée en fréquence ou en amplitude pour atteindre un rapport signal à bruit en sortie donné (choisi arbitrairement à 10 dB), nous avons vu [1] que la limite de détection d’un récepteur DVB-T est de l’ordre de -95 dBm. Deux questions se posent donc : comment la puissance reçue d’un satellite émettant depuis une orbite située à plus de 20000 km du sol se compare avec cette limite de détection, et comment ces deux valeurs se comparent au bruit thermique intégré sur la bande passante de mesure de 2 MHz.

La puissance P_E émise par un satellite GPS est de l’ordre de 25 à 50 W, ou +14 à 17 dBW [9], avec une antenne de G_E=13 dBi de gain. Après avoir dispersé cette puissance sur une fraction de la sphère entourant le satellite de 21000 km de rayon - la distance entre l’émetteur (satellite) et le récepteur (nous, sur la surface de la Terre) - induisant une perte de 182 dB (Free Space Propagation Loss induisant une atténuation de (4π x d/λ)² par distribution uniforme de la puissance incidente sur la sphère de surface 4π·x d² à une distance d de l’émetteur, avec λ la longueur d’onde du signal transmis, 19 cm à 1,57542 GHz), la puissance reçue P_R par le récepteur est -155 dBW ou -125 dBm pour une antenne à peu près isotrope de gain G_R unitaire (Fig. 3).

Comment se compare cette puissance reçue avec la puissance thermique sur une bande passante de B=2 MHz ? La puissance du bruit thermique pour une antenne visant une région de température T est de k_B·x T·x B (avec k_B la constante de Boltzman) ou -174 dBm/Hz+10·log₁₀(2 MHz)=-111 dBm. Le bruit thermique est donc 14 dB au-dessus du signal GPS arrivant au sol, donc détecterons-nous l’information utile ?

Le signal transporté par la porteuse issue du GPS a la propriété de suivre un motif connu et distribué sur toute la bande passante de mesure. Il n’est évidemment pas possible de prétendre rompre avec la limite de capacité de canal de Shannon qui nous informe que la bande passante de transmission, C en bits/s, est liée à la bande passante du canal B en Hz et au rapport signal à bruit SNR (sans unité) de la liaison par C=B x·log₂(1+SNR) : augmenter le nombre d’interlocuteurs nécessite d’augmenter la bande passante, ici en ajoutant au message communiqué (les bits du message de navigation) un identifiant unique à chaque satellite, communiqué beaucoup plus rapidement que le message à transmettre. La recherche de ce motif d’identification nous permettra d’extraire une information pertinente du bruit.

Le code pseudo-aléatoire (voir encart 1 pour une description détaillée de ce point) codé sur 10 bits, généré au rythme de 1,023 MHz, se répète tous les 2¹⁰−1=1023 valeurs, ou 1 kHz. Le signal est donc intégré de façon cohérente sur une durée de 1 ms, induisant un gain de compression de 1,023·x 10⁶ x 10⁻³=1023 ou 30 dB. Le gain de compression a été abordé dans [10] : il s’agit de la propriété de l’intercorrélation d’accumuler de façon cohérente l’énergie du motif connu par lequel nous intercorrélons le signal acquis, tandis que le bruit se moyenne progressivement autour de 0 et n’induit pas cette accumulation cohérente. Intégrer sur 10 périodes consécutives du code recherché fait monter cette valeur à 40 dB. Par conséquent, le rapport signal à bruit après extraction de l’information pertinente est SNR=(−125+30)−(−111)=16 dB. Le signal sera donc clairement visible au-dessus du bruit après traitement.

Par ailleurs, est-ce que ces signaux sont au-dessus du plancher de bruit du récepteur (Fig. 4) ? L’antenne active équipant la majorité des récepteurs GPS intègre un préamplificateur faible bruit de +27 dB de gain. La puissance du signal reçu par le DVB-T est donc de -125+27=-98 dBm. Cette valeur est dangereusement proche du plancher de -95 dBm que nous avions mesuré [1], mais le plancher a été obtenu en vue d’induire un rapport signal à bruit de 10 dB après démodulation AM ou FM. Nous serons donc en mesure d’observer un signal exploitable.

Tableau 1 : Illustration de l’encodage de l’information par un signal pseudo-aléatoire connu.

Nous notons par ailleurs une propriété qui sera au cœur de nos analyses par la suite : la phase de l’intercorrélation reproduit l’encodage de phase sur le signal initial. La démonstration se fait simplement en considérant que l’intercorrélation du signal avec le motif connu est un processus linéaire : si nous nommons x(τ) l’intercorrélation du signal acquis s avec le code c en l’absence de codage de phase, alors

(noter que nous avons pris ici la définition de l’intercorrélation pour des éléments complexes pour laquelle nous utilisons le complexe conjugué d’un des termes - sans ce complexe conjugué nous obtiendrons bien une phase de l’intercorrélation dépendante de la phase du code, mais de signe opposé). Si maintenant nous ajoutons un terme de phase à la copie c' du code introduite dans le signal, tel que c'(τ)=c(t)· exp(jϕ) l’intercorrélation devient

donc nous retrouvons bien ϕ dans la phase de l’intercorrélation. Ainsi, si c -> c·exp(jϕ), alors x -> x·exp(jϕ). Le codage BPSK introduit par GPS sur le code se retrouvera dans la phase de l’intercorrélation.

La somme, version discrète de l’intégrale, de c par s s’accumule de façon cohérente lorsque s présente une séquence de c. Cette opération est répétée pour toutes les translations τ possibles de c dans s : xcorr(c,s)(τ)=∫_−∞^+∞ c(t)·x s*(t−τ)·x dt. Si le signal acquis par le récepteur radiofréquence est la somme de plusieurs contributions dont chacune est encodée par une signature c_k(t), alors il est souhaitable que xcorr(c_i,c_j)≃ 0 si i≠ j et xcorr(c_i,c_i)≃ 1. De cette façon, seule l’intercorrélation du signal s(t) par le « bon » code induira une intercorrélation non-nulle (section 3) : il s’agit de la condition d’orthogonalité qui sera discutée ci-dessous.

Dans le cas particulier du GPS, le code c(t) est généré par deux générateurs pseudo-aléatoires (Linear Feedback Shift Register - LFSR) imbriqués. La description détaillée de ce générateur de code - nommé par la suite PRN - n’a que peu d’intérêt, la seule information utile étant qu’il existe 37 déclinaisons de ces codes vérifiant les propriétés de base orthonormale citée ci-dessus, dont seules les 31 premières sont utilisées pour le moment.

2. Orthogonalité des codes

Le principe de l’extraction du code modulé en phase (BPSK) pour retrouver quel satellite a émis un signal sur la porteuse de fréquence commune à toute la constellation consiste à intercorréler le signal reçu par le récepteur radiofréquence avec chacune des séquences identifiant chaque satellite.

Le terme d’orthogonalité traduit que seule la corrélation d’un code avec une copie de la même séquence de bits génère un résultat non nul. Dans le cas du GPS, le code pseudo-aléatoire qui permet d’identifier la transmission issue de chaque satellite est issu d’un registre à décalage de 10 bits de long. Yann Guidon a déjà longuement expliqué les propriétés des LFSR [11] : GPS étend un peu la complexité en imbriquant deux LFSR. Nous nous sommes contentés, pour la génération des séquences pseudo-aléatoires, d’exploiter la routine Matlab (sous GNU/Octave) disponible à http://www.mathworks.com/matlabcentral/fileexchange/14670-gps-c-a-code-generator. Cependant, nous allons exploiter ce programme pour observer expérimentalement d’une part l’orthogonalité des codes, mais aussi leur robustesse à un décalage des fréquences d’échantillonnage (ou, ce qui revient au même, de la porteuse modulée par ces codes).

Le code est échantillonné à 1,023 MHz, et un LFSR de 10 bits de long présente une répétition tous les 2¹⁰−1=1023 bits. Ainsi, nous obtiendrons, si notre séquence est correcte, un maximum d’intercorrélation chaque fois que notre copie locale de la séquence s’accorde avec le signal reçu, donc toutes les 1 ms. Par ailleurs, l’intercorrélation porte l’information de phase qui a été imprimée par le modulateur sur la porteuse : l’exploitation de la phase de l’intercorrélation chaque fois que cet estimateur est maximum nous permettra de remonter aux symboles transmis. La différenciation des informations issues de chaque satellite s’obtient alors en jouant les diverses séquences des divers satellites visibles dans le ciel à un instant donné.

Une subtilité tient cependant au fait que le code pseudo-aléatoire module une porteuse. Or aucune contrainte ne force l’oscillateur local du récepteur de DVB-T à générer exactement la même fréquence que l’horloge atomique à bord du satellite du GPS. Au contraire, la fréquence émise par le satellite est décalée par effet Doppler - décalage variable en fonction de la position du satellite dans le ciel - et l’oscillateur local dérive avec son environnement (température, contrainte [1]). Ne connaissant pas a priori ce décalage de fréquence, nous devons chercher toutes les valeurs possibles afin d’identifier une première estimation grossière, afin d’affiner par la suite jusqu’à annuler la différence des fréquences par une boucle d’asservissement appropriée. Nous n’en sommes pas encore là : pour le moment, afin de minimiser le temps de calcul, nous devons estimer combien de pas de fréquence seront nécessaires pour la recherche. Plus l’excursion potentielle de l’écart des fréquences est importante, plus il faudra d’étapes de calcul pour un pas donné. D’un autre côté, si l’intercorrélation des codes pseudo-aléatoires est robuste à un certain décalage des fréquences entre oscillateur local et porteuse reçue, alors le pas de fréquence pourra être agrandi, afin de réduire le temps de calcul. La Fig. 5 propose le résultat de l’intercorrélation lorsqu'une des copies du code pseudo-aléatoire est échantillonnée 0,3 et 1% plus rapidement que l’autre. Nous constatons que 1% d’erreur se traduit par un effondrement et un étalement du pic d’intercorrélation - encore visible dans cet exemple idéal, mais inutilisable sur des signaux bruités réels - et un décalage de plus de quelques fractions de pourcent se traduira par un pic d’intercorrélation noyé dans le bruit. Nous avons empiriquement constaté qu’une recherche par pas de 100 à 500 Hz fournit un résultat acceptable. Compte tenu du décalage Doppler de ±5 kHz pour un récepteur fixe [1], cela fait tout de même 31 x 100=3100 calculs d’intercorrélation à effectuer pour trouver les satellites visibles dans le ciel (31 codes pseudo-aléatoires possibles, et 100 pas de fréquence). Ce calcul suppose par ailleurs que l’oscillateur local n’est pas biaisé, ce qui n’est pas le cas des DVB-T faible coût que nous utilisons. Nous avons observé des décalages de ±100 kHz (ou ±60 ppm à 1,5 GHz), ce qui étend considérablement la plage de recherche initiale en multipliant par 20 la durée du calcul. Ce dernier aspect justifie par exemple les recherches en cours sur des oscillateurs compacts d’exactitude améliorée (référencés sur une transition atomique [12]) afin de réduire le temps d’acquisition (nom communément fourni pour cette étape du calcul lorsque le récepteur GPS vient d’être allumé) des signaux GPS.

3. Principe de la démodulation

Nous avions discuté dans un article précédent [1] la nécessité de retrouver la porteuse du signal reçu afin de la compenser et donc extraire l’information codée dans la phase du signal acquis. Pour rappel, un signal s(t) est reçu par un récepteur radiofréquence, contenant une modulation de phase ϕ(t) tel que s(t)=sin(ω_R·x t+ϕ(t)), avec ω_R la pulsation du signal émis. La transposition de fréquence par mélange du signal reçu avec un oscillateur local ω_L suivi d’un filtre passe bas génère le signal numérisé par convertisseur analogique-numérique rapide de la forme S(t)=sin((ω_R−ω_L)·x t+ϕ(t)). Nous ne pourrons extraire ϕ(t) qu’en annulant la contribution de la rotation de phase (ω_R−ω_L)·x t, ce qui se traduit par un asservissement de ω_L sur ω_R. En effet, ω_R varie du fait du décalage Doppler alors que le satellite évolue le long de son orbite, tandis que ω_L varie à cause des fluctuations au sol de la température ou de la contrainte appliquée sur le résonateur à quartz. Finalement, il est classique de considérer un mélange non seulement par sin(ω_L·x t) qui fournit une estimation de cos(ϕ(t)) lorsque ω_R=ω_L, mais aussi par sin(ω_L·x t+π/2)=cos(ω_L·x t) afin de générer, lors du mélange, sin(ϕ(t)) : ces deux composantes sont appelées les composantes I et Q (en phase et en quadrature), utiles pour déduire la phase ϕ(t)=arctan(Q,I).

Tout ceci fonctionne fort bien pour une unique porteuse modulée en phase : nous avions présenté auparavant comment une astuce de traitement du signal - la boucle de Costas - permet de s’affranchir de la modulation et ainsi d’asservir ω_L sur ω_R afin d’extraire ϕ(t). Pour rappel, dans le cas particulier de la modulation à deux états BPSK, ϕ(t)∈[0;π] et 2ϕ=0 [2π] donc la mise au carré du signal reçu élimine la modulation en phase.

Le problème est plus complexe en CDMA lorsque la somme des contributions de plusieurs émetteurs est superposée et acquise. Dans ce cas, nous devons identifier la composante de chaque émetteur et travailler sur une démodulation de chacun de ces signaux individuellement. Le secret de cette stratégie de détection tient dans l’intercorrélation du signal reçu avec une copie locale du code identifiant chaque émetteur (Fig. 5). Grâce à la propriété d’orthogonalité des codes, seul le signal issu du « bon » satellite induit un maximum significatif de la magnitude de l’intercorrélation avec le « bon » code (Fig. 7).

4. Acquisition du signal

Commençons par le commencement : nous acquérons un signal dans la bande de fréquence du GPS, sans aucune information de quel satellite est visible ni de la position de la constellation dans le ciel. Nous devons donc rechercher de façon exhaustive toutes les possibilités pour trouver l’identifiant des satellites visibles, et le décalage Doppler induit par leur position le long de leur orbite.

Lors de la recherche exhaustive des fréquences permettant de compenser le décalage Doppler, l’expérience montre que le pas de balayage en fréquence doit être au plus de 100 Hz, faute de quoi le pic d’intercorrélation n’apparaîtra pas, car la modulation en phase est cachée dans la dérive de phase liée à l’écart de l’oscillateur local à la porteuse du signal reçu. Nous avons vu (section 2) que pour un récepteur fixe, le décalage Doppler est au plus de ±5 kHz, donc il s’agit de balayer 100 points de fréquences et 31 PRN possibles. Cependant, nous devons en plus tenir compte de l’énorme biais des oscillateurs locaux des récepteurs de télévision numérique terrestre : nous avons observé des écarts de ±60 ppm entre la valeur nominale de l’oscillateur et la valeur observée, une valeur médiocre, mais peu surprenante pour de l’électronique grand public. Un tel décalage se traduit par une recherche sur une plage de ±100 kHz, qui multiplie par 20 le temps de recherche lors de l’étalonnage du récepteur. Cette opération n’est cependant nécessaire qu’une fois par récepteur DVB-T, et on se contentera ultérieurement de décaler la plage de fréquences analysées du biais observé (Fig. 8).

Nous illustrons l’identification des satellites sur une longue séquence en effectuant l’acquisition automatique, toutes les 5 minutes, d’une seconde de signal dans la bande de fréquence du GPS suivi d’une recherche des satellites présents. L’automatisation de l’acquisition s’obtient en remplaçant l’enregistrement depuis gnuradio-companion par l’utilisation de rtl-sdr qui prend en argument la fréquence de l’oscillateur local (i.e. de quelle fréquence sera transposé le signal reçu sur l’antenne), la fréquence d’échantillonnage et la durée d’écoute. En pratique, cela donne while true; do nom=`date +%s`;echo $nom; rtl_sdr -f 1575420000 -s 2000000 -g 42 -n 200000 ${nom}.bin ;sleep 300;done.

Le fichier résultant de chaque acquisition contient des données encodées sur 8 bits, signées, interlaçant coefficients I et Q. Sous GNU/Octave, les données sont lues par f=fopen('1426853000.bin');x=fread(f,inf,'uchar');zz=x(1:2:end)-127+i*(x(2:2:end)-127);, le reste du traitement étant le même que celui proposé auparavant. Un film présentant l’évolution temporelle, sur une durée totale de 54 heures, des satellites visibles illustre le défilement puisque l’altitude des orbites est inférieure à l’orbite géostationnaire et la succession des divers satellites visibles depuis un point donné sur Terre. Le volume de données est conséquent : chaque fichier contenant 1 seconde d’enregistrement à 2 MS/s occupe 4 MB, soit un volume total pour les 54 h d’enregistrement de l’ordre de 2 GB. Le film résultant de ce traitement est disponible à http://jmfriedt.free.fr/gps_movie.avi.

Un exemple de traitement systématique de fichiers acquis par la méthode décrite ci-dessus est proposé. Dans un premier temps, nous chargeons tous les fichiers, retirons la valeur moyenne, définissons la gamme de décalage fréquence (dépendant du récepteur - ici un récepteur avec un biais de l’ordre de -95 kHz auquel nous ajoutons les ±10 kHz de décalage Doppler et de dérive thermique) ainsi que l’axe du temps, dont le pas est donné par l’inverse de la fréquence d’échantillonnage fs

Ayant chargé chaque fichier, nous allons générer une copie locale de chaque code PRN numéro m dans a, le décaler de la fréquence freq par multiplication du signal acquis par un signal périodique de phase freq*time, puis effectuons l’intercorrélation pour recherche le motif ainsi synthétisé dans le signal acquis

Finalement, pour chaque fichier un graphique proposant la magnitude de l’intercorrélation en fonction du code (abscisse) et du décalage en fréquence (Doppler + dérive thermique + biais) en ordonnée est affiché :

La séquence d’analyse se résume donc en une recherche exhaustive initiale de tous les codes de satellites possibles pour tous les décalages Doppler possibles : il s’agit de l’étape d’acquisition. Ensuite, connaissant quel satellite est visible, nous traiterons la phase de l’intercorrélation pour extraire la grandeur qui nous intéresse, à savoir la fréquence de la porteuse et le message de navigation. Comme la fréquence ne cesse de se décaler, au rythme de 10 kHz/12 h ou⁶ 0,23 Hz/s, nous devons sans arrêt recaler la fréquence de porteuse pour que l’encodage en phase du message de navigation soit exploitable : cette boucle d’asservissement est au cœur de l’étape de suivi du signal (tracking).

5. Suivi du signal

Ayant identifié grossièrement (nous verrons plus bas que « grossièrement » doit signifier une identification à une cinquantaine de Hz près la fréquence de la porteuse afin de permettre à la boucle de s’asservir) la fréquence de la porteuse et l’identifiant du satellite, nous allons remplacer la recherche exhaustive (f, PRN) par un asservissement sur la fréquence reçue. En effet, nous pourrions nous satisfaire de l’identification grossière de la porteuse et nous dire que nous pourrions maintenant simplement extraire la phase avec ses rotations de π pour obtenir le message (Fig. 9). Il n’en est rien : la phase varie bien trop vite pour permettre un décodage de l’information, puisque la phase varie comme Δ f· t avec Δ f le résiduel de l’écart de fréquence entre porteuse du signal incident et l’oscillateur local. À titre d’illustration, nous avons ajusté à la main la fréquence de l’oscillateur local afin que la phase présente une pente nulle au début ou à la fin d’une acquisition de quelques secondes (Fig. 9, gauche). Le décalage de phase entre le début et la fin de cette séquence est suffisant pour rendre l’extraction de l’information de navigation (au rythme de 50 bps ou 20 points issus de l’intercorrélation du PRN avec le signal reçu) inopérante. Il faut donc absolument asservir la phase de l’oscillateur local sur la phase de la porteuse reçue pour extraire l’information pertinente du GPS, ou exploiter les propriétés de la porteuse.

La dernière subtilité de mise en œuvre tient sur le fait que la phase de la porteuse présente deux rôles : elle porte l’information transmise (saut de π pour marquer la transition d’un état de bit à l’autre), et doit permettre de servir de discriminateur si la fréquence de l’oscillateur local dévie de la fréquence de la porteuse reçue.

La solution à cette double fonctionnalité est d’exploiter un estimateur de phase qui soit insensible aux rotations de π (au lieu de 2π comme habituellement rencontré). Or il se trouve que la fonction arctan(Q/I) a perdu, dans le ratio de Q par I, l’information de signe de I et Q individuellement. Dans cet état, le résultat de arctan se trouve toujours dans la moitié droite du cercle trigonométrique, ou en d’autres termes est insensible aux rotations de 180 degrés. Cela tombe bien, c’est exactement ce que nous recherchons (Fig. 10).

Nous allons donc exploiter deux estimateurs de phase pour extraire l’information des coefficients I et Q : arctan(Q/I) pour asservir la fréquence de la porteuse, car insensible aux sauts de 180°, et atan2(Q,I) (dans la nomenclature Matlab ou GNU/Octave) qui couvre bien tout le cercle trigonométrique. Nous utiliserons cette seconde fonction pour extraire la phase encodant les bits du code de navigation et vérifier que notre traitement du signal est pertinent. D’autres discriminateurs insensibles aux rotations de phase de 180° sont connus [13].

Afin de réaliser une boucle à verrouillage de phase numérique, il faut réaliser une loi de commande qui compense, en variant la fréquence f de l’oscillateur numérique exp(2π j f t) par laquelle est multiplié le signal acquis, la phase aux maxima successifs d’intercorrélation entre le signal reçu et la copie locale du code (maxima qui se répètent toutes les millisecondes). Par ailleurs, comme le générateur de code n’a aucune raison d’être synchronisé avec le générateur du satellite - ne serait-ce que du fait du décalage Doppler qui affecte aussi le code pseudo-aléatoire, même si à moindre mesure que la porteuse (il existe un rapport de 1500 entre les deux fréquences et donc l’effet Doppler qui les affecte) - il faut recaler régulièrement le code sur l’information amenée par la porteuse. Ces deux boucles d’asservissement opèrent simultanément et présentent une difficulté de conception dont la Fig. 14 (gauche) illustre un des échecs avant d’aboutir à un résultat fonctionnel (Fig. 14, milieu et droite).

La première loi de commande contrôle donc la fréquence de sortie d’un oscillateur virtuel commandé en tension (VCO) afin de maintenir la phase de la porteuse reçue constante. La seconde loi, qui ne sera pas détaillée, maintient la position (phase) du code localement généré constante par rapport au signal reçu : cette seconde loi ne sera pas détaillée, car apparaît moins fondamentale au bon fonctionnement du décodage sur les durées de mesures (quelques secondes) qui nous intéresseront ici.

La loi de commande la plus naïve - dite correcteur proportionnel - consiste à affirmer que la fréquence de sortie du VCO est proportionnelle à l’écart entre la fréquence de la porteuse reçue et la fréquence de l’oscillateur local. Un écart de fréquence induit donc une fréquence, les deux termes sont de même unité et la constante les liant est sans unité. Une approche aussi naïve échoue dans la pratique compte tenu du rapport signal à bruit médiocre de nos signaux.

Afin de « lisser » la différence de fréquence observée entre porteuse du signal reçu et oscillateur local, l’approche classique consiste à effectuer une moyenne glissante, ou en d’autres termes une intégration, qui nous amène naturellement au correcteur proportionnel-intégral. L’intégrateur accumule les erreurs successives observées, faisant ainsi lentement dériver la fréquence de l’oscillateur local jusqu’à ce que la fréquence de la porteuse du signal reçu soit égalée. Dans ce cas, l’erreur devient nulle et par la suite, le terme proportionnel prendra le dessus pour maintenir la fréquence du VCO égale à la fréquence de porteuse. Noter que pour un décalage Doppler variant de ±4 kHz en 12 h, la variation de fréquence est de 0,2 Hz/s. Cela signifie que la phase effectue une rotation de 2π toutes les 5 secondes, interdisant tout décodage de la modulation BPSK en l’absence d’un asservissement fiable.

En pratique, la commande du VCO (qui détermine donc sa fréquence de sortie) u_k est déduite des erreurs ε_k−i passées entre fréquence de la porteuse du signal reçu et fréquence du VCO, et des valeurs passées de la commande u_k−i. La Fig. 11 indique l’analyse préliminaire (acquisition) pour un jeu de données pour identifier quel satellite est visible avec quel décalage de fréquence, et l’étude se poursuite sur la Fig. 12 pour illustrer l’influence de l’initialisation du suivi de fréquence à l’issue de la recherche exhaustive, mais grossière lors de l’acquisition. Compte tenu de la dépendance aux conditions initiales, une identification de la fréquence initiale à 50 Hz près semble souhaitable lors de la phase d’acquisition. Nous cherchons à asservir, dans le cas de cette boucle à verrouillage de phase, la phase de l’intercorrélation sur une consigne nulle donc le signal d’erreur ε_k est directement la phase de l’intercorrélation ϕ_k et ses valeurs passées ϕ_k−i.

Une approche classique en automatique consiste à déterminer le point de fonctionnement par une variation lente de commande u_k en la faisant croître ou décroître en fonction du signe et de la valeur de ε_k. On choisit donc u_k comme étant l’intégrale de ε_k, ainsi u_k=α ∑_i=0^k_εi. En écrivant les premiers termes, u₀=αε₀, u₁=αε₀ + αε₁ = u₀+αε₁, u₂=αε₀ + αε₁ + αε₂ = u₁+αε₂, ... la relation de récurrence u_k=u_k−1+αε_k devient évidente. Il reste à déterminer α que l’on choisit en général très petit dans un premier temps (donc la commande varie très lentement au risque, dans le cas présent, de ne pas atteindre l’accrochage de la boucle avant d’atteindre la fin du fichier de données). Ensuite, on augmente α jusqu’à obtenir de bonnes performances en termes de temps d’établissement, en l’occurrence en termes de temps d’accrochage de la boucle à verrouillage de phase (Fig. 13). Quelques tâtonnements permettent de déterminer une valeur de α correcte, en l’occurrence α=0,01. Cependant, on observe sur la Fig. 13 que la valeur moyenne de l’erreur de phase est de -0,29. Ce n’est pas intrinsèquement gênant, mais un autre type de correcteur permet d’éliminer cette erreur et donc de diminuer la probabilité de « décrocher » en cas de fort bruit : le double intégrateur, u_k est alors l’intégrale de l’intégrale de ε_k.

On pose alors

v_k = v_k−1 + βε_k

u_k = u_k−1 + v_k +αε_k

d’où l’on déduit que :

u_k = 2u_k−1 − u_k−2 +(α+β)ε_k − αε_k−1

Le système est alors stable et à erreur nulle en choisissant α petit comme précédemment et 0 <β << α.

La loi de commande s’écrit donc uk = 2u_k−1 + u_k + (α+β)ε_k − αε_k−1 dont il faut identifier les coefficients α, β. Nous avons choisi α=0,0053 et β=0,05α. α et β vont déterminer les propriétés de la loi de commande : un α trop grand va rendre l’asservissement trop rapide et risque de dépasser le point de consigne et donc faire échouer la boucle d’asservissement, tandis que α trop faible induira une convergence trop lente. β détermine la pondération relative des valeurs passées et courantes de la phase de l’intercorrélation. Afin d’analyser les effets de ces divers coefficients, la Fig. 13 fait par exemple varier le coefficient α de part et d’autre de la valeur optimale 0,0053 de 0,0070 à 0,0035.

Le bon fonctionnement de cet algorithme est démontré par l’affichage du message de navigation, issu de la phase calculée par la fonction atan2 (donc capable de distinguer les rotations de phase de π), tel que présenté sur les Figures 12, 13 et 14.

Concrètement, nous chargeons le fichier de données, retirons la valeur moyenne des coefficients I et Q et alternons les octets pour générer le signal zz à analyser :

Cette procédure de chargements d’octets non-signés est valable si l’acquisition s’est faite au moyen de rtl-sdr. Si le fichier de données est issu du File Sink de gnuradio-companion, alors les données sont sous forme de flottants (et non d’entiers codés sur 8 bits - d’où un fichier quatre fois plus volumineux) et se chargent au moyen de read_complex_binary fourni dans le répertoire gr-utils des sources de GNURadio.

Nous allons nous focaliser sur le satellite de code numéro 31 (qui a été identifié comme présent lors de la phase exploratoire précédente) et en considérant que le décalage entre oscillateur local et porteuse du signal reçu est de l’ordre de 93705 (position du maximum d’intercorrélation dans l’étape précédente)

L’axe des temps est défini par pas de l’inverse de la fréquence d’échantillonnage - ici 2 MHz - et divers paramètres (alpha, beta) de la loi de commande sont définis. Les coefficients de la loi de commande sont dépendants de la fréquence d’échantillonnage du signal acquis : après intercorrélation, nous obtenons une mesure toutes les millisecondes, donc la période d’échantillonnage est de 1 ms. La loi de commande est établie sous la forme d’une relation récursive reliant une séquence de coefficients de pondération des valeurs passées de la mesure (phase - e_k) et des coefficients de pondération des valeurs passées de la commande (fréquence de sortie du VCO - u_k).

Après diverses initialisations additionnelles - en particulier GNU/Octave et Matlab adorent connaître à l’avance la taille des matrices manipulées et ne pas allouer dynamiquement de mémoire - le code PRN du satellite que nous avons sélectionné (msatellite) est généré, et des valeurs passées de la commande et de la mesure sont initialisées.

Deux asservissements seront nécessaires lors de l’analyse de la séquence : la fréquence de la porteuse et la position du code pseudo-aléatoire dans la séquence acquise. Le premier point sera pris en charge par une analyse de la phase, modulo π, du signal reçu, tandis que la seconde sera issue d’une comparaison de versions locales du code pseudo-aléatoire en retard (l comme late), en avance (e comme early) ou courante (p comme présent) lors de l’intercorrélation. Le code se décalera petit à petit pour suivre la phase induite par le léger décalage entre notre horloge qui échantillonne le signal et l’horloge du satellite qui génère l’identifiant.

Ayant aligné la position du code, nous sommes désormais prêts à engager la boucle d’asservissement sur la fréquence suite à toutes ces initialisations, boucle qui vise à maintenir la fréquence de l’oscillateur local égale à la fréquence de la porteuse reçue du satellite (identifiée comme la phase, modulo π, du signal reçu), et d’autre part la réplique locale du code alignée sur le signal reçu du satellite. Lorsque ces deux boucles d’asservissement sont fonctionnelles, une estimation de la phase (cette fois modulo 2π et non π) contient l’état des bits du message de navigation - 0 et 1 représentés par les angles de 0 et π, dans la variable yyp :

Finalement, le résultat de cette analyse est affiché sous forme de phase du signal après mélange par le VCO asservi, la fréquence du VCO, et la phase du signal comportant l’information de navigation.

Bien qu’intégriste de l’utilisation des logiciels libres, force est de constater que l’outil propriétaire Matlab obtient le résultat considérablement plus vite que GNU/Octave. En effet, le traitement d’1,4 millions de points (700 ms de durée, puisque chaque ms contient 2000 points lorsque le signal radiofréquence a été échantillonné à 2 MS/s) prend 61 s sous Matlab (R2010) contre 371 s sous GNU/Octave (3.8.2).

Nous concluons cette étude de décodage du message de navigation du GPS en recherchant l’entête de chaque message, la seule partie constante qui se répète toutes les 6 secondes (début de transmission d’un nouveau mot). Le mot de synchronisation, nommé TLM, se compose de l’octet 10001011 [2, p.112], que nous recherchons dans la séquence de bits issue des traitements développés ci-dessus. Comme le code transmis par les satellites (CDMA) se répète toutes les millisecondes et que la transmission numérique se fait à 50 bits/seconde, nous améliorons la capacité d’identification de la nature de chaque bit en moyennant sur 20 valeurs successives après s’être calé sur la première transition nette de 0 à 1. Finalement, le motif TLM est recherché dans la séquence par intercorrélation [2, p.119] : le bon décodage de la séquence ne peut être validé que si le maximum d’intercorrélation se répète toutes les 6 secondes (en plus des occurrences aléatoires du motif dans le contenu du message). En décodant près de 15 secondes de télémétrie, nous retrouvons bien deux maxima d’intercorrélation aux dates 4,4 et 10,4 s, démontrant la validité de notre traitement. Évidemment, un code aussi court que TLM ne donne qu’un pic d’intercorrélation peu résolu par rapport au bruit : le maximum d’intercorrélation est de 2 (tlm=[1 0 0 0 1 0 1 1];sum((tlm-mean(tlm)).^2)) et les valeurs de l’intercorrélation se distribuent sur 8 valeurs, qui se réduisent à 4 lorsque nous considérons la valeur absolue de l’intecorrélation (l’intercorrélation vaut -2 dans le cas où nous avons inversé les états 0 et 1 des bits de navigation, i.e. avons interverti les état 0 et π de la phase du signal).

Le programme GNU/Octave pour obtenir le graphique de la Fig. 15 est détaillé ci-dessous :

tandis que le contenu de la variable sortie est analysé de façon détaillée (ici 3 phrases successives, les deux dernières étant annotées pour aider le lecteur) :

Le début de la séquence est le résidu de la trame précédente. Suivent le mot de synchronisation (10001011), la séquence réservée (16 bits) et les 6 bits de parité du premier mot de 30 bits, puis la séquence HOW dont seule l’ID - qui identifie la trame en cours de transmission (entre 1 et 5) - nous intéresse (notez néanmoins l’incrément d’une unité du compteur HOW-count). Nous avons bien un identifiant égal à 5 (101 : fin de la trame précédente) suivi, dans la trame suivante, de 1 (début d’une nouvelle trame). Les deux derniers bits de HOW (second mot, contenus dans la fin des 6 bits de parité) sont bien à 0 : tout ceci est cohérent avec la norme du message de navigation décrit dans http://www.losangeles.af.mil/shared/media/document/AFD-070803-059.pdf (pp.109 et suivantes). Pour la phrase d’identifiant 1, les 10 premiers bits indiquent la date de transmission modulo 1024. Nous avons 1100010101 soit la semaine 789 qui indique une acquisition dans la semaine entre le 5 et 11 octobre 2014 (cf http://adn.agi.com/GNSSWeb/), une valeur cohérente. Suivent les deux bits de statut des transmissions - 01 indique la nature du code émis sur le second canal sur la porteuse L2, puis la précision sur le positionnement que fournira ce satellite : la meilleure résolution possible est accessible, entre 0 et 2,4 m, puisque tous les systèmes fonctionnent de façon nominale (tous les bits à 0). Cette analyse nous rend donc confiants sur la pertinence du résultat atteint.

6. Extension au RADAR passif

Nous avons extrait le signal reçu d’un satellite GPS en particulier par intercorrélation avec une copie locale du code, information connue et faisant partie intégrante du protocole de communication. Une alternative à l’exploitation d’un code connu émis par un émetteur connu, en vue de réaliser un RADAR passif [14], est de rechercher dans le bruit électromagnétique environnant une information pseudo-aléatoire et effectuer le même traitement que nous venons d’effectuer, mais cette fois entre l’information transmise directement de l’émetteur (supposé pseudo-aléatoire) au récepteur DVB-T, et la même information ayant rebondi sur une cible mobile. Cette approche de RADAR passif - très intéressante, car ne nécessitant pas d’émettre volontairement un signal radiofréquence, mais exploitant de façon opportuniste les sources existantes - permet à l’amateur de créer son système de mesure RADAR sans enfreindre la loi en émettant un signal, ou au militaire de détecter des cibles sans risquer d’attirer l’attention avec un émetteur dédié. Il s’agit donc d’un mode de veille très à la mode puisqu’il s’affranchit d’une source active dédiée, qui est elle-même source d’irradiations peu favorables en milieu urbain, voire de vulnérabilités si la cible de ne veut pas être vue. Cette méthode de travail est à l’origine du tout premier RADAR dans lequel un émetteur de signal radiofréquence de la BBC a été la source du signal réfléchi par un avion pour démontrer la faisabilité du RADAR [15]. Plus tard, les signaux émis par les RADARs anglais ont été utilisés par des récepteurs allemands pendant la Seconde Guerre mondiale pour détecter les vols au-dessus du sud de la mer du Nord [15].

Dans ce cas, deux effets affectent le second signal reçu : un décalage Doppler si la cible est mobile, et un retard puisque la seconde onde a cheminé plus longtemps pour atteindre le récepteur que l’onde directe (Fig. 16). Une analyse recherchant, pour tous les décalages de fréquence Doppler possibles, la position du maximum d’intercorrélation représentatif du retard du motif reçu initialement par R₁ puis, après réflexion par la cible par R₂, est le principe du RADAR passif qui exploite le bruit électromagnétique environnant pour localiser des cibles. Cette fois, le diagramme bidimensionnel que nous avions tracé pour le GPS avec en abscisse le numéro de code du satellite (PRN) et en ordonnée le décalage en fréquence, devient un diagramme avec en ordonnée le décalage Doppler (donc de nouveau un décalage en fréquence), et en abscisse l’intercorrélation dont le maximum est représentatif de la différence de temps de vol (t₁+t₂)−t_r.

La principale difficulté expérimentale dans une telle approche multistatique de RADAR (un émetteur et plusieurs récepteurs, tous séparés spatialement) est de synchroniser les horloges, tel que décrit à http://kaira.sgo.fi/2013/09/16-dual-channel-coherent-digital.html, et plus récemment détaillé sur http://hackaday.com/2015/06/05/building-your-own-sdr-based-passive-radar-on-a-shoestring/. Dans le cas du RADAR passif, R₁ et R₂ doivent être cadencés par la même horloge afin de garantir la cohérence en phase des mesures. Si ces deux récepteurs sont spatialement proches, alors un même oscillateur alimente les deux récepteurs radiofréquences : c’est l’approche sélectionnée dans la majorité des documents décrivant des expériences à ce sujet [16, 17, 18, 19].

Il nous faut cependant quelques ordres de grandeur sur les paramètres de mesure pour voir si l’expérience est envisageable. Les deux paramètres que nous cherchons à observer sont le décalage Doppler dû au mouvement de la cible, et la distance de la cible observée comme retard entre le signal de référence et le signal réfléchi par la cible.

La résolution en vitesse est donnée par le pas de fréquence de l’oscillateur local avec lequel le signal incident est mélangé (afin d’annuler tout écart de fréquence entre l’enregistrement de référence et le signal reçu de la cible mobile). Le décalage Doppler induit, sur une porteuse à f, par une cible mobile de vitesse v, est environ f x·v/c avec c la célérité d’une onde électromagnétique dans le vide. Pour une porteuse à 1 GHz et un véhicule roulant à 100 km/h, le décalage est de 100 Hz (pour éventuellement croître à 1 kHz pour un avion volant à 1000 km/h), parfaitement dans la bande de mesure du DVB-T. Nous vérifions, pour le moment par la simulation (Fig. 17), que la détection de cibles est bien faisable de cette façon, et que les propriétés de distance et de vitesse sont contenues dans l’abscisse et l’ordonnée des graphiques ainsi générés.

La résolution en portée est donnée par l’inverse de la bande passante du signal acquis [20]. La bande passante d’acquisition B est de l’ordre de 2,4 MHz avec un récepteur DVB-T donc la résolution sur la distance de la cible est c/(2B) ou environ 60 m. De façon générale, il peut être utile de mentionner que la résolution spatiale de la mesure est liée au contenu spectral du signal émis, tandis que la résolution en vitesse (i.e. en décalage Doppler) à la durée de la mesure, donc à sa résolution spectrale. GPS s’est mis dans le cas idéal d’un contenu informatif pseudo-aléatoire, comprimant au mieux l’impulsion lors de l’intercorrélation et permettant une identification fine de la position du message, même en présence de bruit. Le cas plus classique d’un signal périodique, ou pseudo-périodique (par exemple de la musique si un signal sur la bande FM est utilisé) se traduira par une moins bonne résolution spatiale puisqu’une sinusoïde légèrement décalée dans le temps (donc dans l’espace) ressemble toujours à une sinusoïde (Fig. 18). Les cas pratiques sont entre ces deux cas extrêmes du signal pseudo-aléatoire et pseudo-périodique : nous prendrons soin, pour les applications de RADAR passif, de choisir des émetteurs puissants, rayonnant vers les cibles recherchées, un signal aussi aléatoire que possible et occupant une bande spectrale au moins aussi importante que la carte d’acquisition. La littérature sur le sujet est très abondante [21].

Conclusion

Malgré un intérêt de longue date pour le décodage de signaux GPS, ce n’est que récemment que les circonstances se sont agencées pour permettre de faire aboutir un projet nécessitant un récepteur radiofréquence performant et peu coûteux, offrant une bande passante suffisante pour décoder un signal transmis au Mb/s, et la puissance de calcul nécessaire à une implémentation logicielle des algorithmes de décodage. Cette conjoncture offre l’opportunité à tout lecteur curieux de reproduire ces expériences. La démonstration initiale de la faisabilité du projet a été introduite par la présentation « Hacking GPS » de P. Boven à la conférence OHM2013, qui y expliquait la stratégie de prendre le carré des coefficients I/Q pour retirer la modulation qui induit l’étalement de spectre, et ainsi accumuler de façon cohérente le signal issu d’un satellite tout en abaissant le bruit thermique. La mise en œuvre n’est pas si triviale que cet auteur a voulu nous faire croire, et en particulier il faut savoir à quel écart de fréquence rechercher le signal d’intérêt aux côtés de la forêt de pics parasites induits par les diverses sources de bruit polluant l’acquisition. Après maintes recherches, nous avons fini par obtenir un signal représentatif de l’émission radiofréquence d’un satellite GPS par cette méthode (Fig. 19). Il s’avère cependant que cette méthode, d’apparence simple initialement, fournit un résultat pauvre et complexe à analyser. Elle a cependant été le point de départ de cette investigation, et méritait d’apparaître comme la conclusion d’une étude démontrant une compréhension complète de la nature des signaux et de leur modulation.

Quelle suite à donner à un tel projet ? De nombreux points restent à élucider. Le premier point évident est que nous avons visualisé et extrait la signification des bits du message de navigation, mais nous n’avons pas poussé la démonstration jusqu’à exploiter les éphémérides des satellites, pour finalement se positionner dans l’espace. Cette étape apparaît finalement nécessaire même pour l’objectif moins ambitieux de reproduire au sol une source de fréquence aussi stable que la combinaison des horloges atomiques embarquées, puisqu’il semble que ce soit la seule stratégie pour s’affranchir du décalage Doppler qui pollue chacune des porteuses observées. Par ailleurs, l’utilisation de la forme du signal de corrélation pour extraire des signaux directs et réfléchis une distance précise entre récepteur et réflecteur (application au RADAR passif) apparaît complexe pour des mesures depuis le sol compte tenu de la bande passante réduite de la mesure. En effet, avec un code se répétant à 1,023 MHz, un intervalle de temps du corrélateur occupe 977 ns pendant lesquelles un signal électromagnétique se propage de 294 m. Le réflecteur doit donc se trouver à plusieurs fois cette distance pour que le pic secondaire d’intercorrélation - dû au signal réfléchi et non au signal direct - soit séparé et aisément détectable. Ceci explique que la majorité des applications exploitant le signal réfléchi du GPS pour une mesure de distance se soit faite avec des récepteurs aéroportés, et de plus au-dessus de la mer afin d’avoir une surface réfléchissante importante et donc améliorer le rapport signal à bruit. De tels domaines d’application dépassent les moyens accessibles à l’amateur. Finalement, nous n’avons pas eu de succès à décoder les informations transmises par les satellites géostationnaires d’aide à la navigation - EGNOS en Europe⁷ et WAAS au-dessus des Amériques - bien que leur code et leur protocole de communication soit censés être compatibles GPS. Les causes de cet échec ne sont pas connues.

Remerciements

Cette étude a été menée en vue des cours de l’École d’Été sur le Temps-Fréquence (EFTS) qui se tiendra à Besançon. Le soutien financier pour l’accueil d’une étudiante ERASMUS - Sarà Martinez Gutierrez - est fourni par le labex First-TF (http://first-tf.com/). Les références bibliographiques qui ne sont pas librement disponibles sur le web sont obtenues auprès de Library Genesis, lib.gen.rus.ec : nos recherches ne seraient pas possibles sans cette ressource documentaire. É. Carry a amélioré cette prose par sa relecture.

J.-M Friedt est enseignant-chercheur à l’Université de Franche-Comté à Besançon, consultant pour une société privée concevant des capteurs passifs interrogeables par liaison radiofréquence s’apparentant aux méthodes de mesure par RADAR. G. Cabodevila est enseignant-chercheur à l’École Nationale Supérieure de Mécanique et Microtechniques de Besançon. Son cours d’automatique, qui formalise en particulier les concepts introduits dans la section 5, est disponible à http://jmfriedt.free.fr/Gonzalo_cours1A.pdf. Tous deux sont hébergés par l’Institut FEMTO-ST.

Références

[1] J.-M Friedt, La réception de signaux venus de l’espace par récepteur de télévision numérique terrestre, Open Silicium n°13 (Déc. 2014/Jan-Fev 2015), disponible à http://jmfriedt.free.fr/sdr2.pdf

[2] K. Borre, D.M. Akos, N. Bertelsen, P. Rinder, & S.H. Jensen, A software-defined GPS and Galileo receiver - A single frequency approach, Springer (2007)

[3] E.D. Kaplan, Understanding GPS : Principles and Applications, 2nd Ed, Artech House Mobile Communications (2005)

[4] Programme et transparents de la Journée GNSS 2015 du CNES, http://geodesie.ign.fr/journee-gnss-science/programme/

[5] R.L. Beard & J.D. White, GPS Application to Time Transfer and Dissemination, GPS Solutions 3 (1), pp.17-25 (1999)

[6] C. Bruyninx, P. Defraigne & J.-M Sleewaegen, Time and Frequency Transfer using GPS Codes and Carrier Phases: Onsite Experiments, GPS Solutions 3 (2), pp.1-10 (1999)

[7] Un excellent ouvrage sur la genèse de la définition du temps sur des transitions atomiques en replacement de la définition astronomique, à lire absolument : T. Jones, Splitting the second - The Story of Atomic Time, Institute of Physics Publishing (2000)

[8] G. Lachapelle, GNSS Solutions: Atomic clocks on satellites and mitigating multipath, Inside GNSS (Sept. 2006), disponible à http://www.insidegnss.com/auto/Sept06GNSS_Solutions_secure.pdf

[9] http://gpsinformation.net/main/gpspower.htm, bien que cette valeur semble plutôt venir d’un calcul inverse puisque la spécification de GPS définit la puissance minimum que doit recevoir un récepteur sur Terre plutôt que la puissance émise par le satellite. Une telle approche est proposée dans http://www.unoosa.org/pdf/icg/2010/ICG5/wgA/05.pdf et donne effectivement environ 27 à 28 W (14 dBW - notez l’erreur d’unité probable dans le transparent 17). http://www.nxp.com/documents/other/75016740.pdf propose une puissance émise de 50 W.

[10] J.-M Friedt, Auto et intercorrélation, recherche de ressemblance dans les signaux : application à l’identification d’images floutées, GNU/Linux Magazine France n°139 (juin 2011), disponible à http://jmfriedt.free.fr/xcorr.pdf

[11] Y. Guidon, Sciences/Théorie de l’information : propriétés et dérivés des CRC et LFSR, GNU/Linux Magazine n°085 (juillet 2006)

[12] R. Boudot et al., A high-performance frequency stability compact CPT clock based on a Cs-Ne microcell, IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 59 (11), pp.2584-7 (2012) ou http://members.femto-st.fr/sites/femto-st.fr.laurent_larger/files/content/Homepage_ll_fichiers/pdf_files/FRGLS13/mo1700boudot.pdf

[13] C. O’Driscoll, M.G. Petovello & G. Lachapelle, Choosing the coherent integration time for Kalman filter-based carrier-phase tracking of GNSS signals, GPS Solutions 15, pp.345-356 (2011)

[14] C. Guyard, Voir sans être vu : le radar passif, La Recherche n°447 (décembre 2010), p.64, disponible à http://www.larecherche.fr/savoirs/technologie/voir-etre-radar-passif-01-12-2010-87441

[15] Y. Blanchard, Le radar, 1904-2004 : Histoire d’un siècle d’innovations techniques et opérationnelles, Ellipses (2004)

[16] A. Capria & al., Ship Detection with DVB-T Software Defined Passive Radar Proc. of the IEEE Gold Remote Sensing Conference (2010), disponible à http://ieee.uniparthenope.it/chapter/_private/proc10/6.pdf

[17] K. Chetty, G.E. Smith, & K. Woodbridge, Through-the-Wall Sensing of Personnel Using Passive Bistatic WiFi Radar at Standoff Distances, IEEE Transactions on Geoscience and Remote Sensing 50 (4), avril 2012

[18] P. Samczynski & al., A Concept of GSM-based Passive Radar for Vehicle Traffic Monitoring, Proc. Microwaves, Radar and Remote Sensing Symposium, (2011), pp.271-274

[19] F. Berizzi & al., USRP technology for multiband passive radar, Proc. IEEE Radar Conference (2010), pp.225-229

[20] J. M. Christiansen, DVB-T based Passive Bistatic Radar Simulated and experimental data analysis of range and Doppler walk, rapport du Norwegian University of Science and Technology (2009), disponible à http://www.diva-portal.org/smash/get/diva2:348777/FULLTEXT01.pdf

[21] P.E. Howland, D. Maksimiuk & G. Reitsma, FM radio based bistatic radar, IEE Proc.-Radar Sonar Navig., 152 (3), juin 2005