A two-stage ensemble Kalman filter based on multiscale model reduction for inverse problems in time fractional diffusion-wave equations
Yuming Ba Lijian Jiang Na Ou






 A two-stage ensemble Kalman filter based on multiscale model reduction  for inverse problems in time fractional diffusion-wave equations
    
Yuming BaCollege of Mathematics and Econometrics, Hunan University, Changsha 410082, China. Email:yumingb@hnu.edu.cn.
Lijian JiangInstitute  of Mathematics, Hunan University, Changsha 410082, China. Email: ljjiang@hnu.edu.cn. Corresponding author
Na OuCollege of Mathematics and Econometrics, Hunan University, Changsha 410082, China. Email: oyoungla@hnu.edu.cn.

    
================================================================================================================================================================================================================================================================================================================================================================================




ABSTRACT




 Ensemble Kalman filter (EnKF) has been widely used in state estimation  and parameter estimation for the dynamic system where observational data is obtained sequentially in time.
 Very burdened  simulations for the forward problem are  needed to update a large number of  EnKF ensemble samples.  This will slow down the  analysis efficiency  of the EnKF for large-scale and
 high dimensional models.  To reduce uncertainty and accelerate posterior inference, a two-stage ensemble Kalman filter is presented to  improve the sequential  analysis of EnKF.
  It is known that the final posterior  ensemble may be concentrated in a small portion of the entire support of the initial prior ensemble.
   It will be much more efficient if we first build  a new prior by some  partial  observations, and construct  a surrogate only over the significant region of the new prior.
  To this end,  we construct a very  coarse model using generalized multiscale finite element method (GMsFEM) and generate a new prior ensemble in the first stage.
   GMsFEM provides a set of hierarchical multiscale basis functions supported in coarse blocks. This gives flexibility and adaptivity  to choosing degree of freedoms
   to construct a reduce model. In  the second stage,  we build an initial surrogate model based on the new prior by using GMsFEM and sparse generalized polynomial chaos (gPC)-based stochastic collocation methods.
  To improve the initial surrogate model,  we dynamically update the surrogate model, which is adapted to the sequential availability of data and the updated  analysis.
  The two-stage EnKF can achieve a better  estimation than  standard EnKF,  and   significantly improve the efficiency to update the ensemble analysis (posterior exploration).
  To enhance the applicability and flexibility in Bayesian inverse problems,  we  extend the two-stage EnKF to non-Gaussian models and hierarchical models.
  In the paper, we focus on the time fractional diffusion-wave models in porous media and investigate their  Bayesian inverse problems using the proposed two-stage EnKF.
  A few numerical examples  are carried out to demonstrate the performance of the two-stage EnKF method by taking account of parameter and structure inversion in permeability fields and source functions.



keywords: GMsFEM, time fractional diffusion-wave equation, two-stage EnKF, inverse problems



§ INTRODUCTION



The  model inputs (parameters, source, domain geometry and system structure, et. al.) in many practical systems are often unknown.  We need to identify or estimate these inputs by partial and noisy  observations to
construct predictive models and calibrate the models.  This results in inverse problems.  In the paper, we consider the inverse problems in anomalous diffusion models.
The anomalous diffusion can be roughly classified into two categories:  subdiffusion (0<γ<1) and superdiffusion (1<γ<2). Here γ denotes the fraction derivative with respect to time.
The anomalous diffusion equations are also called fractional diffusion-wave equations.
 Theoretical results such as  existence, uniqueness, stability and numerical error estimates are presented in <cit.> for some type of anomalous diffusion equations.
  The relationship between  anomalous diffusion equations and regular diffusion equations is discussed in <cit.>.
  The fractional order diffusion-wave equation as a typical fractional partial differential equation <cit.>, is a generalization of the classical diffusion and wave equation and can be used to better characterize anomalous diffusion phenomena in various  fields. The fractional diffusion-wave equations can model porous media applications,   viscoelastic mechanics, power-law phenomenon in fluid and complex network, allometric scaling laws in biology and ecology,  quantum evolution of complex systems and fractional kinetics <cit.>.

In practice, the inputs and  parameters in the anomalous diffusion models are often unknown and need to be identified based on some observation data and prior information.
The problem of identifying unknown inputs in mathematical models has been intensively studied in the framework of inverse problems and various numerical methods have been developed <cit.>. The mathematical model of inverse problem  is featured with quantities which renders  useful simulation prediction obtained by imperfect model equations and measurements.  The inverse problem  is usually ill-posed.
Many methods <cit.> such as regularization or penalty can be used to overcome the ill-posedness.  The unknown inputs (e.g., permeability field) of the anomalous diffusion models in porous media may have multiscale structure, complex geometry  patterns and uncertainty. This significantly increases the challenge of the inverse problems for these models.


Practical models  usually involve  uncertainty. Moreover,  the prior information for unknown parameters and observations are often characterized by random variables. Thus, it is desirable to treat the computational model and its
inverse problem in  statistical perspective.   Once of statistical approach for inverse problems is  Bayesian inference. The Bayesian approach <cit.>  incorporates uncertainties in noisy observations and prior information,  and can derive  the posterior probability density of the parameters, which enables us to quantify the uncertainty in the parameters.
 The popular sampling methods in Bayesian inversion are
 Markov chain Monte Carlo (MCMC) <cit.> method and its variants <cit.>,  which require  costly  computation to achieve convergence and explore the whole state space in  high dimension sample spaces.
 MCMC simulation has to run a long enough chain to give an accurate estimate, and  entails repeated solutions of the forward model.  This leads to great  challenge for solving the Bayesian inverse problem.



The ensemble Kalman filter (EnKF) is another Bayesian method. It can be seen as a reduced-order Kalman filter (KF) or a Monte Carlo implementation of KF <cit.>.
Since its introduction by Evenson in <cit.>,
 EnKF has been applied in many fields such as  oceanography, numerical weather prediction, hydrology and petroleum reservoir history matching.
   EnKF can be used for both data assimilation, where the goal is to estimate model states by incorporating  dynamical observation data,
 and inverse problem, where the objective is estimate unknown parameters appearing in models.
 Ensemble members of EnKF are forecasted in time by solving forward models  and updated by an approximate Kalman filter scheme.
 EnKF has the significant advantage that its inherent recursive process is adapted to the sequential availability of observation data in dynamical systems.
 Unlike MCMC, the ensemble samples are updated independently each other in EnKF and it is not necessary to propose samplers in a very tricky manner. Thus,  EnKF and other filter methods  have attracted much more attention in community of Bayesian inversion <cit.>.
For nonlinear models,
 the linearization of model can be avoided in EnKF using  ensemble covariance as an approximation of the posterior error covariance.
  The EnKF methods provide the first and second moments of random parameter, which are approximated by ensemble mean and ensemble covariance, respectively. Thus, the EnKF algorithm makes Gaussian approximation in a sequential manner. These approximation is accurate  for Gaussian prior models. Some insight analysis of EnKF for inverse problems has been made in recent years <cit.>.
  For the non-Gaussian models, a normal-score ensemble Kalman filter is proposed in <cit.>, where the normal-score transformation is applied to transform unknown non-Gaussian  parameters to Gaussian and make the  parameters follow marginal Gaussian distributions.

   As a Bayesian sampling method,  EnKF   needs to compute  the forward problem repeatedly.
When the forward model is computationally intensive, such as multiscale models, a direct
application of EnKF forecast  with full order model  would be computationally prohibitive. In order to significantly improve the simulation efficiency, seeking more efficient sampling from the previous posterior and  building  surrogates of the forward models <cit.> are necessary to accelerate the EnKF analysis (posterior exploration).
Multiscale models can be solved efficiently and accurately by the numerical multiscale methods in a coarse grid instead of resolving all scales in very fine grid.
 As a numerical multiscale method,   Generalized Multiscale Finite Element Method (GMsFEM) <cit.> can provide a reduce model with a good trade-off between accuracy and
 computation efficiency.  The main idea of GMsFEM is
to use a variational formulation to capture the impact of small-scale features on the coarse-scale by using multiscale basis functions.  The small-scale information is integrated into multiscale basis functions, which can be used repeatedly for  different source terms and boundary conditions of the model <cit.>.
GMsFEM has been developed to solve multiscale models with complex multiscale structures and  its convergence  is independent of the contrastness of the multiscales <cit.>.

In the framework of EnKF, the output of model depends on  random parameters.
 We  use generalized polynomial chaos (gPC)-based stochastic collocation methods to propagate prior uncertainty through the forward model in a sequential manner.
 The gPC stochastic collocation methods require only a few number of uncoupled deterministic simulations with no reformulation of the governing equations
of the forward model. We assume that the model's output is a stochastic field and admits a gPC expansion. Then we select  a set of
collocation nodes and use least-squares methods to determine the coefficients of the gPC basis functions.
 To taking account of the potential sparsity of the gPC expansion, we  use l_1 regularization  to the least-squares problem. This allows
 using much fewer samples to construct a gPC surrogate model.  The idea has been employed in Bayesian inverse problems <cit.>.
 To accelerate computation for the the l_1 regularized  least-squares problem, the lagged diffusivity fixed point method is used.
The gPC surrogate is usually constructed based on a prior density <cit.> . However, the posterior is concentrated in a small portion of the entire prior support in many inference problems <cit.>.
In EnKF, the prior is sequentially updated by incorporating new data information.
 Thus, it may be much more efficient to build surrogates only over the important region of the updated  prior than the initial  prior support.

In this paper, we propose a two-stage EnKF to take care of  the challenges and concerns mentioned above.  In the first stage,  we construct a coarse GMsFEM model with very few multiscale
basis functions, and build a new prior using standard EnKF based on the first partial  measurement data information in time.  The initial ensemble samples are drawn from the new prior for the second stage of EnKF.
By integrating GMsFEM and sparse gPC stochastic collocation method based on the new prior,  we build an initial surrogate model for the second stage.
 Because the ensemble samples are  updated by the new analysis of  EnKF,  this also sequentially updates the prior based on the new ensemble samples.
To improve the initial surrogate model,  we dynamically update the  surrogate model based on the updated prior in each assimilation step.
We note that the surrogates are constructed efficiently in EnKF procedure by using GMsFEM and sparse gPC stochastic collocation method.
By virtue of building  new priors,    we exclude the unimportant region of the posterior.
  We may use some other methods such as  ensemble smoother (ES) <cit.> to build  the new prior. In general, ES is used to estimate the parameters and states  when simulation models are typically stable functions.
   To extend the two-stage EnKF to non-Gaussian models, we integrate the proposed EnKF with normal-score transformation to broaden the applicability.
   The two-stage EnKF is also explored in hierarchical Bayesian inverse problems.  This increases  flexibility  in prior modeling for the Bayesian inference.



The structure of the paper is as follows. We begin with the general framework of EnKF for inverse problems.
In Section <ref>, we focus on the time fractional diffusion-wave models and  the surrogate model construction using GMsFEM and sparse gPC.
 Section <ref> is devoted to presenting  the two-stage EnKF based on the surrogate model.
 In Section <ref>, we present a few numerical examples to illustrate the performance of proposed EnKF  with applications of inverse problems for  time fractional diffusion-wave equations.
 Some conclusions and comments are made finally.












































§ ENSEMBLE KALMAN FILTER  FOR INVERSE PROBLEMS



  Let U be  a Hilbert space  and   𝒩  a generic  forward operator on  U for some physical system.
  We assume that the forward operator  describes the relation of parameter , state u and source term f, i.e.,

    𝒩( u;)=f,

where u∈𝒰, f∈𝒰^*, the dual space  of 𝒩. We assume that the solution of the problem has the form

    𝐮=𝒳(;f).

 Let  H be the observation  operator mapping  the model state u∈U to the observation space Y

    y=H(u()):=H () ∈Y.

Let ε be  additive Gaussian noise for observation.  Then  the observation data can be expressed by

    d= H()+ε.

In the paper, we assume that  ε is  independent of . In practical setting,  the observation data is in a finite dimensional space and can be expressed by
 
    d= H()+ε∈ℝ^n_d,

 where n_d is the dimension of observations.



 §.§ Bayesian inference using EnKF



  The EnKF was introduced by Evensen <cit.> as  a powerful  data assimilation method.
  Kalman filter is used for sequential update for states in  linear dynamical systems and Gaussian distribution. It provides  the mean and covariance information of the posterior distribution. When the prior is Gaussian, the filter gets the posterior Gaussian distribution from the joint Gaussian observation and the parameter. But for nonlinear dynamical system, EnKF has been widely used for data assimilation. In the paper,
  we use EnKF for Bayesian  inverse problems. This is a particular application of EnKF in recent years <cit.>.

 Given some observation data,  we want to estimate the parameter .  In  Bayesian context, both  and d are random variables. Thus Bayes rule gives  the posterior probability density for  by

    p(| d) ∝ p( d|)p(),

where p() is the prior distribution  before the data is observed. The data enter the Bayesian inference through the likelihood function p( d|).


 If the information of observations is incomplete,  the covariance of the noise observation may be unknown.  For this situation, we  need to estimate  the covariance of observation noise.
 Let ε  be independent and identically distributed (i.i.d.) Gaussian random vector with zero mean  and variance  σ^2, i.e.,

    ε∼ N(0,σ^2I),

where I is the n_d× n_d identity matrix. Thus the likelihood function p( d|) obeys the Gaussian distribution.
Let · be the Euclidean norm and ·_ P= P^-1/2· the weighted norm, where P is the prior's covariance matrix.
 If the prior p() is also Gaussian distribution, then

    p(| d) ∝exp(- d- H()^2/2σ^2--^b_ P^2/2),

where ^b is the mean of  prior (background information).   When σ^2 is unknown, σ^2 is a hyperparameter in the hierarchical Bayesian model. Then the
 corresponding posterior

    p(, σ^2| d) ∝ p( d|, σ^2)p()p(σ^2).

The marginal posterior of σ^2 is

    p(σ^2|,  d) ∝ p( d, |σ^2)p(σ^2).

Because the likelihood

    p( d, |σ^2)=1/(2πσ^2)^n_d/2exp(- d-H()^2/2σ^2)

belongs to the inverse-gamma family, the conjugate prior p(σ^2) can be the inverse-gamma distribution

    p(σ^2) ∝ (σ^2)^-(α+1)e^β/σ^2.

From (<ref>) and (<ref>), we get

    σ^2|,  d∼Inv-gamma(α+n_d/2, β+ d-H()^2/2).

As in <cit.>, we choose two numbers   n_s  (n_s often small and between 0.01 and 1 ) and

    σ_s^2= d-H()^2/n_d-n_p,

where ∈ R^n_p, such that α=n_s/2 and β=σ_s^2α.
  Once the posterior distribution of  is inferred, we can extract the posterior mean or the maximum a posteriori (MAP) of the unknown parameter .
  We note that the MAP estimate is equivalent to the solution of a regularized  optimization  problem.  In fact,  maximizing  the right hand side of (<ref>)  is equivalent to the minimization problem

    min_∈ R^n_p(1/2( d- H())^T R^-1( d- H())+1/2(-^b)^T P^-1(-^b)),

where R=σ^2 I. The hyperparameter σ^2 can be drawn from Inv-gamma distribution of (<ref>).
 When the observation operator H is linear,  Kalman filter (KF)  method can be derived from (<ref>) by completing the squares on the variable  and gives the following analysis

    {^a   =^b+ K( d- H(^b)),
     P^a   =( I- K H) P,
    .

where H=H and  the Kalman gain  K is  given by

    K= P H^T( H P H^T+ R)^-1.



Let k be an artificial time for the data assimilation  in a dynamic system.  We denote unknown parameter, hyperparameter and observation  data as _k ∈ℝ^n_p, σ_k^2∈ℝ and d_k∈ℝ^n_d at artificial time step k, respectively.  Then we define the artificial discrete  dynamic system

    {_k    =_k-1, 
     d_k    = H(_k)+ε_k,
    σ_k^2    ∼Inv-gamma(α+n_d/2, β+ d_k- H(_k)^2/2).
      .

In the framework of EnKF, an estimate for  is updated in each data assimilation step. The sequential update needs  forecast steps and analysis steps, which transport information of the current time to the next observation time in the forecast step. At the time step k,   we denote the forecast by ^f_k, the analysis by ^a_k, forecast error covariance matrix by P_k^f  and analysis error covariance matrix by P_k^a. Then we have
 
    {^f_k    =^a_k-1, 
    σ_k^2    ∼Inv-gamma(α+n_d/2, β+ d_k- H(_k^f)^2/2),
     P^f_k    = P^a_k-1.
    .

As in (<ref>),  the posterior is the weighted sum of observations and forecast in the analysis step k, i.e.,

    {^a_k    =^f_k+ K_k( d_k-H(^f_k)), 
     P^a_k    =( I- K_k H) P^f_k,
    .

where  Kalman gain

    K_k=P^f_kH^ T(HP^f_kH^ T+R_k)^-1,  R_k=σ_k^2 I.

 Here  H is the Jacobian matrix of H in Extended Kalman Filter (EKF) when H is nonlinear.

 It is well known that the KF and EKF are numerically scarcely affordable and the storage of a few state vectors is impossible in the high-dimensional systems. To overcome the difficulty,
    EnKF is desirable for  nonlinear data assimilation problems in high-dimensional space. The advantage of EnKF is that we apply a useful approximation to the Kalman filter to avoid propagating the first and second order statistical moments. To this end,  Monte Carlo method is  used to propagate an ensemble of realizations from the prior distribution. In  EnKF, we just update the propagating  ensemble and the Kalman gain matrix is approximated by

    K_k=Cov(Θ^f_k,Z_k)Cov(Z_k,Z_k)^-1,

where Θ^f_k is the forecast ensemble and Z_k is the ensemble of simulated observations. Thus, the forecast error covariance matrix and analysis error covariance matrix are not necessary to compute. The true mean and covariance are approximated by ensemble mean and ensemble covariance, respectively. In the paper, we make use of the stochastic analysis ensemble generation method, where the simulated observations are perturbed by simulated observation error ε^f. The  ε^f is independent of  ε.

 Let {d_1, d_2, ⋯, d_k, ⋯, d_I} be a time series of observations  and d_k∈ℝ^n_d (k=1,⋯ I). We assume that the prior distribution  of  is μ_0 (Gaussian), and ε^f_k∼N(0,σ_k^2 I), where σ_k^2 is unknown.  We initially pick  M  ensemble members for EnKF,  through which we obtain the  analysis ensemble Θ^a_I. Furthermore, the mean and covariance of Θ^a_I can be used to estimate the unknown parameter.
 The pseudo-code of EnKF algorithm is presented in Algorithm <ref>.


 We  denote (d_1^T, d_2^T, ⋯, d_k^T, ⋯, d_I^T) ∈ℝ^(I· n_d)×1  by d^1:I.  We replace H by H^1:I  in (<ref>).
  We can use ensemble smoother (ES) to do  a single global update. Then the analysis in ES is
 
    Θ^a=Θ^f+K(D-Z),

where D=[ d^1:I, ⋯,  d^1:I] ∈ℝ^(I· n_d) × M.




 §.§ EnKF for non-Gaussian model using  normal-score transformation


In general,  EnKF is a Gaussian approximation for the estimated parameter  because it  reproduces the mean and covariance. If the target distribution is Gaussian  and unimodal, EnKF inherently gives an accurate  estimation. However, if the the target distribution is non-Gaussian  or  multimodal, the approximation may  not capture the properties of target distribution. In this situation, we can invoke the normal-score  transformation, which maps non-Gaussian into Gaussian and is invertible <cit.>. We perform the normal-score  transformation  after each forecast step.  Let  F_k be  the normal-score operator at the assimilation  step k and satisfy

    q_k= F_k(_k)∼ N(0,  I).

Because  the support of cumulative distribution function (CDF) is  [0, 1], the transformation can be fulfilled by   CDF.
The normal-score transformation  renders the Gaussian random variables one by one, and the  multivariate properties of  parameter vector are also changed but  not necessary to be multi-Gaussian <cit.>.


We want to incorporate the  normal-score  transformation  into EnKF for non-Gaussian cases.   Let q^f_k and q^a_k be the forecast and analysis after the  normal-score transformation. Then the forecast step of EnKF implies that

    {^f_k    =^a_k-1, 
     q^a_k-1   = F_k-1(^a_k-1), 
     q^f_k   = q^a_k-1. 
    .

The analysis step is followed by

    { q^a_k    = q^f_k-1+ K_k( d_k-H(^f_k)), 
    ^a_k    = F^-1_k( q^a_k),
    .

where the Kalman gain matrix K_k is approximated by

    K_k=Cov(Ξ^f_k,Z_k)Cov(Z_k,Z_k)^-1.

 Here Ξ^f_k is the forecast ensemble after the transformation and Z_k is the ensemble of simulated observations.
 Let μ_0 be a non-Gaussian distribution and ε^f_k∼N(0,σ^2 I). We describe the normal-score EnKF (NS-EnKF) in Algorithm <ref>.




§ SURROGATE MODEL CONSTRUCTION USING GMSFEM AND SPARSE GPC


 For the EnKF methods presented in Algorithm <ref> and  Algorithm <ref>,  we need to repeatedly compute the forward model for all ensemble members.  This computation is very expensive when the  forward model is a complex PDE model and the number of ensemble members is large.  In order to significantly  accelerate the forward model computation, we construct a surrogate model for the forward model using model reduction methods.

 The goal is to approximate a large-scale problem in  a low dimensional space. To this end, the key idea is to choose a set of appropriate basis functions,
 which can span a good approximation space for the solution. If equation   (<ref>) is linear with respect to u,
  we can derive an algebraic system for (<ref>) as follows by applying  suitable discretization method

    K() u= f,

where u∈R^N_h is the numerical solution vector and f∈R^N_h the source vector. The N_h is the number of spatial degree of freedoms and is usually very large if we
straightforwardly solve the equation in fine grid.  We can use a model reduction method and reduce the number of basis functions to improve the efficiency.
Then we can get a reduced  algebraic system for (<ref>),

    K_r() u_r= f_r.

Let  R∈ R^N_h× M_v (M_v≪ N_h) be the matrix comprised of the M_v reduced basis functions. Then a projection reduce method implies

    K_r()=R^T K()R,     f_r=R^T f.


In order to accelerate evaluations of the posterior density for each updated parameter ensemble, we use stochastic response surface methods to construct surrogate. The solution u_r() of the reduced model can be expressed by stochastic basis functions such as polynomial chaos <cit.>, radial basis functions <cit.>, and wavelet basis functions <cit.>.
The surrogate model is constructed through the stochastic collocation method by solving  a  l_1 penalized least-square problem.
We use the lagged diffusivity fixed point method for the l_1 optimization problem and get a sparse representation for u_r() using fewer samples.




 §.§ GMsFEM


In the paper, we consider the following time fractional PDE model

    ^cD_t^γu-∇(k(x)∇ u)=f(x,t),     x∈Ω,    t∈ (0, T]

 where γ∈(0, 1) ∪(1, 2)  and is the fractional order of the derivative with respect to  time. Here we consider  the Caputo fractional derivative defined by

    ^cD_t^γu=1/Γ(m-γ)∫_0^t(t-τ)^m-γ-1∂^m u(x, τ)/∂τ^mdτ,     m-1<γ<m ,

where Γ(·) is the Gamma function and m is a positive integer. The equation (<ref>) has a close form subject to a suitable boundary condition and initial condition.


In the model equation,  k(x) usually refers to a permeability field in porous media applications.  The permeability field has heterogeneous and multiscale structure inherently and
results in a multiscale model.  We will use general multiscale finite element method (GMsFEM) to reduce the model and get a coarse GMsFEM model. This
can achieve a good trade-off between  efficiency and  accuracy for  simulating the forward model.
 We will apply   GMsFEM presented in  <cit.>  to the time-fractional diffusion-wave equation (<ref>).
For GMsFEM, we need to pre-compute a set of multiscale basis functions. To this end, the first step is to construct a snapshot space V^ω_i_snap for multiscale basis by solving  local eigenvalue problem on each coarse block ω_i,

    -div(k(x, _j)∇φ_l,j) =λ_l,j k(x, _j) φ_l,j  in ω_i,
        k(x, _j)∇φ_l,j·n⃗ = 0  on ∂ω_i,

where the samplers {_j}^N_θ_j are drawn from the prior distribution of .  By a finite element method discretization on underlying fine grid, the local eigenvalue problem can be formulated as an algebraic system,

    A(_j)φ_l,j=λ_l,j S(_j)φ_l,j,

where

    [A(_j)]_mn=∫_ω_i k(x, _j)∇ v_n∇ v_m,     
    [S(_j)]_mn=∫_ω_i k(x, _j) v_nv_m,

and v_n are the basis functions in fine grid.  We take the first M^i_snap eigenfunctions corresponding to the dominant eigenvalues for each coarse neighborhood ω_i (see Figure <ref>),
 i=1,2,⋯,N_H,  where N_H is the number of coarse nodes. Hence we  construct the space of snapshots by

    V^ω_i_snap=span{φ_l,j, 1 ≤ j ≤ N_θ, 1 ≤ l ≤ M^i_snap}.

The snapshot functions can be stacked into a matrix as

    R_snap=[φ_1,⋯,φ_M_snap],

where M_snap=N_θ× M^i_snap denotes the total number of snapshots used in the construction.
The second step  is to solve the following  local problems in the snapshot space
  
    -div(k(x, )∇ψ^i_k) =λ_k k(x, ) ψ^i_k  in ω_i,
        k(x, )∇ψ^i_k·n⃗ = 0  on ∂ω_i,

where =1/N_θ∑_j=1^N_θ_j.
We define

    {   [A]_mn=∫_ω_i k(x, )∇φ_n∇φ_m=R_snap^TA̅R_snap,
       [S]_mn=∫_ω_i k(x, ) φ_nφ_m=R_snap^TS̅R_snap,
    .

where A̅ and S̅ denote fine-scale matrices corresponding to the stiffness and mass matrices, respectively,  with the permeability k(x, ).
We choose the smallest M_i eigenvalues of the equation
  
    Aψ^i_k=λ_k Sψ^i_k

and take the corresponding eigenvectors in the snapshot  space  by setting ψ_k^i=∑_j ψ_k,j^iφ_j, for k=1,⋯,M_i, to form the reduced snapshot space, where ψ_k,j^i are the coordinates of the vector ψ^i_k.



   Let {χ_i}_i=1^N_H be a set of partition of unity functions associated with the open cover {ω_i}_i=1^N_H of Ω.
 
 Then we  multiply the partition of unity functions by the eigenfunctions to construct GMsFE space, i.e.,
  
    V_H=span{Ψ_il: Ψ_il=χ_iψ^i_l: 1≤ i≤ N_H   and   1≤ l≤ M_i}.

We can  use  a single index  for  the multiscale basis function set  {Ψ_il} and place them in  the following matrix
  
    R=[Ψ_1, Ψ_2,⋯,Ψ_M_v],

where M_v=∑_i=1^N_H M_i denotes the total number of multiscale basis  functions. We note that once the matrix R is constructed, it can be repeatedly used for simulation.

 Next we present the temporal discretization for the equation (<ref>).   When 0<γ<1,  the equation (<ref>) is the subdiffusion equation. We use the  method in <cit.> to discretize the fractional derivative and have
 
    ∫_0^t(t-τ)^-γ∂ u(x, τ)/∂τdτ    = ∑_k=0^n-1u^k+1-u^k/Δ t∫_t_k^t_k+1dτ/(t_n-τ)^γ+O(Δ t) 
       =Γ(1-γ)/s∑_k=0^n-1b_k(u^k+1-u^k)+O(Δ t),

 where u^k=u(x,t_k), 0=t_0<t_1⋯<t_n=t, t_k=kΔ t and

    {   b_k:=(n-k)^1-γ-(n-k-1)^1-γ,      k=0,1,⋯,n-1,
       s:=Δ t^γΓ(2-γ).
    .

When 1<γ<2,  the (<ref>) is the superdiffusion equation. We can use the method in <cit.> to discretize the fractional  derivative and have
 
    ∫_0^t(t-τ)^1-γ∂^2 u(x,τ)/∂τ^2dτ   =∑_k=0^n-1u^k+2-2u^k+1+u^k/Δ t^2∫_t_k^t_k+1dτ/(t_n-τ)^1-γ
       +O(Δ t)^2
       =Γ(2-γ)/s̃∑_k=0^n-1b̃_k(u^k+2-2u^k+1+u^k)+O(Δ t)^2,

 where

    {   b̃_k:=(n-k)^2-γ-(n-k-1)^2-γ,   k=0,1,⋯,n-1,
       s̃:=Δ t^γΓ(3-γ).
    .


Let U^n be  the solution  at the n-th time level. Then we have the weak formulation for the subdiffusion equation (0<γ<1),

    {1/sã(∑_k=0^n-1(U^k+1-U^k)b_k, v)+a(U^n, v)   =(f(t_n), v),  ∀ v∈ V_H
    
      (U^0,v)   =(u(x,0), v),  ∀ v∈ V_H,
      .

where s and b_k are defined in (<ref>)  and
  
    a(u,v)=∫ k(x,)∇ u∇ vdx,  ã(u,v)=∫ q(x) uvdx.

 The weak formulation for the superdiffusion equation (1<γ<2) reads

    {1/s̃ã(∑_k=0^n-1(U^k+2-2U^k+1+U^k)b̃_k, v)+a(U^n+1, v)   =(f(t_n+1), v)  ∀ v∈ V_H,
    
      (U^0,v)   =(u(x,0), v)   ∀ v∈ V_H.
      .



For subdiffusion case (0<γ<1),  we define

    c_k:=
    
    b_0    k=0,
    
    b_k-b_k-1    1≤ k ≤ n-1.

For superdiffusion case (1<γ<2), we define

    c̃_k=
    
    n^2-γ-(n-1)^2-γ    k=0,
    
    -2n^2-γ+3(n-1)^2-γ-(n-2)^2-γ    k=1,
    
    (n+2-k)^2-γ-3(n+1-k)^2-γ+3(n-k)^2-γ-(n-k-1)^2-γ    2≤ k ≤ n-1,
    
    -3+2^2-γ    k=n,
    
    1    k=n+1.


Then the weak formulation of (<ref>) can be rewritten as
  
    ã(U^n, v)+sa(U^n, v)=ã(∑_k=0^n-1U^kc_k, v)+s((f(t_n), v)).

 We assume that U^n has the expansion
  
    U^n=∑_j=1^M_v u_Hj^n Ψ_j(x),

where {Ψ_j(x) } denote the GMsFEM basis functions. Let
  
    u_H^n=(u_H1^n,u_H2^n,⋯,u_HM_v^n)^T.

Then for k=1,⋯,M_v,
  
    ∑_j=1^M_v u_Hj^nã(Ψ_j, Ψ_k)+s∑_j=1^M_vu_Hj^n a(Ψ_j, Ψ_k)
       =∑_i=1^n∑_j=1^M_v c_i u_Hj^iã(Ψ_j, Ψ_k)+s(f^n, Ψ_k).

 Let B, K and F be the weighted mass, stiffness matrices and load vector using FEM basis functions in fine grid, respectively. Then the equation gives the following algebraic system,
  
    R^TBRu_H^n+sR^TKRu_H^n=∑_i=1^nc_iR^TBRu_H^i+sR^TF.

 If we define
 
    B̃=R^TBR,      K̃=R^TKR,

 
 then u_H^n can be computed  by the iteration
  
    u_H^n=(B̃+sK̃)^-1(∑_i=0^n-1c_iB̃u_H^i+s R^TF).

By using the multiscale basis functions, the solution of (<ref>) in fine grid can be  obtained by downscaling  through  the transformation Ru_H^n.



In a similar way to solving equation (<ref>), the GMsFEM solution of (<ref>) can be computed  by the iteration
  
    u_H^n+1=(B̃+s̃K̃)^-1(∑_i=0^nc̃_iB̃u_H^i+s̃ R^TF).

We note that when GMsFEM is not applied, the full order model solution in fine grid is obtained by the iteration
  
    u_h^n=(B+s K)^-1(∑_i=0^n-1c_iBu_h^i+sF)   for  0<γ<1,

and

    u_h^n+1=(B+s̃K)^-1(∑_i=0^nc̃_iBu_h^i+s̃ F)  for  1<γ<2.



By comparing  the GMsFEM model with the full order model, we see that  the size of K̃ and B̃ are M_v × M_v, but the size of K and B are N_h × N_h  (M_v≪ N_h).
Thus a much smaller system is solved in GMsFEM.  The matrix R for multiscale basis functions is computed overhead and it can be repeatedly used for all the time levels. This significantly improves the efficiency for forward model simulations.


In the sequential data assimilation process,  ensemble members {_j} update and so  of (<ref>)  updates as well.
 We can update the GMsFEM basis matrix R  to improve the GMsFEM model.




 §.§ Stochastic collation method using  l_1 regularized  least-squares

Stochastic collocation method is an efficient approach to approximate the solution of PDEs with random inputs. In this paper, we use stochastic collocation method (SCM) to obtain an expansion
for observation operator H()    and efficiently evaluate the simulated observation ensemble.

  We use  generalized polynomial chaos (gPC) functions to represent H() by  l_1 penalized least-squares method. Let i be a multi-index with |i|=i_1+⋯+i_n_z and N_0 be a nonnegative integer. The N_0th-degree gPC expansion of H() is then approximated  by a linear combination of gPC basis  {Φ_i()}_i=1^P, i.e.,

    H()≈H^N_0():=∑_i^Pc_iΦ_i(),      P=(N_0+n_z)!/N_0!n_z!.

 The coefficients of expansion are obtained by choosing some collocation points and least-squares method. We first take Q realizations {^i}_i=1^Q of  in the support of prior distribution p(). Then for each i=1,⋯,Q, we solve a deterministic problem at the node ^i to obtain H(^i).
   After we obtain all pairs {^i, H(^i)} (i=1,⋯,Q), we are able to construct a approximation of H () such that H^N_0(^i)=H(^i) for all i=1,⋯,Q. Thus, (<ref>) can produce a system of linear equations

    A c= b,

where A∈R^Q× P is the matrix with the entries
  
    A_ij=Φ_j(^i),    i=1,⋯,Q,   j=1,⋯,P.

  and the right term b satisfing
  
    b_i= H(^i),    i=1,⋯,Q.

  If we solve the system in the ordinary least-squares method, the system (<ref>) should be overdetermined, i.e.,  Q should be much larger than P.
  This means that we need solve the forward model for a large number of samples. To reduce the computation burden, we can take much fewer samples  (Q≤ P) and use l_1 regularized least-squares, i.e.,
  
    min c_1    s.t.    Ac=b,

which is equivalent to the optimization  problem
  
    min_ cAc-b_2+α c_1 .

 We use the lagged diffusivity fixed point method <cit.> for the l_1 penalized least-squares problem.
 Due to the nondifferentiability of the l_1 norm, we take an approximation to the penalty x_1 such as  √(| x|^2+β^2), where x is a scalar and β is a small positive parameter. We denote the approximated penalty by
 
    J( c)=∑_i=1^P ψ(|c_i|^2),

where
 
    ψ(t)=2√(t+β^2).

For any v∈R^P,

    d/dτ J(c+τ v )   =   ∑_i=1^P ψ'(|c_i|^2)c_iv_i = ⟨diag(ψ'( c)) c,  v⟩,

where diag (ψ'(𝐜)) denotes the n× n diagonal matrix whose ith diagonal entry is ψ'(|c_i|^2), and ⟨· ,·⟩ denotes the Euclidean inner product on R^P. From this we obtain the gradient

    grad𝒥(𝐜)=ℒ(𝐜)𝐜 ,

where ℒ(𝐜)=diag(ψ'( c)) and  is positive semidefinite.
 For convenience, we present the lagged diffusivity fixed point method for the l_1 regularized least-squares problem (<ref>)  in Algorithm <ref>.




 We note that the forward model is  solved Q times to obtain the sampling vector b.  The accuracy of surrogate model can be ensured  using much fewer samples
 (i.e., Q<P)  when the l_1 regularized least-squares method is used. As the number of  unknown parameters increases, the simulation times for the forward model will increase significantly.
 To treat the challenge,  we use  GMsFEM to  solve the forward model on a coarse grid for each sample to  improve  the simulation efficiency.
 






§ TWO-STAGE ENSEMBLE KALMAN FILTER USING GMSFEM COARSE MODEL




In this section, we present a two-stage EnKF using GMsFEM coarse model to accelerate posterior exploration and improve the sequential assimilation performance.
The goal of inverse problems is to identify an appropriate solution which can minimize the misfit  between the forward model and measurements.
 Here it is equivalent to solving the minimization  problem (<ref>). In this paper, we use EnKF methods to solve the minimization problem.
 EnKF is a sampling  method and avoids expensive  gradient computation for solving  minimization problem.

 Although EnKF method can  avoid  lineralization and  repeated  sampling to  explore  posterior density,  it requires  to compute the forward model many times
 in each forecast step described  in algorithm <ref>.  When the number of unknown  parameters  is large, we need a large number of ensemble members to estimate .
   This implies that the computation of the forecast is very expensive.
    To improve the computation efficiency,  we construct a surrogate model based on sparse gPC  and GMsFEM to approximately represent the full order model.
    However, EnKF has an inherent  constraint assumption  for the prior distribution, which must be a Gaussian distribution.
    As we know, the support of Gaussian distribution is R, but the posterior is often concentrated in a small portion of the entire prior support in many inference problems. Thus, it may be much more efficient to build a surrogate only over the important region of a posterior than the entire prior support. Inspired by the idea,  we proposed a two-stage EnKF method.
    In the first stage, we build a new prior by a very coarse GMsFEM model, where we exclude the unimportant region of the posterior, and the initial ensemble members are drawn from the new prior to enter the EnKF assimilation process.     The second stage is the surrogate model based EnKF.  We update the surrogate model dynamically when a new analysis is obtained.

The objective of the  surrogate model is to construct a representation that quantifies the primary features of the high-fidelity model while providing the computational efficiency required for uncertainty quantification. In stationary Bayesian inference, we may need to build surrogate model only once. However,  EnKF method integrates new measurement data in each   assimilation step.
Thus, the surrogate model need to be updated sequentially.  To construct the current  surrogate model,  we exclude the unimportant region by the previous analysis.




In the first stage, we can apply standard  EnKF method based on the very first few levels of measurement data to construct the new prior.  In the second stage,
we use only a few ensemble members from the previous analysis and l_1 regularized least-squares  to build the current surrogate model, which allows fast forward model evaluations to generate observation ensemble.  The outline of two-stage EnKF is presented in Algorithm <ref>.


The algorithm <ref> can render an  accurate posterior for Gaussian models.  However, for non-Gaussian distributed parameter, we need to construct the surrogate model based on the support of non-Gaussian distribution, while the EnKF only works well for Gaussian prior. For this situation,   we use normal-score EnKF to perform non-Gaussian models.
 We note that the Legendre orthogonal polynomials can be used in surrogate model construction  if the support of unknown parameters is bounded.


In algorithm <ref> and <ref>,  I_2 is the total number of data assimilation steps. In step 3, we can use  ES method to build the new prior  when some of the very first measurement data is uninformative.






§ NUMERICAL EXAMPLES


In this section, we consider the time fractional PDE model (<ref>) and estimate the model's unknown parameters and structures using the proposed two-stage EnKF. A few numerical results
will be presented for the estimation for different unknown sources of the dynamic model.
In Subsection <ref>, we recover   a channel structure in permeability k(x) when the fractional derivative γ is known. In Subsection <ref>, we estimate  the source locations
 when the diffusion type of equation (<ref>) is  unknown. In Subsection <ref>, we will recover a permeability pattern  when measurement noise is  unknown and is treated as a hyperparameter.



For the numerical examples,  we consider a dimensionless  square domain Ω=[0,1]×[0,1] for spatial variable and (0, T] for time, and we set the initial condition as u(x,0)=0 and in addition, for the superdiffusion equation, we set

    ∂ u(x;t)/∂ t|_t=0=0.

Measurement data are generated synthetically by using FEM in a fine grid with time step Δ t=0.001, and the measurement noise is set to be σ=0.01. For any given realizations of , we solve the time fractional diffusion-wave equation using GMsFEM with time step Δ t=0.002. The regularization parameter α is set as 0.01 in the lagged diffusivity fixed point algorithm. For all numerical examples, the number of samplers in constructing the gPC   surrogate is set to be the number of the gPC basis functions for approximation. We parameterize the Gaussian random fields by
 Karhunen-Loève expansion (KLE) with a given covariance function C(x,y)  and truncate the KLE to approximately represent the random fields. For the random filed h(x,ω) by the first N terms can be represented by

    h(x,ω)= E[h(x,ω)]+∑_i=1^N √(λ)_iθ_i(ω)φ_i(x),

where E[·] is the expectation operator, θ_i(ω)∼ N(0,1) and (λ_i, φ_i) are the eigenpairs  of the eigenvalue problem

    ∫_Ω C(x,y)φ_i(y)dy=λ_i φ_i(x).

 We sort the eigenvalues in ascending order, i.e., λ_1≥λ_2⋯, and their corresponding eigenfunctions are also sorted accordingly and {θ_i}_i=1^N are uncorrelated random variables. We will compare the estimation results obtained by using standard EnKF method with the proposed  two-stage EnKF.

 



 §.§ Recover  a channel structure in permeability field


In this subsection, we consider the  subdiffusion model (<ref>) with mixed boundary condition, where Dirichlet boundary conditions is

    u(0,y;t)=1 ,       u(1,y;t)=0,

and there is no flow on the other boundaries. The source term is set as f=10, the end time is T=0.11, and the fractional derivative is given by γ=0.5. The permeability field is unknown here,
and we only  have the prior information of the permeability field, which is structured with a channel that lies between y=0 and y=1. The spatial domain is divided into 3 parts by this channel and the permeability is a constant at each subregion. Thus, we can describe the boundaries of the channel by two curves Γ_1 and Γ_2, which can be expressed as


    Γ_1(x)=∑_i=1^m_1 w_i^1ϕ_i(x),   Γ_2(x)=∑_i=1^m_2 w_i^2ϕ_i(x),

where ϕ_i(x) are interpolation basis functions. In order to reduce dimension of unknown parameters, ϕ_i(x) are set to be Karhunen-Loève expansion (KLE) basis functions, i.e. ϕ_i(x)=√(λ)_iφ_i(x), and E[Γ_i(x)]=0  (i=1,2). The covariance function is set as

    C(x_1,x_2)=ξ^2exp(-x_1-x_2^2/2l^2)

to solve the eigenvalue problem (<ref>), where ξ=1 and l^2=0.1. A  criterion which is adopted for the choice of truncated term N is

    ∑_j=1^Nλ_j/∑_j=1^∞λ_j>99.99%.

We truncate the KLE expansion for functions Γ_i(x),  (i=1,2) by this criterion, and then we get m_1=5 and m_2=5.
In order to constrain the curves in the unit square domain, we make a bijective transformation
 
    Γ̃_i(x)=1/2+1/πarctan(Γ_i(x)),    i=1, 2.

Thus, Γ̃_1(x)∈[0,1] and  Γ̃_2(x)∈[0,1] can be guaranteed, i.e., the curves lie in the physical domain. Then we can construct the level set functions corresponding to the two curves as

    L_i(x,y)=Heaviside(y>Γ̃_i(x)),    i=1,2.

Hence, the random field can be parameterized as

    log k(x,y)=c_1L_1L_2+c_2(1-L_1)L_2+c_3(1-L_1)(1-L_2),

where {c_i}_i=1^i=3 are unknowns.
Thus, we have the unknown parameter vector

    =(c_1,c_2,c_3,w_1^1,⋯,w_m_1^1,w_1^2,⋯,w_m_2^2).

In this example, the reference permeability is generated by setting

    Γ_1(x)    =    0.7+0.1sin(3π x), 
    Γ_2(x)    =    0.4+0.2sin(2π x+0.1),

and the value of each subregion is set as (c_1,c_2,c_3)=(0,4,1), which is shown in Figure <ref>(left). In the example, the number of  artificial time steps for data assimilation is  I_2=9.
 Measurements are taken at time instances 0.012+0.01I:0.002:0.018+0.01I in each data assimilation step, where I∈{1,2,⋯,I_2} and the measurement locations are distributed on the uniform 5× 5 grid of the domain [0.1, 0.9]× [0.1, 0.9] as shown in Figure <ref> (right).


The forward model is defined  on a uniform 80×80 fine grid, and we set the coarse grid 5×5 for GMsFEM simulation. We
construct the local snapshot space with dimension  M^i_snap=20 and  select M_c=10 multiscale basis functions at each block to construct the coarsen reduced order model and solve the optimization problem (<ref>).

To construct a new prior, we use 5× 10^3 ensemble members by the standard EnKF method in the first stage. Then the new prior is constructed  by incorporating data information from the first three data assimilation steps. Seven local multiscale basis functions (M_i=7) are selected in constructing the basis  matrix R to construct coarse model. Then, we construct the surrogate model after obtaining the new prior ensemble in the second stage. Eight local multiscale basis functions ( M_i=8) are selected  in constructing the matrix R to construct the gPC surrogate model in the second stage. When the order of gPC is set as N_0=3, the number of random samplers is 560 for computing vectors b and A in Section <ref>. In this stage, the number of ensemble members is set as 10^4. For a comparison, we also implement the standard EnKF, which directly uses the previous analysis as a prior to the current moment in data assimilation process.
In the standard EnKF, we also use GMsFEM with 8 local GMsFE basis functions to solve forward model for 10^4 ensemble members. But no new prior and no gPC surrogate model are build in the
standard EnKF method.  To complete the assimilation process,  the standard EnKF method takes  about 4.7 hours for CPU time, while the two-stage EnKF just take about  1 hour CPU time.
Thus, the two-stage EnKF is much more efficient than the standard EnKF.






Figure <ref> shows  the  initial prior ensemble (mean) and new prior ensemble (mean) for the  two-stage EnKF method in logarithmic scale.
We can see that the initial prior gives a very rough structure for  the channel, and  the new prior improves the channel pattern.  Then we compare the mean estimate by standard EnKF with the one by two-stage EnKF.  The posterior mean and variance for the two methods are presented   in Figure <ref>, from which we find the the two-stage EnKF  offers  more accurate estimation than standard EnKF.
The large  variance is concentrated around  the two boundaries of the channel because  two boundaries are unknown.






In order to measure the estimate accuracy in the update process, we define the relative errors ε_Γ_1 and ε_Γ_2 corresponding to the two curves after each update by

    ε_Γ_1 := Γ̃_1(x)-Γ_1(x)/Γ_1(x),   ε_Γ_2 := Γ̃_2(x)-Γ_2(x)/Γ_2(x),

where Γ_1(x) and Γ_2(x) are the reference boundary functions. We note that the first three data assimilation steps are used to produce the new prior for the two-stage EnKF.
We list  the relative errors  in table <ref>.
As expected, the relative error gradually decreases when more measurement data is used in the inference.  The table shows that  the two-stage EnKF gives  more accurate estimates  than standard EnKF.



To access the prediction using the posterior model,  we compute  the 95% credible and predictive intervals for model response at u((x, 0.5); t) and u((0.5, y); t)
for the two EnKF methods.  We note that  the realizations of model response of initial prior are  constructed by GMsFE model. This coarsen model error  is added to the estimated error variance to construct prediction intervals. As illustrated in Figure <ref> and <ref>,  both the credible interval and predictive interval become  narrower as  assimilation moves on. This means that as the uncertainty from input  decreases with respect to assimilation step and the uncertainty associated with the model fit and predictions decreases.  The observation/measurment data are almost contained in the predictive intervals when sufficient data is used in the posterior model.





For this example, the measurement noise is the additional Gaussian type. Let ^* denotes the truth  parameter, the discrepancy between the data and forward model is defined by

    [d-H(^*)]_i∼ N(0,σ^2),      i=1,⋯,n_d,

where n_d is the number of measurements and its expectation is given by

    E[d-H(θ^*)^2]=n_dσ^2.

We plot E[d-H(^a)^2] against data assimilation step in Figure <ref>. Then we can see the expectation of the discrepancy tends to n_dσ^2 as assimilation time moves on.
 This implies that  the final ensemble mean is an accurate estimate.



Figure <ref> shows the marginal posterior density against the data assimilation step in the two-stage EnKF. Although there are no truth for these parameters  {ω_i^1}_i=1^5 and {ω_i^2}_i=1^5, the data is sufficiently informative to identify a small range of values for the unknown parameters. The important region of marginal densities becomes narrower as data information gains.








 §.§ Estimate  the source locations and fractional derivative


In this section, we consider a diffusion-wave problem, i.e., the diffusion type is unknown, where the boundary conditions are the same as in Subsection <ref>. The end time is set as T=0.1, and
the source term is given by

    f(x, t)   =s_1/2πτ_1^2exp{-χ_1-x^2/2τ_1^2}[1-H(t-T_m)] 
       +∑_i=2^m s_i/2πτ_i^2exp{-χ_i-x^2/2τ_i^2}H(t-T_i-1),

where H(t) is the heaviside function, s_i is the strength, τ_i is the width and χ_i is the source location. In this example,
both χ_i and γ are unknown. We just have the prior χ_i=(χ_i^1, χ_i^2)∈[0,1]×[0,1], and γ∈ (0,2) and γ≠ 1. The prior is bounded, which can be seen as uniform distribution and γ  lies  in the  interval (0,2) for diffusion-wave equation. But its values are  uncontrollable during the data assimilation steps, the samples may run out of the  interval. Thus, the two-stage EnKF by normal-score method can be used to avoid the issue. To this end, we use a bijective map F:E →R^5, where

    E=(0,2)× (0,1)^4.

Let =(γ,χ_1,χ_2)^ T and F()=q, the map F is set as

    q_i=F_i(θ_i),   i=1,⋯,5.

The dynamic system correspondingly becomes

    {q_k-1   = F_k-1(_k-1), 
    _k    =_k-1,
    q_k    =q_k-1, 
    y_k    =H(_k-1)+ε_k,
      .

where q_k=(q^k_1,⋯,q^k_5)^ T and k is the artificial time for data assimilation, and the Kalman gain is approximated by K_k=Cov(q_k,y_k)Cov(y_k,y_k)^-1.

For the source term, we take m=2, T_1=0.05, T_2=T, and τ_i=0.1 (i=1,2), s_1=3, s_2=1, and the truth source locations are set as χ^tr_1=(0.2,0.6), χ^tr_2=(0.5,0.3) and γ^tr=0.5. Here, the total artificial time steps for data assimilation is set as I_2=8. Measurements are taken from the subdiffusion equation at time instances 0.012+0.01I:0.006:0.018+0.01I in each data assimilation step, where I∈{1,2,⋯,8} and the locations are distributed on the uniform 5× 4 grid of the domain [0.1, 0.9]× [0, 1] as shown in Figure <ref> (right). The high-contrast permeability is known and has the spatial distribution shown in  Figure <ref> (left).





The forward model is defined on 100×100 uniform fine grid, and GMsFEM is implemented on 5×5 coarse grid.
 We choose M_snap=10 eigenfunctions for the local snapshot space  and use  GMsFEM with 9 local multiscale basis functions  to approximate the forward model in solving the optimization problem (<ref>).

In the subsection, we focus on the two-stage EnKF for the inverse problem.  To this end,  we use 3× 10^3  ensemble members to construct a new prior by standard EnKF in the first stage. The new prior is constructed by incorporating data information from the  first data assimilation step. 9 local multiscale basis functions (M_i=9) are chosen  in constructing the matrix R to construct coarse model. Then, we construct the surrogate model after obtaining the new prior ensemble in the second stage. The data assimilation process begins from I=1 in second stage. The same local multiscale basis functions are used  in constructing the matrix R to construct the gPC surrogate model. When the order of gPC is set as N_0=7, the number of random samplers is 792 when computing vectors b and A in Section <ref>. In the second  stage of the proposed EnKF, the number of ensemble members is also set as 3× 10^3.





We only know the problem is a diffusion-wave equation, which may be  the subdiffusion or  superdiffusion. Thus, it is important  to identify  the forward model before estimating   the source locations.
The result of identifying model is shown in Figure <ref>, where we count the number of ensemble members applied to the two diffusion models.
By the figure, we can find that the new prior has actually identified  the forward model, where I=0 denotes the process of generating a new prior. In the second stage of two-stage EnKF, the superdiffusion model is almost not used. We define the relative error ε corresponding to the posterior distribution by

    ε := E()-^tr/^tr,

where ^tr=(γ^tr, χ_1^tr, χ_2^tr). Then we compute  the relative error of the final ensemble mean and get  ε=0.08.
Figure <ref> shows all of the one and two-dimensional posterior marginal of ,  where there exist some correlation  between the source locations and γ, such as χ_1^2 and χ_2^2, but the others appear  uncorrelated and mutually independent based on the shape of their 2-D marginal.



 To make prediction, we plot the 95% credible interval and predictive  for the model response with the ensemble members and realizations  by the surrogate model  in Figure <ref> and Figure <ref>. We note the realizations of model response for initial prior is constructed by GMsFE model. The prediction intervals are constructed by taking account of measurement errors. The credible interval and prediction interval, along with the truth, and observation data, are illustrated for u((x, 0); t) and u((0.9, y); t) in Figure <ref> and Figure <ref>, respectively. We note that  the credible intervals are tight in the final data assimilation step. Then we can observe from these figures that the uncertainty associated with both the model fit and predictions decreases with respect to assimilation step. We also plot the marginal density for  the initial prior, new prior and posteriors together in  Figure <ref>, which clearly shows that the posterior is dynamically  updated. From this figure, we note that the new prior identifies the model belonging to subdiffusion.










 §.§ Estimate  a permeability field in a hierarchical model


In this subsection, we consider the fractional superdiffusion equation with homogeneous Dirichlet boundary condition.
The source term is set as f=20, and the end time is set as T=0.11. In the inversion model, both the permeability filed and observation noise are unknown and we want to estimate them
by the two-stage EnKF.  To parameterize the permeability field,  we set l_x=0.2, l_y=0.3, σ^2=1, and the covariance function is assumed to be
  
    C(x,y)=σ^2exp(-x_1-x_2^2/2l_x^2-y_1-y_2^2/2l_y^2).

 We assume that log k(x,ω)  can be represented by the following truncate the KLE

    h(x,ω):=log k(x,ω)= E[h(x,ω)]+∑_i=1^N √(λ)_iθ_i(ω)φ_i(x),    N=20,

where the E[h(x,ω)] is plotted in Figure <ref>.



The ground true parameter ^* are randomly drawn from the standard multivariate normal distribution and the truth permeability map is depicted in Figure <ref> (left). The total number of  data assimilation steps is set as I_2=9. Measurements are taken at time instances 0.012+0.01I:0.002:0.018+0.01I in each data assimilation step, where I∈{1,2,⋯,9} and the locations are distributed on the uniform 5× 5 grid of the domain [0.1, 0.9]× [0.1, 0.9] as shown in Figure <ref> (right).



The full order forward model is defined  on 100×100 uniform fine grid, and GMsFEM is implemented on 5× 5 coarse grid. We choose the number of eigenfunctions selected in calculating the snapshot space as M^i_snap=20. GMsFEM with 10 offline local multiscale basis functions (i.e., M_c=10) is used to solve optimization problem (<ref>).

In the numerical simulation, we use 3× 10^3  ensemble members to construct a new prior by standard EnKF in the first stage, where the first two levels of data information is used.
Three online  local multiscale basis functions (M_i=3) are selected in constructing the matrix R to construct coarse model to obtain the new prior. Then, we construct the surrogate model after obtaining the new prior ensemble in the second stage. Five online  local multiscale basis functions ( M_i=5) are selected  in constructing the multiscale basis  matrix  to construct the gPC surrogate model. When the order of gPC is set as N_0=3, we use 1771 samplers  when computing vectors b and A in Section <ref>. In this stage, the number of ensemble members is set as 10^4. By Algorithm <ref>, we need to simultaneously update the hyperparameter σ^2. We take n_s=0.05 in in Algorithm <ref>.  The  σ^2 for true measurement noise is set as 10^-4. The  posterior mean and posterior variance via  the data assimilation step are illustrated in Figure <ref>, which shows that the uncertainty mainly lies on high-contrast part and boundary.

 

 



Figure <ref> depicts the discrepancy principle when σ^2 is unknown in the forecast step. As we have mentioned before in Subsection <ref>, the measurement noise is the additional Gaussian type. Then we can see that  the expectation of the discrepancy tends to n_dσ^2 as assimilation time moves on.
To measure the estimate accuracy, we  define the relative errors  corresponding to the posterior distribution  by

    ε_k:= k(x,)-k(x,^*)/k(x,^*),
        ε_σ^2 := σ^2-σ_t^2/σ_t^2,
        ε_θ:= -/,

where  and σ^2 are the mean of final assimilation step by the two-stage EnKF and σ_t^2 is the truth observation noise. We just consider the error of the second stage in the two-stage EnKF method. The relative error is plotted in the Figure <ref>. As expected, the relative error gradually decreases with the assimilation steps.




The credible interval and prediction interval, along with the true, and observation data, are illustrated for u((x,0.5); t) and u((0,5, y); t) in Figure <ref> and Figure <ref>, respectively.  We find  that  the final credible intervals are tight and that the uncertainty associated with both the model fit and predictions decreases with respect to data assimilation step. The corresponding uncertainty decreases as x and y get closer to boundary, which is due to the deterministic Dirichlet boundary condition.
The marginal densities  for the unknown parameters in different assimilation steps are  plotted in Figure <ref>,  which shows the the support of density distribution becomes narrower as more data information is used in the posterior exploration.












§ CONCLUSION


We presented a two-stage EnKF  using coarse GMsFEM models in the paper. In the first stage, we used a very coarse GMsFEM model  to approximate the forward model, and the corresponding misfit-to-observed  problem is solved by using GMsFEM. We then constructed the new prior based on the ensembles by traditional EnKF, and then  used  GMsFEM and gPC to obtain a compact  representation for the model response based on the new prior. The two-stage EnKF  was employed to explore the surrogate posterior density, which was incorporated by the surrogate likelihood and the updated prior. It showed that the proposed method leads to the approximate posterior with better efficiency and accuracy  than traditional EnKF method.

The deterministic and statistical methods were combined  together to solve inverse problems.  We obtained not only the point estimate and  confidence interval but also the statistical properties of the unknowns.    A new prior was constructed for Bayesian inference  using  coarse GMsFEM models.  The new prior  contains the significant  region or support of the posterior and  is incorporated with the likelihood to be explored. For non-Gaussian models, we presented a two-stage EnKF using normal score transformation.



L. Jiang acknowledges the support of Chinese NSF 11471107.


99



anb
A. N. Bondarenko and D. S. Ivaschenko, Numerical methods for solving inverse problems for time fractional diffusion equation with variable coefficient, J. Inverse Ill-Posed Probl., 17 (2009), 419-440.
SD
S. Duane, A. D. Kennedy, B. J. Pendleton and D. Roweth, Hybrid monte carlo, Physics letters B, 195 (1987), pp. 216-222.
efendiev2005efficient
Y. Efendiev, A. Datta-Gupta and V. Ginting, X. Ma and B. Mallick, An efficient two-stage Markov chain Monte Carlo method for dynamic data integration, Water Resources Research, 41 (2005), DOI: 10.1029/2004WR003764.
YE
Y. Efendiev, J. Galvis and T. Y. Hou, Generalized multiscale finite element methods (GMsFEM), J. Comput. Phys., 251 (2013), pp. 116-135.
Em
A. A. Emerick and A. C. Reynolds, Ensemble smoother with multiple data assimilation, Computers & Geosciences, 55 (2013), pp. 3-15.
OGE
O. G. Ernst, B. Sprungk and H. J. Starkloff, Bayesian inverse problems and Kalman filters in Extraction of Quantifiable Information from Complex Systems, Lect. Notes Comput. Sci. Eng.,
Springer, 2014, pp. 133-159.

ogh15
O. G. Ernst, B. Sprungk and H. J. Starkloff, Analysis of the ensemble and polynomial chaos Kalm filters in Bayesian inverse problems, SIAM/ASA J. Uncertainty Quantification, 3 (2015), pp. 823-851.


Ge
G. Evensen, Sequential data assimilation with a nonlinear quasi-geostrophic model ssing Monte-Carlo methods to forecast error statistics, Journal of Geophysical Research: Oceans, 99 (1994), pp. 10143-10162.
AG
A. Gelman, J. B. Carlin, H. S. Stern and D. B. Dunson, Bayesian Data Analysis, Chapman and Hall/CRC, 3rd dition,  2013.
Gil
W. R. Gilks, S. Richardson and D. J. Spiegelhalter, Markov chain Monte Carlo in practice, Chapman & Hall, London, 1996.


Han
P. C. Hansen, Analysis of discrete ill-posed problems by means of the L-curve, SIAM Rev., 34 (1992), pp. 561-580.
hou1997multiscale
T. Y. Hou and X. H. Wu, A multiscale finite element method for elliptic problems in composite materials and porous media, J. Comput. Phys., 134 (1997),  pp. 169-189.

ihs13
M. A. Iglesias, K. Law and A. M. Stuart, Ensemble Kalman methods for inverse problems. Inverse Problems 29 (2013),  045001.

jl17
L. Jiang and Q. Li, Model's sparse representation based on reduced mixed GMsFE basis methods,  J. Comput. Phys.,  338 (2017), pp. 285-312.


jo16
L. Jiang and N. Ou,  Multiscale model reduction method for Bayesian inverse problems of subsurface flow, J. Comput. Appl. Math., 319 (2017), pp. 188-209.


jo18
L. Jiang and N. Ou,  Bayesian inference using intermediate distribution based on coarse multiscale model for time fractional diffusion equations, to appear in Multiscale Modeling and Simulation.




kaipio2006statistical
J. Kaipio and E. Somersalo, Statistical and computational inverse problems, Springer Science & Business Media, 2006.



Lc
C. Li and F.Zeng, Numerical methods for fractional calculus, Chapman & Hall/CRC,  2015.

lz
C. Li, Z. Zhao and Y. Q. Chen, Numerical approximation of nonlinear fractional differential equations with subdiffusion and superdiffusion, Comput. Math. Appl., 62 (2011), pp. 855-875.
GLi
G. Li, W. Gu and X. Jia, Numerical inversions for space-dependent diffusion coefficient in the time fractional diffusion equation,  J. Inverse Ill-Posed Probl., 20 (2012), pp. 339-366.



lm14
J. Li and Y. M. Marzouk, Adapptive construction of surrogates for the Bayesian solutions of inverse problems, SIAM J. Sci. Comput., 36 (2014), A1163-A1186.



JG
J. S. Liu and R. Chen, Sequential Monte Carlo methods for dynamic systems,  J. Amer. Statist. Assoc., 93 (1998), pp. 1032-1044.



XI
X. Luo and I. M. Moroz, Ensemble Kalman filter with the unscented transform, Phys. D, 238 (2009), pp.549-562.
SMA
S. Mallat, A wavelet tour of signal processing: the sparse way, Academic press, 2008.





martin2012stochastic
J. Martin, L. C. Wilcox, C. Burstedde and O. Ghattas, A stochastic Newton MCMC method for large-scale statistical inverse problems with application to seismic inversion, SIAM J. Sci. Comput., 34 (2012), pp. A1460-A1487.
LMi
L. Miller and M. Yamamoto, Coefficient inverse problem for a fractional diffusion equation, Inverse Problems, 29 (2013), 075013.
MKS
K. S. Miller and B. Ross, An introduction to the fractional calculus and fractional differential equations, Wiley-Intersci. Publ., 1st edition,  1993.
OK
K. Oldham and J. Spanier, The fractional calculus theory and applications of differentiation and integration to arbitrary order, Elsevier Science,  1974.

Pres
P. Restrepo, O. Rakovec and S. J. Noh, Advancing data assimilation in operational hydrologic forecasting: progresses, challenges, and emerging opportunities, Hydrology and Earth System Sciences, 16 (2012).


ASV
A. Safdari-Vaighani, A. Heryudono and E. Larsson, A radial basis function partition of unity collocation method for convection-diffusion equations arising in financial applications, J. Sci. Comput., 64 (2015), pp. 341-367.

ss17
C. Schillings and A. M. Stuart, Analysis of the ensemble Kalman filter for inverse problems,  SIAM J. Numer. Anal. 55 (2017), pp.  1264-1290.


AS
A. Seiler, G. Evensen and J. A. Skjervheim, Using the EnKF for history matching and uncertainty quantification of complex reservoir models, Computational Methods for Large Scale Inverse Problems and Quantification of Uncertainty, 2011, pp. 247-271.

stu10
A. M. Stuart. Inverse problems: A Bayesian perspective, Acta Numer., 19 (2010), pp. 451-559.

CR
C. R. Vogel, Computational methods for inverse problems, SIAM, Philadelphia, PA,  2002.
WW
W. Wyss, The fractional diffusion equation, J. Math. Phys., 27 (1986), pp. 2782-2785.


xk02
D. Xiu and G. Karniadakis, The Wiener-Askey polynomial chaos for the stochastic differential equations, SIAM J. Sci. Comput., 24 (2002), pp. 619–644.


yan2015stochastic
L. Yan and L. Guo, Stochastic Collocation Algorithms Using l_1-Minimization for Bayesian Solution of Inverse Problems, SIAM J. Sci. Comput., 37 (2015), pp. A1410-A1435.

HZ
H. Zhou, J. J. Gomez-Hernandez and H. J. H. Franssen, An approach to handling non-Gaussianity of parameters and state variables in ensemble Kalman filtering, Adv. in Water Res., 34 (2011), pp. 844-864.

