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
Download the executable JAR file of Kibrary, anselme branch
Dowload the executable JAR file of TomoTool
Download shbpcat, and psvbpcat
Dowload tish, tipsv, shfp, and psvfp

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

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), 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.

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

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 on which the partials will be computed, and then integrated to the coarser mesh. 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.

DSM information files

java io.github.kensuke1984.kibrary.Property

When prompted to, enter 7 (InformationFileMaker).

Edit the file io.github.kensuke1984.kibrary.dsminformation.InformationFileMaker.properties as follows

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

station.inf

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

threedPartial folder

Run

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

The threedPartial folder contains 4 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

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

cd FPinfo/031704A
mpirun -n 8 mpi-shfp < PREM_SH.inf
mpirun -n 8 mpi-psvfp < PREM_PSV.inf

Partial dataset

java io.github.kensuke1984.kibrary.Property

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

Create the file verticalMapping.inf by running

awk '{print $1,$2,NR-1}' ../perturbationLayer.inf > verticalMapping.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.inf

Run

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

Rename the output files partial... and partialID... to partial_ScS.dat and partialID_ScS.dat.

Do the same for the PcP phase, replacing

timewindowPath ../synthetic/filtered/selectedTimewindow_PcP.dat

in atam.properties. Here you can ignore the newUnknowns... file.

9. Visualizing the kernels (Optional)

To output ASCII files for the kernel for the ScS phase, run

java io.github.kensuke1984.kibrary.quick.KernelVisual partialID_ScS.dat partial_ScS.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.

10. 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, run:

java io.github.kensuke1984.kibrary.waveformdata.addons.CombineParameters partialID_ScS.dat partial_ScS.dat ../unknowns_1.25deg.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 ../unknowns_1.25deg.inf ../horizontalMapping_sampled32.inf

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

11. 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, because 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, 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 unknowns_2.5deg_MG.inf file by running:

cat ../unknowns_2.5deg.inf | sed s/LAMBDA/LAMBDA2MU/ > unknowns_2.5deg_MG.inf

12. 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.dat
waveformPath ../../synthetic/filtered/waveform_ScS.dat
partialIDPath ../partialID_ScS_2.5deg_MG.dat
partialPath ../partial_ScS_2.5deg_MG.dat
unknownParameterListPath ../unknowns_2.5deg_MG.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_ScS.

Do the same for the PcP phase, replacing:

waveformIDPath ../../synthetic/filtered/waveformID_PcP.dat
waveformPath ../../synthetic/filtered/waveform_PcP.dat
partialIDPath ../partialID_PcP_2.5deg_MG.dat
partialPath ../partial_PcP_2.5deg_MG.dat

in lmi.properties.

Launch the inversion by running:

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

Rename the output folder to lmi_PcP.

13. Visualizing the results of the inversion

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

DSM information files

station.inf

threedPartial folder

Partial dataset

9. Visualizing the kernels (Optional)

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

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

12. Inversion

13. Visualizing the results of the inversion

Clone this wiki locally