ČESKÉ VYSOKÉ UČENÍ TECHNICKÉ V PRAZE FAKULTA STAVEBNÍ – OBOR GEODÉZIE A KARTOGRAFIE KATEDRA VYŠŠÍ GEODÉZIE název předmětu
Fyzikální geodézie úloha/zadání
název úlohy
Výpočet lokálního geoidu pro body na území ČR
3/7 školní rok
2010/11
semestr
1
skupina
NG1-88
zpracoval
datum
Jan Dolista
30. 11.
klasifikace
Výpočet lokálního geoidu pro body na území ČR Zadání: K dispozici máte síť bodů na území střední Evropy se známými zeměpisnými souřadnicemi 𝐵, 𝐿, nivelační výškou 𝐻 a měřeným tíhovým zrychlením 𝑔. Síť bodů je uložena v souboru TIHBODY.ASC. Dále vám poskytujeme parametry modelu geoidu EGM96, jmenovitě odlehlost od elipsoidu GRS80 a Fayovu anomálii v pravidelné mřížce bodů na území střední Evropy. Parametry modelu jsou uloženy v souboru EGM96. Vypočítejte odlehlosti zadaných bodů X1,...,X9 od geoidu EGM96 (tzv. lokální složku odlehlosti). Výpočet odlehlosti proveďte pomocí Stokesova integrálu. Jako vstupní hodnoty do integrálu volte rozdíly mezi Fayovou anomálií bodů sítě TIHBODY.ASC a Fayovou anomálií modelu EGM96 z oblasti obklopující zadané body. Tuto oblast (integrační síť) volte čtvercovou o rozměrech 2∘ × 2∘ . Vypočítejte rovněž celkovou odlehlost bodů X1,...,X9 od elipsoidu GRS80. Číselné zadání 7: bod X1 𝜙
bod X2 𝜆
𝜙
bod X3 𝜆
𝜙
bod X4 𝜆
𝜙
bod X5 𝜆
𝜙
𝜆
∘
′
∘
′
∘
′
∘
′
∘
′
∘
′
∘
′
∘
′
∘
′
∘
′
48
54
15
38
48
54
15
48
48
54
15
58
49
04
15
38
49
04
15
48
bod X6 𝜙
bod X7 𝜆
𝜙
bod X8 𝜆
𝜙
bod X9 𝜆
𝜙
𝜆
∘
′
∘
′
∘
′
∘
′
∘
′
∘
′
∘
′
∘
′
49
04
15
58
49
14
15
38
49
14
15
48
49
14
15
58
Vypracování: Veškeré výpočty byly provedeny v programu Octave.
1
Převod 𝜙, 𝜆 výpočetních bodů na 𝐵, 𝐿
Jako první byly převedeny zeměpisné souřadnice 𝜙, 𝜆 výpočetních bodů na geodetické 𝐵, 𝐿. K tomu bylo nutné určit hodnoty tížnicových odchylek ve výpočetních bodech. Ty byly určeny plošnou interpolací z globálního modelu EGM96. 𝐵 =𝜙−𝜉 𝐿=𝜆−
𝜂 , cos 𝜙
kde 𝜉, 𝜂 jsou tížnicové odchylky pro daný bod. Hodnoty tížnicových odchylek byly získány interpolací ze čtyř nejbližších bodů EGM96. Jelikož v modelu EGM96 jsou uvedeny geodetické souřadnice, byly pro interpolaci použity zeměpisné souřadnice výpočetního bodu jako přibližné hodnoty geodetických souřadnic. Vzhledem k malým hodnotám tížnicových odchylek nemá tato nepřesnost na výpočet vliv.
Plošná interpolace Jelikož je plošná interpolace používána ve výpočtu několikrát, byla pro tento účel v programu Octave vytvořena funkce. Vstupem této funkce jsou geodetické souřadnice bodu, v němž má být hodnota určena interpolací a čtveřice souřadnic bodů a funkčních hodnot, z nichž má být hodnota vyinterpolována. Výstupem funkce je pak funkční hodnota v zadaném bodě. Plošná interpolace je váženým průměřem funkčních hodnot ve známých bodech, kde jako váha je použita převrácená hodnota vzdálenosti známého a určovaného bodu. ∑︀
𝑋=
𝑋𝑖 𝑝𝑖 , 𝑝𝑖
kde 𝑋 je funkční hodnota v určovaném bodě, 𝑋𝑖 jsou funkční hodnoty ve známých bodech a 𝑝𝑖 jsou váhy. 1 𝑝𝑖 = , 𝑑𝑖 kde 𝑑𝑖 = 𝑅 · arccos (sin 𝐵 · sin 𝐵𝑖 + cos 𝐵 · cos 𝐵𝑖 · cos ∆𝐿), kde 𝑅 = 6, 371 · 106 𝑚 je poloměr Země, 𝐵 geodetická šířka určovaného bodu, 𝐵𝑖 geodetické šířka známého bodu a ∆𝐿 rozdíl délek známého a určovaného bodu. Kromě funkce pro plošnou interpolaci byla v programu Octave vytvořena funkce pro vyhledání čtyř nejbližších bodů v modelu EGM96 a funkčních hodnot v nich. Pro tu jsou vstupem souřadnice bodu k němuž mají být nejbližší body vyhledány a číslo sloupce ve kterém jsou funkční hodnoty pro hledanou veličinu.
2
Volba integrační oblasti
Integrační oblast byla zvolena 2∘ × 2∘ symetricky vzhledem ke středovému bodu, tedy bodu 𝑋5 . V rámci této oblasti byly zvoleny integrační plošky o velikosti ∆𝐵×∆𝐿, kde ∆𝐵 = 5′ a ∆𝐿 = 7.5′ . Integrační oblast byla zvolena pouze jedna pro všechny určované body, tedy kromě bodu 𝑋5 je výpočetní bod v rámci oblasti umístěn mírně excentricky, avšak výpočet je touto volbou značně urychlen.
3
Redukce měřeného tíhového zrychlení z povrchu na geoid
Soubor měřených tíhových bodů obsahuje měření na povrchu Země, tedy v určité výšce nad geoidem. Proto je nutné měřené tíhové zrychlení redukovat. K tomu je nutné kromě polohy bodu znát i nadmořskou výšku. Měřené tíhové zrychlení je v souboru tíhových bodů uvedeno v mGal/10 a je tedy nutné jej pro další výpočet převést na mGal. Dále je v souboru tíhových bodů použit jiný tíhový systém než u modelu EGM96 a proto je nutné tíhové zrychlení opravit o hodnotu -14mGal.
3.1
Fayova redukce
Jako první je provedena Fayova redukce na volném vzduchu, tedy výpočet tíhového zrychlení na geoidu. 𝑔(𝑃 1) = 𝑔(𝑃 ) + ∆𝐹 , kde ∆𝐹 = 0.3086 · 𝐻, kde 𝑔(𝑃 ) je tíhové zrychlení v bodě na povrchu Země, 𝑔(𝑃 1) je tíhové zrychlení v bodě na geoidu a 𝐻 je nivelovaná výška.
3.2
Fayova anomálie
V druhém kroku je vypočtena tíhová anomálie, tedy rozdíl mezi tíhovým zrychlením na geoidu a na elipsoidu. Jelikož pro redukce tíhového zrychlení na geoid byla použita Fayova redukce, jedná se o Fayovu anomálii. ∆𝑔 = 𝑔(𝑃 1) − 𝛾(𝑃 0), kde
(︁
)︁
𝛾(𝑃 0) = 9.780327 · 105 · 1 + 0.0053024 sin2 𝐵 − 0.0000058 sin2 (2𝐵) , kde 𝛾(𝑃 0) je normální tíhové zrychlení na elipsoidu a 𝐵 geodetická šířka bodu.
4
Tíhová anomálie pro každou z integračních plošek
Pro každou integrační plošku byla určena tíhová anomálie ve středovém bodě. Nejprve byly určeny souřadnice středového bodu a poté souřadnice rohových bodů integrační plošky. Následně byly pro každou plošku procházeny všechny body v souboru tíhových bodů a bylo zjišťováno, zda bod leží uvnitř plošky: ∙ Pokud v plošce neleží žádný bod, je hodnota tíhové anomálie středového bodu rovna nule. ∆𝑔 = 0 ∙ Pokud v plošce leží pouze jeden bod, je tíhová anomálie středového bodu rovna tíhové anomáli právě v tomto bodě. ∆𝑔 = ∆𝑔𝑏𝑜𝑑𝑢 ∙ Pokud v plošce leží více bodů, je hodnota tíhové anomálie středového bodu vypočtena jako vážený průměr ze všech bodů, které v plošce leží. Vahou je převrácená hodnota vzdálenosti od středu plošky. Platí tedy stejné vztahy jako pro plošnou interpolaci a liší se pouze počet známých bodů. ∆𝑔 = ∆𝑔𝑝𝑟𝑢𝑚𝑒𝑟
5
Redukce úplné tíhové anomálie na lokální složku tíhové anomáli
Pro středový bod každé integrační plošky byla z globálního modelu EGM96 určena hodnota globální složky tíhové anomálie ∆𝑔 𝐺 . Použita byla opět funkce pro vyhledání čtyř nejbližších bodů v modelu EGM96 a následně funkce pro plošnou interpolaci. Lokální složka tíhové anomálie pro středový bod integrační plošky je pak rozdílem celkové anomálie a globální anomálie. ∆𝑔 𝐿 = ∆𝑔 − ∆𝑔 𝐺 V případě, že v dané plošce neleží žádný tíhový bod, je i lokální složka tíhové anomálie rovna 0. Tento fakt byl při výpočtu zohledněn podmínkou, bez ní by v těchto ploškách byla hodnota lokální anomálie vypočtena jako −∆𝑔 𝐺 .
6
Lokální složka odlehlosti
Lokální složka odlehlosi se vypočte ze Stokesova integrálu, který je pro konečnou velikost integračních plošek nahrazen sumací přes všechny tyto plošky. 𝑁𝐿 =
𝑅 ∑︁ 𝑆(𝜓) · ∆𝑔 𝐿 · cos 𝐵 · 𝑑𝐵 · 𝑑𝐿, 4𝜋𝛾
kde 𝑆(𝜓) je Stokesova funkce 𝑆=
1 − 4 + 10𝑠2 − (3 − 6𝑠2 ) · ln(𝑠 + 𝑠2 ) − 6𝑠, 𝑠
kde
𝜓 , 2 kde 𝜓 je geocentrický úhel mezi výpočetním bodem 𝐵, 𝐿 a středem integrační plošky 𝐵 ′ , 𝐿′ 𝑠 = sin
𝜓 = arccos sin 𝐵 · sin 𝐵 ′ + cos 𝐵 · cos 𝐵 ′ · cos ∆𝐿 (︀
)︀
Pro dosazení do sumace je nutné velikost plošky 𝑑𝐵 a 𝑑𝐿 uvažovat v radiánech. Poloměr Země je volen 𝑅 = 6, 371 · 106 𝑚, normální tíhová anomáli ve výpočetním bodě je dána vztahem: (︁
)︁
𝛾 = 9.780327 · 105 · 1 + 0.0053024 sin2 𝐵 − 0.0000058 sin2 (2𝐵) , kde 𝐵 je geodetická šířka výpočetního bodu.
7
Úplná odlehlost
Úplná odlehlost geoidu a elipsoidu ve výpočetním bodě je součtem lokální a globální složky odlehlosti. Globální složka odlehlosti 𝑁 𝐺 je získána opět z modelu EGM96 a to plošnou interpolací. Použity jsou opět funkce pro vyhledání čtyř nejbližších bodů a funkce pro interpolaci. 𝑁 = 𝑁𝐿 + 𝑁𝐺
8
Číselné výsledky
bod 𝑁 𝐿 [m] 𝑁 [m]
X1 -0.31 46.84
X2 -0.36 46.58
X3 -0.38 46.28
X4 -0.30 46.96
X5 -0.31 46.83
X6 -0.32 46.60
X7 -0.30 46.90
X8 -0.33 46.82
X9 -0.34 46.63
Závěr: Na základě tíhových měření a globálního modelu geoidu EGM96 byly v zadaných bodech určeny hodnoty lokální odlehlosti geoidu od modelu EGM96 a následně i celková odlehlost od elipsoidu. Použita byla tíhová měření v bodech nerovnoměrně rozmístěných v okolí výpočetních bodů. Pro výpočet byla použita integrační oblast o velikosti 2∘ × 2∘ rozdělená na integrační plošky 5′ × 7.5′ . K výpočtu lokální odlehlosti byl použit Stokesův integrál nahrazený sumací. Pro výpočty vzdálenosti mezi body byla Země nahrazena koulí o poloměru 𝑅 = 6, 371 · 106 𝑚. Pro získání hodnot z globálního modelu pro určitý bod byla používána plošná interpolace ze čtyř nejbližších bodů modelu. Výpočty byly provedeny v programu Octave.
Přílohy: 1. Zdrojový kód a funkce pro výpočet v programu Octave
V Kralupech nad Vltavou 30.11.2010
Jan Dolista (
[email protected])
