Spike Removal and Filtering

Summary

Spikes are single-sample or very narrow excursions that are not formation signal: tool glitches, telemetry errors, cycle skips and noise. A Median filter or a Hampel filter removes them, but any filter also damages real thin beds. This page describes how the two work, how to choose the window and threshold, and how to tell when too much has been removed.

Inputs and outputs

Equations

A median filter of odd window \(W = 2h + 1\) replaces each sample with the median of the window around it:

It changes nearly every sample, however small the change. The Hampel filter replaces a sample only when it is an outlier relative to the local median. With the local median absolute deviation

the sample is flagged when

The factor 1.4826 scales the MAD to a standard deviation for Gaussian noise. The floor \(\sigma_{\min}\), in curve units, is needed because the MAD tends to zero in a flat zone, where every small variation would otherwise be flagged. A flagged sample is replaced by \(\tilde{x}_i\) or set to missing. The usual threshold is \(k = 3\).

A feature of \(n\) samples survives a median filter of window \(W\) only if \(n > W/2\). A narrower feature is removed, whether it is noise or a thin bed.

Single-value calculator

No calculator: this method is a procedure, not a single equation.

Behavior

The synthetic curve has 200 samples at 0.5 ft, five planted single-sample spikes and a real 3-sample (1.5 ft) thin bed. Every setting of the Hampel filter in the table catches all five spikes. A window of 5 with k = 3 flags exactly 5 samples (2.5%) and keeps the thin bed. Windows of 7, 11 and 21 each flag 8 samples (4.0%): the 5 spikes plus the 3 samples of the thin bed, which are removed. A larger k does not protect the thin bed here, because it is the window and not the threshold that defines which features survive. The plain 5-point median alters 155 of the 200 samples, which is why it should not be used when the aim is to remove only spikes.

Parameter guidance

Window. Choose it in feet, not samples: the window length should be less than twice the thickness of the thinnest bed you want to keep, since a feature must occupy more than half the window to survive. At 0.5 ft, a window of 5 keeps features of 1.5 ft and a window of 7 needs 2 ft. Spikes are one or two samples wide, so a window of 5 is the smallest that works.

Threshold k and noise floor. k = 3 is a common starting value. Raise it to flag fewer samples. Set \(\sigma_{\min}\) near the noise level of the curve in a flat interval; it is the curve's own noise and not a universal number. In the example it is 2.0 gAPI for a noise of 1.5.

Curve type. Despike resistivity in logarithms, because its noise is multiplicative. Sonic cycle skips are positive spikes that can be several samples wide; they are better flagged with a rule on slowness than with a narrow filter. Do not despike the caliper before using it as a washout flag, since the filter would hide narrow washouts.

How much is too much. Count the flagged samples. A few percent of a good curve is normal. If more than about 5 to 10% are flagged, either the threshold or window is wrong or the curve is bad and needs repair and not filtering. Also look at the flags: they should be isolated samples and not clustered at bed boundaries. As a further check, the removed part (original minus cleaned) should look like noise, with no formation signal.

Keep the original. Store the flags and the cleaned curve separately and keep the raw curve. Filtering before the thin-bed analysis, or the cutoff count, changes the answer. Where a curve is smoothed for display only, say so.

Worked example

A blocky gamma-ray curve with noise, five planted spikes and one real thin bed, filtered with different windows and thresholds. The table shows spikes caught and whether the thin bed survives:

import numpy as np


def running_median(x, window):
    """Centred running median; the window shrinks at the ends. window must be odd."""
    h = window // 2
    return np.array([np.nanmedian(x[max(0, i - h): i + h + 1]) for i in range(len(x))])


def hampel(x, window=7, k=3.0, min_sigma=2.0):
    """Flag samples further than k scaled MADs from the local median. Returns (cleaned, flags)."""
    h = window // 2
    med = running_median(x, window)
    mad = np.array([np.nanmedian(np.abs(x[max(0, i - h): i + h + 1] - med[i])) for i in range(len(x))])
    sigma = 1.4826 * mad                     # MAD scaled to a standard deviation for Gaussian noise
    # min_sigma is a noise floor in curve units: in a flat zone the MAD tends to zero and everything would flag
    flags = np.abs(x - med) > k * np.maximum(sigma, min_sigma)
    y = x.copy()
    y[flags] = med[flags]
    return y, flags


if __name__ in ("__main__", "worked_example"):
    rng = np.random.default_rng(8)
    step = 0.5                                              # ft
    truth = np.repeat([60.0, 110.0, 45.0, 95.0, 70.0], 40)      # five beds, 20 ft each
    truth[100:103] = 20.0                                   # a real 1.5 ft thin bed (3 samples)
    x = truth + rng.normal(0, 1.5, len(truth))
    spikes = [15, 52, 77, 133, 170]
    x[spikes] += np.array([60, -50, 70, 55, -45])           # five single-sample spikes
    print(f"{len(x)} samples, step {step} ft, {len(spikes)} planted spikes, a 3-sample thin bed at 100-102")
    print(f"{'window':>8} {'k':>4} {'flagged':>8} {'% of curve':>11} {'spikes caught':>14} {'thin bed kept':>14}")
    for window, k in ((5, 3.0), (7, 3.0), (7, 6.0), (11, 3.0), (21, 3.0)):
        y, flags = hampel(x, window, k)
        caught = int(np.sum(flags[spikes]))
        kept = np.mean(y[100:103]) < 40
        print(f"{window:8d} {k:4.1f} {int(flags.sum()):8d} {100 * flags.mean():10.1f}% {caught:11d}/5 {str(kept):>14}")
    print()
    m = running_median(x, 5)
    print("plain 5-point running median: spikes at their samples ->",
          np.round(m[spikes] - truth[spikes], 1), "(error vs truth)")
    print(f"plain median changes {np.sum(np.abs(m - x) > 0)} of {len(x)} samples (it alters almost every sample)")

Output

200 samples, step 0.5 ft, 5 planted spikes, a 3-sample thin bed at 100-102
  window    k  flagged  % of curve  spikes caught  thin bed kept
       5  3.0        5        2.5%           5/5           True
       7  3.0        8        4.0%           5/5          False
       7  6.0        8        4.0%           5/5          False
      11  3.0        8        4.0%           5/5          False
      21  3.0        8        4.0%           5/5          False

plain 5-point running median: spikes at their samples -> [ 2.2 -1.7  0.8  0.1  0.5] (error vs truth)
plain median changes 155 of 200 samples (it alters almost every sample)

Assumptions and limitations

• Spikes are narrower than the window and than any real feature that matters. A wide artefact is not removed.
• The noise is similar over the window, so one scale estimate (MAD) is meaningful. Strong changes of noise level inside a window affect the flags.
• The curve has no nulls. Nulls must be converted first, and the running median ignores them.
• The curve is sampled uniformly, because the window is counted in samples.
• The vertical resolution of the tool is not finer than the window. Otherwise real features are treated as spikes.

QC checks

• The flagged samples are a small fraction of the curve and are mostly isolated.
• The difference between the raw and cleaned curve looks like noise and has no bed-boundary structure.
• Known thin beds (a marker, a coal, a tight streak) survive the filter.
• The cleaned curve compares with other curves that measure related quantities (for example the density with the neutron) with no new disagreement.
• Cross-plots of the cleaned curve have fewer outliers, and the bulk of the data is unchanged.

Going Deeper

The median filter dates from the statistics of robust smoothing and was widely adopted in signal and image processing because it removes impulse noise while keeping edges, which a moving average blurs. The Hampel identifier, usually credited to Hampel (1974), is a robust outlier rule based on the median and MAD. Both fail when the 'spike' is the signal, as in a thin, high-contrast bed, and the distinction between noise and signal is not made by the filter but by the analyst. Where spikes come from a known cause, such as cycle skipping in the sonic or telemetry dropouts, a rule tied to the cause (a slowness limit, a flag from the tool) is better than a generic filter. Smoothing to match the vertical resolution of one tool to another, for example before comparing density with a lower-resolution neutron, is a separate operation from despiking and should be done with a resolution matching filter.

References

1. Hampel, F.R., 1974. The influence curve and its role in robust estimation. Journal of the American Statistical Association, 69(346), 383–393.

Python reference implementation

Related methods

• Null and Invalid Value Handling
• Depth Reference and Sampling Checks
• Log-Quality Flags
• Caliper-Based Flags