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)[source]

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

(solution, dict), str

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) – 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)[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 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)[source]

Bases: object

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

notify_accept()[source]
take_step(x)[source]
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
class seismic.inversion.wavefield_decomp.solvers.HistogramIncremental(bounds, nbins=20)[source]

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[source]

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)[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 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)[source]

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.

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

times()[source]

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