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.

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.

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, we can directly run the Python Script assign_gaussian_noise.py as following:


def assign_noise_to_src_rec_file(in_fname, out_fname, noise_level=0.1):
    sr = SrcRec.read(in_fname)
    sr.add_noise(noise_level)
    sr.write(out_fname)
 
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_noisy.dat" # output source receiver file
    sigma = 0.1 # noise level in seconds
    assign_noise_to_src_rec_file(in_fname, out_fname, noise_level=sigma)

This script will generate the noisy data file OUTPUT_FILES/OUTPUT_FILES_signal/src_rec_file_forward_noisy.dat based on the noise-free data file OUTPUT_FILES/OUTPUT_FILES_signal/src_rec_file_forward.dat