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

seismic.inversion.wavefield_decomp.runners module

Batch execution interfaces for wavefield continuation methods and solvers.

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 (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:

(solution, dict), str

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 (numpy.array) – 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 (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)

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)

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)

Save solution to HDF5 file. In general soln_configs will be an ensemble of per-station 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 dict-like) – 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)

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 (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 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.

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)

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

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

times()

Get array of discrete time values of the time axis

Returns:

Array of evenly spaced time values

Return type:

numpy.array

seismic.inversion.wavefield_decomp.wfd_plot module

Module contents