1 Zagadnienie odwrotne Krzysztof MarkowiczInstytut Geofizyki, Wydział Fizyki Uniwersytet Warszawski Szkoła letnia sieci badawczej Poland-AOD, Warszawa 5-8 lipca 2017r.

2 Wprowadzenie W teledetekcji atmosferycznej zagadnienie odwrotne jest kluczowym problemem w niemal wszystkich metodach pomiarowych aerozoli oraz gazów śladowych. Z całego widma promieniowania analizowane są takie przedziały spektralne, w których promieniowanie elektromagnetyczne oddziałuje z materią (jest rozpraszanie lub absorbowane) W ogólności sygnał S odbierany przez detektor może być zapisany w postaci:

3 gdzie, T jest badanym obiektem, F reprezentuje zaś pewną funkcję. S=F(T), gdzie, T jest badanym obiektem, F reprezentuje zaś pewną funkcję. Funkcja ta opisuje procesy radiacyjne w ośrodku i jest najczęściej funkcją nieliniową. Funkcja odwrotna F-1 opisuje nam badany obiekt zgodnie z relacją: T=F-1(S). W większości przypadków, z jakimi spotykamy się w zagadnieniach atmosferycznych funkcji odwrotnej F-1 nie możemy wyznaczyć. W takim przypadku poszukujemy pewnych parametrów naszego obiektu, które najlepiej odpowiadają rejestrowanemu promieniowaniu. Target (T) Signal (S) S=F(T) T=F-1(S)

4 Podstawowym problemem, jaki napotykamy w metodach odwrotnych jest brak jednoznacznego rozwiązania.Wynika to z faktu, że nasz problem jest najczęściej problemem niedookreślonym, ze względu na większą liczbę parametrów, które chcemy wyznaczać w stosunku do liczny niezależnych obserwacji. Np. w przypadku wyznaczania profilu temperatury zazwyczaj mamy pomiary w kilkunastu czy w kilkudziesięciu kanałach spektralnych, zaś naszą niewiadomą jest funkcja ciągła. Dlatego temperaturę powietrza wyznacza się tylko dla kilku lub kilkunastu warstw powietrza. W innym przypadku jeśli dana warstwa ośrodka składający się z różnych gazów oraz substancji ciekłych i stałych, ich kombinacja może dawać ten sam efekt radiacyjny (w pewnym obszarze spektralnym) dla różnych koncentracji składników atmosfery.

5 Poza niejednoznacznością pojawia się problem stabilności rozwiązania oraz problem uzyskania tego rozwiązania. Niestabilności rozwiązania mogą pojawia się np. ze względu na błędy obserwacyjne lub błędne założenia poczynione na temat własności fizycznych badanego ośrodka. W wielu metodach teledetekcyjnych problem odwrotny sprowadza się do równania Fredholma pierwszego rodzaju gdzie funkcja f(x) może opisywać np. profil pionowy temperatury atmosfery, K(x) jest jądrem, zaś gi wartościami mierzonej radiancji w „i” kanałach spektralnych.

6 Rozważmy ciało doskonale czarne o temperaturze TRozważmy ciało doskonale czarne o temperaturze T. Dokonujmy pomiaru natężenia (radiancji) promieniowania emitowanego przez to ciało w dowolnej odległości. Zakładamy jednak brak atmosfery miedzy detektorem a ciałem. Wyznaczenie temperatury tego ciała (zgodnie ze wzorem Plancka) wymaga pomiaru natężenia promieniowania jedynie dla pojedynczej długości fali. T

7 W przypadku gdy między detektorem a ciałem znajduje się izotermiczna atmosfera o temperaturze TA oraz grubości optycznej τ wówczas promieniowanie docierające do detektora zależy od 3 zmiennych (nie uwzględniając długości fali). Tak, więc musimy mierzyć promieniowanie dla co najmniej trzeych długościach fali, tak aby wyznaczyć niewiadome wielkości. Dodatkowo, własności optyczne dla tych tych długości fali muszą się różnic znacząco. W atmosferze temperatura zmienia się z wysokością więc sytuacja jest znacznie bardziej skomplikowana. T TA τ

8 W tym przypadku interesuje nas wyznaczenie stosunku zmieszania pary wodnej oraz CO2 w atmosferze przy założeniu, że nie zmieniają się one z wysokością. Jeśli wybierzemy przedział spektralny gdzie występuje znaczna absorpcja promieniowania przez parę wodną oraz CO2 oraz stosunek ich współczynnika absorpcji słabo zmienia się z długością fali to łatwo wyobrazić sobie sytuacje, że ten sam efekt radiacyjny mogą powodować różne proporcje zawartości pary wodnej oraz CO2. Problem jest wiec źle uwarunkowany. xH20 =const xCO2=const TA=const Problemu tego możemy uniknąć przez odpowiedni wybór obszaru spektralnego, gdzie stosunek współ. absorpcji obu gazów zmienia się znacząco z długością fali

9 Jeśli uwzględnić niedokładności pomiarowe równania Fredholma sprowadza się do postacigdzie błędy i mogą powodować znacznie zmiany profilu funkcji f(x). Czułość rozwiązania na błędy pomiarowe jest rzeczą bardzo niepożądana. Można ją jednak minimalizować poprzez odpowiedni dobór obszaru spektralnego dla którego wykonujemy pomiary promieniowania.

10 Przykład – czułość na błędy pomiaroweRozważmy przypadek, gdy mamy tylko dwie obserwowane wartości promieniowanie I1 oraz I2. Dyskretyzując równanie Fredholma mamy: Załóżmy, że wagi W mają następujące wartości: W1,1=W1,2 =1 W2,1=2 W2,2 = Zaś wartości promieniowania wynoszą: I1 =2 I2 = Wówczas uzyskujemy poszukiwane wielkości: B1=1, B2 =1 Następnie niech wartości I2 będzie nieco inna ze względu na niepewności pomiarowe i wynosi I2 =4. Uzyskujemy wówczas: B1=2, B2=0.

11 Uwagi do przykładu Problem niestabilności rozwiązania pojawił się ze względu na wartości wag Wij. Wagi W1,1 i W1,2 oraz W2,1 i W2,2 są równe sobie co oznacza, że własności optyczne atmosfery w tych 2 kanałach są identyczne. Dlatego, więc pomiar dla drugiej długości fali nie zawiera dodatkowej informacji, a układ równań dwóch równań jest układem zredukowanym do jednego. W celu wyznaczania informacji atmosferycznych należy wybierać kały spektralne różniące się transmisją atmosfery

12 Zagadnienie odwrotne a asymilacja danych w modelach prognoz pogodyZ matematycznego punktu widzenia problem zagadnienia odwrotnego jest równoznaczny problemowi asymilacji danych w numerycznych prognozach pogody. W obu przypadkach problem jest na ogół źle postawiony gdyż liczba obserwacji jest mniejsza od liczby wyznaczanych wielkości fizycznych.

13 Przez y (y1,y2,…,ym) oznaczmy wektor obserwacji, zaś x (x1,x2,…,xn) wektor wyznaczanych (niewiadomych) wielkości (wektor stanu). Przez  oznaczamy wektor błędów obserwacji. Relacje pomiędzy wektorem obserwacji i wektorem stanu zapisujemy w postaci: gdzie F(x) oznacza model fizyczny (model do przodu – forward model). Używamy terminu model gdyż powyższy związek jest często określony przez skompilowane relacje fizyczne zapisywane w postaci numerycznej.

14 Funkcja wagowa W wielu rozważaniach wygodnie jest rozważać problem liniowy. Dokonujemy linearyzacji modelu fizycznego w otoczeniu pewnego stanu referencyjnego xo. Macierz K (m x n) oznaczamy funkcją wagową. Macierz ta nie koniecznie musi być kwadratowa. W przypadku gdy m m>n mamy nadmiarową liczbę obserwacji.

15 3-wymiarowa analiza wariacyjna : 3D-VARW metodzie 3D-Var poszukujemy wektora analizy xa, który minimalizuje skalarną funkcję kosztu. Zdefiniowana jest ona przez odległość pomiędzy wektorem stanu x a wektorem pierwszego przybliżenia xb mnożoną przez wagę będąca odwrotnością kowariancji błędu i odległość pomiędzy wektorem stanu x, a wektorem obserwacji yo mnożoną przez odwrotność kowariancji błędów obserwacyjnych. W metodzie 3D-Var minimalizacji dokonujemy w przestrzeni wektora stanu.

16 3-wymiarowa analiza wariacyjna : 3D-VARRozważamy funkcję koszu oraz jej gradient w postaci: Minimalizacja wariacyjnej funkcji kosztu (na podstawie 2-wymiarowego modelu). Kwadratura funkcji kosztu ma kształt paraboloidy (w tym przypadku) z wartością minimalna dla optymalnej wartości analizy xa. Algorytm poszukiwania wartości minimalnej sprowadza się do poruszania po krzywej funkcji kosztu w kierunku największego gradientu funkcji.

17 W praktyce punkt startowy minimalizacji zwany pierwszym przybliżeniem (first guess) jest często wybierany na podstawie informacji a priori (background) xb. Nie jest to wybór obowiązkowy jednak należy pamiętać o różnicy pomiędzy informacją a priori, która jest używana w definicji funkcji kosztu od pierwszego przybliżenia, które jest używane do inicjalizacji procedury minimalizacyjnej. Jeśli minimalizacja jest zadowalająca to wynik analizy nie zależy istotnie od wyboru wartości startowej. Jednak zawsze zależy od informacji a priori. Znaczącym problemem analizy 3D-Var jest konieczność znalezienia metody pozwalającej wyznaczyć macierz kowariancji B, która określa błędy informacji a priori dla każdej pary zmiennych modelu. W większości przypadków macierz kowariancji błędu związana z obserwacjami jest przekątna macierzą blokową lub macierzą diagonalą.

18 Łatwo zauważyć, że przekątna macierz blokowa implikuje, iż funkcja kosztu Jo jest sumą N skalarnych funkcji kosztu Jo,i każdej zdefiniowanej dla podmacierzy Ri oraz odpowiadającej Hi oraz yi. Rozbicie funkcji kosztu Jo staje się użytecznym narzędziem do badania zachowania metody 3D-Var ze względu na każdą obserwację (jej wartość i dopasowanie do wektora stanu x) Dodatkowo pozwala to na wymuszenie słabszych więzów (ograniczeń) przez dodanie dodatkowego czynnika w funkcji kosztu Jc. Prowadzi to jednak do warunku wstępnego co utrudnia i komplikuje minimalizację.

19 Teoria Bayesa W podejściu Bayesa używamy pojęcia prawdopodobieństwa do opisu naszej wiedzy na temat wektora stanu oraz obserwacji. Definiujemy: P(x) - gęstość praw-sta (pdf) wektora stanu x. P(x)dx jest prawdopodobieństwem przed wykonaniem obserwacji, że wektor stanu znajduje się w przedziale (x,x+dx). P(y) - pdf obserwacji P(x,y) - pdf złożone x i y. P(x,y)dxdy oznacza prawdopodobieństwo, że wektor x znajduje się w przedziale (x,x+dx) zaś y w przedziale (y.y+dy). P(y|x) - pdf warunkowe wektora y dla danego x. Oznacza, że P(y|x)dy jest prawdopodobieństwem, że wektor obserwacji y znajduje się w przedziale (y,y+dy) gdy wektor stanu x przyjmuje określoną wartość P(x|y) – analogicznie jak powyższej

20 Rodgers, 2000

21 Twierdze Bayesa : opisuje prawdopodobieństwo warunkowe Koncepcyjne przybliżenie problemu odwrotnego: Przed wykonaniem obserwacji mamy wiedzę a priori w postaci rozkładu gęstości prawdopodobieństwa (pdf-u). Proces obserwacyjny jest utożsamiany jako mapowanie wektora stanu w przestrzeni obserwacji przy użyciu modelu (forward model) Teoria Bayesa opisuje formalizm procesu odwrotnego do powyższego mapowania i wyznaczania pdf-u aposteriori poprzez poprawianie pdf-u a priori przez pdf obserwacji.

22 Rozważmy problem liniowyBłędy pomiarowe  mogą być często przybliżane rozkładem Gaussa stąd wyrażenie na P(y|x) ma postać: gdzie c1 jest stałą zaś R jest macierzą kowariancji błędów pomiarowych

23 Ma ono rozkład Gaussa więc może być zapisane w postaci:Podobnie można zdefiniować pdf wektora stanu. Jednak w tym przypadku przybliżenie rozkładem Gaussa jest mnie realistyczne aczkolwiek wygodne do opisu. gdzie xa jest a priori znanym stanem x, zaś B odpowiadającą mu macierzą kowariancji. Podstawiając i wykorzystując twierdzenie Bayesa dostajemy związek na pdf a posteriori Ma ono rozkład Gaussa więc może być zapisane w postaci: gdzie oznacza oczekiwaną wartość

24 Porównując czynniki kwadratowe w x otrzymujemy:Co daje: Analogicznie równanie liniowe w xT: Upraszczając czynnik xT ponieważ równanie musi być spełnione dla każdego x oraz podstawiając za S-1 otrzymujemy:

25 alternatywnie

26 Wyznaczanie rozkładu wielkości aerozolu jako przykład problemu odwrotnegoRozkład wielkości n(r) możemy rozbić na dwie części; wolno f(r) oraz szybko h(r) zmienną, gdzie h(r) ma na przykład postać rozkładu Junge Naszym zadaniem jest wyznaczenie funkcji f(r) zakładając współczynnik refrakcji m. Równanie to sprowadza się do równania Fredholma pierwszego rodzaju jeśli przyjmiemy, że g=() zaś

27 Rozważmy równanie Fredholma w postaci:gdzie K(r) jest funkcją wagową (jądrem). W praktyce, ponieważ mamy tylko skończona ilość mierzonych parametrów gi powyższy problem jest źle postawiony nawet jeśli funkcja wagowa K(r) oraz wartości mierzone gi pozbawione są niepewności. Rozważmy równanie Fredholma w postaci: i=1,2,…,M M jest liczbą obserwacji spektralnych wielkości g Dla wygody poszukiwać będziemy rozwiązanie w postaci: gdzie fj są nieznanymi współczynnikami, zaś wj oznacza funkcje ortogonalne. Podstawiając dostajemy:

28 gdzie i-1,2,…,M Musimy wyznaczyć fj(j=1,…,N) korzystając z obserwacji gi(i=1,…,M). Wprowadzamy oznaczenia: Wyjściowe równanie sprowadza się do równania wektorowego Jak wiadomo macierz odwrotna istnieje tylko wtedy gdy M=N oraz gdy detA0. Jeśli więc M=N to nasze rozwiązanie jest postaci:

29 Z reguły macierz A nie może być odwrócona dlatego używa się metody najmniejszych kwadratów. Różnicę lewej i prawej strony powyższego równania zapisujemy w postaci: i=1,…,M Musimy więc zminimalizować wielkość Poprzez przyrównaniu wszystkich pochodnych względem fk(k=1,2,…,n) do zera

30 lub w formie macierzowejRozwiązanie zagadnienia odwrotnego w powyżej formie jest niestabilne. Niestabilności związane ze źle postawionym problemem wyjściowym nie są jedyne. Istotnym wkład wnoszą błędy kwadratur używane do obliczeń elementów macierzy Aij, błędy obcięcia numerycznego oraz przede wszystkim błędy pomiarowe wielkości gi. W praktyce gi nigdy nie jest znane i dlatego musi być zapisane w postaci: Niepewności i wpływają na niejednoznaczność rozwiązania fi, która może być usunięta poprzez narzucenie dodatkowego warunku pozwalającego wybrać jedno z możliwych fi

31 Rozpatrujemy wyrażeniegdzie  jest współczynnikiem wygładzającym mówiącym jak silnie rozwiązanie fi jest zmuszane do zbiegania do zadanej wartości (w tym przypadku do wartości średniej a ogólnie do wektora informacji a priori). Stosując metodę najmniejszych kwadratów mamy Równoważny zapis macierzowy gdzie macierz H ma postać

32 Postać macierzy H wynika ze wzoruRozwiązanie problemu odwrotnego: Wprowadzony dodatkowy warunek ma na celu zbliżać rozwiązanie do pewnej klasy rozwiązań określonych przez Rozwiązanie na może być wyznaczane na podstawie danych historycznych. gdzie I jest macierzą jednostkową

33 Warunek wygładzania rozwiązania można konstruować również przez wyrażenia:1) (pierwsza pochodna) (druga pochodna) które nie zawierają wyrażenia Stosując przedstawioną po wyżej metodę rozwiązywania problemu odwrotnego dla grubości optycznej aerozolu zakładaliśmy, że znamy współczynnik refrakcji. Założenie to jest bardzo silne i może prowadzić do znacznych błędów, które w tej metodzie wchodzą do wielkości i Wartość współ. refrakcji wpływa na zmienność efektywnego przekroju czynnego na ekstynkcje i tak cześć rzeczywista odpowiada ze przesuwanie kolejnych maksimów w zależności od parametru wielkości x zaś część urojona za wygładzanie oscylacji rezonansowych.

34 Fitowanie rozkładu log-normalnegoZakładamy, że mamy dwu-modowy rozkład log-normalny. Mamy do wyznaczenia 6 parametrów swobodnych (przy założeniu współczynnika refrakcji) N1,N2, rm1, rm2, 1,2 Możemy liczbę niewiadomych zredukować o 1i 2 na podstawie informacji klimatycznych dla odpowiednich modów. Dodatkowo przy dużej licznie kanałów spektralnych AOT możemy fitować współczynnik załamania światła.