Théorie – Datice

Inférence bayésienne

Lors du développement de Datice, Bénédicte Lemieux-Dudon s’est basée sur l’approche et le vocabulaire introduit par Mosegaard et Tarantola, 2002¹K. Mosegaard and A. Tarantola. International Handbook of Earthquake and Engineering Seismology (Part A). Chapter : Probabilistic Approach to Inverse Problems. Academic Press, 2002. Avant d’appliquer le théorème de Bayes à notre problème, il est nécessaire d’introduire différentes définitions.

Tout d’abord, les observations $Y$ et les paramètres du modèle $X$ sont des grandeurs physiques qui peuvent être considérées comme des variables aléatoires et donc définies par une densité de probabilité. On définit une densité de probabilité homogène (indicé $^h$ ) lorsque qu’aucune information a priori ne permet d’établir un lien entre les variables aléatoires $X$ et $Y$ . Elle est notée $\mu _{X,Y}^h (X,Y)$ . Au couple de variables aléatoires $(X,Y)$ est associé une densité de probabilité jointe $\Theta (X,Y)$ qui décrit les corrélations possibles entre $X$ et $Y$ , ainsi que les incertitudes associées au modèle théorique de façon symétrique Tarantola, 2005²Albert Tarantola. Inverse Problem Theory and Methods for Model Parameter Estimation. Society for Industrial and Applied Mathematics, 2005. . Il est possible que les informations a priori ne soient pas complètement indépendantes. Il n’est pas exclu qu’une partie des observations $Y$ ait contribué à produire les paramètres d’ébauche $X$ . On introduit alors la densité de probabilité jointe des informations a priori: $\rho_{X,Y}^b (X,Y)$ où l’indice $^b$ fait référence à l’ébauche (« background » en anglais).

L’objectif de cette modélisation inverse est de trouver une valeur vraie à $X$ en se basant sur les informations a priori à notre disposition. Afin de simplifier l’écriture du théorème de Bayes, nous allons maintenant présenter les hypothèses sur lesquelles repose l’inférence bayésienne:

Dissymétrie des rôles de $X$ et $Y$ dans le modèle théorique direct.
Le qualificatif « direct » signifie que l’on cherche la probabilité d’avoir une mesure $Y$ connaissant la valeur de $X$ ( $\rho (Y\vert X)$ ). Ceci conduit à la réécriture de la densité de probabilité conditionnelle:

(1) $\begin{equation*} \Theta (X,Y) = \rho (Y\vert X) \mu_X^h (X) \end{equation*}$

Indépendance de l’information a priori et des observations.
On introduit alors les densités de probabilité associées à l’ébauche ( $\rho_X^b (X)$ ) et aux observations ( $\rho_Y$ ). Dans ce cas, la densité de probabilité jointe $\rho_{X,Y}^b (X,Y)$ s’écrit:

(2) $\begin{equation*} \rho_{X,Y}^b (X,Y) = \rho_X^b (X) \rho_Y \end{equation*}$

La limite homogène de la densité de probabilité jointe.
En utilisant la limite homogène de la densité de probabilité jointe (équation 2), il est possible de décomposer $\mu _{X,Y}^h (X,Y)$ en deux densités de probabilité homogènes $\mu _X^h (X)$ et $\mu _Y^h (Y)$ :

(3) $\begin{equation*} \mu _{X,Y}^h (X,Y) = \mu _X^h (X) \mu _Y^h (Y) \end{equation*}$

Afin d’obtenir l’état d’information a posteriori $\rho_{X,Y}^a (X,Y)$ , Mosegaard et Tarantola, 2002³K. Mosegaard and A. Tarantola. International Handbook of Earthquake and Engineering Seismology (Part A). Chapter : Probabilistic Approach to Inverse Problems. Academic Press, 2002 s’appuient sur l’opération de conjonction de probabilité simplifiée dans le cadre de l’inférence bayésienne:

(4) $\begin{equation*} \rho_{X,Y}^a (X,Y) = k \frac{\Theta (X,Y) \rho_{X,Y}^b (X,Y)}{\mu_{X,Y}^h (X,Y)} = k \frac{\rho (X\vert Y) \rho_X^b (X) \rho_Y}{\mu _Y^h (Y)} \end{equation*}$

avec $k$ une constante de normalisation.
En intégrant cette équation sur l’espace des observations ( $D$ ), on obtient l’expression de l’information a posteriori sur $X$ :

(5) $\begin{equation*} \rho_{X}^a (X) = \int_D \rho_{X,Y}^a (X,Y) dY = k \rho_X^b (X) \int_D \frac{\rho (X\vert Y) \rho_Y}{\mu _Y^h (Y)} dY \end{equation*}$

Dans le but d’obtenir un contrôle optimal sur $X$ , on introduit la notion de fonction coût, $J(X)$ , qui traduit les contraintes du problème imposées à $X$ . Cette fonction coût mesure la distance entre les observations $Y$ et la prédiction du modèle. L’optimalité est obtenue lorsque le critère de maximum de vraisemblance selon Mosegaard et Tarantola, 2002⁴K. Mosegaard and A. Tarantola. International Handbook of Earthquake and Engineering Seismology (Part A). Chapter : Probabilistic Approach to Inverse Problems. Academic Press, 2002 est atteint. Ce critère d’optimalité peut se traduire par la recherche du minimum de la fonction coût:

(6) $\begin{equation*} J(X) = -ln \frac{\rho_{X}^a (X)}{\mu_X^h} \end{equation*}$

(7) $\begin{equation*} X^a = argmin_X J(X) \end{equation*}$

Définition des âges glace et gaz

Les âges glace et gaz en fonction de la profondeur sont déterminés à partir de trois grandeurs glaciologiques:

le taux d’accumulation $A(z)$ ,
la fonction d’amincissement $T(z)$ ,
la profondeur de piégeage de l’air en équivalent glace (LIDIE) $C(z)$ .

Le taux d’accumulation correspond à la quantité de neige précipitée à la surface de la calotte sur l’année. Il est généralement exprimé en $cm/an$ ou son équivalent en glace $cm_{ie}/an$ . La fonction d’amincissement est définie comme le rapport entre l’épaisseur réelle d’une couche annuelle à la profondeur $z$ et son épaisseur d’origine lors de son dépôt en surface en équivalent glace. La fonction d’amincissement est donc sans unité. On peut directement relier l’épaisseur d’une couche annuelle au taux d’accumulation et à la fonction d’amincissement: $L(z) = A(z) \cdot T(z)$ .

Ainsi, la relation profondeur-âge d’une carotte de glace correspond à l’intégrale de l’inverse de l’épaisseur de couche annuelle. On définit l’âge glace $\Psi (z)$ comme:

(8) $\begin{equation*} \Psi (z) = \int _0^z \frac{D(z')}{L(z')}dz' = \int _0^z \frac{D(z')}{T(z')\cdot A(z')}dz' \end{equation*}$

avec $D(z)$ la densité relative du matériau (neige/glace) par rapport à la glace: $D (z) = \rho (z) / \rho _{glace}$ .

Pour pouvoir calculer l’âge du gaz, $\chi (z)$ , il est nécessaire de connaître la différence de profondeur entre du gaz et de la glace ayant le même âge, appelée $\Delta depth$ . Ce $\Delta depth$ va donc correspondre à la profondeur de piégeage de l’air au moment du dépôt, qui a ensuite subi un amincissement lors de son enfouissement jusqu’à la profondeur $z$ (équation ref{eq_ddepth}). On en déduit donc l’âge du gaz à la profondeur $z$ comme étant l’âge de la glace à la profondeur $(z-\Delta depth)$ (équation 10).

(9) $\begin{equation*} \Delta depth(z) = C(z)\cdot T(z) \end{equation*}$

(10) $\begin{equation*} \chi (z) = \Psi (z - \Delta depth(z)) \end{equation*}$

On définit aussi ce qu’on appelle le $\Delta age$ comme étant la différence d’âge entre la glace et le gaz à une même profondeur:

(11) $\begin{equation*} \Delta age(z) = \Psi (z) -\chi (z) \end{equation*}$

Formulation du modèle d’inversion

Fonctions de correction

Nous avons vu dans la section précédente les différents paramètres glaciologiques nécessaires pour le calcul des âges gaz et glace : le taux d’accumulation, la fonction d’amincissement et la LIDIE. L’objectif de Datice est d’améliorer les ébauches de ces paramètres par modélisation inverse. L’inversion se fait simultanément sur $N$ forages où chaque forage $k$ ( $\in 1,N$ ) possède des ébauches pour le taux d’accumulation ( $A^{b,k}(z)$ ), la fonction d’amincissement ( $T^{b,k}(z)$ ) et la LIDIE ( $C^{b,k}(z)$ ). Datice cherche à obtenir de nouvelles estimations de ces paramètres glaciologiques sous la forme de corrections des ébauches:

(12) $\begin{equation*} A^k (z) = \alpha ^k (z) A^{b,k} (z) \end{equation*}$

(13) $\begin{equation*} T^k (z) = \tau ^k (z) T^{b,k} (z) \end{equation*}$

(14) $\begin{equation*} C^k (z) = \gamma ^k (z) C^{b,k} (z) \end{equation*}$

où les fonctions $\alpha^k (z)$ , $\tau ^k (z)$ et $\gamma ^k (z)$ sont respectivement les fonctions de correction pour l’accumulation, l’amincissement et la LIDIE. Datice cherche alors à déterminer les fonctions de correction optimales simultanément pour chacun des forages. A partir de cela, l’outil est en mesure de calculer les paramètres glaciologiques optimisés $A^{a,k} (z)$ , $T^{a,k} (z)$ et $C^{a,k} (z)$ , également appelés grandeurs analysées. Datice peut ensuite en déduire les chronologies gaz et glace optimales pour chacun des forages.

Ainsi, le vecteur $X$ , qui représente les paramètres du modèle, est caractérisé par les fonctions de correction des trois paramètres glaciologiques pour chacun des forages:

(15) $\begin{equation*} X^k (z) = \left( \begin{array}{c} \alpha ^k (z)\\ \tau ^k (z)\\ \gamma ^k (z)\\ \end{array} \right) \qquad \text{et} \qquad X (z) = \left( \begin{array}{c} X^1 (z)\\ \vdots\\ X^k (z)\\ \vdots\\ X^N (z) \end{array} \right) \end{equation*}$

Le vecteur $X$ est ce qu’on appelle la variable d’état du problème inverse.
Par leur nature, et par les définitions 8 et 10 qui permettent de calculer les âges glace et gaz, le taux d’accumulation, la fonction d’amincissement et la LIDIE ne peuvent prendre que des valeurs strictement positives. De ce fait, les phénomènes tels que l’érosion en surface par le vent, conduisant à une accumulation négative, ou encore le retournement de couches de glace au cours de l’écoulement ne sont pas pris en compte. Ces contraintes sur les grandeurs glaciologiques se propagent aussi sur les fonctions de correction. Ainsi, les composantes du vecteur $X (z)$ sont toutes strictement positives.

Changement de variable

Du fait que les paramètres caractérisant le vecteur $X$ soient strictement positifs, les composantes du vecteur $X$ sont alors mieux décrites par des densités de probabilité lognormales.
Les grandeurs $X^{b,k}$ et $B^k$ caractérisant les ébauches sont donc assimilées à la médiane et à une mesure de dispersion.
Comme il est plus facile de travailler avec des densités de probabilité gaussiennes, on effectue le changement de variable suivant:

(16) $\begin{equation*} \tilde{X} = \ln X \end{equation*}$

(17) $\begin{equation*} \tilde{X}^{b,k} = \ln X^{b,k} \end{equation*}$

Ainsi, le vecteur d’ébauche $\tilde{X}^{b,k}$ correspond à la moyenne/médiane/mode de la distribution et $B^k$ la matrice de covariance d’erreur d’ébauche.

La variable de contrôle du problème est alors $\tilde{X}$ . Le développement de la fonction coût et l’optimisation du problème se font donc à partir de $\tilde{X}$ .

Les paramètres d’ébauches

Il nous est possible de déterminer la valeur du vecteur d’ébauche dans l’espace de contrôle, $\tilde{X}^{b,k}$ , en nous intéressant au vecteur $X^{b,k}$ puisque ces deux vecteurs sont reliés par la relation 17. Au regard de la définition des fonctions de correction, le vecteur $X^{b,k}$ est défini au travers de l’information a priori dont on dispose sur les paramètres glaciologiques. Rappelons que ces grandeurs ( $A^{b,k}$ , $T^{b,k}$ et $C^{b,k}$ ) sont obtenues à l’aide de modèles glaciologiques et correspondent aux valeurs les plus probables des grandeurs $A^k$ , $T^k$ et $C^k$ a priori. En terme de vecteur d’ébauche, ceci se traduit par:

(18) $\begin{equation*} X^{b,k}_i = 1 \qquad \text{et} \qquad \tilde{X}^{b,k}_i = 0 \qquad \text{pour } i = 1,..,n^k \end{equation*}$

La matrice de covariance d’erreur d’ébauche $B^k$ associée à la carotte de glace $k$ mesure les variances et les covariances d’erreur a priori sur le vecteur d’ébauche $\tilde{X}^{b,k}$ . Cette matrice est de la forme:

(19) $\begin{equation*} B^k = \left(\begin{array}{lcr} B_{\alpha}^k & B_{\alpha \gamma}^k & B_{\alpha \tau}^k \\ B_{\alpha \gamma}^k & B_{\gamma}^k & B_{\gamma \tau}^k \\ B_{\alpha \tau}^k & B_{\gamma \tau}^k & B_{\tau}^k \end{array} \right) \end{equation*}$

Les blocs diagonaux correspondent aux covariances d’erreur sur $\tilde{\alpha}^k$ , $\tilde{\gamma}^k$ et $\tilde{\tau}^k$ . Les blocs non-diagonaux rendent compte des covariances d’erreur croisées entre $\tilde{\alpha}^k - \tilde{\gamma}^k$ , $\tilde{\alpha}^k - \tilde{\tau}^k$ et $\tilde{\gamma}^k - \tilde{\tau}^k$ . On se restreint aux cas de matrices $B^k$ inversibles en faisant l’hypothèse que la matrice $B^k$ est définie strictement positive.

Dans le cadre des applications présentées dans cette thèse, les hypothèses suivantes ont été considérées:

pas de corrélation d’erreur entre les scénarios d’ébauche de 2 carottes de glace différentes,
pas de corrélation d’erreur entre $\tilde{\alpha}^{k}$ , $\tilde{\tau}^{k}$ et $\tilde{\gamma}^{k}$ pour une même carotte.

De ce fait, la matrice $B$ est seulement composée des blocs diagonaux $B^k$ associés à chaque forage:

(20) $\begin{equation*} B^k = \left(\begin{array}{lcr} B_{\alpha ^{k}} & 0 & 0 \\ 0 & B_{\gamma ^{k}} & 0 \\ 0 & 0 & B_{\tau ^{k}} \end{array} \right) \end{equation*}$

avec :

(21) $\begin{equation*} [B_{\alpha ^{k}} ]_{ij} = [\sigma _\alpha ^{b,k}]_i [\sigma _\alpha ^{b,k}]_j [\rho _\alpha ^{b,k}]_{ij} \end{equation*}$

(22) $\begin{equation*} [B_{\gamma ^{k}} ]_{ij} = [\sigma _\gamma ^{b,k}]_i [\sigma _\gamma ^{b,k}]_j [\rho _\gamma ^{b,k}]_{ij} \end{equation*}$

(23) $\begin{equation*} [B_{\tau ^{k}} ]_{ij} = [\sigma _\tau ^{b,k}]_i [\sigma _\tau ^{b,k}]_j [\rho _\tau ^{b,k}]_{ij} \end{equation*}$

Les écart-types $\sigma _\alpha ^{b,k}$ , $\sigma _\gamma ^{b,k}$ et $\sigma _\tau ^{b,k}$ sont déduits à partir de l’analyse des erreurs commises sur les paramètres d’ébauche $A^{b,k}$ , $C^{b,k}$ et $T^{b,k}$ respectivement. On rappelle que ces paramètres sont généralement obtenus grâce à des modèles glaciologiques.

Dans le cas de l’accumulation et la LIDIE, l’hypothèse de non-corrélation d’erreur est discutable. En effet, ces deux grandeurs étant liées, leurs fonctions de correction devraient présenter une corrélation entre elles.

Les observations

Datice peut prendre en compte différents types d’observations pour contraindre les chronologies. On peut les classer en deux groupes:

les marqueurs d’âge et de grandeurs qui ne contraignent qu’une seule carotte de glace (marqueurs d’âge glace et gaz, marqueurs de $\Delta depth$ ).
les liens stratigraphiques définis entre deux carottes de glace (gaz et glace).

Ainsi, il est possible de décomposer le vecteur d’observation sous la forme:

(24) $\begin{equation*} Y = \left( \begin{array}{c} Y^*\\ \Delta Y \end{array} \right) \end{equation*}$

Avec $Y^*$ le vecteur regroupant les contraintes spécifiques à un seul forage, et $\Delta Y$ celui regroupant les observations se rapportant à un couple de carottes.

Plusieurs hypothèses sont faites quant aux incertitudes associées à ces observations dans le cadre du problème inverse:

les densités de probabilité des différentes observations sont correctement décrites par des gaussiennes,
les incertitudes sur les marqueurs d’âge et de grandeurs sont indépendantes entre les forages,
les incertitudes sur les liens stratigraphiques sont indépendantes entre les forages,
les incertitudes sur les marqueurs appliqués à un forage et sur un couple de forages sont indépendantes,
les incertitudes des différents marqueurs sont toutes indépendantes les unes des autres,
les incertitudes sur les profondeurs des observations sont négligeables devant les incertitudes sur les mesures.

Toutes ces hypothèses permettent de séparer les densités de probabilité de $Y^*$ et $\Delta Y$ en produit de densités marginales, simplifiant la dérivation de la fonction coût.

Pour la suite, on a besoin de définir l’opérateur d’observation généralisé $h(\tilde{X})$ . Ce vecteur permet, étant donnée une valeur pour le vecteur de contrôle $\tilde{X}$ , de prédire le vecteur d’observation $Y$ moyennant une erreur:

(25) $\begin{equation*} Y = h(\tilde{X}) + \eta \end{equation*}$

Du fait de la décomposition faite précédemment du vecteur $Y$ , l’opérateur d’observation s’écrit:

(26) $\begin{equation*} h(\tilde{X}) = \left( \begin{array}{c} h^*\\ \Delta h \end{array} \right) \qquad \text{et} \qquad \begin{array}{c} Y^* = h^*(\tilde{X}) + \eta\\ \Delta Y = \Delta h(\tilde{X}) + \eta ' \end{array} \end{equation*}$

De plus, on fait l’hypothèse que les incertitudes théoriques associées au modèle directe sont négligeables par rapport aux incertitudes de mesure.

Fonction coût et minimisation

Fonction coût

Comme chaque forage possède sa propre discrétisation en fonction de la profondeur, le vecteur $\tilde{X}$ peut devenir très grand. En conséquence, il a été choisi que la densité de probabilité \textit{a posteriori} serait déterminée à l’aide du critère d’optimalité de maximum de vraisemblance qui fait intervenir la notion de fonction coût $J(\tilde{X})$ .Comme chaque forage possède sa propre discrétisation en fonction de la profondeur, le vecteur $\tilde{X}$ peut devenir très grand. En conséquence, il a été choisi que la densité de probabilité \textit{a posteriori} serait déterminée à l’aide du critère d’optimalité de maximum de vraisemblance qui fait intervenir la notion de fonction coût $J(\tilde{X})$ .
Sous les hypothèses de l’inversion bayésienne, la fonction coût peut se décomposer en un terme d’ébauche et un terme associé aux observations.

(27) $\begin{equation*} J(\tilde{X}) = J^b (\tilde{X}) + J^o (\tilde{X}) \end{equation*}$

La fonction coût détaillée pour toutes les observations prises en compte par Datice s’écrit alors:

(28) $\begin{equation*} J (\tilde{X}) & = \sum ^N_{k=1} \frac{1}{2} \left( \tilde{X}_k - \tilde{X}_k^b \right) ^T [\mathbf{B}]^{-1} \left( \tilde{X}_k - \tilde{X}_k^b \right) + \sum ^N_{k=1} \frac{1}{2} \left( Y^{ii}_k - \mathbf{h}_k^{ii} ( \tilde{X}_k ) \right) ^T [\mathbf{R}^{ii}_k]^{-1} \left( Y^{ii}_k - \mathbf{h}_k^{ii} ( \tilde{X}_k ) \right) \end{equation*}$

Dans l’équation 28, le premier terme correspond à l’ébauche, et le second terme représente de façon abrégé les différents types d’observation (où $^{ii}$ peut correspondre aux marqueurs de $\Delta depth$ , aux marqueurs d’âge glace, aux marqueurs d’âge gaz , aux liens stratigraphiques glace et gaz tous les autres marqueurs intégrés dans Datice).

Minimisation

Afin d’obtenir le vecteur $X^a$ optimisé, le critère de maximum de vraisemblance impose de trouver le minimum de la fonction coût (équation 7). Pour ce faire, l’optimisation repose sur des méthodes numériques dites de descente. On cherche à se diriger vers le minimum de $J$ , par itérations successives suivant un chemin de descente.

Datice utilise le minimiseur m1qn3 (Gilbert et Lemaréchal, 1989)⁵J.-C. Gilbert and C. Lemaréchal. Some numerical experiments with variable-storage quasi-Newton algotithms. Mathematical Programming, 45 :407–435, 1989.. Il est adapté pour résoudre des problèmes de grande échelle peu contraints. Ce minimiseur repose sur la connaissance de la fonction à minimiser ainsi que son gradient afin de pouvoir approximer son hessien (matrice des dérivées secondes de la fonction). A partir de cela, il est possible de déduire la direction de descente de la minimisation pour la prochaine itération.

La minimisation est arrêtée lorsque le test sur le gradient est satisfait, c’est à dire lorsqu’on atteint une certaine précision définie par l’utilisateur (par exemple, la valeur demandée pour Datice est $\epsilon_g = 10^{-8}$ ). Ce critère est basé sur la norme du gradient de la fonction coût :

$\epsilon_i = \frac{\vert g_i \vert}{\vert g_1 \vert} < \epsilon_g$

avec $g_i$ le gradient de la fonction coût de la $i$ ème itération et $g_1$ celui de la première itération.

Nouvelle chronologie optimisée

Une fois le test du gradient satisfait, ou presque satisfait dans certains cas, Datice détermine le vecteur $\tilde{X}^a$ regroupant les fonctions de correction optimales. En appliquant le changement de variable inverse, on obtient le vecteur $X^a$ regroupant les fonctions de correction $\alpha (z)$ , $\tau (z)$ et $\gamma (z)$ . L’outil de datation en déduit les grandeurs glaciologiques optimales puis les chronologies correspondantes, gaz et glace, pour chacun des forages. L’atout majeur de cette méthode de datation est qu’elle calcule directement l’incertitude associée aux chronologies.

Erreur a posteriori sur $\tilde{X}^a$

Suite à la minimisation, Datice est en mesure de calculer la matrice de covariance d’erreur a posteriori sur $\tilde{X}^a$ , en faisant l’hypothèse que l’opérateur d’observation généralisé est faiblement non linéaire:

(29) $\begin{equation*} P^a \simeq \left( \hat{H} \tilde{R}^{-1} \hat{H} + \tilde{B}^{-1} \right) ^{-1} \end{equation*}$

où $\hat{H}$ est l’opérateur linéaire tangent associé à l’opérateur d’observation généralisé, $\tilde{B}$ la matrice de covariance d’erreur d’ébauche généralisée (section ??) et $\tilde{R}$ la matrice de covariance d’erreur d’observation généralisée (section ??).

Nous pouvons noter que les blocs diagonaux de la matrice $P^a$ donnent les écart-types \textit{a posteriori} sur les composantes de $\tilde{X}^a$ :

$\tilde{\sigma}_i^{a,k} \simeq \sqrt{P_{ii}^{a,k}}$

avec $i$ et $k$ les indices associés à chaque carotte de glace.

Erreur a posteriori sur la chronologie

On rappelle que le changement de variable de $X$ à $\tilde{X}$ a permis de passer de distributions lognormales à gaussiennes. La seule quantité conservée entre ces deux distributions est la médiane. Afin de bien comprendre les étapes successives pour en arriver à l’estimation de l’erreur sur les chronologies analysées, nous allons raisonner sous forme d’intervalle de confiance.

Ainsi, en considérant 1–sigma pour l’erreur a posteriori sur $\tilde{X^a}$ , cela correspond à un intervalle de confiance à 68\% . Si on prend l’exemple des grandeurs associées au taux d’accumulation, l’optimisation de $\tilde{\alpha}$ pour la carotte de glace $k$ à une profondeur donnée nous permet de dire que l’on fait une erreur de l’ordre de 32 \% en affirmant :

(30) $\begin{equation*} \tilde{\alpha} \in [ \tilde{\alpha} ^a - \sigma ; \tilde{\alpha} ^a + \sigma ] \end{equation*}$

En appliquant le changement de variable inverse ( $\exp \tilde{X} = X$ ), on en déduit l’intervalle de confiance à 68\% sur $\alpha$ :

(31) $\begin{equation*} \exp (\tilde{\alpha}) \in [ \exp (\tilde{\alpha} ^a - \sigma) ; \exp(\tilde{\alpha} ^a + \sigma) ] \end{equation*}$

(32) $\begin{equation*} \alpha \in [ \alpha ^a \exp (- \sigma) ; \alpha ^a \exp (\sigma) ] \end{equation*}$

Après conversion en accumulation, en utilisant la relation $A = \alpha A^b$ de la section ??, cela devient:

(33) $\begin{equation*} A = A^b \alpha \in [ A^b \alpha ^a \exp (- \sigma) ; A^b \alpha ^a \exp (\sigma) ] \end{equation*}$

(34) $\begin{equation*} A \in [ A^a \exp (- \sigma) ; A^a \exp (\sigma) ] \end{equation*}$

La même démarche s’applique pour la fonction d’amincissement et la LIDIE afin d’en déduire leurs intervalles de confiance à 68\% .

L’incertitude sur l’âge est donc fonction de nombreux facteurs. En considérant que tous ces facteurs sont statistiquement indépendants, on peut faire l’hypothèse que les incertitudes sur les chronologies optimales pour chaque site $\Psi^{a,k}$ sont correctement décrites par des gaussiennes.

Le vecteur de l’âge de la glace en fonction de la profondeur pour chaque forage est obtenu à partir de la relation 8. On définit la matrice de covariance d’erreur \textit{a posteriori} $Q^a$ associée à la chronologie glace $\Psi^a$ comme:

(35) $\begin{equation*} Q^a \simeq \triangledown \Psi^{aT} P^a \triangledown \Psi^a \end{equation*}$

où $\triangledown \Psi^a$ est l’opérateur linéaire tangent du modèle d’âge glace.
La matrice $Q^a$ est symétrique et ses blocs diagonaux donnent les covariances d’erreur associées à chacun des $N$ forages. L’intervalle de confiance de chacune des chronologies glace optimales est directement donné par les éléments diagonaux de la matrice $Q^a$ . Pour le forage $k$ ayant une grille de profondeur $z_{i\in 1..n}$ , l’erreur sur l’âge glace à la profondeur $z_i$ correspond à:

(36) $\begin{equation*} \sigma_i^{\Psi ,k} = [Q_{kk}^a]_{ii}^{1/2} \end{equation*}$

References[+]

References
↑1, ↑3, ↑4	K. Mosegaard and A. Tarantola. International Handbook of Earthquake and Engineering Seismology (Part A). Chapter : Probabilistic Approach to Inverse Problems. Academic Press, 2002
↑2	Albert Tarantola. Inverse Problem Theory and Methods for Model Parameter Estimation. Society for Industrial and Applied Mathematics, 2005.
↑5	J.-C. Gilbert and C. Lemaréchal. Some numerical experiments with variable-storage quasi-Newton algotithms. Mathematical Programming, 45 :407–435, 1989.