Inversion Model-Error QC

Summary

Model-error QC compares the logs predicted by the inversion with the measured logs and expresses the difference in units of the assigned uncertainty. Inversion average error near or below 1 means the model fits the data to within the stated errors; a high value flags intervals where the inversion fails. A good fit is necessary but not sufficient: a missing mineral can be absorbed with little misfit.

Inputs and outputs

Equations

For each active curve, predict the log from the solution, \(\hat d_i = \sum_j A_{ij}\hat V_j\), and form the residual \(r_i = d_i - \hat d_i\). Each residual is divided by the uncertainty of its curve. For the three logs of the calculator:

The average error divides by the number of active curves, here 3, so that it can be compared between intervals with different curves present:

Curves that are missing or switched off are left out of both the sum and the count. The RMS form weights large residuals more and is the square root of the usual reduced misfit.

Single-value calculator

Behavior

The average error is a V-shaped function of the density residual, with its minimum at zero, and is steeper when the density uncertainty is smaller. A density residual of 0.06 g/cm³ with the other curves perfect gives an average error of 1.00, 0.67 and 0.40 for density uncertainties of 0.02, 0.03 and 0.05. The same residual is therefore a warning or noise depending on the uncertainty assigned, so the uncertainties have to be set honestly before the error curve means anything.

Parameter guidance

Uncertainties. Use the same values as in the inversion. The error measure only has meaning relative to them.

Thresholds. As a starting point treat \(E_{avg}\) below 1 as consistent with the stated errors, 1 to 2 as a reason to look, and above 2 as a failure of the model in that interval. Values well below 0.5 over long intervals can mean the uncertainties are too large, or that there are more unknowns than the data can constrain.

Which residuals to map.

• Predicted against measured, one crossplot per curve, with the 1:1 line. Scatter about the line is the noise level; a curved or offset cloud is a systematic problem.
• Residual against depth, in the same track as the logs, with the uncertainty band shaded. Look for intervals where the measured curve leaves the band.
• Residual against caliper and density correction. A residual that follows the caliper is borehole damage, not mineralogy.
• Residual against each component volume. A residual that correlates with one component points to an error in that component's endmember.
• Pattern of signs across curves. A mineral that is missing from the model leaves a characteristic set of residual signs: measured above predicted on U, and on density, for pyrite, for example.

Worked example

Five rocks of known composition, inverted with a model that has only quartz, calcite, clay, kerogen and water. A is a rock the model describes; B has borehole damage added to the logs (density 0.25 g/cm³ low, neutron 0.06 high, sonic 15 µs/ft high); C has gas in the pores; D and E have 5% and 15% pyrite that the model does not have. The columns show each log residual in units of its uncertainty, the average, and the volumes the inversion reports:

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


END["gas"] = [0.25, 0.30, 220.0, 0.0]
END["pyrite"] = [4.99, -0.02, 39.2, 84.0]
base = ["quartz", "calcite", "clay", "kerogen", "water"]
A5, sig5, labels = system(base)           # the model the analyst runs: five components

cases = {
    "A consistent":        ({"quartz": .45, "calcite": .20, "clay": .15, "kerogen": .05, "water": .15}, [0, 0, 0, 0]),
    "B washout":           ({"quartz": .45, "calcite": .20, "clay": .15, "kerogen": .05, "water": .15}, [-0.25, 0.06, 15.0, 0]),
    "C gas in pores":      ({"quartz": .45, "calcite": .20, "clay": .15, "kerogen": .05, "water": .08, "gas": .07}, [0, 0, 0, 0]),
    "D pyrite (5%)":       ({"quartz": .40, "calcite": .20, "clay": .15, "kerogen": .05, "water": .15, "pyrite": .05}, [0, 0, 0, 0]),
    "E pyrite (15%)":      ({"quartz": .30, "calcite": .20, "clay": .15, "kerogen": .05, "water": .15, "pyrite": .15}, [0, 0, 0, 0]),
}
print(f"{'case':16s} {'r_rhob':>7s} {'r_nphi':>7s} {'r_dt':>6s} {'r_u':>6s} {'avg|r|':>7s} {'quartz':>7s} {'calcite':>8s} {'water':>6s}")
for name, (comp, bad) in cases.items():
    full = list(comp)
    A, s, _ = system(full)
    d = A @ np.array([comp[m] for m in full])      # data generated with the real composition
    d[:4] += np.array(bad)                          # borehole damage added to the logs
    d[-1] = 1.0
    v, _ = solve(A5, sig5, d)                       # same seven rows, five unknowns
    r = (d - A5 @ v) / sig5
    print(f"{name:16s} {r[0]:7.2f} {r[1]:7.2f} {r[2]:6.2f} {r[3]:6.2f} {np.mean(np.abs(r[:4])):7.2f} "
          f"{v[0]:7.3f} {v[1]:8.3f} {v[4]:6.3f}")

Output

case              r_rhob  r_nphi   r_dt    r_u  avg|r|  quartz  calcite  water
A consistent       -0.00    0.00   0.00  -0.00    0.00   0.450    0.200  0.150
B washout          -1.41   -2.10  -0.16   0.44    1.03   0.319    0.222  0.276
C gas in pores     -1.06   -1.83   0.21   0.41    0.88   0.488    0.164  0.161
D pyrite (5%)       1.09    0.74   1.30  -0.07    0.80   0.048    0.621  0.114
E pyrite (15%)      2.75    3.96   3.76   6.09    4.14   0.000    0.751  0.010

Assumptions and limitations

• The assigned uncertainties are realistic. The error measure is only as meaningful as they are.
• The residuals are independent between curves. In practice a washout moves density, neutron and sonic together, so the average understates the problem.
• Constraint rows (clay volume, kerogen volume, closure) are left out of the error measure here, as they are in the log-only form. Including them gives a different but useful measure of how far the logs pull the solution from the upstream result.
• A small misfit means the model describes the data, not that it describes the rock. The model can absorb a missing component by changing others.

QC checks

• E_avg is mostly below 1 in good-quality hole and rises in known problem intervals, such as washouts, thin beds and tool-calibration problems.
• Flag intervals where E_avg exceeds the chosen limit and carry the flag to the final results, so that volumes there are not used for cutoffs without review.
• Compare the flags with the caliper, the density correction and the cycle-skip indicator before reading a high misfit as geology.
• A mineral that is geologically expected but not in the model (pyrite, anhydrite, siderite, salt, coal) is a candidate whenever a high misfit extends over a thick interval with a consistent residual sign pattern.
• Test the model by adding the suspected component and checking that the misfit falls and the other volumes move by an amount that makes geological sense.

Going Deeper

The example shows the central limitation. With seven equations for five unknowns, the system has little redundancy, and a damaged or incomplete data set is absorbed by shifting volume between components. Borehole damage of 0.25 g/cm³ in density (more than eight times its uncertainty), 0.06 in neutron and 15 µs/ft in sonic produces an average error of only about 1, and 5% pyrite produces less than 1 while moving quartz from 0.40 to 0.05. For that reason model error is used together with the borehole flags and the conditioning, and not alone. Adding logs raises the redundancy and makes the error measure more sensitive. In practice the most reliable sign of a missing mineral is a persistent residual pattern over a thick interval, not a single high value.

References

1. Quirein, J., Kimminau, S., LaVigne, J., Singer, J. and Wendel, F., 1986. A coherent framework for developing and applying multiple formation evaluation models. Transactions of the SPWLA 27th Annual Logging Symposium, Paper DD.
2. Doveton, J.H., 1994. Geologic Log Analysis Using Computer Methods. AAPG Computer Applications in Geology No. 2, American Association of Petroleum Geologists, Tulsa, OK.

Python reference implementation

Related methods

• Mineral Inversion
• Inversion Uncertainty and Confidence
• Inversion Endmember Parameters
• Gas and Organic-Matter Corrections