Bei der Datenassimilation geht es um die Bestimmung eines Modellzustandes aus Beobachtungen.
Bis zur zweiten Hälfte des 20. Jahrhunderts mussten Meteorologen Wetterkarten manuell auf Grundlage der zur Verfügung stehenden Daten zeichnen, was man als Analyse bezeichnet. Auch die Anfangszustände der ersten operationellen Modellläufe wurden auf diese Weise festgelegt [30]. Dies war zeitaufwendig, ineffizient und fehleranfällig. Daher begann man schnell damit, diesen Prozess zu automatisieren, was man als objective analysis bezeichnet. Die Analyse ist objektiv, da sie nun nicht mehr von der Person abhängt, von der sie erstellt wird. Im Wesentlichen handelt es sich bei solchen Analysen um Interpolationen und Umrechnungen der gemessenen Größen in Modellvariablen.
Ein Nachteil von Interpolationen ist, dass sie sehr ungenau werden können, wenn keine Messdaten in der näheren Umgebung vorliegen. Am Beginn der operationellen numerischen Wettervorhersage in den 1950er-Jahren gab es ausschließlich konventionelle Messungen, die meisten von ihnen in Europa und den USA. Auf der Südhalbkugel gab es kaum Wetterbeobachtungen. Eine naheliegende Lösung ist nun, in Regionen, in denen kaum Messungen vorgenommen werden, einen klimatologischen Hintergrundzustand oder den vorherigen Modelllauf als Ersatz für Messungen zu verwenden. Hierdurch können die Informationen aus den Gebieten mit vielen Beobachtungen mittels der vom Modell simulierten physikalischen Gesetze in die Regionen mit wenigen Beobachtungen propagieren.
Bei der sogenannten optimum interpolation geht es nun darum, Informationen aus den Beobachtungsdaten und einem Hintergrundzustand auf optimale Art und Weise miteinander zu kombinieren. Damit ist gemeint, dass das Resultat (die Analyse) möglichst nahe am realen atmosphärischen Zustand zum Analysezeitpunkt liegen soll.Dies ist nicht zwangsläufig auch der Zustand, der zur besten Vorhersage führt. Hierzu verwendet man eine Kostenfunktion $J$, welche zwei Terme enthält:
Als Formel notiert man
\[ \begin{align} J\left(\mathbf{x}\right) = J_B\left(\mathbf{x}, \mathbf{x}_B\right) + J_y\left(\mathbf{x}, \mathbf{y}\right). \end{align} \]
Zunächst stellt man sich vor, dass man nur eine Beobachtung $y$ hat und das Modell nur aus einer einzigen Zahl $x$ besteht. Als Ansatz $J'$ für $J$ orientiert man sich an der euklidischen Norm und notiert
\[ \begin{align} J'\left(x\right) = \sqrt{J_B'\left(x, x_B\right)^2 + J_y'\left(x, y\right)^2} = \sqrt{\alpha_B\left(x - x_B\right)^2 + \alpha_y\left[x - \newhat{I}\left(y\right)\right]^2}. \end{align} \]
Hierbei sind $\alpha_B$, $\alpha_y$ zwei Gewichtungsfaktoren mit
\[ \begin{align} \alpha_B + \alpha_y = 1. \end{align} \]
Da die Kostenfunktion minimert werden soll und die Wurzel streng monoton steigen ist, vereinfacht man die notation der Kostenfunktion zu
\[ \begin{align} J\left(x\right) = J_B\left(x, x_B\right) + J_y\left(x, y\right) = \frac{1}{2}\alpha_B\left(x - x_B\right)^2 + \frac{1}{2}\alpha_y\left[x - \newhat{I}\left(y\right)\right]^2. \end{align} \]
Der Faktor $\frac{1}{2}$ ist Konvention und ändert das Minimum der Kostenfunktion nicht. Da im Falle von konventionellen Beobachtungen meist mehr Modellfreiheitsgrade als Messwerte vorliegen, führt die Ersetzung
\[ \begin{align} \text{Modell} - \text{interpolierte Beobachtung} &\to \text{Beobachtung} - \text{interpoliertes Modell}\\ \Rightarrow x - \newhat{I}\left(y\right) &\to y - \newhat{h}\left(x\right)\tag{31.6}\label{eq:replace_for_oi} \end{align} \]
zu einem kleineren Problem. Hierbei ist $\newhat{h}$ ein i. A. nichtlinearer Operator, der sogenannte Beobachtungsoperator. Die Kostenfunktion lautet somit
\[ \begin{align} J\left(x\right) = \frac{1}{2}\alpha_B\left(x - x_B\right)^2 + \frac{1}{2}\alpha_y\left[y - \newhat{h}\left(x\right)\right]^2.\tag{31.7}\label{eq:3dvar_cost_one_obs} \end{align} \]
Das Minimum der Kostenfunktion erfüllt
\[ \begin{align} \nabla J = 0.\tag{31.8}\label{eq:3dvar_criterion} \end{align} \]
Diese Gleichung ist auch im mehrdimensionalen Fall notwendig und hinreichend für das Minimum: Da $J$ quadratisch ist, ist $\nabla J$ linear mit konstanten Gliedern, die man durch eine Koordinatentransformation verschwinden lassen kann. Die $\nabla J$ festlegende Matrix ist weiterhin symmetrisch. Da die Einträge auf der Hauptdiagonalen alle ungleich Null sind, ist diese Matrix invertierbar. Es gibt also nur ein $\mathbf{x}$ mit $\nabla J\left(\mathbf{x}\right) = \mathbf{0}$. Wegen $\lim_{\left|\mathbf{x}\right| \to \infty} J = \infty$ ist dies ein Minimum.
Um das Problem mit linearer Algebra bearbeiten zu können, wird $\newhat{h}$ an der Stelle $x_B$ linearisiert, also
\[ \begin{align} \newhat{h}\left(x\right) = \newhat{h}\left(x_B\right) + \newhat{H}\left(x - x_B\right) + \mathcal{O}\left[\left(x - x_B\right)^2\right], \end{align} \]
hierbei ist
\[ \begin{align} \newhat{H}\left(x\right) \coloneqq \newhat{h}'\left(x_B\right) \end{align} \]
die Ableitung von $\newhat{h}$ an der Stelle $x_B$. Setzt man dies in Glg. (31.7) ein, erhält man unter Vernachlässigung eines Ungefährzeichens
\[ \begin{align} J\left(x\right) = \frac{1}{2}\alpha_B\left(x - x_B\right)^2 + \frac{1}{2}\alpha_y\left[y - \newhat{h}\left(x_B\right) - \newhat{H}\left(x - x_B\right)\right]^2.\tag{31.11}\label{eq:3dvar_cost_simple} \end{align} \]
Die Ableitung hiervon lautet
\[ \begin{align} J'\left(x\right) = \alpha_B\left(x - x_B\right) - \alpha_y\newhat{H}\left[y - \newhat{h}\left(x_B\right) - \newhat{H}\left(x - x_B\right)\right].\tag{31.12}\label{eq:3dvar_deriv_1obs} \end{align} \]
Bei nur einer Messung und einer Modellvariable ist $\newhat{H}$ einfach eine Zahl,
\[ \begin{align} \newhat{H} = H. \end{align} \]
Die Kostenfunktion soll dimensionslos sein. Aus einer Dimensionsanalyse erhält man
\[ \begin{align} \left[\alpha_B\right] = \frac{1}{\left[x^2\right]}, & {} & \left[\alpha_y\right] = \frac{1}{\left[y^2\right]}. \end{align} \]
Ein sinnvoller Ansatz für die Gewichtungsfaktoren sind die inversen Varianzen:
\[ \begin{align} \alpha_B = \frac{1}{\sigma_B^2}, & {} & \alpha_y = \frac{1}{\sigma_y^2} \end{align} \]
Hierbei ist $\sigma_B$ die Standardabweichung des Hintergrundzustandes und $\sigma_y$ die Standardabweichung der Messung $y$. Dies kann man auch formaler begründen, indem man annimmt, dass sowohl der Fehler des Hintergrundzustandes als auch der Messung normalverteilt sind, also
\[ \begin{align} P\left(x, y\right) \propto \exp\left[-\frac{\left(x - x_B\right)^2}{2\sigma_B^2}\right]\exp\left[-\frac{\left(y - \newhat{h}\left(x\right)\right)^2}{2\sigma_B^2}\right] = \exp\left[-\frac{\left(x - x_B\right)^2}{2\sigma_B^2} - \frac{\left(y - \newhat{h}\left(x\right)\right)^2}{2\sigma_B^2}\right]. \end{align} \]
$P\left(x, y\right)$ ist die Wahrscheinlichkeitsdichte für $x$ unter der Annahme, dass $y$ eingetreten ist. Dies maximiert man, indem man den negativen Logarithmus minimiert, also
\[ \begin{align} P\left(x, y\right)\text{ maximal} \Rightarrow \frac{\left(x - x_B\right)^2}{2\sigma_B^2} + \frac{\left(y - \newhat{h}\left(x\right)\right)^2}{2\sigma_B^2}\text{ minimal.} \end{align} \]
Aus Glg. (31.12) erhält man mit Glg. (31.8)
\[ \begin{align} \frac{1}{\sigma_B^2}\left(x - x_B\right) - \frac{H}{\sigma_y^2}\left[y - \newhat{h}\left(x_B\right) - H\left(x - x_B\right)\right] &= \frac{1}{\sigma_B^2}\left(x - x_B\right) - \frac{H}{\sigma_y^2}\left[y - \newhat{h}\left(x_B\right)\right] + \frac{H^2}{\sigma_y^2}\left(x - x_B\right) = 0\nonumber\\ \Rightarrow\left(\frac{1}{\sigma_B^2} + \frac{H^2}{\sigma_y^2}\right)\left(x - x_B\right) - \frac{H}{\sigma_y^2}\left[y - \newhat{h}\left(x_B\right)\right] &= 0\nonumber\\ \Rightarrow x = x_B + \frac{H}{\sigma_y^2\left(\frac{1}{\sigma_B^2} + \frac{H^2}{\sigma_y^2}\right)}\left[y - \newhat{h}\left(x_B\right)\right] &= x_B + \frac{H}{\frac{\sigma_y^2}{\sigma_B^2} + H^2}\left[y - \newhat{h}\left(x_B\right)\right]. \end{align} \]
Hieraus folgen zwei interessante Grenzfälle:
\[ \begin{align} x &= x_B + \frac{H}{H^2}\left[y - \newhat{h}\left(x_B\right)\right] = x_B + \frac{y - \newhat{h}\left(x_B\right)}{H}. \end{align} \]
Die Analyse bestimmt sich in diesem Fall nur aus den Messungen.
Nun nimmt man an, dass ein Modellzustand nicht durch eine einzelne Zahl $x$, sondern durch einen Zustandsvektor $\mathbf{x}$ der Länge $N$ festgelegt wird. Für diesen Fall verallgemeinert die Kostenfunktion Glg. (31.11) zu
\[ \begin{align} J\left(\mathbf{x}\right) = \frac{1}{2}\left(\mathbf{x} - \mathbf{x}_B\right)^TB^{-1}\left(\mathbf{x} - \mathbf{x}_B\right) + \frac{1}{2}\frac{1}{\sigma_y^2}\left[y - \newhat{h}\left(\mathbf{x}_B\right) - \mathbf{H}\left(\mathbf{x} - \mathbf{x}_B\right)\right]^2.\tag{31.20}\label{eq:3dvar_cost_1_obs} \end{align} \]
$\mathbf{H}$ ist ein $1 \times N-$Vektor. $B \in \mathbb{R}^{N \times N}$ ist die Kovarianzmatrix der Fehler des Hintergrundzustandes. Diese Matrix hat Diagonalform, auf ihrer Hauptdiagonalen stehen die Varianzen des Hintergrundzustandes:
\[ \begin{align} B = \left(\begin{array}{cccc} \sigma_{0, B}^2 & 0 & \dots & 0\\ 0 & \sigma_{1, B}^2 & 0 & \vdots\\ \vdots & \dots & \ddots & 0\\ 0 & \dots & 0 & \sigma_{N - 1, B}^2 \end{array}\right) \end{align} \]
Hieraus folgt
\[ \begin{align} B^{-1} = \left(\begin{array}{cccc} \frac{1}{\sigma_{0, B}^2} & 0 & \dots & 0\\ 0 & \frac{1}{\sigma_{1, B}^2} & 0 & \vdots\\ \vdots & \dots & \ddots & 0\\ 0 & \dots & 0 & \frac{1}{\sigma_{N - 1, B}^2} \end{array}\right). \end{align} \]
Somit erhält man
\[ \begin{align} \frac{1}{2}\left(\mathbf{x} - \mathbf{x}_B\right)^TB^{-1}\left(\mathbf{x} - \mathbf{x}_B\right) = \sum_{i = 0}^{N - 1}\frac{\left(x_i - x_{i, B}\right)^2}{2\sigma_{i, B}^2}. \end{align} \]
Man erhält hieraus
\[ \begin{align} \left(\nabla J\right)_i = \frac{\left(x_i - x_{i, B}\right)}{\sigma_{i, B}^2} - \frac{H_i}{\sigma_y^2}\left[y - \newhat{h}\left(\mathbf{x}_B\right) - \mathbf{H}\left(\mathbf{x} - \mathbf{x}_B\right)\right]. \end{align} \]
Dies kann man in der kompakten Form
\[ \begin{align} \nabla J = B^{-1}\left(\mathbf{x} - \mathbf{x}_B\right) - \frac{\mathbf{H}^T}{\sigma_y^2}\left[y - \newhat{h}\left(\mathbf{x}_B\right) - \mathbf{H}\left(\mathbf{x} - \mathbf{x}_B\right)\right]. \end{align} \]
notieren. Aus Glg. (31.8) folgt in diesem Fall also
\[ \begin{align} B^{-1}\left(\mathbf{x} - \mathbf{x}_B\right) &= \frac{\mathbf{H}^T}{\sigma_y^2}\left[y - \newhat{h}\left(\mathbf{x}_B\right) - \mathbf{H}\left(\mathbf{x} - \mathbf{x}_B\right)\right]\nonumber\\ \Rightarrow B^{-1}\left(\mathbf{x} - \mathbf{x}_B\right) + \frac{\mathbf{H}^T}{\sigma_y^2}\left[\mathbf{H}\left(\mathbf{x} - \mathbf{x}_B\right)\right] &= \frac{\mathbf{H}^T}{\sigma_y^2}\left[y - \newhat{h}\left(\mathbf{x}_B\right)\right]\nonumber\\ \Rightarrow\left(B^{-1} + \frac{\mathbf{H}^T\mathbf{H}}{\sigma_y^2}\right)\left(\mathbf{x} - \mathbf{x}_B\right) &= \frac{\mathbf{H}^T}{\sigma_y^2}\left[y - \newhat{h}\left(\mathbf{x}_B\right)\right]\nonumber \end{align} \] \[ \begin{align} \Rightarrow\left(1 + B\frac{\mathbf{H}^T\mathbf{H}}{\sigma_y^2}\right)\left(\mathbf{x} - \mathbf{x}_B\right) &= \frac{B\mathbf{H}^T}{\sigma_y^2}\left[y - \newhat{h}\left(\mathbf{x}_B\right)\right]\nonumber\\ \Rightarrow\mathbf{x} - \mathbf{x}_B &= \left(1 + B\frac{\mathbf{H}^T\mathbf{H}}{\sigma_y^2}\right)^{-1}\frac{B\mathbf{H}^T}{\sigma_y^2}\left[y - \newhat{h}\left(\mathbf{x}_B\right)\right]\nonumber \end{align} \]
\[ \begin{align} \Rightarrow\mathbf{x} &= \mathbf{x}_B + \left(B\frac{\mathbf{H}^T\mathbf{H}}{\sigma_y^2} + 1\right)^{-1}\frac{B\mathbf{H}^T}{\sigma_y^2}\left[y - \newhat{h}\left(\mathbf{x}_B\right)\right]. \end{align} \]
Nun nimmt man an, dass nicht eine, sondern $M$ Beobachtungen vorliegen. Die Kostenfunktion Glg. (31.20) verallgemeinert man für diesen Fall zu
\[ \begin{eqnarray} J\left(\mathbf{x}\right) = \frac{1}{2}\left(\mathbf{x} - \mathbf{x}_B\right)^TB^{-1}\left(\mathbf{x} - \mathbf{x}_B\right) + \frac{1}{2}\left(\mathbf{y} - \newhat{h}\left(\mathbf{x}_B\right) - H\left(\mathbf{x} - \mathbf{x}_B\right)\right)^TR^{-1}\left(\mathbf{y} - \newhat{h}\left(\mathbf{x}_B\right) - H\left(\mathbf{x} - \mathbf{x}_B\right)\right).\tag{31.27}\label{eq:3d-var_cost_function} \end{eqnarray} \]
Hierbei ist $R$ die Kovarianzmatrix der Fehler der Beobachtungen. Da es sich hierbei vorwiegend um statistische (und keine systematischen) Fehler handelt, kann man davon ausgehen, dass diese Diagonalform hat. Es gilt
\[ \begin{align} \left(\mathbf{y} - H\mathbf{x}\right)^TR^{-1}\left(\mathbf{y} - H\mathbf{x}\right) &= \mathbf{y}^TR^{-1}\mathbf{y} - \mathbf{y}^TR^{-1}H\mathbf{x} - \mathbf{x}^TH^TR^{-1}\mathbf{y} + \mathbf{x}^TH^TR^{-1}H\mathbf{x}. \end{align} \]
Man rechnet nun
\[ \begin{align} \mathbf{y}^TR^{-1}H\mathbf{x} &= \left(R^{-1}H\mathbf{x}\right)^T\mathbf{y} = \mathbf{x}^TH^TR^{-1}\mathbf{y}. \end{align} \]
Dabei wurde verwendet, dass $R^{-1}$ orthogonal ist. Somit gilt
\[ \begin{align} \left(\mathbf{y} - H\mathbf{x}\right)^TR^{-1}\left(\mathbf{y} - H\mathbf{x}\right) &= \mathbf{y}^TR^{-1}\mathbf{y} - 2\mathbf{x}^TH^TR^{-1}\mathbf{y} + \mathbf{x}^TH^TR^{-1}H\mathbf{x}\nonumber\\ &= \mathbf{y}^TR^{-1}\mathbf{y} - 2\mathbf{x}^TH^TR^{-1}\mathbf{y} + \left(\mathbf{x}^TH^TR^{-1}H\right)^T\cdot\mathbf{x}\nonumber\\ &= \mathbf{y}^TR^{-1}\mathbf{y} - 2\mathbf{x}^TH^TR^{-1}\mathbf{y} + \left(H^TR^{-1}H\mathbf{x}\right)\cdot\mathbf{x}.\tag{31.30}\label{eq:3d-var_deriv} \end{align} \]
Dabei wurde berücksichtigt, dass wegen
\[ \begin{align} \left(H^TR^{-1}H\right)^T &= H^TR^{-1}H \end{align} \]
die Matrix $H^TR^{-1}H$ orthogonal ist. Bildet man den Gradienten von Glg. (31.30) (in Bezug auf $\mathbf{x}$), erhält man
\[ \begin{align} \nabla\left[\left(\mathbf{y} - H\mathbf{x}\right)^TR^{-1}\left(\mathbf{y} - H\mathbf{x}\right)\right] &= -2H^TR^{-1}\mathbf{y} + 2H^TR^{-1}H\mathbf{x} = -2H^TR^{-1}\left(\mathbf{y} - H\mathbf{x}\right). \end{align} \]
Hieraus folgt
\[ \begin{eqnarray} \nabla J = B^{-1}\left(\mathbf{x} - \mathbf{x}_B\right) - H^TR^{-1}\left(\mathbf{y} - \newhat{h}\left(\mathbf{x}_B\right) - H\left(\mathbf{x} - \mathbf{x}_B\right)\right). \end{eqnarray} \]
Aus Glg. (31.8) folgt also
\[ \begin{align} B^{-1}\left(\mathbf{x} - \mathbf{x}_B\right) - H^TR^{-1}\left(\mathbf{y} - \newhat{h}\left(\mathbf{x}_B\right) - H\left(\mathbf{x} - \mathbf{x}_B\right)\right) &= 0\nonumber\\ \Leftrightarrow\left(\mathbf{x} - \mathbf{x}_B\right) - BH^TR^{-1}\left(\mathbf{y} - \newhat{h}\left(\mathbf{x}_B\right) - H\left(\mathbf{x} - \mathbf{x}_B\right)\right) &= 0\nonumber\\ \Leftrightarrow\left(BH^TR^{-1}H + 1\right)\left(\mathbf{x} - \mathbf{x}_B\right) - BH^TR^{-1}\left(\mathbf{y} - \newhat{h}\left(\mathbf{x}_B\right)\right) &= 0\nonumber \end{align} \]
\[ \begin{align} \Leftrightarrow\left(BH^TR^{-1}H + 1\right)\left(\mathbf{x} - \mathbf{x}_B\right) = BH^TR^{-1}\left(\mathbf{y} - \newhat{h}\left(\mathbf{x}_B\right)\right).\tag{31.34}\label{eq:oi_lin} \end{align} \]
Dies ist das für $\mathbf{x}$ geltende lineare Gleichungssystem. Man kann dies weiter umformen zu
\[ \begin{align} \mathbf{x} &= \mathbf{x}_B + \left(BH^TR^{-1}H + 1\right)^{-1}BH^TR^{-1}\left(\mathbf{y} - \newhat{h}\left(\mathbf{x}_B\right)\right).\tag{31.35}\label{eq:oi_result_pre} \end{align} \]
Die Matrix
\[ \begin{align} K \coloneqq \left(BH^TR^{-1}H + 1\right)^{-1}BH^TR^{-1}\tag{31.36}\label{eq:gain_first_version} \end{align} \]
bezeichnet man als Gain-Matrix. Um die rechte Seite dieser Gleichung zu berechnen, muss man eine $N \times N-$Matrix invertieren.
Man kann den Ausdruck für die Gain-Matrix noch etwas vereinfachen. Sei $A \in \mathbb{R}^{N \times M}$. Dann gilt
\[ \begin{align} A &= A\left(HA + 1\right)\left(HA + 1\right)^{-1}\nonumber\\ \Leftrightarrow A &= \left(AHA + A\right)\left(HA + 1\right)^{-1}\nonumber\\ \Leftrightarrow A &= \left(AH + 1\right)A\left(HA + 1\right)^{-1}\nonumber\\ \Leftrightarrow\left(AH + 1\right)^{-1}A &= A\left(HA + 1\right)^{-1}\nonumber. \end{align} \]
Definiert man nun
\[ \begin{align} A \coloneqq BH^TR^{-1}, \end{align} \]
folgt
\[ \begin{align} K = \left(BH^TR^{-1}H + 1\right)^{-1}BH^TR^{-1} = BH^TR^{-1}\left(HBH^TR^{-1} + 1\right)^{-1}. \end{align} \]
Für die Berechnung der rechten Seite dieser Gleichung ist nun eine $M \times M-$Matrix zu invertieren, was für $M < N$ effizienter als Glg. (31.36) ist. Man kann dies weiter zu
\[ \begin{align} K = BH^T\left(HBH^T + R\right)^{-1} \end{align} \]
vereinfachen. Setzt man dies in Glg. (31.35) ein, erhält man
\[ \begin{align} \mathbf{x} &= \mathbf{x}_B + BH^T\left(HBH^T + R\right)^{-1}\left(\mathbf{y} - \newhat{h}\left(\mathbf{x}_B\right)\right).\tag{31.40}\label{eq:oi_result} \end{align} \]
Das Windfeld kann als Teil der OI assimiliert werden. Dies hat bei groben Modellauflösungen jedoch den Nachteil, dass der so erhaltene Wind häufig weiter von der Balancegleichung entfernt ist, als dies auf der synoptischen Skala üblich ist und somit zu starke Schall- und Schwerewellen auslöst. Diese lösen dann Vorhersagefehler aus. Daher kann es sinnvoll sein, das Windfeld nicht auf Grundlage von Messungen zu bestimmen, sondern durch Lösen der Balancegleichung. Dies ist ein Beispiel dafür, dass die OI die Qualität der Analyse und nicht der Vorhersage maximiert.
Geht man von Hydrostatik aus, also von Gültigkeit der hydrostatischen Grundgleichung Glg. (13.122), vereinfacht sich der thermodynamische Zustand der Atmosphäre erheblich: Üblicherweise reichen zwei thermodynamische Größen $q_1$, $q_2$ aus, um diesen festzulegen. Bei Gültigkeit einer Gleichung der Form
\[ \begin{align} \frac{\partial q_1}{\partial z} = f\left(q_2\right) \end{align} \]
kann man jedoch die Festlegung der Größe $q_1$ in der gesamten Atmosphäre durch ihre Festlegung auf einer Fläche ersetzen, zum Beispiel der Erdoberfläche oder einer Druckfläche. Die exakte Gültigkeit einer Aussage über das Resultat einer Datenassimilation zu fordern bezeichnet man als strong constraint.
Die Inversion der $M \times M-$Matrix $HBH^T + R$ ist sehr aufwendig, wenn die Anzahl der Einzelbeobachtungen groß ist. Hierzu mache man sich klar, dass beispielsweise ein Radarbild aus einer sehr großen Zahl von Einzelbeobachtungen besteht. Dies ist ein Nachteil der OI. Diesen könnte man umgehen, indem man die Ersetzung Glg. (31.6) rückgängig macht. Dies ist bei Fernerkundungsdaten nicht direkt möglich, da man beispielsweise eine spektrale Strahldichte an einem Satelliten nicht ohne weiteres in Temperaturen umrechnen kann. Eine Lösung hierfür wird in diesem Abschnitt dargestellt.
Beim 3D-Var-Verfahren geht es darum, Glg. (31.40) iterativ ohne explizite Berechnung der Gain-Matrix zu lösen. Definiere
\[ \begin{align} \delta\mathbf{x} &\coloneqq \mathbf{x} - \mathbf{x}_B,\\ \mathbf{d}_B &\coloneqq \mathbf{y} - \newhat{h}\left(\mathbf{x}_B\right),\\ \end{align} \]
Glg. (31.34) lautet dann
\[ \begin{align} \left(BH^TR^{-1}H + 1\right)\delta\mathbf{x} = BH^TR^{-1}\mathbf{d}_B. \end{align} \]
Multipliziert man mit $B^{-1}$, erhält man
\[ \begin{align} \left(H^TR^{-1}H + B^{-1}\right)\delta\mathbf{x} = H^TR^{-1}\mathbf{d}_B. \end{align} \]
Hier muss also eine $N \times N-$Matrix invertiert werden. Hieraus würde man erwarten, dass diese Gleichung für $N < M$, falls also mehr Beobachtungen als Modellfreiheitsgrade vorliegen, effizienter ist als Glg. (31.40).
Ein weiterer Vorteil von 3D-Var gegenüber OI besteht darin, Nebendingungen ohne großen Aufwand einfordern zu können, indem man sie in die Kostenfunktion aufnimmt. Dies bezeichnet man als weak constraint. Das in Abschn. 31.2.4 dargelegte Verfahren kann man als weak constrain formulieren, indem man einen Term
\[ \begin{align} J \propto \frac{1}{2}\int_A\left(\frac{\partial p}{\partial z} + g\rho\right)^2d^3r \end{align} \]
in die Kostenfunktion aufnimmt.
Der Horizontalwind ist auf der synoptischen Skala aufgrund der näherungsweisen Gültigkeit des geostrophischen Gleichgewichts sowie der Balancegleichung eng mit den thermodynamischen Größen verknüpft. Daher kann man ihn als diagnostische Größe betrachten. Ihn direkt aus der Geostrophie in Form von Glg. (13.187) zu berechnen ist aufgrund des in Abschn. 13.10.1 dargestellten Divergenzproblems in den Tropen nicht möglich. Daher kann man anstattdessen die Balancegleichung (15.137) oder die linearisierte Balancegleichung Glg. (15.138) verwenden. Diese filtert sowohl Schall- als auch Schwerewellen, da beide mit Divergenzen verbunden sind.
Analog zum Fall der Hydrostatik ist es auch hier wieder möglich, die Gültigkeit der Balancegleichung entweder exakt zu fordern, oder, was einfacher ist, sie in die Kostenfunktion aufzunehmen.
4D-Var basiert auf 3D-Var, berücksichtigt jedoch zusätzlich noch die Zeitdimension. Hierzu verallgemeinert man die Kostenfunktion $J$ in der Form
\[ \begin{eqnarray} J = J_{\text{3D-Var}} + \sum_{n = 1}^NJ_n, \end{eqnarray} \]
hierbei ist $N$ die Anzahl der involvierten Zeitschritte (außer dem Analysezeitpunkt selbst) und $J_n$ ist die Kostenfunktion an Schritt $n$. $J_{\text{3D-Var}}$ ist wie in Glg. (31.27) definiert. $\mathbf{Y}_n$ sei der Vektor der Beobachtungen am Zeitschritt $n$. An jedem Zeitschritt kann dieser Vektor eine andere Größe sowie eine andere Bedeutung haben. $\mathbf{X}_n$ sei der Zustandsvektor der Modellatmosphäre bei Zeitschritt $n$, welcher aus dem Analysezustand $\mathbf{X}$ über
\[ \begin{eqnarray} \mathbf{X}_n = M_n\mathbf{X}, \end{eqnarray} \]
ermittelt wird. Hierbei ist $M_n$ die Integration des Modells vom Analysezeitpunkt zu Zeitschritt $n$. Um diese Integration zu linearisieren, berechnet man einen sogenannten tangentialen linearen Modelloperator $M$, welcher eine Matrix ist. Es gilt dann
\[ \begin{eqnarray} M_n\approx M^n, \end{eqnarray} \]
was für zunehmende $n$ immer ungenauer wird. Somit gilt
\[ \begin{eqnarray} J_n = \left(\mathbf{Y}_n - H_nM^n\mathbf{X}\right)^TA_{O, n}\left(\mathbf{Y}_n - H_nM^n\mathbf{X}\right). \end{eqnarray} \]
Liegt der Anfangszustand eines anderen Modells vor (beispielsweise auf einem Länge-Breite-Gitter), ist es durchaus möglich, diesen einfach auf das gewünschte Modellgitter zu interpolieren und in die prognostischen Variablen umzurechnen. Hierbei wird die Analyse des anderen Modells als „ Beobachtung“ aufgefasst und es wird die Annahme gemacht, dass der Fehler der Datenassimilation des anderen Modells vernachlässigbar ist im Vergleich zum Fehler der Vorhersage des eigenen Modells. Diese Annahme ist häufig gerechtfertigt, und selbst wenn sie nicht gerechtfertigt ist, kann es sinnvoll sein, trotzdem diese Methode zu verwenden, um die durch die Einfachheit dieses Verfahrens freiwerdenden Rechenkapazitäten für den Modelllauf zu nutzen.