Signal Smoothing in Seismology

Smoothing is a fundamental preprocessing step in seismic signal processing. It reduces high-frequency noise while preserving the underlying trends or structures in the signal. Smoothing is often used before visualization, detection, or feature extraction.

In ISP, the smoothing function supports three types of smoothing techniques:

• Mean Smoothing
• Gaussian Smoothing
• Teager-Kaiser Energy Operator (TKEO)

1. Mean Smoothing

This technique replaces each sample with the average of its neighbors within a fixed time window. It’s useful for simple noise reduction.

Equation:

Python Snippet:

for i in range(k, n - k - 1):
    filtsig[i] = np.mean(tr.data[i - k:i + k])

Pros:
- Easy to implement
- Effective for simple noise

Cons:
- Can blur signal edges or transients
- Uniform weighting may suppress meaningful features

2. Gaussian Smoothing

Gaussian smoothing weights neighbors with a bell-shaped kernel, giving more importance to samples near the center.

Equation (kernel):

The kernel is normalized to preserve signal energy:

Smoothing Operation:

Python Snippet:

gauswin = np.exp(-(4 * np.log(2) * gtime ** 2) / fwhm ** 2)
gauswin /= np.sum(gauswin)
for i in range(k + 1, n - k - 1):
    filtsigG[i] = np.sum(tr.data[i - k:i + k] * gauswin)

Pros:
- Smooths without sharp cutoffs
- Preserves shape better than mean smoothing

Tips:
- Use FWHM in seconds to control the spread
- Always use tapering/detrending before smoothing

3. Teager-Kaiser Energy Operator (TKEO)

TKEO is a nonlinear operator used to enhance the energy content of rapid oscillations. It's especially useful in event detection and muscle/EMG analysis.

Equation:

Python Snippet:

emgf[1:-1] = emg[1:-1] ** 2 - emg[0:-2] * emg[2:]

Pros:
- Highlights instantaneous energy
- Useful for impulsive signal detection

Cons:
- Sensitive to noise
- Should be followed by smoothing or thresholding

Summary Table

Signal smoothing is a simple yet powerful tool to reduce noise and extract structure. Use with care to preserve signal characteristics needed for further seismic analysis.

4. Spike Removal with the Hampel Filter

The Hampel filter is a robust method for detecting and removing outliers (spikes) in time series. It replaces points that deviate significantly from the local median using a threshold based on standard deviation.

Concept:

For each point \( x[n] \), define a window \( W \) around it and compute: - Median: \( m \) - Median Absolute Deviation (MAD): \( ext{MAD} = ext{median}(|x[i] - m|) \)

If:

then \( x[n] \) is replaced by the local median \( m \).

Python Code:

def hampel_filter(trace, window_size=10, threshold=3):
    x = trace.data
    n = len(x)
    new_x = x.copy()

    k = window_size
    for i in range(k, n - k):
        window = x[i - k:i + k + 1]
        median = np.median(window)
        mad = np.median(np.abs(window - median))
        if np.abs(x[i] - median) > threshold * 1.4826 * mad:
            new_x[i] = median

    trace.data = new_x
    return trace

Pros:
- Excellent for spike removal without affecting underlying trends
- Robust to outliers

Tips:
- Use before smoothing or normalization
- Adjust window_size and threshold depending on noise level

Summary Update