Méthodes variationnelles d’ensemble et optimisation variationnellepour les géosciences

L’assimilation de données est une discipline visant à estimer l’état d’un système physique à l’aide d’observations physiques partielles. En géophysique, ce problème devient l’estimation des paramètres du système terrestre à l’aide de mesures terrestres bruitées et parcellaires. Lorsque les observations deviennent disponibles à mesure que le temps passe, il s’agit d’estimer séquentiellement la trajectoire du système à l’aide d’un historique d’observations. Le contexte géophysique impose un grand nombre de paramètres et d’observations atteignant les limites des ressources numériques disponibles (Evensen, 2009).

L’assimilation de données sans dimension temporelle

Le système géophysique est une abstraction des processus physiques terrestres. Il n’est donc pas la réalité mais une simplification que l’on espère suffisamment précise pour être capable de fournir de bonnes prévisions. Pour cela, l’état du système géophysique est quantifié par un vecteur réel représentant température, pression, humidité, etc en différents points du globe pour un total de plusieurs millions (voir milliards) de paramètres. On supposera donc que notre système est paramétrisé par un vecteur x ∈ Rn avec n suffisamment grand pour ne pas pouvoir négliger les ressources numériques disponibles. De même, les observations physiques sont paramétrisées par un vecteur y ∈ Rd . Les équations régissant le système géophysique sont rassemblées dans l’opérateur d’observation H. L’observation y est donc sensée être déterminée par l’état x avec une faible erreur d’observation e :
y = H (x) + e.   (1.1)

C’est parce que le système n’est qu’une représentation simplifiée de la réalité dont les observations y sont une mesure physique, que cette relation entre état et observation est inexacte. Un objectif de l’assimilation de données consiste à déterminer l’état x à partir des observations y avec une faible erreur d’estimation n :
x = K (y) + n.    (1.2)

Il s’agit d’un problème inverse (Tarantola, 2004). Malheureusement, sans plus d’informations sur l’erreur d’observation e on ne peut déduire de y aucune information sur x. On sait cependant qualitativement que l’opérateur d’observation est efficace : les petites erreurs sont plus probables que les grandes. Il faut alors quantifier l’aléatoire du système pour pouvoir tenir compte de cette information de nature probabiliste.

l’assimilation de données sans dimension temporelle 

Nous avons vu dans cette section des techniques permettant d’analyser l’état x à l’aide de l’assimilation de l’observation y. Dans le cadre bayésien, la solution de cette analyse est la densité a posteriori p (x|y). La règle de Bayes permet d’exprimer cette densité a posteriori en fonction de la vraisemblance et de la densité a priori. Cependant, cette expression est de faible intérêt en pratique, car l’intégration nécessaire au calcul de ses moments n’est pas calculable. L’alternative proposée par l’échantillonnage d’importance est d’approcher ces intégrales par des formules de quadrature. Cela revient à approcher les densités par des sommes de Dirac. Bien qu’exacte dans la limite du nombre d’échantillons, cette méthode requiert un grand nombre d’évaluations de l’opérateur d’observation H. Puisque cet opérateur porte l’essentiel de la complexité numérique, la méthode par échantillonnage d’importance n’est pas immédiatement envisageable en grande dimension.

En revanche, l’opérateur d’observation H est affine si et seulement si la densité jointe p (x, y) est gaussienne. La moyenne et la matrice de covariance a posteriori caractérisent alors la densité a posteriori et peuvent être calculées algébriquement. La moyenne a posteriori est, dans ce cas, le meilleur estimateur de x fonction de y en erreur quadratique moyenne (MMSE), le meilleur estimateur de x fonction affine de y en erreur quadratique moyenne (BLUE), ainsi que la maximum a posteriori (MAP). Dans le cas général, on peut s’inspirer du cas linéaire pour faire une approximation gaussienne de p (x|y) ne retenant que ses deux premiers moments. La moyenne et la matrice de covariance a posteriori ne caractérisent plus la densité a posteriori, elles engendrent néanmoins sa meilleure approximation gaussienne au sens de l’entropie relative. En particulier, la moyenne a posteriori est le meilleur estimateur de x fonction de y au sens de l’erreur quadratique moyenne (MMSE). Le calcul de ces moments a posteriori demande cependant une intégration par rapport à p (x|y). Cette intégration est impossible parce que l’intégrande p (x|y) est justement la densité que l’on cherche à calculer. En revanche, si l’on sait intégrer (x, H (x)) par rapport à la densité a priori p (x), la moyenne et la matrice de covariance de (x, y) sont accessibles et fournissent la meilleure approximation gaussienne de la densité jointe p (x, y) au sens de l’entropie relative. Cette approximation gaussienne de p (x, y) correspond à une approximation affine de H minimisant l’erreur quadratique moyenne avec une approximation gaussienne de son erreur. On obtient alors une approximation de moments a posteriori dont l’approximation de la moyenne a posteriori est le meilleur estimateur de x fonction affine de y en erreur quadratique moyenne.

Finalement, le MAP de p (x|y) fournit un estimateur fonction non linéaire de y calculable grâce aux méthodes itératives d’optimisation. Cette technique ne fournit malheureusement pas naturellement d’estimateur de second ordre. Une alternative consiste alors à donner une portée statistique à l’approximation affine de H au sein de la méthode de Gauss-Newton. Cette approximation se doit d’être précise dans les zones de forte probabilité pour que le calcul des moments soit relativement bon. A leur tour, ces moments peuvent servir à mieux localiser la zone de forte probabilité et ainsi améliorer notre approximation affine de H.

Localisation et inflation 

En géophysique, l’état représente des paramètres du système en chaque point du globe terrestre. On observe alors fréquemment une faible corrélation entre des points éloignés. Cette connaissance a priori peut être utilisée pour corriger l’estimation des matrices de covariance et réduire l’erreur d’échantillonnage. Deux principales méthodes sont employées pour rendre locale la matrice de covariance.  première est la localisation en domaine et consiste à faire plusieurs analyses indépendantes pour chaque domaine du globe. Pour chacune de ces analyses, seules les observations à proximité du domaine sont utilisées. La seconde est la localisation par produit de Schur et consiste à multiplier terme à terme la matrice estimée avec une seconde matrice atténuant les corrélations à longue distance. Pour un traitement plus complet de la localisation on pourra consulter Houtekamer and Mitchell (2001), Hamill et al. (2001), Evensen (2003), Ott et al. (2004) ou Ménétrier et al. (2015) .

Pour estimer l’état x, on utilise une estimation de la moyenne et de la covariance a priori. L’estimation étant inexacte, il en résulte une incertitude sur ces moments. Cependant, l’algorithme n’en tient pas compte. Il sous-estime alors l’incertitude sur x. Puisque la variance a posteriori est une mesure de l’incertitude sur x, l’algorithme la sous-estime aussi. Lors du prochain cycle, l’incertitude sur les moments a priori sera alors d’autant plus importante. On observe donc une amplification de l’incertitude à chaque cycle qui mène nécessairement à une divergence de l’algorithme si elle n’est pas compensée. C’est le rôle de l’inflation. Elle consiste à dilater les matrices de covariances par un facteur multiplicatif afin de corriger l’erreur d’échantillonnage et de mieux estimer l’incertitude sur x.

Table des matières

1 Introduction
1.1 L’assimilation de données sans dimension temporelle
1.1.1 Modélisation des probabilités
1.1.2 Solution bayésienne
1.1.3 Méthode par échantillonnage d’importance
1.1.4 Spécifications gaussiennes a priori
1.1.5 Échantillonner une densité
1.1.6 Solution dans le cadre linéaire et gaussien
1.1.7 Méthode dans un cadre moins linéaire et gaussien
1.1.8 Méthodes variationnelles
1.1.9 Conclusion de l’assimilation de données sans dimension temporelle
1.2 L’assimilation de données avec dimension temporelle
1.2.1 Solution bayésienne
1.2.2 Le filtre particulaire
1.2.3 Solution dans le cadre linéaire : le filtre de Kalman
1.2.4 Solution dans le cadre linéaire : le lisseur de Kalman
1.2.4.1 Remarque sur la méthode aller-retour
1.2.5 Le filtre de Kalman d’ensemble stochastique
1.2.5.1 Localisation et inflation
1.2.6 Le filtre de Kalman d’ensemble déterministe
1.2.7 Méthode variationnelle : 4DVar
1.2.8 Conclusion de l’assimilation de données avec dimension temporelle
1.3 L’assimilation de données ensemble variationnelle fortement contrainte
1.3.1 Méthodes hybrides
1.3.2 Ensemble d’assimilations
1.3.3 Méthodes 4D variationnelles d’ensemble
1.3.4 Le lisseur de Kalman itératif
1.3.5 Conclusion de l’assimilation de données ensemble variationnelle fortement contrainte
2 Assimilation de données ensemble variationnelle quasi-statique
2.1 Introduction
2.1.1 Contexte
2.1.2 Assimilation de données variationnelle quasi-statique
2.1.3 Méthodes variationnelles d’ensemble
2.1.4 Plan
2.2 La fenêtre et les performances d’assimilation
2.2.1 Le 4DVar et l’IEnKS
2.2.2 Performance de l’assimilation
2.2.3 Performance dans le cas linéaire, diagonal et autonome
2.2.4 Performance dans le cas non linéaire et chaotique
2.2.4.1 Des minimums locaux multiples
2.2.4.2 Longueur effective de la fenêtre d’assimilation
2.3 Algorithmes quasi-statiques
2.4 Expériences numériques
2.5 Conclusion
3 Le lisseur de Kalman itératif faiblement contraint
3.1 Introduction
3.1.1 Méthodes 4D ensemble variationnelles faiblement contraintes
3.1.2 Plan
3.2 Assimilation de données séquentielle faiblement contrainte dans l’espace de l’état
3.2.1 Le système d’assimilation de données
3.2.2 Principes d’assimilation de données séquentielle
3.2.3 Comparaison avec le biais d’erreur modèle comme variable de contrôle
3.3 Construction de l’IEnKS-Q
3.3.1 Représentation dans le sous-espace de l’ensemble
3.3.2 Analyse
3.3.2.1 Analyse linéaire de la variable d’ensemble
3.3.2.2 Construction récursive des dérivées
3.3.2.3 Analyse variationnelle de la variable d’ensemble
3.3.2.4 Décomposition de la matrice de déviation et modèle de substitution
3.3.3 Propagation
3.3.3.1 Réduction de la dimension et marginalisation
3.4 Algorithmes
3.5 Résultats numériques
3.5.1 Validation de l’IEnKS-Q
3.5.2 Test du paramètre G
3.5.3 Test avec des modèles de substitution
3.6 Conclusion
4 Conclusion
References

Télécharger le rapport complet