seismic.inversion.wavefield_decomp package¶
Submodules¶
seismic.inversion.wavefield_decomp.call_count_decorator module¶
Decorator to count number of times a function is called.
seismic.inversion.wavefield_decomp.plot_nd_batch module¶
Batch plotting a MCMC solution for batch of stations to a single pdf file.

seismic.inversion.wavefield_decomp.plot_nd_batch.
convert_Vs_to_k
(soln, config)[source]¶ Transform Vs variable into k variable in MCMC solution. Modifies soln inplace.
 Parameters
soln (Customized scipy.optimize.OptimizeResult) – Solution container
config (dict) – Solution configuration
 Returns
None

seismic.inversion.wavefield_decomp.plot_nd_batch.
plot_aux_data
(soln, config, log, scale)[source]¶ Plot auxiliary data such as energy distribution and receiver functions.
 Parameters
soln (Customized scipy.optimize.OptimizeResult) – Solution container
config (dict) – Solution configuration
log (logging.Logger) – Logging instance
scale (float) – Overall image scaling factor
 Returns
Matplotlib figure containing the plotted data
 Return type
seismic.inversion.wavefield_decomp.runners module¶
Batch execution interfaces for wavefield continuation methods and solvers.

seismic.inversion.wavefield_decomp.runners.
curate_seismograms
(data_all, curation_opts, logger, rotate_to_zrt=True)[source]¶ Curation function to remove bad data from streams. Note that this function will modify the input dataset during curation.
 Parameters
data_all (seismic.network_event_dataset.NetworkEventDataset) – NetworkEventDataset containing seismograms to curate. Data will be modified by this function.
curation_opts (dict) – Dict containing curation options.
logger (logging.Logger) – Logger for emitting log messages
rotate_to_zrt (bool) – Whether to automatically rotate to ZRT coords.
 Returns
None, curation operates directly on data_all

seismic.inversion.wavefield_decomp.runners.
load_mcmc_solution
(h5_file, job_timestamp=None, logger=None)[source]¶ Load Monte Carlo Markov Chain solution from HDF5 file.
 Parameters
h5_file (str or pathlib.Path) – File from which to load solution
job_timestamp (str or NoneType) – Timestamp of job whose solution is to be loaded
logger (logging.Logger) – Output logging instance
 Returns
(solution, job configuration), job timestamp
 Return type

seismic.inversion.wavefield_decomp.runners.
mcmc_solver_wrapper
(model, obj_fn, mantle, Vp, rho, flux_window)[source]¶ Wrapper callable for passing to MCMC solver in scipy style, which unpacks inputs into vector variables for solver.
 Parameters
model (numpy.array) – Perlayer model values (the vector being solved for) as flat array of (H, Vs) value pairs ordered by layer.
obj_fn – Callable to WfContinuationSuFluxComputer to compute SU flux.
mantle (seismic.model_properties.LayerProps) – Mantle properties stored in class LayerProps instance.
Vp (numpy.array) – Array of Vp values ordered by layer.
rho (numpy.array) – Array of rho values ordered by layer.
flux_window ((float, float)) – Pair of floats indicating the time window over which to perform SU flux integration
 Returns
Integrated SU flux energy at top of mantle

seismic.inversion.wavefield_decomp.runners.
run_mcmc
(waveform_data, config, logger)[source]¶ Top level runner function for MCMC solver on SU flux minimization for given settings.
 Parameters
waveform_data (Iterable obspy.Stream objects) – Collection of waveform data to process
config (dict) – Dict of job settings. See example files for fields and format of settings.
logger (logging.Logger or NoneType) – Log message receiver. Optional, pass None for no logging.
 Returns
Solution object based on scipy.optimize.OptimizeResult
 Return type
scipy.optimize.OptimizeResult with additional custom attributes

seismic.inversion.wavefield_decomp.runners.
run_station
(config_file, waveform_file, network, station, location, logger)[source]¶ Runner for analysis of single station. For multiple stations, set up config file to run batch job using mpi_job CLI.
The output file is in HDF5 format. The configuration details are added to the output file for traceability.
 Parameters
config_file (dict) – Config filename specifying job settings
waveform_file (str or pathlib.Path) – Event waveform source file for seismograms, generated using extract_event_traces.py script
network (str) – Network code of station to analyse
station (str) – Station code to analyse
location (str) – Location code of station to analyse. Can be ‘’ (empty string) if not set.
logger (logging.Logger) – Output logging instance
 Returns
Pair containing (solution, configuration) containers. Configuration will have additional traceability information.
 Return type
(solution, dict)

seismic.inversion.wavefield_decomp.runners.
save_mcmc_solution
(soln_configs, input_file, output_file, job_timestamp, job_tracking, logger=None)[source]¶ Save solution to HDF5 file. In general soln_configs will be an ensemble of perstation solutions.
 Parameters
soln_configs (list((solution, dict))) – List of (solution, configuration) pairs to save (one per station).
input_file (str) – Name of input file used for the job. Saved to job node for traceability
output_file (str or pathlib.Path) – Name of output file to create
job_timestamp (str) – Job timestamp that will be used to generate the top level job group
job_tracking (dict or dictlike) – Dict containing job identification information for traceability
logger (logging.Logger) – [OPTIONAL] Log message destination
 Returns
None
seismic.inversion.wavefield_decomp.solvers module¶
Objective function minimization solvers.

class
seismic.inversion.wavefield_decomp.solvers.
AdaptiveStepsize
(stepper, accept_rate=0.5, interval=50, factor=0.9, ar_tolerance=0.05)[source]¶ Bases:
object
Class to implement adaptive stepsize. Pulled from scipy.optimize.basinhopping and adapted.

class
seismic.inversion.wavefield_decomp.solvers.
BoundedRandNStepper
(bounds, initial_step=None)[source]¶ Bases:
object
Step one dimensional at a time using normal distribution of steps within defined parameter bounds.

property
stepsize
¶

property

class
seismic.inversion.wavefield_decomp.solvers.
HistogramIncremental
(bounds, nbins=20)[source]¶ Bases:
object
Class to incrementally accumulate Ndimensional histogram stats at runtime.

property
bins
¶

property
dims
¶

property
histograms
¶

property

class
seismic.inversion.wavefield_decomp.solvers.
SolverGlobalMhMcmc
[source]¶ Bases:
object
Dropin custom solver for scipy.optimize.minimize, based on MetrolpolisHastings Monte Carlo Markov Chain random walk with burnin and adaptive acceptance rate, followed by Ndimensional clustering.
Rather than returning one global solution, the solver returns the N best ranked local solutions. It also returns a probability distribution of each unknown based on Monte Carlo statistics.

seismic.inversion.wavefield_decomp.solvers.
optimize_minimize_mhmcmc_cluster
(objective, bounds, args=(), x0=None, T=1, N=3, burnin=100000, maxiter=1000000, target_ar=0.4, ar_tolerance=0.05, cluster_eps=0.05, rnd_seed=None, collect_samples=None, logger=None)[source]¶ Minimize objective function and return up to N local minima solutions.
 Parameters
objective (Callable(*args) > float) – Objective function to minimize. Takes unpacked args as function call arguments and returns a float.
bounds (scipy.optimize.Bounds) – Bounds of the parameter space.
args (tuple or list) – Any additional fixed parameters needed to completely specify the objective function.
x0 (numpy.array with same shape as elements of bounds) – Initial guess. If None, will be selected randomly and uniformly within the parameter bounds.
T (float) – The “temperature” parameter for the accept or reject criterion. To sample the domain well, should be in the order of the typical difference in local minima objective valuations.
N (int) – Maximum number of minima to return
burnin (int) – Number of random steps to discard before starting to accumulate statistics.
maxiter (int) – Maximum number of steps to take (including burnin).
target_ar (float between 0 and 1) – Target acceptance rate of point samples generated by stepping.
ar_tolerance (float) – Tolerance on the acceptance rate before actively adapting the step size.
cluster_eps (float) – Point proximity tolerance for DBSCAN clustering, in normalized bounds coordinates.
rnd_seed (int) – Random seed to force deterministic behaviour
collect_samples (int or NoneType) – If not None and integral type, collect collect_samples at regular intervals and return as part of solution.
logger – Logger instance for outputting log messages.
 Returns
OptimizeResult containing solution(s) and solver data.
 Return type
scipy.optimize.OptimizeResult with additional attributes
seismic.inversion.wavefield_decomp.wavefield_continuation_tao module¶
Class encapsulating algorithm for 1D inversion using wavefield continuation.
Based on reference:
Kai Tao, Tianze Liu, Jieyuan Ning, Fenglin Niu, “Estimating sedimentary and crustal structure using wavefield continuation: theory, techniques and applications”, Geophysical Journal International, Volume 197, Issue 1, April, 2014, Pages 443457, https://doi.org/10.1093/gji/ggt515

class
seismic.inversion.wavefield_decomp.wavefield_continuation_tao.
WfContinuationSuFluxComputer
(station_event_dataset, f_s, time_window, cut_window)[source]¶ Bases:
object
Implements computation of the upwards mean Swave energy flux at the top of the mantle for an ensemble of events for one station.
Class instance can be loaded with a dataset and then evaluated for arbitrary 1D earth models.
Process:
Load data from a station event dataset. Copy of the data is buffered by this class in efficient format for energy flux calculation.
Define 1D earth model and mantle halfspace material properties (in external code).
Call instance with models and receive energy flux results.
Constructor
 Parameters
station_event_dataset – Iterable container of obspy.Stream objects.
f_s – Processing sample rate. Usuually much less than the sampling rate of the input raw seismic traces.
time_window – Time window about onset to use for wave continuation processing.
cut_window – Shorter time segment within time_window to from which to extract primary arrival waveform and its multiples.

grid_search
(mantle_props, layer_props, layer_index, H_vals, k_vals, flux_window= 10, 20, ncpus= 1)[source]¶ Compute SU energy flux over a grid of H and k values for a particular (single) layer. The layer to compute over is indicated by layer_index, which is a zerobased index into layer_props.
 Parameters
mantle_props (LayerProps) – Mantle bulk properties
layer_props (list(seismic.model_properties.LayerProps)) – Layer bulk properties as a list
layer_index (int) – Index of layer to vary.
H_vals (numpy.array) – Set of thickness (H) values to apply to the selected layer.
k_vals (numpy.array) – Set of k values to apply to the selected layer.
flux_window (pair of numeric values) – Time window in which to compute energy flux
ncpus (int) – Number of CPUs to use. Use 1 to use all.
 Returns
tuple(H grid, k grid, energy values)
 Return type
tuple(numpy.array, numpy.array, numpy.array)

propagate_to_base
(layer_props)[source]¶ Given a stack of layers in layer_props, propagate the surface seismograms down to the base of the bottom layer.
 Parameters
layer_props (list(seismic.model_properties.LayerProps)) – List of LayerProps through which to propagate the surface seismograms.
 Returns
(Vr, Vz) velocity components time series per event at the bottom of the stack of layers.
 Return type
numpy.array
seismic.inversion.wavefield_decomp.wfd_plot module¶
Plotting helper functions for Wavefield Decomposition module.

seismic.inversion.wavefield_decomp.wfd_plot.
plot_Esu_space
(H, k, Esu, title=None, savefile_name=None, show=True, c_range=None, None, decorator=None)[source]¶ Plot SU energy as function of Hk.
 Parameters
H (numpy.array) – Grid of H coordinates corresponding to Esu point values.
k (numpy.array) – Grid of k coordinates corresponding to Esu point values.
Esu (numpy.array) – Energy values on Hk grid.
title (str) – Plot title [OPTIONAL]
savefile_name (str) – Output file in which to save plot [OPTIONAL]
show (bool) – If True, display the image and block until it is closed.
c_range (tuple(float, float)) – Custom range of Esu to contour (min, max values)
decorator (Callable) – Callback function to customize plot.
 Returns
None

seismic.inversion.wavefield_decomp.wfd_plot.
plot_Nd
(soln, title='', scale=1.0, vars=None)[source]¶ Plotting routine for Ndimensional solution in grid format. Diagonal contains histograms of solution samples for each variable, and offdiagonal contains pairwise scatter plots of sample points plus solution clusters. The distinct solutions coordinates are colour coded and labelled on the histograms.
This function is intended to be a domainagnostic means of plotting an Ndimensional solution generated by optimize_minimize_mhmcmc_cluster(). It should not have seismologyspecific fields or terminology added into it.
 Parameters
soln (scipy.optimize.OptimizeResult with additional custom fields) – Solution container returned from optimize_minimize_mhmcmc_cluster()
title (str, optional) – Overall plot title, defaults to ‘’ (no title)
scale (float, optional) – Scale size of overall plot, defaults to 1.0. Adjust this depending on dimensionality or desired size.
vars (list(str)) – List of names of variables, should be same length as dimensionality of solution
 Returns
Tuple containing the PairGrid, list of axes for the secondary (histogram) axes in the diagonal of the grid, and list of text items representing solution labels.
 Return type
tuple(seaborn.PairGrid, list(matplotlib.axes.Axes), list(matplotlib.text.Text))