Preparations for inversion

This part prepares necessary files to run this examples, including:

Download input file for source and receiver locations raw link
Build a 1-D isotropic initial model.
Build a checkerboard true model by adding perturbations to the initial model.
Calculate the traveltimes in the true model (checkerboard model) as the observed traveltimes.
Add Gaussian noise on the calculated traveltimes to mimic picking uncertainty.
Add random deviations following the uniform distribution on the earthquake locations and origin times as the initial source information.

Users can directly run the following codes to generate all necessary files automatically.


python prepare_input_files.py
mpirun -n 8 TOMOATT -i 3_input_params/input_params_signal.yaml
python assign_gaussian_noise.py

Here, we demonstrate the preparation step-by-step. First, please load necessary modules in Python Script:


import os
import numpy as np
import sys
from pytomoatt.model import ATTModel
from pytomoatt.checkerboard import Checker
from pytomoatt.src_rec import SrcRec

Prepare the initial model

The initial model is built as an isotropic 1D model, satisfying:

a constant velocity of 5.0 km/s above 0 km depth.
a constant velocity gradient from 5.0 km/s at 0 km to 8.0 km/s at 40 km.
a constant velocity of 8.0 km/s below 40 km depth.

Users can build the initial model by


class BuildInitialModel():
    def __init__(self, par_file="./3_input_params/input_params_signal.yaml", output_dir="2_models"):
        """
        Build initial model for tomography inversion
        """
        self.am = ATTModel(par_file)
        self.output_dir = output_dir
 
    def build_initial_model(self, vel_min=5.0, vel_max=8.0):
        """
        Build initial model for tomography inversion
        """
        self.am.vel[self.am.depths < 0, :, :] = vel_min
        idx = np.where((0 <= self.am.depths) & (self.am.depths < 40.0))[0]
        self.am.vel[idx, :, :] = np.linspace(vel_min, vel_max, idx.size)[::-1][:, np.newaxis, np.newaxis]
        self.am.vel[self.am.depths >= 40.0, :, :] = vel_max
 
output_dir = "2_models"
os.makedirs(output_dir, exist_ok=True)
bim = BuildInitialModel(output_dir=output_dir)
bim.build_initial_model()

The initial model file is generated as 2_models/model_init_N61_61_61.h5

Prepare the true model (checkerboard model)

The true model is built by adding perturbations to the initial model:

2 * 2 * 2 positive and negative 20% velocity perturbations within the range of 0.5-1.5 in longitude, 0.5-1.5 in latitude, and 0-40 in depth.
2 * 2 * 2 azimuthal anisotropy with magnitude of 0.1 and orthogonal fast velocity directions.

Users can build the true model by


def build_ckb_model(output_dir="2_models"):
    cbk = Checker(f'{output_dir}/model_init_N61_61_61.h5', para_fname="./3_input_params/input_params_signal.yaml")
    cbk.checkerboard(
        n_pert_x=2, n_pert_y=2, n_pert_z=2,
        pert_vel=0.2, pert_ani=0.1, ani_dir=60.0,
        lim_x=[0.5, 1.5], lim_y=[0.5, 1.5], lim_z=[0, 40]
    )
    cbk.write(f'{output_dir}/model_ckb_N61_61_61.h5')
 
# build ckb model
build_ckb_model(output_dir)

The true model file is generated as 2_models/model_ckb_N61_61_61.h5

n_pert_x, n_pert_y, n_pert_z indicate the number of checkers (perturbation block) in longitude, latitude, and depth.
pert_vel is the magnitude of assigned velocity perturbation, e.g., pert_vel=0.2 means a 20% velocity perturbation.
pert_ani is the magnitude of assigned azimuthal anisotropy.
lim_x, lim_y, lim_z mean the range of the perturbation block in longitude, latitude, and depth.

Prepare observed traveltime file

The observed traveltimes are generated by calculating the traveltimes in the true model (checkerboard model). Then, random Gaussian noise with a standard deviation of 0.1 s is assigned to the calculated traveltimes to mimic picking uncertainty. In addition, random deviations following the uniform distribution are added on the earthquake locations and origin times as the initial source information.

First, we can download the data file here , by using


url = 'https://zenodo.org/records/14053821/files/src_rec_config.dat'
path = "1_src_rec_files/src_rec_config.dat"
os.makedirs(os.path.dirname(path), exist_ok=True)
if not os.path.exists(path):
    sr = SrcRec.read(url)
    sr.write(path)

Then, we can run TomoATT in bash to calculate the traveltimes in the true model


mpirun -n 8 TOMOATT -i 3_input_params/input_params_signal.yaml

Here we use 8 processors in parallel.

If users want to change the number of used processors, please correspondingly change the parallel setting n_sims in the YAML file 3_input_params/input_params_signal.yaml:


parallel: # parameters for parallel computation
    n_sims: 8 # number of simultaneous runs (parallel the sources)

After calculating the traveltimes, we can find the calculated data file in OUTPUT_FILES/OUTPUT_FILES_signal/src_rec_file_forward.dat

To add Gaussian noise as well as deviate earthquake locations and origin times, we can directly run the Python Script assign_gaussian_noise.py as following:


class AssignNoise:
    def __init__(self, in_fname, out_fname):
        self.in_fname = in_fname
        self.out_fname = out_fname
        self.sr = SrcRec.read(self.in_fname)
 
    def assign_noise_for_tt(self, noise_level=0.1):
        self.sr.add_noise(noise_level)
 
    def assign_noise_for_src(self, lat_pert=0.1, lon_pert=0.1, dep_pert=10, tau_pert=0.5):
        self.sr.add_noise_to_source(lat_pert, lon_pert, dep_pert, tau_pert)
 
if __name__ == "__main__":
    in_fname = "OUTPUT_FILES/OUTPUT_FILES_signal/src_rec_file_forward.dat" # input source receiver file
    out_fname = "OUTPUT_FILES/OUTPUT_FILES_signal/src_rec_file_forward_errloc.dat" # output source receiver file
    sigma = 0.1 # noise level in seconds
    lat_pert = 0.1 # assign noise for latitude in degrees
    lon_pert = 0.1 # assign noise for longitude in degrees
    dep_pert = 10 # assign noise for depth in km
    tau_pert = 0.5 # assign noise for origin time in seconds
    # Initialize the instance
    an = AssignNoise(in_fname, out_fname)
    # Assign noise for travel time
    an.assign_noise_for_tt(sigma)
    # Assign noise for source
    an.assign_noise_for_src(lat_pert, lon_pert, dep_pert, tau_pert)
    # Write the output file
    an.sr.write(out_fname)

This script will generate the noisy data file with deviated locations OUTPUT_FILES/OUTPUT_FILES_signal/src_rec_file_forward_errloc.dat based on the noise-free data file with true locations OUTPUT_FILES/OUTPUT_FILES_signal/src_rec_file_forward.dat

When deviating the earthquake location, only the latitude, longitude, and depths are deviated in the DATA file.

However, when deviating the origin time, the traveltime should be negatively deviated correspondingly, as the factual arrival time (absolute picking time) remain unchanged.