Inversion Uncertainty and Confidence

Summary

The uncertainty of an inversion has two parts: how log errors and endmember errors propagate into the volumes, and how well the model fits the data. The first is controlled by the Condition number of the inversion matrix, the second by the model error QC. Confidence combines both. The calculator gives the standard deviation of the three-component example from the log errors.

Inputs and outputs

Equations

Propagation of log errors. For the linear least-squares problem on the Mineral Inversion page, if the errors are independent with standard deviations \(\sigma_i\) and the model is right, the covariance of the volumes is:

The standard deviation of each volume is the square root of the diagonal. This holds when no non-negativity constraint is active. When one is, the volume at its bound has a smaller spread than this gives.

Conditioning. The Condition number of the weighted matrix is the ratio of its largest to its smallest singular value:

A large \(\kappa\) means that some combination of volumes changes the logs very little, so that a small log error produces a large change in that combination.

Propagation of endmember errors. To first order, an error \(\delta\mathbf A\) in the matrix changes the volumes by:

This error is systematic, not random, and it is not reduced by averaging over depth.

The calculator is the three-component case with errors in the bulk density and in \(U\) only (closure is exact). The coefficients of the volumes with respect to the data are the entries of the inverse matrix, which are the numerators of the formulas on the Mineral Inversion page divided by \(D\). The standard deviations are:

with \(D = \rhoQz(\uCc - \uFl) - \rhoCc(\uQz - \uFl) + \rhoFl(\uQz - \uCc)\) as before. The results do not depend on the measured values, only on the endmembers and the log uncertainties.

Single-value calculator

Behavior

With the default endmembers and uncertainties of 0.03 g/cm³ and 1 barn/cm³, the porosity is determined to 0.019 but quartz and calcite to 0.120 and 0.113. This is the practical meaning of a poorly conditioned system: the sum of the two minerals is known, the split is not. The plot sweeps the calcite \(U\). As it comes down toward the quartz value of 4.8 the minerals become harder to tell apart: the quartz standard deviation is 0.12 at 13.8, 0.21 at 10, 0.34 at 8 and 1.00 at 6 barns/cm³, while the porosity standard deviation hardly moves (0.019, 0.020, 0.023 and 0.041). The density difference between quartz and calcite is only 0.06 g/cm³, so the density log cannot help; the photoelectric log carries the separation, and its accuracy sets the answer.

Parameter guidance

Log uncertainties. Use the repeatability and accuracy of each curve after environmental corrections, not the bit resolution. Values used in practice are of the order 0.02 to 0.03 g/cm³ for density, 0.02 to 0.03 for neutron porosity, 3 to 5 µs/ft for sonic and about 1 barn/cm³ for \(U\). They also act as weights: a curve with a larger uncertainty has less influence on the answer, so raise the uncertainty of curves that are less reliable (for example the sonic in poor hole).

Constraint uncertainties. The uncertainty on the clay volume and kerogen volume rows should reflect the accuracy of the upstream result, about 0.03 and 0.01 here. They may be reduced if the upstream result is calibrated to core.

Confidence classes. A usable confidence flag needs both measures. As a starting point, for the model error \(E_{avg}\) and the propagated standard deviation of the result of interest:

Use the class of the worse of the two. The limits are starting points for review.

Worked example

First, the effect of the choice of logs on a three-component problem (quartz, calcite, water) and, second, a five-component problem solved with and without the clay and kerogen constraints, comparing a Monte Carlo result (300 trials with noise on the logs) with the linear propagation:

import itertools
import numpy as np

# Endmember responses. Illustrative mid-range values, not a recommendation.
#             density  neutron  sonic   U
#             g/cm3    v/v      us/ft   barns/cm3
END = {
    "quartz":  [2.65, -0.02,  55.5,  4.8],
    "calcite": [2.71,  0.00,  47.5, 13.8],
    "clay":    [2.55,  0.35, 110.0,  9.0],
    "kerogen": [1.26,  0.60, 160.0,  0.25],
    "water":   [1.00,  1.00, 189.0,  0.4],
}
LOGS = ["rhob", "nphi", "dt", "u"]
SIGMA = {"rhob": 0.03, "nphi": 0.03, "dt": 5.0, "u": 1.0, "vcl": 0.03, "vker": 0.01, "unity": 0.001}


def system(minerals, curves=LOGS, constraints=("vcl", "vker")):
    """Rows: the log curves, then the constraint rows, then the closure row."""
    rows, sig, labels = [], [], []
    for c in curves:
        rows.append([END[m][LOGS.index(c)] for m in minerals])
        sig.append(SIGMA[c]); labels.append(c)
    for c, name in (("vcl", "clay"), ("vker", "kerogen")):
        if c in constraints:
            rows.append([1.0 if m == name else 0.0 for m in minerals])
            sig.append(SIGMA[c]); labels.append(c)
    rows.append([1.0] * len(minerals)); sig.append(SIGMA["unity"]); labels.append("unity")
    return np.array(rows), np.array(sig), labels


def solve(A, sig, d):
    """Weighted least squares with V >= 0: try every active set, keep the best feasible one."""
    Aw, dw = A / sig[:, None], d / sig
    n = A.shape[1]
    best, best_cost = None, np.inf
    for k in range(1, n + 1):
        for subset in itertools.combinations(range(n), k):
            x = np.linalg.lstsq(Aw[:, subset], dw, rcond=None)[0]
            if (x < -1e-12).any():
                continue
            v = np.zeros(n)
            v[list(subset)] = x
            cost = np.sum((Aw @ v - dw) ** 2)
            if cost < best_cost:
                best, best_cost = v, cost
    return best, best_cost


truth = {"quartz": 0.45, "calcite": 0.20, "clay": 0.15, "kerogen": 0.05, "water": 0.15}

# 1. Three unknowns (quartz, calcite, water): how the choice of logs changes the conditioning.
print("quartz-calcite-water, no constraints, linear error propagation")
print(f"{'logs':14s} {'cond':>7s} {'sd quartz':>10s} {'sd calcite':>11s} {'sd water':>9s}")
for curves in (("rhob", "nphi"), ("rhob", "u"), ("rhob", "nphi", "u")):
    A, sig, _ = system(["quartz", "calcite", "water"], curves, constraints=())
    Aw = A / sig[:, None]
    cov = np.linalg.inv(Aw.T @ Aw)
    sd = np.sqrt(np.diag(cov))
    print(f"{'+'.join(curves):14s} {np.linalg.cond(Aw):7.0f} {sd[0]:10.3f} {sd[1]:11.3f} {sd[2]:9.3f}")

# 2. Five unknowns: Monte Carlo on the log errors against the linear estimate.
minerals = list(truth)
v_true = np.array([truth[m] for m in minerals])
rng = np.random.default_rng(11)
print()
print("five components, standard deviation of the result (volume fraction)")
print(f"{'case':26s}" + "".join(f"{m:>9s}" for m in minerals))
for label, cons in (("logs + closure only", ()), ("logs + VWCL + kerogen", ("vcl", "vker"))):
    A, sig, _ = system(minerals, constraints=cons)
    d0 = A @ v_true
    out = []
    for _ in range(300):
        out.append(solve(A, sig, d0 + rng.normal(0.0, sig))[0])
    out = np.array(out)
    Aw = A / sig[:, None]
    lin = np.sqrt(np.diag(np.linalg.inv(Aw.T @ Aw)))
    print(f"{label + ' (MC)':26s}" + "".join(f"{x:9.3f}" for x in out.std(axis=0)))
    print(f"{label + ' (linear)':26s}" + "".join(f"{x:9.3f}" for x in lin))
    print(f"{'  condition number':26s}{np.linalg.cond(Aw):9.0f}")

Output

quartz-calcite-water, no constraints, linear error propagation
logs              cond  sd quartz  sd calcite  sd water
rhob+nphi         1535      0.632       0.618     0.020
rhob+u             287      0.120       0.113     0.019
rhob+nphi+u        277      0.115       0.111     0.016

five components, standard deviation of the result (volume fraction)
case                         quartz  calcite     clay  kerogen    water
logs + closure only (MC)      0.118    0.116    0.099    0.096    0.084
logs + closure only (linear)    0.174    0.126    0.142    0.231    0.201
  condition number              794
logs + VWCL + kerogen (MC)    0.113    0.106    0.029    0.009    0.018
logs + VWCL + kerogen (linear)    0.115    0.112    0.029    0.010    0.018
  condition number              358

Assumptions and limitations

• Log errors are independent, random, and have the assigned standard deviations. Real errors are correlated across logs (a washout affects density, neutron and sonic together) and systematic.
• The model is right: the right components and the right endmembers. Model error is not propagated, and shows up as misfit.
• The linear propagation ignores the non-negativity bounds. Near a bound the real spread is smaller and skewed, as the Monte Carlo columns in the example show.
• Uncertainties of the constraint rows (clay and kerogen volume) are known and independent of the log errors, although they depend on the same logs in practice.
• The three-component calculator ignores error in the endmembers.

QC checks

• The propagated standard deviation of porosity is a few hundredths or less, and is smaller than the porosity itself in tight rock.
• Quartz and calcite (or any similar pair) standard deviations are reported with the volumes. A pair with deviations far above its volumes is not resolved and should be merged or reported as a sum.
• Adding a log, or a constraint, lowers the standard deviations. If it does not, it carries no new information about the volume of interest.
• The condition number does not change much from well to well unless the component list or logs change.
• Monte Carlo and linear results agree where no bound is active.

Going Deeper

Two measures matter and they are different. The condition number and the propagated standard deviation are a property of the model and the data errors, and can be computed before looking at any logs. The model error is a property of the fit and is computed from the logs afterwards. A system can be well conditioned and fit badly (a missing mineral), or fit perfectly and be badly conditioned (too many similar minerals), and in both cases the volumes are uncertain. In a full Bayesian treatment the same quantities appear as a posterior covariance, and priors (for example on clay volume and kerogen volume) work as extra rows, as they do in the inversion method.

References

1. Aster, R.C., Borchers, B. and Thurber, C.H., 2018. Parameter Estimation and Inverse Problems, 3rd edition. Elsevier.
2. Lawson, C.L. and Hanson, R.J., 1974. Solving Least Squares Problems. Prentice-Hall, Englewood Cliffs, NJ (reprinted by SIAM, 1995).

Python reference implementation

Related methods

• Mineral Inversion
• Inversion Endmember Parameters
• Inversion Model-Error QC