Aktivita A07-03: Teoretické řešení problematiky transformace výšek a určení vybraných parametrů tíhového pole Země. Příloha 1 Popis řešení projektu za rok 2007 Všechny uvedené vzorce pro globální modely předpokládají vstup ve sférických souřadnicích (r, ϕ, λ). Uživatel bude moci zadávat své souřadnice buďto v systému WGS84/ETRS89 nebo v systému JTSK. Souřadnice v systému JTSK budou nejprve transformovány pomocí globálního transformačního klíče pro ČR do ETRS89/WGS84. Následovat bude transformace z geodetických souřadnic do sférických souřadnic (r, θ, λ). Tyto transformace jsou všeobecně známé ([1], [2]). Chyby transformace rovinných souřadnic ze systému JTSK do výše popsaných geocentrických systému nejsou z pohledu naší aplikace podstatné, protože jsou o několik řádů menší než rozlišení použitých globálních i lokálních modelů (5’x5’, resp. 30”x30”). Lokální modely jsou uloženy ve formě rastru v souřadnicovém systému WGS84.
Výpočty z globálních modelů Protože koeficienty globálních modelů potenciálu (GGM) mají obecně jiné parametry GM, a než jsou parametry elipsoidu WGS84, je nutná jejich transformace [1] (1)
∀n, m : C nm = C nm
GM EGM GM WGS 84
⎛ a EGM ⎜⎜ ⎝ aWGS 84
⎞ ⎟⎟ ⎠
n −1
, S nm = S nm
GM EGM GM WGS 84
⎛ a EGM ⎜⎜ ⎝ aWGS 84
⎞ ⎟⎟ ⎠
n −1
kde GM EGM , aEGM jsou parametry GGM a GM WGS 84 , aWGS 84 jsou parametry elipsoidu WGS84. Cnm a Snm jsou (plně normalizované) Stokesovy koeficienty GGM. V následujícím textu, není-li uvedeno jinak, používáme zkratky GM = GM WGS 84 , a = aWGS 84 . Výpočet gravitačního potenciálu V Vstup: souřadnice ϕ, λ, h v systému WGS84 Výstup: hodnota V v jednotkách m2.s-2
Nejprve transformujeme souřadnice do sférické soustavy souřadnic (r, θ, λ). Gravitační potenciál V počítáme z aktuálního globálního modelu geopotenciálu (GGM) pomocí vzorce (2)
n ⎤ GM ⎡ N max ⎛ a ⎞ V (r ,θ ,λ ) = ⎢1 + ∑ ⎜ ⎟ Vn (θ , λ )⎥ , r ⎣⎢ n=2 ⎝ r ⎠ ⎦⎥
kde Nmax je nejvyšší stupeň a řád daného GGM, θ = 90˚ - ϕ je pólová vzdálenost a funkce Tn je funkce definovaná předpisem: (3)
V n (θ , λ ) =
∑ (C n
m =0
nm
cos mλ + S nm sin mλ ) ⋅ Pnm (sin θ )
kde Cnm , S nm jsou (plně normalizované a transformované podle (1)) Stokesovy koeficienty GGM a Pnm jsou plně normalizované přidružené Legendrerovy funkce.
Výpočet potenciálu normálního tíhového pole U Vstup: souřadnice ϕ, h v systému WGS84 (na souřadnici λ výsledek nezávisí) Výstup: hodnota U v jednotkách m2.s-2
Za normální pole bereme pole elipsoidu WGS84. Po transformaci vstupních souřadnic do sférické soustavy souřadnic (r, θ, λ) určíme hodnotu normálního potenciálu jako (4) U = Vnorm + Φ, kde Vnorm je normální gravitační potenciál, který spočteme podle
(5)
2n k ⎤ GM ⎡ ⎛a⎞ Vnorm (r ,θ ) = ⎢1 − J 2 n ∑ ⎜ ⎟ P2 n (sin θ )⎥ , r ⎣⎢ n =1 ⎝ r ⎠ ⎦⎥
a Φ je odstředivý potenciál, který počítáme pomocí vztahu 1 (6) Φ (r ,θ ) = ω 2 r 2 sin 2 θ , RRR5 2 kde k určuje stupeň rozvoje normálního potenciálu (stačí zvolit k=10) a P2n(cos θ) jsou Legendrerovy polynomy. Konstanta J2 je pro elipsoid WGS84 pevně určena a koeficienty J2n se z ní určují pomocí vzorce
(7)
J 2 n = (−1) n +1
(
3e 2 n 1 − n + 5nJ 2 e − 2 (2n + 1)(2n + 3)
)
Plně normalizované koeficienty z nich spočítáme pomocí vztahu J 2n (8) J 2n = 2n + 1 Výpočet tíhového potenciálu W Vstup: souřadnice ϕ, λ, h v systému WGS84 Výstup: hodnota W v jednotkách m2.s-2
Hodnotu W spočítáme snadno pomocí vzorce (9) W(r, θ, λ) = V(r, θ, λ) + Φ(r, θ), kde V spočítáme podle (2) a Φ podle (6). Výpočet poruchového potenciálu T Vstup: souřadnice ϕ, λ, h v systému WGS84 Výstup: hodnota T v jednotkách m2.s-2
Poruchový potenciál T se spočítáme podle vztahu (10) T(ϕ, λ, h) = W(ϕ, λ, h) - U(ϕ, λ, h), kde hodnotu W(ϕ, λ, h) určíme podle (2) a U(ϕ, λ, h) podle (5). Poruchový potenciál můžeme rovněž vyjádřit ve formě řady (11) kde
GM EGM − GM WGS 84 GM + T (r ,θ ,λ ) = aWGS 84 r
N max
∑ n=2
n
⎛a⎞ ⎜ ⎟ Tn (θ , λ ) , ⎝r⎠
(12)
Tn (θ , λ ) =
∑ (C n
m =0
nm
cos mλ + S nm sin mλ ) ⋅ Pnm (sin θ )
Tn (θ , λ ) = (C n 0 + J n ) Pn (sin θ ) + ∑ (C nm cos mλ + S nm sin mλ ) ⋅ Pnm (sin θ )
pro liché n
n
pro sudé n
m =1
Výpočet tíhové anomálie Δg Vstup: souřadnice ϕ, λ, h v systému WGS84 Výstup: hodnota Δg v jednotkách mGal (1 mGal = 10-5 ms-2)
Opět nejprve transformujeme geodetické souřadnice do sférických a hodnotu Δg vyjadřuje vzorec n
⎛a⎞ (13) (n − 1)⎜ ⎟ Tn (θ , λ ) ∑ ⎝r⎠ n=2 Výsledek je nutné vynásobit konstantou k=105, abychom dostali výsledek v mGal. GM Δg (r ,θ , λ ) = 2 r
Nmax
Výpočet radiální derivace tíhového zrychlení gr Vstup: souřadnice ϕ, λ, h v systému WGS84 Výstup: hodnota gr v jednotkách mGal na km (fyzikální rozměr s-2)
Hodnotu gr rozložíme na dvě složky ∂g ∂γ ∂Δ g (14) gr = ≅ + ∂r ∂h ∂r Dílčí složky radiální derivace spočítáme pomocí vztahů (15)
2γ aϖ 2 ∂γ = − e (1 + f + m − 2 f sin2 θ ), m = , γe a ∂h
kde γ e je velikost normálního zrychlení na elipsoidu, kterou spočítáme podle (19) a f jeho zploštění. Druhou složku spočítáme podle vzorce ∂Δg GM = 3 r ∂r
n
⎛a⎞ (n + 2)(n − 1)⎜ ⎟ Tn (θ , λ ) ∑ ⎝r⎠ n=2 Následně vynásobíme řešení konstantou k=108, abychom dostali výsledek v požadovaných jednotkách. (16)
Nmax
Výpočet výškové anomálie ζ Vstup: souřadnice ϕ, λ, h v systému WGS84 Výstup: hodnota ζ v jednotkách m
Výškovou anomálii určíme ze vztahu T (ϕ , λ , h) (17) ζ (ϕ , λ , h) = γ (ϕ , λ , h) Výpočet velikosti vektoru normálního zrychlení γ Vstup: souřadnice ϕ, λ, h v systému WGS84 Výstup: hodnota γ v jednotkách mGal (1 mGal = 10-5 ms-2)
(18)
⎡ ⎣
1+
(19)
2 a
γ (ϕ , h) = γ (ϕ ,0) ⎢1 − (1 + f + m − 2 f sin 2 ϕ )h +
γ (ϕ ,0) = γ a
3 2 ⎤ aϖ 2 ) , = h m ⎥⎦ γe a2
bγ b − aγ a sin 2 ϕ aγ b 1 − f sin 2 ϕ
Konstanty a, b, f, γa, γb, ω jsou parametry zvoleného elipsoidu. Parametry elipsoidu WGS84 jsou (podle [3]) a = 6 378 137 m b = 6 356 752, 3142 m e2 = 6,694 379 990 14 . 10-3 f = 1/298,257 223 563 ω = 7 292 115 . 10-11 rad s-1 GM = 3 986 004, 418 . 108 m3 s-2 γa = 9,780 325 3359 m s-2 = 9,832 184 9378 m s-2 γb Parametry se pro různé modely GGM, ale vždy jsou uvedeny v jeho dokumentaci. Generování přidružených Legendrerových funkcí Pnm Klasický způsob generování přidružených Legendrerových funkcí (viz např. [1]) je založen na rekurentních formulích. S těmi je možné počítat modely do stupně a řádu 360, v nichž máme současný GGM. Nový EGM08 bude ovšem obsahovat koeficienty do stupně a řadu 2160. Pro takto vysoký stupeň a řád již nelze pomocí těchto jednoduchých rekurentních formulí přidružené Legendrerovy funkce počítat. Vzhledem k důležitosti generování těchto funkcí je toto téma v současnosti intenzivně studováno. Jedním z úkolů na příští rok je volba optimálních algoritmů pro generování přidružených Legendrerových funkcí.
Výpočty z lokálních modelů Všechny výpočty z lokálních modelů se dělají pomocí interpolace z předem připraveného rastru. Výhodou proti použití globálních modelů je vyšší rychlost výpočtu, vyšší prostorové rozlišení lokálního modelu i jeho vyšší přesnost. Nevýhodou ovšem je, že veličiny jsou vždy vztaženy pouze k jedné ploše, na níž jsou spočítány (topografie, geoid apod.). Z lokálního modelu lze pak interpolovat hodnoty opět pouze na této ploše. Většinu hodnot počítaných z lokálních modelů tedy budeme určovat prostou interpolací z předpočítaných rastrů. Některé hodnoty však velmi silně závisí na výšce a uživatel musí svou přesnou nadmořskou výšku pro výpočet zadat. Jedná se o velikosti vektoru tíhového zrychlení g a tíhovou anomálii Δg. Výpočet velikosti vektoru tíhového zrychlení g Vstup: souřadnice ϕ, λ v systému WGS84, H nadmořská výška* (normální Moloděnského) Výstup: hodnota g v jednotkách m.s-2
Výpočet velikosti vektoru tíhového zrychlení na povrchu Země se bude počítat podle vzorce g (ϕ , λ , H ) = Δg B (ϕ , λ ) + AB (ϕ , λ ) + AT (ϕ , λ ) − F ( H ) + γ (ϕ , H ) (20) *
Výška H musí být výška na povrchu Země, tímto postupem nelze interpolovat hodnoty nad zemským povrchem. Nelze ji ovšem brát z DMT, protože velikost tíže velmi silně závisí na výšce a proto musí být spočtena přesně pro uživatelem zadanou nadmořskou výšku.
Hodnota Δg B (ϕ , λ ) je Bouguerova tíhová anomálie, která bude vyinterpolována z předpočítaného rastru (výpočet je velmi náročný a nelze jej provádět online). AB (ϕ , λ ) je gravitační zrychlení generované Bouguerovou sférickou slupkou, která bude interpolovaná z předpočítaného rastru. AT (ϕ , λ ) je terénní korekce ve sférické aproximaci, která bude rovněž interpolovaná z předpočítaného rastru. Hodnota F(H) je redukce ve volném vzduchu, která se spočítá podle vztahu (21) F(H) = 0,3086h a γ (ϕ , H ) je velikost normálního zrychlení, které spočteme podle vzorce (18). Výpočet tíhové anomálie Δg Vstup: souřadnice ϕ, λ v systému WGS84, H nadmořská výška* (normální Moloděnského) Výstup: hodnota Δg v jednotkách m.s-2
Hodnotu spočítáme podle vztahu Δg (ϕ , λ , H ) = g (ϕ , λ , H ) − γ (ϕ , H ) (22) Hodnotu g (ϕ , λ , H ) spočítáme pomocí (20) a γ (ϕ , H ) podle (18). Výpočet ostatních parametrů počítaných z lokálních modelů Vstup: souřadnice ϕ, λ v systému WGS84. Výstup: • příčná složka tížnicové odchylky η ve stupních • meridiánová složka tížnicové odchylky ξ ve stupních • Bouguerova anomálie Δg B (ϕ , λ ) v mGal • terénní korekce AT (ϕ , λ ) v mGal • gravitační zrychlení generované Bouguerovou sférickou slupkou AB (ϕ , λ ) • převýšení kvazigeoidu ζ v m • převýšení geoidu N v m
Tyto hodnoty máme předpočítané a uložené ve formě rastru. Jednotlivé rastry jsou předpočítané vždy vztažené k nějaké výšce (ke geoidu, k topografii, ...) a výsledek interpolace je vždy vztažen k odpovídající výšce. Zadávání výšky od uživatele zde tedy nemá význam. Interpolace
Hodnoty z jednotlivých rastrů budeme interpolovat pomocí bilineární interpolace z okolních buněk rastru (23)
⎛ ⎜ f (ϕ , λ ) = (ϕ 2 − ϕ ϕ1 − ϕ )⎜ ⎜ ⎜ ⎝
f (ϕ1 , λ1 ) ΔϕΔλ f (ϕ 2 , λ1 ) ΔϕΔλ
f (ϕ1 , λ 2 ) ⎞ ⎟ ΔϕΔλ ⎟⎛⎜ λ 2 − λ ⎞⎟ , f (ϕ 2 , λ 2 ) ⎟⎜⎝ λ1 − λ ⎟⎠ ΔϕΔλ ⎟⎠
kde ϕi, λi jsou souřadnice okolních 4 buněk rastru, f(ϕi, λi) hodnoty těchto buněk rastru a Δϕ , Δλ prostorové rozlišení rastru. Popis jednotlivých předpočítaných rastrů
Rastry máme uložené v textovém formátu a v databázi GIS GRASS, která nám umožňuje export do celé řady dalších formátu. Konkrétní formát rastrových dat pro aplikaci bude zvolen až ve fázi programování znalostního systému. Podrobný popis těchto rastrů a jejich vzniku
zde nebude pro svůj rozsah uváděn, většinou jsou uvedeny odkazy na články, ve kterých jsou data popsána. Digitální model terénu Digitální model hustoty hornin Tíhové zrychlení na povrchu Země Popis těchto rastrů je uveden v [6]. Příčná složka tížnicové odchylky η Meridiánová složka tížnicové odchylky ξ Převýšení geoidu N Tyto rastry byly spočítány v rámci řešení grantu MŠMT „Řešení přesných modelů geoidu a kvazigeoidu pro oblast střední Evropy“ P. Novákem. Jejich podrobný popis je uveden v [4]. Terénní korekce AT (ϕ , λ ) Výpočet terénní korekce pro území ČR je popsán v [5]. Gravitační zrychlení generované Bouguerovou sférickou slupkou AB (ϕ , λ ) Tento jsme získali z rastru výšek pomocí jednoduchého vztahu (24) AB(ϕ, λ) = 0,1119 H(ϕ, λ) Bouguerova anomálie Δg B (ϕ , λ ) Hodnoty vygenerované na povrchu geoidu. Rastr Bouguerových anomálii jsme spočítali z rastru pozemních tíhových dat ( g (ϕ , λ ) ) a z digitálního modelu terénu (H) podle známého vztahu [3] Δg B (ϕ , λ ) = g (ϕ , λ ) − AB (ϕ , λ ) − AT (ϕ , λ ) + F ( H ) − γ (ϕ , H ) (25) Reference [1] Heiskanen, W. A., Moritz, H.: Physical geodesy. Freeman and Co., San Francisco 1967. [2] Cimbálník, M., Kostelecký, J.: Direct transformation between ETRS-89 and the Czech Cadastral System S-JTSK. Report on the Symposium of the IAG Subcomm. for the EUREF held in Ankara, 1996. Veroeff. der Bayer. Kommission fuer Int. Erdmessung der Bayer. Akad. der Wissenschaften, Heft No. 57, Muenchen 1996, printed 1997, p. 325 - 330. ISSN 03407691, ISBN 3 7696 9620 4 [3] Hoffmann-Wellenhof, B., Moritz, H.: Physical Geodesy. SpringerWienNewYork, Wien, 2005. ISBN 3-211-23584-1 [4] Novák, P.: Evaluation of local gravity field parameters from high resolution gravity and elevation data. Contributions to Geophysics and Geodesy 36: 1-33. 2006. [5] Kadlec, M.: Výpočet topografických oprav tíhových dat pro určení přesného regionálního modelu geoidu. (2007) JUNIORSTAV 2007 - Sborník anotací, plné texty na CD, ISBN 97880-214-3337-3, Vysoké učení technické v Brně, Fakulta stavební, Brno [6] Kadlec, M., Kostelecký J. ml., Novák P.: Databáze pro výpočty parametrů tíhového pole Země pro střední Evropu. (2007) Geodetický a kartografický obzor, č. 12/2007.
