jOas wEb Sci

Latex est revenu !

jOas — Sun, 11 Apr 2010 12:56:00 +0200

J'ai trouvé un nouveau site qui heberge le script permettant de faire apparaître les formules latex !

Il suffit d'insérer :


<script type="text/javascript" src="http://mathcache.s3.amazonaws.com/replacemath.js"></script>
<script type="text/javascript">
replaceMath( document.body );
</script>

Je peux ainsi remettre en ligne tous mes billets !!!

Revue des méthodes d'assimilation : tout en un

jOas — Fri, 08 Feb 2008 22:11:00 +0000

La catastrophe !
Le script que j'utilisais pour faire apparaître mes formules mathématiques n'est plus hébergé.
Mes billets ne ressemblent plus à rien... Je les ai tous mis hors ligne.
Mais...
Mais tout ça est disponible en format postscript sur le site du CERFACS : revue_methodes_assimilation.ps.gz
Ou en format pdf en pièce jointe.

L'assimilation de données

jOas — Thu, 22 Nov 2007 08:23:00 +0100

Chose promise, chose due. Alors, voici une vision générale de ce qui occupe les plus longues de mes heures : l'assimilation de données.

Ce billet regroupe les différents billets publiés sur le contexte général de l'assimilation de données sous forme d'une table des matières hiérarchisée. Peut-être verrez-vous d'ici quelques temps fleurir des billets sur le cadre plus particulier de mon travail : l'assimilation de données océanographique avec une méthode variationnelle. Qui sait ?

Bonne lecture...

I - Préambule

Introduction à l'assimilation de données

II - Introduction

III - Présentation du problème

IV - Interpolation statistique

V - Méthodes d'assimilation

L'interpolation optimale (OI)
Les filtres de Kalman de rang plein
Les filtres de Kalman de rang réduit
Le filtre de Kalman d'ensemble (EnKF)
Méthodes variationnelles
1. 3D-Var
2. 4D-Var

VI - Modèle de covariance d'erreur

Modèle de covariance d'erreur - Modélisation des erreurs

jOas — Fri, 16 Nov 2007 01:24:00 +0100

Comme il a été expliqué précédemment, la matrice de covariances d'erreur d'ébauche \[\mathbf{B}\] complète est trop grande pour être spécifiée explicitement. En général, les variances qui représentent les \[n\] termes de la diagonale de \[\mathbf{B}\] sont tous définis. Les termes non-diagonaux sont plus difficiles à définir. En effet, la matrice \[\mathbf{B}\] doit être définie positive. Les modélisations des termes non-diagonaux doivent donc conserver cette qualité. De plus, ces modélisations doivent, en général, imposer des propriétés physiques qui seront réfléchies dans l'analyse :

Les corrélations doivent être lissées à l'échelle des processus physiques ;
Les corrélations doivent tendre vers zéro pour des grandes distances de séparation car les observations ne doivent avoir qu'un impact local ;
Les corrélation ne doivent pas avoir des variations en fonction des directions ou de la position qui ne puissent être expliquées physiquement ;
Les principales propriétés d'équilibre du système doivent être renforcées ;
Les corrélations ne doivent pas amener des variances d'erreur d'ébauche sur les paramètres observés qui ne soient pas raisonnable.

Ces différentes exigences conduisent à une spécification des covariances d'erreur d'ébauche très complexe et qui peut être comparée à un réglage de paramètres physiques. Les hypothèses reposant sur la physique doivent effectuées et testées avec beaucoup de précaution. La liste ci-dessous cite un certains nombre de techniques populaires :

Les modèles de corrélation peuvent être définis indépendamment des champs de variances à la condition que les échelles de variation des variances soit plus grande que celles des corrélations (sinon la forme des covariances sera très différente de celle des corrélations et les conséquences sur les propriétés d'équilibre seront imprévisibles) ;
Les matrices d'autocorrélations verticales de chaque paramètre sont suffisamment petites pour être spécifiées explicitement ;
Les autocorrélations horizontales ne peuvent pas être spécifiées explicitement, mais elles peuvent être réduites à des matrices creuses sous l'hypothèse d'homogénéité et d'isotropie. Ces matrices creuses sont alors diagonales dans l'espace spectral et comparables à des filtres digitales passe-bas dans l'espace physique ;
Des modèles de corrélations tridimensionnelles et multivariées peuvent être construits en utilisant les hypothèses de séparabilité, d'homogénéité et d'indépendance ;
Les contraintes d'équilibres peuvent être utilisées en transformant les variables du modèle en composantes équilibrées et non-équilibrées. La partie non-équilibrée est supposée avoir une variance d'erreur d'ébauche plus faible que la partie équilibrée, ce qui signifie quelle contribue moins à la construction de l'incrément d'analyse ;

Il existe, bien sûr, beaucoup d'autres techniques permettant de construire des opérateurs modélisant les covariances d'erreur d'ébauche.

Modèle de covariance d'erreur - Estimation des erreurs

jOas — Fri, 16 Nov 2007 01:15:00 +0100

Il est difficile d'estimer les erreurs car elles ne sont jamais observées directement. En effet, l'état vrai n'étant pas accessible, il est impossible d'obtenir des échantillons des erreurs d'ébauche et d'observation. Les données statistiques sont donc difficilement disponibles et très largement insuffisantes pour déterminer tous les éléments. Par ailleurs, les matrices de covariances d'erreur sont très grandes. Pour ces deux raisons, elles doivent être simplifiées et modélisées. De tailles réduites, ces matrices sont manipulables informatiquement parlant et nécessitent moins de d'informations statistiques pour les décrire.

La modélisation des covariances d'erreur est donc un problème difficile et il est nécessaire de faire des hypothèses d'homogénéités. La meilleure source d'information est clairement l'étude de l'innovation (\[\mathbf{d}=\mathbf{y}-H\mathbf{x}^b\]) et peut être utilisée de plusieurs manières différentes. D'autres informations peuvent être obtenues à partir du vecteur d'erreur d'analyse (\[\mathbf{y}-H\mathbf{x}^a\]) ou à partir de la valeur de la fonction coût pour les méthodes variationnelles. D'autres méthodes permettent d'estimer les covariances d'erreur d'ébauche avec des quantités dont les statistiques d'erreur sont équivalentes à celle de l'erreur d'ébauche. Parmi ce type de techniques, la méthode NMC est très connue mais ces bases théoriques sont friables. Une autre possibilité est d'utiliser une méthode d'ensemble de la même manière que pour le filtre de Kalman éponyme. Cette méthode est néanmoins applicable quelque soit la méthode d'assimilation utilisée.

Méthode basée sur l'innovation

Cette méthode a été introduite en météorologie à la fin de années 80 (Hollingworth et Lönnberg, 1986). Elle est basée sur l'utilisation de l'innovation (\[\mathbf{d}=\mathbf{y}^o-H\mathbf{x}^b\]) d'un réseau d'observations suffisamment grand et dense tel qu'il puisse fournir des informations sur les différentes échelles présentes dans la physique du système. Cette méthode permet d'obtenir des statistiques moyennes permettant de construire des matrices de covariances d'erreur statiques. Deux hypothèses importantes sont émises : les erreurs d'ébauches sont indépendantes des erreurs d'observation et les erreurs d'observation ne sont pas corrélées spatialement. Le principe est ensuite assez simple (Fig. innovation)). Il suffit de construire un histogramme représentant les covariances de du vecteur d'innovation en fonction de la distance de séparation. Pour une séparation nulle, l'histogramme fournit une information moyenne sur les variances d'erreurs d'ébauche et d'observation. Pour une distance non-nulle, l'histogramme ne fournit plus qu'une information moyenne sur les corrélations d'erreur d'ébauche.

Fig. innovation : Représentation de la méthode basée sur l'innovation. Les statistiques de covariances de l'innovation (\[\mathbf{y}^o-H\mathbf{x}^b\]) d'un système d'assimilation sont rangées dans un histogramme en fonction de la distance séparant les deux points. L'histogramme à l'origine permet d'estimer les variances d'ébauche et d'observation moyennes.

Soit deux points d'observation \[i\] et \[j\], la covariance d'innovation \[c(i,j)\] s'écrit

\[ c(i,j)=E\left[(\mathbf{y}^o_i-\mathbf{H}_i\mathbf{x}^b)(\mathbf{y}^o_j-\mathbf{H}_j\mathbf{x}^b)^T\right]\]

\[c(i,j)=E\left[\left((\mathbf{y}^o_i-\mathbf{H}_i\mathbf{x}^t)+(\mathbf{H}_i\mathbf{x}^t-\mathbf{H}_i\mathbf{x}^b)\right)\left((\mathbf{y}^o_j-\mathbf{H}_j\mathbf{x}^t)+(\mathbf{H}_j\mathbf{x}^t-\mathbf{H}_j\mathbf{x}^b)\right)^R\right]\]

\[c(i,j)=E\left[(\mathbf{y}^o_i-\mathbf{H}_i\mathbf{x}^t)(\mathbf{y}^o_j-\mathbf{H}_j\mathbf{x}^t)^T\right]+ \mathbf{H}_i E\left[(\mathbf{x}^t-\mathbf{x}^b)(\mathbf{x}^t-\mathbf{x}^b)^T\right]\mathbf{H}_j^T\]
\[ +E\left[(\mathbf{y}^o_i-\mathbf{H}_i\mathbf{x}^t)(\mathbf{x}^t-\mathbf{x}^b)^T\right]\mathbf{H}_j^T+ \mathbf{H}_iE\left[\mathbf{x}^t-\mathbf{x}^b)(\mathbf{y}^o_j-\mathbf{H}_j\mathbf{x}^t)^T\right]\]

(067)

\[c(i,j)=\mathbf{R}_{ij}+\mathbf{H}_i\mathbf{B}\mathbf{H}_j^T\]

En utilisant l'hypothèse que l'erreur d'ébauche n'est pas corrélée à l'erreur d'observation, l'Eq. (067) ne conserve que deux termes : le premier est la covariance d'erreur d'observation entre les points \[i\] et \[j\] ; le second est la matrice de covariances d'erreur d'ébauche interpolée (si l'opérateur d'observation n'agit qu'en tant qu'opérateur d'interpolation) en ces points. Ceci en supposant que ces deux termes sont homogènes sur l'ensemble des observations.

Si les points \[i\] et \[j\] sont identiques (\[i=j\]), alors la corrélation du vecteur d'innovation au point \[i\] est la somme des variances d'erreurs d'ébauche et d'observation (\[c(i,j)=s^o_i+s^b_i\]). Si les points \[i\] et \[j\] sont différents (\[i\neq j\]) et que l'erreur d'observation n'est pas corrélée spatialement, alors la corrélation du vecteur d'innovation entre les points \[i\] et \[j\] est la covariance d'erreur d'ébauche entre ces points (\[c(i,j)=\mathbf{H}_i\mathbf{B}\mathbf{H}_j^T\]). A noter que la décorrélation spatiale d'erreur d'observation est fondamentale, car seule cette hypothèse permet de séparer l'information provenant de la matrice de covariances d'erreur d'observation \[\mathbf{R}\] et d'ébauche \[\mathbf{B}\].

A partir de ces hypothèses, si les points \[i\] et \[j\] sont très proches l'un de l'autre sans jamais être égaux, alors la corrélation du vecteur d'innovation entre les point \[i\] et \[j\] tend vers la variance d'erreur d'ébauche au point \[i\] (\[\lim_{i\to j}c(i,j)=s^b_i\]). En prolongeant la courbe formée par la corrélation du vecteur d'innovation vers une séparation nulle, il est donc possible d'obtenir la variance d'erreur d'ébauche. La variance d'erreur d'observation est alors la différence entre la corrélation du vecteur d'innovation pour une séparation nulle et la variance d'erreur d'ébauche obtenue (\[s^o_i=c(i,j)-s^b_i\]). Il est aussi possible d'obtenir les corrélations d'erreur d'ébauche en fonction de la distance de séparation en prenant le rapport de la corrélation du vecteur d'innovation sur la variance d'erreur d'ébauche (\[c(i,j)/s^b_i\]). Ce résultat n'est possible que si les variances d'erreur d'ébauche sont homogènes sur tout le jeu d'observations.

Si les covariances d'erreur d'ébauche ne tendent pas vers zéro pour unegrande distance de séparation, c'est le signe de la présence d'un biais dans l'ébauche et/ou les observations. Dans ce cas, cette méthode ne fonctionnera pas correctement.

La méthode basée sur l'innovation est la seule méthode directe permettant de diagnostiquer les statistiques d'erreur. Cependant, elle ne fournit des informations que dans l'espace des observations et donc que dans les régions observées. Pour obtenir de bons résultats, il faut un réseau d'observations uniforme et pas trop dense pour ne pas biaiser les statistiques. Cette méthode n'est donc pas toujours très pratique pour spécifier les statistiques des erreurs. De plus, elle ne fournit que des valeurs moyennes ne permettant que de construire des matrices de covariances d'erreur statiques.

Méthode NMC

La méthode NMC (Parrish et Derber, 1992) est ou a été utilisée dans de nombreux centres de prévision météorologique. Elle permet de construire de construire une matrice de covariances d'erreur d'ébauche statique. L'idée est de calculer des différences entre des prévisions valides au même instant mais de durées différentes. A partir de ces différences, il est possible de d'obtenir des statistiques qui peuvent être liées à la matrice de covariances d'erreur d'ébauche \[\mathbf{B}\]. A partir d'un système d'assimilation séquentielle, il est très facile de mettre cette méthode en oeuvre. A la fenêtre d'assimilation démarrant à l'instant \[t_{i-1}\], une prévision d'une durée de deux cycles d'assimilation est effectuée (de \[t_{i-1}\] à \[t_{i+1}\]). A partir du même état à l'instant \[t_{i-1}\], un cycle d'assimilation est effectué et permettent d'obtenir un état à l'instant \[t_i\] à partir duquel une prévision est effectuée de d'une durée d'un cycle (de \[t_i\] à \[t_{i+1}\]). Ces deux prévisions sont donc valides au même instant \[t_{i+1}\]. Il est possible d'envisager ces deux prévisions comme des prévisions d'une durée d'un seul cycle d'assimilation mais dont les conditions initiales à \[t_i\] varient. Ces différences de conditions initiales reflètent l'impact de l'assimilation. En calculant les différences entre ces deux prévisions et en reproduisant ces expériences suffisamment de fois, il est alors possible de calculer des statistiques sur ces différences. Le principe de la méthode NMC est illustré par la Fig. NMC.

Fig. NMC : Méthode NMC. Une prévision démarre à \[t_{i-1}\] et dure jusqu'à \[t_{i+1}\]. Au même instant, une analyse est effectuée entre \[t_{i-1}\] et \[t_i\]. Suite à l'analyse, une prévision est effectuée jusqu'à l'instant \[t_{i+1}\]. Les différences entre les deux prévisions à l'instant \[t_{i+1}\] sont calculées. Ce processus est répété à partir de l'instant \[t_i\] et ainsi de suite. Toutes les différences permettent alors d'estimer la matrice de covariances d'erreur d'ébauche \[\mathbf{B}\].

Comme le met en évidence Berre \etal (2006), l'erreur d'analyse s'écrit

\[ {\epsilon}^a_{i}={\epsilon}^b_{i}+\mathbf{K}({\epsilon}^o_{i}+H{\epsilon}^b_{i})\],

\[ {\epsilon}^a_{i}=(\mathbf{I}-\mathbf{K} H){\epsilon}^b_{i}+\mathbf{K}{\epsilon}^o_{i}\],

et la différence entre les conditions initiales des deux prévisions à l'instant \[t_i\] s'écrit

\[ \delta \mathbf{x}^a_i=\mathbf{x}^a_i-\mathbf{x}^b_i\],

\[ \delta \mathbf{x}^a_i=\mathbf{K}({\epsilon}^o_i-H{\epsilon}^b_i)\],

\[ \delta \mathbf{x}^a_i=-\mathbf{K} H{\epsilon}^b_i+\mathbf{K}{\epsilon}^o_i\].

Les matrices \[\mathbf{I}-\mathbf{K} H\] et \[\mathbf{K}\] représentent les poids des erreurs d'ébauche et d'observation dans l'équation d'analyse. Dans la méthode NMC, le poids de l'erreur d'ébauche est approximé par \[-\mathbf{K} H\]. Cette approximation est raisonnable si \[\mathbf{K} \sim \mathbf{I} / 2\] (Boutier, 1994). Ce cas de figure est décrit par un réseau d'observations très denses (\[H \sim \mathbf{I}\]) et des matrices de covariances d'erreur presque identiques (\[\mathbf{R} \sim H\mathbf{B} H^T \sim \mathbf{B}\]). La deuxième condition signifie que l'erreur d'observation possède la même intensité et les mêmes structures spatiales que l'erreur d'ébauche.

Cependant, dans les régions pauvres en observations ou avec des observations de piètre qualité, l'incrément d'analyse risque d'être faible tandis que l'erreur d'analyse sera grande.

De plus, l'erreur d'observation est généralement moins corrélée spatialement que l'erreur d'ébauche. Comme le montre Daley (1991, section 4.5), l'opérateur \[\mathbf{K} H\] agit comme un filtre passe-bas. Par conséquence, l'opérateur \[\mathbf{I}-\mathbf{K} H\] devrait agir comme un filtre passe-haut. Ce qui signifie que l'incrément d'analyse doit avoir un spectre plus large que l'erreur d'analyse. Les corrélations d'erreur d'analyse risquent donc d'être surestimées avec la méthode NMC.

La perturbation d'analyse à l'instant \[t_i\] est ensuite propagée par le modèle d'évolution jusqu'à \[t_{i+1}\]. Les différences entre les prévisions permettent donc d'estimer l'erreur d'ébauche à la condition que \[\mathbf{K} \sim \mathbf{I} / 2\].

La méthode NMC a de nombreux avantages. Elle permet d'obtenir des statistiques dans l'espace du modèle et donc pour toutes les variables du modèle. De plus, elle est très bon marché. Cependant, elle ne représente pas parfaitement l'erreur d'ébauche car les hypothèses faites ne sont pas respectées. Ainsi, l'estimation de l'erreur d'ébauche est trop faible dans les régions peu ou mal observée.

Méthode d'ensemble

La méthode d'ensemble a d'abord été proposé par Evensen (Evensen, 1994) dans le cadre du filtre de Kalman d'ensemble. Néanmoins, cette méthode peut s'appliquer aux autres méthodes d'assimilation.

L'idée de cette méthode est de construire un ensemble composé d'une série de membres perturbés. Chacun des membres est analysé puis propagé de fenêtre d'assimilation en fenêtre d'assimilation. Ainsi, chaque membre est traité individuellement. Il est alors possible de calculer des différences entre ces membres à n'importe quel instant, puis d'obtenir des statistiques sur ces différences. La figure ens_algo permet d'illustrer l'algorithme.

Fig. ens_algo : Méthode d'ensemble. Un ensemble est constitué de \[n\] membres perturbés qui analysés et propagés indépendamment. Après chaque cycle d'assimilation, les différences entre ces membres permettend d'obtenir des statistiques estimant la matrice de covariances d'erreur d'ébauche \[\mathbf{B}\].

Il existe un lien entre les statistiques obtenues avec les différences entre les membre et l'erreur d'ébauche. En effet, les perturbations ajoutées aux membres de l'ensemble évoluent de manière similaire à l'erreur du système d'assimilation. Ainsi, à condition de bien spécifier les perturbations, les statistiques obtenues avec les différences entre les membres sont une très bonne estimation de l'erreur d'ébauche.

Cependant, il est difficile de bien perturber les membres de l'ensemble, car les perturbations appliquées aux divers champs doivent être similaires aux covariances d'erreur de ces champs. Le problème de la connaissance de la matrice de covariances d'erreur d'ébauche est ainsi déplacé vers la connaissance des matrices de covariances d'erreur des champs perturbés. Néanmoins, ces champs à perturber peuvent être mieux connus ou leurs covariances d'erreur plus accessible.

La méthode d'ensemble est donc une méthode complexe et coûteuse. Elle a cependant des attraits non-négligeables. Elle permet d'estimer réellement les erreurs d'ébauche de toutes les variables du modèle au cours du temps. Il est ainsi possible d'obtenir une matrice de covariance d'erreur d'ébauche \[\mathbf{B}\] dynamique. Cette méthode a cependant un défaut important. Si le système d'analyse est bruité, les statistiques le seront aussi et amplifieront le bruit du système d'analyse. Il est donc nécessaire d'être attentif aux risques de rétroaction.

Modèle de covariance d'erreur - Description des erreurs

jOas — Fri, 16 Nov 2007 01:07:00 +0100

Variances d'erreur d'observation

Les variances d'erreur d'observation sont le plus souvent estimées à l'aide des connaissances sur les caractéristiques techniques des instruments de mesures. Ces caractéristiques peuvent être, par exemple, déterminées par des observations positionnées au même endroit. Ces variances doivent aussi inclure les variances d'erreur de représentativité qui sont loin d'être négligeables tant qu'il existe des phénomènes physiques qui ne sont pas bien représentés dans l'espace du modèle.

D'autre part, il ne faut absolument pas considérer qu'un biais puisse être considéré comme une contribution aux variances d'erreur d'observation. En effet, il occasionnerait un biais dans l'incrément d'analyse. Ainsi, à chaque fois qu'un biais est mis en évidence, il doit être retiré des observations ou de l'état d'ébauche en fonction de son origine supposée. Il est cependant souvent difficile de déterminer son origine.

Covariances d'erreur d'observation

Les covariances d'erreur d'observation sont le plus souvent considérées comme nulles. En d'autres termes, des mesures distinctes sont affectées par des erreurs physiques indépendantes. Cette hypothèse paraît raisonnable pour des observations mesurées par des instruments différents. Elle paraît moins évidente quand un jeu d'observations est obtenu par le même instrument de mesure (mesures satellite, bouées dérivant es...) ou quand une série temporelle de mesures d'une même station, par exemple par une bouée fixe, est utilisée dans un 4D-Var. Dans de tels cas, il apparaît intuitivement que des corrélations d'erreur doivent exister entre ces mesures proches les une des autres.

La présence d'un biais, par exemple, se traduit par une corrélation d'erreur permanente. De plus, le prétraitement des observations peut produire artificiellement des corrélations d'erreur entre les observations transformées. Ces transformations sont assez courantes et permettent, par exemple, de transformer des températures en températures potentielles, la distance altimétrique d'un satellite en anomalie de hauteur de mer ou des radiances obtenues par des satellites en température de surface de la mer.

Quand l'état d'ébauche est utilisé dans le prétraitement des observations, cela peut créer artificiellement des corrélations entre l'erreur d'observation et l'erreur d'ébauche qui sont très difficile à prendre en compte. Par exemple, rapprocher une observation de l'état d'ébauche donne l'impression de réduire les erreurs d'observation et d'ébauche, mais cela réduit de manière irréaliste le poids de l'information apportée par l'observation originale.

Enfin, les erreurs de représentativité sont corrélées par nature. Les d'erreurs d'interpolation sont toujours corrélées quelque soit la densité des observations vis-à-vis de la résolution du modèle. D'autre part, les erreurs dans la définition de l'opérateur d'observation, comme le modèle d'évolution pour le 4D-Var, sont corrélées aux mêmes échelles que le problème modélisé.

La présence de corrélations d'erreur d'observation positives réduit le poids moyen des observations. Ainsi, elle augmente l'importance relative des différences entre les observations, comme les gradients ou les tendances. Cependant, les corrélations d'erreur d'observation sont très difficiles à estimer et peuvent causer des problèmes numériques lors de l'analyse ou dans les algorithmes de contrôle qualité. En pratique, il est souvent plus facile de minimiser leurs impacts en utilisant des méthodes de corrections de biais, en évitant trop de prétraitement des observations, en refusant des observations dans les régions très denses et en améliorant le modèle d'évolution et les opérateurs d'observation. Ainsi, les modèles de covariances d'erreur d'observation sont, en général, diagonaux ou presque.

Variances d'erreur d'ébauche

En général, l'ébauche provient d'une prévision obtenue avec le modèle d'évolution. Les variances d'erreur d'ébauche sont donc les variances d'erreur de la prévision utilisée pour obtenir l'état d'ébauche initial \[\mathbf{x}^b\]. Avec les filtres de Kalman, ces variances d'erreur sont estimées automatiquement à travers le modèle tangent-linéaire. Il n'est donc pas nécessaire de les spécifier. Cependant, le problème n'est que déplacé puisqu'il faut alors spécifier l'erreur modèle \[\mathbf{Q}\] et, pour les filtres réduits, mettre au point les approximations nécessaires à ces algorithmes.

Une première estimation des variances d'erreur d'ébauche peut être obtenue en prenant une fraction des variances climatologiques des champs des variables du vecteur d'état.

Une autre possibilité est d'utiliser des quantités dont les statistiques d'erreur sont équivalentes à celle de l'erreur d'ébauche. Parmi ces méthodes, les plus connues sont la méthode NMC et la méthode d'ensemble décrites par la suite. Une des hypothèses de ces méthodes est que l'analyse soit de bonne qualité. C'est-à-dire, en d'autres termes, que les observations soient nombreuses.

Enfin, la méthode la meilleure est, sans conteste, celle qui utilise l'innovation pour estimer les variances d'erreurs. Cette méthode sera aussi décrite par la suite mais elle repose sur l'hypothèse que les erreurs d'observation ne sont pas corrélées.

Néanmoins, dans la plupart des problèmes, l'erreur d'ébauche est supposée dépendre largement de l'état lui-même de l'ébauche. Il est alors très bénéfique que l'erreur d'ébauche dépende de l'écoulement et prenne en compte les variations temporelles. Cette caractéristique est obtenue avec les filtres de Kalman, dans les fenêtres temporelles des 4D-Var, par l'utilisation de lois empiriques basées les connaissances physiques du système ou par des méthodes d'ensemble (ou équivalentes).

Si les variances d'erreurs d'ébauches sont mal spécifiées, l'incrément d'analyse sera ou trop grand ou trop petit. Avec les algorithmes variationnelles basées sur la méthode des moindres carrés, seul le rapport des variances d'erreur d'ébauche et d'observation est important pour l'analyse. Néanmoins, la valeur absolue peut aussi avoir son importance lors du contrôle de la qualité des observations. En effet, les observations éloignées de l'ébauche peuvent être tout à fait acceptées si l'erreur d'ébauche est importante.

Covariances d'erreur d'ébauche

Les covariances d'erreur d'ébauche sont essentielles pour faire une bonne analyse.

Propagation des informations

Dans les régions pauvres en observations, la forme de l'incrément d'analyse est complètement déterminée par les structures de covariances d'erreur d'ébauche. Ainsi, la forme de l'incrément d'analyse d'une observation esseulée est donnée directement par \[\mathbf{B}\mathbf{H}^T\]. C'est donc les corrélations de \[\mathbf{B}\] qui propagent l'information spatialement autour du point d'observation.

Lissage des informations

Dans les régions riches en observations, le lissage des informations est gouverné par les corrélations de la matrice de covariances d'erreur d'ébauche \[\mathbf{B}\]. Ceci est clairement mis en évidence au regard du gain d'analyse optimal \[\mathbf{K}\] (Eq. (015)) dont le dernier terme à être utilisé, celui le plus à gauche, est \[\mathbf{B}\]. Le lissage de l'incrément d'analyse est très important en ceci qu'il doit permettre à l'analyse d'avoir des échelles statistiquement compatibles avec les propriétés des champs physiques. La spécification des corrélations d'erreur d'ébauche est donc à la fois importante et délicate, car les échelles spatiales des champs physiques sont diverses et variables.

Propriétés d'équilibre

Le nombre de degrés de liberté d'un modèle est souvent supérieur à celui de la réalité. Par exemple, les courants marins sont supposés géostrophique hormis à l'équateur. Cette hypothèse d'équilibre entre gradient de pression et force de Coriolis dans les équations de Navier-Stokes permet d'obtenir une relation directe entre courant géostrophique et dérivée au premier ordre de la hauteur de mer. Ces propriétés d'équilibre peuvent être considérées comme des contraintes gênantes au problème d'analyse et être a posteriori appliquées brutalement. Un autre point de vue est de considérer qu'il existe des propriétés statistiques qui lient les différentes variables du modèle. S'il existe des relations d'équilibre entre les différentes variables du modèle, il doit donc y avoir aussi des relations d'équilibre linéarisés dans la matrice de covariances d'erreur d'ébauche. Ces équilibres sont très intéressants car ils permettent d'apporter des informations sur toutes les variables en équilibre avec celle observée. Ainsi, par exemple, les mesures de température en océanographie permettent de corriger la salinité. Combiné avec le lissage spatial des informations, les propriétés d'équilibre peuvent avoir un impact considérable sur la qualité de l'analyse. Une mesure de température propagée autour du point d'observation peut modifier, en plus de la température, la salinité, la hauteur de mer et les courants. C'est-à-dire toutes les variables utilisées dans un modèle océanique. L'amplitude de toutes ces modifications dépendra de l'estimation de la corrélation entre deux variables différentes et des variances d'erreur de ces variables.

Paramètres additionnels

Avec les méthodes variationnelles, il est possible d'inclure dans le vecteur de contrôle d'autres paramètres additionnels, tels que des paramètres de réglage du modèle ou des estimations du biais. Cette technique indirecte d'estimation des paramètres peut être très efficace à condition qu'il y ait une véritable relation entre ces paramètres et les observations. Cette relation passe, en général, au travers de l'opérateur d'observation ou du modèle dans le cas du 4D-Var. Il n'est généralement ni possible ni prudent de spécifier explicitement les corrélations avec les autres variables de l'état du modèle dans \[\mathbf{B}\]. De plus, les erreurs d'ébauche de tous les paramètres du vecteur de contrôle doivent être spécifiées prudemment, à moins d'être sûr que le problème est surdéterminé par les observations. Une variance d'erreur trop faible évitera, logiquement, de corriger les paramètres additionnels. Une variance d'erreur trop forte pourra, par contre, transformer les paramètres additionnels en source de bruit et créer des variations sans justification physique. D'importants problèmes peuvent alors apparaître car des couplages implicites dans l'analyse sont souvent créés par des dépendances dans l'opérateur d'observation ou dans le modèle en 4D-Var. Ainsi, la spécification des erreurs d'ébauche des paramètres additionnels peut avoir un impact sur les variables d'état du modèle analysées.

Dépendance à l'écoulement

La matrice de covariances d'erreur d'ébauche \[\mathbf{B}\] peut dépendre de l'incertitude d'une précédente analyse ou prévision à condition que la dynamique du problème soit suffisamment connue. Non seulement les variances d'erreur d'ébauche peuvent évoluer au cours du temps, mais aussi les corrélations. En atmosphère comme en océanographie, certaines ondes suivent des motifs spécifiques qui peuvent apparaître dans les erreurs d'ébauche. Par exemple, dans un secteur enclin au développement cyclonique (région de basses pressions), les erreurs d'ébauche le plus probables devraient avoir la forme des structures les plus instables, avec peut-être une inclinaison des ondes baroclines et des anti-corrélations entre les masses d'air chaud et froid. C'est dépendance à l'écoulement est assez équivalent aux propriétés d'équilibre. Ainsi, si de tels informations peuvent être incorporées dans la matrice de covariances d'erreur \[\mathbf{B}\], alors les observations pourront être mieux propagées spatialement et mieux distribuées entre les variables du modèle. Ce type d'information peut être utilisé dans le cadre des filtres de Kalman, en 4D-Var mais aussi en 3D-Var comme le mettra en évidence ce manuscrit.

Modèle de covariance d'erreur

jOas — Fri, 16 Nov 2007 01:05:00 +0100

Une spécification correcte et adaptée des covariances d'erreur d'observation et d'ébauche est primordiale pour la qualité de l'analyse. En effet, ces covariances déterminent comment et à quel point les observations corrigeront l'état d'ébauche. Les variances d'erreur sont les paramètres essentiels. Pour autant, les corrélations sont aussi très importantes car elles déterminent comment les informations apportées par les observations seront lissées dans l'espace du modèle s'il y a un décalage entre la résolution du modèle et la densité d'observations. Dans le cadre des filtres de Kalman ou du 4D-Var à contrainte faible, les covariances d'erreur modèle sont aussi à spécifier.

Méthode d'assimilation - 4DVar

jOas — Fri, 16 Nov 2007 00:14:00 +0100

Le 4D-Var est l'extension temporelle du 3D-Var. Cette méthode ne vise pas à obtenir l'état optimal à un instant donné, mais la trajectoire optimale sur une fenêtre de temps donné. Les observations sont donc prises en compte aussi bien dans leur distribution spatiale que temporelle. Cet aspect est déjà pris en compte par le 3D-Var FGAT présenté dans la section ad-hoc. Néanmoins, le 4D-Var apporte un aspect temporel en plus car il propage l'information apportée par les observations à l'instant initial de la fenêtre d'assimilation. De ce fait, l'analyse obtenue doit permettre au modèle d'évolution d'avoir la trajectoire la plus proche possible de l'ensemble des observations utilisées.

Cette amélioration du 3D-Var permet d'ajouter la connaissance de l'évolution du système comme information pour l'analyse.

De nombreuses applications à des modèles réalistes météorologiques (Thépaut et Courtier, 1991 et Zupanski, 1993) et océanographiques (Moore, 1986 ; Shröter etal., 1993 ; Luong etal., 1998 et Greiner etal., 1998) ont depuis longtemps été effectuées.

L'amélioration ainsi apportée, conjuguée au fort développement des moyens de calculs, font que le 4D-Var est venu remplacer le 3D-Var dans les systèmes de prévision opérationnels atmosphériques du CEPMMT en 1997 et de Météo-France en 2000.

4D-Var classique

Soit \[M_{0 \to i}\] l'opérateur a priori non-linéaire qui permet de propager l'état du système \[\mathbf{x}\] de \[t_0\] à \[t_i\] :

\[ \forall i, \qquad \mathbf{x}(t_i) = M_{0 \to i}(\mathbf{x})\].

En supposant, dans un premier temps, le modèle parfait, la fonction coût \[J\] du 4D-Var se décompose, comme pour le 3D-Var, en un terme \[J^b\] lié à l'ébauche et un autre \[J^o\] lié aux observations.

En tenant compte des instants de mesures \[t_i\], la matrice de covariance d'erreur des observation est notée \[\mathbf{R}_i\] et l'opérateur d'observation non-linéaire \[H_i\]. Les observations \[\mathbf{y}^o_i\] sont donc comparées à leur équivalent modèle \[H_i\mathbf{x}(t_i)\] à chaque instant d'observation. Le calcul du terme \[J^o\] nécessite l'intégration du modèle d'évolution de \[t_0\] à \[t_N\]. Le vecteur d'état \[\mathbf{x}\] est ainsi propagé par le modèle numérique \[M_{0 \to N}\] de \[t_0\] à \[t_N\] où \[N\] représentent le nombre de pas de temps de l'intégration du modèle à l'intérieur d'un cycle d'assimilation.

L'algorithme d'assimilation identifie un état \[\mathbf{x}^a\] de la variable \[\mathbf{x}\] à l'instant \[t_0\] (une condition initiale), qui, intégré par le modèle d'évolution fournit une trajectoire optimale au sens des moindres carrés (la trajectoire analysée) sur l'ensemble de la fenêtre d'assimilation (Les vecteurs \[\mathbf{x}^a(t_0)\] et \[\mathbf{x}^b(t_0)\] sont notés \[\mathbf{x}^a\] et \[\mathbf{x}^b\]. \`A tout autre moment que l'instant initial \[t_0\], les notations \[\mathbf{x}^a(t_i)\] et \[\mathbf{x}^b(t_i)\] seront utilisées.) :

(051)

\[ J^b(\mathbf{x}) = \frac{1}{2}\left(\mathbf{x}-\mathbf{x}^b\right)^T\mathbf{B}^{-1}\left(\mathbf{x}-\mathbf{x}^b\right)\],

\[ J^o(\mathbf{x}) = \frac{1}{2}\sum_{i=0}^N \left(\mathbf{y}^o_i - H_i\mathbf{x}(t_i)\right)^T\mathbf{R}_i^{-1} \left(\mathbf{y}^o_i - H_i\mathbf{x}(t_i)\right)\],

\[J^o(\mathbf{x}) = \frac{1}{2}\sum_{i=0}^N \left(\mathbf{y}^o_i - H_i M_{0\to i}(\mathbf{x})\right)^T\mathbf{R}_i^{-1} \left(\mathbf{y}^o_i - H_i M_{0\to i}(\mathbf{x})\right)\],

(052)

\[ J^o(\mathbf{x}) = \frac{1}{2}\sum_{i=0}^N \left(\mathbf{y}^o_i - G_i\mathbf{x}\right)^T\mathbf{R}_i^{-1} \left(\mathbf{y}^o_i - G_i\mathbf{x}\right)\],

où \[G_i\mathbf{x}=H_i\mathbf{x}(t_i)=H_i M_{0\to i}(\mathbf{x})\]. L'opérateur \[G_i\] est appelé l'opérateur d'observation généralisé pour l'état \[\mathbf{x}\] propagé par le modèle de \[t_0\] à \[t_i\].

L'état optimal \[\mathbf{x}^a\], qui minimise la fonction coût \[J\], est obtenu quand le gradient de cette fonctionnelle est nul. Comme pour le 3D-Var, ce gradient s'obtient simplement :

(053)

\[ \nabla J(\mathbf{x}) = \mathbf{B}^{-1}\left(\mathbf{x}-\mathbf{x}^b\right) - \sum_{i=0}^N \mathbf{M}^T_{0\to i}\mathbf{H}^T_i\mathbf{R}_i^{-1} \left(\mathbf{y}^o_i - G_i\mathbf{x}\right)\].

Les opérateurs \[\mathbf{H}_i\], \[\mathbf{M}_{0\to i}\] et \[\mathbf{G}_i=\mathbf{H}_i \mathbf{M}_{0\to i}\] sont les opérateurs linéarisés de \[H_i\], \[M_{0\to i}\] et \[G_i\] au voisinage de l'ébauche. L'opérateur \[\mathbf{M}^T_{0\to i}\] est l'adjoint de l'opérateur linéarisé \[\mathbf{M}_{0\to i}\].

Équivalence avec le filtre de Kalman

Si les opérateurs \[H\] et \[M\] sont linéaires, alors la fonction coût \[J\] est quadratique. Si de plus le modèle est parfait (hypothèse émise dans la section précédente), alors la solution du 4D-Var à la fin de la fenêtre d'assimilation est identique à celle du filtre de Kalman (Jazwinski, 1970 ; Ghil etal., 1981 et Lorenc, 1986). En météorologie comme en océanographie, \[H\] et \[M\] sont souvent faiblement non-linéaires. Dans ce cas, la minimisation peut être effectuée avec un algorithme adapté aux fonctions coûts non-quadratiques. Généralement, l'opérateur généralisé d'observation linéarisé \[\mathbf{H}_i\] et le modèle linéaire-tangent \[\mathbf{M}_{0\to i}\] sont supposés de bonnes approximations de \[H_i\] et \[M_{0\to i}\] sur la fenêtre temporelle d'assimilation. La validité du linéaire-tangent dépend d'une part de la formulation du modèle numérique et de l'opérateur d'observation considérés mais aussi du contexte de l'assimilation, notamment la durée de la fenêtre d'assimilation, de la physique des phénomènes représentés et de la région d'étude.

Calcul de la fonction coût et de son gradient

Le terme d'ébauche de la fonction coût \[J^b\] est identique à celui décrit pour le 3D-Var. Son évaluation est directe pour tout état de la variable d'état \[\mathbf{x}\]. L'évaluation du terme lié aux observations \[J^o\] est, par contre, plus ardue. Dans le cadre du 3D-Var, \[J^o\] est une combinaison linéaire des écarts entre les observations et l'état du modèle à un instant donné. \`A présent, chaque évaluation de \[J^o\] requiert l'intégration du modèle sur la fenêtre d'assimilation.

L'équation \[\nabla J=0\] ne peut être résolue directement. Une solution minimisant la fonction coût \[J\] par une méthode de descente itérative utilisant la valeur de \[\nabla J\] est envisageable. Généralement, l'état d'ébauche \[\mathbf{x}^b\] est utilisé comme une première estimation de l'état analysé. Le gradient de la fonction coût \[\nabla J\] peut être estimé de manière très efficace par la méthode adjointe (Une description complète des méthodes adjointes est présentée dans Thacket et Long (1988). En fait, l'évaluation du terme \[\nabla J^o\] peut se résumer à une intégration du modèle direct et une intégration du modèle adjoint (Le Dimet et Talagrand, 1986). Le gradient de la fonction coût est évalué par rapport à \[\mathbf{x}\] décrivant le vecteur d'état à l'instant initial de la fenêtre d'assimilation

(054)

\[ \nabla J^o(\mathbf{x}) = \mathbf{x}^*\],

où \[\mathbf{x}^*\] représente l'adjoint de la variable \[\mathbf{x}\].

La méthode adjointe permet une évaluation efficace du gradient de la fonction coût dans l'algorithme du 4D-Var. Néanmoins, cette opération implique qu'une version adjointe du modèle d'évolution linéarisé soit disponible. Cette contrainte est souvent l'étape cruciale à surmonter lors de l'implémentation de l'algorithme d'assimilation 4D-Var. En effet, dans de nombreux domaines, comme la météorologie et l'océanographie, les modèles d'évolution sont souvent complexes et non-linéaires. Deux étapes sont généralement nécessaire : écrire le modèle linéaire-tangent du modèle d'évolution ; puis écrire son adjoint. Ces opérations peuvent être effectuées soit manuellement, soit par des méthodes de linéarisation et d'adjoint automatiques (Giering et Kaminski, 1998). Malheureusement, ces méthodes automatiques ne sont pas capables de présumer les hypothèses de linéarisation qui sont faites manuellement et les modèles linéaire-tangent et adjoint obtenus sont souvent très coûteux. Néanmoins, l'utilisation conjointe d'une méthode automatique et d'une écriture manuelle permet d'obtenir assez rapidement des modèles tangent-linéaire et adjoint efficaces.

Coût de calcul du 4D-Var

L'intégration du modèle direct pour le calcul des innovations, puis pour la propagation de l'analyse, ainsi que l'intégration du modèle adjoint font du 4D-Var une méthode très coûteuse. Elle est néanmoins applicable en météorologie et, dans une moindre mesure, en océanographie. En effet, la plupart des grands centres opérationnels de météorologie utilisent des méthodes 4D-Var. Néanmoins, à la différence des centres océanographiques, leur capacité de calcul est très importante. Ainsi, le 4D-Var a été utilisé dans de nombreuses études en océanographie, notamment pour assimiler des données altimétriques dans des modèles régionaux quasi-géostrophiques (Moore, 1991), dans des modèles à gravité réduite (Weaver et Anderson, 1997), dans des modèles "shallow-water" (Greiner et Périgaud, 1994), dans des modèles aux équations primitives (Stammer etal., 1997 ; Greiner etal., 1998a, b) et dans des modèles couplés (Lee etal., 2000).

Le 4D-Var à contrainte faible

Dans la présentation du 4D-Var, le modèle a été considéré comme parfait. C'est-à-dire que le modèle décrivait exactement le comportement du système. Dans ce cas, le 4D-Var est décrit comme à contrainte forte (Sasaki, 1970). Cependant, malgré l'utilisation de modèles extrêmement sophistiqués, ceux-ci comportent des erreurs qui ne peuvent être négligées et qui ne pourront jamais l'être pour des systèmes aussi complexes que l'atmosphère ou l'océanographie.

Comme pour le filtre de Kalman, l'ajout d'un terme correctif directement dans la fonction coût est possible. En définissant l'erreur du modèle telle que

(055)

\[\mathbf{q}_i=\mathbf{x}(t_i) - M_{i\to i+1}\mathbf{x}(t_{i-1})\],

la fonction coût \[J\] peut alors s'écrire

(056)

\[ J(\mathbf{x}, \mathbf{q}_1,\cdots,\mathbf{q}_N) = \frac{1}{2}\left(\mathbf{x}-\mathbf{x}^b\right)^T\mathbf{B}^{-1}\left(\mathbf{x}-\mathbf{x}^b\right)\]
\[+\frac{1}{2}\sum_{i=0}^N \left(\mathbf{y}^o_i-G_i\mathbf{x}\right)^T\mathbf{R}_i^{-1} \left(\mathbf{y}^o_i-G_i\mathbf{x}\right)\]
\[+\frac{1}{2}\sum_{i=0}^N\mathbf{q}_i^T\mathbf{Q}_i^{-1}\mathbf{q}_i\],

où \[\mathbf{Q}_i\] est la matrice de covariance d'erreur modèle à l'instant \[t_i\].

Cette formulation du 4D-Var est dite à contrainte faible et il est nécessaire de proposer une modélisation de la matrice de covariance d'erreur modèle \[\mathbf{Q}\]. De nombreux travaux par Derber (1989), Stammer etal. (1997), Bennet etal. (1998), Lee et Marotzke (1998) ou Vidard (2001) utilisent cette algorithme du 4D-Var à contrainte faible.

Évolution temporelle de la matrice de covariance d'erreur d'ébauche

L'impact de la prise en compte du caractère temporel est clairement mis en évidence dans la fonction coût liée aux observations \[J^o\]. Cependant, l'aspect temporel de la matrice de covariance d'erreur d'ébauche \[\mathbf{B}\] est moins évident dans un premier temps. En effet, la dynamique du modèle d'évolution est prise en compte implicitement par la matrice de covariance d'erreur d'ébauche \[\mathbf{B}\] sur chaque fenêtre d'assimilation, et donc dans la fonction coût liée à l'ébauche \[J^b\].

L'erreur d'ébauche à l'instant \[t_0\] est l'écart entre l'ébauche et l'état vrai au début de la fenêtre d'assimilation \[{\epsilon}^b=\mathbf{x}^b-\mathbf{x}^t\]. \`A tout instant \[t_i\] du cycle d'assimilation, cette erreur notée \[{\epsilon}^b(t_i)\], est l'écart entre l'état vrai \[\mathbf{x}^t(t_i)\] et l'ébauche \[\mathbf{x}^b(t_i)\] sous l'hypothèse que le modèle est linéaire :

(057)

\[ {\epsilon}^b (t_i)= \mathbf{x}^b(t_i) - \mathbf{x}^t(t_i)\]

\[{\epsilon}^b (t_i)=M_{0\to i}\mathbf{x}^b - M_{0\to i}\mathbf{x}^t\]

\[{\epsilon}^b(t_i)\approx\mathbf{M}_{0\to i}\left(\mathbf{x}^b-\mathbf{x}^t\right)\]

\[{\epsilon}^b(t_i)\approx\mathbf{M}_{0\to i}{\epsilon}^b\].

La matrice de covariance d'erreur à l'instant \[t_i\] est alors donnée par

(058)

\[ E\left[{\epsilon}^b(t_i)({\epsilon}^b(t_i))^T\right] = E\left[\mathbf{M}_{0\to i}{\epsilon}^b({\epsilon}^b)^T\mathbf{M}_{0\to i}^T\right]\]

\[ E\left[{\epsilon}^b(t_i)({\epsilon}^b(t_i))^T\right]=\mathbf{M}_{0\to i}\mathbf{B}\mathbf{M}_{0\to i}^T\].

Cette matrice décrit les erreurs liées à l'ébauche aux instants d'observations. La matrice de covariance d'erreur d'ébauche \[\mathbf{B}\] est donc implicitement propagée en temps par le 4D-Var à travers la dynamique du modèle linéaire-tangent \[\mathbf{M}_{0\to i}\] et son adjoint \[\mathbf{M}_{0\to i}^T\]. Le 4D-Var présente donc une analyse cohérente avec la dynamique du système (Thépaut etal., 1993).

Équivalence avec le 3D-Var FGAT

Il existe une similarité frappante entre le 3D-Var FGAT (Eqs. (048a} et (048b}) et le 4D-Var à contrainte forte (Eqs. (051} et (052}). En effet, il suffit de définir le modèle d'évolution dans le 4D-Var comme l'identité (\[M=\mathbf{I}\]), pour que le l'opérateur généralisé d'observation \[G\] se réduise à l'opérateur d'observation non-linéaire \[H\] et que le 4D-Var à contrainte forte devienne un 3D-Var FGAT.

C'est une caractéristique très intéressante car elle permet très facilement, du point de vue de l'implémentation informatique, de passer d'une méthode à l'autre. Étant donné le coût informatique du 4D-Var, cette approche permet une évolution simple et naturelle du 3D-Var FGAT vers le 4D-Var en fonction de l'évolution des capacités informatiques disponibles.

Formulation incrémentale du 4D-Var

La formulation incrémentale a déjà été abordée rapidement à propos du 3D-Var dans la section ad-hoc, mais elle prend avec le 4D-Var tout son sens.

En effet, l'introduction de l'approche incrémentale en météorologie a été motivée par la réduction de coût qu'elle propose. Dans le cadre du 4D-Var classique décrit précédemment, à chaque itération de la minimisation de la fonction coût (Eqs. (051} et (052}), l'intégration du modèle direct non-linéaire et de l'adjoint du modèle linéarisé est très coûteuse. Les non-linéarités des modèles numériques atmosphériques peuvent conduire à des fonctions coûts complexes. Les minimiseurs utilisés sur ces fonctions n'aboutissent pas forcément à une minimisation fiable et ce, de surcroît, à un coût élevé. Ces non-linéarités compliquent de plus, lourdement la tâche de l'écriture de l'adjoint (Xu, 1996). Dans l'approche incrémentale, la fonction coût est rendue quadratique, ce qui garantit l'identification d'un minimum unique par une méthode de descente pour un coût de calcul notablement inférieur à celui du problème non-linéaire. Le modèle linéaire-tangent est écrit avec une physique simplifiée, ce qui facilite grandement l'écriture de l'adjoint du modèle. De plus, une approche communément choisie en météorologie est d'utiliser un modèle linéaire-tangent à une résolution plus basse que celle du modèle non-linéaire. Le coût de la minimisation du 4D-Var en est significativement réduit. C'est formulation incrémentale de l'approche variationnelle du 4D-Var qui a permis de le rendre applicable de façon opérationnelle pour la prévision météorologique (Courtier etal., 1994 et Rabier etal., 2000).

Dans la formulation incrémentale de l'assimilation variationnelle, l'objectif est de minimiser la fonction coût, non plus par rapport à la variable d'état \[\mathbf{x}\], mais par rapport à un incrément \[\delta\mathbf{x}\] tel que \[\mathbf{x}=\mathbf{x}^b+\delta\mathbf{x}\]. L'hypothèse principale de la formulation incrémentale est d'utiliser un modèle d'évolution et un opérateur d'observation linéarisé pour propager l'incrément mais de conserver le modèle non-linéaire pour la propagation de l'ébauche \[\mathbf{x}^b\]. La solution de la minimisation est l'incrément d'analyse \[\delta\mathbf{x}^a\] à \[t_0\] tel que le vecteur d'analyse \[\mathbf{x}^a\] soit

(059)

\[ \mathbf{x}^a=\mathbf{x}^b+\delta\mathbf{x}^a\].

En supposant que l'ébauche \[\mathbf{x}^b\] est une "bonne" approximation a priori de l'état optimal du système au sens des moindres carrés, l'incrément \[\delta\mathbf{x}\] devrait être petit. Les opérateurs non-linéaires d'observation \[H\] et du modèle \[M\] sont linéarisés au voisinage de l'ébauche de sorte que pour tout état du modèle \[\mathbf{x}\], à chaque instant \[t_i\] de la fenêtre temporelle d'assimilation,

\[ \mathbf{x}(t_{i+1})=M_{0\to i}\mathbf{x}\] \[ \mathbf{x}(t_{i+1})=M_{0\to i}(\mathbf{x}^b+\delta\mathbf{x})\]

(060)

\[ \mathbf{x}(t_{i+1})\approx M_{0\to i}\mathbf{x}^b+\mathbf{M}_{0\to i}\delta\mathbf{x}\].

Par conséquent

(061)

\[H_i\mathbf{x}(t_i)=H_i\mathbf{x}^b(t_i)+\mathbf{H}_i \delta\mathbf{x}+\parallel\delta\mathbf{x}\parallel^2+\cdots\],

(062)

\[G_i\mathbf{x}(t_i)=G_i\mathbf{x}^b(t_i)+\mathbf{G}_i \delta\mathbf{x}+\parallel\delta\mathbf{x}\parallel^2+\cdots \],

où \[\mathbf{M}_{0\to i}\], \[\mathbf{H}_i\] et \[\mathbf{G}_i\] sont les opérateurs linéarisés à \[t_0\] autour de l'ébauche \[x^b\] de \[M_{0\to i}\], \[H_i\] et \[G_i\].

\[ \mathbf{M}_{0\to i}=\left.\frac{\partial M_{0\to i}}{\partial \mathbf{x}}\right|_{\mathbf{x}=\mathbf{x}^b} \],

\[ \mathbf{H}_i=\left.\frac{\partial H_i}{\partial \mathbf{x}}\right|_{\mathbf{x}=\mathbf{x}^b} \],

(063)

\[ \mathbf{G}_i=\left.\frac{\partial G_i}{\partial \mathbf{x}}\right|_{\mathbf{x}=\mathbf{x}^b}\].

En insérant les Eqs. (060}, (061} et (062} dans la formulation de la fonction coût donnée par les Eqs. (051} et (052}, elle peut être reformulée de manière incrémentale

(064)

\[ J^b(\mathbf{x}) = \frac{1}{2}\delta\mathbf{x}^T\mathbf{B}^{-1}\delta\mathbf{x}\],

(065)

\[J^o(\mathbf{x}) = \frac{1}{2}\sum_{i=0}^N \left(\mathbf{d}_i-\mathbf{G}_i \delta \mathbf{x} \right)^T\mathbf{R}_i^{-1} \left(\mathbf{d}_i-\mathbf{G}_i \delta \mathbf{x} \right)\],

où \[\mathbf{d}_i=\mathbf{y}^o_i-G_i\mathbf{x}^b=\mathbf{y}^o_i-H_i\mathbf{x}^b(t_i)\] représente l'innovation au temps \[t_i\], c'est-à-dire l'écart entre les observations et l'équivalent de l'ébauche donné par l'opérateur d'observation généralisé \[G_i\] dans l'espace des observations à chaque temps \[t_i\].

La fonction coût \[J=J^b+J^o\] du 4D-Var incrémental à contrainte forte est quadratique et la minimisation possède une solution unique. Si la linéarisation des opérateurs \[H_i\], \[M_{0\to i}\] et \[G_i\] est exacte, alors la solution est identique à celle obtenue par le filtre de Kalman étendu.

L'incrément d'analyse qui minimise la fonction coût donné par les Eqs. (064} et (065} est

(066)

\[ \delta\mathbf{x}^a= \mathbf{B}\mathbf{G}^T\left(\mathbf{G}\mathbf{B}\mathbf{G}^T+\mathbf{R}\right)^{-1}\mathbf{d}\].

La fonction coût incrémentale (Eqs. (064} et (065}) est minimisée par une méthode itérative de descente. Cette minimisation nécessite le calcul de la fonction coût et de son gradient à chaque itération de la minimisation comme pour l'algorithme classique du 4D-Var. Avant le début de la minimisation, l'état d'ébauche à \[t_0\] noté \[\mathbf{x}^b\] est propagé par le modèle non-linéaire permettant le calcul des innovations \[\mathbf{d}_i\] à chaque instant d'observation \[t_i\].

La figure 4DVar permet de représenter simplement l'utilisation du 4D-Var incrémental. Elle est à comparer aux Figs. 3DVar3 et 3DVar4 qui représentent, selon les mêmes codes, les 3D-Var FGAT incrémental avec ou sans IAU.

Fig. 4DVar : Illustration de la procédure pour cycler le 4D-Var incrémental. Pour chaque cycle \[c\], le modèle d'évolution est intégré de \[t_0\] à \[t_N\] à partir de l'état initial d'ébauche \[\mathbf{x}^b_c(t_0)\] (courbe noir pleine) et le vecteur d'innovation \[\mathbf{d}_i\] est calculé pour les différentes observations \[\mathbf{y}^o_i\] avec \[i=1,\cdots,N\] (ligne fine verticale). L'analyse est effectuée à l'instant \[t_0\] en ramenant les innovations à l'instant \[t_0\] à l'aide du modèle adjoint. Après la minimisation, un incrément est obtenu qui est rajouté à l'état de l'ébauche initial pour obtenir l'état analysé \[\mathbf{x}^a_c(t_0)=\mathbf{x}^b_c(t_0)+\delta\mathbf{x}^a\]. Cet état analysé tient compte de la dynamique du modèle de sorte que la trajectoire analysée (courbe grise pointillée) minimise au mieux l'écart aux observations tout au long du cycle d'assimilation (de \[t_0\] à \[t_N\]). L'état analysé \[\mathbf{x}^a_c(t_N)\] est ensuite utilisé comme état initial d'ébauche pour le cycle suivant.

Prise en compte des non-linéarités

A chaque itération de la minimisation, le terme de la fonction coût lié aux observations \[J^o(\delta\mathbf{x})\] est calculé en propageant l'incrément \[\delta\mathbf{x}\] dans le temps avec le modèle linéaire-tangent \[\mathbf{M}\]. Le calcul du gradient de la fonction coût, notamment de la partie relative aux observations \[\nabla J^o(\delta\mathbf{x})\], nécessite l'intégration de l'adjoint du modèle linéaire-tangent \[\mathbf{M}^T\] sur la fenêtre d'assimilation. A la fin de la minimisation, l'incrément d'analyse est ajouté à l'ébauche \[\mathbf{x}^b\] (Eq. (059}). L'état analysé à l'instant initial de la fenêtre d'assimilation \[\mathbf{x}^a\] est ensuite propagé par le modèle non-linéaire \[M\] jusqu'à la fin de la fenêtre permettant d'obtenir une trajectoire analysée \[\mathbf{x}^a(t_i)\]. En pratique, il est possible de prendre en compte les faibles non-linéarités des opérateurs \[H\] et \[M\] en mettant à jour la trajectoire de référence au cours de la minimisation. Ces mises à jours sont aussi appelée boucles externes. Le modèle linéaire est relinéarisé au voisinage du nouvel état de référence à chaque boucle externe et la fonction coût est ensuite minimisée par une série de boucles internes. Cette méthode permet de conserver la fonction coût quadratique tout en tenant compte, jusqu'à un certain point, des non-linéarités du système.

Méthode d'assimilation - 3DVar

jOas — Thu, 15 Nov 2007 23:51:00 +0100

3D-Var

La méthode d'assimilation variationnelle tri-dimensionnelle, notée 3D-Var pour "3Dimensional VARiational assimilation", consiste à chercher l'état le plus vraisemblable à partir des connaissances disponibles sur les lois de probabilités des erreurs d'observation et d'ébauche.

Comme sont nom l'indique clairement, le 3D-Var traite de problèmes tri-dimensionels. Par abus de langage, cette appellation est aussi utilisée pour des problèmes à une ou deux dimensions afin d'éviter les risques de confusions avec l'extension temporelle de cette méthode. En effet, sur un problème bi-dimensionnel, le 3D-Var s'appellerait 2D-Var, tandis que le 4D-Var se nommerait 3D-Var. Ce qui serait particulièrement ambigu. De ce fait, tous les problèmes ne prenant pas en compte l'aspect temporel sont appelés 3D-Var.

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. Le nouvel état analysé est, en général, utilisé comme point de départ de la prévision suivante.

En reprenant les Eqs. (023) à (029), il est possible d'écrire la fonction coût

\[ J(\mathbf{x}) = \frac{1}{2}(\mathbf{x}-\mathbf{x}^b)^T\mathbf{B}^{-1}(\mathbf{x}-\mathbf{x}^b)+\frac{1}{2}(\mathbf{y}^o-H\mathbf{x})^T\mathbf{R}^{-1}(\mathbf{y}^o-H\mathbf{x})\],

(045)

\[J(\mathbf{x}) = J^b(\mathbf{x})+J^o(\mathbf{x})\].

En général, \[\mathbf{x}^b\] est issu de l'intégration par le modèle d'évolution de l'état analysé à l'étape précédente. L'équilibre entre le terme d'écart aux observations \[J^o\] et celui de l'ébauche \[J^b\] est effectué grâce aux inverses des matrices de covariances d'erreur d'observation et d'ébauche. C'est-à-dire grâce à la confiance portée dans les observations et l'ébauche.

La minimisation se fait à l'aide du gradient égale à

(046)

\[ \nabla J(\mathbf{x}) = \mathbf{B}^{-1}(\mathbf{x}-\mathbf{x}^b)-\mathbf{H}^T\mathbf{R}^{-1}(\mathbf{y}^o-H\mathbf{x})\].

Comme montré dans la section sur le BLUE, si l'opérateur d'observation est linéaire (\[H=\mathbf{H}\]), le 3D-Var est alors équivalent au BLUE à l'optimalité, et donc aussi à l'interpolation optimale (Lorenc, 1986).

3D-Var classique

Dans un cadre théorique, il est possible d'imaginer des observations réparties régulièrement dans le temps. il est alors possible d'effectuer une analyse à chaque pas de temps où celles-ci sont présentes. En pratique, les observations sont diverses et très inégalement réparties. Afin de ne ma multiplier les phases d'analyse, une fenêtre temporelle de taille arbitraire est définie sur laquelle une analyse est effectuée. Dans le 3D-Var classique, toutes les observations sont regroupées (moyennées) à l'instant où est effectuée l'assimilation (Fig. 3DVar1). Dans ce cas, la méthode ne tient absolument pas compte de l'origine temporelle des observations hormis qu'elles doivent être incluses dans cette fenêtre temporelle. En faisant cette approximation, l'Eq. (045) s'écrit alors

\[ J^b(\mathbf{x}) = \frac{1}{2}(\mathbf{x}-\mathbf{x}^b)^T\mathbf{B}^{-1}(\mathbf{x}-\mathbf{x}^b)\],

(047)

\[J^o(\mathbf{x}) = \frac{1}{2}(\overline{\mathbf{y}}^o-\overline{H}\mathbf{x})^T\overline{\mathbf{R}}^{-1}(\overline{\mathbf{y}}^o-\overline{H}\mathbf{x})\].

Fig. 3DVar1 : Répartition des observations dans un 3D-Var classique.

3D-Var FGAT

Le 3D-Var FGAT (First Guess at Appropriate Time) est un 3D-Var pour lequel les observations sont utilisées en tenant compte de l'instant de leurs mesures (Fig. 3DVar2). C'est-à-dire que chaque observation est comparée à l'état du modèle au pas de temps le plus proche. Néanmoins, la correction se fait toujours à un instant donné de la fenêtre (typiquement au début ou au milieu). L'équation (045} peut alors s'écrire

(048a)

\[ J^b(\mathbf{x}) = \frac{1}{2}(\mathbf{x}-\mathbf{x}^b)^T\mathbf{B}^{-1}(\mathbf{x}-\mathbf{x}^b)\],

(048b)

\[ J^o(\mathbf{x}) = \frac{1}{2}\sum_{i=0}^N(\mathbf{y}^o_i-H\mathbf{x}(t_i))^TR^{-1}(\mathbf{y}^o_i-H\mathbf{x}(t_i))\],

où \[\mathbf{y}_i\] représente les \[N\] observations de la fenêtre d'assimilation aux instants \[t_i\] et \[\mathbf{x}(t_i)\] représente l'état du modèle aux instants \[t_i\].

Fig. 3DVar2 : Répartition des observations dans un 3D-Var FGAT.

L'avantage du 3D-Var FGAT sur le 3D-Var classique est qu'il prend réellement en compte l'instant de mesure des observations. De plus, il est très similaire au 4D-Var présenté dans la suite.

3D-Var incrémental

Il est possible d'écrire le 3D-Var de manière incrémentale. L'intérêt de cette formulation sera discutée dans la section 4DVar où il apparaîtra clairement. Cette formulation consiste à considérer comme contrôle l'écart entre l'état du système et l'ébauche, et non plus l'état du système. L'introduction de l'approche incrémentale en météorologie a été motivée par la réduction de coût qu'elle propose en 4D-Var. Soit

\[\delta\mathbf{x} = \mathbf{x} - \mathbf{x}^b \].

La fonction coût devient alors

(049)

\[ J(\delta\mathbf{x}) = \frac{1}{2}\delta\mathbf{x}^T\mathbf{B}^{-1}\delta\mathbf{x} +\frac{1}{2}(\mathbf{d}-\mathbf{H}\delta\mathbf{x})^T\mathbf{R}^{-1}(\mathbf{d}-\mathbf{H}\delta\mathbf{x})\],

où \[\mathbf{d}\] est le vecteur d'innovation \[\mathbf{y}-H\mathbf{x}^b\].

La formulation 3D-Var FGAT (Eqs. (048a} et (048b}) prend alors la forme suivante :

\[ J^b(\delta\mathbf{x}) = \frac{1}{2}\delta\mathbf{x}^T\mathbf{B}^{-1}\delta\mathbf{x}\],

(050)

\[ J^o(\delta\mathbf{x}) = \frac{1}{2}\sum_{i=0}^N(\mathbf{d}_i-\mathbf{H}\delta\mathbf{x})^T R^{-1}(\mathbf{d}_i-\mathbf{H}\delta\mathbf{x})\],

avec \[\mathbf{d}_i=\mathbf{y}^o_i - H\mathbf{x}(t_i)\].

Une illustration schématique (Fig. 3DVar3) du 3D-Var FGAT incrémental permet de comprendre simplement son utilisation. Il apparaît ainsi clairement que chaque analyse introduit un "choc" dans le modèle au moment de son introduction. Pour éviter cet inconvénient, il est possible d'utiliser une technique nommée IAU (Incremental Analysis Updating) et introduite par Bloom etal. (1996) qui consiste à répartir l'incrément d'analyse sur toute le cycle d'assimilation comme s'il s'agissait d'un forçage (Fig. 3DVar4).

Fig. 3DVar3 : Illustration de la procédure pour cycler le 3D-Var FGAT incrémental. Pour chaque cycle \[c\], le modèle d'évolution est intégré de \[t_0\] à \[t_N\] à partir de l'état initial d'ébauche \[\mathbf{x}^b_c(t_0)\] (courbe noir pleine) et le vecteur d'innovation \[\mathbf{d}_i\] est calculé pour les différentes observations \[\mathbf{y}^o_i\] avec \[i=1,\cdots,N\] (ligne fine verticale). L'analyse est effectuée à l'instant \[t\] compris entre \[t_0\] et \[t_N\]. Après l'analyse, un incrément est obtenu et est rajouté à l'état d'ébauche \[\mathbf{x}^b_c(t)\] permettant d'obtenir l'état analysé \[\mathbf{x}^a_c(t)\]. L'état analysé est ensuite propagé jusqu'à le fin du cycle \[c\] (courbe grise pointillée). Cet état analysé \[\mathbf{x}^a_c(t_N)\] est ensuite utilisé comme état initial d'ébauche pour le cycle suivant.

Illustration de la procédure pour cycler le 3D-Var FGAT incrémental avec IAU. Pour chaque cycle \[c\], le modèle d'évolution est intégré de \[t_0\] à \[t_N\] à partir de l'état initial d'ébauche \[\mathbf{x}^b_c(t_0)\] (courbe noir pleine) et le vecteur d'innovation \[\mathbf{d}_i\] est calculé pour les différentes observations \[\mathbf{y}^o_i\] avec \[i=1,\cdots,N\] (ligne fine verticale). L'analyse est effectuée à l'instant \[t\] compris entre \[t_0\] et \[t_N\]. Après l'analyse, un incrément est obtenu qui est rajouté comme un forçage lors de l'intégration du modèle d'évolution de \[t_0\] à \[t_N\] à partir de l'état initial d'ébauche \[\mathbf{x}^b_c(t_0)\] (courbe grise pointillée). Cet état analysé \[\mathbf{x}^a_c(t_N)\] est ensuite utilisé comme état initial d'ébauche pour le cycle suivant.

Méthode d'assimilation - Méthodes variationnelles

jOas — Thu, 15 Nov 2007 23:43:00 +0100

Introduites au début des années cinquante par Sasaki (Sasaki, 1955 et Sasaki, 1958), les méthodes variationnelles sont devenues pendant les années 1990 très populaires. De grands centres de prévisions météorologiques, tels que le NMC (U.S. National Meteorological Center, maintenant appelé National Centers for Environmental Prediction) en 1991 (Parrish et Derber, 1992), le CEPMMT (Centre Européen pour les Prévisions Météorologiques à Moyen Terme, aussi appelé ECMWF) en 1996 (Courtier etal., 1998 et Anderson etal., 1998) ou Météo-France en 1997, ont adopté ce type de méthode.

Cette approche de l'assimilation de données n'est plus basée sur des théories statistiques, mais sur la théorie de l'optimisation. En opposition aux méthodes séquentielles qui ne traitent les observations qu'au fur et à mesure de leur disponibilité sans jamais utiliser des observations futures, l'approche variationnelle traite le problème globalement sous la forme de la minimisation d'une fonctionnelle (fonction objective) mesurant les caractéristiques indésirables de la solution du modèle. Ces caractéristiques peuvent être l'écart aux observations, la présence d'onde de gravité, le non respect de certains équilibres, ou d'autres. Dans la suite, seul l'écart aux observations et l'éloignement à l'ébauche de la condition initiale seront pris en compte.

Si les statistiques ne sont plus les bases de ces méthodes, elles restent indispensables pour les calibrer et définir la fonction à minimiser.

L'approche variationnelle a déjà été abordée dans la section sur le BLUE en mettant en évidence, entre autre, l'équivalence à l'optimalité avec le BLUE.

Méthode d'assimilation - Filtre de Kalman d'ensemble (EnKF)

jOas — Thu, 15 Nov 2007 23:16:00 +0100

Le filtre de Kalman d'ensemble a été proposé par Evensen en 1994, puis corrigé en 1998. Pour une description détaillée, il est possible de se référer à Evensen (2003). Cette méthode a d'abord été présentée comme une alternative stochastique au filtre de Kalman étendu qui est déterministe. L'utilisation d'une méthode de Monte Carlo a été imaginée pour résoudre les deux principaux problèmes du filtre de Kalman étendu dans le cadre de système de grande taille non linéaire qui sont son coût très important et sa mauvaise réponse en cas de forte non-linéarité.

Le filtre de Kalman d'ensemble est très populaire car il est conceptuellement très simple et sa mise en oeuvre est aisée. En effet, il ne nécessite ni dérivation des opérateurs tangent-linéaires et des équations adjointes, ni intégration rétrograde du modèle d'évolution.

Le filtre de Kalman d'ensemble reste un filtre Gaussien et n'est pas un filtre particulaire malgré l'emprunt de la nation de particule (comme le filtre SEIK) car il ne gère les statistiques d'erreur que jusqu'à l'ordre deux. Au lieu de propager une matrice de covariance, les erreurs sont représentées statistiquement par un nuage de points propagés par le modèle d'évolution, sans aucune linéarisation. L'étape d'analyse est ensuite celle d'un filtre de Kalman standard.

Comme la montré Burgers etal. (1998),il est essentiel de perturber les observations pour chacun des membres de l'ensemble avec l'estimation de la matrice de covariance d'erreur d'observation \[\mathbf{R}\]. En effet, comme un échantillon statistique a tendance à s'appauvrir par coalescence des points, l'ajout de bruit dans les observations peut être interprété comme l'adjonction d'une partie stochastique permettant d'enrichir l'échantillon.

L'algorithme du filtre de Kalman d'ensemble peut être décrit de la manière suivante (cf. Fig. 1). A partir d'un ensemble de conditions initiales, un ensemble d'états d'ébauche à l'instant \[t_i\] est construit par de courtes intégrations du modèle d'évolution. La matrice de covariance d'erreur de prévision \[\mathbf{P}^f_i\] est calculée à partir de cet échantillon de telle manière que \[\mathbf{P}^f=E\left[(\mathbf{x}^f-\overline{\mathbf{x}}^f)(\mathbf{x}^f-\overline{\mathbf{x}}^f)^T\right]\]. La matrice de gain \[\mathbf{K}\] peut alors être calculée. Chaque ébauche est utilisée pour effectuer une analyse à l'instant \[t_i\] comme décrit par l'Eq. (038}. Les analyses sont obtenues avec des données bruitées. L'ensemble de ces états analysés est ensuite propagé jusqu'à l'instant \[t_{i+1}\] et permet alors d'estimer la matrice \[\mathbf{P}^f_{i+1}\]. Le rang des matrices ainsi estimées est inférieur ou égal à la taille de l'échantillon stochastique, c'est-à-dire très largement inférieur à la taille du vecteur d'état. Cette déficience de rang signifie que l'utilisation directe de la matrice \[\mathbf{P}^f\] dans l'algorithme d'assimilation contraint les corrections identifiées par l'analyse à être définies dans l'espace des membres de l'échantillon. Afin de palier à ce problème, le vecteur d'état \[\mathbf{x}\] peut être séparé en un vecteur \[\mathbf{x}^p\] projeté sur le sous-espace constitué par les échantillons et un vecteur orthogonal. La matrice estimée \[\mathbf{P}^f\] est alors liée dans l'algorithme à \[\mathbf{x}^p\] et une matrice de covariance d'erreur statique est liée au vecteur orthogonal à \[\mathbf{x}^p\]. Outre la conséquente économie de calcul et de stockage par rapport au filtre de Kalman étendu, l'algorithme du filtre de Kalman d'ensemble présente l'avantage d'être particulièrement adapté aux machines de calcul parallèle puisque chaque membre de l'échantillon d'analyse peut être calculé indépendamment des autres, et donc simultanément.

\input{EnKF.pst}
Fig. 1 : Représentation schématique des différentes étapes du filtre filtre de Kalman d'ensemble lors d'un cycle d'assimilation du temps \[t_i\] au temps \[t_{i+1}\]. L'indice \[k\] variant de \[1\] à \[N\] représente les différents membres de l'ensemble.

Méthode d'assimilation - Coût de Calcul et filtres dégradés

jOas — Thu, 15 Nov 2007 22:07:00 +0100

Les différents filtres de Kalman à rang réduit représentent des approches réalistes permettant l'implémentation du filtre de Kalman à des problèmes complexes et de grandes tailles.

En effet, alors que l'utilisation des filtres de Kalman ou de Kalman étendu nécessitent des ressources informatiques hors de portée pour des problèmes de grandes tailles comme l'océanographie ou la météorologie, le passage à un sous-espace représentatif de taille beaucoup plus petite permet la mise en oeuvre réaliste des méthodes utilisant cette technique.

L'intérêt pour ces différentes formes d'assimilation de données séquentielles varie en fonction des besoins et capacités. Autant l'interpolation optimale est simple à mettre en oeuvre et d'un coût extrêmement réduit, autant les filtres SEEK et SEIK représentent déjà un investissement important dans le cadre de problèmes de grandes tailles. Le filtre SEIK peut être considéré comme une variante du SEEK pour lequel la linéarisation est remplacée par une interpolation linéaire. Il résout partiellement le problème de la non-linéarité du filtre SEEK. Outre sa plus grande stabilité, il est plus simple à implémenter car il ne nécessite par le calcul du gradient du modèle d'évolution. Néanmoins, son coût est quasiment équivalent à celui du filtre SEEK car l'évolution de leurs bases de correction respectives requiert l'intégration du modèle \[r+1\] fois afin de faire évoluer les états interpolés du premier et les vecteurs de la base du deuxième. Ces deux filtres restent donc chers pour des problèmes de grandes tailles comme l'océanographie. En effet, le coût de ces deux filtres est supérieur à \[r+1\] fois le coût du modèle d'évolution.

Il est néanmoins possible de réduire ces coûts en simplifiant l'évolution de leur base de correction. Ce qui est évidement le seul moyen permettant de réduire les coûts.

Brasseur etal. (1999) ont été les premiers à suivre cette idée en proposant de ne garder que la base de correction initiale du filtre SEEK, calculée par une analyse EOFs, fixe dans le temps. Néanmoins, l'expérience montre que le filtre ainsi obtenu présente des instabilités lors du passage du modèle dans une période instable. En effet, les EOFs formant la base de correction ne captent, le plus souvent, que la variabilité globale de l'état du système. Ainsi, comme ce filtre ne fait pas évoluer sa base de correction pour suivre la dynamique du modèle, sa correction ne sera plus efficace lors des instabilités du modèle. Ce filtre se nomme Filtre de Kalman Singulier à base Fixe (SFEK)(Filtre de Kalman Singulier à base fixe) et assume donc l'hypothèse que l'opérateur du modèle d'évolution peut être considérer comme l'identité (\[M=\mathbf{I}\]).

Un autre moyen efficace de réduire les coûts des filtres SEEK et SEIK est d'utiliser la propriété de convergence du filtre de Kalman vers un régime permanent. Ainsi, après un certains temps, la base de correction des filtres SEEK et SEIK reste constante. Il est alors possible de construire un nouveau filtre de Kalman étendu (ou interpolé) singulier à base asymptotique noté SAEK (Filtre de Kalman étendu singulier à base asymptotique) (ou SAIK (Filtre de Kalman interpolé singulier à base asymptotique)) qui opère exactement comme le SEEK (ou le SEIK) durant les premiers temps du filtrage et comme le SFEK ensuite. Néanmoins, la propriété de convergence n'est pas valable dans le cas de modèles non-linéaires comme ceux en océanographie. Ce n'est pas non plus une condition suffisante et seule l'expérience et la comparaison des filtres SEEK et SAEK (ou SEIK et SAIK) peuvent permettre de vérifier la valeur de ces filtres. Enfin, en ce qui concerne le coût, ces filtres sont nettement moins coûteux que des filtres SEEK ou SEIK.

Dans le cas de problèmes faiblement non-linéaires, il est imaginable que la matrice de covariance d'erreur d'analyse tende rapidement vers un régime "semi-fixe" dans lequel elle évolue lentement. Dans ce cas, la base de correction va aussi évoluer lentement. Il alors possible de créer une variante du filtre SEEK (ou SEIK) nommé SIEEK (Filtre de Kalman étendu singulier à base évolutive par intermittence) (ou SIEIK (Filtre de Kalman interpolé singulier à base évolutive par intermittence) qui est capable de faire évoluer cette base selon deux régimes après une période d'initialisation avec le filtre SEEK (ou SEIK). Un régime fixe pour lequel la méthode est la même que pour le SFEK et un régime de rattrapage pour lequel la méthode redevient un SEEK (ou SEIK). Ces nouveaux filtres sont, eux aussi, beaucoup moins coûteux.

Il est aussi possible de faire évoluer les vecteurs de la base de correction sur une grille plus grossière que celle du modèle. Cette méthode réduit les coûts de calcul mais aussi les besoins de stockage.

D'une autre manière, il est aussi possible de ne faire évoluer qu'une partie de la base de correction et de garder l'autre partie fixe dans le temps. Ce filtre se nomme SSEEK (Filtre de Kalman étendu singulier à base semi-évolutive) (ou SSEIK (Filtre de Kalman interpolé singulier à base semi-évolutive)) dans le cas d'une évolution du SEEK (ou du SEIK). Son coût est, bien sûr, dépendant du nombre de vecteurs de base qui évoluent.

Dans tous ces exemples, lorsque le modèle est stable, les formes dégradées du SEEK et SEIK se comportent aussi bien que le SEEK et SEIK mais peuvent être de deux à dix fois plus rapides. Lorsque le modèle est en période instable, les performances de ces filtres se dégradent notablement.

Enfin, il est aussi possible travailler sur les EOFs en utilisant une base d'EOFs locales ou mixte. La thèse de Ibrahim Hoteit (2001) permet d'avoir une vision précise d'un grand nombre de ces variantes.

Méthode d'assimilation - Filtre SEIK

jOas — Thu, 15 Nov 2007 21:57:00 +0100

Le filtre de Kalman étendu peut présenter des instabilités en présence de fortes non-linéarités jusqu'à, parfois, diverger complètement (Evensen, 1992 ; Gauthier etal., 1994 et Kushner, 1967). Une possibilité pour tenter de résoudre cette difficulté est de remplacer la linéarisation dans le filtre de Kalman étendu par un développement de Taylor d'ordre supérieur. Malheureusement, cette approche n'est pas envisageable sur des systèmes de grandes dimensions comme l'océanographie. Une autre approche est possible en utilisant des méthodes stochastiques de type Monte Carlo pour estimer l'évolution de la matrice de covariance d'erreur par un nuage d'états centrés autour de l'état courant et donc la matrice de covariance empirique est celle de la matrice considérée. Cette approche, introduite par Evensen en 1994 avec son filtre de Kalman d'ensemble, est un très bon moyen pour traiter les modèles d'évolution fortement non-linéaires. Cette méthode sera présentée dans la section sur le filtre de Kalman d'ensemble. Néanmoins, cette méthode est limitée par la taille de l'échantillon à considérer.

En 2001, Pham etal. ont proposé une variante du filtre de SEEK, appelé filtre de Kalman Singulier Évolutif Interpolé (SEIK), dans lequel la taille de l'échantillon est, en un certain sens, minimale. En effet, il substitue à la linéarisation opérée dans le filtre de Kalman étendu et dans le SEEK une interpolation sur un échantillon d'états bien choisis propagés dans l'étape de prévision. L'idée du SEIK est donc de faire évoluer la matrice de covariance d'erreur à l'aide d'un nuage de points de taille raisonnable. Dans ce but, Pham a émit l'hypothèse de rang faible \[r\] de la matrice de covariance d'erreur pour réduire la taille du nuage de points à \[r+1\] points exactement. L'autre originalité de ce filtre réside dans le choix des états d'interpolation qui sont tirés "au hasard" à chaque pas de filtrage afin de ne pas privilégier une direction particulière de l'espace d'état. La Fig. 1 permet de mettre en évidence les différentes étapes nécessaire au filtre SEIK.

Fig.1 : Représentation schématique des différentes étapes du filtre SEIK lors d'un cycle d'assimilation du temps \[t_i\] au temps \[t_{i+1}\]. L'indice \[k\] variant de \[1\] à \[r+1\] représente les différents membres du nuage de points.

Méthode d'assimilation - Filtre SEEK

jOas — Thu, 15 Nov 2007 21:45:00 +0100

Le filtre SEEK (Singular Evolutive Extended Kalman filter) a été introduit par Pham etal. en 1998. Il s'agit d'un filtre réduit déduit du filtre de Kalman étendu. Il repose sur la stagnation ou la décroissance du rang des matrices de covariances d'erreur, une propriété avérée ou forcée selon les cas.

Dans le cas d'un filtre sans erreur modèle, il résulte des Eqs. (039) et (041) que le rang de \[\mathbf{P}^f_i\], noté \[r=\mathrm{rang}(\mathbf{P}^f_i)\], est une fonction décroissante de \[t_i\], puisque, \emph{in fine}, la récurrence est de la forme \[\mathbf{P}^f_{i+1}=\mathbf{A}_i\mathbf{P}^f_i\mathbf{B}_i\] . En conséquence, si le rang de la matrice de covariance d'erreur initiale est faible comparée à la dimension \[n\] de l'espace du modèle, il le restera. Il est alors possible de décomposer la matrice \[\mathbf{P}^a_i\] avec une matrice diagonale à coefficient positifs ou nuls et une matrice orthogonale décrivant les \[r\] directions principales d'erreur. A partir de cette décomposition, il est possible de poser le problème dans l'espace des directions principales d'erreur de taille très inférieure à la dimension du système original. L'analyse est alors effectuée dans cet espace réduit et a pour caractéristique de ne pas modifier l'espace engendré par les directions principales des erreurs. Ce qui n'empêche pas, en général, à ces directions de changer.

Dans un cadre plus général, le modèle n'est pas parfait. Il n'est pas possible de négliger \[\mathbf{Q}_i\] . Ainsi, il apparaît que le rang de \[\mathbf{P}^f_{i+1}\] peut être supérieur à celui de \[\mathbf{P}^a_i\] . De plus, il ne peut plus y avoir de réduction du filtre sans approximation. La solution la plus naturelle et la plus simple pour réduire le rang consiste à projeter \[\mathbf{Q}_i\] avec une projection orthogonale. C'est-à-dire que seule la composante de l'erreur modèle agissant dans le sous-espace sur lequel agit \[\mathbf{P}^a_i\] est retenue.

L'idée de filtre SEEK peut donc être résumé à ceci : le système d'évolution \[M\] amplifie les erreurs associées à un sous-espace \[\mathcal{A}_i\] de l'espace tangent à l'espace modèle, tandis que les erreurs associées au complémentaire de cet espace sont atténuées. Si le rang de la matrice de covariance d'erreur est supérieur à la dimension de \[\mathcal{A}_i\], alors il est possible d'espérer que le système d'évolution se chargera d'atténuer toutes les erreurs commises et non corrigées dans l'espace complémentaire de \[\mathcal{A}_i\] .

Méthode d'assimilation - Filtre RRSQRT

jOas — Thu, 15 Nov 2007 21:36:00 +0100

Le filtre RRSQRT est une réponde à ce problème. Il permet d'éviter les différentes difficultés d'implémentation mise en évidence auparavant en représentant les directions principales des matrices d'erreur par des modes réduits. Ainsi, il possible d'utiliser exclusivement les modes au détriment des matrices.

En reprenant l'espace du modèle de dimension \[n\], avec un état du système initial \[\mathbf{x}_0^f\] auquel est associé la matrice de covariance d'erreur \[\mathbf{P}^f_0\], il faut réaliser une décomposition en mode principaux de cette matrice telle que

(041)

\[ \mathbf{P}^f_0 & \simeq & \mathbf{S}^f_0(\mathbf{S}^f_0)^T\],

où \[\mathbf{S}^f_0\] est une matrice de dimension \[(m \times n)\] avec \[m\] représentant les \[m\] premiers modes principaux de \[\mathbf{P}^f_0\] . L'erreur sur l'ébauche a donc été réduite. Il est alors possible de définir un opérateur d'observation \[ {\Psi}=(H_i \mathbf{S}^f_i)^T\] . Le gain de Kalman, calculé dans l'espace d'analyse, peut alors être décrit en fonction de \[{\Psi}^T\] de taille \[(p\times m)\]

(042)

\[ \mathbf{K}^*_i = \mathbf{S}_i {\Psi}\left({\Psi}^T{\Psi} +\mathbf{R}_i\right)^{-1}\].

Et, il est aussi possible d'obtenir la matrice racine de covariance d'analyse sans faire de calcul directement avec les matrices de covariance d'erreur :

(043)

\[ \mathbf{S}^a_i = \mathbf{S}^f_i\left(\mathbf{I}-{\Psi}\left({\Psi}^T{\Psi} +\mathbf{R}_i\right)^{-1}{\Psi}^T\right)^{1/2}\].

Le calcul de la racine de \[\mathbf{S}^a_i\] pourrait être coûteux, mais il n'en est rien puisque \[\mathbf{I}-{\Psi}\left({\Psi}^T{\Psi} +\mathbf{R}_i\right)^{-1}{\Psi}^T\] est de taille \[(m\times m)\]. De plus, la matrice racine est mieux conditionnée, ce qui assure une meilleure précision numérique.

Après l'analyse, la dimension du système est réduit en passant de \[m\] modes à \[m-q\] modes. Pour cela, il suffit de diagonaliser \[(\mathbf{S}^a_i)^T\mathbf{S}^a_i\] et de ne retenir que les \[m-q\] plus grandes valeurs propres de la matrice de passage, puis de réduire \[\mathbf{S}^a_i\] .

A l'étape de prévision, l'état analysé est propagé par le modèle d'évolution et la racine de la matrice réduite \[\tilde{\mathbf{S}}^a_i\] par le modèle linéaire-tangent. Elle est ensuite élargie en ajoutant \[q\] modes imputés à l'erreur modèle :

(044)

\[ \mathbf{S}^f_{i+1} = [\mathbf{M}_{i \to i+1}\tilde{\mathbf{S}}^a_i(\mathbf{M}_{i \to i+1})^T,\mathbf{T}_i]\].

Cette matrice comporte alors \[m\] modes.

Méthode d'assimilation - Les filtre de Kalman réduits

jOas — Thu, 15 Nov 2007 21:25:00 +0100

Depuis R. E. Kalman, les filtres ont été utilisés dans de nombreuses applications. Mais très vite, les aspects limitants de l'implémentation du filtre de Kalman sont apparus. Ainsi, l'assimilation de données n'était pas possible dans des domaines comme la météorologie, ou plus tard, l'océanographie car les dimensions du problème rendaient excessif le coût numérique et, de plus, les statistiques nécessaires au filtre de Kalman ne sont que rarement connues.

Pour résoudre ce problème, une hypothèse peut permettre de le contourner. L'idée est, qu'à un instant donné, la physique du modèle est contrôlée par un nombre ou une combinaison limitée de variables. L'hypothèse est alors que les statistiques d'erreurs significatives sont données par celles portant sur ces variables contrôlant la physique du modèle (les modes réduits). Il est alors nécessaire des les identifier. De plus, il faut aussi être capable d'enrichir stochastiquement le système afin que la base de modes réduits puisse évoluer sans contraintes trop fortes. En effet, le risque est que ces modes, s'ils dégénèrent, ne sous-tendent plus la fraction de l'espace des états dans lequel évolue le système.

Méthode d'assimilation - L'exemple du naufragé

jOas — Thu, 15 Nov 2007 21:11:00 +0100

Revenons aux mésaventures de notre naufragé introduites précédemment. Finalement, ne sachant comment atteindre le rivage, il se résout à évaluer la distance le séparant du rivage toutes les heures. Il dispose ainsi de \[i\] mesures de la distance du canot au rivage (\[v^o_i\]) entre l'instant de son naufrage \[t_0\] et la dernière mesure au temps \[t_i\]. Cette évaluation est supposée sans biais et sa variance, notée comme précédemment \[s^o\], est supposée stationnaire. Les coordonnées réelles du canot sont \[(u_i,v_i)\],tandis que celle issues de l'analyse \[(u_i^a,v_i^a)\] et celles de la prévision \[(u_i^f,v_i^f)\]. A l'instant du naufrage (\[t_0\]), la position du canot est \[(u_0^a,v_0^a)=(0,0)\]. Entre deux mesures aux instants \[t_i\] et \[t_{i+1}\],le canot dérive mais sa direction n'est pas connue. Le naufragé imagine donc un modèle d'évolution comme un modèle de diffusion autour de son point d'origine. Il peut donc écrire le modèle (linéaire en l'occurrence) tel que \[\mathbf{M}_{i \to i+1}=\mathbf{I}\]) et l'erreur modèle, qu'il suppose importante, telle que

\[ \mathbf{Q}_i=\left( \begin{array}{cc} s^m & 0 \\ 0 & s^m \end{array} \right)\],

où \[s^m\] est proportionnel au temps écoulé entre \[t_{k+1}\] et \[t_0\].

Le naufragé fait ensuite la supposition qu'il n'y a aucune corrélation entre les deux coordonnées de la position du canot de sauvetage. Il écrit alors les matrices de covariances d'erreur d'analyse et de prévision sous la forme :

\[\mathbf{P}^a_i=\left( \begin{array}{cc} \lambda_i & 0 \\ 0 & \mu_i \end{array}, \right) \quad \textrm{et} \quad \mathbf{P}^f_i=\left( \begin{array}{cc} \nu_i & 0 \\ 0 & \rho_i \end{array}, \right)\].

Analyse

En raisonnant à partir de l'instant \[t_i\], le naufragé dispose du vecteur résultant de sa dernière analyse \[(u_i^f,v_i^f)^T\]. Il effectue alors une nouvelle mesure \[v_i^o\] et effectue une analyse optimale en utilisant le gain de Kalman (Eq. (032)) :

\[ \mathbf{K}^*_i=\left( \begin{array}{cc}\nu_i & 0 \\ 0 & \rho_i \end{array} \right)\left( \begin{array}{c} 0 \\ 1 \end{array} \right)\left(s^o+ \left(0,1\right) \left( \begin{array}{cc} \nu_i & 0 \\ 0 & \rho_i \end{array} \right) \left( \begin{array}{c} 0 \\ 1 \end{array} \right) \right)^{-1} \]

\[\mathbf{K}^*_i=\frac{\rho_i}{s^o+\rho_i} \left( \begin{array}{c} 0 \\ 1 \end{array} \right)\].

Il peut alors en déduire les coordonnées de la position du canot (Eq. (033))

\[ \left( \begin{array}{c} u_i^a \\ v_i^a \end{array} \right) = \left( \begin{array}{c} u_i^f \\ v_i^f \end{array} \right) +\mathbf{K}^*_i\left(v_i^o-\left(0,1\right) \left( \begin{array}{c} u_i^f \\ v_i^f \end{array} \right) \right)\]

\[ \left( \begin{array}{c} u_i^a \\ v_i^a \end{array} \right)=\left( \begin{array}{c} u_i^f \\ v_i^f \end{array} \right) +\mathbf{K}^*_i\left(v_i^o-v_i^f\right)\]

\[ \left( \begin{array}{c} u_i^a \\ v_i^a \end{array} \right)=\left( \begin{array}{c} u_i^f \\ v_i^f+\frac{\rho_i}{s^o+\rho_i}\left(v_i^o-v_i^f\right) \end{array} \right)\],

et la matrice de covariance d'erreur d'analyse (Eq. (034))

\[ \left( \begin{array}{cc} \lambda_i & 0 \\ 0 & \mu_i \end{array} \right)=\left(\left( \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right)-\frac{\rho_i}{s^o+\rho_i}\left( \begin{array}{c} 0 \\ 1 \end{array} \right)\left(0,1\right)\right)\left( \begin{array}{cc} \nu_i & 0 \\ 0 & \rho_i\end{array} \right)\]

\[ \left( \begin{array}{cc} \lambda_i & 0 \\ 0 & \mu_i \end{array} \right)=\left( \begin{array}{cc} \nu_i & 0 \\ 0 & \frac{s^o \rho_i}{s^o+\rho_i} \end{array} \right)\].

Il en déduit alors que

\[\lambda_i=\nu_i \quad \textrm{et} \quad \frac{1}{\mu_i}=\frac{1}{s^o}+\frac{1}{\rho_i} \].

Prévision

Le naufragé passe maintenant à l'étape de prévision du filtre de Kalman et applique son modèle d'évolution \[\mathbf{I}\] (Eq. (035)) pour estimer sa prochaine position

\[ \left( \begin{array}{c} u^f_{i+1} \\ v^f_{i+1} \end{array} \right)=\left( \begin{array}{c} u^a_{i} \\ v^a_{i} \end{array} \right)\].

Ensuite, il calcule sa prochaine matrice de covariance d'erreur de prévision (Eq. (036)) facilement et trouve que \[\mathbf{P}^f_{i+1}=\mathbf{P}^a_i+s^m\mathbf{I}\], ce qu'il écrit aussi sous la forme

\[ \nu_{i+1}=\lambda_i+s^m \quad \textrm{et} \quad \rho_{i+1}=\mu_i+s^m \].

Il résulte alors de la succession de l'analyse et de la prévision que

\[ \left\{ \begin{array}{l} u^f_{i+1} = u^f_i, \\ v^f_{i+1} = v^f_i + \frac{\rho_i}{s^o+\rho_i}\left(v^o_i-v^f_i\right),\end{array} \right.\]

\[ \left\{ \begin{array}{l} \nu_{i+1} = \nu_i + s^m, \\ \frac{1}{\rho_{i+1}-s^m} = \frac{1}{s^o}+\frac{1}{\rho_i}. \end{array} \right.\]

Comme le naufragé sait que \[u_0^f=0\],il en déduit que \[u_k^f=0\]. Sans mesure, il n'apprend rien sur sa position le long de la côte. Cependant, l'incertitude croît linéairement avec le temps puisque $\nu_k=k s^m\].

De l'incertitude sur la coordonnée \[v\],le naufragé en déduit que, dans un premier temps, la confiance sur l'analyse est la somme de la confiance sur l'observation et sur la prévision, et, dans un second temps, que l'erreur de la prévision est la somme de l'erreur modèle et d'analyse. C'est donc le résultat d'un compromis entre la réduction de l'incertitude par l'assimilation et son accroissement due à la dérive non contrôlée de l'embarcation.

Méthode d'assimilation - Le coût de calcul

jOas — Thu, 15 Nov 2007 20:52:00 +0100

L'algorithme du filtre de Kalman complète le système d'équations lié à la détermination de l'état analysé et à sa propagation dans le temps avec deux équations de calcul et de propagation de la matrice de covariance d'erreur d'analyse. Le coût numérique du filtre de Kalman est donc la somme du coût du traitement du vecteur d'état et des covariances d'erreur. Pour les systèmes de grande taille tels que l'océan ou l'atmosphère, le coût de calcul principal provient du traitement des covariances d'erreur d'analyse. La première étape coûteuse est l'inversion de la matrice \[\left(\mathbf{H}_i\mathbf{P}^f_i\mathbf{H}_i^T+\mathbf{R}_i \right)\]. La propagation par les équations de la dynamique du modèle linéaire-tangent de \[\mathbf{P}^a\] requiert ensuite la multiplication par la matrice \[\mathbf{M}\] par chaque colonne (chaque ligne pour \[ \mathbf{M}^T \]) de \[ \mathbf{P}^a \] (autour de \[ 10^7 \times 10^7 \] opérations). Au delà du coup de calcul exorbitant de ces opérations, il est impossible de stocker entre chaque étape d'analyse de telles matrices malgré les capacités déjà importantes disponibles. Pour ces raisons, l'algorithme du filtre de Kalman ne peut être appliqué qu'à des systèmes de taille réduite.

Il doit donc être simplifié pour permettre son application aux systèmes océaniques et atmosphériques. Plusieurs études visent notamment à réduire le nombre d'intégration du modèle linéaire-tangent en ne propageant pas la matrice de covariance d'erreur que suivant certaines directions (Fukomori etal., 1995 ; Evensen, 1994 ; Fisher, 1998 et Evensen, 2003). Il faut tout d'abord identifier un sous-espace de dimension réduite. Ensuite, seule la projection de la matrice de covariance dans ce sous-espace, et non la matrice complète, est propagée. Divers de ces filtres seront présentés dans la section suivante.

Méthode d'assimilation - Filtre de Kalman étendu (EKF)

jOas — Thu, 15 Nov 2007 20:43:00 +0100

Dans le filtre de Kalman classique, le modèle d'évolution et l'opérateur d'observation sont supposés linéaires. Cependant, il arrive souvent que l'hypothèse de linéarité ne soit pas valide. Dans ce cas, il est possible de généraliser le filtre de Kalman en utilisant des formes linéarisées de l'opérateur d'observation et du modèle d'évolution pour les Eqs. (032), (034) et (036) et la forme non-linéaire pour les Eqs. (033) et (035). Ce filtre est appelé filtre de Kalman étendu (EKF).

Les cinq étapes de l'analyse peuvent alors s'écrire ainsi :

Calcul de la matrice de gain \[\mathbf{K}\] au temps \[t_i\] :

(037)
\[ \mathbf{K}_i = \mathbf{P}^f_i\mathbf{H}_i^T\left(\mathbf{H}_i\mathbf{P}^f_i\mathbf{H}_i^T+\mathbf{R}_i\right)^{-1} \].
Analyse au temps \[t_i\] à l'aide le l'opérateur d'observation non-linéaire :

(038)
\[ \mathbf{x}^a_i = \mathbf{x}^f_i +\mathbf{K}_i\left(\mathbf{y}^o_i - H_i \mathbf{x}^f_i \right)\].
Calcul de la matrice de covariance d'erreur d'analyse au temps \[t_i\] :

(039)
\[ \mathbf{P}^a_i = \left( \mathbf{I} - \mathbf{K}_i\mathbf{H}_i \right) \mathbf{P}^f_i\].
Prévision au temps \[t_{i+1}\] par propagation de l'analyse de \[t_i\] à $t_{i+1}\] par le modèle non-linéaire d'évolution :

(040)
\[ \mathbf{x}^f_{i+1} = M_{i \to i+1}(\mathbf{x}^a_i)\].
Calcul de la matrice de covariance d'erreur de prévision au temps \[t_{i+1}\] par propagation de la matrice de covariance d'erreur d'analyse de \[t_i\] à \[t_{i+1}\] par le modèle linéaire d'évolution :

(041)
\[ \mathbf{P}^f_{i+1} = \mathbf{M}_{i \to i+1} \mathbf{P}^a_i \mathbf{M}^T_{i \to i+1} + \mathbf{Q}_i\].

A noter que, bien que le filtre de Kalman soit une analyse optimale, le filtre de Kalman étendu perd cette qualité (il ne fournit pas la solution de variance minimale). Néanmoins, l'utilisation du filtre de Kalman étendu dans un cadre faiblement non-linéaire permet d'obtenir de bonnes analyses. De plus, la linéarisation du modèle d'évolution \[M\] peut interagir avec les erreurs modèle de manière assez compliquée.

Méthode d'assimilation - Le filtre de Kalman (KF)

jOas — Thu, 15 Nov 2007 17:01:00 +0100

Filtre de Kalman - KF

En 1960-61, Kalman et Bucy ont décrit une solution récursive pour des problèmes de filtrage linéaire de données discrète. Cette solution est depuis nommée filtre de Kalman. Ce filtre peut être appréhendé comme une extension du BLUE pour laquelle l'état analysé pour une étape donnée définit l'ébauche à l'étape d'analyse suivante. Outre ceci, le filtre de Kalman incorpore un modèle d'évolution de l'état du système entre deux instants \[t_i\] et \[t_{i+1\].

Pour décrire cette méthode, les notations usuellement utilisées seront reprises. En particulier, la matrice des covariances d'erreur d'analyse, jusqu'ici notée \[\mathbf{A}\], se nommera maintenant \[\mathbf{P}^a\] et celle d'ébauche, jusqu'ici notés \[\mathbf{B}\], s'appellera \[\mathbf{P}^f\] de manière à mettre en évidence que l'ébauche (b comme background) est maintenant une prévision (f comme forecast). Le modèle d'évolution non-linéaire sera noté \[M\], linéarisé il se nommera \[\mathbf{M}\] et son adjoint \[\mathbf{M}^T\]. L'utilisation d'un modèle d'évolution entre deux instant \[t_i\] et \[t_{i+1}\] entraîne un nouveau type d'erreur nommée l'erreur modèle. Elle est supposée non-biaisée et est décrite par la matrice de covariance d'erreur du modèle à chaque instant \[ t_i\] : \[\mathbf{Q}_i\]. De plus, les erreurs d'analyse et modèle sont supposées non-corrélées.

Le schéma d'assimilation peut être décrit de la manière suivante : à partir d'une prévision à l'instant \[t_i\] et de sa matrice de covariance d'erreur de prévision \[\mathbf{P}^f_i\], une analyse est effectuée permettant d'obtenir un état analysé et une matrice de covariances d'erreur d'analyse \[\mathbf{P}^a_i\] à l'instant \[t_i\]. Ensuite, une prévision du temps \[t_i\] à \[t_{i+1}\] est effectuée en partant de l'état analysé. De manière similaire, la matrice de covariance d'erreur d'analyse est propagée par le modèle d'évolution linéaire afin d'obtenir la matrice de covariance d'erreur de prévision \[\mathbf{P}^f_{i+1}\] à l'instant \[t_{i+1}\]. Il suffit ensuite de répéter cette opération.

La deuxième étape, durant laquelle l'état analysé et la matrice de covariances d'erreur d'analyse est propagée jusqu'au temps d'observation suivant, est clairement la plus coûteuse.

De manière plus formelle, l'algorithme du filtre de Kalman entre les instants d'observation \[t_i\] et \[t_{i+1}\] peut être décrit les Eqs. (032), (033), (034), (035) et (036).

Calcul de la matrice de gain \[\mathbf{K}\] au temps \[t_i\] :
(032)
\[\mathbf{K}_i = \mathbf{P}^f_i\mathbf{H}_i^T\left(\mathbf{H}_i\mathbf{P}^f_i\mathbf{H}_i^T+\mathbf{R}_i\right)^{-1}\].
Analyse au temps \[t_i\] :

(033)
\[\mathbf{x}^a_i = \mathbf{x}^f_i +\mathbf{K}_i\left(\mathbf{y}^o_i - \mathbf{H}_i \mathbf{x}^f_i \right)\].
Calcul de la matrice de covariance d'erreur d'analyse au temps \[t_i\] :

(034)
\[\mathbf{P}^a_i = \left( \mathbf{I} - \mathbf{K}_i\mathbf{H}_i \right) \mathbf{P}^f_i\].
Prévision au temps \[t_{i+1}\] par propagation de l'analyse de \[t_i\] à \[t_{i+1}\] par le modèle linéaire d'évolution :

(035)
\[ \mathbf{x}^f_{i+1} = \mathbf{M}_{i \to i+1}(\mathbf{x}^a_i)\].
Calcul de la matrice de covariance d'erreur de prévision au temps \[t_{i+1}\] par propagation de la matrice de covariance d'erreur d'analyse de \[t_i\] à \[t_{i+1}\] par le modèle linéaire d'évolution :

(036)
\[\mathbf{P}^f_{i+1} = \mathbf{M}_{i \to i+1} \mathbf{P}^a_i \mathbf{M}^T_{i \to i+1} + \mathbf{Q}_i\].