Jak znaleźć pierścień zasięgu satelity GPS na elipsoidzie WGS-84?

14

Biorąc pod uwagę następujące kwestie:

  1. Czas, t
  2. Zbiór danych efemeryd IS-200, E, satelity GPS odpowiadający czasowi t
  3. Pozycja ECEF satelity GPS, P = (x, y, z), wyprowadzona z czasu i efemeryd, (t, E).
  4. Załóżmy, że Ziemia to tylko elipsoida WGS-84.
  5. Wszystkie punkty na WGS-84 mają kąt maski, m.

Znajdź następujące:

  1. pierścień zasięgu, R, na WGS-84 satelity GPS. tj. granica, która odróżnia, które punkty WGS-84 są widoczne dla satelity w punkcie P = (x, y, z) i które punkty WGS-84 nie są widoczne

Koncepcyjna ilustracja problemu.  P jest czerwonym punktem, PRN12;  a czarny pierścień jest „pierścieniem zasięgu”

Dopuszczalne rozwiązania:

  1. Splajn nad WGS-84 zbliżony do R.
  2. Wielokąt nad WGS-84 zbliżony do R.
  3. Lub wzór (formuły), które dają mi R.

Co próbowałem do tej pory:

  • Niech e ^ 2 = 0,0066943799901264; ekscentryczność do kwadratu

Posiadamy pozycję ECEF WGS-84 według szerokości geograficznej phi i długości geograficznej lambda:

r = 1 / (sqrt (1-e ^ 2 sin ^ 2 (phi))) * (cos (phi) * cos (lambda), cos (phi) * sin (lambda), (1-e ^ 2) * sin (phi))

Następnie przekształcam ECEF w ramkę geograficzną wschód-północ w górę (ENU) za pomocą phi i lambda za pomocą macierzy:

     (-sin(lambda)                  cos(lambda)                  0       )
C=   (-cos(lambda)*sin(phi)        -sin(lambda)*sin(phi)         cos(phi))
     ( cos(lambda)*cos(phi)         sin(lambda)*cos(phi)         sin(phi))
  • Niech G = C (P - r)
  • Weźmy składnik Z G. Jeśli składnik Z G jest większy niż sin (m), to wiem, że punkt r jest widoczny. Ale to nie wystarczy uzyskać rozwiązanie, którego szukam. Mógłbym po prostu znaleźć kilka punktów, które są widoczne, i wziąć wypukły kadłub tych punktów, ale to wcale nie jest skuteczne.
torrho
źródło
1
Cześć @torrho, witamy w GIS.stackexchange. Bardziej prawdopodobne jest uzyskanie pomocy, jeśli pokażesz swoją pracę - to, czego próbowałeś do tej pory i co (konkretnie!) Sprawia ci kłopot.
Simbamangu,
@Simbamangu Jak korzystać ze znaczników lateksowych w GIS.stackexchange? czy mogę tylko powiedzieć $$ \ pi $$?
torrho
1
@tomfumb Nie, to nie jest praca domowa. Pomyślałem, że nie jestem jedynym, który napotkał ten problem, więc pomyślałem, że zapytam społeczność, która może mieć.
torrho
1
Widzę, że ktoś na [ meta.gis.stackexchange.com/questions/3423/… uważa, że ​​to zadanie domowe. To nie jest praca domowa, wyczerpująco przejrzałem ten temat i nie znalazłem nic na ten temat.
torrho
Niestety nie mogę znaleźć sposobu na użycie LaTeX na tej stronie! Możesz umieścić równania w tekście najlepiej, jak to możliwe, lub link do zrzutów ekranu równań LaTeX w innym miejscu (np. Folder Dropbox; nie możesz dodawać obrazów, dopóki nie uzyskasz wyższej reputacji). Poinformuj nas o kontekście tego problemu (dlaczego to robisz) oraz o tym, który konkretny komponent GIS sprawia mu kłopot, oraz jakie inne badania lub zasoby wykorzystałeś.
Simbamangu,

Odpowiedzi:

17

Rozwiązanie elipsoidy jest dość niechlujne - ma nieregularny kształt, a nie okrąg - i najlepiej jest je obliczać numerycznie niż za pomocą wzoru.

Na mapie świata różnica między rozwiązaniem WGS84 a rozwiązaniem czysto sferycznym będzie ledwo zauważalna (około jednego piksela na ekranie). Tę samą różnicę można uzyskać, zmieniając kąt maski o około 0,2 stopnia lub stosując przybliżenie wielokąta. Jeśli te błędy są akceptowalnie małe, możesz wykorzystać symetrię kuli, aby uzyskać prostą formułę.

Postać

Ta mapa (przy użyciu rzutu w kształcie prostokąta) pokazuje zasięg satelity w odległości 22 164 km (od środka ziemi) przy kącie maski m = 15 stopni na sferoidie WGS84. Ponowne obliczenie zasięgu kuli nie zmienia widocznie tej mapy.

Na kuli zasięg będzie naprawdę okręgiem wyśrodkowanym w miejscu satelity, więc musimy tylko określić jego promień, który jest kątem. Nazwij to t . W przekroju znajduje się trójkąt OSP utworzony przez środek ziemi (O), satelitę (S) i dowolny punkt (P) na okręgu:

  • Bok OP promień Ziemi, R .

  • Boczny system operacyjny to wysokość satelity (powyżej środka ziemi). Nazwij to h .

  • Kąt OPS wynosi 90 + m .

  • Kąt SOP to t , który chcemy znaleźć.

  • Ponieważ trzy kąty trójkąta sumują się do 180 stopni, trzeci kąt OSP musi wynosić 90 - ( m + t ).

Rozwiązanie jest teraz kwestią elementarnej trygonometrii. Potwierdza to (płaskie) prawo sinusów

sin(90 - (m+t)) / r = sin(90 + m) / h.

Rozwiązaniem jest

t = ArcCos(cos(m) / (h/r)) - m.

W ramach kontroli rozważ niektóre ekstremalne przypadki:

  1. Gdy m = 0, t = ArcCos (r / h), co można zweryfikować za pomocą podstawowej geometrii euklidesowej.

  2. Gdy h = r (satelita nie wystartował), t = ArcCos (cos (m) / 1) - m = m - m = 0.

  3. Gdy m = 90 stopni, t = ArcCos (0) - 90 = 90 - 90 = 0, jak powinno być.

Zmniejsza to problem rysowania koła na kuli, które można rozwiązać na wiele sposobów. Na przykład, możesz buforować lokalizację satelity za pomocą t * R * pi / 180, używając projekcji w równej odległości, wyśrodkowanej na satelicie. Techniki pracy z okręgami bezpośrednio na kuli są zilustrowane na /gis//a/53323/664 .


Edytować

FWIW, w przypadku satelitów GPS i małych kątów maski (mniejszych niż 20 stopni), to nie trygonometryczne przybliżenie jest dokładne (do kilku dziesiątych stopnia i mniej niż kilka setnych stopnia, gdy kąt maski jest mniejszy niż 10 stopni ):

t (degrees) = -0.0000152198628163333 * (-5.93410042925107*10^6 + 
              3.88800000000000*10^6 r/h + 65703.6145507725 m + 
              9.86960440108936 m^2 - 631.654681669719 r/h m^2)

Na przykład przy kącie maski m = 10 stopni i satelicie w odległości 26 559,7 km powyżej centrum Ziemi (co jest nominalną odległością satelity GPS ), to przybliżenie daje 66,32159 ... natomiast wartość (poprawna dla kuli ) to 66.32023 ...

(Przybliżenie oparte jest na rozszerzeniu szeregu Taylora wokół m = 0, r / h = 1/4.)

Whuber
źródło