pysmo.tools.noise #

Generate realistic synthetic noise.

This module provides support for calculating random synthetic noise that matches the naturally observed amplitude spectrum.

Further Reading

Peterson, J., 1993. Observations and modelling of background seismic noise. Open-file report 93-322, U. S. Geological Survey, Albuquerque, New Mexico.

Examples:

Given the spectral amplitude in observed seismic noise on Earth is not flat (i.e. white noise), it makes sense to calculate more realistic noise for things like resolution tests with synthetic data.

In this example, random noise seismograms are generated from three different noise models. These are Peterson's NLNM, NHNM, and an interpolated model that lies between the two.

peterson

Example source code

peterson.py

#!/usr/bin/env python

# Example script for pysmo.tools.noise
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from pysmo.tools.noise import generate_noise, peterson
from datetime import timedelta

npts = 200000
delta = timedelta(seconds=0.1)
sampling_frequency = 1 / delta.total_seconds()
nperseg = npts / 4
nfft = npts / 2
time_in_seconds = np.linspace(0, npts * delta.total_seconds(), npts)


def calc_power(data, sampling_frequency):
    """Calculuate power and drop first element (f=0Hz) to avoid dividing by 0"""
    freqs, psd = signal.welch(
        data, sampling_frequency, nperseg=nperseg, nfft=nfft, scaling="density"
    )
    return freqs[1:], psd[1:]


# Calculate noise models
low_noise_model = peterson(noise_level=0)
mid_noise_model = peterson(noise_level=0.5)
high_noise_model = peterson(noise_level=1)

# Generate random noise seismograms
low_noise_seismogram = generate_noise(npts=npts, model=low_noise_model, delta=delta)
mid_noise_seismogram = generate_noise(npts=npts, model=mid_noise_model, delta=delta)
high_noise_seismogram = generate_noise(npts=npts, model=high_noise_model, delta=delta)

# Calculuate power spectral density
f_low, Pxx_dens_low = calc_power(low_noise_seismogram.data, sampling_frequency)
f_mid, Pxx_dens_mid = calc_power(mid_noise_seismogram.data, sampling_frequency)
f_high, Pxx_dens_high = calc_power(high_noise_seismogram.data, sampling_frequency)

fig = plt.figure(figsize=(20, 10))
# Plot random high and low noise
plt.subplot(321)
plt.plot(time_in_seconds, high_noise_seismogram.data, "r", linewidth=0.2)
plt.ylabel("Ground Accelaration")

plt.subplot(323)
plt.plot(time_in_seconds, mid_noise_seismogram.data, "g", linewidth=0.2)
plt.ylabel("Ground Accelaration")

plt.subplot(325)
plt.plot(time_in_seconds, low_noise_seismogram.data, "b", linewidth=0.2)
plt.ylabel("Ground Accelaration")
plt.xlabel("Time [s]")

# Plot PSD of noise
plt.subplot(122)
plt.plot(
    1 / f_high,
    10 * np.log10(Pxx_dens_high),
    "r",
    linewidth=0.5,
    label="generated high noise",
)
plt.plot(
    1 / f_mid,
    10 * np.log10(Pxx_dens_mid),
    "g",
    linewidth=0.5,
    label="generated medium noise",
)
plt.plot(
    1 / f_low,
    10 * np.log10(Pxx_dens_low),
    "b",
    linewidth=0.5,
    label="generated low noise",
)
plt.plot(
    high_noise_model.T,
    high_noise_model.psd,
    "k",
    linewidth=1,
    linestyle="dotted",
    label="NHNM",
)
plt.plot(
    mid_noise_model.T,
    mid_noise_model.psd,
    "k",
    linewidth=1,
    linestyle="dashdot",
    label="Interpolated noise model",
)
plt.plot(
    low_noise_model.T,
    low_noise_model.psd,
    "k",
    linewidth=1,
    linestyle="dashed",
    label="NLNM",
)
plt.gca().set_xscale("log")
plt.xlabel("Period [s]")
plt.ylabel("Power Spectral Density (dB/Hz)")
plt.legend()
plt.savefig("peterson.png")
plt.show()

Classes:

NoiseModel –

Class to store seismic noise models.

Functions:

generate_noise –

Generate a random seismogram from a noise model. Realistic seismic noise is
peterson –

Generate a noise model by interpolating between Peterson's

NoiseModel `dataclass` #

Class to store seismic noise models.

Attributes:

psd (NDArray) –

Power spectral density of ground acceleration [dB].
T (NDArray) –

Period [seconds].

Source code in pysmo/tools/noise.py

@dataclass(frozen=True)
class NoiseModel:
    """Class to store seismic noise models.

    Attributes:
        psd: Power spectral density of ground acceleration [dB].
        T: Period [seconds].
    """

    psd: npt.NDArray = field(default_factory=lambda: np.array([]))
    T: npt.NDArray = field(default_factory=lambda: np.array([]))

    def __post_init__(self) -> None:
        if np.size(self.psd) != np.size(self.T):
            raise ValueError(
                "dB and T arrays are not of same size",
                f"({np.size(self.psd)} != {np.size(self.T)}",
            )
        self.psd.flags.writeable = False
        self.T.flags.writeable = False

generate_noise #

generate_noise(model: NoiseModel, npts: int, delta: timedelta = delta.value, begin_time: datetime = begin_time.value, return_velocity: bool = False, seed: int | None = None) -> MiniSeismogram

Generate a random seismogram from a noise model. Realistic seismic noise is generated by assigning random phases to a fixed amplitude spectrum prescribed by a noise model.

Parameters:

model (NoiseModel) –

Noise model used to compute seismic noise.
npts (int) –

Number of points of generated noise
delta (timedelta, default: delta.value ) –

Sampling interval of generated noise
begin_time (datetime, default: begin_time.value ) –

Seismogram begin time
return_velocity (bool, default: False ) –

Return velocity instead of acceleration.
seed (int | None, default: None ) –

set random seed manually (usually for testing purposes).

Returns:

MiniSeismogram –

Seismogram with random seismic noise as data.

Source code in pysmo/tools/noise.py

def generate_noise(
    model: NoiseModel,
    npts: int,
    delta: timedelta = SEISMOGRAM_DEFAULTS.delta.value,
    begin_time: datetime = SEISMOGRAM_DEFAULTS.begin_time.value,
    return_velocity: bool = False,
    seed: int | None = None,
) -> MiniSeismogram:
    """Generate a random seismogram from a noise model. Realistic seismic noise is
    generated by assigning random phases to a fixed amplitude spectrum prescribed by a
    noise model.

    Parameters:
        model: Noise model used to compute seismic noise.
        npts: Number of points of generated noise
        delta: Sampling interval of generated noise
        begin_time: Seismogram begin time
        return_velocity: Return velocity instead of acceleration.
        seed: set random seed manually (usually for testing purposes).

    Returns:
        Seismogram with random seismic noise as data.
    """
    # Sampling frequency
    Fs = 1 / delta.total_seconds()

    # Nyquist frequency
    Fnyq = 0.5 / delta.total_seconds()

    # get next power of 2 of the nunmber of points and calculate frequencies from
    # Fs/NPTS to Fnyq (we skip a frequency of 0 for now to avoid dividing by 0)
    NPTS = int(2 ** np.ceil(np.log2(npts)))
    freqs = np.linspace(Fs / NPTS, Fnyq, NPTS - 1)

    # interpolate psd and recreate amplitude spectrum with the first
    # term=0 (i.e. mean=0).
    Pxx = np.interp(1 / freqs, model.T, model.psd)
    spectrum = np.append(
        np.array([0]), np.sqrt(10 ** (Pxx / 10) * NPTS / delta.total_seconds() * 2)
    )

    # phase is randomly generated
    rng = np.random.default_rng(seed)
    phase = (rng.random(NPTS) - 0.5) * np.pi * 2

    # combine amplitude with phase and perform ifft
    NewX = spectrum * (np.cos(phase) + 1j * np.sin(phase))
    acceleration = np.fft.irfft(NewX)

    start = int((NPTS - npts) / 2)
    end = start + npts
    if return_velocity:
        velocity = cumulative_trapezoid(acceleration, dx=delta.total_seconds())
        velocity = velocity[start:end]
        return MiniSeismogram(begin_time=begin_time, delta=delta, data=velocity)
    acceleration = acceleration[start:end]
    return MiniSeismogram(begin_time=begin_time, delta=delta, data=acceleration)

peterson #

peterson(noise_level: float) -> NoiseModel

Generate a noise model by interpolating between Peterson's New Low Noise Model (NLNM) and New High Noice Model (NHNM).

Parameters:

noise_level (float) –

Determines the noise level of the generated noise model. A noise level of 0 returns the NLNM, 1 returns the NHNM, and anything > 0 and < 1 returns an interpolated model that lies between the NLNM and NHNM.

Returns:

NoiseModel –

Noise model.

Source code in pysmo/tools/noise.py

def peterson(noise_level: float) -> NoiseModel:
    """Generate a noise model by interpolating between Peterson's
    New Low Noise Model (NLNM) and New High Noice Model (NHNM).

    Parameters:
        noise_level: Determines the noise level of the generated noise model.
                     A noise level of 0 returns the NLNM, 1 returns the NHNM,
                     and anything > 0 and < 1 returns an interpolated model
                     that lies between the NLNM and NHNM.

    Returns:
        Noise model.
    """
    # check for valid input
    if not 0 <= noise_level <= 1:
        raise ValueError(
            f"Parameter noise_level={noise_level} is not within 0-1 range."
        )

    # calculate noise model
    if noise_level == 0:
        return NLNM
    elif noise_level == 1:
        return NHNM
    else:
        T_common = np.unique(np.concatenate((NLNM.T, NHNM.T)))
        NLNM_interp = np.interp(T_common, NLNM.T, NLNM.psd)
        NHNM_interp = np.interp(T_common, NHNM.T, NHNM.psd)
        dB = NLNM_interp + (NHNM_interp - NLNM_interp) * noise_level
        return NoiseModel(psd=dB, T=T_common)

pysmo.tools.noise #

NoiseModel dataclass #

generate_noise #

peterson #

NoiseModel `dataclass` #