12 Wyznaczanie orientacji: czujniki, TRIAD i filtr · Matematyka i fizyka od podstaw

0 Prolog: Hubble na jednym żyroskopie

Kosmiczny Teleskop Hubble’a musi utrzymać kierunek z dokładnością do kilku tysięcznych sekundy kątowej przez wiele minut ekspozycji. Robi to, łącząc żyroskopy (mierzą obrót między pomiarami) z czujnikami prowadzenia, które śledzą gwiazdy w polu widzenia. Żyroskopy się zużywają: teleskop ma ich sześć, misje serwisowe je wymieniały, a od 2024 roku NASA prowadzi obserwacje w trybie jednego żyroskopu, żeby zachować sprawne w rezerwie. Teleskop dalej działa, ale wolniej się przestawia i nie może obserwować części nieba w danym momencie.

Ta historia zawiera całą lekcję w skrócie. Czujniki kierunku (gwiazdy, Słońce, pole magnetyczne) mówią, jak satelita jest ustawiony, ale rzadko, z szumem i nie zawsze. Żyroskopy mówią, jak się obraca, często i gładko, ale z dryfem. Dobra orientacja to połączenie jednych z drugimi. Ten sam problem rozwiązywaliśmy dla samolotu (dynamika lotu, lekcja 3 Fuzja czujników w 3D: estymacja orientacji na kwaternionach); w kosmosie nie ma pionu z akcelerometru, więc zaczniemy od tego, co jest zamiast niego.

1 Czujniki orientacji

Czujniki dzielą się na dwie grupy. Czujniki wektorowe mierzą kierunek czegoś znanego w osiach satelity: jednostkowy wektor $\mathbf b$ w osiach ciała, który odpowiada wektorowi $\mathbf r$ znanemu w układzie odniesienia (np. inercjalnym). Orientacja to macierz $A$ , dla której $\mathbf b = A\mathbf r$ . Czujniki prędkości kątowej (żyroskopy) mierzą $\boldsymbol\omega$ i pozwalają przenieść znaną orientację w czasie.

Typowe dokładności (rzędy wielkości, $1\sigma$ )
Czujnik	Co mierzy	Dokładność	Uwagi
Czujnik Słońca, zgrubny	kierunek do Słońca	ok. 1°	fotodiody na ściankach; nie działa w cieniu Ziemi
Czujnik Słońca, dokładny	kierunek do Słońca	0,01–0,1°	maska i matryca
Czujnik horyzontu	kierunek do środka Ziemi	0,05–0,3°	podczerwień; zakłócenia od atmosfery
Magnetometr	wektor pola magnetycznego	0,5–3° (z modelem pola)	ograniczony błędem modelu pola i polem własnym satelity
Tracker gwiazd	pełna orientacja	1–10″ w poprzek osi, kilka razy gorzej wokół osi	oślepia go Słońce, Ziemia, Księżyc
Żyroskop	prędkość kątowa	bias od 0,001 do kilku °/h	dryf: trzeba korygować

Magnetometr mierzy pole bardzo dokładnie, ale do orientacji potrzebny jest też kierunek pola w układzie odniesienia, a ten bierze się z modelu (IGRF). Błąd modelu rzędu 100 nT przy polu ok. 30 000 nT to $0{,}2^\circ$ ; do tego dochodzi pole wytwarzane przez sam satelitę. Tracker gwiazd to kamera z katalogiem gwiazd: rozpoznaje wzór gwiazd w polu widzenia i podaje od razu całą orientację, z dokładnością rzędu sekund kątowych.

Zadania · część 1

3 zadania · ok. 20 min

Model pola magnetycznego ma błąd 100 nT, a pole na orbicie ma 30 000 nT. Jaki błąd kierunku to daje w najgorszym przypadku?

Odpowiedź
Błąd prostopadły do pola obraca kierunek o $100/30000 = 3{,}3\cdot10^{-3}$ rad $= 0{,}19^\circ$ .
Tracker ma pole widzenia $20^\circ$ i matrycę 1024 piksele. Ile sekund kątowych przypada na piksel? Jak tracker osiąga lepszą dokładność?

Odpowiedź
$20\cdot 3600/1024 = 70$ ″ na piksel. Obraz gwiazdy jest celowo rozmyty na kilka pikseli, a jego środek (centroid) wyznacza się z dokładnością ok. 1/10 piksela (7″). Uśrednienie kilkunastu gwiazd poprawia to jeszcze ok. $\sqrt N$ razy.
Żyroskop ma bias $0{,}01^\circ$ /h, a satelita liczy orientację z samego żyroskopu przez godzinę. Jaki błąd narasta? A przy biasie $1^\circ$ /h (tani żyroskop MEMS)?

Odpowiedź
$0{,}01^\circ = 36$ ″ w godzinę; przy $1^\circ$ /h cały stopień. Dlatego bias trzeba szacować na bieżąco (część 4).

2 Orientacja z dwóch kierunków: TRIAD

Jeden kierunek nie wystarcza: obrót wokół samego kierunku do Słońca nie zmienia tego, co widzi czujnik Słońca. Dwa nierównoległe kierunki wyznaczają orientację jednoznacznie, a nawet z nadmiarem (cztery liczby na trzy niewiadome, bo kąt między wektorami jest taki sam w obu układach). Najprostszy algorytm, TRIAD (używany od lat 60.), buduje z obu par wektorów dwie ortonormalne trójki i składa je ze sobą.

TRIAD

\mathbf t_1 = \mathbf b_1,\quad \mathbf t_2 = \frac{\mathbf b_1\times\mathbf b_2}{|\mathbf b_1\times\mathbf b_2|},\quad \mathbf t_3 = \mathbf t_1\times\mathbf t_2;\qquad \text{tak samo } \mathbf s_1, \mathbf s_2, \mathbf s_3 \text{ z } \mathbf r_1, \mathbf r_2,

A = [\mathbf t_1\ \mathbf t_2\ \mathbf t_3]\,[\mathbf s_1\ \mathbf s_2\ \mathbf s_3]^{\mathsf T}.

Pierwszy wektor jest odtworzony dokładnie ( $A\mathbf r_1 = \mathbf b_1$ ); z drugiego bierze się tylko płaszczyzna. Dlatego pierwszym wektorem musi być ten dokładniejszy.

Kolejność ma znaczenie. Czujnik Słońca ( $0{,}05^\circ$ ) i magnetometr ( $1^\circ$ ), wektory pod kątem $90^\circ$ : ze Słońcem jako pierwszym wektorem RMS błędu wynosi $1{,}02^\circ$ , z magnetometrem jako pierwszym $1{,}41^\circ$ . Ponadto obrót wokół kierunku Słońca widzi tylko magnetometr, więc błąd wokół tej osi wynosi $\sigma_2/\sin\alpha$ , gdzie $\alpha$ to kąt między wektorami. Przy $20^\circ$ to już $2{,}9^\circ$ .

Zadania · część 2

3 zadania · ok. 30 min

W układzie odniesienia $\mathbf r_1 = [1;0;0]$ , $\mathbf r_2 = [0;1;0]$ ; w osiach ciała zmierzono $\mathbf b_1 = [0;1;0]$ , $\mathbf b_2 = [-1;0;0]$ . Policz $A$ algorytmem TRIAD i sprawdź wynik.

Odpowiedź
$\mathbf t = ([0;1;0],\ [0;0;1],\ [1;0;0])$ , $\mathbf s = ([1;0;0],\ [0;0;1],\ [0;-1;0])$ . $A = \sum_i\mathbf t_i\mathbf s_i^{\mathsf T} = \begin{bmatrix}0&-1&0\\1&0&0\\0&0&1\end{bmatrix}$ : obrót o $90^\circ$ wokół $z$ . Sprawdzenie: $A\mathbf r_1 = [0;1;0] = \mathbf b_1$ , $A\mathbf r_2 = [-1;0;0] = \mathbf b_2$ .
Dlaczego błąd obrotu wokół kierunku Słońca wynosi $\sigma_2/\sin\alpha$ ? Policz go dla $\sigma_2 = 1^\circ$ i $\alpha = 60^\circ$ .

Odpowiedź
Obrót o $\delta\phi$ wokół $\mathbf r_1$ przesuwa wektor $\mathbf r_2$ o $\delta\phi\sin\alpha$ (tylko składowa prostopadła do osi obrotu się porusza). Błąd pomiaru $\sigma_2$ odpowiada więc błędowi kąta $\sigma_2/\sin\alpha = 1{,}155^\circ$ .
W cieniu Ziemi czujnik Słońca nie działa. Co zostaje satelicie z magnetometrem i czujnikiem Słońca?

Odpowiedź
Jeden kierunek (pole magnetyczne) nie wyznacza orientacji. Trzeba propagować orientację żyroskopami albo modelem dynamiki do wyjścia z cienia, albo dodać trzeci czujnik (np. horyzontu). Pole magnetyczne zmienia jednak kierunek wzdłuż orbity, więc filtr z modelem ruchu wyciąga z niego informację o wszystkich trzech osiach w dłuższym czasie.

3 Problem Wahby i metoda q

TRIAD nie umie wykorzystać trzech wektorów ani zważyć dwóch o podobnej dokładności. Grace Wahba postawiła w 1965 roku problem ogólnie: znaleźć macierz obrotu, która najlepiej przekształca wszystkie wektory odniesienia w zmierzone, z wagami $a_i$ .

Problem Wahby i metoda q Davenporta

\min_A\ L(A) = \tfrac12\sum_i a_i\,|\mathbf b_i - A\mathbf r_i|^2 = \sum_i a_i - \operatorname{tr}(AB^{\mathsf T}),\qquad B = \sum_i a_i\mathbf b_i\mathbf r_i^{\mathsf T}.

Z kwaternionem $q = (q_0, \mathbf q_v)$ : $\operatorname{tr}(AB^{\mathsf T}) = q^{\mathsf T}Kq$ , gdzie

K = \begin{bmatrix}\operatorname{tr}B & \mathbf z^{\mathsf T}\\ \mathbf z & B + B^{\mathsf T} - \operatorname{tr}B\,I\end{bmatrix},\qquad \mathbf z = \sum_i a_i\,\mathbf b_i\times\mathbf r_i .

Tu $q$ to kwaternion obrotu ciało → odniesienie w konwencji z lekcji 1 Orientacja w 3D: macierz obrotu, kąty Eulera, oś i kąt, kwaterniony (Hamilton, skalar na początku), a $A = R(q)^{\mathsf T}$ ; w podręczniku Markleya i Crassidisa konwencja jest inna i znak $\mathbf z$ może się różnić. Optymalny kwaternion to wektor własny $K$ dla największej wartości własnej $\lambda_{\max}$ ; $L_{\min} = \sum a_i - \lambda_{\max}$ . Wagi: $a_i = 1/\sigma_i^2$ .

Szukanie maksimum formy kwadratowej na sferze jednostkowej to dokładnie problem wartości własnych (mnożnik Lagrange’a dla więzu $|q| = 1$ ). Algorytm QUEST (Shuster, 1979–1981) przyspiesza go: wie, że $\lambda_{\max}$ jest bliskie $\sum a_i$ , więc szuka go kilkoma krokami metody Newtona na wielomianie charakterystycznym zamiast pełnego rozkładu. Na komputerach pokładowych lat 80. to była różnica decydująca.

Czego lepiej się nie da? Dla małych błędów kowariancja błędu orientacji optymalnego estymatora to

P = \Big(\sum_i \frac{I - \mathbf b_i\mathbf b_i^{\mathsf T}}{\sigma_i^2}\Big)^{-1}.

Każdy wektor daje informację tylko o obrotach prostopadłych do siebie (macierz $I - \mathbf b\mathbf b^{\mathsf T}$ zeruje obrót wokół $\mathbf b$ ). Przy dwóch pomiarach po $0{,}5^\circ$ i kącie $60^\circ$ : granica z kowariancji to $0{,}890^\circ$ , metoda q osiąga $0{,}888^\circ$ (Monte Carlo, 6000 prób), TRIAD $0{,}955^\circ$ .

Zadania · część 3

3 zadania · ok. 30 min

Pokaż, że $L(A) = \sum a_i - \operatorname{tr}(AB^{\mathsf T})$ dla wektorów jednostkowych.

Odpowiedź
$|\mathbf b - A\mathbf r|^2 = 1 + 1 - 2\mathbf b^{\mathsf T}A\mathbf r$ (obrót zachowuje długość). Zatem $L = \sum a_i - \sum a_i\mathbf b_i^{\mathsf T}A\mathbf r_i$ , a $\mathbf b^{\mathsf T}A\mathbf r = \operatorname{tr}(A\mathbf r\mathbf b^{\mathsf T}) = \operatorname{tr}(A(\mathbf b\mathbf r^{\mathsf T})^{\mathsf T})$ . Suma daje $\operatorname{tr}(AB^{\mathsf T})$ .
Dwa prostopadłe wektory: $\mathbf b_1 = \hat x$ z $\sigma_1 = 0{,}05^\circ$ i $\mathbf b_2 = \hat y$ z $\sigma_2 = 1^\circ$ . Policz $P$ . Wokół której osi błąd jest największy?

Odpowiedź
$\sum (I - \mathbf b\mathbf b^{\mathsf T})/\sigma^2 = \operatorname{diag}(1/\sigma_2^2,\ 1/\sigma_1^2,\ 1/\sigma_1^2 + 1/\sigma_2^2)$ . Odchylenia: wokół $x$ (Słońce) $\sigma_2 = 1^\circ$ , wokół $y$ $\sigma_1 = 0{,}05^\circ$ , wokół $z$ $0{,}0499^\circ$ . Największy wokół kierunku dokładnego wektora, bo ten obrót widzi tylko gorszy czujnik.
Bez szumu, dla trzech pomiarów z wagami 1, jaka jest $\lambda_{\max}$ ? Co oznacza $\lambda_{\max}$ wyraźnie mniejsze od $\sum a_i$ ?

Odpowiedź
$\lambda_{\max} = 3$ ( $L_{\min} = 0$ ). Duża różnica oznacza, że pomiary są niespójne: zły pomiar, błąd modelu pola, pomylona gwiazda. To prosty test wiarygodności.

4 Tracker gwiazd i żyroskop: filtr

Metody statyczne biorą pomiary z jednej chwili. Filtr łączy je w czasie: między pomiarami trackera orientację niesie żyroskop, a każdy pomiar trackera poprawia orientację i bias żyroskopu. W satelitach standardem jest multiplikatywny filtr Kalmana (MEKF), który znasz z dynamiki lotu (lekcja 3 Fuzja czujników w 3D: estymacja orientacji na kwaternionach): stanem jest mały błąd orientacji i błąd biasu, sześć liczb. Dla małych błędów trzy osie są prawie niezależne, więc wiele daje analiza jednej osi.

Filtr jednej osi: stan (kąt, bias)

Żyroskop: $\omega_m = \omega + b + \eta_v$ , bias błądzi: $\dot b = \eta_u$ . Szum kątowy (ARW) $\sigma_v$ [rad/√s], błądzenie biasu (RRW) $\sigma_u$ [rad/s^3/2], tracker $\sigma_n$ co $\Delta t$ .

W stanie ustalonym, gdy $\sigma_v^2\Delta t \ll \sigma_n^2$ :

\sigma_\theta \approx \sqrt{\sigma_n\sigma_v}\;\Delta t^{1/4}.

Bez pomiarów (sam żyroskop): $\sigma_\theta^2(t) = \sigma_0^2 + \sigma_b^2t^2 + \sigma_v^2t + \sigma_u^2t^3/3$ .

Liczby (wartości z podręcznika Markleya i Crassidisa: $\sigma_v = \sqrt{10}\cdot10^{-7}$ rad/√s, czyli $0{,}0011^\circ/\sqrt{\text{h}}$ , $\sigma_u = \sqrt{10}\cdot10^{-10}$ rad/s^3/2, tracker $\sigma_n = 17$ µrad $= 3{,}5$ ″ co 1 s): filtr daje $0{,}49$ ″, siedem razy lepiej niż pojedynczy pomiar trackera. Przybliżenie daje $0{,}48$ ″. Przy pomiarach co 10 s: $0{,}90$ ″ (przybliżenie $0{,}85$ ″), co 60 s: $1{,}5$ ″. Jeśli tracker zgaśnie (Słońce w polu widzenia), po 10 min błąd rośnie do 2,2″, po godzinie do 12″, po trzech do 49″; na długich przerwach dominuje błądzenie biasu, które rośnie jak $t^{3/2}$ .

Zadania · część 4

3 zadania · ok. 30 min

Policz $\sigma_\theta$ z przybliżenia dla $\sigma_n = 17$ µrad, $\sigma_v = 3{,}16\cdot10^{-7}$ rad/√s, $\Delta t = 1$ s i $\Delta t = 10$ s.

Odpowiedź
$\sqrt{17\cdot10^{-6}\cdot 3{,}16\cdot10^{-7}} = 2{,}32\cdot10^{-6}$ rad $= 0{,}48$ ″; dla 10 s razy $10^{1/4} = 1{,}78$ : $0{,}85$ ″. Dokładna wartość z równania Riccatiego: $0{,}49$ ″ i $0{,}90$ ″.
Skąd przybliżenie? Rozważ sam kąt (bez biasu): między pomiarami wariancja rośnie o $q = \sigma_v^2\Delta t$ , pomiar ma wariancję $r = \sigma_n^2$ . Pokaż, że dla $q \ll r$ wariancja w stanie ustalonym to ok. $\sqrt{qr}$ .

Odpowiedź
Stan ustalony przed pomiarem: $P^- = P^+ + q$ , $P^+ = P^-r/(P^- + r)$ . Stąd $P^-(P^- + r) - P^-r = q(P^- + r)$ , czyli $(P^-)^2 = q(P^- + r)$ . Dla $P^- \ll r$ : $(P^-)^2 \approx qr$ , $P^- \approx \sqrt{qr} = \sigma_v\sigma_n\sqrt{\Delta t}$ , a $\sigma = \sqrt{P}$ daje wzór z ramki.
W godzinnej przerwie błąd rośnie do 11,8″. Z czego się składa, jeśli bias oszacowano z dokładnością $0{,}0021^\circ$ /h?

Odpowiedź
Bias: $0{,}0021^\circ/\text{h}\cdot 1\,\text{h} = 7{,}6$ ″; ARW: $\sigma_v\sqrt{3600} = 3{,}9$ ″; RRW: $\sigma_u\sqrt{3600^3/3} = 8{,}1$ ″. Suma kwadratów: $\sqrt{7{,}6^2 + 3{,}9^2 + 8{,}1^2 + 0{,}5^2} = 11{,}8$ ″.

5 Od czujników do systemu

Kilka rzeczy, których nie widać we wzorach. Tracker po włączeniu nie wie nic: rozpoznaje gwiazdy, porównując kąty między nimi (trójkąty gwiazd) z katalogiem, i podaje orientację w kilka sekund bez żadnej wiedzy wstępnej (tryb „zagubiony w kosmosie”). Tracker oślepia Słońce, a także jasna Ziemia i Księżyc, w promieniu kilkudziesięciu stopni od osi, więc satelity mają zwykle dwa albo trzy trackery patrzące w różne strony. Dwa trackery prostopadłe do siebie dają też dobrą dokładność we wszystkich osiach: obrót wokół osi jednego to obrót w poprzek drugiego. Wreszcie czujniki są zamontowane z błędami ustawienia rzędu minut kątowych, więc ich kalibracja na orbicie (np. przez dodatkowe stany w filtrze) jest częścią projektu.

Zadania · część 5

2 zadania · ok. 20 min

Dlaczego tracker mierzy obrót wokół swojej osi kilka razy gorzej niż w poprzek?

Odpowiedź
Obrót w poprzek osi przesuwa wszystkie gwiazdy w polu widzenia o ten sam kąt. Obrót wokół osi przesuwa gwiazdę oddaloną od osi o $\rho$ tylko o $\delta\phi\cdot\rho$ , a $\rho$ jest co najwyżej połową pola widzenia (np. $10^\circ = 0{,}17$ rad). Ten sam błąd położenia gwiazdy daje więc kilka razy większy błąd kąta wokół osi.
Satelita ma jeden tracker o $30^\circ$ kącie wykluczenia dla Słońca. Ile czasu na orbicie może on być oślepiony, jeśli patrzy prostopadle do płaszczyzny orbity, a Słońce leży w tej płaszczyźnie? A jeśli patrzy w zenit?

Odpowiedź
Prostopadle do orbity: kąt do Słońca stale ok. $90^\circ$ , tracker nie jest oślepiany (dlatego tak się je ustawia). W zenit: oś obraca się raz na obieg i przechodzi przez Słońce; kąt mniejszy niż $30^\circ$ przez $60/360$ obiegu, ok. 16 min z 98, zawsze po stronie dziennej orbity.

6 Metoda i przykłady rozwiązane

Metoda: budżet dokładności orientacji

Wymaganie w każdej osi i czujniki, które je mogą spełnić (tabela z części 1).
Statycznie: kowariancja $P = (\sum (I - \mathbf b\mathbf b^{\mathsf T})/\sigma^2)^{-1}$ ; uwaga na mały kąt między wektorami.
Z żyroskopem: $\sigma \approx \sqrt{\sigma_n\sigma_v}\Delta t^{1/4}$ w stanie ustalonym, wzrost w przerwach.
Geometria: kiedy czujniki są oślepione albo w cieniu; ile trwają przerwy.

Przykład 6.1 Mały satelita: Słońce i magnetometr

CubeSat ma czujnik Słońca i magnetometr, oba z dokładnością $0{,}5^\circ$ ; kąt między wektorami $60^\circ$ . Jaka jest dokładność orientacji z TRIAD i z metody q?

Granica z kowariancji: $0{,}89^\circ$ RMS całkowity.
Metoda q z wagami $1/\sigma^2$ : $0{,}888^\circ$ (Monte Carlo). TRIAD: $0{,}955^\circ$ , ok. 7% gorzej, bo zachowuje pierwszy wektor dokładnie i odrzuca część informacji z drugiego.

Wynik: ok. $1^\circ$ ; przy równych dokładnościach metoda q ma niewielką przewagę.

Sprawdzenie: przy $\sigma_1 \ll \sigma_2$ TRIAD z dokładnym wektorem jako pierwszym i metoda q dają to samo (1,13° i 1,13° przy $0{,}05^\circ$ / $1^\circ$ ) ✓.

Przykład 6.2 Satelita obserwacyjny: tracker i żyroskop

Wymaganie: $1$ ″ ( $1\sigma$ ) na osi. Tracker 3,5″ co 1 s, żyroskop ARW $0{,}0011^\circ/\sqrt{\text h}$ . Czy wystarczy? A z przerwą 10 min?

Stan ustalony: $0{,}49$ ″, spełnione z zapasem.
Przerwa 10 min: $2{,}2$ ″, niespełnione; przerwa 1 min: $0{,}7$ ″, spełnione.

Wynik: trzeba drugiego trackera albo ograniczyć przerwy do minut.

Sprawdzenie: w Symulacji 2 przy tych parametrach koperta $3\sigma$ w stanie ustalonym ma ok. 1,5″ ✓.

Przykład 6.3 Zły pomiar

Metoda q z trzema wektorami (wagi znormalizowane do sumy 1) daje $\lambda_{\max} = 0{,}9990$ . Typowy szum to ok. $0{,}5^\circ$ . Czy to normalne?

$L_{\min} = 1 - 0{,}999 = 10^{-3}$ . Każdy pomiar ma dwie składowe błędu o wariancji $\sigma^2$ , więc przed dopasowaniem $L \approx \frac12\sum a_i\cdot 2\sigma^2 = \sigma^2$ . Dopasowanie trzech kątów „zjada” trzy z sześciu składowych, więc oczekiwane $L_{\min} \approx \frac{6 - 3}{6}\sigma^2 = \frac12(0{,}0087)^2 = 3{,}8\cdot10^{-5}$ (Monte Carlo: $3{,}86\cdot10^{-5}$ ).
$L_{\min}$ jest ok. 26 razy większe od oczekiwanego: jeden z pomiarów jest najpewniej błędny (np. magnetometr zakłócony przez prąd w satelicie).

Wynik: odrzucić pomiar o największym residuum i policzyć ponownie.

Sprawdzenie: $L_{\min}$ jest sumą kwadratów residuów, więc przy poprawnych pomiarach ma rozkład zbliżony do $\chi^2$ (statystyka, matematyka) ✓.

7 Błędne intuicje, których trzeba się oduczyć

„Dwa dokładne czujniki zawsze dają dokładną orientację.”: Nie, jeśli wektory są prawie równoległe: obrót wokół wspólnego kierunku widać słabo, błąd rośnie jak $1/\sin\alpha$ .
„W TRIAD kolejność wektorów jest obojętna.”: Pierwszy wektor jest zachowany dokładnie; musi to być ten dokładniejszy.
„Filtr nie może być dokładniejszy niż jego najlepszy czujnik.”: Może: żyroskop wiąże kolejne pomiary trackera i jego szum się uśrednia (0,49″ z pomiarów po 3,5″).
„Magnetometr mierzy pole dokładnie, więc kierunek też.”: Kierunek w układzie odniesienia pochodzi z modelu pola, a ten ma błąd rzędu 0,1–0,5°; do tego pole własne satelity.
„Kwaternion można filtrować jak zwykły wektor czterech liczb.”: Więz $|q| = 1$ czyni taką kowariancję osobliwą; poprawny jest filtr błędu (MEKF), jak w dynamice lotu.

8 Zadania do przemyślenia

Te zadania wymagają więcej myślenia niż liczenia.

Zadania do przemyślenia

5 zadań · ok. 1,5 godziny

Skąd wagi $1/\sigma^2$ ? Uzasadnij, że minimalizacja $L(A)$ z $a_i = 1/\sigma_i^2$ to estymacja największej wiarygodności przy gaussowskim szumie prostopadłym do wektora.
Rozwiązanie
1. Model: $\mathbf b_i = A\mathbf r_i + \boldsymbol\nu_i$ , $\boldsymbol\nu_i$ prostopadły do $\mathbf b_i$ , gaussowski z wariancją $\sigma_i^2$ na składową.
2. Logarytm wiarygodności to $-\sum |\mathbf b_i - A\mathbf r_i|^2/(2\sigma_i^2) + \text{stała}$ ; jego maksimum to minimum $L$ z wagami $1/\sigma_i^2$ (Shuster nazwał to modelem pomiaru QUEST).
Sam magnetometr. Czy satelita z samym magnetometrem i żyroskopami może wyznaczyć pełną orientację?
Rozwiązanie
1. W jednej chwili nie: obrót wokół kierunku pola jest niewidoczny.
2. W czasie tak: kierunek pola w układzie inercjalnym zmienia się wzdłuż orbity (na orbicie polarnej obraca się dwa razy na obieg), a żyroskopy łączą pomiary z różnych chwil. Filtr widzi kilka różnych kierunków, czyli dostaje informację jak od kilku czujników. Dokładność jest ograniczona błędem modelu pola.
Bez żyroskopów. Coraz więcej małych satelitów nie ma żyroskopów. Skąd biorą prędkość kątową i jaka jest cena?
Rozwiązanie
1. Z różnic kolejnych orientacji trackera (różniczkowanie) albo z modelu dynamiki w filtrze (równania Eulera ze znanymi momentami sił od kół).
2. Różniczkowanie wzmacnia szum ( $\sigma_\omega \sim \sqrt2\sigma_n/\Delta t$ ); model dynamiki wymaga dobrej znajomości momentów bezwładności i zakłóceń. Przerwa w pomiarach trackera jest groźniejsza, bo nie ma czym jej przetrwać.
Zbyt pewny filtr. Co się stanie, jeśli w filtrze ustawimy szum żyroskopu dziesięć razy za mały?
Rozwiązanie
1. Kowariancja filtru będzie za mała, wzmocnienie Kalmana też: filtr za bardzo ufa propagacji żyroskopem i za mało poprawia pomiarami.
2. Prawdziwy błąd wyjdzie poza kopertę $3\sigma$ ; przy dużych kątach linearyzacja przestaje działać i filtr może się rozbiec. Test: znormalizowane innowacje powinny mieć wariancję 1; za duża oznacza za mały szum modelu.
TRIAD a metoda q. Kiedy TRIAD jest równie dobry jak metoda q, a kiedy wyraźnie gorszy?
Rozwiązanie
1. Prawie równoważne, gdy jeden czujnik jest dużo dokładniejszy i ustawiony jako pierwszy: cała informacja o dwóch osiach pochodzi z niego, a o trzeciej z drugiego.
2. Gorszy przy podobnych dokładnościach (kilka procent), przy złej kolejności (dziesiątki procent) i przy więcej niż dwóch wektorach, których TRIAD nie umie użyć.

9 Sprawdź się

Odpowiedz w myślach, zanim rozwiniesz odpowiedź.

1 Dlaczego jeden kierunek nie wystarcza?

Obrót wokół tego kierunku go nie zmienia; potrzeba dwóch nierównoległych.

2 Który wektor dać jako pierwszy w TRIAD?

Dokładniejszy: pierwszy jest odtwarzany dokładnie.

3 Jak metoda q znajduje orientację?

Wektor własny macierzy $K$ dla największej wartości własnej to optymalny kwaternion.

4 Jaka jest najlepsza możliwa kowariancja z wektorów?

$P = (\sum (I - \mathbf b\mathbf b^{\mathsf T})/\sigma^2)^{-1}$ .

5 Jak dokładność filtru zależy od odstępu pomiarów?

Jak $\Delta t^{1/4}$ : $\sigma \approx \sqrt{\sigma_n\sigma_v}\Delta t^{1/4}$ .

6 Co dominuje w długiej przerwie w pomiarach?

Błąd biasu (liniowo w czasie) i błądzenie biasu (jak $t^{3/2}$ ).

10 Karta wzorów i co dalej

Metody statyczne

TRIAD: $A = [\mathbf t_1\ \mathbf t_2\ \mathbf t_3][\mathbf s_1\ \mathbf s_2\ \mathbf s_3]^{\mathsf T}$

Wahba: $L = \sum a_i - \operatorname{tr}(AB^{\mathsf T})$ , $B = \sum a_i\mathbf b_i\mathbf r_i^{\mathsf T}$

metoda q: $Kq = \lambda_{\max}q$

$P = (\sum (I - \mathbf b\mathbf b^{\mathsf T})/\sigma^2)^{-1}$

Filtr

$\omega_m = \omega + b + \eta_v$ , $\dot b = \eta_u$

$\sigma_\theta \approx \sqrt{\sigma_n\sigma_v}\,\Delta t^{1/4}$

przerwa: $\sigma^2 = \sigma_0^2 + \sigma_b^2t^2 + \sigma_v^2t + \sigma_u^2t^3/3$

tracker: wokół osi kilka razy gorzej

Do przeczytania

Obowiązkowo: F. L. Markley, J. L. Crassidis, Fundamentals of Spacecraft Attitude Determination and Control (Springer, 2014), rozdz. 4 (czujniki), 5 (statyczne wyznaczanie orientacji: TRIAD, metoda q, QUEST), 6 (filtracja, MEKF).
Dla chętnych: M. D. Shuster, S. D. Oh, „Three-axis attitude determination from vector observations”, Journal of Guidance and Control (1981): QUEST i model pomiaru w oryginale.

Satelita wie już, jak jest ustawiony, z dokładnością od stopnia do ułamka sekundy kątowej. Ostatnia lekcja kursu: jak go obrócić tam, gdzie trzeba, kołami reakcyjnymi, cewkami magnetycznymi i silniczkami, i co zrobić z momentem pędu, który koła zbierają od zakłóceń.