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.call_count_decorator.call_counter(func)

Decorator to count calls to a function. The number of calls can be queryied from func.counter.

Parameters

func – Function whose calls to count

Returns

func wrapper

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)

Transform Vs variable into k variable in MCMC solution. Modifies soln in-place.

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)

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)

Curation function to remove bad data from streams. Note that this function will modify the input dataset during curation.

Parameters

• data_all – NetworkEventDataset containing seismograms to curate.

• curation_opts – Dict containing curation options.

• logger – Logger for emitting log messages

• rotate_to_zrt – 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)

Load Monte Carlo Markov Chain solution from HDF5 file.

Parameters

• h5_file – File from which to load solution

• job_timestamp – Timestamp of job whose solution is to be loaded

• logger – Output logging instance

Returns

(solution, job configuration), job timestamp

seismic.inversion.wavefield_decomp.runners.mcmc_solver_wrapper(model, obj_fn, mantle, Vp, rho, flux_window)

Wrapper callable for passing to MCMC solver in scipy style, which unpacks inputs into vector variables for solver.

Parameters

• model – Per-layer 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 – Mantle properties stored in class LayerProps instance.

• Vp – Array of Vp values ordered by layer.

• rho – Arroy of rho values ordered by layer.

• flux_window – 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)

Top level runner function for MCMC solver on SU flux minimization for given settings.

Parameters

• waveform_data – Iterable container of obspy.Stream objects.

• config – Dict of job settings. See example files for fields and format of settings.

• logger – Log message receiver. Optional, pass None for no logging.

Returns

Solution object based on scipy.optimize.OptimizeResult

seismic.inversion.wavefield_decomp.runners.run_station(config_file, waveform_file, network, station, location, logger)

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 – Config filename specifying job settings

• waveform_file – Event waveform source file for seismograms, generated using extract_event_traces.py script

• network – Network code of station to analyse

• station – Station code to analyse

• location – Location code of station to analyse. Can be ‘’ (empty string) if not set.

• logger – Output logging instance

Returns

Pair containing (solution, configuration) containers. Configuration will have additional traceability information.

seismic.inversion.wavefield_decomp.runners.save_mcmc_solution(soln_configs, input_file, output_file, job_timestamp, job_tracking, logger=None)

Save solution to HDF5 file.

Parameters

• soln_configs – List of (solution, configuration) pairs to save (one per station).

• input_file – Name of input file. Saved to job node for traceability

• output_file – Name of output file to create

• job_timestamp – Job timestamp that will be used to generate the top level job group

• job_tracking – Dict containing job identification information for traceability

• 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)

Bases: object

Class to implement adaptive stepsize. Pulled from scipy.optimize.basinhopping and adapted.

notify_accept()

take_step(x)

class seismic.inversion.wavefield_decomp.solvers.BoundedRandNStepper(bounds, initial_step=None)

Bases: object

Step one dimensional at a time using normal distribution of steps within defined parameter bounds.

property stepsize

class seismic.inversion.wavefield_decomp.solvers.HistogramIncremental(bounds, nbins=20)

Bases: object

Class to incrementally accumulate N-dimensional histogram stats at runtime.

property bins

property dims

property histograms

class seismic.inversion.wavefield_decomp.solvers.SolverGlobalMhMcmc

Bases: object

Drop-in custom solver for scipy.optimize.minimize, based on Metrolpolis-Hastings Monte Carlo Markov Chain random walk with burn-in and adaptive acceptance rate, followed by N-dimensional 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)

Minimize objective function and return up to N local minima solutions.

Parameters

• objective – Objective function to minimize

• bounds – Bounds of the parameter space.

• args – Any additional fixed parameters needed to completely specify the objective function.

• x0 – Initial guess. If None, will be selected randomly and uniformly within the parameter bounds.

• T – 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 – Maximum number of minima to return

• burnin – Number of random steps to discard before starting to accumulate statistics.

• maxiter – Maximum number of steps to take (including burnin).

• target_ar – Target acceptance rate of point samples generated by stepping.

• ar_tolerance – Tolerance on the acceptance rate before actively adapting the step size.

• cluster_eps – Point proximity tolerance for DBSCAN clustering, in normalized bounds coordinates.

• rnd_seed – Random seed to force deterministic behaviour

• collect_samples – If not None and integral type, collect collect_samples at regular intervals and return in solution.

• logger – Logger instance for outputting log messages.

Returns

OptimizeResult containing solution(s) and solver data.

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 443-457, 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)

Bases: object

Implements computation of the upwards mean S-wave 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:

1. Load data from a station event dataset. Copy of the data is buffered by this

   class in efficient format for energy flux calculation.

2. Define 1D earth model and mantle half-space material properties (in external code).

3. Call instance with models and receive energy flux results.

grid_search(mantle_props, layer_props, layer_index, H_vals, k_vals, flux_window=- 10, 20, ncpus=- 1)

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 zero-based index into layer_props.

Parameters

• mantle_props (LayerProps) – Mantle bulk properties

• layer_props (list of 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)

Given a stack of layers in layer_props, propagate the surface seismograms down to the base of the bottom layer.

Parameters

layer_props – List of LayerProps through which to propagate the surface seismograms.

Returns

(Vr, Vz) time series per event at the bottom of the stack of layers.

times()

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)

Plot SU energy as function of H-k.

Parameters

• H – Grid of H coordinates corresponding to Esu point values.

• k – Grid of k coordinates corresponding to Esu point values.

• Esu – Energy values on H-k grid.

• title – Plot title [OPTIONAL]

• savefile_name – Output file in which to save plot [OPTIONAL]

• show – If True, display the image and block until it is closed.

• c_range – Custom range of Esu to contour (min, max values)

• decorator – Callback function to customize plot.

Returns

None

seismic.inversion.wavefield_decomp.wfd_plot.plot_Nd(soln, title='', scale=1.0, vars=None)

Plotting routine for N-dimensional solution in grid format. Diagonal contains histograms of solution samples for each variable, and off-diagonal contains pair-wise 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 domain-agnostic means of plotting an N-dimensional solution generated by optimize_minimize_mhmcmc_cluster(). It should not have seismology-specific 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))

Module contents