Reservoir characterization

Sparse History Matching: Nonlinear-Orthogonal-Matching-Pursuit Algorithm

A nonlinear orthogonal-matching pursuit (NOMP) for sparse calibration of reservoir models has been proposed. Sparse calibration is a challenging problem because the unknowns are the nonzero components of the solution and their associated weights.

jpt-2013-04-hmffocushero.jpg
Source: Getty Images.

A nonlinear orthogonal-matching pursuit (NOMP) for sparse calibration of reservoir models has been proposed. Sparse calibration is a challenging problem because the unknowns are the nonzero components of the solution and their associated weights. NOMP is a greedy algorithm that, at each iteration, discovers components of the basis functions that are most correlated with the residual. The proposed algorithm relies on approximate-gradient estimation by use of an iterative-stochastic-ensemble method (ISEM). ISEM uses an ensemble of directional derivatives to approximate gradients efficiently.

Introduction

Subsurface-flow models rely on many parameters that cannot be measured directly. Instead, a sparse set of measurements may exist at well locations. The complete distributions of these unknown fields commonly are inferred by a model-calibration process that takes into account historical records of the input/output of the model. However, the amount of data available to constrain the models often is limited in quantity and quality. The result is an ill-posed inverse problem that might allow many different solutions. ­Parameter-estimation techniques that can be applied to this problem can be classified into Bayesian methods based on ­Markov-chain Monte Carlo methods, gradient-based-optimization methods, and ensemble-Kalman-filter methods.

An important step in the automatic-calibration process is to define a proper parameterization of these unknown fields. Most of the parameterization methods depend on a previous model’s assumptions that define the spatial correlations of the unknown fields implicitly with a parameter-covariance matrix. Karhunen-Loève expansion (KLE) can be used for parameterizing spatially distributed fields. KLE, also known as proper orthogonal decomposition or principal-component analysis in the ­finite-dimensional case, is widely used for parameterizing the permeability field in subsurface-flow models. KLE is an effective method that is simple to implement; however, it preserves only the second-order moments of the distribution. For complex continuous geological structures such as a channelized domain, KLE fails to preserve higher-order moments.

Sparse calibration and compressed sensing are active research areas in the signal-processing community. Standard reconstruction methods rely on defining a set of basis functions that are orthogonal, as in KLE methods, and then an attempt is made to find the optimal set of weights to reconstruct the measurements. This reconstruction problem is an ill-posed problem, and regularization techniques (i.e., Tikhonov regularization) that constrain the ℓ2 norm of the solution are applied commonly. The quality of the solution depends on the class of basis functions used to parameterize the search space. In sparse-­calibration methods, a large collection of basis functions is included in a dictionary, and the solution process consists of picking the best basis functions for accurate reconstruction of the unknown field and finding the associated weights.

The authors used sparse-­dictionary-learning-based parameterization. Given a set of realizations of the unknown field (e.g., permeability field), the dictionary-learning problem was formulated as an optimization problem to find the best basis functions such that each realization could be represented as a linear combination of only a few basis functions. These dictionaries were overcomplete and had a certain amount of redundancy. This redundancy is desirable because it provides a robust representation. Building the optimal dictionary that approximates a signal with a minimum error is a nondeterministic polynomial-time hard problem, and approximate algorithms can be used. The authors used the K-SVD algorithm for parameterizing the unknown subsurface fields.

Once the dictionary is defined, the sparse calibration can proceed in two different directions. The first direction is to solve an optimization problem that penalizes the solution in the ℓ1 norm and minimizes the reconstruction error. The second class of algorithms is greedy algorithms, which iteratively find and remove elements from the dictionary that are maximally correlated with the residuals. This study used an iteratively-­reweighted-least-squares (IRLS) algorithm to identify the important dictionary elements (solution support) and their associated weights by minimizing a sparsely regularized objective function and then used an adjoint code to estimate sensitivities for solving the nonlinear-­parameter-estimation problem.

The authors developed an ­ensemble-based method for solving the sparse-­calibration problem given a dictionary built by use of the K-SVD algorithm. ­Ensemble-based methods have proved to be an effective tool for subsurface-­model calibration. The proposed algorithm enables the use of sparse-­calibration techniques for computer models when adjoint codes are not available. For the sparse-­calibration problem, a new algorithm based on the orthogonal-­matching-pursuit (OMP) algorithm was proposed. The proposed algorithm is of the class of greedy algorithms for sparse recovery and extends the standard OMP algorithm (limited to linear reconstruction problems) to nonlinear parameter-­estimation problems.

Parameter-Estimation Algorithm

Calibration of subsurface-flow models given a dictionary of basis functions for the unknown fields is a nonlinear ­parameter-estimation problem. Here, an ensemble-based method is used for parameter estimation. The proposed ­parameter-estimation method relies on the Gauss-Newton method and stochastic estimation of the derivatives by use of an ensemble of directional derivatives.

ISEM. Directional derivatives are used in a stochastic ensemble method for parameter estimation. The authors use an ensemble of perturbations to approximate the standard derivative (gradient) from an ensemble of directional derivatives. Eq. 15 in the complete paper is the main update equation of the proposed ISEM. This update equation was used in the nonlinear OMP algorithm for this study.

Linear Sparse Reconstruction

The calibration process is converted into a sequence of linear problems formulated by Eq. 15 in the complete paper. However, distributed parameter fields (e.g., permeability or porosity) commonly are parameterized to obtain efficient calibration methods. The authors used a parameterization that builds a large dictionary of basis functions. These large dictionaries have the advantage of dealing with non-Gaussian models and a mix of models. The sparse-reconstruction problem is applicable for linear problems and is applicable for nonlinear parameter-­estimation problems at the linearized iteration level as formulated by Eq. 15 in the complete paper.

OMP Algorithm. OMP is an iterative greedy algorithm. The algorithm is an extension of basis-pursuit algorithms. The OMP algorithm has been applied to sparse-signal recovery in many studies. The OMP algorithm pseudocode is detailed in Algorithm 1 in the complete paper.

Dictionary Learning. Learning, or building, a dictionary aims to provide a pool of basis functions in which a few basis functions can be combined linearly to approximate a novel signal or field.

Sparse-Nonlinear-Parameter-Estimation Algorithm. The authors present NOMP for sparse calibration of nonlinear models. To solve this problem, a natural extension of the OMP algorithm, NOMP, was used as a greedy algorithm for nonlinear problems by storing the discovered solution support between the subsequent nonlinear iterations. This is consistent with the logic of the OMP as a greedy algorithm. Once an atom of the dictionary is included in the support, it then is carried over all subsequent update iterations. The pseudocode of NOMP for sparse calibration combined with the ISEM is described in Algorithm 2 in the complete paper. Note two major changes of NOMP from the standard OMP algorithm. First, the solution support is carried between the nonlinear iterations. Second, once the solution support is identified, ℓ2 regularization is used for calculating the residuals but the solution space is limited to the identified support. The ℓ2 regularization is needed because the estimated sensitivity matrix is rank deficient and may contain sampling errors.

Problem Formulation

A two-phase immiscible flow in a heterogeneous porous subsurface region is considered. To simplify this example, gravity and capillary effects were neglected. However, the proposed model-­calibration algorithm is independent of the selected physical mechanisms. The two phases were water and oil. This ­subsurface-flow problem is described by the mass-conservation equation and Darcy’s law. The pore space was assumed to be filled with fluid, and thus the sum of the fluid saturations should add up to unity. Therefore, only the water-­saturation equations need to be solved.

K-SVD Parameterization. The reference permeability fields for the test problem represent channelized models. Different realizations of channelized models were generated from the training image by use of geostatistical modeling software. The training image was from a similar published example. A total of 1,000 realizations were generated and used as input to the K-SVD algorithm to produce a sparse parameterization of the search space. Twelve basis functions were selected randomly from a dictionary of 500 basis functions built with the K-SVD algorithm with target sparsity of 20 elements.

In numerical testing, the discretized model used a 2D regular grid of 50×50 blocks in the x- and y-directions, with each gridblock being 10 m. The z-direction had unit thickness. The porosity was assumed to be constant in all gridblocks (0.2). The water viscosity was set to 0.3 cp, and the oil viscosity was set to 3 cp. The irreducible water saturation and irreducible oil saturation were set at 0.2, and the simulations were run until 1 pore volume was injected. Two injection/production patterns were used for the test problem. Fig. 1 shows the injection wells as black dots and the production wells as white dots. The first pattern has one injection well and four production wells arranged in the inverted-five-spot pattern, as shown in Fig. 1a. Pattern 2, shown in Fig. 1b, had nine production wells distributed around four injection wells. For parameter estimation, the production curves at the production wells were used to define the misfit function and guide the inverse-problem solution. Each water-cut curve was sampled at 50 points, and these samples were used for calculating the errors and the update equation. The observation data (water-cut values) were perturbed with uncorrelated white noise with a small standard deviation of 10−6 to enable performing a convergence study for different ensemble sizes.

jpt-2013-04-sparsehistoryf1.jpg
Fig. 1: Injection/production patterns (injector as black dots and producer as white dots), (a) Pattern 1 and (b) Pattern 2.

 

Waterflooding Test Case. For comparison purposes, all runs were initialized using a uninformed prior of a uniform permeability k field with log (k)=0. For injection/production Pattern 1, the optimized permeability fields shown in Fig. 2 resulted from different ensemble sizes of 5, 10, and 20 members. The smallest ensemble, 5 members, managed to reproduce the locations of the two channels running along the model. However, this is not expected from every run of the algorithm because of the approximate nature of the estimated derivatives. Also, the problem is ill posed and may admit different solutions. Fig. 3 shows the corresponding weights of the basis functions selected from the dictionary for ensembles of 5, 10, and 20 members. For an ensemble of 5 members, the inferred support has 73 nonzero bases, as shown in Fig. 3a. The initial root-mean-square error (RMSE) for the initial permeability field of log (k)=0 is 0.0691. The number of nonlinear iterations is set to 40, 20, and 10 for the ensemble size of 5, 10, and 20, respectively. This corresponds to the same number of 200 forward runs. It was observed that smaller ensembles produced solutions with larger support because of the increased number of nonlinear iterations. This was reflected in the final RMSE that is smaller for smaller ensembles compared with larger ensembles after 200 forward runs. The optimized permeability fields for injection Pattern 2 are shown in Fig. 4 for ensembles of 5, 10, and 20 members. The stem graph of the discovered weights is shown in Fig. 5. The initial RMSE for the initial permeability field of log (k)=0 was 0.0729. The number of nonlinear iterations was set to 40, 20, and 10 for the ensemble size of 5, 10, and 20, respectively. Again, it was observed that smaller ensembles were more effective in matching the data because the number of nonlinear iterations was larger for the same number of forward runs.

jpt-2013-04-sparsehistoryf2.jpg
Fig. 2: Calibrated log-permeability field for different ensemble sizes for the waterflooding test case under injection/production Pattern 1. (a) n=5, (b) n=10, and (c) n=20.

 

jpt-2013-04-sparsehistoryf3.jpg
Fig. 3: Stem graphs of the weights of the calibrated permeability fields for the waterflooding test case under injection/production Pattern 1. (a) n=5, (b) n=10, and (c) n=20.

 

jpt-2013-04-sparsehistoryf4.jpg
Fig. 4: Calibrated log-permeability field for different ensemble sizes for the waterflooding test case under injection/production Pattern 2. (a) n=5, (b) n=10, and (c) n=20.

 

jpt-2013-04-sparsehistoryf5.jpg
Fig. 5: Stem graphs of the weights of the calibrated permeability fields for the waterflooding test case under injection/production Pattern 2. (a) n=5, (b) n=10, and (c) n=20.

 

Convergence Study. A complete convergence study was performed for the test cases. The stochastic nature of the estimated gradients resulted in different solution paths for each run. The average of 50 runs showed the ensemble-size effect on the error-reduction rates. Fig. 6 shows the average RMSE in water cut vs. the total number of forward runs under injection/production Pattern 1 (left) and Pattern 2 (right). Smaller ensembles were more effective for sparse calibration in terms of error reduction and support detection for the same number of forward runs. On average, smaller ensembles outperformed larger ensembles in all the numerical test cases because of the increased number of nonlinear iterations for the same number of forward runs. At each nonlinear iteration, the support was updated and the major search directions were detected. Thus, applying more nonlinear iterations had a positive effect on the error reduction.

jpt-2013-04-sparsehistoryf6.jpg
Fig. 6: Average RMSE in water-cut curve vs. the total number of forward runs for different ensemble sizes with different injection/production patterns, (a) injection/production Pattern 1 and (b) injection/production Pattern 2.

Conclusions

The solution of the nonlinear-sparse-­calibration problem is challenging. The algorithm must find the optimal weights to reproduce the measured values, and it must select the basis functions that are included in the solution support. A complete combinatorial exploration by running standard parameter-­estimation algorithms on a subset of the basis functions led to a huge combinatorial problem that was impossible to solve. In the linear setting, different algorithms for sparse reconstruction can be used. These algorithms can be classified simply as forward-stagewise-selection algorithms (such as the OMP algorithm) and as ­optimization-based algorithms (such as the IRLS algorithm).

The calibrated models using the NOMP algorithm did not show extreme values in the inferred permeability fields. This was attributed to applying ℓ2 regularization at each iteration once the solution support was discovered. The ℓ2 regularization has the advantage of penalizing large weights that produce realizations with extreme permeability values. This was evident from the stem plots showing the weights of the different dictionary atoms. This is a clear advantage of NOMP over different sparse-reconstruction algorithms that penalize only the ℓ1-norm of the solution. Another advantage of the proposed algorithm is in the efficient use of an ­ensemble-based approximate derivative using ISEM. The proposed algorithm combining ISEM and NOMP facilitates sparse calibration for numerical-­simulation packages when the adjoint code is not available.

This article, written by Senior Technology Editor Dennis Denney, contains highlights of paper SPE 163582, “An Ensemble-Based Nonlinear-Orthogonal-Matching-Pursuit Algorithm for Sparse History Matching of Reservoir Models,” by Ahmed H. Fsheikh, University of Texas at Austin, King Abdullah University of Science and Technology; Mary F. Wheeler, SPE, University of Texas at Austin; and Ibrahim Hoteit, King Abdullah University of Science and Technology, prepared for the 2013 SPE Reservoir Simulation Symposium, The Woodlands, Texas, 18–20 February. The paper has not been peer reviewed.