Samouczek 1 — Moja pierwsza orbita: gdzie są satelity GPS

Kind: orbit · Scenariusz: scenarios/orbit-sgp4-gps.toml (kopia dydaktyczna: scenarios/orbit.toml) · Poziom trudności: początkujący · ~15 min

Na końcu będziesz mieć propagowaną prawdziwą operacyjną konstelację GPS, odczytaną dostępność i geometrię oraz wyeksportowaną tabelę tego, gdzie każdy satelita faktycznie jest w przestrzeni — i będziesz wiedzieć, że ta tabela jest sprawdzona względem zewnętrznego wzorca, a nie względem samej siebie.

Czym jest ten scenariusz#

Satelity GPS mówią użytkownikowi, gdzie jest i która jest godzina. Zanim to nastąpi, trzeba wiedzieć, gdzie są satelity i ile z nich użytkownik widzi. To geometria, a geometria wyznacza dolną granicę dokładności wyznaczenia pozycji (współczynnik DOP, rozmycie dokładności). Ten samouczek:

bierze prawdziwą migawkę gps-ops z Celestrak konstelacji GPS (2021-07-28, 30 satelitów, prawdziwe elementy dwuliniowe, TLE),
propaguje każdego satelitę sprawdzonym względem zewnętrznego wzorca modelem SGP4/SDP4 (Simplified General Perturbations 4 / Simplified Deep-space Perturbations 4),
w każdym kroku czasu ustala, które satelity widzi użytkownik na niskiej orbicie okołoziemskiej (linia widzenia, zasłonięcie przez Ziemię i maska elewacji 5°),
zamienia tę geometrię na PDOP (rozmycie dokładności pozycji) i dokładność pozycji oraz
eksportuje faktyczne położenia satelitów w układzie związanym z Ziemią o środku w Ziemi (ECEF) do pliku SP3 (Standard Product 3, format precyzyjnych orbit Międzynarodowej Służby GNSS, IGS; GNSS to globalny system nawigacji satelitarnej).

Kręgosłup jest taki: TLE → propagacja SGP4 → ECEF dla każdej epoki → widoczność → liczba widocznych satelitów → PDOP → sigma pozycji = PDOP · σ_UERE, gdzie σ_UERE to równoważny błąd zasięgu użytkownika (błąd pomiaru odległości 1-σ dla jednego satelity).

Uruchom#

cargo run -- scenarios/orbit-sgp4-gps.toml

Python:

import json, kshana
result = json.loads(kshana.run(open("scenarios/orbit-sgp4-gps.toml").read()))
print(result["geometry"]["best_pdop"], result["geometry"]["best_position_sigma_m"])

Albo otwórz Studio w przeglądarce i wybierz scenariusz orbity.

Przeczytaj jednowierszowe podsumowanie#

Przebieg deterministyczny (seed = 17) drukuje:

scenario 4c51512369a6 | 345/361 samples GNSS-nominal | best PDOP 1.07 pos 1.07m | quantum holdover 0s p95 2.59e-5ns integrity 1.000 security n/a (no attack) | classical holdover 0s p95 6.7ns integrity 1.000 security n/a (no attack)

Pole po polu:

scenario 4c51512369a6 — 12-znakowy skrót scenariusza. Jest odciskiem dokładnych danych wejściowych (ziarno, progi, parametry modelu, blok TLE). Zmień dowolne wejście, a skrót się zmieni; zachowaj je bez zmian, a odtworzysz przebieg bit w bit. Ten sam skrót jest umieszczony w stopce wykresu i w pliku JSON (JavaScript Object Notation) wyniku.
345/361 samples GNSS-nominal — na 361 kroków czasu użytkownik ma użyteczne rozwiązanie w 345 z nich (95,6 %). Przebieg to przelot 12 h (43 200 s) z krokiem 120 s.
best PDOP 1.07 — najlepsze (najniższe) rozmycie dokładności pozycji w czasie przelotu. PDOP bliski 1 to doskonała geometria; im niższy, tym lepiej.
pos 1.07m — najlepsza dokładność pozycji: sigma pozycji = PDOP · σ_UERE. W tym scenariuszu sigma_uere_m = 1.0, więc te dwie liczby są tu liczbowo równe — ale reguła brzmi PDOP × UERE, a nie „pozycja = PDOP”. Przy UERE 3 m sigma pozycji wynosiłaby ~3,2 m. (Zob. pułapkę na końcu.)
Wartości zegara (holdover, p95, integrity, security) pojawiają się, ponieważ pakiet orbity niesie też zegar w czasie przelotu; są tematem Samouczka 2, a nie tego. (security pokazuje n/a (no attack), ponieważ nie skonfigurowano żadnego ataku.)

Blok geometrii w JSON zawiera liczby źródłowe:

"geometry": {
  "samples_total": 361,
  "samples_with_fix": 345,
  "sigma_uere_m": 1.0,
  "best_pdop": 1.0713624151805556,
  "median_pdop": 2.318376881311817,
  "best_position_sigma_m": 1.0713624151805556,
  "median_position_sigma_m": 2.318376881311817
}

Gdzie są satelity? (nagroda)#

Geometria jest dobra, ponieważ satelity naprawdę są w powłoce GPS. Udowodnij to — wyeksportuj ich położenia w układzie związanym z Ziemią:

cargo run -- scenarios/orbit-sgp4-gps.toml --export-sp3 gps.sp3

Pierwsze trzy satelity w pierwszym kroku próbkowania, ECEF w km (plik SP3 niesie czwartą kolumnę, zegara, zapisaną jako znacznik „brak wartości” 999999,999999):

PG01   9771.575585  24612.128830      0.003744 999999.999999
PG02  20449.643959 -11875.958365 -12379.457079 999999.999999
PG03  26467.367408     45.510223   4010.411755 999999.999999

Weź geocentryczny promień PG01:

|r(PG01)| = sqrt(9771.575585^2 + 24612.128830^2 + 0.003744^2) = 26480.9 km

Niecykliczna wyrocznia: czy to naprawdę orbita GPS?#

Samouczek, który sprawdza silnik względem niego samego, niczego nie dowodzi. Oto zewnętrzna, niecykliczna kotwica dla każdego powyższego twierdzenia.

1. Satelity leżą na prawdziwej powłoce średniej orbity okołoziemskiej (MEO) GPS. Nominalna półoś wielka GPS wynosi a = 26 560 km, wysokość ≈ 20 180 km nad średnim promieniem Ziemi 6 378 km (IS-GPS-200, specyfikacja interfejsu GPS; GPS Standard Positioning Service (SPS) Performance Standard, Departament Obrony USA (DoD); Misra & Enge, Global Positioning System, wyd. 2). Chwilowy promień PG01 równy 26 480,9 km mieści się w kilkuset km od a — dokładnie to, co daje rzeczywisty mimośród GPS (~0,005–0,02). Test tests/tutorials.rs::tutorial1_satellites_are_in_the_gps_meo_shell sprawdza 26000 km < |r(PG01)| < 27200 km względem tej opublikowanej wartości a, która jest zewnętrzna wobec Kshana.

2. Długość przelotu to jeden obrót GPS. Trzecie prawo Keplera daje okres T = 2π√(a³/μ) przy a = 26,560 km i μ = 398,600.4418 km³/s² (World Geodetic System 1984, WGS-84; raport techniczny National Imagery and Mapping Agency NIMA TR8350.2): T = 43,078 s = 11.967 h, pół doby gwiazdowej. Czas trwania scenariusza 43 200 s („~jeden obrót GPS”) zgadza się z tym z dokładnością lepszą niż 0,3 %, a ślad naziemny GPS powtarza się co dobę gwiazdową (dwa obroty), zgodnie z IS-GPS-200.

3. Propagowane położenia są sprawdzone względem zewnętrznego wzorca, a nie tylko wewnętrznie spójne. Propagator SGP4/SDP4, który wytworzył tę tabelę ECEF, zgadza się ze wszystkimi 666 wektorami weryfikacyjnymi AIAA (American Institute of Aeronautics and Astronautics) 2006-6753 z najgorszym przypadkiem 4,12 mm (tests/sgp4_verification.rs, docs/SGP4-VALIDATION.md) i zgadza się z niezależną biblioteką sgp4 poniżej mikrometra. Wyeksportowana tabela jest więc sprawdzona względem zewnętrznej implementacji wzorcowej, a nie względem własnej arytmetyki Kshana.

Uwaga o datach SP3. Nagłówek epoki SP3 pokazuje zastępczą datę kalendarzową (2000-01-01); propagacja faktycznie rusza od własnej epoki każdego TLE (2021). Czytaj tabelę jako „położenia ECEF w pierwszym kroku próbkowania”, a nie „w dniu 1 stycznia 2000”.

Co przypina test#

tests/tutorials.rs zamienia każdą powyższą liczbę w kontrakt CI:

tutorial1_orbit_headline_holds — best_pdop ≈ 1.071, best_position_sigma_m == best_pdop (ponieważ sigma_uere_m = 1.0), samples_total == 361 i podsumowanie zawiera 345/361.
tutorial1_satellites_are_in_the_gps_meo_shell — promień PG01 mieści się w paśmie powłoki GPS, sprawdzony względem opublikowanego a = 26,560 km.

Jeśli silnik kiedykolwiek przestanie wytwarzać prawdziwą powłokę GPS, kompilacja zrobi się czerwona.

Pułapki i jednostki#

Sigma pozycji to PDOP × UERE, a nie PDOP. Obie są równe tylko dlatego, że tu sigma_uere_m = 1.0. Zmień UERE, a się rozdzielą.
Wartości synchronizacji czasu są w nanosekundach, wartość pozycji w metrach. Nie mieszaj ich.
holdover_s jest skwantowane siatką (dolna granica przy kroku siatki czasu). W tym przelocie użytkownik praktycznie zawsze ma rozwiązanie, więc holdover wynosi 0 s (nie ma przerwy do przetrwania) — to dobry przypadek.

Co dalej#

Zacieśnij lub rozluźnij maskę elewacji (mask_deg) i zobacz, jak zmienia się dostępność — to ćwiczenie poziomu 2.
Podmień na świeżą migawkę TLE za pomocą scripts/fetch_tles.sh (poziom 2, „użyj świeżych danych”).
Przejdź do Samouczka 2 — Holdover zegara, gdzie sygnał GNSS jest odebrany, a zegar pokładowy musi dryfować sam.