Adaptive Strategies For Efficient Model Reduction In High-Dimensional Inverse Problems

TL;DR

This paper introduces an adaptive, differentiable PCA-based method for efficient model reduction in high-dimensional inverse problems, enhancing optimization processes like history matching.

Contribution

It presents a novel adaptive PCA-based parametrization technique that adjusts based on objective function behavior, improving model reduction in large-scale inverse problems.

Findings

Enables efficient model reduction in high-dimensional inverse problems.

Compatible with gradient-based optimization algorithms.

Improves accuracy of history matching processes.

Abstract

This work explores a novel approach for adaptive, differentiable parametrization of large-scale non-stationary random fields. Coupled with any gradient-based algorithm, the method can be applied to variety of optimization problems, including history matching. The developed technique is based on principal component analysis (PCA), but, in contrast to other PCA-based methods, allows to amend parametrization process regarding objective function behaviour.

Figures

Equations

$$\mathbf{m^{*}}=\operatorname*{argmin}_{\mathbf{m}}\big{[}S(\mathbf{m})=(g(\mathbf{m})-\mathbf{d_{obs}})^{T}\mathbf{C}_{D}^{-1}(g(\mathbf{m})-\mathbf{d_{obs}})\big{]}$$

$$S(\mathbf{m})=\beta (g(\mathbf{m})-{\mathbf{d}}_{\mathbf{obs}}{)}^T$$

$$+(1-\beta )(\mathbf{m}-{\mathbf{m}}_{\mathbf{prior}}{)}^T{\mathbf{C}}_M^{-1}(\mathbf{m}-{\mathbf{m}}_{\mathbf{prior}})$$

$$\displaystyle\min_{\mathbf{m}}\Big{[}S(\mathbf{m})=\sum_{n=0}^{N-1}L^{n}(\mathbf{x}^{n+1},\mathbf{m})\Big{]}\qquad\forall n\in(0,...,N-1)$$

$$g^n({\mathbf{x}}^{n+1},{\mathbf{x}}^n,\mathbf{m})=0,\forall n\in (0,...,N-1)$$

$$\displaystyle\mathbf{\lambda}^{n}=-\Big{[}\frac{\partial L^{n-1}}{\partial\mathbf{x}^{n}}+\mathbf{\lambda}^{n+1}\frac{\partial g^{n}}{\partial\mathbf{x}^{n}}\Big{]}\Big{[}\frac{\partial g^{n-1}}{\partial\mathbf{x}^{n}}\Big{]}^{-1}\qquad\forall n\in(0,...,N-1)$$

$$\displaystyle\mathbf{\lambda}^{N}=-\Big{[}\frac{\partial L^{N-1}}{\partial\mathbf{x}^{N}}\Big{]}\Big{[}\frac{\partial g^{N-1}}{\partial\mathbf{x}^{N}}\Big{]}^{-1}$$

$$\frac{dS}{d\mathbf{m}}=\sum_{n=0}^{N-1}\Big{[}\frac{\partial L^{n}}{\partial\mathbf{m}}+\mathbf{\lambda}^{n+1}\frac{\partial g^{n}}{\partial\mathbf{m}}\Big{]}$$

$$\mathbf{m}=\sum _{k=1}^{\infty }{A}_k{\mathbf{\phi }}_{\mathbf{k}}$$

$$\tilde{\mathbf{m}}=\sum _{k=1}^N{A}_k{\mathbf{\phi }}_{\mathbf{k}}$$

$$\underset{{\mathbf{\phi }}_{\mathbf{k}}}{\min }\text{E}{\int }_r(\delta \mathbf{m}(\mathbf{r}){)}^{2}d\mathbf{r}$$

$$\delta \mathbf{m}(\mathbf{r})=\mathbf{m}(\mathbf{r})-\tilde{\mathbf{m}}(\mathbf{r})=\sum _{k=1}^{\infty }{A}_k{\mathbf{\phi }}_{\mathbf{k}}(\mathbf{r})-\sum _{k=1}^N{A}_k{\mathbf{\phi }}_{\mathbf{k}}(\mathbf{r})=\sum _{k=N+1}^{\infty }{A}_k{\mathbf{\phi }}_{\mathbf{k}}(\mathbf{r})$$

$$A_k={\int }_{\mathbf{r}}\mathbf{m}(\mathbf{r}){\mathbf{\phi }}_{\mathbf{k}}(\mathbf{r})d\mathbf{r}$$

$$\text{E}{\int }_r(\delta \mathbf{m}(\mathbf{r}){)}^{2}d\mathbf{r}=$$

$$=$$

$$\mathbf{K}(\mathbf{q},\mathbf{m})=\text{E}[\mathbf{m}(\mathbf{s}),\mathbf{m}(\mathbf{q})]$$

$${\int }_q\mathbf{K}(\mathbf{s},\mathbf{q}){\phi }_{\mathbf{i}}(\mathbf{q})\mathbf{dq}={\beta }_i{\mathbf{\phi }}_{\mathbf{i}}(\mathbf{s})$$

$$\mathbf{m}\to \mathbf{m}(\xi )=\bar{\mathbf{m}}+\sum _{k=1}^N{\mathbf{\phi }}_{\mathbf{k}}(\mathbf{x}){\mathbf{\xi }}_{\mathbf{i}}=\bar{\mathbf{m}}+W\Sigma \mathbf{\xi }$$

$$W\Sigma V^{\ast }=\text{SVD}(K)$$

$${\text{Energy}}_i=\frac{{\beta }_i}{{\sum }_i{\beta }_i}$$

$$\frac{dS}{d\boldsymbol{\xi}}=\frac{dS}{d\mathbf{m}}\frac{d\mathbf{m}}{d\boldsymbol{\xi}}=\Bigg{[}\sum_{n=0}^{N-1}\Big{[}\frac{\partial L^{n}}{\partial\mathbf{m}}+\boldsymbol{\lambda}^{n+1}\frac{\partial g^{n}}{\partial\mathbf{m}}\Big{]}\Bigg{]}(W\Sigma)^{T}$$

$$\min_{\varphi_{k}}\textbf{E}\Big{[}(\delta S)^{2}\Big{]}$$

$$\delta S(m)$$

$$=S(\tilde{m}+\delta m)-S(\tilde{m})$$

$$\delta S=\frac{\partial S}{\partial m(r)}{|}_{\tilde{m}}\delta m$$

$$\delta S=\sum _{k=N+1}^{\infty }{A}_k{\int }_r\frac{\partial S}{\partial m(r)}{\phi }_k(r)dr$$

$$(\delta S{)}^2=\sum _{i=N+1}^{\infty }\sum _{j=N+1}^{\infty }{A}_i(\omega )A_j(\omega ){\int }_{r_1}{\int }_{r_2}J(r_1)J(r_2){\phi }_i(r_1){\phi }_j(r_2)d{r}_{1}d{r}_2$$

$$(\delta S{)}^2=\sum _{i=N+1}^{\infty }\sum _{j=N+1}^{\infty }$$

$$\times {\int }_{r_1}{\int }_{r_2}J(r_1)J(r_2){\phi }_i(r_1){\phi }_j(r_2)d{r}_{1}d{r}_2$$

$$\displaystyle\textbf{E}\big{[}(\delta S)^{2}\big{]}=\sum_{i=N+1}^{\infty}\sum_{j=N+1}^{\infty}\int_{s}$$

$$\times {\int }_{r_1}{\int }_{r_2}J(r_1)J(r_2){\phi }_i(r_1){\phi }_j(r_2)d{r}_{1}d{r}_2$$

$$J(r)=\sum _{k=1}^{\infty }{c}_k{\phi }_k(r)\to c_k={\int }_{r}J(r){\phi }_k(r)dr$$

$$\underset{\phi }{\min }\sum _{i,j}{c}_i{c}_j{\int }_s{\int }_{q}K(s,q){\phi }_i(s){\phi }_j(q)dsdq$$

Full text

To define an efficient basis update, Adaptive Strategies PCA (AS-PCA) uses gradients of objective function with respect to model parameters. Gradients are already available in gradient-based optimization process combined with adjoint method. Optimality, correctness and low computational cost of the new parametrization procedure is guaranteed by theoretical derivation of basis with stationary perturbation theory from quantum mechanics. Three modifications of method are presented. One of them not only improves quality of parametrization, but also extends applicability of method for uncertainty quantification.

The AS-PCA is then applied to synthetic history matching. Forward problem is represented by non-linear parabolic equation. Through simple cases, we demonstrate an improve in matching performance in comparison with standard linear PCA for misfit minimization and overall field consistency.

1 Introduction

Inverse problems appear in many areas of comparative research, where the problem of defining inner structure of uncertain object having a set of its lifecycle observations is considered. Although correct solution for a lot of practical inverse problems requires an efficient parametrization algorithm, our specific interest here is in the field of history matching problem. The purpose of this procedure is to generate detailed reservoir description that is consistent with prior information and match production data to within some tolerance. This is usually done using two types of data, namely static and dynamic. Static data is mostly constant over time, e.g. geological concept of formation, well logs and petrophysical data and is commonly given as prior information. Dynamic data is time-dependent and represents properties that change during production process, e.g. pressure and flow rates, flow responses. Initial reservoir model is constructed with static data and then has to be adjusted to correctly match dynamic data. This is done by iterative change of model parameters (e.g. facies or permeability) untill model flow simulation is in well agreement with production data. A relatively recent review on history matching problem can be found in [6].

A common approach is to perform history matching either within optimization formulation, or as data assimilation problem. For the latter, ensemble-based methods, such as ensemble Kalman filters (EnKF) [1], recently have gained popularity. Such methods are non-invasive regarding forward flow simulator and provide multiple matched models. These features simplify process of uncertainty quantification and allow to work with black-box simulator. However, among limitations of ensemble-based methods is problem of ensemble collapse. This leads to need of large number of realizations within ensemble and, in hence, high computational costs. Also, these methods are incapable to represent models with complex geological structures.

In optimization context, the history matching problem is often addressed by stochastic methods such as genetic algorithm [8], particle swarm optimization [2], evolutionary algorithms [12] and others. Although these methods perform a global search and allow to use forward simulator as black-box, their use for complex models has a relatively high computational cost and can be applicable only via decent computational clusters.

A gradient-based optimization approach is in the main focus of this paper. These methods are invasive w.r.t the forward simulator and provide a local search, but they are sufficiently faster than methods described above. Adjoint-based techniques are usually applied in history matching procedure in this context, since they provide neccesary gradients at cost of one additional simulation. These techniques are investigated for partial differential equations (PDE) based systems, e.g. with applications in closed-loop reservoir management [11] and, recently, for integro-differential equation (IDE) systems with a memory effect [4].

Since history matching often has to be performed on the real large-scale fields, parametrization techniques are very useful as they substantially reduce the number of parameters that have to be determined. Also, it allows for maintaining geological consistency of result, to some extent. Approaches like discrete cosine transform (DCT) [3] or discrete wavelet transform (DWT) [9] allow to represent model in term of relatively few parameters, but their performance on complex models is not suitable since information about prior covariance of geological model is not considered in basis construction process. Another choice is to use PCA-based algorithms. Classic linear PCA, often referred to as Karhunen-Loeve expansion in this context, has been succesfully applied to history matching problem [11]. Due to the fact that linear PCA considers only covariance matrix of model, it preserves only two-point geostatistics, and number of approaches were developed to use PCA for complex non-Gaussian fields. Among them are kernel PCA (kPCA) [10], optimization-based PCA (O-PCA) [13], regularized kPCA (R-kPCA) [14], convolutional neural network PCA (CNN-PCA) [5].

One of the missing key in existing methods is the lack of ability to adaptively include information about objective function structure into parametrization process. Sensitivity of data is considered in a number of parametrization techniques [15], [7]. Though, to our knowledge, these methods were investigated only for cases when forward problem is linear and seem to be hardly applicable for complex flow modelling. A novel technique, Adaptive Strategies PCA, was developed to address these challenges. Since linear PCA represents fields in term of truncated eigen-decomposition of its covariance matrix, AS-PCA offers two ways to adaptively adjust this representation during optimization process. This can be done either by performing rotation of the truncated eigenvectors of covariance matrix or by swaping them by more sensitive ones. To calculate optimal update, gradients of objective functions are required and using of adjoint-procedures within gradient-based optimization is proposed. Under these conditions, the method provides fast and effective reparametrization with relatively low computational cost. AS-PCA can be extended for use with other state of the art PCA-based techniques, though this question is not considered in the paper.

This paper proceed as follows. After describing the basics of gradient-based history matching in Section 2.1 and short introduction into PCA parametrization in Section 2.2, we provide theoretical and practical description of AS-PCA in Section 2.3. In Section 3, detailed results for method performance on synthetic problem are presented. We conclude with a summary and highlights for future research.

2 Gradient-based history matching

2.1 History matching as inverse problem

The inverse modelling as discussed above require finding the set of unknown model parameters m. In history matching, the m typically represented by permeability in grid blocks of reservoir model. Let $\mathbf{d_{obs}}$ represent the production data, e.g. pressure and flow rates. Then, the history matching problem can be written as follows:

[TABLE]

where $(g(\mathbf{m})$ represent forward flow simulation, $\mathbf{C}_{D}^{-1}$ is the covariance of the observation errors.

For real reservoirs the model m is typically large-scale (about $10^{6}-10^{9}$ grid cells) and exceeds amount of observed data. The parameter-estimation problem is ill-posed and results in non-unique estimates, often geologically unrealistic. A common approach to make problem a well-posed one is to regularize objective function $S(\mathbf{m})$ using prior knowledge about model [6]:

[TABLE]

where $\mathbf{m_{prior}}$ is prior model, $\mathbf{C}_{M}^{-1}$ is the covariance matrix of the prior parameter errors. $\beta$ is a scale constant.

2.2 Adjoint-based history matching

To solve the problem (1) using gradient-based methods one need an efficient procedure to calculate the gradient of objective function $S(m)$ with respect to model $m$. Because $S(m)$ consists of $g(m)$ which involves solving system of PDE, the most efficient way to obtain these gradients is to use adjoint method. In [11] adjoint system was defined for PDE-constrained problem. For similarity with referenced work, the problem (1) can be rewritten as:

[TABLE]

where $\mathbf{x}^{n}$ is the PDE solution at n-th timestep, $L^{n}(\mathbf{x}^{n+1},\mathbf{m})$ is the same least squares misfit as in (1).

Following [11], adjoint model can be derived as:

[TABLE]

And finally, gradient of the objective w.r.t. the model can be calculated as:

[TABLE]

Then it can be used with any gradients-based optimization algorithm. In this paper, a nonlinear conjugate gradient (CG) method was applied.

2.3 PCA-based model parametrization

The models, generated in history matching process, have to preserve geological features of formations, i.e. to be geologically consistent. Also, performing history matching with the large scale reservoir models is a high cost process. Principal component analysis parametrization technique was developed to adress this problem. It allows to efficiently reduce model size while preserving key model features known from prior information. PCA utilize idea of stochastic process representation as an finite linear combination of orthogonal functions.

It is well known that any stochastic process can be represented as linear combination of infinite number of orthogonal functions:

[TABLE]

where $\boldsymbol{\varphi_{k}}$ is k-th orthogonal function, $A_{k}$ - k-th coefficient of decomposition

The idea of data decomposition is to define such type of orthogonal basis functions, that model $\hat{m}$ given by the truncated basis approximate initial model with minimal error.

[TABLE]

PCA or Karhunen-Loève expansion yields the best such basis in the sense that it minimizes the total mean squared error:

[TABLE]

where $\delta\mathbf{m}$ is defined as approximation misfit:

[TABLE]

Using the fact that

[TABLE]

Total mean square error from (10) can be expressed as:

[TABLE]

where $K(q,m)$ - is a covariance function defined as:

[TABLE]

Optimal orthogonal basis vector will be given by solving the following problem:

[TABLE]

Equation (15) states that optimal choice for basis vectors is eigenvectors of prior models covariance matrix.

In practice we are interested in generating new models with the structure of initial model. Using KL expansion it can be performed in following manner:

[TABLE]

where $\mathbf{\bar{m}}$ - dataset mean model, $\boldsymbol{\xi}$ - vector from $N(0,1)$, $W$ - transition matrix with eigenvectors in columns, $\Sigma$ - diagonal matrix containing eigenvalues

Practically, $W$ and $\Sigma$ can be defined from SVD decomposition of dataset covariance matrix:

[TABLE]

where $K$ is a covariance matrix

Number of truncated eigenvectors is determined by energy retained in them:

[TABLE]

where $\beta_{i}$ is eigenvalue of corresponding eigenvector

Structure of parametrization by PCA given by (16) allows easily get gradient of objective function with respect to coefficients of decomposition. Combining (7) and (16):

[TABLE]

3 Adaptive parametrization using AS-PCA

3.1 Construction of a new eigenvalue problem

The described above principal component analysis parametrization doesn’t consider sensitivity of objective function w.r.t basis components. Therefore choice of including eigenvectors into optimization is limited by those with highest retained energy (18). This fact leads to an incorrect model parametrization when important geological structure is underrepresented in prior models. Due to high uncertainty of data and inner reservoir structure, the problem can occur in most of real cases. In this section we propose a new technique which allows to easily include information about sensitivity into basis construction process. Our Adaptive Strategies Principal Component Analysis (AS-PCA) is much faster than sensitivity analysis process, since it uses only a few linear transforms on initial basis and doesn’t require any additional run of forward model.

The basis of classic PCA is given as truncated Karhunen-Loeve expansion which provides the best approximation of the original process in the sense that it reduces the total mean-square error resulting of its truncation as described above (10), (15).

Our suggestion is to find a new basis that will minimize the misfit between objective function value of original and approximated model.

[TABLE]

[TABLE]

From (3.1):

[TABLE]

Substituting definition of $\delta x$ from (11):

[TABLE]

For the sake of convenience let’s define $\frac{\partial S}{\partial m(r)}$ as $J(r)$

Let’s then calculate total mean square error:

[TABLE]

Using the same workflow as for classic KL expansion:

[TABLE]

[TABLE]

Then, because gradient $J(r)$ has the same shape as $r$, we can decompose it by the same basis:

[TABLE]

Therefore, we can define minimization problem as:

[TABLE]

To find optimum basis $\varphi_{k}$ we will define augmented with Lagrange multipliers $\beta_{i}$ function:

[TABLE]

[TABLE]

In order to find optimal $\varphi_{k}$ we need to calculate $\frac{\partial L_{A}}{\partial\varphi_{k}}$ and set it to zero:

[TABLE]

[TABLE]

Finally, we have:

[TABLE]

Dividing whole equation by $c_{i}^{2}$ and substituting $\frac{c_{j}}{c_{i}}$ by $\alpha_{ij}$ and $\frac{\beta_{i}}{c_{i}}$ by $\beta_{i}$ we get:

[TABLE]

Without second addend in left side, equation (34) is just the same eigenproblem as classical PCA gives (15) and solution is already known. However, the key difference here is the presence of relatively small addend, form of which doesn’t allow us to easily solve the new eigenvalue problem.

3.2 Stationary perturbation theory

One can notice that equation (34) looks similar to quantum mechanics problem of solving Schrodinger equation, with complex Hamiltonian $\hat{H}$ which doesn’t allow to directly solve problem. Classic approach here is to divide Hamiltonian $\hat{H}=\hat{H}_{0}+\hat{V}$, in such ways that eigenvectors $\varphi_{k}^{0}$ and eigenvalues $E^{0}$ of Hamiltonian $\hat{H}_{0}$ are known, and $\hat{V}$ is a small corrective of $\hat{H}_{0}$ or as said, perturbation and then apply Stationary Perturbation Theory (SPT) approach.

We can rewrite eigenvalue problem (34) in notation that is commonly used in quantum mechanics:

[TABLE]

Known eigenvalues and eigenvectors for classic PCA case we will define as $\beta^{0}$ and $\varphi^{0}_{k}$.

As SPT suggests, let’s decompose solution of (35) by the basis of standard PCA:

[TABLE]

Then, we can rewrite equation (34) as:

[TABLE]

Let’s multiply (37) by $\varphi_{n}^{0}$ and integrate by $ds$. Then we can use that

[TABLE]

And get:

[TABLE]

Now we’re going to use following representation:

[TABLE]

where $c_{k}^{i}\rightarrow c_{i.k}$ for the sake of convenience.

From which we get:

[TABLE]

Because $c_{i,n}^{0}=\delta_{i,n}$, then:

[TABLE]

Finally, we get coefficients of decomposition (36):

[TABLE]

Eigenvectors of the problem (34) are given by:

[TABLE]

Where $\gamma$ is manually setted coefficient that will keep $V$ small relatively to $H$, as requirement for usage of SPT.

[TABLE]

3.3 AS-PCA modifications

A new parametrization basis (44) is given by performing small rotation of the initial basis from classic KL expansion by using information of truncated components. Orthogonality is saved due to the smallness of addend and it allows us to efficiently embed this method into existing PCA + gradient-based optimization code just by substituting old basis by new without significant changes in workflow. Due to the fact that technique performs rotation of original PCA basis, the method was called rotation strategy of AS-PCA.

We find interest in two heuristic modifications of idea explained above. Having total input of each (unused in the old basis) component, we can sort them by impact and then do one of two possible variants:

Extend old basis by adding some number of vectors with the highest impact, or 2. 2.

Substitute weak component of basis by them

The first technique is called extension strategy of AS-PCA

The procedure of substituting vectors has been named ’swap’ and the whole technique is called swap strategy of AS-PCA

4 Results

4.1 Synthetic problem description

In this section the comparative analysis of the classic PCA and proposed method is presented. A 1D non-linear diffusion problem with second-order boundary conditions was chosen as synthetic problem to test the methods:

[TABLE]

By previously used nomenclature:

• •

As a model $m$ - spatial distribution of diffusion coefficient $D(x)$;

• •

As observed data $d_{obs}$ - concentration fields $u_{obs}$, given by solving (4.1) with known $\hat{D}(x)$

Model size is 100 grid cells. Our goal is to define $\hat{D}(x)$ having only $u_{obs}$

The problem (4.1) was discretized by fully implicit scheme using finite volumes method and then solved using Newton-Raphson method.

A prior information was represented by dataset of possible model realizations. Firstly, the true model was defined as:

[TABLE]

Then, dataset was generated as true realization with random smooth perturbations:

The AS-PCA algorithm was tested on two cases. In first, the true model was represented by (47). The second case was constructed by adding low-frequency noise into (47)

The solutions of the problems are presented in figure below. All three strategies of AS-PCA were compared with PCA parametrization and with full model opimization without reduction. Since every basis adaptation provides new solution, we presented all of them in plots corresponding to methods.

Our general conclusions from results described below are summarized as follows:

Rotation Strategy. A small change of basis allows to efficiently investigate area of local minima. Our proposition is to use the strategy when all structures of solution a presumed to be well-known. 2. 2.

Swap Strategy. Because every swap significantly change form of basis, a big amount of relatively different realizations can be obtained. Despite the fact that use of strategy doesn’t always leads to a lower values of objective functions, it provides a number of realization with significantly different structures. Our assumption is to use the strategy for uncertainty quantification and forecasting, though this is topic of further research. 3. 3.

Extension Strategy. On presented synthetic problem and dataset, significant impact of extension strategy wasn’t observed. Also, method increase computational cost of problem solving, due to increasing number of parameters that have to be adjusted.

5 Discussion

A novel algorithm for adaptive parametrization of random fields was developed. This PCA-based method is expanded to include information about objective function structure into basis construction process. This can be done using one of three suggested strategies. One of them provide significant improve in local minima investigation (3), second allows multiple different models obtaining. Although algorithm has been tested on synthetic examples with artificially created model complexities, applicability on more practical models is expected, since real data always have uncertain structures. Also, method provides a relatively fast parametrization, since only a few linear transforms of initial basis are required altogether with already given gradients.

Although, in this paper, use of the method is limited by gradient-based optimization, its implementation simplicity (1, 3) and relatively low computational cost allow to easily improve overall quality of inverse problem solution. Algorithm development for effective usage with non-gradient methods is one of possible directions of further research.

Despite that AS-PCA was considered on the case of history matching problem, it is applicable for a wide class of other inverse problems, such as tomography, computer vision, etc. and can be extended for usage with a variety of parametrization techniques.