## 2 Problem Statement and Related Work  
Regarding ocean remote sensing data, spatio-temporal interpolation issues (also referred to as data assimilation in geoscience) can be regarded as the reconstruction of some hidden states from partial and/or noisy observation series [31]. Data assimilation techniques usually involve a state-space evolution model [31]:  
\[x_{t+1} =\mathcal{F}(x_{t})+\eta_{t} \tag{1}\] \[y_{t+1} =\mathcal{H}(x_{t+1})+e_{t} \tag{2}\]  
where \(t\in\{0,...,T\}\) represents the temporal resolution of our time series and \(\mathcal{F}\) the dynamical model describing the temporal evolution of the physical variables \(x\). The observation model \(\mathcal{H}\) links the observation \(y\) to the physical variable \(x\). \(\eta_{t}\) and \(e_{t}\) are random processes accounting for the uncertaintiesin the dynamical and observation models. They are usually defined as centered Gaussian processes with covariances \(Q_{l}\) and \(R_{t}\) respectively. From a probabilistic point of view, the spatio-temporal interpolation problem can be seen as a Bayesian filtering problem where the main goal is to evaluate the conditional probabilities \(p(x_{t+1}|y_{1},...,y_{t})\) (prediction distribution of the state \(x_{t+1}\) given observations up to time \(t\)) and \(p(x_{t+1}|y_{1},...,y_{t},y_{t+1})\) (posterior distribution of \(x_{t+1}\) given observations up to time \(t+1\)). Under certain assumptions over the state space model (the dynamical and observation models are linear with Gaussian uncertainties), the prediction and posterior distributions are also Gaussian and can be written as:  
\[p(x_{t+1}|y_{1},...,y_{t}) =\mathcal{N}(x_{t+1}^{-},\Sigma_{t+1}^{-}) \tag{3}\] \[p(x_{t+1}|y_{1},...,y_{t+1}) =\mathcal{N}(x_{t+1}^{+},\Sigma_{t+1}^{+}) \tag{4}\]  
with the means and covariances computed for each time \(t\) using the well known Kalman recursion  
\[x_{t+1}^{-} =Fx_{t}^{+} \tag{5}\] \[\Sigma_{t+1}^{-} =F\Sigma_{t}^{+}F^{T}+Qt\] (6) \[x_{t+1}^{+} =x_{t+1}^{-}+K_{t+1}[y_{t+1}-H_{t+1}x_{t+1}^{-}]\] (7) \[\Sigma_{t+1}^{+} =\Sigma_{t+1}^{-}-K_{t+1}H_{t+1}\Sigma_{t+1}^{-} \tag{8}\]  
with  
\[K_{t+1}=\Sigma_{t+1}^{-}H_{t+1}^{T}[H_{t+1}\Sigma_{t+1}^{-}H_{t+1}^{T}+R_{t}] ^{-1}. \tag{9}\]  
Here \(F\) and \(H_{t+1}\) corresponds respectively to some linear dynamical and observation models. The superscript (-) refers to the forecasting of the mean of the state variable \(x_{t+1}^{-}\) and of its covariance matrix \(\Sigma_{t+1}^{-}\) given observations up to time \(t\) but without the new observation at time \(t+1\). The superscript (+) refers in the other hand to the mean of the state variable \(x_{t+1}^{+}\) and of the covariance matrix \(\Sigma_{t+1}^{+}\) given all observations up to time \(t+1\). They are referred to as the assimilated mean and covariance. \(K_{t+1}\) is the Kalman gain. Kalman filters provide a sequential formulation of the Optimal Interpolation (OI) [15] which may also be solved directly knowing the space-time covariance of processes \(x\) and \(y\). For non-linear and high-dimensional dynamical systems, the Probability Density Functions (PDFs) are not Gaussian anymore and the above Kalman recursion does define their means and covariances. Ensemble Kalman methods have been proposed to address these issues. The ensemble Kalman filter and smoother [31] are the first sequential filtering techniques used reliably in the reconstruction of geophysical fields. The key idea here is to approximate the forecasting mean \(x_{t+1}^{-}\) and covariance \(\Sigma_{t+1}^{-}\) by a sample mean and covariance matrix computed by propagating an ensemble of \(M\) members, \(\{x_{t+1}^{i-}\}_{i=1}^{M}\), using the dynamical model \(\mathcal{F}\). \[x_{t+1}^{i-} =\mathcal{F}(x_{t}^{i+},i\in\{0,...,N\})\] (10) \[y_{t+1} =\mathcal{H}(x_{t+1})+e_{t} \tag{2}\]0 (11) \[Besides all its advantages, EnKF techniques do not escape the curse of dimensionality. High-dimensional systems require using large ensemble sizes \(M\) which may lead to very high-computational complexity.