|
Algorytmy zamiany współrzędnych kartezjańskich na elipsoidalne |
|
W pracy przedstawiono dwa nowe sposoby bardzo precyzyjnego
przeliczenia współrzędnych prostokątnych na elipsoidalne (geodezyjne). Oba
polegają na rozwiązaniu równania 2 sin(ψ
– Ω) =
c sin2ψ,
gdzie Ω = arctg [bz/(ar)], c = (a2
– b2)/[(ar)2 + (bz)2]1/2,
a i b są półosiami elipsoidy odniesienia, r i z — składową
rownikową i biegunową geocentrycznego wektora, a
ψ — tzw. parametryczną szerokością.
W ścisłym podejściu podstawienie
tg(π/4 –
ψ/2) ≡
t prowadzi do wielomianu: t4 + 2Et3 + 2Ft
– 1 = 0,
gdzie E = tgΩ
– c/cosΩ,
zaś F = tgΩ
+ c/cosΩ,
do którego można zastosować metodę Ferrari. Szerokość oraz
wysokość geodezyjną oblicza się ze wzorów φ
= arctg[a(1 – t2)/(2bt)] oraz
h = (r – at)cosφ +
(z – b)sinφ. Zarówno
to rozwiązanie ścisłe jak i rozwiązanie przybliżone (metodą
Newtona) dają porównywalnie wysokie dokładności — o kilka
rzędów wielkości lepiej niż w typowych algorytmach. |
|
Zagadnienie transformacji współrzędnych pomiędzy układem geocentrycznym a geodezyjnym występuje w codziennej praktyce w geodezji i astrometrii, a współczesne pomiary astrogeodezyjne wymagają zgodnej pary transformat, których dokładności utrzymują się w dużym przedziale wysokości na dowolnej szerokości geograficznej. O ile zamiana elipsoidalnych współrzędnych geodezyjnych, szerokości φ, długości λ i wysokości nad elipsoidą odniesienia h, na geocentryczne współrzędne kartezjańskie (x,y,z) lub biegunowe (szerokość, długość i promień wodzący) nie stanowi żadnego problemu, o tyle proces odwrotny sprawia pewne trudności. Odwrotne przekształcenie wykonuje się zwykle w sposób przybliżony w procesie iteracyjnym (np. [3], [8], [9]), [17]), albo używając przybliżeń w postaci szeregów potęgowych (np. [10], [12], [13], [16]) lub innych zwartych wyrażeń (np. [15], [2]). Dzieje się tak mimo, że istnieją ścisłe rozwiązania na zamianę współrzędnych kartezjańskich na geodezyjne ([6], [7], [14], [18]), o czym zdaje się nie wiedzieć wielu autorów poważnych publikacji (np. [1], [2], [8], [10], [16], [19]).
Rozwiązania ścisłe są w zasadzie wolne od błędów, jednak w praktyce błędy powstają wskutek zaokrągleń wyników pośrednich w stosunkowo złożonych algorytmach. Niemniej, w pewnych zastosowaniach te ścisłe algorytmy mogą być preferowane ze względu na ich dokładność i kompletność (zakresy współrzędnych).
W tym raporcie przedstawię dwa nowe rozwiązania, z których jedno jest bardzo prostym algorytmem przybliżającym rozwiązanie ścisłe tak dobrze, że w praktyce mu wcale nie ustępuje. Druga propozycja zawiera prawdziwie ścisłe rozwiązanie niezależne od znanych mi z literatury, a przewyższające je głównie swą prostotą. Prostota wynika z faktu, że rozwiązaniu poddaje się równanie czwartego stopnia, w którym tylko dwa spośrod pięciu współczynników różnią się od zera lub jedynki (nieliczni inni autorzy rozwiązują ogólną postać takiego równania). Podam także wyniki porównań moich algorytmów z kilku innymi.
Na początek przypomnijmy mniej lub bardziej znane wzory
zamiany współrzędnych geodezyjnych na kartezjańskie:
| (1a) |
| (1b) |
|
(2) |
| Rys. 1. Przykład określenia współrzędnych geodezyjnych (φ,h) dla elipsoidy o dużym spłaszczeniu [(a – b)/a = 0,23; parametr ten dla elipsoidy ziemskiej wynosi zaledwie ok. 0,003]. Współrzędne prostokątne r i z biegną od środka elipsy wzdłuż osi a i b, odpowiednio. D jest ujemne wewnątrz asteroidy (jej połowę zakreskowano w sposób zgodny z kierunkami szerokości geodezyjnych — tak, by pokazać że ich obwiednia wyznacza krzywą D = 0) i występują tu 4 rozwiązania rzeczywiste na (φ,h) każdego punktu (w kierunkach prostopadłych do elipsy: jeden w danej ćwiartce, trzy skierowane na przeciwległą półkulę). W innych miejscach (D ≥ 0) są dwa takie rozwiązania, z których jedno ma |φ| > 90° (tu: linia przerywana wzdłuż ujemnej wysokości). |
Częściej wzory (1) przedstawiane są równoważnie następująco:
|
(3a) |
| (3b) |
gdzie N = a/(1 – e2sin2 φ)1/2, zaś e2 = (a2 – b2)/a2.
Problemem jest znalezienie φ i h, gdy dane są r i z. Geometrycznie oznacza to wyznaczenie kierunku prostej prostopadłej do elipsoidy z danego punktu (r,z). Kąt, pod którym prosta taka przebija płaszczyznę równika, jest szerokością geodezyjną, φ, konwencjonalnie liczoną dodatnio na północ od równika, a ujemnie — na południe. Odległość punktu od elipsoidy jest wysokością, h (dodatnią na zewnątrz elipsoidy). Chociaż nie jest to oczywiste, ale Czytelnik używając elementarnych przekształceń może w wolnej chwili sprawdzić, że problem ten sprowadza się do rozwiązania równania ar tgψ = bz + (a2 – b2)sinψ albo:
| (4) |
gdzie Ω = arctg[bz/(ar)], zaś c = (a2 – b2)/[(ar)2 + (bz)2]1/2.
Taka reprezentacja omawianego problemu wydaje się być całkiem oryginalna i głównie dzięki przedstawieniu (4) nasze rozwiązanie okaże się o wiele dokładniejsze od innych. Wartość ψ uzyskaną z rozwiązania równania (4), na mocy równiania (2), można użyć do obliczenia szerokości geodezyjnej:
| (5) |
Ten wynik z kolei pozwala znaleźć wysokość nad elipsoidą ziemską z jednego albo obu równań (1), np.:
| (6) |
(W swoich obliczeniach, które przedstawiam niżej, stosowałem inne wyrażenie na h — to, które podałem w angielskim streszczeniu tej pracy).
Rozwiązanie równania typu (4) nie powinno w zasadzie sprawiać żadnych kłopotów (dziś nawet podręczne kalkulatory posiadają wbudowaną operację znajdowania zera funkcji przestępnych) niemniej, dla wygody Czytelnika, podam tutaj gotowy algorytm sprawdzony w praktyce. Rozwiązanie zaczynamy od obliczenia wygodnej wartości ψ′ — początkowej dla następnych przybliżeń ψ. Możemy w to miejsce wziąć np. arctg[az/(br)], które — jak to widać z (1) — jest ścisłe na powierzchni elipsoidy (dla h = 0), a więc będzie z pewnością optymalnym wyborem dla większości zastosowań praktycznych. Innego rodzaju wymierne korzyści (uproszczenie pierwszych wyrażeń w ciągu iteracyjnym) osiągniemy wybierając mniej optymalnie ψ′ = Ω. Zastosowanie metody Newtona na znajdowanie zera lewej strony równości (4) w każdym przypadku prowadzi do bardzo szybko zbieżnego ciągu przybliżeń. W metodzie tej korzysta się z pochodnej badanej funkcji względem danej zmiennej (u nas ψ). W naszym przypadku pochodna ta ma postać:
|
(7) |
Ocenę następnego przybliżenia dostaje się z:
| (8) |
gdzie ψ′ oznacza wartość początkową albo poprzednie przybliżenie. Z testów numerycznych wynika, że dwie iteracje dają dokładności o całe rzędy wielkości przewyższające potrzeby nawet bardzo wymagającego użytkownika. A są one rzędu 10–9 rad w szerokości geodezyjnej dla wszystkich punktów położonych dalej niż ok. 1000 km od środka Ziemi. W praktyce jest to jakieś sześć rzędów wielkości lepiej niż w przypadku algorytmów znanych mi z literatury (oczywiście wyłączając rozwiązania ścisłe).
W tabeli 1 zebrałem błędy pozycji obliczonych za pomocą 9 różnych algorytmów. Błędy są tam zdefiniowane jako odległości prawdziwych położeń od tych obliczonych. Pierwsze trzy algorytmy (oparte na pracach [8], [1] oraz [3], w takiej kolejności) są iteracyjnymi, a wyniki podane przed pochyłą kreską (/) dotyczą dwóch iteracji, zaś trzech iteracji — po tej kresce. Dalsze algorytmy zaczerpnąłem z prac [10], [12] (wzory 8 i 15 w tej pracy), [2] i [15], odpowiednio. Wzór Baranova i in. ([2]) poprawiłem na niewątpliwe przekłamanie w druku (parametr a w wyrażeniu na \tgí został w cytowanej pracy pomyłkowo podniesiony do kwadratu). Z tego zestawienia jasno wynika, że z naszym względnie prostym algorytmem, którego błędy zawiera przedostatnia kolumna tabeli, mogą konkurować jedynie rozwiązania ścisłe (ostatnia kolumna tabeli odpowiada mojemu rozwiązaniu przedstawianemu wcześniej w [4] oraz niżej w tej pracy). Chociaż nasze przybliżone rozwiązanie jest w istocie iteracyjnym to jednak, dzięki wyjątkowo szybkiej zbieżności, cały proces można ograniczyć do zaledwie dwóch wzorów, które odpowiadają pierwszemu przybliżeniu: ψ′ – Ω + 0,5 c sin2Ω/[1 – c cos2Ω] (wzór 8 z ψ′ = Ω) i drugiemu (wzór 8). Taką właśnie wersję użyłem do obliczeń odpowiednich składników tabeli 1.
|
Tabela 1 Błędy pozycji [mm] i czasy wykonania [0,1 ms] dziewięciu algorytmów. Ukośna kreska (/) oddziela wyniki dla dwóch i trzech iteracji, a gwiazdkami zaznaczono wartości większe od 10 000 mm. Oznaczenia algorytmów są zgodne z odpowiednimi publikacjami ze spisu literatury z wyjątkiem ostatnich dwóch kolumn, odnoszących się do niniejszej pracy (też [4] i [5]), z których pierwsza odpowiada algorytmowi rozwiązania przybliżonego, a druga — ścisłego. W algorytmie [3] rozwiązania dla h = 0 są nieokreślone. |
| |||||||||||||
Przedstawię teraz rozwiązanie ścisłe. Równanie (4) na
kilka sposobów można przekształcić do wielomianu czwartego
stopnia. Jeśli wyrazimy je w tg(π/4
– ψ/2) ≡ t,
to dostaniemy wyjątkowo prostą postać takiego wielomianu:
| (9) |
|
(10) |
|
(11) |
Istnieją standartowe rozwiązania równań czwartego stopnia (np. [11]), z których w naszym przypadku najodpowiedniejszym zdaje się być tzw. rozwiązanie Ferrariego. W już zredukowanej formie jest ono następujące:
| (12) |
| (13) |
|
(14a) |
| (15) |
| (16) |
| (17) |
Dla przypadku D < 0, aby uniknąć arytmetyki liczb zespolonych, można użyć równoważnego ale wygodniejszego od (14a) wyrażenia:
| (14b) |
Każde rozwiązanie rzeczywiste z czterech podanych (odpowiadających czterem kombinacjom znaków +/– w równaniach 12 i 13) reprezentuje inną parę współrzędnych geodezyjnych, a odpowiadające im kierunki są stycznymi do właściwych gałęzi ewoluty. Współrzędne te można odzyskać z t w następujący sposób:
| (18) |
| (19) |
Oczywiście, w praktyce nie są potrzebne wszystkie cztery rozwiązania (wszystkie one są rzeczywiste dla D < 0, a tylko dwa są takie w przeciwnym przypadku). Aby uzyskać wygodny do zaprogramowania algorytm dający to pojedyncze pożądane rozwiązanie wystarczy opuścić podwójne znaki (±) we wzorach (12) i (13) pozostawiając tam dodatnie wartości pirwiastków. Taki algorytm będzie poprawny tylko dla a > b (co w praktyce jest zawsze spełnione) i φ > 0. Także to ostatnie ograniczenie niewiele ujmuje z ogólności rozwiązania, gdyż na południowej półkuli mamy takie same wyniki, tyle że z przeciwnym znakiem szerokości. Ponadto, w naszym algorytmie istnieje prosty sposób automatycznego przypisania właściwego znaku wynikowi: nadanie małej półosi b znaku współrzędnej z przed procesem rozwiązywania problemu czyni podany algorytm (ten uproszczony do dodatnich pierwiastków) poprawnym na obu półkulach — bez konieczności wyszukiwania właściwego rozwiązania pośród czterech.
Przetestowałem ten algorytm (wzory od 10 do 19) w szerokim zakresie płaszczyzny (r,z) nie napotykając żadnych problemów. Również błędy zaokrągleń okazują się nieznaczące dla praktycznych zastosowań (używałem fortranowskiej opcji DOUBLE PRECISION), chociaż są one wyraźnie większe w pobliżu osi r i z niż gdzie indziej. Głównym źródłem tych błędów jest obliczanie wyrażenia (14a), które zawiera różnicę pierwiastków sześciennych (zaprogramowanych jako eksponent z trzeciej części logarytmu naturalnego). Problem ten obszedłem przez poprawienie wyniku otrzymanego z (14a), mówmy v′, korzystając z resolwenty równania (9) (której jednym z trzech pierwiastków jest właśnie v bądź v′):
| (20) |
W celu porównania efektywności omawianych algorytmów zmierzyłem czas potrzebny na ich wykonanie na IBM PC (z koprocesorem) z kompilatorem FORTRANu firmy Lahey Computer Systems. Wyniki, w jednostkach 0,1 ms podane w ostatnim wierszu tabeli 1, wskazują że wszystkie czasy są porównywalne mimo dość znacznych różnic w złożoności algorytmów. Wydaje się, że jest tak dlatego, że główny wkład do czasu wykonania pochodzi od obliczania funkcji trygonometrycznych, których w każdym programie jest po kilka. Na przykład, obliczenie funkcji kosinus zajmuje około 0,6 ms albo 6 jednostek tabeli 1. Ponieważ czasy obliczeń okazały się tak mało zróżnicowane, a obliczanie w sumie tanie, nie próbowałem optymalizować w tym względzie żadnych z implementowanych algorytmów. Dotyczy to także objętości poszczególnych programów, które są porównywalne (zajmują 7 – 15 linii FORTRANu) z wyjątkiem algorytmu [15], ale ten i tak zdaje się być najgorszy w każdym względzie. Wiedzieć jednak warto, że nasz algorytm ścisły będzie znacząco szybszy, gdy równanie (14a) zastąpimy przez równoważne:
| (14c) |
subroutine GEOD(r,z,fi,h)
c Program to transform Cartesian to geodetic coordinates
c Input: r, z = equatorial [m] and polar [m] components
c Output: fi, h = geodetic coord's (latitude [rad], height [m])
implicit real*8(a-h,o-z)
c IAU ellipsoid: semimajor axis and inverse flattening
data a,frec /6378140.d0,298.257d0/
b=dsign(a-a/frec,z)
if(r.eq.0d0) return
E=((z+b)*b/a-a)/r
F=((z-b)*b/a+a)/r
P=(E*F+1.)*4d0/3.
Q=(E*E-F*F)*2.
D=P*P*P+Q*Q
if(D.ge.0d0) then
s=dsqrt(D)+Q
s=dsign(dexp(dlog(dabs(s))/3d0),s)
v=P/s-s
v=-(Q+Q+v*v*v)/(3*P)
else
v=2.*dsqrt(-P)*dcos(dacos(Q/P/dsqrt(-P))/3.)
endif
G=.5*(E+dsqrt(E*E+v))
t=dsqrt(G*G+(F-v*G)/(G+G-E))-G
fi=datan((1.-t*t)*a/(2*b*t))
h=(r-a*t)*dcos(fi)+(z-b)*dsin(fi)
end
Test:
r z fi h
4000000. 6000000. 0.985526645027216 847786.688189974
4000.000 -6000.000 -1.48883906081174 -6350591.52477262
|
Podsumowując można stwierdzić, że algorytmy iteracyjne, takie jak [8], albo publikowane w roczniku astronomicznym [1], są najprostszymi do zaprogramowania i dają przyzwoite dokładności po 3 – 4 iteracjach. Z jeszcze jedną lub dwiema dodatkowymi iteracjami można ich bezpiecznie używać także, gdy jest wymagana wyższa dokładność. Jednak, z powodu błędów maszynowej reprezentacji liczb, próby z więcej niż 6 iteracjami mogą nie poprawiać już dokładności, która może być niekiedy nieco gorsza niż przy użyciu algorytmów opartych o rozwiązania ścisłe. W takich przypadkach bardziej właściwe wydają sie dwa algorytmy przedstawione tutaj albo algorytm z pracy [7]. Gwarantują one dokładności kilka rzędów wielkości poniżej 1 mm na całym praktycznie użytecznym obszarze płaszczyzny (r,z) — od kilkudziesięciu kilometrów od centrum Ziemi do wielu promieni ziemskich ponad jej powierzchnię. W tych zastosowaniach nasza poprawka (20) do matematycznie ścisłego rozwiązania może być bezpiecznie pominięta. W przypadku poszukiwania bardziej ogólnych rozwiązań na transformację współrzędnych prostokątnych na geodezyjne, jedynym gotowym algorytmem jest nasze ścisłe rozwiązanie (oryginalnie przedstawione pełniej, lecz z przykrym przekłamaniem, w pracy [4]).
Po złożeniu tego raportu w Redakcji ukazała się praca Laskowskiego [20],
który zwrócił uwagę na całkiem prosty algorytm iteracyjny Bowringa [21]
i stwierdził, że jest on równie dokładny jak nasz przybliżony, ale jest
nieco odeń szybszy. Ponadto w rozwiązaniu ścisłym odkryłem osobliwość
przy EF = –1, która występuje w odległości 41 – 46 km od środka Ziemi, jest więc niegroźna dla praktycznych zastosowań.
| LITERATURA |
[1] Astronomical Almanac for the Year 1989, USNO and RGO, Washington and London, str. K11.
[2] Baranov V.N. i in., Kosmicheskaya geodeziya, Nedra, Moskva, 1986, str. 27.
[3] Bartelme N., Meissl P., Ein einfaches, rasches und numerisch stabiles Verfahren zur Bestimmung des kürzesten Abstandes eines Punktes von einen sphäroidischen Rotationsellipsoid, AVN, Heft 12, 1975, 436 – 439 (Wichmann, Karlsruhe).
[4] Borkowski K.M., Transformation of Geocentric to Geodetic Coordinates without Approximations, Astrophys. Space Sci., 139, 1987, 1 – 4 (Erratum: 146, 1988, 201).
[5] Borkowski K.M., Accurate Algorithms to Transform Geocentric to Geodetic Coordinates, Bull. Géod., 63, 1989, 50 – 56.
[6] Hedgley D.R., An Exact Transformation from Geocentric to Geodetic Coordinates for Nonzero Altitudes, NASA Techn. Rep. R-458 (1976).
[7] Heikkinen M., Geschlossene Formeln zur Berechnung raumlicher geod„tischer Koordinaten aus rechtwinkligen Koordinaten, Zeitschrift Vermess., 107, 1982, 207 – 211.
[8] Heiskanen W.A. i Moritz H., Physical Geodesy, W.H. Freeman and Co., San Francisco, 1967, str. 181.
[9] Hlibowicki R. (red.), Geodezja wyższa i astronomia geodezyjna, PWN, Warszawa, 1981, str. 273.
[10] Izotov A.A. i in., Osnovy sputnikovoj geodezii, Nedra, Moskva, 1974, str. 34.
[11] Korn G.A. i Korn T.M., Mathematical Handbook for Scientists and Engineers, McGraw-Hill, New York, 1961, str. 24 (tłum. polskie: PWN, Warszawa, 1983, str. 36).
[12] Long S.A.T., General Altitude Transformation between Geocentric and Geodetic Coordinates, Celestial Mech., 12, 1975, 225 – 230.
[13] Morrison J. i Pines S., The Reduction from Geocentric to Geodetic Coordinates, Astron. J., 66, 1961, 15 – 16.
[14] Paul M.K., A Note on Computation of Geodetic Coordinates from Geocentric (Cartesian) Coordinates, Bull. Géod., Nouv. Ser., No 108, 1973, 135 – 139.
[15] Pick M., Closed Formulae for Transformation of the Cartesian Coordinate System into a System of Geodetic Coordinates, Studia geoph. et geod., 29, 1985, 112 – 119.
[16] Taff L.G., Celestial Mechanics. A Computational Guide for the Practitioner, J. Wiley and Sons, New York, 1985, rozdz. 3.
[17] Urmaev M.S., Orbital'nye metody kosmicheskoj geodezii, Nedra, Moskva, 1981, str. 14.
[18] Vanček P. i Krakiwski E.J., Geodesy: The Concepts, North Holland Publ. Co., Amsterdam, 1982, str. 324.
[19] Woolard E.W. i Clemence G.M., Spherical Astronomy, Academic Press, New York, 1966, rozdz. 3.
[20] Laskowski P., Is Newton's Iteration Faster than Simple Iteration for Transformation between Geocentric and Geodetic Coordinates?, Bull. Géod., 65, 1991, 14 – 17.
[21] Bowring B.R., Transformation from Spatial to Geographical Coordinates, Survey Review, XXIII, 181, 1976, 323 – 327.
Praca wpłynęła do Redakcji 1991.04.08
Akceptowana do druku 1992.10.20
|
Kazimierz M. Borkowski
Two very accurate closed solutions are proposed of which one is
approximate in nature and the other is exact. They are shown to be
superior to others, found in literature and in practice, in both or
either accuracy and/or simplicity. In the approximation algorithm
one determines the eccentric latitude: |