-
Notifications
You must be signed in to change notification settings - Fork 234
08_doing_tomography
Table of Contents
- Doing tomography, i.e., updating the model based on the sensitivity kernels obtained
See also the next chapter (Chapter [cha:FWI]) about how to perform full waveform inversion (FWI) or source inversions.
STILL UNDER CONSTRUCTION ... New content will be added soon.
One of the fundamental reasons for computing sensitivity kernels (Section [sec:Adjoint-simulation-finite]) is to use them within a tomographic inversion, i.e., perform imaging. In other words, use recorded seismograms, make measurements with synthetic seismograms, and use the misfit between the two to iteratively improve the model described by (at least)
As explained for instance in (Monteiller et al. 2015), full waveform inversion means that one considers the observed seismograms (possibly filtered) as the basic observables that one wants to fit. One thus searches for the model that minimizes the mean squared difference between observed and synthetic seismograms. In other words, the goal is to find a structural model that can explain a larger portion of seismological records, and not simply the phase of a few seismic arrivals.
Whatever misfit function you use for the tomographic inversion (you can find several examples in Tromp, Komatitsch, and Liu (2008) and Tromp, Tape, and Liu (2005)), you will weigh the sensitivity kernels with measurements. Furthermore, you will use as many measurements (stations, components, time windows) as possible per event; hence, we call these composite kernels “event kernels”, which are volumetric fields representing the gradient of the misfit function with respect to one of the variables (
We want to minimize the classical waveform misfit function:
A direct method to compute the gradient is to take the derivative of ([misfit]) with respect to model parameters:
Let us give an outline of the theory to compute the gradient based on the adjoint method, and refer the reader to e.g. (Tromp, Tape, and Liu 2005) and (Tromp, Komatitsch, and Liu 2008) for further details. The perturbation of the misfit function can be expressed as:
As can be seen from these expressions, the principle of the adjoint-state method is to correlate two wavefields: the direct (i.e. forward) field that propagates from the source to the receivers, and the adjoint field that propagates from all the receivers backward in time. The same approach can be followed for any type of seismic observable (phase, amplitude, envelope, time series...), provided the appropriate adjoint source is used (Tromp, Tape, and Liu 2005; Tromp, Komatitsch, and Liu 2008). For example, for the cross-correlated traveltime of a seismic phase, the adjoint source is defined as the velocity of that synthetic phase weighted by the travel-time residual.
Computing the gradient based on the adjoint-state method requires performing two simulations per source (forward and adjoint fields) regardless of the type of observable. However, to define the adjoint field one must know the adjoint source, and that source is computed from the results of the forward simulation. One must therefore perform the forward simulation before the adjoint simulation. A straightforward solution for time-domain methods would be to store the whole forward field to disk at each time step during the forward run and then read it back during the adjoint simulation to calculate the interaction of these two fields. In 2-D this is feasible but in the 3-D case for very short seismic periods and without lossy compression, downsampling, or a large amount of disk or memory checkpointing (e.g., Fichtner et al. 2009; Rubio Dalmau et al. 2014; Cyr, Shadid, and Wildey 2015) the amount of disk storage required would currently be too large.
However let us note again that this situation will change in the future. In the mean time, a standard possible solution is to perform three simulations per source (Tromp, Komatitsch, and Liu 2008; Peter et al. 2011), i.e., perform the forward calculation twice, once to compute the adjoint sources and once again at the same time as the adjoint simulation to correlate the two fields and sum their interaction on the fly over all the time steps. Doing so for an elastic Earth, one only needs a small amount of disk storage to store the last time step of the forward run, which is then used as an initial condition to redo the forward run backwards, as well as the field on the outer edges of the mesh for each time steps in order to be able to undo the absorbing boundary conditions.
Besides the ability to compute adjoint sensitivity kernels, the SPECFEM3D package provides a number of tools useful for post-processing kernels, gradients and conducting tomographic model updates. The tomographic and post-processing tools can be compiled by typing in the main directory:
make tomography
make postprocess
You can modify specific settings affecting, e.g., the (isotropic) parameterization or directory setup structure in file setup/constants_tomography.h
. The compiled binaries can be found in directory bin/
.
The iterative inversion workflow may comprise three main steps for which tools are provided:
-
Summing event kernels, e.g.
$K_{\alpha,\beta,\rho}$ , to build the misfit gradient$\boldg(\boldm)$ -
Smoothing and post-processing of the gradient
$\boldg(\boldm)$ -
Updating the model
$\boldm_{i+1} = \boldm_i + \delta\boldm$ .
Please note that the tomographic routines are very similar for both SPECFEM3D_Cartesian and SPECFEM3D_GLOBE versions. Some differences occur, e.g., when reading in mesh files and when the maximum of the gradient is taken for the update step length. Also note that currently only the binary database format is supported, no ADIOS file format support is implemented yet. In future, we plan to extend and merge these into the same set of tools within a common SPECFEM function library.
For summing different event kernels, the executables xsum_kernels
and xsum_preconditioned_kernels
can be used. The binaries sum up either a set of isotropic or transversely isotropic kernel files, depending on the setting in file setup/constants_tomography.h
. The following parameterizations can be chosen:
-
isotropic kernels
$K_{\alpha,\beta,\rho}$ or -
isotropic bulk kernels
$K_{c,\beta,\rho}$ , or -
transversely isotropic bulk kernels
$K_{c,\beta_v,\beta_h,\eta}$ .
The tools for summing event kernels can be used in parallel, using the same number of processes (NPROC
) as the forward/kernel simulations. The summation tools use the following input/output format (see Fig. 1.1 for the default directory structure):
-
Input: Input is provided by a list of event directories, given in file
kernels_list.txt
(default setting forKERNEL_FILE_LIST
inconstants_tomography.h
). The file lists all event kernel directories which should be summed together. All event directories have to be located in directoryINPUT_KERNELS/
, which can be setup as links to all individual event kernel directories. -
Output: The summed kernels will be stored in corresponding kernel files in output directory
OUTPUT_SUM/
.
Kernel summation:
As example with a total of 4 MPI processes
mpirun -np 4 ./bin/xsum_kernels
adds sensitivity kernels kernel_list.txt
.
Kernel Summation with Preconditioning:
As example with a total of 4 MPI processes
mpirun -np 4 ./bin/xsum_preconditioned_kernels
adds sensitivity kernels
Kernel names: The summation tool will look for kernels with following names:
-
for an isotropic model parameterization, kernel names are:
[TABLE]
-
for an isotropic bulk model parameterization, kernel names are:
[TABLE]
-
for a transversely isotropic model parameterization, kernel names are:
[TABLE]
Note that these event kernels are stored after a kernel simulation (
SIMULATION_TYPE = 3
) in theLOCAL_PATH
directory, which by default points to directoryDATABASES_MPI/
. Isotropic kernels will be created by default. To create transversely isotropic kernels, you setANISOTROPIC_KL
andSAVE_TRANSVERSE_KL
to.true.
in the parameter fileDATA/Par_file
. -
for preconditioning, the approximate Hessian kernels taken for summation are:
[TABLE]
These Hessian kernels approximate the diagonal elements of the Hessian matrix and can be used as preconditioner in a gradient optimization scheme. To create these kernels, you have to set
APPROXIMATE_HESS_KL
to.true.
in the parameter fileDATA/Par_file
.
For SPECFEM3D_GLOBE, by default only kernels in the crust/mantle region ("reg1") will be considered for summation. You can change the default by setting parameter REG
in file constants_tomography.h
to the region of interest.
Note that although we provide the preconditioned summation here, we recommend to smooth first both, the summed kernels and summed approximate Hessians, before inverting the summed Hessian and applying it as preconditioner to the gradient. From experience, we prefer doing $$\boldg = \frac{1}{\mathcal{F}{smooth} \left( \Sum^{N}i \tilde{H}^{i} \right) } ; \mathcal{F}{smooth} \left( \Sum_i K^{i} \right) \nonumber$$ rather than $$\boldg = \mathcal{F}{smooth} \left( \frac{1}{\Sum^{N}_i \tilde{H}^{i}} ; \Sum_i K^{i} \right) : . \nonumber$$ Due to the large sensitivity kernel values close to source and receiver locations, the inverse of the thresholded summed Hessian becomes better balanced in the former.
For additional smoothing and post-processing of kernels, gradients, or models, we provide tools useful for gradient and model preparation. These additional tools for post-processing sensitivity kernels are provided together with the main package. For compilation of these tools, you type in the main directory:
make postprocess
The following tools are provided:
-
xsmooth_sem
is used for smoothing with a gaussian function. -
xclip_sem
is used to threshold a kernel, gradient or model, clipping to a minimum and maximum value. -
xcombine_sem
is used to combine different event kernels, gradient or model files, summing all individually specified files together (similar toxsum_kernels
but more flexible).
The post-processing tools are primarily intented to be used to process kernel files. They can be used though on any scalar field of dimension (NGLLX,NGLLY,NGLLZ,NSPEC)
. The tools are parallel programs – they must be invoked with mpirun or other appropriate utility. Operations are performed in embarrassingly-parallel fashion.
Smoothing:
To smooth a kernel, gradient or model file, use
mpirun -np NPROC bin/xsmooth_sem SIGMA_H SIGMA_V KERNEL_NAME INPUT_DIR OUPUT_DIR
with command line arguments: SIGMA_H
horizontal smoothing radius, SIGMA_V
vertical smoothing radius, KERNEL_NAME
kernel name, e.g. alpha_kernel
, INPUT_DIR
directory from which kernels are read, OUTPUT_DIR
directory to which smoothed kernels are written.
This smooths kernels by convolution with a Gaussian, and writes the resulting smoothed kernels to OUTPUT_DIR
.
Files written to OUTPUT_DIR
have the suffix ’smooth’ appended, e.g. proc***alpha_kernel.bin
becomes proc***alpha_kernel_smooth.bin
.
Clipping:
Values in a kernel, gradient or model binary file can be clipped with a minimum/maximum threshold value using
mpirun -np NPROC bin/xclip_sem MIN_VAL MAX_VAL KERNEL_NAMES INPUT_FILE OUTPUT_DIR
with command line arguments: MIN_VAL
threshold below which array values are clipped, MAX_VAL
threshold above which array values are clipped, KERNEL_NAMES
one or more kernel names separated by commas, INPUT_DIR
directory from which arrays are read, OUTPUT_DIR
directory to which clipped array are written.
For each name in KERNEL_NAMES
, reads kernels from INPUT_DIR
, applies thresholds, and writes the resulting clipped kernels to OUTPUT_DIR
.
KERNEL_NAMES
is a comma-delimited list of kernel names, e.g. alphav_kernel,alphah_kernel
.
-
For example, in this case type
mpirun -np 4 ./bin/xclip_sem -0.1 0.1 alphav_kernel,alphah_kernel DIR1/ DIR1/
to clip the corresponding kernels in directory
DIR1
to be in range$[-0.1,0.1]$ .
Files written to OUTPUT_DIR
have the suffix ’clip’ appended, e.g. proc***alpha_kernel.bin
becomes proc***alpha_kernel_clip.bin
Combining:
To combine kernel values from different kernel directories, you can use
mpirun -np NPROC bin/xcombine_sem KERNEL_NAMES INPUT_FILE OUTPUT_DIR
with command line arguments: KERNEL_NAMES
one or more kernel names separated by commas, INPUT_FILE
text file containing list of kernel directories, OUTPUT_PATH
directory to which summed kernels are written.
For each name in KERNEL_NAMES
, sums kernels from directories specified in INPUT_FILE
. Writes the resulting sums to OUTPUT_DIR
.
INPUT_FILE
is a text file containing a list of absolute or relative paths to kernel directories, one directory per line.
KERNEL_NAMES
is a comma-delimited list of kernel names, e.g. alpha_kernel,beta_kernel,rho_kernel
.
We provide (simple) tools for updating the current model
-
xadd_model_iso
can be used to update isotropic model files with a (summed & smoothed) gradient. The gradient files are given for isotropic parameters or isotropic bulk parameters.
Note that instead of using the density kernels USE_RHO_SCALING
to .true.
will ignore the density kernel and use density perturbations scaled from isotropic
Isotropic model update:
The program xadd_model_iso
can be used to update isotropic model files with (smoothed & summed) event kernels:
mpirun -np 4 bin/xadd_model_iso step_factor [INPUT-KERNELS-DIR/] [OUTPUT-MODEL-DIR/]
with command line arguments: step_factor
the step length to scale the gradient, e.g. INPUT-KERNELS-DIR/
(optional) directory which holds summed kernels (e.g. proc***alpha_kernel.bin
,.., by default directory INPUT_GRADIENT/
is used), OUTPUT-MODEL-DIR/
(optional) directory which will hold new model files (e.g. proc***vp_new.bin
,.., by default directory OUTPUT_MODEL/
is used).
The gradients are given for isotropic parameters (alpha,beta,rho)
or (bulk_c,beta,rho)
.
The algorithm uses a steepest descent method with a step length determined by the given maximum update percentage.
By default, the directory and file setup used is:
-
directory
INPUT_MODEL/
contains:proc***_vs.bin proc***_vp.bin proc***_rho.bin
-
directory
INPUT_GRADIENT/
contains:proc***_bulk_c_kernel_smooth.bin proc***_bulk_beta_kernel_smooth.bin proc***_rho_kernel_smooth.bin
or
proc***_alpha_kernel_smooth.bin proc***_beta_kernel_smooth.bin proc***_rho_kernel_smooth.bin
depending on the model parameterization.
-
directory
topo/
contains:proc***_external_mesh.bin
for the Cartesian version, and
proc***_solver_data.bin
for the GLOBE version.
The new model files are stored in directory OUTPUT_MODEL/
as
proc***_vp_new.bin
proc***_vs_new.bin
proc***_rho_new.bin
Note that additional post-processing and model update tools are provided for the SPECFEM3D_GLOBE version and will be provided for the Cartesian version as well in future.
There are many other versions of gradient-based inversion algorithms that could alternatively be used (see e.g. (Virieux and Operto 2009; Monteiller et al. 2015) for a list). The tomographic inversion of Tape et al. (2009, 2010) used SPECFEM3D_Cartesian as well as several additional components which are also stored on the CIG svn server, described next.
The directory containing external utilities for tomographic inversion using SPECFEM3D Cartesian (or other packages that evaluate misfit functions and gradients) is in directory utils/ADJOINT_TOMOGRAPHY_TOOLS/
:
flexwin/ -- FLEXWIN algorithm for automated picking of time windows
measure_adj/ -- reads FLEXWIN output file and makes measurements,
with the option for computing adjoint sources
iterate_adj/ -- various tools for iterative inversion
(requires pre-computed "event kernels")
This directory also contains a brief README
file indicating the role of the three subdirectories, flexwin
(Maggi et al. 2009), measure_adj
, and iterate_adj
. The components for making the model update are there; however, there are no explicit rules for computing the model update, just as with any optimization problem. There are options for computing a conjugate gradient step, as well as a source subspace projection step.
The best single file to read is probably: ADJOINT_TOMO/iterate_adj/cluster/README
.
Aravkin, Aleksandr, Tristan van Leeuwen, and Felix Herrmann. 2011. “Robust Full-Waveform Inversion Using the Student’s
Bamberger, A., G. Chavent, Ch. Hemons, and P. Lailly. 1982. “Inversion of Normal Incidence Seismograms.” Geophysics 47 (5): 757–70.
Bamberger, A., G. Chavent, and P. Lailly. 1977. “Une Application de La Théorie Du Contrôle à Un Problème Inverse de Sismique.” Annales de Géophysique 33 (1/2): 183–200.
Brossier, R., S. Operto, and J. Virieux. 2010. “Which Data Residual Norm for Robust Elastic Frequency-Domain Full Waveform Inversion?” Geophysics 75 (3): R37–46.
Bube, Kenneth P., and Robert T. Langan. 1997. “Hybrid
Chavent, G. 1974. “Identification of Function Parameters in Partial Differential Equations.” In Identification of Parameter Distributed Systems, edited by R. E. Goodson and M. Polis, 31–48. American Society of Mechanical Engineers.
Crase, E., A. Pica, M. Noble, J. McDonald, and A. Tarantola. 1990. “Robust Elastic Non-Linear Waveform Inversion: Application to Real Data.” Geophysics 55: 527–38.
Cyr, E. C., J. N. Shadid, and T. Wildey. 2015. “Towards Efficient Backward-in-Time Adjoint Computations Using Data Compression Techniques.” Comput. Methods Appl. Mech. Eng. 288: 24–44. https://doi.org/10.1016/j.cma.2014.12.001.
Fichtner, A., H.-P. Bunge, and H. Igel. 2006. “The Adjoint Method in Seismology: I. Theory.” Phys. Earth Planet. Inter. 157 (1-2): 86–104. https://doi.org/10.1016/j.pepi.2006.03.016.
Fichtner, A., B. L. N. Kennett, H. Igel, and H. P. Bunge. 2009. “Full Seismic Waveform Tomography for Upper-Mantle Structure in the Australasian Region Using Adjoint Methods.” Geophys. J. Int. 179 (3): 1703–25.
Ha, Taeyoung, Wookeen Chung, and Changsoo Shin. 2009. “Waveform Inversion Using a Back-Propagation Algorithm and a Huber Function Norm.” Geophysics 74 (3): R15–24. https://doi.org/10.1190/1.3112572.
Hudson, J. A. 1977. “Scattered Waves in the Coda of
Jeong, Woodon, Minji Kang, Shinwoong Kim, Dong-Joo Min, and Won-Ki Kim. 2015. “Full Waveform Inversion Using Student’s
Lailly, P. 1983. “The Seismic Inverse Problem as a Sequence of Before-Stack Migrations.” In Proceedings of the Conference on Inverse Scattering, Theory and Application Expanded Abstracts, edited by J. B. Bednar, R. Redner, E. Robinson, and A. Weglein, 206–20. Philadelphia, PA, USA: Society of Industrial; Applied Mathematics.
Maggi, A., C. Tape, M. Chen, D. Chao, and J. Tromp. 2009. “An Automated Time-Window Selection Algorithm for Seismic Tomography.” Geophys. J. Int. 178: 257–81.
Monteiller, Vadim, Sébastien Chevrot, Dimitri Komatitsch, and Yi Wang. 2015. “Three-Dimensional Full Waveform Inversion of Short-Period Teleseismic Wavefields Based Upon the SEM-DSM Hybrid Method.” Geophys. J. Int. 202 (2): 811–27. https://doi.org/10.1093/gji/ggv189.
Peter, Daniel, Dimitri Komatitsch, Yang Luo, Roland Martin, Nicolas Le Goff, Emanuele Casarotti, Pieyre Le Loher, et al. 2011. “Forward and Adjoint Simulations of Seismic Wave Propagation on Fully Unstructured Hexahedral Meshes.” Geophys. J. Int. 186 (2): 721–39. https://doi.org/10.1111/j.1365-246X.2011.05044.x.
Plessix, R. E. 2006. “A Review of the Adjoint-State Method for Computing the Gradient of a Functional with Geophysical Applications.” Geophys. J. Int. 167 (2): 495–503.
Rubio Dalmau, F., M. Hanzich, J. de la Puente, and N. Gutiérrez. 2014. “Lossy Data Compression with DCT Transforms.” In Proceedings of the EAGE Workshop on High Performance Computing for Upstream, HPC30. Chania, Crete, Greece. https://doi.org/10.3997/2214-4609.20141939.
Tape, Carl, Qinya Liu, Alessia Maggi, and Jeroen Tromp. 2009. “Adjoint Tomography of the Southern California Crust.” Science 325: 988–92.
———. 2010. “Seismic Tomography of the Southern California Crust Based on Spectral-Element and Adjoint Methods.” Geophys. J. Int. 180: 433–62.
Tape, Carl, Qinya Liu, and Jeroen Tromp. 2007. “Finite-Frequency Tomography Using Adjoint Methods - Methodology and Examples Using Membrane Surface Waves.” Geophys. J. Int. 168 (3): 1105–29. https://doi.org/10.1111/j.1365-246X.2006.03191.x.
Tarantola, A. 1984. “Inversion of Seismic Reflection Data in the Acoustic Approximation.” Geophysics 49: 1259–66.
Tromp, Jeroen, Dimitri Komatitsch, and Qinya Liu. 2008. “Spectral-Element and Adjoint Methods in Seismology.” Communications in Computational Physics 3 (1): 1–32.
Tromp, Jeroen, Carl Tape, and Qinya Liu. 2005. “Seismic Tomography, Adjoint Methods, Time Reversal and Banana-Doughnut Kernels.” Geophys. J. Int. 160 (1): 195–216. https://doi.org/10.1111/j.1365-246X.2004.02453.x.
Virieux, J., and S. Operto. 2009. “An Overview of Full-Waveform Inversion in Exploration Geophysics.” Geophysics 74 (6): WCC1–26. https://doi.org/10.1190/1.3238367.
This documentation has been automatically generated by pandoc based on the User manual (LaTeX version) in folder doc/USER_MANUAL/ (Nov 15, 2023)
Development wiki for SPECFEM3D
Development wiki
- Home
- Development plan
- Best Practices
- Using Git for SPECFEM
- Advanced Git Topics
- Versioning Conventions
- Merging Development Branches into Master
User manual
- 01_introduction
- 02_getting_started
- 03_mesh_generation
- 04_creating_databases
- 05_running_the_solver
- 06_fault_sources
- 07_adjoint_simulations
- 08_doing_tomography
- 08b_performing_full_waveform_inversion_FWI_or_source_inversions
- 09_noise_simulations
- 10_gravity_calculations
- 11_graphics
- 12_running_scheduler
- 13_changing_the_model
- 14_post_processing
- 15_informations_for_developers
- A_reference_frame
- B_channel_codes
- C_troubleshooting
- D_license
- authors
- bug_reports
- copyright_and_version
- features
- manual_SPECFEM3D_Cartesian
- notes_and_acknowledgement
- sponsors