Unsupervised Clustering

Summary

Unsupervised clustering groups depth samples by how close they are on the standardised logs, with no core description. k-means clustering is the standard choice: it splits the samples into k groups, each described by a centre. The number of clusters is chosen with the Silhouette width and the geology, and the clusters are given facies names only afterwards, by reading their centres. Use it to explore a well, to find groups a rule table missed, or when no labelled data exist.

Inputs and outputs

Equations

Each sample is a vector \(x_i\) of standardised logs: each curve has its mean subtracted and is divided by its standard deviation, and resistivity enters as its logarithm. k-means finds \(k\) centres \(\mu_1, \ldots, \mu_k\) and assignments \(c_i\) that minimise the within-cluster sum of squares (the inertia):

It alternates two steps until the assignments stop changing: assign each sample to the nearest centre, then move each centre to the mean of its samples. The result depends on the starting centres, so the run is repeated from several seeds (k-means++ seeding picks well-spread starting centres) and the lowest inertia is kept.

The inertia always falls as \(k\) grows, so it cannot choose \(k\) by itself. The silhouette width of sample \(i\) compares the mean distance \(a_i\) to the other samples of its own cluster with the mean distance \(b_i\) to the samples of the nearest other cluster:

and the mean of \(s_i\) over all samples is the score of a clustering, between -1 and 1. A Gaussian mixture is the soft version: it fits \(k\) Gaussian components with their own covariance and gives each sample a probability for each component, and the number of components is chosen with an information criterion in place of the silhouette. Neither method uses depth, so the result is a label per sample, which is smoothed afterwards. A majority filter of width \(w\) replaces each label by the commonest label in a window of \(w\) samples.

Single-value calculator

No calculator: this method is a procedure, not a single equation. The worked example below runs the whole procedure on synthetic logs.

Behavior

On the synthetic well of 687 samples and four true facies, the inertia falls steadily from 1282 at k = 2 to 393 at k = 7, with the largest drop to k = 3 and a smaller one after k = 4, but the silhouette peaks at 0.561 for k = 3 and is 0.496 at k = 4. The silhouette therefore prefers three groups: the limestone and the dolomite overlap on the density and are merged at k = 3. The geologist's four facies are recovered at k = 4 anyway, with a purity of 0.988 against the known facies, and the five-cluster silhouette drops to 0.346 because a real group is split. Choosing k by one score alone would have merged a facies, so the score is a guide and not the decision. The 5-sample majority filter lowers the number of label changes down the well from 38 to 27, which equals the 27 changes of the true beds, and raises the agreement from 0.988 to 0.997.

Parameter guidance

Curves. Use curves that respond to different rock properties, such as gamma ray, bulk density, neutron porosity and the logarithm of resistivity; see Curve Selection and Missing Data. Standardisation. Without it, the curve with the largest numbers (gamma ray in API units, resistivity in ohm·m) dominates the distance. Standardise with statistics of the key well, and apply the same values to the other wells after normalization. Number of clusters. Compute the inertia and the silhouette for k from 2 to about 8, look for the elbow and the silhouette peak, and then check each choice against the geology: a cluster that cannot be named is a candidate to merge, and a cluster that mixes two known facies is a candidate to split. Restarts. Use at least five seeds and keep the lowest inertia. Labelling. After clustering, name each cluster from its centre in physical units (gamma ray, density, neutron, resistivity), a crossplot, or core. Do the naming once for the key well and keep the names with the saved model, because a cluster number has no meaning from one run to the next. Depth continuity. Apply a majority filter with a width near the thinnest bed that matters.

Worked example

A k-means written in numpy, with k-means++ seeding, five restarts, a silhouette function written out in full and a majority filter. It is run on the synthetic logs for k from 2 to 7. The four-cluster result is labelled from its centres, compared with the known facies in a contingency table, and filtered for depth continuity:

import numpy as np

NAMES = ["Shale", "Sandstone", "Limestone", "Dolomite"]
# facies means: GR (API), RHOB (g/cm3), NPHI (v/v, limestone scale), Rt (ohm.m). Illustrative values.
MEAN = np.array([[105.0, 2.55, 0.30, 3.0],
                 [45.0, 2.34, 0.17, 12.0],
                 [28.0, 2.60, 0.09, 35.0],
                 [28.0, 2.72, 0.07, 45.0]])


def synthetic_well(n_beds=40, seed=7, probs=(0.35, 0.30, 0.25, 0.10), bed_scale=1.0):
    """Beds of 8 to 25 samples. Each bed has its own offset, shared by all its samples, and each
    sample adds independent noise. bed_scale multiplies the bed offsets. Resistivity is simulated in log10."""
    rng = np.random.default_rng(seed)
    truth, logs = [], []
    for _ in range(n_beds):
        f = rng.choice(4, p=probs)
        n = rng.integers(8, 26)
        bed = rng.normal(0, np.array([8.0, 0.04, 0.025, 0.20]) * bed_scale)
        gr = MEAN[f, 0] + bed[0] + rng.normal(0, 12.0, n)
        rhob = MEAN[f, 1] + bed[1] + rng.normal(0, 0.04, n)
        nphi = MEAN[f, 2] + bed[2] + rng.normal(0, 0.025, n)
        rt = 10 ** (np.log10(MEAN[f, 3]) + bed[3] + rng.normal(0, 0.15, n))
        truth += [f] * n
        logs.append(np.column_stack([gr, rhob, nphi, rt]))
    return np.vstack(logs), np.array(truth)


def standardise(x):
    """Resistivity enters as log10; every column then gets zero mean and unit standard deviation."""
    x = x.copy()
    x[:, 3] = np.log10(x[:, 3])
    return (x - x.mean(axis=0)) / x.std(axis=0), x


def kmeans(z, k, seed=0, n_iter=100):
    """Plain k-means with k-means++ seeding. Returns labels, centroids and the inertia."""
    rng = np.random.default_rng(seed)
    centres = [z[rng.integers(len(z))]]
    for _ in range(k - 1):
        d2 = np.min([((z - c) ** 2).sum(axis=1) for c in centres], axis=0)
        centres.append(z[rng.choice(len(z), p=d2 / d2.sum())])
    c = np.array(centres)
    for _ in range(n_iter):
        d = ((z[:, None, :] - c[None, :, :]) ** 2).sum(axis=2)
        lab = d.argmin(axis=1)
        new = np.array([z[lab == j].mean(axis=0) if np.any(lab == j) else c[j] for j in range(k)])
        if np.allclose(new, c):
            break
        c = new
    d = ((z[:, None, :] - c[None, :, :]) ** 2).sum(axis=2)
    lab = d.argmin(axis=1)
    return lab, c, d[np.arange(len(z)), lab].sum()


def silhouette(z, lab):
    """Mean silhouette width: (b - a) / max(a, b), with a the mean distance inside the cluster
    and b the mean distance to the nearest other cluster."""
    dist = np.sqrt(((z[:, None, :] - z[None, :, :]) ** 2).sum(axis=2))
    s = np.zeros(len(z))
    for i in range(len(z)):
        same = lab == lab[i]
        if same.sum() < 2:
            continue
        a = dist[i, same].sum() / (same.sum() - 1)
        b = min(dist[i, lab == j].mean() for j in np.unique(lab) if j != lab[i])
        s[i] = (b - a) / max(a, b)
    return s.mean()


def majority_filter(lab, width=5):
    """Replace each label by the most common label in a centred window (depth continuity)."""
    half = width // 2
    out = lab.copy()
    for i in range(len(lab)):
        w = lab[max(0, i - half): i + half + 1]
        out[i] = np.bincount(w).argmax()
    return out


logs, truth = synthetic_well()
z, phys = standardise(logs)
print("choosing k (best of 5 starts for each k):")
print("   k   inertia  silhouette")
runs = {}
for k in range(2, 8):
    best = min((kmeans(z, k, seed=s) for s in range(5)), key=lambda r: r[2])
    runs[k] = best
    print(f"  {k:2d}  {best[2]:8.1f}  {silhouette(z, best[0]):9.3f}")

k = 4
lab, cent, _ = runs[k]
# physical-unit centroids, then label each cluster after the fact by its position on the logs
print("\ncluster centroids in log units and the label given afterwards:")
names = {}
for j in range(k):
    m = phys[lab == j].mean(axis=0)
    ndsep = m[2] - (2.71 - m[1]) / 1.71
    if m[0] >= 75:
        names[j] = "Shale"
    elif ndsep < -0.01:
        names[j] = "Sandstone"
    elif m[1] >= 2.66:
        names[j] = "Dolomite"
    else:
        names[j] = "Limestone"
    print(f"  cluster {j}: GR {m[0]:6.1f}  RHOB {m[1]:.3f}  NPHI {m[2]:.3f}  Rt {10 ** m[3]:5.1f}  n={np.sum(lab == j):4d}  -> {names[j]}")

# agreement with the (known, synthetic) facies: contingency table, clusters in rows
tab = np.zeros((k, 4), dtype=int)
np.add.at(tab, (lab, truth), 1)
print("\ncontingency table (rows clusters, columns true facies):")
print("               " + " ".join(f"{n[:5]:>6s}" for n in NAMES))
for j in range(k):
    print(f"  cl {j} {names[j][:5]:5s}  " + " ".join(f"{v:6d}" for v in tab[j]))
print(f"purity {tab.max(axis=1).sum() / tab.sum():.3f}")

# depth continuity: count label changes down the well before and after a 5-sample majority filter
smooth = majority_filter(lab, 5)
to_code = np.array([NAMES.index(names[j]) for j in range(k)])
print(f"\nlabel changes down the well: raw {np.sum(lab[1:] != lab[:-1])}, filtered {np.sum(smooth[1:] != smooth[:-1])}, true beds {np.sum(truth[1:] != truth[:-1])}")
print(f"agreement with true facies: raw {np.mean(to_code[lab] == truth):.3f}, filtered {np.mean(to_code[smooth] == truth):.3f}")

Output

choosing k (best of 5 starts for each k):
   k   inertia  silhouette
   2    1282.5      0.521
   3     580.1      0.561
   4     515.3      0.496
   5     462.6      0.346
   6     425.2      0.239
   7     393.3      0.234

cluster centroids in log units and the label given afterwards:
  cluster 0: GR   44.4  RHOB 2.327  NPHI 0.168  Rt  11.7  n= 318  -> Sandstone
  cluster 1: GR   34.6  RHOB 2.564  NPHI 0.098  Rt  36.0  n= 124  -> Limestone
  cluster 2: GR  100.2  RHOB 2.538  NPHI 0.294  Rt   2.8  n= 191  -> Shale
  cluster 3: GR   32.9  RHOB 2.749  NPHI 0.075  Rt  37.6  n=  54  -> Dolomite

contingency table (rows clusters, columns true facies):
                Shale  Sands  Limes  Dolom
  cl 0 Sands       0    318      0      0
  cl 1 Limes       0      0    119      5
  cl 2 Shale     191      0      0      0
  cl 3 Dolom       0      0      3     51
purity 0.988

label changes down the well: raw 38, filtered 27, true beds 27
agreement with true facies: raw 0.988, filtered 0.997

Assumptions and limitations

• Clusters in log space correspond to facies. Rock with a gradual change of properties, such as a sand with changing porosity, forms a continuum that k-means cuts into arbitrary slices.
• k-means assumes compact, similarly sized, roughly spherical clusters on the standardised axes. Elongated or very unequal clusters are split or merged badly; a Gaussian mixture with full covariance handles elongated clusters better.
• The curves are normalized, repaired and complete. A shifted or missing curve changes the distance between samples and therefore the groups.
• The number of clusters is known or can be chosen. The silhouette finds the best separation, which is not always the geologically useful split.
• Samples are independent. Depth is not used, so a thin bed of a rare facies can be absorbed into a larger cluster, and a noisy log gives a flickering track.

QC checks

• Every cluster can be named, and the names are consistent with the geology and the core. A cluster with no name is a reason to try a smaller k.
• The result is stable: a different seed or a resample of the well gives the same clusters, up to the order of the numbers.
• The silhouette of each cluster is positive and none is a tiny cluster of a few samples that could be a log artefact, such as a washout or a casing effect.
• Label changes down the well are close to the number of beds that the geologist would pick. A track that changes every few samples needs the depth filter.
• The same clustering model, and the same cluster names, are applied to the other wells and give a facies proportion that the geology supports.

Going Deeper

The k-means algorithm is attributed to Lloyd, whose least-squares quantization was written in 1957 and published in 1982, and to MacQueen, who named it in 1967. The k-means++ seeding of Arthur and Vassilvitskii, in 2007, made the result less dependent on the start. In petrophysics the method is used for electrofacies, a term for groups defined by their log response and not by core description. Self-organising maps and hierarchical clustering are older alternatives. A fixed model fitted on a key well and applied to other wells, which is how clustering is used in a field study, is only meaningful after normalization, since a cluster centre fixed in absolute log values is moved by a tool shift. A Gaussian mixture gives a probability for each facies, which allows a result to be flagged as uncertain where two probabilities are close, and is preferred when a probability, not just a label, is needed downstream.

References

1. Lloyd, S.P., 1982. Least squares quantization in PCM. IEEE Transactions on Information Theory, 28(2), 129–137.
2. Arthur, D. and Vassilvitskii, S., 2007. k-means++: the advantages of careful seeding. Proceedings of the 18th ACM-SIAM Symposium on Discrete Algorithms (SODA), 1027–1035.
3. Rousseeuw, P.J., 1987. Silhouettes: a graphical aid to the interpretation and validation of cluster analysis. Journal of Computational and Applied Mathematics, 20, 53–65.

Python reference implementation

Related methods

• Deterministic Rules from Crossplots
• Supervised Classification
• Curve Selection and Missing Data