Solving Multiphysics, Multiparameter, Multimodal Inverse Problems: An Application to NMR Relaxation in Porous Media

A general and robust Bayesian optimization framework for the extraction of intrinsic physical properties from an integration of pore-scale forward modeling and experimental measurements of macroscopic system responses is developed. The efficiency of the scheme, which utilizes Gaussian process regression, enables the simultaneous extraction of multiple intrinsic physical properties with a minimal number of function evaluations. Here it is applied to nuclear magnetic resonance (NMR) relaxation responses, paving the way for inverse problem approaches to digital rock physics given its general nature. NMR relaxation responses of fluids in porous media may be described by sums of multiexponential decays resulting in a relaxation time distribution. The shape of this distribution is dependent on intrinsic physical system properties, but also effects like diffusion coupling between different relaxation regimes in heterogeneous porous materials. Forward models based on high-resolution images are employed to naturally incorporate structural heterogeneity and diffusive motion without limiting assumptions. Extracting the required multiple intrinsic parameters of the system poses an ill-conditioned multiphysics multiparameter inverse problem where multiple scales are covered by the underlying microstructure. Exploration of the multidimensional search space given an expensive cost function makes an efficient solution strategy mandatory. We propose a workflow to match experimental measurements with simulations via Bayesian optimization, with special attention paid to the multimodal nature of the topography of the objective function using solution space partitioning. A multimodal search strategy using state-of-the-art evolutionary algorithms and gradient-based optimization algorithms guarantees that the multimodal nature is captured. The workflow is demonstrated on $T_2$ relaxation responses of Bentheimer sandstone, extracting three physical parameters simultaneously: the surface relaxivity of quartz grains, the effective transverse relaxation time, and the effective diffusion coefficient in clay regions. Multiple mathematically sound and physically plausible solutions corresponding to global minimum and multiple local minima of the objective function are identified within a limited number of function evaluations. Importantly, the shape of the experimental $T_2$ distribution is recovered almost perfectly, enabling the use of classical interpretation techniques and local analysis of responses based on numerical simulation.

Figure 1

The detailed workflow to solve an inverse problem of finding feasible combinations of physical properties of saturated sandstone resulting in NMR $T_2$ distributions identical to the experiment.

Figure 2

Cross section through a ${1120}^3$ voxel subdomain of Bentheimer sandstone; (a) tomogram and (b) phase segmented image with resolution $ϵ=2.89\mu \mathrm{m}$ (blue represents the pore space, red represents the quartz, green represents the clay region).

Figure 3

Illustration of GPR (main axis) and EI (secondary axis) in a minimization problem. EI is high where the predictive value is low (exploitation) and where the uncertainty is high (exploration). The four peaks on the EI plot indicate the four local maxima of EI. How EI is calculated at the four locations is also displayed: the horizontal red line indicates the minimum value at the current iteration and the purple line indicates the probability density distribution of the predictive mean given $\mathbf{x}$, $\mu$, and $\sigma$. The shaded purple region indicates the distribution of the possible improvement at that location, the integral of which is precisely EI.

Figure 4

Illustration of the subroutine searching for promising candidates by partitioning the solution space. GSP and LSP indicate the global solution partition and local solution partition, respectively, UP indicates the unique partition, _PSO indicates the solution is obtained using SL PSO, _R indicates the result is refined by L-BFGS-B.

Figure 5

Illustration of two consecutive steps using the ISW to find global and local minima of the designed one-dimensional objective function. GP mean and GP variance are calculated using Eq. (18). At each step a new candidate is proposed by maximization of EI (green line, right axis).

Figure 6

In step 3, three UPs and the corresponding minima, i.e., local minimum 1 (LM 1) in UP 1 (LSP 1), global minimum (GM) in UP 2 (GSP), and local minimum 3 (LM 3) in UP 3 (LSP 2), colored orange, green, and brown, respectively, are identified using SSP. Of all observations, LMs are colored in denser colors. LM 2 will be proposed for evaluation in the next iteration. Hyperparameters and the GPR proxy are updated whenever a new candidate is evaluated.

Figure 7

Isosurfaces of the parameter space. (a) predictive GP mean, step 67; (b) predictive GP variance, step 67; (c) GP integrated EI, step 67; (d) predictive GP mean, step 114; (e) predictive GP variance, step 114; (f) GP integrated EI, step 114. Five levels of isosurfaces are displayed, corresponding to the minimum (mean and variance) or maximum (integrated EI) 0.01%, 0.1%, 1%, 3%, and 5% of the values of that scalar field. The absolute values corresponding to these isosurface levels are shown on the colorbar. The correlations among physical parameters evaluated at the best candidate are shown in the ${\rho }_q$-${\log }_{10}T_{2e,c}$, ${\rho }_q$-${\log }_{10}{D}_{e,c}$, and ${\log }_{10}T_{2e,c}$-${\log }_{10}{D}_{e,c}$ planes. The contour plot and the isosurface plot share the same colorbar.

Figure 8

Trace of the fitness value and minimum fitness value for the decay of interest at $\mathrm{SNR}=100$. Candidates with low fitness values are identified at step 67 (rankings first) and step 114 (rankings fourth). The horizontal axis starts at 13 since the first 12 candidates are preevaluated without using the ISW to avoid a cold start.

Figure 9

The comparison of the correlation plots between each pair of hyperparameters at the $67\mathrm{th}$ and $114\mathrm{th}$ steps. The panels along the diagonal show histograms of 20 000 samples for the six hyperparameters $\mathit{\theta }$ drawn using Bayesian MCMC after a burn-in period, and the off-diagonal panels display correlations between each pair of hyperparameters.

Figure 10

Normal distribution fitted to the marginal distribution of all six hyperparameters inferred at step 120. The three length scales are plotted together for easier comparison of the kurtosis. Either horizontal axes or vertical axes are displayed in the same scale and the histograms are normalized so that the integral of the bar area is 1. The mean, standard deviation (SD), and expectation for each hyperparameter are given in Table 2.

Figure 11

Distribution of the two identified UPs at step 120 with the top three solutions from each UP. Five isosurface levels are displayed, corresponding to the minimum 0.01%, 0.1%, 1%, 3%, and 5% of the GP mean. Colored circles and crosses shown in the $\rho$-$T_2$, $\rho$-$D_e$, and $T_2$-$D_e$ planes are projections of the candidates. Details of plotted candidates are shown in Table 3.

Figure 12

The comparison of fit for the two LMs from the two UPs (see also Table 3).

Figure 13

The comparison of fit of the four LMs out of 2431 candidates (see also Table 4).

Figure 14

Illustration of (a) the mean $T_2$ distributions at seven SNRs and (b) an enlarged view of the short-time peak. Within each SNR group, the solid line and shaded area represent the mean value and confidence interval, respectively, averaged over 2 to 16 distributions, indicated in Table 5.

Figure 15

Illustration of the representative fits in each SNR group. Simulated $T_2$ distributions are colored according to which generalized unique partition (GUP) they belong to. Fitness values and identified parameter values corresponding to each $T_2$ distribution are listed in Table 6.

Figure 16

Illustration of the decays on both the linear scale and log scale for low-SNR (50) and high-SNR (400) groups. Decays are colored according to which GUP they belong to.

Figure 17

(a) Illustration of the spatial distribution of the 139 UPs and LMs identified from 63 decays followed by categorization of all UPs into six GUPs. Potential correlations between pairs of physical parameters are shown as projections. (b) The strongest correlation is found between relaxation and diffusion in the clay effective phase. Categorization of the six GUPs is based on the value of ${\log }_{10}{D}_{e,c}$, which can be found in Table 7. Here $\mathrm{SNR}=50$ ($\bullet$), 71 ($\star$), 100 ($\ast$), 141 (+), 200 ($\times$), 282 ($◂$), and 400 ($\diamond$) are collectively demonstrated. The 139 LMs are colored according to which GUP they belong to.

Figure 18

Minimum fitness values across seven SNR groups. Confidence intervals for the first three SNR groups with 16 decays are displayed as shaded areas. The first 12 steps are not displayed (same fixed preevaluated candidates without using the ISW). The zig-zag shapes are due to a limited number of realizations of decay in high-SNR groups.

Figure 19

The results of comparisons of fitness values and identified physical parameters of the four local minima at seven SNR values. Each plot shows four groups of bars, corresponding to the four GUPs that passed the physical validity test (GUP 1 and GUP 6 are not shown). Mean and standard deviations are shown for SNR groups with 16 decays.

Figure 20

(a) Overfit and (b) underfit regression for the same dataset as Fig. 5. The hyperparameter settings are shown in Table 9.

Figure 21

The mean fitness values for five ISW variants. Means and standard deviations of the inverse problem using the ISW are displayed, repeated 128 times but using different random streams.

Figure 22

The comparisons of fitness values and identified physical parameter values in the four physically valid GUPs for five ISW variants. Each plot shows four groups of bars (one for each GUP). Means and standard deviations for the 128 repetitions are displayed.