pysmo.tools.noise

Generate realistic synthetic noise.

matches the naturally observed amplitude spectrum.

Examples:

Given the spectral amplitude in observed seismic noise on Earth is not flat (i.e. not consisting of 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 NHNM (red), NLNM (blue), and an interpolated model that lies between the two (green).

peterson example

Example source code

peterson.py

"""
Example script for pysmo.tools.noise
"""

#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from pysmo.tools.noise import generate_noise, peterson
from pandas import Timedelta


def calc_power(
    data: np.ndarray, sampling_frequency: float, nperseg: int, nfft: int
) -> tuple[np.ndarray, np.ndarray]:
    """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:]


def main() -> None:
    # Set parameters
    npts: int = 200000  # multiple of 4
    delta = Timedelta(seconds=0.1)
    sampling_frequency = 1 / delta.total_seconds()
    nperseg = int(npts / 4)
    nfft = int(npts / 2)
    time_in_seconds = np.linspace(0, npts * delta.total_seconds(), npts)

    # 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, nperseg, nfft
    )
    f_mid, Pxx_dens_mid = calc_power(
        mid_noise_seismogram.data, sampling_frequency, nperseg, nfft
    )
    f_high, Pxx_dens_high = calc_power(
        high_noise_seismogram.data, sampling_frequency, nperseg, nfft
    )

    _ = plt.figure(figsize=(13, 9), layout="tight")

    # Plot random high and low noise seismograms
    ax1 = plt.subplot2grid((4, 3), (0, 0))
    plt.plot(time_in_seconds, high_noise_seismogram.data, "r", linewidth=0.2)
    plt.ylabel("Ground Accelaration")
    plt.xlabel("Time [s]")
    plt.xlim(time_in_seconds[0], time_in_seconds[-1])
    xticks = np.arange(
        tmin := time_in_seconds[0], (tmax := time_in_seconds[-1]) + 1, (tmax - tmin) / 5
    )
    ax1.set_xticks(xticks)

    ax2 = plt.subplot2grid((4, 3), (0, 1))
    ax2.set_xticks(xticks)
    plt.plot(time_in_seconds, mid_noise_seismogram.data, "g", linewidth=0.2)
    plt.xlabel("Time [s]")
    plt.xlim(time_in_seconds[0], time_in_seconds[-1])

    ax3 = plt.subplot2grid((4, 3), (0, 2))
    ax3.set_xticks(xticks)
    plt.plot(time_in_seconds, low_noise_seismogram.data, "b", linewidth=0.2)
    plt.xlabel("Time [s]")
    plt.xlim(time_in_seconds[0], time_in_seconds[-1])

    # Plot PSD of noise
    _ = plt.subplot2grid((4, 3), (1, 0), rowspan=3, colspan=3)
    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.xlim(low_noise_model.T[0], low_noise_model.T[-1])
    plt.xlabel("Period [s]")
    plt.ylabel("Power Spectral Density (dB/Hz)")
    plt.legend()
    plt.savefig("peterson.png", transparent=True)
    plt.show()


if __name__ == "__main__":
    main()

Classes:

Name	Description
`NoiseModel`	Class to store seismic noise models.

Functions:

Name	Description
`generate_noise`	Generate a random seismogram from a noise model. Realistic seismic noise is
`peterson`	Generate a noise model by interpolating between Peterson's[^1]

NoiseModel `dataclass`

Class to store seismic noise models.

Parameters:

Name	Type	Description	Default
`psd`	`ndarray`	Power spectral density of ground acceleration [dB].	`(lambda: array([]))()`
`T`	`ndarray`	Period [seconds].	`(lambda: array([]))()`

Source code in src/pysmo/tools/noise.py

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

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

    psd: np.ndarray = field(default_factory=lambda: np.array([]))
    T: np.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: Timestamp = 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:

Name	Type	Description	Default
`model`	`NoiseModel`	Noise model used to compute seismic noise.	required
`npts`	`int`	Number of points of generated noise	required
`delta`	`Timedelta`	Sampling interval of generated noise	`delta.value`
`begin_time`	`Timestamp`	Seismogram begin time	`begin_time.value`
`return_velocity`	`bool`	Return velocity instead of acceleration.	`False`
`seed`	`int \| None`	set random seed manually (usually for testing purposes).	`None`

Returns:

Type	Description
`MiniSeismogram`	Seismogram with random seismic noise as data.

Source code in src/pysmo/tools/noise.py

def generate_noise(
    model: NoiseModel,
    npts: int,
    delta: Timedelta = SEISMOGRAM_DEFAULTS.delta.value,
    begin_time: Timestamp = 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.

    Args:
        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).

Peterson, Jon R. Observations and Modeling of Seismic Background Noise. Report, 93–322, 1993, https://doi.org/10.3133/ofr93322. USGS Publications Warehouse. ↩

Parameters:

Name	Type	Description	Default
`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.	required

Returns:

Type	Description
`NoiseModel`	Noise model.

Source code in src/pysmo/tools/noise.py

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

    [^1]: Peterson, Jon R. Observations and Modeling of Seismic Background
        Noise. Report, 93–322, 1993, https://doi.org/10.3133/ofr93322. USGS
        Publications Warehouse.

    Args:
        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`