Assimilação de dados

A previsão numérica da evolução da atmosfera depende grandemente das condições iniciais que a ele estão proporcionadas. Um bem que é difícil de determinar, a um dado instante, o estado do atmosfera, isto é o conjunto das variáveis atmosféricas (pressão, temperatura, humidade, etc.) no conjunto do volume, com uma boa resolução e uma boa precisão.

As únicas informações disponíveis a um dado momento são as observações meteorológicas de diferente natureza (rádio sondas, Estações meteorológicas, bóia oceânica, etc.). Mas estas informações não são suficientes. Efectivamente o modelo atmosférico requer da ordem de ${\displaystyle 10^{7}}$  valores (para todos os campos físicos considerados, em todos os pontos do modelo). O bem atual das observações são da ordem de ${\displaystyle 10^{6}}$ . Uma mera interpolação não basta nestas condições. Ter então recursos a um método chamado "assimilação de dados".

A assimilação de dados é um método "predição/correcção". Uma previsão, calculada ao passo de tempo precedente e válido ao instante considerado, é utilizada como predição. As observações disponíveis permitem corrigir este esboço para estimar ao melhor o estado real da atmosfera.

Deseja-se conhecer a temperatura numa peça que dispõe de uma fonte de calor pontual de 20°C.

 

A fonte está presa. Ao tempo ${\displaystyle t_{0}}$  , faz 15°C em toda a peça. A fonte de calor activa-se, e é observada na peça.

A previsão consiste em dizer o que fará, ao cabo de um verdadeiro tempo ${\displaystyle dt}$  20°C no ponto de aplicativo da fonte, depois a cada vez mais frio em se afastando : trata-se aqui da previsão válida no conjunto da peça.

O observador volta 3 horas depois. Um termômetro fixado na peça indica 17°C num ponto bastante afastado da fonte onde está suposto fazer 18°C. A assimilação parte da ideia que esta informação vai corrigir a previsão precedente. Por exemplo supondo que localmente, uma reação do ar faz baixar esta temperatura. Ou ainda que a descida da temperatura para além da fonte de calor faz-se mais rapidamente. Obtemos assim uma análise da situação.

O termômetro não é muito preciso, por exemplo um erro de +/- 0,5°C. O conhecimento do erro de observação devida à carência de precisão do termômetro reduzirá o impacto desta observação durante a nossa assimilação. O conhecimento do erro de previsão (por exemplo a carência de informação na isolação exacta da peça), vai jogar no outro sentido. Estes diferentes aspectos estarão expostos mais longe após a formulação matemática.

Querer conhecer o estado de um sistema que não evolui no tempo representado por um vetor ${\displaystyle x}$ .(com frequência de dimensão infinita). Se supor-mos ${\displaystyle x}$  espacialmente com ajuda de um operador ${\displaystyle \Pi }$  que da o vetor ${\displaystyle x_{t}}$ de dimensão ${\displaystyle n}$ .

Fazer medidas a certos lugares no sistema. Reunir estas informações no vetor das observações ${\displaystyle y}$ de dimensão p. Associa-se um vetor de erro de observação ${\displaystyle e_{0}}$ cuja se não conhece, em geral, que a esperança e a variante (ou melhor a sua estimativa). Poder pois construir uma matriz de covariantes ${\displaystyle R_{i,j}=E[e_{0}^{T}.e_{0}]}$ . Comparar as observações ao estado real do sistema com ajuda da função de observação ${\displaystyle H}$ (que pode ser linear ou não):

${\displaystyle y=Hx_{t}+e_{0}}$ 

Utilizar um modelo que nós permite avaliar o estado do sistema ${\displaystyle x_{t}}$ de maneira aproximada. A estimativa de ${\displaystyle x_{t}}$ pelo modelo está chamado vetor de esboço anotado ${\displaystyle x_{b}}$ . Associar este vetor com um erro ${\displaystyle e_{b}}$ e uma matriz de covariante ${\displaystyle P_{i,j}^{b}=E[e_{b}^{T}.e_{b}]}$ .

A assimilação de dados consiste em dar a melhor aproximação do estado do sistema ${\displaystyle x_{t}}$ a retirar do esboço e das observações. O vetor resultante é o vetor de análise ${\displaystyle x_{a}}$ . Procurar também a matriz de covariante de erro ${\displaystyle P_{i,j}^{a}}$ .

Considera-se num primeiro momento o operador ${\displaystyle H}$  linear sem prejuízo do linearizar se não o é . Procurar minimizar o erro cometido a posteriori de ${\displaystyle e_{a}}$ minimizando ${\displaystyle Tr(P_{i,j}^{a})}$ .

Procurar a solução com ajuda de uma regreção linear (ver método do mínimo quadrado), chamado um Ansatz.^[1]em assimilação de dado:

${\displaystyle x_{a}=Lx_{b}+Ky}$ 

Supor que os erros de observação e do esboço estão sem distorção sem prejuízo de retorcer a distorção se existe. Se queremos que o erro da análise se apresente sem distorção, se tem ${\displaystyle L=I-KH}$ . Obtemos então :

${\displaystyle x_{a}=x_{b}+K(y-Hx_{b})}$ onde ${\displaystyle y-Hx_{b}}$ se chama o vetor inovação

${\displaystyle P^{a}=(I-KH)P^{b}(I-KH)^{T}+KRK^{T}}$ 

Procurar agora o ganho óptimo ${\displaystyle K^{*}}$  para minimizar ${\displaystyle Tr(P_{i,j}^{a})}$ . A análise BLUE Best linear unbiased estimator permite obter o ganho óptimo.

${\displaystyle K^{*}=P^{b}H^{T}(R+HP^{b}H^{T})^{-1}}$ 

Supomos agora que o estado do sistema evolui no tempo. Desejar efectuar uma sucessão de análise a todos os instantes possíveis. Temos previsões procedentes do modelo às datas${\displaystyle t_{0},t_{1},\cdots ,t_{k},\cdots ,t_{n}}$  e das observações a várias datas cujas ${\displaystyle t_{k}}$ . Anotar o vetor de previsão ${\displaystyle x_{k}^{f}}$ (que corresponde a ${\displaystyle x_{b}}$  no parágrafo que precede), o vetor das observações ${\displaystyle y_{k}}$ e o vecteur de análise ${\displaystyle x_{k}^{a}}$ .

Pode-se primeiramente resolver este problema com ajuda de método dito sequenciais. Neste tipo de método, há primeiramente a etapa de previsão onde se obtém ${\displaystyle x_{k}^{f}}$ , depois a etapa de análise onde se combina a informação das observações e da previsão para ter ${\displaystyle x_{k}^{a}}$ . Poder resumir este problema baixo o jogo de equação que segue:

$${\displaystyle {\begin{cases}x_{k+1}^{f}&=M_{k+1}(x_{k})+v_{k+1}\\y_{k}&=H_{k}(x_{k})+e_{k}^{0}\end{cases}}}$$ 

Aqui ${\displaystyle v_{k+1}}$ é o erro do modelo da passagem do tempo ${\displaystyle t_{k}}$ a ${\displaystyle t_{k+1}}$  ao modelo. ${\displaystyle e_{k+1}^{f}}$ é o erro de previsão acumulada durante a sucessão das etapas. Associar a ${\displaystyle e_{k+1}^{f}}$ a matriz de covariante ${\displaystyle Q_{k}}$ .

Supõe-se para este filtro que o operador ${\displaystyle H_{k}}$  e ${\displaystyle M_{k}}$  são lineares e que os erros de observação e de previsão são sem distorção. Poder demonstrar que os erros de análises estão mostrados sem distorções.

Aqui está o algoritmo do filtro de Kalman no marco da assimilação de dados.

1. Initialização 
         Estimar ${\displaystyle x_{0}^{f}}$ 
         Estimar a matris de covariante ${\displaystyle P_{0}^{f}}$ 

2. Volta sobre as diferentes datas da observação ${\displaystyle t_{k}}$ 
         a. Análise

                 Calcular o ganho com o método BLUE
                 ${\displaystyle K_{k}=P_{k}^{f}H_{k}^{T}(H_{k}P_{k}^{f}H_{k}^{T}+R_{k})^{-1}}$ 

                 Estimativa de ${\displaystyle x_{0}^{a}}$ 
                 ${\displaystyle x_{k}^{a}=x_{k}^{f}+K_{k}(y_{k}-H_{k}x_{k}^{f})}$ 

                 Calculo da matriz da covariante ${\displaystyle P_{k}^{a}}$ 
                 ${\displaystyle P_{k}^{a}=(I-K_{k}H_{k})P_{k}^{f}}$ 

         b. Previsão

                 Calcular a nova previsão ${\displaystyle x_{k+1}^{f}}$ 
                 ${\displaystyle x_{k+1}^{f}=M_{k+1}x_{k}^{a}}$ 

                 Calcular a matriz de covariante ${\displaystyle P_{k}^{f}}$ 
                 ${\displaystyle P_{k+1}^{f}=M_{k+1}P_{k}^{a}M_{k+1}^{T}+Q_{k}}$ 

O filtro de Kalman extenso retoma exactamente o mesmo princípio que o filtro de Kalman. É apenas necessário de linearizar os operadores ${\displaystyle H_{k}}$ e ${\displaystyle M_{k}}$  em torno do estado ${\displaystyle x_{k}}$ . Aplicar logo exactamente o mesmo algoritmo que precedente. Este filtro funciona bem se a amostragem das observações é elevada bastante ou se as não linearizar do modelo não são demasiado grande.

Neste caso, não procuramos as matrizes ${\displaystyle P_{k}^{f}}$ e ${\displaystyle P_{k}^{a}}$  mas a densidade de probabilidade de ${\displaystyle x_{k}^{a}}$ . Falha primeiramente propor este problema sob esta forma chamada filtro bayesiana.

Anotar ${\displaystyle Y_{k}}$ , o conjunto das observações ${\displaystyle y_{0},\cdots ,y_{k}}$  passadas entre os instantes ${\displaystyle t_{0}}$  e ${\displaystyle t_{k}}$ . Considerar agora que o operador de observação ${\displaystyle H_{k}}$  não é necessariamente linear e depende também do erro ${\displaystyle y_{k}=H_{k}(x_{k}^{f},e_{k}^{0})}$ . A priori, conhecemos ${\displaystyle p_{Y_{k}|X_{k}^{f}}(y_{k}|x_{k}^{f})}$ e ${\displaystyle p_{X_{k}^{f}}(x_{k}^{f})}$ . Em realidade, ${\displaystyle p_{Y_{k}|X_{k}^{f}}(y_{k}|x_{k}^{f})}$  corresponde a ${\displaystyle p_{E^{0}}(e_{k}^{0})}$ .

A ideia do filtro de partículas é de calcular as distribuições de probabilidade com ajuda de uma amostragem do espaço do estado do sistema. Criar partículas a marchar dos pontos escolhidos para a amostragem e seu estado vai evoluir à ajuda do modelo.

Aqui está o algoritmo do filtro de partículas bootstrap.

1. Initialização 
         Amostragem ${\displaystyle x_{k}^{f}}$ à ajuda de ${\displaystyle N}$ partículas ${\displaystyle {x_{k}^{1},\cdots ,x_{k}^{N}}}$ 

         Atribuir um peso idêntico ${\displaystyle w_{k}^{i}={\frac {1}{M}}}$ às diferentes partículas ${\displaystyle x_{k}^{i}}$ 

2. Voltar sobre as diferentes datas de observação ${\displaystyle t_{k}}$ 
         a. Previsão
                 Propagar as partículas à ajuda do model
                 ${\displaystyle x_{k+1}^{i}=M_{k+1}x_{k}^{i}}$ 


         b. Analisar

                 Calcular os novos pesos das partículas
                 ${\displaystyle w_{k+1,\ a}^{i}=w_{k+1}^{i}p(y_{k+1}|x_{k+1}^{i})}$ 

                 Normalizar os pesos para obter a distribuição de ${\displaystyle x_{k}^{a}}$ 

         c. Re-amostragem
                 O filtro vai privilegiar uma partícula se se não o amostrado  (fenómeno chamado degenarativo).
                 Se re-amostrar ${\displaystyle xf_{k}}$ com do peso idêntico.   

Em geral este método é eficaz para modelos fortemente não lineares mas se a dimensão do estado do sistema é demasiado grande então o filtro já não funciona (em geral maior que 8). Poder também encontrar variantes onde se re-amostra só as partículas que têm um peso demasiado elevado.

O filtro de grupos utiliza também a noção de partícula mas só gerará os momentos de ordem 1 e 2 da análise. A análise é a mesma que o filtro de Kalman mas das partículas estão criar para propagar os erros da observação.
Este filtro funciona com um modelo não linear mas faz falta linearizar a função de observação para calcular o ganho.

Aqui está o algoritmo:

1. Initialização 
         Estimar ${\displaystyle x_{0}^{f}}$ 
         Estimar a matriz de covariante ${\displaystyle P_{0}^{f}}$ 
         Criar N partículas que estimam ${\displaystyle x_{0}^{f}}$ à ajuda a matriz de covariante ${\displaystyle P_{0}^{f}}$ 

2. Volta nas diferentes datas de observação ${\displaystyle t_{k}}$ 

         a. Observação

                 Criar um jogo de observação ${\displaystyle {y_{k}^{1},\cdots ,y_{k}^{N}}}$  de distorção nula em torno do valor observado ${\displaystyle y_{k}}$ 

                 Calcular a matriz de covariante ${\displaystyle R_{k}}$ associada
         b. Análise

                 Cálculo do ganho com o método BLUE
                 ${\displaystyle K_{k}=P_{k}^{f}H_{k}^{T}(H_{k}P_{k}^{f}H_{k}^{T}+R_{k})^{-1}}$ 
                 Aqui ${\displaystyle H_{k}}$ lineariza

                 Estimativa de ${\displaystyle x_{k,i}^{a}}$ 
                 ${\displaystyle x_{k,i}^{a}=x_{k,i}^{f}+K_{k}(y_{k}^{i}-H_{k}(x_{k,i}^{f}))}$ 
                 Aqui ${\displaystyle H_{k}}$ não lineariza

                 Calcular a média ${\displaystyle x_{k,i}^{a}}$ 

                 Cálculo da matriz de covariante ${\displaystyle P_{k}^{a}}$ 
                 ${\displaystyle P_{k}^{a}={\frac {1}{N-1}}\sum _{j=1}^{N}(x_{k,j}^{a}-{\bar {x_{k}^{a}}})(x_{k,j}^{a}-{\bar {x_{k}^{a}}})^{T}}$ 

         c. Previsão

                 Calcular as novas previsões ${\displaystyle x_{k+1,i}^{f}}$ 
                 ${\displaystyle x_{k+1,i}^{f}=M_{k+1}x_{k,i}^{a}}$ 

                 Calcular a matriz de covariante ${\displaystyle P_{k}^{f}}$ 
                 ${\displaystyle P_{k}^{f}={\frac {1}{N-1}}\sum _{j=1}^{N}(x_{k,j}^{f}-{\bar {x_{k}^{f}}})(x_{k,j}^{f}-{\bar {x_{k}^{f}}})^{T}}$ 

É possível de associar filtros para reduzir a dimensionalidade do sistema. Há vários filtros como o filtro RRSQRT,^[2] SEEK^[3] ou ainda SEIK.^[4]

 Ver artigo principal: Cálculo variacional

O método de assimilação variacional é utilizada para obter os valores aos pontos de grade do modelo os mais para perto da realidade. Implica encontrar um conjunto de pontos do modelo cuja descrição por uma função se acerca o mais perto dos valores aos pontos observados sem introduzir instabilidade no modelo numérico. Consiste pois de procurar o estado o mais verdadeiro a par dos conhecimentos disponíveis nas leis de probabilidades dos erros de observação.

Isto se faz minimizando por iteração a função perda, com mais frequência a soma dos mínimos quadradas dos desvios entre a análise e a observação ponderada pela qualidade destas últimas. Este processo pode ser feito em 3 ou 4 dimensões.

O método a três dimensões, comunemente chamada 3D-Var, faz-se a um passo de tempo fixo nas três dimensões cartesianas X, Y e Z. Como para o filtro de Kalman, o 3D-Var consiste em minimizar a distância ao sentido dos mínimos quadrados entre o estado estimado e as diferentes fontes de informações tais que a previsão precedente e as observações ao tempo inicial. O novo estado analisado é, em geral, utilizado como ponto de saída da previsão seguinte.

A função perda expressa-se como^[5]: ${\displaystyle J(\mathbf {x} )=(\mathbf {x} -\mathbf {x} _{b})^{\mathrm {T} }\mathbf {B} ^{-1}(\mathbf {x} -\mathbf {x} _{b})+(\mathbf {y} -{\mathit {H}}[\mathbf {x} ])^{\mathrm {T} }\mathbf {R} ^{-1}(\mathbf {y} -{\mathit {H}}[\mathbf {x} ]),}$ 

Onde :

• ${\displaystyle \mathbf {B} }$  é a matriz de covariante do erro de ruído de fundo ;
• ${\displaystyle \mathbf {R} }$  é a matriz de covariante do erro de observação.

A quatro dimensões, a análise faz-se a vários passos de tempos entre o tempo inicial e um tempo futuro da previsão. Trata-se pois de uma extensão do método 3D-Var que não aponta por obter o estado óptimo a um dado instante, mas a trajectória óptima numa janela de dado tempo. As observações são pois tomadas em conta também em sua distribuição espacial e temporária e o 4D-Var propaga pois a informação contribuída pelas observações ao instante inicial da janela de assimilação.^[6]

Esta melhoria do 3D-Var permite acrescentar o conhecimento da evolução do sistema como informação para a análise. Ainda que pede um bem maior de grande potência de cálculo que o método precedente, está feito para ser a mais utilizada nos sistemas de previsão operativa atmosféricas do CEPMMT em 1997, da Météo-France em 2000, e de numeroso outros centros meteorológicos internacionais.^[6]

As técnicas variacionais são mais eficazes para encontrar uma boa análise e as técnicas sequenciais permitem uma caracterização dos erros. Assim, os novos métodos estão inventadas para combinar estes dois aspectos.

• Marc Bocquet, [1], notas de cursos do Groupe des écoles nationales supérieures de techniques avancées e a École nationale des ponts et chaussées.

•   Portal da meteorologia