MSc
BMEEOAFML01
Fizikai geodézia és gravimetria / 18. A FIZIKAI GEODÉZIÁBAN ALKALMAZOTT SZOFTVEREK ÁTTEKINTÉSE. A fizikai geodézia szakterületén sok olyan eljárás létezik, ami igen számításigényes. Ilyenek a különböző nagy felbontású geopotenciális modellekkel végzett számítások, amelyekkel a nehézségi erőtér különböző összetevői (potenciálzavar, geoidunduláció, nehézségi rendellenességek, függővonal elhajlások, Eötvös-tenzor különböző elemei) határozhatók meg. A spektrális eljárásokon alapuló szoftverek széles körben alkalmazást nyertek a geoidszámításban és a különböző felszínmodellekkel végzett hatékony terephatás számítás tekintetében. A legkisebb négyzetek módszerén alapuló kollokáció is igen számításigényes lehet a rendelkezésre álló adatok mennyisége függvényében. A terephatás számításához sokszor közvetlen numerikus integrálási eljárásokat alkalmaznak. Az árapály korrekciók numerikus előállítása is elengedhetetlen a gravimetriai adatok redukciójához, vagy akár a GNSS mérések korrekciójához – különösen a precíz pontmeghatározási, PPP feldolgozási technológia – esetében. Az űrgravimetriai és űrgradiometriai mérések feldolgozása pedig legújabban került előtérbe, és a felhasználók munkáját olyan programok segítik, amelyek nagyban megkönnyítik az ilyen jellegű mérési, számítási eredmények felhasználását és feldolgozását. A továbbiakban megismerünk néhány olyan, a fizikai geodézia területén sokak által használt szoftvert és feldolgozási eljárást, amelyekkel akár a saját adatainkkal mi magunk is különböző számításokat végezhetünk. Néhány esetben az ilyen számításokat konkrét példákkal is illusztráljuk.
Számítások geopotenciális modellekkel A geopotenciális modellben használható gömbfüggvény-sorfejtés maximális lehetséges fokszáma szerint megkülönböztethetünk alacsony fokszámú (nmax = 720) modellel végzett számításra alkalmas programokat, illetve kifejezetten magas fokszámú (nmax = 2190) EGM2008 modellel végzett számításra alkalmas programokat. A különbségtétel azért is indokolt, mert a magas fokszámú Legendre függvények számításához speciális eljárás szükséges, és erre az egyszerűbb programok nem képesek. A számítható mennyiségek szerint a legtöbb szoftverben a geoidundulációk illetve a nehézségi rendellenességek mind pontbeli mind rács értékekre előállíthatók. A függővonal elhajlások és gravitációs gradiensek számítására viszont már csak viszonylag kevés program alkalmas. Az ilyen programokra példaként megemlítjük a http://www.geod.bme.hu/tantargyak/korszmat/Euler-egyeb.zip állományban található geopot.e eljáráskönyvtárban levő szubrutinokat, melyekkel koordinátákkal megadott pontokban, illetve szabályos rácspontokban történő számítás is lehetséges. A spektrális módszereknél már megismert FFT-vel való gyorsítás is lehetséges az említett eljárások esetében, és a számított mennyiségek lehetnek geoidmagasságok, nehéz-
1
ségi rendellenességek, az Eötvös-tenzor összes eleme, az Eötvös-tenzor sugárirányú 1. és 2. deriváltjai, valamint a fentiek hibái a gömbfüggvény-együtthatók hibái alapján.
Az EGM2008 geopotenciális modell Az EGM2008 geopotenciális modell Az NGA (National GeospatialIntalligence Agency) legújabb fejlesztésű modellje, amely 2160 fokig és rendig tartalmazza (további 2190 fokig a 2160-nál nem nagyobb rendű tagok tekintetében) a gömbfüggvény együtthatókat, összesen 4 802 666 Cnm és Snm együtthatót. A modell a nehézségi erőtér hosszú hullámú összetevőit a GRACE műholdak mérései alapján tartalmazza. A modell megalkotásához felhasználtak egy 5’x5’ felbontású nehézségi adatokat tartalmazó adatbázist (lásd az alábbi ábrát), az SRTM (Shuttle Radar Topographic Mission) 30’x30’ felbontású topográfiai adatbázisát az északi szélesség 60°-tól a déli szélesség 58°-ig, az északi és a déli sarkok környezetében az ICESat méréseit, valamint a közepes tengerszintre (MSS és DOT) az altiméteres méréseket a nehézségi rendellenességek számításához.
1. ábra. Az 5’x5’ nehézségi rendellenességek adatnyerési forrásai az EGM2008-as modell számára
2. ábra. Az EGM2008 modellből számítható nehézségi rendellenességek 2
Az EGM2008 modell esetében a geoidmagasságok gömbfüggvényszoból történő számítása igen időigényes. Ezért a modell megalkotói a felhasználók számára többek között elkészítették a geoidmagasságok különböző felbontású (10’, 2.5’, 1’) rács adatbázisait. Ez viszonylag gyors és egyszerű geoidszámítást tesz lehetővé. Megemlítjük ezzel kapcsolatban az AllTrans EGM2008 Calculator nevű programot (Hans-Gerd Duenck-Kerst munkája), ami pontonkénti és rácsra történő számítást tesz lehetővé az említett különböző felbontású rács adatbázisok és különböző interpolációs eljárások segítségével.
3. ábra. Az AllTrans EGM2008 Calculator (Hans-Gerd Duenck-Kerst)
4. ábra. A programmal számított geoidkép Magyarországra (1’ x 1’-es rács adatbázis alapján):
3
0
5
0 5
00
5 7.5
7. 5
10
.5 -2
-7.5 -5
2.5
0
-2.5
10
2.5
5
5
0
10
-5
0
2.5
2.5
10
0
-5 -2.5 0
0
7.5
2.5
7.5
-5
5
7.5 5
2.5
0
-2.5
0 .5 2
7.5
5
5
0
5
0
7.5
2.5
10 2.5
5 7.5
7.5
10
7.5
5
7.5 5
0
2.5
5 2.
7.5 5
5 7.
7.5
7.5
0
-5
5
5
2.5
5 2.
5
5
0
10
5 7.
2.5
5
5 2.
10 5
5
5
2.5
7.5
2.5
0
7.5
7.5
5 7.5
10
7.5
2.5
7. 5
10
10
0 7.55
5
Az EGM2008-as modellből számított geoid illeszkedését megvizsgáltuk az Országos GPS hálózat (OGPSH) szintezett pontjain. Az eltéréseket a pontszám függvényében az alábbi ábra bal oldalán láthatjuk (Szűcs E. munkája). A jobb oldali ábrán az eltérések területi eloszlását figyelhetjük meg. Szembetűnő, hogy az ország délkeleti részén negatív értékeket kapunk.
-5
5. ábra. Az EGM2008 modell illeszkedése 340 OGPSH pontban 1. táblázat. Az EGM2008 modellből számított geoidundulációk eltéréseinek statisztikai jellemzői az OGPSH 340 szintezett pontjában, illetve kihagyva a legnagyobb eltérést mutató 34 pontot átlagérték min. érték max. érték szórás (m) (m) (m) (cm) 340 OGPSH pont 306 OGPSH pont
0.051
-0,235
0,386
7,81
0.042
-0,115
0,121
6,32
10 17 .5
7.5 10 12.5 15
5
20
7.5
10
7.5
2.5 2.5
5
10
0
7. 5
-2.5 5
6. ábra. Az EGM2008 modell illeszkedése OGPSH pontban
4
0
10
2.5
7.5
7.5
5
5
0
10
10
-2.5
7.5
10
-5
2.5
10
7.5
12.5
10
7.5
5 7.
7.5 10
2.5
10
5 7.
5 7.
10
5
12.5
10
5
12.5
5
5
12 .5
Az EGM2008-as modellből számított geoid illeszkedését megvizsgáltuk OGPSH 95 újonnan szintezett pontjain is. Az eltéréseket az alábbi ábrákon láthatjuk. Ismét szembetűnő, hogy az ország délkeleti részén negatív értékeket kapunk. Ezek az eltérések talán kapcsolatban vannak az EOMA újramérése kapcsán tapasztalt jelentős magasságváltozásokkal.
2. táblázat. Az EGM2008 modellből számított geoidundulációk eltéréseinek statisztikai jellemzői az OGPSH 95 újonnan szintezett pontjában átlagérték min. érték max. érték szórás (m) (m) (m) (cm) 95 OGPSH pont
0.061
-0,055
0,221
5,3
A speciális alapozású szintezési ún. K-pontok magasságváltozásait a KeletMagyarország északi részét lefedő 8, 9, 10-es poligon 1. és 2. epochában végzett mérései (1. epocha: 1975-1978, 2. epocha: 2007-2009) alapján, Busics Gy. munkája nyomán az alábbi ábrán láthatjuk.
7. ábra. A K-pontok magasságváltozása a KMO 1. és 2. mérése között. Forrás: Busics Gy (2010): Az EOMA újramérésének előzetes eredményei az első három poligonban. Geomatikai Közlemények XIII/2, 141-148. o.
GMT (Generic Mapping Tools) és Mirone A rendkívül rugalmasan használható, Wessel és Smith által fejlesztett és a földtudományokban sokszor használt Generic Mapping Tools (GMT) program nem csak az adatok ábrázolására alkalmas, hanem a programrendszerbe sokféle feldolgozási eljárást is beépítettek. A GMT-hez Matlab-ban fejlesztett egyik felhasználóbarát grafikus felület a GMT-hez a Mirone (http://w3.ualg.pt/~jluis/mirone/ ), amelynek futtatásához nem szükséges a Matlab-ot telepíteni. Ez a program további képességekkel is rendelkezik azokon túl, amelyek a GMT-be is be lettek építve. Többek között a Mirone többféle rács adatformátumot képes kezelni. Ilyenek a GMT/Netcdf, SURFER 6/7, Encom, Arc/Info, ENVI raster, Erdas, ESRI, Geosoft, GeoTIFF, JPEG2000, ENVISAT, DTED, SRTM, USGS DEM, stb. formátumok. A Mirone tud a Google Earth által használt .kmz fájlba is menteni, így az adataink a Google Earth-ben is megjeleníthetők. A Mirone adatfeldolgozási képességei között szerepel például az FFT spektrum számítás, digitális szűrés, különböző képfeldolgozási eljárások alkal-
5
mazása, georeferálás, rajzeszközök, lemeztektonikai kalkulátor, szeizmológia (fészekmechanizmusok), cunami terjedés, rugalmas deformáció számítás, domborzatelemzés, vetületek, stb. A Mirone munkafelületét az alábbi ábrán láthatjuk.
8. ábra. A Mirone program kezelőfelülete
A GRAVSOFT programrendszer A GRAVSOFT 49 önálló FORTRAN nyelven írt program együttese, amelyeket 1973-tól kezdve C.C. Tscherning, R. Forsberg, P. Knudsen (Dánia), D. Arabelos (Görögország) fejlesztettek és még jelenleg is folyamatosan fejlesztenek. Az adatok kezelésére különböző interpolációs (pontbeli és rács) eljárásokat tartalmaz. A GRAVSOFT-tal lehetséges nagy fokszámú gömbfüggvény sorfejtés segítségével különböző nehézségi erőtér paraméterek számítását elvégezni, így képes az EGM2008-as modellel is dolgozni. A programrendszerben szerepel még a terepi korrekció számítása, a Stokes integrál numerikus megoldása, különböző spektrális eljárások alkalmazása. Külön programok szolgálnak a legkisebb négyzetes kollokáció gömbi és sík közelítésben történő alkalmazására, valamint altiméteres mérések feldolgozására. Megemlítjük a GRAVSOFT programhoz Python programozási nyelven készített grafikus felületet, a PyGravsoft-ot, amely a GRAVSOFT 25 önálló programjának közös grafikus felülete.
6
9. ábra. A PyGravsoft grafikus kezelőfelülete Ez a grafikus felület kényelmesebbé teszi a programok futtatását, a számítási paraméterek megadását a gyakrabban előforduló számítási feladatok esetében.
A HGTUB2000B geoidmeghatározás lépései Az alábbi ábrák (amelyek a Mirone-nal készültek) illusztrálják azokat a konkrét lépéseket, amelyeket elvégeztünk a HGTUB2000B magyarországi geoidmegoldás számításához. Zárójelben megadjuk az egyes lépésekben alkalmazott programokat, illetve az előállított adatállományok nevét. 1. EGM96 geopotenciális modell (GRAVSOFT/harmexp): • geoidundulációk meghatározása (egm96.n) • nehézségi rendellenességek számítása (egm96.g)
10a. ábra. EGM96 geoidundulációk (m) (egm96.n)
10b. ábra. EGM96 nehézségi rendellenességek (mGal) (egm96.g)
7
2. Faye anomáliák redukciója • geopotenciális modellel (dgfreet.grd) • terepi korrekcióval
11a. ábra. Faye nehézségi rendellenességek (mGal)
11b. ábra. EGM96 redukált Faye nehézségi rendellenességek (mGal) (dgfreet.grd)
3. Maradék geoidunduláció FFT-vel (FFTGEOID) • Indirekt hatás hozzáadása • EGM96 geoidundulációk hozzáadása
12a. ábra. Terepi korrekciók (mGal)
12b. ábra. Terepi korrekciókkal redukált Δg (mGal)
Összehasonlítás céljából a maradék geoidundulációkat a kanadai Calgary egyetem PhD hallgatója, YeCai Li (1994) által fejlesztett FFTGEOID programmal is kiszámítottuk. Ez a program 1D FFT-vel számol, hasonlóan a GRAVSOFT programrendszer GEOFOUR programjához. A két program közötti alapvető különbség az, hogy az FFTGEOID program a transzformáció előtt 100% zérus feltöltést alkalmaz a nehézségi rendellenességek tömbjére, és ezáltal elkerüli a körkonvolúciónak nevezett jelenséget, amely a számítási terület szélén meghamisítja az adatokat. Ez jól látható az alábbi ábrákon, ahol a két programmal számított maradék geoidundulációkat, illetve azok eltéréseit mutatjuk be. Ha kiküszöböljük a körkonvolúciót, akkor (ahogy azt M. Sideris is kimutatta), az FFT-vel történő számítás pontosan ugyanazt az eredményt fogja szolgáltatni mint az adatok numerikus integrálása, csak sokkal gyorsabban. 8
13a. ábra. Maradék geoidunduláció FFTvel (m) – FFTGEOID
13b. ábra. Maradék geoidunduláció FFT-vel (m) – GEOFOUR
14a. ábra. Maradék geoidundulációk különbsége (m) az FFTGEOID és GEOFOUR programmal
14b. ábra. Indirekt hatás (m)
15. ábra. A kész geoidmegoldás (HGTUB200B) A szintvonalak értékköze 0.5 m
9
Terephatás számítása A GRAVSOFT programrendszer részei a TC, illetve TCFOUR programok, amelyekkel a terephatás számítása végezhető el mindenféle gravimetriai mennyiségre (nehézségi rendellenességek, nehézségi zavar, függővonal-elhajlás, magassági rendellenesség, vertikális gravitációs gradiens, Eötvös-tenzor összes eleme). A TC program prizma integrálást hajt végre a számítási pont környezetében az alábbi ábra szerint besűrítve a magassági adatok rácsát. Képes arra is, hogy egy belső és egy külső (kisebb felbontású) rács kombinációját alkalmazza a számítás felgyorsítása érdekében.
16. ábra. Prizma integrálás a TC programmal A TCFOUR program a terephatás számítását gyors Fourier transzformációval végzi.
Terephatás számítása poliéder modellből A Fortran 90 programnyelven írt PolyGrav program (http://www.geod.bme.hu/gtoth/PolyGrav.html) képes egy tetszőleges homogén sűrűségeloszlású poliéder (síklapokkal határolt test) gravitációs hatásának számítására a Holstein (2003) által kidolgozott összefüggések alapján.
P
bij r2ij
ni · ri
r1ij
ni
j él hij · rij
tij
i lap
poliéder
17. ábra. 10
hij = tij × ni
Az összefüggések egy (i-edik lapon a j-edik) élhez kapcsolt ortonormális ( hij , tij , ni ) vektorhármasra épülnek (lásd a 17. ábrát). Az ni vektor az i-edik lap kifelé mutató normális egységvektora. A tij vektor az i-edik lapon a j-edik él ni normálvektor irányából nézve óramutató járásával ellenkező irányba mutató értelmű egységvektora. A hij vektor pedig a hij = tij × ni vektorszorzatból adódik, és jobbsodrású rendszerré egészíti ki a vektorhármast. A számítási összefüggésekben alapvető szerepe van még a bij vektornak, amely a hij és ni vektorok síkjába esik, tehát azok lineáris kombinációja az alábbiak szerint: b ij = h ij ln
1 + Λ ij 1 − Λ ij
− n i 2sign(n i ⋅ ri ) arctan(λ ij )
Az összefüggésben a sign( ) az előjel függvény, ri a P számítási pontból a poliéder iedik lapjának tetszőleges pontjába mutató vektor, továbbá ni · ri a két vektor skaláris szorzatát jelöli. A dimenziótlan Λij és λij számokat pedig az alábbi képletek definiálják: Λ ij =
sij r1ij + r2ij
λ ij =
,
1 2
(r
1ij
(h
ij
⋅ rij )Λ ij
+ r2ij − sij Λ ij ) + ni ⋅ ri
.
A fenti összefüggésekben r1ij és r2ij a P pontnak a j-edik él két végpontjától mért távolsága illetve sij a j-edik él hossza. A P pontban test V gravitációs potenciálja, a gravitációs erő g vektora illetve az E Eötvös tenzor (gravitációs gradiens tenzor) a test összes lapjára és az adott lap öszszes élére vett alábbi összegekkel állítható elő (a ⊗ jel a vektorok ún. külső, diadikus szorzatát jelenti és ez egy tenzort eredményez; G valamint ρ a Newton-féle tömegvonzási állandót illetve a sűrűséget jelöli): V (P ) = 12 Gρ ri ⋅ n i b ij ⋅ rij i
j
g (P ) = Gρ n i b ij ⋅ rij i
.
j
E (P ) = Gρ n i ⊗ b ij i
j
A fenti képleteket az alábbi Fortran 90 nyelven programozott számítási eljárással értékeljük ki: subroutine GrEl(Corner, Face, rho, Points, dV, dG, dE) A szubrutin (eljárás) paraméterei a következők: Corner – a test csúcspontjainak koordinátái Face – az egyes lapok sarokpontjainak indexei a Corner tömbben rho – a test sűrűsége (kg/m3) Points – a számítási pont(ok) koordinátái dV, dG, dE – a számított eredmények az adott pontokban (1, 3 ill. 6 elem) Megjegyzés: Mivel az Eötvös-tenzor szimmetrikus, ezért van csak 6 független eleme: Exx, Eyy, Ezz, Exy, Exz, Eyz. A szubrutin az eredményeket J/kg, μGal illetve E (Eötvös, 10-9 s-2) mértékegységekben adja meg.
11
A számításhoz szükség van még a poliéder adatainak megadására, illetve – az eredmények ábrázolása céljából – azok megfelelő formátumú kiíratására. Erre a célra a számítást vezérlő főprogram (PolyGrav) a következő szubrutinokat hívja meg:
PolyGrav
Read_OFF
VTK_Write_Grav
18. ábra. A számítást vezérlő PolyGrav főprogram Az input adatokat a program (a Read_OFF eljáráson keresztül) egy ún. ’object file format’ (.off) formátumban képes fogadni.
Példa: terephatás számítása a Makádi teszt területen Az alábbi ábrán láthatjuk az 5 m × 5 m-es felbontású domborzatmodellel poliéder tömegmodell alapján a felszín felett 1 m-es magasságban számított vizszintes gravitációs gradienseket a Csepel-sziget déli részén, Makád közelében található teszt területen. Megfigyelhetjük, hogy még közel sík terepen is jelentős, 20 E-t meghaladó nagyságú gradienseket kaptunk. 194800
9
8.6
10
23.2
11
4.6
12
0.9
13
4.5
14
3.7
15
3.2
19
0.4
20
16.5
21
0.1
22
13.6
23
9
24
10.1
25
11.1
30
5.5
31
8.9
32
10
33
18.4
34
18.2
35
9.9
40
5.7
41
2.1
42
3.8
43
4
44
12.8
45
6
50
8.6
51
7.9
52
4.9
53
6.1
54
3.3
55
15
60
3.9
61
7
62
16.6
63
21
64
26.7
65
22.5
194600
194400
EOV-X
194200
194000
193800
193600 Horizontális gradiens [Eötvös] 193400
10.0 20.0
193200 641200 641400 641600 641800 642000 642200 642400 642600 642800 643000 643200 EOV-Y
19. ábra. A terephatás a vízszintes gradiensekre (E)
12
Árapály hatás számítása A mérési eredményeink feldolgozásához sokszor szükségünk van az árapály hatás számítására. A számításhoz gyakran ún. árapály katalógusokat használnak. megemlítjük ezek közül Tamura (1987) árapály katalógusát, amely 1200 tagot tartalmaz, illetve Hartmann-Wenzel (1995) árapály katalógusát, amelyben már 12935 tag szerepel. Legújabban a KSM03 katalógust használják, amelyben 28806 tag szerepel (Kudryavtsev, 2004). A szakterületen leginkább a néhai H-G. Wenzel által írt ETERNA 3.4 árapály számító program használata terjedt el. A program képes a közvetlen árapályhatás számítása mellett azt a hatást is számítani, amely az óceáni árapály teher földkéregre gyakorolt hatása miatt jelentkezik a mért mennyiségekben.
Szilárd földkéreg árapálya A GNSS mérések, különösen pedig a GNSS ún. PPP (precíz pontmeghatározási) eljárása esetében nem elhanyagolható a szilárd földkéreg árapálya. Ezt például a D. Milbert által írt solid programmal számíthatjuk ki (http://home.comcast.net/~dmilbert/softs/solid.htm). Az alábbi ábrán Budapesten egy adott napra láthatjuk a program által a helyi derékszögű koordináta-rendszerben számított árapály-komponenseket. A függőleges összetevő közel 30 cm-t változik a nap folyamán, és csak azért nem észleljük a hagyományos GNSS-mérésekkel, mert általában relatív helymeghatározást végzünk, vagyis egy (vagy több) bázisállomáshoz képest vett relatív koordináta-különbségeket határozunk meg.
20. ábra. A földkéreg árapálya Budapesten 2011.11.29.-én
13
Űrgravimetria (GOCE) Említettük, hogy az űrgravimetriai és űrgradiometriai mérések feldolgozásában a felhasználók munkáját olyan programok segítik, amelyek nagyban megkönnyítik az ilyen jellegű mérési, számítási eredmények elemzését. Az egyik ilyen eszköztár, a szabadon elérhető GUT (GOCE User Toolbox) a GOCE ESA mesterséges hold ún. Level 2 szintű adatainak feldolgozásához készült (http://earth.esa.int/gut/). Alapja a mérési eredmények gömbfüggvény-sora, és képes a geoid, a nehézségi rendellenességek, illetve a függővonal-elhajlások számítására. Figyelembe tudja venni az árapály hatását, megadhatók a feldolgozásnál figyelembe veendő maximális gömbfüggvény együttható fokszám, továbbá a referencia ellipszoid paraméterei. Képes különböző szűrők alkalmazására, illetve az eredmények elemzésére, összehasonlítására is.
14
