Assimilation de données

On considère dans un premier temps l'opérateur ${\displaystyle H}$ linéaire quitte à le linéariser s'il ne l'est pas. On cherche à minimiser l'erreur commise a posteriori de ${\displaystyle e_{a}}$ en minimisant ${\displaystyle Tr(P_{i,j}^{a})}$.

On cherche la solution à l'aide d'une régression linéaire (voir Méthode des moindres carrés), appelée un Ansatz[1] en assimilation de données :

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

On suppose que les erreurs d'observation et de l'ébauche sont sans biais, quitte à retrancher le biais s'il existe. Si nous voulons que l'erreur de l'analyse reste sans biais, on a ${\displaystyle L=I-KH}$. On obtient alors :

${\displaystyle x_{a}=x_{b}+K(y-Hx_{b})}$

où ${\displaystyle y-Hx_{b}}$ s'appelle le vecteur innovation

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

On cherche maintenant le gain optimal ${\displaystyle K^{*}}$ pour minimiser ${\displaystyle Tr(P_{i,j}^{a})}$. L'analyse BLUE Best linear unbiased estimator permet d'obtenir le gain optimal

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

Supposons maintenant que l'état du système évolue dans le temps. On souhaite effectuer une succession d'analyses à tous les instants possibles. Nous avons des prévisions provenant du modèle aux dates ${\displaystyle t_{0},t_{1},\cdots ,t_{k},\cdots ,t_{n}}$ et des observations à plusieurs dates dont ${\displaystyle t_{k}}$. On note le vecteur de prévision ${\displaystyle x_{k}^{f}}$ (correspondant à ${\displaystyle x_{b}}$ dans le paragraphe précédent), le vecteur des observations ${\displaystyle y_{k}}$ et le vecteur d'analyse ${\displaystyle x_{k}^{a}}$.

On peut d'abord résoudre ce problème à l'aide de méthode dites séquentielles. Dans ce type de méthode, il y a d'abord l'étape de prévision où l'on obtient ${\displaystyle x_{k}^{f}}$, puis l'étape d'analyse où l'on combine l'information des observations et de la prévision pour avoir ${\displaystyle x_{k}^{a}}$. On peut résumer ce problème sous le jeu d'équation suivant :

$${\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}}}$$

Ici ${\displaystyle v_{k+1}}$est l'erreur modèle du passage au temps ${\displaystyle t_{k}}$à ${\displaystyle t_{k+1}}$du au modèle. ${\displaystyle e_{k+1}^{f}}$est l'erreur de prévision accumulée lors de la succession des étapes. On associe à ${\displaystyle e_{k+1}^{f}}$la matrice de covariance ${\displaystyle Q_{k}}$.

On suppose pour ce filtre que les opérateurs ${\displaystyle H_{k}}$ et ${\displaystyle M_{k}}$ sont linéaires et que les erreurs d'observation et de prévision sont sans biais. On peut démontrer que les erreurs d'analyse sont alors sans biais.

Voici l'algorithme du filtre de Kalman dans le cadre de l'assimilation de données.

1. Initialisation  
         Estimer ${\displaystyle x_{0}^{f}}$
         Estimer la matrice de covariance ${\displaystyle P_{0}^{f}}$

2. Boucle sur les différentes dates d'observation ${\displaystyle t_{k}}$
         a. Analyse

                 Calcul du gain avec la méthode BLUE
                 ${\displaystyle K_{k}=P_{k}^{f}H_{k}^{T}(H_{k}P_{k}^{f}H_{k}^{T}+R_{k})^{-1}}$

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

                 Calcul de la matrice de covariance ${\displaystyle P_{k}^{a}}$
                 ${\displaystyle P_{k}^{a}=(I-K_{k}H_{k})P_{k}^{f}}$

         b. Prévision

                 Calculer la nouvelle prévision ${\displaystyle x_{k+1}^{f}}$
                 ${\displaystyle x_{k+1}^{f}=M_{k+1}x_{k}^{a}}$

                 Calculer la matrice de covariance ${\displaystyle P_{k}^{f}}$
                 ${\displaystyle P_{k+1}^{f}=M_{k+1}P_{k}^{a}M_{k+1}^{T}+Q_{k}}$

Le filtre de Kalman étendu reprend exactement le même principe que le filtre de Kalman. Il est juste nécessaire de linéariser les opérateurs ${\displaystyle H_{k}}$et ${\displaystyle M_{k}}$ autour de l'état ${\displaystyle x_{k}}$. On applique ensuite exactement le même algorithme que précédemment. Ce filtre fonctionne bien si l'échantillonnage des observations est assez élevé ou si les non linéarités du modèle ne sont pas trop grandes.

Dans ce cas, nous ne cherchons pas les matrices ${\displaystyle P_{k}^{f}}$et ${\displaystyle P_{k}^{a}}$ mais la densité de probabilité de ${\displaystyle x_{k}^{a}}$. Il faut d'abord poser ce problème sous cette forme appelée filtre bayésien.

On notera ${\displaystyle Y_{k}}$, l'ensemble des observations ${\displaystyle y_{0},\cdots ,y_{k}}$ passées entre les instants ${\displaystyle t_{0}}$ et ${\displaystyle t_{k}}$. On considère maintenant que l'opérateur d'observation ${\displaystyle H_{k}}$ n'est pas nécessairement linéaire et dépend aussi de l'erreur ${\displaystyle y_{k}=H_{k}(x_{k}^{f},e_{k}^{0})}$. A priori, nous connaissons ${\displaystyle p_{Y_{k}|X_{k}^{f}}(y_{k}|x_{k}^{f})}$et ${\displaystyle p_{X_{k}^{f}}(x_{k}^{f})}$. En réalité, ${\displaystyle p_{Y_{k}|X_{k}^{f}}(y_{k}|x_{k}^{f})}$ correspond à ${\displaystyle p_{E^{0}}(e_{k}^{0})}$.

L'idée du filtre particulaire est de calculer les distributions de probabilité à l'aide d'un échantillonnage de l'espace de l'état du système. On crée des particules à partir des points choisis pour l’échantillonnage et leur état va évoluer à l'aide du modèle.

Voici l'algorithme du filtre particulaire bootstrap.

1. Initialisation  
         Échantillonner ${\displaystyle x_{k}^{f}}$à l'aide de ${\displaystyle N}$particules ${\displaystyle {x_{k}^{1},\cdots ,x_{k}^{N}}}$

         Assigner un poids identique ${\displaystyle w_{k}^{i}={\frac {1}{M}}}$aux différentes particules ${\displaystyle x_{k}^{i}}$

2. Boucle sur les différentes dates d'observation ${\displaystyle t_{k}}$
         a. Prévision
                 Propager les particules à l'aide du modèle
                 ${\displaystyle x_{k+1}^{i}=M_{k+1}x_{k}^{i}}$


         b. Analyse

                 Calculer les nouveaux poids des particules
                 ${\displaystyle w_{k+1,\ a}^{i}=w_{k+1}^{i}p(y_{k+1}|x_{k+1}^{i})}$

                 Normaliser les poids pour obtenir la distribution de ${\displaystyle x_{k}^{a}}$
 
         c. Re-échantillonnage

                 Le filtre va privilégier une particule si on ne le ré-échantillonne pas (phénomène appelé dégénérescence). 
                 On ré-échantillonne ${\displaystyle x_{k}^{f}}$avec des poids identiques.    

En général cette méthode est efficace pour des modèles fortement non linéaires mais si la dimension de l'état du système est trop grande alors le filtre ne fonctionne plus (en général plus grand que 8). On peut aussi trouver des variantes où l'on ré-échantillonne seulement les particules qui ont un poids trop élevé.

Le filtre d'ensemble utilise lui aussi la notion de particule mais il ne générera que les moments d'ordre 1 et 2 de l'analyse. L'analyse est la même que le filtre de Kalman mais des particules sont créées pour propager les erreurs dues à l'observation.
Ce filtre fonctionne avec un modèle non linéaire mais il faut linéariser la fonction d'observation pour calculer le gain.

Voici l'algorithme:

1. Initialisation  
         Estimer ${\displaystyle x_{0}^{f}}$
         Estimer la matrice de covariance ${\displaystyle P_{0}^{f}}$
         Créer N particules estimant ${\displaystyle x_{0}^{f}}$à l'aide la matrice de covariance ${\displaystyle P_{0}^{f}}$

2. Boucle sur les différentes dates d'observation ${\displaystyle t_{k}}$

         a. Observation

                 Créer un jeu d'observation ${\displaystyle {y_{k}^{1},\cdots ,y_{k}^{N}}}$ de biais nulle autour de la valeur observée ${\displaystyle y_{k}}$

                 Calculer la matrice de covariance ${\displaystyle R_{k}}$associée
         b. Analyse

                 Calcul du gain avec la méthode BLUE
                 ${\displaystyle K_{k}=P_{k}^{f}H_{k}^{T}(H_{k}P_{k}^{f}H_{k}^{T}+R_{k})^{-1}}$
                 Ici ${\displaystyle H_{k}}$linéarisé

                 Estimation 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}))}$
                 Ici ${\displaystyle H_{k}}$non linéarisé

                 Calculer la moyenne ${\displaystyle x_{k,i}^{a}}$

                 Calcul de la matrice de covariance ${\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. Prévision

                 Calculer les nouvelles prévisions ${\displaystyle x_{k+1,i}^{f}}$
                 ${\displaystyle x_{k+1,i}^{f}=M_{k+1}x_{k,i}^{a}}$

                 Calculer la matrice de covariance ${\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}}$

Il est possible d'associer des filtres pour réduire la dimensionnalité du système. Il existe plusieurs filtres comme le filtre RRSQRT[2], SEEK[3] ou encore SEIK[4].

La méthode d'assimilation variationnelle est utilisée pour obtenir les valeurs aux points de grille du modèle les plus près de la réalité. Elle implique de trouver un ensemble de points du modèle dont la description par une fonction se rapproche le plus des valeurs aux points observés sans introduire d'instabilité dans le modèle numérique. Elle consiste donc à chercher l'état le plus vraisemblable à partir des connaissances disponibles sur les lois de probabilités des erreurs d'observation.

Ceci se fait en minimisant par itération la fonction coût, le plus souvent la somme des moindres carrés des déviations entre l'analyse et l'observation pondérée par la qualité de ces dernières. Ce processus peut être fait en 3 ou 4 dimensions.

La méthode à trois dimensions, communément appelée 3D-Var, se fait à un pas de temps fixe dans les trois dimensions cartésiennes X, Y et Z. Comme pour le filtre de Kalman, le 3D-Var consiste à minimiser la distance au sens des moindres carrés entre l'état estimé et les différentes sources d'informations telles que la prévision précédente et les observations au temps initial. Le nouvel état analysé est, en général, utilisé comme point de départ de la prévision suivante.

La fonction coût s'exprime comme[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} ]),}$

Où :

• ${\displaystyle \mathbf {B} }$ est la matrice de covariance de l'erreur de bruit de fond ;
• ${\displaystyle \mathbf {R} }$ est la matrice de covariance de l'erreur d'observation.

À quatre dimensions, l'analyse se fait à plusieurs pas temps entre le temps initial et un temps futur de prévision. Il s'agit donc d'une extension de la méthode 3D-Var qui ne vise pas à obtenir l'état optimal à un instant donné, mais la trajectoire optimale sur une fenêtre de temps donnée. Les observations sont donc prises en compte aussi bien dans leur distribution spatiale que temporelle et le 4D-Var propage donc l'information apportée par les observations à l'instant initial de la fenêtre d'assimilation[6].

Cette amélioration du 3D-Var permet d'ajouter la connaissance de l'évolution du système comme information pour l'analyse. Bien qu'elle demande une beaucoup plus grande puissance de calcul que la méthode précédente, elle est devenue la plus utilisée dans les systèmes de prévision opérationnels atmosphériques du CEPMMT en 1997, de Météo-France en 2000, et de nombreux autres centres météorologiques internationaux[6].

Les techniques variationnelles sont plus efficaces pour trouver une bonne analyse et les techniques séquentielles permettent une caractérisation des erreurs. Ainsi, de nouvelles méthodes sont inventées pour combiner ces deux aspects.