II.2.  Adatasszimiláció és inicializáció

Bölöni Gergely

Horányi András

II.2.1. Adatasszimiláció

Az adatasszimiláció célja

A légköri mozgásokat és állapotváltozásokat leíró prognosztikai egyenletrendszer vegyes feladat, azaz olyan parciális differenciálegyenlet-rendszer, amelynek a megoldásához kezdeti- és határfeltételek megadása szükséges. Az adatasszimiláció célja a HTER megoldásához szükséges kezdeti feltételek (analízis) megadása az alkalmazott numerikus modell rácshálózatán. A HTER-ben előforduló nem-lineáris advekciós tagok miatt a megoldás (a numerikus előrejelzés) rendkívüli érzékenységet mutat a kezdeti feltételekre ( II.6. ábra), ezért kulcsfontosságú ezek minél pontosabb megadása az adatasszimiláció során.

II.6. ábra. Fent: a „valóságot” reprezentáló verifikációs analízis (tengerszinti légnyomás), amely egy mély ciklont rajzol ki Nyugat-Európa fölött. Balra, lent: a verifikációs időpontra vonatkozó előrejelzés (tengerszinti légnyomás) egy „jó” kezdeti feltételből indítva. Jobbra, lent: a verifikációs időpontra vonatkozó előrejelzés (tengerszinti légnyomás) egy „rossz” kezdeti feltételből indítva. Az alsó ábrákhoz tartozó kezdeti feltételek ugyanazon megfigyelések felhasználásával készültek, azonban különböző adatasszimilációs módszerekkel (bal: variációs asszimiláció, jobb: optimális interpoláció, a módszerek magyarázatát lásd később). Amíg a két kiindulási mező csak csekély mértékben tér el egymástól, addig az abból indított előrejelzések már komoly különbségeket mutatnak.

A modell kezdeti feltételeinek minél pontosabb meghatározása érdekében az adatasszimiláció során célunk, hogy minden rendelkezésünkre álló információt (adatforrást) figyelembe vegyünk. Ezek a következők:

  1. A légköri megfigyelések térben és időben szabálytalanul elhelyezkedő sokasága (ezek általában megfelelő megbízhatóságú adatok, de nem állnak kellő számban rendelkezésre, azaz az adatasszimilációs probléma csupán a megfigyelések figyelembevételével erőteljesen alulhatározott lenne);

  2. Az ún. háttérmező, ami a kívánt kezdeti időpontra vonatkozó rövidtávú numerikus modell előrejelzés (ez a modell rácspontjaiban áll rendelkezésre, ugyanakkor kevésbé megbízható, mint a megfigyelések, lévén, hogy előrejelzésről van szó);

  3. A légköri egyensúlyi folyamatokról alkotott ismereteink (pl. a tömeg- és a szélmező geosztrofikus igazodása).

Az első két információforrást az adatasszimilációs eljárások képesek közvetlenül figyelembe venni, míg a harmadikból csak közvetett módon profitálnak. A légköri egyensúlyi folyamatokkal kapcsolatos ismereteink felhasználására (a 3. információ-forrás) a következő, az inicializációról szóló II.2.2. alfejezetben világítunk rá. Az adatasszimiláció során célunk a fenti információ sokaság matematikailag optimális felhasználása a kezdeti feltételek meghatározásához. Az „optimalitás” szó itt azt jelenti, hogy az egyes információk megbízhatóságuk függvényében járulnak hozzá a kezdeti mezőhöz (pl. pontosabb mérőeszközből származó információ nagyobb, míg a kevésbé pontos műszerekből származó mérések kisebb súlyt kapnak). Az egyes információk (megfigyelések és háttér) megbízhatósága (vagy hibája) sajnos nem ismert egzaktul, ezért ezeket becsülni kényszerülünk. Az adatasszimiláció egyik fontos részterülete a megfigyelések és a háttér megbízhatóságának becslése, a jegyzetben azonban erre csak érintőlegesen térünk ki, és az egyszerűség kedvéért alapvetően ismertnek tekintjük a felhasznált információk megbízhatóságát (hibáját).

Elméleti alapok

Az adatasszimiláció feladatát a matematikai statisztikára támaszkodva a következőképpen lehet egyszerűen megfogalmazni. Tekintsünk egyetlen térbeli pontot, amelyben az x változó xt valódi értékére szeretnénk xa becslést (analízist) adni a rendelkezésünkre álló két megfigyelés, y1 és y2 alapján:

.

(II.7.)

Az adatasszimilációs módszerekkel lényegében arra az egyszerű kérdésre keressük a választ, hogy mi legyen az az F függvény, amely a lehető legjobb becslést adja xt-re a fentiek ismeretében. Az F függvény meghatározására a két leggyakrabban alkalmazott közelítést mutatjuk be. Ezek 1.) a legkisebb négyzetek módszere és 2.) a maximum likelihood módszer.

1.) Legkisebb négyzetek módszere

Tudjuk, hogy az y1 és y2 megfigyelések hibával terheltek: y1 = xt + ε 1  illetve y2 = xt + ε 2 , ahol a mérési hibákat ε 1-gyel, illetve ε 2-vel jelöljük. Feltesszük, hogy

  • a mérések torzítatlanok (nincs szisztematikus hibájuk, azaz hosszabb idő átlagában sem felülbecslést és sem alulbecslést nem mutatnak): E(ε 1) = E(ε 2) = 0 ;

  • ismerjük a mérési hibák szórásnégyzetét (azaz a hibák nagyságának mérőszámát): σ 12 = E(ε 12), σ 22 = E(ε 22) ;

  • a mérési hibák korrelálatlanok (azaz a mérések függetlenek egymástól): E(ε 1ε 2) = 0  , ahol E a várható értéket jelöli.

Keressük az xa becslést a megfigyelések lineáris kombinációjaként!

   

(II.8.)

azaz ebben az esetben a feladat tulajdonképpen a k1, k2 együtthatók meghatározása. További elvárásunk az xa becsléssel szemben, hogy torzítatlan legyen (azaz a mérésekhez hasonlóan ne legyen szisztematikus hibája):

    

(II.9.)

illetve hogy a négyzetes hibája minimális legyen:

   

(II.10.)

A becslés torzítatlanságából és a megfigyelési hibákra tett fenti feltételekből következik, hogy az együtthatók összege egy: k1 + k2 = 1. Ennek az összefüggésnek és a megfigyelésekre tett feltételezések (torzítatlanság és korrelálatlanság) kihasználásával könnyen kifejezhetők a keresett együtthatók a megfigyelések ismert hiba szórásnégyzeteinek függvényében:

  

(II.11.)

A k1, k2 súlyokat visszahelyettesítve a becslési hiba szórásnégyzetét megadó fenti egyenletbe adódik, hogy:

    

(II.12.)

Nevezzük el a megfigyelések és az analízis hiba szórásnégyzeteinek reciprokait a fenti egyenletben megbízhatóságnak (hiszen minél nagyobb a megfigyelési hiba, annál kisebb a mérés megbízhatósága). Így a fenti egyenlet rendkívül szemléletesen világít rá az adatasszimiláció lényegére, azaz arra, hogy az analízis megbízhatóságát minden egyes megfigyelés növeli (mivel a megfigyelések megbízhatóságai mind 1-nél kisebb pozitív számok és az analízis megbízhatósága az egyes mérések megbízhatóságainak az összege).

2.) Maximum likelihood módszer

Tekintsük az xt-t ismert eloszlású valószínűségi változónak, ismert f sűrűségfüggvénnyel. Tekintsük továbbá az y1, y2 megfigyeléseket ugyanilyen eloszlású statisztikai mintának és xa -t (xt becslését) egy ismeretlen becsülni kívánt paraméternek. Az xa becslést ekkor az ún. likelihood függvény (minta együttes sűrűségfüggvénye) maximalizálásával kapjuk meg (lásd pl. Dévényi és Gulyás, 1988):

  

(II.13.)

Tegyük fel, hogy xt és az y1, y2 minta normális eloszlású és hogy a megfigyelésekhez tartózó hiba szórásnégyzetek (σ 12 és σ 22) most is ismertek. Ekkor a likelihood függvényt az alábbi alakban írhatjuk fel:

(II.14.)

Az exponenciális függvény negatív kitevője miatt a fenti likelihood függvény értéke akkor lesz maximális, ha a kitevő értéke minimális, azaz xa értékét az alábbi feladat megoldásával kapjuk:

   

(II.15.)

Másodfokú egyenlet lévén, a fenti összefüggés minimumát úgy kapjuk meg, hogy xa szerinti deriváltját nullával tesszük egyenlővé:

    

(II.16.)

Világosan látszik tehát, hogy a maximum likelihood módszerrel is ugyanarra az eredményre jutunk mint a legkisebb négyzetek módszerével. Természetesen a maximum likelihood módszer esetében is igaz tehát, hogy a becslés megbízhatóságát az egyes megfigyelések megbízhatóságának összege adja.

Fontos megjegyezni, hogy a fenti módszereknél tett feltételezések (megfigyelések torzítatlansága, korrelálatlansága és normális eloszlása) nem feltétlenül teljesülnek minden megfigyelési típus és meteorológiai változó esetén (például a műholdas megfigyelések adatai nem függetlenek egymástól). A fenti feltételezésekben foglaltak az adatasszimilációs módszerek korlátaiként tekinthetők, és amennyiben a valóság eltér a feltételezésekben foglaltaktól, a statisztikai becslés nem tekinthető megbízhatónak. A gyakorlatban ezért rendkívül fontos a fenti feltételezések vizsgálata, és amennyiben szükséges, a megfigyelések szisztematikus hibáinak szűrése az asszimilációs feladat elvégzése előtt. Szintén fontos, hogy az asszimilációban résztvevő változókat normális eloszlású változókká kell transzformálni (amennyiben az eredeti változók nem tekinthetőek normális eloszlásúnak). Természetesen az asszimiláció elvégzése után az így kapott mesterséges változók visszatranszformálandóak az eredeti változókba. A megfigyelések szisztematikus hibájából adódó problémákat a II.7. ábra szemlélteti.

II.7. ábra. Szél megfigyelések szisztematikus hibájának hatása 12 órás (bal) és 24 órás (jobb) szélsebesség előrejelzések (ms-1) minőségére. A szisztematikus hibával terhelt és a kontroll előrejelzések négyzetes hibáinak különbsége a földrajzi szélesség–magasság keresztmetszeten (a földrajzi hosszúság szerint átlagolt értékek). Az ábrán jól látható, hogy a megfigyelések szisztematikus hibái nagymértékben rontják az előrejelzések minőségét (sárgás-pirosas árnyalatok; a keresztek a négyzetes hiba különbségének szignifikanciáját jelölik).

Adatasszimilációs módszerek alkalmazása a gyakorlatban

A gyakorlatban alkalmazott numerikus modellek esetében a kezdeti (analízis) mezőt 3-dimenziós rácson kell megadnunk, amely a II.1. fejezetben adott becslésnek megfelelően 107 nagyságrendű pontból áll. A gyakorlatban a megfigyelések száma 105 nagyságrendű (például felszíni, rádiószondás, repülőgépes, műholdas, radar, GPS, wind-profiler megfigyelések). A megfigyelések relatív ritkasága miatt (a ritkaság a meghatározni kívánt modell rácspontokhoz képest értendő), nélkülözhetetlen a háttérmező (a numerikus modell korábbi időpontból induló előrejelzése az analízis időpontra vonatkozóan) felhasználása az asszimiláció elvégzéséhez. A gyakorlatban használt adatasszimilációs rendszerek esetében tehát az analízis mező a háttérmező megfigyelésekkel való korrekciója révén áll elő. Mint látni fogjuk, a korrekció mértékét a háttérmező és az aktuális megfigyelések megbízhatósága alapján határozzuk meg. Fontos megjegyezni, hogy a megfigyelések általában nem rácspontokban, hanem térben szabálytalanul helyezkednek el, ezért a háttérmező és a megfigyelések értékeinek összehasonlításához legalább egy térbeli interpolációt kell végrehajtanunk, amelyet az ún. megfigyelési operátor hivatott elvégezni. Amennyiben a megfigyelés, amelyet asszimilálni kívánunk nem közvetlenül a modell valamely állapothatározójára vonatkozik, a térbeli interpoláción felül a megfigyelési operátor hivatott leírni a megfigyelt és a modell-változó között fennálló fizikai kapcsolatot is. Ez a kapcsolat gyakran bonyolult és nem-lineáris (pl. a műholdas megfigyelések esetében a sugárzás-átviteli egyenlet, a radar reflektivitás megfigyelések esetében pedig a radar egyenlet írja le a megfigyelt és a modellváltozók közötti kapcsolatot). A fentiek alapján a gyakorlatban elvégzendő adatasszimilációs feladat matematikai leírásához a következő jelöléseket célszerű bevezetnünk:

  • xa : analízis mező (107 méretű vektor);

  • xb : háttérmező (107 méretű vektor);

  • : megfigyelések (105 méretű vektor);

  • : megfigyelési operátor, amely nem-lineáris (a 107 méretű xb vektort 105 méretű vektorrá alakítja át, amelyeket a megfigyelési pontokban értelmezünk);

  • : a H megfigyelési operátor lineáris közelítése (107 x 105-es mátrix), azaz H(x) – H(xb) = H(x – xb) ;

  • : háttérmező hiba kovariancia mátrixa (107 x 107 méretű mátrix);

  • : megfigyelések hiba kovariancia mátrixa (105 x 105 méretű mátrix, amely korrelálatlan megfigyelési hibák esetében diagonális).

1.) Optimális interpoláció

Megmutatható, hogy a fenti jelöléseket felhasználva a legkisebb négyzetek módszerével adott becslés több dimenzióban és háttérmezőt is felhasználva az alábbi formában írható fel (lásd pl. Gandin, 1963, Bouttier és Courtier, 1999):

    

(II.17.)

ahol K az ún. Kalman súlymátrix (Kálmán Emil Rudolf, 1930-ban született magyar matematikus és mérnök után elnevezve; Kalman, 1960), amely az alábbi alakban írható föl:

(II.18.)

A fenti első egyenlet szemléletesen azt fejezi ki, hogy az analízist a háttérmező pontosításával kapjuk, méghozzá a megfigyelésekből jövő korrekciók (inkrementumok) lineáris kombinációinak hozzáadásával. A lineáris kombináció súlytényezői a K (10x 105 méretű mátrix) elemeiből adódnak. A fenti becslést gyakran nevezik BLUE (best linear unbiased estimate)-nak az irodalomban, jelezvén, hogy a K mátrix fenti összetételéből kapjuk a legjobb torzítatlan becslést. A szemléletesség kedvéért érdemes egyetlen térbeli pont esetére felírni a K mátrixot. Ekkor H = I, R = σ o2 (megfigyelési hiba szórásnégyzete), B = σ b2  (háttér hiba szórásnégyzete), azaz K = σ b2 / (σ b2 + σ o2) , amely visszaadja az előző fejezetben kapott k1, k2 súlyok alakját (a különbség az, hogy ott két megfigyelés alapján becsültünk, a jelen egyszerűsített példa pedig egy megfigyelt és egy háttér információból való becslést mutat be). Az optimális interpoláció gyakorlati megvalósításánál problémaként merül fel, hogy a nagy dimenziójú, 10x 105 méretű (HBHT + R)-1  inverz mátrix kiszámítása meglehetősen költséges. A dimenzió csökkentése céljából az optimális interpolációt a gyakorlatban rácspontonként szokás elvégezni, úgy, hogy csak meghatározott távolságon belüli megfigyeléseket használunk fel az adott rácsponti analízis érték meghatározásához. Így a fent említett invertálandó mátrix ~10 x 10-es dimenziójúra csökken. A feladat rácspontonkénti megoldásának hátránya, hogy a szomszédos pontok konzisztenciája sérülhet az analízisben (ugyanis az egyes rácspontokban más és más megfigyeléseket használunk az adatasszimiláció során).

2.) Variációs asszimiláció

Megmutatható, hogy a maximum likelihood becslés több dimenzióban és háttérmezőt is felhasználva az alábbi skalár függvény minimalizálására vezethető vissza (lásd pl. Bouttier és Courtier, 1999, Kalnay, 2003):

(II.19.)

A fenti függvényt variációs veszteségfüggvénynek nevezzük. A veszteségfüggvény x szerinti minimumában kapjuk meg az xa analízist. A veszteségfüggvény első tagja a Jb háttér tag, amely szemléletesen a keresett analízis háttérmezőtől vett távolságát méri súlyozva a háttérmező hiba kovariancia mátrixával. A veszteségfüggvény második tagja a Jo megfigyelési tag, amely szemléletesen a keresett analízis mező megfigyelésektől vett távolságát méri súlyozva a megfigyelési hiba kovariancia mátrixszal.

Mint ahogy az előző fejezet egyszerű példájában is megmutattuk a legkisebb négyzetek módszerének és a maximum likelihood módszer ekvivalenciáját, a jelenleg tárgyalt nagy dimenziószámú esetben is megmutatható, hogy lineáris megfigyelési operátor esetén [azaz ha H(x) = H(x)] , a fenti veszteségfüggvény minimalizálása pontosan ugyanazt az eredményt adja, mint az optimális interpoláció, azaz a (II.17–18.) BLUE egyenletet. Az ekvivalencia a fenti veszteségfüggvény xa szerinti gradiensének számításával és annak zérussá tételével látható be formálisan. A gyakorlatban a nagy dimenziószám miatt a veszteségfüggvény zérus-helyének keresése nem kapható meg triviálisan egy gradiens számításból, hanem azt numerikusan közelítjük minimum-kereső algoritmusok alkalmazásával.

A minimum-kereső algoritmusok iteratív módon a gradiens értéke és iránya alapján javasolnak olyan újabb és újabb x vektorokat, amelyekre a J(x) veszteségfüggvény egészen addig csökken, amíg a gradiens értéke közel nulla nem lesz. A fenti veszteségfüggvény alakjából közvetlenül nem látszik, de belátható, hogy a gradiens számításához szükség van a H operátor linearizáltjára (H) és adjungáltjára (HT) is. A veszteségfüggvény értékének számításakor problémaként merül fel a nagy dimenziójú B-1  és R-1  mátrix inverzek számítása is, hasonlóan az optimális interpolációhoz. Az R-1  számítása a megfigyelések korrelálatlanságának feltétele mellett egyszerűvé válik, hiszen ebben az esetben R  diagonális mátrix a σ oi2 szórásnégyzetekkel a főátlójában, amelyeknek a reciprokát véve kapjuk az inverz mátrixot. A B-1  explicit számolását a χ = B-1 / 2 (x – xb) vektor transzformációval szokás áthidalni, amelynek eredményeképpen a háttér tag az egyszerű χTχ alakot ölti. A B-1 / 2 mátrixszal való szorzás pedig a gyakorlatban helyettesíthető az ismert háttér hiba szórásnégyzetek és kovarianciák négyzetgyökeit tartalmazó vektorszorzásokkal.

Ha a megfigyelések az analízis időpontjára vonatkoznak, a variációs veszteség függvény minimalizálása egy három-dimenziós analízis mezőt ad meg, ezért a variációs módszert három-dimenziósnak nevezzük (3DVAR). Abban az esetben azonban, ha a megfigyeléseket egy időintervallumon belül szeretnénk figyelembe venni, a veszteségfüggvény minimalizálása a modell egy időbeli trajektóriáját fogja megadni, azaz egy négy-dimenziós mezőt kapunk (4DVAR). Ekkor az  y – H(x) inkrementum számításához a numerikus modell időbeli integrálására is szükség van, azaz a H operátor magában foglalja az időbeli fejlődést leíró modell operátorát is (M) , valamint a gradiens számításához szükség van a numerikus modell linearizáltjára (M) és adjungáltjára (MT). A modell linearizáltjának és adjungáltjának létrehozása jelentős kódfejlesztést igényel, ugyanis a nem-lineáris modell kódjának soronkénti deriválását, illetve annak transzponálását igényli. Ez a kódfejlesztés igen munkaigényes, hiszen a nemlineáris modell folyamatos fejlesztése mellett a lineáris és adjungált változatokat is folyamatosan aktualizálni kell. Végül de nem utolsósorban érdemes azt is megjegyezni, hogy minden adott numerikus modell és használandó asszimilációs időablak esetén külön vizsgálni kell azt, hogy teljesül-e a linearitás feltétele (érvényes-e a tangens-lineáris közelítés), azaz hogy a H(x) – H(xb) = H(x – xb) összefüggés jó közelítéssel fennáll-e.

3.) Kalman filter

A Kalman filter (vagy Kálmán-szűrő) jelentőségének megértéséhez fontos hangsúlyozni, hogy az asszimilációban a B háttér hiba kovariancia mátrix szerepe rendkívül nagy, mégpedig azért mert a benne foglalt térbeli és változók közötti kovarianciák teszik lehetővé a megfigyelésekben rejlő információk kiterjesztését a modell rácsra ( II.8. ábra) (Fisher, 2001).

II.8. ábra. Az analízis és a háttérmező különbsége (analízis inkrementum) a hőmérsékletre (K) a kb. 500 hPa-os szinten egy megfigyelés felhasználásával, amely 1 K-nel tér el a háttértől a megfigyelési pontban. Az ábra jól érzékelteti, hogy a Ljubljana fölötti hőmérséklet megfigyelés hatása horizontálisan kiterjed az analízisben a B háttér hiba kovariancia mátrix hatására. A kísérlet az AROME nem hidrosztatikus modellel készült 2,5 km-es horizontális rácstávolságot használva.

A fent tárgyalt optimális interpoláció és variációs asszimiláció esetében B-t időben konstans ismert mátrixnak tekintettük. A valóságban azonban a háttér hibák (rövidtávú előrejelzési hibák) mértéke és struktúrája meglehetősen függ az időjárási helyzettől (hiszen a modellek bizonyos helyzetekben jobb, más helyzetekben rosszabb előrejelzéseket produkálnak). A Kalman filter elmélet a háttér hiba kovariancia mátrix időfüggésének leírását célozza meg a (II.17–18.) BLUE egyenleteket kiegészítve a B mátrix időbeli fejlődését leíró egyenlettel (lásd pl. Kalman, 1960; Kalnay, 2003). Legyen az analízis hiba kovariancia mátrixa A és a modell hiba kovariancia mátrixa Q (modell hiba alatt HTER diszkretizációjából és a fizikai parametrizációk bizonytalanságából adódó hibát értjük). Jelöljük az analízis időbeli indexét k val. A Kalman filter egyenleteket ekkor a következőképpen írhatjuk fel:

  • az analízis számítása a k időpontban: xak = xbk + K[yk – H(xbk)]  a Kk = BkHT(HBkHT + R)-1 súlymátrixot használva;

  • az analízis hiba kovariancia mátrix számítása: Ak = (I – KkH)  ;

  • a numerikus modell integrálása k-ból + 1 időpontba: xbk+1 = M(xak) ;

  • a háttér hiba kovariancia mátrix időbeli fejlődésének (k-ból + 1 időpontba) leírása Bk+1 = MAkMT +  .

A fenti egyenletek megoldása után a + 1-edik időpontban ismert megfigyelések mellett minden információnk rendelkezésre áll a BLUE becslés kiszámításához, beleértve a + 1-edik időpontra jellemző háttér hiba kovariancia mátrixot is. Fontos azonban megjegyezni, hogy a Kalman filter elmélet alkalmazása az operatív gyakorlatban nem lehetséges, mivel az MAkMT számítása 10x 107 integrálást igényel a numerikus modell linearizáltjával, illetve adjungáltjával, amelynek gyors (operatív gyakorlatban igényelt) számítására nincs lehetőség a jelenleg alkalmazott számítógépek kapacitása mellett. Mindazonáltal, az ensemble módszerek alkalmazásával lehetőség nyílik arra, hogy a B mátrix időbeli fejlődését leíró egyenletet ne explicit módon, hanem egyszerűsítve, dimenziócsökkentéssel oldjuk meg. Ezeket az egyszerűsített módszereket ensemble Kalman filtereknek nevezzük [erre a témakörre nem áll módunkban részletesebben kitérni a jegyzet keretein belül, de az érdeklődők például Evensen (2007) könyvéből tájékozódhatnak].

Adatasszimilációs ciklus

Az analízisek és időbeli modellintegrálások (háttér előrejelzések) egymásutánját adatasszimilációs ciklusnak nevezzük. Az adatasszimilációs ciklus úgy is felfogható, mint a modell folytonos trajektóriájának a megfigyelésekhez való illesztése, amelyből tetszőleges időpontban (de leginkább azokban a szinoptikus időpontokban, amikor a legtöbb megfigyelés rendelkezésre áll) indíthatunk hosszabb távú, a gyakorlatban is felhasználható előrejelzéseket. Az adatasszimilációs ciklus lényegét az alábbi II.9. ábra szemlélteti sematikusan.

II.9. ábra. Az adatasszimilációs ciklus illusztrációja. Az analízis mindig a megfigyelés és a háttér közé esik. Az analízisből rövid-távú modell előrejelzésekkel kapjuk a háttérbecslést a következő analízis időpontra.

Az ábrán látszik, hogy a háttér előrejelzések az analízisből indulnak, amely mindig a megfigyelés és az előző háttér közé esik, hiszen mint azt az előző fejezetekben láttuk, az adatasszimilációs módszerek a megfigyelésekkel korrigálják a háttérmezőt. Az adatasszimilációs ciklus jelentősége abban áll, hogy segítségével a háttérmező egy adott analízis időpontban az összes korábbi, a ciklusban figyelembe vett, megfigyelés információtartalmával is bírni fog, azaz a ciklus ismerni fogja a légkör korábbi „történetét” (vagy legalábbis annak becslését), illetve más szavakkal azt is mondhatjuk, hogy minden egyes múltbeli megfigyelés hatással van minden egyes jövőre vonatkozó előrejelzésre.

II.2.2. Inicializáció

Az inicializáció célja

Az időjárás előrejelzésben alkalmazott numerikus modellek megalkotásakor általában valamilyen egyszerűsítő közelítéssel élünk a HTER-re nézve (pl. a ma még gyakran használt hidrosztatikus közelítés, vagy a múltban előszeretettel használt kvázi-geosztrofikus és sekélyvíz közelítések). Ha magukat a prognosztikai egyenleteket nem is egyszerűsítjük le drasztikusan (nem-hidrosztatikus modellek), a numerikus megoldáshoz szükséges diszkretizáció során mindenképpen információ veszteséggel kell számolnunk. Más szóval képtelenek vagyunk a valós légkör minden egyes folyamatát a modell rácson leírni. Mint ahogy azt az előző fejezetben megfogalmaztuk, az adatasszimiláció során megfigyelt információkkal tápláljuk a numerikus modellt a kezdeti feltételeken keresztül. Fontos tény, hogy a meteorológiai megfigyelési rendszerek „belemérhetnek” olyan légköri folyamatokba is, amelyeket a fent említett egyszerűsítések miatt az alkalmazott numerikus modell nem képes leírni (az egyenletek egyszerűsítése vagy a térbeli felbontás elégtelensége miatt). Amennyiben ilyen megfigyelést asszimilálunk, az zajhoz vezet a kezdeti mezőben, illetve a modell-egyenletek integrálása során.

Az, hogy mely megfigyelések vezetnek zajhoz és melyek hordoznak hasznos információt (jelet), az adott numerikus modell egyszerűsítéseitől (az alkalmazott közelítésektől) és a rácsfelbontástól függ. Egy ~10 km-es horizontális rácsfelbontású hidrosztatikus modellben például egy zivatarban elvégzett szél megfigyelés zajként fog jelentkezni, hiszen ebben az esetben sem az egyenletek, sem a rácsfelbontás nem engedi meg a zivatarcellák explicit leírását a modellben, azaz a méréshez hasonló nagy vertikális sebességeket. Ebben a példában a kezdeti mezőben jelentkező zaj főként a szél- és tömeg-mező kiegyensúlyozatlanságából adódik az adott hidrosztatikus közelítés mellett. A légköri hullám-mozgásokat (Práger, 1992) alapul véve úgy is fogalmazhatunk, hogy a fenti esetben olyan a zivatarokra jellemző nagy frekvenciájú gravitációs hullámokat is figyelembe veszünk a kezdeti mezőben, amelyek leírására az adott numerikus modell nem képes. Az ilyen hullámok csak zajként lehetnek tehát jelen az integrálás során, torzítva a megoldást.

Az asszimilált megfigyelések gyakran a modellben parametrizált fizikai (termodinamikai) változókkal sem konzisztensek. A fenti példát alapul véve, egy hidrosztatikus modellben, amelyben a mély-konvekciót parametrizáljuk, egy a kezdeti mezőben előforduló nagy vertikális sebesség irreálisan nagy kondenzációt és csapadékképződést okozhat a kezdeti időlépcsőkben, ezzel elrontva az előrejelzést. A kezdeti zajnak betudható, a fizikai parametrizációkból származó instabilitásokat „spin-up”-nak vagy felpörgésnek nevezzük. Spin-upot okozhat a modellekben az is, ha a prognosztikus csapadékelemek hiányoznak a kezdeti mezőből (ezek asszimilációja igen nehézkes és ezért a legtöbb modell kezdeti feltételeiben ezek nem adottak), hiszen ekkor néhány időlépcsőnyi integrálásra van szükség, hogy a parametrizációs séma létrehozza ezeket.

Korlátos tartományú modellek esetén (vagyis amikor a modell nem az egész Földre, hanem csak egy jól behatárolt kisebb területre fókuszál) gyakori, hogy a kezdeti mezőt nem adatasszimiláció útján, hanem egy gyengébb felbontású meghajtó modell kezdeti feltételeinek interpolációjával határozzuk meg. Az interpoláció a meghajtó modell felbontása alatt óhatatlanul zajt generál, hiszen ez alatt a felbontás alatti térskálákról a meghajtó modellünk semmiféle információval nem rendelkezik.

Az inicializáció célja, hogy a kezdeti mezőben valamilyen módon megkülönböztesse a modell számára hasznos információt (a jelet), illetve a zajt, majd az utóbbit kiszűrje. Ez a zajszűrés meglehetősen fontos, hiszen szélsőséges esetben a kezdeti zaj hatására a modell numerikus megoldása instabillá válhat, azaz a valóságban elő nem forduló értékeket jelezhetünk előre. Általában a hatékony numerikus sémák alkalmazásának köszönhetően a numerikus instabilitások nem végzetesek inicializáció alkalmazása nélkül sem (a modell lefut), viszont az eredmény torzulhat, főként az integrálás kezdeti szakaszában (II.10. ábra).

A ma alkalmazott hidrosztatikus modellekben az inicializációs feladat tulajdonképpen a gyorsan terjedő gravitációs hullámok kiszűrésére szűkül le, hiszen ennek a hullám-mozgásnak a leírására a hidrosztatikus közelítés nem ad módot. Az egyre nagyobb felbontású nem-hidrosztatikus modellek esetében az inicializációra egyre kevésbé van szükség, hiszen ezek a modellek a gravitációs hullám-mozgásokat is képesek leírni, azaz ha van is zaj a kezdeti mezőben, azt egyre nehezebb elválasztani az értékes információtól. Ezekben a modellekben az inicializáció fő célja a spin-up kiküszöbölése vagy lerövidítése, ami annál is inkább elengedhetetlen, mert ezek a nagy felbontású modellek legtöbbször az ultra-rövidtávú előrejelzési időszakra koncentrálnak, azaz kulcskérdés, hogy az előrejelzés eredményes legyen már a modellintegrálás első óráiban is.

II.10. ábra. Az inicializáció hatása. A teljes modell tartományra kiátlagolt felszíni nyomás tendencia jó indikátora a gravitációs hullámok által keltett zajnak. Az inicializáció szűri ezt a zajt, amely leginkább a modellintegrálás kezdeti szakaszára jellemző. A fenti ábra az ALADIN/CHAPEAU modell hidrosztatikus változatának futtatásával készült 8 km-es horizontális rácstávolsággal és interpolált kezdeti feltételeket használva az ECMWF (European Centre for Medium-Range Weather Forecasts) globális modelljéből.

Alkalmazott inicializációs módszerek

Az alkalmazott inicializációs módszerek közül az 1970-es években elterjedt normál módus inicializációt (NMI) (Leith, 1980) és az 1990-es években teret hódító digitális filter inicializációt (DFI) (Lynch és Huang, 1992) mutatjuk be tömören.

1.) Normál módus inicializáció

A Normál módus inicializáció az alkalmazott modell fő hullám-megoldásainak (normál módusainak) lassú (kis frekvenciájú) és gyors (nagy frekvenciájú) hullámokra való felosztásán alapszik, amelyekből természetesen csak a lassú, az adott modell-egyenletek által leírt hullámokat kívánjuk megtartani a kezdeti mezőben. A normál módusok és frekvenciájuk kiszámításához az adott egyenlet-rendszer linearizálása szükséges, hiszen a normál módusok által reprezentált tiszta hullám-mozgások csak periodikus, lineáris rendszerekben vannak jelen (ezek vizsgálatával reméljük azonban, hogy a teljes, nem-lineáris egyenletrendszert is tudjuk jellemezni, feltételezve, hogy a nem-lineáris mozgásformák a linearizált rendszerben előforduló hullám-mozgások szuperpozíciói) (Práger, 1992). A linearizálást a kis perturbációk módszerével végezhetjük el, és a linearizált egyenletrendszert a következő alakban írhatjuk fel:

(II.20.)

ahol x’ az állapotvektor kis perturbációja (egy feltételezett alapállapottól vett kis eltérése), L pedig egy konstans mátrix, amely a linearizált modell operátora. Feltesszük, hogy a normál módusok az alábbi hullám-függvény alakban írhatóak fel:

   

(II.21.)

ahol X az amplitúdó és ω a frekvencia. A fenti normál módus alakot komponensenként visszahelyettesítve a linearizált egyenletbe és kiírva az L mátrix elemeit, egy lineáris egyenletrendszerhez jutunk, amelyet megoldhatunk az ω frekvenciára (Temperton, 1987; Kalnay, 2003). A normál módus inicializáció alkalmazásával így kapjuk meg a diszperziós relációt, vagyis az adott modell normál módusaihoz tartozó légköri hullám-mozgások frekvenciáját (pl. Rossby- és gravitációs-hullámok frekvenciái). Az inicializáció a gyors, az adott modellben túl nagy frekvenciájú módusok amplitúdójának lenullázásával történik meg a kezdeti mezőben.

2.) Digitális filter inicializáció

A digitális filter alkalmazásánál a zajt valamilyen határértéknél nagyobb frekvenciájú légköri hullámokkal azonosítjuk. A zaj szűrése ilyen módon egy alul-áteresztő (alacsony frekvenciákat megtartó és nagy frekvenciákat elimináló) digitális szűrő alkalmazásával valósulhat meg. A digitális szűrő alkalmazásakor a modell-változók időbeli sorozatát képezzük, majd a sorozat elemeit súlyozzuk (lineáris kombinációjukat képezzük):

 

(II.22.)

A változók időbeli sorozatát (xn) a modell hátra-, majd előre-integrálásával állítjuk elő úgy, hogy a sorozat a kezdeti időpontra nézve centrált legyen. A hátrafele-integrálásnál (t = –NΔt-ig) a fizikai parametrizációkat nem tudjuk alkalmazni a termo-dinamikai folyamatok irreverzibilitása miatt, azonban az előre-integrálás a teljes fizikai parametrizációs csomag alkalmazásával történhet (t = +NΔt-ig). A kezdeti időpontra vonatkozó inicializált mezőt (x*) a hn együtthatók megfelelő megválasztásával nyerjük. Az együtthatók meghatározására Huang és Lynch (1992) valamint Lynch és Huang (1992) a következő formulát javasolja:

 

(II.23.)

ahol Θ a frekvencia és Θ c az alkalmazott frekvencia határérték, amely fölött szűrünk. A fenti integrálást elvégezve kapjuk a hn súlyokra a következő egyszerű alakot:

  

(II.24.)

A hn együtthatók Fourier-sorát képezve kapjuk meg az ún. válaszfüggvényt:

.

(II.25.)

amely visszaadja az idealizált alul-áteresztő szűrő fent megfogalmazott tulajdonságát, miszerint Θ c alatt az értéke 1, a felett pedig 0. A gyakorlatban a Fourier-sor számításakor csonkítást alkalmazunk (csak véges intervallumban számoljuk ki a fenti összeget, mint ahogyan adott véges hosszúságú időablakban hozzuk létre az xn adatsort is, amelyen a szűrést végezzük). A csonkítás következménye, hogy az alul-áteresztő szűrőnk nem lesz többé ideális, azaz a nagy frekvenciák felé haladva a Θ c környezetében, egy sima átmenettel veszi fel a válaszfüggvény a 0 értéket az 1 helyett ().

II.11. ábra. A digitális szűrő (DFI) válasz függvénye idealizált, illetve valós esetben. A Lanczos ablak (Lynch és Huang, 1992) alkalmazásával jobb eredmény érhető el a DFI alkalmazásakor, ennek részleteire azonban a jegyzetben nem térünk ki.

II.2.3. Feladatok

Feladat 1: Egyetlen térbeli pontban szeretnénk becsülni az xt paramétert két megfigyelés, y1 és y2 alapján, amelyek normális eloszlású mintát alkotnak. Tudjuk, hogy a megfigyelések hibái torzítatlanok és korrelálatlanok, valamint ismerjük szórásnégyzetüket: σ 12 és σ 22. Adjuk meg a becslés alakját a legkisebb négyzetek módszerével vagy a maximum likelihood becsléssel. Mutassa meg, hogyan jutott el a becslés alakjához.

Megoldás 1:

Feladat 2: Egyetlen térbeli pontban szeretnénk becsülni az xt paramétert két megfigyelés, y1 = 5 és y2 = 10 alapján, amelyek normális eloszlású mintát alkotnak. Tudjuk, hogy a megfigyelések hibái torzítatlanok és korrelálatlanok, valamint ismerjük szórásnégyzetüket:σ 12 = 1 és σ 22 = 2. Adjuk meg a becslést, illetve hibájának szórásnégyzetét. Ábrázoljuk a megfigyelésekhez tartozó sűrűségfüggvényeket, illetve a becslés sűrűségfüggvényét.

Megoldás 2: xa = 6, σ a2 = 0,8 ; az ábrázoláshoz olyan programot kell írni (legegyszerűbb GnuPlot-ban), amely ábrázolja a normális eloszláshoz tartozó sűrűségfüggvényt:

.

A sűrűségfüggvény alakjában várható érték (m) helyére az adott megfigyelés (vagy az analízis) értékét, míg a szórás és szórásnégyzet helyére az adott megfigyelés (vagy az analízis) szórását, illetve szórásnégyzetét kell behelyettesíteni. Az eredmény:

Feladat 3: Vezessük le a Kalman filter elméletben megadott Bk+1 = MAkMT + Q egyenletet a háttér hiba kovariancia mátrix időfejlődésére, az alábbi definíciókat és feltételezéseket felhasználva:

  • legyen M lineáris modell, amelyre xbk+1 = xak;

  • háttér hiba definíciója: ε bk+1 = xtk+1  – xbk+1  ;

  • analízis hiba definíciója: ε ak = xtk  – xak  ;

  • modell hiba definíciója: ε M k+1 = xtk+1  – xtk  ;

  • a modell hiba és az analízis hiba korrelálatlanok (a fenti definíciókból adódik);

  • háttér hiba kovariancia mátrix definíciója: Bk+1 = [ε bk+1, (ε bk+1)T] ;

  • analízis hiba kovariancia mátrix definíciója Ak = [ε ak, (ε ak)T] ;

  • a modell hiba kovariancia mátrix definíciója: Qk+1 = [ε M k+1, (ε M k+1)T] .

Megoldás 3: Fejezzük ki a háttér hibát, a modell hibát (ε M k+1 = xtk+1  + M xtk ) leíró és az időbeli fejlődést (xbk+1  = M xak) leíró egyenletekkel. Ebből adódik, hogy: ε bk+1 = ε ak + ε M k+1, azaz a háttér hiba az analízis hiba időbeli fejlődéséből és a modellhibából tevődik össze. A háttér hiba ezen alakját behelyettesítve a háttér hiba kovariancia mátrix definíciójába, megkapjuk a keresett egyenletet.