Uncertainty quantification for inverse problems with weak partial-differential-equation constraints

In statistical inverse problems, the objective is a complete statistical description of unknown parameters from noisy observations to quantify uncertainties in unknown parameters. We consider inverse problems with partial-differential-equation (PDE) constraints, which are applicable to many seismic problems. Bayesian inference is one of the most widely used approaches to precisely quantify statistics through a posterior distribution, incorporating uncertainties in observed data, modeling kernel, and prior knowledge of parameters. Typically when formulating the posterior distribution, the PDE constraints are required to be exactly satisfied, resulting in a highly nonlinear forward map and a posterior distribution with many local maxima. These drawbacks make it difficult to find an appropriate approximation for the posterior distribution. Another complicating factor is that traditional Markov chain Monte Carlo (MCMC) methods are known to converge slowly for realistically sized problems. To overcome these drawbacks, we relax the PDE constraints by introducing an auxiliary variable, which allows for Gaussian errors in the PDE and yields a bilinear posterior distribution with weak PDE constraints that is more amenable to uncertainty quantification because of its special structure. We determine that for a particular range of variance choices for the PDE misfit term, the new posterior distribution has fewer modes and can be well-approximated by a Gaussian distribution, which can then be sampled in a straightforward manner. Because it is prohibitively expensive to explicitly construct the dense covariance matrix of the Gaussian approximation for problems with more than 105 unknowns, we have developed a method to implicitly construct it, which enables efficient sampling. We apply this framework to 2D seismic inverse problems with 1800 and 92,455 unknown parameters. The results illustrate that our framework can produce comparable statistical quantities with those produced by conventional MCMC-type methods while requiring far fewer PDE solves, which are the main computational bottlenecks in these problems.