Metody teledetekcyjne w badaniach atmosfery i oceanów. Wykład 7.

Metody teledetekcyjne w

badaniach atmosfery i oceanów.

Wykład 7.

Krzysztof Markowicz kmark@igf.fuw.edu.pl

2

Zagadnienie odwrotne

Z 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

odzyskiwanych parametrów (czy parametrów wektora stanu modelu w numerycznych prognozach pogody).

• Przez y (y₁,y₂,…,y_m) oznaczmy wektor obserwacji, zaś x (x₁,x₂,…,x_n) wektor odzyskiwanych 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ż związek ten może być tylko przybliżeniem lub oparty jest na

teorii fizycznej nie do końca jeszcze poznanej.





 F(x) y

4

Funkcja wagowa













 

 

 (x x ) K(x x )

x ) x ( ) F

x ( F

y _{o o o}

W wielu rozważaniach wygodnie jest rozważać problem liniowy. Dokonujemy linearyzacji modelu fizycznego w otoczeniu pewnego stanu referencyjnego x_o.

Macierz K (m x n) oznaczamy funkcją wagową. Macierz ta nie koniecznie musi być kwadratowa. W przypadku gdy

m<n problem jest niedookreślony (źle postawiony) m>n mamy nadmiarową liczbę obserwacji.

• Rozkład macierzy wagowej K według wartości osobliwych (dekompozycja na wartości singularne, SVD)

• Każdą macierz rzeczywistą K można przedstawić w postaci rozkładu SVD:

• U i V - macierze ortonormalne (U^-1 =U^T , V^-1 = V^T

  - macierz diagonalna, taka że = diag(σi), gdzie σi - nieujemne wartości szczególne (osobliwe) macierzy K, zwyczajowo uporządkowane nierosnąco.

V U

K  

6

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 przed jej wykonaniem

• 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 warunkowy 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

Rodgers, 2000

8

• Twierdze Bayesa :

opisuje prawdopodobieństwo warunkowe

Koncepcyjne przybliżenie problemu odwrotnego:

• Przed wykonaniem obserwacji mamy wiedzę a priori w postaci 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.

Zauważmy, że teoria Bayesa nie opisuje metody odwrotnej, która może być wykorzystana do uzyskania rozwiązania ale metodę

połączenia wszystkich metod odwrotnych w celu scharakteryzowania klasy możliwych rozwiązań i wyznaczenia pdf-u dla każdego z nich.

) y ( P

) x ( P ) x

| y ( ) P

y

| x (

P 

Rozważmy problem liniowy



  

  (y y) S^ (y y)

2 exp 1

| S

| ) 2 ( ) 1 y (

P _1/2 ^T _y¹

y 2 / n

 ^  ^  ^

 _{i i j j}

ij y y y y

S











 F(x) Kx y

Błędy pomiarowy  mogą być często przybliżane rozkładem Gaussa stąd wyrażenie na P(y|x) ma postać:

 

1

TS (y Kx) c

) Kx y

( )

x

| y ( P ln

2    

 ^_

gdzie c₁ jest stałą zaś S_ jest macierzą kowariancji błędów pomiarowych

10

• Podobnie można zdefiniować pdf wektora stanu. Jednak w tym przypadku przybliżenie rozkładem Gaussa jest mnie realistyczne aczkolwiek wygodne do opisu.

 ^  ^  ^

 _{a a} ^T

a x x x x

S





) x ( P ln

2    

 ^

gdzie x_a jest a priori znanym stanem x, zaś S_a odpowiadającą mu macierzą kowariancji.

Podstawiając i wykorzystując twierdzenie Bayesa dostajemy związek na pdf a posteriori

  



1 a T a 1

TS (y Kx) (x x ) S (x x ) c

) Kx y

( )

y

| x ( P ln

2       

 _^{ }

Ma ono rozkład Gaussa więc może być zapisane w postaci:

 

1

TSˆ (x xˆ) c )

xˆ x

( )

y

| x ( P ln

2    

 ^ gdzie oznacza

oczekiwaną wartośćxˆ

Porównując czynniki kwadratowe w x otrzymujemy:

x Sˆ x K

S x Kx

S K

x^{T T }_¹  ^T _a^1  ^{T 1}

1 a 1

1 KS K S

Sˆ^  _^  ^

Co daje:

Analogicznie równanie liniowe w x^T:

) xˆ ( Sˆ x )

x ( S x )

y ( S ) Kx

( ^T _^1  ^T _a^1  _a  ^{T 1} 

Upraszczając czynnik x^T ponieważ równanie musi być spełnione dla każdego x oraz podstawiając za S^-1 otrzymujemy:

xˆ ) S K

Sˆ K ( x

S y

S

K^T _^1  _a^1 _a  ^{T }_¹  _a^1

12

) Kx y

( Sˆ K )

S K

Sˆ K ( x

) x S

y S

K ( ) S K

Sˆ K ( xˆ

a 1

T 1

1 a 1

T a

a 1 a 1

T 1

1 a 1

T





































) Kx y

( ) S K

Sˆ K ( K S x

xˆ  _a  _a ^T _a ^T  _ ^1  _a

alternatywnie

• Rysunek obrazuje relacje pomiędzy kowariancją a priori obserwacji oraz kowariancją a posteriori w przypadku 3D wektora stanu oraz 2D wektora obserwacji. Duża elipsoida centrowana w xa opisuje kontur kowariancji a priori. Cylinder opisuje przestrzeń zgodności wektora stanu i obserwacji.

• Małą elipsoida przedstawia obszar zgodności informacji a priori oraz

obserwacji. Jej środek x nie pokrywa się z osią obrotu cylindra co świadczy,

14

Liczba stopni swobody

• Rozważmy przypadek gdy mamy p niezależnych informacji (p

pomiarów) nie obarczonych błędami (gdy dopuścimy błędy pomiarowe oznaczać to może duże błędy zmniejszają liczbę niezależnych

informacji).

• Rozważmy przypadek, gdy mamy dwu elementowy wektor stanu (x₁,x₂) oraz dwa pomiary (y₁, y₂) i prosty model do przodu:



 







 



 







 



 



 





2 1 2

1 2

1

x x 01

. 1 99 . 0

99 . 0 01 . 1 y

y

gdzie błędy są niezależne o wariancji ². Jest to równoznaczne z pomiarem ortogonalnej kombinacji z₁ oraz z₂.

2 1

2 1

2 1

2

2 1

2 1

2 1

1

) x x

( 02 . 0 )

y y

( z

) x x

( 2 )

y y

( z































 _{Zmienna z}

2 ma znacznie mniejszą wartość niż z₁ a więc nie zawiera użytecznej informacji na temat różnicy x₂ – x₁.

15

• Ponieważ, macierze kowariancji mogą posiadać niezerowe elementy poza diagonalą (będące odzwierciedleniem korelacji pomiędzy poszczególnymi elementami) transformujemy

macierz do nowej bazy w której wszystkie wartości pozadiagonalne są zerowe.

) x x

( S

x~  _a^1/2  _a y~  S^_^1/2y











 _^{ }_ K~x~ ~ S

x~

KS S

y~ ^{1/2 1}_a^{/2 1/2}

gdzie: K~ S_^1/2KS¹_a^/2

 

Liczba niezależnych obserwacji jest równa liczbie wartości osobliwych macierzy:

2 / 1

a 2

/

1 KS

S K~ _

 

16

Jest to równoznaczne z liczbą wartości własnych macierzy

większych od jedności. K~ T

K~

17

Analiza błędów





 f (x,b) y

) c , x , bˆ , y ( R

xˆ  _a

Zapiszmy wektor obserwacji w postaci:

gdzie b oznacza wektor parametrów nie wchodzących w

skład wektora stanu (np. natężenie linii widmowej, zależność poszerzenia linii widomych od temperatury itd.), zaś f jest

„forward function” opisującą fizykę pomiaru uwzględniającą np. transfer promieniowania, czy pełny opis aparatury

pomiarowej.

Wektor odzyskiwanych parametrów może być umownie zapisany w postaci:

gdzie R oznacza umownie metodę odwrotną, oznacza

najlepsze oszacowanie parametrów funkcji do przodu f, zaś c jest wektorem parametrów nie występujących podobnie jak wektor informacji a priori x_a w funkcji f, które jednak mogą wpływać na wartości odzyskiwanych parametrów np. przez

bˆ

18

Podstawiając otrzymujemy:

) b , b , x ( f )

b , x (

F  ^'

) c , x , bˆ , )

b , x ( f ( R

xˆ    _a

Dokonujemy linearyzacji modelu do przodu F (y=F(x)+)

gdzie wektor b został podzielony na b i b’ zaś b’ opisuje te parametry funkcji do przodu f, które zostały zignorowane przy konstrukcji modelu do przodu F.

Wyznaczany wektor możemy przepisać do postaci:

) c , x , bˆ , )

' b , b , x ( f )

b , x ( F ( R

xˆ      _a

gdzie f jest błędem modelu do przodu związanym z niepoprawnym opisem fizycznym

) b , x ( F )

' b , b , x ( f

f  



Dokonujemy linearyzacji modelu F w otoczeniu

otrzymujemy _{b bˆ}

x x _a





) c , x , bˆ , )

' b , b , x ( f )

bˆ b

( K )

x x

( K )

bˆ , x ( F ( R

xˆ  _a  _x  _a  _b      _a

gdzie

b K F

x K F

b x



 



 

Obecnie linearyzujemy operator R względem wektora y:

] )

' b , b , x ( f )

bˆ b

( K )

x x

( K [ G ]

c , x , bˆ ), bˆ , x ( F [ R

xˆ  _{a a}  _{x x}  _a  _b     

x G_x R



 

20 y

y

a

a a

a a

G

) x x

( A

x ] c , x , bˆ ), bˆ , x ( F [ R x

xˆ















Ostatecznie różnica pomiędzy wektorem odzyskany a wektorem informacji a priori wynosi:

x K xˆ

G

A _{y x}



 

 _y  K_b(b bˆ)  f (x,b,b')  

bias

wygładzanie

błąd metody odwrotnej

gdzie





















y y

b y

a

G

) ' b , b , x ( f G

) bˆ b ( K G

) x x )(

I A ( x xˆ

Ostatecznie różnica pomiędzy wektorem odzyskany a wektorem informacji a priori wynosi:

błąd wygładzania

błąd parametrów modelu błąd modelu do przodu szum metody odwrotnej