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!