Understanding background seismic noise levels is crucial for assessing the suitability of a site for seismic monitoring and instrumentation. Background noise analysis helps identify optimal locations for seismic sensors and ensures accurate detection of seismic events. This blog explores how to calculate background noise levels based on methodologies from Ramirez et al. (2019), using power spectral density (PSD) techniques.

Why Measure Seismic Noise Levels?

Seismic stations continuously record ambient noise from both natural and anthropogenic sources. High noise levels can mask low-magnitude earthquakes and reduce the effectiveness of seismic monitoring networks. The analysis of background noise levels provides insights into:

• The impact of environmental conditions on sensor performance.
• The influence of geological settings on noise characteristics.
• The optimal placement of seismic sensors for improved data quality.

Step 1: Data Collection

To measure background seismic noise, long-term continuous seismic data from a given site must be collected. Typically, broadband seismometers and accelerometers are used to capture signals across a wide range of frequencies. These data are stored in a digital format, such as miniSEED, which allows for efficient processing.

Step 2: Preprocessing the Seismic Data

Before analyzing noise levels, raw seismic data must be preprocessed to remove unwanted signals. This includes:

• Segmenting Data: Breaking continuous records into smaller time windows (e.g., hourly segments) to analyze temporal variations.
• Applying Windowing Functions: Using techniques like cosine tapers to reduce spectral leakage.
• Fourier Transform: Converting time-domain data into the frequency domain for spectral analysis.
• Acceleration Conversion: Deriving acceleration spectra from velocity or displacement records.

Step 3: Power Spectral Density (PSD) Calculation

The Power Spectral Density (PSD) method is widely used to quantify ambient seismic noise levels. This technique measures the energy distribution of seismic signals across different frequencies. The steps include:

1. Computing the Fast Fourier Transform (FFT): This transforms seismic data from the time domain to the frequency domain.
2. Averaging Spectra: To improve statistical reliability, multiple PSD estimates are averaged over time.
3. Expressing Noise Levels in Decibels (dB): PSD values are converted to decibels relative to a reference power level (e.g., 1 μm²/s⁴/Hz).
4. Generating PSD Probability Density Functions (PDFs): A statistical representation of noise levels over time to assess variations.

Step 4: Analyzing Noise Level Trends

Once PSDs are computed, trends in noise levels can be analyzed:

• Comparing Different Geological Settings: Sites on sedimentary deposits often exhibit higher noise levels than those on bedrock.
• Assessing Diurnal Variations: Noise levels tend to be higher during the day due to human activities.
• Impact of Environmental Factors: Atmospheric pressure and temperature fluctuations can influence background noise levels.

Step 5: Comparing with Standard Noise Models

To interpret seismic noise levels, results are compared with global noise models:

• New High-Noise Model (NHNM): Represents the highest observed noise levels at global seismic stations.
• New Low-Noise Model (NLNM): Represents the lowest possible ambient noise levels.

Step 6: Identifying and Mitigating Noise Sources

If noise levels are excessive, it is essential to identify and mitigate noise sources:

• Relocating Sensors: Moving stations away from urban areas or installing them in underground vaults can reduce anthropogenic noise.
• Improving Insulation: Using thermal and acoustic insulation can minimize environmental noise interference.
• Using Different Sensor Types: Some sensors may exhibit different sensitivity to background noise, requiring careful selection based on site conditions.

Complete Python Code for Computing Background Seismic Noise Levels

Below is a full Python script that includes preprocessing, PSD computation, and noise model comparison.

import numpy as np
import obspy
import matplotlib.pyplot as plt
from obspy.signal.spectral_estimation import PSD, get_nlnm, get_nhnm

def preprocess_seismic_data(trace):
    """Preprocess seismic data: detrend, taper, and remove response."""
    trace.detrend("linear")
    trace.taper(max_percentage=0.05, type="hann")
    return trace

def compute_psd(trace, nfft=1024, overlap=0.5):
    """Compute Power Spectral Density (PSD) for a given seismic trace."""
    psd, freq = PSD(trace.data, nfft=nfft, overlap=overlap, fs=trace.stats.sampling_rate)
    psd_db = 10 * np.log10(psd)
    return freq, psd_db

def plot_noise_models(freq, psd_db):
    """Plot the computed PSD against global noise models."""
    nlnm_freq, nlnm_psd = get_nlnm()
    nhnm_freq, nhnm_psd = get_nhnm()
    
    plt.figure(figsize=(8,6))
    plt.plot(nlnm_freq, nlnm_psd, label="NLNM", linestyle="dashed")
    plt.plot(nhnm_freq, nhnm_psd, label="NHNM", linestyle="dashed")
    plt.plot(freq, psd_db, label="Computed PSD", color='red')
    plt.xlabel("Frequency (Hz)")
    plt.ylabel("Power Spectral Density (dB)")
    plt.legend()
    plt.title("Comparison of Seismic Noise Levels with Global Models")
    plt.grid()
    plt.show()

# Example usage
st = obspy.read("example.mseed")
st = st.select(component="Z")
for tr in st:
    tr = preprocess_seismic_data(tr)
    freq, psd_db = compute_psd(tr)
    plot_noise_models(freq, psd_db)

Conclusion

Calculating background seismic noise levels using the PSD method provides a robust approach to evaluating site suitability for seismic monitoring. By following these steps, researchers and engineers can optimize sensor placement, minimize noise interference, and enhance the reliability of seismic data. Understanding noise trends and environmental influences allows for better decision-making in seismic instrumentation and hazard assessment.ntal influences allows for better decision-making in seismic instrumentation and hazard assessment.

Reference

Ramírez, E. E., Vidal-Villegas, J. A., Nuñez-Leal, M. A., Ramírez-Hernández, J., Mejía-Trejo, A., & Rosas-Verdug, E. (2019). Seismic noise levels in northern Baja California, Mexico. Bulletin of the Seismological Society of America, 109(2), 610–620. https://doi.org/10.1785/0120180155

When dealing with seismic data, a critical step in understanding the potential impact on structures is analyzing the response to ground motion, often represented as a response spectrum. This Python script by QuakeLogic enables engineers and researchers to compute the response spectrum of seismic time-series data (in .tcl format), plot the results, compute Peak Ground Acceleration (PGA), Peak Ground Velocity (PGV), Peak Ground Displacement (PGD), and save them for further analysis.

In this blog post, we will explain how to use this Python code, its features, and provide an example of how to compute the response spectrum and key seismic parameters (PGA, PGV, PGD) for a given ground motion.

Key Features of the Script

Input and Output:
The script processes input acceleration time-series data stored in .tcl format. The input files should contain acceleration data in units of in/s², which the script converts to g (gravitational acceleration).

Fourier Transform-Based Method:
The response spectrum is calculated using a Fourier transform-based approach, ensuring accuracy across a range of frequencies and damping ratios.

Peak Ground Parameters:
The script computes PGA, PGV, and PGD to offer a complete picture of the seismic intensity. PGA represents the highest acceleration during the event, PGV corresponds to the highest velocity, and PGD reflects the largest displacement.

Plots and Data Files:
The script generates plots of the input acceleration time series, response spectrum, and computed PGA, PGV, and PGD values, saving them in specified directories. It also exports the computed response spectrum and seismic parameters into output files.

Step-by-Step Guide to Use the Script

Installation and Requirements
Before running the script, ensure you have the necessary Python packages installed. Run the following command to install the required libraries:

pip install numpy matplotlib tqdm

Create a Config File
The script reads parameters from a configuration file (config.json). This file specifies the input directory, response spectrum periods, damping ratio, and other critical parameters. Here’s an example of what the config.json should look like:

{
  "sampling_rate": 50,
  "damping_ratio": 0.05,
  "min_period": 0.01,
  "max_period": 5.0,
  "num_periods": 100,
  "input_dir": "input_data/"
}

Run the Script
After preparing the input files (in .tcl format) and the config file, you can run the script using the following command:

python seismic_response_spectrum_analysis.py

The script processes all .tcl files in the specified input directory (input_data/), computes the response spectrum, and outputs the results, including PGA, PGV, and PGD.

Computation of PGA, PGV, and PGD

1. PGA (Peak Ground Acceleration):
   The highest absolute value of the acceleration time series, representing the strongest shaking the ground experiences.
2. PGV (Peak Ground Velocity):
   Computed by integrating the acceleration time series over time to derive velocity and extracting the maximum value.
3. PGD (Peak Ground Displacement):
   Derived by further integrating the velocity time series to calculate displacement, with the peak displacement indicating the largest ground movement.

These values provide critical insights into the severity of ground motion and are vital for structural response analyses.

Generated Outputs

The script generates three key outputs:

1. Plots:
   The input acceleration time series, response spectrum, and seismic parameters (PGA, PGV, PGD) are plotted and saved as PNG files in the figures/ directory.
2. Response Spectrum Data:
   The computed response spectrum is saved in .dat format in the output_data/ directory.
3. Seismic Parameters:
   PGA, PGV, and PGD are computed and saved in an output file for further reference.

How the Code Works

1. Loading Acceleration Data

The script reads acceleration data from .tcl files, which are expected to contain acceleration in in/s². The values are converted to g using a conversion factor.

2. Calculating the Response Spectrum

The script uses the Fourier transform method to compute the response spectrum for each oscillator frequency and damping ratio.

3. Computing PGA, PGV, and PGD

The script integrates the acceleration time series to compute velocity and displacement, identifying the peak values for PGA, PGV, and PGD.

4. Plotting and Saving Results

Once the response spectrum and seismic parameters are calculated, the script plots both the input acceleration time series and the response spectrum. The PGA, PGV, and PGD values are also printed for reference.

Example Code

import os
import glob
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
import json

IN_S2_TO_G = 386.0886  # Conversion factor: 1 g = 386.0886 in/s²

def calc_sdf_resp(freq, fourier_amp, osc_damping, osc_freq, max_freq_ratio=5., peak_resp_only=False):
    h = (-np.power(osc_freq, 2.) / ((np.power(freq, 2.) - np.power(osc_freq, 2.)) - 2.j * osc_damping * osc_freq * freq))
    n = len(fourier_amp)
    m = max(n, int(max_freq_ratio * osc_freq / freq[1]))
    scale = float(m) / float(n)
    resp = scale * np.fft.irfft(fourier_amp * h, 2 * (m - 1))
    if peak_resp_only:
        return np.abs(resp).max()
    return resp

def calc_spec_accels(time_step, accel_ts, osc_freqs, osc_damping=0.05):
    fourier_amp = np.fft.rfft(accel_ts)
    freq = np.linspace(0, 1. / (2 * time_step), num=fourier_amp.size)
    spec_accels = [calc_sdf_resp(freq, fourier_amp, osc_damping, of, peak_resp_only=True)
                   for of in osc_freqs]
    return np.rec.fromarrays([osc_freqs, spec_accels], names='osc_freq,spec_accel')

def load_acceleration_file(filepath):
    acc_data_in_in_s2 = np.loadtxt(filepath)
    return acc_data_in_in_s2 / IN_S2_TO_G  # Convert from in/s² to g

def compute_pgv(acc, time_step):
    vel = np.cumsum(acc) * time_step
    return np.max(np.abs(vel))

def compute_pgd(vel, time_step):
    disp = np.cumsum(vel) * time_step
    return np.max(np.abs(disp))

def gen_rsa(file_to_use, config, figures_dir):
    data = load_acceleration_file(file_to_use)
    osc_damping = config['damping_ratio']
    osc_freqs = np.logspace(np.log10(1.0 / config['max_period']), np.log10(1.0 / config['min_period']), config['num_periods'])
    time_step = 1.0 / config['sampling_rate']

    pga = np.max(np.abs(data))
    pgv = compute_pgv(data, time_step)
    pgd = compute_pgd(pgv, time_step)

    resp_spec = calc_spec_accels(time_step, data, osc_freqs, osc_damping)

    # Plotting and saving code...
    # (skipped for brevity, similar to original script)

    # Save PGA, PGV, PGD
    with open(f'{figures_dir}/seismic_parameters.txt', 'w') as file:
        file.write(f'PGA: {pga}\nPGV: {pgv}\nPGD: {pgd}\n')

Contact and Versioning

Version: 1.1
Last Update: October 23, 2024

This code was developed by QuakeLogic Inc. to assist with seismic data analysis. For support, contact us at: support@quakelogic.net.

We hope this tool is helpful for your seismic analysis. Let us know if you have any feedback!