Changeset 174 for altifloat/doc

Timestamp:

09/21/15 10:31:16 (9 years ago)

Author:

leila_ocean

Message:

method 80%

Location:

altifloat/doc/ocean_modelling

Files:

• Draft1.tex (modified) (4 diffs)
• fig/Dts_win24_f14_tf72.fig (modified) (previous)
• fig/Dts_win24_f14_tf72.pdf (modified) (previous)

altifloat/doc/ocean_modelling/Draft1.tex

@@ -155 +155 @@
 
 \section{Data}
-\subsection {altimetric data}
+\subsection {Altimetry data}
 \subsection{Drifter data}
 \subsection {Model data}
@@ -166 +166 @@
 \subsection{Statement of the problem}
 
-We consider $N_f$ Lagrangian drifters released at the time $t=0$ at various locations. These drifters provide their positions each $\Delta t$, over a period $[0,t_f]$, where $t_f=N \Delta t.$ 
-Our objective is to determine an estimate of the Eulerian two-dimensional velocity field  \begin{equation}\notag
+We consider $N_f$ Lagrangian drifters released at time $t=0$ at various locations. These drifters provide their positions each $\Delta t$, over a period $[0,t_f]$.
+Our objective is to determine an estimate of the two-dimensional Eulerian velocity field  \begin{equation}\notag
 \mathbf{u}(x,y,t)=(u(x,y,t),v(x,y,t)) 
 \end{equation}
-on a specified grid, in a time frame $[0, t_f],$ given the observations 
+characterized by a length scale $R$ [Refs], given observations of the drifters' positions 
 \begin{equation}
-\bo{r}^{obs}_i(n\Delta t), \,\,\, i=1,2, \cdots, N_f, \,\,\, n=1, 2, \cdots N
+\bo{r}^{obs}_i(n\Delta t), \,\,\, i=1,2, \cdots, N_f, \,\,\, n=1, 2, \cdots N, \,\,\, \text{where}\,\,\, N \Delta t= t_f.
 \end{equation}
-Based on a variational assimilation approach (REFS), where we 
-match observations  
-to a model $\mathbf{r}(t)=(x(t),y(t))=\mathcal{M} (\bo{u}, \bo{r}$ that simulates drifters' trajectories. 
-The proposed improvements are:
- \begin{enumerate}
- \item Accounting for the wind: we shall consider that $\bo{u}=\bo{u_{geo}}+\bo{u_{wind}}$ 
- \item Physical constraint of the geostrophic part to be divergence free.
- using a first guess or background made of AVISO and wind fields. details in section.
-different is: all at once, ,moving window, time varying, wind and divergence. added to the objective function
-\item smoothed correction with the moving window.
- \end{enumerate}
-
-
-
-\subsection{Model for Lagrangian data}
-
-
-The position of a specific drifter is the solution of the non-linear advection equation
+The velocity shall be estimated on a specified grid with resolution $\Delta x,$ in the time frame $[0, T_f].$ 
+
+The estimation is done following a variational assimilation approach [Refs], whereby the first guessed velocity, or background $\bo{u_b}$, is corrected by matching the observations with a model that simulates the drifters' trajectories. This correction is done using a sliding time window of size $T_w$, where we assume $\Delta t<T_w<T_L,$ and where $T_L$ is the Lagrangian time scale associated with the drifters in the concerned region [Refs]. The background field is considered to be the sum of a geostrophic component (provided by altimetry) on which we impose a divergence free constraint, and a velocity due the wind () [Refs]. The details of this procedure are given in Section 3.5.
+
+
+
+
+\subsection{Linearized model for Lagrangian data}
+
+
+
+The position of a specific drifter $\mathbf{r}(t)=(x(t),y(t))$ is the solution of the non-linear advection equation
 \begin{equation} \label{advection}
  \odif{\mathbf{r}}{t}=\mathbf{u}(\mathbf{r}(t),t), \,\,\,\,\bo{r}(0)=\bo{r}_0, \mathbf{u}(x,y,0)=\bo{u}_0.
  \end{equation}
-This equation \eqref{advection} is integrated numerically using an Euler scheme for example. Also, since the drifters may not coincide the Eulerian velocity's grid, a spatial interpolation of $\mathbf{u}$ to these positions is needed. The numerical model is therefore:
-
+This equation is integrated numerically using an Euler scheme for example. Since the drifters positions may not coincide with the Eulerian velocity's grid points, a spatial interpolation of $\mathbf{u}$ to these positions is also needed. 
+
+
+The observation operator, denoted it schematically by  $\bo{r}=\mathcal{M} (\bo{u}, \bo{r}),$ consists then of numerical advection and interpolation, and it is given by  
 \begin{equation} \label{euler_advection}
-\bo{r}(k\delta t)=\bo{r}((k-1)\delta t)+2\delta t \, interp(\bo{u}, \bo{r}((k-1)\delta t)), \,\,\,\,\, k=1,2, \cdots
+\bo{r}(k\delta t)=\bo{r}((k-1)\delta t)+\delta t \, interp(\bo{u}((k-1)\delta t), \bo{r}((k-1)\delta t)), \,\,\,\,\, k=1,2, \cdots
 \end{equation}
 %Let us denote by $\bo{r}_k = (x_k ,y_k)$ the position of the drifter at time $t_k=k\delta t$ , $\mathbf{u}_k$ the Eulerian velocity at position $\mathbf{r}_k$ and time $t_k$, and
-where $\delta t$ the time step of the scheme, typically a fraction of $\Delta t$. We can choose bilinear interpolation 
+where $\delta t$ the time step of the scheme, typically a fraction of $\Delta t$. We choose bilinear interpolation 
 \begin{align}
 interp(\mathbf{u}, (x , y )) &= \mathbf{u}_1 + (\mathbf{u}_2 -\mathbf{u}_1 )\frac{(x - x_1)}{\Delta x} + (\mathbf{u}_3 -\mathbf{u}_1 )\frac{(y -y_1 )}{ \Delta y} \\ \notag
-&+ (\mathbf{u}_1-\mathbf{u}_2-\mathbf{u}_3 + \mathbf{u}_4 )\frac{(x-x_1 )(y -y_1 )}{\Delta x \, \Delta y} 
+&+ (\mathbf{u}_1-\mathbf{u}_2-\mathbf{u}_3 + \mathbf{u}_4 )\frac{(x-x_1 )(y -y_1 )}{\Delta x \, \Delta y}, 
 \end{align}
-
+where
 \begin{align} \notag
 \mathbf{u}_1 &= \mathbf{u} (x_1 , y_1 ), \\  \notag
@@ -213 +209 @@
  \mathbf{u}_4 &= \mathbf{u}(x_1 + 1, y_1 + 1). \notag
 \end{align}
-
-
-\subsection{The Linear tangent model}
-%We write $\mathbf{r}=\mathbf{r}^b+\delta \bo{r}$ and $\mathbf{u}=\mathbf{u}^b+\delta \mathbf{u}$ and obtain
-%\begin{align} \label{Evolve}
-%\bo{r}^b_k&=\bo{r}^b_{k-2}+2\delta t \, interp(\mathbf{u}^b, \bo{r}^b_{k-1})  \\ 
-%\delta \bo{r}_k &= \delta \bo{r}_{k-2} + 2 \delta t \, \{ interp(\delta \bo{u},\bo{r}^b_{k-1})+ \delta \bo{r}_{k-1} \cdot \partial _{(x,y)} interp (\bo{u}^b,\bo{r}^b_{k-1})\} \notag
-%\end{align}
-%where $  \partial _{(x,y)} interp$ is the derivative of the ÔinterpÕ function with respect to $(x,y)$ (not differentiable at $(x_1,y_1)$). 
-%We initialize this with $\bo{r}^b_{0,1}=\bo{r}^{obs}_{0,1}$ and $\delta \bo{r}_{0,1}=\bo{0}.$
-%
-%
-%\textbf{Remark}: Technically here, $\delta \bo{u}$ and $\bo{u}^b$ are time independent. 
-%If we were to include the full physical model (denote it by $\mathcal{M}$), we would have \[\mathcal{M}(t_i,t_0)(\bo{u_0})=\mathcal{M}(t_i,t_0)(\bo{u^b+\delta \bo{u}})=\bo{u}_i\approx \mathcal{M}(t_i,t_0)(\bo{u}^b)+\mathcal{M}'(t_i,t_o) \delta \bo{u},\]   that is we would have time dependent $\bo{u}^b$ and  $\delta \bo{u}$.
-%
+Here, $x_1=\left \lfloor{x}\right \rfloor $ is the floor function and $(x_1, y_1), (x_1 + 1, y_1), (x_1, y_1 + 1)$ and $(x_1 + 1, y_1 + 1)$ are the grid points which are nearest to $(x, y)$. 
+
+
+
+%This has the advantage that the cost function becomes quadratical; it has a unique minimum and this minimum is assumed to be close to that of the full non-quadratical cost function (below). 
+
+Using the incremental approach [Refs], the nonlinear observation operator $\mathcal{M}$ is linearized around a reference state. In a specific time window, we consider time independent perturbations $\delta \bo{u}$ on top of the background velocity field, that is
+\begin{align}\label{totalR}
+\bo{r}&=\bo{r^b}+\delta \bo{r} \\ \notag
+\bo{u}&=\bo{u^b}+\delta \bo{u}. 
+\end{align}
+The linearized equations become 
+\begin{align} \label{REquations}
+&\bo{r^b}(k\delta t)=\bo{r^b}((k-1)\delta t)+\delta t \, interp\bigl(   \bo{u^b}((k-1)\delta t)), \bo{r^b}((k-1)\delta t     \bigr) ,\,\,\,\,\, \text{background}  \\  \notag
+&\bo{\delta r}(k\delta t) = \bo{\delta r}((k-1)\delta t) + \delta t \, \{ interp(\bo{\delta u},\bo{r^b}((k-1)\delta t)) \\ \notag 
+&+ \bo{\delta r}((k-1)\delta t) \cdot \partial _{(x,y)} interp \bigl(\bo{u^b}((k-1)\delta t),\bo{r^b}((k-1)\delta t)\bigr)\},  \,\,\, \text{tangent} 
+\end{align}
+where the drifters' positions are initialized with observations, and where $k=1,2,3, \cdots    \left \lfloor{T_w/\delta t}\right \rfloor .$ 
+Here, $ \partial _{(x,y)} interp$ is the derivative of the interpolation function with respect to $(x,y)$. 
+
+
+The background velocity used in the advection of the drifters is the aggregate of a geostrophic component $\bo{u}_{geo}$ provided by altimetry and a component driven by the wind $\bo{u}_{wind}$, which is parametrized by two parameters as described in PPM. So we have  
+\[\bo{u}^b=\bo{u}_{geo}+\bo{u}_{wind}
+\]
+The data for both velocities are provided daily as described in the Data Section. This means that an interpolation in time of these velocities may be needed.
+
+
+
 %\subsection{Discussion on the choice of $\delta t$ }
 %The observations $\bo{r}^{obs}$ are available to us for each $\Delta t$, that is we have $\bo{r}^{obs}(i \Delta t), \,\,  i=0,1,2, \cdots $. Let us denote by $\Delta x$ the typical distance traveled by a buoys during $\Delta t$ and let $\Delta s$ denote the spatial resolution of the correction $\delta \bo{u}$. We have two scenarios: 
@@ -234 +243 @@
 %\item If $\Delta s << \Delta x$, we need to iterate many $\delta t$ to arrive at an observation. We take $N \delta t=\Delta t$
 %\end{enumerate}
-%
-%
-%\subsection{Transformation}
-%We can transform the state vector using $\delta \bo{v}=\bo{B}^{-1/2}\delta \bo{u}$. We can define this matrix directly using smoothing operators \cite{weaver2006correlation}. This formulation has the advantage that $\delta \bo{v}$ has now the identity matrix as a covariance matrix. 
-%
-%
-%\subsection{SLIDES}
-%Minimize objective function with respect to the state $\vec{U_0}$ 
-%\[
-%\mathcal{J}(\vec{U_0})= \frac{1}{2}
-% \left \{
-% \sum_{i,m}
-% \vectornorm{ \mathcal{G} \mathcal{M}(\vec{U_0})-\vec{r}_i^{\,o}(m\Delta t) } ^2+ \vectornorm{ \vec{U_0}-\vec{U}_b}^2_{\mathbf{B}}
-% \right \}\]
-%\vspace{2mm}
-%where $\vectornorm{\psi}^2_{\bo{B}} \equiv \psi^T \mathbf{B}^{-1} \psi.$
-%
-%\begin{itemize}
-%\item $\vec{U}^b$ is background state of  $\vec{U_0}$, with error covariance matrix $\mathbf{B}^{-1}$
-%\small{
-%\begin{itemize}
-%\item  dual purpose: regularization and information spreading (smoothing)
-%%\item in extremely idealized assimilation contexts, for example with floaters "everywhere", can use identity
-%\item use the diffusion filter method [Weaver and Courtier, 2001]
-%%need R 
-%\end{itemize}
-%}
-%\item If $\mathcal{J}$ is not quadratic in $\vec{U}_0$, optimization is very costly and may not converge: work with linear tangents (around the background) of model and advection 
-%\item  Incremental approach takes into account weak nonlinearities, as models are updated several times.
-%\end{itemize}
-%
-%
-%%%%%%---
-%
-%
-%
-
-\subsection{Optimization}
-
-\begin{itemize} 
-
-\item $\vec{U_0}=\vec{U}_{b}+\vec{\delta U},$ $\vec{U}_b=\vec{U}_{geo}+\vec{U}_{wind}$ and $\vec{r}=\vec{r}^{\,b}+\vec{\delta r}$
-\item Model and Linear tangent (forward Euler for the numerical integration, and bilinear Lagrange interpolation in space for instance):  
-\small{
-\begin{align} \notag
-&\vec{r}^{\,b}_i(k\delta t)=\vec{r}^{\,b}_i((k-1)\delta t)+2\delta t \, interp(\vec{U}^b, \vec{r}^{\,b}_i((k-1)\delta t)),\,\,\,\,\, \text{model}  \\  \notag
-&\vec{\delta r}_i(k\delta t) = \vec{\delta r}_i((k-1)\delta t) + 2 \delta t \, \{ interp(\vec{\delta U},\vec{r}^{\,b}_i((k-1)\delta t)) \\ \notag 
-&+ \vec{\delta {r}}_i((k-1)\delta t) \cdot \partial _{(x,y)} interp (\vec{U}^b,\vec{r}^{\,b}_i((k-1)\delta t))\}  \,\,\, \text{tangent} \notag
-\end{align}
-initialized with observations, for $k=1,2,3, \cdots $ and $i=1, \cdots N_{drift}$ 
-}
-%\vec{r}(0)= \vec{r}^{o}(0)$
-\item Perform sequences of optimization
-\[\mathcal{J}(\delta \vec{U})=  \sum _i \sum_{m=1}^{P} \vectornorm{\vec{r}^{\,b}_i+\delta \vec{r}_i(\delta \vec{U}) -\vec{r}_i^{\,o}(m\Delta t) }^2 +\alpha_1 \,\vectornorm{ \vec{\delta U} }^2_{\bo{B}} +\alpha_2 \,div (\vec{U}_{geo}) \]
-$P$ depends on $T_{L}/\Delta t$. $T_L$ characteristic of autocorrelation of drifters (typical 1-3 days) 
-
-%+\sum_{i,m=1}^{m=P} \vectornorm{\vec{r}^{\,b}+\delta \vec{r}(\textcolor{green}{\delta \vec{u}}) -\vec{r}_i^{\,o}(m\Delta t) }^2   \right]   \]
-\end{itemize}
-
-
-
-
-
-
-
-
+
+\subsection{Algorithm for assimilation}
+
+
+
+We perform sequences of optimizations, where we minimize the following objective function with respect to the time independent control correction $\delta \bo{u},$ in a specific time window $[0,T_w].$
+\begin{equation}
+\mathcal{J}(\delta \bo{u})=  \sum _{i=1}^{N_f} \sum_{m=1}^{\left \lfloor{T_w/\Delta t}\right \rfloor} \vectornorm{\bo{r}^{\,b}_{i}(\bo{u^b})+\delta \bo{r}_i(\delta \bo{u}) -\bo{r}_i^{\,obs}(m\Delta t) }^2 +\alpha_1 \vectornorm{ \bo{\delta u} }^2_{\bo{B}} +\alpha_2 \,div (\bo{u}_{geo}) 
+\end{equation}
+
+The first component of the objective function quantifies the misfit between the model 
+obtained by iterations of \eqref{REquations}, to observations $\bo{r}^{\,obs}(m\Delta t)$. 
+We highlight the dependence of $\bo{r}^b$ on the background velocity only, whereas $\bo{\delta r}$ depends on both background and correction.
+The second component states that the corrected field is required to stay close to the background velocity. Here the $B$-norm is defined as $\vectornorm{\psi}^2_{\bo{B}} \equiv \psi^T \mathbf{B}^{-1} \psi,$ where the error covariance matrix  is $\mathbf{B}^{-1}.$ This term serves the dual purpose of regularization and information spreading or smoothing. We use the diffusion filter method of [Weaver and Courtier, 2001], where the Eulerian scale $R$ enters in ..
+The parameter $\alpha_1$ represents the relative weight of this regularization term with respect to the other terms, and it must be chosen carefully.
+The last component is a constraint on the geostrophic part of the velocity, required to stay divergence free. More here?.The total velocity may have a divergent component due to the wind. 
+
+
+Inside a time window $[0, T_w]$, we assimilate a whole trajectory of drifters, which gives one fixed correction $\bo{\delta u}.$ We want to reconstruct a smooth time dependent velocity inside $[0, T_f]$: we achieve this by sliding the window by a time shift $\sigma$ to obtain other corrections $\bo{\delta u}_k$ in $[k\sigma, k\sigma +T_w], \, \, k=0, 1, 2 \cdots$. The reconstructed velocity is then \[\bo{u_{corrected}}(t_i)=\bo{u^{b}}(t_i)+\sum_{k=0} w_k \bo{\delta u}_k.\]  A correction at at specific instant $t_i$ only takes into accounts windows sliding through $t_i$ and weight $w_i$ given to correction $\bo{\delta u}_i$ is inversely proportional to the ``distance" between time $t_i$ and the window's position. Here $i=1,2, \cdots T_f/\delta t$
+
+
+
+
+DIAGRAM IS USEFUL!
+
+\subsection{Implementation of the assimilation method in YAO}
 
 \section{Twin Experiments}