Tutorial for waveform inversion using Kibrary

Note: This tutorial uses the "anselme" branch of the Kibrary project.
The data, property files (.properties), and scripts for this tutorial can be found in Kibrary/tutorial in the Kibrary project on Github, anselme branch.

1. Prerequisite

Check you have Java 8 installed on your system
Follow the installation instructions for Kibrary (anselme branch). Alternatively, you can download the executable JAR file of Kibrary.
Follow the installation instructions for TomoTool. Alternatively, you can download the executable JAR file of TomoTool
Download shbpcat, and psvbpcat
Dowload tish, tipsv, shfp, and psvfp
Download the seismic data (SAC files for event 031704A) here (or from the tutorial folder). Create the folder tutorial_inversion, and move the data folder 031704A you just downloaded to tutorial_inversion. tutorial_inversion is the root folder for this tutorial.

Optional - for visualization

Install anaconda
Create a conda environment tuto and install the following python packages by running

conda create -n tuto python=3.8
conda install -n tuto numpy matplotlib cartopy -y

Before running the python scripts in this tutorial, activate the tuto environment using

conda activate tuto

2. Pre-tests (optional but recommended)

SHBPCAT

In folder shbpcat:

make
cd examples
cp -r AK135 AK135_TEST
mpirun -n 3 ../mpi-shbpcat < AK135_SH_5.inf

The run should complete in less than 1 minute. You can compare the results with pre-computed files:

head -n 5 XY01.340A_TA.PB...SH.spc.txt <(java io.github.kensuke1984.kibrary.quick.LookAtBPspc XY01.340A_TA.PB...SH.spc)

Note: There might be small differences, but they should be within ~10^-8 %

PSVBPCAT

In folder psvbpcat:

make
cd examples
cp -r AK135 AK135_TEST
mpirun -n 1 ../mpi-psvbpcat < AK135_PSV.inf

The run should complete in less than 2 minute. You can compare the results with pre-computed files:

head -n 5 XY01.340A_TA.PB...SH.spc.txt <(java io.github.kensuke1984.kibrary.quick.LookAtBPspc XY01.340A_TA.PB...SH.spc)

Note: There might be small differences, but they should be within ~10^-8 %

3. Computation of synthetics

In folder tutorial_inversion, run:

java io.github.kensuke1984.kibrary.Property

When prompted to, enter 19 (SyntheticDSMInformationFileMaker). Rename the output file to sdsmifm.properties.

structureFile PREM
tlen 3276.8
np 512

Run

java io.github.kensuke1984.kibrary.Operation sdsmifm.properties

Rename the output folder to synthetic

cd synthetic/031704A
mpirun -n 8 mpi-tipsv < PREM_PSV
mpirun -n 8 mpi-tish < PREM_SH

java io.github.kensuke1984.kibrary.Property

When prompted to, enter 17 (SPC_SAC). Rename the output file to ss.properties.

components Z T
psvPath .
shPath .
sourceTimeFunction 2

Run

java io.github.kensuke1984.kibrary.Operation ss.properties

Rename the output folder to sac

java io.github.kensuke1984.kibrary.Property

When prompted to, enter 4 (FilterDivider) Rename the output file to fd.properties

components Z T
obsPath ..
synPath sac
lowFreq 0.005
highFreq 0.08
filter bandpass
np 4
backward false

Run

java io.github.kensuke1984.kibrary.Operation fd.properties

Rename the output folder to filtered

4. Time windows

In folder filtered, run

java io.github.kensuke1984.kibrary.timewindow.addons.MakeWindowPcPandScS

This will create 6 time windows binary files

timewindow_P.dat used for computing static corrections for the PcP phase
timewindow_S.dat same, but for the ScS phase
selectedTimewindow_PcP.dat contains selected windows for the PcP phase
selectedTimewindow_ScS.dat same, but for the ScS phase
timewindow_PcP.dat and timewindow_ScS.dat contain windows before selection, and are not used in this tutorial

The time window are selected based on the following criteria:

(zero-lag) correlation coefficient > -0.5
variance < 2.5
0.4 < amplitude ratio < 2.5
signal-to-noise ratio > 2

Currently, these criteria are hard-coded. Alternatively, you can use timewindow_PcP.dat and timewindow_ScS.dat together with io.github.kensuke1984.kibrary.selection.DataSelection.

5. Static corrections

In folder filtered, run

java io.github.kensuke1984.kibrary.Property

When prompted to, enter 6 (FujiStaticCorrection). Rename the output file to fsc.properties

To compute static corrections for the ScS phase, edit fsc.properties as follows:

components Z T
timewindowInformationPath timewindow_S.dat
convolute true

and run

java io.github.kensuke1984.kibrary.Operation fsc.properties

Rename the output file to staticCorrection_S.dat.

Do the same for the PcP phase, replacing

timewindowInformationPath timewindow_P.dat

in the fsc.properties file. Rename the output file to staticCorrection_P.dat

6. Mantle corrections

Corrections for the 3-D mantle above D" (the target region) are computed using TomoTool. To compute corrections for PcP-P differential travel times for the SEMUCB-WM1 model (French and Romanowicz, 2014), in folder filtered, run

java raytheory.ComputeCorrection selectedTimewindow_PcP.dat semucb pcp

This creates a binary correction file mantleCorrection_P-PcP.dat with the same structure as the one computed using FujiStaticCorrection (section 5)

For ScS-S differential travel times, run

java raytheory.ComputeCorrection selectedTimewindow_ScS.dat semucb scs

This creates the file mantleCorrection_S-ScS.dat.

7. Residual vector

In folder filtered, run

java io.github.kensuke1984.kibrary.Property

When prompted to, enter 9 (ObservedSyntheticDatasetMaker). Rename the output file to osdm.properties

components Z T
timeCorrection true
timewindowPath selectedTimewindow_ScS.dat
staticCorrectionPath staticCorrection_S.dat
correctMantle true
mantleCorrectionPath mantleCorrection_S-ScS.dat

Run

java io.github.kensuke1984.kibrary.Operation osdm.properties

Rename the output files to waveform_ScS.dat and waveformID_ScS.dat

Do the same for the PcP phase, replacing

timewindowPath selectedTimewindow_ScS.dat
staticCorrectionPath staticCorrection_P.dat
mantleCorrectionPath mantleCorrection_P-PcP.dat

in osdm.properties. Rename the output files to waveform_PcP.dat and waveformID_PcP.dat.

8. Partial derivatives

Voxel mesh

We first define the mesh of voxels (3-D pixels) to parametrize the target region in the inversion model. Download the files horizontalPoints_5deg.inf and radii.inf from the tutorial/INF_FILES folder into the tutorial_inversion folder:

horizontalPoints_5deg.inf: each line gives the latitude and longitude (in degree) of the center of a 5 by 5 degree pixel at the Earth's surface
radii.inf: each line gives the radius of the center of a depth layer

In tutorial_inversion, create the files perturbationPoint_5deg.inf and perturbationLayer.inf by running:

while read p; do while read r; do echo "$p $r" >> perturbationPoint_5deg.inf; done < radii.inf; done < horizontalPoints_5deg.inf
awk '{print $1,"50"}' radii.inf > perturbationLayer.inf

To create the unknown parameter file (i.e. model parameters) for the mesh perturbationPoint_5deg.inf for partial types MU and LAMBDA, run:

java io.github.kensuke1984.kibrary.inversion.addons.MakeUnknownParameterFile perturbationPoint_5deg.inf perturbationLayer.inf 5 MU LAMBDA

Rename the output file unknowns... to unknowns_5deg.inf

In order to achieve higher accuracy, we define a new finer mesh (1.25 by 1.25 degree) on which the partials will be computed, and then later integrated to a coarser mesh (2.5 by 2.5 degree). To make the finer mesh, run

java io.github.kensuke1984.kibrary.inversion.addons.ResampleGrid unknowns_5deg.inf

This creates the files

```newHorizontalPoints_sampled2.inf``: horizontal point file with voxel size of 2.5 by 2.5 degree
newHorizontalPoints_sampled3.inf: same, but with voxel size of 1.25 by 1.25 degree
horizontalMapping_sampled2.inf: file that maps voxels from newHorizontalPoints_sampled2.inf```` to horizontalPoints_5deg.inf```
horizontalMapping_sampled3.inf: file that maps voxels from newHorizontalPoints_sampled3.inf```` to horizontalPoints_5deg.inf```
horizontalMapping_sampled3.inf: file that maps voxels from newHorizontalPoints_sampled3.inf```` to newHorizontalPoints_sampled2.inf````

Creates the file pointLocations_1.25deg.inf and perturbationPoint_1.25deg.inf by running:

awk '{printf "%s ", $1}' radii.inf > tmp; echo "" >> tmp; cat newHorizontalPoints_sampled3.inf >> tmp; mv tmp pointLocations_1.25deg.inf
while read p; do while read r; do echo "$p $r" >> perturbationPoint_1.25deg.inf; done < radii.inf; done < newHorizontalPoints_sampled3.inf

Create the new unknown parameter file for the 1.25 degree mesh perturbationPoint_1.25deg.inf by running:

java io.github.kensuke1984.kibrary.inversion.addons.MakeUnknownParameterFile perturbationPoint_1.25deg.inf perturbationLayer.inf 1.25 MU LAMBDA

Rename the output file unknowns... to unknowns_1.25deg.inf.

Create the file perturbationPoint_2.5deg.inf by running:

while read p; do while read r; do echo "$p $r" >> perturbationPoint_2.5deg.inf; done < radii.inf; done < newHorizontalPoints_sampled2.inf

Create the new unknown parameter file for the 2.5 degree mesh perturbationPoint_2.5deg.inf by running:

java io.github.kensuke1984.kibrary.inversion.addons.MakeUnknownParameterFile perturbationPoint_2.5deg.inf perturbationLayer.inf 2.5 MU LAMBDA

Rename the output file unknowns... to unknowns_2.5deg.inf.

station.inf

In tutorial_inversion, run:

java io.github.kensuke1984.kibrary.inversion.StationInformationFile .

station.inf gives the station code, the network code, the latitude, and longitude. The 5 first lines are:

A06 XN 60.6312 -123.4932
GGN CN 45.1184 -66.842
A07 XN 60.9368 -123.0947
ICQ CN 49.5217 -67.2719
A04 XN 59.9811 -122.9164

DSM information files

In folder tutorial_inversion, run:

java io.github.kensuke1984.kibrary.Property

When prompted to, enter 7 (InformationFileMaker). Rename the output file to ifm.properties, and edit it as follows:

locationsPath pointLocations_1.25deg.inf
stationInformationPath station.inf
np 512
tlen 3276.8
catalogue True
thetaRange 0.01 110 0.01

To create the DSM information files to compute partials, run:

java io.github.kensuke1984.kibrary.Operation ifm.properties

This will create the threedPartial folder. The threedPartial folder contains 4 sub-folders with input files for computation of partial derivatives using DSM:

BPcat: for catalog of backward-propagated wavefield using shbpcat and psvbpcat
FPinfo: for the forward-propagated wavefield using shfp and psvfp
FPcat and BPinfo: not used in this tutorial

Move to the threedPartial/BPcat folder, and run:

mpirun -n 8 mpi-shbpcat < PREM_SH.inf
mpirun -n 8 mpi-psvbpcat < PREM_PSV.inf

Move to the threedPartial/FPinfo/031704A folder, and run:

mpirun -n 8 mpi-shfp < PREM_SH.inf
mpirun -n 8 mpi-psvfp < PREM_PSV.inf

Partial dataset

In the threedPartial folder, run:

java io.github.kensuke1984.kibrary.Property

When prompted to, enter 22 (AtAMaker). Rename the output file to atam.properties.

Create the file verticalMapping_4layers.inf by running

awk '{print $1,$2,int((NR-1)/2)}' ../perturbationLayer.inf > verticalMapping_4layers.inf

Edit atam.properties as follows:

nproc 8
components Z T
rootWaveformPath ../synthetic/filtered
waveformIDPath waveformID.dat
waveformPath waveform.dat
correctionTypes S
mode BOTH
thetaInfo 0.01 110. 0.01
sourceTimeFunction 2
timewindowPath ../synthetic/filtered/selectedTimewindow_ScS.dat
partialTypes MU LAMBDA
unknownParameterPath ../unknowns_1.25deg.inf
minFreq 0.005
maxFreq 0.08
filterNp 4
backward false
finalSamplingHz 1
resamplingRate 1
verticalMappingFile verticalMapping_4layers.inf

To compute the partial derivatives from the backward- and forward-propagated wavefield computed using DSM, run:

java io.github.kensuke1984.kibrary.Operation atam.properties

Rename the output files partial... and partialID... to partial_ScS.dat and partialID_ScS.dat. Rename the output file newPerturbationLayers... to newPerturbationLayers_4layers.inf, and newUnknowns... to newUnknowns_1.25deg_4layers.inf.

Do the same for the PcP phase, replacing

timewindowPath ../synthetic/filtered/selectedTimewindow_PcP.dat

in atam.properties.

Note: partials computed using AtAMaker are already multiplied by the volume of each voxel (in contrast to partials computed using PartialDatasetMaker.

9. Integrate the partials to a coarser mesh (2.5 by 2.5 degree)

Voxels of size 1.25 by 1.25 degrees result in too many model parameters to be handled in the inversion (i.e., to compute the AtA matrix). In order to reduce the number of model parameters, we integrate the mesh into a coarser mesh with voxels of size 2.5 by 2.5 degree. To integrate the partials for the ScS phase, in threedPartial, run:

java io.github.kensuke1984.kibrary.waveformdata.addons.CombineParameters partialID_ScS.dat partial_ScS.dat newUnknowns_1.25deg_4layers.inf ../horizontalMapping_sampled32.inf

Rename the output files to partialID_ScS_2.5deg.dat and partial_ScS_2.5deg.dat.

For the PcP phase, run:

java io.github.kensuke1984.kibrary.waveformdata.addons.CombineParameters partialID_PcP.dat partial_PcP.dat newUnknowns_1.25deg_4layers.inf ../horizontalMapping_sampled32.inf

Rename the output files to partialID_PcP_2.5deg.dat and partial_PcP_2.5deg.dat.

10. Convert partials from LAMBDA and MU to (LAMBDA+2MU) and MU

We change the parametrization from LAMBDA and MU to the P- and S-wave moduli M=LAMBDA+2MU and G=MU, since the former parametrization introduce an apparent trade-off between S- and P-velocities. The elastic tensor for an isotropic medium is written equivalently as:

$C_{ijpq}=\lambda\delta_{ij}\delta_{pq}+\mu\left(\delta_{ip}\delta_{jq}+\delta_{iq}\delta_{jp}\right)$
$\qquad =M\delta_{ij}\delta_{pq}+G\left(\delta_{ip}\delta_{jq}+\delta_{iq}\delta_{jp}-2\delta_{ij}\delta_{pq}\right)$

That is, the kernels for LAMBDA and M are identical, but the LAMBDA kernel is subtracted from the MU kernel to obtain the G kernel. This makes the M kernel sensitive only to P-velocity, and the G kernel only sensitive to the S-velocity, reducing trade-off between model parameters.

To compute the kernels for M and G for the ScS phase for the 2.5 by 2.5 degree mesh, in threedPartial, run

java io.github.kensuke1984.kibrary.waveformdata.convert.LambdaMuToLambda2MuMu partialID_ScS_2.5deg.dat partial_ScS_2.5deg.dat

Rename the output files to partialID_ScS_2.5deg_MG.dat and partial_ScS_2.5deg_MG.dat.

For the PcP phase, run

java io.github.kensuke1984.kibrary.waveformdata.convert.LambdaMuToLambda2MuMu partialID_PcP_2.5deg.dat partial_PcP_2.5deg.dat

Rename the output files to partialID_PcP_2.5deg_MG.dat and partial_PcP_2.5deg_MG.dat.

The new partial types are named LAMBDA2MU and MU in Kibrary. Create a new newUnknowns_1.25deg_4layers.inf file by running:

cat newUnknowns_1.25deg_4layers.inf | sed s/LAMBDA/LAMBDA2MU/ > newUnknowns_1.25deg_4layers.inf

11. Visualizing the kernels (Optional)

To output ASCII files for the kernel for the PcP phase, in threedPartial, run

java io.github.kensuke1984.kibrary.quick.KernelVisual partialID_PcP_2.5deg_MG.dat partial_PcP_2.5deg_MG.dat

This will create the folder KernelTemporalVisual with ASCII text files containing kernel values. A python script kernel.py to make an .mpg movie from the kernel ASCII files can be found in the tutorial folder. Movies for PcP and ScS kernels can be found here.

12. Combine data and partials for ScS and PcP

In threedPartial, run:

java io.github.kensuke1984.kibrary.waveformdata.addons.PartialIDMerge partialID_ScS_2.5deg_MG.dat partial_ScS_2.5deg_MG.dat partialID_PcP_2.5deg_MG.dat partial_PcP_2.5deg_MG.dat

Rename the output files to partialID_ScS_PcP_2.5deg_MG.dat and partial_ScS_PcP_2.5deg_MG.dat.

In tutorial_inversion/synthetic/filtered, run:

java io.github.kensuke1984.kibrary.waveformdata.addons.BasicIDMerge waveformID_ScS.dat waveform_ScS.dat waveformID_PcP.dat waveform_PcP.dat

Rename the output files to waveformID_ScS_PcP.dat and waveform_ScS_PcP.dat.

13. Inversion

Create the folder inversion (in threedPartial). In folder inversion, run:

java io.github.kensuke1984.kibrary.Property

When prompted to, enter 8 (LetMeInvert). Rename the output file to lmi.properties.

Edit lmi.properties as follows:

waveformIDPath ../../synthetic/filtered/waveformID_ScS_PcP.dat
waveformPath ../../synthetic/filtered/waveform_ScS_PcP.dat
partialIDPath ../partialID_ScS_PcP_2.5deg_MG.dat
partialPath ../partial_ScS_PcP_2.5deg_MG.dat
unknownParameterListPath ../newUnknowns_1.25deg_4layers.inf
stationInformationPath ../../station.inf
inverseMethods CG
modelCovariance true
cm0 1.
cmH 1.5
cmV 0.

modelCovariance true adds an a priori smoothness constraint to the model parameters using Gaussian smoothing. cm0 controls the strength of the constraint (larger leads to a more constrained model), and cmH and cmV give the horizontal, and vertical correlation lengths (in degrees). Note that we set cmV 0 (i.e., no vertical smoothness constraints) because we use 50 km thick layers, which are thick enough to be completely resolved. You will have to play with values for cm0 and cmH to find suitable regularization parameters for your dataset. For details on the definition of cm0 and cmH, see eqs. (4) and (5) of Borgeaud et al. (2017).

Launch the inversion by running:

java io.github.kensuke1984.kibrary.Operation lmi.properties

Rename the output folder to lmi.

14. Visualizing the results of the inversion

In lmi, run:

cp ../../newPerturbationLayers_4layers.inf perturbationLayer.inf
java io.github.kensuke1984.kibrary.inversion.addons.VelocityField3D 2.5

This will creates ASCII text files CG/velocityCGi.txt, where i is the number of conjugate gradient (CG) vectors used. Each line gives for a given voxel: latitude longitude depth dVs dVp dVb where dVb is the perturbation in bulk velocity.

To visualize the velocity file using i CG vectors, you can use the python script inversion.py in the tutorial/SCRIPTS folder, running:

python inversion.py i

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Tutorial for waveform inversion using Kibrary

1. Prerequisite

Optional - for visualization

2. Pre-tests (optional but recommended)

SHBPCAT

PSVBPCAT

3. Computation of synthetics

4. Time windows

5. Static corrections

6. Mantle corrections

7. Residual vector

8. Partial derivatives

Voxel mesh

station.inf

DSM information files

Partial dataset

9. Integrate the partials to a coarser mesh (2.5 by 2.5 degree)

10. Convert partials from LAMBDA and MU to (LAMBDA+2MU) and MU

11. Visualizing the kernels (Optional)

12. Combine data and partials for ScS and PcP

13. Inversion

14. Visualizing the results of the inversion

Clone this wiki locally