Understanding uncertainty in CSEM

Michael E. Glinsky and James Gunning, Commonwealth Scientific and Industrial Research Organization (CSIRO), Australia

Controlled-source electromagnetic (CSEM) remote sensing of petroleum reservoirs has come into increasing use. An understanding of the uncertainty and resolution limitations of this data is important for practical integration into the E&P decision-making process. We have developed tools for examining these issues based on Bayesian statistical methods. These methods have been applied to a range of synthetic models. For 1D common midpoint (CMP) inversions, CSEM data can usually be explained with rather simple models of 10-layer order. The limits of the resolution in depth, resistivity and thickness of those layers were also determined, along with the sometimes significant correlations in the properties.

We identified and overcame major technical hurdles in the inversion problem. The results were analyzed with a novel cluster identification method that identifies the resistive anomalies and separates them from a background trend. We also studied the effects of systematic noise on the result and successfully implemented methods to account for this effect in the inversion. The result is a practical tool for CSEM data analysis that gives a reliable uncertainty for use in business decisions.

INVERSION PROBLEM

CSEM is naturally suited to detecting resistive anomalies in the marine subsurface due to the presence of hydrocarbon deposits, so regional geological knowledge or other remote sensing techniques must be sufficient to exclude the possibility of highly resistive rocks such as evaporites, volcanics or carbonates. Taken in conjunction with seismic data for geological and structural delineation, the tool is potentially a powerful discriminator between high and low gas saturation. The technique is a valuable complement to seismic methods, since these are well known to have difficulty in detection of gas saturation in amplitude-versus-offset (AVO) applications.

Many articles have been published outlining the general nature of the CSEM acquisition framework, Fig. 1.¹  CSEM inversion is a challenging problem,^2,3 not just because the forward energy propagation is heavily dispersive and, thus, destroys resolution, but also because the subsurface response is poorly modeled by a weak scattering (Born) approximation when the conceivable range of resistivity variation of submarine rocks spans several decades. Subsurface resistors behave almost as waveguides, so the overall field response in the presence of resistors is qualitatively very different from typical background response caused by conducting shales—hence the failure of weak scattering theory. The imaging problem is thus controlled by both dispersive physics and strong nonlinearity, which, when combined with expensive forward models, makes the problem qualitatively similar to the difficult history matching problem in petroleum reservoir characterization.

Fig. 1. The data acquisition of CSEM data. Courtesy of Scripps Institute of Oceanography.

One of the main consequences of the waveguide-like energy flow along subsurface resistors, if present, is that received signals are largely controlled by the resistivity-thickness product (RTP) of these resistors. This characteristic can be shown both numerically and analytically in 1D.⁴ Imaging for structure—that is, extracting thickness and depth—is, thus, very challenging in the absence of complementary information, such as seismic delineation of resistor boundaries.

These issues, being germane to the forward physics, affect both regularizing and Bayesian approaches to inversion. But Bayesian approaches enable a probabilistically coherent statement of what inversion degeneracies or multimodalities imply about subsurface structure, and are able to embrace multiple data types with much more logical consistency. This is our preferred approach.

Furthermore, an important goal of modern Bayesian inversion methods is an assessment of the range of uncertainty of important reservoir parameters as a means of informing risk management decisions. Typically, this approach will involve quantities such as gas saturation, reservoir thickness or, more ambitiously, total hydrocarbon in place. This is another level of difficulty over and above trying to compute “maximum-likelihood”-type inversion images for several reasons. First, the strong nonlinearity means the posterior distributions are rarely approximately quadratic (“bowl shaped”) in parameter space, so local-linearization approaches to uncertainty mapping are not useful. This means sampling methods such as Markov-Chain Monte Carlo or parametric bootstrap-type techniques must be used. These are computationally demanding, since very many forward models (MCMC) or many inversions (bootstrap) need to be computed. These methods have been successfully demonstrated for 1D problems.^3,5

A second challenge is that uncertainty estimates are dependent on the noise model to leading order. Simple maximum-likelihood parameter estimates in inverse problems that are not too nonlinear or degenerate are roughly independent of estimated noise levels. Statisticians speak of such parameters as being robustly estimated or “insensitive to noise mis-specification.” But the uncertainty of model parameters is always dependent on the noise model to leading order, even in kind problems. In hard problems such as CSEM, the nonlinearity and near degeneracy make parameter uncertainties rather sensitive to the choice of the noise model.

A pervasive difficulty is that forward models in geophysics virtually always make assumptions that render faster computation possible, and errors resulting from these assumptions are expected to be absorbed in the “effective” noise model. Examples in the CSEM context might be use of isotropic resistivities, neglect of bathymetry or approximation of 3D responses by 2D or 1D models. These kinds of modeling approximations usually lead to strongly correlated modeling errors, and if the modeling error dominates the overall instrumental and environmental noise, the overall “effective” noise process can then be expected to be fairly spatially correlated. The challenge is to provide estimates of parameter uncertainty that are useful in decision making and that are not hopelessly corrupted by inadequacies in noise modeling.

BAYESIAN METHODOLOGY

Bayesian methods ^6,7 start with Bayes’ theorem, which is a relation for conditional probabilities based on two ways of expressing the joint probability of events A and B:

where P(A|B) is the probability of A given B, P(A) is the probability of A and P(A^B) is the probability of A and B. For convenience, we identify A as the geologic model m that we wish to estimate based on the observed data d, identified as B. We can then rewrite this relationship as:

We want to know P(m|d), the probability of the model given the data. P(m) is the declared prior probability of the model, provided by inversions of other types of data or by prior knowledge about the geology. P(d|m) is the probability of the data given the model, or likelihood, which involves the forward computational model and some noise specification. The analogy with Boltzmann statistical physics comes from writing the probability of the data given the model as:

where F(m) is the forward synthetic data of the model. Note the similarity of the mismatch term to the statistical mechanical Hamiltonian H divided by the temperature T. The set of models, {mi} distributed with probability P(m|d) is analogous to a statistical ensemble of realizations with a Hamiltonian dependent on the data mismatch divided by a temperature determined by the noise (error) in the data. This is a classical problem in the estimation of equations of state. There is rich literature on methods to sample this ensemble.8 From such a set of samples, for example, any average property can be calculated as:

The “best” inversion, or maximum a posteriori (MAP) estimate, coincides with the zero-temperature configuration of the model, or global minimum of ||d-F(m)||2. This is what statisticians call a “point estimate.” But there may be several competitive local minima, or the minima may be quite diffuse. It is much better to base decisions on the range of solutions supported by the posterior distribution, which is why generating samples from the posterior is so helpful. Metropolis methods enable this sampling to be done even when F(m) is complicated or exists only as a forward computer model.

TECHNICAL HURDLES

The major challenge of CSEM data is the badly scaled, nonlinear and multi-modal character of the forward model’s null space. The large dynamic range of resistivity in EM problems makes the Born approximation largely useless, unlike problems such as seismic reflection imaging. As an illustration, in Fig. 2 a simple model of an objective layer near resolution limits is split in two. For a seismic reflection problem in impedances, the null space is linear, so a Gaussian prior on each layer’s impedance yields a unimodal posterior.

Fig. 2. a) Model used to understand the character of the solution: a near-resolution target layer buried under about a kilometer of sea and a kilometer of rock. b) Characteristic null space and posterior probability contours for reflection seismic. c) Characteristic null space and posterior probability contours for CSEM. d) Possible slow diffusion problem with Metropolis CSEM sampling. e) Parametric bootstrap methods quickly explore the posterior by repeated optimizations.

For the CSEM case, conditions are quite different. The null space is very nonlinear, and the posterior surface has very steep-sided curving valleys, joining multiple maxima. This creates several problems in the solution. How are the multiple optima identified? How is the step size for the Metropolis method chosen? Multiple modes are found using multiple Gauss-Newton minimizations, with strategic starting points. The sampling calibration uses the curvature scales and estimated topology to estimate good jumping choices. An alternative to the Metropolis approach is a “Bayesianized parametric bootstrap,” where, instead of trying to do a random walk around a very contorted and tight geometry,⁹ the data errors are “resampled,” and multiple MAP inversions computed for each resampled data set. Together, these methods dramatically increase the sampling acceptance rate, while ensuring that the sample space is adequately explored.

Another challenge is incorporating bound constraints on resistivity without slowing down the optimization scheme or biasing the posterior distribution sampling. The bounds are needed because unconstrained inversions often either stagnate in unwanted basins of attraction with unphysically low resistivities, or converge successfully, but with very low resistivity values that are nonphysical. From effective-media arguments using seawater and typical high-porosity sediments, a lower bound on resistivities of about 0.6 Ω-m is reasonable. The distribution is expected to be reasonably log-normal. In badly scaled problems, optimization methods that require hard bounds have to be coded carefully to prevent stagnation when the inversion trajectory “crashes” into a constraint.¹⁰ Similarly, sampling methods have to be implemented carefully to make sure they produce the right truncated distributions near hard boundaries.

TECHNICAL APPROACHES

Three different approaches were developed and applied to understand the resolution and uncertainty in the CSEM inversion. The first makes an assumption that the resistivity model is smooth in depth. The characteristic length of the smoothing is an additional parameter that is estimated in the inversion. This is done by parametrizing the smoothing using an exponential variogram structure with unknown smoothing/correlation length. The smoothing length is then determined by the noise in the CSEM data. The larger the noise, the greater the smoothing in the solution. We call this Bayesian smoothing.

The second method assumes the opposite: The Earth is made up of a set of discrete layers with constant properties discontinuous at the boundaries. We perform a hierarchical search through a “family tree” of models, constructed by recursive splitting or merging of layers in a starting “parent” model. Model parameter inversions are performed for each model geometry, and the final model geometry is selected using Bayesian model-selection techniques.

The third uses the physics of the problem to construct a fixed logarithmic depth grid consistent with the expected depth resolution of the CSEM method. That geometry is kept fixed, no smoothing is imposed, and constant properties are assumed within a layer with discontinuity at the layer boundaries.³

RESOLUTION AND UNCERTAINTY

The first model is shown in Fig. 3: a thin resistor of 10 Ω-m, 640 m below a sea layer 780 m in depth. The background resistivity increases from 1 Ω-m to 1.6 Ω-m in three steps. Figure 3a shows the result using Bayesian smoothing. Various noise levels were assumed increasing from 2% to 10%. Note that the greater the noise level, the greater the smoothing length, hence the lower the resolution—our first lesson. Figure 3b shows the result of the second type of inversion: Bayesian model selection from splitting/merging. The model in this case is a 100-m-thick, 100-Ω-m resistor buried 1 km under deep sea. Note that the optimal model has four layers with an obvious identification of the resistive anomaly. This behavior is generic, and constitutes our second lesson: CSEM data can be fit with quite parsimonious models, typically of order 10 parameters per common midpoint.

Fig. 3. Two types of inversion: a) Bayesian smoothing where the “truth case” model is shown as well as a series of inversion where the noise level was varied between 2% and 10%; b) Bayesian model selection combined with splitting and merging.

To understand the resolution in more detail, a wedge model was constructed. Two types of inversion were done: Bayesian smoothing (Fig. 4b) and fixed logarithmic gridding with no smoothing, Fig. 4c. Figure 4d shows the results of the estimate of the resistivity-thickness product (RTP) and the thickness with uncertainty. Note that the RTP is determined much better than the thickness (scaled uncertainty in RTP is more than an order of magnitude smaller). This is an explicit manifestation of poorly scaled inversion. There is no apparent bias in the estimate of either quantity. Even though the thickness is not as well determined as the RTP, there is distinct correlation between the estimate and the true value with no evident bias in the mean and a standard deviation from the true value consistent with the estimated standard deviation.

Fig. 4. a) Truth-case wedge model. b) Most likely inversion image using Bayesian smoothing on a regular 50-m grid. c) Most likely inversion image using a logarithmic grid without smoothing. d) Three-layer inversions for depth, thickness and resistivity, showing marginal-distribution 95% error bars for thickness and RTP.

Not only does the inversion estimate the uncertainty in quantities, but it can also can be used to investigate the correlation in the estimates of different quantities. Figure 5 shows the correlation of two such properties: the thickness and the depth of the resistor for a specific point along the wedge model.

Fig. 5. a) Contour plot of model probability, or marginal model likelihood (MML), of the three-layer model fit to common-midpoint-5 data in the wedge model, as a function of reservoir depth and thickness. Contours are at unit spacings of log10(MML), so three contours is approximately the conceivable span of model support in the data. Marginal distributions of b) thickness and c) depth of the reservoir layer associated with the MML measure are shown with the original truth-case model values from which the data were generated and the most likely estimates of parameters.

Finally, the reliability of the uncertainty estimates formed from Metropolis sampling is examined using a more complex “bird” model. Figure 6 compares the model to the results of log grid inversion with no smoothing. The model lies within one standard deviation of the truth case, as determined by examining the layer-by-layer statistics for the resistivity. What it does not capture is the correlation between the resistivities of different layers.  In order to examine this correlation, a multivariate cluster analysis was performed on the ensemble of realizations produced by the inversion.

Fig. 6. Statistical resistivity quantiles for the “bird” model.

IDENTIFYING ANOMALIES

The ensemble of realizations can be visualized by displaying the resistivity profile for each realization as a stacked grayscale image sequence, Fig. 7. All of the models fit the data to within the error. Note the uncertainty in the depth of the two resistive anomalies. There is obviously a correlation in the resistivities near these anomalies. Only one layer will have a larger resistivity. This image also explains the “jumpiness” seen in the resistor depths of weakly smoothed deterministic inversions. If the resistivity profile of the Earth were flat, this image would resemble a deterministic inversion’s cross-section through the Earth, with resistor location jumping around as the noise is resampled. This characteristic reflects the uncertainty in depth. To render a more stable image, the smoothing length, or regularization, of the deterministic inversion could be increased, but the continuity of the resistor would be increased at the cost of decreasing the peak resistivity. Without careful interpretation, this could increase the risk of predicting a lower-resistivity objective that is not a hydrocarbon feature.

Fig. 7. Multiple realizations of a resistivity profile that is all consistent with the CSEM data within a reasonable noise level. The horizontal axis is the realization number.

To better understand the correlation and to identify the resistive anomalies, a cluster analysis was applied to the set of realizations. A separate data point was created for each layer of each realization. Each point was assigned three attributes: resistivity, depth (a random value between the layer’s top and bottom) and the second derivative of the resistivity in depth normalized by the resistivity divided by the square of width. A cluster analysis11 was run on these data points, and those clusters with an average normalized second derivative less than 0.5 were identified as anomalies, Fig. 8. For each anomaly, the distribution of resistivity and depth can be calculated. The model value falls well within these distributions.

Fig. 8. Scatter plot of all realizations showing the background trend (black) and two anomalies (blue and red). The distributions of the two anomalies in resistivity and depth are also shown, with black arrows indicating the true value in the “bird” model.

SYSTEMATIC NOISE EFFECTS

As discussed in the introduction, modeling assumptions can introduce strongly correlated modeling error. This violates the assumption of Gaussian noise and can lead to bias in the result of the inversion. Correlated measurements are worth less than independent ones, so we expect uncertainties to increase if systematic effects apply. We have two approaches to compensating for this systematic noise, both assuming that the ratio of systematic to random noise is estimable.12 In both cases, the uncertainty in the result is increased as a result of the systematic noise.

The first approach regresses systematic noise trends directly as a smooth function of the offset, a common modeling artifact that can easily be confused with a real model response. The size of the additional term shown in Fig. 9a is related to the ratio of the systematic to random noise power. The greater this ratio, the larger the uncertainty in the result. A second approach is reducing the data precision in proportion to the relative ratio of the systematic noise to total noise (Fig. 9b), leading again to an increase in the uncertainty of the result. To demonstrate this effect, data were generated from a 250-m resistor buried 850 m below the mudline in 1-km water. Random and systematic error of about 5% was added to the data. If no systematic error is assumed in the inversion, the thickness is underestimated by more than a standard deviation. If the correct amount of systematic error is assumed and corrected for by either approach, the standard deviation of the inversion result is increased sufficiently to cover the correct thickness of the resistor.

Fig. 9. a) Inferred thickness marginal distribution of the resistivity target as a function of relative tolerance of systematic errors in inline |E| data, using quadratic systematics with offset. Curves are colored according to the assumed ratio of systematic to random noise power. b) Very similar effects are obtained through data set reduction by a factor of 2.

BUSINESS VALUE

Oil and gas operators must make decisions in an uncertain world. Should a wildcat well be drilled? How much flexibility should be in the development plan, given that there is uncertainty in the reservoir? This decision making is complicated by the ambiguity in the data caused by uncertainty in the acquisition and natural ambiguity in the physical response of the Earth to the measurement. Only when the uncertainty in the inversion is known can these questions be answered.

Even more important is the ability to combine the knowledge of multiple measurements, such as CSEM, seismic and time lapse. The Bayesian approach and noise model determines how these measurements should be combined. This significantly reduces uncertainty, as each measurement provides complementary information. Future work is expected to consider joint and time-lapse inversion of seismic and CSEM data, and the Bayesian inversion of CSEM data in 2D and 3D.  

LITERATURE CITED

 ¹ “Special section: CSEM,” Geophysics, No. 72 (2), 2007.
 ² Constable, S. C., Parker, R. L. and C.G. Constable, “Occam’s inversion: A practical algorithm for generating smooth models from electromagnetic sounding data,” Geophysics, No. 52, 1987, pp. 289.
 ³ Gunning, J., Glinsky, M. and J. Hedditch, “Resolution and uncertainty in 1D CSEM inversion: A Bayesian approach and open-source implementation,” Geophysics, 75, No. 6, 2010.
 ⁴ Loseth, L., “Modeling of controlled source electromagnetic data, ”PhD thesis, NTNU, 2007.
 ⁵ Chen, J., Hoversten, G. M., Vasco, D., Rubin, Y. and Z. Hou, “A Bayesian model for gas saturation estimation using marine seismic AVA and CSEM data,” Geophysics, No. 72, 2007, WA85.
 ⁶ Sambridge, M., Gallagher, K., Jackson, A. and P. Rickwood, “Trans-dimensional inverse problems, model comparison and the evidence,” Geophysical Journal International, No. 167, 2006, pp. 528.
 ⁷ Evans, S. and P. Stark, “Inverse problems as statistics,” Inverse Problems 2002, No. 18, R55.
 ⁸ Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H. and E. Teller, “Equation of state calculations by fast computing machines,” Journal of Chemical Physics, No. 21, 1953, pp. 1087.
 ⁹ Oliver, D., He, N., and A Reynolds, “Conditioning permeability fields to pressure data,” Proceedings of the 5th European Conference on the Mathematics of Oil Recovery, 1996; P. K. Kitanidis, “Quasi-linear geostatistical theory for inversing,” Water Resources Research, No. 31, 1995, pp. 2411.
 ¹⁰ Bersekas, D. P., “Projected Newton methods for optimization problems with simple constraints,” MIT Technical Report, 1981.
 ¹¹ Fraley, C. and A. E. Raftery, “Enhanced model-based clustering, density estimation, and discriminant analysis software: MCLUST, ”Journal of Classification, No. 20, 2003, pp. 263.
 ¹² Gunning, J. and M.E. Glinsky, “Error modeling in Bayesian CSEM Inversion,” Extended Abstracts, 72nd EAGE Conference & Exhibition, Barcelona, 2010.