Metody Przetwarzania Danych Meteorologicznych Wykład 5

1. Metody Przetwarzania Danych MeteorologicznychWykład 5 Krzysztof Markowicz Instytut Geofizyki UW kmark@igf.fuw.edu.pl

2. Asymilacja Danych Meteorologicznych • To proces znajdywania reprezentacji modelu „najbardziej zgodnej” z obserwacjami meteorologicznymi. • Asymilacja danych jest procesem polegającym na uzupełnieniu niekompletnych danych w pewnym modelu opisującym system dynamiczny • Głównym problemem asymilacji danych są zróżnicowane źródła danych meteorologicznych: pomiary in-situ, pomiary satelitarnych, radarowe, lidarowych wykonywane w różnym czasie i w różnych miejscach na kuli ziemskiej. Np. Jak asymilować dane o odbiciowości radarowej do modeli numerycznych prognoz pogody? • Asymilacja danych jest bardzo złożonym procesem i można ją rozdzielić na dwa procesy; • (a) pierwszy etap polega na sprawdzeniu jakości danych - to jest etap skomplikowany, ale prosty do zrozumienia. Chodzi o wyeliminowanie oczywistych błędów pomiarowych; • (b) drugim i ważniejszym etapem jest uzgodnienie danych i wykorzystanie informacji dostępnej z poprzednich godzin czy z poprzednich dni. Jest to skomplikowany proces polegający na całościowej analizie danych dość często z wykorzystaniem modelu prognoz numerycznych.

13. Zarys historii • Pierwsze metody asymilacji danych nosiły nazwą analizy obiektywnej (np. algorytm Cressmana, 1959). Było to podejście przeciwne do metod subiektywnych opierających się na analizie pól meteorologicznych przez synoptyka • Metody obiektywne wykorzystywały przybliżenia oparte na prostych technikach interpolacyjnych. Są to metody 2 lub 3D. • Podobne metody w 4D (z czasem) noszą nazwę "nudging" (np. w modelu MM5) Bazują one na idei relaksacji Newtona, która ma za zadanie dodanie we właściwych członach równań dynamicznych modelu różnicy pomiędzy obliczonymi zmiennymi meteorologicznymi a wartościami obserwacyjnymi. Człony te maja ujemny znak co pozwala utrzymywać obliczone wartości zmiennych (wektora stanu) w sąsiedztwie wartości obserwowanych (wektora obserwacji).

14. Przełomowe stało się wprowadzenie do asymilacji danych statystycznej (optymalnej) interpolacji przez L. Gandin (1963). Metoda jego bazowała na idei Kolmogorov. Jest ona pewnym rodzajem analizy regresyjnej, która wykorzystuje informacje o przestrzennym rozkładzie funkcji kowariancyjnych błędów pierwszego przybliżenia (wcześniejszej prognozy) oraz rzeczywistego pola meteorologicznego. Chociaż funkcje te nigdy jednak nie są znane to stosuje się szereg ich przybliżeń. Optymalna interpolacja odpowiada zredukowanej wersji filtru Kalmana, gdy macierze kowariancji nie są liczone z równań dynamiki a wyznaczane na podstawie kolejnych przybliżeń pół. • Wprowadzenie filtru Kalmana do asymilacji 4D jest bardzo trudnym zadaniem i wymaga rozwiązania dodatkowo ogromnej ilości równań: ~N*N~1012, gdzie N=Nx*Ny*Nzjest rozmiarem wektora stanu, typowy rozmiar domeny obliczeniowej : Nx~100, Ny~100, Nz~100

15. W celu uniknięcia tego problemu stosuje się specjalne wersje filtru Kalmana np Reduced-Rank Kalman filters (RRSQRT) Kolejnym rozszerzeniem metod 4D-Var jest wykorzystanie przybliżenia wariacyjnego (Le Dimet and Talagrand, 1986 wcześniej opracowane przez G. Marchuk) • Uwzględnia ono fakt, iż pola meteorologiczne spełniają równania dynamiki zawarte w model prognoz pogody. Tym samym minimalizują funkcjonał opisujący ich różnice w stosunku do obserwacji. • Jak pokazał Lorenz, 1986 wszystkie wspomniane powyżej metody 4D są równoważne (przy pewnych założeniach ) minimalizacji pewnej funkcji kosztu. • W praktyce założenia te nigdy nie są jednak spełnione.

16. Gwałtowny rozwój metody asymilacji danych meteorologicznych do numerycznych prognoz pogody związany jest z dwoma zagadnieniami: • wykorzystaniem obecnych obserwacji do poprawienia jakości prognozy dla różnych skal przestrzennych (od skali planetarnej do skali ulic miasta) oraz czasowych • wykorzystanie równych typów dostępnych obserwacji (sodary, radary,lidary), które gwałtownie rozwijają się • Ciągle aktualnym pytaniem pozostaje: czy asymilacja danych pozwoli na przezwyciężenie trudności z prognozowaniem stanów atmosfery?

17. Do czego używamy asymilacji danych • Globalne i lokalne mapy synoptyczne(Primary -Constrained Product) • Niezmierzone wielkości(Primary - Derived Product) • Wiatr ageostroficzny, pionowe profile, • Wyznaczane wielkości: • Ruch pionowy / dywergencja,cyrkulacja residualna, diabatyczne oraz radiacyjne własności, ozon troposferyzny • Inicjalizacja prognozy • Poprawki radiacyjne w metodach teledetekcyjnych “Background,” (a prioriprofile dla metod teledetekcyjnych) • Monitoring • Kalibracja przyrządów • Ocena jakości obserwacji • Walidacja i ocena modeli

18. Główne strategie asymilacji danych • Wykorzystuje się dwa główne podejścia: • asymilacja sekwencyjna, w której wykorzystuje się obserwacje wykonane przed rozpoczęciem analizy, która jest częścią „real-time assymilation systems”. • asymilacja niesekwencyjna w której również obserwacje z „przyszłości są wykorzystywane. Ma to miejcie w tzw. re-analize. • Podział ze względu na metody • nieciągła (intermittent), w której obserwacje asymilowane są małymi porcjami co jest technicznie wygodną techniką • ciągła (continuous), w której asymilowane są dane obserwacyjne z znacznie dłuższego okresu czasu. Pozawala to na korekcje wektora stanu i jego wygładzanie co jest fizycznie bardziej realistyczne.

21. Definicje podstawowych wielkości • Wektor stanu (x) opisujący stan układu. Jego związek z rzeczywistym stanem układu zależy od dyskretyzacji co z matematycznego punktu widzenia związane jest z wyborem bazy. • Wyróżniamy więc: • xt – prawdziwy wektor stanu (true state vector) • xb- wektor informacji a priori lub „background” będącym oszacowaniem rzeczywistości przez wykonaniem analizy • xa– poszukiwany wektor (analiza) • Problem analizy sprowadza się do znalezienia poprawek do stanu podstawowego xb

22. Terminologia: informacje a priori / a posteriori • Informacje a priori – zawierają dane o systemie przed wykonaniem obserwacji. Są to najczęściej dane klimatologiczne (background) lub wynik wcześniejszej prognozy. • Informacje a posteriori (po fakcie)– określają naszą znajomość systemu (wektora stanu) po wykonaniu obserwacji.

23. Obserwacje • Wektor obserwacji: y • Operator obserwacji H, pozwala dokonać transformacji od wektora stanu (zdefiniowanego w modelu) do wektora obserwacji. W praktyce jest to operator związany interpolacją od dyskretnych punktów siatki modelu do nieregularnej sieci obserwacyjnej. • H(x) – są wartościami obserwacyjnymi jakie uzyskalibyśmy gdyby wektor stanu był idealny i model pozbawiony był błędów. • Głównym zadaniem asymilacji danych jest minimalizacja różnicy pomiędzy obserwacjami a wektorem stanu: y-H(x) Różnice tą możemy liczyć dla x=xb oraz dla x=xa (mówimy wówczas o analizie residualnej)

24. Błędy i niepewności • Ze względu na błędy zawarte w obserwacjach i w polu pierwszego przybliżenia (np. dane klimatologiczne) musimy założyć pewien model błędów. • Do tego celu wykorzystamy funkcji gęstości prawdopodobieństwa (pdf). • Dla danego pierwszego przybliżenia xb tuż przed wykonaniem asymilacji mamy jeden wektor błędu, który oddziela to pole od pola rzeczywistego Jeśli moglibyśmy wykonać to bardzo wiele razy w tych samych warunkach ale z różna realizacją błędów moglibyśmy obliczyć statystyki (średnia, wariancje, histogram εb). Dla bardzo dużej liczby realizacji moglibyśmy wyznaczyć gęstość prawdopodobieństwa pdf i z niej wyznaczać wszystkie statystyki błędów.

25. Błędy zerowego przybliżenia (background errors): • Związana z nim kowariancja: Nie zawierają one błędów dyskretyzacji • Błędy obserwacyjne: Zawiera błędy powstałe w czasie wykonywania obserwacji, ale również związanie z konstrukcja operatora H a więc zawiera błędy dyskretyzacji. H(xt) nie jest perfekcyjnym obrazem prawdziwego stanu. • Błędy analizy: • Błąd wektora stanu dany jest przez ślad macierzy kowariancji. • Naszym zadaniem jest minimalizacja tego błędu • Średni błąd (obciążenie) określa błąd systematyczny. Niezerowa wartość wskazuje na problemy w systemie asymilacji danych i może być związanym z dryfem modelu czy błędami systematycznym obserwacji.

26. Macierz kowariancji • Jeśli wektora stanu ma wymiar n wówczas macierz kowariancji ma wymiar n x n, elementy diagonalne są wariancjami dla poszczególnych zmiennych. Elementy pozadiagonalne są kowariancją poszczególnych elementów wektora. • Marzcież jest dodatnio określona: xTAx>0 dla x>0, oraz ma dodatnie własności własne. • W przypadku gdy dokonujemy liniowej transformacji P wektora stanu macierz kowariancji B po transformacji ma postać: PBPT.

27. Praktyczne wyznaczanie błędów • Statystyki błędów są funkcjami procesów fizycznych rządzących sytuacjami meteorologicznymi oraz własnościami sieci obserwacyjnej. • Zależą one również od naszej znajomości a priori błędów. • Generalnie mamy tylko jedna możliwości oszacowania statystyki błędów. Musimy założyć stacjonarności w czasie i jednorodność w przestrzeni statystyki błędów. Dzięki czemu dostajemy wiele realizacji błędów i możemy wyznaczamy empiryczne statystyki. Podejście takie ma sens klimatologiczny.

28. Analiza Cressman’a • To jedno z najprostszych podejść asymilacja danych, w którym zmienne modelu są ustawiane na podstawie obserwacji meteorologicznych w najbliższym ich otoczeniu. Po za tym obszarem wektor stanu modelu ustawiany jest na podstawie danych klimatologicznych lub wcześniejszej prognozy modelu.

29. Zakładamy, że składowe wektor stanu modelu opisywane są przez zmienne skalarne określone w punktach węzłowych modelu. xb – wektor stanu określony na podstawie klimatologii (background) lub wcześniejszej prognozy, xb(i) – jest tym samym wektorem prze- interpolowanym do punktu i, y(i) – wektor obserwacji (i=1,2,…,n), xa - wektor modelu określony w punktach siatki j, d – jest odległością pomiędzy punktami i oraz j, funkcja wagowa w(i,j) wynosząca 1 dla punktu siatki modelu (i=j) oraz malejąca z odległością osiągającą wartość zero poza tzw. promieniem wpływu (di,j >R). Istnieje wiele odmian metody Cressman’a, w których definiuje się różne postacie funkcji wagowej w.

30. W metodzie „successive correction” funkcja wagowa może mieć wartość mniejszą od jedności w punkcie siatki modelu (i=j) co oznacza, że zarówno wartość klimatyczna jak i obserwacyjna ma wkład do wartości osiąganej w tym punkcie siatki. Słabe strony metod Cressman’a • Jeśli posiadamy wcześniejszy wektor stanu modelu o wysokiej jakości i nie chcemy modyfikować przez słabej jakości dane obserwacyjne • Nie jest oczywiste jak oddalając się od punktu obserwacyjnego dokonać relaksacji do danych klimatycznych • Analiza powinna uwzględniać znane własności rzeczywistości (zależności pomiędzy zmiennymi, równowagę hydrostatyczna) a metoda ta nie uwzględnia tego • Błędy obserwacyjne mogą generować niefizyczne stany modelu.

31. Nasze oczekiwania • Asymilacje powinniśmy zacząć od stanu o wysokiej jakości (opartym np. na wcześniejszej prognozie) zwanym pierwszym (startowym) przybliżeniem (first guess) • Jeśli sieć obserwacji jest gęsta wówczas zakładamy, że prawdziwy stan znajduje się „blisko” średniej wartości z tych obserwacji. Musimy dokonać kompromisu pomiędzy pierwszy przybliżeniem a wartościami pochodzącymi z obserwacji. • Analiza powinna wygładzać nasze pole, gdyż wiemy, że taka jest cecha pól meteorologicznych. Gdy odchodzimy od punktu obserwacyjnego analiza powinna gładko przejść do pierwszego przybliżenia. • Analiza powinna uwzględniać znane własności fizyczne opisujące stan atmosfery. • Chcemy minimalizować różnice pomiędzy analizą a prawdziwym stanem.

32. Interpolacja statystyczna – metoda najmniejszych kwadratów Zakładamy: • zmienność operatora obserwacji H w otoczeniu pierwszego przybliżenia jest liniowa: H(x)-H(xb)=H(x-xb), H jest operatorem linowym • Nietrywialne postacie macierzy kowariancji B i R • Średnie błędy są zerowe:

33. Błędy nie są skorelowane: • Analiza liniowa: poszukujemy poprawek do pierwszego przybliżenia, które zależą liniowo od różnicy pierwszego przybliżenia i obserwacji. • Analiza optymalna: poszukujemy wektora stanu, który w sensie odchylenia średnio-kwadratowego jest najbliżej stanu rzeczywistego. • Z metody najmniejszych kwadratów otrzymujemy: K jest macierzą wagową Macierz kowariancji błędu w ogólnym przypadku dana jest wzorem: Dla metody najmniejszych kwadratów ma postać:

34. Jest ona równoważna metodzie optymizacyjno - wariacyjnej gdzie J jest funkcją kosztu analizy, Jbjest czynnikiem związanym z pierwszym przybliżeniem zaś Joz obserwacjami. Jeśli funkcje gęstości prawdopodobieństwa błędów pierwszego przybliżenie oraz obserwacji są gaussowskie wówczas xa jest estymatorem rzeczywistego stanu xt w sensie maksymalnego prawdopodobieństwa.

35. Dowód poprawności wzorów metody najmniejszych kwadratów • Minimalizacja funkcji kosztu odpowiada zerowej pochodnej funkcji kosztów dla optymalnego wektora stanu xa. Łatwo można pokazać iż postać ta jest identyczna ze wzorem pokazanym w metodzie najmniejszych kwadratów, gdyż

36. Realizacja metody najmniejszych kwadratów • W obecnych modelach wektor stanu x jest rzędu n=107 • Liczba obserwacji p=105 dla każdej analizy. • Dlatego problem z punktu matematycznego jest niedookreślony.

37. Uwagi do założeń • Założenie dodatnio określoności macierzy kowariancji B i R jest spełnione w „dobrze” postawionych problemach asymilacji. Jeśli B nie jest dodatnio określona transformujemy ją do bazy ortogonalnej. Co oznacza, że pierwsze przybliżenie jest idealne. Jeśli R nie jest dodatnio określona to macierz K jest dobrze określona a analiza będzie równa obserwacji w punkcie siatki. Jednak metoda wariacyjna w tym przypadku nie może być używana. • Średnie błędy przeważnie nie są zerowe. Jednak jeśli są znane mogą być odjęte od pierwszego przybliżenia oraz pola obserwacji. Jeśli nie są znane analiza nie będzie optymalna. Dlatego istotne staje się monitorowanie średniego odchylenia przybliżenia zerowego w czasie asymilacji. • Założenie, iż błędy nie są skorelowane jest najczęściej spełnione, ponieważ błędy pierwszego przybliżenia oraz błędy obserwacji są zupełnie niezależne. Jednak w przypadku metod odwrotnych używanych np. w obserwacjach satelitarnych może istnieć niezerowa korelacja ze względu na fakt iż w metodach odwrotnych wykorzystuje się informacje z pierwszego przybliżenia.

38. Uwagi do liniowości operatora H • Założenie liniowości jest potrzebne do wyprowadzenia wyrażenia na macierz K. W praktyce H może nie być liniowa ale możemy dokonać linearyzacji w sąsiedztwie wektora przybliżenia zerowego. • Bardziej ogólnie, możemy dokonać rozwiniecie w szereg Taylora • Operator H zwany jest stycznym • W przypadku metody najmniejszy kwadratów wymagamy aby: Problem nieliniowości operatora H nie jest związany z błędami obserwacyjnymi ale z błędami pierwszego przybliżenia, które w asymilacji sekwencyjnej są błędami wcześniejszej prognozy i zależą od zasięgu prognozy i jakości modelu.

39. Teoria Bayesa • W podejściu Bayesa używamy pojęcia prawdopodobieństwa do opisu naszej wiedzy na temat wektora stanu oraz obserwacji. • Twierdze Bayesa : opisuje prawdopodobieństwo warunkowe • Jeśli A jest zdarzeniem x=xt, B jest zdarzeniem y=yo wówczas rozkład prawdopodobieństwa a posteriori (po fakcie) wektora stanu dla nowej obserwacji yo wynosi:

40. Bayesowskie oszacowanie wektora stanu odpowiada maksymalnemu prawdopodobieństwu a posteriori zgodnie z poprzednim wzorem. Celem naszej analizy jest wyznaczenie stanu o maksymalnym prawdopodobieństwie a posteriori znając rozkład prawdopodobieństwa dla pierwszego przybliżenia (tła) oraz dla obserwacji. • Zakładamy, że rzeczywisty wektor stanu jest realizacją procesu losowego zdefiniowanego przez wielowymiarowy rozkład Gaussa.

41. Korzystając z tw. Bayesa otrzymujemy: Szukamy więc wektora stanu, dla którego prawdopodobieństwo warunkowe osiąga maksimum co odpowiada minimalnej wartości funkcji kosztu J.

42. Podsumowanie Mamy dwie możliwości zdefiniowania analizy statystycznej: (1) Gdy znamy macierze kowariancji błędów pierwszego przybliżenia oraz obserwacji i wyprowadzamy równania analizy wymagając aby całkowita wariancja błędów analizy była minimalna. (2) Gdy zakładamy gaussowskie rozkłady gęstości prawdopodobieństwa dla zerowego przybliżenia i wyprowadzamy równania analizy opisujący wektor stanu o maksymalnym prawdopodobieństwie. Oba przybliżenia prowadzą do matematycznie równoważnych algorytmów Z punktu widzenia numerycznego mają one różne własności.

43. Przykład – ilustracja metody najmniejszych kwadratów – przypadek skalarny • Chcemy oszacować temperaturę powietrza w pokoju na podstawie wskazań termometru o znanej dokładności σo(odchylenie standardowe). W wyniku pomiaru otrzymaliśmy wartość To. • Jeśli nie mamy żadnych innych dodatkowych informacji najlepsze oszacowanie temperatur powietrza wynosi oczywiście To z dokładnością σo . • Załóżmy, że posiadamy dokładne pomiary z dnia ubiegłego, które możemy traktować jako informację a priori (pierwsze przybliżenie) Tb , σb. • Nie trudno domyśleć się, że kombinacja liniowa wartości Tooraz Tb pozwoli nam na lepsze oszacowanie nieznanej temperatury rzeczywistej Tt. • Rozpatrzmy temperaturę będąca średnią ważona: oraz wariancja gdzie założyliśmy że błędy nie są ze sobą skorelowane

44. Jest to równoważne minimalizacji funkcji kosztu Zakładamy wartość k minimalizując wartość błędu zgodnie ze wzorem • W przypadku malej dokładności pomiaru (σo>>σb), k=0 • W przypadku dużej dokładności pomiaru • (σo <<σb), k=1 • Gdy (σo =σb), k=0.5 • W pozostałych przypadkach wartość analizy będzie średnią ważoną pomiędzy obserwacja a informacją a priori.

45. Wariancja analizy wyraża się wzorem błąd analizy jest zawsze mniejszy niż błędy obserwacji i informacji a priori razem wzięte.

46. Wyznaczanie macierzy kowariancji błędów • Poprawne wyznaczenie macierzy kowariancji błędów obserwacyjnych oraz zerowego przybliżenia jest kluczowe dla procesu asymilacji danych. • Poza wariancjami (wyrazy diagonalne macierzy korelacji) również współczynniki korelacji są istotne gdyż decydują o tym jak dane obserwacyjne będą wygładzane w przestrzeni modelu gdy istnieje niedopasowanie rozdzielczości modelu oraz gęstości sieci obserwacyjnej.

47. Wariancja błędów danych obserwacyjnych • Często zakłada się, że błędy wielkości pomiarowych nie są ze sobą skorelowane. Założenie to jest często racjonalne jednak w przypadku takich pomiarów jak: radiosondażowe czy satelitarne może nie być spełnione. • Powinno unikać się sytuacji gdy wartości obserwacyjne zawierają błędy statystyczne. • Generalnie jednak wyznaczanie macierzy kowariancji błędów obserwacyjnych R jest trudnym zadaniem. Dlatego w większości modeli macierz R jest diagonalna.

48. Wariancja błędów informacji a priori • Błędne oszacowanie wariacji błędów pierwszego przybliżenia prowadzi do zbyt małych lub zbyt dużych poprawek (innowacji) w procesie asymilacji danych w kolejnych (analysis increment). • W przypadku metody najmniejszych kwadratów jedynie względna wartość wariancji błędów obserwacyjnych i pierwszego przybliżenia jest istotna. • Jednak bezwzględne wartości wariancji mogą być istotne gdy dokonujemy kontroli jakości danych obserwacyjnych.

49. Korelacje błędów informacji a priori są istotne ze względu na: • Rozkład stacji obserwacyjnych W obszarach o rzadkiej gęstości sieci stacji pomiarowych kolejne poprawki analizy są określone przez strukturę macierzy kowariancji (dla pojedynczej obserwacji są one dane przez wielkość BHT). Tak wiec współczynniki korelacji macierzy B mówią jak informacje pochodzące ze stacji pomiarowej są propagowane na ich otoczenie. • Wygładzanie informacji W przypadku gęstych sieci obserwacyjnych istotne staje się wygładzanie informacji, które jednak zależy od samego pola meteorologicznego. Inaczej wygładzać powinno się obszary frontowe a inaczej gdzie mamy antycykloniczny charakter cyrkulacji.

50. (3) Różne typy równowag występujących w atmosferze W modelu zwykle mamy znacznie więcej stopni swobody niż w rzeczywistości. Np. w dużej skali przeważnie mamy równowagę hydrostatyczną, zaś w obszarach poza tropikalnych odchylenie wiatru od równowagi geostroficznej jest niewielkie. Tak, więc jedna zmienna obserwacyjna w modelu zwiera informacje o pozostałych, które są z nią powiązane. Np. pole przypowierzchniowego wiatru pozwala skorygować pole ciśnienia przy założeniu, że wiatr w pewnym obszarze jest geostroficzny.