MRAR-Cp
Č. úlohy 4
Transformace geodetických souřadnic
ZADÁNÍ 4.1
Sestavte v Matlabu funkci pro stanovení výšky geoidu WGS84.
4.2
Sestavte v Matlabu funkci pro přepočet geodetických souřadnic na kartézské souřadnice pro geodetický systém WGS84 s využitím funkce pro stanovení výšky geoidu z předchozího bodu.
4.3
Sestavte v Matlabu funkci pro přepočet kartézských souřadnic na geodetické souřadnice v systému WGS84 s aplikací prosté iterační metody a s využitím funkce pro stanovení výšky geoidu z bodu 4.1.
4.4
Otestujte funkčnost sestavených funkcí na příkladech.
ROZBOR Cílem cvičení je připravit si funkce pro výpočty určení polohy z GPS signálu, které jsou orientovány na přepočty souřadných systémů. Při určování polohy ze signálů GPS družic je výhodné (a z hlediska komplexnosti výpočtu prakticky nejsnadnější) pracovat v kartézském souřadném systému. Až po určení polohy v kartézském systému pak provést transformaci na všeobecně užívané geodetické souřadnice a připravit si funkce pro tyto transformace je náplní tohoto zadání. Výška geoidu pro daný geodetický referenční systém je stěžejní pro kalibraci určení skutečné nadmořské výšky určované polohy. Geoid je fyzikální model povrchu Země při střední hladině světových oceánů a je definován jako ekvipotenciální plocha vůči gravitaci, tzn., že je to plocha se stejnou úrovní tíhového potenciálu, na kterou je vektor tíhového zrychlení kolmý. Geoid se vůči referenčnímu zemskému elipsoidu může lišit až o ± 100 m. Známe-li výšku geoidu N pro daný bod měření, pak nadmořská výška HASL daného bodu je:
H ASL = H − N ,
(4.1)
přičemž H reprezentuje stanovenou výšku na referenčním elipsoidem. Výška geoidu je vždy stanovena k příslušnému geodetickému referenčnímu elipsoidu jako funkce závislá na geodetické šířce a délce. V tomto úkolu budeme používat systém WGS84. Při popisu geodetických přesných zeměpisných souřadnic vycházíme z definice referenčního elipsoidu a definujeme: •
ZEMĚPISNOU GEODETICKOU ŠÍŘKU ϕ – úhel svírající rovina rovníku s normálou k ploše elipsoidu (kladná na severní polovině zemského elipsoidu)
•
ZEMĚPISNOU GEODETICKOU DÉLKU λ – úhel svírající rovina místního poledníku s rovinou základního poledníku (kladná východním směrem)
•
ELIPSOIDICKOU VÝŠKU H – vzdálenost od elipsoidu, měřená po normále (kladná vně elipsoidu), kterou při aplikaci výšky geoidu přepočteme na nadmořskou výšku (viz vztah 4.1)
Na obrázku 4.1 je prezentována graficky definice geodetický souřadnic. Jednotlivé geodetické souřadné systémy se liší typem použitého referenčního elipsoidu, který je definován velkou poloosou a a zploštěním f, které je dáno výrazem
f =
a −b a ,
(4.2)
b je malá poloosa. Pro námi uvažovaný systém WGS84 je a = 6378137 m a 1/f 298,257223563.
=
Obr. 4.1 Geodetické souřadnice a definice kartézských souřadnic. Mezi pravoúhlými a geodetickými souřadnicemi platí vztahy:
p y = (ρ + H )cos ϕ sin λ p x = (ρ + H )cos ϕ cos λ
[(
)
]
p z = 1 − e ρ + H sin ϕ 2
,
kde e je excentricita (výstřednost) elipsoidu a ρ je příčný poloměr křivosti elipsoidu:
(4.3)
b2 e = 1− 2 a
,
(4.4)
a
ρ=
1 − e 2 sin 2 ϕ
.
(4.5)
Přepočet kartézských na geodetické souřadnice je o poznání komplikovanější. Je-li bod p (z obr. 4.1) kolmým průmětem bodu P do roviny geodetického rovníku pak jeho vzdálenost od počátku dp bude dána:
dp =
p x2 + p y2
,
(4.6)
a pro geodetickou délku platí vztahy:
sin λ =
py
dp p cos λ = x dp
,
(4.7)
z nichž můžeme jednoznačně určit geodetickou délku v celém jejím intervalu:
⎛
⎞ ⎟ ⎟. ⎝ px + d p ⎠ py
λ = 2arctg ⎜⎜
(4.8)
Eliminací geodetické délky získáme pro šířku a výšku soustavu dvou transcendentních rovnic:
⎛ ⎞ a ⎜ dp = + H ⎟ cos ϕ ⎜ 1 − e 2 sin 2 ϕ ⎟ ⎝ ⎠ ⎛ a 1 − e2 ⎞ . ⎜ pz = + H ⎟ sin ϕ ⎜ 1 − e 2 sin 2 ϕ ⎟ ⎝ ⎠
(
)
(4.9)
Jejich řešení je komplikované a nabízí se několik způsobů řešení. Jednou z možností je zavést substituci t = tg(ϕ), přičemž získáme rovnici:
pz ae 2
t= dp −
(
,
)
(4.10)
1+ 1− e2 ⋅ t 2
kterou řešíme analyticky nebo numericky. Analytické řešení není příliš vhodné pro počítačové systémy (např. mikroprocesor v GPS přijímači) a navíc vede na několik řešení, proto se zaměříme na numerické řešení s využitím metody prosté iterace, kde nové řešení parametru t získáme jako:
pz e2a
ti = dp −
(
.
)
1 + 1 − e 2 ⋅ t i2−1
(4.11)
x
Počáteční hodnotu lze definovat následovně:
t0 =
pz (1 − e 2 ) ⋅ d p
.
(4.12)
Geodetickou šířku pak určíme ze vztahu:
ϕ = arctg (t ) .
(4.13)
A elipsoidickou výšku určíme z rovnice:
⎛ a H = 1+ t ⎜d p − ⎜ 1 + 1 − e2 ⋅ t 2 ⎝ 2
(
)
⎞ ⎟ ⎟. ⎠
(4.14)
POSTUP ŘEŠENÍ Ad 4.1) V matlabovském skriptu funkce GeoidHeightWGS84(Lat, Lon) je definována tabulka výšky geoidu WGS84_GH_Table pro systém WGS84 pro geodetické ířky a délky po 10°. Jednotlivé řádky odpovídají geodetickým šířkám od +90° až po -90°. Sloupce pak geodetickým délkám v rozsahu od -180° po 180°. To je pochopitelně pro přesné určení nadmořské výšky příliš hrubé rozlišení. Doplňte proto skript této funkce o aproximační algoritmus, který odhadne výšku geoidu i mimo definovanou síť. Spolehlivé a relativně
jednoduché je použití postupného lineárního odhadu (viz obr. 2). Úkolem je určit výšku geoidu v bodě G. Nejprve provedeme odhad výšek v bodě E a F na hranách sítě, například pro geodetickou délku, ve druhém kroku pak již odhad výšky v cílové poloze G na základě znalosti výšek v polohách E a F, tedy aproximace v geodetické šířce.
Obr. 4.2 Grafické znázornění odhadu výšky ve 2D síti. Pro otestování použijte následující matlabovský skript: lon lat las los
= = = =
-180:1:+179; 90:-1:-89; length(lat) length(lon)
for lai=1:las, for loi=1:los, N(lai,loi)=GeoidHeightWGS84(lat(lai),lon(loi)); end end figure(1) surf(lon,lat,N) Skript je rovněž k dispozici jako GeoidWGS84plot.m. Výsledkem by měl být výstup v podobě aproximovaného geoidu WGS84 v rozlišení po 1°, který je na obrázku 4.3.
Obr. 4.3 Model aproximovaného geoidu pro WGS84. ad 4.2) Pro funkci přepočtu geodetických na kartézské souřadnice aplikujte rovnice (4.1) až (4.5), hlavičku definujte následovně: function [X Y Z]=GeoWGS842Cart(Lat, Lon, Alt) Výška nechť je udávána jako nadmořská, proto musí funkce volat sestavenou funkci z části 4.1 pro určení výšky geoidu. ad 4.3) Pro funkci přepočtu kartézských na geodetické souřadnice použijte postup prosté iterace podle rozboru. Jako kritérium pro ukončení iteračního procesu si zvolte podmínku, kdy rozdíl v hodnotách parametru t pro nové a předchozí řešení je menší než 1·10-6. Hlavičku definujte následovně: function [Lat Lon Alt]=Cart2GeoWGS84(X, Y, Z) Výška je opět udávána jako nadmořská. Na závěr proveďte několik kontrolních testů pro různé souřadnice s využitím zpětného výpočtu. Inspirujte se také přednáškou, kde jsou uvedeny příklady.
LITERATURA [4.1]
PANY, T. Navigation Signal Processing for GNSS Software Receivers. 1st ed. Boston: Artech House, 2010.
[4.2]
KAPLAN, E.D., HEGARTY, Ch.J. Understanding GPS. Principles and Applications. 2nd ed. Boston: Artech House, 2006.
[4.3]
ZIEDAN, N.I. GNSS Receivers for Weak Signals. 1st ed. Boston: Artech House, 2006
